The impact of local stellar radiation on the HI CDDF

# The impact of local stellar radiation on the HI column density distribution

Alireza Rahmati, Joop Schaye, Andreas H. Pawlik, Milan Raičevic̀
Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, 85748 Garching, Germany
rahmati@strw.leidenuniv.nl
###### Abstract

It is often assumed that local sources of ionizing radiation have little impact on the distribution of neutral hydrogen in the post-reionization Universe. While this is a good assumption for the intergalactic medium, analytic arguments suggest that local sources may typically be more important than the meta-galactic background radiation for high column density absorbers (). We post-process cosmological, hydrodynamical simulations with accurate radiation transport to investigate the impact of local stellar sources on the column density distribution function of neutral hydrogen. We demonstrate that the limited numerical resolution and the simplified treatment of the interstellar medium (ISM) that are typical of the current generation of cosmological simulations provide significant challenges, but that many of the problems can be overcome by taking two steps. First, using ISM particles rather than stellar particles as sources results in a much better sampling of the source distribution, effectively mimicking higher-resolution simulations. Second, by rescaling the source luminosities so that the amount of radiation escaping into the intergalactic medium agrees with that required to produce the observed background radiation, many of the results become insensitive to errors in the predicted fraction of the radiation that escapes the immediate vicinity of the sources. By adopting this strategy and by drastically varying the assumptions about the structure of the unresolved ISM, we conclude that we can robustly estimate the effect of local sources for column densities . However, neither the escape fraction of ionizing radiation nor the effect of local sources on the abundance of systems can be predicted with confidence. We find that local stellar radiation is unimportant for , but that it can affect Lyman Limit and Damped Ly systems. For the impact of local sources increases with redshift. At the abundance of absorbers is substantially reduced for , but at the effect only becomes significant for .

###### keywords:
radiative transfer – methods: numerical – galaxies: evolution – galaxies: formation – galaxies: high-redshift – intergalactic medium

## 1 Introduction

After the reionization of the Universe at , hydrogen residing in the intergalactic medium (IGM) is kept highly ionized primarily by the meta-galactic UV background (UVB). The UVB is the integrated radiation that has been able to escape from sources into the IGM. Because the mean free path of ionizing photons is large compared to the scale at which the sources of ionizing radiation cluster, the UVB is expected to be close to uniform in the IGM. However, close to galaxies, the radiation field is dominated by local sources and hence more inhomogeneous.

Observations of neutral hydrogen (HI) in the Ly forest mostly probe the low-density IGM which is typically far from star-forming regions. The statistical properties of the Ly forest are therefore insensitive to the small-scale fluctuations in the UVB (e.g., Zuo, 1992; Croft, 2004). On the other hand, the neutral hydrogen in Damped Ly (DLA; i.e., ) and Lyman Limit systems (LL; i.e., ), which are thought to originate inside or close to galaxies, might be substantially affected by radiation from local sources that are stronger than the ambient UVB (Gnedin, 2010). As a result, the abundances of the high HI column densities may also change significantly by locally produced radiation.

Indeed, Schaye (2006) and Miralda-Escudé (2005) have used analytic arguments to show that the impact of local radiation may be substantial for LL and (sub-)DLA systems, but should generally be very small at lower column densities. However, relatively little has been done to go beyond idealized analytic arguments and to simulate the effect of local radiation by taking into account the inhomogeneous distribution of sources and gas in and around galaxies. One of the main reasons for this is the computational expense of radiative transfer (RT) calculations in simulations with large numbers of sources. In addition, high resolution is required to capture the distribution of gas on small scales accurately. Because of these difficulties, most simulations of the cosmological HI distribution have ignored the impact of local radiation and focused only on the effect of the UVB (e.g., Katz et al., 1996; Gardner et al., 1997; Haehnelt et al., 1998; Cen et al., 2003; Nagamine et al., 2004; Razoumov et al., 2006; Pontzen et al., 2008; Altay et al., 2011; McQuinn et al., 2011; van de Voort et al., 2012; Rahmati et al., 2012; Bird et al., 2013). A few studies have taken into account local stellar radiation but their results are inconclusive: while Nagamine et al. (2010) and Yajima et al. (2012) found that local stellar radiation has a negligible impact on the distribution of HI, Fumagalli et al. (2011) found that the HI column density distribution above the Lyman limit is reduced by dex due to local stellar radiation.

In this paper, we investigate the impact of local stellar radiation on the HI distribution by combining cosmological hydrodynamical simulations with accurate RT. We find that the inclusion of local sources can dramatically change the predicted abundance of strong DLA systems and, depending on redshift, also those of LL and weak DLA systems. On the other hand, lower HI column densities are hardly affected. We also show that resolution effects have a major impact. For instance, the resolution accessible to current cosmological simulations is insufficient to resolve the interstellar medium (ISM) on the scales relevant for the propagation of ionizing photons. The limited resolution also affects the source distribution which can change the resulting HI distribution, especially in low-mass galaxies. On top of that, assumptions about the structure of the unresolved multiphase ISM significantly affect the escape of stellar radiation into the IGM. Therefore, any attempt to use cosmological simulations to investigate the impact of local stellar radiation on the HI distribution may suffer from serious numerical artifacts.

Some of these difficulties can be circumvented by tuning the luminosities of the sources such that the escaped radiation can account for the observed UVB. Then the interaction between the radiation that reaches the IGM and the intervening gas can be studied on scales that are properly resolved in the simulations. We adopt this procedure to generate the observed UVB for various ISM models but find that our fiducial simulation reproduces the observed UVB without any tuning.

Among the known sources of radiation, quasars and massive stars are the most efficient producers of hydrogen ionizing photons and are therefore thought to be the main contributors to the UVB. Star-forming galaxies however are thought to be the dominant producers of the UVB at (e.g., Haehnelt et al., 2001; Bolton et al., 2005; Faucher-Giguère et al., 2008a). Therefore, we only account for local radiation that is produced by star formation in our simulations. Moreover, we assume that the ionizing emissivity of baryons strictly follows the star formation rate. We use star-forming gas particles rather than young stellar particles as ionizing sources. Since the gas consumption time scale in the ISM is much larger than the lifetime of massive stars ( yr vs. yr), there are many more star-forming gas particles than there are stellar particles young enough to efficiently emit ionizing radiation. Therefore, using star-forming gas particles as ionizing sources allows us to sample the source distribution better and hence to reduce the impact of the limited resolution of cosmological simulations.

The structure of the paper is as follows. in 2 we use the observed relation between the star formation rate and gas surface densities to provide an analytic estimate for the photoionization rate that is produced by young stars in a uniform and optically thick ISM that is in good agreement with our simulation results. In 3 we discuss the details of our numerical simulations and RT calculations. The simulation results are presented in 4. We show how the local stellar radiation can generate the observed UVB. We also discuss the effect of unresolved ISM on the escape fraction of ionizing radiation and we investigate the impact of local stellar radiation on the HI column density distribution. Finally, we conclude in 5.

## 2 Photoionization rate in star-forming regions

A strong correlation between gas surface density and star formation rate has been observed in low-redshift galaxies (e.g., Kennicutt, 1998; Bigiel et al., 2008). In principle, such a relation can be combined with simplified assumptions to derive a typical photoionization rate that is expected from stellar radiation. In this section, we use this approach to estimate the average ionization rate that is produced by star formation as a function of gas (surface) density on galactic scales.

Assuming a Chabrier (2003) initial mass function (IMF), the observed Kennicutt-Schmidt law provides a relation between gas surface density, , and star formation rate surface density, , on kilo-parsec scales (Kennicutt, 1998):

 ˙Σ⋆≈1.5×10−4 M⊙yr−1kpc−2(Σgas1M⊙pc−2)1.4. (1)

Equation (1) can be used to derive a relation between ionizing emissivity (per unit area), , and the gas surface density. As we will discuss in 3.4, stellar population synthesis models indicate that the typical number of ionizing photons produced per unit time by a constant star formation rate is

 ˙Qγ∼2×1053 s−1 (SFR1 M⊙ yr−1). (2)

