Dark-Ages Reionisation & Galaxy Formation Simulation XVI: The Thermal Memory of Reionisation

TL;DR

This study uses thermal statistics of the intergalactic medium to constrain the epoch of reionisation, revealing how temperature evolution retains memory of reionisation timing and spectral properties.

Contribution

It introduces a method to constrain reionisation timing and spectral slope using thermal history and temperature-density relations from galaxy formation models.

Findings

Reionisation likely completed at redshift ~6.8.

Thermal memory diminishes as gas cools, forming a power-law temperature-density relation.

Future observations can break degeneracies between reionisation timing and spectral properties.

Abstract

Intergalactic medium temperature is a powerful probe of the epoch of reionisation, as information is retained long after reionisation itself. However, mean temperatures are highly degenerate with the timing of reionisation, with the amount heat injected during the epoch, and with the subsequent cooling rates. We post-process a suite of semi-analytic galaxy formation models to characterise how different thermal statistics of the intergalactic medium can be used to constrain reionisation. Temperature is highly correlated with redshift of reionisation for a period of time after the gas is heated. However as the gas cools, thermal memory of reionisation is lost, and a power-law temperature-density relation is formed, $T = T_0(1+\delta)^{1-\gamma}$ with $\gamma \approx 1.5$. Constraining our model against observations of electron optical depth and temperature at mean density, we find that reionisation likely finished at $z_{\rm{reion}} = 6.8 ^{+ 0.5} _{-0.8}$ with a soft spectral slope of $\alpha = 2.8 ^{+ 1.2} _{-1.0}$. By restricting spectral slope to the range $[0.5,2.5]$ motivated by population II synthesis models, reionisation timing is further constrained to $z_{\rm{reion}} = 6.9 ^{+ 0.4} _{-0.5}$. We find that, in the future, the degeneracies between reionisation timing and background spectrum can be broken using the scatter in temperatures and integrated thermal history.

Equations

$$\bar{\epsilon }=\frac{f_{esc}{N}_{\gamma }}{\frac{4}{3}\pi r^3(1-0.75Y_{\mathrm{He}})m_{\mathrm{p}}}{\dot{m}}_{\ast }(r),$$

$${\bar{J}}_{21}=\frac{(1+z{)}^2}{4\pi }{\lambda }_{\mathrm{mfp}}{h}_{\mathrm{p}}\alpha \bar{\epsilon },$$

$$J(\nu )={\bar{J}}_{21}{\left(\frac{\nu }{{\nu }_{\mathrm{HI}}}\right)}^{-\alpha },$$

$${\Gamma }_i={\int }_{{\nu }_i}^{\infty }\frac{4\pi J_{\nu }}{h_{\mathrm{p}}\nu }{\sigma }_i(\nu )d\nu ,$$

$$g_i={\int }_{{\nu }_i}^{\infty }\frac{4\pi J_{\nu }}{\nu }(\nu -{\nu }_i){\sigma }_i(\nu )d\nu ,$$

$$\frac{d{\tilde{X}}_i}{dt}=-{\Gamma }_i{\tilde{X}}_i+\sum _{j,k}{\tilde{X}}_j{\tilde{X}}_k{R}_{jk}\frac{{\bar{\rho }}_b(1+\delta )}{m_{\mathrm{p}}}$$

$$\begin{array}{rl}\frac{dT}{dt}=& \frac{2}{3k_B{\sum }_i{\tilde{X}}_i}[G(t)-\Lambda (t,{\tilde{X}}_i)]-2HT\\ +& \frac{2T}{3(1+\delta )}\frac{d\delta }{dt}-\frac{T}{{\sum }_i{\tilde{X}}_i}\frac{d{\sum }_i{\tilde{X}}_i}{dt},\end{array}$$

$$f_{esc}=\text{min}\left[f_5{\left(\frac{1+z}{6}\right)}^{\beta },1.0\right]$$

$$T_{\text{reion}}\approx \frac{1}{3k_b}⟨E_{\gamma }⟩.$$

$$E_{\text{excess}}=\frac{{\int }_{{\nu }_i}^{\infty }\frac{4\pi J_{\nu }}{h_{\mathrm{p}}\nu }{f}_{i}d\nu }{{\int }_{{\nu }_i}^{\infty }\frac{4\pi J_{\nu }}{h_{\mathrm{p}}\nu }{h}_{\mathrm{p}}(\nu -{\nu }_i)f_{i}d\nu }$$

$$f_i=\frac{n_i{\sigma }_i}{\sum n_j{\sigma }_j}$$

Full text

††pubyear: 2019††pagerange: Dark-Ages Reionisation & Galaxy Formation Simulation XVI: The Thermal Memory of Reionisation–B

1 Introduction

Over 13 billion years ago, when the universe was approximately 400,000 years old, it consisted mostly of a neutral atomic gas of hydrogen and helium. Over time the gas began to cool and collapse into the first stars and galaxies. Some of the radiation from these sources was energetic enough to strip electrons from the surrounding atomic gas, ionising it. This period of time, from the birth of the first stars until almost all of the atomic gas in the universe had been reionised, is referred to as the Epoch of Reionisation (EoR). The EoR is the last large-scale cosmic event to be studied in detail, and is of great interest to cosmology as it contains information about the formation processes behind the first galaxies in our universe via their effect on the intergalactic medium (IGM). There are many unanswered questions concerning the EoR, including details of its structure, duration, and effect on subsequent galaxies.

However, observations of the EoR at optical to near infrared wavelengths are made difficult by the absorption of Lyman alpha photons by neutral hydrogen, which is optically thick to this wavelength even at low concentrations. Since all wavelengths below the Ly$\alpha$ line will eventually redshift to that wavelength, measuring the presence of neutral hydrogen using Ly$\alpha$ optical depth probes the tail end of reionisation (Greig & Mesinger, 2017), since the absorption saturates around neutral fractions of $x_{HI}\lesssim 10^{-4}$ (Fan et al., 2006). While saturation varies between sightlines, the strength of this absorption line prohibits direct observations deeper into the EoR. In addition, we can infer the number of free electrons, and hence ionised gas, along a sightline by measuring the Thomson scattering of CMB photons. However this is an integrated measure that cannot distinguish between reionisation histories of different durations.

The reionisation of the IGM is accompanied by a large increase in temperature, to $\sim 2\times 10^{4}$ K, followed by cooling on a cosmological timescale (Miralda-Escudé & Rees, 1994; Hui & Gnedin, 1997; Schaye et al., 2000; Hui & Haiman, 2003; Furlanetto & Oh, 2009; Upton Sanderbeck et al., 2016; Oñorbe et al., 2017a; Keating et al., 2018; Puchwein et al., 2018; Oñorbe et al., 2018; Gaikwad et al., 2018; Wu et al., 2019). As a result, the thermal imprint of reionisation will last much longer than reionisation itself. IGM temperature measurements are therefore a potentially powerful way to probe the EoR, as they contain information about the reionisation history of a region that lasts long after it reionises. Modeling the temperature evolution during the EoR allows us to relate various parameters of ionisation history to IGM temperature. Comparing observations of temperature at various redshifts to the model will place constraints on the nature of the EoR and the sources driving it (Theuns et al., 2002; Raskutti et al., 2012; Lidz & Malloy, 2014; Oñorbe et al., 2017b; Boera et al., 2019).

