Simulations of galaxy formation with radiative transfer

Simulations of galaxy formation with radiative transfer: Hydrogen reionisation and radiative feedback

Margarita Petkova and Volker Springel
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, 85748 Garching, Germany
Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany
Centre for Astronomy, Heidelberg University, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
E-mail: mpetkova@mpa-garching.mpg.deE-mail:
Accepted 2010 October 29; Received 2010 October 18

We carry out hydrodynamical simulations of galaxy formation that simultaneously follow radiative transfer of hydrogen-ionising photons, based on the optically-thin variable Eddington tensor approximation as implemented in the GADGET code. We consider only star-forming galaxies as sources and examine to what extent they can yield a reasonable reionisation history and thermal state of the intergalactic medium at redshifts around . This serves as an important benchmark for our self-consistent methodology to simulate galaxy formation and reionisation, and for future improvements through accounting of other sources and other wavelength ranges. We find that star formation alone is sufficient for reionising the Universe by redshift . For a suitable choice of the escape fraction and the heating efficiency, our models are approximately able to account at the same time for the one-point function and the power spectrum of the Lyman- forest. The radiation field has an important impact on the star formation rate density in our simulations and significantly lowers the gaseous and stellar fractions in low-mass dark matter halos. Our results thus directly demonstrate the importance of radiative feedback for galaxy formation. The spatial and temporal importance of this effect can be studied accurately with the modelling technique explored here, allowing more faithful simulations of galaxy formation.

galaxy formation – cosmic reionisation – intergalactic medium
pagerange: Simulations of galaxy formation with radiative transfer: Hydrogen reionisation and radiative feedbackReferencespubyear: 2010

1 Introduction

In the standard cosmological model, the primordial gas recombines and becomes neutral around redshift . However, the absence of Gunn-Peterson troughs in the spectra of high redshift quasars up to (White et al., 2003; Fan et al., 2006) suggests that hydrogen is highly ionised at low redshift. Thus, there must be a period in the history of the Universe when hydrogen became ionised again, but it is still an open question when the process of this cosmic reionisation started, how it proceeded in detail, and which sources of radiation were primarily responsible for it.

An important observational clue about reionisation is provided by the total electron-scattering optical depth to the last scattering surface of the cosmic microwave background (CMB), found to be by the latest WMAP7 data release (Larson et al., 2010). This points to an early start and possibly extended period of reionisation. Further observational information about the history of hydrogen reionisation can be inferred through different astrophysical phenomena. The Gunn-Peterson troughs in quasar spectra are sensitive to small trace amounts of neutral hydrogen during the late stages of reionisation (); Lyman- emitting galaxies probe the intermediate stages of reionisation (); the 21-cm line background and gamma ray bursts (GRBs) probe the early stages of reionisation, when the Universe was mostly neutral (); and finally, the cosmic microwave background (CMB) polarisation provides important data on the free electron column density integrated over a large range of redshifts (Alvarez et al., 2006). In the future, upcoming observations from new radio telescopes such as LOFAR promise to be able to map out the epoch of cosmic reionisation in unprecedented detail, making theoretical studies of this process especially timely.

Population III stars are commonly believed to be the first significant sources of photons for reionisation. They are very massive, short-lived stars and form around redshift , or even earlier (Gao et al., 2007). Another, in principle plausible candidate for causing reionisation, are quasars. Madau et al. (1999) compare luminosity and spacial density functions of quasars and star-forming galaxies and show that quasars alone are unlikely to be the dominant sources of photons for reionising the Universe. On the other hand other groups (Trac & Gnedin, 2009; Thomas et al., 2009) predict in recent work that quasars can produce enough ionising photons to be able to reionise the Universe.

In many scenarios, reionisation is assumed to be driven mainly by high-mass objects with mass (Ciardi et al., 2000). However, the dim but abundant low mass sources () may still strongly influence the way in which the ionised regions grow (Sokasian et al., 2003; Choudhury & Ferrara, 2007). These small galaxies are preferentially found in low-density environments along the cosmic web and may contain many metal-free Pop-III stars in the early Universe. It has been estimated that these low-mass halos account for about 80% of the total ionising power at (Choudhury & Ferrara, 2007). Disregarding them may lead to an overestimate of the number of photons required to reionise the Universe (Sokasian et al., 2003).

Numerous theoretical studies have begun to investigate in detail the characteristic scales and the topology of the reionisation process through the use of numerical simulations (e.g. Gnedin & Ostriker, 1997; Miralda-Escudé et al., 2000; Gnedin, 2000; Gnedin & Abel, 2001; Sokasian et al., 2001; Ciardi et al., 2003; Sokasian et al., 2004; Iliev et al., 2006; Zahn et al., 2007; Croft & Altay, 2008; Shin et al., 2008; Wise & Abel, 2008). However, due to the high computational cost and complexity of the radiative transfer problem, most simulations, with very few exceptions (Gnedin & Ostriker, 1997; Kohler et al., 2007; Shin et al., 2008; Wise & Abel, 2008), have treated reionisation through post-processing applied to static or separately evolved gas density fields. This neglects the fact that the radiation field may exert important feedback effects on galaxy formation itself (e.g. Iliev et al., 2005; Yoshida et al., 2007; Croft & Altay, 2008). It is therefore an important task to develop more accurate theoretical models based on self-consistent simulations where the ionisation field is evolved simultaneously with the growth of cosmic structures, a topic that we address in this work.

