Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

`progenax.imf`

San Diego State University

progenax.imf

Auto-generated by scripts/build_api_reference.py from progenax source. Re-run after public-API changes; the script is idempotent.

Module path: progenax/imf/

Public symbols: 46

Contents

imf.IMFProtocol

protocol

IMFProtocol(*args, **kwargs)

Protocol for all IMF implementations.

Any class with these attributes and methods can be used as an IMF, enabling TruncatedIMF to wrap any compatible IMF.

Source: src/progenax/imf/base.py#L25

imf.BaseIMF

class

BaseIMF(m_min: float = 0.0, m_max: float = inf) -> None

Abstract base for all IMFs with shared PPF solver and sampling modes.

Provides automatic normalization, Newton PPF solver with custom gradients, and four sampling modes (N, M_total, M_total_packed, fixed-N).

Subclasses must implement: - _logpdf_unnorm(m): Unnormalized log-PDF (shape function) - _cdf_unnorm(m): Unnormalized CDF (integral of unnormalized PDF)

The base class handles normalization automatically via _log_norm property.

Attributes

ParameterDescription
m_minMinimum mass [M_sun] (default: 0.0)
m_maxMaximum mass [M_sun] (default: inf)

Source: src/progenax/imf/base.py#L92

imf._ppf_newton

function

_ppf_newton(imf: 'BaseIMF', u: Float[Array, '...']) -> Float[Array, '...']

Inverse CDF via fixed Newton iteration.

Thin wrapper over jaxstro’s generic newton_ppf solver: supplies the IMF’s CDF, a log-space initial guess, and the [m_min, m_max] bounds. The Newton loop (fixed iterations, safe division, clip-to-bounds) lives in jaxstro.numerics.rootfinding.newton_ppf — JIT-safe, no convergence loops. Initial guess uses linear interpolation in log-mass space.

Gradients flow through all iterations via automatic differentiation, enabling differentiation w.r.t. both u and IMF parameters.

Args

ParameterDescription
imfIMF instance
uUniform samples in [0, 1]

Returns: Mass values m such that CDF(m) ≈ u

Source: src/progenax/imf/base.py#L48

imf.TruncatedIMF

class

📇 model card

TruncatedIMF(inner: progenax.imf.base.IMFProtocol, m_min: float, m_max: float)

Wrapper for hard mass truncation of any IMF.

Takes an existing IMF and enforces strict bounds [m_min, m_max]. Renormalizes the PDF over the truncated domain.

Attributes

ParameterDescription
innerThe wrapped IMF
m_minMinimum mass [M_sun]
m_maxMaximum mass [M_sun]

Example. >>> from progenax.imf import ChabrierIMF, TruncatedIMF

chabrier = ChabrierIMF() truncated = TruncatedIMF(chabrier, m_min=0.08, m_max=150.0) masses = truncated.sample(key, 1000)

Source: src/progenax/imf/truncated.py#L15

imf.PowerLawIMF

class

📇 model card · ∇ gradient-verified — 9 audit cases

PowerLawIMF(exponents: 'List[float]', breakpoints: 'List[float]', m_min: 'float' = 0.01, m_max: 'float' = 100.0)

N-segment piecewise power-law IMF.

The PDF is: ξ(m) ∝ m^(-α_i) for m ∈ [b_{i-1}, b_i]

where α_i are the exponents and b_i are the breakpoints. Continuity is enforced at breakpoints.

Attributes

ParameterDescription
m_minMinimum mass [M_sun]
m_maxMaximum mass [M_sun]
breakpointsMass breakpoints (excluding m_min, m_max)
exponentsPower-law exponents for each segment
_continuity_factorsMultiplicative factors for PDF continuity
_segment_integralsIntegral of PDF over each segment
_segment_cumprobCumulative probability at segment boundaries

Example. >>> imf = PowerLawIMF.kroupa() # Kroupa (2001)

masses = imf.sample(key, 1000)

Source: src/progenax/imf/power_law.py#L20

imf.prepare_imf_samples

function

prepare_imf_samples(N: 'int', key: 'PRNGKeyArray') -> "Float[Array, 'N']"

Prepare uniform samples for reparameterized IMF sampling.

Args

ParameterDescription
NNumber of samples
keyJAX random key

Returns: Uniform samples in [0, 1]

Source: src/progenax/imf/power_law.py#L287

imf.estimate_N_max_for_M_total

function

estimate_N_max_for_M_total(m_total: 'float', imf: 'PowerLawIMF', safety_factor: 'float' = 2.0) -> 'int'

Estimate N_max for M_total mode sampling.

Source: src/progenax/imf/power_law.py#L302

imf.estimate_pool_size

function

estimate_pool_size(m_total: 'float', imf: 'PowerLawIMF') -> 'int'

Alias for estimate_N_max_for_M_total.

Source: src/progenax/imf/power_law.py#L312

imf.Maschberger

class

📇 model card · ∇ gradient-verified — 3 audit cases

Maschberger(m_min: float = 0.01, m_max: float = 300.0, mu: float = 0.2, alpha: float = 2.3, beta: float = 1.4) -> None

Maschberger (2013) smooth IMF.

Single formula bridging lognormal turnover at low mass and power-law tail at high mass.

PDF: f(m) ∝ (m/μ)^(-α) * (1 + (m/μ)^(1-α))^(-β)

Default parameters from Maschberger (2013) Table 1 (canonical single-star IMF): mu = 0.2 M_sun (scale parameter; the pdf peak is near here) alpha = 2.3 (high-mass slope; Kroupa/Chabrier canonical, not Salpeter’s 2.35) beta = 1.4 (low-mass turnover; gives effective low-mass slope γ=α+β(1-α)=0.48)

