Lyman-\alpha from Cold Accretion

# Lyα Cooling Emission from Galaxy Formation

Claude-André Faucher-Giguère11affiliation: Department of Astronomy, Harvard University, Cambridge, MA 02138, USA. 22affiliation: Department of Astronomy and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720-3411, USA. $\star$$\staraffiliation: Miller Fellow; cgiguere@berkeley.edu , Dušan Kereš11affiliation: Department of Astronomy, Harvard University, Cambridge, MA 02138, USA. 22affiliation: Department of Astronomy and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720-3411, USA. \dagger$$\dagger$affiliation: Hubble Fellow , Mark Dijkstra11affiliation: Department of Astronomy, Harvard University, Cambridge, MA 02138, USA. 33affiliation: Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, 85748 Garching, Germany. , Lars Hernquist11affiliation: Department of Astronomy, Harvard University, Cambridge, MA 02138, USA. , Matias Zaldarriaga44affiliation: School of Natural Sciences, Institute for Advanced Study, Princeton, NJ 08540, USA.
###### Abstract

Recent numerical and analytical studies have shown that galaxies accrete most of their baryons via the cold mode, from streams with temperatures K. At these temperatures, the streams should radiate primarily in the Ly line and have therefore been proposed as a model to power the extended, high-redshift objects known as Ly blobs, and may also be relevant for powering a range of less luminous Ly sources. We introduce a new Ly radiative transfer code, , and calculate the transport of the Ly emission from cold accretion in cosmological hydrodynamical simulations. In this paper, we describe our methodology, and address physical and numerical issues that are critical to making accurate predictions for the cooling luminosity, but that have been mostly neglected or treated simplistically so far. In particular, we highlight the importance of self-shielding and of properly treating sub-resolution models in numerical simulations. Most existing simulations do not self-consistently incorporate these effects, which can lead to order-of-magnitude errors in the predicted cooling luminosity. Using a combination of post-processing ionizing radiative transfer and re-simulation techniques, we develop an approximation to the consistent evolution of the self-shielded gas. We quantify the dependence of the Ly cooling luminosity on halo mass at for the simplified problem of pure gas accretion embedded in the cosmic radiation background and without feedback, and present radiative transfer results for a particular system. While pure cooling in massive halos (without additional energy input from star formation and AGN) is in principle sufficient to produce erg s blobs, this requires including energy released in gas of density sufficient to form stars, but which is kept 100% gaseous in our optimistic estimates. Excluding emission from such dense gas yields lower luminosities by up to one to two orders of magnitude at high masses, making it difficult to explain the observed Ly blobs with pure cooling. Resonant scattering produces diffuse Ly halos, even for centrally concentrated emission, and broad double peaked line profiles. In particular, the emergent line widths are in general not representative of the velocity dispersion within galactic halos and cannot be directly used to infer host halo masses.

Galaxies: formation, evolution, high-redshift – cooling flows – radiative transfer
slugcomment: Draft

## 1. Introduction

Steidel et al. (2000) discovered two extremely bright, large, and diffuse Ly-emitting ‘blobs’ in a narrowband survey of the SSA 22 proto-cluster region at (for the early detection of extended Ly emission at , see also Francis et al., 1996, 2001; Keel et al., 1999). These two blobs, labelled LAB1 and LAB2, have physical extent kpc, are more luminous ( erg s) than typical line emitters at the same redshift by a factor , and unlike similar halos around radio galaxies have no detectable radio continuum. Since their discovery, these blobs have become some of the most spectacular of a new class of sources that now counts several tens of members (e.g., Matsuda et al., 2004; Palunas et al., 2004; Francis et al., 2004; Dey et al., 2005; Nilsson et al., 2006; Saito et al., 2006; Smith & Jarvis, 2007; Prescott et al., 2008, 2009; Ouchi et al., 2009; Yang et al., 2009, 2010) and whose nature remains unclear. Because detecting spatially extended Ly emission usually requires deep narrowband imaging, which covers thin redshift slices, only a small volume the Universe has been effectively surveyed for them to date; the blobs and their fainter analogues are therefore likely to be numerous and cosmologically significant. In fact, existing observations and theoretical models indicate that they may be the sites of massive galaxy formation and also may display signatures of associated active galactic nuclei (AGN) and/or supernova feedback. The Ly blobs thus provide a unique opportunity to probe the processes driving and regulating galaxy formation, and may be intimately related to phenomena including proto-clusters, mergers, and submillimeter galaxies.

Although the most extreme Ly blobs have received the most attention, there in fact exists a wide continuum of spatially extended Ly sources at high redshift, for example the ones discovered by Saito et al. (2006) with line luminosities erg s and those discovered by Rauch et al. (2008), with line luminosities as low as erg s. Understanding the nature of these fainter but more numerous sources is equally important to develop a physical picture of galaxy formation. In fact, the fainter sources likely probe different (perhaps earlier) stages of galaxy assembly.

A central puzzle for the Ly blobs is that many of them do not appear to have a central source energetic enough to power their entire Ly emission (e.g., Matsuda et al., 2004; Nilsson et al., 2006). Even for the blobs which do have energetically-sufficient counterparts (e.g., Geach et al., 2005, 2007, 2009; Webb et al., 2009), it is unclear whether that energy can actually couple effectively to the Ly emission. Submillimeter starbursts and obscured AGN imply the presence of large quantities of dust, which acts to destroy Ly photons particularly efficiently (Neufeld, 1990). It is therefore uncertain whether Ly photons produced by such dust-enshrouded sources can escape in significant amounts. Nevertheless, several different mechanisms been proposed to power the blobs and can be broadly divided into three categories:

Embedded star formation or AGN, possibly obscured from direct view by dust, could photoionize the surrounding hydrogen nebula (Moller & Warren, 1998; Haiman & Rees, 2001; Weidinger et al., 2004, 2005; Laursen & Sommer-Larsen, 2007).

Superwinds driven by starburst supernovae could explain the observed sizes and kinematics of the blobs, with the Ly emission being generated in the swept up material (Taniguchi & Shioya, 2000; Taniguchi et al., 2001; Ohyama et al., 2003; Mori et al., 2004; Wilman et al., 2005; Geach et al., 2005). AGNs could also drive similar winds (e.g., Murray et al., 2005; Springel et al., 2005; Di Matteo et al., 2005; Hopkins & Hernquist, 2006).

Cooling radiation, emitted as gas accretes onto forming galaxies, could produce luminous and extended structures. If a large fraction of the accreting gas has a temperature K, then most of the cooling radiation could be Ly (Fabian & Nulsen, 1977; Katz & Gunn, 1991; Hu, 1992; Haiman et al., 2000; Fardal et al., 2001; Birnboim & Dekel, 2003; Kereš et al., 2005; Furlanetto et al., 2005; Dijkstra et al., 2006; Yang et al., 2006; Dijkstra & Loeb, 2009; Goerdt et al., 2010; Dayal et al., 2010).

In this work, motivated by recent progress on our understanding of how galaxies get their gas, we focus on the third possibility. The methods developed could however be applied to the first two classes of models as well, and we plan to extend our calculations to model those processes in the future.

In the classic sketch of galaxy formation (Rees & Ostriker, 1977; Silk, 1977; White & Rees, 1978), gas falling into dark matter halos is shocked and heated to the virial temperature. For a galaxy with a mass similar to that of the Milky Way, the shocked gas attains a temperature K. In the dense inner regions of the halos, this gas efficiently radiates its thermal energy, loses its pressure support, and settles into compact discs where it can form stars. A wealth of recent work however suggests that this picture requires an important modification: most of the gas is never strongly shocked as it flows toward the central forming galaxy, but rather accretes in a “cold mode”, maintaining a temperature K. Moreover, this cold accretion proceeds through dense filaments rather than in a spherically symmetric fashion.

Although the importance of cold accretion in galaxy formation has only recently been demonstrated in high-resolution three-dimensional hydrodynamical simulations (e.g., Fardal et al., 2001; Katz et al., 2003; Kereš et al., 2005, 2009b; Ocvirk et al., 2008; Brooks et al., 2009; Dekel et al., 2009), Binney (1977) had argued on the basis of analytic models of protogalaxy collapse that the amount of shock heating could be small for plausible physical conditions. Moreover, already in the first simulations of forming galaxies, most of the gas never heated above K (Katz & Gunn, 1991) and the importance of filamentary structures was recognized by Katz & White (1993) and Katz et al. (1994). Birnboim & Dekel (2003) carried out a stability analysis, supported by one-dimensional hydrodynamical simulations (see also Birnboim et al., 2007; Dekel & Birnboim, 2006), and found that when the radiative cooling is efficient compared with the infall rate, the post-shock gas becomes unstable and cannot support the shock. When applied to cosmology, their results agree well with those of three-dimensional simulations.

The Ly emission from cold accretion has already been the subject of some studies. Fardal et al. (2001) first evaluated the Ly cooling luminosity from hydrodynamical simulations and suggested that it could account for the Ly blobs discovered by Steidel et al. (2000). Haiman et al. (2000) reached a similar conclusion using simplified analytic arguments. Also using simulations, Furlanetto et al. (2005) found that the Ly cooling radiation from structure formation could account for some, but not all, of the luminosity of the Ly blobs. Yang et al. (2006) studied both the hydrogen and helium cooling radiation using simulations, but found the hydrogen Ly luminosity to be strongly dependent on the self-shielding correction applied. Recently, Dijkstra & Loeb (2009) developed an analytic model and suggested that cooling radiation from the cold mode could account for all the Ly blobs under reasonable assumptions. Using adaptive mesh refinement (AMR) simulations, Goerdt et al. (2010) provided supporting evidence for this picture. In our discussion (§5), we will contrast our main results with those of Goerdt et al. (2010), concluding that the differences with theirs most likely originate from the treatments of self-shielding and sub-resolution modeling, which are a focus of our study.

No study focusing specifically on cooling emission has however combined realistic Ly radiative transfer with hydrodynamical simulations before.111Other authors have included a cooling component in radiative transfer calculations of Ly-emitting galaxies (e.g., Tasitsiomi, 2006a; Laursen & Sommer-Larsen, 2007; Laursen et al., 2009a, a), but have not explicitly separated out the signatures of pure cooling or investigated the important uncertainties in detail.Because Ly photons resonantly scatter, the resultant morphologies, spatial extents, and spectra are strongly modified by radiative transfer effects (e.g., Dijkstra et al., 2006). Since the Ly photons tend to follow paths of least resistance in space and frequency (§4), the spatial geometry and bulk velocity fields play critical roles in determining the radiation transport (for observational evidence of these effects, see e.g. Kunth et al., 1998; Mas-Hesse et al., 2003). Fully three-dimensional calculations are therefore necessary to make realistic predictions. Moreover, both the existing analytical and numerical studies have limitations that could induce important errors in quantities as basic as the integrated Ly luminosity of the cold streams. One such uncertainty arises from the exponential dependence of the Ly emissivity on the gas temperature for the temperatures K that are characteristic of cold accretion (§3). At present, most galaxy formation simulations do not self-consistently predict the temperature distribution within the streams. In fact, existing simulations usually do not follow the transport of the ultra-violet (UV) radiation that ionizes and heats dense gas. This seriously limits the predictive power of these calculations, since small errors in the temperatures can result in large errors in the Ly cooling emission, and potentially grossly violate energy conservation. As we will demonstrate, models of sub-resolution physics in hydrodynamical simulations can also introduce large errors if not properly taken into account. Analytical studies based on energetic considerations are not as sensitive to the temperature of the cold streams (e.g., Dijkstra & Loeb, 2009), but are not immune of uncertainties either, since they rely on assumptions regarding the efficiency of Ly emission. Moreover, their simplified nature does not lend itself to detailed radiative transfer predictions. Resolving these issues is critical to relating the Ly emission from cold accretion to observations.