In the literature, predictions based on simulations for the onset of reionisation range from redshifts (Iliev et al., 2007; Wise & Abel, 2008) to (Norman et al., 1998; Abel et al., 2002). Reionisation then probably proceeds in an inhomogeneous and patchy fashion (Lidz et al., 2007; Iliev et al., 2007), reflecting the inhomogeneous density distribution of the large-scale structure. At first, many isolated ionised bubbles are formed. They then grow in size from during the early stages of reionisation up to in the late phases (Furlanetto et al., 2006; Iliev et al., 2006). Around redshift (Iliev et al., 2006), (Lee et al., 2008), or (Gnedin & Fan, 2006), the ionised regions overlap. These values are very uncertain and depend strongly on the modelling details of reionisation and the parameters of the underlying galaxy formation simulations. However, it is plausible that the future use of more self-consistent simulation techniques should be able to reduce the systematic modelling uncertainties.

In this work, we therefore use the new radiative transfer algorithm we developed in a previous study (Petkova & Springel, 2009) to carry out high-resolution simulations of cosmic structure growth in the proper cosmological context. The approximation to radiative transfer employed, the optically-thin variable Eddington tensor approach (Gnedin & Abel, 2001), is fast enough to allow coupled radiative-hydrodynamic simulations of the galaxy formation process. At the same time, the employed moment-based approximation to the radiative transfer problem can be expected to be still reasonably accurate for the reionisation problem. In particular, thanks to the photon-conserving character of our implementation of radiative transfer and of the chemical network, ionisation fronts are bound to propagate with the right speed. The Lagrangian smoothed particle approach (SPH) we use automatically adapts to the large dynamic range in density developing in the galaxy formation problem. Combined with the fully adaptive gravitational force solver implemented in GADGET, this yields a numerical scheme that is particularly well suited for the cosmic structure formation problem.

For a first assessment of our new approach we study simulations of the standard CDM cosmology and treat star formation and supernova feedback with the ISM sub-resolution model developed in Springel & Hernquist (2003a). For simplicity, we shall here only consider ordinary star-formation regions as sources of ionising radiation. We are especially interested in whether the star formation predicted by the simulations results in a plausible reionisation history of the Universe, and whether it at the same time yields a thermal and ionisation state of the IGM at intermediate redshifts that is consistent with that probed by observations of the Lyman- forest. Finally, we are interested in possible differences induced in galaxy formation due to the spatially varying radiative feedback in the radiative transfer simulations, especially in comparison with the much simpler and so far widely adopted treatment where a spatially homogeneous UV background is externally imposed.

This paper is structured as follows. We start in Section 2 with a brief summary our methods for simulating hydrogen reionisation. We then present our results in Section 3, focusing on the history of reionisation in Section 3.1, the Lyman- forest in Section 3.2 and the feedback from reionisation in Section 3.3. We end with a discussion and our conclusions in Section 4.

2 Simulating hydrogen reionisation

2.1 Radiative transfer modelling

For our work we use an updated version of the cosmological simulation code GADGET (Springel et al., 2001; Springel, 2005), combined with the radiative transfer (RT) implementation of Petkova & Springel (2009). The RT equation is solved using a moment-based approach similar to the one proposed by Gnedin & Abel (2001). The resulting partial differential equation essentially describes an anisotropic diffusion of the photon density field ,111Differently from Petkova & Springel (2009) we solve equation (1) for the photon density , rather than the photon over-density, with respect to the hydrogen density .


where is the speed of light, is the absorption coefficient, is the Eddington tensor and is the source function. The closure relation for this particular moment-based method is obtained by approximating the Eddington tensor as


where is the radiation pressure tensor


This estimate of the Eddington tensors is carried out in the optically thin regime, giving the method its name.

Figure 1: Ionised fraction in a slice through the simulated volume at evolution time of the cosmological density field test described in Iliev et al. (2006)

The abundances of hydrogen and helium species are updated by accounting for the processes of photo-ionisation, collisional ionisation and recombination:


where all abundances of hydrogen species (, ) are expressed in dimensionless form relative to the total number density of hydrogen, all helium abundances are fractions with respect to the helium number density, and the electron abundance is expressed as . Furthermore, denotes the collisional ionisation coefficient, is the recombination coefficient and is the photo-ionisation cross-section for hydrogen at . These rate equations are solved using a semi-implicit scheme.

Figure 2: Volume averaged neutral fraction as a function of redshift for the low resolution simulations. The different colours represent simulation results for different values of the efficiency parameter (or ‘ISM escape fraction’) . Reionisation happens earlier for higher efficiency since then more photons become available for ionising the gas. In all cases, the final phase of reionisation proceeds rapidly; over a small range of redshift, the neutral volume fraction drops from 10% to negligibly small values.
Figure 3: Ionising background for the low resolution simulation as a function of redshift for different efficiency . As expected, the background becomes higher for a larger efficiency since more of the photons emitted by the sources are made available to build up the ionising background. The thin solid line shows the background computed by Haardt & Madau (1996) for quasars, which provides for an interesting comparison. Our results for the time evolution of the neutral volume fraction agree quite well with previous studies (e.g. Gnedin, 2000).

The photo-heating of the gas by the ionising part of the spectrum is described by the heating rate




is a frequency averaged photo-ionisation cross-section and