Note on m_max: Maschberger (2013) Table 1 adopts the fiducial upper limit m_u = 150 M_sun; these limits “are only needed for the normalization” (Table 1 caption). progenax defaults to m_max = 300 M_sun to admit very massive stars; since the limit only sets the normalization constant, this is a convention choice, not a change to the IMF shape. Pass m_max=150.0 to reproduce the paper exactly.

Reference: Maschberger, T. (2013), MNRAS, 429, 1725 “On the function describing the stellar initial mass function”

Source: src/progenax/imf/smooth.py#L51

imf.TaperedPowerLaw

class

📇 model card

TaperedPowerLaw(m_min: float = 0.01, m_max: float = 300.0, alpha: float = 2.3, m_peak: float = 0.3, beta: float = 2.0) -> None

Tapered Power Law IMF.

PDF: f(m) ∝ m^(-α) * (1 - exp(-(m/m_peak)^β))

The exponential taper suppresses low masses below m_peak, creating a smooth turnover. Functional form: the tapered power law of Parravano, McKee & Hollenbach (2011) (cf. Maschberger 2013, Eq. 12).

No closed-form CDF inverse: the CDF uses the shared cumulative-trapezoid grid and ppf is the inherited BaseIMF fixed-iteration Newton solver (unlike Maschberger, which has an analytic quantile). Differentiable and JIT-safe.

Attributes

ParameterDescription
alphaPower-law slope (default: 2.3, canonical high-mass)
m_peakTurnover/peak mass [M_sun]
betaTaper sharpness (higher = sharper cutoff)

Source: src/progenax/imf/smooth.py#L180

imf.Schechter

class

📇 model card · ∇ gradient-verified — 1 audit case

Schechter(m_min: float = 0.01, m_max: float = 300.0, alpha: float = 2.3, m_star: float = 100.0) -> None

Schechter function with exponential high-mass cutoff.

PDF: f(m) ∝ m^(-α) * exp(-m/m_star)

Originally developed for galaxy luminosity functions (Schechter 1976), also used for IMFs in extreme environments or IGIMF theory.

No closed-form CDF inverse: the CDF uses the shared cumulative-trapezoid grid and ppf is the inherited BaseIMF fixed-iteration Newton solver. Differentiable/JIT-safe.

WARNING: Default parameters (α=2.3, m_star=100) give a function that is essentially just a power-law since the cutoff is far above where the power-law already suppresses high masses.

For stellar IMFs with meaningful high-mass cutoff, consider: - α = 1.0-1.5 (flatter slope) - m_star = 10-50 M_sun (relevant cutoff scale)

Example. >>> # IGIMF-style with cluster-dependent upper mass limit

imf = Schechter(alpha=1.35, m_star=m_max_from_cluster_mass(M_cluster))

Attributes

ParameterDescription
alphaPower-law slope (default: 2.3, Salpeter-like)
m_starExponential cutoff mass [M_sun] (default: 100)

Source: src/progenax/imf/smooth.py#L228

imf.ChabrierIMF

class

📇 model card · ∇ gradient-verified — 4 audit cases

ChabrierIMF(m_min: float = 0.08, m_max: float = 100.0, m_c: float = 0.08, sigma: float = 0.69, alpha: float = 2.3, m_trans: float = 1.0, A_ln: float = 0.158) -> None

Chabrier (2003) lognormal + power-law IMF.

The Chabrier IMF has two components:

Uses standard Chabrier (2003) parameters with A_pl computed for continuity at the transition mass m_trans = 1 M☉.

Implements IMFProtocol for compatibility with BaseIMF framework.

Attributes

ParameterDescription
m_minMinimum mass [M☉] (default: 0.08 - hydrogen burning limit)
m_maxMaximum mass [M☉] (default: 100)
m_cCharacteristic mass for lognormal [M☉] (default: 0.08, Chabrier 2003)
sigmaWidth of lognormal in log-space (default: 0.69)
alphaHigh-mass slope for ξ(m)=dN/dm ∝ m^(-α) (default: 2.3, Chabrier 2003 Table 1 high-mass tail x=1.3 ⇒ α=2.3; the original Salpeter slope is 2.35, available via PowerLawIMF.salpeter())
m_transTransition mass between lognormal and power-law (default: 1.0)
A_lnLognormal coefficient (default: 0.158, Chabrier 2003 single-star disk IMF)
A_plPower-law coefficient (computed for continuity at m_trans)

Examples. >>> imf = ChabrierIMF() # Default Chabrier (2003) single-star (disk) IMF

key = jax.random.PRNGKey(42) masses = imf.sample(key, 1000) print(f"Mean mass: {imf.mean_mass():.3f} M☉") # ~0.61 M☉ over [0.08, 100]

References. Chabrier (2003), PASP, 115, 763 - Table 1: single-star disk IMF coefficients Chabrier (2005), ASSL, 327, 41 - Review of IMF determinations

Source: src/progenax/imf/chabrier.py#L29

imf.MassRatioProtocol

protocol

MassRatioProtocol(*args, **kwargs)

Protocol for mass-ratio distributions.

All mass-ratio distributions must implement: - pdf(q): Probability density at mass ratio q - cdf(q): Cumulative distribution up to q - ppf(u): Inverse CDF (percent point function) - sample(key, n): Generate n samples

Note: q = M_secondary / M_primary ∈ [q_min, 1] by convention.

Source: src/progenax/imf/binary/mass_ratio.py#L22

imf.FlatMassRatio

class

FlatMassRatio(q_min: 'float' = 0.1) -> None

Uniform mass-ratio distribution q ∈ [q_min, 1].

Reference: Raghavan et al. (2010) ApJS 190, 1 - Fig. 16 Wide solar-type binaries show approximately flat q distribution.

