The fate of baryons in counterfactual universes

The simplest anthropic ensemble is a subset within the multiverse in which only the value of $\Lambda$ varies between the universes. All other key parameters, in the form of dimensionless numbers such as the photon-to-baryon ratio, are assumed to remain unchanged between members of the ensemble. In effect, we can imagine making multiple copies of the universe at high redshifts when $\Lambda$ is so small that its effect is non-existent; we can then scale $\Lambda$ up or down in these universes, which will have negligible effect until $\Lambda$ dominates at some later time. The amount of scaling will modify the cosmological parameters at $z=0$ – which we define for convenience as the time when the temperature of the CMB is $2.725\,{\rm K}$.

The altered $z=0$ cosmological parameters are of interest as they are the practical means by which we will carry out simulations of the counterfactual universes. The dependence on $\alpha$, the scaling factor for $\Lambda$, can be deduced as follows. First note that all present-day physical densities of the different cosmological constituents are unchanged:

where a tilde denotes the modified value of the parameter, and the suffix $j$ can refer to radiation, CDM, or baryons (we assume that the neutrino mass is negligibly small). This arises because $\rho_{j}\propto\Omega_{j}h^{2}$, because defining $z=0$ at fixed temperature fixes $\rho_{r}$, and because we fix the baryon fraction and the number of photons per baryon. In order obtain $\tilde{h}$ and $\tilde{\Omega}_{j}$ separately, we make appropriate use of the Friedmann equation, using the requirement that the expansion rate is not changed at early times where $\Lambda$ is negligible. Suppose we scale $\Lambda$ by $\alpha$ at some early time when the scale factor is $a_{i}$, while maintaining spatial flatness. The subsequent evolution of the Hubble parameter is

where $\epsilon=\rho_{v}/\rho_{\rm tot}$ and $\rho_{v}$ is the vacuum density at $a_{i}$. At early times, $\epsilon$ is small enough that we can set $(1-\alpha\epsilon)$ to unity. Comparing with the reference $\Lambda$CDM model, $H^{2}=H_{0}^{2}(\Omega_{m}a^{-3}+1-\Omega_{m})$, and setting ${\tilde{H}}=H$ at early times, this yields $H_{i}^{2}a_{i}^{3}=\Omega_{m}H_{0}^{2}$. Finally, the fiducial vacuum fraction at $a_{i}$ is

Thus $\epsilon=(1-\Omega_{m})H_{0}^{2}/H_{i}^{2}$ and we can express $\tilde{H}(a)$ in terms of our reference parameters,

which gives the $z=0$ Hubble parameter as

This relation can be obtained directly from the Friedmann equation by noting that we do not change the matter density at $z=0$, but only scale the vacuum density. Since $\rho\propto\Omega H_{0}^{2}$, the result follows. Because the total density remains critical, $\tilde{\Omega}_{m}$ is the matter density divided by the sum of matter and vacuum densities at $a=1$: $\Omega_{m}/[\Omega_{m}+\alpha(1-\Omega_{m})]$, or

This recasting of the cosmological parameters is reminiscent of the ‘separate universe’ technique common in numerical simulations (e.g. Wagner et al. 2015), but with the critical difference that in the latter case the curvature is allowed to vary. Here the assumption is that whatever generates the multiverse ensemble, such as inflation, has a mechanism that guarantees flatness in all cases.

Finally, we also need the normalization $\tilde{\sigma}_{8}$ in this modified cosmology. If the comoving filter length in the reference $\Lambda$CDM model is $R=8\,h^{-1}\,$cMpc, then $\sigma(R)$ needs to have the same value at high redshift in all models. However, $\sigma(R)$ at $z=0$ when the temperature of the CMB is 2.725 K is now altered because the modification of the value of $\Omega_{m}$ will change the linear growth factor (Linder & Cahn, 2007). Furthermore, $R$ is no longer $8\,h^{-1}{\rm cMpc}$ in terms of $\tilde{h}$, but we wish by convention to specify the normalization in terms of $\tilde{\sigma}_{8}$, which is the $z=0$ value of $\sigma$ with a filtering radius of $8\,\tilde{h}^{-1}{\rm cMpc}$. We therefore need to scale the $z=0$ $\sigma(R)$ by a factor that depends on the shape of the power spectrum in order to specify a normalization that keeps the early universe physically identical in all cases. The results of this exercise, together with the dependence of $\tilde{\Omega}_{m}$ and $\tilde{h}$ on $\alpha$, are shown in Figure 1.

