Inhomogeneous Reionization

Inhomogeneous Reionization Models in Cosmological Hydrodynamical Simulations

Jose Oñorbe, F. B. Davies, Z. Lukić, J. F. Hennawi and D. Sorini
Institute for Astronomy, University of Edinburgh, Blackford Hill, EH9 3HJ, Edinburgh, United Kingdom
Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
Department of Physics, University of California, Santa Barbara, CA 93106-9530, USA
Lawrence Berkeley National Laboratory, CA 94720-8139, USA
Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD)
E-mail: onorbe@roe.ac.uk
Accepted XXX. Received YYY; in original form ZZZ
Abstract

In this work we present a new hybrid method to simulate the thermal effects of the reionization in cosmological hydrodynamical simulations. The method improves upon the standard approach used in simulations of the intergalactic medium (IGM) and galaxy formation without a significant increase of the computational cost allowing for efficient exploration of the parameter space. The method uses a small set of phenomenological input parameters and combines a semi-numerical reionization model to solve for the topology of reionization and an approximate model of how reionization heats the IGM, with the massively parallel Nyx hydrodynamics code, specifically designed to solve for the structure of diffuse IGM gas. We have produced several large-scale high resolution cosmological hydrodynamical simulations (, Mpc/h) with different instantaneous and inhomogeneous reionization models that use this new methodology. We study the IGM thermal properties of these models and find that large scale temperature fluctuations extend well beyond the end of reionization. Analyzing the 1D flux power spectrum of these models, we find up to differences in the large scale properties (low modes, s/km) of the post-reionization power spectrum due to the thermal fluctuations. We show that these differences could allow one to distinguish between different reionization scenarios already with existing Ly forest measurements. Finally, we explore the differences in the small-scale cutoff of the power spectrum and we find that, for the same heat input, models show very good agreement provided that the reionization redshift of the instantaneous reionization model happens at the midpoint of the inhomogeneous model.

keywords:
methods: numerical, early universe, intergalactic medium, large-scale structure of universe, reionization
pubyear: 2018pagerange: Inhomogeneous Reionization Models in Cosmological Hydrodynamical SimulationsA

1 Introduction

The reionization of neutral hydrogen () in the Universe is believed to be driven by the ultraviolet radiation emitted by the first galaxies and quasars. Our current knowledge about reionization is mainly informed by two observational probes: an integral constraint from the electron scattering optical depth of the cosmic microwave background (CMB) and the measured transmission of the Lyman- (Ly) forest (McQuinn, 2016).

The most recent Planck-satellite measurements of the CMB Thomson scattering optical depth, , imply that hydrogen was reionized by ( limits from Planck Collaboration et al., 2018). However, complete Gunn-Peterson absorption (Gunn & Peterson, 1965) observed in the spectra of many of the highest quasars, along with the steep rise of both the Ly optical depth and its sightline-to-sightline scatter with redshift, has led to the consensus that we are witnessing the end of reionization only at (Fan et al., 2006; Becker et al., 2015; Bosman et al., 2018; Eilers et al., 2018). At higher redshifts , the Ly opacity can only set lower limits on the redshift of reionization, but measurements of the damping wing of in the intergalactic medium (IGM) in the spectra of quasars have recently started to provide interesting constraints of the hydrogen neutral fraction at these redshifts consistent with CMB observations (Mortlock et al., 2011; Greig et al., 2017; Davies et al., 2018b; Bañados et al., 2018). Therefore the exact details about hydrogen reionization, e.g. main sources responsible of the ionizing photons, are still not known and in fact several large observational projects within the astrophysics community are underway, with the goal of finally fully understanding this fundamental cosmological process (e.g. JWST, SKA, HERA).

A key aspect of the reionization of is that it plays a crucial role in the thermal history of the Universe. During reionization, ionization fronts around the first sources propagate supersonically and the IGM gas that has been cooling since the Big Bang is photoheated. Studies have suggested a post-reionization temperature of order Kelvin depending on the spectral shape of ionizing sources and the speed of ionization fronts (Miralda-Escudé & Rees, 1994; Abel & Haehnelt, 1999; Tittley & Meiksin, 2007; McQuinn, 2012; Davies et al., 2016; Park et al., 2016; Finlator et al., 2018; D’Aloisio et al., 2018a). After reionization, hydrodynamical simulations show that the IGM expands and cools, but because cooling times in the low-density IGM are long, the thermal memory of this heat injection can persist for hundreds of Myr () (Trac et al., 2008; Furlanetto & Oh, 2009). Thus if the temperature of the IGM can be measured just after reionization at , before the thermal state has fully relaxed, the nature and timing of reionization can be constrained (Cen et al., 2009; Lidz & Malloy, 2014; Nasir et al., 2016; Oñorbe et al., 2017a, b; Boera et al., 2018). However the temperature of the IGM is not the entire story. Although baryons in the IGM trace dark matter fluctuations on Mpc scales, on smaller scales ( kpc) the gas is pressure supported against gravitational collapse by its finite temperature ( K, Gnedin & Hui, 1998a; Kulkarni et al., 2015; Oñorbe et al., 2017a, b). Analogous to the classic Jeans argument, baryonic fluctuations are suppressed relative to the pressureless dark matter (which can collapse). The pressure smoothing scale of the IGM thus provides an integrated record of the thermal history of the IGM, and is sensitive to the timing of and heat injection by reionization events (Gnedin & Hui, 1998a; Kulkarni et al., 2015; Nasir et al., 2016; Oñorbe et al., 2017a, b).

Recent measurements of the Ly optical depth at high redshift have found enhanced scatter at that exceeds what can be attributed to density fluctuations alone (Fan et al., 2006; Becker et al., 2015; Bosman et al., 2018; Eilers et al., 2018). However the origin of these fluctuations has been hotly debated. It has been argued that they are driven by fluctuations in the radiation field (Davies & Furlanetto, 2016; Gnedin et al., 2017; D’Aloisio et al., 2018b) or the temperature field (D’Aloisio et al., 2015; Davies et al., 2018a), both of which may result from a patchy, extended, and late-ending reionization process. Still others have interpreted these fluctuations as evidence that reionization was driven by rare AGN (Chardin et al. 2015, 2017, but see D’Aloisio et al. 2017), rather than the faintest galaxies (Robertson et al., 2015), intriguing given that an abundant population of faint AGN may have just been discovered (Giallongo et al., 2015, see however Weigel et al. 2015).

Thus, previous work has demonstrated that there is not a sharp transition from the epoch of reionization to the simple uniform UV background and tight temperature-density relation that govern the Ly forest at lower redshift . Instead, spatially coherent radiation and temperature fluctuations are expected to persist even after reionization is complete, leaving observable artifacts in the post-reionization IGM. As a result, study of the statistical properties of the IGM, where the transmitted flux is still significant, can address fundamental questions about reionization. In fact, current Ly forest measurements at (Viel et al., 2013a; Iršič et al., 2017; D’Aloisio et al., 2018b) provide constraints on the power at the relevant scales ( s/km, h/Mpc) and therefore may already tell us something about the effects of inhomogeneous reionization on the IGM at these redshifts. Moreover the past several years have seen a dramatic fivefold increase in the number of quasars from deep wide-field optical/IR surveys (Willott et al., 2010; Bañados et al., 2014; Venemans et al., 2015b, a; Wang et al., 2017). The ongoing spectroscopic follow-up of these objects (e.g. Eilers et al., 2018) ensures a significant precision boost of Ly forest measurements at these redshifts in the near future.

Despite the apparent simplicity, accurate simulations of the IGM are hard because of two spatial scales. High spatial resolution ( kpc) is required to resolve the density structure of the IGM, while a relatively large box size ( Mpc) is required both to capture the large-scale power and to obtain a fair sample of the universe to capture the reionization fluctuation scale (see e.g. Iliev et al., 2014; Lukić et al., 2015; Oñorbe et al., 2017b; Bolton et al., 2017). Ideally one would run coupled radiative transfer hydrodynamical simulations that include extra physics governing the sources of ionizing photons (stars, quasars, etc.). Although significant progress on this front has been made (e.g. Meiksin & Tittley, 2012; Wise et al., 2014; So et al., 2014; Gnedin, 2014; Pawlik et al., 2015; Norman et al., 2015; Ocvirk et al., 2016; Kaurov, 2016; La Plante et al., 2017), these simulations are still too costly for exploration of the parameter space. Post-processing radiative transfer approaches have also been used in the field (e.g. Trac & Cen, 2007; Iliev et al., 2014; Bauer et al., 2015; Chardin et al., 2017; Keating et al., 2018; Kulkarni et al., 2018) as they provide a cheaper alternative and a coarser resolution can be used to solve the radiative transfer, however they do not capture the full hydrodynamical evolution. For this reason, semi-numerical and analytical methods of reionization modeling continue to remain attractive for efficient and flexible exploration of the parameter space (e.g. Furlanetto et al., 2004; Mesinger et al., 2011; Alvarez et al., 2012; Kaurov, 2016; Kim et al., 2016; Hassan et al., 2016).

The dominant approach, implemented in the vast majority of hydrodynamical codes that aim for moderate to high spatial resolution, is to assume that all gas elements are optically thin to ionizing photons. In this way their ionization state can be fully described by a uniform and isotropic UV+X-ray background radiation field (e.g. Katz et al., 1996). Thus, the radiation field is encapsulated by a set of photoionization and photoheating rates that evolve with redshift for each relevant ion (Haardt & Madau, 1996, 2001; Faucher-Giguère et al., 2009; Haardt & Madau, 2012). However, these UVB synthesis models surely break down during reionization events, and therefore the inhomogeneous nature of these processes is not captured (see discussion e.g. Oñorbe et al., 2017a; Puchwein et al., 2018).

The necessary next step to improve the predictions of numerical models is therefore to drop the assumption of a uniform and isotropic background in self consistent hydrodynamical simulations. Motivated by this, we present in this paper a hybrid scheme for implementing inhomogeneous reionization physics into large scale cosmological hydrodynamical simulations based on a small set of phenomenological input parameters. These models combine a semi-numerical reionization model to solve for the morphology of reionization (Mesinger et al., 2011) and an approximate model of how reionization heats the IGM, with the massively parallel Nyx hydrodynamics code (Almgren et al., 2013), specifically designed to solve for the structure of diffuse IGM gas.

The structure of this paper is as follows. In Section 2 we describe our new method to model reionization in hydrodynamical simulations and compare it with other methods used in the literature. Section 3 describes the characteristics of all the hydrodynamical simulations studied in this work. Section 4 presents the results of simulations using one inhomogeneous and instantaneous reionization models but different heat input during reionization. We compare the outcome of different inhomogeneous reionization models in Section 5. In Section 6 we further investigate the differences between our simulations studying the properties of their temperature and density fields. Section 7 shows the effects that fluctuations of the UV background could have on the Ly forest in both our flash and inhomogeneous reionization models. We discuss in Section 8 the limitations of our new approach and compare it with previous work. We conclude by presenting a summary of our results in Section 9. Finally in Appendix A we present a set of tests in which we run different variations of the new reionization modelling method introduced in this paper.

2 Modeling Reionization in Hydrodynamical Simulations

In this paper we will focus on non-radiative transfer methods to implement reionization and the UVB in hydrodynamical simulations that are suited to create a large grid of different models and/or easily used in combination with expensive subgrid star formation and feedback prescriptions. In this context we define three different approaches depending on how reionization is modeled in the simulations: homogeneous reionization, inhomogeneous reionization, and flash reionization.

