General relativity
${\displaystyle G_{\mu \nu }+\Lambda g_{\mu \nu }={\kappa }T_{\mu \nu }}$
Introduction History Timeline Tests Mathematical formulation
Fundamental concepts Equivalence principle Special relativity World line Pseudo-Riemannian manifold
Phenomena Kepler problem Gravitational lensing Gravitational waves Frame-dragging Geodetic effect Event horizon Singularity Black hole Spacetime Spacetime diagrams Minkowski spacetime Einstein–Rosen bridge
Equations Formalisms Equations Linearized gravity Einstein field equations Friedmann Geodesics Mathisson–Papapetrou–Dixon Hamilton–Jacobi–Einstein Formalisms ADM BSSN Post-Newtonian Advanced theory Kaluza–Klein theory Quantum gravity
Solutions Schwarzschild (interior) Reissner–Nordström Gödel Kerr Kerr–Newman Kerr–Newman–de Sitter Kasner Lemaître–Tolman Taub–NUT Milne Robertson–Walker Oppenheimer–Snyder pp-wave van Stockum dust Weyl−Lewis−Papapetrou
Scientists Einstein Lorentz Hilbert Poincaré Schwarzschild de Sitter Reissner Nordström Weyl Eddington Friedman Milne Zwicky Lemaître Oppenheimer Gödel Wheeler Robertson Bardeen Walker Kerr Chandrasekhar Ehlers Penrose Hawking Raychaudhuri Taylor Hulse van Stockum Taub Newman Yau Thorne others
Physics portal  Category

In general relativity, Schwarzschild geodesics describe the motion of test particles in the gravitational field of a central fixed mass ${\textstyle M,}$ that is, motion in the Schwarzschild metric. Schwarzschild geodesics have been pivotal in the validation of Einstein's theory of general relativity. For example, they provide accurate predictions of the anomalous precession of the planets in the Solar System and of the deflection of light by gravity.

Schwarzschild geodesics pertain only to the motion of particles of masses so small they contribute little to the gravitational field. However, they are highly accurate in many astrophysical scenarios provided that ${\textstyle m}$ is many-fold smaller than the central mass ${\textstyle M}$, e.g., for planets orbiting their star. Schwarzschild geodesics are also a good approximation to the relative motion of two bodies of arbitrary mass, provided that the Schwarzschild mass ${\textstyle M}$ is set equal to the sum of the two individual masses ${\textstyle m_{1}}$ and ${\textstyle m_{2}}$. This is important in predicting the motion of binary stars in general relativity.

Historical context

The Schwarzschild metric is named in honour of its discoverer Karl Schwarzschild, who found the solution in 1915, only about a month after the publication of Einstein's theory of general relativity. It was the first exact solution of the Einstein field equations other than the trivial flat space solution.

In 1931, Yusuke Hagihara published a paper showing that the trajectory of a test particle in the Schwarzschild metric can be expressed in terms of elliptic functions.^[1]

Samuil Kaplan in 1949 has shown that there is a minimum radius for the circular orbit to be stable in Schwarzschild metric.^[2]

Schwarzschild metric

An exact solution to the Einstein field equations is the Schwarzschild metric, which corresponds to the external gravitational field of an uncharged, non-rotating, spherically symmetric body of mass ${\textstyle M}$. The Schwarzschild solution can be written as^[3]

${\displaystyle c^{2}{d\tau }^{2}=\left(1-{\frac {r_{\text{s}}}{r}}\right)c^{2}dt^{2}-{\frac {dr^{2}}{1-{\frac {r_{\text{s}}}{r}}}}-r^{2}d\theta ^{2}-r^{2}\sin ^{2}\theta \,d\varphi ^{2}}$

where

${\displaystyle \tau }$, in the case of a test particle of small positive mass, is the proper time (time measured by a clock moving with the particle) in seconds,

${\displaystyle c}$ is the speed of light in meters per second,

${\displaystyle t}$ is, for ${\displaystyle r>r_{\text{s}}}$, the time coordinate (time measured by a stationary clock at infinity) in seconds,

${\displaystyle r}$ is, for ${\displaystyle r>r_{\text{s}}}$, the radial coordinate (circumference of a circle centered at the star divided by ${\displaystyle 2\pi }$) in meters,

${\displaystyle \theta }$ is the colatitude (angle from North) in radians,

${\displaystyle \varphi }$ is the longitude in radians, and

${\displaystyle r_{\text{s}}}$ is the Schwarzschild radius of the massive body (in meters), which is related to its mass ${\textstyle M}$ by

${\displaystyle r_{\text{s}}={\frac {2GM}{c^{2}}},}$