For any value of $\alpha$, we can thus obtain an altered set of $z=0$ cosmological parameters, which can then be used to generate cosmological simulations with an unaltered $\Lambda$CDM simulation code. The only complication is the sign of $\alpha$: the scaling formulae are general, and will yield correct results if the high-$z$ cosmological constant is replaced by a negative value, but complications can arise in the subsequent evolution. With $\Lambda<0$, the cosmological expansion will eventually reach a maximum and then undergo recollapse. A standard $N$-body code is likely to fail at this point, since the relation between time and scale factor is normally assumed to be monotonic. Moreover, if $|\Lambda|$ is sufficiently large, the recollapse will occur very early: the minimum temperature will be above $2.725\,$K, so that $z=0$ will not be reached. We can see this from the above expressions, in which $\tilde{h}$ is undefined for $\alpha<-1/\Omega_{m}+1$. These complications are not too difficult to deal with, but we do not address them in this paper, which is concerned purely with the astrophysical impact of a large positive $\Lambda$ in freezing out structure growth and subsequent star formation.

As an alternative to these scaled-$\Lambda$ models, we also conduct a set of simulations with a simpler parameter variation in which we alter only the normalization $\sigma_{8}$, while retaining unchanged the values of the other $\Lambda{\rm CDM}$ parameters. The rationale for this is that a principal effect of raising the value of $\Lambda$ is to cause earlier freezeout and hence yield a lower value of $\sigma_{8}$. This reduced normalization is the principal reason for the operation of Weinberg’s anthropic argument, since it causes a corresponding reduction in the abundance of galaxy-scale dark-matter haloes. It is therefore of interest to look at the effect of a direct change in normalization, as opposed to a change induced as a consequence of changing $\Lambda$. For a given value of $\sigma_{8}$ in the far future, to what degree will the results depend on the evolutionary track by which it was attained? This question motivates us to scale $\sigma_{8}$ in order to obtain an equivalent combination of the cosmological parameters to one where $\Lambda$ is changed. We know that $d\ln\delta/d\ln a\approx\Omega_{m}^{0.55}$ (Linder & Cahn, 2007) and growth freezes when $\Omega_{m}=\Omega_{\Lambda}=0.5$. Suppose we take a universe where we increase $\rho_{\Lambda}$ by a factor $\alpha$: the scale factor where $\rho_{\Lambda}$ = $\rho_{m}$ is then reduced by a factor $\alpha^{1/3}$. Since fluctuations grow linearly with $a(t)$ until this point, we therefore reduce the value of $\sigma_{8}$ by the same factor of $\alpha^{1/3}$ relative to the reference $\Lambda\rm CDM$ universe.

The scaling relations from Section 2.1 allow us to produce sets of cosmological parameters for models where we change either the value of $\Lambda$ or $\sigma_{8}$ from a reference $\Lambda\rm CDM$ universe. The chosen values are summarised in Table 1. The reference $\Lambda{\rm CDM}$ cosmological parameters are adopted from WMAP-9 (Bennett et al., 2013):

where the parameters have their usual definitions. We generate the initial conditions for these universes with MUlti-Scale Initial Conditions (MUSIC) for cosmological simulations (Hahn & Abel, 2011) and evolve them with Enzo (Bryan et al., 2014), using the hydrodynamic solver ZEUS (Stone & Norman, 1992), and an N-body adaptive particle-mesh gravity solver (Efstathiou et al., 1985). We also couple the simulation with the Grackle gas-phase chemistry package (Smith et al., 2017). We used this simulation apparatus in paper 1: Oh et al. (2020), a study of cosmic star formation focusing on galaxies similar to the Milky Way, where we found that reliable results to $z=0$ for a fiducial LCDM simulation could be achieved using a base mesh of $128^{3}$ and four levels of AMR refinement, yielding a spatial resolution of 35 ckpc and a mass resolution of $6.8\times 10^{9}\,M_{\odot}$. In paper 2: Oh et al. (2021), we extended this calculation into the future; this required a number of modifications to the standard code, which we include here.

