The differentiable lowered-model family (Engine A)

The idea¶

The classical truncated-isothermal models — Woolley (no energy truncation), King (1966) (a linear lowered Maxwellian), and Wilson (a quadratic lowered form) — are members of one family distinguished by how sharply the distribution function is truncated near the escape energy. Gieles & Zocchi (2015) made this explicit with a single continuous truncation parameter $g$: their LIMEPY DF reads

↗ model card

with $g=0$ recovering the Woolley cutoff, $g=1$ the King model, and $g=2$ the Wilson model, interpolating smoothly in between. The same framework extends naturally to multiple mass components (each mass group its own $\sigma_j$, sharing one potential) and to radial anisotropy (an Osipkov–Merritt / Michie term).

The coupled multi-component equilibrium¶

progenax generalizes the single-mass solve to $n$ components sharing one dimensionless potential $\psi(\xi)$, $\xi = r/r_c$. The model-defining insight is that a component’s single free Engine-A scale is its velocity-scale ratio

so component $j$ sees the shared potential at the rescaled depth $W_j(\xi) = \texttt{rescale}_j\,\psi(\xi)$ and contributes a density $\hat\rho_j \propto E_\gamma\!\big(g+\tfrac32,\, \texttt{rescale}_j\,\psi\big)$ (the “lowered exponential” of Gieles & Zocchi (2015), Appendix B), giving the coupled Poisson equation

↗ model card

with $\alpha_j$ the central density fractions ($\sum_j\alpha_j = 1$). A colder component ($w_j < 1$) feels a deeper effective well and concentrates — as a genuine equilibrium, not a post-hoc radial reshuffling. Mass segregation is then the equipartition convenience

↗ model card

(Gieles & Zocchi (2015), Eqs. 24–26): $\delta = 1/2$ is the standard partial-equipartition choice, $\delta = 0$ collapses every component to the single-mass model exactly (the cleanest oracle). The representative stellar masses $m_j$ are otherwise decoupled labels — the structure depends only on $(\alpha_j, w_j, \texttt{ra\_hat}_j, W_0, g)$, so equal-mass populations of different concentration (GC 1G/2G, halo+core, binaries-vs-singles) set $w_j$ directly.

Per-component radial anisotropy enters as an Osipkov–Merritt / Michie term with its own anisotropy radius $\hat r_{a,j} = r_{a,j}/r_c$ (Michie (1963); Merritt (1985)); the anisotropic density is evaluated by an exact, numerically stable 1-D quadrature (the bounded $T(\beta)$ Poisson sum) because JAX’s hyp1f1 is NaN for arguments $\gtrsim 100$.

Entry points¶

The family is exposed through three MultiComponentCluster constructors (see Multi-component populations for the engine map and Two-component populations — worked examples for a worked example):

from progenax import MultiComponentCluster

# General: per-component velocity-scale ratios (GC 1G/2G, halo+core, ...)
model = MultiComponentCluster.from_components(
    alpha_j=[0.5, 0.5], w_j=[0.7, 1.0], m_j=[1.0, 1.0], W0=7.0, g=1.0)

# Mass segregation as equipartition: w_j = mu_j^(-delta)
model = MultiComponentCluster.from_mass_segregation(
    alpha_j=[0.6, 0.4], m_j=[0.3, 3.0], W0=7.0, g=1.0, delta=0.5)

# IMF-driven: bin the IMF, eigenvalue-solve for the alpha_j that
# reproduce the per-bin mass budget
model = MultiComponentCluster.from_imf(imf, n_comp=8, W0=7.0, g=1.0, delta=0.5)

ic = model.sample_cluster(key, n_stars=50_000)   # ICResult with component_id