is a frequency averaged photon excess energy (Spitzer, 1998). For our calculations, we usually assume a black body spectrum with , which leads to and . However, given the uncertainties in accurately calculating the net energy input during the photo-ionisation transition, which involves non-equilibrium processes, we will later on vary in order to investigate the influence of different heating efficiencies.

We have tested our radiative transfer implementation on the static cosmological density field that was used in the radiative transfer code comparison study by Iliev et al. (2006). In Figure 1 we show the ionised fraction in a slice through the simulated volume at evolution time . Reassuringly, our result is in good agreement with the ones obtained by other radiative codes in the comparison study.

2.2 Treatment of star formation and radiative cooling

We follow the radiative cooling of helium and hydrogen with a non-equilibrium treatment, with the rates as cited in Petkova & Springel (2009). Winds and star formation in the dense, cold gas is modelled in a sub-resolution fashion (Springel & Hernquist, 2003a), where the interstellar medium is pictured as being composed of cold clouds that are embedded at pressure equilibrium in a hot tenuous phase that is heated by supernova explosions. Through the evaporation of clouds, this establishes a tight self-regulation of the star formation rate. Previous work has shown that this model converges well with numerical resolution and yields star formation rates that are consistent with the Kennicutt relation (Kennicutt, 1998) observed at low redshift. We shall here assume that the same star formation law also holds at high redshift.

We use the star-forming regions in all simulated galaxies as sources of ionising photons in our radiative transfer model. We adopt an ionising source luminosity of (Madau et al., 1999), which relates the number of emitted photons to the star formation rate in units of . This source luminosity is released individually by every star-forming gas particle, hence the number of numerically represented sources is a non-negligible fraction of all simulation particles. Fortunately, the speed of our radiative transfer algorithm is almost insensitive to the total number of sources, because the Eddington tensor calculation can be carried out with a tree algorithm similar to the gravity calculation. This insensitivity of the computational cost to the number of sources is a significant advantage of the method used here, and is not shared by most alternatives for treating radiative transfer.

Figure 4: Slices through the middle of the simulated volume for the low resolution simulation realisation with and . The maps are showing the density (left) with contours at ionised fraction , and overlaid, ionised fraction (middle), and temperature (right) of the gas. The snapshots show a time evolution from top to bottom, with individual redshifts , , , and , respectively.

To account for the uncertain absorption that occurs in reality in the spatially unresolved multi-phase structure of our simulations, we impose a phenomenological efficiency factor on the source luminosities. In our simulation set we explore values in the range to get a feeling for the sensitivity of our results to this uncertainty. We note however that should not be confused with what is usually called galaxy escape fraction, which has a slightly different meaning. Our is meant to be just an ‘interstellar medium escape fraction’ whereas photon losses in the gaseous halos of galaxies will be taken into account self-consistently in our simulations.

Figure 5: Scatter plot of the neutral fraction as a function of over-density for two different redshift of the low resolution simulation realisation with and . The thick dashed line is the median neutral fraction. The neutral fraction shows a clear dependence with density, where high density regions are more ionised than low density regions. Very high density gas is highly ionised as well, because it is in the star forming stage.
Figure 6: Lyman- flux probability (left) and power spectrum (right) from the low resolution simulation with and at redshift . The result is compared with observational data from McDonald et al. (2000) and Kim et al. (2007). The flux probability agrees reasonably well with observations, whereas the power spectrum deviates at the high end. We discuss this problem in the text.
Figure 7: Lyman- flux probability (left) and power spectrum (right) for the low resolution simulation with efficiency , for four different values of the averaged excess photon energy . The simulated data is compared to the observational result from McDonald et al. (2000) and Kim et al. (2007).

2.3 Simulation set

All our simulations assume a CDM universe with cosmological parameters , , , and . In order to have sufficiently high mass resolution, we follow a comparatively small region in a periodic box of comoving size on a side. This region is still sufficiently large to give a representative account of the Lyman- forest at redshift , which is the final time of our runs.

Our primary simulation set has dark matter and gas particles, giving a mass resolution of and in dark matter and gas, respectively. A selected model was also carried out at the much higher resolution of , giving a gas mass resolution of . The gravitational softening was chosen as of the mean particle spacing in each case, corresponding to and in the two resolution sets. In order to save computational time, most runs were restarted from with different radiative transfer treatments, such that the higher redshift evolution did not have to be repeated. We have also systematically varied the ‘escape fraction’ and the heating efficiency , considering the values and for the low resolution set. For the high resolution run we chose the parameters and . In this way we could systematically determine the settings that give the most promising agreement with the observations.

3 Results

3.1 Hydrogen reionisation history

When star formation starts at around redshift in our simulations, the process of reionisation begins through photo-ionisation of the gas around the star-formation sites. However, as the dense gas has a high recombination rate, the progress in the reionisation is sensitively determined by a competition between the luminosity of the sources, the rate with which they turn on, and the density of the neutral gas they are embedded in. Using our simulation set, we first investigate the global reionisation history and its dependence on the source efficiency parameter .