where ${\textstyle G}$ is the gravitational constant. The classical Newtonian theory of gravity is recovered in the limit as the ratio ${\textstyle {\frac {r_{\text{s}}}{r}}}$ goes to zero. In that limit, the metric returns to that defined by special relativity.

In practice, this ratio is almost always extremely small. For example, the Schwarzschild radius ${\textstyle r_{\text{s}}}$ of the Earth is roughly 9 mm (3⁄8 inch); at the surface of the Earth, the corrections to Newtonian gravity are only one part in a billion. The Schwarzschild radius of the Sun is much larger, roughly 2953 meters, but at its surface, the ratio ${\textstyle {\frac {r_{\text{s}}}{r}}}$ is roughly 4 parts in a million. A white dwarf star is much denser, but even here the ratio at its surface is roughly 250 parts in a million. The ratio only becomes large close to ultra-dense objects such as neutron stars (where the ratio is roughly 50%) and black holes.

Orbits of test particles

Comparison between the orbit of a test particle in Newtonian (left) and Schwarzschild (right) spacetime; note the apsidal precession on the right.

We may simplify the problem by using symmetry to eliminate one variable from consideration. Since the Schwarzschild metric is symmetrical about ${\textstyle \theta ={\frac {\pi }{2}}}$, any geodesic that begins moving in that plane will remain in that plane indefinitely (the plane is totally geodesic). Therefore, we orient the coordinate system so that the orbit of the particle lies in that plane, and fix the ${\textstyle \theta }$ coordinate to be ${\textstyle {\frac {\pi }{2}}}$ so that the metric (of this plane) simplifies to

${\displaystyle c^{2}d\tau ^{2}=\left(1-{\frac {r_{\text{s}}}{r}}\right)c^{2}dt^{2}-{\frac {dr^{2}}{1-{\frac {r_{\text{s}}}{r}}}}-r^{2}d\varphi ^{2}.}$

Two constants of motion (values that do not change over proper time ${\displaystyle \tau }$) can be identified (cf. the derivation given below). One is the total energy ${\textstyle E}$:

${\displaystyle \left(1-{\frac {r_{\text{s}}}{r}}\right){\frac {dt}{d\tau }}={\frac {E}{mc^{2}}}.}$

and the other is the specific angular momentum:

${\displaystyle h={\frac {L}{\mu }}=r^{2}{\frac {d\varphi }{d\tau }},}$

where ${\textstyle L}$ is the total angular momentum of the two bodies, and ${\textstyle \mu }$ is the reduced mass. When ${\textstyle M\gg m}$, the reduced mass is approximately equal to ${\textstyle m}$. Sometimes it is assumed that ${\textstyle m=\mu }$. In the case of the planet Mercury this simplification introduces an error more than twice as large as the relativistic effect. When discussing geodesics, ${\textstyle m}$ can be considered fictitious, and what matters are the constants ${\textstyle {\frac {E}{m}}}$ and ${\textstyle h}$. In order to cover all possible geodesics, we need to consider cases in which ${\textstyle {\frac {E}{m}}}$ is infinite (giving trajectories of photons) or imaginary (for tachyonic geodesics). For the photonic case, we also need to specify a number corresponding to the ratio of the two constants, namely ${\textstyle {\frac {mh}{E}}}$, which may be zero or a non-zero real number.

Substituting these constants into the definition of the Schwarzschild metric

${\displaystyle c^{2}=\left(1-{\frac {r_{\text{s}}}{r}}\right)c^{2}\left({\frac {dt}{d\tau }}\right)^{2}-{\frac {1}{1-{\frac {r_{\text{s}}}{r}}}}\left({\frac {dr}{d\tau }}\right)^{2}-r^{2}\left({\frac {d\varphi }{d\tau }}\right)^{2},}$

yields an equation of motion for the radius as a function of the proper time ${\textstyle \tau }$:

${\displaystyle \left({\frac {dr}{d\tau }}\right)^{2}={\frac {E^{2}}{m^{2}c^{2}}}-\left(1-{\frac {r_{\text{s}}}{r}}\right)\left(c^{2}+{\frac {h^{2}}{r^{2}}}\right).}$

The formal solution to this is

${\displaystyle \tau =\int {\frac {dr}{\pm {\sqrt {{\frac {E^{2}}{m^{2}c^{2}}}-\left(1-{\frac {r_{\text{s}}}{r}}\right)\left(c^{2}+{\frac {h^{2}}{r^{2}}}\right)}}}}.}$

Note that the square root will be imaginary for tachyonic geodesics.

Using the relation higher up between ${\textstyle {\frac {dt}{d\tau }}}$ and ${\textstyle E}$, we can also write

