3.1. Halo Occupation Fraction

The halo occupation fraction presents the distribution of dark matter halos with chemically enriched starparticles at a particular redshift, as a function of the halo mass. We present the halo occupation fraction in the fiducial model Ch 18 512 MDW in Figure 1, including mass bins of log(M _h/M_ʘ) ₌ 0.1. It is worth noting that we only take into account halos above the mass resolution M _h,min = 1.48 x 10⁹M_ʘ.

Fig. 1 Halo occupation fraction in the simulation Ch 18 512 MDW (our fiducial run). The distributions show the percentage of dark matter halos containing formed stars at z = 8, 7, and 6 (light, navy, and dark blue lines, respectively). The color figure can be viewed online. v_fid = 600 km/s, and match the observations at the time when ^{García et al. (2017b)} was published. Hence, the SFR (z = 8) and the stellar mass function are not predictions of the model. 

The distribution in Figure 1 reaches its maximum occupation at M _h ≈ 1 ≈ 10¹⁰M_ʘ for z = 8, 7, and 6, which is consistent with the hierarchical model of structure formation. Dark matter halos with masses below 10⁹M_ʘ contain less than 470 dark matter particles, then they cannot be considered virialized, and the star formation in such regions is disfavoured for two reasons: first, they are unresolved, and they do not efficiently experience atomic cooling; hence, gas is collisionally excited and less likely to form structure in the dark matter wells until it cools down. On the contrary, halos with large masses (≥ 10¹⁰M_ʘ) preferentially form stars, but they are rare on the simulations at high redshift. Although scarce, dark matter halos with large masses present a non-negligible occupation of chemically enriched stellar populations.

At z = 6 (dark blue line), when the Universe has evolved for ≈ 360 Myr from the start of the simulations (z = 8), the highest probability of finding fully occupied halos occurs at M _h ≈ 1 x 10¹⁰M_ʘ. The latter result indicates that halos in the simulation have grown in mass during this period and, consequently, also the cosmic star formation rate.

3.2. Stellar-to-Halo Mass Function

Another observable that we check in our models is known as the galaxy stellar-to-halo virial mass function, presented in Figure 2 for the fiducial run Ch 18 512 MDW, at z = 8 and 6.

Fig. 2 Galaxy stellar-to-halo virial mass function in our reference model Ch 18 512 MDW. The graph displays the distribution of masses above the threshold for this simulation box-size at redshifts z = 8 and 6 (corresponding to light and dark blue stars, respectively). For comparison, the predicted trends from UniverseMachine (^{Behroozi et al. 2020}) are shown at z = 8 (dark red circles) and z = 6 (red dots). Grey dashed lines indicate constant stellar functions f _∗ = M _∗/M _h, ranging from 10⁻⁴ (bottom) to 10⁻¹ (top). The color figure can be viewed online. 

Figure 2 presents the distribution of the stellar-tohalo mass for each galaxy in the fiducial simulation above the mass resolution threshold (> 10⁹M_ʘ or equivalently, more than 470 dark matter particles). Simulated galaxies follow a trend of 10⁻², and at the high mass end, f _∗ moves to 10⁻¹. f _∗ grows with redshift, with more galaxies with large stellar masses at z = 6. This result is also seen in numerical simulations such as Astraeus (^{Hutter et al. 2021}), with M _{h =} 10^9.5−10M_ʘ corresponding to stellar masses ≈ 10^7.5−8M_ʘ during the epoch of reionization. Moreover, ^{O’Shea et al. (2015)} display this distribution (Figure 2, left), but they stack all the galaxies in their realizations up to their final redshifts. The Renaissance simulations describe galaxies that are evolving during the progress of reionization; therefore, they account for halo masses down to 7 × 10⁶M_ʘ (halos with formed stars). Their plot is not directly comparable with Figure 2 because the mass range covered by their simulations differs from ours. Still, galaxies around 10⁹Min their work show the same difference of two orders of magnitude in M _h/M _s. Recent observations by ^{Stefanon et al. (2021)} reveal the same ratio of two orders of magnitude for stellar-to-halo mass ratios (see Table 7 in their work). They emphasize that there is no significant evolution in the observed stellar-to-halo mass function for galaxies in the first Gyr of cosmic time. The latter conclusion is consistent with Figure 2. Instead, predictions by ^{Behroozi et al. (2020)} are an order of magnitude below the trend of our fiducial simulation, at z = 8 and 6. However, we cannot provide a clear explanation for this discrepancy since UniverseMachine is calibrated to resolve virialized halos up to z ≈ 15, but the ratio predicted by their simulation is off compared with other simulations at z = 6 - 8 (including ours).