Simulating the thermal history of the IGM is not a new idea, with many authors having studied IGM temperature under various assumptions of the background density and ionisation history. Instantaneous reionisation models (Theuns et al., 2002; Hui & Haiman, 2003; Upton Sanderbeck et al., 2016), radiative transfer (Puchwein et al., 2018; Keating et al., 2018), and inhomogeneous reionisation simulations (Furlanetto & Oh, 2009; D’Aloisio et al., 2019; Raskutti et al., 2012) have all been utilised as a basis for IGM temperature models. Once the density field and ionising background are known, most temperature modelling follows a similar process, where photo-heating, adiabatic cooling due to structure growth and the Hubble flow, along with various cooling processes in the IGM are followed over time. The process used in our model is outlined in 2.2.

This work uses the DRAGONS simulation suite, using the density grids from the N-body simulation Tiamat (Poole et al., 2016) and the ionising flux grids from the semi-analytic galaxy formation model Meraxes (Mutch et al., 2016) to determine temperature evolution over time. The ionising flux model in Meraxes captures the inhomogeneous nature of reionisation, allowing us to study the entire thermal structure in our simulation. Importantly, Meraxes directly couples high-redshift galaxy formation to hydrogen reionisation, allowing us to constrain the properties of ionising sources using their effect on the IGM. Using the temperature to probe the EoR involves comparing our simulated temperature distribution to observations of thermal history, placing constraints on the timing of reionisation, and the sources driving it.

Previous works studying IGM temperature have shown the relationships between the thermal state of the IGM and reionisation. Many of these works connected IGM temperature to statistics of the Lyman alpha forest and offered some measurements and constraints on the nautre of the EoR (Becker et al., 2011; Becker & Bolton, 2013; Walther et al., 2018; Boera et al., 2019; Oñorbe et al., 2018; Wu et al., 2019). The hydrodynamic and radiative transfer simulations used in this manner have made it possible to make measurements of the IGM temperature, and draw connections between thermal variables and the EoR. However these simulations are extremely computationally expensive. In order to compare these variables with a wide range of reionisation histories, a faster model is required. Using the DRAGONS simulation suite, we compare these observations to a wide range of reionisation scenarios in order to statistically constrain the global nature of the EoR.

The paper is structured as follows. The DRAGONS simulations will be briefly described and the post-processing temperature model will be laid out in section 2. Overview of the model outputs is given in section 3. Our results, including an investigation of the constraints the IGM thermal history can offer and constraints from current observations, are in sections 4 and 5 respectively, before concluding in section 6. The cosmology utilised throughout this paper is the flat standard $\Lambda$CDM from Planck Collaboration et al. (2016) with $\{\Omega_{\rm m},\Omega_{\rm b},\Omega_{\Lambda},h,\sigma_{8},n_{\rm s}\}=\{0.308,0.0484,0.692,0.678,0.815,0.968\}$.

2 Methodology

2.1 Model Inputs

Meraxes couples early galaxy formation and reionisation in a spatially and temporally dependent way, tracking gas cooling, star formation, and feedback between galaxies and the IGM amongst other processes (see Mutch et al. (2016) for more details). This paper utilises the 100 cMpc Tiamat simulation box (Poole et al., 2016) containing $2160^{3}$ particles of mass $2.64\times 10^{6}\rm{M_{\odot}}$. We use GBPTREES merger trees (Poole et al., 2017) from redshifts $z=35$ to $z=2$. We use the fiducial parameter balues presented in Qin et al. (2017), apart from the variations listed in section 2.3. Meraxes includes a modified version of the excursion-set algorithm 21cmFAST (Mesinger & Furlanetto, 2007) to track the progress of inhomogeneous reionisation in the simulated volume. Emissivity within an ionised bubble of radius $r$ in Meraxes is calculated from the star formation rate within the bubble, $\dot{m}_{*}(r)$

[TABLE]

where $f_{esc}$ is the fraction of ionising photons that escape the host galaxy, $N_{\gamma}$ is the number of ionising photons produced per stellar baryon, and $m_{\rm p}$ is the proton mass. The specific intensity at the hydrogen ionisation threshold, $J_{21}$, within the region is then computed from the emissivity,

[TABLE]

where $h_{\rm p}$ is the Planck constant. The comoving mean-free path, $\lambda_{\rm{mfp}}$, is equal to the ionised bubble radius during reionisation, and limited to $30\rm{cMpc}$ throughout the simulation, due to the maximum scale of the excursion-set algorithm111The value of $30\rm{cMpc}$ was chosen to roughly correspond to the mean-free path of ionising photons through an ionisied IGM at $z\sim 6$ (Sobacchi & Mesinger, 2013).. $\alpha$ is the assumed spectral power-law slope, a free parameter in our model where a small value corresponds to a harder UV background, such that

[TABLE]

where $\nu_{\rm{HI}}$ is the ionisation threshold of hydrogen.

Meraxes includes the quasar model detailed in Qin et al. (2017), where radiation from quasars is included when calculating reionisation structure and feedback, using equations analogous to 1 and 2, with a spectral slope of $\alpha_{q}=1.57$. Grids of the specific intensity from both galaxies and quasars, as well as density grids, are output from Meraxes. As stated in Qin et al. (2017) quasars have a sub-dominant effect on hydrogen reionisation in our model, due to the low number density of these luminous objects.

DRAGONS combines a mass resolution small enough to capture low-mass galaxy formation, a volume large enough to study the structure of reionisation, and an inhomogeneous reionisation model based on galaxy physics. Meraxes can simultaneously reproduce the observed stellar mass function, as well as Thomson scattering optical depth and ionising emissivity measurements with certain parameter choices (Mutch et al., 2016).

2.2 Temperature Model

In this paper we introduce an IGM temperature model in order to better constrain the EoR within DRAGONS. The model is largely based on the semi-numerical approaches of Raskutti (2011) and Hui & Gnedin (1997). Using the ionising background and density grids from the DRAGONS semi-analytic framework, we calculate the temperature and ionisiation of the IGM within the simulated volume. Non-equilibrium photo-ionisation rates, $\Gamma_{i}$ and photo-heating rates, $g_{i}$, are calculated post-ionisation from the specific ionising intensity in Meraxes, $J_{21}$, assuming an optically thin IGM:

[TABLE]

[TABLE]