${\displaystyle t=\int {\frac {dr}{\pm c\left(1-{\frac {r_{\text{s}}}{r}}\right){\sqrt {1-\left(1-{\frac {r_{\text{s}}}{r}}\right)\left(c^{2}+{\frac {h^{2}}{r^{2}}}\right){\frac {m^{2}c^{2}}{E^{2}}}}}}}.}$

Since asymptotically the integrand is inversely proportional to ${\textstyle r-r_{\text{s}}}$, this shows that in the ${\textstyle r,\theta ,\varphi ,t}$ frame of reference if ${\textstyle r}$ approaches ${\textstyle r_{\text{s}}}$ it does so exponentially without ever reaching it. However, as a function of ${\textstyle \tau }$, ${\textstyle r}$ does reach ${\textstyle r_{\text{s}}}$.

The above solutions are valid while the integrand is finite, but a total solution may involve two or an infinity of pieces, each described by the integral but with alternating signs for the square root.

When ${\textstyle E=mc^{2}}$ and ${\textstyle h=0}$, we can solve for ${\textstyle t}$ and ${\textstyle \tau }$ explicitly:

${\displaystyle {\begin{aligned}t&={\text{constant}}\pm {\frac {r_{\text{s}}}{c}}\left({\frac {2}{3}}\left({\frac {r}{r_{\text{s}}}}\right)^{\frac {3}{2}}+2{\sqrt {\frac {r}{r_{\text{s}}}}}+\ln {\frac {\left|{\sqrt {\frac {r}{r_{\text{s}}}}}-1\right|}{{\sqrt {\frac {r}{r_{\text{s}}}}}+1}}\right)\\\tau &={\text{constant}}\pm {\frac {2}{3}}{\frac {r_{\text{s}}}{c}}\left({\frac {r}{r_{\text{s}}}}\right)^{\frac {3}{2}}\end{aligned}}}$

and for photonic geodesics (${\textstyle m=0}$) with zero angular momentum

${\displaystyle {\begin{aligned}t&={\text{constant}}\pm {\frac {1}{c}}\left(r+r_{\text{s}}\ln \left|{\frac {r}{r_{\text{s}}}}-1\right|\right)\\\tau &={\text{constant}}.\end{aligned}}}$

(Although the proper time is trivial in the photonic case, one can define an affine parameter ${\textstyle \lambda }$, and then the solution to the geodesic equation is ${\textstyle r=c_{1}\lambda +c_{2}}$.)

Another solvable case is that in which ${\textstyle E=0}$ and ${\textstyle t}$ and ${\textstyle \varphi }$ are constant. In the volume where ${\textstyle r<r_{\text{s}}}$ this gives for the proper time

${\displaystyle \tau ={\text{constant}}\pm {\frac {r_{\text{s}}}{c}}\left(\arcsin {\sqrt {\frac {r}{r_{\text{s}}}}}-{\sqrt {{\frac {r}{r_{\text{s}}}}\left(1-{\frac {r}{r_{\text{s}}}}\right)}}\right).}$

This is close to solutions with ${\textstyle {\frac {E^{2}}{m^{2}}}}$ small and positive. Outside of ${\textstyle r_{\text{s}}}$ the ${\textstyle E=0}$ solution is tachyonic and the "proper time" is space-like:

${\displaystyle \tau ={\text{constant}}\pm i{\frac {r_{\text{s}}}{c}}\left(\ln \left({\sqrt {\frac {r}{r_{\text{s}}}}}+{\sqrt {{\frac {r}{r_{\text{s}}}}-1}}\right)+{\sqrt {{\frac {r}{r_{\text{s}}}}\left({\frac {r}{r_{\text{s}}}}-1\right)}}\right).}$

This is close to other tachyonic solutions with ${\textstyle {\frac {E^{2}}{m^{2}}}}$ small and negative. The constant ${\textstyle t}$ tachyonic geodesic outside ${\textstyle r_{\text{s}}}$ is not continued by a constant ${\textstyle t}$ geodesic inside ${\textstyle r_{\text{s}}}$, but rather continues into a "parallel exterior region" (see Kruskal–Szekeres coordinates). Other tachyonic solutions can enter a black hole and re-exit into the parallel exterior region. The constant ${\textstyle t}$ solution inside the event horizon (${\textstyle r_{\text{s}}}$) is continued by a constant ${\textstyle t}$ solution in a white hole.

When the angular momentum is not zero we can replace the dependence on proper time by a dependence on the angle ${\textstyle \varphi }$ using the definition of ${\textstyle h}$