2.1 Homogeneous Reionization

A standard approach to model reionization and the UV background in hydrodynamical simulations is to assume that all resolution elements are optically thin to radiation. Thus, radiative feedback is accounted for via a spatially uniform, time-varying ultraviolet background (UVB) radiation field, input to the code as a list of photoionization and photoheating rates that vary with redshift (e.g. Katz et al., 1996). These rates are obtained from one-dimensional cosmological radiative transfer calculations (e.g. Faucher-Giguère et al., 2009; Haardt & Madau, 2012; Khaire & Srianand, 2018; Puchwein et al., 2018), also known as UVB synthesis models, and give the evolution of both the photoionization, and the photoheating rates with redshift for, at least, , and . This method is used in several hydrodynamical simulations that study galaxy formation and evolution (e.g. Schaye et al., 2015; Pillepich et al., 2018; Hopkins et al., 2018) and the circumgalactic medium (CGM) or the IGM (e.g. Lukić et al., 2015; Bolton et al., 2017; Oñorbe et al., 2017b). However, these UVB synthesis models surely break down during inhomogeneous reionization events and therefore their use in optically thin simulations has some intrinsic limitations. One of the most important issues is that the heat injection in this approach does not behave like realistic ionization fronts and it is done too gradually and therefore is usually underestimated (see discussion e.g. Oñorbe et al., 2017a; Puchwein et al., 2018).

2.2 Inhomogeneous Reionization

We introduce here a new method to account for inhomogeneous reionization in cosmological hydrodynamical simulations. A similar but more simple approach was introduced by Feng et al. (2016) based on semi-numerical work by Battaglia et al. (2013). This method assigns a specific reionization redshift to each resolution element of the simulation and works as follows. When the reionization redshift of a resolution element is smaller than the current redshift of the simulation, the UV background seen by the resolution element is assumed to be zero. During the first time step of the simulation in which the reionization redshift is crossed we inject a fixed amount of energy to account for the heating produced by the ionization front (Abel & Haehnelt, 1999). We will first describe how we have implemented this heating in the simulation, and below we will discuss how we have generated the different reionization models.

The heating due to reionization is a complicated process that not only depends on the shape of the ionizing spectrum but also on the local density field and how fast the ionization front travels (Miralda-Escudé & Rees, 1994; Abel & Haehnelt, 1999; McQuinn, 2012; Davies et al., 2016; D’Aloisio et al., 2018a). These studies have suggested a post-reionization temperature of order K, depending on the ionizing photon spectrum and the velocity of the ionization front111Results from 3D radiative transfer cosmological simulations are still not that homogeneous in this regard. The treatment of multifrequency radiative transfer and particularly of the hard tail of the photon spectrum varies significantly among the methods and yields a correspondingly large range of temperatures after the reionization front (see e.g. Iliev et al., 2009).. In this work we choose as a first step to simply parametrize our ignorance of the details of reionization using heat injection as a free parameter, . The temperature of the resolution element after reionization, , is then calculated as follows:

(1)

where and are the temperature and the hydrogen ionization state just before reionization and is again a free parameter in our code. In this way, represents the maximum temperature that can be reached after reionization222We discuss in Appendix A different variations of this approach.. If the resolution element was already ionized the heating will be attenuated proportionally. This avoids any heat injection to regions that may have been already ionized by other physical processes included in the simulations, e.g., collisional ionization. Finally, we want to emphasize that we inject the necessary energy to ensure that the final temperature of the resolution element will be taking into account its final ionization state, with the UVB already turned on.

While the density dependence of the heat injection during reionization may be the easiest further improvement that could be added to the model (see e.g. Kaurov, 2016; Hirata, 2018), its exact solution still depends on the assumed velocity of the ionization front. Therefore we have decided that in order to start exploring the parameter space, the best approach is to use the simplest model and assume that the reionization temperature is independent of density, which is common in studies of the thermal history of the IGM. We will discuss in detail about the implications of these results for our methodology in Section 8.

After its reionization redshift, the resolution element is assumed to be optically thin and is exposed to the (uniform) tabulated photoionization rates given to the code. However, in contrast with the full homogeneous approach typically used in cosmological hydrodynamical simulations, in this method the UV background is only applied to the resolution elements that have already been reionized. Therefore, we need to provide the evolution of the average UV background in ionized regions and not just the evolution of the average UV background in the Universe. Of course once reionization is finished both quantities will be the same but during reionization they will be very different (see e.g. Lidz et al., 2007). For this reason we have used the photoionization and photoheating rates of Haardt & Madau (2012)333with the minor corrections suggested by (Oñorbe et al., 2017b) to provide a better fit to the observed and mean flux evolution for . At higher redshifts, , we have assumed that the photoionization and photoheating in ionized regions are the same as at which guarantees that all ionized regions have a ionized fraction of . While this approach does not directly include UVB fluctuations they can be modeled straightforwardly in post-processing. We will present results of adding UVB fluctuations in this way to our simulations in Section 7

We now describe how we assign a reionization redshift to each resolution element. We use a semi-numerical method based on the excursion set formalism, applied to the initial conditions of the simulation, to pre-compute a reionization redshift field that is read at the beginning of the simulation. We describe in detail in §2.2.1 how we have done this. Alternatively, the reionization redshift could be computed while the simulation is running using the halo mass function computed directly on the fly, or some more advanced metrics if star formation is also modeled in the code (see e.g. Poole et al., 2016). We are already exploring this possibility and we will present these results in the near future. In any case, we want to emphasize that the heat injection method described above is independent of how the reionization redshift is assigned.

2.2.1 Using Semi-numerical Methods to Generate Reionization Fields

Figure 1: Left panel: A kpc/h slice of the reionization redshift field generated using the method described in Section 2.2.1 (model A). The evolution of the ionization fraction for this model can be found in Figure 3. The exact parameters used to generate this field are summarized in Table 1. Middle panel: The same slice as in the left panel, now showing the total density field at . Right panel: The same slice presented in the top left panel this time showing the total density field at to which we have applied a gaussian smoothed filter of 1 Mpc/h.
Figure 2: Left panel: The 2D histogram of model A reionization redshift and the density fields of a cosmological hydrodynamical simulation using this model. No strong correlation between the reionization redshift of one resolution element and its precise density value is found. Middle and right panel: The 2D histogram of model A reionization redshift and the smoothed density field with a gaussian smoothed filter of 1 Mpc/h (middle) and 5 Mpc/h (right). This shows that the reionization redshift of each resolution element is correlated with the density field smoothed on large scales, as expected for reionization fields generated using the excursion set formalism. See text for more details.

Our method relies on having a unique reionization redshift for each resolution element in the simulation. In this paper we use a semi-numerical approach to generate the reionization redshift field based on the analytical method introduced by Furlanetto et al. (2004) to describe the morphology of reionization on large scales. In this model each point at position has assigned a one-dimensional function, trajectory , corresponding to the average matter density within spheres of radius centered in . The physical motivation is that the function can provide an estimate for the fraction of matter collapsed into halos (Press & Schechter, 1974). Then, assuming a rate of ionizing photon production in halos, one can derive the total number of ionizing photons in the spherical regions centered at the point of interest. Therefore as the density field is evolved a ionization history can be obtained.

This approach is therefore physically based on the excursion set framework which is widely used in structure formation theory (Bond et al., 1991). It has become very popular in recent years because of its high computational efficiency. Various improvements can be made to this model, which allow to include more sophisticated physics (e.g. Sobacchi & Mesinger, 2014). We have used the publicly available code 21cmFAST (Mesinger et al., 2011) which is fundamentally based on the idea described above. There, the fraction of material in collapsed objects is computed analytically from an initial density field, which in our case corresponds to a particular realization of the high- initial conditions of the hydrodynamical simulation. To generate the reionization models used in this paper we have just freely vary two parameters: the minimum halo mass for producing ionizing photons, , and the ionizing efficiency, . We adopt the central-pixel-flagging algorithm instead of the slower bubble-painting algorithm (see Zahn et al., 2011, for a discussion of the relative accuracy of these methods) and a k-space top-hat filter when smoothing the density field. In this work we initialize 21cmFAST using the same grid as in the Nyx simulations, 2048, but generate evolved density and ionization fields on a coarser grid of 256 (models A, B, D) or 128 (model C). We generate ionization field using temporal resolution of going from to . From the combination of these ionization fields we construct a reionization redshift field that can be plug in to the hydrodynamical simulation.

The left panel of Figure 1 shows a random slice of a reionization redshift field generated using the method described above. To generate this model we have used the total density field of a CDM cosmological box of size Mpc/, and . The exact values of these parameters are not relevant here and we will discuss them in detail below (See Section 3). This model produces a reionization history consistent with current CMB constraints (Planck Collaboration et al., 2018). The middle panel presents the same slice but now showing the total density field at obtained running a cosmological hydrodynamical simulation. The main cosmological structures in the slice can be identified by eye and one can see that as expected there is a correlation between dense structures and early values of the reionization redshift seen in the left panel. However the right panel shows a gaussian smoothed version of the original density field at with a filter of Mpc/h. In this case, the correlation between the density values and the reionization map is even more clear.

The left panel of Figure 2 demonstrates that the relation between density and reionization redshift is not simple. It displays the 2D histogram of the reionization redshift and the density field. It is clear that there is no strong correlation between the reionization redshift of one resolution element and its precise density value at (or any other redshift). The middle and right panel of Figure 2 show the same slice again but this time for a gaussian smoothed version of the total density field at with a filter of Mpc/h and Mpc/h (respectively). In this case there is a clear, although noisy, correlation between the reionization redshift of a resolution element and its specific smoothed density value. A similar correlation is found regarding of the redshift at which we pick the density field. This shows that the reionization redshift of each resolution element is correlated with the density field smoothed on large scales, as expected for reionization fields generated using the excursion set formalism (see e.g. Battaglia et al., 2013).

In this work we have created four different inhomogeneous reionization models using 21CMFAST for various minimum halo masses: , and (models A, B and C respectively). We adjusted the ionization efficiency in these models so that all of them had the median reionization redshift as close as possible to while still having reionization finish at consistent with the Ly forest. We also generated an extra model with with a faster reionization (model D, cyan solid line) using an ionizing efficiency evolving with redshift, . The upper panel of Figure 3 shows the ionization fraction evolution for each of these models (green, pink and orange solid lines). The lower panel shows the evolution of the CMB electron scattering optical depth444We used the results of the hydrodynamical simulations to computed the density weighted ionization history for this model, which is the relevant ionization history for computing (see e.g. Liu et al., 2016), , for each of these models. The gray band shows the most recent constraints on from CMB data (Planck Collaboration et al., 2018). Table 1 contains a list of different parameters characterizing these reionization models: the median reionization redshift, , the width of reionization, , Thompson scattering optical depth, and the end of reionization redshift, .

Figure 3: Reionization models studied in this work. Upper panel: evolution of the ionization fraction for the different reionization models considered in this work. Notice that the models were built to have a similar median reionization redshift, . Lower panel: integrated electron scattering optical depth, , computed from the above models using the density weighted ionization fraction in the simulation.. The gray band stand for last constrains on coming from Planck Collaboration et al. (2018) data.

2.3 Flash Reionization

We consider one more method, “flash”, to simulate the reionization process in hydrodynamical simulations which relies on the same methodology to inject heat, but has fixed reionization redshift for all resolution elements in the simulation. Therefore these simulations will have only two parameters to describe the reionization process: the redshift of reionization, , and the heat injection, .

