Plummer velocity distribution functions

The isotropic equilibrium DF for the Plummer profile is one of the oldest closed-form distribution functions in stellar dynamics. It falls out of Plummer (1911)’s original analysis as a direct consequence of the smooth $(1 + r^2/a^2)^{-5/2}$ density profile, and provides the canonical example of how the Eddington inversion machinery turns a 3D density profile into a 6D phase-space DF for spherical isotropic systems.

This chapter derives $f(E) \propto (-E)^{7/2}$ from the Plummer density profile, lists the closed-form moments, walks through progenax’s sampling implementation, and connects to the anisotropy and rotation extensions in Anisotropy and rotation.

Eddington inversion in one paragraph¶

For a spherical, isotropic system, the DF depends only on the binding energy $\mathcal{E} \equiv -E = -\tfrac{1}{2}v^2 + \Psi(r)$, where $\Psi \equiv -\Phi$ is the relative potential (positive for a bound system, $\Psi \to 0$ at infinity). The Eddington formula

↗ model card

inverts the spatial-density profile $\rho(\Psi)$ — the density expressed as a function of the relative potential — into the DF $f(\mathcal{E})$. (An earlier revision of this page wrote the integration variable as $\Phi$; the integral runs over the positive relative potential, so $\Psi$ is the honest symbol.) The boundary term vanishes for systems with $\rho \to 0$ as $\Psi \to 0$ (i.e. cluster density vanishing at infinity, which Plummer satisfies). For Plummer the full inversion gives a closed-form DF in $\mathcal{E}$ alone.

The Plummer DF¶

Substituting Plummer’s

into (1). The boundary term vanishes ($\rho \propto \Psi^5 \to 0$ as $\Psi \to 0$), and the second derivative is a pure power: $\mathrm{d}^2\rho/\mathrm{d}\Psi^2 = \tfrac{3M}{4\pi a^3}\, (a/GM)^5\, 20\,\Psi^3$. The Abel integral is then a Beta function,

and collecting constants ($\tfrac{1}{\sqrt8 \pi^2} \cdot \tfrac{3}{4\pi} \cdot 20 \cdot \tfrac{32}{35} = \tfrac{24\sqrt2}{7\pi^3}$) yields

↗ model card

with $\mathcal{E} = -E = GM/\sqrt{r^2 + a^2} - \tfrac{1}{2}v^2$. The $\mathcal{E}^{7/2}$ exponent is the signature of the Plummer profile — it follows directly from $\rho \propto \Phi^5$ and the Abel-integral $\frac{1}{2}$-power kernel of (1). progenax’s sampler uses this exact form, with no truncation or smoothing.

Closed-form moments¶

The first three velocity moments of the Plummer DF are

with $\langle v^2 \rangle = \sigma_r^2 + 2\sigma_t^2 = 3\sigma_r^2$. The ratio $\langle v^2 \rangle / |\Phi(r)| = \tfrac{1}{2}$ at all radii — this is the local statement of the virial theorem and the proximate reason a Plummer cluster sampled from (4) lands at $Q_{\mathrm{vir}} = 0.5$ at the DF level — realised as 0.5 up to the finite-$N$ Monte-Carlo fluctuation, with no rescale (What is an initial condition?).

The central velocity dispersion is

which (combined with (3)) gives a one-line conversion between cluster mass, scale radius, and central velocity dispersion — useful for sanity-checking observational fits.

Sampling: speed and rejection¶

progenax samples velocities from (4) via inverse-CDF on the speed distribution at each particle’s position. At fixed $r$, the DF integrated over angles gives

↗ model card

with $v \in [0, v_{\mathrm{esc}}(r)]$. The distribution is bounded, unimodal, and smooth, so progenax uses a fixed-iteration inverse-CDF lookup table rather than rejection sampling. The lookup table is constructed once per call (vectorised over $r$) and inverted via jnp.interp, so the cost per particle is dominated by the single bisection step — no Python loop, no rejection retries, fully JIT-compatible.

Implementation in progenax¶

from progenax.kinematics import PlummerVelocityDF
from jaxstro.units import STELLAR
import jax, jax.numpy as jnp

df = PlummerVelocityDF(r_h=1.0)        # Half-mass radius in pc
masses = jnp.ones(1000)
key = jax.random.PRNGKey(42)
key_pos, key_vel = jax.random.split(key)   # never reuse a key