As listed in Table 1, we consider scaling factors for $\Lambda$ of $\alpha=0,10,50,100$, yielding simulations labelled as EDS, 10x and 100x respectively. The value of $\alpha=0$ implies a flat universe fully dominated by matter, the Einstein–de Sitter (EdS) model (Einstein & de Sitter, 1932). We also reduce $\sigma_{8}$ by a factor of $\beta^{1/3}$, where $\beta=5,10,20,100$. These simulations are named s$\beta$x with a ‘s’ prefix to distinguish them from the simulations with scaled $\Lambda$.

In addition to the cosmological parameters, Table 1 includes two additional columns, giving the comoving box size and starting redshift used to generate the initial conditions for the counterfactual universes. For LCDM, the mass and maximum spatial resolutions are identical to simulation NL in Oh et al. (2021). The comoving box size is specified in units of $\,h^{-1}{\rm cMpc}$, but we change this value according to $h$ in different cosmologies so that it has a fixed value in cMpc, in order to maintain a consistent spatial resolution and an identical particle mass across all simulations. This issue does not arise for simulations with scaled $\sigma_{8}$, since $h$ is the same for all these models; but for EDS, 10x and 100x, $h$ differs significantly between the universes.

The starting redshift of the simulations is also adjusted such that the universes are at the same starting age. We can convert the starting redshift of 99 in the fiducial model to a corresponding value for the counterfactual universes by using the relation

where $\mathrm{H_{0}}$ and $\Omega_{m}$ is the Hubble parameter and matter density parameter at $z=0$ (Peebles, 1993). Again, we do not need to change the starting redshifts for simulations with scaled $\sigma_{8}$. For EDS, 10x and 100x, the starting redshifts are 25.69, 99.03 and 98.75 respectively.

Our simulations adopt the Cen & Ostriker (1992) and Cen & Ostriker (2006) models for star formation and feedback respectively, with modifications made by Smith et al. (2011). An extensive parameter study and calibration of these algorithms within Enzo has been carried out by Oh et al. (2020). We provide here a short summary of the parameters involved in the modelling.

Out of the two setups explored in Oh et al. (2020), we focus on the application of ‘Setup 2’, which employs a timestep-independent conversion of gas into stars. According to the calibrated results, we set $\epsilon=3.0\times 10^{-5}$, $r\_s=1\_1$ and $f_{*}=0.9$, where $\epsilon$ is a factor that relates the amount of feedback energy to the rest mass energy of the star forming gas (see Equation 6 in Oh et al. 2020); $r\_s$ specifies the volume into which the feedback energy is deposited (6 adjacent cells to the star particle for $r\_s=1\_1$); and $f_{*}$ indicates the conversion efficiency of gas mass in a cell into stellar mass (see sections 2.1 and 2.2 in Oh et al. 2020 for details). This set of parameters successfully reproduced the baryonic makeup of a Milky-way sized halo at $z=0$, as quantified in the analysis by McGaugh et al. (2010).

Extending the study of the impact of feedback physics on the evolution of baryonic properties, Oh et al. (2021) evolved the same setup beyond $z=0$. Such a setup in a standard $\Lambda$CDM universe forms a similar total stellar mass by $z=-0.995$ (approximately $96\,{\rm Gyr}$) to that predicted by extrapolating the analytical SFRD fit from Madau & Dickinson (2014). This agreement is also achieved regardless of the resolution of the simulation.

The starting assumption of the current work is that the above modelling correctly represents the operation of the full baryonic physics that applies on scales below those that can be resolved in the simulation, and that this same small-scale modelling can be used in universes beyond the $\Lambda$CDM case that was used for calibration. To the extent that the parameters of the modelling refer to small-scale local quantities, this should be a defensible assumption. Furthermore, we consider only simple counterfactual cosmologies in which the densities of baryons and dark matter are always in the same ratio, so that the internal dynamics of a halo should be followed self-consistently by Enzo in a manner that is independent of the large-scale cosmological context. In any case, it is certainly of interest to see how the modelling functions in universes very different from our own. In the absence of observations, the robustness of such predictions can only be tested by comparison with alternative calculations that use different subgrid approaches. But here we can do no more than explore the predictions of the existing Enzo model, taking care to ensure that the results are at least converged within our resolution limits.