Furthermore, as both models and observations suggest, the escape fraction of ionizing photons from galaxies is (e.g., Shapley et al., 2006; Gnedin et al., 2008; Vanzella et al., 2010; Yajima et al., 2012; Paardekooper et al., 2011; Kim et al., 2012). This allows us to assume that most of the ionizing photons that are produced by star-forming gas are absorbed on scales kpc. Therefore, the hydrogen photoionization rate, on galactic scales, can be computed:

 Γ⋆=˙ΣγNH=˙Qγ ˙Σ⋆NH, (3)

where is the hydrogen column density which can be obtained from the gas surface density:

 NH≈9.4×1019cm−2 (Σgas1M⊙pc−2) (X0.75), (4)

where is the hydrogen mass fraction. After assuming and substituting equations (1), (2) and (4) in equation (3), one gets

 Γ⋆∼8.5×10−14 s−1(NH1021cm−2)0.4. (5)

If the scale height of the disk is similar to the local Jeans scale, the hydrogen column density can be computed as a function of the hydrogen number density (Schaye, 2001, 2004; Schaye & Dalla Vecchia, 2008)

 NH∼2.8×1021cm−2 (nH1cm−3)1/2T1/24(fgfth)1/2, (6)

where , is the total mass fraction in gas and is the fraction of the pressure that is thermal (i.e., ). In addition, for deriving equation (6) we have adopted an adiabatic index, , mean particle mass, . After eliminating between equation (5) and equation (6), the photoionization rate can be written as a function of the hydrogen number density

 Γ⋆∼1.3×10−13 s−1 (nH1cm−3)0.2(fgfth)0.2T0.24. (7)

Based on equation (7), the photoionization rate produced by star formation is only weakly sensitive to the temperature, total gas fraction, , and the fraction of the pressure that is thermal, . The photoionization rate in equation (7) is also very weakly dependent on the gas density. As we will discuss in 4.1, our simulations show the same trend and are in excellent agreement with this analytic estimate. However, one should note that due to the simplified assumptions we adopted to derive equation (7), it provides only an order-of-magnitude estimate for the typical photoionization rate that is produced by young stars in the ISM, on kilo-parsec scales. This relation does not capture the inhomogeneity of the ISM that may result in large fluctuations in the radiation field on smaller scales. We note that the relation we found here between the gas density and stellar photoionization rate is a direct consequence of the underlying star formation law and is independent of the total amount of star formation in a galaxy.

## 3 Simulation techniques

In this section we describe different parts of our simulations. We start with discussing the details of the hydrodynamical simulations that are post-processed with RT calculations. Then we explain our RT that includes the radiation from local stellar radiation, the UVB and recombination radiation.

### 3.1 Hydrodynamical simulations

We use a set of cosmological simulations that are performed using a modified and extended version of the smoothed particle hydrodynamics (SPH) code GADGET-3 (last described in Springel, 2005). The physical processes that are included in the simulations are identical to what has been used in the reference simulation of the Overwhelmingly Large Simulations (OWLS) described in Schaye et al. (2010). Briefly, we use a subgrid pressure-dependent star formation prescription of Schaye & Dalla Vecchia (2008) which reproduces the observed Kennicutt-Schmidt law. We use the chemodynamics model of Wiersma et al. (2009b) which follows the abundances of eleven elements assuming a Chabrier (2003) IMF. These abundances are used for calculating radiative cooling/heating rates, element-by-element and in the presence of the uniform cosmic microwave background and the Haardt & Madau (2001) UVB model (Wiersma et al., 2009a). We note that the Haardt & Madau (2001) UVB model has been shown to be consistent with observations of HI (Altay et al., 2011; Rahmati et al., 2012) and metal absorption lines (Aguirre et al., 2008). We model galactic winds driven by star formation using a kinetic feedback recipe that assumes of the kinetic energy generated by Type II SNe is injected as outflows with initial velocity of and with a mass loading parameter (Dalla Vecchia & Schaye, 2008).

We adopt cosmological parameters consistent with the most recent WMAP year 7 results: (Komatsu et al., 2011). Our reference simulation has a periodic box of comoving and contains dark matter particles with mass and an equal number of baryons with initial mass . The Plummer equivalent gravitational softening length is set to and is limited to a minimum physical scale of . We also use a simulation with a box size identical to our reference simulation but with 8 (2) times better mass (spatial) resolution to assess the effect of resolution on our findings (see Appendix B.2).

In our hydrodynamical simulations, ISM gas particles (which all have densities ) follow a polytropic equation of state that defines their temperatures. These temperatures are not physical and only measure the imposed pressure (Schaye & Dalla Vecchia, 2008). Therefore, when calculating recombination and collisional ionization rates, we set the temperature of ISM particles to which is the typical temperature of the warm-neutral phase of the ISM. Furthermore, we simplify our RT calculations by assuming that helium and hydrogen absorb the same amount of ionizing radiation per unit mass and by ignoring dust absorption and the possibility that some hydrogen may be molecular (see Rahmati et al., 2012 for more discussion).

For the RT calculations we use TRAPHIC (see Pawlik & Schaye, 2008, 2011; Rahmati et al., 2012). TRAPHIC is an explicitly photon-conserving RT method designed to exploit the full spatial resolution of SPH simulations by transporting radiation directly on the irregular distribution of SPH particles.

The RT calculation in TRAPHIC starts with source particles emitting photon packets to their neighbors. This is done using a set of tessellating emission cones, each subtending a solid angle of . The propagation directions of the photon packets are initially parallel to the central axes of the emission cones. In order to improve the angular sampling of the RT, the orientations of these emission cones are randomly rotated between emission time steps. We adopt for computational efficiency. However, we note that our results are insensitive to the precise value of , thanks to the random rotations.

After emission, photon packets travel along their propagation direction from one SPH particle to its neighbors. Only neighbors that are inside transmission cones can receive photons. Transmission cones are defined as regular cones with opening solid angle , and are centered on the propagation direction. The parameter sets the angular resolution of the RT and we adopt which produces converged results (see Appendix C1 of Rahmati et al., 2012). To guarantee the independence of the RT from the distribution of SPH particles, additional virtual particles (ViP) are introduced. This is done wherever a transmission cone contains no neighboring SPH particle. ViPs do not affect the underlying SPH simulation and are deleted after the photon packets are transferred.

Furthermore, arriving photon packets are merged into a discrete number of reception cones, . This makes the computational cost of RT calculations with TRAPHIC independent of the number of radiation sources. This feature is particularly important for the purpose of the present work as it enables the RT calculations in cosmological density fields with large numbers of sources.

Reception cones are also used for emission from gas particles (e.g., ionizing radiation from star-forming gas particles, recombination radiation). This reduces computational expenses while yielding accurate results. Photon packets are isotropically emitted into reception cones which are already in place at each SPH particle. This obviates the need for constructing any additional emission cone tessellations. Using this recipe, the emission of radiation by star-forming gas particles is identical to that of the emission of diffuse radiation by recombining gas particles, and is described in more detail in Raicevic et al. (in prep.).

Photon packets emitted by the three main radiation sources we consider in our study (i.e., UVB, RR and stellar radiation) are channeled in separate frequency bins and are not merged with each other. This enables us to compute the contribution of each component to the total photoionization rate. The total amount of the absorbed radiation, summed over all frequency bins, determines the ionization state of the absorbing SPH particles. The time-dependent differential equations that control the evolution of the ionization states of different species (e.g., H, He) are solved using a sub-cycling scheme that allows us to choose the RT time step independent of the photoionization and recombination times.

In our RT calculations we use a set of numerical parameters identical to that used in Rahmati et al. (2012). These parameters produce converged results. In addition to the parameters mentioned above, we use 48 neighbors for SPH particles, 5 neighbors for ViPs and an RT time step of . Ionizing photons are propagated at the speed of light, inside the simulation box with absorbing boundaries, until the equilibrium solution for the hydrogen neutral fractions is reached. To facilitate this process, the RT calculation is starting from an initial neutral state, except for the gas at low densities (i.e., ) and high temperatures (i.e., K) which is assumed to be in equilibrium with the UVB photoionization rate and its collisional ionization rate. Typically, the average neutral fraction in the simulation box does not evolve after 2-3 light-crossing times (the light-crossing time for the comoving is Myr at z = 3).

We continue this section by briefly describing the implementation of the UVB, diffuse recombination radiation and stellar radiation.