Our ultimate goal is a systematic investigation of the Ly emission from galaxy formation that is both detailed in its predictions, and robust. By detailed, we envision predictions that can be directly compared with observations, and therefore require both 3D hydrodynamical simulations and realistic radiative transfer. By robust, we mean that the predictions should be free of assumptions that introduce the kind of large uncertainties that existing studies are subject to. Due to the complexity of the problem, this ultimate goal is likely to require a long-term effort. The present paper is dedicated to laying down some of the foundations for this research program. Specifically, we present a new Ly radiative transfer code, named , and describe its application to the cooling radiation in cosmological simulations of galaxy formation (for other applications of Ly radiative transfer codes to hydrodynamical simulations, see e.g. Cantalupo et al., 2005; Tasitsiomi, 2006a; Laursen & Sommer-Larsen, 2007; Laursen et al., 2009a; Kollmeier et al., 2010; Zheng et al., 2010a). We pay particular attention to clarifying the physical and numerical uncertainties of these calculations, in particular with respect to the predicted Ly luminosities, and illustrate the importance of radiative transfer effects. To do so, aside for calculating the Ly luminosities of a sample of halos from a cosmological volume, we focus our radiative transfer calculations on a particular system at and explore variations in both the emission and radiative transfer physics. We limit ourselves to the most basic physical problem of accreting halos embedded in a cosmic ionizing background and neglect feedback processes. Follow up studies will build on the results obtained here and investigate the properties of the Ly emission as a function of halo mass and redshift, and will extend them by incorporating additional physics, including feedback (Faucher-Giguère et al.,, in prep.).

We begin by describing our hydrodynamical simulations in §2. We address the emission of Ly photons in §3 and explicitly demonstrate the sensitivity of the predicted Ly cooling luminosity on assumptions regarding the thermal state of the self-shielded gas. Using a combination of post-processing ionizing radiative transfer and re-simulation techniques, we develop an approximation to consistently model the evolution of the dense gas in the hydrodynamical simulations, resulting in the most robust numerical predictions to date. In §4, we present the results of Ly radiative transfer calculations for a particular system of total mass M at and highlight the role of both the bulk velocity flows and of the resonant scatters in shaping the emergent morphology and spectrum. Finally, we discuss our results and conclude in §5. The Appendices document the radiative transfer code introduced in this work, the ionizing radiative transfer method, and relevant analytical estimates.

Throughout, we assume a cosmology with , as inferred from the Wilkinson Microwave Anisotropy Probe (WMAP) five-year data in combination with baryon acoustic oscillations and supernovae (Komatsu et al., 2009). While some of our hydrodynamical simulations were run with slightly different parameters, none of our conclusions are sensitive to the details of the cosmology. We assume hydrogen and helium mass fractions of and (e.g., Burles et al., 2001), the collisional ionization coefficients given in Katz et al. (1996), the Ly collisional excitation coefficient and average number of Ly photons produced per recombination from Osterbrock & Ferland (2006), and the recombination coefficients in the appendix of Hui & Gnedin (1997). For convenience, some symbols used in this work are defined in Table 1.

## 2. Simulations

### 2.1. Code Details

We compute the hydrodynamics of forming galaxies in a CDM universe using a modified version of the GADGET cosmological simulation code (Springel, 2005). The calculation of the gravitational force uses a combination of the particle mesh algorithm (e.g., Hockney & Eastwood, 1988) for large separations and the hierarchical tree algorithm (e.g., Barnes & Hut, 1986; Hernquist, 1987) at small distances. The gas dynamics is calculated using a smoothed particle hydrodynamics (SPH) algorithm (e.g., Lucy, 1977; Gingold & Monaghan, 1977) that conserves both energy and entropy (Springel & Hernquist, 2002). The modifications with respect to the public version of the code include the treatment of cooling, the effects an uniform ultra-violet background (UVB), and a multiphase star formation algorithm as in Springel & Hernquist (2003). Star formation is implemented by the stochastic spawning of collisionless star particles by the gas particles. In practice, star formation in the multiphase model occurs above a density threshold of cm and is calibrated to the observed Kennicutt (1998) law, although it plays only a tangential role in this work, which focuses the cooling emission. The thermal and ionization properties of the gas are calculated including all relevant processes in a plasma with primordial abundances of hydrogen and helium following Katz et al. (1996).

### 2.2. Simulation Parameters and Halo Identification

We use two types of simulations. To achieve high resolution, we ‘zoom in’ on individual halos within a larger simulation box and only follow the local gas dynamics at the refined resolution. This is done by first running a dark matter only simulation and selecting halos of interests. The simulation is then rerun including gas particles, with 8 times the original mass resolution, in a Lagrangian volume surrounding the halo of interest (e.g., Katz & White, 1993). In this work, we focus on zoom in simulations of an individual halo (labeled A1) selected from a volume of side length 10 comoving Mpc. The A1 halo has a total mass M at . As it is also important to understand the trends and variance between different halos, we simulated an entire cosmological volume consisting of a cubical box with a side length of 40 comoving Mpc. While the resolution in this volume is more limited by computational constraints, it provides us with a large number of halos of different masses. Table 2 lists the simulations used in this work and their parameters. The minimum gas smoothing length is set to 0.1 of the gravitational softening in all the simulations. Given the importance of self-shielding (§3), we rerun our simulations with exactly the same parameters, but with the UVB artificially turned off in regions exceeding a certain density (suffix ssUV), to be discussed in §3.2. All the simulations are also rerun with the star formation model turned off. The simulations without star formation are identified by the additional suffix noSF. The A1 zoom in simulations assume a variant of the Haardt & Madau (1996) model of the UVB, while our cosmological simulations use the more recent model of Faucher-Giguère et al. (2009).

A friends-of-friends (FoF) algorithm (e.g., Davis et al., 1985) with linking length set to in units of the mean interparticle separation is used to identify the dark matter halos in the simulations. The total mass of the particles within each FoF group, , corresponds approximately to the mass in a sphere of mean interior density 180 times the background matter density, (White, 2002). We therefore define the virial radius of a halo as the radius of a sphere containing , , where and is the critical density at . In what follows, we use as a shorthand for . The center of each halo is determined as the point deepest in the gravitational potential. Since we run the FoF algorithm on the dark matter particles only, we multiply the returned masses by to account for the baryons.

As we want to isolate the cold accretion cooling radiation, the simulations studied in this work do no include galactic winds or AGN feedback. It is of course likely that the results would be somewhat modified if these processes were included. We plan to investigate the effects of feedback and their observational manifestations in future work.

### 2.3. Ionization and Thermal Structure

Our basic SPH simulations, like most cosmological simulations to date, assume a UVB that is spatially uniform throughout the simulation volume and calculate the local ionization state of the gas assuming photoionization equilibrium with this background. This approach misses the effects of self-shielding: where the optical depth to ionizing photons from the background is of order unity or more, the gas would in reality be exposed to an attenuated ionizing field. This has consequences for both the ionization state and the thermal properties of the gas. Indirectly, the gas dynamics is also affected by the modification of pressure forces.

The omission of self-shielding implies that (in the absence of local sources) dense gas in the simulations sees an ionizing flux stronger than in reality. As a result, the gas tends to be overionized. This is important for the Ly radiative transfer problem for two reasons. First, the Ly emission mechanisms which seed the Ly photons depend not on the total gas density, but on the number densities of ions (§3.1). Second, the transport of Ly photons depends on the neutral hydrogen distribution, as only this ion provides scattering opacity (§4).

Self-shielding also has important effects on the thermal evolution of the gas. As Fardal et al. (2001) pointed out, neglecting self-shielding in a simulation with a prescribed uniform UVB introduces an artificially high rate of photoheating in dense regions and thus results in overestimated temperatures in these regions. Moreover, the presence of a penetrating ionizing background suppresses the cooling function (e.g., Efstathiou, 1992; Katz et al., 1996; Weinberg et al., 1997; Kereš et al., 2005; Wiersma et al., 2009) and these effects could be amplified by the lack of a dynamical response of the gas to cooling, which would tend to make it denser and hence to cool even more rapidly. Although these effects are critical to accurately predicting the Ly cooling luminosity (§3), previous Ly studies have either neglected them, made simplifying but not necessarily correct assumptions (Fardal et al., 2001; Goerdt et al., 2010), or have explored a range of prescriptions (e.g., Furlanetto et al., 2005; Yang et al., 2006).

We improve significantly over previous work by performing ionizing radiative transfer in post-processing to identify the self-shielded gas, rather than relying on simplified criteria (§3.2). Since knowing the distribution of neutral hydrogen is a fundamental component of Ly radiative transfer, post-processing ionizing radiative transfer has previously been used by other groups for related problems (e.g., Cantalupo et al., 2005; Laursen et al., 2009a; Kollmeier et al., 2010; Zheng et al., 2010a), but never before in focused studies of cooling emission. As we will show, the predicted Ly cooling luminosity is very sensitive to the state of the self-shielded gas. In order to obtain more robust predictions, we rerun our hydrodynamical simulations with the ionizing background turned off in regions above a certain density threshold (informed by our ionizing radiative transfer calculations) as an approximation to the consistent treatment of self-shielding.

For the moment, we pause to discuss the star-forming gas, which also requires special treatment.

### 2.4. The Multiphase ISM

The star-forming gas particles in the multiphase model carry effective ionic densities and temperatures that are mass-weighted averages of the hot and cold components (Springel & Hernquist, 2003). While most of the mass in this model is in the K cold component,222The terminology with respect to temperature is somewhat discrepant in the contexts of galaxy formation and of the interstellar medium (ISM). While the K gas is termed cold in galaxy formation and throughout most of this paper, it is usually qualified as warm in the context of the ISM to distinguish it from the much cooler, star-forming molecular gas. the hot component is in general much hotter ( K). This results in high effective temperatures with simultaneously large neutral fractions, and therefore in high collisional excitation rates (§3.1.2). As we will show, using the effective multiphase temperatures and ionic densities for the star-forming particles yields artificially high cooling luminosities, owing to the non-linearity of the emissivity function with respect to density and temperature. Moreover, in the multiphase model, supernovae are responsible for pressurizing the ISM and are therefore an additional source of energy.

