Population Synthesis

Population Synthesis#

Starburst99 Reader#

class pyathena.util.SB99(basedir='/projects/EOSTRIKE/SB99/Z014_M1E6_GenevaV00_dt01', verbose=True)[source]#

Bases: object

Interface for reading and plotting Starburst99 (SB99) output.

Reads the text output files produced by the Starburst99 stellar population synthesis code and computes derived quantities (band-integrated luminosities, dust cross sections, photoionization cross sections, momentum injection rates).

Parameters

basedir (str) – Directory containing the SB99 output files (e.g. *.spectrum1, *.power1, *.snr1). All files in the directory must share a common prefix, which is determined automatically.
verbose (bool, optional) – If True, print cluster mass and star-formation parameters on initialization. Default is True.

basedir#

Directory passed at construction time.

Type: str

prefix#

Common filename prefix discovered from the output directory.

Type: str

files#

Mapping of file-type key (e.g. 'spectrum', 'power', 'snr') to the corresponding file path.

Type: dict

logM#

log10 of the instantaneous-burst cluster mass in solar masses (set to 0 for continuous star-formation runs).

Type: float

SFR#

Star-formation rate in M_sun yr-1 (continuous SF runs only).

Type: float

tmax_Myr#

Maximum simulation time in Myr read from the .output1 file.

Type: float

cont_SF#

True if the run assumes continuous star formation.

Type: bool

Examples

>>> import pyathena as pa
>>> sb = pa.SB99('/path/to/sb99/Z014_M1E6_GenevaV00_logdt')
>>> rr = sb.read_rad()        # radiation output
>>> rw = sb.read_wind()       # stellar wind output
>>> rs = sb.read_sn()         # supernova output

Parameters

basedir (str) – directory name in which SB99 simulation results are located
verbose (bool) – Print some useful parameters

_find_files()[source]#

_read_params()[source]#: Read some important parameters Full information contained in self.par

static plt_Edot_pdot_evol_cumul(rr, rw, rs, plt_sn=True, lw=2, normed=False, fig=None, axes=None)[source]#

Plot a 2×2 summary figure of luminosity and momentum (evolution + cumulative).

Combines plt_lum_evol(), plt_pdot_evol(), plt_lum_cumul(), and plt_pdot_cumul() into a single figure.

Parameters

rr (dict) – Output of read_rad().
rw (pandas.DataFrame) – Output of read_wind().
rs (pandas.DataFrame) – Output of read_sn().
plt_sn (bool, optional) – If True, include SN contributions. Default is True.
lw (float, optional) – Line width. Default is 2.
normed (bool, optional) – Normalize cumulative panels. Default is False.
fig (matplotlib.figure.Figure, optional) – Existing figure to draw into. Created if None.
axes (array of matplotlib.axes.Axes, optional) – Four axes to draw into. Created if None.

Returns

fig (matplotlib.figure.Figure)
axes (numpy.ndarray of matplotlib.axes.Axes)

static plt_lum_cumul(ax, rr, rw, rs, normed=True, plt_sn=False, lw=2)[source]#

Plot cumulative energy by band vs age.

Parameters

ax (matplotlib.axes.Axes) –
rr (dict) – Output of read_rad().
rw (pandas.DataFrame) – Output of read_wind().
rs (pandas.DataFrame) – Output of read_sn().
normed (bool, optional) – If True, normalize by the total bolometric cumulative energy. Default is True.
plt_sn (bool, optional) – If True, also plot cumulative SN energy. Default is False.
lw (float, optional) – Line width. Default is 2.

Returns

ax

Return type

matplotlib.axes.Axes

static plt_lum_evol(ax, rr, rw, rs, lw=2, plt_sn=False)[source]#

Plot luminosity evolution by band vs age.

Parameters

ax (matplotlib.axes.Axes) –
rr (dict) – Output of read_rad().
rw (pandas.DataFrame) – Output of read_wind().
rs (pandas.DataFrame) – Output of read_sn().
lw (float, optional) – Line width. Default is 2.
plt_sn (bool, optional) – If True, also plot the SN luminosity. Default is False.

