Skip to content

Bending of light by a mass

A light ray that passes a static mass and is deflected towards it.

  • Name: light_bending
  • Field: physics/celestial_mechanics
  • Formalism: ordinary differential equations
  • References: Schwichtenberg 2020; Satz 2022; Einstein 1915; Dyson, Eddington, and Davidson 1920

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('light_bending')
entry = describe('light_bending')
trajectory = create_simulation('light_bending').prepare().run().trajectory_

Factory

@register_system(
    field='physics/celestial_mechanics',
    tags=('general relativity', 'ode', 'second order', 'gravity', 'geodesics', 'light'),
    references=('Schwichtenberg 2020', 'Satz 2022', 'Einstein 1915', 'Dyson, Eddington, and Davidson 1920'),
    defaults={'until': 4000.0, 'dt': 4.0, 'solver': GEODESIC_SOLVER},
)
def light_bending(
    impact_parameter: float = 20.0,
    distance: float = 2000.0,
    *,
    mass: float = 1.0,
) -> System:
    """Return a light ray that passes a static mass and is deflected towards it.

    The ray follows a null geodesic of the Schwarzschild spacetime `SchwarzschildSpace`, in units with `G = c = 1`, with
    the affine parameter as the time of the system. It starts at the distance `D` along the axis behind the mass,
    offset by the impact parameter `b`, aimed parallel to the axis, with the energy `E = 1` and the angular momentum
    `L = b`. Far from the mass the ray bends by `4 M / b`, twice the Newtonian deflection of a particle at the speed of
    light, as the eclipse of 1919 measured for starlight that grazes the Sun.

    Args:
        impact_parameter:
            The impact parameter `b`, the distance of the unbent ray from the mass.

        distance:
            The distance `D` along the axis from which the ray starts.

        mass:
            The mass `M` of the central body.

    Returns:
        system:
            The system, with the position `photon.x = (t, r, theta, phi)` and the tangent `photon.x_d1`.

    Raises:
        ValidationError:
            The ray would be captured by the mass, its impact parameter being at most `3 √3 M`, or it starts closer
            than its impact parameter.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import light_bending
        >>> round(light_bending().elements[0].x_d1[0], 4)
        1.001
    """
    if impact_parameter <= 3 * math.sqrt(3) * mass:
        error_msg = (
            f'The impact parameter {impact_parameter} is at most 3 √3 M = {3 * math.sqrt(3) * mass}, so the mass '
            'captures the ray. Give a larger impact parameter.'
        )
        raise ValidationError(error_msg)
    if distance <= impact_parameter:
        error_msg = f'The distance {distance} is not larger than the impact parameter {impact_parameter}.'
        raise ValidationError(error_msg)
    radius = math.hypot(distance, impact_parameter)
    lapse = 1 - 2 * mass / radius
    radial = -math.sqrt(1 - lapse * impact_parameter**2 / radius**2)
    photon = GeodesicParticleElement(
        'photon',
        x=[0.0, radius, EQUATOR, math.atan2(impact_parameter, -distance)],
        x_d1=[1 / lapse, radial, 0.0, -impact_parameter / radius**2],
    )
    space = SchwarzschildSpace(mass=mass)
    geodesic = GeodesicInteraction('geodesic', target='photon', space=space)
    return System(elements=[photon], interactions=[geodesic], space=space)

Components

Elements:

  • GeodesicParticleElement, from Element

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

photon.x[0]: diff(x_photon_0(t), t, t) == -2*diff(x_photon_0(t), t)*diff(x_photon_1(t), t)/((x_photon_1(t) - 2)*x_photon_1(t))
photon.x[1]: diff(x_photon_1(t), t, t) == (x_photon_1(t) - 2)*sin(x_photon_2(t))^2*diff(x_photon_3(t), t)^2 + (x_photon_1(t) - 2)*diff(x_photon_2(t), t)^2 - (x_photon_1(t) - 2)*diff(x_photon_0(t), t)^2/x_photon_1(t)^3 + (x_photon_1(t) - 2)*diff(x_photon_1(t), t)^2/((2/x_photon_1(t) - 1)^2*x_photon_1(t)^3)
photon.x[2]: diff(x_photon_2(t), t, t) == cos(x_photon_2(t))*sin(x_photon_2(t))*diff(x_photon_3(t), t)^2 - 2*diff(x_photon_1(t), t)*diff(x_photon_2(t), t)/x_photon_1(t)
photon.x[3]: diff(x_photon_3(t), t, t) == -2*cos(x_photon_2(t))*diff(x_photon_2(t), t)*diff(x_photon_3(t), t)/sin(x_photon_2(t)) - 2*diff(x_photon_1(t), t)*diff(x_photon_3(t), t)/x_photon_1(t)