You’d calculate & store dBessel(…)/dorder (generally, dF(x)/dx). This can be done in [autograd.Function’s] forward() or postponed to backward(). Irrespective of how you obtained a gradient, in backward you do multiplication: dLoss/dx = dLoss/dF(x) * dF(x)/dx = grad_out * local_grad  and return that for each function input.
nn.Module is orthogonal to this, it would only be an additional wrapper around a function (or rather, a call to Function.apply). That is mainly useful for adding local nn.Parameters or if you use nn.Sequential declarative pipelines.