# Is there something wrong with this numerical simulation of Schwarzschild photon orbits?

This is a follow-up question to the answer given at What is the exact gravitational force between two masses including relativistic effects?. Unfortunately the author hasn't been online for a few years and therefore does not answer comments any more.

In the answer given there the differential equation of motion in Schwarzschild coordinates was

$$\ddot{r} = -\frac{G m}{r^2} + r\dot{\theta}^2 - \frac{ {\color{red} 3} G m}{c^2}\dot{\theta}^2$$

$$\ddot{\theta} = -\frac{2}{r}\dot{r}\dot{\theta}$$

for the angular acceleration. When I plot the path for an object near the speed of light, with this formula I get a stable orbit at $r_0=2 r_s$:

But shouldn't that be at $r_0=1.5 r_s$, the photon sphere? With that formula the orbiting particle would fall into the black hole very quickly, for example, with $v_0=0.999c$ at $r_0=1.6 r_s$:

When I replace the term 3Gm/c² with 2GM/c² so that

$$\ddot{r} = -\frac{G m}{r^2} + r\dot{\theta}^2 - \frac{ {\color{red} 2} G m}{c^2}\dot{\theta}^2$$

I get the expected result with a stable orbit right at the photon sphere (initial velocity again $v_0=0.999c$):

So my question is: is the formula wrong and the factor 3 needs to be replaced with a factor of 2, or are there different minimum-distances for stable orbits, one for particles, and one for photons? Or did I miss something else? Wikipedia says:

The radius of the photon sphere, which is also the lower bound for any stable orbit, is $1.5 r_s$

so I would expect that also particles with mass should stay in orbit if they are close to the speed of light and slightly above the photon sphere.

For reproduction of the problem the Mathematica-code as I believe it to be correct is available (with the factor 2 instead of 3)

There seem to be several confusions here. Massive and massless particles behave qualitatively differently, even if the massive particle is traveling very fast.

• The minimum radius for a stable orbit for a massive particle is $3 r_s$. Circular orbits above this radius are all stable.
• Massless particles only have circular orbits at the photon sphere, $(3/2) r_s$. These orbits are not stable, Wikipedia is wrong. Massless particles also obey a different equation of motion.

The other confusion is that what your simulations are showing has nothing to do with stability. Your particles are falling into the center because you're not giving them the right initial velocity. It's just like classical mechanics: if you suddenly took away half the Earth's speed, it would start falling inward. In order to initialize them to the right initial velocity, you need to solve for $\dot{\theta}$ so that $\ddot{r} = 0$.

This is in contrast to the massless case, where the initial velocity is already determined for you (i.e. it's the speed of light).

• So the factor of 3 in $\ddot{r} = -\frac{G m}{r^2} + r\dot{\theta}^2 - ({\color{red} 3} G m/c^2)\dot{\theta}^2$ is correct for massive particles? – Симон Тыран May 25 '16 at 1:54
• @СимонТыран Yes. – knzhou May 25 '16 at 2:01
• Thanks and +1. Yet I still don't understand why I get my stable orbit for particles at r0=2rs with v0=0.999c and not 3rs, since the initial angular velocity I give is v0/r0 and this should be proper length and time and as I understand it therefore also proper velocity, but I hope I'll figure that out too... – Симон Тыран May 25 '16 at 2:09
• Very good except for one thing; we test stability by perturbing v and seeing what happens. – Joshua May 25 '16 at 15:34
• Your answer really helped, but are you sure that "The minimum radius for a stable orbit for a massive particle is 3 rs"? When I compensate for all contractions it is still 1.5 rs like Wikipedia says. Do you have any quote for the 3 rs, or is my calculation (see below) correct now? – Симон Тыран May 26 '16 at 4:42

Thanks to the hint given by knzhou I figured out that if one wants to give the particle a proper initial velocity of $v_0$, the initial velocity in terms of Schwarzschild coordinates $v_i$ would then be

$$\dot{\theta}(0)\cdot r(0) = \frac{{v\perp}_0}{ \color{green}{\sqrt{ 1-v_0^2/c^2}}}$$

for the transversal component, and

$$\dot{r}(0) = \frac{{v\parallel}_0\cdot \color{blue}{\sqrt{1-r_s/r_0}}}{ \color{green}{\sqrt{ 1-v_0^2/c^2}}}$$

for the radial component since one has to compensate for the gravitational length contraction (blue) and the length contraction due to the particle's velocity (green) with respect to the particle's proper time.

Now I get the expected results: a circular orbit with transversal initial velocity around the photon sphere, and a stationary particle with outwards initial velocity around the event horizon when v0 is set close to c.