Figure 2 shows the reionisation histories of our low resolution simulated box for several choices of . As expected, the results for the reionisation redshift depend strongly on . For larger efficiencies, the Universe gets reionised earlier, since more photons become available to ionise the hydrogen. This effect always overwhelms the reduction in star formation and hence source luminosity that an increased efficiency parameter also induces. In all cases the last phase of the Epoch of Reionisation (EoR) is rather short, i.e. the final 10% of the volume transition very rapidly from being neutral to being essentially fully ionised. The redshift of reionisation when this occurs varies between and for and , respectively. In contrast, the EoR starts in general at much higher redshift. For example, the highest efficiency model has already ionised nearly 30% of the volume by redshift . Interestingly, there is a systematic variation of the time it takes to complete the last phase of the EoR. This is much more rapid for the high efficiency model than for the lower efficiency ones. Since typical bubble sizes at redshift are up to , we note that the size of our simulation box is really too small to draw any definitive conclusions about the global duration of reionisation in the Universe. A number of authors (Ciardi et al., 2003; Furlanetto & Oh, 2005; Furlanetto & Mesinger, 2009) have pointed out that the local reionisation history depends on the environment, so that, for example, the evolution of the neutral fraction in a field or void region is different than in a proto-cluster region. As a result, only very large simulation volumes of or more on a side can be expected to yield a truly representative account of cosmic reionization. Since at the same time an equal or better mass resolution as we use here is required, such calculations are very expensive and beyond the scope of this work, but the should become feasible in the future on the next generation of high-performance computers.

In Figure 3 we compare the ionising background for the same simulations, which gives further interesting clues about the reionisation histories of the different models. The background is computed as the volume-averaged intensity in the simulation box. The ionising background is compared with an analytical estimate by Haardt & Madau (1996) for the meta-galactic ionising flux from quasars rather than stellar sources, which is an interesting comparison point for the expected value of the background. Clearly, lowering the efficiency decreases the ionising background as well, consistent with the findings above. Interestingly, the rapid rise of the mean background intensity ends when reionisation is complete. From this point on, the background shows only a weak residual evolution.

Figure 8: Volume averaged neutral fraction as a function of redshift. The set of low resolution simulations with efficiency and average photon excess energy is compared with the high resolution run with and . When a higher energy per ionising event is injected, the universe gets ionised slightly earlier since higher temperatures help to maintain higher ionised fractions.

The time evolution of the temperature, ionised fraction and density fields in a representative simulation model ( and ) is illustrated in Figure 4. The different panels correspond to slices through the middle of the simulation box, between redshifts and , from the top row to the bottom row. It is seen that the ionised regions start to grow first around high density peaks, where the star forming regions are concentrated. Then the radiation diffuses into the inter-cluster medium. The filaments remain less ionised than the voids for a while, since their density is much higher. Initially the photons heat up the gas in the ionised region to temperatures slightly above . As the ionising background declines due to the expansion of the Universe and the drop of star formation, the heating becomes less effective and the temperature of the highly ionised gas in the voids drops somewhat as a result of the expansion cooling.

Figure 9: Ionising background as a function of redshift. The set of low resolution simulations with efficiency and averaged photon excess energy is compared with the high resolution run with and . For a higher injected energy per ionisation event, the background also increases as a result of the higher temperature, which leaves more photons unabsorbed so that they can contribute to a higher level of the ionising background.

In Figure 5, we present a contour plot of the neutral fraction versus the over-density of the gas in the representative model for two different redshifts, at before reionisation is completed, and at after reionisation is completed. There is a clear dependence of the neutral fraction on over-density in both cases. The high density regions around star-forming matter are ionised very quickly. The average density regions, e.g. filaments, tend to be less ionised and get ionised after the low density regions. However, note that in the star-forming tale of the diagram all the gas is ionised. This is here due to the star formation scheme adopted in GADGET, where star-forming gas particles are assigned a mean mass-weighted temperature which is so high that all this gas is formally collisionally ionised. We note that these results are very similar to the ones reported by Gnedin (2000), where a similar relation between neutral fraction and density is found. However, we are able to probe somewhat higher densities thanks to better spatial resolution of our simulations.

3.2 Lyman- forest

We next turn to an analysis of the thermal state of the intergalactic medium left behind at by our self-consistent reionisation simulations. To this end we use the gas density, gas temperature, gas velocity and ionisation state of the gas in the simulation box and compute Lyman- absorption spectra for random lines of sight. We then compare the statistics of these artificial absorption spectra with observational data on the Lyman- forest at this epoch, as given by McDonald et al. (2000) and Kim et al. (2007). Here we assume that the main source of ionising sources is stellar and thus discard any contribution from a quasar-type spectra, in agreement with Madau et al. (1999). We also note that some discrepancies are possible due to the fact that helium is not photo-ionised in our simulations, but only collisionally ionised.

Our simulations have the necessary gas mass resolution at redshift required to reproduce realistic Lyman- absorption in the low density regions (Bolton & Becker, 2009). However, we can not match the required simulation volume of size to properly sample the largest voids. This can have a significant effect on our predictions of the Lyman- flux probability distribution and power spectra.

Figure 6 shows the flux probability distribution function (PDF) and flux power spectrum for our simulation model, where an efficiency parameter of and averaged photon excess energy were adopted. For this choice, we achieve the best fit to the flux PDF. However, for all the other models the power spectrum is over-predicted at high wave numbers. We suggest that this overestimation of the power spectrum is due to the insufficient heating of the gas in low density regions, causing an excess of small-scale structure in the Lyman- forest.

