Binary-aware IMF recovery

The stellar initial mass function (IMF) governs chemical enrichment, supernova rates, and the integrated light of galaxies. Measuring its shape — and whether it varies with environment — is one of the central problems in stellar astrophysics Salpeter, 1955Kroupa, 2001Chabrier, 2003Marks et al., 2012Jeřábková et al., 2018. Most stars do not form alone: roughly 50% of solar-type stars and >90% of O-type stars have at least one bound companion Sana et al., 2012Moe & Di Stefano, 2017. When those companions are unresolved, the IMF inferred from observed system masses is systematically biased toward shallower slopes.

This chapter is the single source of truth for the binary-aware IMF recovery pipeline implemented in progenax.imf.BinaryIMF and exercised in scripts/validate_binaries.py (fast MLE) and validation/imf/validate_binary_aware_recovery.py (full numpyro NUTS). It derives the system mass function under the Moe & Di Stefano (2017) joint binary statistics, formulates a likelihood that marginalises over latent binary status and mass ratio, and quantifies the regime — sample size $N \gtrsim 10^4$ — where the naive single-star likelihood becomes “confidently wrong”: its 95% credible interval shrinks below the bias and excludes the true slope.

The Maschberger (2013) functional form¶

progenax uses the Maschberger (2013) “L3 distribution” as the single-star IMF backbone:

with default parameters $\alpha = 2.3$ (Salpeter high-mass slope), $\beta = 1.4$, $\mu = 0.2\,\Msun$. The high-mass limit recovers $f \propto m^{-\alpha}$; the second factor produces the observed turnover below $\sim 0.3\,\Msun$. Crucially, the cumulative distribution has a closed-form inverse, enabling exact analytical sampling — no numerical root-finding, fully jax.lax.scan-compatible, fully differentiable. See `progenax.imf` for the implementation in progenax.imf.Maschberger.

The binary contamination problem¶

For an unresolved binary system, the observer measures the system mass

where $q = m_2/m_1 \in [q_{\min}, 1]$ and $f_b(m_1)$ is the mass-dependent binary fraction. Binaries transfer probability mass from intermediate to high system masses — a primary $m_1 = 5\,\Msun$ with $q = 0.5$ appears as a single $7.5\,\Msun$ system. Consequently the observed system mass function is shallower than the true single-star IMF.

The naive practice of fitting a single-star IMF to the observed system masses — likelihood $\mathcal{L}_{\mathrm{naive}}(\alpha) = \prod_i \xi(m_{\mathrm{sys},i} \mid \alpha)$ — converges as $N\to\infty$ to the wrong value of $\alpha$. Kroupa (2001) quantified the bias at $\Delta\alpha \sim 0.05$–0.10 depending on mass range and binary fraction. Section Paragraph below shows that the bias does not shrink with $N$, and at LSST-scale samples ($N \sim 10^{4}$–10⁵) the naive 95% credible interval becomes narrower than the bias and excludes the truth.

The Moe & Di Stefano (2017) joint binary statistics¶

Moe & Di Stefano (2017) compiled the most comprehensive census of binary star properties to date, combining spectroscopy, eclipsing binaries, long-baseline interferometry, adaptive optics, and common proper motion surveys. After correcting each sample for its selection function, they derived the intrinsic joint distribution $f(M_1, q, P, e)$ of primary mass, mass ratio, orbital period, and eccentricity. Two findings shape progenax’s implementation:

1. Binary properties are not universal. Binary fraction, mass-ratio distribution, and period distribution all vary systematically with $M_1$. Treating them as constants introduces correlated errors.

2. The joint distribution is not separable: $f(M_1, q, P, e) \neq f(M_1)\,f(q)\,f(P)\,f(e)$. The mass-ratio distribution depends on both $M_1$ and $P$; eccentricity depends on $P$; companion frequency depends on both.

Binary fraction $f_b(M_1)$¶