Parameters: q_min: Minimum mass ratio (default: 0.1) Below ~0.1, companions become brown dwarfs.

Source: src/progenax/imf/binary/mass_ratio.py#L59

imf.PowerLawMassRatio

class

PowerLawMassRatio(gamma: 'float' = 0.0, q_min: 'float' = 0.1) -> None

Power-law mass-ratio distribution p(q) ∝ q^γ.

Reference: Sana et al. (2012) Science 337, 444 - main text & Fig. 1 (κ ≈ -0.1; Science Report, no numbered eqs) O-star binaries: γ ≈ -0.1 (nearly flat, slight preference for unequal)

Moe & Di Stefano (2017) ApJS 230, 15 - Section 9.1
Solar-type stars: γ ≈ 0 to +0.3 depending on period

Parameters: gamma: Power-law exponent (default: 0.0 = flat) γ < 0: Prefers unequal masses γ > 0: Prefers equal masses γ = 0: Flat distribution q_min: Minimum mass ratio (default: 0.1)

Note. For γ = -1, the distribution is singular at q = 0. Use q_min > 0 to avoid singularity.

Source: src/progenax/imf/binary/mass_ratio.py#L116

imf.TwinPeakedMassRatio

class

TwinPeakedMassRatio(gamma: 'float' = 0.0, f_twin: 'float' = 0.1, sigma_twin: 'float' = 0.03, q_min: 'float' = 0.1) -> None

Mass-ratio distribution with excess of twins (q ≈ 1).

Reference: Moe & Di Stefano (2017) ApJS 230, 15 - Section 9.2 “For solar-type primaries at intermediate periods (P = 10-200 days), we measure an excess fraction of twins F_twin ≈ 0.1 above the baseline power-law distribution.”

Lucy (2006) A&A 457, 629 - systematic statistical study of the twin
excess at q≈1 (the strong twin hypothesis itself traces to Lucy & Ricco 1979)

Model: p(q) = (1 - f_twin) × q^γ / Z_pl + f_twin × N(q | μ=1, σ_twin)

where Z_pl is the power-law normalization.

Parameters: gamma: Power-law exponent for baseline (default: 0.0) f_twin: Fraction of systems that are twins (default: 0.1) sigma_twin: Width of twin peak (default: 0.03) Moe+17 suggest σ ≈ 0.02-0.05 q_min: Minimum mass ratio (default: 0.1)

Source: src/progenax/imf/binary/mass_ratio.py#L275

imf.MoeDiStefano2017

class

MoeDiStefano2017(q_min: 'float' = 0.1, sigma_twin: 'float' = 0.03) -> None

Mass-dependent mass-ratio distribution from Moe & Di Stefano (2017).

Reference: Moe & Di Stefano (2017) ApJS 230, 15, Table 13.

APPROXIMATION — single-slope, period-averaged. Moe & Di Stefano’s actual mass-ratio distribution is a THREE-parameter, PERIOD-DEPENDENT form (Table 13): two power-law slopes γ_smallq (0.1<q<0.3) and γ_largeq (0.3<q<1.0) plus a twin excess F_twin, with all three tabulated as functions of BOTH primary mass AND orbital period (at logP=1,3,5,7). This class collapses that to a SINGLE period-averaged slope γ(M1) and a period-averaged F_twin(M1) — it captures the qualitative trend (low-mass companions favour equal q, OB companions favour small q) but is not a verbatim Table row. A faithful two-slope, period-dependent implementation is tracked for a future release.

Period-averaged single-slope reduction (γ from the qualitative γ_smallq/γ_largeq trend; f_twin period-averaged from Table 13’s F_twin, which falls from ~0.1–0.3 at logP=1 to <0.03 at long P): - M1 < 0.8 Msun: γ ≈ 0.4, f_twin ≈ 0.05 (M-dwarfs) - 0.8 < M1 < 1.2 Msun: γ ≈ 0.3, f_twin ≈ 0.10 (Solar-type; Table 13 F_twin period-avg) - 1.2 < M1 < 3.5 Msun: γ ≈ 0.0, f_twin ≈ 0.08 (A/F stars) - M1 > 3.5 Msun: γ ≈ -0.5, f_twin ≈ 0.03 (OB stars; Table 13 γ_largeq(logP=1)=-0.5)

Parameters: q_min: Minimum mass ratio (default: 0.1) sigma_twin: Width of twin peak (default: 0.03)

Source: src/progenax/imf/binary/moe_di_stefano.py#L20

imf.MoeDiStefano2017Full

class

📇 model card

MoeDiStefano2017Full(q_min: 'float' = 0.1, q_break: 'float' = 0.3, n_grid: 'int' = 512) -> None

Faithful two-slope, period-dependent Moe & Di Stefano (2017) mass-ratio model.

The mass-ratio pdf is a two-slope power law plus a twin excess block, jointly normalized over [q_min, 1.0]:

p_q(q | M1, P) ∝ p_2slope(q) + [twin block on [0.95, 1.0]]

where p_2slope ∝ q^γsmallq on [q_min, 0.3] then q^γlargeq on [0.3, 1.0] (continuous at q=0.3). F_twin is the EXCESS-twin fraction of the q > 0.3 population (MD17 p.5, Fig. 2 — NOT the fraction of all q > q_min companions, and NOT the total q > 0.95 fraction). To realize that convention the twin block carries unnormalized mass ft/(1-ft)·I_B (I_B = power-law mass on [0.3, 1]), so twin/(twin + I_B) = ft exactly (see _components). The pre-fix code mixed ft against the whole [q_min, 1] population, over-weighting twins by ~22% at solar logP=1 (realized 0.367 vs Table-13 0.300) — audit R3. γsmallq, γlargeq, F_twin are bilinearly interpolated (clamped) over Table 13 (verified against the PDF p.52). η(P,M1) for eccentricity is handled by progenax.binaries.MoeEccentricity (which reproduces Table 13’s η).