where $\sigma_{i}(\nu)$ is the frequency dependent cross-section taken from Verner et al. (1996)222We integrate over the assumed power-law spectrum with 100,000 frequency bins between 1 and 4 Ryd. The results are not sensitive to the number of frequency bins, nor is it a limiting factor computationally, as the spectral slope is homogeneous, so this integral only needs to be performed once per model. The ionisation state of the IGM is then governed by the following differential equation

[TABLE]

for each species $i,j,k$ $\in$ {HI,HII,HeI,HeII} where $R_{jk}$ is the recombination rate of species $j$ and $k$ resulting in $i$, including recombination and collisional ionisation. $\tilde{X}_{i}$ is defined as $\tilde{X}_{i}=\frac{n_{i}m_{\rm p}}{(1+\delta)\bar{\rho_{b}}}$, for the local overdensity, $\delta$, and cosmic mean density, $\rho_{b}$.

The thermal state of the IGM is governed by the balance between photo-heating, adiabatic cooling under the Hubble flow, recombination cooling, and inverse Compton cooling, as well as changes in local overdensity according to Hui & Gnedin (1997)

[TABLE]

where $G$ is the total photo-heating rate of all species, $G=\sum_{i}g_{i}\tilde{X}_{i}$, $k_{b}$ is the Boltzmann constant, and $H$ is the Hubble parameter. The cooling rate, $\Lambda$, takes into account recombinations, collisions, bremsstrahlung, and inverse Compton cooling. We take the rates for these processes from Lukić et al. (2015). Photo-ionisation and heating rates are calculated separately for stellar and quasar sources using the optically thin equations 4 and 5, then added together when solving equations 6 and 7.

Following Raskutti (2011), the coupled equations 6 and 7 are solved recursively, without assuming ionisation equilibrium, using a first order implicit integration scheme (Anninos et al., 1997; Bolton & Haehnelt, 2007) until they converge to a solution with $|\tilde{X}_{e,k}-\tilde{X}_{e,k+1}|<10^{-6}$ at a given attempt $k$, where we use the electron abundance as our convergence statistic. To improve the efficiency of our code, we adopt a variable timestep, where the timestep length is doubled for the next timestep each time a solution is found, or halved if a convergent solution cannot be found within 100 attempts.

In the same manner, we follow the integrated thermal history via $u$, the total energy injected into the IGM per unit mass via photo-heating. This can be observed in Lyman alpha power spectra, distinct from temperature (Nasir et al., 2016; Boera et al., 2019) and can be used to simultaneously measure reionisation timing and amount of photo-heating that exists when only considering mean temperature. The injected energy $u$ is followed by simply integrating the photoheating rate over time. The value of $u$ has been related to the small scale Jeans smoothing of the IGM, as it is dependent on the integrated thermal history throughout the EoR. We cannot calculate the Jeans scale directly in post-processing, but $u$ provides a similar probe into the integrated thermal history throughout the EoR.

In order to achieve the computational speeds required to run the model many times, we track the temperature within a $128^{3}$ grid with cell side length $\approx 800ckpc$ of the Tiamat $100$ cMpc box; from the redshift of reionisation of each voxel, until $z=4$. The results for this paper use the same 10,000 (approximately 0.5%) randomly selected voxels in each box, unless otherwise stated, as a sample of the entire volume. As this is a post-processing model, the thermal state of each voxel is treated independently, although their ionising flux intensities and densities are already coupled within Meraxes and Tiamat.

We set the specific intensity above the helium ionisation threshold $\nu_{HeII}=54.4$ eV to zero, so that there is no reionisation of HeII to HeIII. This is because Meraxes only traces the size of HII bubbles, meaning it does not predict the mean-free path of photons above $\nu_{HeII}$. This will restrict our temperature model to times earlier than HeII reionisation, thought to complete around $z\sim 3$ (Furlanetto & Oh, 2008). We include HeI reionisation, since we expect helium to be singly ionised at the same time as hydrogen (Wyithe & Loeb, 2003). When constraining the EoR, we also ignore the outputs of our model below $z=4$, to minimise confusion with the effects of HeII reionisation.

2.3 Parametrising Reionisation

In order to produce and test a large number of thermal and ionisation histories, we vary three parameters within Meraxes and the post-processing temperature model. Escape fraction normalisation and redshift-scaling, as well as background ionising spectral slope.

The temperature of the IGM is sensitive to the timing of reionisation, and the timing of reionisation is heavily dependent on the escape fraction of photons from galaxies. We utilise a redshift-dependent, uniform escape fraction for ionising photons, which was shown by Mutch et al. (2016) to allow the model to match electron optical depth and ionising emissivity observations simultaneously.

[TABLE]

We vary $f_{5}$ between 0.03 and 0.12 and $\beta$ between 0 and 2.5 to vary the timing and duration of reionisation in the model. These values were chosen to bracket constraints from electron optical depth measurements, producing reionisation histories that finish between redshifts 5 and 10.

The spectral shape of the ionising background sets the energy injected per photoionisation, which affects the reionisation temperature and the subsequent cooling rate. We model the stellar spectrum between 13.6 and 54.4 eV as a power law (equation 2), with a slope, $\alpha$, between 0.2 and 5. The quasar spectral slope in the same frequency range is fixed at 1.57 and the quasar escape fraction is fixed at 1, as in Qin et al. (2017).

The range of spectral slopes considered is both broader and softer than those often used in temperature modelling (Upton Sanderbeck et al., 2016; D’Aloisio et al., 2019), which are based on Population II stellar synthesis models. This range was chosen to produce at least one thermal history that is consistent with observations for each reionisation history. We investigate these scenarios when studying the correlations between heat injection, reionisation timing and temperature. When placing constraints on the EoR, however, we restrict $\alpha$ to be more consistent with these population synthesis models.

2.4 Initial Conditions

For our fiducial model, we start with a 99 per cent ionised (HII and HeII) IGM, with the initial temperature calculated from the UV spectral slope at ionisation and the speed of the ionisation front in Meraxes using fits to radiative transfer simulations, performed by D’Aloisio et al. (2019). If the ionisation front passes through the gas very quickly, the reionisation temperature is decided entirely by the average energy of the ionising photons $\langle E_{\gamma}\rangle$. (Keating et al., 2018; Hui & Gnedin, 1997), yielding

[TABLE]

However, D’Aloisio et al. (2019) found using one dimensional radiative transfer simulations that the ionisation front can pass through slowly enough for collisional cooling within the hot, partially neutral gas to have a large effect on the reionisation temperaure. It was found that the speed of the ionisation front provided the best estimation for reionisation temperature, as the faster the ionisation front passes through the gas, the less time it spends in the hot, semi-neutral state where collisional cooling is efficient.