In principle, local sources of ionization inside the simulation box should be able to generate the UVB. However, the box size of our simulations is smaller than the mean free path of ionizing photons which makes the simulated volume too small to generate the observed UVB intensity (see 4.3). In addition, a considerable fraction of the UVB at is produced by quasars which are not included in our simulations. Therefore, we impose an additional UVB in our simulation box.

The implementation of the UVB is identical to that of Rahmati et al. (2012). The ionizing background radiation is simulated as plain-parallel radiation entering the simulation box from its sides and the injection rate of the UVB ionizing photons is normalized to the desired photoionization rate in the absence of any absorption (i.e., the so-called optically thin limit).

We set the effective photoionization rate and spectral shape of the UVB radiation at different redshifts based on the UVB model of Haardt & Madau (2001) for quasars and galaxies. The same UVB model has been used to calculate heating/cooling in our hydrodynamical simulations and has been shown to be consistent with observations of HI (Altay et al., 2011; Rahmati et al., 2012) and metal absorption lines (Aguirre et al., 2008) at . We adopt the gray approximation instead of an explicit multifrequency treatment for the UVB radiation. Our tests show that a multifrequency treatment of the UVB radiation does not significantly change the resulting hydrogen neutral fractions (see Appendix D1 of Rahmati et al., 2012)

In addition, hydrogen recombination radiation (RR) is simulated by making all SPH particles isotropic radiation sources with emissivities based on their recombination rates. We do not account for the redshifting of the recombination photons and assume that they are monochromatic with energy 13.6 eV (see Raicevic et al. in prep.).

The ionizing photon production rate of star-forming galaxies is dominated by young and massive stars which have relatively short life times. In cosmological simulations with limited mass resolutions, the spatial distribution of newly formed stellar particles (e.g., with ages less than a few tens of Myr) may sample the locally imposed star formation rates relatively poorly. As we show in 4.2, the distribution of star formation is better sampled by star-forming SPH particles and this is particularly important for low and intermediate halo masses. For this reason, we use star-forming particles as sources of stellar ionizing radiation.

For a constant star formation rate, the production rate of ionizing photons reaches an equilibrium value within (i.e., the typical life-time of massive stars). This equilibrium photon production rate per unit star formation rate can be calculated using stellar population synthesis models. We used STARBURST99 (Leitherer et al., 1999) to calculate this emissivity for the Kroupa (2001) and Salpeter (1955) IMF and a metallicity consistent with the typical metallicities of star-forming particles in our simulations (the median metallicity of star-forming particles in our simulations are within the range at redshifts 0 - 5). We found that for a constant star formation rate of and after , the photon production rate of ionizing photons converges to , which is the value we used in equation (2). As shown in Figure 1, the ionizing photon production rate varies only by if the metallicity changes by . Also, reasonable variations in the IMF (e.g., the Chabrier IMF) do not significantly change our adopted value for the photon production rate. For the spectral shape of the stellar radiation, we adopt a blackbody spectrum with temperature which is consistent with the spectrum of massive young stars.

## 4 Results and discussion

In this section, we report our findings based on RT simulations that include the UVB, recombination radiation and local stellar radiation. The RT calculations are performed by post-processing a hydrodynamical simulation with SPH particles in a 6.25 comoving box and at redshifts . In 4.1, we compare different photoionizing components and assess the impact of local stellar radiation on the HI distribution. In 4.2 we illustrate the importance of using star-forming SPH particles instead of using young stellar particles as ionizing sources. In 4.3 we show that in our simulations, the predicted intensity of stellar radiation in the IGM is consistent with the observed intensity of the UVB and calculate the implied average escape fraction. In 4.4 we show that assumptions about the unresolved properties of the ISM are very important and we discuss the impact of local stellar radiation on the HI column density distribution.

### 4.1 The role of local stellar radiation in hydrogen ionization

As we showed in Rahmati et al. (2012), the UVB radiation and collisional ionization are the dominant sources of ionization at low densities where the gas is shock heated and optically thin. However, self-shielding prevents the UVB radiation from penetrating high-density gas. In addition, much of the ISM has temperatures that are too low for collisional ionization to be efficient. Therefore, because UVB photoionization and collisional ionization are both inefficient in those regions, the ionizing radiation from young stars becomes the primary source of ionization.

Figure 2 shows the distribution of neutral hydrogen in and around a galaxy in our reference simulation at redshifts (top row) and 0 (bottom row) and illustrates that the distribution of HI in galaxies may change significantly as a result of local stellar radiation. The total mass of the halo hosting this galaxy at and 0 is and , respectively. Comparing the panels in the left and middle columns, we see that while collisional ionization and photoionization by the UVB ionize the low-density gas around the galaxy, they do not change the high HI column densities in the inner regions. Comparing the yellow regions in the right and middle panels of Figure 2 demonstrates that local stellar radiation significantly changes the distribution of the HI at .

In the left panel of Figure 3 the contributions of the photoionization rates from different radiation components are plotted against the total hydrogen number density at . The purple solid line shows the total photoionization rate. The blue dashed, green dotted and red dot-dashed lines show respectively the contribution of the UVB, recombination radiation (RR) and local stellar radiation (LSR). The resulting hydrogen neutral fraction as a function of gas density is shown by the purple solid curve in the right panel of Figure 3 which illustrates the significant impact of LSR at high densities. We note that there is a sharp feature in the hydrogen neutral fraction at in the presence of local stellar radiation. This feature corresponds to a sharp drop-off in the photoionization rate from local stellar radiation at the same densities (see the red dot-dashed curve in the left panel of Figure 3). The left panel of Figure 4 shows photoionization rate profiles in spherical shells around haloes with at . The right panel of Figure 4 shows the resulting neutral hydrogen fraction profile (purple solid curve) which is compared with the simulation that does not include local stellar radiation (green dashed curve). We note that the trends in the photoionization rate and hydrogen neutral fraction profiles do not strongly depend on the chosen mass bin as long as .

As Figures 3 and 4 show, at low densities (i.e., at large distances from the centers of the halos) the gas is highly ionized by a combination of the UVB and collisional ionization and this optically thin gas does not absorb a significant fraction of the ambient ionizing radiation. Consequently, at , which corresponds to typical distances comoving Mpc from the centers of the haloes, the UVB photoionization rate does not change with density. By increasing the density, or equivalently by decreasing the distance to the center of the haloes, the optical depth increases and eventually the gas becomes self-shielded against the UVB. This causes a sharp drop in the UVB photoionization rate at the self-shielding density (see the blue dashed curves in the left panels of Figures 3 and 4). As shown in Rahmati et al. (2012), the UVB photoionization rate shows a similar behavior in the absence of local stellar radiation. However, local stellar radiation increases the photoionization rate at high and intermediate densities. This decreases the hydrogen neutral fractions around the self-shielding densities, allowing the UVB radiation to penetrate to higher densities. Consequently, the local stellar radiation increases the effective self-shielding threshold against the UVB radiation slightly (by dex) compared to the simulation without the radiation from local stars (not shown).

As shown by the red dot-dashed curve in the left panel of Figure 3, at densities the median photoionization rate due to local stellar radiation increases with decreasing density. The main reason for this is the superposition of radiation from multiple sources (i.e., star-forming SPH particles) as the mean free path of ionizing photons increases with decreasing density (see Figure 11). This is also seen in the left panel of Figure 4 as increasing stellar photoionization rate with increasing the distance from the center of halos for comoving kpc. On the other hand, at densities lower than the star formation threshold (i.e., ), the gas is typically at larger distances from the star-forming regions. Therefore, the photoionization rate of local stellar radiation drops rapidly with decreasing density. The star formation rate density averaged on larger and larger scales becomes increasingly more uniform. This causes the photoionization from galaxies that are emitting ionizing radiation to produce a more uniform photoionization rate at the lowest densities. Note that at the highest densities, the photoionization rate from local stellar radiation agrees well with the analytic estimate presented in 2.