In this section, we test the sensitivity of the Halo Mass Function (HMF) to cosmology, using halo catalogues obtained via the ROCKSTAR halo finder (Behroozi et al., 2013), in which each halo is resolved by a minimum of 20 particles. In Figure 3 we display the evolution of the HMF in our different universes. We choose several snapshots at $t\approx 7,\ 14,\ 26,\ 51\,{\rm Gyr}$, which correspond roughly to $0.5\ t_{\rm H},\ t_{\rm H},\ 2\ t_{\rm H}$ and $4\ t_{\rm H}$ with $t_{\rm H}\approx 13.7\,{\rm Gyr}$ in panels (a) to (d) respectively. Each coloured line refers to the same cosmology consistently across the panels as indicated in the legend and caption.

Since structure formation in the CDM paradigm happens in a bottom-up manner, low mass haloes form first from the primordial density fluctuations before merging to create more massive haloes. Thus the HMF ‘enters’ the plot in Figure 3 from the low mass end. The characteristic cutoff in the HMF then moves towards higher mass as time increases, and is expected to saturate at an asymptotic form once $\Lambda$ dominates. However, the continuing expansion of the universe induces a deterioration of the proper force resolution in our simulations. Once this resolution scale exceeds the virial radius of a halo, the particles in the halo experience too small a gravitational force to remain bound, causing the haloes to dissolve and fade into the background. This eventual loss of haloes is apparent on all mass scales: see e.g. the red line in Figure 3, where the mass function for 10x is seen ‘exiting’ the plot to the left at $t=26$ Gyr, being strongly suppressed with respect to its value at $t=14$ Gyr. It may seem puzzling that the mass function is not eroded first at the lowest masses: these have the smallest virial radii and hence should arguably be the first to experience the effects of insufficient resolution. However, all haloes have a central spike in density where the velocity dispersion is highest, and declining proper resolution will lead to particles being able to escape from this region, so that even the most massive haloes dissolve, in an inside-out manner.

The same general pattern of evolution is seen for different cosmologies, although the unphysical turnaround of the HMF occurs at different times, depending on $\alpha$. With a larger $\alpha$, $\Lambda$ dominates earlier in the universe, leading to an earlier onset of exponential expansion and freezeout of growth. As a consequence, the resolution-induced decline in the number density of haloes across the mass range sets in earlier when $\alpha$ is increased. This difference is obvious when we compare 10x (red), 100x (green) and LCDM (orange) in Figure 3. In panel (a) of that figure, we see that EDS, LCDM, and 10x are all very similar in their mass functions, whereas 100x is displaced to low masses, reflecting the fact that this model is already heavily $\Lambda$-dominated even at $t=7$ Gyr. In turn, 10x is separated from LCDM in panel (b), because this model has already frozen out by $t=14$ Gyr; and by $t=26$ Gyr this model experiences resolution-induced loss of haloes. However, even by the latest time shown (51 Gyr), there is no sign of such problems in LCDM or EDS. For the latter model, the slower expansion means that the force resolution remains adequate at all simulated times, and the lack of freezout means that the most massive halo in the simulation continues to gain mass as structure formation proceeds. For LCDM, we expect that there will also in due course be a loss of haloes as exponential expansion fully sets in, but at $t=51$ Gyr the expansion since $t=14$ Gyr is a factor of 9.9 or 2.4 for LCDM or EDS respectively, so the proper resolutions of these models are not yet vastly different.

