Using the Hamiltonian for a test particle in Kerr spacetime, we arrive at the following equations for generalized position and momenta (in natural units, $G = c = M = 1$):

\begin{align} \dot{r} &= \frac{\Delta}{\Sigma}p_r\\ \dot{\theta} &= \frac{1}{\Sigma}p_\theta\\ \dot{\phi} &= \frac{1}{\Sigma\Delta}\left(2arE + (\Sigma - 2r)L\csc^2\theta\right)\\ \dot{p}_r &= \frac{1}{\Sigma\Delta}\left[(r - 1)((r^2 + a^2)\mu - \kappa + r\Delta\mu + 2r(r^2 + a^2)E^2 - 2aEL\right]-\frac{2p_r^2(r-1)}{\Sigma}\\ \dot{p}_\theta &= \frac{\sin\theta\cos\theta}{\Sigma}\left[L^2\csc^4\theta - a^2(E^2 + \mu)\right]\\ \dot{p}_\phi &= 0 \end{align} where, $\kappa = p_\theta^2 + L^2\csc^2\theta + a^2 - (E^2\sin^2\theta + \mu)$, $\Sigma = r^2 + a^2\cos^2\theta$, $\Delta = r^2 - 2r + a^2$; and $\mu$, $E = -p_t$ and $L = p_\phi$ are the test particle mass, energy and angular momentum, respectively. Also, $\mu = -1$ for massive particles, while for massless particles, $\mu = 0$, and the derivatives are with respect to an affine parameter, say $\lambda$.

I want to implement a symplectic leapfrog/verlet solver to solve this system. But I don't understand, how this set of equations should be discretized. Since $p_i$ are the momenta and not velocities, does the usual Taylor series approach make sense here? If not, could someone help me understand, how a verlet/leapfrog scheme can be used here?

  • 1
    $\begingroup$ Is there a reason you can't just plug these equations into a symplectic integrator. DifferentialEquations.jl contains quite a few that are highly performant and easy to use diffeq.sciml.ai/stable/solvers/dynamical_solve/… $\endgroup$
    – whpowell96
    Aug 25, 2020 at 7:12
  • $\begingroup$ I want to understand the discretization process for this system, as I am looking to solve this in Python. I have tried HamiltonianProblem and DynamicalProblem in Julia, and they work beautifully, but they use Forward Mode Automatic Differentiation to obtain the equations of motion, which abstracts the discretization. @whpowell96 $\endgroup$
    – P_0
    Aug 25, 2020 at 7:20
  • $\begingroup$ @P_0 If instead of this relatively large system you had dynamic equations for 1D simple harmonic oscillator in the Hamiltonian form, would you know how to discretize those for symplectic time integration? It is structurally similar but easier to discuss. $\endgroup$ Aug 27, 2020 at 2:13
  • $\begingroup$ @MaximUmansky The dynamical equations for 1D SHM are $\dot{p} = -k q$ and $\dot{q} = \frac{p}{M}$, where the latter is just "rescaled" velocity or contravariant momentum. In my case, I have covariant momenta in all my equations. but that can be thought of as an index-lowering operation performed on contravariant momenta, and so, the equations in the question already account for that. Given all this, should this leapfrog scheme not work? $\endgroup$
    – P_0
    Sep 3, 2020 at 10:29
  • $\begingroup$ @MaximUmansky I ask, because in Python, some of the momentum derivatives blow up in the first step itself. Should I share my code in the question body? $\endgroup$
    – P_0
    Sep 3, 2020 at 10:32


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.