The purple solid curve in the right panel of Figure 3 shows the hydrogen neutral fractions in the presence of the UVB, recombination radiation and local stellar radiation. For comparison, the hydrogen neutral fractions in the absence of local stellar radiation, and for the optically thin gas that is photoionized only by the UVB, are also shown (with the green dashed and brown dotted curves respectively). Hydrogen at densities is self-shielded and mostly neutral if the UVB and recombination radiation are the only sources of photoionization. However, local stellar radiation significantly ionizes the gas at intermediate and high densities. As mentioned in 2, the typical photoionization rate that is produced by star-forming gas is , which is comparable to the photoionization rate of the unattenuated UVB at . Consequently, the hydrogen neutral fractions in the presence of local stellar radiation are much closer to the optically thin case. It is also worth noting that the scatter around the median hydrogen neutral fractions is largest for intermediate densities (i.e., . This is closely related to the large scatter in the photoionization rate produced by local stellar radiation (see the left panel in Figure 3). At these densities, the large scatter in the distances to the nearest sources and RT effects like shadowing produce a large range of photoionization rates. For the gas at very high and very low densities on the other hand, the scatter becomes smaller because the relative distribution of sources with respect to the absorbing gas becomes more uniform.

The trends discussed so far are qualitatively similar at redshifts other than : Although the star formation rates evolve at high densities, the photoionization rate due to local stars is set by the underlying star formation law (see 2) which does not change with time. The peak of the stellar photoionization rate at , produced by the superposition of multiple sources, exists at all redshifts but it moves to slightly higher densities at lower redshifts. At densities immediately below the adopted star formation threshold (i.e., ), the photoionization rate produced by local sources drops rapidly with decreasing density. The distribution of star formation in the simulation box is almost uniform on large scales. Therefore, at the lowest densities, local stellar radiation produces a photoionization rate which is not changing strongly with density. However, at low redshifts (e.g., ), the star formation density decreases significantly and becomes highly non-uniform on the scales probed by our small simulation box. Consequently, the resulting photoionization rate due to young stars does not converge to a constant value at low densities at these redshifts. For the same reason, at low densities the scatter around the median stellar photoionization rate increases with decreasing redshift (not shown). As we will discuss in 4.3, if we correct for the small size of our simulation box, the asymptotic photoionization rate due to stellar radiation that has reached the IGM, is consistent with the intensity of the observed UVB at the same redshift.

### 4.2 Star-forming particles versus stellar particles

In cosmological simulations, the integrated instantaneous star formation rate of the simulation box closely matches the total amount of mass converted into stellar particles. However, due to limited resolution, the spatial distribution of young stellar particles (e.g., those with ages Myr) may not sample the spatial distribution of star formation in individual haloes very well.

This issue is illustrated in the left panel of Figure 5, where the blue crosses show the total instantaneous star formation rates inside galaxies in our reference simulation at . These star formation rates are calculated by the summation of the star formation rates of all SPH particles in a given galaxy. The star formation rates can also be estimated by measuring the average rate by which stellar particles are formed in a given galaxy over some small time interval. The star formation rate calculated using this method (i.e., measuring the rate by which stellar particles are formed during the last 20 Myr), is indicted in the left panel of Figure 5 by the red circles. For massive galaxies, with high star formation rates, the formation rate of stellar particles agrees reasonably well with the star formation rates computed from the gas distribution. However, most haloes with SPH particles (i.e., ) do not contain any young stellar particles despite having non-zero instantaneous star formation rates. Consequently, for the few low-mass galaxies that by chance contain one or more young stellar particles, the implied star formation rates are much higher than the instantaneous rate that corresponds to the gas distribution. Moreover, as the right panel of Figure 5 shows, the number of star-forming particles in a given simulated galaxy is times larger than that of young stellar particles. This ratio could be understood by noting that the observed star formation law implies a gas consumption time scale in the ISM that is times longer than the life times of massive stars. Hence, if star formation is implemented by the stochastic conversion of gas particles into stellar particles, as is the case here, the number of star-forming (i.e., ISM) particles is expected to be times larger than the number of young stellar particles. This ratio will be somewhat smaller in starbursts or if gas particles are allowed to spawn multiple stellar particles, but under realistic conditions, star formation will be sampled substantially better by gas particles than by young stellar particles.

As a result of the above mentioned sampling effects, using stellar particles as ionizing sources (as was done in all previous work, e.g., Nagamine et al., 2010; Fumagalli et al., 2011; Yajima et al., 2012) would underestimate the impact of local stellar radiation on the HI distribution for a large fraction of simulated galaxies. We therefore use star-forming SPH particles instead of young stellar particles as local sources of radiation.

Figure 6 illustrates the difference between the photoionization rate profiles produced by star-forming SPH particles (blue solid curves) and young stellar particles (red dotted curves) at . The left panel of Figure 6 shows this for a halo with (see the top panels of Figure 2), while the right panel illustrates the results for haloes with . For this comparison we imposed the same total number of emitted photons in the simulation box in both cases. This was done by setting the photon production rates of individual gas particles proportional to their star formation rates (see equation 2) and setting the photon production rate of individual young stellar particles proportional to their masses.

For the massive halo, simulating local stellar radiation using SPH particles results in photoionization rate profiles similar to those produced by using stellar particles. This similarity is mainly due to the agreement between the total star formation rate and the rate by which gas is converted into young stellar particles (see the left panel of Figure 5). However, this galaxy has times fewer young stellar particles than star-forming SPH particles (see the right panel of Figure 5) and as the shaded areas in the left panel of Figure 6 show, the scatter in the photoionization rate profile is much larger when stellar particles are used.

For haloes with slightly lower masses, shown in the right panel of Figure 6, the rates by which gas is converted into young stellar particles are similar to the total star formation rates (see the left panel of Figure 5). Despite this agreement, the median photoionization rates for these two cases differ dramatically, being much higher when star-forming gas particles are used as sources. We conclude that using stellar particles as sources results in the underestimation of the ionization impact of local stellar radiation on most of the galaxies in the simulation. The intensity of the UVB produced by the simulation that uses young stellar particles as sources does not agree with the observed UVB intensity, while using star-forming particles resolves this issue, after correcting for the box size (see 4.3).

Although in this section we have shown that using star-forming particles as sources of ionizing radiation helps to resolve the above mentioned sampling issues in post-processing RT simulations, we note that doing the same may not be a good solution for simulations with sufficient resolution to model the cold, interstellar gas phase. Such simulations can capture the effects of the relative motions of stars and gas, e.g., the effect of stars moving out of, or destroying, their parent molecular clouds.

### 4.3 Stellar ionizing radiation, its escape fraction and the buildup of the UVB

Observations show that the luminosity function of bright quasars drops sharply at redshifts (e.g., Hopkins et al., 2007) and models for the cosmic ionizing background indicate that star-forming galaxies dominate the production of hydrogen ionizing photons at (e.g., Haehnelt et al., 2001; Bolton et al., 2005; Faucher-Giguère et al., 2008a). In the simulations, it should therefore be possible to build up the background radiation from the radiation produced by star-forming galaxies at . In this section, we quantify the contribution of the ionizing photons that are produced in stars to the build-up of the UVB radiation in our simulations. In addition, we calculate the implied average escape fraction that is required to generate the observed UVB.

#### 4.3.1 Generating the UVB

Our simulation box is smaller than the typical mean free path of ionizing photons, (see Table 1). Therefore, the IGM gas can receive stellar ionizing radiation from a region that is larger than the simulation box. We note that using periodic boundaries to propagate the stellar ionizing photons is not a good solution to resolve this issue. In addition to increasing the computational expense, using periodic boundaries for RT would require us to account for the cosmological redshifting of ionizing photons. Moreover, because of the small size of our box, stellar photons traveling along paths that are nearly parallel to a side of the box may never intersect an optically thick absorber, if periodic boundaries are being used for RT. However, we can correct for the small size of the simulation box by requiring consistency between the simulated UV intensities that local sources produce in the IGM and existing observations/models of the UVB. In order to do that, we assume an isotropic and homogeneous universe that is in photoionization equilibrium, and , so that we can ignore evolution and redshifting during the travel time of the photons. Then, the volume-weighted mean ionizing flux, , from stars is given by:

 F⋆=∫∞0u⋆ e−rλmfpdr=¯u⋆ λmfp, (8)

where is the photon production rate per unit comoving volume. Moreover, one can express the mean ionizing flux in equation (8) in terms of the volume-weighted mean ionizing flux produced by source that are inside the simulation box,

 Fin⋆=∫αLbox0u⋆ e−rλmfpdr, (9)

where is a geometrical factor. In the absence of any absorption and for uniform and isotropic distributions of gas and sources, is set by the average distance between two random points inside a cube. For a cube with unit length this would yield (Robbins, 1978). Based on the average distance between the low-density gas (e.g., SPH particles with ) and sources (i.e., star-forming particles) in our simulations, we find that the average value of is which varies mildly with redshift from at to at .

As Table 1 shows, our simulation box is much smaller than the mean free path of ionizing photons. Therefore, we can use in equation (9) and get:

 Fin⋆≈¯u⋆ αLbox, (10)

Using equation (8) and (10), the total stellar UVB radiation flux can thus be written as:

 F⋆≈λmfpα Lbox Fin⋆. (11)

Therefore, the volume-weighted photoionization rate due to radiation produced by stars in the simulation, , which is close to the median photoionization rate in low-density gas, can be used to calculate the implied UVB photoionization rate after correcting for the small box size of the simulations:

 ΓUVB,⋆≈λmfpα Lbox Γin⋆. (12)

As mentioned earlier, our simulations show that the photoionization rate from local stellar radiation approaches a density independent rate at low densities. We take this photoionization rate as . In addition, we use a compilation of available Lyman-limit mean free path measurements at different redshifts from the literature (i.e., from Bolton & Haehnelt, 2007; Faucher-Giguère et al., 2008a; Prochaska et al., 2009; Songaila & Cowie, 2010; see Table 1). After converting the Lyman-limit mean free paths into the typical mean free path of ionizing photons with our assumed stellar spectrum222Because the effective hydrogen ionization cross section of stellar ionizing photons that we use in this work, , their typical mean free paths are times longer than the mean free path of Lyman-limit photons which have hydrogen ionizing cross sections ., we use equation (12) to derive the implied UVB photoionization rate.

Figure 7 shows the predicted contribution of stellar radiation to the UVB (filled circles). The error bars reflect the quoted error in the mean free path measurements. For comparison, the observational measurement of the UVB from the Ly effective opacity by Bolton & Haehnelt (2007) and the modeled Haardt & Madau (2001) UVB photoionization rates are also shown by orange diamonds and green dashed curve, respectively. Both the observational measurements and modeled UVB intensities are in good agreement with our simulation results for . However, their UVB intensity is slightly higher than ours at . The reason for this could be the absence of radiation from quasars in our simulations.

#### 4.3.2 Average escape fractions

The average star formation rate density in our simulations is in good agreement with the observed cosmic star formation rate (Schaye et al., 2010). Therefore, the good agreement between the UVB intensities in our simulation and the ones inferred from the observations suggests that at the average escape fractions of stellar ionizing photons in our simulations are also reasonable. However, as we will discuss in Appendix B.2, the structure of the ISM is unresolved in our simulations. Therefore, the fact that the produced escape fractions are reasonable may be coincidental.

Since in our simulations both the intensity of stellar radiation in the IGM and the photon production rate are known, we can measure the mean star formation rate-weighted escape fraction of ionizing photons from galaxies into the IGM (see Appendix C):

 fesc ∼ 10−2 (Γin⋆10−14 s−1)(¯σ⋆2.9×10−18cm2)−1 (13) ×  (˙ρ⋆0.15 M⊙ yr−1cMpc−3)−1 ×  (α0.7)−1(Lbox10 cMpc)−1(1+z4)−2.

The implied escape fractions are shown in the right panel of Figure 7. The simulated escape fractions are between and 5 and they decrease with decreasing redshift below . This result is consistent with previous observational and theoretical studies (e.g., Shapley et al., 2006; Schaye, 2006; Gnedin et al., 2008; Kuhlen & Faucher-Giguère, 2012).

### 4.4 The impact of local stellar radiation on the HI column density distribution

The observed distribution of neutral hydrogen is often quantified by measuring the distribution of HI absorbers with different strengths in the spectra of background quasars (e.g., Kim et al., 2002; Péroux et al., 2005; O’Meara et al., 2007; Noterdaeme et al., 2009; Prochaska et al., 2009; Prochaska & Wolfe, 2009; O’Meara et al., 2012; Noterdaeme et al., 2012). The HI column density distribution function (CDDF) is defined as the number of systems at a given column density, per unit column density, per unit absorption length, :

 f(NHI,z)≡d2ndNHIdX≡d2ndNHIdzH(z)H01(1+z)2. (14)

In order to study the effect of local stellar radiation on the HI CDDF, we project our simulation boxes on a two-dimensional grid and use this to calculate the column densities (see Rahmati et al., 2012 for more details).

#### 4.4.1 Local stellar radiation and the HI CDDF at z=3

The simulated HI CDDF at is shown in the left panel of Figure 8 for different ionizing sources. The blue solid curve shows the HI CDDF in our reference simulation which includes local stellar radiation, the UVB and recombination radiation. For comparison, the HI CDDF without local stellar radiation (i.e., only including the UVB and recombination radiation) is shown with the green dashed curve. As the ratio between these two HI CDDFs in the top section of the left panel in Figure 8 illustrates, the effect of local stellar radiation increases with HI column density and reaches a dex reduction at . For HI column densities lower than on the other hand, the HI CDDF is insensitive to local stellar radiation. These trends are consistent with previous analytic arguments (Miralda-Escudé, 2005; Schaye, 2006) and numerical simulations performed by Fumagalli et al. (2011).

As we discussed in 4.2, using star-forming SPH particles as ionizing sources (i.e., our reference model) results in a better sampling than using young stellar particles as sources. We showed that if one uses young stellar particles as ionizing sources, the impact of local stellar radiation will be under-estimated in a large fraction of low-mass haloes in our simulations (see Figure 6 and 5). To illustrate this difference, the HI CDDF for the simulation in which young stellar particles are used, is shown with the red dotted curve in the left panel of Figure 8. Indeed, the impact of local sources on the HI CDDF is much weaker if we use young stellar particles as ionizing sources. This may partly explain why Yajima et al. (2012) found that local sources did not affect the HI CDDF significantly.

#### 4.4.2 The impact of the unresolved ISM

In our simulations, the ISM is modeled by enforcing a polytropic equation of state on SPH particles with densities . This means that our simulations do not include a cold ( K) interstellar gas phase. This simple ISM modeling introduces uncertainties in the hydrogen neutral fraction calculations of the dense regions in our simulations. In addition, for gas with the mean free path of ionizing photons is unresolved (see Appendix B.1).

To estimate the impact of the structure of the ISM on the HI CDDF and the effect of local sources, we introduce a porous ISM model in which we assume that molecular hydrogen is confined to clouds with such small covering fractions that we can ignore them when performing the RT. Following Rahmati et al. (2012), we use the observationally inferred pressure law of Blitz & Rosolowsky (2006) to compute the molecular fraction (see Appendix A), but now we not only subtract the fraction when projecting the HI distribution to compute the HI CDDF, we also subtract it before doing the RT calculation (note that in the fiducial case we assumed fractions to be zero during the RT calculation). The conversion of diffuse atomic hydrogen into compact molecular clouds which do not absorb ionizing photons should facilitate the propagation of photons from star-forming regions into the IGM. Indeed, we find that replacing our reference uniform ISM model with the porous ISM model increases the resulting UVB photoionization rate by a factor of . Therefore, we decreased the photon production rate of local stellar radiation by a factor of 3, such that the model with a porous ISM yields the same UVB intensity that is generated by our reference simulation and is in agreement with the observed UVB. This factor of 3 reduction could be interpreted as reflecting absorptions in molecular clouds that are hosting young stars333Note that this is not equivalent to . The radiation that leaves star-forming regions is still subject to significant absorption before it reaches the IGM.. One should note that despite the enforced agreement between the generated UVBs, the stellar photoionization rates at high and intermediate densities are different in the two simulations.

The right panel of Figure 8 shows that for the HI CDDF in the simulation with a porous ISM (orchid dot-dashed curve) is almost identical to the HI CDDF in our reference simulation (blue solid curve). This suggests that for these column densities the predicted impact of local stellar radiation on the HI CDDF is robust to uncertainties regarding the small scale structure of the ISM. Instead, its impact is controlled by the amount of radiation that is propagating through the properly resolved intermediate densities towards the IGM, which is constrained by the observed UVB intensity. Note, however, that the impact of local sources on the high HI column density part of the distribution (i.e., ) is in fact highly dependent on the assumptions that are made about the physics of the ISM.

It is also possible that the intense stellar ionizing radiation within the ISM will effectively increase the HI column densities by dissociating hydrogen molecules. However, we implicitly neglected this effect in our simulation with a porous ISM by assuming that molecular clouds are not affected by local stellar radiation. To put an upper limit on the impact of this effect, one can assume that all molecular clouds are completely dissociated by the absorption of stellar radiation. The result of this exercise is shown by the red long-dashed curve in Figure 8. This shows that accounting for dissociation could reduce (or even reverse) the impact of local stellar radiation only at the very high column density end of the HI CDDF (i.e., ).

We stress that none of our ISM models are realistic. However, by considering very different models, we can nevertheless get an idea of the possible impact of our simplified treatment of the ISM and we conclude that our results are relatively robust for .

#### 4.4.3 Evolution

The evolution of the impact of local stellar radiation on the HI CDDF is illustrated in Figure 9 for our fiducial model. Each curve in this figure shows the ratio between the HI CDDF with and without local stellar radiation at a given redshift. To avoid the uncertainties about the conversion of atomic gas into , the HI CDDFs are computed assuming that the neutral gas is fully atomic (i.e., no ). For all redshifts local stellar radiation has only a very small impact on the HI CDDF for but significantly reduces the abundance of systems with . The impact of local stellar radiation increases with redshift for . While local sources significantly reduce the HI CDDF for at , their effect only becomes significant for at . This might be attributed to decrease in the proper sizes of galaxies with redshift.

In Rahmati et al. (2012) we used simulations that include only the UVB and recombination radiation to show that for , the HI CDDF does not evolve at and increases with increasing redshift for . Figure 10 shows that if we also include local stellar radiation, the weak evolution of the HI CDDF in the Lyman Limit range extends to even higher redshifts (i.e., ).

The predicted HI CDDFs are compared to observations in Figure 10. We note that the HI CDDF predicted by our reference simulation is not converged with respect to the size of the simulation box (see Appendix B in Rahmati et al., 2012). Since the photoionization caused by local stellar radiation only exceeds the UVB photoionization rate close to galaxies (see Figure 4), we can assume that the impact of local stellar radiation on the HI CDDF is independent of the size of the simulation box. Therefore, for all the curves shown in Figure 10 we have corrected the box size effect by multiplying the CDDF of the fiducial box by the ratio of the CDDF in a converged simulation (with ) and in our reference simulations with , both in the absence of local stellar radiation (i.e., with the UVB and recombination radiation).

As we showed in Rahmati et al. (2012), the simulated HI CDDFs in the presence of the UVB and recombination radiation is in a reasonable agreement with observations (see also Altay et al., 2011) with a small deviation of dex. But as the top section of Figure 10 illustrates, the addition of local stellar radiation increases this deviation to dex. While the difference between the predicted HI CDDFs and observations is larger at the highest HI column densities (i.e., ), we have seen that in this regime the effect of local sources is sensitive to the complex physics of the ISM, which our simulations do not capture. These very high HI column densities are also sensitive to the strength and the details of different feedback mechanisms (see Altay et al. in prep.). The small but significant discrepancy in the LL and weak DLA regime is therefore more interesting. It is important to confirm this discrepancy with larger simulations, so that a correction for box size will no longer be necessary, but this requires more computing power than is presently available to us.

## 5 Discussion and conclusions

The column density distribution function (CDDF) of neutral hydrogen inferred from observations of quasar absorption lines is the most accurately determined observable of the distribution of gas in the high-redshift Universe. Moreover, the high column density absorbers () arise in gas that is either already part of the ISM or will soon accrete onto a galaxy (van de Voort et al., 2012) and hence these systems directly probe the fuel for star formation.

To predict the distribution of high column density absorbers, it is necessary to combine cosmological hydrodynamical simulations with accurate radiative transfer (RT) of ionizing radiation. Because of the cost and complexity associated with RT calculations that include many sources, nearly all studies have only considered ionization by the ultraviolet background radiation (UVB), whose intensity can be inferred from observations of the Ly forest. However, analytic arguments suggest that the radiation field to which high column density absorbers are exposed is typically dominated by local stellar sources (Miralda-Escudé, 2005; Schaye, 2006). It is therefore important to investigate whether the remarkable success of simulations that consider only the UVB, such as the agreement with the observed CDDF over 10 orders of magnitude in column density at (Altay et al., 2011) as well as its evolution down to (Rahmati et al., 2012), is compromised when local sources are included. Here we addressed this question by repeating some of the simulations of Rahmati et al. (2012) with our RT code TRAPHIC (Pawlik & Schaye, 2008, 2011; Rahmati et al., 2012), but this time including not only the UVB and diffuse recombination radiation, but also local stellar sources.

In agreement with the analytic predictions of Miralda-Escudé (2005) and Schaye (2006), we found that local stellar radiation is unimportant for and dominates over the UVB for high column density absorbers. For all redshifts considered here (i.e., ), local sources strongly reduce the abundance of systems with . The impact of local sources increases with redshift for . At the CDDF is substantially reduced for , but at the effect only becomes significant for . As a result, the remarkable lack of evolution in the CDDF that we found in Rahmati et al. (2012) for , and which is also observed, extends to if local sources are taken into account. On the other hand, the agreement with the observed CDDF is not quite as good as before, with the simulations underpredicting the rates of incidence of absorbers by factors of a few. However, because of the large corrections that we had to make because of the small size of the simulation box used to study the effect of local sources (), this discrepancy will have to be confirmed with larger simulations. Moreover, we did not account for possible hydrodynamical effects that might be caused by extra heating due to ionizing radiation from local sources. This process might change the distribution of gas that is affected by local stellar radiation and requires further investigation.

We found that the average photoionization rate due to young stars in high-density gas is weakly dependent on the gas density and is . We showed analytically that this rate follows directly from the imposed (and observed) Kennicutt-Schmidt star formation law if we assume that most of the ionizing photons that are produced by star-forming gas are absorbed on scales kpc. However, in reality we expect the photoionization rate in the ISM to fluctuate more strongly than predicted by simulations like ours, which lack the resolution required to model the cold, interstellar phase.

Indeed, the spatial resolution that is required for accurate RT of ionizing radiation through the ISM is several orders of magnitude higher than the smallest scales accessible in current cosmological simulations. This makes tackling this problem in the near future hardly feasible and poses a difficult challenge for studying the impact of local stellar radiation on the distribution of HI in and around galaxies. Fortunately, one can circumvent part of the problem by tuning the production rate of ionizing photons (which is equivalent to adjusting the escape fraction of ionizing photons from the unresolved ISM) such that the models reproduce the observed mean photoionization rate in the IGM (after subtracting the contribution from quasars). In other words, if one knows the amount of ionizing radiation that is required to reach the IGM, its ionization impact on the intervening gas can be determined even if we cannot predict what fraction escapes from the immediate vicinity of the young stars.

We adopted this approach but found that tuning was unnecessary for our reference simulation at . Moreover, since our simulations also yield star formation histories that are in good agreement with observations (Schaye et al., 2010), we used them to constrain the implied star formation rate weighted mean escape fraction that relates the predicted star formation rate density to the intensity of the UVB. We found that the average escape fraction in our simulations is at , which agrees with previous constraints on the escape fraction from observations (e.g., Shapley et al., 2006) and theoretical work (e.g., Schaye, 2006; Gnedin et al., 2008; Kuhlen & Faucher-Giguère, 2012).

The limited spatial resolution of cosmological simulations mandates the use of simplified models for the structure of the ISM. To estimate the impact of such subgrid models on the CDDF and on the effect of local sources, we varied some of the underlying assumptions. In particular, we considered a porous ISM model which assumes that molecular hydrogen is confined to clouds with such small covering fractions that we can ignore them when performing the RT. We also considered a model which assumes that all molecular clouds are completely dissociated, but not ionised, by the absorption of stellar radiation. Although none of these models are realistic, we used them to estimate the potential impact of our simplified treatment of the ISM. We found that provided that we rescale the source luminosities so that the different models all reproduce the observed background radiation, the models predict nearly the same HI CDDFs in the regime where the absorbers are well resolved. We therefore concluded that our results on the effect of local sources are relatively robust for , but that their predicted impact is highly sensitive to the assumptions about the ISM for .

Different studies have found qualitatively different results for the impact of local stellar radiation on the CDDF. For instance, Fumagalli et al. (2011) used relatively high-resolution zoomed simulations of individual objects to demonstrate that local stellar radiation significantly reduces the abundance of high column density absorbers. On the other hand, some studies using cosmological simulations similar to those presented here (with roughly the same resolution) found that local stellar radiation has a negligible impact on the CDDF (Nagamine et al., 2010; Yajima et al., 2012).

We found that difference in the resolutions of the simulations that were used in these previous studies may explain their inconsistent findings. Using star particles as sources, as was done by (Nagamine et al., 2010; Yajima et al., 2012), we also found that local sources have a negligible impact on the abundance of strong HI systems. However, we demonstrated that it is possible to dramatically improve the sampling of the distribution of ionizing sources by using star-forming gas particles (i.e., gas with densities at least as high as those typical of the warm ISM), thus effectively increasing the resolution without modifying the time-averaged production rate of ionizing photons. We adopted this strategy in our fiducial models and found, as summarized above, that the radiation from local sources significantly affects the high column density end of the CDDF. This result is in agreement with Fumagalli et al. (2011), Gnedin (2010) and analytic estimates of Miralda-Escudé (2005) and Schaye (2006), and confirms that poor sampling of the distribution of ionizing sources can lead to an under-estimation of the impact of local stellar radiation.

Further progress will require higher resolution simulations and, most importantly, more realistic models for the ISM. In the near future it will remain unfeasible to accomplish this in a cosmologically representative volume. Until this challenge is met, predictions for the escape fractions of ionizing radiation averaged over galaxy populations should be considered highly approximate. Predictions for the abundances of LL and weak DLA systems based on models that neglect local sources of stellar radiation should be interpreted with care, particularly for . Predictions for the CDDF in the strong DLA regime () must be considered highly approximate at all redshifts.

## Acknowledgments

We thank the anonymous referee for a helpful report. We also would like to thank Dušan Kereš, J. Xavier Prochaska and Tom Theuns for valuable discussions. The simulations presented here were run on the Cosmology Machine at the Institute for Computational Cosmology in Durham (which is part of the DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and Durham University) as part of the Virgo Consortium research programme. This work was sponsored with financial support from the Netherlands Organization for Scientific Research (NWO), also through a VIDI grant and an NWO open competition grant. We also benefited from funding from NOVA, from the European Research Council under the European UnionÕs Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement 278594-GasAroundGalaxies and from the Marie Curie Training Network CosmoComp (PITN-GA-2009-238356). AHP receives funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) under grant agreement number 301096-proFeSsOR.

