zeus21.bursty_sfh

Compute Star Formation Histories with Burstiness.

Authors: zeus21 v2 collaboration - June 2026: Emily Bregou ; Hector Afonso G. Cruz ; Sarah Libanore ; Julian B. Muñoz ; Yonny Sklansky ; Emilie Thélie ; Alessandra Venditti

arXiv:2302.08506, arXiv:2306.09403, arXiv:2407.18294, Sklansky et al. (in prep)

Classes

SFH_class

Star formation histories (SFHs) with stochastic burstiness for Pop II galaxies.

Module Contents

class zeus21.bursty_sfh.SFH_class(UserParams, CosmoParams, AstroParams, HMFinterp, tage, zobs, z_Init=None, SFRD_Init=None)[source]

Star formation histories (SFHs) with stochastic burstiness for Pop II galaxies.

Computes the SFR(t) of a galaxy in a halo of mass Mh, including a power-spectral-density (PSD) model for log-normal SFR fluctuations (damped random walk in ln SFR). See 2601.07912 for the burstiness model.

Only Pop II is currently supported; Pop III raises a ValueError.

Parameters:

UserParams (User_Parameters) – Global run settings.
CosmoParams (Cosmo_Parameters) – Cosmological parameters, including tageofzMyr and zfoftageMyr.
AstroParams (Astro_Parameters) –
Astrophysical parameters. Key attributes used here:
- USE_POPIII : bool — if True, raises ValueError (not implemented).
- FLAG_COMPARE_BAGPIPES : bool — use a toy exponential SFH for comparison with BAGPIPES fits instead of the full model.
- FLAG_RENORMALIZE_AVG_SFH : bool — if True, boosts mean SFR by exp(σ²/2) to account for log-normal stochasticity (<e^x>).
- sigmaPSD, dsigmaPSDdlog10Mh : PSD amplitude and its mass slope.
- tauPSD, dlog10tauPSDdlog10Mh : PSD timescale (Myr) and mass slope.
- _minsigmaPSD, _maxsigmaPSD, _mintauPSD, _maxtauPSD : clamp limits.
- _omegamin, _omegamax : frequency integration range (1/Myr).
- _dt_FFT, _N_FFT : time resolution (Myr) and number of points for FFT.
- epsstar, dlog10epsstardz, _zpivot, Mc, alphastar, betastar, fstarmax, mean_SFR_normalization : star-formation efficiency params.
HMFinterp (HMF_interpolator) – Halo mass function interpolator; must expose Mhtab.
tage (float or array-like) – Lookback time(s) in Myr at which to evaluate the SFH.
zobs (float) – Observed redshift.
z_Init (float, optional) – Initialisation redshift for the SFRD. Computed via Z_init if not given.
SFRD_Init (SFRD_class, optional) – Pre-computed SFRD object. Constructed internally if not given.

SFH_II[source]

Pop II star formation rate in M☉ yr⁻¹ as a function of halo mass and lookback time.

Type:: ndarray, shape (NMh, Ntage)

SFH_II[source]

SFH(CosmoParams, AstroParams, HMFinterp, SFRD_Init, tage, zobs, pop)[source]

Compute the star formation history at a set of lookback times.

Reconstructs halo mass history assuming exponential accretion, then multiplies the mass accretion rate by the stellar efficiency fstar(z, Mh). Optionally renormalises the mean SFR to account for log-normal scatter (see FLAG_RENORMALIZE_AVG_SFH).

Parameters:

CosmoParams (Cosmo_Parameters)
AstroParams (Astro_Parameters)
HMFinterp (HMF_interpolator)
SFRD_Init (SFRD_class)
tage (float or array-like) – Lookback times in Myr relative to zobs.
zobs (float) – Observed redshift.
pop ({2}) – Stellar population. Only Pop II (2) is implemented.

Returns:

SFR in M☉ yr⁻¹ at each (halo mass, lookback time).

Return type:

ndarray, shape (NMh, Ntage)

fstarofz_scaled_Mz(AstroParams, CosmoParams, SFRD_Init, z, Mhlist, pop)[source]

Factorized approximation to fstar(z, Mh) over a 2-D (Mh, z) grid.

Computes fstar(Mh, z) ≈ fstar(Mh, z0) · (Mh / Mh0)^α, separating the mass and redshift dependences. This avoids evaluating fstar_ofz on the full NMh × Nz array; it is exact for Mh < Mc and any residual offset is absorbed into β*.

Only valid under exponential accretion. Use with care if the accretion history deviates strongly from exponential.

Parameters:

AstroParams (Astro_Parameters)
CosmoParams (Cosmo_Parameters)
SFRD_Init (SFRD_class)
z (array-like, shape (Nz,)) – Redshift array (z[0] is used as the reference epoch).
Mhlist (ndarray, shape (NMh, Nz)) – Halo mass grid at each redshift step.
pop ({2}) – Stellar population.