Returns

ax

Return type

matplotlib.axes.Axes

static plt_pdot_cumul(ax, rr, rw, rs, normed=False, plt_sn=False, lw=2)[source]#

Plot cumulative momentum injection vs age.

Parameters

ax (matplotlib.axes.Axes) –
rr (dict) – Output of read_rad().
rw (pandas.DataFrame) – Output of read_wind().
rs (pandas.DataFrame) – Output of read_sn().
normed (bool, optional) – If True, normalize by the total cumulative momentum. Default is False.
plt_sn (bool, optional) – Unused (reserved). Default is False.
lw (float, optional) – Line width. Default is 2.

Returns

ax

Return type

matplotlib.axes.Axes

static plt_pdot_evol(ax, rr, rw, rs, lw=2)[source]#

Plot radiation momentum injection rate evolution vs age.

Parameters

ax (matplotlib.axes.Axes) –
rr (dict) – Output of read_rad().
rw (pandas.DataFrame) – Output of read_wind().
rs (pandas.DataFrame) – Output of read_sn().
lw (float, optional) – Line width. Default is 2.

Returns

ax

Return type

matplotlib.axes.Axes

static plt_spec_sigmad(rr, lambda_Llambda=False, plt_isrf=False, tmax=50.0, nstride=10)[source]#

Plot the SB99 spectrum evolution with dust opacity.

Creates a figure with two or three rows: (1) dust opacity vs wavelength, (2) spectral luminosity at different ages (color-coded by age), and optionally (3) the ISRF at the midplane of a plane-parallel disk.

Parameters

rr (dict) – Output of read_rad().
lambda_Llambda (bool, optional) – If True, plot lambda * L_lambda instead of L_lambda. Default is False.
plt_isrf (bool, optional) – If True, add a third panel showing the ISRF estimate from plt_nuJnu_mid_plane_parallel(). Requires a separate continuous-SF SB99 run; disabled by default.
tmax (float, optional) – Maximum age to show in the spectrum panel [Myr]. Default is 50.
nstride (int, optional) – Plot every nstride-th time step. Default is 10.

Returns

fig

Return type

matplotlib.figure.Figure

read_quanta()[source]#

Read the .quanta1 ionizing-photon rate output file.

Returns: df – DataFrame with columns time_yr, time_Myr, ionizing photon rates Q_HI, Q_HeI, Q_HeII (in s-1 M_sun-1), luminosity fractions Lfrac_HI, Lfrac_HeI, Lfrac_HeII, and bolometric luminosity logL (log10 erg s-1 M_sun-1).
Return type: pandas.DataFrame

read_rad()[source]#

Read the .spectrum1 file and compute derived radiation quantities.

Integrates the SB99 spectrum over standard photometric bands (LyC, LW, PE, OPT, IR, FUV) and computes luminosity-weighted mean dust cross sections, photoionization cross sections, and photon energies using the Draine (2003) dust model (R_V = 3.1).

Band definitions (wavelength in Angstrom):

Band	Wavelength range
LyC	< 912
LW	912 – 1108
PE	1108 – 2068
FUV	912 – 2068
OPT	2068 – 10000
IR	10000 – 200000

Returns

r – Dictionary with the following keys:

Spectral data

time_Myr, time_yr : 1-D array (ntime,) — age of the stellar population in Myr and yr.
wav : 1-D array (nwav,) — wavelength grid [Angstrom].
logf : 2-D array (ntime, nwav) — log10 of the spectral luminosity \(L_\lambda\) [erg s-1 Å-1] for the full instantaneous-burst cluster of mass \(10^{\rm logM}\,M_\odot\). To obtain the specific luminosity per unit stellar mass use 10.0 ** (logf - logM) [erg s-1 Å-1 Msun-1].
logM : float — log10 cluster mass [\(M_\odot\)]; equals self.logM (0 for continuous-SF runs).
df : pandas.DataFrame — raw per-time-step spectral table with columns time_Myr, time_yr, wav, logf, logfstar (stellar continuum), logfneb (nebular continuum). Useful for plotting the full SED at individual ages.