It is worth mentioning that we only calculate the stellar-to-halo mass function for the fiducial model Ch 18 512 MDW because there is minimal variation in the range of masses resolved by our simulations due to our small boxes (12, 18, and 25 cMpc/h). Thus, the trend shown in Figure 2 is barely affected by changes in the box size.

3.3. Galaxy Stellar Mass Function

One can also count the number density of galaxies formed inside the virialized halos per unit volume _V per stellar mass bin ∆M (^{Weigel et al. 2016}). This observable is known as the galaxy stellar mass function, and it is given by:

The analytical form that describes the galaxy stellar mass function is commonly described using a Schechter function (^{Schechter 1976}), as follows:

The exponential term in the expression above shows the evolution for the high- and low-mass and the fore-most right term a power-law behavior as a function of the stellar mass.

We present the best fit parameters and corresponding errors for the galaxy stellar mass functions at z = 8, 7, and 6 for three of our models with same configuration (MDW and no molecular cooling) with box sizes of 12, 18, and 25 cMpc/h, in Table 2. We derive these parameters with an adapted version from python routine emcee. Table 2 also shows the Schechter function parameters from ^{Duncan et al. (2014)}, ^{Grazian et al. (2015)}, ^{Bhatawdekar et al. (2019)} and ^{Stefanon et al. (2021)}.

One highlight from Table 2 is that the slope of the galaxy stellar mass function at z = 6 - 8 remains constant, and in all the cases presented, is close to the value −2. These findings are consistent with the observational constraints also shown in Table 2.

The evolution of the galaxy stellar mass function from z ₌ 8 to 6 is displayed in Figure 3. Remarkably, the galaxy stellar mass function Φ at z = 8 was set to match observations by ^{Song et al. (2016)} and to calibrate the mass loading factor of the simulations v_fid. However, the theoretical trends agree well with the best fit at z = 8 by ^{Bhatawdekar et al. (2019)} and with the most recent Spitzer/IRAC observations by ^{Stefanon et al. (2021)}, which is quite reassuring since the observational detections came after our simulations.

Fig. 3 Simulated galaxy stellar mass function at z = 8, 7 and 6 in the top, middle and bottom panels. Theoretical predictions from our models are compared with observations by ^{Stefanon et al. (2021)} as black diamonds, ^{Song et al. (2016)} as orange circles and ^{González et al. (2011)} as grey circles, and the best fits proposed ^{Bhatawdekar et al. (2019)}, ^{Grazian et al. (2015)} and ^{Duncan et al. (2014)} as light blue, olive and black dotted lines, respectively. The parameters of the Schechter function (4) for each model are presented in Table 2. The color figure can be viewed online. 

The predicted galaxy stellar mass functions at z = 8, 7, and 6 are compatible with the observational data at high z. Nonetheless, the simulations slightly differ from the galaxy stellar mass function reported by ^{Song et al.(2016)} and ^{González et al. (2011)} in the high mass end at z = 6, mainly because massive galactic halos are scarce in the simulations at these redshifts; thus, high-mass galaxies are rare. Larger synthetic boxes could alleviate this mass bias (for instance, Flares or UniverseMachine).

