Velocity distribution functions

A spatial density profile alone does not produce a complete IC: every particle needs a velocity. The velocity distribution function (DF) specifies how velocities are sampled given particle positions and masses. progenax implements the canonical equilibrium DFs for its spatial profiles (Plummer, King, EFF), the self-consistent radially anisotropic Michie–King DF, plus extensions for Osipkov–Merritt radial anisotropy (Merritt (1985)) and rigid/differential rotation.

What “equilibrium” means¶

A spatial profile $\rho(r)$ paired with the wrong velocity DF produces a non-equilibrium IC: the cluster will dynamically relax towards a new equilibrium state via violent relaxation on the first crossing time. For most production work this is undesirable — you wanted a Plummer cluster at $t = 0$, not a Plummer-like distribution heading toward something else. The fix is to use the matched equilibrium DF.

with $\beta = 1 - \sigma_t^2/\sigma_r^2$ the anisotropy parameter (zero for isotropic). progenax’s validation suite (Validation methodology) checks this directly for sampled positions+velocities.

```{figure} ../figures/phase_space_hexbin.webp
:label: fig-phase-space-hexbin
:width: 100%

The DF as a density on the classic $(r, v_r)$ phase-space plane
($4\times10^5$ stars per model, seed 13), with the analytic escape envelope
$\pm v_{\rm esc}(r) = \pm\sqrt{2\Psi(r)}$ computed three different ways
(closed form / ODE $\psi$-grid / direct density quadrature). Plummer's
envelope never closes — and its stars stop visibly *inside* it, because
$f \propto \mathcal{E}^{7/2}$ vanishes toward the escape energy; King hugs
its envelope and pinches shut at $r_t$; EFF fills its power-law halo out to
the hard edge. No star crosses an envelope: bound-ness, visually. Regenerate: `python -m laboratory.icviz --only phase-space-hexbin`.

When you actually want a non-equilibrium IC¶

Three deliberate non-equilibrium configurations are useful enough to deserve mention:

1. Subvirial / cold cluster, $Q_{\mathrm{vir}} \approx 0.3$ — the Allison et al. (2009) cool-fractal setup. Cluster contracts on a crossing time and produces dynamical mass segregation within $\sim 1$ Myr.

2. Supervirial / hot cluster, $Q_{\mathrm{vir}} \approx 0.75$ — post-gas-expulsion. Cluster expands, may unbind partially.

3. Mismatched-DF cluster — Plummer profile with a King DF, or vice versa. Useful for testing the rate of violent relaxation without contaminating the experiment with non-equilibrium energy ratio.

In all three cases, the spatial profile and base DF are sampled independently, then the velocity field is rescaled to the target $Q_{\mathrm{vir}}$ via Virial Q convention (Q = T/|V|). Each non-equilibrium configuration is a deliberate scientific choice; progenax exposes the parameters so the configuration can be made explicit at IC construction time.

Common API contract¶

Every velocity DF satisfies the VelocityDF protocol, whose single required method is sample_velocities:

class VelocityDF(Protocol):
    def sample_velocities(
        self,
        positions: Float[Array, "N 3"],
        masses: Float[Array, "N"],
        key: PRNGKey,
        G: float,
    ) -> Float[Array, "N 3"]:
        """Draw N velocities consistent with the DF given positions+masses."""
        ...

That is the whole contract — there is no velocity_dispersion method on the protocol or on the concrete DF classes; the dispersion profiles shown in the per-DF chapters are derived analytically (Plummer) or built internally by the sampler (King, EFF, Michie) and are not part of the public surface. sample_velocities is differentiable in the model parameters via the chain (e.g. $r_h \to \sigma_r(r) \to \mathbf{v}$), and every DF is JIT-compatible and vectorisable via jax.vmap.

Multiple components in one potential¶

Everything above is single-component: one DF in its own self-consistent potential. For clusters with multiple populations sharing one potential — GC 1G/2G, halo+core, multi-mass equipartition — the right tool is MultiComponentCluster (Multi-component populations), which provides two equilibrium engines: the DF-defined lowered-isothermal family (Engine A) and density-defined Eddington inversion in a shared potential (Engine B). Per-component equilibrium ($Q_j = 0.5$) emerges from the DFs with no external rescale.

Implementation, validation & references¶

• In code: the velocity DFs live under src/progenax/kinematics/ (plummer_df.py, king_df.py, eff_df.py, michie_df.py, plus rotation.py and the Osipkov–Merritt anisotropy options). See the kinematics API; the per-DF chapters (Plummer, King, Michie–King, anisotropy & rotation) carry the exact module paths.

• Validated in: Plummer equilibrium, King profile, EFF profile, Michie anisotropy, and OM / rotation anisotropy.

• Primary sources: the Eddington-inversion machinery (Plummer (1911), King (1966)) is standard textbook material; the anisotropy extensions are Merritt (1985) (Osipkov–Merritt) and Michie (1963) (self-consistent King), and Gieles & Zocchi (2015) covers the multi-mass lowered-model family. Full notes in the bibliography.