A quantification of the effects using different stellar population synthesis models for epoch of reionization

The dark matter simulation and merger trees adopted are from the Millennium-II simulation (MS-II, Boylan-Kolchin et al., 2009), which was run with an updated version of the GADGET code (Springel, 2005) i.e. GADGET-3. The original box size of MS-II is $100\,{\rm Mpc}/h$ with $2160^{3}$ dark matter particles, and each particle is with mass of $6.89\times 10^{6}\,h^{-1}\,{\rm M_{\odot}}$. The halos were identified using the Friends-Of-Friends (FOF) algorithm (More et al., 2011), then the SUBFIND algorithm was applied to identify the self-bound substructures within each FOF group. The halos are with at least 20 particles, i.e. the minimal halo mass $1.38\times 10^{8}\,h^{-1}\,{\rm M_{\odot}}$. The simulations are scaled to the Planck cosmology (Planck Collaboration et al., 2020) following the procedure of Angulo & Hilbert (2015), with the cosmological parameters $\Omega_{m}=0.315$, $\Omega_{b}=0.049$, $\Omega_{\Lambda}=0.685$, $h=0.673$, $\sigma_{8}=0.826$ and $n_{s}=0.965$.

MS-II simulation has outputs of 68 snapshots from $z=127$ to 0, while 22 of them are at $z\geq 6$. The merger trees were constructed with the sub-halos found in these snapshots, in this case the merger trees of MS-II simulation include $\sim$ 590 million sub-halos in total (Boylan-Kolchin et al., 2009). The more or less snapshots should not obviously change the conclusions presented in this paper. We test it with the low resolution simulations described in Ma et al. (2023).

We adopt the public SAM model L-Galaxies 2020 (Henriques et al., 2020) to evolve the galaxies with the merger trees from MS-II simulation. This model has included almost all the physical processes related to galaxy formation, such as gas cooling, star formation, galaxy merger, supernovae and AGN feedback. The application of L-Galaxies 2020 on the high-$z$ galaxies and EoR has been explored in Ma et al. (2023). The L-Galaxies 2020 model is the updated version of Henriques et al. (2015), which has a few differences with respect to the latter. In the new model, the galactic discs are spatially resolved by dividing the discs into $12$ concentric rings with radii $r_{i}=0.01\times 2^{i}h^{-1}\,{\rm kpc},\,i=0,...,11$ (Fu et al., 2013). All the properties and physical processes of discs, such as star formation, chemical enrichment and gas ejection, are evolved for each ring separately. The star formation is linked to the H₂ abundance of each ring that depending on the metallicity of gas (Krumholz et al., 2009; McKee & Krumholz, 2010).

We adopt the default values for the free parameters in the L-Galaxies 2020 code, which are the best fit values with the observations of galaxies at low-$z$ and the merger trees from MS-II simulation (Henriques et al., 2020). We only modify the input SPS model adopted in the L-Galaxies 2020 code.

We adopt four popular SPS models in this paper, i.e. BC03 (Bruzual & Charlot, 2003), M05 (Maraston, 2005), YEPS (Zhang et al., 2005) and BPASS (Stanway & Eldridge, 2018). The parameters of these SPS models, including the wavelength range, age and metallicity, are summarized in Tab. 1. To make consistent comparison, we employ the results of these models with the same Salpeter stellar initial mass function (IMF, Salpeter, 1955) with a power index $\alpha$ of 1.35 and maximum initial mass 100 $\rm M_{\odot}$. For the YEPS and BPASS models, we consider both the ones including binary stars, which are named as YEPS BS and BPASS BS respectively, and the ones with only single stars, which are named as YEPS SS and BPASS SS respectively.