Returns:

Approximate fstar values over the full (Mh, z) grid.

Return type:

ndarray, shape (NMh, Nz)

sigmaPSD_at_Mh(AstroParams, Mh)[source]

PSD amplitude σ(Mh) for the ln SFR damped random walk.

Linear in log10(Mh), clamped to [_minsigmaPSD, _maxsigmaPSD]. To convert to log10(SFR) units, multiply by log(10).

Parameters:

AstroParams (Astro_Parameters)
Mh (float or array-like) – Halo mass(es) in M☉.

Returns:

σ in units of ln SFR, same shape as Mh.

Return type:

ndarray

tauPSD_at_Mh(AstroParams, Mh)[source]

PSD correlation timescale τ(Mh) for the ln SFR damped random walk.

Power-law in Mh, clamped to [_mintauPSD, _maxtauPSD].

Parameters:

AstroParams (Astro_Parameters)
Mh (float or array-like) – Halo mass(es) in M☉.

Returns:

τ in Myr, same shape as Mh.

Return type:

ndarray

PowerlnSFR(AstroParams, omega, Mh)[source]

Power spectrum of ln SFR under a damped (α=2) random walk model.

P(ω, Mh) = σ²(Mh) · τ(Mh) / [1 + (ω τ)²]

As in 2601.07912. Note that compared to other refs (eg 2410.21409) we do not have an implicit factor of 1 Myr in the amplitude; σ here is in ln SFR units; can convert by sigma->sigma*np.log(10) Index for random walk us alpha = 2 by default, can be enhanced to modify.

Parameters:

AstroParams (Astro_Parameters)
omega (float or array-like, shape (Nω,)) – Angular frequency in Myr⁻¹.
Mh (float or array-like, shape (NMh,)) – Halo mass(es) in M☉.

Returns:

Power spectrum in (ln SFR)² Myr.

Return type:

ndarray, shape (NMh, Nω)

Wink_TH(omega, T)[source]

Top-hat temporal window function in Fourier space.

W(ω, T) = sinc(ωT/2) = sin(ωT/2) / (ωT/2)

Used to compute the variance of ln SFR averaged over a timescale T. A small offset (1e-16) is added to the argument to avoid 0/0 at ω=0.

Parameters:

omega (array-like) – Angular frequency in Myr⁻¹.
T (float) – Averaging timescale in Myr.

Returns:

Dimensionless window, same shape as omega. Equals 1 at ω→0.

Return type:

ndarray

Variance_of_lnSFR(AstroParams, T, Mh)[source]

Variance of ln SFR averaged over a timescale T.

Var[ln SFR]_T = (2/2π) ∫ P(ω) |W(ω,T)|² dω

The factor 2/(2π) converts the one-sided integral over positive ω to the full two-sided variance.

Parameters:

AstroParams (Astro_Parameters)
T (float) – Averaging timescale in Myr. Use T=0 for the unsmoothed variance (integrates over all frequencies).
Mh (float or array-like) – Halo mass(es) in M☉.

Returns:

Variance in (ln SFR)², same shape as Mh.

Return type:

ndarray

WindowFourier(CosmoParams, AstroParams, HMFinterp, SFRD_Init, GreensFunction, zobs, tage, pop)[source]

Fourier transform of the convolution G(t) * SFH(t).

Computes W̃(ω, Mh) = FFT[ G(t) · SFH(t) ], where G is a Green’s function (e.g. Greens_function_LUV) and SFH(t) is the star formation history. The result enters the luminosity power spectrum as P_L(ω) = |W̃(ω)|² · P_SFR(ω).

Parameters:

CosmoParams (Cosmo_Parameters)
AstroParams (Astro_Parameters) – Must expose: _dt_FFT, _N_FFT.
HMFinterp (HMF_interpolator)
SFRD_Init (SFRD_class)
GreensFunction (callable) – A function with signature G(AstroParams, tage_Myr, Mhtab) returning an array of shape (NMh, Nt) in erg s⁻¹ M☉⁻¹. Typically one of the Greens_function_L* functions from sed.py.
zobs (float) – Observed redshift.
tage (float) – Maximum lookback time in Myr (sets the FFT window).
pop ({2}) – Stellar population.

Returns:

omegalist (ndarray, shape (Nfft//2 + 1,)) – Angular frequencies in Myr⁻¹.
windowFourier (ndarray, shape (NMh, Nfft//2 + 1)) – Complex Fourier transform of G·SFH in erg s⁻¹ M☉⁻¹ Myr.

Notes

SFH is converted from M☉ yr⁻¹ to M☉ Myr⁻¹ (×10⁶) before the FFT so that the output is in consistent Myr-based units.