On the other hand, results with Flares (^{Lovell et al. 2021}; ^{Vijayan et al. 2021}, supported by a vast cosmological box, that at z = 5 and 6, reach stellar masses up to 10¹¹M_ʘ), claim that the galaxy stellar mass function must be fitted with a double-slope Schechter function with a knee at M_{s =} 10¹⁰M_ʘ. They back the latter argument with recent observational constrains from ^{Stefanon et al. (2021)}. Nevertheless, our small boxes do not allow us to test this range of mass (M_s,max ≈ 10^9.8M_ʘ in our largest realization; hence, we keep a single-slope fit). One of the motivations for the ^{Lovell et al. (2021)} work was to extend the range of stellar masses and the number of resolved galaxies reached by the Eagle simulation. Our cosmological runs have similar modules, box sizes, and an analogous configuration as Eagle; therefore, tests related to galaxy observables must be done with Eagle -or equivalent hydrodynamical simulations. Instead, this particular conclusion derived by ^{Lovell et al. (2021)} is out of the reach of our simulations.

3.4. Halo Mass Function

In order to characterize the galaxies in the simulations, the halo mass functions at redshifts z = 8, 7, and 6 are presented in Figure 4, in light, navy, and dark blue lines, respectively. This quantity is computed with systems above the mass resolution limit (a resolved halo in the simulation contains ≈ 470 dark matter particles, or equivalently, a minimum mass M _h,min = 1.48 × 10⁹M_ʘ for boxes of 18 Mpc/h side length). Whenever a galactic halo is below this mass threshold, the object is considered unresolved and is excluded from the statistics.

Fig. 4 Halo mass functions at z = 6, 7, and 8 in dark, navy, and light blue. We include shaded regions corresponding to Poisson errors for the fiducial run Ch 18 512 MDW. As a reference, we show a constant power-law slope α = −2 for galaxies at high redshift, as a dotted black line. The color figure can be viewed online. 

The evolution of the mass function for dark matter halos is computed only with the fiducial model Ch 18 512 MDW since the number of resolved galaxies is almost independent of the feedback mechanisms or the cooling processes implemented in the simulations. As a reference, a dotted black line is presented on top of our predictions in Figure 4, indicating a constant faint-end slope α = −2.

Theoretically, the number density of halos follows the relations d(log⁡N)d(log⁡Mh)=-1, as shown in Figure 4, and N∝Mhα+1, leading to a faint-end slope α = −2². Although the set of simulations presented in this work is unable to provide a direct prediction of the power slope due to the narrow range of halo masses and the small box sizes of the simulations, the curves in Figure 4 show a trend consistent with a power-law slope of −2 at z ≈ 6-8.

We support the latter claim based on the results described in previous sections. The stellar-to-halo mass ratio (Figure 2) shows little evolution of the mass ratio at z ₌ 8 to 6 (regardless of the increasing number of halos that form galaxies with time). On the other hand, the galaxy stellar mass function (Table 2 and Figure 3) indicates that the slope is close to −2 during the time frame described by the simulations. Since the halo masses, M _h are two orders of magnitude larger than the stellar masses M _s -this ratio stays constant at the tail of the reionizationand the slope for the stellar mass function is −2, with almost no variation in time, the value of the slope of the halo mass function is consistent with −2.

3.5. Star Formation Rate Density

The star formation rate is the mass of the new stars in the simulation, measured in the total volume per year. It is commonly assessed by galaxy surveys or derived from studies with luminosity functions. There are two ways to compute the cosmic star formation rate in the numerical runs: (i) adding up the star formation of each gas particle, per comoving volume V; or (ii) recovering the SFR estimate for galaxy groups from the FoF catalog.

