physicskit.relativity

physicskit.relativity hub · History · Examples · API reference

physicskit.relativity: fast, visually captivating numerical General Relativity.

physicskit.relativity brings curved spacetime geometry, black hole physics, ray-traced gravitational lensing, and gravitational wave dynamics into an interactive, computational Python framework, in geometrized units (\(G = c = 1\); see physicskit.relativity.utils.constants):

• Black hole physics (physicskit.relativity.chapters.schwarzschild, physicskit.relativity.chapters.kerr): timelike orbits, perihelion precession, the innermost stable circular orbit (ISCO), the photon sphere, light deflection and the Shapiro delay, disk redshift, tidal forces, frame dragging, the ergosphere, and the Penrose process.

• Spacetime geometry (physicskit.relativity.core.tensors): a numerical differential geometry engine computing Christoffel symbols and curvature tensors for arbitrary metrics (Schwarzschild, Kerr, Reissner-Nordstrom, FLRW, the Alcubierre warp drive).

• Gravitational lensing (physicskit.relativity.core.raytracer, physicskit.relativity.core.kerr_raytracer, physicskit.relativity.visualizers.shadow_render, physicskit.relativity.chapters.lensing): Numba-accelerated backward photon ray-tracers (Schwarzschild and Kerr) rendering black hole shadows and redshift-shaded lensed accretion disks, plus point-source Einstein rings, multiple images, and microlensing magnification.

• Gravitational waves (physicskit.relativity.chapters.gw_merger): binary black hole inspiral chirps and ringdown, and Hulse-Taylor-style orbital decay.

• Cosmology & compact objects (physicskit.relativity.chapters.cosmology, physicskit.relativity.chapters.neutron_star): the FLRW Friedmann equations for cosmic expansion, and the Tolman-Oppenheimer-Volkoff equations for neutron star structure and the maximum neutron star mass.

• Relativity in everyday technology (physicskit.relativity.chapters.timekeeping): the GPS satellite relativistic clock correction.

Examples

>>> from physicskit.relativity.chapters import SchwarzschildBlackHole
>>> bh = SchwarzschildBlackHole(M=1.0)
>>> bh.isco_radius
6.0
>>> bh.photon_sphere_radius
3.0

class physicskit.relativity.BinaryMerger(m1, m2, distance, inclination=0.0)

Bases: object

A compact binary and its leading-order inspiral-merger-ringdown (IMR) waveform.

Parameters:

• m1 (float) – Component masses, in geometrized length units. Use physicskit.relativity.utils.constants.solar_masses_to_geometrized() to convert from solar masses. Unlike SchwarzschildBlackHole or KerrBlackHole, there is no normalized M=1 convention available here: the radiated strain genuinely depends on the absolute physical scale, so an unconverted mass doesn’t just rescale the answer – it silently produces a wrong waveform. This constructor does not itself warn on an implausible value (many legitimate uses – unit tests of a formula’s mathematical properties, mainly – pass small “toy” numbers on purpose, with no claim of physical realism); use check_geometrized_mass() explicitly wherever real-world realism is actually the point (as the visualizers in physicskit.relativity.visualizers.wave_plots do).

• m2 (float) – Component masses, in geometrized length units. Use physicskit.relativity.utils.constants.solar_masses_to_geometrized() to convert from solar masses. Unlike SchwarzschildBlackHole or KerrBlackHole, there is no normalized M=1 convention available here: the radiated strain genuinely depends on the absolute physical scale, so an unconverted mass doesn’t just rescale the answer – it silently produces a wrong waveform. This constructor does not itself warn on an implausible value (many legitimate uses – unit tests of a formula’s mathematical properties, mainly – pass small “toy” numbers on purpose, with no claim of physical realism); use check_geometrized_mass() explicitly wherever real-world realism is actually the point (as the visualizers in physicskit.relativity.visualizers.wave_plots do).

• distance (float) – Distance to the source, in geometrized length units. Same caveat as m1/m2.

• inclination (float, default=0.0) – Inclination angle between the orbital angular momentum and the line of sight, in radians (0 = face-on, pi/2 = edge-on).

property chirp_mass

Chirp mass, \(\mathcal{M} = (m_1 m_2)^{3/5} / (m_1+m_2)^{1/5}\).

full_waveform(t, t_merger)

A continuous inspiral-merger-ringdown toy waveform over an arbitrary time array.

Uses inspiral_strain() for t < t_merger and ringdown_strain() for t >= t_merger, with the ringdown amplitude anchored to the inspiral amplitude evaluated just before merger for continuity.

Parameters:

• t (array_like) – Time array, spanning both before and after t_merger.

• t_merger (float) – Coalescence time; should be strictly greater than every pre-merger sample (the inspiral formula diverges at t = t_merger).

Returns:

• h_plus (ndarray)

• h_cross (ndarray)

Examples

>>> merger = BinaryMerger(m1=30.0, m2=25.0, distance=1000.0)
>>> t = np.linspace(-1.0, 0.02, 500)
>>> hp, hc = merger.full_waveform(t, t_merger=0.0)
>>> hp.shape
(500,)