To obtain realistic results uncontaminated by feedback energy, it is necessary to exclude the star-forming particles from the cooling luminosity calculations. We explore two ways of doing this. First, we use the simulations with star formation, but ignore all the multiphase particles in the luminosity calculation. For this case, we assume that the star-forming particles are effectively optically thin for the purpose of transporting the Ly photons. In reality, dust may destroy a portion of those photons, but we do not model this effect for this extreme prescription. This case approximates a multiphase medium in which cold neutral clumps embedded in a hot medium are either so compact that their covering factor is negligible, or in which they simply reflect the Ly photons and contribute only a small effective Ly optical depth (e.g., Neufeld, 1991; Hansen & Oh, 2006). A potential worry with this approach, however, is that it misses the cooling that occurs in particles with density above the star formation threshold. In our second approach, we make sure to capture all the cooling by using identical simulations but with the star formation model turned off, in which case we can simply sum over all the particles. The gas is allowed to become arbitrarily dense (as permitted by the resolution) and the absence of ionizing and mechnical feedback from embedded stars makes the density peaks very optically thick to Ly photons. By comparison with the case of optically thin star-forming gas, this prescription therefore allows us to also quantify the effects of the opacity provided by star-forming particles on the transport problem.

We intentionally do not model the Ly photons produced by stars in this work in order to separate out the properties of pure cooling emission; we briefly discuss their importance in the discussion (§5) and in Appendix A.

## 3. Lyα Emission

### 3.1. Emission Processes

#### 3.1.1 Recombination

Ionizing radiation (either from the cosmic background, local star formation, or an AGN) can photoionize gas that recombines and produces Ly photons. Collisions in gas of sufficient density and temperature can also ionize hydrogen and be followed by the reemission of Ly photons via recombination. We group these two processes in ‘recombination emission.’ Recombination emission produces Ly photons at a rate (in units of ph s cm)

 ϵph,recα=fα,recαBHI(T)nHIIne, (1)

where is the hydrogen case B recombination coefficient and is the average number of Ly photons produced per case B recombination. For gas at K that is optically thick to Lyman series transitions (so that higher-order Lyman series recombination photons can ultimately be degraded into a Ly photon), (Osterbrock & Ferland, 2006). This fraction is only weakly dependent on temperature and so we assume this constant value throughout.

#### 3.1.2 Collisional Excitation

Collisions can also excite the Ly line without ionizing hydrogen. The Ly emissivity from this process is given by

 ϵph,collα=CLyα(T)nHIne, (2)

where is the Ly collisional excitation coefficient in units of ph cm s. Note that the Ly recombination and collisional excitation emissivities scale differently with temperature. Moreover, whereas the recombination term is proportional to the HII number density, the collisional excitation term is proportional to the HI number density. The relative importance of the two processes will therefore depend on the local temperature and ionization state of the gas.

#### 3.1.3 Limiting Equilibrium Cases

We assume ionization equilibrium, which is generally valid since the time scale to reach equilibrium is small compared to the dynamical time scale in both optically thin and self-shielded gas. The statistical equilibrium equation for hydrogen is

 ΓHInHI+ΓHI,c(T)nenHI=αAHI(T)nenHII, (3)

where is the photoionization rate, is the collisional ionization coefficient, and is the case A recombination coefficient.

Physical intuition can be gained by considering the two limiting cases of photoionization equilibrium and of pure collisional ionization equilibrium (CIE; ). As we will show, the two cases are directly relevant to our problem: While most of the cosmic volume is well approximated by the photoionization equilibrium regime, the dense cold gas (including the cold streams of interest) can self-shield from the external ionizing radiation and is then more accurately described by the the pure CIE case. Figure 1 shows from both recombination and collisional excitation. The left panel shows the case of gas in ionization equilibrium with a photoionization rate s (approximately the magnitude of the cosmic ionizing background at , e.g. Faucher-Giguère et al., 2008a, b), and the right panel shows the case of gas in collisional ionization equilibrium. The most important point to note is that the collisional excitation contribution is exponentially sensitive to temperature, at the temperatures K characteristic of cold accretion streams, in both cases. As collisional excitation is the dominant Ly cooling mechanism, it is critical to accurately capture the thermal state of the emitting gas. The curves shown in Figure 1 assume that all the helium is in the form of HeIII for simplicity. The emissivity from both processes is only weakly sensitive to the helium ionized fractions, except for temperatures K for the CIE case and at very high densities in the photoionization case, when most of the hydrogen is neutral and helium dominates the free electrons. Since little Ly emission originates from these regimes, our results are broadly insensitive to the helium ionization state, although we do solve for the correct helium ionization fraction in self-shielded regions in our simulations with the on-the-fly self-shielding approximation (§3.2.2).

In our simulations, the Ly luminosity is evaluated directly from the SPH particles. Specifically, each particle is assigned a Ly luminosity , where the emissivity terms are evaluated using its density, ionization state, and temperature, and is its volume, defined as the ratio of its mass to its density. This approach is desirable as it accurately takes into account the clumping of the gas on small scales, relevant to calculate the emission from the density-squared processes, and because it avoids artificial mixing that could occur if hot and cold phases were averaged in a gridding procedure. Such artificial mixing could boost the predicted Ly luminosity by a large factor owing to the non-linearity of the emissivity function.

### 3.2. Self-Shielding

Since the emission processes scale with density squared (eqs 1-2), the emissivity peaks in the densest regions, which are the most likely to self-shield. Because our hydrodynamical simulations lack proper ionizing radiative transfer, they do not correctly capture self-shielding (§2.3). As explained in §2.4, naively integrating over multiphase SPH particles could also induce large errors in the predicted cooling luminosity. As we will show, it is necessary to both exclude multiphase particles from the calculation and to model self-shielding to accurately predict the cooling luminosity. The rest of this section is dedicated to demonstrating the importance of each potential source of error and to developing a consistent approximation to the Ly cooling luminosity.

We follow the following steps:

1. Naively calculate the Ly luminosity from a simulation with standard UV background and star formation treatments.

2. Using post-processing ionizing radiative transfer, identify the self-shielded gas in step 1.

3. Using the post-processed output, illustrate how the predicted Ly luminosity depends on the assumed thermal state of the self-shielded gas.

4. Rerun a hydrodynamical simulation with the same initial conditions, but with the ionizing background turned off in self-shielded regions on the fly as an approximation to the self-consistent effects of self-shielding.

5. By rerunning identical simulations with star formation turned off, separate the effects of incorrectly including multiphase particles from those of ignoring self-shielding.

#### 3.2.1 Post-Processing Self-Shielding

A technical description of our ionizing radiative transfer code is provided in Appendix B. Briefly, the hydrogen photoionization rate of the cosmic background, , is specified and taken as the boundary condition at the faces of the cubical radiative transfer volume. For the radiative transfer calculations, the simulation outputs are interpolated onto a Cartesian grid taking into account the smoothing kernels, with grid points along each dimension. We employ the fiducial choice and a radiative transfer volume of (1 comoving Mpc/h), centered around each halo considered, which convergence tests suggest is sufficient (§4.3). Rays normal to each of the six faces are then sent inward and the optical depth to ionizing photons is calculated along each ray. Given the attenuated photoionization rate at each point, the ionization equilibrium is updated taking into account photoionization, collisional ionization, and recombination. The procedure is iterated until the ionized fraction has converged in all the cells. In solving for the equilibrium ionization balance, the gas temperatures used are those provided by the hydrodynamical simulation. These should be accurate in the optically thin regions and therefore our scheme should accurately capture the onset of self-shielding. In this post-processing treatment, the temperature structure will however be inaccurate in the self-shielded regions, since the modifications of the heating and cooling functions are not properly modeled in the hydrodynamical calculation. As outlined above, we address this in two ways: first, we explore a range of prescriptions for the self-shielded gas, illustrating the sensitivity of the predictions to these prescriptions; we then subsequently improve the accuracy of our calculations by approximating the self-consistent thermal evolution of the gas with simulations in which the ionizing background is switched off in dense regions. The self-shielded cells are defined as those that see an attenuated ionizing background, with , where is the angle-averaged attenuation factor in the cell after post-processing. Since the optical depth rapidly increases within a self-shielded region, the results are weakly sensitive to the choice of the self-shielding threshold.

To investigate how the Ly luminosity depends on assumptions regarding the self-shielded gas (steps 1, 2, and 3 above), we start with hydrodynamical simulations with standard treatments of the UV background and of star formation. Halos covering a broad range of masses are selected from our cosmological volume simulation gdm40n512 at and the luminosity of each is defined as the sum of the luminosities of the gas particles contained within its virial radius. The same prescriptions are also applied to the A1 halo at , which will be used for the Ly radiative transfer calculations (§4). We explore three cases that are illustrative of the range of possibilities (for similar prescriptions, see Furlanetto et al., 2005):

Naively sum all the particles. In this simplistic prescription, we simply sum all the SPH particles within the virial radius. This includes dense particles that would in reality self-shield but that are artificially illuminated by a uniform ionizing background, and star-forming particles that carry effective multiphase values for their temperature and ionization state, which will introduce luminosity errors as described in §2.4.

No Ly emission from self-shielded gas. Self-shielded gas may become neutral and cool substantially. If the gas is not photoionized and cools below K, its Ly emissivity is severely suppressed (Fig. 1). To illustrate how this case might differ from the naive calculation above, we model it by the extreme assumption that self-shielded gas does not produce Ly photons at all. To transport the Ly photons, we assume that this gas has a temperature K.

Collisional ionization equilibrium. Somewhat intermediate between the two above cases, self-shielded gas could settle to collisional ionization equilibrium (CIE) with a temperature K if gravitational heating is sufficiently efficient. In this case, the neutral hydrogen fraction is given by

 xHI=11+ΓHI,c(TCIE)/αAHI(TCIE) (4)

and the Ly emissivity is a well defined function of temperature, as shown in the right panel of Figure 1. We explore how the luminosity depends on the prescribed CIE temperature, for K and 15,000 K. As we will see in our simulations that approximate the effects of self-shielding, the self-shielded gas cools very effectively, so that higher temperatures are not expected. Furthermore, the proximity of 15,000 K to the peak of the cooling curve implies that this case is already quite optimistic.

For the latter two cases, in which we assume either no emission from self-shielded gas or CIE, we exclude emission from the star-forming particles that might fall outside of the self-shielded regions, in order to avoid potential confusion with artificially high luminosities from multiphase particles. Table 3 summarizes the different prescriptions explored for calculating the Ly cooling luminosity; the above cases are labeled 14.

