Emergence of the temperature–density relation in the low-density intergalactic medium Open Access

ABSTRACT

We examine the evolution of the phase diagram of the low-density intergalactic medium during the Epoch of Reionization in simulation boxes with varying reionization histories from the Cosmic Reionization on Computers project. The probability density function (PDF) of gas temperature at fixed density exhibits two clear modes: a warm and a cold temperature mode, corresponding to the gas inside and outside of ionized bubbles. We find that the transition between the two modes is ‘universal’ in the sense that its timing is accurately parametrized by the value of the volume-weighted neutral fraction for any reionization history. This ‘universality’ is more complex than just a reflection of the fact that ionized gas is warm and neutral gas is cold: it holds for the transition at a fixed value of gas density, and gas at different densities transitions from the cold to the warm mode at different values of the neutral fraction, reflecting a non-trivial relationship between the ionization history and the evolving gas density PDF. Furthermore, the ‘emergence’ of the tight temperature–density relation in the warm mode is also approximately ‘universally’ controlled by the volume-weighted neutral fraction for any reionization history. In particular, the ‘emergence’ of the temperature–density relation (as quantified by the rapid decrease in its width) occurs when the neutral fraction is 10⁻⁴ ≲ X_H i ≲ 10⁻³ for any reionization history. Our results indicate that the neutral fraction is a primary quantity controlling the various properties of the temperature–density relation, regardless of reionization history.

1 INTRODUCTION

A second after the big bang, the Universe was filled with a plasma of free protons and electrons. As the Universe expanded over time, it also cooled enough that these protons and electrons combined to form neutral hydrogen atoms, decoupling from the primordial blackbody radiation known today as the cosmic microwave background radiation (Penzias & Wilson 1965). Immediately after this epoch, there were no other sources of light as large-scale structures had not yet been formed. This waiting period, often denoted the ‘Dark Ages’, lasted until the first stars began to form and emit light some 100–200 Myr after the big bang.

This emitted ultraviolet (UV) radiation was then able to ionize nearby regions, or ‘bubbles’, of hydrogen gas in what is known as the Epoch of Reionization (for a recent review, see Wise 2019). This process occurred approximately from z = 15 to z = 6. These bubbles continued to grow as more neutral hydrogen became reionized and the bubbles eventually began to overlap, until nearly the entire intergalactic medium (IGM) was ionized.

It is known that the IGM at z < 6 consists mostly of ionized hydrogen because of the observed absorption spectra of distant quasars. As (UV) photons from a quasar pass through a cloud of neutral hydrogen gas, an absorption line corresponding to the Lyman α transition appears. From passing through multiple gas clouds at different redshifts, a series of absorption lines results called the Lyman α forest. The fact that this absorption is not complete and that flux from these high-redshift objects is transmitted is proof that the IGM hydrogen must be highly ionized, as fully neutral (or even neutral at the ∼10⁻³ level) IGM would absorb all Lyman α photons (Gunn & Peterson 1965).

Simulations of the IGM show that the mean temperature of the IGM is approximately related to its density by a power-law relation, T(δ) = T₀(1 + δ)^{γ − 1}, where T₀ is the temperature at the mean density and γ is the power-law slope (Hui & Gnedin 1997; Schaye et al. 1999; Ricotti, Gnedin & Shull 2000; McQuinn & Upton Sanderbeck 2016).

Measurements of the Lyman α forest can thus be used to constrain the parameters of this temperature–density relation (TDR), T₀ and γ (e.g. Boera et al. 2014, 2016; Nasir, Bolton & Becker 2016; Rorai et al. 2018; Walther et al. 2019; Müller, Behrens & Marsh 2021). Both T₀ and γ are theoretically expected (e.g. Hui & Gnedin 1997; Ricotti et al. 2000; Puchwein et al. 2012, 2023; McQuinn & Upton Sanderbeck 2016; Nasir et al. 2016) and have been observationally measured to evolve (e.g. Chang, Broderick & Pfrommer 2012; Boera et al. 2014; Hiss et al. 2018; Walther et al. 2019).