${\displaystyle \left({\frac {dr}{d\varphi }}\right)^{2}=\left({\frac {dr}{d\tau }}\right)^{2}\left({\frac {d\tau }{d\varphi }}\right)^{2}=\left({\frac {dr}{d\tau }}\right)^{2}\left({\frac {r^{2}}{h}}\right)^{2},}$

which yields the equation for the orbit

${\displaystyle \left({\frac {dr}{d\varphi }}\right)^{2}={\frac {r^{4}}{b^{2}}}-\left(1-{\frac {r_{\text{s}}}{r}}\right)\left({\frac {r^{4}}{a^{2}}}+r^{2}\right)}$

where, for brevity, two length-scales, ${\textstyle a}$ and ${\textstyle b}$, have been defined by

${\displaystyle {\begin{aligned}a&={\frac {h}{c}},\\b&={\frac {cL}{E}}={\frac {hmc}{E}}.\end{aligned}}}$

Note that in the tachyonic case, ${\textstyle a}$ will be imaginary and ${\textstyle b}$ real or infinite.

The same equation can also be derived using a Lagrangian approach^[4] or the Hamilton–Jacobi equation^[5] (see below). The solution of the orbit equation is

${\displaystyle \varphi =\int {\frac {dr}{\pm r^{2}{\sqrt {{\frac {1}{b^{2}}}-\left(1-{\frac {r_{\text{s}}}{r}}\right)\left({\frac {1}{a^{2}}}+{\frac {1}{r^{2}}}\right)}}}}.}$

This can be expressed in terms of the Weierstrass elliptic function ${\textstyle \wp }$.^[6]

Local and delayed velocities

Unlike in classical mechanics, in Schwarzschild coordinates ${\textstyle {\frac {{\rm {d}}r}{{\rm {d}}\tau }}}$ and ${\textstyle r\ {\frac {{\rm {d}}\varphi }{{\rm {d}}\tau }}}$ are not the radial ${\textstyle v_{\parallel }}$ and transverse ${\textstyle v_{\perp }}$ components of the local velocity ${\textstyle v}$ (relative to a stationary observer), instead they give the components for the celerity which are related to ${\textstyle v}$ by

${\displaystyle {\frac {{\rm {d}}r}{{\rm {d}}\tau }}=v_{\parallel }{\sqrt {1-{\frac {r_{\text{s}}}{r}}}}\ \gamma }$

for the radial and

${\displaystyle {\frac {{\rm {d}}\varphi }{{\rm {d}}\tau }}={\frac {v_{\perp }}{r}}\ \gamma }$

for the transverse component of motion, with ${\textstyle v^{2}=v_{\parallel }^{2}+v_{\perp }^{2}}$. The coordinate bookkeeper far away from the scene observes the shapiro-delayed velocity ${\textstyle {\hat {v}}}$, which is given by the relation

${\displaystyle {\hat {v}}_{\perp }=v_{\perp }{\sqrt {1-{\frac {r_{\text{s}}}{r}}}}}$ and ${\displaystyle {\hat {v}}_{\parallel }=v_{\parallel }\left(1-{\frac {r_{\text{s}}}{r}}\right)}$.

The time dilation factor between the bookkeeper and the moving test-particle can also be put into the form

${\displaystyle {\frac {{\rm {d}}\tau }{{\rm {d}}t}}={\frac {\sqrt {1-{\frac {r_{\text{s}}}{r}}}}{\gamma }}}$

where the numerator is the gravitational, and the denominator is the kinematic component of the time dilation. For a particle falling in from infinity the left factor equals the right factor, since the in-falling velocity ${\textstyle v}$ matches the escape velocity ${\textstyle c{\sqrt {\frac {r_{\text{s}}}{r}}}}$ in this case.

The two constants angular momentum ${\textstyle L}$ and total energy ${\textstyle E}$ of a test-particle with mass ${\textstyle m}$ are in terms of ${\textstyle v}$

${\displaystyle L=m\ v_{\perp }\ r\ \gamma }$

and

${\displaystyle E=mc^{2}\ {\sqrt {1-{\frac {r_{\text{s}}}{r}}}}\ \gamma }$

where

${\displaystyle E=E_{\rm {rest}}+E_{\rm {kin}}+E_{\rm {pot}}}$

and

${\displaystyle E_{\rm {rest}}=mc^{2}\ ,\ \ E_{\rm {kin}}=(\gamma -1)mc^{2}\ ,\ \ E_{\rm {pot}}=\left({\sqrt {1-{\frac {r_{\text{s}}}{r}}}}-1\right)\ \gamma \ mc^{2}}$