Figure 10: Star formation rate density as a function of redshift for the low resolution simulation set at different efficiencies of . The results are compared to the SFR history of a low resolution simulation with instantaneous reionisation at and photo-heating by a Haardt & Madau (1996) ionising background (thin black line), and a simulation with neither reionisation nor photo-heating (thick black line). The photo-heating from stellar sources decreases star formation, as suggested by Pawlik & Schaye (2009). As the escape efficiency gets higher, this effect becomes progressively stronger.

In order to examine this effect further, we vary the photon excess energy used in the photo-heating and examine the influence this has on the flux probability and power spectrum. There are two possible reasons why our simulations underestimate the photo heating. First, we expect that some non-equilibrium effects in the photo-heating are treated inaccurately due to our implicit treatment of the radiation transport and chemistry (e.g. Bolton & Becker, 2009). Second, photo-heating is different in optically thin and optically thick regions. For example, in an optically thick region the average photon excess energy obtained from Eqn. (11) is . It is however likely that our approximative radiative transfer scheme leads to inaccurcies in the effective heating rates of regions of different optical depths, due to the varying accuracy of the scheme in different regimes. Part of these inaccuracies can be absorbed into a suitably modified value of the effective heating rate . To explore the full range of plausible values, we therefore vary the values for as follows: and . We aim to bracket what can be expected when non-equilibrium effects are fully taken into account in future treatments, and want to identify the case that provides the best representation of the Universe at redshift .

Figure 6 shows the flux PDF and power spectra for these different heating values. Clearly, the high wave number region of the flux power spectrum is strongly influenced by the amount of injected heat energy into the gas, and the increase of the temperature also affects the flux probability distribution. For the low efficiency of , there is a substantial mismatch already in the flux PDF, simply because there is too little ionisation overall so that the mean transmission ends up being too low. However, as the adopted photo-heating energy increases, the gas is getting hotter and is able to stay ionised longer due to the higher temperatures, yielding a better fit to the flux PDF. At the same time, small-scale structure in the flux power spectrum is erased due to thermal broadening, bringing the simulations into agreement with the observation. This shows the power of detailed Lyman- data to constrain simulations of the reionisation process. In our current models we need to adopt a quite extreme heating efficiency of combined with a low ‘escape fraction’ of to achieve a good match to the data.

Figure 11: Star formation rate density as a function of redshift. The set of low resolution simulations with efficiency and averaged photon excess energy are compared to the high resolution run with and . The results are compared to the SFR history of a low resolution simulation with instantaneous reionisation at and photo-heating by a Haardt & Madau (1996) ionising background (thin black line), and a simulation with neither reionisation nor photo-heating (thick black line). The star formation decreases with increasing heating energy, as expected. For the low resolution run with , the result of the self-consistent radiative transfer calculation matches the simulation with instantaneous reionisation. The high resolution simulation SFR is higher at higher redshift due to better resolution and agrees well with the other results at redshifts less than .
Figure 12: Lyman- flux probability (left) and power spectrum (right) for the high resolution simulation with efficiency and averaged excess photon energy , compared to observational results from McDonald et al. (2000) and Kim et al. (2007).

In Figures 8 and 9 we compare the impact of the different photo-heating efficiencies on the evolution of the neutral volume fraction and the ionising background. We also show for comparison the results from our high resolution simulation, which is discussed below in the text. As expected, an increase in the heating energy leads to a slightly earlier reionisation and to a slightly elevated ionising background flux. Both of these effects can be readily understood from the higher gas temperature produced in the ionised gas when the higher heating efficiency is adopted. However, the effect is quite weak, and very much smaller than the changes resulting from a different choice of .

We have also measured the Thomson electron scattering optical depth in our high resolution simulation and found it to be , which is smaller than the WMAP7 value (Komatsu et al., 2010). This discrepancy is, however, not critical since the simulated volume is too small to obtain a realistic value and we have also not included photo-ionisation of helium.

3.3 Feedback from reionisation

Reionisation can in principle exert a strong feedback effect on the gas through the temperature increase induced by photo-heating. As the gas temperature increases, the gas densities will be lowered through pressure effects. The gas will then cool and collapse more slowly, such that the star formation rate is ultimately reduced. Especially small dark matter halos should be sensitive to this effect. In the extreme case of halos that have virial temperatures comparable to or only slightly larger than the temperature reached by the gas through reionisation, the UV radiation may even completely suppress atomic cooling and efficient star formation. This effect is often invoked to explain why so many of the dark matter satellites expected in CDM in the halos of ordinary galaxies are apparently largely devoid of stars.

In order to highlight the radiative transfer effects on the star formation in galaxies, we compare our simulations with two fiducial models where no radiative transfer is used. The first is a simulation simply without any photo-heating of the gas, while the second one is a simulation where reionisation is induced by an externally imposed, spatially homogeneous UV background based on a modified Haardt & Madau (1996) model that causes reionisation and an associated photo-heating of the gas at (for details see Davé et al., 1999). The latter model corresponds to the standard approach applied in many previous hydrodynamical simulation models of galaxy formation (e.g. Tornatore et al., 2003; Wadepuhl & Springel, 2010).