Figure 5 shows the cosmic star formation rate in our simulations at 4 < z < 8. The left panel shows the total SFR (including contributions from all the collapsed objects inside the box). Conversely, the right panel displays the same observable, but this time, applying a cut in mass; thus, only the most luminous galaxies are taken into account in the calculation. This mass threshold responds to the resolution achieved by our telescopes that only detect the most luminous galaxies (in particular, at high redshift). Current instruments do not detect the faintest objects; therefore, their SFR cannot be inferred with that method. Consequently, there is an excess in the star formation rate predicted by the simulations -on the left- to observational data. The discrepancy between the calculated and the observed SFR is corrected in the right panel by imposing a luminosity cut M_UV < −17 (corresponding to a minimum SFR > 0.331 M_ʘ/yr and the absolute magnitude set by Hubble observations). When we impose the latter criterion to the simulations, the predicted SFR agrees well with the observations to-date (except for data points measured by ^{Steidel et al. (1999)} and ^{Ouchi et al. (2004)} that do not account for dust corrections).

Fig. 5 Cosmic star formation rate density in the redshift range of 4 < z < 8. The predictions from the simulations are compared with observations from ^{Bhatawdekar et al. (2019)} shown as brown diamonts, ^{Driver et al. (2018)} as yellow squares, ^{Bouwens et al. (2015)} as orange circles (with dust corrections), ^{Cucciati et al. (2012)} as an olive pentagon, ^{Hildebrandt et al. (2010)} as a green inverted triangle, ^{Bouwens et al. (2009)} as a pink square, ^{Ouchi et al. (2004}) as cyan triangles and ^{Steidel et al. (1999)} as a grey diamond. In the left panel, the SFR in the simulations is computed including all objects in the box per unit volume. On the right, the observable is limited to masses with the luminosity cut of M_UV < −17, equivalent to a minimum SFR of > 0.331M_ʘ/yr, following the observational constraints of our current telescopes. The color figure can be viewed online. 

Interestingly, data from ^{Driver et al. (2018)} and ^{Bhatawdekar et al. (2019)} had not been published when our simulations were run, but most of the models are in agreement with these observations.

On the other hand, the cosmic SFR reported by ^{Finkelstein (2016)} -with a corresponding comparison with ^{Madau & Dickinson (2014)}- shows an increment of 1 dex in their reference model, consistent with findings from this work with the mass cut M_UV < −17 (right panel of Figure 5).

Figure 6 shows a compilation of theoretical predictions for the cosmic star formation rate by UniverseMachine(^{Behroozietal.2020}), L-Galaxies 2020(^{Henriques et al. 2020}; ^{Yates et al. 2021a}), Flares (^{Lovell et al. 2021}; ^{Vijayan et al. 2021}), TNG100 (^{Nelson et al. 2018}; ^{Pillepich et al. 2018}; ^{Naiman et al. 2018}; ^{Marinacci et al. 2018}; ^{Springel et al. 2018}), Eagle (^{Crain et al. 2015}; ^{Schaye et al. 2015}) and our models.

Fig. 6 Cosmic star formation rate density in the redshift range of 4 < z < 8. The predictions from our models are compared with the calculated SFR with the UniverseMachine (^{Behroozi et al. 2020}), L-Galaxies 2020 (^{Henriques et al. 2020}; ^{Yates et al. 2021a}), Flares (^{Lovell et al. 2021}; ^{Vijayan et al. 2021}), TNG100 (^{Nelson et al. 2018}; ^{Pillepich et al. 2018}; ^{Naiman et al. 2018}; ^{Marinacci et al. 2018}; ^{Springel et al. 2018}), and the Eagle simulations (^{Crain et al. 2015}; ^{Schaye et al. 2015}). The observable is limited to observed galaxies with absolute magnitude cut of M_UV < −17, equivalent to a minimum SFR of > 0.331M_ʘ/yr. The color figure can be viewed online. 

Except for Ch 18 512 EDW, all our simulations show an excess of the SFR when compared with other theoretical models. As mentioned in a previous section, the simulation with the closest configuration to ours is Eagle, and this explains why their calculated SFR is compatible with our runs with the energy-driven winds prescription for the supernova outflows. This is also true with TNG100, which agrees well with EDW runs, but it is always below the prediction with MDW realizations. This outcome from our simulations is promising, since the winds implemented in the IllustrisTNG project have much more complex dynamics than ours: the velocity of the galactic winds v _w also depends on z (suppressing the efficiency of winds with the Hubble factor, ^{Pillepich et al. 2018}), but scales in the same way as our winds with the virial halo mass. Another remarkable difference in the IllustrisTNG winds is that the outflow mass loading is a non-monotonic function of the galaxy stellar mass (^{Nelson et al. 2018}). We do not account for such dependence in our models. Finally, TNG introduces an improved mechanism for the AGN feedback, even for a low accretion rate, whereas this work does not account for AGN feedback.