Band-integrated quantities (each is a 1-D array over time)

L : dict — specific luminosity per unit stellar mass [\(L_\odot\,M_\odot^{-1}\)]; keys are 'tot', 'LyC', 'LW', 'PE', 'FUV', 'OPT', 'IR', 'UV'.
pdot : dict — specific radiation momentum injection rate [\(M_\odot\,{\rm km\,s}^{-1}\,{\rm Myr}^{-1}\,M_\odot^{-1}\)]; same keys as L.
Q : dict — specific ionizing-photon emission rate [\({\rm s}^{-1}\,M_\odot^{-1}\)]; same keys as L.

Dust and photoionization cross sections (1-D arrays over time)

Cext, Cabs, Crpr : dict — luminosity-weighted mean dust cross sections per H nucleon [\({\rm cm}^2\,H^{-1}\)]; keyed by band.
hnu : dict — luminosity-weighted mean photon energy [eV]; keyed by band.
sigma_pi_H, sigma_pi_H2 : array — luminosity-weighted mean photoionization cross sections for H and H2 [cm2].
dhnu_H_LyC, dhnu_H2_LyC : array — luminosity-weighted mean excess photon energy above the ionization threshold for H and H2 [eV].

Timescales (each is a dict keyed by band)

tdecay_lum : luminosity-weighted mean emission time [Myr].
tcumul_lum_50, tcumul_lum_90 : age at which 50%/90% of the cumulative band luminosity has been emitted [Myr].

Cumulative time-averaged quantities

L_avg, Q_avg : time-averaged L and Q.
hnu_Lavg, Cabs_Lavg, Cext_Lavg, Crpr_Lavg : L-weighted running averages.
hnu_Qavg, Cabs_Qavg, Cext_Qavg, Crpr_Qavg : Q-weighted running averages.
sigma_pi_H_Qavg, sigma_pi_H2_Qavg, dhnu_H_LyC_Qavg, dhnu_H2_LyC_Qavg : Q-weighted running averages of the LyC cross sections and excess energies.

Return type

dict

read_sn()[source]#

Read the .snr1 supernova rate output file.

Returns: df – DataFrame with columns time_Myr, time_yr, SN rate, energy injection rates (Edot_SN, Edot_SN_IB, Edot_tot in erg s-1 M_sun-1), and cumulative energy (Einj_*), all normalized by the cluster mass.
Return type: pandas.DataFrame

read_wind()[source]#

Read the .power1 stellar wind output file.

Returns: df – DataFrame indexed by time with columns for energy injection rates (Edot_all/OB/RSG/LBV/WR in erg s-1 M_sun-1), momentum injection rates (pdot_all/... in M_sun km s-1 Myr-1 M_sun-1), mass-loss rates (Mdot_all/... in M_sun Myr-1 M_sun-1), and terminal wind velocities (Vw_all/... in km s-1). Also includes time-averaged quantities (Edot_*_avg, pdot_*_avg).
Return type: pandas.DataFrame

read_yield()[source]#

Read the .yield1 chemical yield output file.

Returns: df – DataFrame with columns time_yr, time_Myr, and mass-loss rates (H, He, C, N, O, Mg, Si, S, Fe, wind, SN, combined) in M_sun yr-1 M_sun-1, normalized by the cluster mass.
Return type: pandas.DataFrame

pyathena.util.sb99.get_ISRF_SB99_plane_parallel(rr=None, rr_long=None, Sigma_gas=<Quantity 10. solMass / pc2>, Sigma_SFR=<Quantity 0.0025 solMass / (yr kpc2)>, sfh=None, t_max_Myr=10000.0, Z_dust=1.0, dust_kind='Rv31', verbose=True)[source]#

Compute the plane-parallel ISRF by convolving SB99 SSP spectra with a star-formation history.

Estimates the angle-averaged mean intensity J_lambda at the disk midplane following Ostriker et al. (2010). For a constant SFR the effective spectrum is