Ionisation front speeds are calculated in each voxel by finding the distance between ionisation boundaries (where adjacent voxels have different ionisation snapshots) at successive snapshots, and assuming that fronts travel at a constant speed within each 11Myr snapshot. We take the distance between random points in each voxel, representing our uncertainty at the grid resolution. The reionisation temperature is calculated from the front speeds and the spectral slope of the background, using fits provided by D’Aloisio et al. (2019). The distribution of reionisation temperatures in our box is presented in section 3.2, this approach introduces a correlation between ionising flux amplitude, gas density, and reionisation temperature, which further complicates the picture of patchy reionisation.

Since we explore softer spectra than D’Aloisio et al. (2019), initial temperatures for models with $3<\alpha\leq 5$ are given by the lowest of our two upper limits; from 1) the temperature at the speed of the ionisation front with $\alpha=3$ and 2) the maximum temperature given by the spectral slope in equation 9. This will slightly overestimate initial temperatures for the slower moving fronts, however softer spectra approach their maximum reionisation temperature at slower speeds, so this effect will be small.

The total photo-heating energy, $u$, is initialised to the mean excess energy of ionising photons, $E_{\textrm{excess}}$ (Puchwein et al., 2018). Assuming total local absorption of ionising photons,

[TABLE]

where

[TABLE]

is a factor that accounts for the preferential absorption of different frequencies by different species, and $j$ sums over HI, HeI and HeII. Regardless of how long it takes the front to move through the gas, the total amount of energy imparted to the gas will depend only on the background spectrum, assuming the number of recombinations and collisional ionisations that occur as the front passes is small, and that all photons with energies between 1 and 4 Ryd are absorbed within the front.

We note that, in contrast to other works, the reionisation temperature is partly decoupled from the average photon energy. This occurs via collisional excitation cooling described above, from the models in D’Aloisio et al. (2019), whereas most previous works assume reionisation occurs very quickly in each region with no cooling in the ionisation front. Using these models lowers Treion without changing the initial $u_{0}$, as $u_{0}$ only takes into account energy changes via photo-heating, as defined by Nasir et al. (2016) and Boera et al. (2019). As a result, we set the initial $u_{0}$ to the mean excess energy of the ionising background, described above.

2.5 Sub-grid Clumping

Gas on scales below our grid resolution is not homogeneous, and clumping on unresolved scales will increase recombination rates. Increasing the recombination rate will increase temperatures after reionisation by shifting ionisation equilibrium to a more neutral state, allowing more photoionisations to occur, which overcomes the increased cooling rates. Large scale clumping is taken into account via our spatial grid, of voxel length $\sim 800$ckpc. We also expect reionisation to erase much of the small scale clumping due to Jeans smoothing. This still leaves some room for subgrid clumping to affect our results, on scales between our grid resolution and the Jeans length.

The increased temperatures due to clumping, while reflective of the total heat inside a voxel, will not represent the wide distribution of temperatures that can exist within the voxel, as it will be dominated by the higher density clumps. In order to compare with higher resolution simulations (see appendix A), and produce a temperature independent of grid resolution, we set the clumping factor to 1. This implies that we are following the gas at the mean density of each voxel throughout the simulation, rather than each voxel as a whole.

We note that there is an effect due to the subgrid structure on the background spectrum due to spectral filtering. This is where intermediary absorbers harden the ionising background by preferentially absorbing lower frequency photons, as the HI ionisation cross-section scales approximately as $\nu^{-3}$ near the ionisation edge. The amount of hardening is dependent on the structure and clumping of the IGM, as well as the structure or the ionising background. Theoretically this could harden the spectral slope by 3 (see Faucher-Giguère et al. (2009) appendix D) but measurements of the column density distribution of absorbers suggest a hardening of the spectral slope by approximately 1 (Songaila & Cowie, 2010). However modeling this filtering, as well as the spatially dependent spectrum of the ionising background, is beyond the scope of this work.

2.6 Model Caveats

Two approximations are made when evaluating temperature evolution that should be noted. First, equations 6 and 7 are applied to the Meraxes grids, rather than individual parcels of gas. Second, we calculate temperatures in post-processing, this assumes independence of each voxel with regards to the temperatures of other voxels, and that any gas influx is at the same temperature as the gas within the voxel. These assumptions can cause some spurious heating or cooling as gas moves through the box, since this will be treated in the same way as structure growth (via the third term of equation 7). Since neighboring voxels tend to have similar ionisation histories and densities, combined with dark matter velocities that are fairly low compared to the voxel sizes, we do not expect this to have a large effect on temperatures.

We tested the effect of gas diffusion on the model by running a model where an extra term was added to equation 7 to account for the gas of differing temperature entering the voxel. We keep track of the bulk flow of matter via the Tiamat velocity grids, and conservatively assumed that the temperatures of adjacent regions scaled with the maximum temperature density slope found in our other models, $T\propto(1+\delta)^{0.6}$. Under these assumptions, heating rate differences of order $\sim 10\%$ were observed in high density regions, and changes of order $\sim 1\%$ were observed at mean density. Considering that this model overestimates the heating changes due to gas diffusion (nearby voxels are likely to be at similar temperatures shortly after reionisation), we believe that the independence of voxels is a safe approximation for the purposes of this work.

The independence of voxels within our model also excludes a treatment of recombination emission, since inhomogeneous recombinations are not yet included in Meraxes. The effects of recombination emission are detailed in Faucher-Giguère et al. (2009) where they find it can contribute up to $10\%$ of the ionising background for hydrogen. The effect of a $10\%$ increase in the photo-ionisation rate will not greatly affect our results (see Appendix B). However, recombination emission could soften the background spectrum by $\approx 1$, since the photons from hydrogen recombination will tend to have lower energies than those in the rest of the ionising background. While this would cause a decrease in temperatures, the effect is degenerate with our free parameter for the background spectral slope, so modelling this effect is outside the scope of this work.

We assume that there is negligible heating by $\rm{HeII}$ reionisation before $z=4$ (Upton Sanderbeck et al., 2016). If Helium reionisation is a highly extended process, then we would be underestimating the temperature, and overestimating the IGM cooling rate for any given reionisation history. With our parameterisation, this would bias our models to earlier reionisation scenarios (due to the flatter temperature gradient), with harder spectral slopes (due to the higher temperatures).

The reionisation model in Meraxes assumes local absorption of ionising photons. Photons that would redshift below the ionisation threshold in a full radiative transfer model can still ionise in Meraxes. Since Meraxes is tuned to match the ionising emissivity measurements of Becker & Bolton (2013), which were calculated based on radiative transfer models, $J_{21}$ could be overestimated by a small amount. This will not have a large effect during Hydrogen reionisation, as the size of HII regions are small enough to make local absorption a good approximation. Furthermore, we find that temperature at all redshifts is insensitive to the value of $J_{21}$ after the region is reionised, as long as it is large enough to maintain a highly ionised IGM (Appendix B).