The left panel of Figure 2 shows how the Ly cooling luminosity varies with halo mass at the fiducial redshift for these different prescriptions. The dashed curve shows an analytic estimate of the average maximum cooling luminosity achievable from the release of gravitational potential energy as a function of halo mass (see Appendix A; also Goerdt et al. 2010 for similar ideas). Briefly, the CDM cosmology predicts the average mass accretion rate onto dark matter halos as a function of mass and redshift (e.g., Neistein & Dekel, 2008; McBride et al., 2009; Fakhouri et al., 2010), as well as the shape of the dark matter halo potential wells (e.g., Hernquist, 1990; Navarro et al., 1997). The product of the halo potential well depth with the gas mass accretion rate provides an estimate of the rate at which gravitational potential energy that can be radiated is “injected” into the halo. Figure 2 assumes an optimistically high efficiency factor (eq. A3). As we discuss in Appendix 2, this gravitational power estimate is not strictly an upper bound for the total amount of cooling achievable even in the case of pure accretion, without feedback from stars and AGN, since the accretion streams are embedded in the cosmic ionizing background, which can transfer further energy to them. This contribution is however included in our simulations and circumstantial evidence suggests that it does not dominate. In fact, in our most realistic simulations with a consistent self-shielding approximation – shown in the right panel of Figure 2 and to be discussed below – the gravitational power upper bound is never systematically violated, and a simulation of our A1 halo with the ionizing background completely turned off yields a nearly equal cooling luminosity (within 30%) at as our prescription 9 when excluding the multiphase star forming regions in the same manner.

The more sophisticated analytic model of Dijkstra & Loeb (2009), shown by the dashed lines for their fiducial efficiency parameter , includes a factor accounting for the decreasing fraction of cold gas in massive halos (e.g., Kereš et al., 2005, 2009b). This model therefore predicts lower cooling luminosities than the above upper bound, with a shallower mass dependence at large masses that is in better agreement with the simulation data points. It is important to note here that the cooling luminosities shown in Figure 2 are “theoretical” or “intrinsic”, meaning that they include all the photons emitted within the virial radii of the halos. These luminosities will in general be higher than the observationally inferred luminosities, which only include the emission above a certain surface brightness threshold determined by the observation. Furthermore, a certain fraction of the emitted photons are in practice absorbed by the intervening intergalactic medium (IGM). The Dijkstra & Loeb (2009) data points plotted here (which were computed for a 50% IGM transmission factor by these authors) have been multiplied by a factor of 2 for a fair comparison with our simulation data points (which assume 100% transmission). In future work, we will quantify how the predicted theoretical luminosities translate into observational ones.

There are two main points to take away from the left panel of Figure 2. First, the predicted Ly cooling luminosity is extremely sensitive to the treatment of the dense gas, with the range exceeding three orders of magnitude at M. Second, some prescriptions actually lead to unphysical results, as comparison with the analytic upper bound indicates that they emit more Ly power than is available from the release of gravitational energy by orders of magnitude (in Appendix A, we show that photoionization from the cosmic background cannot physically account for such large luminosities either). This is the case, in particular, for a standard simulation with optically thin ionizing balance and a multiphase star formation model in which all the gas particles are naively summed over (red x’s). Unsurprisingly, the cases of no Ly emission from self-shielded gas (magenta squares) and of CIE with K (blue circles) yield nearly identical luminosities, since the cooling curve is already strongly suppressed at this temperature. The case of CIE with (cyan diamonds) yields more optimistic Ly luminosities, although these push the limit of the power that can be provided by the release of gravitational potential energy alone (the dashed line in Fig. 2) at low masses. Prescription number 5 (green +’s), which uses a simulation in which the multiphase star formation model was turned off, but with a uniform ionizing background penetrating deep into the dense gas, illustrates that artificial photoheating alone can boost the cooling luminosities by orders of magnitude.

The critical question is therefore: What is the correct Ly cooling luminosity of the cold streams? To address this, we develop a simple approximation to the self-consistent evolution of the gas properties with self-shielding. When self-shielding is properly modeled, it will also be possible to show that naively integrating over multiphase particles alone can also artificially boost the cooling luminosity by a large factor.

#### 3.2.2 On-the-Fly Self-Shielding

In Figure 3 we plot the hydrogen neutral fraction as a function of total proper hydrogen number density for the gas around the A1 halo at . The panel on the left shows this distribution for the standard simulation with a uniform ionizing background. The panel on the right shows exactly the same quantity after the post-processing ionizing radiative transfer. The effect of self-shielding is clear. Roughly, it generates a vertical “plume” above cm, indicating the fact that the gas becomes mostly neutral above this density. This motivates our approximation to the self-consistent evolution of self-shielded gas in the hydrodynamical simulations. Namely, we rerun simulations with exactly the same initial conditions and other physical parameters, but set the ionizing background to zero in regions where the density exceeds the fiducial threshold cm (for an alternative scheme in which the UVB is turned off where the gas is optically thick to ionizing photons on a scale kpc, see Sommer-Larsen, 2006; Laursen & Sommer-Larsen, 2007; Laursen et al., 2009b). By turning the ionizing background off in dense regions on the fly, their thermal and dynamical properties are consistently evolved with the modified cooling and heating functions, and the corresponding dynamical response. The “CDB” simulation analyzed by Goerdt et al. (2010) employed an analogous scheme, but with a density threshold 10 higher, cm; in §5, we argue that this difference likely explains much of the discrepancy with our results. Some uncertainty is introduced by our choice of a fixed density threshold for self-shielding, and in the future it would be useful to improve the methodology by performing proper ionizing radiative transfer on the fly, which our codes do not allow us to do at present. There are however reasons to believe that this choice is a good one, which we outline next.

Figure 4 summarizes the hydrodynamical properties (total gas distribution, neutral gas distribution, and temperature structure) for the A1 system at for the different treatments of self-shielding: standard uniform ionizing background, post-processing ionizing radiative transfer, and the on-the-fly self-shielding approximation. When the ionizing radiative transfer is taken into account, the neutral hydrogen column density of the cold streams can be greatly enhanced, especially in the higher density regions close to the central and satellite galaxies, indicating the fact that they self-shield (some of the cold gas at larger radii however remains optically thin). The ionization structure obtained with the on-the-fly self-shielding approximation is furthermore remarkably similar to the one obtained with the post-processing ray tracing scheme, supporting the validity of using our simple density criterion during the course of the hydrodynamical simulation. The simulation with on-the-fly self-shielding is the most accurate as it consistently captures the thermal evolution and dynamical response of the self-shielded gas. Figure 5 illustrates the effects of self-shielding on the temperature structure of the gas more explicitly: the self-shielded gas with cm is generally cooler (with K) when its evolution is consistently modeled. This simply results from the suppression of artificial photoheating and the enhancement of the cooling function in CIE. At these temperatures, the gas radiates very inefficiently in Ly, which provides further evidence that the simple self-shielding density threshold is not introducing large errors: the Ly cooling luminosity versus halo mass predicted using prescription 9 (discussed below) is quite close to what is obtained by effectively suppressing the Ly emission from all the self-shielded gas, as in prescriptions 2 and 3 in which the self-shielded gas is identified using a ray tracing method and does not assume a particular density threshold. We have also run a simulation of the A1 halo with the ionizing background completely turned off, so that all the cooling in this case originates from gravitational energy and requires no self-shielding correction. The cooling luminosity for this simulation equals the one obtained with the simulation with on-the-fly self-shielding within 30%, when the star-forming regions are identically excised. We are therefore confident that our simple density threshold for self-shielding yields relatively accurate results.

The prescriptions for calculating the Ly cooling luminosity from the simulations with on-the-fly self-shielding approximation are labeled 69 in Table 3 and the corresponding results are shown in the right panel of Figure 2. For the last three (most consistent) cases, the filled symbols show the results for the A1 halo in our zoom in simulation with 8 better mass resolution. Prescription 6 (red x’s), in which self-shielding is modeled but in which we naively sum over multiphase particles, demonstrates how the multiphase particles alone can produce artificially high cooling luminosities; these should therefore always be excluded, or treated separately. As an aside, comparison of prescriptions 1, 5, and 6 indicates that having the high density regions turn into multiphase particles limits the amount of artificial photoheating by effectively shielding the very dense gas. Three prescriptions (7, 8, and 9) correspond to physically plausible cases: one with the on-the-fly self-shielding approximation and star formation, but with multiphase star-forming particles excluded from the Ly luminosity sum (9; blue circles); one also with the on-the-fly self-shielding approximation, but with the star formation model turned off, summed over all the particles (7; magenta squares); and the intermediate case of summing only the particles with cm at which the gas would have become multiphase if the star formation model had been on (cyan diamonds; 8). All three are physically realistic in the sense that they are uncontaminated by either artificial photoheating or by the sub-resolution multiphase model. The only difference between the three cases is in how the gas with cm, the density at which the multiphase model becomes active if on, is treated. Since the multiphase model was calibrated to match the observed Kennicutt-Schmidt relation (Springel & Hernquist, 2003), this threshold density corresponds approximately to the density above which stars should start forming (although the exact value depends on some physical assumptions and may depend on redshift; e.g., Schaye, 2004; Wolfe & Chen, 2006).

In principle, the second approach (prescription 7) might seem more accurate since it captures the entire cooling. However, as is apparent in the radiative transfer results of §4, the extra cooling luminosity is concentrated around the accreting galaxy. It is unclear whether this central cooling emission would be observable in reality for at least two reasons. First, the density of the medium and the immediate proximity of the galaxy imply that locally produced Ly photons could be efficiently destructed by dust (the escape fraction of Ly photons from Lyman break galaxies (LBGs), for example, covers the entire range ; e.g., Kornei et al., 2010). In itself, this is not necessarily an issue for this study in which we focus on a simplified dust-free problem, and would be a well posed problem for follow up studies in which dust would be included. Since the extra cooling luminosity occurs in ISM gas, it should however be accompanied by stellar emission which would most likely swamp it locally (see §5), and in that case is not really cooling luminosity from the cold streams. Second, turning off the multiphase model removes ISM pressurization, which can lead to catastrophic collapse of the gas rich discs and deepen the potential wells, allowing extra energy release. Gravitational interactions with dark matter clumps might also artificially transfer energy to unstable gas discs. Since this prescription assumes that all the baryons are in the gaseous component, while in the central galaxies of actual galaxies in massive halos a large fraction would be locked in stars, it is likely an upper limit. Prescriptions 8 and 9 are more conservative as they exclude all cooling emission occurring at densities cm. We expect these three cases to bracket the true cooling luminosity. A better understanding of the energy release at disc interfaces and within galaxies will likely be required to make more definite predictions and should be addressed in future work.

We conclude by noting that some caution is in order when interpreting the quantitative details of the cooling luminosity predictions for the halos at the low end of the mass range in Figure 2, since they contain relatively few SPH particles ( for M) and may not be well converged. Note, however, that for the realistic prescriptions 79, the slope of the numerically predicted relation agrees well with the analytic expectation based on energy conservation (dashed lines in the Figure), in this regime where the cold mode dominates and where the scaling should apply.

## 4. Lyα Radiative Transfer

Having described our hydrodynamical simulations and the modeling of Ly photon production, we proceed to the Ly radiative transport problem. We use a new three-dimensional Ly radiative transfer code, , described in more detail in Appendix C. To summarize, the fields defining the physical state of the gas from a simulation are interpolated onto a Cartesian grid placed around a halo of interest, as for the post-processing ionizing radiative transfer (§3.2). Monte Carlo Ly photons are then seeded throughout the gridded volume, with a number proportional to the local Ly emissivity (§3). The multiple resonant scatters of each Monte Carlo photon are simulated until escape. As the photons propagate, the Ly image and corresponding spectrum in each pixel on the sky are constructed as seen by an observer on Earth, taking into account cosmological surface brightness dimming. We do not however model the effects of IGM filtering in this work (see §4.2), so that the results are more properly interpreted as the redshifted emission as it emerges from the galactic halos.