Of course this method misses by construction the inhomogeneities associated with reionization. However, it allows us to easily compare them with the inhomogeneous approach and study in more detail the effect that these inhomogeneities can have on different scales. We show in Figure 3 a flash reionization evolution model that happens at (black line). Note that this value is very close the median reionization redshift of the four inhomogeneous runs, .

3 Simulations

Sim method
(K) ()
IR-A Inhomogeneous (Model A) ()
IR-Acold Inhomogeneous (Model A) ()
IR-Ahot Inhomogeneous (Model A) ()
IR-B Inhomogeneous (Model B) ()
IR-C Inhomogeneous (Model C) ()
IR-D Inhomogeneous (Model D) () evol.
FR Flash
FRcold Flash
FRhot Flash
Column 1: Simulation code.
Column 2: Method used to simulate reionization. See text for more details.
Column 3: reionization redshift.
Column 4: Width of reionization. .
Column 5: Thompson scattering optical depth density weighted (volume weighted).
Column 6: Total heating assumed for reionization.
Column 7: Parameter excursion set model 1. Ionizing efficiency.
Column 8: Parameter excursion set model 2. Minimum mass.
All simulations have a box size of length, and resolution elements.
Table 1: Summary of Simulations.

The simulations used in this work were performed with the Nyx code (Almgren et al., 2013). Nyx follows the evolution of dark matter simulated as self-gravitating Lagrangian particles, and baryons modeled as an ideal gas on a uniform Cartesian grid. The Eulerian gas dynamics equations are solved using a second-order accurate piecewise parabolic method (PPM) to accurately capture shock waves. We do not make use of Adaptive Mesh Refinement (AMR) capabilities of Nyx in the current work, as the Ly forest signal spans nearly the entire simulation domain rather than isolated concentrations of matter where AMR is more effective. For more details of these numerical methods and scaling behavior tests, see Almgren et al. (2013).

Besides solving for gravity and the Euler equations, we also include the main physical processes fundamental to model the Ly forest. First we consider the chemistry of the gas as having a primordial composition with hydrogen and helium mass abundances of , and , respectively. In addition, we include inverse Compton cooling off the microwave background and keep track of the net loss of thermal energy resulting from atomic collisional processes. We used the updated recombination, collision ionization, dielectric recombination rates, and cooling rates given in Lukić et al. (2015).

In order to generate the initial conditions, we have used the music code (Hahn & Abel, 2011), and the camb code (Lewis et al., 2000; Howlett et al., 2012) to generate the transfer function. All simulations started at to be sure that non-linear evolution is not compromised. Throughout this paper we assumed a CDM cosmology with the following fundamental parameters: , , , , and . These values are within one sigma agreement with the latest cosmological constrains from the CMB (Planck Collaboration et al., 2018). The choice of hydrogen and helium mass abundances ( and , therefore ) is in agreement with the recent CMB observations and Big Bang nucleosynthesis (Coc et al., 2013). All simulations were run down to , saving 32 snapshots555For all simulations we saved snapshots at the following redshifts: , , , , , , , , , , , , , , , , , , , , , , , , , , . Unless otherwise stated all simulations presented here have a box size of length, and resolution elements.

We have run hydrodynamical simulations with the four inhomogeneous reionization models (IR-A, IR-B, IR-C, IR-D named after the reionization models used in them, see Figure 3), and the flash reionization model (FR) described in Section 2.2. In all these runs we assumed a total heat injection of K which is the standard value obtained in galaxy driven reionization models (see above). Quasar or Population III driven scenarios may inject more heat, as much as K (but see D’Aloisio et al., 2018a, and discussion in Section 8). So in order to study the effects of different heat injection during reionization we have also run the flash reionization model, , and model A inhomogeneous reionization ( ) with two different values for the heat injected during reionization: a colder model with K and a hotter model with K (hereafter referred as IR-Acold and IR-Ahot for the inhomogeneous runs and FR-cold and FR-hot for the flash runs respectively).

4 Comparison of Flash and Inhomogeneous Reionization Simulations

We now present the results of the inhomogeneous and flash reionization simulations. We will first focus on the differences between these two methods. Figure 4 shows the same slice of the temperature field of our fiducial flash reionization simulation (FR, left column) and inhomogeneous reionization simulation (IR-A, right column) at , , and (from top to bottom respectively). These two simulations have the same heat injection, , and the reionization redshift of the flash simulations is equal to the median reionization redshift of the inhomogeneous reionization, . Therefore at reionization has just finished in the flash reionization and we can see an overall homogeneous high temperature (top left panel) while at the same redshift the inhomogeneous reionization (top right panel) is still ongoing and we can clearly see the difference in temperatures between regions that have been reionized and those which remain neutral. At the temperature field in the inhomogeneous approach (bottom middle panel) shows lingering signatures of the temperature fluctuations on large scales ( Mpc/h). In this run different regions of the Universe are reionized and heated at different times, and so they asymptote to the temperature set by the balance between photoheating and the adiabatic and Compton cooling governing the temperature-density relation at different times (McQuinn & Upton Sanderbeck, 2016). In this particular case we see that the high density regions which reionized at high redshift have had time to cool and show lower temperatures than the low density regions, which were reionized much later (D’Aloisio et al., 2015; Davies et al., 2018a). This will have relevant implications in the properties of the Ly at these large scales that we will discuss in detail in Section 4.1. The high temperature filamentary structure observed in all panels of Figure 4 reflects the high density regions where collapsed objects have heated the gas to their virial temperatures. This is of course independent of the reionization method employed666Note that, contrary to some other models as, e.g., Trac et al. (2008), our simulations do not include any kind of thermal feedback in the dense regions from star formation..

It is also interesting to study the differences in the temperature-density relation between these two runs. Figure 5 shows the volume weighted temperature-density histogram for the same flash (FR, upper row) and inhomogeneous (IR-A, lower row) simulations at , , and (left, middle and right columns). At reionization has just happened in the flash simulations () and therefore most of the gas is still at high temperatures, leading to a flat temperature-density relation. As time evolves the adiabatic and Compton cooling is efficient enough to give rise to a steeper temperature-density relation despite the photoheating background. For the inhomogeneous run however the temperature-density relation is not well defined until reionization is finished ( in this model). Before this time, the temperature-density distribution is bimodal, reflecting regions that have experienced reionization versus the regions that still have not. This is clearly seen in the temperature-density distribution (lower left panel of Figure 5). Even when reionization is finished the inhomogeneous simulation shows a much larger scatter than the flash simulation, especially at lower densities (see Trac et al., 2008; Lidz & Malloy, 2014; Keating et al., 2018, for similar findings). This is produced because, as shown in Section 2, for a fixed density there is a wide range of reionization redshifts and we can see the effect of the different timing of the heating and cooling process for all resolution elements. Towards lower redshift the scatter is reduced, but even at we still see a significant difference between the inhomogeneous run and the flash run (see also Trac et al., 2008, who found similar effects using radiative transfer simulations).

We have also measured at each snapshot the volume weighted temperature at mean density, and the slope of the power-law temperature-density relation, , where . The parameters are obtained by fitting the temperature-density relation with linear least squares in and using only cells that satisfy the following criteria: and . The panels in Figure 5 show as a black solid line the results for both simulations. Changing these thresholds for the flash simulation within reasonable IGM densities () and temperatures () produce differences just at a few per cent level at any redshift after reionization. We have tested that we obtain similar conclusions if we change the fitting method and obtain the thermal parameters finding the median of the gas temperature in bins of width dex at and . We find no relevant changes if we modify these two densities slightly (again as long as we do not include too high densities, ). We find however that the fitting procedure used seems to have a more important role in the value of the thermal parameters in the inhomogeneous runs. This is not only true before reionization is complete but also at later times. At our inhomogeneous reionization is not yet complete so one could argue that in order to define a clear temperature-density relation only ionized gas should be used (top right panel of Figure 5). Once reionization is complete we find that different fitting methods and parameters result in differences of and at and differences at . Notice however that these differences are systematics for all simulations. We find that the parameters with the largest effect on the temperature-density relation fit are the density thresholds. Although this effect has no relevant consequences for the results and conclusions presented in this work it will have to be considered in future studies of the IGM thermal parameters that use simulations of inhomogeneous reionization.

Figure 4: The same slice of the temperature field for our fiducial flash (FR, left column) and inhomogeneous (IR-A, right column) reionization simulations at redshifts , , and (from top to bottom respectively). Thermal fluctuations are clearly seen in the inhomogeneous run well after reionization is finished, with hotter (colder) regions tracing late (early) reionization regions.

In Figure 6 we present the thermal histories for a set of flash (dotted lines) and inhomogeneous reionization simulations (solid lines). All flash simulations share the same reionization redshift but they differ in the heat injected during reionization: K, K, K (blue, green and red respectively). We show the analogous inhomogeneous reionization simulations with different heat injection during reionization and the same reionization evolution in which . For the inhomogeneous runs we only show the evolution once reionization is finished () as before that time such global parameters are not very meaningful. The first and second panels of Figure 6 show the evolution of the and thermal parameters (respectively) governing the density-temperature relation determined by fitting the distribution of densities and temperatures in the simulation discussed above. The color bands shown for each inhomogeneous simulation illustrate the maximum variation that we found in these parameters changing the fitting method. As a reference we also show in the first panel the median value of temperature (black dots) and the scatter (black solid errorbars) that we found for the IR-A inhomogeneous model at , and .

The bottom panel of Figure 6 shows the evolution of the pressure smoothing scale, , as defined by Kulkarni et al. (2015). These authors define a pseudo real-space Ly flux field, which is the computed similarly to the true Ly forest flux but without redshift space effects such as peculiar velocities and thermal Doppler broadening. This field naturally suppresses the dense gas that would otherwise dominate the baryon power spectrum, making it robust against the poorly understood physics of galaxy formation and revealing the pressure smoothing in the diffuse IGM.

Figure 5: The temperature-density relation for our fiducial flash (FR, left column) and inhomogeneous (IR-A, right column) reionization simulations at redshifts , , and (from top to bottom respectively). At reionization has just happened in the flash run () while the inhomogeneous run still shows a bimodal distribution. Even after reionization is finished in the inhomogeneous run ()a much wider temperature-density distribution is found in the inhomogeneous simulation.

Inspection of Figure 6 shows that, regardless of the reionization method implemented in the simulation, models with more heat injection give rise to larger temperatures and slightly lower values. These models also produce larger pressure smoothing scales, and these differences persist even at lower redshifts long after the other thermal parameters, and , have relaxed to their asymptotic values. Once reionization is finished all these models share the same photoionization and photoheating rates and therefore and asymptote to the same values much faster than the pressure smoothing scale, which retains memory of the thermal history (Kulkarni et al., 2015; Nasir et al., 2016; Oñorbe et al., 2017b). While is very similar between the analogous flash and inhomogeneous models, this is not the case for and inhomogeneous simulations show systematically lower values. This is expected from the study of the temperature-density relation presented above (Figure 5) as we have seen that the slope of the density-temperature relation is more sensitive to the larger scatter in this relation that appears in inhomogeneous reionization models. However the evolution of the pressure smoothing scale, , unlike the other thermal parameters, is very different between the flash and inhomogeneous models that share the same amount of heat injection during reionization. In the next section we consider how the differences in thermal parameters between these simulations translate into Ly statistics. We return to a more detailed discussion of the statistical properties of the temperature and density field in Section 6.