We only vary the escape fraction (equation 8) to produce different reionisation scenarios, while keeping the same source model. This means that we ignore any degeneracy between our escape fraction parameters and source modelling with regards to the IGM thermal state. In particular the constraints we find due to the scatter in temperatures are likely optimistic, as the patchiness of reionisation will have a significant effect on the range of temperatures observed afterwards.

Shock heating of high density regions as they collapse is also not modeled in this work, as we are primarily interested in the diffuse IGM. We only track the increase in temperatures from structure growth on the scale of our grid. This will exclude the $\sim 10^{5}K$ gas that exists in hydrodynamic simulations. However, studies of IGM temperature have so far focused primarily on gas at or below the critical density, excluding shock heated gas in their analyses (Oñorbe et al., 2018).

3 IGM Thermal History

Using the temperature evolution model described in section 2, we calculated the thermal and ionisation states of the same 10,000 randomly chosen voxels for 2750 Meraxes realisations, in order to obtain a wide range of density and ionising flux histories per model, as well as a wide range of global reionisation and thermal histories. We have also computed the evolution of one full box, to examine the topology of the temperature field. The model parameters used in the illustrative examples within this section are $\{f_{5},\beta,\alpha\}=\{0.08,1,2\}$ on the full $128^{3}$ Meraxes box. This run was chosen as all the correlations between reionisation timing, density and temperature of reionisation are clearly shown in its results. However, this is not our highest likelihood model, based on observational data.

A few examples of thermal histories output by our model are given in Figure 1. Each voxel will experience a sudden rise in temperature when it is ionised, dominated by photo-ionisation heating. Afterwards, the voxel cools on a cosmological timescale at a rate determined mainly by the ratio of photo-heating and recombination cooling at ionisation equilibrium, as well as adiabatic cooling and inverse compton cooling off the CMB. Density evolution will modulate the temperature, but the previous density history of a voxel will have little effect on its thermal asymptote (McQuinn & Upton Sanderbeck, 2016). However, the previous density history will have a substantial effect on the integrated photo-heating, $u$. Each voxel approches its thermal asymptote within $\Delta z\approx 1-3$ of its reionisation, creating a distribution in temperatures in the whole box dependent on the inhomogeneous reionisation history.

3.1 Topology

Figure 2 shows a $128^{2}$ voxel slice of the temperature model at selected redshifts between 6 and 4. It can be seen that shortly after reionisation, at $z=6$, the high density regions reach temperatures of $T\approx 2\times 10^{4}K$, as they have recombination rates high enough to allow continuous photo-heating of the gas. The larger voids also reach similar temperatures, as they reionise last due to the “inside-out" nature of the EoR in our models, and hence have the least time to cool. The coolest regions are the low density regions in close proximity to high density regions, which are reionised early by stars in the nearby dense filament, and cool over time as their recombination rates are not high enough to maintain high levels of photo-heating. Long after reionisation, at $z=4$, the hotter low density regions cool so that temperature is closely correlated with density. The ionisation topology matches that of previous works modeling inhomogeneous reionisation (Trac et al., 2008; Raskutti et al., 2012; Keating et al., 2018; Lidz & Malloy, 2014; Oñorbe et al., 2018). Comparing Figure 2 with the left panel of Figure 3, which show the redshifts of reionisation and temperatures within the same region, we see that shortly after reionisation there is a notable anti-correlation between temperature and redshift of reionisation. Over time this correlation diminishes, and the correlation between temperature and density becomes dominant.

3.2 Reionisation Temperature

In this section we present the reionisation temperatures in our fiducial model, using fits for the reionisation temperature provided by D’Aloisio et al. (2019), to convert ionisation front velocities to reionisation temperatures. Using these fits allows us to include the correlations between temperature, ionisation history and the photon background in more detail; this reduces the number of free parameters we need to include.

The temperature to which a region of the IGM is reionised to is inversely correlated with its redshift of reionisation, as shown in Figure 3. This anti-correlation results from the fact that the ionisation fronts speed up as they enter the low-density regions towards the end of the EoR, resulting in higher temperatures because there is less time for collisional cooling to take effect. Differences in simulated reionisation history, and to a lesser extent grid resolution, snapshot cadence, and algorithm to find ionisation front speed result in changes to the distribution. Figure 4 shows the distribution of reionisation temperatures, compared to redshift of reionisation. As shown in the top panel, reionisation in this model occurs in the redshift range $6\lesssim z\lesssim 12$. The middle panel shows the full distribution of temperatures to which regions ionise. We find a similar range of reionisation temperatures as D’Aloisio et al. (2019), between $(15-25)\times 10^{4}$ K. The bottom panel splits up this distribution into regions that ionised at different times, showing that the temperature of the ionisation fronts increases as reionisation progresses. Lower spatial resolution, as well as our discrete timestep, likely smooths out the extreme ends of the distribution compared to D’Aloisio et al. (2019), and results in a broader distribution of temperatures at any given redshift. However it is difficult to directly compare the distributions due to differences in reionisation history.

The excursion set formalism used to predict reionisation history in Meraxes is likely too crude to predict accurate ionisation front speeds on scales of individual voxels, as the spherical symmetry assumed by 21cmFAST effectively averages over many cells during the later stages of reionisation. However, since we recover a similar range of speeds and reionisation temperatures as D’Aloisio et al. (2019), and the reionisation temperature increases as reionisation progresses, we consider this approach to be a valid approximation for testing the statistics of IGM temperature and more accurate than assuming a constant reionisation temperature for all voxels.

3.3 Temperature-Density Relation

Figure 5 shows the temperature-density relation at various redshifts in our fiducial model. The temperature-density relation (TDR) is a powerful probe of the conditions in the IGM. The shape and scatter in the relation at various densities can reveal information about ionisation history and structure of the IGM during the EoR. A power law fit $T=T_{0}(1+\delta)^{\gamma-1}$ is commonly used when characterising the thermal state of the IGM. However, this is only accurate for regions that ionised homogeneously or long before the measurement is taken. Restricting the study of temperature to a power law misses much of the information in the temperature density relation that can be used to constrain a patchy EoR (Oñorbe et al., 2018; Wu et al., 2019).

There is a wide distribution of temperatures in low density regions shortly after reionisation, as the temperature at low density is highly dependent on the redshift of ionisation at this time. The large scatter lasts long after reionisation itself, halving approximately within $\Delta z\approx 1$ after reionisation in most of our models. Long after reionisation, regions of all densities cool to their asymptotic temperature, and the temperature density relation is well described by a power law $T=T_{0}(\delta+1)^{\gamma-1}$, with a slope approaching $\gamma-1\approx 0.6$ by redshift $z=4$. This is consistent with other inhomogeneous reionisation models and analytic calculations of the thermal asymptote (Hui & Gnedin, 1997; McQuinn & Upton Sanderbeck, 2016). High density regions do not show much variance in their temperatures, as their reionisation temperatures tend to be much closer to their early-reionisation asymptote. As a result high density regions settle much more quickly into the late time power-law. Furthermore, dense regions tend to ionise earlier, so any scatter at the high end of the temperature density relation will likely have disappeared by the time of measurement.