A slightly different scenario is drawn with the two versions of L-Galaxies 2020. Both configurations match the observed SFR at low and intermediate redshifts because their semi-analytical models were built to follow the chemical enrichment at late times, not during reionization. Besides, their models heavily rely on observations from damped Lyman systems (DLAs) that cannot be extended to the redshifts of the EoR (^{García et al. 2017b}). The gap in the calculated SFR grows among our runs and large cosmological simulations, particularly towards z → 8. The UniverseMachine and Flares theoretical models aim to correct the UV luminosity function and to provide a forecast for future wide-field surveys, as the Nancy Roman (previously known as Wfirst), Euclid or JWST. Findings from these large volume boxes, with broader redshift ranges, are poorly constrained by periodic hydrodynamical simulations due to their limited volume, and consequently, reduced number of massive galaxies.

It is worth mentioning that the different supernova feedback mechanisms play a dominant role in the evolution of the star formation rate. Figures 5 and 6 show lower values for the SFR with energy-driven winds than momentum-driven winds (EDW and MDW, respectively), indicating that the former mechanism is more efficient at quenching the SFR because it prevents the overcooling present in the latter case, which leads to an excess of the number of stars that would form during a time interval of ≈ 1 billion years (z = 8 to 4). Notably, once many star formation events occur in the simulation, the stochastic SFR converges to its continuous history, and the galaxies grow in size and mass through this physical scheme.

Finally, it is interesting to study the ratio between the observed M_UV < −17 and the total cosmic star formation rates in the different realizations considered in Figure 5. Although, this is not an observable, it reflects how the mass cut affects the overall SFR in the synthetic realizations.

Figure 7 shows the predicted ratio from the simulations and UniverseMachine. Most of our configurations show a flat trend over the entire redshift range (except for the Ch 25 512 MDW mol, which is, in fact, the simulation with the lowest resolution). This result leads to the conclusion that ≈ 2 out of 5 simulated galaxies are about the luminosity cut of M_UV < −17, and this ratio does not evolve from z = 8 to 4. Conversely, ^{Behroozi et al. (2020)} show a rapid increment in the SFRD(M_UV < −17) / Total SFRD, from 0.5 to 0.9 for z = 8 to 4, respectively. The latter is consistent with a change of 1 dex in Figure 5 -right panel-, indicating that the vast majority of the stars formed during a time frame of 900 Myr are above the luminosity cut M_UV < −17. Beyond the percentage inferred from our simulations, observations show an increasing number of bright galaxies at the tail of reionization and indirectly confirm the predicted ratio by ^{Behroozi et al. (2020)}.

Fig. 7 Predicted ratio between the observed M_UV < −17 and total cosmic star formation rates from the simulations, and comparison with UniverseMachine ratio. The blue band indicates the error, and the dashed black line the mean value from the set of simulations. The color figure can be viewed online. 

3.6. Chemical Enrichment in the Simulations

One of the strengths of this model is the self-consistent chemical enrichment implementation. The metals’ production, spread, and mixing to the cincum- and intergalactic medium come from the assumed stellar lifetime function, stellar yields, and initial mass function.

It is worth noting that there are no measurements of the cosmic mass densities for any metal. The only tight constraint is that elements except for H and He should account for about 1% of the baryonic content in the Universe. This issue becomes even more challenging at high redshift when indirect methods are less precise to quantify the amount of any element. However, astronomers can estimate lower limits for the percentage of individual metals in the total census by measuring the total mass density for metal ions in the IGM (see work from ^{García et al. (2017b)} with CII and CIV, and the detections by ^{Codoreanu et al. (2018)} on SiII and SiIV). Also, an approximative evaluation of the relative metallicity to hydrogen with damped Lyα systems is possible, but this method does not provide any observational constraints for O or Si (among other metals).