\[\frac{L_\lambda}{\rm SFR} = \int_0^{t_{\rm max}} \frac{L_\lambda^{\rm SSP}(t)}{M_*}\,dt\]

where \(L_\lambda^{\rm SSP}(t)/M_*\) is the SSP spectral luminosity per unit cluster mass at stellar age \(t\).

For a time-varying SFH supplied via the sfh parameter the luminosity per unit area is computed directly as

\[\Sigma_\lambda = \int_0^{t_{\rm max}} \frac{L_\lambda^{\rm SSP}(t)}{M_*}\,\Sigma_{\rm SFR}(t)\,dt\]

where \(\Sigma_{\rm SFR}(t)\) is the SFR surface density at lookback time \(t\) returned by the callable sfh.

By default the function stitches together the two bundled SSP datasets:

Z014_M1E6_GenevaV00_dt02 — 251 uniform steps at dt = 0.2 Myr, 0–50 Myr (fine short-timescale coverage for integration accuracy).
Z014_M1E6_GenevaV00_logdt_10Gyr — 77 log-spaced steps, 2 Myr–10 Gyr (long-timescale tail, slow evolution well captured by log spacing).

The stitch point is at the end of rr (i.e. rr['time_Myr'][-1]); rr_long contributes only time steps beyond that.

Note

For higher accuracy at short timescales, pass rr from the original fine-dt run (dt = 0.1 Myr, 1000 steps) if available locally.

Parameters

rr (dict, optional) – Output of SB99.read_rad() for a short-timescale SSP run (e.g. 0.1–100 Myr). If None, the bundled Z014_M1E6_GenevaV00_logdt dataset is loaded automatically.
rr_long (dict, optional) – Output of SB99.read_rad() for a long-timescale SSP run (e.g. 2 Myr–10 Gyr). If None, the bundled Z014_M1E6_GenevaV00_logdt_10Gyr dataset is loaded automatically.
Sigma_gas (astropy.units.Quantity, optional) – Gas surface density used to compute the dust optical depth. Default is 10 M_sun pc-2.
Sigma_SFR (astropy.units.Quantity, optional) – Star-formation rate surface density used when sfh is None (constant SFH). Default is 2.5e-3 M_sun kpc-2 yr-1.
sfh (callable, optional) –
Time-varying star-formation history. When provided, Sigma_SFR is ignored. The callable must accept a numpy array of lookback times in Myr and return an ~astropy.units.Quantity array of SFR surface density in units compatible with M_sun yr-1 kpc-2. Example (Zari et al. 2022 History 1):
```
def sfh_zari1(t_Myr):
    return np.where(t_Myr < 10.0, 4e-3, 2e-3) * au.M_sun/au.kpc**2/au.yr
```
t_max_Myr (float, optional) – Upper integration limit [Myr]. Default is 10000 (10 Gyr).
Z_dust (float, optional) – Dust metallicity relative to solar. Default is 1.0.
dust_kind (str, optional) – Dust model key for DustDraine. One of 'Rv31', 'Rv40', 'Rv55'. Default is 'Rv31'.
verbose (bool, optional) – Print summary statistics. Default is True.

Returns

r – Dictionary with keys Jlambda, Jlambda_unatt, tau_perp, J_FUV, J_FUV_unatt, Sigma_FUV, L_FUV_per_SFR, J, J_unatt, L_per_area, L_per_SFR, w_micron, w_angstrom, Sigma_gas, Sigma_SFR, Z_dust, t_max_Myr.

Return type

dict

Mass-to-luminosity Converter#

class pyathena.util.mass_to_lum(model='Padova')[source]#

Simple class for calculating mass-to-light conversion factor

Initialize a mass_to_lum object.

Parameters: model (string) – Name of stellar evolutionary track. Padova - Power-law approximation to stellar luminosity, MS lifetime based on Padova evolutionary track (Bruzual & Charlot 2003). Data taken from Table 1 in Parravano et al. (2003).

Population Synthesis

Contents

Population Synthesis#

Starburst99 Reader#

Mass-to-luminosity Converter#