4.1 Lyman- Forest Statistics: The 1D Power Spectrum

In order to further explore the differences between the inhomogeneous and flash reionization simulations presented above we now investigate the properties of their respective Ly forests. To compute the Ly forest spectra from the simulation, we compute the optical depth at a fixed redshift, which can then be easily converted into a transmitted flux fraction, . That is, we do not extract outputs along the light cone when we cast rays in the simulation; we use the gas state at a single cosmic time (i.e. a single hydrodynamical simulation snapshot). The simulated spectra are not meant to look like observational ones, but simply recover the statistics of the flux in a small redshift window. Our calculation of the spectra accounts for Doppler shifts due to bulk flows of the gas, as well as for thermal broadening of the Ly line. We refer to Lukić et al. (2015) for specific details of these calculations. This procedure results in the Ly flux as a function of wavelength or equivalently time or distance. Following the standard approach, we then rescale the UV background intensity so that the mean flux of all the extracted spectra from the simulation matches the observed mean flux at the respective redshift We have neglected noise and metal contamination in our skewers so far, but this will not be relevant for the results presented in this paper.

Figure 6: The evolution of thermal parameters with redshift for a set of flash (dashed lines) and inhomogeneous models (solid lines). Models with the same color have the same heat injection during reionization: K (blue), K (green) and K (red). From top to bottom it shows the temperature at mean density, , the slope of the density-temperature relation, and the pressure smoothing scale, , as defined in Kulkarni et al. (2015). In the first panel, the black circles and error bars show the median value of the temperature at mean density and the scatter found for the IR-A run at , and . The color bands illustrate the scatter in an that we find using slightly different fitting methods for the inhomogeneous models. See text for more details.

We have calculated the 1D flux power spectrum, , for each simulation at different redshifts which is sensitive to the parameters governing the thermal state of the IGM. Pressure smoothing, , damps out small-scale fluctuations in the gas, while random thermal motions (sensitive to temperature, or and ) Doppler broadens Ly forest lines, further reducing the amount of small-scale structure. Both of these effects combine to produce a prominent small-scale (high-k) cutoff in the flux power spectrum (Zaldarriaga et al., 2001; Peeples et al., 2010; Walther et al., 2018a). Therefore, by carefully studying this cutoff it is possible to constrain not only the thermal state of the gas but also its full thermal history (Nasir et al., 2016; Oñorbe et al., 2017b; Boera et al., 2018). The large scale properties of Ly forest (low-k) can provide information about possible fluctuations of the UVB or the temperature field that could provide very relevant clues about possible sources responsible for reionization (e.g. D’Aloisio et al., 2018b). We have computed the power spectrum, , of the fractional contrast, , at each redshift defined as =. A total of skewers with a length equal to the size of the box have been used at each redshift. We computed the power spectrum of each skewer and then calculated the average value at each mode, .

In this work we assumed the following mean flux values at each redshift to be consistent (within ) with recent observational constraints. This level of accuracy is more than enough for this work and none of the results presented in this paper depend on the exact values assumed for the mean flux. At , and we have used and respectively. This is in good agreement (within ) with observations (Becker & Bolton, 2013; Viel et al., 2013a; Eilers et al., 2018). At and we have used and respectively consistent with recent measurements (D’Aloisio et al., 2018b; Bosman et al., 2018) and just above the limit of Eilers et al. (2018). Finally, for we have assumed a mean flux of . This is just slightly above the range measured by Bosman et al. (2018, ) and Eilers et al. (2018, ).

Figure 7 presents the dimensionless 1D flux power spectrum, , computed at , , and (left, middle and right panels respectively) for the different flash (dashed lines; FR, FRcold, FRhot) and inhomogeneous (solid lines; IR-A, IR-Acold, IR-Ahot) simulations. The color of each line indicates the heat injected during reionization: K (blue), K (green), and K (red). The first thing to notice between the three panels is that the overall power level increases with redshift, which reflects the fact that as the average mean flux decreases toward higher , density fluctuations are exponentially amplified (e.g. Viel et al., 2004; Palanque-Delabrouille et al., 2013; Viel et al., 2013b; Walther et al., 2018a). At each redshift the lower panels show the percentage difference between each flash simulation and its analogous inhomogeneous run i.e., the one that has the same heat input during reionization, . At low- modes (large scales) there is a systematic difference between the reionization approaches and the inhomogeneous runs show more power at these scales than their flash counterparts. This is due to the large-scale coherence of temperature fluctuations which can be linked to the specific details of the reionization model and the sources responsible of it. We can see that both the overall difference in power and the mode at which this effect is relevant decrease with redshift. This decrease shows how the temperature fluctuations are being attenuated in the simulation as time evolves. As expected, inhomogeneous models with a larger total heat injection during reionization show larger differences with the flash reionization model because the temperature fluctuations in these models are also larger. Note that the difference between the flash and inhomogeneous models gradually increases towards lower modes. At high-, , the difference is as large as for the lowest modes that we can study in these simulations ( s/km). We expect these differences to continue increasing as we go to lower modes (e.g. D’Aloisio et al., 2018b). However, note that the overall power decreases as we go to larger scales (lower ), and Ly forest observations are also limited on the minimum scale that can be measured accurately which is around s/km (Lee et al., 2012; Palanque-Delabrouille et al., 2013; D’Aloisio et al., 2018b). We will compare current high- large scale measurements with our models in the next section.

Figure 7: The 1D flux power spectrum at , and (from left to right respectively) for a set of different flash (dashed lines) and inhomogeneous models (solid lines). Color of each line stands for the heat injected during reionization: K, K, K (blue, green and red respectively). The reionization redshift of the flash reionization simulations agrees with the median reionization redshift of the inhomogeneous runs. The lower panels of each panel show the percentage difference between each flash simulation and its analogous inhomogeneous run i.e., the one with the same heat input during reionization. Notice that inhomogeneous and flash reionization differ at low- modes (large scales) due to temperature fluctuations but still agree at high- scales (small scales).

For high- modes (small scales) the power spectrum shows a larger scale cut-off (i.e. toward lower ) for simulations with a higher heat input during reionization regardless of the reionization methodology used. The different heat injection leads to different levels of pressure smoothing and thermal broadening (see Fig. 6), changing the shape of the small-scale (high-) cutoffs in the power spectra (Oñorbe et al., 2017b). Interestingly, the power spectrum at s/km is very similar for inhomogeneous and flash models when using the same heat input. At the largest difference is between the models with the largest heat input during reionization and it decreases towards lower redshift. We have found that this similarity occurs when the inhomogeneous and flash models have the same heat input due to reionization, , and the median redshift of reionization of the inhomogeneous run is equal to the reionization redshift of the flash simulation, i.e., . This has very important implications as it suggests that the thermal fluctuations associated with the inhomogeneous reionization models do not significantly change the small scale correlation properties of the Ly forest. It shows that once reionization is finished the small scale properties of the Ly forest in an inhomogeneous reionization model can be reasonably well-described with a flash model that shares the same average heat injection and reionization time. It is likely that the lowest (i.e., the largest scale) at which both types of models agree depends on the exact morphology of the reionization redshift field and the relevant scale of temperature fluctuations which arises from it. Below we will further explore this interesting result using different inhomogeneous reionization models.

5 Different inhomogeneous models

We now turn to the results of our simulations using different inhomogeneous reionization models: IR-A, IR-B, IR-C and IR-D in which we have varied the dominant halo masses responsible for reionization (, and respectively) and adjusted the ionizing efficiency such that reionization is half complete () at a very similar redshift (i.e. all have a very similar , see Figure 3 and Section 2 for further details). We run all these models with the same heat injection during reionization, our fiducial value K. To illustrate the differences between these runs the upper row of Figure 8 shows the same random slice of the reionization redshift field for the four simulations. These slices reveal the different characteristic morphologies for these models where we can see a clear trend with the characteristic minimum halo mass assumed to create them. As expected, larger minimum halo masses result in larger characteristic ionized bubbles. Figure 8 shows the same slices of the (middle row) and temperature (bottom row) fields for these runs at , when the ionization fraction between the models is very similar. In the lower row we can see that as expected, the morphology in the field translates directly to the temperature field at these redshifts, which produce temperature fluctuations persisting to lower redshifts. We now explore whether we can use the high- Ly forest to discriminate between these different reionization histories and morphologies. First we will check the evolution of the thermal parameters of these models.

(fast)

IR-C IR-A IR-D IR-B

Figure 8: Top row: slices of the reionization redshift field for some of the inhomogeneous reionization simulations studied in this work ordered from left to right by the minimum halo mass responsible of reionization: IR-C, IR-A, IR-D and IR-B. In these models we have varied the dominant halo masses responsible for reionization ( , , (fast) and respectively) and adjust the ionization efficiency so that they have a similar global reionization history (see Figure 3). As expected reionization models dominated by smaller halo masses show morphologies with smaller ionized bubbles. Middle row: the same slices as in the top row, now showing the hydrogen ionized fraction for these simulations at when they all have very similar overall ionization fraction between the models is very close. Bottom row: the same slices as above showing the temperature field at . We can clearly see how the typical ionized bubble size of the models directly translate into a different temperature scale.
Figure 9: Thermal history for the inhomogeneous reionization simulations presented in this work: IR-A, IR-B, IR-C and IR-D. From top to bottom it shows the temperature at mean density,, the slope of the density-temperature relation, and the pressure smoothing scale, , as defined in Kulkarni et al. (2015).

Figure 9 shows the evolution of the thermal parameters for the different inhomogeneous reionization simulations: IR-A, IR-B, IR-C and IR-D (green, violet, orange and cyan lines respectively). We also show our fiducial flash reionization simulation for comparison (FR, dashed green lines). The top two panels show the evolution of and parameters describing the temperature-density relation of the IGM (from top to bottom). We can see that while all models share the same evolution, they have slightly different values at fixed redshift, with the largest value corresponding to the flash model. We have checked that this is due to a larger scatter in temperature for a fixed density in the inhomogeneous models. Models that finished reionization at slightly later times (IR-B, IR-C, see Figure 3) have larger scatter which translates into smaller values fitted for the temperature-density slope, (see discussion above). The third panel of Figure 9 presents the evolution of the pressure smoothing scale, , as defined by Kulkarni et al. (2015). We see here a similar trend as the one found for the temperature-density slope, . Models that finish reionization at earlier times show a slightly higher value of . As discussed above, this results from the dependence of the IGM pressure-smoothing scale on the full thermal history (Hui & Haiman, 2003; Kulkarni et al., 2015) and not just on the instantaneous temperature, and it could have important consequences for the statistics of the Ly forest (Nasir et al., 2016; Oñorbe et al., 2017b).

Figure 10: The 1D flux power spectrum at redshifts , and for a set of different inhomogeneous reionization simulations: IR-A (green solid line), IR-B (magenta solid line), IR-C (orange solid line) and IR-D (cyan solid line) with the same heat input during reionization, K. In these models we have varied the dominant halo masses responsible for reionization (, and respectively) and adjust the ionization efficiency so that they have a similar global reionization history (see Figure 3). We also show the 1D flux power spectrum for our FR flash reionization simulation, , that also has the same heat input during reionization as the inhomogeneous runs. The bottom panels at each redshift show the percentage difference between each model with our fiducial inhomogeneous simulation (IR-A).