This is period-dependent (sample takes periods AND masses) — it does NOT fit the unconditional MassRatioProtocol; use it via MoeJointOrbit. The single-slope period-averaged MoeDiStefano2017 remains the fast approximation / BinaryIMF default.

Reference: Moe & Di Stefano (2017) ApJS 230, 15, §9.1 Eqs. 2-3, Table 13 (p.52), Fig. 2.

Parameters: q_min: Minimum mass ratio for the normalized range (default: 0.1).

Source: src/progenax/imf/binary/moe_di_stefano.py#L240

imf.MoePeriod

class

MoePeriod(logP_min: 'float' = 0.2, logP_max: 'float' = 8.0, n_grid: 'int' = 257) -> None

Mass-dependent orbital-period distribution from Moe & Di Stefano (2017) Table 13.

The companion-frequency anchors f_logP;q>0.1(M1) at logP={1,3,5,7} (bilinear in M1) define the period-distribution shape; we sample logP by inverting the normalized cumulative of the piecewise-linear density over [logP_min, logP_max] (clamped at the edges). Differentiable wrt M1 via jnp.interp / cumsum.

Solar-type companion frequency peaks at long periods (logP~5); early-type is flatter and weighted to shorter periods — reproducing Moe’s period trend.

Reference: Moe & Di Stefano (2017) ApJS 230, 15, Table 13, Figure 37.

Source: src/progenax/imf/binary/moe_di_stefano.py#L360

imf.MoeJointOrbit

class

MoeJointOrbit(period: 'MoePeriod', massratio: 'MoeDiStefano2017Full', eccentricity: 'eqx.Module') -> None

Faithful joint (P, q, e) sampler — the Moe & Di Stefano (2017) interrelation.

Given a primary mass M1, samples the correlated orbital parameters: logP ~ MoePeriod(M1); q ~ MoeDiStefano2017Full(M1, P); e ~ MoeEccentricity(P, M1). The P–q–e coupling (short-P -> larger q / twins / circular; long-P -> small q approaching random IMF pairings) is the paper’s central result (“Mind your Ps and Qs”).

Construct via MoeJointOrbit.default() for the standard components, or pass your own. The eccentricity sampler is duck-typed (any .sample(key, periods, masses)) so imf need not hard-import binaries.

Reference: Moe & Di Stefano (2017) ApJS 230, 15 (full joint distribution).

Source: src/progenax/imf/binary/moe_di_stefano.py#L400

imf.ConstantBinaryFraction

class

ConstantBinaryFraction(f_bin: 'float' = 0.5) -> None

Constant binary fraction independent of mass.

Reference: Raghavan et al. (2010) ApJS 190, 1 Overall multiplicity fraction ~46% for solar-type stars.

Parameters: f_bin: Binary fraction (default: 0.5)

Source: src/progenax/imf/binary/binary_fraction.py#L14

imf.MassDependentBinaryFraction

class

MassDependentBinaryFraction() -> None

Mass-dependent MULTIPLICITY fraction (probability a primary has ≥1 companion).

These are the multiplicity fractions f = 1 − (single-star fraction). For M ≥ 0.8 Msun they are derived from Moe & Di Stefano (2017) ApJS 230, 15, Table 13 (the single-star fraction F_{n=0;q>0.1} row): e.g. solar 1−0.60=0.40, B-star 1−0.41≈0.59, O-star 1−0.06≈0.94. Below 0.8 Msun the values come from M-dwarf surveys (Raghavan et al. 2010; Duchêne & Kraus 2013). They are NOT the “close binary fraction” (Table 13’s f_{logP<3.7}) nor the total multiplicity frequency f_mult (which exceeds 1 for massive stars because of triples/quadruples); they are the fraction of primaries with a companion.

Model (multiplicity fraction vs primary mass): - M < 0.1 Msun: f ≈ 0.22 (VLM/brown dwarfs; M-dwarf surveys) - 0.1 < M < 0.5 Msun: f ≈ 0.26 (M-dwarfs; surveys) - 0.5 < M < 1.0 Msun: f ≈ 0.44 (K/G-dwarfs; Raghavan 2010 ≈0.44, Moe solar 1−0.60=0.40) - 1.0 < M < 2.0 Msun: f ≈ 0.50 (F/A-stars) - 2.0 < M < 5.0 Msun: f ≈ 0.60 (B-stars; Moe A/late-B 1−0.41≈0.59) - 5.0 < M < 10 Msun: f ≈ 0.80 (early B; Moe mid/early-B 1−{0.24,0.16}≈0.76–0.84) - M > 10 Msun: f ≈ 0.90 (O-stars; Moe O 1−0.06≈0.94)

Note: pairing one companion per “binary” system; higher-order multiples (triples, etc.) are common at high masses but are not modelled here (use the full Moe model for that).

Source: src/progenax/imf/binary/binary_fraction.py#L36

imf.BinaryIMF

class

BinaryIMF(primary_imf: 'BaseIMF', q_distribution: 'Union[MassRatioProtocol, MoeDiStefano2017, MassRatioSamplerCallable, None]' = None, binary_fraction: 'Union[ConstantBinaryFraction, MassDependentBinaryFraction, BinaryFractionCallable, float, None]' = None) -> None

Initial Mass Function for binary star populations.

Combines a primary-mass IMF with a mass-ratio distribution and binary fraction model to generate complete binary populations.

Default Configuration (Moe & Di Stefano 2017): - Mass-ratio distribution: MoeDiStefano2017 (mass-dependent γ and twin excess) - Binary fraction: MassDependentBinaryFraction (increases with primary mass)

