Cosmology

Cosmology: the expansion of a homogeneous universe and the gravitational collapse of matter.

DustShellElement(label, visual=None, *, R=1.0, R_d1=0.0, energy=0.0)

Bases: Element

Shell of a spherical ball of dust, with the areal radius R of order 2 and the energy per unit mass energy.

Parameters:

Name	Type	Description	Default
label	str	The label of the shell.	required
visual	VisualSpec | None	How the shell is drawn, or None for none.	None
R	Annotated[float, StateSpec(domain=open(0, oo), order=2)]	The initial areal radius.	1.0
R_d1	float	The initial rate of change of the areal radius.	0.0
energy	Annotated[float, ParameterSpec(domain=open(-0.5, oo))]	The energy per unit mass, negative for a bound shell.	0.0

Examples:

>>> from skcomplex.systems.physics.cosmology import DustShellElement
>>> DustShellElement('shell', R=2.0).R
2.0

Source code in src/skcomplex/systems/physics/cosmology/_metric_cosmologies.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    R: Annotated[float, StateSpec(domain=RealSet.open(0, oo), order=2)] = 1.0,
    R_d1: float = 0.0,
    energy: Annotated[float, ParameterSpec(domain=RealSet.open(-0.5, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.R = R
    self.R_d1 = R_d1
    self.energy = energy

HomogeneousUniverseElement(label, visual=None, *, a=1.0, a_d1=1.0, H0=1.0, omega_matter=0.3, omega_radiation=0.0, omega_lambda=0.7)

Bases: Element

Homogeneous universe of matter, radiation, and vacuum energy, with the scale factor a of order 2.

The densities are measured by their parameters Ω today, at a = 1, relative to the critical density 3 H0² / κ: matter dilutes as a⁻³, radiation as a⁻⁴, and the vacuum energy of a cosmological constant stays.

Parameters:

Name	Type	Description	Default
label	str	The label of the universe.	required
visual	VisualSpec | None	How the universe is drawn, or None for none.	None
a	Annotated[float, StateSpec(domain=open(0, oo), order=2)]	The initial scale factor.	1.0
a_d1	float	The initial rate of change of the scale factor.	1.0
H0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Hubble constant today.	1.0
omega_matter	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The density parameter of matter.	0.3
omega_radiation	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The density parameter of radiation.	0.0
omega_lambda	Annotated[float, ParameterSpec()]	The density parameter of the cosmological constant.	0.7

Examples:

>>> from skcomplex.systems.physics.cosmology import HomogeneousUniverseElement
>>> HomogeneousUniverseElement('universe', omega_radiation=1e-4).omega_radiation
0.0001

Source code in src/skcomplex/systems/physics/cosmology/_metric_cosmologies.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a: Annotated[float, StateSpec(domain=RealSet.open(0, oo), order=2)] = 1.0,
    a_d1: float = 1.0,
    H0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    omega_matter: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    omega_radiation: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    omega_lambda: Annotated[float, ParameterSpec()] = 0.7,
) -> None:
    super().__init__(label, visual=visual)
    self.a = a
    self.a_d1 = a_d1
    self.H0 = H0
    self.omega_matter = omega_matter
    self.omega_radiation = omega_radiation
    self.omega_lambda = omega_lambda

LemaitreTolmanBondiSpace(*, label='space')

Bases: PseudoRiemannianSpace

Spherical spacetime of dust in comoving coordinates, the Lemaître-Tolman-Bondi ansatz.

The line element in (t, r, theta, phi) is ds² = -dt² + (∂R/∂r)² dr² / (1 + 2 E(r)) + R(t, r)² (dθ² + sin²θ dφ²), where the areal radius R(t, r) of the shell of comoving radius r is the unknown and E(r) is the energy of the shell per unit mass. They are the functions R and energy, a state and a parameter of the target of EinsteinInteraction.

Parameters:

Name	Type	Description	Default
label	str	The label of the space in its system.	'space'

Examples:

>>> from skcomplex.systems.physics.cosmology import LemaitreTolmanBondiSpace
>>> space = LemaitreTolmanBondiSpace().prepare()
>>> space.metric_[2, 2]
R(t, r)^2

Source code in src/skcomplex/systems/physics/cosmology/_metric_cosmologies.py

def __init__(self, *, label: str = 'space') -> None:
    super().__init__(curvature=True, label=label)

RobertsonWalkerSpace(k=0.0, *, label='space')

Bases: PseudoRiemannianSpace

Homogeneous and isotropic spacetime of the Robertson-Walker ansatz with the unknown scale factor a(t).

The line element in the comoving coordinates (t, r, theta, phi) is ds² = -dt² + a(t)² (dr² / (1 - k r²) + r² (dθ² + sin²θ dφ²)), with the spatial curvature k. The scale factor is the function a that EinsteinInteraction solves for, the state a of its target.

Parameters:

Name	Type	Description	Default
k	float	The spatial curvature, negative for an open, zero for a flat, and positive for a closed universe.	0.0
label	str	The label of the space in its system.	'space'

Examples:

>>> from skcomplex.systems.physics.cosmology import RobertsonWalkerSpace
>>> space = RobertsonWalkerSpace(k=0.0).prepare()
>>> (space.inverse_metric_ * space.einstein_)[0, 0]
-3*diff(a(t), t)^2/a(t)^2

Source code in src/skcomplex/systems/physics/cosmology/_metric_cosmologies.py

def __init__(self, k: float = 0.0, *, label: str = 'space') -> None:
    super().__init__(curvature=True, label=label)
    self.k = k

ScaleFactorElement(label, visual=None, *, a=1.0, a_d1=0.0, H0=1.0, omega_matter=0.3, omega_lambda=0.7)

Bases: Element

Homogeneous universe whose size is the scale factor a, a state of order 2 with the expansion rate a_d1.

Parameters:

Name	Type	Description	Default
label	str	The label of the universe.	required
visual	VisualSpec | None	How the universe is drawn, or None for none.	None
a	Annotated[float, StateSpec(domain=open(0, oo), order=2)]	The initial scale factor.	1.0
a_d1	float	The initial rate of change of the scale factor.	0.0
H0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Hubble constant.	1.0
omega_matter	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The density parameter of matter.	0.3
omega_lambda	Annotated[float, ParameterSpec()]	The density parameter of dark energy.	0.7

Examples:

>>> from skcomplex.systems.physics.cosmology import ScaleFactorElement
>>> ScaleFactorElement('universe', H0=2.0).H0
2.0

Source code in src/skcomplex/systems/physics/cosmology/_expansion.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a: Annotated[float, StateSpec(domain=RealSet.open(0, oo), order=2)] = 1.0,
    a_d1: float = 0.0,
    H0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    omega_matter: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    omega_lambda: Annotated[float, ParameterSpec()] = 0.7,
) -> None:
    super().__init__(label, visual=visual)
    self.a = a
    self.a_d1 = a_d1
    self.H0 = H0
    self.omega_matter = omega_matter
    self.omega_lambda = omega_lambda

flrw_cosmology(a=0.3, *, H0=1.0, omega_matter=0.3, omega_radiation=0.0, omega_lambda=0.7)

Return a homogeneous universe whose expansion follows from the Einstein equations of the Robertson-Walker metric.

The metric of RobertsonWalkerSpace gives, through its Einstein tensor, the two Friedmann equations of a perfect fluid of density rho and pressure p. The space-space component G^r_r = -(2 a a'' + a'² + k) / a² = κ p is the law of the scale factor, and the time-time component G^t_t = -3 (a'² + k) / a² = -κ rho is the constraint (a'/a)² = κ rho / 3 - k / a², which sets the initial rate and which the evolution keeps. The fluid holds matter of pressure 0, radiation of pressure rho / 3, and the vacuum energy of a cosmological constant of pressure -rho, with the density parameters today, at a = 1, and the curvature k = -Ω_k H0² of the remainder Ω_k = 1 - Ω_m - Ω_r - Ω_Λ. Units have G = c = 1.

Parameters:

Name	Type	Description	Default
a	float	The initial scale factor.	0.3
H0	float	The Hubble constant, the expansion rate at a = 1.	1.0
omega_matter	float	The density parameter of matter.	0.3
omega_radiation	float	The density parameter of radiation.	0.0
omega_lambda	float	The density parameter of the cosmological constant.	0.7

Returns:

Name	Type	Description
system	System	The system, with the scale factor universe.a of order 2 and its rate universe.a_d1.

Raises:

Type	Description
ValidationError	The universe does not expand at the initial scale factor.

Examples:

>>> from skcomplex.systems.physics.cosmology import flrw_cosmology
>>> flrw_cosmology(1.0, omega_matter=1.0, omega_lambda=0.0).elements[0].a_d1
1.0

Source code in src/skcomplex/systems/physics/cosmology/_metric_cosmologies.py

@register_system(
    field='physics/cosmology',
    tags=('cosmology', 'general relativity', 'ode', 'second order', 'gravity', 'dark energy'),
    references=('Schwichtenberg 2020', 'Satz 2022', 'Friedmann 1922', 'Lemaître 1927', 'Carroll 2004'),
    defaults={'until': 2.0, 'dt': 0.02, 'solver': ScipyODESolver(method='DOP853', rtol=1e-10, atol=1e-12)},
)
def flrw_cosmology(
    a: float = 0.3,
    *,
    H0: float = 1.0,
    omega_matter: float = 0.3,
    omega_radiation: float = 0.0,
    omega_lambda: float = 0.7,
) -> System:
    """Return a homogeneous universe whose expansion follows from the Einstein equations of the Robertson-Walker metric.

    The metric of `RobertsonWalkerSpace` gives, through its Einstein tensor, the two Friedmann equations of a perfect
    fluid of density `rho` and pressure `p`. The space-space component `G^r_r = -(2 a a'' + a'² + k) / a² = κ p` is the
    law of the scale factor, and the time-time component `G^t_t = -3 (a'² + k) / a² = -κ rho` is the constraint `(a'/a)²
    = κ rho / 3 - k / a²`, which sets the initial rate and which the evolution keeps. The fluid holds matter of pressure
    0, radiation of pressure `rho / 3`, and the vacuum energy of a cosmological constant of pressure `-rho`, with the
    density parameters today, at `a = 1`, and the curvature `k = -Ω_k H0²` of the remainder `Ω_k = 1 - Ω_m - Ω_r - Ω_Λ`.
    Units have `G = c = 1`.

    Args:
        a:
            The initial scale factor.

        H0:
            The Hubble constant, the expansion rate at `a = 1`.

        omega_matter:
            The density parameter of matter.

        omega_radiation:
            The density parameter of radiation.

        omega_lambda:
            The density parameter of the cosmological constant.

    Returns:
        system:
            The system, with the scale factor `universe.a` of order 2 and its rate `universe.a_d1`.

    Raises:
        ValidationError:
            The universe does not expand at the initial scale factor.

    Examples:
        >>> from skcomplex.systems.physics.cosmology import flrw_cosmology
        >>> flrw_cosmology(1.0, omega_matter=1.0, omega_lambda=0.0).elements[0].a_d1
        1.0
    """
    curvature = 1 - omega_matter - omega_radiation - omega_lambda
    square = H0**2 * (omega_matter / a + omega_radiation / a**2 + omega_lambda * a**2 + curvature)
    if square < 0:
        error_msg = (
            f'The Friedmann constraint gives a negative squared expansion rate {square} at a = {a}. Choose a and the '
            'density parameters for which it is positive.'
        )
        raise ValidationError(error_msg)
    universe = HomogeneousUniverseElement(
        'universe',
        a=a,
        a_d1=math.sqrt(square),
        H0=H0,
        omega_matter=omega_matter,
        omega_radiation=omega_radiation,
        omega_lambda=omega_lambda,
    )
    einstein = EinsteinInteraction(
        'einstein',
        target='universe',
        states='a',
        space=RobertsonWalkerSpace(k=-curvature * H0**2),
        components=[(1, 1)],
        source=_get_perfect_fluid,
        coupling=COUPLING,
    )
    return System(elements=[universe], interactions=[einstein])

friedmann_expansion(a=1.0, *, H0=1.0, omega_matter=0.3, omega_lambda=0.7)

Return the scale factor of a homogeneous universe of matter and dark energy, after Friedmann and Lemaître.

The acceleration equation a² a'' = -4 pi G M / 3 + Lambda a³ / 3 is written with the density parameters of matter Ω_m = 8 pi G M / (3 H0²) and of the cosmological constant Ω_Λ = Lambda / (3 H0²) in units of the Hubble constant H0 at a = 1, so a² a'' = -Ω_m H0² / 2 + Ω_Λ H0² a³. The initial expansion rate follows from the Friedmann equation a'² = H0² (Ω_m / a + Ω_Λ a² + 1 - Ω_m - Ω_Λ), which stays a first integral of the motion. Matter decelerates the expansion and dark energy accelerates it once a exceeds (Ω_m / (2 Ω_Λ))^(1/3).

Parameters:

Name	Type	Description	Default
a	float	The initial scale factor.	1.0
H0	float	The Hubble constant, the expansion rate at a = 1.	1.0
omega_matter	float	The density parameter of matter.	0.3
omega_lambda	float	The density parameter of the cosmological constant.	0.7

Returns:

Name	Type	Description
system	System	The system, with the scale factor universe.a of order 2 and its rate universe.a_d1.

Raises:

Type	Description
ValidationError	The universe does not expand at the initial scale factor.

Examples:

>>> from skcomplex.systems.physics.cosmology import friedmann_expansion
>>> friedmann_expansion(omega_matter=1.0, omega_lambda=0.0).elements[0].a_d1
1.0

Source code in src/skcomplex/systems/physics/cosmology/_expansion.py

@register_system(
    field='physics/cosmology',
    tags=('cosmology', 'ode', 'second order', 'gravity', 'dark energy'),
    references=('Friedmann 1922', 'Lemaître 1927', 'Satz 2022'),
    defaults={'until': 2.0, 'dt': 0.02},
)
def friedmann_expansion(
    a: float = 1.0,
    *,
    H0: float = 1.0,
    omega_matter: float = 0.3,
    omega_lambda: float = 0.7,
) -> System:
    """Return the scale factor of a homogeneous universe of matter and dark energy, after Friedmann and Lemaître.

    The acceleration equation `a² a'' = -4 pi G M / 3 + Lambda a³ / 3` is written with the density parameters of
    matter `Ω_m = 8 pi G M / (3 H0²)` and of the cosmological constant `Ω_Λ = Lambda / (3 H0²)` in units of the Hubble
    constant `H0` at `a = 1`, so `a² a'' = -Ω_m H0² / 2 + Ω_Λ H0² a³`. The initial expansion rate follows from the
    Friedmann equation `a'² = H0² (Ω_m / a + Ω_Λ a² + 1 - Ω_m - Ω_Λ)`, which stays a first integral of the motion.
    Matter decelerates the expansion and dark energy accelerates it once `a` exceeds `(Ω_m / (2 Ω_Λ))^(1/3)`.

    Args:
        a:
            The initial scale factor.

        H0:
            The Hubble constant, the expansion rate at `a = 1`.

        omega_matter:
            The density parameter of matter.

        omega_lambda:
            The density parameter of the cosmological constant.

    Returns:
        system:
            The system, with the scale factor `universe.a` of order 2 and its rate `universe.a_d1`.

    Raises:
        ValidationError:
            The universe does not expand at the initial scale factor.

    Examples:
        >>> from skcomplex.systems.physics.cosmology import friedmann_expansion
        >>> friedmann_expansion(omega_matter=1.0, omega_lambda=0.0).elements[0].a_d1
        1.0
    """
    square = H0**2 * (omega_matter / a + omega_lambda * a**2 + 1 - omega_matter - omega_lambda)
    if square < 0:
        error_msg = (
            f'The Friedmann equation gives a negative squared expansion rate {square} at a = {a}. Choose a, '
            'omega_matter, and omega_lambda for which it is positive.'
        )
        raise ValidationError(error_msg)
    universe = ScaleFactorElement(
        'universe',
        a=a,
        a_d1=math.sqrt(square),
        H0=H0,
        omega_matter=omega_matter,
        omega_lambda=omega_lambda,
    )
    acceleration = GradientInteraction(
        'acceleration',
        target='universe',
        states='a',
        potential=_get_expansion_potential,
    )
    return System(elements=[universe], interactions=[acceleration])

oppenheimer_snyder_collapse(n_shells=32, *, radius=10.0, mass=1.0)

Return a ball of dust at rest that collapses under its own gravity into a black hole, as Oppenheimer and Snyder.

The ball is a field of comoving shells on a radial lattice, the Lemaître-Tolman-Bondi ansatz of LemaitreTolmanBondiSpace in units with G = c = 1. Dust has no pressure, so the radial component of the Einstein equations, G^r_r = -(R'² + 2 R R'' - 2 E) / R² = 0 with ' the derivative by proper time, is the law of the areal radius R of every shell, and its first integral R'² = 2 M(r) / R + 2 E(r), from G^t_t, holds the mass M(r) inside the shell. A uniform ball of mass M and radius R₀ at rest has M(r) = M (r / R₀)³ and E(r) = -M(r) / r, so every shell follows the same cycloid R = r (1 + cos η) / 2 with t = (η + sin η) / (2 √k) and k = 2 M / R₀³, as the interior of Oppenheimer and Snyder. The surface crosses its Schwarzschild radius 2 M and the ball reaches zero size at the proper time π / (2 √k).

Parameters:

Name	Type	Description	Default
n_shells	int	The number of shells, which split the comoving radius R₀ evenly, the outermost at the surface.	32
radius	float	The initial radius R₀ of the ball.	10.0
mass	float	The mass M of the ball.	1.0

Returns:

Name	Type	Description
system	System	The system, with the areal radius ball.R of every shell and its rate ball.R_d1.

Raises:

Type	Description
ValidationError	The ball starts inside its Schwarzschild radius.

Examples:

>>> from skcomplex.systems.physics.cosmology import oppenheimer_snyder_collapse
>>> round(float(oppenheimer_snyder_collapse(4).elements[0].values['energy'][-1]), 3)
-0.1

Source code in src/skcomplex/systems/physics/cosmology/_metric_cosmologies.py

@register_system(
    field='physics/cosmology',
    tags=('cosmology', 'general relativity', 'pde', 'second order', 'gravity', 'black holes'),
    references=('Schwichtenberg 2020', 'Satz 2022', 'Oppenheimer and Snyder 1939', 'Tolman 1934', 'Bondi 1947'),
    defaults={'until': 34.5, 'dt': 0.1, 'solver': ScipyODESolver(method='DOP853', rtol=1e-10, atol=1e-12)},
)
def oppenheimer_snyder_collapse(
    n_shells: int = 32,
    *,
    radius: float = 10.0,
    mass: float = 1.0,
) -> System:
    """Return a ball of dust at rest that collapses under its own gravity into a black hole, as Oppenheimer and Snyder.

    The ball is a field of comoving shells on a radial lattice, the Lemaître-Tolman-Bondi ansatz of
    `LemaitreTolmanBondiSpace` in units with `G = c = 1`. Dust has no pressure, so the radial component of the Einstein
    equations, `G^r_r = -(R'² + 2 R R'' - 2 E) / R² = 0` with `'` the derivative by proper time, is the law of the
    areal radius `R` of every shell, and its first integral `R'² = 2 M(r) / R + 2 E(r)`, from `G^t_t`, holds the mass
    `M(r)` inside the shell. A uniform ball of mass `M` and radius `R₀` at rest has `M(r) = M (r / R₀)³` and
    `E(r) = -M(r) / r`, so every shell follows the same cycloid `R = r (1 + cos η) / 2` with
    `t = (η + sin η) / (2 √k)` and `k = 2 M / R₀³`, as the interior of Oppenheimer and Snyder. The surface crosses its
    Schwarzschild radius `2 M` and the ball reaches zero size at the proper time `π / (2 √k)`.

    Args:
        n_shells:
            The number of shells, which split the comoving radius `R₀` evenly, the outermost at the surface.

        radius:
            The initial radius `R₀` of the ball.

        mass:
            The mass `M` of the ball.

    Returns:
        system:
            The system, with the areal radius `ball.R` of every shell and its rate `ball.R_d1`.

    Raises:
        ValidationError:
            The ball starts inside its Schwarzschild radius.

    Examples:
        >>> from skcomplex.systems.physics.cosmology import oppenheimer_snyder_collapse
        >>> round(float(oppenheimer_snyder_collapse(4).elements[0].values['energy'][-1]), 3)
        -0.1
    """
    if radius <= 2 * mass:
        error_msg = f'The radius {radius} is inside the Schwarzschild radius 2 M = {2 * mass}. Give a larger radius.'
        raise ValidationError(error_msg)
    radii = radius * np.arange(1, n_shells + 1) / n_shells
    energies = -mass * radii**2 / radius**3
    ball = FieldElement(
        'ball',
        element=DustShellElement,
        space=LatticeSpace((n_shells,), spacing=radius / n_shells),
        values={'R': radii, 'energy': energies},
    )
    einstein = EinsteinInteraction(
        'einstein',
        target='ball',
        states='R',
        space=LemaitreTolmanBondiSpace(),
        components=[(1, 1)],
        coupling=COUPLING,
        axes={1: 0},
    )
    return System(elements=[ball], interactions=[einstein])