profiles.density.PlummerDensityProfile.compute_gravitational_potential#

Numerically compute the gravitational potential \(\Phi(r)\).

The potential of a spherical mass distribution can be expressed as

\[\Phi(r) = -G \left[ \frac{M(r)}{r} \;+\; 4 \pi \int_r^\infty \rho(r') \, r' \, dr' \right],\]

where \(M(r)\) is the enclosed mass within radius \(r\) and the second term accounts for contributions from shells at \(r' > r\).

This formulation ensures convergence even for profiles with infinite total mass (e.g., NFW), provided the density falls off sufficiently fast at large radii.

Parameters:

r (unyt_quantity or unyt_array) – Radius or array of radii at which to compute the potential. Must carry length units.
units (str or Unit, optional) – Desired output units for the potential. Defaults to cm**2/s**2 if not specified.
G (unyt_quantity, optional) – Gravitational constant to use. Defaults to unyt.physical_constants.gravitational_constant.
**kwargs – Additional arguments passed to compute_enclosed_mass() and the integration routine.

Returns:

phi – Gravitational potential at each radius.

Return type:

unyt_quantity or unyt_array

Notes

The potential is defined to vanish at infinity.
The boundary term \(M(r)/r\) is explicitly included.
Uses the internal pisces.math_utils.integration.integrate() utility for stable quadrature.

Examples

import numpy as np, unyt
from pisces.profiles.density import (
    EinastoDensityProfile,
)

profile = EinastoDensityProfile(
    rho_0=0.01 * unyt.Msun / unyt.kpc**3,
    r_s=20 * unyt.kpc,
    alpha=0.18,
)

radii = np.logspace(0, 2, 50) * unyt.kpc
phi = profile.compute_gravitational_potential(
    radii, units="km**2/s**2"
)
print(phi)

profiles.density.PlummerDensityProfile.compute_gravitational_potential#

This Page