Detections of absorption lines from ^{Codoreanu et al. (2018)}, ^{Meyer et al. (2019)}, and ^{Cooper et al. (2019)} show that we can set a lower limit for the mass density of silicon and oxygen through the reconstruction of the cosmic mass densities of their corresponding metal ions. Besides that, the low-to-high ionization ratio is an independent proxy for the end of reionization and reveals the gas state in the IGM. Both oxygen and silicon have observable ionization states that are exhibited in the spectra of high redshift quasars (OI, SiII and SiIV), redward from the Lyman α emission, and could provide complementary constraints to the metal enrichment, apart from carbon.

Figure 8 shows the cosmic evolution of the oxygen and silicon mass densities, from z = 8 to z = 4 when the chemical pollution has been occurring for about a billion years in the Universe from stars and supernovae. The cosmic density Ω as a function of _z is obtained by summing the amount of each metal in all gas particles inside the simulated box. Finally, this calculation is divided by the comoving volume.

Fig. 8 Theoretical predictions for the total cosmological mass densities for oxygen (Ω_O) and silicon (Ω_Si) in the left and right panels, respectively. The color figure can be viewed online. 

In addition to the metal mass densities predicted from our models, a comparison with L-Galaxies 2020 for O and Si is presented for their default and modified model (DM and MM, respectively; ^{Yates et al. 2021b}). Their normalization is similar because the overall amount of cosmic star formation is slightly higher in their modified model than the default setup, despite their distinctive mechanisms introduced to enrich the CGM/IGM; thus, the amount of each element produced overall is similar. For further details on their chemical enrichment modeling, see ^{Yates et al. (2013)}. Both models differ by around one order of magnitude from our mass densities due to three main differences among their models and ours: (i) their predicted SFRD are lower at all redshifts. Therefore, it is expected that chemical pollution is less effective in this time frame. (ii) ^{Yates et al. (2021b)} assume different metal yields than the ones imposed on our set of simulations. The former is around 0.03 -in order to match late metallicities measured with DLAs-. Instead, the metal yield in all our models has a fixed value of 0.02. (iii) Their models include metal outflows released and spread by SNe-II, SNe-Ia, and AGB stars. In our simulations, neither AGB nor AGN are predominant in the feedback mechanisms. Thus fewer processes prevent the outburst of material.

Our predictions are also contrasted with the cosmic mass density calculated with TNG100. Their trends show a much faster evolution than the prediction from L-Galaxies 2020 and ours. These results are consistent with the SFRD exhibited by TNG100, on top of a sophisticated set of stellar yield tables (Table 2 from ^{Pillepich et al. 2018}).

From the observational side, a compelling test arises from the evolution of the mean metallicity in the Universe Figure 14 in ^{Madau & Dickinson (2014)}, and Figure 8. In ^{Madau & Dickinson (2014)} an increment of 1 dex to the solar metallicity is present (under the assumption that the mass of heavy metals per baryon density produced over the cosmic history with a given SFR model and an IMF-averaged yield of y = 0.02), consistent with the predictions from the simulations at z = 4 - 7.

Finally, a correlation between the SFR in the simulations (mainly driven by the feedback processes of the gas) and the cosmic mass densities of these elements is found. Both metals in Figure 8 show a slight boost at all redshifts when the feedback mechanism is MDW. As mentioned above, the latter feedback prescription is more effective in producing an overcooling of the gas in the CGM, leading to a larger star-formation in the simulations. Hence, more metals are generated and expelled from the galaxies.

It is fair to conclude that the metal pollution scheme that occurs inside the synthetic realizations agrees with current limits obtained with metal ions of oxygen and silicon at high redshift ^{Codoreanu et al. (2018)}, despite the limited number of absorption lines detected to date.