We can contrast this behaviour with what is seen when we simply scale the normalization $\sigma_{8}$ within $\Lambda{\rm CDM}$, with otherwise fiducial parameter values. If we neglect issues of numerical resolution, the main factor is a timing issue of when freezeout occurs in each counterfactual universe. For universes with scaled $\Lambda$, e.g., 10x and 100x, freezeout of the HMF into an unevolving asymptotic form will occur at increasingly early times, the more $\Lambda$ is increased. For universes with unchanged $\Lambda$, we can scale $\sigma_{8}$ appropriately to yield this same asymptotic HMF. We would then expect the differences between $\alpha$x, s$\beta$x, and LCDM to be small once both models have frozen out, until resolution effects come to dominate. This is broadly the impression given by Figure 3. For example, in panel (b) of that plot, s10x has a similar HMF to that of 10x, at a time when the latter model has completely frozen out and the former very nearly so. But as discussed earlier, a larger $\alpha$ accelerates the expansion of the universe leading to a deterioration of the proper force resolution. Consequently, by $t=51.00\,{\rm Gyr}$, panel (d) displays HMFs only for s$\beta$x models in addition to LCDM and EDS, and the halo populations of models with increased $\Lambda$ have not been retained to this time.

This inability to follow the halo population into the indefinite future is a potential concern for our principal aim, which is to look at the long-term efficiency of star formation in counterfactual universes. This issue is discussed in the following Section, where we find that in all cases there is an era of peak star formation, followed by a decline in activity that shows good evidence for a convergence in the total production of stars. However, it is important to be convinced that this convergence is not an artefact caused by the disappearance of haloes within the simulation. In order to probe this issue, we chose to repeat the most strongly affected simulation, which is the universe with the largest $\Lambda$ ($100x$). Here, we were able to achieve five additional AMR levels beyond the default of four, which was motivated by our previous results on the evolution to $z=0$ in LCDM. Results from this run are contrasted with the base resolution in Figure 3, where the loss of haloes apparent in $100x$ at 7 Gyr is now removed. We consider below the impact of this variation in resolution on the star-formation history.

We now investigate how the long-term SFRD in these counterfactual models is influenced by cosmology. The maximum time into the future that can be investigated is limited by a number of factors: computational cost; onset of resolution systematics; and on occasion simple crashes of the code in extreme parameter regimes, as discussed in papers 1 & 2. These limits are different for each simulation, as discussed below. But we believe that the star formation in each model shows good evidence of convergence by the point of its maximum reliable time.

Figure 4 shows how the SFRDs evolve in our different simulated counterfactual universes, and a wide diversity of behaviour is apparent. Each cosmology displays a different peak in the SFR, both in terms of the time at which it occurs and the value at the peak. For a universe with higher $\alpha$ scaling applied to $\Lambda$, the peak in the SFRD happens earlier and reaches a lower maximum value. It is interesting to consider the main qualitative drivers of this behaviour. With a higher value of $\alpha$, $\Lambda$ dominates earlier in the universe, so that the HMF freezes out at a lower characteristic halo mass (See Section 3.1). In this exponentially expanding de Sitter phase, we expect that accretion of gas onto haloes will cease, leading ultimately to a reduction in star formation.

But there is also a peak in the SFRD for EDS, in which there is no freezeout, consistent with the findings of Salcido et al. (2018). This behaviour is a reflection of gravitational heating of the IGM by continuous infall of matter and orbiting sub-structures, which convert gravitational potential energy into heat (Khochfar & Ostriker, 2008; Dekel & Birnboim, 2008). Thus baryonic physics must also be playing an important role in the SFRD peak. At early times, there is general agreement that star formation is driven by the infall of new cold material from the cosmic web, rather than coming from the pre-existing baryons in haloes, which are driven into a hot gaseous halo (Dekel et al., 2009a, b; Somerville & Davé, 2015). The question is then what happens to this material in the very long term, when it will eventually be able to cool and potentially generate further star formation. This is in effect the issue that is addressed directly by our simulations.

For counterfactual universes with reduced $\sigma_{8}$, structure formation and its associated star formation is suppressed with respect to LCDM in all cases. However, the detailed effect is very different from when we scale $\Lambda$. As demonstrated by Figure 4, the peak in the SFRD for e.g. s10x occurs much later than for the matched 10x simulation, and the peak is much lower, so that the total amount of star formation in 10x is higher than in s10x, even though the activity in the latter case occurs over a larger range of time.

