Perihelion precession in the Schwarzschild spacetime
A massive particle on a bound orbit around a static mass, whose perihelion advances every turn.
- Name:
schwarzschild_precession - Field:
physics/celestial_mechanics - Formalism: ordinary differential equations
- References: Schwichtenberg 2020; Satz 2022; Schwarzschild 1916; Einstein 1915; Carroll 2004
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('schwarzschild_precession')
entry = describe('schwarzschild_precession')
trajectory = create_simulation('schwarzschild_precession').prepare().run().trajectory_
Factory
@register_system(
field='physics/celestial_mechanics',
tags=('general relativity', 'ode', 'second order', 'gravity', 'geodesics'),
references=('Schwichtenberg 2020', 'Satz 2022', 'Schwarzschild 1916', 'Einstein 1915', 'Carroll 2004'),
defaults={'until': 20000.0, 'dt': 10.0, 'solver': GEODESIC_SOLVER},
)
def schwarzschild_precession(
semi_major_axis: float = 100.0,
eccentricity: float = 0.5,
*,
mass: float = 1.0,
) -> System:
"""Return a massive particle on a bound orbit around a static mass, whose perihelion advances every turn.
The particle follows a timelike geodesic of the Schwarzschild spacetime `SchwarzschildSpace`, in units with
`G = c = 1` and with its proper time as the time of the system. It starts at the periapsis
`r₁ = a (1 - e)` of an orbit whose apoapsis is `r₂ = a (1 + e)`, with the energy `E` and the angular momentum `L`
for which both are turning points of `E² = (1 - 2M/r) (1 + L²/r²) + (dr/dτ)²`. Unlike a Kepler ellipse, the orbit
does not close: the periapsis advances by `6 π M / (a (1 - e²))` per turn when the orbit is far from the
horizon, or `6 π G M / (c² a (1 - e²))` in other units, which accounts for the anomalous 43 arcseconds per century
of Mercury.
Args:
semi_major_axis:
The semi-major axis `a` of the orbit.
eccentricity:
The eccentricity `e` of the orbit, in `[0, 1)`.
mass:
The mass `M` of the central body.
Returns:
system:
The system, with the position `particle.x = (t, r, theta, phi)` and the four-velocity `particle.x_d1`.
Raises:
ValidationError:
The eccentricity is outside `[0, 1)`, or the orbit comes closer than the innermost stable circular orbit.
Examples:
>>> from skcomplex.systems.physics.celestial_mechanics import schwarzschild_precession
>>> round(schwarzschild_precession().elements[0].x[1], 6)
50.0
"""
if not 0 <= eccentricity < 1:
error_msg = f'The eccentricity is {eccentricity!r}. Give a value in [0, 1) for a bound orbit.'
raise ValidationError(error_msg)
periapsis, apoapsis = semi_major_axis * (1 - eccentricity), semi_major_axis * (1 + eccentricity)
if periapsis <= 6 * mass:
error_msg = (
f'The periapsis {periapsis} lies inside the innermost stable circular orbit at 6 M = {6 * mass}. Give a '
'larger semi-major axis or a smaller eccentricity.'
)
raise ValidationError(error_msg)
lapses = [1 - 2 * mass / radius for radius in (periapsis, apoapsis)]
if eccentricity == 0:
momentum_squared = mass * periapsis**2 / (periapsis - 3 * mass)
else:
momentum_squared = (lapses[1] - lapses[0]) / (lapses[0] / periapsis**2 - lapses[1] / apoapsis**2)
energy = math.sqrt(lapses[0] * (1 + momentum_squared / periapsis**2))
particle = GeodesicParticleElement(
'particle',
x=[0.0, periapsis, EQUATOR, 0.0],
x_d1=[energy / lapses[0], 0.0, 0.0, math.sqrt(momentum_squared) / periapsis**2],
)
space = SchwarzschildSpace(mass=mass)
geodesic = GeodesicInteraction('geodesic', target='particle', space=space)
return System(elements=[particle], interactions=[geodesic], space=space)
Components
Elements:
GeodesicParticleElement, fromElement
Interactions:
GeodesicInteraction
class GeodesicParticleElement(Element):
"""Test particle at the spacetime position `x = (t, r, theta, phi)` around a central mass, with its four-velocity.
The position is a vector state of order 2 in the Schwarzschild or Boyer-Lindquist coordinates of the space, and its
companion `x_d1` is the four-velocity, the derivative of the position by the proper time of a massive particle or by
the affine parameter of a light ray. The particle is drawn in the equatorial plane, at `(r cos phi, r sin phi)`.
Args:
label:
The label of the particle.
x:
The initial position `(t, r, theta, phi)`.
x_d1:
The initial four-velocity.
visual:
How the particle is drawn, or `None` for its class visual.
Examples:
>>> from skcomplex.systems.physics.celestial_mechanics import GeodesicParticleElement
>>> GeodesicParticleElement('particle', x=[0.0, 10.0, 1.5708, 0.0]).x
[0.0, 10.0, 1.5708, 0.0]
"""
visual = VisualSpec('bodies', {'positions': _equatorial_position, 'trail': None, 'marks': _central_mass})
def __init__(
self,
label: str,
x: Annotated[npt.ArrayLike, StateSpec(vector=True, order=2)] = 0.0,
x_d1: npt.ArrayLike = 0.0,
visual: VisualSpec | None = None,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.x_d1 = x_d1
class GeodesicInteraction(Interaction):
"""Geodesic equation `x''^mu = -Γ^mu_(alpha beta)(x) x'^alpha x'^beta` of a point of a pseudo-Riemannian space.
The state is a vector of order 2 with one component per coordinate of the space, the position `x^mu` of a test
particle, and its companion holds the tangent `x'^mu`. The time of the system is the affine parameter of the
geodesic, the proper time of a timelike geodesic when the tangent starts with `g_(mu nu) x'^mu x'^nu = -1`. A
tangent that starts null, with `g_(mu nu) x'^mu x'^nu = 0`, follows a light ray. Either norm is conserved along the
geodesic, so it measures the error of a run. The Christoffel symbols are those the space derives from its metric.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state is the position.
state:
The name of the vector state of order 2 that is the position.
space:
The pseudo-Riemannian space whose geodesics the point follows, usually the space of the system too.
Attributes:
space_:
The prepared space, whose Christoffel symbols the rate reads.
Raises:
ValidationError:
The state is not a name.
SpaceError:
The space is not a pseudo-Riemannian space, raised by `prepare`.
Examples:
>>> from sage.all__sagemath_symbolics import SR, matrix
>>> from skcomplex.core.components.elements import PointElement
>>> from skcomplex.core.components.interactions import GeodesicInteraction
>>> from skcomplex.core.components.spaces import PseudoRiemannianSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> r, phi = SR.var('r phi')
>>> plane = PseudoRiemannianSpace(matrix(SR, [[1, 0], [0, r**2]]), (r, phi))
>>> line = GeodesicInteraction('line', target='particle', space=plane)
>>> particle = PointElement('particle', x=[1.0, 0.0], x_d1=[0.0, 1.0])
>>> system = System(elements=[particle], interactions=[line], space=plane)
>>> trajectory = Simulation(system, dt=0.001).prepare().run(until=1.0).trajectory_
>>> [round(float(value), 3) for value in trajectory.values('particle.x')[-1]]
[1.414, 0.785]
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state',))],
state: str = 'x',
space: PseudoRiemannianSpace | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.space = space
def _check(self) -> None:
"""Raise an error for a state that is not a name.
Raises:
ValidationError:
The type declares states, or the state is not a name.
"""
super()._check()
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} moves the state {self.state!r}. Give the name of a state.'
)
raise ValidationError(error_msg)
def _prepare(self, context: BaseSpace | None = None) -> None:
"""Set `space_`, the prepared pseudo-Riemannian space.
Args:
context:
The space the interaction inherits from its system.
Raises:
SpaceError:
The space is not a pseudo-Riemannian space.
"""
super()._prepare(context)
if not isinstance(self.space, PseudoRiemannianSpace):
error_msg = (
f'{type(self).__name__} {self.label!r} moves a point along the geodesics of {self.space!r}. Give a '
'`PseudoRiemannianSpace`.'
)
raise SpaceError(error_msg)
self.space_ = self.space.prepare()
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the second derivative of the position, minus the connection on the tangent.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the target.
Returns:
laws:
The rate law of the position.
Raises:
ValidationError:
The position has not one component per coordinate of the space.
"""
position = getattr(self._get_member(namespace, 'target'), self.state)
components = list(position) if isinstance(position, FreeModuleElement) else [position]
coordinates = self.space_.coordinates_
if len(components) != len(coordinates):
error_msg = (
f'The position {self.state!r} of {self.label!r} has {len(components)} components for the '
f'{len(coordinates)} coordinates of the space. Give a vector state of the dimension of the space.'
)
raise ValidationError(error_msg)
at = dict(zip(coordinates, components, strict=True))
tangent = [derivative(component, 1) for component in components]
indices = range(len(components))
rate = [
-sum(
(
SR(self.space_.christoffel_[mu][alpha][beta]).subs(at) * tangent[alpha] * tangent[beta]
for alpha in indices
for beta in indices
if not SR(self.space_.christoffel_[mu][alpha][beta]).is_zero()
),
SR(0),
)
for mu in indices
]
return [RateLaw(position, vector(SR, rate) if isinstance(position, FreeModuleElement) else rate[0])]
Equations
particle.x[0]: diff(x_particle_0(t), t, t) == -2*diff(x_particle_0(t), t)*diff(x_particle_1(t), t)/((x_particle_1(t) - 2)*x_particle_1(t))
particle.x[1]: diff(x_particle_1(t), t, t) == (x_particle_1(t) - 2)*sin(x_particle_2(t))^2*diff(x_particle_3(t), t)^2 + (x_particle_1(t) - 2)*diff(x_particle_2(t), t)^2 - (x_particle_1(t) - 2)*diff(x_particle_0(t), t)^2/x_particle_1(t)^3 + (x_particle_1(t) - 2)*diff(x_particle_1(t), t)^2/((2/x_particle_1(t) - 1)^2*x_particle_1(t)^3)
particle.x[2]: diff(x_particle_2(t), t, t) == cos(x_particle_2(t))*sin(x_particle_2(t))*diff(x_particle_3(t), t)^2 - 2*diff(x_particle_1(t), t)*diff(x_particle_2(t), t)/x_particle_1(t)
particle.x[3]: diff(x_particle_3(t), t, t) == -2*cos(x_particle_2(t))*diff(x_particle_2(t), t)*diff(x_particle_3(t), t)/sin(x_particle_2(t)) - 2*diff(x_particle_1(t), t)*diff(x_particle_3(t), t)/x_particle_1(t)