In Figure 10, we compare our results for the cosmic star formation rate density evolution as a function of the adopted efficiency parameter (for the low resolution simulation set). We also include the two fiducial comparison models as limiting cases. As we increase the escape efficiency of the ionising radiation, the star formation drops, as expected, since this makes more photons available to photo-heat the gas. We note that our results for the SFR are always lower than the fiducial simulation where no photo-heating is included at all, consistent with findings by Pawlik & Schaye (2009). Towards lower redshift, the reduction of the SFR due to the radiation field becomes progressively larger. The run with quite closely corresponds to the simulation with the imposed reionisation epoch, but starts to slightly differ at redshifts .

We also carry out a corresponding comparison for a low resolution simulation set with constant efficiency but different values for the photon excess energy. In Figure 11 we show the results for the SFR, again including the two fiducial models as limiting cases for comparison. The results confirm the expectation that an increase of the photon excess energy decreases the star formation rate density. Interestingly, the model that best reproduced the Lyman- power spectrum observations, the one with , quite closely follows the star formation rate density obtained for the fiducial model where reionisation is imposed at .

For our high resolution run we chose to repeat the simulation with averaged photon excess energy and adopt a higher escape fraction rather than . In this way we make sure we account for the trapping of photons in high density peaks, which were not present in the low resolution runs. In Figure 11 we show how the star formation rate history compares to the ones from the low resolution runs. They are in good agreement, except for the higher redshift, where the high resolution captures more star formation, as expected (Springel & Hernquist, 2003b). The Lyman- forest flux probability and power spectrum at for this simulation are shown in Figure 12. While the simulated data is in reasonable qualitative agreement with the observational results from McDonald et al. (2000) and Kim et al. (2007), it does not provide in this case a detailed fit within the error bars, again highlighting that simultaneously accounting for the cosmic star formation history, cosmic reionisation and the state of the IGM at intermediate redshifts provides a powerful constraint on self-consistent simulations of galaxy formation and reionisation.

Figure 13: Evolution of the mass averaged temperature at three different over-densities for the low resolution simulation with and , with and without photo-heating. The strongest effect is observed in the low density gas, which is heated by the photons much more than the higher density gas. At all densities, however, photo-heating increases the temperature, as expected.
Figure 14: Median temperature of the gas as a function of over-densities at redshifts and for the high resolution simulation with and . At redshift , shortly after reionisation is completed, the temperature at low densities is clearly higher than that at higher densities (apart from the gas in the star-forming phase). This can be interpreted as an inverted equation of state. At lower redshift the relation reverts again to normal form as the gas in the low density regions cools down adiabatically due to the expansion of the Universe.
Figure 15: Mean stellar and gas masses as a function of the DM halo mass at in the high resolution simulation. The black vertical corresponds to a mass of 100 DM particles, which can be taken as an (optimistic) resolution limit of the simulation. Photo-heating slows down the collapse of gas in halos, which in turn also decreases their stellar and gas masses. The effect becomes stronger for low mass DM halos.
Figure 16: Baryon fraction as a function of the DM halo mass at in the higher resolution simulation. The black vertical corresponds to a mass of 100 DM particles.
Figure 17: Baryon fraction as a function of DM halo mass at in the lower resolution simulations with efficiency , and for different averaged photon excess energy. The dashed line shows the baryon fraction when no photo-heating has taken place. The black vertical line corresponds to a mass of 100 DM particles. Photo-heating does not affect the baryon fraction as strongly as in the high resolution simulation. We also observe that differences in the excess photon energy do not have a large effect on the baryon fraction.

In Figure 13, we explore the temperature evolution of the gas at different characteristic densities, corresponding to under-dense gas by a factor of 10, gas at the mean density, and over-dense gas by a factor of 10 relative to the mean. We compare our default simulation with radiative transfer and photo-heating to the fiducial simulation where no such heating is included at all. Clearly, the effect of photo-heating is most prominent in the lowest density gas. This gas is only weakly heated by structure formation shock waves when photo-heating is not included. In contrast, when reionisation is accounted for, the temperature of this gas reaches a high value of at the end of the epoch of cosmic reionisation, and even before that, the mean temperature of this gas is raised considerably as a result of the patchy and temporally extended reionisation transition in our radiative transfer simulations. Interestingly, after reionisation is complete, the mean temperature of the under-dense gas starts to slowly decline again, while already for the mean density gas structure formation shocks can provide for a slow further increase of the temperature.

We also analysed the median temperature of the gas as a function of over-density. As shown in Figure 14, after reionisation has been completed, the low density gas ends up with a higher median temperature than the higher density gas (except for the gas in the star-forming phase). This points towards an ‘inverted equation of state’, as observed by Bolton et al. (2008), Trac et al. (2008) and Furlanetto & Oh (2009). At later times, the equation of state reverts again to a normal positive slope, when the low density gas cools down due to the adiabatic expansion of the Universe.

Finally, in Figure 15, Figure 16 and Figure 17 we explore the impact of the ionising radiation field on the gas and stellar mass content of individual dark matter halos. To this end we run a group finder on our simulations and simply determine the average gas mass, stellar mass and baryon fraction of halos as a function of their dark matter mass. We compare the results of our higher resolution radiative transfer simulation with the simulation where photo-heating is completely ignored. Interestingly, we find a reduction of the gas and stellar mass for all halo masses when radiative transfer is included. The effect is quite weak for large halos but becomes progressively larger for small halos. At dark matter halo masses of the suppression in baryonic content is approximately 60%, while at it drops to only a few percent. This shows clearly the important impact of the ionising radiation field on small dwarf galaxies, in particular. While an externally imposed UV background can perhaps account for the mean effect of this radiative feedback process (Hoeft et al., 2006; Okamoto et al., 2008), only a spatially resolved treatment of radiative transfer can account for effects of proximity that may well play an important role in shaping, e.g., the satellite luminosity function (Muñoz et al., 2009; Busha et al., 2010; Iliev et al., 2010).