sample_cluster draws every star from its component’s lowered DF at the rescaled potential $W_j(r)$ and velocity scale $s_j = s\,w_j$ — no external virial rescale anywhere: each component is individually virial ($Q_j = 0.5$) and the sampled cluster is globally virial for any mass spectrum. The unequal-mass two-population regression measures theory $Q_j = 0.5 \pm 2\times 10^{-3}$, sampled global $Q$ within $\pm 0.04$ and per-component $Q_j$ within $\pm 0.07$ at test resolution. The single-mass corners were validated against the released models directly: $g=1$ isotropic $\equiv$ King, $g=1$ anisotropic $\equiv$ Michie, both with $Q = 0.5$ unscaled.

Check yourself¶

The DF-table performance layer¶

The anisotropic density quadrature is the 86% hotspot of the coupled solve, and per-star speed draws dominated sampling. Phase 1.5 replaced both with three differentiable table primitives in profiles/limepy_tables.py — with the exact quadrature retained everywhere as a selectable oracle (aniso_method="quadrature"), and every approximation budget asserted in tests:

Measured speedups (warm, two-component, $W_0 = 7$): anisotropic construction 5.6× faster (957 → 170 ms); sampling at $N = 10^5$: isotropic 67× (0.48 µs/star), anisotropic 21.7× (2.9 µs/star). The anisotropic table build dominates at small $N$ — break-even is ~3k stars. Crucially, the equilibrium oracle component_virial_ratios is deliberately quadrature-only (it must not share the approximation it checks): a table-built model proves $Q_j = 0.5001 \pm 1.5\times 10^{-4}$ against it.

Why progenax reimplemented it (rather than depend on LIMEPY)¶

Gieles & Zocchi (2015) is the reference formalism, and the published limepy code is the standard numpy/scipy implementation. progenax does not wrap or depend on it, for one decisive reason: limepy is not differentiable. The entire progenax thesis is JAX-native, end-to-end-differentiable initial conditions, so that structural parameters can be inferred from data by gradient descent or HMC (see Differentiability rules).

The implemented family delivers exactly this: a single continuous parameter vector — $(W_0, g, \{w_j\}, \delta, r_a)$ — can be fit jointly to an observed cluster, with $g$ itself a fitted quantity that selects the truncation sharpness the data prefer (King-vs-Wilson as a posterior, not a modelling choice). $\partial(\text{model})/\partial g$ flows through the coupled Poisson solve and through the DF tables: AD agrees with finite differences to $2.15\times 10^{-4}$, and the table-backed gradient agrees with the quadrature-backed gradient to $2.1\times 10^{-4}$. The anisotropic and mass-segregation parameters $(r_a, \eta, \delta)$ are likewise differentiable through construction and sampling.

Delivered scope¶

Relationship to the other roadmap item¶

This is independent of, but complementary to, the now-resolved differentiable tidal radius $\partial r_t/\partial W_0$ (Differentiable King tidal radius ($\partial r_t/\partial W_0$)): the scalar $r_t$ already carries a finite, exact gradient through the unclamped $\psi=0$ crossing. In the unified family $r_t$ is a function of $(g, W_0)$; the same implicit-function-theorem treatment of the $\psi=0$ crossing carries over to the $g$ generalization unchanged.

Implementation, validation & references¶

• In code: src/progenax/cluster/multicomponent.py (MultiComponentCluster and the coupled-Poisson core solve_multicomponent_limepy), with the lowered-density kernel in src/progenax/profiles/limepy_multimass.py and the DF-table layer in src/progenax/profiles/limepy_tables.py. See the MultiComponentCluster API.

• Validated in: multimass equilibrium — per-component $Q_j$, the single-mass King/Michie corners, and the AD-vs-FD gradient checks for $(W_0, g, w_j, \delta, r_a)$.

• Primary sources: Gieles & Zocchi (2015) (the LIMEPY lowered-model formalism); single-model members are King (1966) (with Woolley 1954 / Wilson 1975 for the $g = 0, 2$ endpoints), and anisotropy follows Michie (1963) and Merritt (1985). Full notes in the bibliography: Gieles & Zocchi 2015 and King 1966.