The mass-dependent binary fraction $f_b(M_1) = 1 - \mathcal{F}_{n=0}$ — the multiplicity fraction ($\le 1$, the probability of $\ge 1$ companion) — is tabulated, with its full per-mass values and provenance (Moe & Di Stefano (2017) Table 13 single-star row for $M_1 \ge 0.8\,\Msun$, M-dwarf surveys below), in the canonical multiplicity-statistics chapter. The trend is steep: O-type stars are almost always in multiples ($f_b \approx 0.9$), while M-dwarfs are mostly single ($f_b \approx 0.2$–0.3) — a factor-of-4 dynamic range that is the leading-order driver of the binary-IMF bias.

This $f_b$ is the multiplicity fraction, not the companion frequency $f_{\mathrm{mult}}$, which exceeds 1.0 ($\sim 2.1$ per O-type primary, Moe & Di Stefano (2017) Table 13) once hierarchical triples are counted — see Limitations below for the single-companion approximation progenax adopts.

Mass-ratio distribution $g(q \mid M_1)$¶

For a primary of mass $M_1$, Moe & Di Stefano (2017) parameterise $g(q \mid M_1)$ as a power law plus a narrow “twin excess” peak near $q = 1$:

where $Z_{\mathrm{pl}} = \int_{q_{\min}}^{1} q^\gamma\,\mathrm{d}q$ and $\sigma_{\mathrm{twin}} \approx 0.03$. progenax uses a period-averaged single-slope $\gamma(M_1)$ — a reduction of Moe & Di Stefano (2017) Table 13’s period-dependent $\gamma_{\mathrm{smallq}}/\gamma_{\mathrm{largeq}}$ (averaged over $\log P$) — appropriate for total mass-function analyses, but losing the small-$q$ structure that matters for surveys sensitive to specific period ranges (see Limitations). The per-mass $\gamma(M_1)$ and $f_{\mathrm{twin}}(M_1)$ tables live in the canonical Mass-ratio distributions.

Qualitatively, the slope trend reverses at high masses: massive stars preferentially have low-q companions ($\gamma < 0$), while low-mass stars prefer near-equal-mass companions ($\gamma > 0$). The twin fraction $f_{\mathrm{twin}}$ peaks at solar-type primaries ($f_{\mathrm{twin}} = 0.10$ for 0.8–$1.2\,\Msun$) and falls to 0.03 for $M_1 > 3.5\,\Msun$.

The system mass function¶

Given the forward generative model — draw $m_1$ from the IMF, decide binary or not, add $m_2$ if binary — the probability density of the observed system mass $M$ is

The binary term integrates over every primary mass $m_1$ that could produce the observed system mass $M$, weighted by the IMF, the binary fraction, and the mass-ratio distribution. The integration limits come from $q_{\min} \le q \le 1$: $q \ge q_{\min}$ implies $m_1 \le M/(1+q_{\min})$; $q \le 1$ implies $m_1 \ge M/2$. The Jacobian $1/m_1$ arises from the change of variables $q = M/m_1 - 1$.

This integral has no closed form because $f_b(m_1)$ and $g(q \mid m_1)$ are piecewise functions of $M_1$. progenax evaluates it via 128-point Gauss-Legendre quadrature, which gives ~14 digits of accuracy for smooth integrands.

The binary-aware likelihood¶

For an observed system mass $M$:

where $\bar{f}_b$ is evaluated at the candidate primary mass $m_1 = M/(1+q)$ inside the integral. The integrand is smooth in $q$, so Gauss-Legendre quadrature on $[q_{\min}, 1]$ converges fast.

Computational cost. For $N$ observed systems and 128 quadrature points, each likelihood evaluation requires $N \times 128$ IMF + mass-ratio evaluations. At $N = 30{,}000$ (LSST-scale) this is $\sim 3.8 \times 10^6$ evaluations per NUTS step, $\sim 6 \times 10^9$ flops per chain. On a MacBook Pro CPU one chain takes $\sim 35$ minutes; on an A100 the same $N \times 128$ tile is embarrassingly parallel via jax.vmap, bringing the per-chain time down to $\sim 30$ s — an estimated $\sim 70\times$ speedup (see Binary-aware likelihood for the per-step breakdown).