For all these star formation histories, it is convenient to have a smooth fit to the time-dependent SFRD:

where $a=0.015,b=2.7,c=2.9,d=5.6$ for the analytic fit to observations (Madau & Dickinson, 2014). It turns out that the SFRD from our simulations can also be fitted accurately with this double power-law in $(1+z)$. We illustrate the goodness of the fit both directly and via its residuals in Figure 5. However, we do see cases where the computed SFRD rises sharply above the declining fit at late times. This upturn was discussed in detail in Oh et al. (2021), where we established that it marks the maximum time at which the SFRD calculation can be treated as numerically reliable on resolution grounds. The model fits to the SFRD data therefore exclude the data beyond these times.

We also note a drop in the EdS SFRD above about 100 Gyr. This drop coincides with the halo mass function ceasing to evolve in our simulation. This is in contrast with expectations for the EdS case, which predicts continual growth of haloes. This lack of evolution reflects our failure to resolve haloes at late times, which eventually limits all our simulations. In the EdS case, this appears to bias the star formation to low values; but by this point the total production of stars has converged.

The values of the fitted parameters of Equation 9 are summarised in Table 2. On top of these parameters, we also include a measure of the precision of the fit, $R^{2}$:

where $\overline{\rm SFRD}_{\rm sim}$ denotes the mean of the SFRD from the simulation and the rest of the symbols have the same meanings as before. This quantity provides a quantitative complement to the residuals shown in Figure 5. The value of $R^{2}$ varies between zero and one with higher values corresponding to a better fit. As we see in Table 2, the SFRD in counterfactual universes with scaled $\Lambda$ can be fitted in this way to high precision, about as well as our observed universe.

The main potential concern with these results is the impact of the inability of the simulations to follow the halo population into the indefinite future, as discussed above in Section 3. Provided the critical time at which haloes are lost is beyond the peak in the SFRD, we can have some confidence that the peak in star formation is physical and that the apparent convergence in the total production of stars is robust. The model that comes closest to violating this criterion is $100x$, where Figure 3 shows that haloes are being lost by $t=7$ Gyr, whereas the peak in the SFRD is at about 2 Gyr. This particular simulation is the one that has the strongest impact on our understanding of the suppression of asymptotic star-formation efficiency at high $\Lambda$, so we felt it was important to repeat this calculation with a higher resolution than our fiducial choice. As discussed in Section 3, this remedies the loss of haloes at 7 Gyr, and so we can consider the associated impact on the SFRD. The runs with two different resolutions are contrasted in Figure 4, where it can be seen that the results are very close, both in terms of the peak in SFRD and the subsequent decline. In summary, the haloes are dissolving at a time when the star forming activity is already low. We therefore conclude that our results are a fair representation of the predictions of the Enzo model, and are not affected by resolution artefacts.

A significant element in the Enzo model for star formation is allowance for the effect of a UV background on the thermal evolution of the IGM. This presents a problem of self-consistency, since the level of the UV background depends on the star-formation history that we are aiming to calculate. We can only approach this iteratively: assume a background, compute the SRFD, and then estimate a revised background from this. We make the simple assumption that the UV background is dominated by young massive stars, whose short lifetime means that their abundance is directly proportional to the SFRD. Therefore, we scale the UV background with the ratio of the SFRDs in each counterfactual universe to LCDM at identical times:

where $\Gamma$ is the photoheating rate from the UV background, subscripts ‘scaled’ and ‘LCDM’ represent counterfactual and LCDM cosmology respectively. In practice, we implement the UV scaling via the smooth fits to the SFRD discussed above.

We thus ignore any contribution from quasars, which are the most probable alternative source of the UV background. However, the emissivity of quasars appears to be sub-dominant in comparison to galaxies unless there is a significant number of very low luminosity AGN (Meiksin, 2005; Bolton & Haehnelt, 2007; Srbinovsky & Wyithe, 2007). We assume that this smaller relative contribution applies in all our models. We also neglect more exotic alternative contributions to the UV background, such as an early generation of black holes (Ricotti & Ostriker, 2004; Ricotti et al., 2005; Venkatesan et al., 2001). While possible in principle, there is currently little observational support for significant contributions from these mechanisms (Meiksin, 2009).