While a number of studies have found γ > 1 between 1 ≲ z ≲ 5 (after reionization) so that temperature increases with increasing density (Boera et al. 2014; Rorai et al. 2018; Walther et al. 2019; Müller et al. 2021), there has also been observational evidence for an ‘inverted’ TDR with γ < 1 at low density and low redshifts (e.g. Bolton et al. 2008; Chang et al. 2012; Puchwein et al. 2012), or a broken power law where the TDR is inverted for densities lower than about twice the mean density (Rorai et al. 2017). Such an inverted TDR can only be produced by non-radiative heating sources, for example, ultra-high-energy cosmic rays from blazars driving plasma instabilities in the IGM (e.g. Chang et al. 2012; Puchwein et al. 2012).

Hence, further theoretical understanding of how the TDR of the low-density IGM emerges during the Epoch of Reionization is needed to better interpret measurements of the thermal state of the IGM post-reionization. Recently, simulations that track the post-reionization thermal state of the IGM have been performed, which both do (Puchwein et al. 2023) and do not (Villasenor et al. 2021) explicitly model the reionization process, but the actual process of ‘emergence’ of the TDR during reionization is not yet well understood.

In this paper, we use existing Epoch of Reionization simulations to study the timing of the thermal evolution of the low-density IGM, including when the TDR emerges. First, we introduce the simulation project used in this study and show how the phase diagrams demonstrate the IGM evolution. We then describe the separation of two thermal modes in this evolution, analyse the timing of the transition between these modes, and look at how the evolution of the neutral fraction is related to this timing. We examine the TDR that emerges post-reionization and see how its slope and full width at half-maximum (FWHM) change with time and across reionization histories.

2 METHODOLOGY

2.1 CROC simulations

In this study, we use numerical simulations from the ‘Cosmic Reionization on Computers’ (CROC) project (to learn more about the CROC project, see Gnedin 2014; Gnedin & Kaurov 2014). These simulations model self-consistently all physical processes generally considered to be relevant to cosmic reionization, such as radiative transfer, star formation and feedback, and ionizing radiation from stars and (approximately) from quasars. We use four 40 h⁻¹ comoving Mpc boxes that sample the full range of reionization histories modelled by the CROC project. Fig. 1 shows the reionization histories of the different simulation boxes. These four models differ only by the value of the rms density fluctuation at the box scale (the so-called DC mode) but otherwise include the same physical modelling and are simulated with the same numerical parameters.

2.2 IGM evolution

Before the Epoch of Reionization, the IGM consists of mostly neutral hydrogen (as well as neutral helium) after the recombination of protons and electrons in the Universe. However, as gas clouds begin to condense and form stars and galaxies, these objects emit radiation that reionizes the IGM over time. The excess energy of ionizing photons is deposited as heat in the IGM, raising the temperature to a few tens of thousands of kelvins, a typical temperature of the photoionized gas. This thermal evolution of the IGM is shown through the temperature–density phase diagram in Fig. 2, represented at three different scale factors.

The quantity on the x-axis of Fig. 2 is

where δ represents the overdensity, calculated as the local density ρ minus the mean density of the universe |$\overline{\rho }$|⁠, divided by the mean density of the universe, as shown by equation (1).

The diagram at three scale factors in Fig. 2, which has an x-axis density bin width of log(1 + δ) = 0.01 and a y-axis temperature bin width of log(T/K) = 0.01, showcases how the gas is being heated from the ionizing radiation and is moving to a higher temperature, settling between T = 10⁴ K and T = 2 × 10⁴ K in the last panel. This warm gas is expected to be within the bubbles of ionized IGM, while the cold, neutral gas is expected to be outside of those ionized bubbles.