inspiral_frequency(t, t_merger)

Gravitational-wave frequency during the inspiral, at leading (Newtonian quadrupole) order.

\[f(t) = \frac{1}{\pi}\left[\frac{5}{256(t_{\text{merger}} - t)}\right]^{3/8} \mathcal{M}^{-5/8}\]

Parameters:

• t (float or array_like) – Time(s), with t < t_merger.

• t_merger (float) – Coalescence time.

Returns:

float or ndarray – GW frequency (twice the orbital frequency).

inspiral_phase(t, t_merger, phase_ref=0.0)

Gravitational-wave phase during the inspiral.

\[\Phi(t) = \Phi_{\text{ref}} - 2\left[\frac{t_{\text{merger}} - t} {5\mathcal{M}}\right]^{5/8}\]

The antiderivative of \(d\Phi/dt = 2\pi f(t)\) with inspiral_frequency().

Parameters:

• t (float or array_like) – Time(s), with t < t_merger.

• t_merger (float) – Coalescence time.

• phase_ref (float, default=0.0) – Reference phase (phase at merger).

Returns:

float or ndarray

inspiral_strain(t, t_merger, phase_ref=0.0)

Plus and cross gravitational-wave strain during the inspiral.

\[h_+(t) = \frac{4}{R} \mathcal{M}^{5/3} (\pi f(t))^{2/3} \frac{1 + \cos^2\iota}{2} \cos\Phi(t), \qquad h_\times(t) = \frac{4}{R} \mathcal{M}^{5/3} (\pi f(t))^{2/3} \cos\iota \, \sin\Phi(t)\]

Parameters:

• t (array_like) – Times, with t < t_merger.

• t_merger (float) – Coalescence time.

• phase_ref (float, default=0.0) – Reference phase at merger, forwarded to inspiral_phase().

Returns:

• h_plus (ndarray)

• h_cross (ndarray)

Examples

>>> merger = BinaryMerger(m1=30.0, m2=25.0, distance=1000.0)
>>> t = np.linspace(-2.0, -0.1, 100)
>>> hp, hc = merger.inspiral_strain(t, t_merger=0.0)
>>> hp.shape
(100,)

orbital_period(semi_major_axis)

Newtonian (Keplerian) orbital period at a given semi-major axis.

\[P = 2\pi \sqrt{\frac{a^3}{M}}\]

Valid for the wide, slowly-decaying orbits this method’s counterparts describe (unlike inspiral_frequency(), which covers the final relativistic plunge).

Parameters:

semi_major_axis (float) – Orbital semi-major axis \(a\).

Returns:

float

period_decay_rate(semi_major_axis, eccentricity)

Rate of orbital period decay, \(dP/dt\), from gravitational-wave emission.

Obtained from semi_major_axis_decay_rate() via the chain rule through Kepler’s third law, \(dP/dt = (dP/da)(da/dt)\) with \(dP/da = 3\pi\sqrt{a/M}\).

Parameters:

• semi_major_axis (float) – Orbital semi-major axis \(a\).

• eccentricity (float) – Orbital eccentricity \(e \in [0, 1)\).

Returns:

float – \(dP/dt\) (dimensionless: a rate of change of time per unit time; negative, since the period shrinks).

Examples

>>> import physicskit.relativity.utils.constants as const
>>> m1 = const.solar_masses_to_geometrized(1.4398)
>>> m2 = const.solar_masses_to_geometrized(1.3886)
>>> psr_b1913 = BinaryMerger(m1=m1, m2=m2, distance=1.0)
>>> a = (psr_b1913.total_mass * (27906.98 * const.C_SI) ** 2 / (4 * np.pi**2)) ** (1 / 3)
>>> dPdt = psr_b1913.period_decay_rate(a, eccentricity=0.6171338)
>>> bool(-3e-12 < dPdt < -2e-12)  # observed: -2.4e-12 (dimensionless)
True

qnm_frequency_damping(M_f=None, a_f=None)

Dominant (l=2, m=2, n=0) quasinormal mode frequency and damping time.

Fitting formulas of Berti, Cardoso & Will (2006):

\[f_{\text{QNM}} = \frac{1}{2\pi M_f}\left[1.5251 - 1.1568(1-a_*)^{0.1292}\right]\]

\[Q = 0.7000 + 1.4187(1-a_*)^{-0.4990}, \qquad \tau_{\text{damp}} = \frac{Q}{\pi f_{\text{QNM}}}\]

where \(a_* = a_f/M_f\).

Parameters:

• M_f (float, optional) – Remnant mass. Defaults to remnant_estimate().

• a_f (float, optional) – Remnant spin parameter. Defaults to remnant_estimate().

Returns:

• f_qnm (float) – Ringdown frequency.

• tau_damp (float) – Damping time.

remnant_estimate()

A simple, illustrative estimate of the merger remnant’s mass and spin.