In order to explore the possibility of distinguishing between these scenarios using Ly forest statistics we have also computed the 1D flux power spectrum for these simulations following the same methodology described in Section 4. Figure 10 shows the 1D flux power spectrum at redshifts , , and . The models differ at large scales (low modes) with the differences decreasing at lower redshifts. The bottom panels show the percentage difference between each model and our fiducial inhomogeneous simulation (IR-A). The largest mode below which these differences are relevant is also decreasing as we go to lower redshifts. This is the same behavior that we found when we compared flash and inhomogeneous models (see Figure 7). These differences are due to the temperature fluctuations resulting from inhomogeneous reionization. In fact, we find that the differences between the inhomogeneous models are directly related to the duration of reionization, . The models with a more extended reionization (e.g. IR-B, see Figure 3) have more power at large scales. This is because the difference in temperature between the regions that reionized early and the ones that reionized last is larger (D’Aloisio et al., 2015). We found differences up to % in the power spectrum at s/km between flash models and the most extended inhomogeneous reionization model (IR-B). While we are limited by the box size of our simulations ( Mpc/h) these results clearly indicate that using the 1D flux power spectrum at high redshift could provide valuable constraints on the extent of reionization.

The top left panel of Figure 10 also shows the recent measurements of the 1D flux power spectrum computed by D’Aloisio et al. (2018b) at using the 21 unique quasar spectra presented in (McGreer et al., 2015). A precise fit to this data is beyond the scope of this paper, as different degeneracies should be considered at the same time (e.g mean flux, cosmology, etc). However we think that by comparing with our different models we can illustrate that based on the current observational precision and the differences that we find between our models we should be able to distinguish between inhomogeneous reionization scenarios, especially considering global fits over a range of redshifts. Note that, contrary to the need for accurate measurements of the cutoff of the 1D flux power spectrum (high- modes, e.g. Walther et al., 2018b), high resolution measurements are not critical and moderate resolution spectra could provide the more precise measurement of the relevant low modes. In this regard, it is also worth mentioning that the data presented here is only using a small subset of all data currently available.

Very interestingly, at small scales (high modes) all the models show a very good agreement and are within of each other. This indicates that the cutoff of the 1D flux power spectrum is mainly sensitive to the heat input during reionization, , and the median reionization redshift, . This is consistent with our findings comparing flash and inhomogeneous reionization models with different heat injection (see Figure 7). This similarity suggests that the heat input by  reionization may be constrained from the Ly forest with relatively simple flash models. Finally, we have not found any relevant features that indicate that the 1D flux power spectrum at is particularly sensitive to the different specific morphologies of the reionization models with different source populations, but this may change with larger simulation volumes.

6 Temperature and Density fields at high- in Reionization Simulations

In order to further investigate the relation between the large scale properties of the Ly forest and the thermal properties of the IGM we have studied in more detail the properties of the density and temperature fields of our simulations and the relation between them.

Our first objective is to exclude the possibility that the large scale power observed in the inhomogeneous runs is due to the scatter of the temperature-density relation and is in fact due to large-scale (i.e. coherent) temperature fluctuations. For this purpose we start with the temperature field of the IR-A inhomogeneous reionization run at . We then generated a new temperature field by randomly shuffling the temperature values777We exclude cells with temperatures higher than K and densities higher than from this analysis as we want to focus on the IGM properties and exclude shock heated cells. In any case the results presented here do not change if we include also these cells in narrow density bins (width in of ). As a result, this new painted temperature field that we obtain, shares the global thermal properties of the original one but it has lost any spatial correlations with the large-scale density field. The two top panels of Figure 11 illustrate this with a random slice of the original temperature field (left panel) and the same slice for the shuffled temperature field (right panel) at . The temperature-density relation between these two models is indistinguishable by construction and therefore and parameters are identical. Since both models share also the exact same density field, the gas pressure effects are also identical888Notice however that is not identical between runs, since in its definition the field is used and this will not be the same in both runs. We will further discuss this below..

Figure 11: Painting the temperature field. Left panel shows a temperature slice of the inhomogeneous run (IR-A) at . Right panel shows an slice of the new painted temperature field constructed to reproduce the same temperature-density relation as in the inhomogeneous run by shuffling its temperature values at fixed density. Although the two temperature-density relations are indistinguishable by construction, the spatial correlation between the two fields is very different for each case. See text for more details.

We now compute the 1D flux power spectrum for the new model with the shuffled temperature field. We generate Ly forest skewers using the new temperature field and the original density and velocity fields of the inhomogeneous run as well as the same UVB. We then renormalize the flux field to the observed mean flux at the respective redshift, as we have done for the original simulation (see §4.1). The left panel of Figure 12 shows the 1D flux power spectrum at of this new model (blue solid line) compared with the original inhomogeneous run (solid green line) and the flash run (dashed green line). The differences at large-scale between the flash and inhomogeneous runs were already discussed in Section 4 (see also Figure 7) where we argued that they are due to the temperature fluctuations in the inhomogeneous reionization models. The 1D flux power spectrum of the model with the painted temperature field provides very interesting insights into the physical origin of the most relevant features of the 1D flux power spectrum, both at large and small scales. First, on large scales (low- modes) its shape seems to be more in agreement with the flash reionization model. The lack of large scale power corroborates the idea that the extra power found at large scales is due to coherent temperature fluctuations. On small scales the model with the shuffled temperatures shows significant extra power compared to both the flash and inhomogeneous models.

Figure 12: Temperature fields in reionization simulations. Left panel: the 1D flux power spectrum of our fiducial flash reionization simulation (FR, dashed green line), the inhomogeneous reionization simulation (IR-A, solid green line) and the painted temperature field (solid blue line). Right panel: The 3D power spectrum of the temperature field fluctuations (i.e., ) for the same simulations.

Since the overall normalization of the flux field obscures the effect of the different temperature fields we have also computed the 3D power spectrum of the temperature field fluctuations (i.e., ) at for these three simulations: the fiducial flash (FR) and inhomogeneous (IR-A) reionization simulations and the new model with the shuffled temperature. The right panel of Figure 12 shows these results. The first thing to notice is that the power spectrum obtained from the model in which we painted the temperature field has a very different overall power as the other two models. The power signal is reduced at all relevant scales and only increases at very large values ( s/km). Regarding the power spectrum at large scales ( s/km), only the fiducial inhomogeneous run, IR-A, presents a significant signal. This additional power comes from the temperature fluctuations due to inhomogeneous reionization and it is missed in both the flash model and the model in which we shuffled the temperature field. It is also interesting to focus on the small scale differences between the models. The flash and inhomogeneous models agree quite well for ( s/km) confirming our conclusions from the analysis of the 1D flux power spectrum of these two type of models: for a similar heat injection during reionization the small scale thermal properties are very similar as long as . On top of the overall normalization mentioned above, the shape of the 3D temperature power spectrum at small scale ( s/km) for the shuffled temperature model seems to be steeper as we go to lower than the one we find for the inhomogeneous or flash runs. Since we expect structure at these scales to be naturally suppressed, it is not surprising that by shuffling the temperature values at a fixed density from the inhomogeneous run we will introduce some extra power at small scale.

We now want to further study how the density field is affected by the reionization process. Following the analysis that we have done with the temperature field we want to compute the 3D power spectrum of the gas density fluctuations field (i.e., ). However the gas power spectrum at small scales in cosmological simulations is dominated by high density regions (see for example Kulkarni et al., 2015) and since we are interested in the properties of the IGM we have computed the power spectrum clipping the density at different thresholds. This allows us to explore the properties of the gas at lower densities corresponding to the IGM. We set all densities above the threshold to this exact value and then compute the 3D power spectrum of the gas density fluctuations of the resulting field (therefore is different for each threshold). The left panel of Figure 13 shows the 3D power spectrum of the gas density field fluctuations for the flash (FR) and inhomogeneous (IR-A) simulations at . We find that the gas density field of both simulations have almost identical properties. This confirms that the gas pressure effects in both simulations is on average very similar despite the different reionization histories.

The resulting 3D clipped gas power spectrum shows a clear cutoff which, we argue, describes the scale at which the pressure smoothing is relevant for that density. As we decrease the value of , this cutoff moves to lower modes (i.e., the gas pressure scale increases) until we reach when the pressure scale seems to decrease again. This cutoff can be fitted by the function suggested in Kulkarni et al. (2015) for the pseudo real-space Ly flux field:

(2)

which has three parameters: , , and . We define the scale associated with this cutoff as . This naturally brings out the pressure smoothing scale - density relation present in the density field (e.g. Gnedin & Hui, 1998b). The right panel of Figure 13 shows the relation between this cutoff scale and the clipping overdensity measured at for the flash (FR, FRcold, FRhot, dashed lines) and inhomogeneous reionization (IR-A, IR-Acold, IR-Ahot, solid lines) simulations that we discuss in Section 4. As expected models with a higher heat injection during reionization have a larger pressure scale, especially at densities below . Interestingly, the analogous flash and inhomogeneous simulations, i.e. the ones that share the same heat injection during reionization, have similar values at mean density. We think that the power spectrum of the smoothed gas density field could be a complementary and more detailed description of the gas pressure support in hydrodynamical simulations as it directly connects with the gas density field. It has the advantage that it is well defined at any redshift and it does not have any second order dependence with the UV background or the temperature fields. A full detailed study of this possibility is beyond the scope of this paper but deserves theoretical consideration for future work as it could help characterize in a much more accurate way the thermal properties of the IGM in hydrodynamical simulations.

Figure 13: Left: The 3D power spectrum of the gas density field fluctuations (i.e., ) at clipped at different density values for the FR flash reionization simulation (dashed lines) and the IR-A inhomogeneous reionization simulation (solid lines). Bottom panel show the difference between the flash and inhomogeneous runs for each different clipped density value. Right: the cutoff scale of the different clipped density fields, , versus the density value used to clip the field. Our tests show that this relation could help quantifying the gas pressure scale - density relation in hydrodynamical simulations. See text for more details.

7 UVB Fluctuations

Until now in our models we have considered a uniform UVB that evolves with redshift999In the inhomogeneous reionization simulations, if a resolution element is still not reionized the UVB is considered to be zero, but all resolution elements see the same UVB once they are reionized.. As we discuss above, recent observations of the optical depth at redshift indicate that the photoionization and photoheating rates can not be well described by an uniform field, high- studies (see e.g. Becker et al., 2018). In this section we try to explore the effect that UVB fluctuations can have in the context of our new simulations.

We model fluctuations of the ionizing background in post-processing following the approach of Davies & Furlanetto (2016). In this approach we consider spatial variations of the mean free path, , from modulations in the ionization state of optically thick absorbers assuming that , where is the local matter density. We refer the reader to this paper for technical details.

Using this approach we have created one model with a mean free path of cMpc. This is roughly a factor of 4 lower than the extrapolation of the measurements in Worseck et al. (2016) to this redshift and in this sense it can be considered an extreme model101010The box size of our simulations is too small to fully simulate strong UVB fluctuations, but we can approximate it by shrinking the mean free path.. However our goal in this work is to explore the possible effect that UVB fluctuations can have on the 1D flux power spectrum in the context of our new inhomogeneous and flash reionization models. We compute the fluctuating photoionization rate, , at on a uniform grid of , which resolves the large-scale fluctuations in the radiation field, and scale this field up to a (i.e., the same size of the hydrodynamical simulation) via linear interpolation in log space. The left panel of Figure 14 shows a slice of this field and shows that for this model the difference in the photoionization rates between regions can be as large as an order of magnitude. From this fluctuating UVB we compute a new neutral hydrogen density, , for our fiducial flash (FR) and inhomogeneous reionization (IR-A) simulations using their original temperature and density fields for both simulations.

