Lyman-alpha catalogs

Note: The Subhalo_NoSFR HDF5 group contains luminosities based on the exact star particle ages, whereas Subhalo uses a noise reduction technique to reinterpret ages to be more consistent with the star-formation rate (SFR). The simulations employ stochastic sampling to determine when new star particles are created, with comparable mass resolutions for star and gas particles. The age discretization leads to characteristically bursty star-formation histories (SFHs) in marginally resolved haloes. To extend the reliability of halo-by-halo predictions for galaxies at the faint-end of the luminosity function we employ the following SFH smoothing procedure. We first calculate the total mass of young stars (\(<5\,{\rm Myr}\)) in each subhalo and combine this with the instantaneous SFR to define a duration over which we expect these stars to have formed: \(\Delta t_{\rm SFH} = M_{\rm stars}(<5\,{\rm Myr})/{\rm SFR}\). In the event that this duration is longer than \(5\,{\rm Myr}\) we reassign the young stars as having all been formed with a constant SFR with the expected ages. This is done with 100 equal age bins in the interval \([0, \Delta t_{\rm SFH}]\), retaining the mass and metallicity distribution of each stellar population. This prescription helps us to smooth out an artificial bump around \(M_{1500} \sim -15\) or \(L_\alpha \sim 10^{41}\,{\rm erg\,s}^{-1}\), corresponding to newly spawned star particles in low mass haloes, while the rest of the luminosities are essentially unaltered.

import numpy as np
import h5py

c = 2.99792458e10          # Speed of light [cm/s]
h = 6.626069573e-27        # Planck's constant [erg/s]
pc = 3.085677581467192e18  # Units: 1 pc  = 3e18 cm
lambda_1500 = 1500. * angstrom # Continuum wavelength [cm]
R_10pc = 10. * pc              # Reference distance for continuum [cm]
fnu_1500_fac = lambda_1500**2 / (4. * np.pi * c * R_10pc**2 * angstrom)
E_AGN = 5.29e-11 # Mean photon energy [erg]
E_Lya = h * c / lambda_1216 # Lyman-alpha energy [erg]

def read_Lya(snap=70, sim_dir='Thesan-1'):
    Lya_filename = f'{sim_dir}/postprocessing/Lya/Lya_{snap:03d}.hdf5'
    with h5py.File(Lya_filename, 'r') as f:
        g = f['Header']
        s = {'snap': snap}
        for key in ['BoxSize', 'EscapeFraction', 'Omega0', 'OmegaBaryon', 'OmegaLambda', 'HubbleParam', 'Time', 'Redshift']:
            s[key] = g.attrs[key]
        h = g.attrs['HubbleParam']
        a = g.attrs['Time']
        z = g.attrs['Redshift']
        UnitLength_in_cm = g.attrs['UnitLength_in_cm']
        UnitMass_in_g = g.attrs['UnitMass_in_g']
        UnitVelocity_in_cm_per_s = g.attrs['UnitVelocity_in_cm_per_s']
        UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s
        UnitEnergy_in_cgs = UnitMass_in_g * UnitVelocity_in_cm_per_s * UnitVelocity_in_cm_per_s
        UnitLum_in_cgs = UnitEnergy_in_cgs / UnitTime_in_s
        length_to_cgs = a * UnitLength_in_cm / h
        s['BoxSize_cMpc'] = s['BoxSize'] * length_to_cgs / (a * Mpc) # Box size [cMpc]
        s['V_box_cMpc3'] = s['BoxSize_cMpc']**3 # Box volume [cMpc^3]
        s['Lya'] = UnitLum_in_cgs * f['Subhalo']['LyaLum'][:].astype(np.float64) # Lya = LyaCol + LyaRec + LyaStars
        s['L1500'] = UnitLum_in_cgs * f['Subhalo']['1500LumStars'][:].astype(np.float64) # erg/s/Angstrom
        s['M1500'] = -2.5 * np.log10(fnu_1500_fac * s['L1500']) - 48.6 # Continuum absolute magnitude
        s['NdotAGN'] = UnitLum_in_cgs * f['Subhalo']['IonLumAGN'][:].astype(np.float64) / E_AGN # AGN ionizing photon rate [photons/s]
        s['NdotStars'] = UnitLum_in_cgs * f['Subhalo']['LyaLumStars'][:].astype(np.float64) / (0.68 * E_Lya * (1. - s['EscapeFraction'])) # Star ionizing photon rate [photons/s]
    return s