# Positions from the matched profile (any sampled positions work)
from progenax.profiles import PlummerProfile
positions = PlummerProfile(r_h=1.0).sample_positions(masses, key_pos)

# Sample velocities consistent with f(E) ∝ E^(7/2)
velocities = df.sample_velocities(positions, masses, key_vel, G=STELLAR.G)

# The DF is an exact equilibrium, so Q_vir -> 0.5 with no rescale
# (finite-N: 0.5 +/- ~5e-3 at N=1e4, the Monte-Carlo fluctuation)

PlummerVelocityDF is an Equinox module. Its single Python parameter is r_h; everything else (the central potential depth, the inverse-CDF lookup table) is derived analytically from $r_h$ inside the class. Differentiability flows through $r_h \to \Phi_0 \to f(\mathcal{E}) \to \mathbf{v}$.

Sanity checks¶

The Plummer DF satisfies four checks that progenax’s validation suite (Plummer equilibrium) verifies on every release:

The first row implements the full Jeans-equation check: Plummer positions sampled from PlummerProfile plus velocities from PlummerVelocityDF produce $Q_{\mathrm{vir}} = 0.5$ to 0.5%, with the residual being the finite-$N$ fluctuation of the Monte Carlo energy estimate. No virial rescaling is needed — equilibrium is exact at the DF level.

Check yourself¶

Anisotropy and rotation extensions¶

The Plummer DF above is isotropic ($\beta = 0$) and non-rotating. progenax extends it with two compositions:

• Osipkov-Merritt anisotropy: introduces a free anisotropy radius $r_a$ such that $\beta(r) = r^2/(r^2 + r_a^2)$. The DF is no longer a function of $\mathcal{E}$ alone but of $Q = \mathcal{E} - L^2/(2 r_a^2)$, the augmented integral of motion. progenax implements this as an intrinsic DF option, PlummerVelocityDF(r_h=..., anisotropy_radius=r_a), using the analytic OM Plummer DF (Merritt 1985, Eq. 45), which is non-negative only for $r_a \ge 0.75\,a$ (Merritt 1985, Eq. 46). See Anisotropy and rotation.

• Solid-body / differential rotation: adds a tangential velocity component $\mathbf{v}_\phi(\mathbf{r}) = \mathbf{\Omega}(r) \times \mathbf{r}$ to the isotropic Plummer velocities. Implemented via apply_solid_body_rotation(df, omega) and apply_differential_rotation(df, omega_profile).

Both extensions are documented in detail at Anisotropy and rotation. They compose freely with each other and with non-Plummer profiles (e.g. King + Osipkov-Merritt is a common combination for old globular clusters with observed radial anisotropy).

Domain of validity and limitations¶

1. Single-mass equilibrium only. The Eddington DF assumes one stellar mass species. Multi-mass equilibrium DFs require simultaneous treatment of all mass groups — the multi-mass generalisation is formalized by Gieles & Zocchi (2015), which progenax plans to implement natively as its own differentiable generalization (see The differentiable lowered-model family (Engine A)).

2. Spherical and isotropic by default. Anisotropy and rotation extensions exist (above) but are layered on top of the isotropic DF; they do not change the underlying $f(\mathcal{E})$.

3. Untruncated. The DF extends to $r \to \infty$. For tidally truncated Plummer-like systems use King profile instead, or apply post-sampling truncation via apply_tidal_truncation (Tidal physics).

Implementation, validation & references¶

• In code: src/progenax/kinematics/plummer_df.py (the isotropic $f(\mathcal{E})$, the inverse-CDF speed sampler, and the optional Osipkov–Merritt anisotropy) — see the PlummerVelocityDF API.

• Validated in: Plummer equilibrium — the regression suite that locks the velocity moments, isotropy $|\beta| < 0.02$, and the unscaled $Q_{\mathrm{vir}} = 0.5$.

• Primary sources: Plummer (1911) (the model and its DF); modern textbook treatments follow Binney & Tremaine Galactic Dynamics §4, and Aarseth et al. (1974) is a clean N-body initialisation reference — full notes in the bibliography. The mass-first sampling contract follows the IC redesign at IC redesign spec (2026-02-12).