4 Discussion and conclusions

We have presented the first application of our new implementation (Petkova & Springel, 2009) of radiation hydrodynamics in the cosmological simulation code GADGET. We focused on the problem of cosmic reionisation, aiming in particular at a first test on whether the default star formation model in the code combined with our radiative transfer modelling can yield a plausible reionisation history of the Universe and a reasonable thermal state of the intermediate redshift intergalactic medium. For simplicity, we have here only studied star-forming galaxies as ionising sources, and restricted the analysis to hydrogen reionisation alone. Based on the encouraging results collected here, it is clearly worthwhile to extend the model further in future work.

Since the level of internal absorption in the interstellar medium is uncertain and cannot be resolved by our simulations, we have examined models with different effective source efficiencies . Likewise, as we have not included a detailed spectral treatment and the time evolution of non-equilibrium in the chemistry may be inaccurate, we have parametrised the heat input per ionisation event in terms of a parameter .

We find that our simulated universes can get reionised by star formation in ordinary galaxies alone, with the epoch of reionisation ending between redshifts to , depending on the assumed escape efficiency. The final phase transition is always quite rapid in this our setup, but sizable fractions of the volume begin to be reionised much earlier. The heating efficiency has only a weak influence on the reionisation history, but a stronger one on the cosmic star formation rate density. In fact, we have shown that photo-heating plays an important role in the evolution of the baryonic gas. As a result of the associated heating, it changes baryonic structure formation by slowing down the collapse of gas in dark matter halos, thereby delaying star formation. This effect is strongest for the lowest mass halos, where the DM potential well is not deep enough to easily overcome the thermal pressure from the effects of photo-ionisation.

Our simple models of a self-consistent treatment of galaxy formation and radiative transfer are not only able to produce a plausible history of reionisation, but they also manage to approximately match the basic statistics of the Lyman- forest, at least for an appropriate choice of the parameters and . This suggests that the low-redshift IGM data can be a powerful additional constraint on future reionisation modelling in structure formation simulations.

Despite these encouraging results it is also clear that our simulation results are likely still affected by numerical resolution effects, because the resolution in the lowest mass halos is still too coarse to yield fully converged results. Ideally, we would like to resolve the full range of star-forming halos with enough particles to achieve fully converged results. While this is unlikely to qualitatively change any of the results presented here, future precision work will require such calculations. Another important caveat that will require further study are uncertainties due to the radiative transfer approximation itself. This is probably best addressed by comparing the results with a completely different approach to radiative transfer. We are presently finalising such a scheme in the moving-mesh code AREPO. It avoids the use of the diffusion approximation and can cast sharp shadows, hence a direct one-to-one comparison with the optically-thin variable Eddington tensor approach applied here will be particularly interesting.