## References

• Aguirre et al. (2008) Aguirre, A., Dow-Hygelund, C., Schaye, J., & Theuns, T. 2008, ApJ , 689, 851
• Altay et al. (2011) Altay, G., Theuns, T., Schaye, J., Crighton, N. H. M., & Dalla Vecchia, C. 2011, ApJL, 737, L37
• Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846
• Bird et al. (2013) Bird, S., Vogelsberger, M., Sijacki, D., et al. 2013, MNRAS , 514
• Blitz & Rosolowsky (2006) Blitz, L., & Rosolowsky, E. 2006, ApJ, 650, 933
• Bolton et al. (2005) Bolton, J. S., Haehnelt, M. G., Viel, M., & Springel, V. 2005, MNRAS , 357, 1178
• Bolton & Haehnelt (2007) Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS , 382, 325
• Cen et al. (2003) Cen, R., Ostriker, J. P., Prochaska, J. X., & Wolfe, A. M. 2003, ApJ , 598, 741
• Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
• Croft (2004) Croft, R. A. C. 2004, ApJ , 610, 642
• Dalla Vecchia & Schaye (2008) Dalla Vecchia, C., & Schaye, J. 2008, MNRAS, 387, 1431
• Faucher-Giguère et al. (2008a) Faucher-Giguère, C.-A., Lidz, A., Hernquist, L., & Zaldarriaga, M. 2008, ApJ , 688, 85
• Fumagalli et al. (2011) Fumagalli, M., Prochaska, J. X., Kasen, D., et al. 2011, MNRAS, 418, 1796
• Gardner et al. (1997) Gardner, J. P., Katz, N., Hernquist, L., & Weinberg, D. H. 1997, ApJ , 484, 31
• Gnedin et al. (2008) Gnedin, N. Y., Kravtsov, A. V., & Chen, H.-W. 2008, ApJ , 672, 765
• Gnedin (2010) Gnedin, N. Y. 2010, ApJL , 721, L79
• Haardt & Madau (2001) Haardt F., Madau P., 2001, in Clusters of Galaxies and the High Redshift Universe Observed in X-rays, Neumann D. M., Tran J. T. V., eds.
• Haehnelt et al. (1998) Haehnelt, M. G., Steinmetz, M., & Rauch, M. 1998, ApJ , 495, 647
• Haehnelt et al. (2001) Haehnelt, M. G., Madau, P., Kudritzki, R., & Haardt, F. 2001, ApJL , 549, L151
• Hopkins et al. (2007) Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ , 654, 731
• Katz et al. (1996) Katz, N., Weinberg, D. H., Hernquist, L., & Miralda-Escude, J. 1996, ApJL , 457, L57
• Kennicutt (1998) Kennicutt, R. C., Jr. 1998, ApJ , 498, 541
• Kim et al. (2012) Kim, J.-h., Krumholz, M. R., Wise, J. H., et al. 2012, arXiv:1210.3361
• Kim et al. (2002) Kim, T.-S., Carswell, R. F., Cristiani, S., D’Odorico, S., & Giallongo, E. 2002, MNRAS , 335, 555
• Komatsu et al. (2011) Komatsu, E., et al. 2011, ApJS, 192, 18
• Kroupa (2001) Kroupa, P. 2001, MNRAS , 322, 231
• Kuhlen & Faucher-Giguère (2012) Kuhlen, M., & Faucher-Giguère, C.-A. 2012, MNRAS , 423, 862
• Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS , 123, 3
• McQuinn et al. (2011) McQuinn, M., Oh, S. P., & Faucher-Giguère, C.-A. 2011, ApJ, 743, 82
• Miralda-Escudé (2005) Miralda-Escudé, J. 2005, ApJL , 620, L91
• Nagamine et al. (2004) Nagamine, K., Springel, V., & Hernquist, L. 2004, MNRAS , 348, 421
• Nagamine et al. (2010) Nagamine, K., Choi, J.-H., & Yajima, H. 2010, ApJL , 725, L219
• Noterdaeme et al. (2009) Noterdaeme, P., Petitjean, P., Ledoux, C., & Srianand, R. 2009, A&A, 505, 1087
• Noterdaeme et al. (2012) Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, arXiv:1210.1213
• O’Meara et al. (2007) O’Meara, J. M., Prochaska, J. X., Burles, S., et al. 2007, ApJ , 656, 666
• O’Meara et al. (2012) O’Meara, J. M., Prochaska, J. X., Worseck, G., Chen, H.-W., & Madau, P. 2012, arXiv:1204.3093
• Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei, 2nd. ed. by D.E. Osterbrock and G.J. Ferland. Sausalito, CA: University Science Books, 2006,
• Paardekooper et al. (2011) Paardekooper, J.-P., Pelupessy, F. I., Altay, G., & Kruip, C. J. H. 2011, A&A, 530, A87
• Pawlik & Schaye (2008) Pawlik, A. H., & Schaye, J. 2008, MNRAS , 389, 651
• Pawlik & Schaye (2011) Pawlik, A. H., & Schaye, J. 2011, MNRAS , 412, 1943
• Péroux et al. (2005) Péroux, C., Dessauges-Zavadsky, M., D’Odorico, S., Sun Kim, T., & McMahon, R. G. 2005, MNRAS , 363, 479
• Pontzen et al. (2008) Pontzen, A., Governato, F., Pettini, M., et al. 2008, MNRAS , 390, 1349
• Prochaska & Wolfe (2009) Prochaska, J. X., & Wolfe, A. M. 2009, ApJ , 696, 1543
• Prochaska et al. (2009) Prochaska, J. X., Worseck, G., & O’Meara, J. M. 2009, ApJL, 705, L113
• Rahmati et al. (2012) Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2012, arXiv:1210.7808, MNRAS in press.
• Razoumov et al. (2006) Razoumov, A. O., Norman, M. L., Prochaska, J. X., & Wolfe, A. M. 2006, ApJ , 645, 55
• Robbins (1978) Robbins, D. 1978, Amer. Math. Monthly 85, 278
• Salpeter (1955) Salpeter, E. E. 1955, ApJ , 121, 161
• Schaye (2001) Schaye, J. 2001, ApJ, 559, 507
• Schaye (2004) Schaye, J. 2004, ApJ , 609, 667
• Schaye (2006) Schaye, J. 2006, ApJ , 643, 59
• Schaye & Dalla Vecchia (2008) Schaye, J., & Dalla Vecchia, C. 2008, MNRAS, 383, 1210
• Schaye et al. (2010) Schaye, J., Dalla Vecchia, C., Booth, C. M., et al. 2010, MNRAS, 402, 1536
• Shapley et al. (2006) Shapley, A. E., Steidel, C. C., Pettini, M., Adelberger, K. L., & Erb, D. K. 2006, ApJ , 651, 688
• Songaila & Cowie (2010) Songaila, A., & Cowie, L. L. 2010, ApJ , 721, 1448
• Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
• Vanzella et al. (2010) Vanzella, E., Giavalisco, M., Inoue, A. K., et al. 2010, ApJ , 725, 1011
• van de Voort et al. (2012) van de Voort, F., Schaye, J., Altay, G., & Theuns, T. 2012, MNRAS , 421, 2809
• Wiersma et al. (2009a) Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009a, MNRAS, 399, 574
• Wiersma et al. (2009b) Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009b, MNRAS, 393, 99
• Yajima et al. (2012) Yajima, H., Choi, J.-H., & Nagamine, K. 2012, MNRAS , 427, 2889
• Zuo (1992) Zuo, L. 1992, MNRAS , 258, 36
• Zwaan et al. (2005) Zwaan, M. A., van der Hulst, J. M., Briggs, F. H., Verheijen, M. A. W., & Ryan-Weber, E. V. 2005, MNRAS , 364, 1467