For massive testparticles ${\textstyle \gamma }$ is the Lorentz factor ${\textstyle \gamma =1/{\sqrt {1-v^{2}/c^{2}}}}$ and ${\textstyle \tau }$ is the proper time, while for massless particles like photons ${\textstyle \gamma }$ is set to ${\textstyle 1}$ and ${\textstyle \tau }$ takes the role of an affine parameter. If the particle is massless ${\textstyle E_{\rm {rest}}}$ is replaced with ${\textstyle E_{\rm {kin}}}$ and ${\textstyle mc^{2}}$ with ${\textstyle hf}$, where ${\textstyle h}$ is the Planck constant and ${\textstyle f}$ the locally observed frequency.

Exact solution using elliptic functions

The fundamental equation of the orbit is easier to solve^{[note 1]} if it is expressed in terms of the inverse radius ${\textstyle u={\frac {1}{r}}}$

${\displaystyle \left({\frac {du}{d\varphi }}\right)^{2}={\frac {1}{b^{2}}}-\left(1-ur_{\text{s}}\right)\left({\frac {1}{a^{2}}}+u^{2}\right)}$

The right-hand side of this equation is a cubic polynomial, which has three roots, denoted here as ${\textstyle u_{1}}$, ${\textstyle u_{2}}$, and ${\textstyle u_{3}}$

${\displaystyle \left({\frac {du}{d\varphi }}\right)^{2}=r_{\text{s}}\left(u-u_{1}\right)\left(u-u_{2}\right)\left(u-u_{3}\right)}$

The sum of the three roots equals the coefficient of the ${\textstyle u^{2}}$ term

${\displaystyle u_{1}+u_{2}+u_{3}={\frac {1}{r_{\text{s}}}}}$

A cubic polynomial with real coefficients can either have three real roots, or one real root and two complex conjugate roots. If all three roots are real numbers, the roots are labeled so that ${\textstyle u_{1}<u_{2}<u_{3}}$. If instead there is only one real root, then that is denoted as ${\textstyle u_{3}}$; the complex conjugate roots are labeled ${\textstyle u_{1}}$ and ${\textstyle u_{2}}$. Using Descartes' rule of signs, there can be at most one negative root; ${\textstyle u_{1}}$ is negative if and only if ${\textstyle b<a}$. As discussed below, the roots are useful in determining the types of possible orbits.

Given this labeling of the roots, the solution of the fundamental orbital equation is

${\displaystyle u=u_{1}+\left(u_{2}-u_{1}\right)\,\mathrm {sn} ^{2}\left({\frac {1}{2}}\varphi {\sqrt {r_{\text{s}}\left(u_{3}-u_{1}\right)}}+\delta \right)}$

where ${\textstyle \mathrm {sn} }$ represents the sinus amplitudinus function (one of the Jacobi elliptic functions) and ${\textstyle \delta }$ is a constant of integration reflecting the initial position. The elliptic modulus ${\textstyle k}$ of this elliptic function is given by the formula

${\displaystyle k={\sqrt {\frac {u_{2}-u_{1}}{u_{3}-u_{1}}}}}$

Precession of orbits

The function sn and its square sn² have periods of 4K and 2K, respectively, where K is defined by the equation^{[note 2]}

${\displaystyle K=\int _{0}^{1}{\frac {dy}{\sqrt {\left(1-y^{2}\right)\left(1-k^{2}y^{2}\right)}}}}$

Therefore, the change in φ over one oscillation of ${\textstyle u}$ (or, equivalently, one oscillation of ${\textstyle r}$) equals^[7]

${\displaystyle \Delta \varphi ={\frac {4K}{\sqrt {r_{\text{s}}\left(u_{3}-u_{1}\right)}}}}$

In the classical limit, u₃ approaches ${\textstyle {\frac {1}{r_{\text{s}}}}}$ and is much larger than ${\textstyle u_{1}}$ or ${\textstyle u_{2}}$. Hence, ${\textstyle k^{2}}$ is approximately

${\displaystyle k^{2}={\frac {u_{2}-u_{1}}{u_{3}-u_{1}}}\approx r_{\text{s}}\left(u_{2}-u_{1}\right)\ll 1}$

For the same reasons, the denominator of Δφ is approximately

${\displaystyle {\frac {1}{\sqrt {r_{\text{s}}\left(u_{3}-u_{1}\right)}}}={\frac {1}{\sqrt {1-r_{\text{s}}\left(2u_{1}+u_{2}\right)}}}\approx 1+{\frac {1}{2}}r_{\text{s}}\left(2u_{1}+u_{2}\right)}$