The authors would like to thank Benedetta Ciardi for the helpful discussion and comments. This research was supported by the DFG cluster of excellence “Origin and Structure of the Universe” (


  • Abel et al. (2002) Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93
  • Alvarez et al. (2006) Alvarez M. A., Komatsu E., Doré O., Shapiro P. R., 2006, ApJ, 647, 840
  • Bolton & Becker (2009) Bolton J. S., Becker G. D., 2009, MNRAS, 398, L26
  • Bolton et al. (2008) Bolton J. S., Viel M., Kim T., Haehnelt M. G., Carswell R. F., 2008, MNRAS, 386, 1131
  • Busha et al. (2010) Busha M. T., et al., 2010, ApJ, 710, 408
  • Choudhury & Ferrara (2007) Choudhury T. R., Ferrara A., 2007, MNRAS, 380, L6
  • Ciardi et al. (2000) Ciardi B., Ferrara A., Governato F., Jenkins A., 2000, MNRAS, 314, 611
  • Ciardi et al. (2003) Ciardi B., Stoehr F., White S. D. M., 2003, MNRAS, 343, 1101
  • Croft & Altay (2008) Croft R. A. C., Altay G., 2008, MNRAS, 388, 1501
  • Davé et al. (1999) Davé R., Hernquist L., Katz N., Weinberg D. H., 1999, ApJ, 511, 521
  • Fan et al. (2006) Fan X., et al., 2006, AJ, 131, 1203
  • Furlanetto et al. (2006) Furlanetto S. R., McQuinn M., Hernquist L., 2006, MNRAS, 365, 115
  • Furlanetto & Mesinger (2009) Furlanetto S. R., Mesinger A., 2009, MNRAS, 394, 1667
  • Furlanetto & Oh (2005) Furlanetto S. R., Oh S. P., 2005, MNRAS, 363, 1031
  • Furlanetto & Oh (2009) Furlanetto S. R., Oh S. P., 2009, ApJ, 701, 94
  • Gao et al. (2007) Gao L., Yoshida N., Abel T., Frenk C. S., Jenkins A., Springel V., 2007, MNRAS, 378, 449
  • Gnedin (2000) Gnedin N. Y., 2000, ApJ, 535, 530
  • Gnedin & Abel (2001) Gnedin N. Y., Abel T., 2001, New Astronomy, 6, 437
  • Gnedin & Fan (2006) Gnedin N. Y., Fan X., 2006, ApJ, 648, 1
  • Gnedin & Ostriker (1997) Gnedin N. Y., Ostriker J. P., 1997, ApJ, 486, 581
  • Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
  • Hoeft et al. (2006) Hoeft M., Yepes G., Gottlöber S., Springel V., 2006, MNRAS, 371, 401
  • Iliev et al. (2006) Iliev I. T., Ciardi B., Alvarez M. A., Maselli A., Ferrara A., Gnedin N. Y., Mellema G., Nakamoto T. an d Norman M. L., Razoumov A. O., Rijkhorst E.-J., Ritzerveld J., Shapiro P. R., Susa H., Umemura M., Whalen D. J., 2006, MNRAS, 371, 1057
  • 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. (2010) Iliev I. T., Moore B., Gottloeber S., Yepes G., Hoffman Y., Mellema G., 2010, ArXiv e-prints: 1005.3139
  • Iliev et al. (2007) Iliev I. T., Shapiro P. R., Mellema G., Pen U.-L., McDonald P., Alvarez M. A., 2007, ArXiv e-prints: 0710.2451
  • Iliev et al. (2005) Iliev I. T., Shapiro P. R., Raga A. C., 2005, MNRAS, 361, 405
  • Kennicutt (1998) Kennicutt Jr. R. C., 1998, ApJ, 498, 541
  • Kim et al. (2007) Kim T., Bolton J. S., Viel M., Haehnelt M. G., Carswell R. F., 2007, MNRAS, 382, 1657
  • Kohler et al. (2007) Kohler K., Gnedin N. Y., Hamilton A. J. S., 2007, ApJ, 657, 15
  • Komatsu et al. (2010) Komatsu E., et al., 2010, ArXiv e-prints
  • Larson et al. (2010) Larson D., et al., 2010, ArXiv e-prints
  • Lee et al. (2008) Lee K., Cen R., Gott III J. R., Trac H., 2008, ApJ, 675, 8
  • Lidz et al. (2007) Lidz A., McQuinn M., Zaldarriaga M., Hernquist L., Dutta S., 2007, ApJ, 670, 39
  • Madau et al. (1999) Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648
  • McDonald et al. (2000) McDonald P., et al., 2000, ApJ, 543, 1
  • Miralda-Escudé et al. (2000) Miralda-Escudé J., Haehnelt M., Rees M. J., 2000, ApJ, 530, 1
  • Muñoz et al. (2009) Muñoz J. A., Madau P., Loeb A., Diemand J., 2009, MNRAS, 400, 1593
  • Norman et al. (1998) Norman M. L., Paschos P., Abel T., 1998, Memorie della Societa Astronomica Italiana, 69, 455
  • Okamoto et al. (2008) Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920
  • Pawlik & Schaye (2009) Pawlik A. H., Schaye J., 2009, MNRAS, 396, L46
  • Petkova & Springel (2009) Petkova M., Springel V., 2009, MNRAS, 396, 1383
  • Shin et al. (2008) Shin M.-S., Trac H., Cen R., 2008, ApJ, 681, 756
  • Sokasian et al. (2003) Sokasian A., Abel T., Hernquist L., Springel V., 2003, MNRAS, 344, 607
  • Sokasian et al. (2001) Sokasian A., Abel T., Hernquist L. E., 2001, New Astronomy, 6, 359
  • Sokasian et al. (2004) Sokasian A., Yoshida N., Abel T., Hernquist L., Springel V., 2004, MNRAS, 350, 47
  • Spitzer (1998) Spitzer L., 1998, Physical Processes in the Interstellar Medium
  • Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  • Springel & Hernquist (2003a) Springel V., Hernquist L., 2003a, MNRAS, 339, 289
  • Springel & Hernquist (2003b) Springel V., Hernquist L., 2003b, MNRAS, 339, 312
  • Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astronomy, 6, 79
  • Thomas et al. (2009) Thomas R. M., Zaroubi S., Ciardi B., Pawlik A. H., Labropoulos P., Jelić V., Bernardi G., Brentjens M. A., de Bruyn A. G., Harker G. J. A., Koopmans L. V. E., Mellema G., Pandey V. N., Schaye J., Yatawatta S., 2009, MNRAS, 393, 32
  • Tornatore et al. (2003) Tornatore L., Borgani S., Springel V., Matteucci F., Menci N., Murante G., 2003, MNRAS, 342, 1025
  • Trac et al. (2008) Trac H., Cen R., Loeb A., 2008, ApJL, 689, L81
  • Trac & Gnedin (2009) Trac H., Gnedin N. Y., 2009, ArXiv e-prints: 0906.4348
  • Wadepuhl & Springel (2010) Wadepuhl M., Springel V., 2010, MNRAS, pp 1538–+
  • White et al. (2003) White R. L., Becker R. H., Fan X., Strauss M. A., 2003, AJ, 126, 1
  • Wise & Abel (2008) Wise J. H., Abel T., 2008, ApJ, 684, 1
  • Yoshida et al. (2007) Yoshida N., Oh S. P., Kitayama T., Hernquist L., 2007, ApJ, 663, 687
  • Zahn et al. (2007) Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L., Zaldarriaga M., Furlanetto S. R., 2007, ApJ, 654, 12
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
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

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 description