In agreement with other inhomogeneous reionisation models (Trac et al., 2008; Raskutti et al., 2012; Keating et al., 2018), we find that the large scatter in low density regions shortly after reionisation is highly correlated with the redshift of ionisation of the region (shown in the colours in Figure 5). The size of this scatter, and its correlation with $z_{r}$, is what causes the changes in the $T-\delta$ slope. The scatter in temperatures can be a powerful probe of the patchiness of reionisation, as it shows differences in reionisation redshift between different regions of the low-density IGM, allowing us to quantify the patchiness of reionisation, as well as estimate how long ago reionisation occurred as the scatter diminishes (see section 4.1). Once every region has reached its thermal asymptote, and the TDR settles into a power law, the thermal memory of reionisation has been lost, as regions that ionised at different times approach the same temperature, based on their density.

Trac et al. (2008); Furlanetto & Oh (2009) and Raskutti et al. (2012) note an inversion of the temperature-density relation, with $\gamma-1\sim-0.2$, at low densities shortly after reionisation, meaning the lowest density regions are hotter on average than mean density regions at the end of the EoR. This inversion is due to the lowest density regions ionising last, hence having less time to cool, as well as ionising to higher temperatures due to faster ionisation front speeds. In the above model, we find this inversion towards the end of reionisation, lasting until $z\approx 5$ as there are many hot low-density regions that have recently ionised, and have yet to cool towards the asymptotic power-law relation. The average slope of the TDR is highly dependent on reionisation history and spectral slope, as they alter the strength of the correlations between density, redshift of reionisation, and temperature. As a result, inversion in the TDR is not observed in all of our models, although the slope will always be at a minimum towards the end of the EoR.

4 Constraining the EoR Using Temperature

In this section we investigate how the distribution of temperatures in the IGM can be used to probe reionisation. We have run a suite of realisations of our model, with different escape fractions and background spectral slopes (see section 2.3), controlling the timing and duration of reionisation, as well as the temperature of reionisation and subsequent cooling rates.

In order to simulate earlier and later reionisation scenarios, we vary the escape fraction normalisation and redshift scaling (equation 8) in Meraxes. A higher (lower) escape fraction results in more (less) ionising photons escaping galaxies, and therefore an earlier (later) reionisation, and a colder (hotter) IGM at a fixed redshift, while thermal memory of reionisation still exists. Models with different escape fractions will tend towards the same temperature at late times, as all regions approach their thermal asymptote. The effects of the escape fraction parameters on global redshift of reionisation (defined as when the global neutral fraction is 10 percent) can be seen in Figure 6. Reionisation history also has an effect on reionisation temperatures, by changing the speed of ionisation fronts in the IGM.

Figure 7 shows the temperature at mean density, it’s standard deviation, and the energy injected from photoionisations at $z=4.5$ versus the model’s global redshift of reionisation (where the mass-weighted neutral fraction falls below 0.1) for varying spectral slopes. As noted in section 3.3, both the mean and spread of temperature are maximised shortly after the bulk of reionisation occurs. The mean temperature will decrease towards an asymptote dependent on the background spectrum, and the scatter will decrease towards zero. A harder spectral slope will impart more energy to the IGM on average per ionising photon, increasing the ratio of the photoheating rate to the ionisation rate. As a result, the thermal asymptote of each voxel becomes hotter. The reionisation temperature also increases, however this can be limited by collisional cooling if the ionisation front passes through the IGM slowly (see section 3.2 and D’Aloisio et al. (2019)).

Changes in temperature due to reionisation timing and those due to background spectrum are difficult to differentiate using temperature measurements alone. However, mean temperature and scatter in temperature have different correlations with the timing of reionisation and background spectrum, and will evolve at different rates over time. As a result, observations of mean and scatter in IGM temperature at multiple redshifts can be used to break the degeneracy between the timing of reionisation and the background spectrum, offering tighter constraints on the redshift of reionisation.

The integrated photo-ionisation energy at mean density, $u_{0}$, can be used to further tighten constraints. Unlike temperature, an earlier reionisation increases the photo-ionisation energy at mean density, $u_{0}$, since residual photo-heating that occurs after reionisation will have been occurring for longer. Because injected energy and temperature have opposite correlations with redshift of reionisation, their observations can be used together to simultaneously constrain the background spectral slope and the timing of reionisation (Nasir et al., 2016; Boera et al., 2019).

4.1 Mock Observations

In order to explore how temperature observations constrain the EoR, we perform mock observations on our simulation. We create mock observations on the model discussed in section 3 of temperature at mean density $T_{0}$, scatter of temperature at mean density $\sigma(T_{0})$, and injected photo-ionisation energy at mean density $u_{0}$ at redshifts 4.5, 5.0 and 5.5, similar to the redshifts of the temperature measurements from Boera et al. (2019), which are 4.2,4.6, and 5.0. $T_{0}$ and $u_{0}$ are derived from $T$ and $u$ by binning their respective distributions around mean density $\Delta=0\pm 0.1$ and averaging. We use 0.125 dex as a $1\sigma$ measurement error, of similar size to the largest uncertainties given in Boera et al. (2019). By comparing the mock observations to our series of models, we then estimate the timing of reionisation and background spectral slope, to see how well the input parameters can be recovered. We estimate the likelihood of each model assuming Gaussian errors, based on the true values from one model. The likelihoods of each of our models, compared to these mock observations from our model (with parameters $\{f_{5},\beta,\alpha\}=\{0.08,1.00,2.00\}$, and $z_{\rm{reion}}\approx 6.84$), are shown in Figure 8. It is important to note that we do not create mock Lyman alpha spectra, due to the low resolution of our simulations. Here, “mock observations" refers to the summary statistics of the IGM thermal state, that we generate from one of our models with some assumed error margin. These statistics correspond to the temperature and energy measurements attained from the analysis of Lyman alpha spectra, for example those published in Becker et al. (2011); Walther et al. (2018) and Boera et al. (2019).

We recover a strong peak in $\alpha$, and as seen in the contours of the $\beta$ versus $f_{5}$ panel, the highest likelihood models are those with the same $z_{\rm{reion}}$ as the true model. However, we cannot recover our two escape fraction parameters independently with this sample, as our observables are more sensitive to the timing of reionisation than the duration of it. We will require more precise observations, across a wider range of redshifts to begin to distinguish these scenarios using the scatter in temperature. A more extended reionisation will have a wider distribution of temperatures, from a wider distribution of reionisation times.