While the underlying physical processes of the TDR are well understood, the resulting quantitative trends of the evolution are not. This is particularly because the Lyman α forest observations used are difficult to obtain at higher redshifts corresponding to the Epoch of Reionization given the higher neutral fractions, and therefore these observations do not constrain the TDR during the Epoch of Reionization. We can investigate this evolution by defining the two gas ‘modes’ (warm and cold) and finding relations between the timing of the transition from the cold to the warm mode and the timing of the overall reionization process.

2.3 Mode separation methods

The two modes of the IGM can easily be seen in a one-dimensional histogram for gas temperature at a fixed density. An example of such a histogram is shown in Fig. 3. The leftmost peak of each panel of the figure represents the colder mode of neutral hydrogen, whereas the rightmost peak of each panel represents the warmer mode of reionized hydrogen. The modes are almost always distinct but are not completely separated. Because of this, we define two different methods for separating all or most of the cells in the simulation into the two modes.

The left panel of Fig. 3 displays the mode separation by use of the cut-off method, where the modes are first separated by the |${\sim }10^{4} \, \mathrm{K}$| threshold for ionized hydrogen. For a peak above |${\sim }10^{4} \, \mathrm{K}$|⁠, we label the resulting mode the ‘warm mode’, and for a peak below |${\sim }10^{4} \, \mathrm{K}$|⁠, we label it the ‘cold mode’. Then, we only count the bins that lie above a threshold (such as 10 per cent of the peak height, as shown in the left panel in Fig. 3). This is done to establish a more constrained mode definition that focuses on the peak values and a more limited range of their distributions. The right panel of Fig. 3 displays the mode separation by use of the minimum method, where the modes are simply separated by the lowest histogram value between the two peaks. This method allows for a wider distribution in the mode definition, and it does not require the explicit |${\sim }10^{4} \, \mathrm{K}$| threshold. However, the temperature at which this minimum occurs is always near this threshold. Further analysis of these modes throughout the text suggests that the two mode separation methods yield indistinguishable results. See Appendix  A to view the insignificant differences when comparing the two mode separation methods.

3 RESULTS

3.1 Timing of the transition

By assigning all or most of the cells in the simulation to one of the modes, we can calculate the cold mode fraction with respect to the sum of both the warm and cold modes. Fig. 4 displays this cold mode fraction as well as the volume-weighted neutral fraction (the same as in Fig. 1) over time for each box. The cold mode fraction has a similar evolution as the neutral fraction, but it transitions slightly earlier than the neutral fraction. This is because higher energy photons of the ionizing radiation have a longer mean free path (MFP) and will heat the neutral gas before the more numerous lower energy ionizing photons can fully ionize it. Because the reionization process is dominated by these low-energy photons and heating is more sensitive to higher energy photons, any particular place in the IGM is heated somewhat before being fully ionized.

From Fig. 4 we can see that the cold mode fraction equals 0.50 at some scale factor, which is the point where the cold and warm mode fractions are equal. This signifies the moment when the warm mode begins overtaking the cold mode in volume. We call the scale factor at which this occurs the ‘critical scale factor’, and this value can be used as a quantitative measure of the timing of the cold-to-warm mode transition in each density bin. Fig. 5 displays the critical scale factor across the density bins for each simulation box.

The heating and the ionization of the IGM are two aspects of the same process, and hence one may wonder whether the connection between them is ‘universal’, i.e. independent of the details of the actual reionization history. In order to explore that, we show in Fig. 6 a variation of Fig. 5 with values for critical scale factors replaced by the values of the volume-weighted neutral fraction at that moment for that simulation box.

Different simulation boxes (with their differing ionization histories) at each density have near-universal neutral fractions at the critical scale factors for the cold modes. In other words, when the thermal history is parametrized by the value of volume-weighted neutral fraction, the transition from cold to warm mode becomes almost universal. Note that this conclusion also holds if we separate the cold and warm modes via the cut-off method (see Appendix  A). We elaborate on this universality in Section 4.

3.2 Emergence of a temperature–density relation

