# Stability criterion for leapfrog in relativistic physics

I am doing a 2D MD simulations of charge carriers in graphene using the Leapfrog algorithm. I am trying to prove that, in some specific cases (when distance between particles is small), the method is unstable.

The Hamiltonian is given as $$H = |\vec{p_{1}}| + |\vec{p_2}| - \frac{\alpha}{r_{12}}$$. With $$p$$ the momentums, $$\alpha$$ some numerical constant and $$r_{12}$$ some constants.

The leapfrog iteration would be:

$$\vec{p_i} \to \vec{p_i} - 0.5 \Delta t \frac{\alpha}{r_{ij}^2} \hat{r_{ij}}\\ \vec{x_i} \to \vec{x_i} + \Delta t \hat{p_i} \\ \vec{p_i} \to \vec{p_i} - 0.5 \Delta t \frac{\alpha}{r_{ij}^2} \hat{r_{ij}}$$

where the vectors with hat's are the direction vectors. I expect the stability to break down when $$r_{ij}$$ gets small. Yet, I tried using analysis from this (page 7 onward) source and it didn't really help.

Also, a different way to derive leapfrog is to define coordinate shift operators $$L_p$$ and $$L_r$$ that move $$p \to p +\Delta p$$ and $$x \to x + \Delta x$$. Then the differential equation can be written as: $$f(p(t), r(t)) = \exp\left(t (L_p + L_r) \right) f(p(0), r(0)).$$ Using Trotter's identity we can write: $$\exp\left( \Delta t (A + B) \right) \approx \exp\left( 0.5\Delta t A \right) \exp\left( \Delta t B \right) \exp\left(0.5 \Delta t A \right)$$ and this is exactly the Leapfrog algorithm. I tried figuring out if things breakdown in 1D case of two particles going straight at each other. To do this I worked out the commutator term $$[L_p,L_r]$$ that comes up in Zassenhaus formula. This also didn't help.

Anyone has a fresh look or has some derivations that go beyond basic 1D cases where things go nicely? Any help is appreciated.

• For completness sake: I have an ongoing bounty for this exact question on the math exchange. I don't expect much contributions there - it has been ignored by the community there so far. Jan 21, 2019 at 11:12