We next present the constraints from each observable separately. The different correlations of temperature at mean density, its scatter, and injected energy with redshift of reionisation and background spectral slope allow us to constrain these reionisation parameters simultaneously. A similar result was demonstrated by Boera et al. (2019), using the integrated thermal history to break the degeneracy between reionisation timing and initial temperatures. Figure 9 illustrates this for the $\{f_{5},\beta,\alpha\}=\{0.08,1.00,2.00\}$ model, showing the constraints from mock observations of mean temperature, scatter and photo-heating separately on the redshift of reionisation and background spectrum. While each observable alone is degenerate between the timing of reionisation and the background spectrum, their degeneracies differ in magnitude and direction, allowing us to perform this analysis.

We have applied this analysis to each of our models and present the precision achievable for recovery of reionisation redshift and spectral slope in Figure 10. We can recover the redshift of reionisation in our model within $\Delta z\approx 1$ and spectral slopes within $\Delta\alpha\approx 2$ to 95% confidence, using an error margin of 0.125 dex for temperature and injected energy observations.

4.2 Future Observations

As discussed in section 4.1, our mock observations with errors corresponding to existing datasets were unable to independently recover escape fraction normalisation and scaling. We now examine a much more optimistic dataset from possible future observations, where temperatures at mean density, scatters, and injected energies are known within a $1\sigma$ error of 0.05 dex at 9 evenly spaced redshifts between 4 and 6. We show the results for this case in Figure 11 for the parameter set $\{f_{5},\beta,\alpha\}=\{0.08,1.00,2.00\}$ and Figure 12 for all models. The recovery of parameters is much more precise in this scenario. With such an extensive dataset we can recover the redshift of reionisation within $\Delta z\approx 0.5$ and background spectral slope within $\Delta\alpha\approx 0.5also$ to 95% confidence333We note that our model grid is not sufficiently fine to resolve details of the 2d parameter space, however the likelihood values illustrate the improved constraints. This example is also closer to recovering the escape fraction normalisation and scaling independently. These mock measurements of the scatter in temperatures disfavour reionisation scenarios that are too extended or too sudden, however there is still a significant degeneracy between these two recovered parameters. Considering that this is a very optimistic dataset, recovering the duration of reionisation from its thermal state will likely require a more detailed analysis.

5 Constraints Using Existing Observations

In the last section, we demonstrated how IGM thermal state observables constrain the nature of the EoR using mock observations of our suite of models. In this section, we apply this process to recent observations of the thermal and ionisation history of the universe, in order to constrain our model parameters. We constrain against measurements of temperature at mean density from Boera et al. (2019), and electron optical depth measurements from Planck Collaboration et al. (2018). Temperature measurements from Becker et al. (2011) and Walther et al. (2018) were also consdidered, but model likelihoods on these datasets will not be presented in this paper. The former dataset produces similar but looser constraints on reionisation, whereas we are unable to match temperatures from the latter dataset due to the sudden increase in temperatures at $z=5$. Figure 13 shows the observations utilised, and maximum likelihood models based on each temperature dataset. Electron Optical depth measurements place constraints on the redshift of reionisation444Electron optical depth can also be used to constrain both the timing and duration of reionisation, though not independently (Greig & Mesinger, 2017). However, we only constrain $z_{\textrm{reion}}$ in this work, while temperature measurements constrain both the timing of reionisation and the background slope.

We restrict our attention to temperature at mean density, $T_{0}$, and electron optical depth, $\tau_{e}$. No measurements currently exist for scatter in temperatures, so we are unable to use this observable for our EoR constraints, as in our mock examples. Regarding $u_{0}$, we are unable to reliably relate the integrated thermal history in the optically thin models that produce this measurement in Boera et al. (2019), with the thermal histories in our patchy reionisation models. A patchy reionisation effectively has a much harder spectrum within ionisation fronts, reducing the effective spectral slope by approximately 3 due to all ionising photons being absorbed within the front555Assuming the number of recombinations in the ionisation front is small, the photoheating energy at ionisation is set by $E_{\rm{excess}}$ (equation 10, see section 3.2) regardless of front speed. In a optically thin reionisation, both integrands would have a factor of the ionisation cross-section $\sigma_{\rm{HI}}$, which scales approximately as $\nu^{-3}$, softening the effective UV spectrum. This creates a hotter reionisation, but most importantly, permanently offsets the injected energy by a certain amount compared to an optically thin model, because $u_{0}$ is a time-integrated statistic. Since we cannot model how $u_{0}$ affects the small scale Lyman alpha power spectrum in our patchy reionisation models, we do not include $u_{0}$ in our fiducial constraints.

We have ignored temperature measurements at redshifts $z<4$ to minimise confusion resulting from the beginning of HeII reionisation. If a substantial amount of HeII ionisation has occurred at the redshift of measurement, this could be confused with a hotter post HI reionisation IGM with a flatter or positive evolution. This would result in a bias towards models where the overall cooling rate is suppressed, either due to a harder ionising spectrum or earlier reionisation, where the gas is closer to its thermal asymptote and would show a flatter evolution.

Given the strong degeneracy between escape fraction parameters, we instead show constraints on $\alpha$ and the redshift of reionisation, $z_{\rm{reion}}$. Figure 14 shows our constraints on these parameters based on the measurements from Boera et al. (2019) and Planck Collaboration et al. (2018). We find a strong degeneracy between $z_{\rm{reion}}$ and $\alpha$ where $T_{0}$ alone is considered. Using $\tau_{e}$ to break the degeneracy results in a late reionisation $z_{\rm{reion}}=6.8^{+0.5}_{-0.8}$ and a soft ionising background $\alpha=2.8^{+1.2}_{-1.0}$.

In order to match temperature observations at $z>4$ as well as observations of electron optical depth, our simulations favour a softer UV spectrum than assumed in other works. This may result from tension between the observations, with temperature observations favouring an earlier reionisation than electron optical depth measurements. This is certainly the case when we restrict the background spectral slope to ranges supported by stellar population synthesis models, with $0.5<\alpha<2.5$, given by the red shaded region in figure 14 (see Upton Sanderbeck et al. (2016) and D’Aloisio et al. (2019) for discussions of possible background spectral slopes) which tightens our reionisation timing constraints to $z_{\rm{reion}}=6.9^{+0.4}_{-0.5}$. However many models without this restriction, with $\alpha<2.5$ and a slightly earlier reionisation $z_{\rm{reion}}\approx 7$ are within $1\sigma$ uncertainty of our results. Since the duration of reionisation also has an effect on temperatures via the ionisation front speed (see section 3.2, D’Aloisio et al. (2019)), models with harder spectra require slower reionisation histories in order to agree with both temperature and CMB measurements (Boera et al., 2019). As a result constraining reionisation duration will result in tighter constraints on spectral slope. We also note that spectral slopes $\alpha<0.5$ are ruled out at $2\sigma$, as these spectra produce temperatures that are too high compared to observations, even when reionisation occurs relatively early.

6 Conclusions