Using the new neutral hydrogen density we compute the optical depth and the 1D flux power spectrum for each of these simulations. The right panel of Figure 14 compares the 1D flux power spectrum of these two new models with the original non-fluctuating ones. All models were normalized to the same mean flux. For the flash simulation the effect seems to be an overall rescaling of the power spectrum at all scales. There is a small bump in power at low- ( s/km), however this scale is very close to the maximum scale measured by the simulation. The increase suggests that the effects of the UVB fluctuations could be more significant at even lower (larger scales), but we will need a larger simulation volume to further study this effect (D’Aloisio et al., 2018b). In the case of the modified inhomogeneous reionization simulation however the effects are more clear, with a change of the overall shape of the power spectrum. There is clear drop in the power at large scales (lower ) compared with the original model indicating that the effect of thermal fluctuations in the power is canceled by the photoionization fluctuations. This cancellation is a result of the opposite direction in which both effects act to increase or decrease Ly forest opacity on large scales (Davies & Furlanetto, 2016; Davies et al., 2018a). UVB fluctuations tend to decrease the opacity in overdense regions with more sources of ionizing photons, and increase the opacity in underdense regions that are distant from sources. Temperature fluctuations act in the opposite direction: overdense regions reionize early and are cold and opaque, while underdense regions reionize late and are hot and transparent. It is interesting to note that the UVB fluctuations do not seem to produce significant changes of the 1D flux power spectrum at small scales ( s/km) for the inhomogeneous run. The effect in this range of scales is larger when the UVB are applied to the flash reionization. The rise in power in the flash runs could be related to more dense/evolved parts of the universe dominating the transmission. However the similar power between the inhomogeneous runs at intermediate and small scales is not clear. This could be linked with the dispersion of the temperature-density relation, indicating that the small scale structure of the Ly forest in models with very tight relation is much more sensitive to UVB fluctuations than in models which have more scatter about their temperature-density relations.

Figure 14: Adding UVB fluctuations in post-processing to our fiducial flash (FR) and inhomogeneous (IR-A) reionization models. Left panel: A slice showing the Mpc model of UVB fluctuation field used to post-process the hydrodynamical simulations. Right panel: The 1D flux power spectrum at for the FR flash reionization model with and without UVB fluctuations (dashed green and solid red lines respectively), and the the IR-A inhomogeneous reionization model also with and without UVB fluctuations (dashed blue and solid orange lines respectively). See text for more details.

8 Discussion

In this work we have combined excursion set semi-numerical models with cosmological hydrodynamical simulations in order to study the properties of the IGM and the Ly forest to learn about reionization. The combination of box size (L Mpc/h) and spatial resolution ( kpc/h) chosen in our test runs have been mainly driven by the convergence constraints of the small scale Ly forest properties. We will first discuss how this can limit the methodology presented in this work as well as affect the results.

The first relevant effect to be considered is the convergence of the reionization histories and field properties produced by our semi-numerical approach that we used as input in our hydrodynamical simulations. Results from 3D radiative transfer cosmological simulations argued that box sizes of L Mpc/h or larger are needed to yield convergent reionization histories (e.g. Iliev et al., 2014). This is of course going to be related with the typical size of ionization bubbles in the reionization model, which is related with the sources responsible for reionization. The box size constraints cited above assume that galaxies in M halos are the main drivers of reionization. From a purely computational perspective the goal of reaching L Mpc/ simulations while maintaining the necessary spatial resolution to resolve the Ly forest is already doable. However if larger halos were responsible for reionization, as in quasar dominated scenarios, the box size needed for convergence will be larger. This is in fact the main limitation to currently use the method introduced in this work to simulate inhomogeneous reionization models in cosmological hydrodynamical simulations. reionization is driven by rare, luminous quasars (see e.g. Madau & Meiksin, 1994) and therefore a a much larger box is needed to obtain converged reionization models, L Mpc/ (see e.g. Dixon et al., 2014; Davies et al., 2017; D’Aloisio et al., 2017). A short term solution for this problem could be to use large-scale zoom-in cosmological simulations that can generate accurate and inhomogeneous reionization models and also Ly forest predictions in a smaller sub-volume.

In order to explore the effect of box size on the results from the hydrodynamical simulations we ran the same flash and inhomogeneous reionization models presented in this work in a smaller box (L Mpc/h) but using the same spatial resolution as in the larger boxes ( kpc/h). We did not find any major difference in the main results shown in this work. We found convergence at the level for the 1D flux power spectrum at 111111Notice that this quoted convergence level is the worst one at the highest mode studied, s/km, but it decreases as we move to lower values. which agree with similar recent tests using the Nyx code (Oñorbe et al., 2017b; Walther et al., 2018a). We refer to Lukić et al. (2015) for a more detailed discussion on resolution and box effects in the Ly forest properties using Nyx. Similar results at these redshifts but for simulations using the Gadget code can be found in Bolton & Becker (2009) and Bolton et al. (2017).

There have been several studies to test the accuracy of semi-numerical reionization models compared with radiative transfer simulations. These comparisons have shown that despite the initial general idea that semi-numerical methods will not be accurate to recreate the ionization history, since they do not chronologically follow the state of ionization at individual grid cells, they produce ionization histories that agree on the percentage level when compared with radiative transfer codes (Zahn et al., 2011; Majumdar et al., 2014; Hutter, 2018). The agreement of the morphology of reionization between both methods is however more complicated. Zahn et al. (2011) compared Mpc/h box size 3D radiative transfer cosmological simulations with excursion set simulations and found a good agreement with the cross-correlation coefficient of the ionization fields between /Mpc. The averaging of the density and ionizing emissivity fields in semi-numerical methods results in more connected ionized regions, and strengthens the inside-out character of the ionization topology relative to radiative transfer simulations. More recent comparisons (Majumdar et al., 2014; Hutter, 2018) between radiative transfer codes and new more updated semi-numerical methods, which used pre-run halo catalogs, have confirmed these results. However using higher resolution runs Hutter (2018) found that the correlation between both methods decreases towards smaller scales /Mpc mainly for two reasons. First, the more fractal shape of the ionization boundaries in the radiative transfer simulations leads to more small-scale structures in the ionization fields and second, as reionization proceeds the ionized regions in both methods merge at different points in space and time, and consequently develop increasingly different shapes with time. It is therefore necessary to perform a detailed comparison of the Ly forest properties between 3D cosmological radiative transfer simulations with high enough spatial resolution and the method introduced in this work by using the reionization redshift field produced by the radiative transfer runs. We consider this to be beyond the scope of this paper and we plan to do such a comparison in future work.

It is also interesting to consider the time scale of the reionization process in the context of our method as the spatial and time step resolution used in both generating the reionization redshift fields and in the hydrodynamical simulations will set up an upper limit on how well the exact morphology of the reionization models can be tracked in the simulation. Therefore we will now discuss the expected velocity of ionization fronts from previous works and then put that in context of our results. First analytical estimates of the typical velocity of ionization fronts in the IGM during reionization give values of h ckpc/Myr ( km/s, see e.g. Strömgren, 1939). This is in good agreement cosmological simulations with 3D radiative transfer that obtain velocities from h ckpc/Myr ( km/s) around individual sources at first but accelerate up to h ckpc/Myr ( km/s) as larger bubbles form (see e.g. Iliev et al., 2006; McQuinn et al., 2007; Trac & Cen, 2007; Hirata, 2018). Recently, Kaurov (2016) presented a detailed analysis of the time scale of the reionization process using the results of a 3D radiative transfer cosmological simulations. These authors computed the time delay between the moments when each of the cells of the simulation (with a size of ckpc/h) reaches the 10% and 90% ionization. They found a distribution of values that peaks at Myr (i.e., implying an ionization front velocity of h ckpc/Myr, km/s) with most of the cells ranging between and Myr but with an extended tail to values as high as Myr121212The time intervals between the snapshots limited the lowest value that could be measured to Myr but this does not affect the main results.. While to generate the reionization redshift fields from the excursion set we used a very small time step integration, , the spatial resolution, ckpc, is the current limiting factor (see Section 2.2). Regarding the hydrodynamical simulations, the spatial resolution of all runs presented in this work is ckpc/h and the typical timestep of these simulations during reionization is Myr131313This is assuming Courant factors of .. Therefore the lower limit for the ionization front velocity in our simulations is h ckpc/Myr ( km/s). This is only a factor of larger than the fastest ionization front velocities expected and therefore should be enough to capture the thermal structure of reionization.

Recently D’Aloisio et al. (2018a) have studied in detail the peak gas temperatures behind the ionization fronts, T, using 1D high-resolution radiative transfer simulations. They have found that T is only mildly sensitive to the spectrum of incident radiation over most of the parameter space, with temperatures set primarily by the ionization front speeds. They measured the velocity of ionization fronts in 3D cosmological radiative transfer SCORCH simulations (Doussot et al., 2017) finding also a broad range of velocities during the early phases of reionization, from to h ckpc/Myr ( to km/s). However the velocity increases to h ckpc/Myr ( km/s) and the distribution narrows by the end of reionization. Mapping these velocities to temperatures yields T K during the first half of reionization, but hotter temperatures of T K can be reached during the very last phases of reionization. These findings are extremely interesting in the context of this work and they indicate that assuming just one value of T during reionization, as we have done in our simulations, could be too simplistic. However these findings are easy to implement in the context of the method introduced here and therefore also allow for a more accurate, but yet cheap, modeling of the process in cosmological hydrodynamical simulations. We definitely plan to explore this issue in detail in the near future.

Finally, in order to further study the effects on small scale Ly statistics between our different reionization models we have also calculated the curvature flux statistics (Becker et al., 2011). We have confirmed that all results presented in this work for the shape of the 1D flux power spectrum cut-off regarding the comparison between different reionization models hold also for the curvature statistics.

9 Conclusions

In this paper we have introduced a new method to model inhomogeneous reionization scenarios in optically thin cosmological hydrodynamical simulations which does not increase their cost. In this method each resolution element in the simulation is assigned its own reionization redshift. To account for the heating produced by the ionization front during the first time step of the simulation in which the reionization redshift is crossed a fixed amount of energy is injected, , a free parameter in the model. From the first time step after crossing the reionization redshift and onwards, the resolution element is assumed to be optically thin and can see the homogeneous UV background.

Using this method we have produced a suite of and Mpc/ Nyx cosmological hydrodynamical simulations with different reionization models parameterized by the mid point of reionization, , its duration, and the heat injection, . To generate different redshift reionization models we have used semi-numerical methods that combine the excursion set formalism with the initial conditions of the simulation and approximate radiative transfer methods (e.g. 21cmFAST, Mesinger et al., 2011), which allow us to efficiently explore the parameter space. In order to compare with the inhomogeneous models we also run a set of instantaneous or flash reionization models in which we fixed all resolution elements in the simulation to same the reionization redshift.

We have first analyzed the temperature and density fields and their relationship for all these simulations during and after reionization. As expected, inhomogeneous reionization models enhance the temperature fluctuations found in the IGM once reionization is completed, as regions that reionized early have had time to cool while regions that reionize at late times are still warmer. For this reason reionization models with the same heat injection but that occurred over a longer period of time produce larger temperature fluctuations (D’Aloisio et al., 2018b). Temperature fluctuations alter the temperature-density relation in the IGM after reionization, significantly increasing its scatter compared with flash reionization models.