We can now investigate the self-consistency of the SFRDs in counterfactual universes by implementing the scaled UV background. The initial computations assumed the standard UV background from the LCDM run, and these were then repeated with the UV background based on this initial SFRD calculation, adopting the scaling from equation 11. If we refer to Figure 5a, the SFRDs of counterfactual universes relative to LCDM suggest that the impact of the scaled UV background on the SFRDs is small. At early times, the SFRDs of all counterfactual universes are similar to LCDM or are lower, meaning that deviations in the photoheating rates in these universes will not be as extreme as switching the UV background off completely. At late times, even though the SFRD of the EDS cosmology is much higher than LCDM, the photoheating rates at these times are already low. Therefore, any effect from a scaled UV background will be less extreme than removing the UV background entirely. In practice, therefore, for computations of the SFRD it seems sufficient to use the Haardt & Madau (2012) UV background model, even though in principle this should vary with cosmology.

Even though the history of star formation turns out not to have a strong dependence on the assumed UV background, this hides some of the underlying complexity of the situation, and we show in this Section that the properties of the IGM do certainly depend on the assumed level of the UV radiation field. Following Oh et al. (2021), we focus on the phase distribution of the IGM material, defined as having an overdensity less than $10^{3}$ (Davé et al., 2001). We adopt the power-law model of Hui & Gnedin (1997), in which the IGM temperature is given by

where $T_{0}$ is the temperature at the mean density, $\delta$ is the gas overdensity and $\gamma$ is a sensitivity parameter for the equation of state. We show the differences in these IGM parameters between simulations with the default Haardt & Madau (2012) UV background (HM) and the scaled UV background (UVB) in Figure 6, assuming that the SFRD remains exactly the same in each universe regardless of the UV background. We break down the analysis into universes with a scaling of $\Lambda$ (top) and $\sigma_{8}$ (bottom).

In Figure 6a, we see that the scaled UV background does not affect the evolution of $T_{0}$. The solid lines (HM: standard Haardt–Madau background) and dashed-dotted lines (UVB: scaled background) are indistinguishable from each other. But the equation of state is affected: at intermediate times, the HM $\gamma$ is consistently lower than the UVB $\gamma$, indicating that the temperature of the IGM is more sensitive to its overdensity in the simulations with a scaled UV background. The reduced heating from the scaled UV background makes cooling the dominant process, which is sensitive to density.

On the other hand, the IGM in the EDS simulation shows an almost identical evolution regardless of the UV background. The only difference occurs in the value of $T_{0}$ before $20\,{\rm Gyr}$, when the SFRD reaches a peak. Before this time, the photoheating rate in the UVB run is lower than the HM case because the SFRD of EDS is lower than that of LCDM. Therefore, the resulting value of $T_{0}$ in the UVB case is slightly smaller than for HM, although this slight difference is not sufficient to affect $\gamma$. Beyond this time, the HM photoheating rates are low, reducing the impact of any scaling of the SFRDs, so that the subsequent evolution of $T_{0}$ and $\gamma$ is independent of the UV background.

In contrast to these results, $T_{0}$ for the scaled $\sigma_{8}$ simulations varies significantly according to changes in the UV background, as shown in in Figure 6b. $T_{0}$ values for simulations with a scaled UV background are consistently lower than for their counterparts with an HM UV background. If we compare the SFRDs of s5x, s10x, s20x and s100x with LCDM, they are significantly lower, translating into a weaker UV background. This reduced heating rate also causes $T_{0}$ to reach the simulation’s assumed temperature floor of $1\,{\rm K}$, which correspondingly causes $\gamma$ to approach a value of unity more rapidly. The opposite is true for EDS, for which the SFRD is higher than for LCDM. This difference leads to the increase in $\gamma$ at late times seen in Figure 6a. Simulations with a scaled UV background thus experience a slightly different evolution from the HM calculation. We will therefore henceforth consider only the more self-consistent analyses of the IGM that use simulations with a scaled UV background.