Using an inhomogeneous reionisation model, we have probed the full distribution of IGM temperatures, and its correlations with structure growth, the ionising background and patchy reionisation. We have begun to show how these correlations can be used to characterise the EoR, placing simultaneous constraints on the timing of reionisation and the background spectrum. We recover thermal behaviours from other inhomogeneous reionisation models (Raskutti et al., 2012; Trac et al., 2008; Oñorbe et al., 2017b; Keating et al., 2018; Furlanetto & Oh, 2009; Oñorbe et al., 2018), where shortly after reionisation, there is a large scatter in the low-density end of the temperature-density relation, and a strong correlation between temperature and redshift of ionisation. By redshift z = 4, the temperature correlation shifts towards density, creating a power law temperature density relation, leaving no memory of the redshift of reionisation in the low-density IGM. Our main findings are as follows:

• •

The large initial scatter in the temperature density relation is highly correlated with the redshift of reionisation of the region, differentiating between regions near the cosmic filaments, which reionise early, and those in large voids, ionising later. Temperature measurements taken while the low-density scatter still exists contain information concerning when a region reionised, and the mean and scatter of multiple measurements can be used to constrain the reionisation history. However interpretations of these measurements will be highly dependent on the nature of the ionising sources, and the structure of the gas in the region.

• •

Performing mock observations of the mean and scatter of the temperature density relation at different redshifts illustrates the potential for constraints on the EoR. Using uncertainties similar to those for recent temperature measurements (Boera et al., 2019; Walther et al., 2018), we can recover the redshift of reionisation from our mock observations to $1\sigma$ within $\Delta z_{\rm{reion}}\approx 0.5$ and background spectral slopes within $\Delta\alpha\approx 1$.

• •

Comparing our suite of models to electron optical depth and temperature measurements, our modeling favours a reionisation history that finishes around $z_{\rm{reion}}=6.8^{+0.5}_{-0.8}$ and a UV background with a power-law spectral index of $\alpha=2.8^{+1.2}_{-1.0}$ ($1\sigma$ uncertainties) between 912 Å and 228 Å. If we restrict our models to those with spectral slopes consistent with population II stellar sysnthesis models ($0.5<\alpha<2.5$), the redshift of reionisation is restricted further to $z_{\rm{reion}}=6.9^{+0.4}_{-0.5}$.

Knowledge of the distribution and history of temperatures within the IGM in addition to temperature at mean density will greatly improve our understanding of the EoR. while the timing of reionisation is probed by mean temperature, as discussed in section 4.1, an estimate of the scatter in temperatures between different lines of sight from multiple quasar spectra, or by the large scale features of the Lyman alpha forest (Oñorbe et al., 2018; Wu et al., 2019), would begin to offer information about the duration and patchiness of reionisation. This is because the distribution of temperatures shortly after reionisation is closely related to the distribution of reionisation redshifts within a volume (see Figure 5). This is especially true for low density gas, where the distribution of temperatures is widest.

One way to further probe the temperature density relation could involve the cross-correlation of temperature measurements with the galaxy field. In this manner, information could be gained on the correlation of temperatures with density at various scales throughout reionisation. From the temperature-density plots in this work, we would expect to see the correlation on small scales flattening, or even inverting, shortly after reionisation.

Acknowledgements

JD would like to thank James Bolton for helpful comments on the draft version of this paper. JD is supported by the University of Melbourne under the Melbroune Research Scholarship (MRS). This research was also supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. This work was performed on the OzSTAR national facility at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). Parts of this work were performed on the gSTAR national facility at Swinburne University of Technology. gSTAR is funded by Swinburne and the Australian Government’s Education Investment Fund. Parameter constraints plots were generated using the corner python package (Foreman-Mackey, 2016).

Appendix A Numerical convergence

This section explores the effects of numerical factors in our model. The grid resolution of Meraxes will be investigated, as well as the parameters used in our DE solver, namely the convergence threshold that we use to determine the timestep throughout the model.

For computational reasons, all of model runs were performed using $128^{3}$ grids in Meraxes. However to ensure inhomogeneities below this scale do not greatly affect our results, we examine results from one run on a $256^{3}$ grid. As shown in Figure 15, while the overall clumping factor has increased, the temperatures and photo-heating energy at mean density are largely unaffected. This is due to the fact that we follow the gas at the mean density of each voxel, rather than the voxel as a whole. This will produce a very similar temperature-density relation regardless of grid resolution; The distribution of densities will change, but the temperature at any particular density will remain the same.

We have altered the convergence conditions in our differential equation solver to test for convergence. We change the conversion threshold between $|\tilde{X}_{e,k+1}-\tilde{X}_{e,k}|<10^{-8}$ and $10^{-4}$. The results of these changes for an individual voxel that ionised at $z=10$ are presented in Figure 16. The weakest threshold gives a maximum difference in temperatures of $\approx 2\%$, and our fiducial threshold of $10^{-6}$ is negligibly different from the strictest thresholds. We are therefore confident that our choice of differential equation solver parameters does not affect our results.

Appendix B Variation with Ionising Flux Amplitude

Previous simulations have found little dependence of the long term thermal state on the amplitude of the ionising background, as long as it is strong enough to maintain an ionised IGM (Hui & Gnedin, 1997; Furlanetto & Oh, 2009). This result is verified in our model, as we can see no change in the long term thermal state when we artificially increase or decrease the amplitude of the ionising background post-reionisation, keeping our other parameters constant. The ionisation rate in a voxel is directly proportional to the amplitude, however the number of ionisations that actually occur is limited by the recombination rate, which is negligibly changed for all highly ionised states. As a result the long term cooling is largely unaffected at most densities. Figure 17 shows the mean temperature at mean density, a third of mean density, and three times mean density, for the three amplitudes tested.

It is important to note that the ionising flux amplitude has little effect on temperature only when there is enough ionising flux to maintain the ionised state in a region. If the amount of ionising flux is low enough such that it is comparable to or lower than the recombination rate, there will be a significant drop in temperature as the gas recombines. Furthermore, the long term temperatures will decrease, due to the ionisation rate falling below the recombination rate for a fully ionised state; this causes the ionisation equilibrium state to shift to a more neutral state, lowering the rate of ionisations and recombinations, and therefore the photo-heating rate near equilibrium, as the cooling rate is now limited by ionisations, rather than recombinations.

In our simulation, there is enough ionising flux with the fiducial values from Meraxes to maintain a highly ionised state in the vast majority of the simulated volume, so a significant drop in temperature is observed only for highly dense regions. Voxels at mean density will only show a 5% decrease in temperatures when the ionising flux is decreased by a factor of 10. Differences in temperature will also diminish over time, as the gas cools to its thermal asymptote.

Regions near early galaxies that have their star formation suppressed or move between voxels can also show this behaviour, however in this case the temperature drop is also usually very short-lived, as nearby HII bubbles expand to re-heat the voxel, washing out any memory of previous ionisation events.