This default follows the comprehensive observational constraints from Moe & Di Stefano (2017) ApJS 230, 15, which represents the most complete characterization of binary statistics to date.

Reference: Kroupa (1995) MNRAS 277, 1491 - IMF-consistent binary populations Moe & Di Stefano (2017) ApJS 230, 15 - Modern comprehensive model Raghavan et al. (2010) ApJS 190, 1 - Solar-type binary statistics Sana et al. (2012) Science 337, 444 - O-star binary statistics

Parameters: primary_imf: IMF for primary stars (must implement BaseIMF interface) q_distribution: Mass-ratio distribution. Options: - MoeDiStefano2017() [DEFAULT] - Mass-dependent from Moe+17 - FlatMassRatio() - Uniform q ∈ [q_min, 1] - PowerLawMassRatio(gamma) - p(q) ∝ q^γ - TwinPeakedMassRatio() - Flat + Gaussian twin peak - Custom callable: f(key, m1) -> q array binary_fraction: Binary fraction model. Options: - MassDependentBinaryFraction() [DEFAULT] - From Moe+17 Table 13 - ConstantBinaryFraction(f_bin) - Constant fraction - float - Shorthand for ConstantBinaryFraction - Custom callable: f(m) -> f_bin array

Examples. >>> from progenax.imf import PowerLawIMF

from progenax.imf.binary import BinaryIMF

Default: full Moe+17 mass-dependent model

imf = BinaryIMF(primary_imf=PowerLawIMF.kroupa())

Equivalent to:

imf = BinaryIMF(

primary_imf=PowerLawIMF.kroupa(),

q_distribution=MoeDiStefano2017(),

binary_fraction=MassDependentBinaryFraction(),

)

Simple constant binary fraction