Confidently wrong: the naive likelihood at large N¶

The naive likelihood $\mathcal{L}_{\mathrm{naive}}(\alpha) = \prod_i \xi (M_i \mid \alpha)$ is a misspecified model: the data-generating process includes binaries; the likelihood does not. As $N\to\infty$ both naive and binary-aware posteriors shrink as $\sigma \propto 1/\sqrt{N}$ (Bernstein–von Mises), but they shrink around different values:

• The binary-aware posterior shrinks around $\alpha_{\mathrm{true}}$.

• The naive posterior shrinks around $\alpha_{\mathrm{true}} + \Delta\alpha$, where $|\Delta\alpha| \sim 0.05$–0.10 is set by the binary statistics and is independent of $N$.

The 95% credible interval shrinks; the bias does not. At some $N$ the CI becomes narrower than the bias, and the posterior excludes the true value:

This regime — narrow posterior centred on the wrong answer — is the central danger of ignoring binaries with large datasets. progenax’s binary-aware likelihood eliminates the bias by construction; the ${\sim}128\times$ cost of the inner $q$-quadrature is the price.

Observation operators: what does the telescope measure?¶

The likelihood (5) uses the dynamical-mass operator $m_{\mathrm{sys}} = m_1 + m_2$. This is exact for eclipsing binaries, SB2 spectroscopic binaries, and astrometric binaries — observations that recover the true system mass. It overstates the distortion for photometric surveys, where the telescope measures combined flux rather than combined mass.

For a main-sequence luminosity scaling $L \propto m^s$ with $s \approx 3$–4, the inferred photometric mass of an unresolved binary is

Because $s \gtrsim 3$, the secondary contributes appreciably to the flux only when $q \to 1$:

The mass-addition operator $m_1 + m_2$ therefore gives an upper bound on the binary distortion. Photometric masses are much less affected, especially at moderate $q$ — but the bias does not vanish, and the naive model still becomes confidently wrong at sufficiently large $N$. Full LSST-realistic inference replaces $m_1 + m_2$ with isochrone-mediated flux addition; see the Add a binary population recipe for the planned extension.

Limitations¶

Implementation in progenax¶

The pipeline lives in three modules:

• progenax.imf.BinaryIMF — the Moe & Di Stefano (2017)-conditioned forward model and likelihood. Differentiable via jax.grad, JIT-compatible, vectorisable via jax.vmap.

• progenax.binaries.MoeEccentricity and friends — the joint binary statistics tables used by BinaryIMF (see Binary period distributions).

• scripts/validate_binaries.py — the fast (MLE, no-MCMC) figure suite that regenerates the “$N$ vs. CI vs. bias” result; validation/imf/validate_binary_aware_recovery.py is the full numpyro NUTS cross-check.

See `progenax.imf` for full API signatures and Add a binary population for a recipe building a binary-corrected mock catalog end-to-end.

Implementation, validation & references¶

• In code: src/progenax/imf/binary/imf.py (BinaryIMF — the forward model and sampling helpers), with the joint statistics in src/progenax/imf/binary/moe_di_stefano.py and the single-star backbone in src/progenax/imf/smooth.py (Maschberger) — see the IMF API. The recovery figures are regenerated by scripts/validate_binaries.py (fast MLE) and validation/imf/validate_binary_aware_recovery.py (full numpyro NUTS).

• Validated in: binary-aware recovery — the regression suite that demonstrates unbiased $\alpha$ recovery at $N = 30{,}000$ where the naive fit is confidently wrong.

• Primary sources: the single-star backbone is Maschberger (2013); the binary-contamination statistics are Sana et al. (2012)Moe & Di Stefano (2017)Moe et al. (2019); the “confidently wrong” framing is original to progenax. Full notes in the bibliography. For the environment-dependent IMF that consumes $\alpha$ from this chapter, see Environment-dependent IMFs.