The tight TDR is perhaps the most recognizable feature of the Lyman α forest at intermediate redshifts. Fig. 2 shows that such a relation is not present during the Epoch of Reionization, and hence it must ‘emerge’ shortly after reionization is completed. In order to determine the emerging TDR T(δ) = T₀(1 + δ)^{γ − 1}, we extract the most likely temperature at each density in the warm mode – i.e. the peak of the temperature histogram in each density bin, shown with vertical lines within the rightmost (warm) modes in Fig. 3.

The parameters of the TDR, the temperature at the mean density T₀ and the power-law slope γ, can be measured by fitting a power law to a range of densities around the mean cosmic density. This procedure is not well defined, though, as it depends on the specific density range chosen. Instead, we noticed that the full TDR is well fitted by a third-degree polynomial in log–log space, which we can then use to determine T₀ and the slope γ (at the mean density) accurately. The ‘emergence’ of the TDR is not quantified by its parameters, however, but rather by how tight it is. In order to quantify that, we also measure the  FWHM of the warm mode around the peak temperature. Fig. 7 shows the most likely temperature of the warm mode as well as its FWHM at three different scale factors.

From Fig. 7, we can see that the FWHM of the warm mode in the mean density bin 0 < log(1 + δ) < 0.01 is decreasing as the scale factor increases, illustrating the emergence of the TDR over time. We see for each scale factor shown that the FWHM is greater at lower densities than near the mean density, confirming that as the low-density gas is transitioning from cold to warm later than gas near the mean density (see Fig. 5), the TDR also emerges later in the low-density gas.

The behaviour of the slope and the FWHM at the mean density (which we use to characterize the ‘emergence’ of the TDR) can also be investigated over time or with respect to the neutral fraction, shown in Figs 8 and 9, respectively. From the differences in the simulation boxes, we can see that in Fig. 8, earlier reionization leads to a consistently steeper slope and smaller FWHM than later reionization does. From the evolution of the slope and FWHM across neutral fractions shown in Fig. 9, we again observe an approximate universality among the different simulation boxes, i.e. among the different reionization histories. From this figure, we can see that this universality is not exact and is not as precise as the universality for the timing of the cold-to-warm mode transition (Fig. 6). Additionally, the uncertainty in these lines has not been quantified, so we cannot say whether the differences between the reionization histories are larger than the inherent uncertainty for the individual histories (due to the inhomogeneous nature within a single simulation box). Multiple factors can affect this slight non-universality, from actual physical non-universality, to numerical artefacts in the simulations, to particular numerical choices in our analysis (such as binning). At this point, it does not seem worthwhile to explore these factors – the fact that the (near-) universality exists is much more important than the reasons for it not being perfect. In particular, Fig. 9 shows that the emergence of the TDR primarily occurs when the neutral fraction is 10⁻⁴ ≲ X_H i ≲ 10⁻³ irrespective of other details of reionization.

4 DISCUSSION AND CONCLUSIONS

In this paper, we investigate the evolution of the thermal state of the IGM during cosmic reionization. Our primary goal is to better understand how the reionization history impacts the emergence and the evolution of the tight TDR, so ubiquitous in the lower redshift Lyman α forest.

Our main finding is that the emergence of the TDR is nearly ‘universal’, in the sense that when expressed as a function of the volume-weighted neutral fraction instead of the actual time itself, both timing and the shape of the TDR become only weakly dependent on the details of ionization history. In particular, the timing of the transition from the cold neutral IGM into the warm ionized IGM (Fig. 6) is almost perfectly universal at all densities. This fact is not necessarily surprising, since ionized gas is warm, but it is not trivial either. First, not all warm gas is ionized as Fig. 4 shows (the heating history precedes the ionization history). Secondly, and most importantly, Fig. 6 relates the cold-to-warm transition at a given density to the global volume-weighted neutral fraction, and hence depends on the overall density distribution. Hence, the (near-) universality at a fixed density is not necessarily assured. The mean neutral fraction depends on the full ranges of densities. Hence, the cold-to-warm transition at a fixed density does not have to be correlated with the mean neutral fraction. For example, if the density PDF evolves faster at higher densities than at lower densities, the mean neutral fraction would evolve faster than the cold-to-warm transition at lower densities. The fact that this is not the case and the universality is highly accurate implies that the timing of reionization is in some way ‘synchronized’ with the evolution of the density PDF – an idea that has not been clearly stated or explored previously. More than that, we observe the universality between the four distinct 40 h⁻¹ Mpc boxes, implying that the ‘synchronization’ between reionization and the evolution of the density PDF occurs on scales at least as small as several tens of cMpc.