Since the modulus ${\textstyle k}$ is close to zero, the period K can be expanded in powers of ${\textstyle k}$; to lowest order, this expansion yields

${\displaystyle K\approx \int _{0}^{1}{\frac {dy}{\sqrt {1-y^{2}}}}\left(1+{\frac {1}{2}}k^{2}y^{2}\right)={\frac {\pi }{2}}\left(1+{\frac {k^{2}}{4}}\right)}$

Substituting these approximations into the formula for Δφ yields a formula for angular advance per radial oscillation

${\displaystyle \delta \varphi =\Delta \varphi -2\pi \approx {\frac {3}{2}}\pi r_{\text{s}}\left(u_{1}+u_{2}\right)}$

For an elliptical orbit, ${\textstyle u_{1}}$ and ${\textstyle u_{2}}$ represent the inverses of the longest and shortest distances, respectively. These can be expressed in terms of the ellipse's semi-major axis ${\textstyle A}$ and its orbital eccentricity ${\textstyle e}$,

${\displaystyle {\begin{aligned}r_{\text{max}}&={\frac {1}{u_{1}}}=A(1+e)\\r_{\text{min}}&={\frac {1}{u_{2}}}=A(1-e)\end{aligned}}}$

giving

${\displaystyle u_{1}+u_{2}={\frac {2}{A\left(1-e^{2}\right)}}}$

Substituting the definition of ${\textstyle r_{\text{s}}}$ gives the final equation

${\displaystyle \delta \varphi \approx {\frac {6\pi GM}{c^{2}A\left(1-e^{2}\right)}}}$

Bending of light by gravity

Diagram of gravitational lensing by a compact body

In the limit as the particle mass m goes to zero (or, equivalently if the light is heading directly toward the central mass, as the length-scale a goes to infinity), the equation for the orbit becomes

${\displaystyle \varphi =\int {\frac {dr}{r^{2}{\sqrt {{\frac {1}{b^{2}}}-\left(1-{\frac {r_{\text{s}}}{r}}\right){\frac {1}{r^{2}}}}}}}}$

Expanding in powers of ${\textstyle {\frac {r_{\text{s}}}{r}}}$, the leading order term in this formula gives the approximate angular deflection δφ for a massless particle coming in from infinity and going back out to infinity:

${\displaystyle \delta \varphi \approx {\frac {2r_{\text{s}}}{b}}={\frac {4GM}{c^{2}b}}.}$

Here, ${\textstyle b}$ is the impact parameter, somewhat greater than the distance of closest approach, ${\textstyle r_{3}}$:^[8]

${\displaystyle b=r_{3}{\sqrt {\frac {r_{3}}{r_{3}-r_{\text{s}}}}}}$

Although this formula is approximate, it is accurate for most measurements of gravitational lensing, due to the smallness of the ratio ${\textstyle {\frac {r_{\text{s}}}{r}}}$. For light grazing the surface of the sun, the approximate angular deflection is roughly 1.75 arcseconds, roughly one millionth part of a circle.

