# 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.

• Are you sure about your Hamiltonian? There are squares and a factor 1/2 missing in the kinetic part. Jan 14, 2019 at 17:52
• Yes. If velocities are relativistic ($v \to c$) then kinetic energy no longer equals $\frac{p^2}{2m}$. This is the case of $m = 0$. Jan 14, 2019 at 18:41
• @PiotrBenedysiuk This question should be asked on the physics-site. Jan 21, 2019 at 8:35
• Good point. I mirrored it to over there Jan 21, 2019 at 10:02
• The exact claim, as can be found in the works of Hairer et al., is that the Verlet method in any of its implementation preserves to $O(Δt^4)$ an $O(Δt^2)$ modification of the Hamiltonian. This modification includes derivatives of the Hamiltonian, so that close to singular points of the Hamiltonian the $O(Δt^2)$ modification need not be small at all. // Of course this does not help if the Hamiltonian used is not even a first integral of the system. Jan 23, 2019 at 13:58