## Appendix A Hydrogen molecular fraction

To account for the effect of molecular hydrogen, we adopt the same conversion relation that Altay et al. (2011) used to successfully reproduce the HI CDDF high-end cut-off. We follow Blitz & Rosolowsky (2006) and adopt an observationally inferred scaling relation between the gas pressure and the ratio between molecular and total hydrogen surface densities:

 Rmol=(PextP0)α, (15)

where , is the galactic mid-plane pressure, and . Furthermore, if we assume gives also the local mass ratio between molecular and atomic hydrogen, the fraction of gas mass which is in molecular form can be written as

 fH2=MH2MTOT=MH2MH2+MHI=11+R−1mol. (16)

The last equation also assumes that at high densities ionization is not dominant and the gas is either molecular or neutral. In our simulations we model the multiphase ISM by imposing an effective equation of state with pressure for densities , where which is normalized to at the threshold. Therefore we have:

 Pext=P⋆(nHn⋆H)γeff. (17)

To make the Jeans mass and the ratio of the Jeans length to the SPH kernel independent of the density, we use (Schaye & Dalla Vecchia, 2008). Consequently, after combining equations (15), (16) & (17) the fraction of gas mass which is converted into , can be written as:

 fH2=⎛⎝1+A(nHn⋆H)−β⎞⎠−1, (18)