More generally, the geodesics of a photon emitted from a light source located at a radial coordinate ${\displaystyle r={1 \over u}\in [r_{s},\infty [}$ can be calculated as follows, by applying the equation

${\displaystyle {\left({du \over d\varphi }\right)^{2}}=r_{s}\ u^{3}-u^{2}+{\frac {1}{b^{2}}}}$

Geodesic of a photon emitted from a light source located on the event horizon of a black hole and back to it, with an impact parameter ${\displaystyle b>b_{crit}={\frac {3{\sqrt {3}}}{2}}r_{s}}$.

Geodesic of a photon emitted from a light source located on the event horizon of a black hole, with an impact parameter ${\displaystyle b=b_{crit}={\frac {3{\sqrt {3}}}{2}}r_{s}}$ and then moving to the unstable orbit ${\displaystyle r_{crit}={3 \over 2}r_{s}}$. If ${\displaystyle b<b_{crit}}$ the photon is released towards infinity.

The equation can be derived as

${\displaystyle 2\ {\frac {du}{d\phi }}{\frac {d^{2}u}{d\phi ^{2}}}=3\ r_{s}u^{2}{\frac {du}{d\phi }}-2\ u{\frac {du}{d\phi }}}$

which leads to

${\displaystyle {d^{2}u \over d\varphi ^{2}}={\frac {3}{2}}\ r_{s}\ u^{2}-\ u}$

This equation with second derivative can be numerically integrated as follows by a 4^th order Runge-Kutta method, considering a step size ${\displaystyle \Delta \varphi }$ and with:

${\displaystyle k_{1}={d^{2}u \over d\varphi ^{2}}(u)}$,

${\displaystyle k_{2}={d^{2}u \over d\varphi ^{2}}{\bigl (}u+{\frac {\Delta \varphi }{2}}{du \over d\varphi }{\bigr )}}$,

${\displaystyle k_{3}={d^{2}u \over d\varphi ^{2}}{\Bigl (}u+{\frac {\Delta \varphi }{2}}{du \over d\varphi }+{\frac {\Delta \varphi ^{2}}{4}}k_{1}{\Bigr )}}$ and

${\displaystyle k_{4}={d^{2}u \over d\varphi ^{2}}{\Bigl (}u+\Delta \varphi {du \over d\varphi }+{\frac {\Delta \varphi ^{2}}{2}}k_{2}{\Bigr )}}$.

The value at the next step ${\displaystyle {du \over d\varphi }(\varphi +\Delta \varphi )}$ is

${\displaystyle {du \over d\varphi }(\varphi )+{\frac {\Delta \varphi }{6}}(k_{1}+2k_{2}+2k_{3}+k_{4})}$

and the value at the next step ${\displaystyle u(\varphi +\Delta \varphi )}$ is

${\displaystyle u(\varphi )+\Delta \varphi {du \over d\varphi }(\varphi )+{\frac {\Delta \varphi ^{2}}{6}}(k_{1}+k_{2}+k_{3})}$

The step ${\displaystyle \Delta \varphi }$ can be chosen to be constant or adaptive, depending on the accuracy required on ${\displaystyle r={1 \over u}}$.

Relation to Newtonian physics

Circular orbits and their stability

Effective radial potential for various angular momenta. At small radii, the energy drops precipitously, causing the particle to be pulled inexorably inwards to ${\textstyle r=0}$. However, when the normalized angular momentum ${\textstyle {\frac {a}{r_{\text{s}}}}={\frac {L}{mcr_{\text{s}}}}}$ equals the square root of three, a metastable circular orbit is possible at the radius highlighted with a green circle. At higher angular momenta, there is a significant centrifugal barrier (orange curve) and an unstable inner radius, highlighted in red.

The effective potential ${\textstyle V}$ can be re-written in terms of the length ${\textstyle a={\frac {h}{c}}}$.

${\displaystyle V(r)={\frac {\mu c^{2}}{2}}\left[-{\frac {r_{\text{s}}}{r}}+{\frac {a^{2}}{r^{2}}}-{\frac {r_{\text{s}}a^{2}}{r^{3}}}\right]}$

Circular orbits are possible when the effective force is zero

${\displaystyle F=-{\frac {dV}{dr}}=-{\frac {\mu c^{2}}{2r^{4}}}\left[r_{\text{s}}r^{2}-2a^{2}r+3r_{\text{s}}a^{2}\right]=0}$

i.e., when the two attractive forces — Newtonian gravity (first term) and the attraction unique to general relativity (third term) — are exactly balanced by the repulsive centrifugal force (second term). There are two radii at which this balancing can occur, denoted here as r_inner and r_outer

${\displaystyle {\begin{aligned}r_{\text{outer}}&={\frac {a^{2}}{r_{\text{s}}}}\left(1+{\sqrt {1-{\frac {3r_{\text{s}}^{2}}{a^{2}}}}}\right)\\[3pt]r_{\text{inner}}&={\frac {a^{2}}{r_{\text{s}}}}\left(1-{\sqrt {1-{\frac {3r_{\text{s}}^{2}}{a^{2}}}}}\right)={\frac {3a^{2}}{r_{\text{outer}}}}\end{aligned}}}$

which are obtained using the quadratic formula. The inner radius r_inner is unstable, because the attractive third force strengthens much faster than the other two forces when r becomes small; if the particle slips slightly inwards from r_inner (where all three forces are in balance), the third force dominates the other two and draws the particle inexorably inwards to r = 0. At the outer radius, however, the circular orbits are stable; the third term is less important and the system behaves more like the non-relativistic Kepler problem.

When ${\textstyle a}$ is much greater than ${\textstyle r_{\text{s}}}$ (the classical case), these formulae become approximately

${\displaystyle {\begin{aligned}r_{\text{outer}}&\approx {\frac {2a^{2}}{r_{\text{s}}}}\\[3pt]r_{\text{inner}}&\approx {\frac {3}{2}}r_{\text{s}}\end{aligned}}}$

The stable and unstable radii are plotted versus the normalized angular momentum ${\textstyle {\frac {a}{r_{\text{s}}}}={\frac {L}{mcr_{\text{s}}}}}$ in blue and red, respectively. These curves meet at a unique circular orbit (green circle) when the normalized angular momentum equals the square root of three. For comparison, the classical radius predicted from the centripetal acceleration and Newton's law of gravity is plotted in black.

Substituting the definitions of ${\textstyle a}$ and r_s into r_outer yields the classical formula for a particle of mass ${\textstyle m}$ orbiting a body of mass ${\textstyle M}$.

${\displaystyle r_{\mathrm {outer} }^{3}={\frac {G(M+m)}{\omega _{\varphi }^{2}}}}$

where ω_φ is the orbital angular speed of the particle. This formula is obtained in non-relativistic mechanics by setting the centrifugal force equal to the Newtonian gravitational force:

${\displaystyle {\frac {GMm}{r^{2}}}=\mu \omega _{\varphi }^{2}r}$

Where ${\textstyle \mu }$ is the reduced mass.

In our notation, the classical orbital angular speed equals

${\displaystyle \omega _{\varphi }^{2}\approx {\frac {GM}{r_{\mathrm {outer} }^{3}}}=\left({\frac {r_{\text{s}}c^{2}}{2r_{\mathrm {outer} }^{3}}}\right)=\left({\frac {r_{\text{s}}c^{2}}{2}}\right)\left({\frac {r_{\text{s}}^{3}}{8a^{6}}}\right)={\frac {c^{2}r_{\text{s}}^{4}}{16a^{6}}}}$

At the other extreme, when a² approaches 3r_s² from above, the two radii converge to a single value

${\displaystyle r_{\mathrm {outer} }\approx r_{\mathrm {inner} }\approx 3r_{\text{s}}}$

The quadratic solutions above ensure that r_outer is always greater than 3r_s, whereas r_inner lies between 3⁄2 r_s and 3r_s. Circular orbits smaller than 3⁄2 r_s are not possible. For massless particles, a goes to infinity, implying that there is a circular orbit for photons at r_inner = 3⁄2 r_s. The sphere of this radius is sometimes known as the photon sphere.

Precession of elliptical orbits

In the non-relativistic Kepler problem, a particle follows the same perfect ellipse (red orbit) eternally. General relativity introduces a third force that attracts the particle slightly more strongly than Newtonian gravity, especially at small radii. This third force causes the particle's elliptical orbit to precess (cyan orbit) in the direction of its rotation; this effect has been measured in Mercury, Venus and Earth. The yellow dot within the orbits represents the center of attraction, such as the Sun.

The orbital precession rate may be derived using this radial effective potential V. A small radial deviation from a circular orbit of radius r_outer will oscillate stably with an angular frequency

${\displaystyle \omega _{r}^{2}={\frac {1}{m}}\left[{\frac {d^{2}V}{dr^{2}}}\right]_{r=r_{\mathrm {outer} }}}$

which equals

${\displaystyle \omega _{r}^{2}=\left({\frac {c^{2}r_{\text{s}}}{2r_{\mathrm {outer} }^{4}}}\right)\left(r_{\mathrm {outer} }-r_{\mathrm {inner} }\right)=\omega _{\varphi }^{2}{\sqrt {1-{\frac {3r_{\text{s}}^{2}}{a^{2}}}}}}$

Taking the square root of both sides and performing a Taylor series expansion yields

${\displaystyle \omega _{r}=\omega _{\varphi }\left[1-{\frac {3r_{\text{s}}^{2}}{4a^{2}}}+{\mathcal {O}}\left({\frac {r_{\text{s}}^{4}}{a^{4}}}\right)\right]}$

Multiplying by the period T of one revolution gives the precession of the orbit per revolution

${\displaystyle \delta \varphi =T\left(\omega _{\varphi }-\omega _{r}\right)\approx 2\pi \left({\frac {3r_{\text{s}}^{2}}{4a^{2}}}\right)={\frac {3\pi m^{2}c^{2}}{2L^{2}}}r_{\text{s}}^{2}}$

where we have used ω_φT = 2п and the definition of the length-scale a. Substituting the definition of the Schwarzschild radius r_s gives

${\displaystyle \delta \varphi \approx {\frac {3\pi m^{2}c^{2}}{2L^{2}}}\left({\frac {4G^{2}M^{2}}{c^{4}}}\right)={\frac {6\pi G^{2}M^{2}m^{2}}{c^{2}L^{2}}}}$

This may be simplified using the elliptical orbit's semiaxis A and eccentricity e related by the formula

${\displaystyle {\frac {h^{2}}{G(M+m)}}=A\left(1-e^{2}\right)}$

to give the precession angle

${\displaystyle \delta \varphi \approx {\frac {6\pi G(M+m)}{c^{2}A\left(1-e^{2}\right)}}}$

Mathematical derivations of the orbital equation