With this restriction to a single prescription for the UV background, Figure 7 superimposes the evolution of the IGM associated with all the various counterfactual universes into a single panel. At the start of the simulations, $T_{0}$ in 10x and 100x is higher than in LCDM for a short period of time. If we compare the SFRDs of these simulations in Figure 5a, we find that the initial SFRDs in 10x and 100x are slightly higher than in LCDM. This difference raises the initial heating rates, resulting in a higher $T_{0}$. The IGM then cools rapidly to the temperature floor, because of a combination of a much lower SFRD and a faster expansion as compared to LCDM. As a result the equation of state also rapidly flattens to $\gamma=1$. The evolution of the IGM that takes place over a period of approximately $100\ {\rm Gyr}$ in LCDM is essentially compressed into $40\ {\rm Gyr}$ for 10x and even further into $20\ {\rm Gyr}$ for 100x. We also observe a similar trend for simulations with scaled $\sigma_{8}$. A lower $\sigma_{8}$ leads to less star formation, resulting in a lower heating rate from the UV background. Therefore, simulations with larger downward scaling of $\sigma_{8}$ experience a faster decline of $T_{0}$ towards zero and a similarly rapid flattening of $\gamma$.

In summary, parameter scaling determines the speed at which cosmic evolution causes the characteristic IGM temperature $T_{0}$ to approach the temperature floor and the equation of state slope $\gamma$ to flatten, but the behaviour is different for universes with scaled $\Lambda$ and $\sigma_{8}$. For the former, the histories of the IGM in all models follow a common track at high $z$, departing from this at the point where $\Lambda$ dominates. In the simulations where $\sigma_{8}$ is lowered, the initial state is already different, with LCDM always having a higher $T_{0}$ than the s$\beta$x simulations. The reduced clustering also leads to both $T_{0}$ dropping to the temperature floor more rapidly, accompanied by a faster flattening of $\gamma$.

Given these substantial changes to the IGM, it is stiking that there is little to no directly resulting effect on the SFRD. At late times when the universe is dominated by $\Lambda$ and the UV background is insignificant, the gas in the IGM cools to the temperature floor regardless. At early times, we find that the absence of the UV background affects gas of low overdensity most significantly, and this material lies below the threshold for star formation.

An interesting aspect of the long-term star formation history is its associated chemical evolution. Will stars continue to form from ever more polluted gas in their vicinity, or will there be time for pristine gas from the IGM to mix? We investigate this by looking at how the average metallicity of young stars ($<500\ {\rm Myr}$) evolves as a function of time in our counterfactual universes. As described in Section 2 of Oh et al. (2020), Enzo tracks metals as part of its feedback prescription, following the injection of metals into the IGM from star-forming activity, and assigning an appropriate metallicity to each newly-formed star particle. We plot the resulting average metallicity as a ratio to the solar metallicity against time in Figure 8. Except for EDS and s100x, all other simulations exhibit a similar evolution in the average metallicity of their young stars: this increases up to a peak, declines, and finally starts to increase again at very late times. This final increment in metallicity might plausibly be attributed to the isolated evolution of haloes as a result of freezeout: the haloes could then function as closed boxes, so that the IGM becomes ever more polluted, raising the metallicity of stars formed at late times. However, we note that this feature seems to occur at the same point as the spurious star formation associated with worsening resolution in these models, and we therefore assume that this feature in the metallicity history is similarly not physical.

The general behaviour is then that the stellar metallicity scales with SFRD in the simulation. With the peak metallicity of stars coinciding with the peak of SFRD, the decline in SFRD is also associated with a drop in metallicity. When the simulation is not actively forming stars, there is time for the fuel store of potential future star formation to mix with the inflow of fresh pristine gas from the IGM. Because of the difference in the gas cooling and star formation time scale discussed earlier, the metallicity of the young stars decreases. But in any case, the peak metallicity is lower in the highly scaled models where the total star-forming activity is suppressed, reflecting the lack of metal pollution of the IGM. Conversely, with the active star formation seen in EDS, the gas is continuously injected with new metals associated with feedback. Therefore, after the initial rise in metallicity, stars in EDS remain highly enriched due to the consistently higher SFR shown in Figure 4.