To study the effects of these different reionization models on the properties of the Ly forest, we compute the Ly 1D flux power spectrum at , i.e. after reionization is finished, for our simulation ensemble. We found that thermal fluctuations can produce significant increase of power (up to ) at large scales, s/km. Reionization models with larger thermal injection have more power at these scales with the differences decreasing with redshift as the fluctuations decrease. We quantify these differences for our suite of simulations and show that using available observational data could already provide interesting constraints on the duration of reionization.

We have confirmed that the small scale differences in the 1D power spectrum (high modes) are driven by the thermal history of the IGM and therefore by studying them we can constrain when reionization happened and how much heat was injected (Nasir et al., 2016; Oñorbe et al., 2017b; Boera et al., 2018). Moreover, we find that once reionization is finished, the small scale properties of the Ly forest in inhomogeneous and flash reionization are very similar provided they share the same heat injection during reionization, , and the reionization in the flash models happens at the median reionization redshift of the extended model, i.e, . This result indicates that there is a scale below which the correlations due to inhomogeneous reionization average out and that the the Ly forest properties are well described by an average of the density dependent thermal properties on these scales. We have shown that this scale could be as large as s/km at and therefore relevant for available observations and constraining not only reionization but other physical processes relevant at these scales (e.g. nature of dark matter).

To further investigate these findings we have compared the properties of the density and temperature fields of the different reionization models. The spatial correlations of these fields at large and small scales show a clear correlation with the Ly forest properties. We confirmed that the difference in the large scale power between the inhomogeneous and flash models is indeed due to these correlations and not the larger scatter in the temperature-density relation that naturally arises in the inhomogeneous models. Moreover we have shown that the spatial correlations of the density and temperature fields and their relation are crucial to properly model the small scales of the Ly forest. In fact we find that the classical paradigm of describing the thermal state of the IGM with three parameters, two for the temperature-density relationship, plus one describing the pressure scale of the gas, does not fully capture the differences between models once inhomogeneous reionization models are considered with high enough precision.

We also studied the effects of UVB fluctuations in our simulations to see how they could affect the Ly forest. We construct an extreme model ( Mpc) that we applied in post-processing to our fiducial flash and inhomogeneous reionization simulations and compare the 1D flux power spectrum between all these models. We find that UVB fluctuations alter the overall normalization of the power spectrum by changing the density - optical depth mapping, producing effects on all scales. It is interesting however that the UVB fluctuations seem to cancel the large-scale effects that appear in the inhomogeneous reionization model due to thermal fluctuations. In any case these results indicate that this effect should be also considered when using the 1D flux power spectrum at high- to constrain reionization, cosmological parameters and/or the nature of dark matter. It is worth noting that thermal fluctuations will last way past the end of reionization while UVB fluctuations will decay faster, therefore in principle, these effects could be disentangled.

Finally, we believe that the new method introduced in this work is an important step forward to accurately model the and reionization processes in cosmological hydrodynamical simulations. The method is easy to implement and adds no significant extra computing cost compared with the standard approach used in IGM and galaxy formation and evolution studies which assumes a simple homogeneous reionization model. We plan to further explore the limitations of this approach both by implementing more accurate thermal heating models and comparing with high resolution 3D radiative transfer cosmological simulations.

Acknowledgements

We thank the members of the ENIGMA group at the Max Planck Institute for Astronomy (MPIA) for helpful discussions. Calculations presented in this paper used the hydra and draco cluster of the Max Planck Computing and Data Facility (MPCDF, formerly known as RZG) MPCDF is a competence center of the Max Planck Society located in Garching (Germany). We also used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The ASCR Leadership Computing Challenge (ALCC) program has provided NERSC allocation under the project: “Cosmic Frontier Computational End-Station”. This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org.

References

  • Abel & Haehnelt (1999) Abel T., Haehnelt M. G., 1999, ApJ, 520, L13
  • Almgren et al. (2013) Almgren A. S., Bell J. B., Lijewski M. J., Lukić Z., Van Andel E., 2013, ApJ, 765, 39
  • Alvarez et al. (2012) Alvarez M. A., Finlator K., Trenti M., 2012, ApJ, 759, L38
  • Bañados et al. (2014) Bañados E., et al., 2014, AJ, 148, 14
  • Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  • Battaglia et al. (2013) Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81
  • Bauer et al. (2015) Bauer A., Springel V., Vogelsberger M., Genel S., Torrey P., Sijacki D., Nelson D., Hernquist L., 2015, MNRAS, 453, 3593
  • Becker & Bolton (2013) Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023
  • Becker et al. (2011) Becker G. D., Bolton J. S., Haehnelt M. G., Sargent W. L. W., 2011, MNRAS, 410, 1096
  • Becker et al. (2015) Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Venemans B. P., 2015, MNRAS, 447, 3402
  • Becker et al. (2018) Becker G. D., Davies F. B., Furlanetto S. R., Malkan M. A., Boera E., Douglass C., 2018, ApJ, 863, 92
  • Boera et al. (2018) Boera E., Becker G. D., Bolton J. S., Nasir F., 2018, preprint, (arXiv:1809.06980)
  • Bolton & Becker (2009) Bolton J. S., Becker G. D., 2009, MNRAS, 398, L26
  • Bolton et al. (2017) Bolton J. S., Puchwein E., Sijacki D., Haehnelt M. G., Kim T.-S., Meiksin A., Regan J. A., Viel M., 2017, MNRAS, 464, 897
  • Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  • Bosman et al. (2018) Bosman S. E. I., Fan X., Jiang L., Reed S., Matsuoka Y., Becker G., Haehnelt M., 2018, MNRAS, 479, 1055
  • Cen et al. (2009) Cen R., McDonald P., Trac H., Loeb A., 2009, ApJ, 706, L164
  • Chardin et al. (2015) Chardin J., Haehnelt M. G., Aubert D., Puchwein E., 2015, MNRAS, 453, 2943
  • Chardin et al. (2017) Chardin J., Puchwein E., Haehnelt M. G., 2017, MNRAS, 465, 3429
  • Coc et al. (2013) Coc A., Uzan J.-P., Vangioni E., 2013, preprint, (arXiv:1307.6955)
  • D’Aloisio et al. (2015) D’Aloisio A., McQuinn M., Trac H., 2015, ApJ, 813, L38
  • D’Aloisio et al. (2017) D’Aloisio A., Upton Sanderbeck P. R., McQuinn M., Trac H., Shapiro P. R., 2017, MNRAS, 468, 4691
  • D’Aloisio et al. (2018a) D’Aloisio A., McQuinn M., Maupin O., Davies F. B., Trac H., Fuller S., Upton Sanderbeck P. R., 2018a, preprint, (arXiv:1807.09282)
  • D’Aloisio et al. (2018b) D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018b, MNRAS, 473, 560
  • Davies & Furlanetto (2016) Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328
  • Davies et al. (2016) Davies F. B., Furlanetto S. R., McQuinn M., 2016, MNRAS, 457, 3006
  • Davies et al. (2017) Davies F. B., Furlanetto S. R., Dixon K. L., 2017, MNRAS, 465, 2886
  • Davies et al. (2018a) Davies F. B., Becker G. D., Furlanetto S. R., 2018a, ApJ, 860, 155
  • Davies et al. (2018b) Davies F. B., et al., 2018b, ApJ, 864, 142
  • Dixon et al. (2014) Dixon K. L., Furlanetto S. R., Mesinger A., 2014, MNRAS, 440, 987
  • Doussot et al. (2017) Doussot A., Trac H., Cen R., 2017, preprint, (arXiv:1712.04464)
  • Eilers et al. (2018) Eilers A.-C., Davies F. B., Hennawi J. F., 2018, ApJ, 864, 53
  • Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
  • Faucher-Giguère et al. (2009) Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
  • Feng et al. (2016) Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, MNRAS, 455, 2778
  • Finlator et al. (2018) Finlator K., Keating L., Oppenheimer B. D., Davé R., Zackrisson E., 2018, MNRAS, 480, 2628
  • Furlanetto & Oh (2009) Furlanetto S. R., Oh S. P., 2009, ApJ, 701, 94
  • Furlanetto et al. (2004) Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1
  • Giallongo et al. (2015) Giallongo E., et al., 2015, A&A, 578, A83
  • Gnedin (2014) Gnedin N. Y., 2014, ApJ, 793, 29
  • Gnedin & Hui (1998a) Gnedin N. Y., Hui L., 1998a, MNRAS, 296, 44
  • Gnedin & Hui (1998b) Gnedin N. Y., Hui L., 1998b, MNRAS, 296, 44
  • Gnedin et al. (2017) Gnedin N. Y., Becker G. D., Fan X., 2017, ApJ, 841, 26
  • Greig et al. (2017) Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239
  • Gunn & Peterson (1965) Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633
  • Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
  • Haardt & Madau (2001) Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays. p. 64 (arXiv:astro-ph/0106018)
  • Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  • Hahn & Abel (2011) Hahn O., Abel T., 2011, MNRAS, 415, 2101
  • Hassan et al. (2016) Hassan S., Davé R., Finlator K., Santos M. G., 2016, MNRAS, 457, 1550
  • Hirata (2018) Hirata C. M., 2018, MNRAS, 474, 2173
  • Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  • Howlett et al. (2012) Howlett C., Lewis A., Hall A., Challinor A., 2012, Journal of Cosmology and Astroparticle Physics, 4, 27
  • Hui & Haiman (2003) Hui L., Haiman Z., 2003, ApJ, 596, 9
  • Hutter (2018) Hutter A., 2018, MNRAS, 477, 1549
  • Iliev et al. (2006) Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625
  • Iliev et al. (2009) Iliev I. T., et al., 2009, MNRAS, 400, 1283
  • Iliev et al. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725
  • Iršič et al. (2017) Iršič V., et al., 2017, Phys. Rev. D, 96, 023522
  • Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
  • Kaurov (2016) Kaurov A. A., 2016, ApJ, 831, 198
  • Keating et al. (2018) Keating L. C., Puchwein E., Haehnelt M. G., 2018, MNRAS, 477, 5501
  • Khaire & Srianand (2018) Khaire V., Srianand R., 2018, preprint, (arXiv:1801.09693)
  • Kim et al. (2016) Kim H.-S., Wyithe J. S. B., Park J., Poole G. B., Lacey C. G., Baugh C. M., 2016, MNRAS, 455, 4498
  • Kulkarni et al. (2015) Kulkarni G., Hennawi J. F., Oñorbe J., Rorai A., Springel V., 2015, ApJ, 812, 30
  • Kulkarni et al. (2018) Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I., Puchwein E., Chardin J., Aubert D., 2018, preprint, (arXiv:1809.06374)
  • La Plante et al. (2017) La Plante P., Trac H., Croft R., Cen R., 2017, ApJ, 841, 87
  • Lee et al. (2012) Lee K.-G., Suzuki N., Spergel D. N., 2012, AJ, 143, 51
  • Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  • Lidz & Malloy (2014) Lidz A., Malloy M., 2014, ApJ, 788, 175
  • Lidz et al. (2007) Lidz A., McQuinn M., Zaldarriaga M., Hernquist L., Dutta S., 2007, ApJ, 670, 39
  • Liu et al. (2016) Liu A., Pritchard J. R., Allison R., Parsons A. R., Seljak U., Sherwin B. D., 2016, Phys. Rev. D, 93, 043013
  • Lukić et al. (2015) Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697
  • Madau & Meiksin (1994) Madau P., Meiksin A., 1994, ApJ, 433, L53
  • Majumdar et al. (2014) Majumdar S., Mellema G., Datta K. K., Jensen H., Choudhury T. R., Bharadwaj S., Friedrich M. M., 2014, MNRAS, 443, 2843
  • McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499
  • McQuinn (2012) McQuinn M., 2012, MNRAS, 426, 1349
  • McQuinn (2016) McQuinn M., 2016, ARA&A, 54, 313
  • McQuinn & Upton Sanderbeck (2016) McQuinn M., Upton Sanderbeck P. R., 2016, MNRAS, 456, 47
  • McQuinn et al. (2007) McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007, MNRAS, 377, 1043
  • Meiksin & Tittley (2012) Meiksin A., Tittley E. R., 2012, MNRAS, 423, 7
  • Mesinger et al. (2011) Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955
  • Miralda-Escudé & Rees (1994) Miralda-Escudé J., Rees M. J., 1994, MNRAS, 266, 343
  • Mortlock et al. (2011) Mortlock D. J., et al., 2011, Nature, 474, 616
  • Nasir et al. (2016) Nasir F., Bolton J. S., Becker G. D., 2016, MNRAS, 463, 2335
  • Norman et al. (2015) Norman M. L., Reynolds D. R., So G. C., Harkness R. P., Wise J. H., 2015, ApJS, 216, 16
  • Oñorbe et al. (2017a) Oñorbe J., Hennawi J. F., Lukić Z., 2017a, ApJ, 837, 106
  • Oñorbe et al. (2017b) Oñorbe J., Hennawi J. F., Lukić Z., Walther M., 2017b, ApJ, 847, 63
  • Ocvirk et al. (2016) Ocvirk P., et al., 2016, MNRAS, 463, 1462
  • Palanque-Delabrouille et al. (2013) Palanque-Delabrouille N., et al., 2013, A&A, 559, A85
  • Park et al. (2016) Park H., Shapiro P. R., Choi J.-h., Yoshida N., Hirano S., Ahn K., 2016, ApJ, 831, 86
  • Pawlik et al. (2015) Pawlik A. H., Schaye J., Dalla Vecchia C., 2015, MNRAS, 451, 1586
  • Peeples et al. (2010) Peeples M. S., Weinberg D. H., Davé R., Fardal M. A., Katz N., 2010, MNRAS, 404, 1281
  • Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  • Planck Collaboration et al. (2018) Planck Collaboration et al., 2018, preprint, (arXiv:1807.06205)
  • Poole et al. (2016) Poole G. B., Angel P. W., Mutch S. J., Power C., Duffy A. R., Geil P. M., Mesinger A., Wyithe S. B., 2016, MNRAS, 459, 3025
  • Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  • Puchwein et al. (2018) Puchwein E., Haardt F., Haehnelt M. G., Madau P., 2018, preprint, (arXiv:1801.04931)
  • Robertson et al. (2015) Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19
  • Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  • So et al. (2014) So G. C., Norman M. L., Reynolds D. R., Wise J. H., 2014, ApJ, 789, 149
  • Sobacchi & Mesinger (2014) Sobacchi E., Mesinger A., 2014, MNRAS, 440, 1662
  • Strömgren (1939) Strömgren B., 1939, ApJ, 89, 526
  • Tittley & Meiksin (2007) Tittley E. R., Meiksin A., 2007, MNRAS, 380, 1369
  • Trac & Cen (2007) Trac H., Cen R., 2007, ApJ, 671, 1
  • Trac et al. (2008) Trac H., Cen R., Loeb A., 2008, ApJ, 689, L81
  • Venemans et al. (2015a) Venemans B. P., et al., 2015a, MNRAS, 453, 2259
  • Venemans et al. (2015b) Venemans B. P., et al., 2015b, ApJ, 801, L11
  • Viel et al. (2004) Viel M., Haehnelt M. G., Springel V., 2004, MNRAS, 354, 684
  • Viel et al. (2013a) Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013a, Phys. Rev. D, 88, 043502
  • Viel et al. (2013b) Viel M., Schaye J., Booth C. M., 2013b, MNRAS, 429, 1734
  • Walther et al. (2018a) Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2018a, preprint, (arXiv:1808.04367)
  • Walther et al. (2018b) Walther M., Hennawi J. F., Hiss H., Oñorbe J., Lee K.-G., Rorai A., O’Meara J., 2018b, ApJ, 852, 22
  • Wang et al. (2017) Wang F., et al., 2017, ApJ, 839, 27
  • Weigel et al. (2015) Weigel A. K., Schawinski K., Treister E., Urry C. M., Koss M., Trakhtenbrot B., 2015, MNRAS, 448, 3167
  • Willott et al. (2010) Willott C. J., et al., 2010, AJ, 139, 906
  • Wise et al. (2014) Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560
  • Worseck et al. (2016) Worseck G., Prochaska J. X., Hennawi J. F., McQuinn M., 2016, ApJ, 825, 144
  • Zahn et al. (2011) Zahn O., Mesinger A., McQuinn M., Trac H., Cen R., Hernquist L. E., 2011, MNRAS, 414, 727
  • Zaldarriaga et al. (2001) Zaldarriaga M., Hui L., Tegmark M., 2001, ApJ, 557, 519