In typical astrophysical situations, the optical depth at the center of the Ly resonance is very large, easily or more (eq. C1). As a result, a photon propagates only a short distance before being absorbed by and exciting another neutral hydrogen atom to its state. Because of the high Einstein A coefficient for the transition, another Ly photon is quickly reemitted ( s). Since the reemission is in general in a different direction than the incident photon, the propagation can be effectively pictured as a single photon being scattered and undergoing a random walk. Even if the scattering is coherent in the frame of the scattering atom, the motion of the atom combined with the redirection of the photon in general results in a small shift in the frequency as viewed by an external observer. This results in a random walk in frequency space as well. The frequency-space random walk plays a crucial role in shaping the emergent line profile as a photon tends to escape the medium when it finds itself sufficiently far from the line center that the optical depth it sees is reduced to (unless the medium is so optically thick that the photon spatially diffuses out first). The transport of Ly radiation is thus very different than ordinary lines and a simplified optically thin treatment would lead to fundamental errors both in the theoretical predictions and in interpreting observations.

In this section, we present our basic results on Ly cooling emission radiative transfer. As the focus of this work is to understand the theoretical and numerical uncertainties, and the relevant physical effects, we limit ourselves to examining a particular example, our A1 halo at .

### 4.1. Emission Physics

In §3, we demonstrated how the predicted Ly cooling luminosity depends on the treatment of self-shielding and of sub-resolution physics. We illustrate how the different possible assumptions manifest themselves in other observational properties by performing Ly radiative transfer calculations for some of the prescriptions studied above (Table 3). The results are shown in Figures 6 (prescriptions 1, 3, and 4) and 7 (prescriptions 5, 7, and 8). These fiducial calculations use Monte Carlo photons. For the radiative transfer plots, we subtracted the peculiar motion of the central galaxy with respect to the simulation box (40 km s toward the observer), so that corresponds to the galaxy rest frame. As we had found from our study of a sample of halos in §3.2, incorrectly treating self-shielding (either by assuming an uniform ionizing background, or post-processing ionizing radiative transfer with default simulation temperatures) or including gas particles contaminated by the effective sub-resolution multiphase model can lead to order-of-magnitude overestimates of the cooling luminosity. However, it is not only the total luminosity of a system that can be incorrectly predicted, but also its morphology and line spectrum.

It is easy to understand how different prescriptions for seeding the Ly photons can affect the observed morphology and spectrum. In terms of morphology, different prescriptions seed the photons not only in different amounts, but also in different places. For instance, including artificially luminous multiphase gas particles produces a disproportionally large Ly luminosity concentrated in the densest central parts of the system, approximating a bright point source rather than spatially extended emission from the cold streams. Seeding the photons in different places also has implications for the emergent spectrum: photons seeded deep within dense self-shielded regions have a much larger HI column density to traverse before escape, and so result in more widely separated double peak profiles.

Some specific points are worth noting. As indicated by the surface brightness maps and the circularly averaged surface brightness profiles, the resulting objects are spatially extended by tens of proper kpc, comparable to the virial radius of the halo (see also Dijkstra & Loeb, 2009) and to many of the observed Ly blobs (e.g., Matsuda et al., 2004; Saito et al., 2008). This spatial extent results from a combination of some of the cooling radiation being emitted in the accretion streams far from the galaxy, and of spatial diffusion owing to resonant scattering. However, the spatial extent is strongly dependent on the surface brightness threshold and consequently on the luminosity prescription, so that even an intrinsically diffuse source could appear relatively compact in observations For instance, for our optimistic prescription 7 for the A1 halo at shown in the second row of Figure 7, only the central few kpc would stick out above the surface brightness threshold of erg s cm arcsec of the narrowband images of Matsuda et al. (2004). In this case, the more diffuse cooling halo would be completely missed, but would show up over a larger area in deep long slit spectra sensitive to lower surface brightnesses (e.g., Rauch et al., 2008). The predicted line spectra are double peaked, a common characteristic of Ly radiative transfer reflecting the fact that the photons can escape the medium either on the blue side or on the red side of the line center, where the opacity is typically too extreme (e.g., Neufeld, 1990; Zheng & Miralda-Escudé, 2002a; Dijkstra et al., 2006; Verhamme et al., 2006). In most cases (but not all, reflecting the effects of the complex geometry in the different prescriptions), the blue peak is slightly more pronounced than the red peak. This is a signature of systematic infall in the problem at hand, in which the velocity gradients tend to smear the line opacity on the red side of the Ly line, making it easier for the photons to escape on the blue side. Almost all the star-formation powered Ly emitters (Shapley et al., 2003; Steidel et al., 2010), as well as many Ly blobs and fainter analogues (Matsuda et al., 2006; Saito et al., 2008; Rauch et al., 2008), instead show dominant red peaks indicative of outflows. We do not see this phenomenon here simply because we have not modeled galactic winds in our simulations to simplify the physical problem. Unlike the thin accretion streams (with small covering factor) that produce only slightly stronger blue peaks, outflows (with order unity covering factor) are expected to boost the red peak more drastically (e.g., Dijkstra et al., 2006; Verhamme et al., 2006). Intergalactic absorption, neglected here, would also tend to preferentially suppress the blue peak. In future work, we will include outflows and IGM filtering, which should provide a better match to the observational data.

In §3, we had identified three physically plausible prescriptions for calculating the Ly cooling luminosity (7, 8, and 9); the radiative transfer results for 7 and 9 are shown in the second and third rows of Figure 7. The surface brightness maps now make it clear that the extra luminosity obtained using prescription 7, i.e. when using the on-the-fly self-shielding approximation but turning off the star formation model and summing over the all the particles, is concentrated in the densest central parts of the system. This also manifests itself in the predicted line spectrum, which shows more widely separated peaks as a result of the higher column densities of the gas through which the bulk of the photons must propagate to escape, as discussed in the following section.

### 4.2. Radiative Transfer Physics

Since radiative transfer effects play a crucial role in shaping the observational properties, we pause to illustrate exactly how important each piece of physics is. To do so, we repeat the same fiducial calculations but sequentially turn off different physical effects. We only repeat the calculations for the plausible prescription 9, which suffices to illustrate the physics. The results are presented in Figure 8.

First, we keep resonant scatters but artificially set the gas velocities to zero. In this case, the most important change is in the line spectrum, which becomes a nearly symmetric double peak. This is easily understood in the context of the plane-parallel analytic solution derived by Neufeld (1990) (for an adaptation to spherical geometry, see Dijkstra et al., 2006) for a monochromatic source in an extremely optically thick, static medium. Assuming that the source is located at the center of the slab and that is the line center optical depth from the source to the surface, Neufeld (1990) showed that the emergent spectrum is a symmetric double peak profile, with each peak offset by

 |Δvp|≈191 km s−1(T104 K)1/6(NHI1020 cm−2)1/3 (5)

from the center. In dimensionless units in which is the frequency offset in Doppler-broadening units, . For convenience, the corresponding offset in observed wavelength units can be obtained from

 Δλobs≈8.1 \AA(1+z4)(Δv500 km s−1). (6)

In particular, the emergent line profile and corresponding effective line width are in this case determined by the physics of Ly radiative transfer and of self-shielding (which sets the HI column densities) and have little to do with the global properties of the host halo.

Next, we turn velocities back on but artificially turn off resonant scattering. This is equivalent to assuming that the Ly line were an ordinary, optically thin line. In this case, there are two important changes. First, the morphology on the sky is slightly sharper since there is no longer a spatial diffusion effect arising from photon random walks. The spatial diffusion effect appears subtle here as the cooling emission is intrinsically diffuse. It is better illustrated in Figure 9, in which we have kept the same physical set up but seeded the Ly photons in proportion to the local star formation rate, which produces a much more intrinsically compact source. To facilitate visual comparison of the morphology with the pure cooling cases, we have normalized the star formation Ly emission to the cooling luminosity of the halo. This halo is however forming stars at a rate M/yr, which could result in up to erg s of star formation powered Ly emission if the escape fraction were unity (Leitherer et al., 1999); in this case, it would dominate over the cooling luminosity of the halo. Second, the line spectrum is very different and has essentially lost its double peaked nature. Most importantly, the spectrum (particularly its width) is determined by completely different physics, since it is now directly representative of the velocity dispersion within the host halo through the Doppler effect.

An important implication of the radiative transfer effects on the observed Ly line width concerns the way the masses of extended Ly sources like the Ly blobs are estimated. In fact, it is often assumed that the Ly line width can be associated with random motion within the halo (e.g., Bower et al., 2004; Matsuda et al., 2006). As could be anticipated from more general Ly radiative transfer studies (e.g., Zheng & Miralda-Escudé, 2002a; Dijkstra et al., 2006; Verhamme et al., 2006), our results show that neglecting radiative transfer effects will tend to overestimate the mass, since they alone broaden the Ly line even for completely static emitting media. However, the situation is in reality more complex since the Ly lines that one actually observe on Earth are further filtered by the intervening IGM, which can significantly attenuate the Ly line (e.g., Zheng et al., 2010a, b; Laursen et al., 2010). At , for example, the diffuse IGM transmits only % of photons immediately blueward of Ly and this fraction decreases rapidly with increasing redshift (e.g., Faucher-Giguère et al., 2008d). Moreover, local matter overdensities and systematic infall around massive halos (e.g., Faucher-Giguère et al., 2008c) can amplify and shift the absorption, whereas galactic winds act to “let through” more radiation by redshifting it away from the Ly line center (e.g., Santos, 2004; Dijkstra et al., 2007; Dijkstra & Wyithe, 2010). It is therefore in general difficult to relate the observed Ly line alone to the properties of the halo producing it, and it is not even clear in any given case whether the observed line width over or underestimates the halo velocity dispersion. More systematic studies quantifying the relation between the observed Ly line properties and the host halo mass, extending the type of calculations presented in this work, are certainly warranted. We here simply caution against taking halo masses estimated from the Ly line too literally. This remark applies even if the Ly blobs are not predominantly powered by cooling radiation since radiative transfer effects will necessarily be at play in any Ly source.

### 4.3. Convergence

Most of our radiative transfer calculations employed physical fields stored on a Cartesian grid with cells in a volume of (1 comoving Mpc/h). Since the Ly photon trajectories are not tied to the grid points (see Appendix C), their increments were however finer than that this by a factor of 5. Furthermore, the cell luminosities were always calculated directly from the SPH particles so that they captured the full gas clumping and were not subject to resolution degradation. To verify that our radiative calculations are robust, we have repeated many of them with grid points in the same volume instead. As shown in Figure 10 for prescription 9 applied to the A1 halo at , the calculations in fact appear well converged both in terms of Ly morphology and spectrum.