As a reference, in Fig. 1 we present the initial SEDs from different SPS models at four metallicities ($Z=10^{-4}$, $10^{-3}$, $4\times 10^{-3}$ and 0.02) and four ages (1 Myr, 10 Myr, 100 Myr and 1 Gyr). Note that the SEDs in Fig. 1 are only in the frequency range $[6.2,200]\,\rm eV$ (i.e. UV band), while the Fig. 6 in Appendix A shows one example of the full SEDs. As showed in Fig. 1, the initial SEDs from all four SPS models are not very sensitive to the metallicities, while the amplitudes reduce obviously with the increasing ages. The different SPS models present roughly consistent SEDs, especially at $h_{\rm P}\nu<13.6\,\rm eV$, while some differences are still visible. The inclusion of binary stars in the YEPS and BPASS models (i.e. YEPS BS and BPASS BS) leads to harder SEDs at $h_{\rm P}\nu>13.6\,\rm eV$ compared to those with single stars (i.e. YEPS SS and BPASS SS), and thus expects to increase the emission of ionizing photons. In the following, we will briefly describe these SPS models. For more comparisons about the SPS models, one can also refer to e.g. Chen et al. (2010) and Han & Han (2019).

The SEDs of galaxies from L-Galaxies 2020 are computed with the SPS results (i.e. the Fig. 1) and the history of star formation and metal enrichment within each galaxy. More specifically, with the information of star formation and metal enrichment history produced, the L-Galaxies 2020 post-processes the SEDs of each galaxy by linearly interpolated the input SEDs from different SPS models at specific ages and metallicities, and then multiply the star mass. The final SED of galaxies is the sum of the results at all ages. Fig. 2 shows the average rest-frame SEDs of high-$z$ galaxies within the same halo mass range obtained from L-Galaxies 2020 with four SPS models at different $z$s. To compare with the input SEDs showed in Fig. 1, the results of Fig. 2 are normalized by the stellar mass of galaxies. By comparing with the Fig. 1, we can see that the SEDs of high-$z$ galaxies within UV band are dominated by the young stars, i.e. those with age $<$ 100 Myr. Due to the same reason, the SEDs of high-$z$ galaxies do not evolve too much with the decreasing redshift. The SEDs of massive halos are slightly lower than the less massive ones, since the latter ones have higher ratio of star formation rate over stellar mass (Henriques et al., 2020).

Four SPS models predict similar SEDs of high-$z$ galaxies at $h_{\rm p}\nu<13.6\,\rm eV$, while some differences are obvious at higher energy band. The BPASS model shows obviously higher amplitudes of SEDs than other models at $h_{\rm p}\nu=13.6-50\,\rm eV$. Within the same band, the BC03 model is slightly lower than the BPASS model, while the YEPS model has the lowest amplitude, and the M05 model is between BC03 and YEPS models. Differently, the YEPS model has the highest SED of galaxies at $h_{\rm p}\nu>50\,\rm eV$ and $z\geq 7$, while becomes similar to the BPASS model at $z=6$. Within the same band, the SEDs of galaxies with M05 model are roughly consistent with the BPASS model at $z\geq 7$, while obviously lower than the latter at $z=6$. The SEDs of galaxies with BC03 model at $z\geq 7$ are similar to those with M05 and BPASS models at $M_{\rm vir}<10^{10}\,\rm M_{\odot}$, while become similar to those with YEPS models with the increasing halo mass e.g. at $M_{\rm vir}=10^{11}-10^{12}\,\rm M_{\odot}$. The inclusion of binary stars (BPASS BS) does not obviously change the SEDs of BPASS model at $h_{\rm p}\nu<50\,\rm eV$, while increases the high energy radiation at $h_{\rm p}\nu>50\,\rm eV$. Such effect is weaker on the SEDs of massive halos (e.g. $M_{\rm vir}=10^{11}-10^{12}\,\rm M_{\odot}$) than the less massive ones, due to the increasing contributions of old stars (age $>$ 100 Myr). With the same reason, the SEDs of galaxies at $z=6$ has no significant features of binary stars. Instead, the inclusion of binary stars in YEPS model (YEPS BS) has no obvious effects on the SEDs of high-$z$ galaxies, except slightly higher amplitude at $h_{\rm p}\nu>40\,\rm eV$ and $z=6$.