Not a precision numerical-relativity surrogate: it linearly interpolates between the extreme-mass-ratio limit (negligible radiated energy and spin-up, \(\eta \to 0\)) and the well-measured equal-mass, non-spinning result (\(M_f \approx 0.95 M\), \(a_f/M_f \approx 0.69\) at \(\eta = 1/4\)), which is sufficient to demonstrate the qualitative ringdown physics.

Returns:

• M_f (float) – Estimated remnant mass.

• a_f (float) – Estimated remnant spin parameter.

ringdown_strain(t, t_merger, amplitude, M_f=None, a_f=None, phase_ref=0.0)

Damped-sinusoid ringdown strain for \(t \ge t_{\text{merger}}\).

\[h_+(t) = A \, e^{-(t-t_{\text{merger}})/\tau_{\text{damp}}} \cos\left(2\pi f_{\text{QNM}} (t - t_{\text{merger}}) + \Phi_{\text{ref}}\right)\]

with \(h_\times\) the corresponding sine (quarter-cycle phase-shifted) quadrature, modulated by the same inclination dependence as inspiral_strain().

Parameters:

• t (array_like) – Times, with t >= t_merger.

• t_merger (float) – Merger (ringdown onset) time.

• amplitude (float) – Ringdown strain amplitude at \(t = t_{\text{merger}}\), typically chosen to match the inspiral amplitude just before merger for a continuous (if not perfectly smooth) toy waveform.

• M_f (float, optional) – Remnant mass and spin, forwarded to qnm_frequency_damping().

• a_f (float, optional) – Remnant mass and spin, forwarded to qnm_frequency_damping().

• phase_ref (float, default=0.0) – Reference phase at \(t = t_{\text{merger}}\).

Returns:

• h_plus (ndarray)

• h_cross (ndarray)

semi_major_axis_decay_rate(semi_major_axis, eccentricity)

Rate of orbital shrinkage from gravitational-wave emission (Peters 1964).

\[\frac{da}{dt} = -\frac{64}{5} \frac{m_1 m_2 (m_1+m_2)}{a^3} \frac{1 + \frac{73}{24}e^2 + \frac{37}{96}e^4}{(1-e^2)^{7/2}}\]

This is the formula behind the first (indirect) detection of gravitational waves: Hulse and Taylor’s 1974 discovery of the binary pulsar PSR B1913+16, whose orbital period was found to decay at precisely the rate this equation predicts (1993 Nobel Prize in Physics).

Parameters:

• semi_major_axis (float) – Orbital semi-major axis \(a\).

• eccentricity (float) – Orbital eccentricity \(e \in [0, 1)\).

Returns:

float – \(da/dt\) (negative: the orbit shrinks).

property symmetric_mass_ratio

Symmetric mass ratio, \(\eta = m_1 m_2 / M^2\).

property total_mass

Total mass, \(M = m_1 + m_2\).

class physicskit.relativity.FLRWCosmology(H0=70.0, Omega_m=0.3, Omega_r=0.0001, Omega_Lambda=0.7)

Bases: object

A flat or curved FLRW universe with radiation, matter, curvature, and dark energy.

Parameters:

• H0 (float, default=70.0) – Hubble constant today, in km/s/Mpc.

• Omega_m (float, default=0.3) – Matter density parameter today.

• Omega_r (float, default=1e-4) – Radiation density parameter today.

• Omega_Lambda (float, default=0.7) – Dark energy (cosmological constant) density parameter today.

Variables:

Omega_k (float) – Curvature density parameter, \(1 - \Omega_m - \Omega_r - \Omega_\Lambda\) (computed so all density parameters sum to 1, as required by the Friedmann equation at \(a=1\)).

E(a)

Dimensionless expansion rate \(E(a) = H(a)/H_0\).

\[E(a) = \sqrt{\Omega_r a^{-4} + \Omega_m a^{-3} + \Omega_k a^{-2} + \Omega_\Lambda}\]

Parameters:

a (float or array_like) – Scale factor (\(a=1\) today).

Returns:

float or ndarray

age_gyr(a=1.0)

Cosmic age at scale factor a, in Gyr.

