Binary-inflated dynamical mass (B12)

A star cluster’s dynamical mass is read off its velocity dispersion, $M \propto \sigma^2 r_h / G$. But an unresolved binary contributes its orbital motion to the measured velocity, widening the dispersion and biasing the mass high. In low-dispersion systems — ultra-faint dwarfs, low-mass globular clusters — the binary motion is a large fraction of $\sigma$ itself, so the bias is large and has driven real debates about mass-to-light ratios and dark-matter content.

This demo (a) quantifies the bias, (b) shows a dispersion-only analysis cannot remove it — the $(\sigma_{\rm true}, f_b)$ problem is rank-1 degenerate — and (c) demonstrates a differentiable joint recovery from the non-Gaussian wings of the velocity distribution that returns an unbiased dynamical mass, with a Fisher/CRLB forecast vs sample size $N$ and RV precision $\epsilon$.

It is the kinematic companion to B4: B4 recovers the binary fraction $f_b$ photometrically (from the unresolved-binary mass function); B12 measures the dynamical mass kinematically in the presence of the same binaries. Both reuse the Moe & Di Stefano (2017) $P$–$q$–$e$ statistics Moe & Di Stefano, 2017 and the Tout et al. (1996) ZAMS mass–luminosity relation Tout et al., 1996.

The forward model¶

Each star has a cluster centre-of-mass line-of-sight velocity $v_{\rm COM}\sim\mathcal N(0,\sigma_{\rm true}^2)$. A fraction $f_b$ are unresolved binaries: you measure one blended velocity, not two. The blend sits at $v_{\rm obs}=v_{\rm COM}+\Delta$, where $\Delta$ is the internal orbital motion projected onto the line of sight and flux-weighted by the two components’ luminosities,

with $L_{\rm ZAMS}$ the Tout (1996) ZAMS relation (in-package via progenax.stellar). Two limits make $\Delta$ intuitive:

• $q\to1$ (equal mass): $L_1=L_2$ and $v_1=-v_2$, so the blend cancels, $\Delta\to0$. Twins are nearly invisible contaminants.

• $q\to0$ (extreme ratio): the faint secondary adds no light, so you see the primary’s reflex, $\Delta\to v_{1,\rm los}$ (an SB1).

The orbits $(P, q, e)$ are drawn from the Moe coupling; periods become semimajor axes via Kepler’s third law, and the orbit is sampled at a uniform mean anomaly (uniform in time, so eccentric orbits are correctly weighted toward slow apocenter passages). The distribution of $\Delta$ over a large binary pool is the contamination kernel $K_{\rm orb}$ — crucially independent of $\sigma_{\rm true}$, so the observed velocity density factorizes,

where $\epsilon$ is the per-star RV precision (added in quadrature). Because $K_{\rm orb}$ does not depend on $\sigma_{\rm true}$, it is precomputed once. The expected binned counts $\mu_k(\sigma_{\rm true},f_b)=N\!\int_{{\rm bin}\,k}p$ are differentiable in both parameters, and $(\sigma_{\rm true},f_b)$ are recovered by a per-bin Poisson MLE (scripts/_demo_inference.py), with the Fisher from the Poisson information.

Inputs and assumptions¶

The fit recovers only two parameters, $(\sigma_{\rm true}, f_b)$; every other quantity is an assumed-known input. Making that split explicit is essential — the demo’s optimism lives entirely in what it treats as known.

Result — freshly run, ALL GATES PASS¶

UFD-like regime: $\sigma_{\rm true}=5$ km/s, $f_b=0.5$, $N=1500$ RV stars, $\epsilon=1$ km/s, metal-poor $Z=10^{-3}$, $r_h=30$ pc (exit 0; wall $\approx 90$ s).

The bias, and why dispersion-only cannot fix it¶

The dispersion grows by the variance budget $\sigma_{\rm obs}^2 = \sigma_{\rm true}^2 + f_b\,\mathrm{Var}(K_{\rm orb}) + \epsilon^2$, so the naive virial mass is biased high by exactly $(\sigma_{\rm obs}/\sigma_{\rm true})^2$ — $1.28\times$ at $f_b=0.5$ (Figure 2). Crucially, this is a systematic: averaging more stars at fixed $f_b$ measures the wrong mass more precisely, it does not remove the bias.

A dispersion-only analysis sees a single number $\sigma_{\rm obs}$, so any $(\sigma_{\rm true}, f_b)$ on the curve $\sigma_{\rm true}^2 + f_b\,\mathrm{Var}(K_{\rm orb}) = \text{const}$ fits equally well — an infinite ridge (Figure 3, orange). Its $2\times2$ Fisher is the outer product of one gradient: exactly rank-1 (one zero eigenvalue, cond $=2.3\times10^{16}$). One number cannot separate two unknowns.

How the wings break the degeneracy¶

The fix is to fit the whole shape of $p(v)$, not just its width. The Gaussian core constrains one combination of $(\sigma_{\rm true},f_b)$; the non-Gaussian wings — the heavy $K_{\rm orb}$ tail no single-star Gaussian can mimic (Figure 1) — constrain a different combination. Two independent constraints turn the open ridge into a closed ellipse (Figure 3, blue): the full-distribution Fisher is full-rank (cond $=878$, a 10¹³ improvement). The recovered mass is unbiased, $\sigma=4.97\pm0.13$ km/s, $M_{\rm ratio}=0.99$ vs the naive 1.28 (verified unbiased over 24 realizations: $\langle\sigma\rangle=5.00$, $\langle f_b\rangle=0.50$).

The RV-precision floor¶

Because $\epsilon$ is known and folded into (2), the recovered mass stays unbiased at every $\epsilon$ (Figure 4) — even at $\epsilon=5$ km/s $\approx\sigma_{\rm obs}$. What degrades is precision: $\sigma(\sigma_{\rm true})$ doubles as $\epsilon$ washes out the wing signature, and the detectable binary fraction $f_b\,P(|\Delta|>\epsilon)$ collapses from 0.21 to 0.02. The honest scope: you can only correct for the binaries you can detect, and that detectable fraction shrinks as the spectrograph gets noisier — but the correction you do make is unbiased. The Fisher forecast (Figure 5) shows both uncertainties fall as $N^{-1/2}$; $f_b$ (and hence the bias correction) is the harder parameter, since the wings it lives in are sparsely populated.

Figures¶

Caveats¶

How to run¶

env -u VIRTUAL_ENV uv run --no-sync python scripts/demo_binary_dynamical_mass.py

References¶

The binary $P$–$q$–$e$ coupling is Moe & Di Stefano (2017); the ZAMS mass–luminosity relation is Tout et al. (1996). The Moe and companion models are documented on the binary statistics theory pages; the photometric sibling is B4.