## 5. Discussion and Conclusion

Motivated by observations of a variety of extended Ly sources at high redshift and by recent studies suggesting that cold mode accretion could account for the majority of them (e.g., Dijkstra & Loeb, 2009), we have used hydrodynamical simulations to predict the Ly cooling emission of forming galaxies. We have introduced a new Ly radiative transfer code, , and for the first time applied such a code to the particular problem of cooling radiation in cosmological hydrodynamical simulations. In this work, which will serve as the foundation for follow up studies of the astrophysical implications, we have quantified the theoretical and numerical uncertainties in predicting the properties of Ly cooling radiation. To do so, we have considered the simplest physical problem of galaxy assembly embedded in the cosmic UV background but without feedback processes. We have shown that the Ly cooling luminosities, morphologies, and spectra of the cold streams are strongly dependent on the treatment of the ionization and thermal state of the emitting gas, and that naive assumptions can produce artificially high cooling luminosities. Previous studies have usually relied on such assumptions, because most existing hydrodynamical simulations do not self-consistently follow the transfer of ionizing radiation (but see Sommer-Larsen, 2006; Laursen & Sommer-Larsen, 2007; Laursen et al., 2009b), whereas the cold streams in reality self-shield. The predicted cooling luminosity is so sensitive to the thermal state of the gas principally because of the exponential dependence of the collisional excitation coefficient on temperature in the range K characteristic of the cold streams. We have also demonstrated that subresolution physics models, when not carefully taken into account, can lead to large errors as well.

Having explicitly illustrated the range of uncertainties, we have made a systematic attempt to converge on the correct answer. Our strategy consisted of first post-processing our hydrodynamical simulations with ionizing radiative transfer to identify the self-shielded regions. By inspection, we found that there is a fairly well defined total hydrogen density cm above which the gas self-shields, at least at the redshifts characteristic of the existing observations. We have then rerun simulations with exactly the same initial conditions and other physical parameters, but with the ionizing background turned off in regions exceeding the fiducial density threshold cm. By turning the ionizing background off in dense regions on the fly, this dense gas is consistently evolved with modified heating and cooling functions, as well as the corresponding dynamical response. These simulations provide us with our best estimates for the actual Ly cooling luminosity of the cold streams, in good agreement with energetic considerations.

For our Ly radiative transfer calculations, we focused on a particular halo (A1) of total mass M at . We have shown that it is not only the integrated Ly luminosity of the system that depends strongly on assumptions regarding the thermal state of the gas, but also its apparent morphology and spectrum. To isolate the role of different physical effects in shaping the observational properties, we have sequentially turned off separate pieces of physics. Both the Ly resonant scatters as well as the bulk velocity structure of the system are critical in shaping the emergent spectrum. The resultant effective line width in general differs strongly from the optically thin expectation owing to these radiative transfer effects, and we therefore caution against measuring the masses of systems based on the Ly line alone. The intrinsically extended nature of the cooling emission and the spatial diffusion owing to resonant scattering combine to produce objects that are spatially extended on the sky, but with an apparent size that depends sensitively on the observational surface brightness threshold. As a result, only the central few kpc of the A1 halo Ly cooling emission at sticks out above the surface brightness threshold erg s cm arcsec typically achieved to date (Fig. 7). If the characteristic size of the Ly emission scales as , then the surface brightness should scale as . For our consistent prescriptions 79, Figure 2 shows that the high-mass luminosity slope varies from about 3/4 to 1 in the most optimistic case, i.e. , so that even the higher-mass systems are unlikely to show up as sources with extents kpc in existing observations from cooling luminosity alone. Fainter sources may however well be detectable in deeper observations, such as the 100-hr long-slit spectrum reported by Rauch et al. (2008).

Our analysis of a sample of halos from a cosmological volume at confirms that (in agreement with previous studies; Haiman et al., 2000; Fardal et al., 2001; Dijkstra & Loeb, 2009; Goerdt et al., 2010) cooling emission alone can in principle produce luminosities erg s sufficient to account for luminous Ly blobs. This however requires quite optimistic assumptions under which most of the energy is released close to, or within, the accreting discs, at densities sufficient to form stars according to the observed Kennicutt-Schmidt relation. In those optimistic predictions based on prescription 7, star formation was artificially turned off (yielding purely gaseous discs), so that it is unlikely that this entire cooling emission can be realized in reality without being overwhelmed by stellar emission.333In their simulations of LBGs at , Laursen et al. 2009a find that only % of the total Ly emission comes from cooling, although we show in Appendix A that this ratio should depend on halo mass and redshift. Our findings in this respect are at odds with the recent simulations of Goerdt et al. (2010), who argue that pure cooling radiation can explain the observed giant Ly blobs. Our investigation of the sources of error in the Ly luminosity predicted from simulations (§3) allows us to understand the differences. First, we have demonstrated that the predicted Ly luminosity is very sensitive to the treatment of self-shielded gas and that at , the characteristic density above which gas self-shields from the UVB is cm. Goerdt et al. (2010) did not explicitly perform ionizing radiative transfer and instead assumed a higher self-shielding threshold of cm. This higher density threshold translates into an overestimate of the Ly luminosity owing to artificial photoheating of gas in the density range cm. Examination of the luminosity-weighted histogram in their Figure 7 in fact indicates that the bulk of the Ly luminosity they predict originates from this density regime (including in clumps that would correspond to satellite galaxies, even when the central galaxy is excised), whereas Ly emission is strongly suppressed at these densities in our approximation to the consistent thermal evolution of self-shielded gas, owing to rapid cooling below K (Fig. 5). While there is some uncertainty in the precise value of the self-shielding threshold, our radiative transfer calculations show that it is clearly below 0.1 cm at the redshifts under consideration (Fig. 3), and the value cm is further supported by recent radiative transfer calculations by other groups (Nagamine et al., 2010; Aubert & Teyssier, 2010). The cleanest contrast between our results and those of Goerdt et al. (2010) is provided by their predictions of the Ly luminosity originating from the “streams” alone, in which they excluded all emission from within 20% of the virial radius. Those can be compared, at best, with our pessimistic prescriptions 8 and 9, in which only gas with cm (above the subresolution star formation threshold) is excluded; our predictions are in those cases a factor lower than theirs. It would not be fair to compare our most optimistic prescription 7 with the Goerdt et al. (2010) results with the central 20% of the virial radius excluded, since most of the extra emission under this prescription occurs very close to, or inside, the central galaxy. Finally, the simulations of Goerdt et al. (2010) included supernova feedback, which provides an additional source of energy and could also contribute significantly to boosting the cooling luminosity, even in the absence of ionizing photons from associated local sources (e.g., Taniguchi & Shioya, 2000; Taniguchi et al., 2001; Ohyama et al., 2003; Mori et al., 2004).