where and .

## Appendix B Resolution effects

### b.1 Limited spatial resolution at high densities

At high densities, ionizing photons are typically absorbed within a short distance from the sources. The propagation of ionizing photons in these regions is therefore controlled by distribution of gas on scales that are comparable to the very short mean free path of ionizing photons. In Figure 11, the mean free path of stellar ionizing photons, 444 is the spectrally averaged hydrogen photoionization cross section for stellar ionizing photons. For the spectral shape of the stellar radiation, we adopt a blackbody spectrum with K., is plotted as a function of density for our reference simulation at . The blue dashed curve in Figure 11 shows the result when only the UVB and recombination radiation are included. The effect of adding local stellar radiation is shown by the red dot-dashed curve. For comparison, also the result for a constant UVB radiation (i.e., optically thin gas) is illustrated with the green dotted curve. The ionization by local stellar radiation decreases the hydrogen neutral fraction at which in turn increases the mean free path of ionizing photons making it comparable to the optically thin case.

The simulations cannot provide any reliable information about the spatial distribution of gas and ionizing sources on scales smaller than the resolution limit (i.e., the typical distance between SPH particles) which is shown by the black dotted line in Figure 11. At densities for which the mean free path of ionizing photons is shorter than the spatial resolution, the RT results may not be accurate since all the photons that are emitted by sources are absorbed by their immediate neighbors. This happens at densities without local stellar radiation (blue dashed curve in Figure 11) and at densities with local stellar radiation (red dot-dashed curve in Figure 11). On the other hand, Rt is irrelevant if the gas is highly neutral, which is predicted to be the case at slightly higher densities (see Figure 3).