\[t(a) = \int_0^a \frac{da'}{a' H(a')}\]

Parameters:

a (float, default=1.0) – Scale factor at which to evaluate the age (a=1 gives the current age of the universe).

Returns:

float

Examples

>>> cosmo = FLRWCosmology(H0=70.0, Omega_m=0.3, Omega_r=0.0, Omega_Lambda=0.7)
>>> 12.0 < cosmo.age_gyr() < 15.0
True

comoving_distance_mpc(z)

Comoving distance to redshift z, in Mpc.

\[D_C(z) = c \int_0^z \frac{dz'}{H(z')}\]

Parameters:

z (float) – Redshift.

Returns:

float

hubble_parameter(a)

Hubble parameter \(H(a) = H_0 E(a)\), in km/s/Mpc.

Parameters:

a (float or array_like) – Scale factor.

Returns:

float or ndarray

luminosity_distance_mpc(z)

Luminosity distance to redshift z, in Mpc.

\[\begin{split}D_L = (1+z) D_M, \qquad D_M = \begin{cases} \frac{D_H}{\sqrt{\Omega_k}} \sinh\left(\sqrt{\Omega_k}\, D_C/D_H\right), & \Omega_k > 0, \\ D_C, & \Omega_k = 0, \\ \frac{D_H}{\sqrt{|\Omega_k|}} \sin\left(\sqrt{|\Omega_k|}\, D_C/D_H\right), & \Omega_k < 0, \end{cases}\end{split}\]

where \(D_H = c/H_0\) is the Hubble distance and \(D_M\) the transverse comoving distance, which differs from the line-of-sight comoving_distance_mpc() \(D_C\) only in a curved universe (Hogg 1999, astro-ph/9905116, eqs. 16 and 21).

Parameters:

z (float) – Redshift.

Returns:

float

Examples

An open, matter-only universe matches Mattig’s closed form, \(D_L = \frac{2 D_H}{\Omega_m^2}\left[\Omega_m z + (\Omega_m - 2)(\sqrt{1+\Omega_m z} - 1)\right]\):

>>> cosmo = FLRWCosmology(H0=70.0, Omega_m=0.3, Omega_r=0.0, Omega_Lambda=0.0)
>>> round(cosmo.luminosity_distance_mpc(1.0), 1)
5872.3

redshift_to_scale_factor(z)

Convert redshift to scale factor, \(a = 1/(1+z)\).

scale_factor_history(n_points=200, a_min=0.001, a_max=3.0)

Tabulate the age of the universe at a range of scale factors, for plotting \(a(t)\).

Parameters:

• n_points (int, default=200) – Number of scale-factor grid points.

• a_min (float) – Range of scale factors to tabulate.

• a_max (float) – Range of scale factors to tabulate.

Returns:

• t_gyr (ndarray of shape (n_points,)) – Cosmic age at each scale factor, in Gyr.

• a (ndarray of shape (n_points,)) – The scale factor grid.

scale_factor_to_redshift(a)

Convert scale factor to redshift, \(z = 1/a - 1\).

class physicskit.relativity.KerrBlackHole(M=1.0, a=0.0)

Bases: object

A rotating black hole of mass M and spin parameter a, in geometrized units.

Parameters:

• M (float, default=1.0) – Mass, in geometrized length units.

• a (float, default=0.0) – Spin parameter \(a = J/M\), with \(0 \le a < M\). a=0 recovers Schwarzschild.

circular_orbit_conserved_quantities(r0, prograde=True)

Conserved specific energy and angular momentum of a circular equatorial orbit.

Standard Bardeen-Press-Teukolsky (1972) result:

\[E = \frac{r^{3/2} - 2Mr^{1/2} + a M^{1/2}} {r^{3/4}\sqrt{r^{3/2} - 3Mr^{1/2} + 2aM^{1/2}}}, \qquad L = \frac{M^{1/2}\left(r^2 - 2aM^{1/2}r^{1/2} + a^2\right)} {r^{3/4}\sqrt{r^{3/2} - 3Mr^{1/2} + 2aM^{1/2}}}\]

for a prograde orbit. For a retrograde orbit, substitute \(a \to -a\) and flip the sign of \(L\): it counter-rotates, so \(L<0\) (d(phi)/d(tau) < 0) in the hole’s frame, which is the sign integrate_equatorial_geodesic() needs to reproduce the circular orbit.

Parameters:

• r0 (float) – Orbital radius; should exceed the appropriate isco_radius() for the orbit to be stable.

• prograde (bool, default=True) – Whether the orbit co-rotates with the black hole.

Returns:

• E (float) – Specific energy.

• L (float) – Specific angular momentum (negative for a retrograde orbit).

property dimensionless_spin

Dimensionless spin \(a_* = a/M \in [0, 1)\).

disk_redshift_factor(r, impact_parameter, prograde=True)

Combined gravitational + Doppler redshift factor for a circularly orbiting disk emitter.

Same construction as physicskit.relativity.chapters.schwarzschild.SchwarzschildBlackHole.disk_redshift_factor(), \(g = 1/[u^t(1 - b\Omega)]\), using the Kerr equatorial metric and keplerian_angular_velocity() in place of their Schwarzschild counterparts.

Parameters:

• r (float or array_like) – Disk radius at emission.

• impact_parameter (float or array_like) – The received photon’s impact parameter \(b = L_{\text{ph}}/E_{\text{ph}}\).

• prograde (bool, default=True) – Whether the disk co-rotates with the black hole.

Returns:

float or ndarray – The redshift factor \(g\).

ergosphere_radius(theta)

Static-limit (ergosphere outer boundary) radius at polar angle theta.

\[r_E(\theta) = M + \sqrt{M^2 - a^2\cos^2\theta}\]

The ergosphere is the region \(r_+ < r < r_E(\theta)\), where the t-Killing vector becomes spacelike: no observer can remain at fixed \((r, \theta, \phi)\) there, no matter how powerful their rocket – they are inevitably dragged around in \(\phi\) by the rotating spacetime itself.

Parameters:

theta (float or array_like) – Polar angle (pi/2 is the equator, where the ergosphere is widest, \(r_E = 2M\)).

Returns:

float or ndarray

frame_dragging_angular_velocity(r)

Angular velocity of a zero-angular-momentum observer (ZAMO) in the equatorial plane.

\[\omega(r) = -\frac{g_{t\phi}}{g_{\phi\phi}} = \frac{2Ma}{r^3 + a^2 r + 2Ma^2}\]

The rate at which locally non-rotating observers are swept around in \(\phi\) purely by the geometry – frame dragging. At large \(r\) this falls off as \(2Ma/r^3\), the classic Lense-Thirring precession rate; it rises sharply approaching the ergosphere, where even a photon fired against the rotation is dragged forward.

Parameters:

r (float or array_like) – Equatorial radial coordinate, \(r > r_+\).

Returns:

float or ndarray

property inner_horizon_radius

Inner (Cauchy) horizon radius, \(r_- = M - \sqrt{M^2 - a^2}\).

integrate_equatorial_geodesic(r0, E, L, mu2, dtau, n_steps)

Integrate an equatorial Kerr geodesic of conserved energy E and momentum L.

Parameters:

• r0 (float) – Initial radius.

• E (float) – Conserved specific energy.

• L (float) – Conserved specific angular momentum.

• mu2 (float) – 1.0 for a timelike geodesic, 0.0 for a null (photon) geodesic.

• dtau (float) – Affine-parameter (or proper-time) step size.

• n_steps (int) – Maximum number of steps.

Returns:

dict of str -> ndarray – Keys "t", "r", "phi", "ur", plus "norm": the four-velocity normalization \(g_{\mu\nu}u^\mu u^\nu\) at every step (should stay close to \(-\mu^2\); its spread measures integration error).

irreducible_mass()

Irreducible mass, \(M_{\text{irr}} = \sqrt{M r_+ / 2}\).

The irreducible mass can never decrease in any classical process (Hawking’s area theorem: it is proportional to \(\sqrt{\text{horizon area}}\)), which is what caps how much rotational energy the Penrose process can extract.

Returns:

float

is_inside_ergosphere(r, theta)

Whether a point \((r, \theta)\) lies inside the ergosphere.

Parameters:

• r (float) – Radial coordinate.

• theta (float) – Polar angle.

Returns:

bool

isco_radius(prograde=True)

Innermost stable circular equatorial orbit radius (Bardeen-Press-Teukolsky 1972).

\[Z_1 = 1 + (1-a_*^2)^{1/3}\left[(1+a_*)^{1/3} + (1-a_*)^{1/3}\right], \quad Z_2 = \sqrt{3a_*^2 + Z_1^2}\]

\[\frac{r_{\text{ISCO}}}{M} = 3 + Z_2 \mp \sqrt{(3-Z_1)(3+Z_1+2Z_2)}\]

with the \(-\) sign for a prograde (co-rotating) orbit and \(+\) for retrograde. Prograde ISCOs shrink toward the horizon as spin increases (\(r_{\text{ISCO}} \to M\) as \(a \to M\)); retrograde ISCOs grow toward \(9M\).

Parameters:

prograde (bool, default=True) – Whether the orbit co-rotates with the black hole.

Returns:

float

keplerian_angular_velocity(r, prograde=True)

Coordinate angular velocity of a circular equatorial geodesic orbit.

\[\Omega(r) = \pm\frac{\sqrt{M}}{r^{3/2} \pm a\sqrt{M}}\]

with the upper sign for a prograde (co-rotating) orbit. Reduces to the Keplerian \(\sqrt{M/r^3}\) at \(a=0\).

Parameters:

• r (float or array_like) – Orbital radius.

• prograde (bool, default=True) – Whether the orbit co-rotates with the black hole.

Returns:

float or ndarray

max_penrose_efficiency()

Maximum fraction of the black hole’s mass extractable via the Penrose process.

\[\eta_{\max} = 1 - \frac{M_{\text{irr}}}{M} = 1 - \sqrt{\frac{r_+}{2M}}\]

Rises from 0 at \(a=0\) to \(1 - 1/\sqrt{2} \approx 29.3\%\) for an extremal (\(a \to M\)) black hole.

Returns:

float

property outer_horizon_radius

Event horizon radius, \(r_+ = M + \sqrt{M^2 - a^2}\).

penrose_energy_gain(initial_energy, fragment_energy_infalling)

Energy gained by the escaping fragment in a Penrose-process split.

A particle of specific energy initial_energy (measured at infinity) splits, inside the ergosphere, into two fragments. If one fragment falls into the horizon with negative specific energy fragment_energy_infalling (possible only inside the ergosphere, where the t-Killing vector is spacelike), energy conservation at the split point requires the escaping fragment to carry away

\[E_{\text{out}} = E_{\text{in}} - E_{\text{infalling}} > E_{\text{in}}\]

Parameters:

• initial_energy (float) – Specific energy of the infalling particle before it splits.

• fragment_energy_infalling (float) – Specific energy of the fragment that falls into the horizon; must be negative for energy to be extracted.

Returns:

float – Specific energy of the escaping fragment.

Examples

>>> bh = KerrBlackHole(M=1.0, a=0.9)
>>> e_out = bh.penrose_energy_gain(initial_energy=1.0, fragment_energy_infalling=-0.05)
>>> e_out > 1.0
True

class physicskit.relativity.NeutronStar(K=100.0, Gamma=2.0)

Bases: object

A neutron star modeled with a polytropic equation of state, solved via the TOV equations.

The equation of state is the polytrope \(P = K \rho_0^\Gamma\) relating pressure to rest-mass density; the total energy density (rest mass plus internal energy) used as the gravitational source is \(\varepsilon = \rho_0 + P/(\Gamma - 1)\). All quantities are in geometrized units; K and Gamma are illustrative “toy” parameters unless calibrated to a specific realistic equation of state.

Parameters:

• K (float, default=100.0) – Polytropic constant.

• Gamma (float, default=2.0) – Polytropic index (\(\Gamma = 2\) is a common toy choice loosely representative of nuclear matter).

energy_density(pressure)

Total energy density \(\varepsilon(P) = (P/K)^{1/\Gamma} + P/(\Gamma-1)\).

Parameters:

pressure (float or array_like) – Pressure.

Returns:

float or ndarray

mass_radius_curve(central_densities)

Trace the mass-radius relation over a range of central densities.

The existence of a maximum mass (a turning point in \(M\) as a function of central density, beyond which increasing the central density decreases the stable mass) is the TOV equation’s central physical prediction.

Parameters:

central_densities (array_like) – Central densities to sample.

Returns:

• masses (ndarray)

• radii (ndarray)

solve(central_density, r_max=50.0, r_start=1e-06)

Integrate the TOV equations outward from the center to the stellar surface.

Parameters:

• central_density (float) – Central rest-mass density \(\rho_0(0)\).

• r_max (float, default=50.0) – Outer integration bound (should exceed any physically expected stellar radius for the given equation of state).

• r_start (float, default=1e-6) – Small starting radius (avoids the coordinate singularity at exactly \(r=0\)).

Returns:

• mass (float) – Total gravitational mass \(M\), read off where the pressure drops to (numerically) zero – the stellar surface.

• radius (float) – Stellar radius \(R\).

• solution (scipy.integrate.OdeSolution) – The full solve_ivp result, for inspecting the interior profile.

Examples

>>> star = NeutronStar(K=100.0, Gamma=2.0)
>>> M, R, sol = star.solve(central_density=1.28e-3)
>>> bool(M > 0.0 and R > 0.0)
True

class physicskit.relativity.PointMassLens(M, D_L, D_S)

Bases: object

A point-mass (weak-field) gravitational lens: Einstein ring, images, and magnification.

Uses the thin-lens approximation with the leading-order deflection angle \(\hat\alpha(b) = 4M/b\), valid whenever the source is not aligned so closely that its image angle approaches the photon sphere (see exact_deflection_angle() for the strong-field case).

Parameters:

• M (float) – Lens mass, in geometrized units.

• D_L (float) – Observer-lens distance, in the same geometrized length units as M.

• D_S (float) – Observer-source distance (\(D_S > D_L\)), in the same units as M.

Notes

Every formula here is a ratio of M, D_L, D_S (angles and magnifications, not an absolute strain), so it is scale-invariant: a normalized convention (M=1, distances “in units of M”) gives exactly the same angles as real SI-consistent geometrized values would. Unlike BinaryMerger, there is no absolute-scale trap here – any self-consistent unit choice works.

Variables:

D_LS (float) – Lens-source distance, \(D_S - D_L\) (flat-space subtraction; this is a local, non-cosmological thin-lens treatment).

einstein_angle()

Einstein ring angular radius, \(\theta_E = \sqrt{4M D_{LS}/(D_L D_S)}\).

The angular radius at which a perfectly aligned (\(\beta=0\)) point source appears as a complete ring.

Returns:

float

image_angles(beta)

Angular positions of the two images of a source at angular offset beta.

Solves the thin-lens equation \(\beta = \theta - \theta_E^2/\theta\) (equivalent to a quadratic in \(\theta\)) exactly:

\[\theta_\pm = \frac{\beta \pm \sqrt{\beta^2 + 4\theta_E^2}}{2}\]

Parameters:

beta (float) – True (unlensed) angular offset of the source from the lens.

Returns:

• theta_plus (float) – Position of the primary image, on the same side as the source (\(\theta_+ > 0\)).

• theta_minus (float) – Position of the secondary image, on the opposite side (\(\theta_- < 0\)), closer to the lens and fainter.

Examples

>>> lens = PointMassLens(M=1.0, D_L=1000.0, D_S=2000.0)
>>> theta_plus, theta_minus = lens.image_angles(beta=0.0)
>>> bool(np.isclose(theta_plus, lens.einstein_angle()))
True
>>> bool(np.isclose(theta_minus, -lens.einstein_angle()))
True

magnification(beta)

Magnification of each image, and their sum (total observed brightening).

Standard point-lens (Paczynski 1986) result, with \(u = \beta/\theta_E\):

\[\mu_\pm = \frac{1}{2}\left[\frac{u^2+2}{u\sqrt{u^2+4}} \pm 1\right]\]

Parameters:

beta (float) – True angular offset of the source from the lens (beta != 0).

Returns:

• mu_plus (float) – Magnification of the primary image (always \(> 1\)).

• mu_minus (float) – Magnification of the secondary image (its absolute value is \(< 1\); it has opposite parity, conventionally signed negative).

• total (float) – Total observed magnification, \(|\mu_+| + |\mu_-|\).

class physicskit.relativity.SchwarzschildBlackHole(M=1.0)

Bases: object

A non-rotating black hole (or any spherical mass) of mass M, in geometrized units.

Parameters:

M (float, default=1.0) – Mass, in geometrized length units (\(GM/c^2\)). Use physicskit.relativity.utils.constants.solar_masses_to_geometrized() to convert from a physical mass in solar masses.

circular_orbit_angular_velocity(r)

Coordinate angular velocity \(\Omega = d\phi/dt\) of a circular timelike orbit.

Exactly Keplerian in Schwarzschild coordinates: \(\Omega = \sqrt{M/r^3}\).

Parameters:

r (float) – Orbital radius.

Returns:

float

circular_orbit_initial_state(r0)

Initial 8-component state for a circular, equatorial, timelike orbit at radius r0.

Parameters:

r0 (float) – Orbital radius. Should exceed isco_radius for the orbit to be stable.

Returns:

ndarray of shape (8,) – \((t, r, \theta, \phi, u^t, u^r, u^\theta, u^\phi)\) at \(\tau = 0\).

property critical_impact_parameter

Critical photon impact parameter, \(b_c = 3\sqrt{3}\,M\).

disk_redshift_factor(r, impact_parameter)

Combined gravitational + Doppler redshift factor for a circularly orbiting disk emitter.

For a photon with impact parameter \(b = L_{\text{ph}}/E_{\text{ph}}\) received from gas in a circular Keplerian orbit at radius \(r\), the ratio of observed to emitted photon energy is

\[g = \frac{E_{\text{obs}}}{E_{\text{emit}}} = \frac{1}{u^t(r)\left[1 - b\,\Omega(r)\right]}\]

where \(u^t(r) = 1/\sqrt{f(r) - r^2\Omega(r)^2}\) is the emitter’s time dilation factor and \(\Omega(r) = \sqrt{M/r^3}\) its orbital angular velocity (see circular_orbit_angular_velocity()). \(g < 1\) is redshifted (gravity, or the far/receding side of the disk); \(g > 1\) is blueshifted (the near/approaching side, when Doppler beaming wins).

Parameters:

• r (float or array_like) – Disk radius at emission.

• impact_parameter (float or array_like) – The received photon’s impact parameter \(b\).

Returns:

float or ndarray – The redshift factor \(g\).

eccentric_orbit_initial_state(r0, eccentricity_boost=0.0)

Initial state for a bound, generally eccentric equatorial timelike orbit.

Starts at radius r0 (taken as the initial apoapsis) with purely tangential velocity reduced below the circular value by eccentricity_boost, producing a radially oscillating (rosette, precessing) orbit for eccentricity_boost > 0.

Parameters:

• r0 (float) – Starting radius (apoapsis).

• eccentricity_boost (float, default=0.0) – Fractional reduction of the tangential velocity relative to the circular orbit at r0; 0 gives a circular orbit, larger values give more eccentric orbits.

Returns:

ndarray of shape (8,) – Initial state, as in circular_orbit_initial_state().

effective_potential(r, L)

Effective potential for equatorial timelike orbits, \(V_{\text{eff}}(r) = (1-2M/r)(1+L^2/r^2)\).

Radial motion obeys \((dr/d\tau)^2 = E^2 - V_{\text{eff}}(r)\).

Parameters:

• r (float or array_like) – Radial coordinate(s).

• L (float) – Specific angular momentum.

Returns:

float or ndarray

property horizon_radius

Event horizon radius, \(r_s = 2M\).

integrate_geodesic(y0, dtau, n_steps)

Integrate a geodesic from initial state y0 via fixed-step RK4.

Parameters:

• y0 (ndarray of shape (8,)) – Initial state, e.g. from circular_orbit_initial_state(), eccentric_orbit_initial_state(), or null_geodesic_initial_state().

• dtau (float) – Affine-parameter (or proper-time) step size.

• n_steps (int) – Maximum number of steps.

Returns:

dict of str -> ndarray – Keys "t", "r", "theta", "phi", "ut", "ur", "utheta", "uphi", each a 1D array over the (trimmed) trajectory, plus "energy" and "angular_momentum": the conserved \(E = -u_t = (1-2M/r) u^t\) and \(L = u_\phi = r^2 \sin^2\theta \, u^\phi\) evaluated at every step (constant along an exact geodesic; their spread measures integration error).

Examples

>>> bh = SchwarzschildBlackHole(M=1.0)
>>> y0 = bh.circular_orbit_initial_state(r0=10.0)
>>> traj = bh.integrate_geodesic(y0, dtau=0.01, n_steps=500)
>>> float(np.std(traj["energy"])) < 1e-6
True

property isco_radius

Innermost stable circular orbit radius, \(r_{\text{ISCO}} = 6M\).

light_deflection_angle(impact_parameter)

Leading-order light deflection angle for a photon passing at a given impact parameter.

\[\delta\phi \approx \frac{4M}{b}\]

The weak-field (\(b \gg M\)) limit of gravitational light bending; famously confirmed for starlight grazing the Sun during the 1919 solar eclipse.

Parameters:

impact_parameter (float) – Impact parameter \(b\).

Returns:

float – Deflection angle, in radians.

null_geodesic_initial_state(r0, impact_parameter, ingoing=True)

Initial state for an equatorial null (photon) geodesic with a given impact parameter.

Parameters:

• r0 (float) – Starting radius.

• impact_parameter (float) – Impact parameter \(b = L/E\).

• ingoing (bool, default=True) – Whether the photon starts moving toward decreasing \(r\).

Returns:

ndarray of shape (8,) – Initial state with \(E = 1\) (an affine-parameter normalization choice; only the ratio \(b\) is physical).

perihelion_precession(trajectory)

Measure the apsidal (perihelion) precession per orbit from an integrated trajectory.

Locates successive periapsis passages (local minima of \(r(\lambda)\)) and returns the excess azimuthal angle swept between them beyond \(2\pi\).

Parameters:

trajectory (dict of str -> ndarray) – Output of integrate_geodesic() for an eccentric orbit.

Returns:

ndarray – Precession (in radians) measured between each pair of consecutive periapsis passages found in the trajectory.

property photon_sphere_radius

Photon sphere radius, \(r_{ph} = 3M\).

shapiro_delay(r1, r2, impact_parameter)

Shapiro (gravitational) time delay for light passing between two radii.

\[\Delta t = 2M \ln\left[ \frac{\left(r_1 + \sqrt{r_1^2 - b^2}\right) \left(r_2 + \sqrt{r_2^2 - b^2}\right)}{b^2} \right]\]

The excess light-travel time (beyond flat-spacetime expectations) caused by spacetime curvature, first measured via radar ranging to Venus and the Viking Mars landers.

Parameters:

• r1 (float) – Radial coordinates of the two endpoints (e.g. Earth and a reflecting planet or spacecraft).

• r2 (float) – Radial coordinates of the two endpoints (e.g. Earth and a reflecting planet or spacecraft).

• impact_parameter (float) – Closest approach distance of the light path to the mass.

Returns:

float – Time delay, in geometrized time units (length; divide by physicskit.relativity.utils.constants.C_SI for seconds).

tidal_acceleration(r, proper_separation)

Tidal (geodesic deviation) acceleration on a radially infalling object.

In the local orthonormal frame of a freely-falling observer, the Riemann tensor’s nonzero components give a stretching acceleration along the radial direction and an equal compression along each transverse direction:

\[a_{\text{radial}} = \frac{2M}{r^3}\,\ell, \qquad a_{\text{transverse}} = -\frac{M}{r^3}\,\ell\]

for two points separated by a small proper length \(\ell\). This is “spaghettification”: an object falling feet-first is stretched head-to-toe and squeezed side-to-side, with both effects growing as \(1/r^3\) and formally diverging at the singularity.

Parameters:

• r (float or array_like) – Radial coordinate.

• proper_separation (float) – Proper length between the two points being torn apart (e.g. the height of an infalling astronaut).

Returns:

• radial (float or ndarray) – Stretching acceleration along the radial direction (positive: pulls the two points apart).

• transverse (float or ndarray) – Compressing acceleration along each transverse direction (negative: pushes together).

Examples

>>> bh = SchwarzschildBlackHole(M=1.0)
>>> radial, transverse = bh.tidal_acceleration(r=6.0, proper_separation=1.0)
>>> bool(radial > 0.0 and transverse < 0.0)
True

weak_field_precession_per_orbit(semi_major_axis, eccentricity)

Leading-order post-Newtonian perihelion precession per orbit.

\[\Delta\phi = \frac{6\pi M}{a (1 - e^2)}\]

Parameters:

• semi_major_axis (float) – Orbital semi-major axis \(a\).

• eccentricity (float) – Orbital eccentricity \(e \in [0, 1)\).

Returns:

float – Precession per orbit, in radians.