In order to gain some intuition on how such a ‘synchronization’ may develop, let us consider a toy model of a non-evolving density field. Let us also imagine, merely for the sake of the argument, that reionization proceeds totally randomly. Then, half of the gas at a given density would become warm when half of the universe is ionized; i.e. the timing would be universal, but the dependence on the mean neutral fraction would be trivial, and all lines in Fig. 6 would be strictly horizontal. Notice that accounting for the finite width of ionization fronts would not solve the difference, but rather exacerbate it: since the MFP is longer in lower density gas, the lower density gas would get heated at slightly higher values of the mean neutral fractions than the denser gas; i.e. in the random ionization scenario, lines in Fig. 6 would slope slightly downward.

A much more realistic model for reionization is a barrier-crossing formalism of Furlanetto, Zaldarriaga & Hernquist (2004). In that model, the lower density gas gets ionized later than the higher density gas. Hence, there is indeed a positive correlation between the neutral fraction at the time of the cold-to-warm transition at fixed density and the value of density. However, since different 40 h⁻¹ Mpc boxes have different density PDFs, the dependence shown in Fig. 6 would not be universal among our four simulation boxes. In this paper, we do not offer an alternative scenario that explains Fig. 6 – understanding the universality will require a significantly larger effort, and a focus significantly different from the subject of this paper, and hence we leave such effort for future work.

The near-universality in the evolution of the TDR of the warm mode (the only mode where the well-defined TDR actually exists) is less striking than the near-universality of the timing of the cold-to-warm transition, but it is still significant (Fig. 9). This allows us to make generic conclusions about the emergence of the tight TDR in reionization history-independent terms. For example, the emergence primarily occurs when the volume-weighted neutral fraction decreases from 10⁻³ to 10⁻⁴.

The (near-) universality is a prediction that is easily testable by observations: since different regions of the universe reionize at different times, observations of absorption spectra of different quasars should show this near-universality – or demonstrate a major shortcoming of the CROC project simulations.

ACKNOWLEDGEMENTS

AW would like to thank the National Science Foundation for funding the University of Michigan Research Experiences for Undergraduates programme (grant no. #2149884). This work was supported in part by the National Aeronautics and Space Administration (NASA) Theoretical and Computational Astrophysics Network grant 80NSSC21K0271. This paper has been co-authored by Fermi Research Alliance, LLC under contract no. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. This work used resources of the Argonne Leadership Computing Facility, which is a Department of Energy Office of Science User Facility supported under contract no. DE-AC02-06CH11357. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment programme. This research is also part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana–Champaign and its National Center for Supercomputing Applications. This work was completed in part with resources provided by the University of Chicago’s Research Computing Center. DR acknowledges funding provided by the Leinweber Graduate Fellowship at the University of Michigan. We would lastly like to thank the anonymous referee for constructive and helpful comments.

DATA AVAILABILITY

The phase diagram and reionization history data for the four simulation boxes used in this paper can be made available upon request.

The code used to display all plots in this work can be found within the GitHub repository: https://github.com/wellsalexandra/IGM-Phase.

References

APPENDIX A: INDISTINGUISHABLE RESULTS FROM MODE SEPARATION METHODS

Fig. A1 is a version of Fig. 6 that shows both mode separation methods. Because the two methods do not yield notable differences in their results, we resolved to only show this plot for one method in Fig. 6.