### b.2 The impact of a higher resolution on the RT

In Rahmati et al. (2012) we showed that the self-shielding limit is not very sensitive to the resolution of the simulation. On the other hand, the small scale structure of the ISM may significantly change the propagation of stellar radiation and their impact on the HI distribution. While this suggests that the impact of local stellar radiation might be sensitive to the resolution of the underlying simulation, our approach of tuning the escaped stellar radiation such that it can generate the desired UVB, circumvent most resolution effects.

To study the impact of resolution on the propagation of stellar ionizing photons in the ISM and their escape to the IGM, we performed a simulation similar to our reference simulation but using 8 times more dark matter and SPH particles (i.e., using dark matter particles and the same number of SPH particles whose masses are 8 times lower than in our reference simulation). Figure 12 shows that at , the photoionization rate of stellar radiation in the higher resolution simulation (red dashed curve) is qualitatively similar to that obtained from our reference simulation (blue solid curve). As expected, the stellar photoionization rate at the highest densities is in agreement with our analytic estimate in 2. However, because of the shorter inter-particle distances in the higher resolution simulation, the stellar photoionization rate peaks at slightly higher densities. The biggest difference between the stellar photoionization rates of the high-resolution and our reference simulation occurs at the lowest densities. Despite having a times higher total star formation rate, the high-resolution simulation results in a UVB which has a times lower photoionization rate. This means that the effective escape fraction in the high-resolution simulation is times lower than for the reference simulation which has .

Finally, we emphasize that due to our simplified treatment of the ISM, it is not even clear whether the results become more realistic with increasing resolution.

## Appendix C Calculation of the escape fraction

Using the equilibrium photon production rate per unit star formation rate from equation (2), the total production rate of ionizing photons in the simulation box, , can be calculated as a function of the comoving star formation rate density within the simulation box, , :

 ˙Nγ,⋆=2×1053 (˙ρ⋆ L3box1 M⊙ yr−1). (19)

Noting that the typical time the ionizing photons spend inside the simulation box is , and letting be the star formation rate-weighted mean fraction of ionizing photons that escape into the IGM, one can calculate the comoving equilibrium photon number density in the IGM:

 ¯nγ∼fesc α ˙Nγ,⋆L2box c (1+z), (20)

where we ignored absorptions outside the host galaxies because .

The comoving number density of ionizing photons in the simulation box is related to the effective hydrogen photoionization rate they produce in the IGM, :

 ¯nγ=Γin⋆c ¯σ⋆(1+z)3, (21)

where is the effective hydrogen ionization cross section for stellar photons. Combining equations (20) and (21) yields the effective escape fraction of ionizing photons from stars into the IGM:

 fesc∼Γin⋆ L2box¯σ⋆ α ˙Nγ,⋆(1+z)2. (22)

After putting numbers in equation (22) and using equation (19), the escape fraction becomes:

 fesc ∼ 10−2 (Γin⋆10−14 s−1)(¯σ⋆2.9×10−18cm2)−1 (23) ×  (˙ρ⋆0.15 M⊙ yr−1cMpc−3)−1 ×  (α0.7)−1(Lbox10 cMpc)−1(1+z4)−2.
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