It is further notable that many Ly blobs show spectral signatures of outflows (dominant red peak), at odds with the expectation of a more prominent blue peak in the case of pure cooling of infalling gas.444There are exceptions, for example certain regions of the LAB2 blob at (Wilman et al., 2005). Nonetheless, the release of gravitational potential energy through cooling radiation at observationally interesting levels is an inevitable prediction of galaxy formation in CDM and the growing body of theoretical work supporting the cold mode scenario implies that a significant fraction would come out in Ly. Such sources are therefore poised to be routinely detected in future surveys, and likely account for at least a subpopulation of the more modest extended Ly sources, such as the ones detected by Rauch et al. (2008). In fact, at least three of the Rauch et al. (2008) sources (#15, 36, and 37) show spectral signatures suggestive of infall (with a more prominent blue peak). Interestingly, GALEX non-detections indicate that the bright Ly blobs are much rarer at than at (Keel et al., 2009). This redshift evolution is reminiscent of the gradual disappearance of the cold mode at low redshift predicted by theoretical studies (e.g., Birnboim & Dekel, 2003; Dekel & Birnboim, 2006; Kereš et al., 2005, 2009b). It is important to note that the cold mode could play an important role in the existence and properties of the Ly blobs even if they are not energetically dominated by cooling emission. In fact, the presence of cold neutral gas in galactic halos enhances the conversion stellar or AGN power into Ly photons, and scattering off such gas may be necessary to produce the morphological and spectral properties observed.

Although we have made important progress in accurately predicting the Ly properties of cold accretion, much work remains to be done. For instance, we have focused on the simplest physical problem of accretion embedded in the cosmic UV background and neglected feedback processes. In the actual Universe, galactic winds are observed to be ubiquitous at high redshift (e.g., Shapley et al., 2003; Steidel et al., 2010) and simulations suggest that they are in fact needed to reproduce the observed baryonic mass function of galaxies (e.g., Kereš et al., 2009a; Oppenheimer et al., 2010). Furthermore, AGN can inject large amounts of energy in the surrounding gas (Springel et al., 2005; Hopkins & Hernquist, 2006; DeBuhr et al., 2009), which could radiate even after the nucleus has shut off. Such processes will modify the kinematic and thermal properties of the circumgalactic medium, and therefore its Ly signatures. Moreover, they provide additional sources of energy that could enhance the total emission. The ionizing radiation produced by embedded star formation or AGN (e.g., Kollmeier et al., 2010), the effects of metals on the cooling function (e.g., Cantalupo, 2010), and the destruction of Ly photons by dust (as shown by Laursen et al., 2009b) are also likely to be important to reproduce the observed sources, but have not been modeled here for simplicity. Since cold accretion likely fuels star formation in halo centers, at rates similar to the cold gas accretion rates (e.g., Kereš et al., 2005, 2009b; Faucher-Giguère et al., in prep.), gas accretion and star formation are intrinsically linked and their roles in Ly emission are difficult to decouple. Follow up studies will build upon the technical foundation established in this work and include some of these effects. Ultimately, cosmological statistics like the luminosity and correlation functions of the Ly sources will be predicted. It will also be interesting to consider complementary observational probes of the cold streams, including their Ly polarization (e.g., Rybicki & Loeb, 1999; Dijkstra & Loeb, 2008), their absorption signatures, and other emission lines (including from metals).

On a more basic level, our understanding of the physics of the cold mode is incomplete. In particular, it is not exactly clear how the cold streams release their gravitational energy. At least some simulations indicate that they maintain a roughly constant velocity, rather than accelerate, as they fall into galactic halos (Kereš et al., 2005). This implies that the streams would continuously radiate the gravitational work done on them, perhaps by undergoing a series of weak shocks. However, these weak shocks have not been explicitly identified in existing simulations. Birnboim & Dekel (2003) instead envision a scenario in which the cold streams free fall into the halos and release little energy before hitting the central disc in a strong shock, which may be partially supported by our results which suggest that a large fraction of the energy could be released near the halo center. Infalling cold gas might also release some energy through interactions with halo sub-structure and with the lower density, hot halo gas. These interactions could involve hydrodynamic instabilities below the resolution of current generation cosmological hydrodynamic simulations, both SPH and AMR. While we have taken the pragmatic point of view of predicting the observational signatures implied by our current galaxy formation models, it is conceivable that the existing simulations are subject to numerical limitations. Work is underway to elucidate some of these hydrodynamical issues using the new shock-capturing, moving mesh code AREPO (Springel, 2010). When the capability becomes available, it would also be useful to perform truly self-consistent ionizing radiative transfer on the fly in order to definitively capture the thermal evolution of the gas.

We are grateful to Volker Springel for making the expanded GADGET code available to us, and to Yuval Birnboim, Juna Kollmeier, Adam Lidz, Avi Loeb, Matt McQuinn, Eliot Quataert, Chuck Steidel, and David Weinberg for useful discussions. We also thank the referee, Jesper Sommer-Larsen, for a detailed and constructive review. CAFG is supported by a fellowship from the Miller Institute for Basic Research in Science, and received further support from the Harvard Merit Fellowship and FQRNT during the course of this work. DK is supported by NASA through Hubble Fellowship grant HST-HF-51276.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under contract NAS 5-26555. Partial funding was also provided by NSF grants ACI 96-19019, AST 00-71019, AST 02-06299, AST 03-07690, AST 05-06556, AST 05-06556, AST 09-07969, and PHY 08-55425, and NASA ATP grants NAG5-12140, NAG5-13292, NAG5-13381, and NNG-05GJ40G. Further support was provided by the David and Lucile Packard Foundation, the Alfred P. Sloan Foundation, the John D. and Catherine T. MacArthur Foundation, and Harvard University funds. The computations in this paper were run on the Odyssey cluster supported by the FAS Sciences Division Research Computing Group at Harvard University.

## Appendix A A. Analytic Considerations

In this section, we outline analytic arguments that motivate the expectation of significant Ly emission from galactic gas accretion, at a level comparable to observed Ly blobs. Similar arguments were made by Haiman et al. (2000), Dijkstra & Loeb (2009), and Goerdt et al. (2010); their essence is summarized here as a basis to understand our numerical results. We also provide an analytic estimate of the amount of power that could be contributed by photoionization from the cosmic background.

We consider gas accreting from the IGM to the bottom of a dark matter halo potential well and ask how much energy is available to be radiated. Averaged over the accretion time scale, the gravitational power available can be expressed as , where is the gas mass accretion rate, is the potential difference from the IGM to the radius from the bottom of the potential well at which the gas is assumed to settle, and is an efficiency factor. The efficiency factor accounts the fraction of the gravitational energy that remains in bulk kinetic or thermal form. In addition, only a certain fraction of the energy that is radiated comes out in the Ly line.

Using the fitting formula derived by Neistein & Dekel (2008) for the average halo mass accretion rate, we can write (for the cosmology)

 ˙Mgas≈210 M⊙ yr−1(M1012 M⊙)1.4(1+z4)2.5(fgas0.165), (A1)

where is the fraction of the accreted mass that is gaseous (see also McBride et al., 2009; Fakhouri et al., 2010). For a Navarro et al. (1997) halo profile,

 |ΔΦgrav(rmin)|≈GM200crminln[1+(rmin/r200c)c]ln(1+c)−c/(1+c)≈GM200cr200ccln(1+c)−c/(1+c), (A2)

where is the halo mass defined as that exceeding 200 times the critical density, is the corresponding radius, and is the halo concentration parameter. In the last step, we have assumed to simplify the expression. In general, (; see §2.2), but since may exceed by as much of 75% on cluster scales (White, 2002), can be lower than by a factor of as much as 5 at . For the present crude estimate, we will ignore the distinction between and , which is much smaller at the high redshifts of interest. Combining equations A1 and A2, we find

 ⟨˙Egrav⟩≈3.8×1043 erg s−1feff(M1012 M⊙)1.8(1+z4)3.5(fgas0.165)(c/[ln(1+c)−c/(1+c)]5.2). (A3)

This expression is valid at (where the cosmological constant can be neglected) and fiducially assumes , approximately the median concentration of all massive halos at (e.g., Zhao et al., 2003).

As an estimate of the sensitivity of this analytic prediction to the shape of the dark matter halos, the calculation can be repeated for a Hernquist (1990) profile,

 Φgrav(r)=−GM200cr+a, (A4)

using the relation (e.g., Springel et al., 2005). We then obtain the ratio

 ⟨˙Egrav⟩Hernquist⟨˙Egrav⟩NFW=√ln(1+c)−c/(1+c)2. (A5)

For , this ratio takes the values , respectively. Since these differences, at the tens of percent level, are subdominant compared to the other sources of uncertainty, we will only plot the prediction for the NFW model in this work.

In the absence of additional sources heating (see below, where we consider the importance of the ionizing background), the actual Ly cooling luminosity emergent from a given halo will in general be somewhat lower than predicted by equation A3, since some of the gravitational power will be converted to kinetic and thermal forms, rather than radiated, and the radiated fraction will not come out entirely in Ly. For halos dominated by the cold mode, however, the is expected to be significant, perhaps up to % if half of the gravitational energy is radiated and most of this radiation comes out in Ly. At higher masses, the efficiency will be suppressed as the hot mode becomes more important and a higher fraction of the baryons are accreted in the form of stars.

The analytic model of Dijkstra & Loeb (2009) is more sophisticated than the simple energetic argument above, as it accounts for the fact that the fraction of cold gas in halos diminishes with increasing mass. Since only the cold gas radiates efficiently in Ly, this yields a shallower slope for the relation at high masses, which is in fact in better agreement with the simulation results (c.f. Fig. 2). Equation A3 is thus really an upper bound that we do not expect to be saturated in this regime.

In addition to gravitational potential energy, there is an additional source of power that can be converted into Ly photons even in the case of pure accretion, without feedback from stars or AGN: the cosmic ionizing background. An upper bound for how much power the ionizing background can contribute to Ly cooling within a dark matter halo can be estimated as the total inward flux of energy from ionizing photons across a sphere of a virial radius. A more realistic estimate is obtained by multiplying this quantity by a factor accounting for the fact that only a small fraction of these ionizing photons are actually absorbed within the virial shell:

 ˙Eion=4πr2virFionfcov, (A6)

where and is the specific intensity of the ionizing background, assumed homogeneous and isotropic, which is a good assumption just above the hydrogen ionization edge at (e.g., Faucher-Giguère et al., 2009). Since the background spectrum has a strong absorption edge at 4 owing to intergalactic HeII absorption, we can optimistically take for , and beyond . Then:

 ˙Eion≈12π2r2virνHIJνHIfcov, (A7)

For the hydrogen photoionization rate s measured from the Ly forest at (Faucher-Giguère et al., 2008b), this spectrum implies erg s cm Hz. Expressing in terms of halo mass, this gives

 ˙Eion≈1.9×1041 erg s−1(Mh1012 M⊙)2/3(41+z)2(JνHI1.5×10−22 erg s−1 cm−2 Hz−1)(fcov0.05). (A8)

Comparison of equations A3 and A8 suggests that photoionization from the cosmic background can account for only about a percent of the maximum power available from the release of gravitational potential energy at the fiducial halo mass M. Owing to the different mass dependences, and because our numerical results indicate that the gravitational power upper bound can be far from saturated in some plausible cases (§3.2), we however cannot exclude that the ionizing background contributes a significant fraction of the Ly cooling emission in some regimes. This contribution is included in our simulations, though, and while it is difficult to disentangle from gravitational energy, circumstantial evidence suggests that it does not dominate. Indeed, in our most realistic simulations with a consistent self-shielding approximation the gravitational power upper bound is never systematically violated (see the right panel of Fig. 2), and a simulation of our A1 halo with the ionizing background completely turned off yields a nearly equal cooling luminosity (within 30%) at as our prescription 9 when excluding the multiphase star forming regions in the same manner.

Although we intentionally ignore Ly photons produced by star formation throughout most of this work, it is interesting to analytically estimate the photoionization stellar Ly emission relative to pure cooling under simplified assumptions. Numerical simulations indicate that, in the absence of strong feedback, the star formation rate in high-redshift halos closely follows the cold gas accretion rate (Kereš et al., 2005, 2009b; Faucher-Giguère et al., in prep.), i.e. , except for modest . Stellar Ly emission results from the conversion of the ionizing radiation from young stars via absorption and recombination by the surrounding medium, and therefore scales with the star formation rate, with a proportionality constant that depends on the initial mass function. For standard assumptions, erg s [SFR/(M yr)] (e.g., Leitherer et al., 1999). If we assume that , 555This assumption is increasingly violated at high halo masses, where hot halo gas becomes more prevalent, but since the cooling Ly emission also scales with instead of to first order, the errors made should approximately cancel when taking the ratio of stellar to cooling Ly below. then we can again make use of the average fitting formula in equation A1 to derive the average stellar Ly luminosity as a function of halo mass and redshift:

 ⟨LSFα⟩≈2.1×1044 erg s−1fα,esc(M1012 M⊙)1.4(1+z4)2.5(fgas0.165), (A9)

where the factor quantifies the fraction of Ly photons that avoid destruction by dust. We can then combine this with equation A3 to find the characteristic ratio of stellar to cooling Ly:

 ⟨LSFα⟩⟨˙Egrav⟩≈5.53(fα,escfeff)(M1012 M⊙)−0.4(1+z4)−15.2c/[ln(1+c)−c/(1+c)]. (A10)

The main point to take away is that it is a priori difficult to predict whether stellar Ly emission should dominate over cooling Ly, since the ratio depends on a number of uncertain parameters, including and . Furthermore, the ratio depends on both halo mass and redshift, so that theoretical predictions based on model galaxies at fixed and (e.g., Laursen et al., 2009a, b) cannot be universally applied. The actual ratio of stellar ionization to cooling Ly is likely to be somewhat suppressed relative to the simple scaling in equation A10 because galactic winds can remove a large fraction of the gas from the star-forming reservoir (e.g., Murray et al., 2005; Oppenheimer et al., 2010), therefore suppressing the Ly emission from star formation ionization, without significantly affecting the accretion of cold material and hence the cooling Ly.

## Appendix B B. Ionizing Radiative Transfer

The opacity in the Ly line (eq. C1) and the Ly emissivity (eq. 1 and 2) depend on the ionization state of the hydrogen atoms. The SPH code GADGET calculates the ionization of gas elements assuming photoionization equilibrium with a uniform UV background (e.g., Haardt & Madau, 1996; Faucher-Giguère et al., 2009) (see §2.3). This optically thin treatment misses the effects of self-shielding in dense regions (§2.3). To correct for this, we use a post-processing method to solve the radiative transfer problem on the gas distribution provided by GADGET to obtain more realistic ionization fractions and to identify the regions that self-shield from external radiation. The basic assumption is that the gas dynamics does not significantly depend on the radiative transfer effects missed in the SPH calculation. However, re-simulations are also used to approximate the dynamical and thermodynamical response to self-shielding (§3.2.2).

We use a ray tracing algorithm on the same grid as the Ly radiative transfer (Appendix C). Specifically, we send a normal ray from each cell on the surface of the radiative transfer volume along the inward direction (, , , , , or ) and iteratively solve for the equilibrium ionization structure. For this calculation, we take the gas temperatures provided by the SPH calculation and assume that all the helium is in the form of HeIII. These are good approximations for the gas outside of self-shielded regions and they therefore correctly capture the onset of self-shielding. In principle, the distribution of helium ionization affects the predicted Ly luminosity through the abundance of free electrons. The exact ionization state of helium can however matter significantly only inside self-shielded gas, since the free electron number density is dominated by hydrogen elsewhere. Assuming that all the helium is in HeIII overestimates the number of free electrons and hence the recombination and collisional excitation rates. However, the overestimation is only significant for temperatures K, at which the self-shielded gas is found to contribute negligibly to the total Ly emission (§3.2). During the re-simulations, turning off the ionizing background in dense gas results in most of the helium being in the form of HeI in self-shielded gas, as should be the case since HeI traces HI for realistic ionizing spectra. Errors associated with the helium ionization state are therefore expected to be small in all cases. The rays are assumed to originate from a diffuse ionizing background with hydrogen photoionization rate . For simplicity, we only update the hydrogen ionization fractions and only keep track of the photoionization rate along each ray, rather than the full ionizing spectrum. This is a fair approximation because of the small thickness of the self-shielding layer, in which the gas transitions from almost fully ionized to almost completely neutral. Indeed, the thickness of this layer is approximately equal to the mean free path of ionizing photons within it:

 Δlss∼Δlmfp∼1nHIσHI(νHI)≈5 pc(nHI0.01 cm−3)−1, (B1)

where is the hydrogen photoionization cross section and is its ionization frequency, corresponding to 13.6 eV. This mean free path is generally much smaller than the spatial resolution of both our SPH simulations (Table 2) and our radiative transfer grids. One circumstance in which the self-shielding layer could be significantly thicker would be the case in which the illumination is dominated by a nearby quasar (e.g., Cantalupo et al., 2008; Kollmeier et al., 2010); we do not attempt to model this here, although it would be interesting to explore in the future.

Along each ray , the optical depth from the background is calculated assuming that all the ionizing photons effectively have a frequency , , and the photoionization rate is correspondingly attenuated as . The equilibrium solution is obtained iteratively as follows. Trial values for the ionized fraction of hydrogen in each cell are first arbitrarily chosen. The optical depth from the background to each cell, along each normal direction, is then computed. The effective photoionization rate within the cell is taken to be the average of the attenuated value from each direction, . Using this photoionization rate, an updated hydrogen fraction is calculated using the balance equation 3. For this step, the electron number density is taken equal to its value after the previous update. The procedure is repeated until the ionization fraction has converged in each cell. The convergence criterion is set to one part in .

In reality, the gas temperatures provided by the basic hydrodynamical simulations are in error within self-shielded regions. We thus store the mean attenuation from the background, , in each cell as an indicator of self-shielding. Different prescriptions for the self-shielded gas and re-simulation procedures are explored in §3.2.

## Appendix C C. Lyα Radiative Transfer: αRt

In this Appendix, we describe the new three-dimensional Ly radiative transfer code, called “”, developed for this project.

### c.1. Basic Physics and Algorithm

The basic algorithm is the standard Monte Carlo method, which follows the scatters of a prescribed number, , of Ly photons and statistically estimates the emergent morphology and line profile from them (e.g., Zheng & Miralda-Escudé, 2002a; Cantalupo et al., 2005; Dijkstra et al., 2006; Verhamme et al., 2006; Tasitsiomi, 2006a; Laursen & Sommer-Larsen, 2007; Laursen et al., 2009a). The principal improvements over many of these codes are that it is fully three-dimensional and its implementation is parallel to take advantage of distributed computing resources. The code therefore scales well to large problems with complex geometries and so can be run on the outputs of large-scale hydrodynamical simulations. It produces Ly images, spectra, as well as integral field images (spectra as a function of position on the sky). As it was independently developed, its results provide useful checks of existing numerical calculations. It has been demonstrated (Tasitsiomi, 2006a; Laursen et al., 2009a) that the Monte Carlo algorithm can be combined with mesh refinement techniques relatively straightforwardly, which is particularly valuable for problems of large dynamic range, for instance when attempting to capture the details of scattering within the clumpy ISM of galaxies at the same time as the large scale cosmological radiative transport. We plan to build on the infrastructure developed in this work and implement adaptive refinement in future versions of the code. Here, we briefly outline the basic physics and methodology of the current implementation. Our notation follows that of Dijkstra et al. (2006).

There are three relevant reference frames: the frame of the observer, the fluid frame, and the frame of a particular atom, which may have thermal motion with respect to the fluid frame. The observer is assumed to be located at and to be stationary with respect to the simulation volume, apart from a possible Hubble flow that redshifts all the photons from the simulation volume by the same factor. We assume that the simulation volume is sufficiently small that its volume is essentially at a single redshift. A quantity in the observer’s frame is denoted in the fluid frame and in the frame of the atom.

The Monte Carlo photons are injected with a probability per unit volume proportional to the local Ly emissivity (§3). The photons are assumed to have energy exactly at the line center in the frame moving with the fluid and the emission is taken to be isotropic. The natural and thermal broadening of the emission line have a negligible impact on the results, as the resonant scatters rapidly erase their memory. The Ly photons are then each transported by simulating their resonant scatters.

The Ly optical depth through a HI column density for a photon at the line center is given by

 τ0=NHIσLyα(ν0)=NHIf12√πe2mecΔνD≈8.3×106(NHI2×1020 cm−2)(T2×104 K)−1/2, (C1)

where is the Ly cross section, is the Ly line center, is the oscillator strength of the transition, is the charge of the electron, is its mass, is the speed of light, and , with . It is convenient to locally define a dimensionless frequency . The optical depth for photons of arbitrary energy is then

 τx=τ0H(a, x)=τ0aπ∫∞−∞e−y2dy(y−x)2+a2, (C2)

where is the Voigt parameter and is the Ly Einstein A coefficient.

In the absence of perturbations, the absorption of a Ly photon is quickly followed (after s) by the reemission of another Ly photon of the same energy. Owing to the motion of the atom, the scatter is however not coherent in the observer’s frame. Defining to be the net velocity of the scattering atom, the frequency of a photon of incident frequency after scattering is

 xo=xi+va⋅(ki−ko)vth+g(ki⋅ko−1)+O((vthc)2), (C3)

where and are unit vectors in the directions of the incident and outgoing photons, respectively (e.g., Dijkstra et al., 2006). The Hubble expansion is modeled by including a component , where is the displacement vector from the center of the system and is the Hubble parameter, in the bulk velocity term. The “recoil” term accounts for the average transfer of momentum from the photon to the atom during the scatter and can be written as (Field, 1959). The effect of recoil is usually small (Adams, 1971) and is ignored in our calculations.

Because the optical depth near the Ly line center is very large in astrophysical conditions, a Ly photon typically travels only a short distance before being scattered. The numerous scatters result in a random walk in space. At the same time, the scatters cause the photon to oscillate in frequency and it usually escapes the medium during an excursion far into a damping wing of the line, where the optical depth is greatly reduced (e.g., Zanstra, 1949; Unno, 1952; Field, 1959; Osterbrock, 1962; Adams, 1972; Harrington, 1973; Neufeld, 1990; Gould & Weinberg, 1996). As a result, the emergent Ly line profile is heavily affected by the medium through which it propagates and therefore provides a probe of the latter.

To simulate the scatters of the Monte Carlo photons, we at each step randomly pick the optical depth before the next scatter. By definition, has a PDF . Given the frequency and propagation direction of the photon, we integrate through the medium until an optical depth has been reached; this defines the position of the next scatter. The frequency seen by a fluid element along the way in general differs from that seen by the observer and is given by the Doppler shift formula . The frequency of the outgoing photon is then calculated using equation C3, which requires picking and .

The direction of the outgoing photon is picked from a dipolar phase function, . Owing to quantum mechanical effects, this Rayleigh phase function is strictly valid only for scatters with (see Dijkstra & Loeb, 2008, and references therein), but the multiple scatters act to quickly randomize propagation directions so that the precise phase function is unimportant. An isotropic phase function, , gives identical results in practically all astrophysical conditions, aside for polarization, which we neglect here.

The first part of the scattering atom velocity is simply the local bulk velocity of the fluid, . To this, a thermal velocity must be added. In a frame moving with the fluid, but oriented such that the axis is parallel to , the thermal velocity can be expressed as ), where is the normalized component parallel to the direction of propagation, and are perpendicular components. The probability of scattering off an atom with parallel component is the product of the one-dimensional thermal distribution and the absorption cross section,

 P(u′||)=aπe−u′2||(x′i−u′||)2+a2H−1(a, x′i), (C4)