Note that, due to the abundant neutral hydrogen during EoR, it will be hard to measure the SEDs at $>13.6\,\rm eV$, although the uncertainties of different SPS models are mostly at such band. The measurements and comparisons at rest-frame UV band e.g. the UV luminosity function ($\phi$ showed in Fig. 3) can exclude the uncertainties of physical processes except the SPS models, while the results of SED fitting, e.g. the ionizing photon production efficiency $\zeta_{\rm ion}$, can help to distinguish the different SPS models (Seeyave et al., 2023).

Fig. 3 shows the UV luminosity function $\phi$ at the rest-frame wavelength $\lambda=1600$ Å from four SPS models at different $z$s. The absolute magnitude of galaxy luminosity at $\lambda=1600$ Å is computed by:

where $F_{1600}$ is the brightness of galaxy at $\lambda=1600$Å, and $R=10\,\rm pc$. As a comparison, we also present some of recent observations of $\phi$ from HST (Bouwens et al., 2021) and JWST (Donnan et al., 2023; Adams et al., 2023; Finkelstein et al., 2023; Donnan et al., 2024) telescopes. Note that the results of Bouwens et al. (2021) are at $\lambda=1600$Å, while others that observed by the JWST telescope are at $\lambda=1500$ Å. Such differences are not significant on the UV luminosity functions. Table 2 shows the $\chi^{2}$ of $\phi$ from different SPS models compared to the observations at three $z$s, which is defined as:

where $\phi_{\rm obs}$ is the observational $\phi$, $\phi_{\rm sim}$ is the $\phi$ from simulations, and $\sigma_{\rm obs}$ is the 1-$\sigma$ error of $\phi_{\rm obs}$. Due to the lack of bright galaxies from MS-II simulation, the $\chi^{2}$ is computed only for the galaxies with $M_{1600,\rm AB}>-21$.

The UV luminosities of galaxies can be reduced by the extinction models of dust and molecular clouds. We test the dust model within L-Galaxies 2020, which does not change too much on the Fig. 3, especially at $z>7$. This is different to previous studies, e.g. the results of L-Galaxies 2015 (Clay et al., 2015). Although the dust model is not included to precisely match the results from the MS-II simulation and L-Galaxies 2020 code with the observations, it indeed affects the UV luminosity of the bright galaxies at low redshifts (Yung et al., 2020a; Bhagwat et al., 2023), but not too much on the ones at high redshifts and the faint ones.

The UV luminosity function $\phi$ from four SPS models are roughly consistent with the observations at six $z$s. The differences of four SPS models on $\phi$ are $\lesssim 0.5$ dex, smaller than the uncertainties (i.e. error bars) of the current measurements of $\phi$. Specifically, the BPASS model has higher amplitudes of $\phi$ than other models. The BC03 model is globally similar to the M05 model, which $\phi$s have amplitudes lower than the YEPS and BPASS models. The YEPS model is similar to the BPASS model at $M_{1600,\rm AB}>-18$, while closer to the BC03 and M05 models at $M_{1600,\rm AB}<-18$. The inclusion of binary stars has negligible effects on $\phi$, both for the YEPS and BPASS models, consistent with the SEDs of high-$z$ galaxies showed in Fig. 2. With the $\chi^{2}$ showed in Table 2, at $z=10$ the YEPS BS model fits better with the observations, then is the BPASS SS models. At $z=9$, the BPASS SS model is the best fit one. At $z=8$, the best fit model is the BPASS BS model, then is the BPASS SS model.

To properly compute the number of ionizing photons from high-$z$ galaxies, we rerun the L-Galaxies 2020 SAM simulations with the integrated SED (iSED) over the age of stars (i.e. the time from the birth of stars to the output $z$s) for four SPS models. The iSED can easily include the effects of star formation history (Ma et al., 2023). Meanwhile, as showed in Ma et al. (2023), after normalized by the stellar mass the iSED of high-$z$ galaxies is sensitive neither to the galaxy formation models nor to the output $z$s. The number of ionizing photons ($n_{\rm ion}$) from galaxies is then calculated by:

where $L_{\nu}$ is the iSED of galaxies computed by L-Galaxies 2020, $h_{\rm P}$ is the Planck constant, and $\nu$ is the frequency of photons. The integration is done with the full iSED at $h_{\rm p}\nu>13.6\,\rm eV$. The cosmic volume averaged ionizing photons $N_{\rm ion}$ (i.e. number density of ionizing photon) is expressed as:

where the sum $\sum$ is for all the selected galaxies, and $V_{\rm box}$ is the comoving volume of MS-II simulation.

Fig. 4 shows the $N_{\rm ion}$ of galaxies within different halo mass $M_{\rm vir}$ range from four SPS models as functions of $z$s. As a comparison, we also present the minimal budget of ionizing photon number density to fully ionize neutral hydrogen and first ionizing of helium by assuming 75% of baryon is hydrogen and 25% is helium, and show as the horizontal gray lines. From Fig. 4, the ionizing photons are mostly from the galaxies with $M_{\rm vir}>10^{9}\,\rm M_{\odot}$, and $\sim$ half of them are from the galaxies with $M_{\rm vir}>10^{10}\,\rm M_{\odot}$. The galaxies with $M_{\rm vir}<10^{9}\,\rm M_{\odot}$ only slightly increase the $N_{\rm ion}$. This is partly due to the incomplete sample of halos with $M_{\rm vir}<10^{9}\,\rm M_{\odot}$ from MS-II simulation. Meanwhile, the supernovae and radiation feedback also suppresses the star formation on such low mass halos (Hutter et al., 2021; Legrand et al., 2023). In Fig. 7 of Appendix B, we show one sample of the distributions of ionizing photons at $z=7$ versus halo mass.

Since the different SPS models and the binary stars mostly affect the radiation of galaxies at high energy band ($h_{\rm p}\nu>13.6\,\rm eV$) as showed in Fig. 2, the four SPS models have much larger differences on the budget of ionizing photons (i.e. $N_{\rm ion}$) than that on the UV luminosity function showed in Fig. 3. Specifically, the BPASS model has the highest $N_{\rm ion}$, which is $\sim 2$ times that of YEPS model. Comparing with the minimal budget of ionizing photon number density to fully ionize neutral hydrogen and first ionizing of helium, this can lead to redshift difference $\delta z\sim 1$ on the end redshift of EoR. Note that, the precise end redshift of EoR should be computed with the radiative transfer simulations, i.e. the redshift difference estimated might be different due to the ionizing and recombination models. The M05 model has $N_{\rm ion}$ similar to that of YEPS model, which is slightly lower than that of BC03 model. With the inclusion of binary stars, the BPASS model (BPASS BS) can have $\sim 40\%$ more budget of $N_{\rm ion}$ than that with single stars (BPASS SS), consistent with the results of e.g. Ma et al. (2022), while the YEPS model (YEPS BS) has $\sim 28\%$ higher $N_{\rm ion}$ than YEPS SS. The differences on $N_{\rm ion}$ caused by the binary stars are much smaller than that due to different SPS models. Such results are consistent with the SEDs of galaxies showed in Fig. 2.

To compare with the observations, in Fig. 5 we show the ionizing photon production efficiency $\zeta_{\rm ion}$ from four SPS models, which is defined as:

where $\dot{n}_{\rm ion}$ is the ionizing photon emissivity of galaxies. Our results are roughly consistent with the previous studies e.g. Wilkins et al. (2016); Yung et al. (2020b); Seeyave et al. (2023). We only present some results of recent measurements about $\zeta_{\rm ion}$, while for more ones one can refer to the recent paper Simmonds et al. (2024). From Fig. 5, the $\zeta_{\rm ion}$ from our simulations are lower than the measurements by Tang et al. (2023); Whitler et al. (2024), while roughly consistent with other observations.