Sim method
(K) ()
IR-A Inhomogeneous (Model A) ()
IR-A-noTthres Inhomogeneous (Model A-no T threshold) ()
IR-A-dz Inhomogeneous (Model A-) ()
IR-A-li Inhomogeneous (Model A-lin-interp) ()
IR-A-gs Inhomogeneous (Model A-gauss-smo) ()
Column 1: Simulation code.
Column 2: Method used to simulate reionization. See text for more details.
Column 3: reionization redshift.
Column 4: Width of reionization. .
Column 5: Thompson scattering optical depth density weighted (volume weighted).
Column 6: Total heating assumed for reionization.
Column 7: Parameter excursion set model 1. Efficiency.
Column 8: Parameter excursion set model 2. Minimum mass.
The simulations used for these tests have a box size of length, and resolution elements.
Table 2: Summary of simulations used to test the reionization field.

Appendix A Inhomogenous Reionization Fields and Heat Injection: Tests

In this Section we want to address how different assumptions in the creation of the initial reionization field and how the heat due to reionization is injected can affect the results presented in this paper as well how exactly the heat injection is implemented in the hydrodynamical simulation. All simulations discussed here have a box size length of Mpc/h and resolution elements. Therefore, although they cover a smaller volume than the fiducial runs, they do have the same spatial resolution.

We first want to modify how the heat due to reionization is injected in the simulation. In our default implementation we have considered that there is a maximum temperature that any resolution element can be reached after reionization which corresponds to the free parameter, . We have run a simulation (IR-A-noTthres) modifying how the final temperature of each resolution element, , is computed in the simulation by removing the condition of setting up an upper limit of the final temperature equal to . This basically just translates in changing Equation 1 to:

(3)

where is the temperature of the resolution element just before reionization and is its neutral fraction also before reionization. Therefore now resolution elements can in fact be above the temperature injected due to reionization and a slightly larger scatter in temperatures is produced. The exact parameters of this simulation can be found in Table 2. The top panel of Figure 15 shows a slice of the temperature field at of a simulation run using our fiducial model (IR-A). The left panel in the second row shows the same slice without imposing any temperature threshold. Both slices are very similar, however there are some cells which are reionized at very late times in the simulation that in the fiducial simulation have lower temperature values (e.g. in the lower left part of the plot). However, as the left panel of Figure 16 shows for , when the 1D flux power spectrum for both models are compared they do not show any significant differences at scales reached by observations ( s/km). The right panel of Figure 16 shows the 3D power spectrum of the fluctuations of the temperature field () for these runs. Note that despite having slightly different mean temperature values both simulations show a very similar 3D temperature power spectrum at all scales.

In order to address the effects of the time resolution used when generating the reionization models we create new model (IR-A-dz) in which we use a higher time resolution than our default runs (). The exact parameters of this simulation can be found in Table 2. We first confirmed that using a different time resolution does not change the reionization history of the model. Notice that as in the previous test, in this simulation we used Equation-3 to calculate the temperature after reionization (i.e. we did not consider a temperature threshold). The right panel in the second row of Figure 15 shows a slice of the temperature field at of the new model with higher resolution time. While the effect is very subtle the temperature map is smoother in the run with a higher time resolution. Figure 16 shows the 1D flux power spectrum at for this model (left panel) and the 3D temperature power spectrum (right panel) and demonstrates that the differences between both simulations are very small and only relevant above scales of k s/km, above the modes in which we are interested in this work. The differences in the temperature field are in agreement with the test run in which we only modified the reionization heating (see above) and in any case the.

Finally, the bottom row panels of Figure 15 show the effects of using different smoothing algorithms to generate the high resolution reionization map used as input in the simulations from the one obtained using our semi-numerical approach. Our fiducial model (top panel) does not post-process this field but we explore the effects of applying a a simple linear interpolation scheme (IR-A-li, bottom row left panel) and a gaussian smoothing approach (IR-A-gs, bottom row right panel). Notice that in these two simulations we also used Equation-3 to calculate the temperature after reionization (i.e. we did not consider a temperature threshold). The temperature slices show very clearly the effect of smoothing the reionization field generating a smoother temperature field by .

Figure 16 shows the 1D flux power spectrum (left panel) and the 3D temperature power spectrum (right panel) at for the two smoothed models (dotted orange and red lines for the linear and gaussian interpolation respectively). While the flux power spectrum shows that the differences with our fiducial run are not relevant at the scales studied in this work (k s/km) they make a substantial effect in smaller scales. This is even more clear when we look the 3D temperature power spectrum as at small scales the smoothing is clearly changing its overall shape. Therefore we decide not to smooth any reionization field in our fiducial runs. Although this result is not relevant for the results presented in this work, it raises some doubts about the accuracy at very small scales s/km. It is however still not clear if these differences are just due to the spatial resolution of the hydrodynamical simulation or some intrinsic limitations of the reionization models used in this work. Recent works have shown that pushing observations of the 1D flux power spectrum beyond the classical limit of s/km could significantly improve the different constraints obtained from this observable on reionization and probably other physics (e.g. nature of dark matter) (Nasir et al., 2016). Therefore a detailed analysis of these high resolution effects will be crucial for accurate models of the 1D flux power spectrum.

Figure 15: Slice of the temperature field at for different simulation tests on the small scale structure of the reionization field. Top row: method used in our fiducial runs. Second row left panel: reionization heating in the simulation uses no temperature threshold. Second row righ panel: input reionization field generated with a higher time resolution. Bottom row left panel: input reionization field smoothed using linear interpolation. Bottom row right panel: input reionization field smoothed using gaussian interpolation. See text for more details.
Figure 16: The 1D flux power spectrum at (left panel) and the 3D temperature power spectrum at for set of different simulation tests on our new method to simulate inhomogeneous reionization models in hydrodynamical simulations. Both plots show the following models: our fiducial method (cyan solid line, IR-A), reionization heating in the simulation uses no temperature threshold (dashed green line, IR-A-noTthres), input reionization field generated with a higher time resolution (dotted-dash blue line, IR-A-dz). input reionization field smoothed using linear interpolation (dotted orange line, IR-A-li). input reionization field smoothed using gaussian interpolation (red dotted line, IR-A-gs). We do not find any significant changes in the main results presented in this work for any of these variations. See text for more details.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
316036
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description