while are simply thermally distributed since the probability of absorption is independent of these normal components.

The procedure is repeated for each scatter of each Monte Carlo until it escapes the computational volume. In this work, we assume that all photons escape the medium but the algorithm can be straightforwardly extended to model photon destruction (e.g., by dust) by specifying a destruction probability at each integration step or scattering event.

### c.2. Numerical Aspects

Having described the essential aspects of the algorithm, we now elaborate on numerical aspects of its implementation relevant to its accuracy and performance.

The generation of random numbers is a key element of the Monte Carlo method. In most instances, these are straightforwardly generated using either a transformation method or a rejection method (e.g., Press et al., 1992). An important exception is the generation of , whose PDF (eq. C4) is not analytically integrable and for which there is not an obvious efficient bounding function to use in the rejection method. Since a variate from this distribution is generated at each scatter, it is important to have a fast implementation. We generalize the rejection method described in the appendix of Zheng & Miralda-Escudé (2002b) to use three critical velocities (in their notation, and in addition to ) in order to minimize the number of rejections.

Although the distance between neighboring grid points is , photon trajectories are traced exactly in the sense that the photon can be located at any coordinate in the box along its ray. This can greatly improve the accuracy of the calculations possible with limited computational resources. For instance, a large and/or very optically thick medium may be well described by a relatively coarse grid if it is fairly homogeneous. However, because of the short mean free path of the Ly photons, the radiative transfer would be poorly captured if the position of the photon were restricted to widely separated grid points. The integration along rays is still done in finite steps of length , but this parameter can be specified arbitrarily and in particular can be . Since the photon position is in general between the grid points on which physical properties are stored, these must be interpolated onto the ray. Our code can use either nearest grid point (NGP) or cloud-in-cell (CIC) interpolation (e.g., Hockney & Eastwood, 1988). While CIC is more accurate, NGP produces similar results but requires a factor of several fewer operations per integration step. We use NGP interpolation in the calculations presented in this work.

Finally, we use an “accelerated scheme” (e.g., Ahn et al., 2002) that greatly speeds up radiative transfer calculations by skipping scatters in the core of the Ly line, during which the photon is nearly stationary in space. Specifically, we define a critical frequency and say that a photon with is in the core of the line. Photons in the line core are then forced into the wings by allowing them to scatter only off high-velocity atoms. This is achieved by truncating the thermal Gaussian distribution from which and are drawn to zero for