imf = BinaryIMF( ... primary_imf=PowerLawIMF.kroupa(), ... binary_fraction=0.5, # 50% binaries ... )

Custom mass-ratio sampler

def my_q_sampler(key, m1): ... # Custom logic: q depends on m1 ... return jax.random.uniform(key, m1.shape, minval=0.3, maxval=1.0) imf = BinaryIMF( ... primary_imf=PowerLawIMF.kroupa(), ... q_distribution=my_q_sampler, ... )

Custom binary fraction function

def my_f_bin(m): ... # Custom: 40% for low mass, 80% for high mass ... return jnp.where(m < 1.0, 0.4, 0.8) imf = BinaryIMF( ... primary_imf=PowerLawIMF.kroupa(), ... binary_fraction=my_f_bin, ... )

Source: src/progenax/imf/binary/imf.py#L37

imf.IMFParams

class

∇ gradient-verified — 2 audit cases

IMFParams(alpha0: Float[Array, ''], alpha1: Float[Array, ''], alpha2: Float[Array, ''], alpha3: Float[Array, ''], m_break0: float = 0.08, m_break1: float = 0.5, m_break2: float = 1.0, m_min: float = 0.01, m_max: float = 150.0) -> None

4-segment piecewise power-law IMF parameters.

Matches Marks+2012 Eq. 2 convention: ξ(m) ∝ m^(-α₀) for m ∈ [0.01, 0.08] M☉ (sub-stellar, typically FIXED) ξ(m) ∝ m^(-α₁) for m ∈ [0.08, 0.50] M☉ (low-mass) ξ(m) ∝ m^(-α₂) for m ∈ [0.50, 1.00] M☉ (intermediate) ξ(m) ∝ m^(-α₃) for m ∈ [1.00, m_max] M☉ (high-mass, main target)

Primary inference targets: - α₃: Environment-dependent high-mass slope (main science driver) - α₁, α₂: Metallicity-dependent (optional refinement via Marks Eq. 12) - α₀: Fixed at 0.3 (brown dwarf regime, rarely constrained)

Attributes

ParameterDescription
alpha0Slope for [0.01, 0.08] M☉ (canonical 0.3, typically FIXED)
alpha1Slope for [0.08, 0.50] M☉ (canonical 1.3)
alpha2Slope for [0.50, 1.00] M☉ (canonical 2.3)
alpha3Slope for [1.00, m_max] M☉ (canonical 2.3, environment-dependent)
m_break0First break at 0.08 M☉ (static)
m_break1Second break at 0.50 M☉ (static)
m_break2Third break at 1.00 M☉ (static)
m_minMinimum mass 0.01 M☉ (static)
m_maxMaximum mass 150.0 M☉ (static)

Example. >>> params = IMFParams.kroupa()

print(f"High-mass slope: {params.alpha3}") High-mass slope: 2.3

Top-heavy IMF (e.g., in dense starburst)

params = IMFParams( ... alpha0=jnp.array(0.3), ... alpha1=jnp.array(1.3), ... alpha2=jnp.array(2.3), ... alpha3=jnp.array(1.8), # Top-heavy ... )

Source: src/progenax/imf/params.py#L23

imf.log_prob_masses

function

log_prob_masses(masses: Float[Array, 'N'], params: progenax.imf.params.IMFParams) -> Float[Array, 'N']

Compute log probability of each mass under the 4-segment IMF.

Evaluates the normalized piecewise power-law PDF: log p(m | params) = log(ξ(m)) - log(normalization)

where ξ(m) is the unnormalized IMF with 4 segments.

Args

ParameterDescription
massesStellar masses [M☉], shape (N,)
paramsIMF parameters

Returns: Log probability for each mass, shape (N,)

Example. >>> params = IMFParams.kroupa()

masses = jnp.array([0.5, 1.0, 10.0]) log_probs = log_prob_masses(masses, params)

Source: src/progenax/imf/differentiable.py#L74

imf.sample_masses_from_params

function

sample_masses_from_params(params: progenax.imf.params.IMFParams, u: Float[Array, 'N']) -> Float[Array, 'N']

Sample masses via inverse CDF - fully differentiable.

Uses the reparameterization trick: masses = F⁻¹(u; params) where F is the cumulative distribution function.

Gradients flow through params, not through the random samples u.

Args

ParameterDescription
paramsIMF parameters
uUniform samples in [0, 1], shape (N,)

Returns: Sampled masses [M☉], shape (N,)

Example. >>> params = IMFParams.kroupa()

key = jax.random.PRNGKey(42) u = jax.random.uniform(key, (1000,)) masses = sample_masses_from_params(params, u)

Source: src/progenax/imf/differentiable.py#L176

imf.individual_mass_nll

function

individual_mass_nll(masses: Float[Array, 'N'], params: progenax.imf.params.IMFParams) -> Float[Array, '']

Compute negative log-likelihood for individual resolved masses.

NLL = -Σᵢ log p(mᵢ | params)

This is the simplest likelihood for gradient-based IMF inference.

Args

ParameterDescription
massesObserved stellar masses [M☉], shape (N,)
paramsIMF parameters to evaluate

Returns: Negative log-likelihood (scalar)

Example. >>> params = IMFParams.kroupa()

masses = jnp.array([0.5, 1.0, 10.0, 50.0]) nll = individual_mass_nll(masses, params)

Minimize NLL to fit params to data

Source: src/progenax/imf/differentiable.py#L265

imf.BirthEnvironment

class

📇 model card

BirthEnvironment(metallicity: "Float[Array, ''] | None" = None, log_mecl: "Float[Array, ''] | None" = None, sfe: "Float[Array, ''] | None" = None, log_rho_cl: "Float[Array, ''] | None" = None)

Physical birth environment with SFE parameter.

Represents conditions during star formation that influence the IMF. All fields are JAX arrays for gradient-based inference.

Attributes

ParameterDescription
metallicity[Fe/H], default 0.0 (solar). Calibrated range [-2.5, +0.5]
log_mecllog₁₀(M_ecl / M☉). Calibrated range [3, 8]
sfeStar formation efficiency ε = M_ecl / M_cl, default 0.33
log_rho_clOptional override for log₁₀(ρ_cl / 10⁶ M☉ pc⁻³)

Example. >>> # Solar neighborhood cluster

env = BirthEnvironment.solar() params = env_to_imf_params(env) print(f"α₃ = {float(params.alpha3):.2f}") # ~2.3

Dense, metal-poor globular cluster with low SFE

env = BirthEnvironment.from_cluster_mass( ... M_ecl=1e6, FeH=-1.5, sfe=0.1 ... ) params = env_to_imf_params(env) print(f"α₃ = {float(params.alpha3):.2f}") # Top-heavy!

Source: src/progenax/imf/environment/birth_environment.py#L21

imf.env_to_imf_params

function

📇 model card

env_to_imf_params(env: 'BirthEnvironment', model: 'str' = 'jerabkova_generalized', include_lowmass_variation: 'bool' = False, smooth_alpha3: 'bool' = False, smooth_width: 'float' = 0.2, clamp_domain: 'bool' = True) -> 'IMFParams'

Convert birth environment to 4-segment IMF parameters.

SINGLE API for all environment-dependent IMF models. Fully differentiable - gradients flow from IMFParams back to BirthEnvironment.

Models: - “kroupa”: Standard Kroupa (2001), ignores environment - “jerabkova_generalized”: Eq. 9 with explicit ε (RECOMMENDED) - “jerabkova_mecl”: Eq. 9 (mass-based, assumes ε = 0.33) - “jerabkova_rho”: Eq. 7 (density-based, requires log_rho_cl) - “marks_plane”: Fundamental Plane Eq. 14-15 (requires log_rho_cl) - “marks_mcl”: Table 3 (cloud mass) - “marks_mecl”: Table 3 (stellar mass) - “marks_rho”: Table 3 (density, requires log_rho_cl) - “marks_feh”: Table 3 (metallicity only)

Args

ParameterDescription
envBirthEnvironment with metallicity, log_mecl, sfe
modelModel name (see above)
include_lowmass_variationApply Marks Eq. 12 to α₁, α₂ (default False)
smooth_alpha3Use tanh smoothing at threshold (for gradients)
smooth_widthWidth of tanh transition
clamp_domainClamp inputs to calibrated ranges

Returns: 4-segment IMFParams (α₀, α₁, α₂, α₃)

Example. >>> env = BirthEnvironment.from_cluster_mass(M_ecl=1e6, FeH=-1.5)

params = env_to_imf_params(env) print(f"α₃ = {float(params.alpha3):.2f}")

With explicit SFE

env = BirthEnvironment.from_cluster_mass(M_ecl=1e6, FeH=-1.5, sfe=0.1) params = env_to_imf_params(env, model=“jerabkova_generalized”)

Source: src/progenax/imf/environment/mapping.py#L382

imf.compute_r_half

function

compute_r_half(M_ecl: "Float[Array, '...']") -> "Float[Array, '...']"

Radius-mass relation from Marks & Kroupa (2012).

r_h [pc] = 0.1 × (M_ecl / M☉)^0.13

Args

ParameterDescription
M_eclStellar mass of embedded cluster [M☉]

Returns: Half-mass radius [pc]

Source: src/progenax/imf/environment/density.py#L13

imf.compute_rho_ecl

function

compute_rho_ecl(M_ecl: "Float[Array, '...']") -> "Float[Array, '...']"

Half-mass density (Marks & Kroupa 2012 definition).

ρ_ecl = 3 × M_ecl / (8π × r_h³)

The 8π factor arises because half the cluster mass is within r_h: ρ = (M_ecl/2) / (4π/3 × r_h³) = 3M_ecl / (8π × r_h³)

This is the authoritative definition in Marks & Kroupa (2012), A&A 543, A8 (p. 2, “ρ_ecl = 3 M_ecl/8π r_h³”), and reproduces Marks+2012 (MNRAS 422, 2246) Table 1 densities exactly (e.g. NGC 104: 3·9.40e6/(8π·0.49³)=9.54e6 = the tabulated ρ_cl). NOTE: Jerabkova+2018 Eq. 8 writes 4π, but that is internally inconsistent with its own ρ_ecl=0.61·logM+2.08 relation (which is 8π); progenax follows the 8π convention that matches the actual α₃–ρ calibration data.

Args

ParameterDescription
M_eclStellar mass of embedded cluster [M☉]

Returns: Stellar half-mass density [M☉ pc⁻³]

Source: src/progenax/imf/environment/density.py#L27

imf.compute_rho_cl

function

compute_rho_cl(M_ecl: "Float[Array, '...']", sfe: "Float[Array, '...']") -> "Float[Array, '...']"

Cloud-core density from stellar mass and SFE.

ρ_cl = ρ_ecl / ε

Args

ParameterDescription
M_eclStellar mass of embedded cluster [M☉]
sfeStar formation efficiency ε = M_ecl / M_cl

Returns: Cloud density [M☉ pc⁻³]

Source: src/progenax/imf/environment/density.py#L52

imf.compute_log_rho_cl_6

function

compute_log_rho_cl_6(M_ecl: "Float[Array, '...']", sfe: "Float[Array, '...']") -> "Float[Array, '...']"

log₁₀(ρ_cl / 10⁶ M☉ pc⁻³) from cluster mass and SFE.

Args

ParameterDescription
M_eclStellar mass of embedded cluster [M☉]
sfeStar formation efficiency

Returns: log₁₀(ρ_cl / 10⁶)

Source: src/progenax/imf/environment/density.py#L70

imf.x_jerabkova_generalized

function

x_jerabkova_generalized(FeH: "Float[Array, '...']", M_ecl: "Float[Array, '...']", sfe: "Float[Array, '...']") -> "Float[Array, '...']"

Mass-based x with explicit ε dependence.

x = -0.14 × [Fe/H] + 0.6039 × log₁₀(M_ecl/10⁶) + 0.2161 - 0.99 × log₁₀(ε/0.33)

Derived from Jeřábková Eq. 7 using Marks & Kroupa r_h-M_ecl relation and 8π half-mass density convention. See JERABKOVA_COEFFICIENTS for derivation.

Args

ParameterDescription
FeHMetallicity [Fe/H]
M_eclCluster stellar mass [M☉]
sfeStar formation efficiency ε

Returns: x parameter for α₃ calculation

Source: src/progenax/imf/environment/mapping.py#L29

imf.x_jerabkova_rho

function

x_jerabkova_rho(FeH: "Float[Array, '...']", log_rho_6: "Float[Array, '...']") -> "Float[Array, '...']"

Jerabkova+2018 Eq. 7: x from density.

x = -0.14 × [Fe/H] + 0.99 × log₁₀(ρ_cl/10⁶)

NOTE: This density-based formula has NO constant term (unlike the mass-based Eq. 9).

Args

ParameterDescription
FeHMetallicity [Fe/H]
log_rho_6log₁₀(ρ_cl / 10⁶ M☉ pc⁻³)

Returns: x parameter for α₃ calculation

Source: src/progenax/imf/environment/mapping.py#L60

imf.x_hat_marks_plane

function

x_hat_marks_plane(FeH: "Float[Array, '...']", log_rho_6: "Float[Array, '...']") -> "Float[Array, '...']"

Marks+2012 Fundamental Plane coordinate (Eq. 14).

x̂ = cos(θ) × [Fe/H] + sin(θ) × log₁₀(ρ_cl/10⁶) = -0.139 × [Fe/H] + 0.990 × log₁₀(ρ_cl/10⁶)

Args

ParameterDescription
FeHMetallicity [Fe/H]
log_rho_6log₁₀(ρ_cl / 10⁶ M☉ pc⁻³)

Returns: x̂ coordinate on Fundamental Plane

Source: src/progenax/imf/environment/mapping.py#L82

imf.alpha3_jerabkova_generalized

function

alpha3_jerabkova_generalized(FeH: "Float[Array, '...']", M_ecl: "Float[Array, '...']", sfe: "Float[Array, '...']", smooth: 'bool' = False, smooth_width: 'float' = 0.2) -> "Float[Array, '...']"

α₃ from generalized Jerabkova with explicit ε (RECOMMENDED).

Uses x_jerabkova_generalized(), the mass-based x built self-consistently from Jerabkova Eq. 7 + the Marks r_h–M_ecl relation + the 8π half-mass density (constant 0.2161; see coefficients.py). NOTE: this is NOT the constant (2.83) printed in Jerabkova Eq. 9 — that printed value is internally inconsistent with her own Eq. 7+8 density relation (which reconstructs to ~0.50, not 2.83). We deliberately use the density-consistent 8π constant so the mass- and density- based paths agree; it does not reproduce Eq. 9 as literally printed.

Args

ParameterDescription
FeHMetallicity [Fe/H]
M_eclCluster stellar mass [M☉]
sfeStar formation efficiency
smoothUse tanh smoothing
smooth_widthSmoothing width

Returns: α₃, clipped to [0.5, 2.3]

Source: src/progenax/imf/environment/mapping.py#L148

imf.alpha3_jerabkova_mecl

function

alpha3_jerabkova_mecl(log_mecl_6: "Float[Array, '...']", FeH: "Float[Array, '...']", smooth: 'bool' = False, smooth_width: 'float' = 0.2) -> "Float[Array, '...']"

Mass-based α₃ from cluster mass (assumes ε = 0.33).

x = -0.14 × [Fe/H] + 0.6039 × log₁₀(M_ecl/10⁶) + 0.2161

Args

ParameterDescription
log_mecl_6log₁₀(M_ecl / 10⁶ M☉)
FeHMetallicity [Fe/H]
smoothUse tanh smoothing
smooth_widthSmoothing width

Returns: α₃, clipped to [0.5, 2.3]

Source: src/progenax/imf/environment/mapping.py#L188

imf.alpha3_jerabkova_rho

function

alpha3_jerabkova_rho(log_rho_6: "Float[Array, '...']", FeH: "Float[Array, '...']", smooth: 'bool' = False, smooth_width: 'float' = 0.2) -> "Float[Array, '...']"

Jerabkova+2018 Eq. 7: α₃ from density.

Args

ParameterDescription
log_rho_6log₁₀(ρ_cl / 10⁶ M☉ pc⁻³)
FeHMetallicity [Fe/H]
smoothUse tanh smoothing
smooth_widthSmoothing width

Returns: α₃, clipped to [0.5, 2.3]

Source: src/progenax/imf/environment/mapping.py#L220

imf.alpha3_marks_plane

function

alpha3_marks_plane(log_rho_6: "Float[Array, '...']", FeH: "Float[Array, '...']", smooth: 'bool' = False, smooth_width: 'float' = 0.2) -> "Float[Array, '...']"

Marks+2012 Fundamental Plane (Eq. 14-15), with the 2014 erratum applied.

x̂ = -0.139 × [Fe/H] + 0.990 × log₁₀(ρ_cl/10⁶) α₃ = -0.4072·x̂ + 1.9383 for x̂ ≥ -0.87 (else canonical 2.3)

The threshold is -0.87 (NEGATIVE): the originally printed Marks+2012 Eq.14/15 had a missing-minus-sign typo (“x̂ ≥ 0.87”), corrected by the 2014 erratum (Marks et al. 2014, MNRAS 442, 3315) and visible in Marks+2012 Fig.6, where the canonical-plateau knee sits at x̂ ≈ -0.87 (the line meets 2.3 continuously there). With this correction the Marks plane coincides with the Jerabkova (2018) IGIMF density relation alpha3_jerabkova_rho (which adopts the same erratum-corrected form) to within the -0.4072-vs-(-0.41) rounding (~0.01).

Args

ParameterDescription
log_rho_6log₁₀(ρ_cl / 10⁶ M☉ pc⁻³)
FeHMetallicity [Fe/H]
smoothUse tanh smoothing
smooth_widthSmoothing width

Returns: α₃, clipped to [0.5, 2.3]

Source: src/progenax/imf/environment/mapping.py#L250

imf.alpha3_marks_table3

function

alpha3_marks_table3(lambda_param: "Float[Array, '...']", relation: 'str', smooth: 'bool' = False, smooth_width: 'float' = 0.2) -> "Float[Array, '...']"

Marks+2012 Table 3: 1D relations for α₃.

α₃(λ) = p × λ + q, if λ ≷ λ_lim

Args

ParameterDescription
lambda_paramParameter value (log₁₀(M/10⁶), log₁₀(ρ/10⁶), or [Fe/H])
relationOne of “mcl”, “mecl”, “rho”, “feh”
smoothUse tanh smoothing
smooth_widthSmoothing width

Returns: α₃, clipped to [0.5, 2.3]

Source: src/progenax/imf/environment/mapping.py#L291

imf.lowmass_slopes_metallicity

function

lowmass_slopes_metallicity(FeH: "Float[Array, '...']", clamp_FeH: 'bool' = True) -> "tuple[Float[Array, '...'], Float[Array, '...']]"

Marks+2012 Eq. 12: Low-mass slopes from metallicity.

For 4-segment IMF convention: α₁([Fe/H]) = 1.3 + 0.5 × [Fe/H] (0.08-0.50 M☉) α₂([Fe/H]) = 2.3 + 0.5 × [Fe/H] (0.50-1.00 M☉)

TENTATIVE: Extrapolation for [Fe/H] < -0.5 is uncertain.

Args

ParameterDescription
FeHMetallicity [Fe/H]
clamp_FeHIf True, clamp to calibrated range [-2.5, +0.5]

Returns: (α₁, α₂) - slopes for segments 1 and 2

Source: src/progenax/imf/environment/mapping.py#L345

imf.JERABKOVA_COEFFICIENTS

value

dict() -> new empty dictionary dict(mapping) -> new dictionary initialized from a mapping object’s (key, value) pairs dict(iterable) -> new dictionary initialized as if via: d = {} for k, v in iterable: d[k] = v dict(**kwargs) -> new dictionary initialized with the name=value pairs in the keyword argument list. For example: dict(one=1, two=2)

imf.MARKS_COEFFICIENTS

value

dict() -> new empty dictionary dict(mapping) -> new dictionary initialized from a mapping object’s (key, value) pairs dict(iterable) -> new dictionary initialized as if via: d = {} for k, v in iterable: d[k] = v dict(**kwargs) -> new dictionary initialized with the name=value pairs in the keyword argument list. For example: dict(one=1, two=2)

imf.MARKS_TABLE3_COEFFICIENTS

value

dict() -> new empty dictionary dict(mapping) -> new dictionary initialized from a mapping object’s (key, value) pairs dict(iterable) -> new dictionary initialized as if via: d = {} for k, v in iterable: d[k] = v dict(**kwargs) -> new dictionary initialized with the name=value pairs in the keyword argument list. For example: dict(one=1, two=2)

imf.DEFAULT_SFE

value

Convert a string or number to a floating-point number, if possible.