The Aurora radiation-hydrodynamical simulations of reionization: calibration and first results.
We introduce a new suite of radiation-hydrodynamical simulations of galaxy formation and reionization called Aurora. The Aurora simulations make use of a spatially adaptive radiative transfer technique that lets us accurately capture the small-scale structure in the gas at the resolution of the hydrodynamics, in cosmological volumes. In addition to ionizing radiation, Aurora includes galactic winds driven by star formation and the enrichment of the universe with metals synthesized in the stars. Our reference simulation uses dark matter and gas particles in a box of size with a force softening scale of at most . It is accompanied by simulations in larger and smaller boxes and at higher and lower resolution, employing up to particles, to investigate numerical convergence. All simulations are calibrated to yield simulated star formation rate (SFR) functions in close agreement with observational constraints at redshift and to achieve reionization at , which is consistent with the observed optical depth to reionization. We focus on the design and calibration of the simulations and present some first results. The median stellar metallicities of low-mass galaxies at are consistent with the metallicities of dwarf galaxies in the Local Group, which are believed to have formed most of their stars at high redshifts. After reionization, the mean photoionization rate decreases systematically with increasing resolution. This coincides with a systematic increase in the abundance of neutral hydrogen absorbers in the IGM.
keywords:cosmology: reionization – methods: numerical – radiative transfer – galaxies: high-redshift – intergalactic medium – HII regions
Reionization is the transformation of the initially cold and neutral gas which resulted from recombination after the Big Bang into the hot and ionized plasma that we observe between galaxies today. Measurements of the electron scattering optical depth from observations of the cosmic microwave background (CMB; e.g., Planck Collaboration et al., 2015) and the detection of extended Gunn-Peterson absorption troughs in the spectra of high-redshift quasars (e.g., Becker et al., 2001; Fan, Carilli, & Keating, 2006) have firmly established that reionization took place in the first billion years. This is a key epoch in the history of the Universe that has also seen the birth of the first stars, the initial enrichment of the cosmic gas with metals synthesized in the stars, and the formation of galaxies, including the progenitors of more massive galaxies like our own galaxy, the Milky Way (for reviews see, e.g., Barkana & Loeb, 2001; Ricotti, 2010; Bromm & Yoshida, 2011; Zaroubi, 2013; Lidz, 2015, McQuinn, 2015).
Understanding the physics of reionization is a primary goal of modern astrophysics. One of the main open questions that drives research in the field is whether the ionizing radiation from stars was sufficient to accomplish reionization, or if there were other, perhaps yet to be discovered, sources at play (e.g., Volonteri & Gnedin, 2009; Loeb, 2009; Madau & Haardt, 2015). Estimates from observations of the early Universe with the Hubble Space Telescope are consistent with a picture in which reionization is driven mainly by star-forming galaxies (e.g., Pawlik, Schaye, & van Scherpenzeel, 2009; Fontanot, Cristiani, & Vanzella, 2012; Finkelstein et al., 2015; Sharma et al., 2015). However, such estimates still require assumptions about, e.g., the escape fraction of ionizing radiation (e.g., Wise & Cen, 2009; Yajima, Choi, & Nagamine, 2011; Paardekooper, Khochfar, & Dalla Vecchia, 2015), absorption of ionizing photons at small scales (e.g., Finlator et al., 2012; Kaurov & Gnedin, 2014; Pawlik, Schaye, & Dalla Vecchia, 2015, hereafter PSD15) and the abundance of faint galaxies below the current detection threshold, all of which are difficult to test observationally (for a recent review see Bouwens, 2015).
Upcoming observations of galaxies with, e.g., ALMA (e.g., Ota et al., 2014), MUSE (e.g., Wisotzki et al., 2015) and JWST (e.g., Finkelstein et al., 2015), and of the intergalactic medium (IGM) with, e.g., LOFAR (e.g., Zaroubi et al., 2012), MWA (e.g., Lidz et al., 2008) and SKA (e.g., Mellema et al., 2013), promise to shed further light on these issues. Theoretical models of the early Universe will be critical to help interpret the data that these telescopes will collect. Cosmological radiative transfer (RT) simulations of reionization have emerged as some of the most powerful approaches to building such models. These simulations track the gravitational growth of density fluctuations in the expanding Universe and the hydrodynamical evolution of the cosmic gas, include recipes for star formation and the associated feedback, and also follow the propagation of ionizing radiation (e.g., Ciardi, Stoehr, & White, 2003; Iliev et al., 2006; McQuinn et al., 2007; Petkova & Springel, 2011; Hasegawa & Semelin, 2013; Gnedin, 2014; PSD15; Aubert, Deparis, & Ocvirk, 2015). The main drawback of these simulations is the extreme computational challenge that they pose, in particular the accurate transport of ionizing photons from multiple sources across a large range of scales in cosmological volumes (for reviews see, e.g., Trac & Gnedin, 2011; Finlator, 2012; Dale, 2015).
In this work we introduce Aurora111Aurora is the Latin word for dawn, and the goddess of dawn in Roman mythology and Latin poetry., a new suite of cosmological radiation-hydrodynamical Smoothed Particle Hydrodynamics (SPH) simulations of reionization. The Aurora simulations track the formation of stars and galaxies at high resolution across cosmological scales, and also follow the hydrodynamically coupled transport of ionizing radiation from stars, the explosion of stars as supernovae (SNe) and the enrichment of the universe with metals.
Aurora addresses the computational challenges imposed by the RT by using the spatially adaptive radiation transport technique TRAPHIC (Pawlik & Schaye, 2008) that solves multi-scale problems efficiently and has a computational cost that is independent of the number of radiation sources. Ionizing radiation is transported at the same high spatially adaptive resolution at which the gas is evolved. The ionizing radiation is coupled hydrodynamically to the evolution of the gas, which lets us accurately account for feedback from photoheating.
For comparison, in most previous reionization simulations the radiation transport was carried out on a uniform grid superimposed on the cosmological simulation (e.g., Iliev et al., 2006; Finlator, Davé, Özel, 2011; Ciardi et al., 2012; Iliev et al., 2014). Because computational resources are finite, such a grid typically implies a substantial reduction in dynamic range in comparison with the underlying spatially adaptive gas simulation and impedes the coupling of the radiation and the dynamics of the gas. In fact, the gas dynamics has often been ignored altogether, assuming the gas traces the dark matter, thus preventing an accurate treatment of the feedback from photoionization heating. Only recently have cosmological simulations started to approach the scales and the high resolution needed to capture the relevant physics (e.g., Gnedin, 2014; Norman et al., 2015; PSD15; Bauer et al., 2015; Ocvirk et al., 2015).
Simulating the internal structure of galaxies in cosmological volumes from first principles, is currently not computationally feasible. This holds even for substantially less expensive simulations in which the accurate transfer of ionizing radiation is ignored. For this reason, cosmological simulations like ours must make use of subresolution models to incorporate processes that are not resolved but which are known to impact the dynamics at resolved scales. Because the coupling between unresolved and resolved physics is often not sufficiently understood, this introduces free parameters. Two parameters that are especially relevant to simulations of galaxy formation during reionization are the efficiency of stellar feedback due to the injection of energy by stellar winds and by SN explosions, and the subresolution escape fraction of ionizing radiation.
In Aurora, we calibrate the efficiency of stellar feedback and the subresolution escape fraction to match current observational constraints on the SFR function at and the redshift of reionization. The calibration is carried out for a set of simulations that cover a range of box sizes and resolutions. This lets us investigate the numerical convergence of our simulations at fixed reionization redshift, which removes a major variable that may complicate the interpretation of convergence tests. Quantities that have not been calibrated, such as the slope of the SFR function at SFRs that are lower than currently observed, the time-dependent distribution of neutral gas, and the strength of the photoionizing background that fills space after reionization, are direct predictions coming out of Aurora. Calibration is key to embedding these predictions in observationally supported models of reionization
Our approach follows that of Schaye et al. (2015), who introduced the notions of strong and weak convergence. Simulations converge strongly if the simulation outcome does not change significantly when the resolution is increased while holding all other parameters fixed. Strong convergence is computationally expensive to achieve and might still say little about the quality of the simulations. This is because subresolution models are typically only applicable within a limited range of resolutions. Weak convergence acknowledges these difficulties by demanding convergence only after an appropriate adaption of the subresolution model. In other words, model parameters that are critical for the outcome but whose values cannot be predicted from first principles should be calibrated, using observations if possible.
Gnedin (2016) also investigated the weak and strong convergence of the SFR in simulations of reionization. Our approach differs from that work because we calibrate subresolution model parameters to match the observed SFR function rather than the theoretically estimated SFR function, which Gnedin (2016) determines by extrapolating the results from a series of simulations of increasing resolution. In addition, unlike Gnedin (2016), we calibrate also the subresolution escape fraction, which is needed for achieving weak convergence in the reionization redshift, the other major observable of the epoch of reionization that currently cannot be predicted from first principles.
This paper focuses on the design of the Aurora simulations and discusses results from the initial analysis of a subset of the simulations carried out in the Aurora programme. A more comprehensive and detailed analysis of the simulations will be presented in future work. There we will also discuss additional simulations, for which we have varied the numerical parameters away from their calibrated values to help gain further insight in the simulation results (similar to what we did in PSD15).
The Aurora simulations start from the same initial conditions as the spatially adaptive radiation-hydrodynamical reionization simulations presented in PSD15. They complement the simulations in PSD15 by means of calibration, which we use to remove the dependence on resolution of the comparison between simulated and observed SFR functions and of the redshift of reionization, two key observables of the early Universe. They improve on them by the increased physical realism, additionally including, e.g., absorption of ionizing photons by helium and chemical enrichment, and by the use of a larger range of box sizes and more particles.
The structure of this paper is as follows. In Section 2 we describe our numerical techniques and the set of simulations that we investigate here. In Section 3 we describe the calibration of the efficiency of the stellar feedback and of the escape fraction of ionizing photons. In Section 4 we describe a set of initial results, including the star formation and reionization history, the build-up of the ionizing background and metal enrichment. We also discuss the properties of the first galaxies. In Section 5 we conclude with a brief summary. Lengths are expressed in physical units, unless noted otherwise.
We use a modified version of the N-body/TreePM Smoothed Particle Hydrodynamics (SPH) code GADGET (last described in Springel, 2005) to perform a suite of cosmological radiation-hydrodynamical simulations of galaxy formation down to redshift (see Table 1 and Figure 1). The simulations presented here start at from the same initial conditions as the simulations published in PSD15. However, the two suites of simulations differ substantially in a number of points.
A major difference is that in the current work we calibrate the stellar feedback strength and the unresolved escape fraction of ionizing radiation separately for each of the simulations. Consequently, all simulations of the current suite yield nearly identical SFR functions and reionization histories, independent of box size and numerical resolution and consistent with observations. This is complementary to our previous work, where the parameters were identical for all simulations and the match between simulated and observed SFR functions and the redshift of reionization both depended on resolution.
Other differences include that in this work we follow the chemistry and cooling of helium in addition to that of hydrogen, including ionization of helium by stellar radiation. New is also that we track the enrichment with metals of the IGM. Finally, the current work makes use of an improved implementation of stellar feedback and also puts stricter limits on the size of the time steps over which the dynamics of the SPH particles are integrated, which improves the accuracy at which the hydrodynamics equations are solved.
In the following summary of our methodology we will focus on the differences with respect to PSD15, to which we refer the reader for details. We adopt the CDM cosmological model with parameters , and , , , and (Komatsu et al., 2011), consistent with the most recent constraints from observations of the CMB by the Planck satellite (Planck Collaboration et al., 2015).
As in PSD15, we adopt the entropy-conserving formulation of SPH (Springel, 2005) and average over 48 SPH neighbor particles. New is that we use the time step limiter described in Durier & Dalla Vecchia (2012) to keep the ratio of time steps of neighboring SPH particles . This ensures an accurate energy-conserving integration of the equations of hydrodynamics (see also Saitoh & Makino, 2009; Springel, 2010).
2.2 Chemistry and cooling
We follow the non-equilibrium chemistry and radiative cooling of the gas, assuming that the gas is of primordial composition with a hydrogen mass fraction and a helium mass fraction . This is a substantial improvement on PSD15, in which we neglected helium to speed up the simulations. We achieve this improvement by replacing the explicit chemistry solver described in PSD15 with the implicit solver described in Pawlik, Milosavljević, & Bromm (2013), which is faster. As in PSD15, we impose a temperature floor above densities to prevent spurious fragmentation in dense underresolved gas, following Schaye & Dalla Vecchia (2008).
2.3 Star formation
Our implementation of star formation is identical to that used in PSD15. This means that we use the pressure-dependent stochastic recipe of Schaye & Dalla Vecchia (2008) to sample from a Chabrier (2003) initial mass function (IMF) in the range and turn gas particles with densities above into star particles, at a rate that reproduces the observed star formation law in simulations of isolated disk galaxies (see Schaye & Dalla Vecchia, 2008 for a detailed discussion).
2.4 Chemical enrichment
Another major improvement on PSD15 is that we account for the age-dependent release of hydrogen, helium, and metals by Type II core-collapse and Type Ia SNe, stellar winds from massive stars, and AGB stars, following Schaye et al. (2010) and Wiersma et al. (2009). We track both the total ejected metal mass, and, separately, the ejected mass of 11 individual elements, including H, He, C, N, O, Ne, Mg, Si, and Fe. Unlike Schaye et al. (2010), we use particle metallicities, not SPH smoothed metallicities, and we neglect metal-line cooling (see, e.g., Schaye et al. (2015) for a discussion of the impact on metal line cooling on the calibration of SFR functions). Chemical enrichment thus impacts our simulations solely by setting the metallicity-dependent luminosities of the star particles (Schaerer, 2003) and the duration of the various phases of stellar evolution, determining the mass loss from stellar winds and the timing of SN explosions (Wiersma et al., 2009).
2.5 Ionizing luminosities
We use the Schaerer (2003) models to compute ionizing luminosities for the star particles in our simulations and adopt a Chabrier IMF in the range , consistent with our implementation of star formation. Schaerer (2003) does not tabulate models for the adopted IMF, but only for a Salpeter (1955) IMF in the range . We therefore correct the tabulated luminosities by a factor 0.4 to account for the difference in the mass range and by a factor 1.7 to account for the difference in the shape between the tabulated and target IMFs.
New with respect to PSD15 is that we compute metallicity-dependent luminosities, where we use the nearest available model metallicity larger than or equal to the metallicity of the gas particle from which the star particle was born. We approximately account for the pre-enrichment with metals by unresolved minihalos by adopting a metallicity floor (only) for the computation of the ionizing luminosities, (e.g., Bromm & Yoshida, 2011).
As in PSD15, we multiply the luminosities of star particles by the subresolution escape fraction to account for the removal of photons due to absorption in the unresolved interstellar medium (ISM) and for uncertainties in the intrinsic luminosities of the stars determined by the population synthesis model. New is that we calibrate the subresolution escape fraction separately for each simulation. The calibration ensures that all simulations achieve a volume-weighted hydrogen ionized fraction at nearly the same redshift (see Table 1). Note that the total escape fraction is the product , where is the fraction of photons absorbed by the gas particles in the galaxy hosting the star particles. The latter escape fraction is an outcome of the RT and is thus not calibrated. The calibration of the subresolution escape fraction will be discussed in more detail in Section 3 below.
2.6 Ionizing radiative transfer
We use the RT code TRAPHIC (Pawlik & Schaye, 2008; Pawlik & Schaye, 2011; Pawlik, Milosavljević, & Bromm, 2013; Raičević et al., 2014) to transport hydrogen and helium ionizing photons. TRAPHIC transports radiation directly on the spatially adaptive, unstructured computational grid traced by the SPH particles and hence exploits the full dynamic range of the hydrodynamic simulation. Figure 1 shows that this lets us capture small-scale structure in the ionized fraction in cosmological volumes.
TRAPHIC further employs directional averaging of the radiation field to render the computational cost of the RT independent of the number of sources. This lets us cope efficiently with the large number of ionizing sources typical of reionization simulations.
Photons are transported in the grey approximation using a single frequency bin. We compute grey photoionization cross sections , , using the fits to the frequency-dependent photoionization cross sections by Verner et al. (1996). The kinetic energy of the electron freed upon photoionizations of an ion of HI, HeI and HeII is , , and . All other parameters of the RT are as described in PSD15. In particular, photons leaving the box re-enter the box on the opposite side, that is, we are adopting periodic boundary conditions for the radiation. This lets us track the build-up of the UV background in a cosmological setting.
2.7 Stellar feedback
We account for the injection of mass and energy by both core-collapse and Type Ia SNe as well as the mass loss due to winds from massive stars and asymptotic giant branch (AGB) stars. This gives rise to a collective feedback on the properties of galaxies to which we refer as stellar feedback. Most of that feedback is due to the energy released by massive stars.
For each core-collapse SN that occurs, of thermal energy is stochastically injected among a subset of the neighboring SPH particles, using the method of Dalla Vecchia & Schaye (2012). The parameter accounts for radiative energy losses that our simulations may under- or overestimate because they lack both the resolution and the physics to accurately model the ISM. Our simulations also lack the resolution to distinguish between different types of energy feedback, such as SNe and radiatively driven winds. The factor can hence also be thought to account for the injection of energy by sources other than SNe. We improve on PSD15 by calibrating separately for each simulation to bring simulated SFR functions in even better agreement with each other and with observational constraints. The feedback efficiency parameter therefore depends on resolution (see Table 1). The calibration of the stellar feedback will be described below in Section 3 in more detail. New is also that those gas particles that have been heated by core-collapse SNe are immediately activated for the time integration so that they can react promptly to the change in energy, as described in Durier & Dalla Vecchia (2012).
Our numerical implementation of the age-dependent injection of energy by Type Ia SNe follows the deterministic method described in Schaye et al. (2010). For each Type Ia SN that occurs, of thermal energy is injected and distributed among all neighboring gas particles. Unlike for core-collapse SNe, we do not activate neighboring particles for the prompt response to the change in their energy. The less accurate treatment of Type Ia SNe is for computational convenience and is justified because they contribute little to the total energy released by a stellar population above , the range of redshifts relevant in this work (e.g., Wiersma et al., 2009). In future simulations, especially when simulating down to lower redshifts, it may make sense to adopt an implementation of Type Ia SNe that is more similar to that of Type II SNe.
Finally, the implementation of mass loss follows Wiersma et al. (2009).
2.8 Identification of galaxies
We calibrate the two main parameters of our simulations, the stellar feedback efficiency and the subresolution escape fraction of ionizing radiation , to yield agreement between simulated and observed SFR functions at independently of the resolution and to achieve a reionization redshift that is independent of the resolution and box size and in agreement with observational constraints on the optical depth for electron scattering from measurements of the CMB. The reionization redshift is thereby defined as the redshift at which the mean volume-weighted neutral hydrogen fraction reaches 0.5. Calibration is required, because cosmological simulations like ours lack both the physics and the resolution to predict SFRs and reionization histories from first principles (e.g., Rahmati et al., 2013b). A similar calibration was carried out in Gnedin (2014) and Schaye et al. (2015).
Each of these two numerical parameters can impact both the SFR functions and the redshift of reionization. This could potentially complicate the calibration, which might need to proceed iteratively, calibrating each parameter at a time, possibly converging slowly. However, reionization does not significantly impact the SFR function at the SFRs at which it is observationally constrained, it only affects the SFR function at much lower SFRs. We can therefore adopt the following two-step calibration procedure. First, we choose for each simulation a stellar feedback efficiency that yields agreement between observed and simulated SFRs. Second, we choose for each simulation a subresolution escape fraction such that the redshift of reionization, defined as the redshift at which half of the simulation volume is ionized, equals . The latter implies an optical depth to electron scattering in agreement with constraints from observations of the CMB (; Planck Collaboration et al., 2015).
The calibrated values of the feedback efficiency and the subresolution escape fraction will generally depend on box size and resolution, for several reasons. First, the radiative losses that are simulated explicitly may change if we better resolve the structure of the ISM. Keeping the feedback efficiency parameter fixed while changing the resolution may therefore affect the match between simulated and observed SFR function. Second, the reionization redshift is sensitive to resolution because resolution impacts the absorption of photons in the simulation, and it is sensitive to box size unless the box is much larger than the characteristic size of individual ionized regions. Third, a change in resolution impacts the SFR function at small, currently unobservable SFRs, at which the SFR function is not fixed by calibration. This, in turn, affects the redshift of reionization, also requiring us to adjust the subresolution escape fraction.
In the following, we discuss the calibration of each of the parameters in more detail. In addition to the calibrated simulations, we have carried out simulations in which we have varied the stellar feedback efficiency and the subresolution escape fraction away from their calibrated values, including simulations in which we have turned off SN explosions and photoheating altogether. While this work focuses on the calibrated set of simulations, these additional simulations can be used to gain further insight into the role of these parameters and the roles of galactic winds and photoheating, similar to our approach in PSD15.
3.1 Stellar feedback efficiency
Figure 2 shows the SFR functions at redshifts , 7, 8, and 9. We calibrated the stellar feedback efficiency at such that all simulations yield nearly identical SFR functions when considering galaxies with SFRs that are accessible by current observations using the HST, i.e., galaxies with SFRs . Where necessary, before carrying out the calibration, we have divided the observed SFRs, which assume a Salpeter IMF, by 1.7 to enable comparisons with our simulations, which assume a Chabrier IMF.
Figure 2 demonstrates generally good agreement between simulated and observationally inferred SFR functions above z = 6. By construction, our results are very close to the observational constraints by Smit et al. (2012) at z = 7, which we used for our calibration. However, differences with respect to more recent estimates of the SFR function and at other redshifts are minor.
The good match was facilitated by our calibration of the stellar feedback efficiency parameter, since stellar feedback dominates over feedback from photoheating in the range of SFRs at which the SFR function has been observationally inferred (PSD15). Achieving this match required us to reduce the stellar feedback efficiency parameter from at the lowest resolution to at the highest resolution (see Table 1). This implies that the cooling losses are smaller at higher resolution, which is consistent with Dalla Vecchia & Schaye (2012). The calibration does not depend on box size because for our set of simulations, box size impacts the SFR functions only at very large SFRs at which stellar feedback is inefficient. At these large SFRs, additional processes such as, e.g., dust corrections and black hole feedback may need to be invoked to reduce any remaining mismatch with the observations.
3.2 Escape fraction
The left panel of Figure 3 shows the evolution of the mean volume-weighted hydrogen neutral fraction. Because of our calibration of the subresolution escape fraction, , all simulations reach a neutral hydrogen fraction of 0.5 at , independently of the resolution (except for the dotted curve which shows a simulation with a different calibration using a lower subresolution escape fraction of 0.3, which is explained in Section 4.3). Furthermore, as the right panel of Figure 3 shows, all simulated reionization histories are in excellent agreement with the observational estimate of the electron scattering optical depth (Planck Collaboration et al., 2015).
Achieving this agreement required us to decrease the subresolution escape fraction with increasing resolution, from at the lowest resolution to at the highest resolution (see Table 1). This points at an increased contribution of low-mass galaxies as well as a more leaky gas distribution, facilitating the escape of photons, at increased resolution. The former can be deducted from the increasing star formation rate density with increasing resolution (see Figure 5). However, box size does not impact the choice of escape fraction, because in our simulations, the reionization histories do not change with box size above (please see Section 3.1 in PSD15 for a more detailed discussion).
Simulations of reionization can also be calibrated against the optical depth (e.g., Ciardi et al., 2012), but that does not necessarily guarantee the same redshift of reionization. Direct calibration against the redshift of reionization, as done here, still results in an excellent match between simulated and observationally inferred optical depths, while also rendering feedback from reionization insensitive to changes in resolution and box size. Because of this, a fixed redshift of reionization facilitates the analysis of the impact of resolution and box size on our simulations.
Figure 3 shows that in our calibrated set of simulations, reionization occurs earlier than inferred from observations of Lyman- emitters and quasar damping wings at . However, as we discuss in Section 4.3 below, the observational constraints are strongly model-dependent, and therefore matching these constraints was not our primary aim. Instead, the goal of the calibration was to yield a reionization redshift that is approximately independent of resolution, that implies an electron optical depth consistent with observations, and that is sufficiently high such that we can investigate the impact of the feedback from photoheating at redshifts .
The subresolution escape fraction merely multiplies the stellar luminosities and thus does not fully specify the escape fraction of galaxies in our simulations. The fraction of ionizing photons that escape a galaxy and are available to reionize the IGM is the product , where is the fraction of photons absorbed by simulation gas particles on their way to the IGM. Because the subresolution escape fraction is not a physical quantity, we are free to set it to values larger than unity, as we did in our low resolution simulations, as long as the number of ionizing photons reaching the IGM remains consistent with expectations from population synthesis models. A value then implies that the RT simulation overestimates the fraction of photons that are absorbed in the ISM, for example because it does not resolve holes in the ISM, or that we need to compensate for missing radiation from unresolved low-mass galaxies. Finally, we note that the calibrated subresolution escape fraction is chosen independently of the properties of the galaxies and time, while in reality the escape fraction may differ strongly between galaxies and change with time (e.g., Wise & Cen, 2009; Paardekooper, Khochfar, & Dalla Vecchia, 2015).
Here we discuss results from a first analysis of the simulations in Table 1. Because simulation L100N1024 was stopped at , we will show the SFR function at derived from this simulation, but otherwise omit it from our analysis.
Figure 4 shows images of the neutral hydrogen fractions, gas temperatures and gas metallicities in our reference simulation L025N0512 at three characteristic redshifts , 8, and 6, corresponding to volume-weighted ionized fractions of , 0.6, and 1, respectively. Stellar ionizing radiation creates individual ionized regions that grow and overlap. While the mean ionized fraction of these regions is very high, the neutral fraction can locally reach high values as gas is shadowed or self-shields from the ionizing radiation (see also Figure 1).
The reionization of the gas is accompanied by an increase in its temperature to about , the characteristic temperature of gas photoionized by stellar radiation. Locally the gas reaches higher temperatures, up to , where it is heated by structure formation shocks and explosive stellar feedback. Galactic winds enrich the gas with metals synthesized in the stars.
Photoheating raises the pressure in the galaxies and the surrounding reionized regions, which impedes the gravitational collapse of the gas in galaxies and in the IGM. As a consequence, reionization is associated with a decrease in the rates at which stars are formed and thus provides a negative feedback on reionization (e.g., Shapiro, Giroux, & Babul, 1994; Gnedin, 2000; Pawlik & Schaye, 2009; Noh & McQuinn, 2014; PSD15). In the IGM, however, the dynamical impact from photoheating reduces the recombination rate, and this provides a positive feedback on reionization (e.g., Wise & Abel, 2005; Pawlik, Schaye, & van Scherpenzeel, 2009; Sobacchi & Mesinger, 2014). Both types of feedback, negative and positive, are captured self-consistently by our radiation-hydrodynamical simulations.
In the following, we will provide a more quantitative discussion of our simulations, and investigate the dependence on box size and resolution.
4.2 Star formation history
Figure 2 shows that thanks in part to the calibration of stellar feedback near , the simulated and observed SFR functions at match very well near SFRs . At lower and higher SFRs, the calibration is not as effective and the simulated SFRs depend still strongly on resolution and box size.
Generally, higher resolution enables the formation of lower-mass galaxies with lower SFRs, and also renders stellar feedback more effective. As a result, at low SFRs, despite calibration, the simulations overpredict the SFR function, and this effect is strongest at the lowest simulated redshift and for the simulations with the lowest resolution. On the other hand, a larger box size reduces the bias on the estimate of the abundance of massive galaxies that results from missing large-scale power in the initial conditions. This impacts the comparison at high SFRs, where the abundance of star-forming galaxies is also affected by cosmic variance. The comparison of the simulated and observed SFR functions at large SFRs may also be affected by uncertain dust corrections and by feedback from massive black holes (e.g., Feng et al., 2016), both of which our simulations ignore.
Negative feedback from photoheating reduces star formation in low-mass galaxies, as discussed above. One may expect this to imprint an observable change in the slope of the SFR function at low SFRs, which correspond to low-mass halos (e.g., Barkana & Loeb, 2006; Muñoz & Loeb, 2011; Mashian & Loeb, 2013; Wise et al., 2014). However, strong stellar feedback may mask the impact of photoheating on the SFR function, because it also strongly reduces the normalization of the SFR function (PSD15; see Gnedin & Kaurov, 2014 for an alternative mechanism that reduces the signature of photoheating). In our simulations, the turnover in the slope of the SFR function at low SFRs is mostly a numerical artefact caused by the lack of resolution and low-temperature coolants, such as molecular hydrogen.
Figure 5 shows that the agreement between simulated and observed SFR functions carries over to the comparison between the evolution of the simulated SFR density contributed by observationally accessible galaxies with SFRs and the corresponding observational constraints. The figure also shows explicitly that our simulations predict a substantial population of faint star-forming galaxies that are currently below the detection threshold. These will be prime targets for upcoming observations with, e.g., the JWST (e.g., Johnson et al., 2009; Pawlik, Milosavljević, & Bromm, 2011; Zackrisson et al., 2011; Wise et al., 2014).
4.3 Reionization history
The left panel of Figure 3 shows that all simulations yield mean neutral hydrogen fractions in agreement with estimates from observations of the Lyman- forest at (filled triangles). Our calibrations of the SFR function and of the subresolution escape fraction implies that the gas is, on average, ionized earlier than inferred from observations of Lyman- emitters (e.g., Tilvi et al., 2014, Pentericci et al., 2014, Schenker et al., 2014; see Robertson et al., 2015 for an overview) and quasar damping wings (Schroeder, Mesinger, & Haiman, 2013) at . However, the interpretation of these observations is strongly model-dependent which makes their constraining power unclear (e.g., Dijkstra, 2014; Gnedin & Kaurov, 2014). Note that at higher resolution, reionization is more extended. This is most likely because at higher resolution, the abundance of low-luminosity sources, as well as that of absorbers of ionizing radiation, is increased. The right panel of Figure 6 shows that the adopted calibration of the subresolution escape fraction implies a reionization optical depth in excellent agreement with observational constraints (Planck Collaboration et al., 2015).
One should keep in mind that our simulations do not resolve any ionization by an early population of star-forming minihaloes (e.g., Abel, Bryan, & Norman, 2002; Bromm, Coppi, & Larson, 2002; Greif et al., 2008; Wise et al., 2014) and also do not account for ionization by non-stellar sources, such as, e.g., accreting black holes (e.g., Shull & Venkatesan, 2008; Haiman, 2011; Jeon et al., 2014a). A significant contribution from these sources to the total amount of ionizing radiation during reionization will require us to reduce the subresolution escape fraction in order to retain the match between simulated and observed reionization optical depths (e.g., Ahn et al., 2012; Wise et al., 2014). This might delay reionization as these very low-mass sources are unlikely to continue to contribute substantially during the later stages of reionization, and thereby improve the match between simulated and observed neutral hydrogen fractions at . For reference, a simulation at the resolution of our reference simulation but with subresolution escape fraction matches most observational constraints on the neutral fraction, while retaining consistency with the observed electron scattering optical depth (dotted curve in Figure 3; for computational efficiency, this simulation was carried out in a smaller box of size ).
4.4 Build-up of the ionizing background
The left panel of Figure 6 shows the build-up of the background of hydrogen ionizing radiation that accompanies the reionization of the gas. Immediately apparent is that the simulated volume-weighted mean photoionization rates at are substantially larger than the photoionization rate inferred from current observations (e.g., Wyithe & Bolton, 2011; Calverley et al., 2011). This could indicate that our adopted subresolution escape fractions are too high and result in reionization occuring too early. Indeed, a simulation at the resolution of our reference simulation but with a lower subresolution escape fraction yields a mean photoionization rate in agreement with the observations, while also matching the constraints on the neutral fraction and the optical depth (dashed curves in Figure 3).
However, additional factors may contribute to produce a high mean photoionization rate, and these would need to be considered upon any attempts to recalibrate the subresolution escape fraction. For instance, in our simulations, there is a large scatter around the mean photoionization rate at fixed density. This scatter causes the mean photoionization rate to be much higher than the median photoionization rate (see the long-dashed curve in the left panel of Figure 6), and this median might be more appropriate for comparison with the observational estimates (Rahmati et al. in prep.). Additionally, Figure 6 reveals a marked trend of decreasing mean photoionization rates at increasing resolution near the end of reionization, which is not yet fully converged in our reference simulation. If this trend continues to higher resolution, this may further help bring simulated and observed photoionization rates in agreement.
The trend of decreasing mean photoionization rates with increasing resolution is interesting in itself. Thanks to the calibration, the redshift of reionization in our simulations is independent of resolution. If the redshift of reionization is mostly set by the total amount of ionizing radiation escaping into the IGM, which our calibration controls, the decrease with increasing resolution in the mean photoionization rate, which scales with the total amount of ionizing radiation escaping into the IGM and the mean free path, may be caused primarily by the decrease in the mean free path of ionizing photons. The mean free path is set by the abundance of neutral hydrogen absorption systems with column densities (e.g., Prochaska, Worseck, & O’Meara, 2009). The right panel of figure 6 shows the column density distribution of neutral hydrogen computed as described in Rahmati et al. (2013a). The abundance of absorbers with column densities grows significantly with increasing resolution between our low-resolution and reference simulation. However, other factors may play a role, such as, e.g., the increased duration of reionization at higher resolution seen in the left panel of Figure 3. We will analyze this in more detail in subsequent work (Rahmati et al. in prep.).
4.5 Chemical enrichment history
The galaxies that reionize the IGM also enrich the cosmic gas with metals due to mass loss from stellar winds and explosions as SNe. Constraints on the metal enrichment history may therefore help in understanding how galaxies reionized the universe (e.g., Oh, 2002; Keating et al., 2014; Finlator et al., 2015; Ferrara, 2015). Note that the outflow of metals into the IGM can be facilitated by a prior episode of photoheating that reduces the densities near the sites at which the metals are produced (e.g., Pawlik & Schaye, 2009; Finlator, Davé, Özel, 2011; Jeon et al., 2014b; Walch & Naab, 2015).
Figure 7 shows the mean mass-weighted metallicity, defined as the ratio of metal mass and total mass, of all gas and of gas in the ISM only (i.e., gas with ). The mean total gas metallicity is substantially smaller than the metallicity of the ISM, consistent with the small fraction of the volume filled with metals (Figure 4). The mean metallicity of the star particles (not shown here) is very similar to that of the gas in the ISM. At , the mean mass-weighted gas metallicity is nearly converged with resolution in our reference simulation, reaching values , consistent with previous simulations (e.g., Wiersma et al., 2009; Salvaterra, Ferrara, & Dayal, 2011; Pallottini et al., 2014).
4.6 Galaxy properties
Figure 8 shows select properties of the population of galaxies at and explores their dependence on resolution in our calibrated set of reionization simulations. The baryon fraction (top left panel) is strongly reduced by photoheating and stellar feedback. Stellar feedback suppresses star formation significantly in galaxies as massive as , consistent with previous works (e.g., Finlator, Davé, Özel, 2011; see Fig. 6 in PSD15 for an explicit demonstration). At the reference resolution, stellar feedback is efficient only in galaxies more massive than . Because photoheating only removes substantial amounts of gas from galaxies , consistent with previous works (e.g., Okamoto, Gao, & Theuns, 2008; Pawlik & Schaye, 2009; Finlator, Davé, Özel, 2011; Gnedin & Kaurov, 2014), the baryon fraction near is artificially increased (see Figure 6 in PSD15 for an explicit demonstration). In the high resolution simulation, stellar feedback is efficient down to smaller masses, and the baryon fraction is strongly suppressed across the range of masses.
The SFRs of the simulated galaxies (top right and bottom left panels) are consistent with current observational estimates (McLure et al., 2011). The simulations predict that the currently observable galaxies with SFRs are hosted by halos with virial masses . Stellar feedback determines the overall normalization of the relation between SFR and virial mass in halos with virial masses , while photoheating and the lack of resolution and low temperature physics (such as molecular cooling) are mainly responsible for the strong suppression of the SFRs in halos with virial masses (see Figure 6 in PSD15 for an explicit demonstration).
By , the most massive galaxies with stellar masses have metallicities (bottom right panel). The relationship between stellar metallicities and masses of simulated galaxies is consistent with the stellar mass-metallicity relation of Local Group dwarf galaxies (Kirby et al., 2013). This is interesting because those galaxies are believed to have formed most of their stars at high redshifts (Weisz et al., 2014).
We have introduced Aurora, a new suite of cosmological radiation-hydrodynamical simulations of galaxy formation and reionization. Aurora uses accurate radiative transfer (RT) at the native, spatially adaptive resolution of the hydrodynamic simulation. The Aurora simulations track the stellar feedback from the explosion of stars as supernovae (SNe) and also account for the enrichment of the universe with metals synthesized in the stars. The reference simulation, which makes use of a box of size and contains dark matter and baryonic particles and therefore resolves atomically cooling haloes with virial temperatures with DM particles, is accompanied by simulations in larger and smaller boxes, with, respectively, lower and higher resolution to investigate numerical convergence. While our most demanding simulation, which uses particles in a box, was stopped at , all other simulations were run down to .
The Aurora simulations are designed to yield star formation rate (SFR) functions in agreement with observations of galaxies at and to complete reionization at , consistent with measurements of the electron scattering optical depth, independently of the box size and the resolution. We accomplished this by calibrating the stellar feedback efficiency, which is the main subresolution parameter impacting the SFR function over the range of observed SFRs, and the subresolution escape fraction, which is the main subresolution parameter impacting the redshift of reionization. The calibration is motivated by the fact that current cosmological simulations lack both the physics and the resolution to reliably simulate the structure of the ISM in galaxies. The latter determines how much of the energy injected by massive stars drives galactic outflows reducing the SFR, and how many of the emitted ionizing photons escape into the IGM.
While the main aim of this work was to describe the design and calibration of the simulations, we have also reported results from an initial analysis. We found that among Aurora simulations and after reionization, there is a strong decrease of the mean photoionization rate with increasing resolution. This coincides with an increase in the abundance of small-scale optically thick HI absorbers in the IGM. Given the large scatter in the photoionization rates in the Aurora simulations, the predicted IGM photoionization rates are consistent with the existing observational constraints at . The fact that the radiative transfer technique used to carry out Aurora is spatially adaptive, was critical for enabling us resolve the absorption of ionizing radiation on small scales in radiation-hydrodynamical simulations of cosmological volumes. In future works, we will exploit this unique feature of the Aurora simulations to study different properties of small scale absorbers and radiation field, and their evolutions.
We are grateful to Volker Springel for letting us use GADGET and the halo finder Subfind. We are grateful to Volker Bromm for letting us use his chemistry solver. We thank Milan Raicevic for his contributions in an early, preparatory phase of this project, and we thank Gunjan Bansal for helpful discussions. AHP thanks the Max Planck Institute for Astrophysics (MPA) for its hospitality during the work on this paper. Computer resources for this project have been provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre under grant:pr83le (PI Pawlik). We further acknowledge PRACE for awarding us access to resource SuperMUC based in Germany at LRZ Garching (proposal number 2013091919, PI Pawlik). Some of the simulations presented here were run on Odin at the Rechenzentrum Garching (RZG) and the Max Planck Institute for Astrophysics (MPA) and on Hydra at the Leibniz-Rechenzentrum (LRZ) in Garching. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s613 (PI Rahmati). This work was sponsored with financial support from the Netherlands Organization for Scientific Research (NWO). We also benefited from funding from the European Research Council under the European Unions Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement 278594-GasAroundGalaxies. AHP received funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) under grant agreement number 301096-proFeSsoR.
- Abel, Bryan, & Norman (2002) Abel T., Bryan G. L., Norman M. L., 2002, Sci, 295, 93
- Ahn et al. (2012) Ahn K., Iliev I. T., Shapiro P. R., Mellema G., Koda J., Mao Y., 2012, ApJ, 756, L16
- Alvarez, Finlator, & Trenti (2012) Alvarez M. A., Finlator K., Trenti M., 2012, ApJ, 759, L38
- Aubert, Deparis, & Ocvirk (2015) Aubert D., Deparis N., Ocvirk P., 2015, MNRAS, 454, 1012
- Barkana & Loeb (2001) Barkana R., Loeb A., 2001, PhR, 349, 125
- Barkana & Loeb (2006) Barkana R., Loeb A., 2006, MNRAS, 371, 395
- 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 et al. (2001) Becker R. H., et al., 2001, AJ, 122, 2850
- Bouwens et al. (2014) Bouwens R. J., et al., 2014, ApJ, 795, 126
- Bouwens et al. (2015) Bouwens R. J., et al., 2015, ApJ, 803, 34
- Bouwens (2015) Bouwens R. J., 2015, arXiv, arXiv:1511.01133
- Bromm, Coppi, & Larson (2002) Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23
- Bromm & Yoshida (2011) Bromm V., Yoshida N., 2011, ARA&A, 49, 373
- Calverley et al. (2011) Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543
- Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
- Ciardi, Stoehr, & White (2003) Ciardi B., Stoehr F., White S. D. M., 2003, MNRAS, 343, 1101
- Ciardi et al. (2012) Ciardi B., Bolton J. S., Maselli A., Graziani L., 2012, MNRAS, 423, 558
- Coe et al. (2013) Coe D., et al., 2013, ApJ, 762, 32
- Dale (2015) Dale J. E., 2015, NewAR, 68, 1
- Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
- Dijkstra (2014) Dijkstra M., 2014, PASA, 31, e040
- Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
- Duncan et al. (2014) Duncan K., et al., 2014, MNRAS, 444, 2960
- Durier & Dalla Vecchia (2012) Durier F., Dalla Vecchia C., 2012, MNRAS, 419, 465
- Ellis et al. (2013) Ellis R. S., et al., 2013, ApJ, 763, L7
- Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
- Fan, Carilli, & Keating (2006) Fan X., Carilli C. L., Keating B., 2006, ARA&A, 44, 415
- Feng et al. (2016) Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, MNRAS, 455, 2778
- Ferrara (2015) Ferrara A., 2015, arXiv, arXiv:1511.01120
- Finkelstein et al. (2015) Finkelstein S. L., et al., 2015, ApJ, 810, 71
- Finkelstein et al. (2015) Finkelstein S. L., Dunlop J., Le Fevre O., Wilkins S., 2015, arXiv, arXiv:1512.04530
- Finlator, Davé, Özel (2011) Finlator K., Davé R., Özel F., 2011, ApJ, 743, 169
- Finlator (2012) Finlator K., 2012, arXiv, arXiv:1203.4862
- Finlator et al. (2012) Finlator K., Oh S. P., Özel F., Davé R., 2012, MNRAS, 427, 2464
- Finlator et al. (2015) Finlator K., Thompson R., Huang S., Davé R., Zackrisson E., Oppenheimer B. D., 2015, MNRAS, 447, 2526
- Fontanot, Cristiani, & Vanzella (2012) Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413
- Gnedin (2000) Gnedin N. Y., 2000, ApJ, 542, 535
- Gnedin (2014) Gnedin N. Y., 2014, ApJ, 793, 29
- Gnedin & Kaurov (2014) Gnedin N. Y., Kaurov A. A., 2014, ApJ, 793, 30
- Gnedin (2016) Gnedin N. Y., 2016, arXiv, arXiv:1601.05802
- Greif et al. (2008) Greif T. H., Johnson J. L., Klessen R. S., Bromm V., 2008, MNRAS, 387, 1021
- Haiman (2011) Haiman Z., 2011, Natur, 472, 47
- Hasegawa & Semelin (2013) Hasegawa K., Semelin B., 2013, MNRAS, 428, 154
- 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. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725
- Jeon et al. (2014a) Jeon M., Pawlik A. H., Bromm V., Milosavljević M., 2014, MNRAS, 440, 3778
- Jeon et al. (2014b) Jeon M., Pawlik A. H., Bromm V., Milosavljevic M., 2014, 2014, MNRAS, 444,
- Johnson et al. (2009) Johnson J. L., Greif T. H., Bromm V., Klessen R. S., Ippolito J., 2009, MNRAS, 399, 37
- Kaurov & Gnedin (2014) Kaurov A. A., Gnedin N. Y., 2014, ApJ, 787, 146
- Keating et al. (2014) Keating L. C., Haehnelt M. G., Becker G. D., Bolton J. S., 2014, MNRAS, 438, 1820
- Khakhaleva-Li & Gnedin (2016) Khakhaleva-Li Z., Gnedin N. Y., 2016, arXiv, arXiv:1601.00641
- Kennicutt (1998) Kennicutt R. C., Jr., 1998, ApJ, 498, 541
- Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102
- Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
- Kuhlen & Faucher-Giguère (2012) Kuhlen M., Faucher-Giguère C.-A., 2012, MNRAS, 423, 862
- Lidz et al. (2008) Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Hernquist L., 2008, ApJ, 680, 962
- Lidz (2015) Lidz A., 2015, arXiv, arXiv:1511.01188
- Loeb (2009) Loeb A., 2009, JCAP, 3, 22
- Loeb (2010) Loeb A., 2010, hdfs.book,
- Madau & Haardt (2015) Madau P., Haardt F., 2015, ApJ, 813, L8
- Maio et al. (2010) Maio U., Ciardi B., Dolag K., Tornatore L., Khochfar S., 2010, MNRAS, 407, 1003
- Mashian & Loeb (2013) Mashian N., Loeb A., 2013, JCAP, 12, 17
- Mashian, Oesch, & Loeb (2016) Mashian N., Oesch P. A., Loeb A., 2016, MNRAS, 455, 2101
- McGreer, Mesinger, & D’Odorico (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499
- McLure et al. (2011) McLure R. J., et al., 2011, MNRAS, 418, 2074
- McQuinn et al. (2007) McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007, MNRAS, 377, 1043
- McQuinn et al. (2008) McQuinn M., Lidz A., Zaldarriaga M., Hernquist L., Dutta S., 2008, MNRAS, 388, 1101
- McQuinn (2015) McQuinn M., 2015, arXiv, arXiv:1512.00086
- Mellema et al. (2013) Mellema G., et al., 2013, ExA, 36, 235
- Muñoz & Loeb (2011) Muñoz J. A., Loeb A., 2011, ApJ, 729, 99
- Noh & McQuinn (2014) Noh Y., McQuinn M., 2014, MNRAS, 444, 503
- Norman et al. (2015) Norman M. L., Reynolds D. R., So G. C., Harkness R. P., Wise J. H., 2015, ApJS, 216, 16
- Ocvirk et al. (2015) Ocvirk P., et al., 2015, arXiv, arXiv:1511.00011
- Oesch et al. (2013) Oesch P. A., et al., 2013, ApJ, 773, 75
- Oesch et al. (2014) Oesch P. A., et al., 2014, ApJ, 786, 108
- Okamoto, Gao, & Theuns (2008) Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920
- Ono et al. (2012) Ono Y., et al., 2012, ApJ, 744, 83
- Ouchi et al. (2010) Ouchi M., et al., 2010, ApJ, 723, 869
- Oh (2002) Oh S. P., 2002, MNRAS, 336, 1021
- Ota et al. (2008) Ota K., et al., 2008, ApJ, 677, 12
- Ota et al. (2014) Ota K., et al., 2014, ApJ, 792, 34
- Paardekooper, Khochfar, & Dalla Vecchia (2015) Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544
- Pallottini et al. (2014) Pallottini A., Ferrara A., Gallerani S., Salvadori S., D’Odorico V., 2014, MNRAS, 440, 2498
- Pawlik & Schaye (2008) Pawlik A. H., Schaye J., 2008, MNRAS, 389, 651
- Pawlik, Schaye, & van Scherpenzeel (2009) Pawlik A. H., Schaye J., van Scherpenzeel E., 2009, MNRAS, 394, 1812
- Pawlik & Schaye (2009) Pawlik A. H., Schaye J., 2009, MNRAS, 396, L46
- Pawlik & Schaye (2011) Pawlik A. H., Schaye J., 2011, MNRAS, 412, 1943
- Pawlik, Milosavljević, & Bromm (2011) Pawlik A. H., Milosavljević M., Bromm V., 2011, ApJ, 731, 54
- Pawlik, Milosavljević, & Bromm (2013) Pawlik A. H., Milosavljević M., Bromm V., 2013, ApJ, 767, 59
- Pawlik, Schaye, & Dalla Vecchia (2015) Pawlik A. H., Schaye J., Dalla Vecchia C., 2015, MNRAS, 451, 1586 (PSD15)
- Pentericci et al. (2014) Pentericci L., et al., 2014, ApJ, 793, 113
- Petkova & Springel (2011) Petkova M., Springel V., 2011, MNRAS, 412, 935
- Planck Collaboration et al. (2015) Planck Collaboration, 2015, arXiv, arXiv:1502.01589
- Prochaska, Worseck, & O’Meara (2009) Prochaska J. X., Worseck G., O’Meara J. M., 2009, ApJ, 705, L113
- Rahmati et al. (2013a) Rahmati A., Pawlik A. H., Raičevic M., Schaye J., 2013, MNRAS, 430, 2427
- Rahmati et al. (2013b) Rahmati A., Schaye J., Pawlik A. H., Raičevic M., 2013, MNRAS, 431, 2261
- Raičević et al. (2014) Raičević M., Pawlik A. H., Schaye J., Rahmati A., 2014, MNRAS, 437, 2816
- Ricotti (2010) Ricotti M., 2010, AdAst, 2010
- Robertson et al. (2013) Robertson B. E., et al., 2013, ApJ, 768, 71
- Robertson et al. (2015) Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19
- Saitoh & Makino (2009) Saitoh T. R., Makino J., 2009, ApJ, 697, L99
- Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
- Salvaterra, Ferrara, & Dayal (2011) Salvaterra R., Ferrara A., Dayal P., 2011, MNRAS, 414, 847
- Scannapieco et al. (2012) Scannapieco C., et al., 2012, MNRAS, 423, 1726
- Schaerer (2003) Schaerer D., 2003, A&A, 397, 527
- Schaye (2004) Schaye J., 2004, ApJ, 609, 667
- Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
- Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
- Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
- Schenker et al. (2014) Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014, ApJ, 795, 20
- Schroeder, Mesinger, & Haiman (2013) Schroeder J., Mesinger A., Haiman Z., 2013, MNRAS, 428, 3058
- Shapiro, Giroux, & Babul (1994) Shapiro P. R., Giroux M. L., Babul A., 1994, ApJ, 427, 25
- Sharma et al. (2015) Sharma M., Theuns T., Frenk C. S., Bower R. G., Crain R. A., Schaller M., Schaye J., 2015, arXiv, arXiv:1512.04537
- Shull & Venkatesan (2008) Shull J. M., Venkatesan A., 2008, ApJ, 685, 1
- Smit et al. (2012) Smit R., Bouwens R. J., Franx M., Illingworth G. D., Labbé I., Oesch P. A., van Dokkum P. G., 2012, ApJ, 756, 14
- Sobacchi & Mesinger (2014) Sobacchi E., Mesinger A., 2014, MNRAS, 440, 1662
- Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Springel (2010) Springel V., 2010, MNRAS, 401, 791
- Tilvi et al. (2014) Tilvi V., et al., 2014, ApJ, 794, 5
- Trac & Gnedin (2011) Trac H. Y., Gnedin N. Y., 2011, ASL, 4, 228
- Verner et al. (1996) Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487
- Volonteri & Gnedin (2009) Volonteri M., Gnedin N. Y., 2009, ApJ, 703, 2113
- Walch & Naab (2015) Walch S., Naab T., 2015, MNRAS, 451, 2757
- Weisz et al. (2014) Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 147
- Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574
- Wise & Abel (2005) Wise J. H., Abel T., 2005, ApJ, 629, 615
- Wise & Cen (2009) Wise J. H., Cen R., 2009, ApJ, 693, 984
- Wise et al. (2012) Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50
- 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
- Wisotzki et al. (2015) Wisotzki L., et al., 2015, arXiv, arXiv:1509.05143
- Wyithe & Bolton (2011) Wyithe J. S. B., Bolton J. S., 2011, MNRAS, 412, 1926
- Yajima, Choi, & Nagamine (2011) Yajima H., Choi J.-H., Nagamine K., 2011, MNRAS, 412, 411
- Zackrisson et al. (2011) Zackrisson E., Rydberg C.-E., Schaerer D., Östlin G., Tuli M., 2011, ApJ, 740, 13
- Zaroubi et al. (2012) Zaroubi S., et al., 2012, MNRAS, 425, 2964
- Zaroubi (2013) Zaroubi S., 2013, ASSL, 396, 45