Lyman Alpha Radiative Transfer in a Multi-Phase Medium

Hydrogen Ly-alpha is our primary emission-line window into high redshift galaxies. Surprisingly, despite an extensive literature, Ly-alpha radiative transfer in the most realistic case of a dusty, multi-phase medium has not received detailed theoretical attention. We investigate resonant scattering through an ensemble of dusty, moving, optically thick gas clumps. We treat each clump as a scattering particle and use Monte Carlo simulations of surface scattering to quantify continuum and Ly-alpha surface scattering angles, absorption probabilities, and frequency redistribution, as a function of the gas dust content. This atomistic approach speeds up the simulations by many orders of magnitude, making possible calculations which are otherwise intractable. With these surface scattering results, we develop an analytic framework for estimating escape fractions and line widths as a function of gas geometry, motion, and dust content. We show that the key geometric parameter is the average number of surface scatters for escape in the absence of absorption. We consider two interesting applications: (i) Equivalent widths. Ly-alpha can preferentially escape from a dusty multi-phase ISM if most of the dust lies in cold neutral clouds, possibly explaining anomalously high EWs seen in many high redshift/submm sources. (ii) Multi-phase galactic outflows. We show the characteristic profile is asymmetric with a broad red tail, and relate the profile features to the outflow speed and gas geometry. Many future applications are envisaged. [Abridged]


INTRODUCTION
The hydrogen Lyα line is our primary emission line window on the high-redshift universe. It is almost invariably crucial in securing redshift-identifications for the highest redshift galaxies (e.g., Hu et al. (2002a,b); Ajiki et al. (2003); Kodaira et al. (2003); Rhoads et al. (2003); Santos et al (2004)). Besides yielding redshifts, the shape of the line profile, equivalent width, and offset from other emission/absorption lines also encode information about the geometry, kinematics and underlying stellar population of the host galaxy. For instance, features in Lyα emission have been used to suggest strong galactic outflows (Kunth et al. 1998), as a signature of strong accretion shocks (Barkana & Loeb 2003), and as evidence for an unusually strong ionizing continuum, perhaps due to Pop III stars (Malhotra & Rhoads 2002). Even after escaping the environs of the host galaxy, Lyα photons undergo processing in the surrounding intergalactic medium (IGM), and the presence or absence of observed Lyα emission can be used to place constraints on the epoch of reionization (Fan et al 2002;Haiman 2002;Santos 2004). Because of the numerous factors which con-tribute to Lyα radiative transfer, the interpretation of such features is fraught with complexity. For instance, in a comprehensive set of ∼ 1000 Lyman Break galaxies at z ∼ 3, a plethora of Lyα strengths and line shapes were seen, ranging from pure damped absorption, to emission plus absorption, to pure strong emission (Shapley et al 2003). Because of the tremendous potential returns for interpreting some rich data-sets, it is crucial to strive for a more detailed theoretical understanding of Lyα emission line features.
A very important factor in Lyα transmission is the presence of dust. Since the massive stars which produce metals evolve on a short timescale, and indeed supersolar metallicities (Pentericci et al. 2002) and CO emission (Bertoldi et al. 2003) have been observed in the highestredshift quasars at z ∼ 6, dust is likely to be present in the ISM of even high-redshift galaxies. Because of their long scattering path-lengths, Lyα photons are extremely vulnerable to dust attenuation (Neufeld 1990;Charlot & Fall 1991), and it was thought that this could account for low observed Lyα equivalent widths compared to that expected from optical Balmer emission lines (Meier & Terlevich 1981;Hartmann et al. 1988), as well as c 0000 RAS early failures to detect high-redshift galaxies in blank sky surveys. However, further work has shown that dust content is not strongly correlated with Lyα equivalent width (where dust content can be inferred from metallicity or submillimeter emission). For instance, some dust-rich galaxies have significantly higher Lyα photon escape fractions than less dusty counterparts (Kunth et al. 1998(Kunth et al. , 2003. Indeed, Giavalisco, Koratkar, & Calzetti (1996) found a lack of correlation between the equivalent width of Lyα and the UV continuum slope β, which measures continuum extinction. They interpreted this as evidence for decoupling of the extinction of continuum and resonant line photons.
Such decoupling could take place if the interstellar medium is clumpy. Neufeld (1991) and Charlot & Fall (1993) emphasized the importance of the geometry and multi-phase nature of the ISM in affecting the observed Lyα line. In particular, Neufeld (1991) showed that in a clumpy, dusty ISM, the emergent Lyα emission could have a higher equivalent width than the unprocessed spectrum of the underlying stellar population. For instance, if the dust survives primarily in cold neutral clouds, Lyα photons scatter off the clouds and spend most of their time in the intercloud medium, whereas continuum photons propagate unhindered into the clouds and suffer greater extinction. Observationally, the ISM of our Galaxy is known to be clumpy down to small scales (Marscher, Moore, & Bania 1993;Stutzki & Guesten 1990) , with a power-law cloud mass spectrum based on CO (Sanders, Scoville, & Solomon 1985) and 21cm (Dickey & Garwood 1989) emission data. From IRAS 100µm, CO and 21cm data, there is evidence for a multi-scale fractal structure for both the diffuse HI clouds (Bazell & Desert 1988) as well as the molecular component (Elmegreen & Falgarone 1996). The clumpiness of the ISM is well-established and it must be taken into account in radiative transfer calculations.
Surprisingly, there has not been any detailed, quantitative, three-dimensional study of the effects of a dusty, clumpy, ISM on Lyα radiative transfer. The pioneering work of Neufeld (1991) was a semi-analytic calculation for a planeparallel slab: many issues, such as the detailed line profile and the effect of geometry, cannot be addressed with such an approach. Recently, Richling (2003) made a first attempt at a quantitative, three-dimensional calculation, but the slow convergence of the numerical technique employed restricted the study to line center optical depths of τ 100, corresponding to neutral hydrogen column densities of N 10 16 cm −2 for velocities v ∼ 100 km s −1 : orders of magnitude too low to be applicable to high-redshift galaxies. There have been many studies of the radiative transfer of UV continuum photons in a clumpy, dusty ISM, using a variety of techniques (e.g. Witt & Gordon (1996); Vársoi & Dwek (1999); Gordon et al (2001)), but none with extensions to resonance line photons. Conversely, while there have been Monte-Carlo radiative transfer studies of Lyα photons in both static media (Ahn, Lee, & Lee 2001 and expanding supershells (Ahn, Lee, & Lee 2003), all have only considered a uniform medium. This paper therefore represents a first attempt at numerically investigating Lyα radia-tive transfer incorporating both the effects of dust and gas clumping.
A key motivation is understanding recent puzzling observations of anomalous equivalent widths in high-redshift galaxies. For instance, high redshift z = 4.5, 5.7 sources observed by Rhoads et al. (2003) in the Large Area Lyman Alpha survey (LALA) show anomolously large Lyα equivalent widths of EW 150Å (rest-frame), many far in excess of any known nearby stellar population. An AGN origin is unlikely, as the observed upper limit on the X-ray to Lyα ratio is about 4-24 times lower than the ratio for known type II quasars (Malhotra & Rhoads 2002). The radiative transfer effects studied in this paper can produce an anomolously large Lyα equivalent width from a standard stellar population. Such an effect could also be at work in the mysterious Lyman-alpha emitters observed at z ∼ 3.1 by Steidel et al (2000), which have enormous Lyα fluxes of ∼ 10 −15 erg s −1 cm −2 (a factor ∼ 20 − 40 times large than typical line emitters at the same redshift), but no observed continuum. Finally, our calculation could be of particular interest in interpreting the large (∼ 1000) sample of Ly-break galaxy spectra (Shapley et al 2003), as well as understanding the spectra of galactic starbursts with winds.
The outline of this paper is as follows. In §2, we derive the basic multiphase Lyα scaling relations. We then consider radiative transfer off opaque gas surfaces in §3, describing the Monte Carlo simulations and obtaining fitting formulae for the absorption probability, angular and frequency redistribution functions for both continuum and resonant scattering. With these surface scattering formulae in hand, we then develop a framework for multi-phase radiative transfer in §4, where we derive escape fractions and Lyα line widths, discuss the role of the gas geometry, analyze several geometries of astrophysical relevance, and discuss the effects of dilute gas in between the opaque clumps. The surface scattering formulae substantially reduce the computational cost of simulating Lyα transfer, making otherwise intractable calculations feasible. We also develop a simple analytic model with a single geometric parameter that shows good agreement with the full Monte-Carlo simulations. In §5, we discuss some applications of our formalism. We show how preferential absorption of continuum photons can lead to strong enhancement of the Lyα equivalent width. We also consider the typical Lyα line profiles resulting from outflows/inflows of multiphase gas, and relate the profile characteristics to the outflow/inflow speed and the gas geometry.

SCALING RELATIONS FOR Lyα ABSORPTION
In this section we build some physical intuition, by making simple order-of-magnitude estimates for the absorption of Lyα photons in both homogeneous and multi-phase media, and summarizing some of the most important results from §3 and §4. We shall see that for conditions prevailing in most galaxies, Lyα photons cannot escape unless the medium is multi-phase. Before beginning, it is useful to define some terms. Let ν dop = (V dop /c)ν0 be the line Doppler width, where ν0 is T 4 ≡ T /10 4 K, xHI is the hydrogen neutral fraction, ν 0 = 2.48 × 10 15 s −1 is the Lyα line center frequency, ν dop = (V dop /c)ν 0 is the doppler frequency, ν L = 4.03 × 10 −8 ν 0 is the width of the quantum broadening Lorentzian profile, and θscat is the angle between the incident and outgoing photon directions.
the Lyα line center frequency, and V dop = (2kBT /mp) 1/2 is the characteristic atomic velocity dispersion times √ 2. We evaluate the frequency shift from line center in Doppler units, x ≡ (ν − νo)/ν dop 1 . The Lyα scattering cross section is σ(x) = σ0Φ(x), where Φ(x) is the Voigt function, which is characterized by a Gaussian Doppler core, and Lorentzian damping wings due to quantum broadening. For frequencies in the line wing, |x| > 3, the Voigt function is dominated by the Lorentzian: Φ(x) ≈ a/( √ πx 2 ), where a ≡ νL/2ν dop = 4.72 × 10 −4 T4 −1/2 , νL = 4.03 × 10 −8 ν0 is the width of the Lorenztian profile, and T4 ≡ T /10 4 K. For ease of reference, we have listed the most common radiative transfer parameters used in this paper in Table 1.

Homogeneous Slab
We begin by reviewing the physics of radiative transfer of Lyα photons through an optically thick slab, a problem that was first correctly solved by Adams (1972), and subsequently verified and explored in much greater detail (Harrington 1973;Hummer & Kunasz 1980;Bonilha et al. 1979;Frisch 1980;Neufeld 1990;. We use these classical results to test our Monte-Carlo code in Appendix A.

Homogenous Slab: Dust-Free
Consider a Lyα photon escaping from a dust-free slab of pure HI with line center optical depth τ0. When the photon is the Doppler core, its mean free path is very short, and it barely diffuses spatially. It is always scattered by atoms with the same velocity along its direction of motion as the atom that emitted it. On rare occasions, it will encounter a fast moving atom in the tail of the Maxwellian velocity distribution, with large velocities perpendicular to the photon's direction. When this photon is re-emitted, it will be far from line-center, where the slab is optically thin. For a line-center optical depth of τ0 = 10 3 , a frequency shift of x ≈ 2.6 is sufficient to render the slab optically thin, τ ≈ τ0e −x 2 ≈ 1, and the photon can escape. So escape from the medium is dominated by rare scattering events. However, if the medium is sufficiently optically thick, τ0a > 10 3 , the non-negligible optical depth due to the damping wings still prevents escape. In this case, the photon will suffer repeated scatterings in the Lorentzian wings of typical atoms, and diffuse slowly in space and frequency, executing a random walk. Each scatter induces an r.m.s. Doppler shift of order x ∼ 1, has a mean Doppler shift per scatter of −1/|x| (with a bias to return to line center, due to the large probability for photons to scatter in the core ; Osterbrock (1962)), and transverses an optical depth τ0Φ(x)∆τ ∼ 1, or a mean free path which is ∆τ ∼ 1/Φ(x) line-center optical-depths. Hence, a photon at frequency |x| ≫ 1 returns toward line center after Ne ∼ x 2 scatterings, having travelled an r.m.s. optical depth τrms ∼ √ Ne∆τ ∼ |x|/Φ(x). If on its single longest excursion, the photon diffuses an r.m.s. distance of order the system size, τrms ∼ |x|/Φ(x) ∼ τ0, then the photon can escape. Since Φ(x) ∼ a/x 2 , this implies a critical escape frequency: or almost ∼ 400N 1/3 21 T4 1/6 km s −1 away from line center, where N21 ≡ NHI /(10 21 cm −2 ). This displacement of photons away from line center can be seen in our Monte-Carlo simulations in Fig. A1.

Homogeneous Slab: Dusty
Now let the gas contain dust, with an total (scattering + absorption) interaction cross section per hydrogen atom of σ d , and an absorption probability per dust interaction ǫ d . The average absorption probability per interaction with either dust or hydrogen is: where xHI is the hydrogen neutral (HI) fraction (which must be introduced because σ d is the cross-section per hydrogen nuclei). In the third step we used the fact that, except very far from line center, HI scattering dominates, xHIΦ(x)σ0 ≫ σ d , and defined the absorption parameter where σ a ≡ ǫ d σ d . For the diffuse HI phase of the Milky Way, σ a −21 ≡ σ a /10 −21 cm 2 /H ≈ 1, ǫ d ≈ 0.5, and xHI ≈ 1 (Draine & Lee 1984;Whittet 2003;Draine 2003), and so β ≈ 10 −8 . The fourth step in eqn. (2) uses the wing photon Under what conditions can the Lyα photon escape from such a dusty medium? While this has been the subject of detailed analytic and numerical work (e.g., Hummer & Kunasz (1980);Frisch (1980); Neufeld (1990)), we can understand the basic scaling laws quite easily. The probability that a photon will be absorbed at a given frequency x is simply the number of scatterings at that frequency times the probability of absorption per scattering: Thus, P slab abs (x) ∼ 1 for This implies that a photon will be absorbed before escape if it has to diffuse far into the line wings in order to escape from the slab, or if xe > x abs , with xe and x abs given by eqns. (1) and (5)  photons cannot escape from the medium. This simple criterion is borne out by more detailed calculations (e.g., see Fig. 5 of Ahn, Lee, & Lee (2000), and references therein). In terms of the HI column density N21, Lyα photons cannot escape from a homogeneous dusty slab once Since typical HI column densities in the Milky Way and other galaxies is N21 ∼ 1, Lyα photons could not escape if most of the HI is in a homogeneous dusty slab.
In the next section, we see that if the gas is instead inhomogeneous/multi-phase, Lyα photons can escape much more easily.
2 Note that at x abs , the absorption probability per interaction is still small, ǫ ≪ 1, and scatterings still strongly predominate. The scattering and absorption cross-sections are comparable only at a much larger frequency, x(ǫ ∼ 0.5) ∼ (a/β) 1/2 ∼ x 2 absorb ≫ x absorb , by which time all photons have been absorbed.

Multiphase Gas
We shall now estimate the absorption criteria for a multiphase dusty HI distribution. It is worth first noting that a medium is always more transparent when it is clumpy, for fairly generic and model independent reasons. The effective optical depth in an inhomogeneous medium is τ clumpy = −ln ( exp [−τ ] ), where the average is over all lines of sight. However, for a uniform medium, τ = constant along all lines of sight, so that τ uniform = τ . From the standard triangle inequality, and applying the negative logarithm to both sides, we see that τ clumpy τ uniform . Thus, for instance, flux transmission in quasar absorption spectra is increased for an inhomogeneous intergalactic medium (IGM), where transmission is dominated by underdense voids (e.g., Fan et al (2002)). This effect is strongly exacerbated if most of the absorbing material lies in dense clumps which are optically thick to scattering. In this case, most of the photons scatter off the cloud surfaces without penetrating the clouds, which effectively shields the absorbing material. This situation naturally arises in a multi-phase ISM, when most of the dust lies in dense molecular/atomic clouds. For now, let us assume that the inter-cloud medium is highly ionized and relatively dust-free, so that all of the dust and HI lies in dense clouds.
In §3.3, we show that we can calculate analytically the escape probability of Lyα photons in a multi-phase medium quite accurately, given just two parameters: N0 and ǫc. We define N0 as the mean number of cloud surfaces a photon would encounter before escape in the absence of absorption. It only depends on the geometry of the multi-phase medium (the trajectory of photons is independent of frequency, provided clouds are very optically thick). For most of the cases we will consider, N0 ∼ 1 − 30, with typical values N0 ∼ 5. The cloud albedo ǫc is the probability of absorption upon hitting a cloud surface. In §3.3, we show that: where ǫ is the absorption probability per interaction given by eqn.
(2). Eqn. (9) is easily understood in the case where ǫ is constant (e.g., for coherent scattering). The effective absorption optical depth of a medium with scattering is τ * ≈ τa(τa + τs) (e.g., Rybicki & Lightman, 1979, p. 38), where τa and τs are the absorption and scattering optical depths, respectively. Hence, the albedo is ǫc = τ * /(τs +τa) ≈ τa/(τa + τs) = √ ǫ. In §3.3, we show this scaling still holds for Lyα photons, despite the fact that ǫ(x) changes as the photon random walks in frequency whilst scattering within the cloud. We find that the typical frequency shift after scattering off an optically thick surface is ∆x ∼ 1.5, with most of the redistribution in a symmetric profile about the incident frequency xi. With ǫ(x) given by eqn.
Under the above approximations, the probability that a Lyα photon at frequency x will be absorbed is: which should be compared against eqn. (4) for a homogenous slab. Notice the much weaker scaling with frequency: P multi−phase abs ∝ x, instead of P slab abs ∝ x 4 . From eqn. (10), P multi−phase abs ∼ 1 for frequencies Comparing against eqn. (5), the cut-off absorption frequencies for the slab and multi-phase case are actually comparable, modulo the value of N0. However, there is an important difference: Lyα photons have to diffuse far into the wings in order to diffuse spatially out of an optically thick slab. Since escape requires xe > x absorb for a very optically thick slab, the photons will inevitably be absorbed. By contrast, there is generally much less diffusion into the line wings when scattering off surfaces in a multi-phase medium. Photons typically only penetrate small optical depths, τ ∼ 1 − 10 in the cloud surfaces before escaping, and the number of scatterings is much less. Thus, the majority of photons need not necessarily stray far from line center. Clearly, the crucial parameter which determines if Lyα photons can escape from a multi-phase medium is xe, the characteristic escape frequency; this must be small, xe < x absorb , for photons to escape. Lyα photons in a multiphase medium acquire Doppler frequency shifts in two ways: through the thermal motions of HI atoms, as before, and also through the bulk motions of the clouds/scattering surfaces. For this reason, it is useful to rewrite x abs in units of velocity: where V2 ≡ V /100 km s −1 . In §4.3, we show that atomic motions cause a net r.m.s. frequency shift of ∼ 1.5N0V dop after N0 surface scatterings, and consequently result in an escape frequency of: By contrast, we find that cloud motions (either random motions or bulk inflows/outflows) with characteristic velocity V c causes a net r.m.s. frequency shift of: Note that V e,atomic ∝ N0, while V e,clouds ∝ √ N0, which we discuss in §4.3. When both the frequency redistribution and surface motion are combined, we find that the typical escape velocity is simply given by the sum (rather than the sum in quadrature), For absorption to be important, we require V e > V absorb , which constrains N0 to be If N0 satisfies this inequality then the Lyα photons will be significantly absorbed. This approximate constraint has the correct limits when either V e,atomic = 0 or V e,cloud = 0, but is off by a factor of ∼ 2 when V e,atomic ∼ V e,cloud . The geometric parameter N0 thus plays a key role in determining if Lyα photons can escape in a multi-phase medium, and plays an analogous role to the column density NHI in a homogeneous slab. In §4.2, we provide formulas for N0 for various different basic types of multi-phase geometries.

SURFACE SCATTERING & ABSORPTION
Multi-phase radiative transfer typically involves photon propagation through an optically thin inter-cloud medium, and repeated scattering off optically thick clouds. In a fullblown Monte-Carlo simulation, the latter consumes by far the lion's share of computational time. This is extremely inefficient: the same scattering/absorption problem off cloud surfaces is being solved over and over again for each photon. A better approach is to consider each cloud as a scattering/absorbing particle with its own radiative transfer properties (for other applications of this viewpoint, see Neufeld (1991); Hobson & Padman (1993); Vársoi & Dwek (1999)). We characterize these cloud scattering properties in this section. For Lyα photons, clouds are extremely optically thick and have essentially the same radiative transfer properties as a semi-infinite slab. This eliminates detailed dependence on the geometry of the cloud: all that matters is its dust content and the initial photon frequency. Surprisingly, the radiative transfer properties of a dusty semi-infinite slab to Lyα photon scattering has not been characterized in detail. We do so in this section. We derive formulas for the net absorption probability (the "cloud albedo") ǫc, the exiting photon angular distribution D(θ), and the exiting photon frequency redistribution R(xi, x), as a function of the initial photon frequency xi and the gas composition. With these surface transfer formulae, radiative transfer through regions containing opaque gas clouds can be quickly estimated and/or simulated without performing any scattering calculations within the individual gas clouds. This allows for both vast speed-ups of Monte-Carlo simulations (outlined in §4.6) and a tractable analytic multiphase radiative transfer analysis ( §4). If a photon typically scatters N times before exiting a cloud (where N ∼ 10 5−7 for incident frequencies xi ∼ 5 and σ a −21 ∼ 1-see §3.2.4), then this allows a speedup of order ∼ N , making tractable multi-phase calculations which would otherwise be prohibitively expensive.
Our approach is to find fitting formulas to Monte-Carlo simulations of an ensemble of incident photons. Whenever possible, we base the fits on known analytic formulas for simple cases, extending the analytic formulas to encompass the more general cases that we simulate. We begin by describing the Monte-Carlo algorithm we use. Surface radiative transfer of continuum photons, where scattering by dust is effectively coherent, is discussed next. We then consider the more complex case of Lyα surface scattering, where resonant frequency redistribution effects must be dealt with. Lastly, we discuss two kinematic aspects of surface scattering: the average scattering angle, and the frequency shift due to a bulk surface velocity.

Monte-Carlo Code Description
The Monte Carlo algorithm we use is similar to the code used by Ahn, Lee, & Lee (2001 and Zheng & Miralda-Escudé (2002). These papers provide a fuller description of the algorithm than that which we give here. An ensemble of photons is run through a medium with neutral HI and dust, and statistics are gathered. Each photon is tracked until it either escapes the medium or is absorbed, at which point the photon's flight is terminated. The optical depth τ between each interaction is drawn from the distribution exp(−τ ), i.e. τ = − ln u where u ∈ [0, 1] is a random variable drawn from a uniform distribution (hereafter "univariate"), and the photon's position is updated. The photon then interacts with either the HI or the dust, resulting in either HI resonant scattering, dust scattering, or dust absorption, all of which we describe next.
We model the dust as particles which can either absorb or coherently scatter photons. Although dust scattering is not necessarily coherent, in practice ignoring frequency redistribution due to dust is an excellent approximation. The trajectory of a continuum photon is unaffected by small deviations from coherent scattering, since the dust albedo ǫ d only varies weakly with frequency. However, the Lyα absorption probability is a strong function of frequency (see eqn. (2), so for resonant scattering the effects of frequency redistribution must be taken carefully into account.
A continuum photon only interacts with dust, with an absorption probability ǫ d per interaction. We determine if the photon is absorbed during a given interaction by drawing a random variable u ∈ [0, 1] from a uniform distribution; if u ǫ d , the photon is deemed to be absorbed. If the continuum photon is not absorbed, then its direction is changed by the scattering angle θscat off the incident direction, with a random azimuthal angle. We use the Henyey-Greenstein scattering angle distribution 3 (e.g., (Witt 1977)) which is parameterized by the dust scattering asymmetry parameter g d ∈ [−1, 1], and where we use the normalization 1 = π 0 dθ P (θ). The scattering asymmetry parameter is defined as g ≡ cos θscat . To approximate dust absorption and scattering in the Milky Way (Draine & Lee 1984; 3 Draine (2003) shows that the Henyey-Greenstain distribution is inaccurate for wavelengths λ 4700Å. However, since the surface scattering problem we consider is essentially planar, the details of the dust scattering distribution should significantly affect the results. Witt & Gordon 1996;Whittet 2003;Draine 2003) at wavelengths near 1216Å, we use ǫ d = 0.5 and g d = 0.5, unless otherwise noted.
If the photon is a Lyα photon, then the interaction can either be with dust or with neutral HI. The probability of a dust interaction is σ d /(σ d +xHIΦ(x)σ0); if a random univariate u is less than this, then the interaction is identical to the "continuum" dust interaction described above. Otherwise, the Lyα photon scatters resonantly off neutral hydrogen. In all our simulations we take the hydrogen in the cold phase to be completely neutral, and so adopt xHI = 1 unless otherwise noted. Although the velocity distribution of the hydrogen is Maxwellian, the velocity distribution of atoms that scatter photons depends upon the frequency of the photon. Letn be the direction of the photon before scattering. In the two directions perpendicular ton, the scattering atom's velocity distribution is Gaussian, where w ⊥ = v ⊥ /V dop and v ⊥ is a velocity component in one of the two transverse directions ton. In the direction parallel ton the velocity distribution of scattering atoms is a Gaussian weighted by a Lorentzian: the Gaussian is due to the thermal motion of the gas, while the Lorentzian is due to the increased probability for scattering in the wings from quantum mechanical broadening. The velocity of the atom vz along the direction of the incident photon is determined by drawing a random variable from the distribution: where wz = vz/V dop (see Zheng & Miralda-Escudé (2002) for a rapid algorithm for generating random numbers with this distribution) and xi is the photon's incident frequency.
In the rest frame of the atom, the frequency of the outgoing photon is the same as the incident frequency (strictly speaking it differs slightly due to the recoil effect (Field 1959), but for our purposes this is negligible). The new directionn ′ is given by a dipole distribution, with the symmetry axis defined by the incident directionn: where θ is the polar angle off the directionn. Although resonant scattering can result in either isotropic or dipole scattering angle distributions, depending upon the intermediate excited quantum state (Stenflo 1980;, the difference is immaterial for calculating spectra and escape fractions; a more careful treatment would be required, for instance, to accurately simulate Lyα polarization. Given the new photon directionn ′ and the scattering atom Doppler velocity w, the new photon frequency x ′ is given by

Avoiding Core Scatters
Lyα photons spend most of their scatters in the line core, where spatial diffusion is typically negligible. Essentially, each time a photon enters the line core it scatters in place until it is scattered by a high speed atom which moves the frequency out of the core. The frequency at the core-wing boundary, xc ≈ 3, is defined by: Hence, it typically takes exp(xc 2 ) ∼ exp(9) ∼ 10 4 scatters to scatter out of the core. Note that xc depends only logarithmically on the gas temperature, through the damping parameter a. For a photon that starts out at frequency xi in the line wing, the photon typically returns to the core fairly quickly, after ∼ xi 2 = 25(xi/5) 2 scatters (see discussion in §2.1). Consequently, most of the simulation time is spent calculating core scatters. By circumventing the core scatters, the simulation can be greatly sped up. We have adopted a scheme to do so that is similar to that used by , with the addition that we also consider absorption. One might think that absorption whilst scattering in the line core is negligible, due to the small physical path lengths traversed whilst scattering in the core. We confirm this quantitatively below.
For a photon with an initial core frequency xi where |xi| < xc, letnc be the average number of scatters for the photon to leave the core,xc ≡ |x| be the average frequency (absolute magnitude) while in the core, andxw ≡ |xw| be the average frequency (absolute magnitude) of the first scatter that leaves the line core, xw > xc. We ran simulations for gas at 10 4 K and adopted xc = 3. We find that for any initial frequency in the core,nc ≈ 2.9 × 10 4 ,xc ≈ 0.57, and xw ≈ 3.3 with an equal probability of leaving the core at x = 3.3 and x = −3.3. Thus, the probability of absorption during core scatters can be approximated by (see eqn. (2)): Since the typical photon scatters N ∼ 10 times before escaping the surface ( §3.3.4) and it takes N core ∼ 9(xi/3) 2 scatters to reach the core, a typical photon injected in the line wing will visit the core perhaps once. The probability of a Lyα photon being absorbed in the core during the surface scattering is, therefore, ∼ 10 −4 σ a −21 (ǫ d /0.5), which is negligible.
We therefore devised the following acceleration scheme 4 : 1) If a photon starts off in the line core, we do the exact core scattering, and only employ the approximation scheme on subsequent visits to the core. This is computationally cheap, since an incident core photon does not penetrate deep into the surface, and typically leaves after a few scatters. 2) If a photon enters the line core from the wing, the probability of absorption is ǫcore. 3) If the photon is not absorbed, then it is given a wing frequency x = ±xw, with an equal probability for plus and minus. The spatial position of the photon is exactly the same as where it entered the line core, and the new angular direction is randomly drawn from an isotropic distribution. In §A2 we compare this accelerated scheme to exact simulations of surface scattering. In practice, it gives accurate results, and gives a vast speed-up of the simulations, typically of order ∼ 10 5 .

Surface Scattering of Continuum Photons
In this section, we study the properties of coherent surface scattering of continuum photons. It is very useful to understand the properties of coherent scattering surface transfer in order to have a baseline for comparison with Lyα surface transfer. As such, in this section we do not use the Henyey-Greenstein scattering angle distribution, eqn. (18), but instead use the same distribution that we use for Lyα scattering, which is the dipole distribution, eqn. (21) 5 . We begin by calculating how thick a slab of gas must be before the surface scattering approximations apply. We then describe fits for the cloud albedo ǫc, the exiting photon scattering angle distribution D(θ), and the typical number of scatters N .

The Surface Approximation
For a slab of material with a finite optical thickness τ , the radiative transfer of photons incident on a surface will be approximately the same as for a semi-infinite slab if the fraction of photons that are transmitted through the slab, f T , is small. In this limit, the surface is not translucent but acts as an absorbing mirror, and all photons are either reflected or absorbed. We define the penetration column density N pt such that when N > N pt the transmitted fraction is less than 10%, f T < 0.1. From a series of Monte-Carlo simulations, we find that a decent fitting formula for f t is as shown in Figure 1. The transmission will be negligible (f T 0.1) when either τ ǫ 1/2 > ∼ 3 or τ > ∼ 9. The corresponding penetration column density is (26)

Surface Absorption
To derive a formula for the cloud albedo ǫc, we ran simulations for an isotropic surface source and averaged the absorption over this ensemble. For one dimensional radiative transfer, an exact formula for ǫc can be derived for photons incident on a semi-infinite line of material, where each scatter is front-back symmetric (g = 0) 6 . As shown by Figure 2, eqn. (27) provides a very good fit for the three dimensional, semi-infinite plane case. When ǫ ≪ 1, we find ǫc ≈ 2 √ ǫ, which is similar to the power law found by Neufeld (1991).

Escape Angles
To find a fit for the distribution of exiting angles for photons that escape, we ran simulations for various incident angles θi relative to the surface normal. We find that for nearly all θi, and for both isotropic and dipole single particle scattering, the distribution of exiting angles θ is well fit by as shown in Figure 3. This distribution can be understood as the combination of two effects. First, if photons scatter multiple times before escape, the photons effectively lose all "memory" of the incident angle. This effect leads to a random exiting angle distribution, Dss(θ) ∝ sin θ. Second, photons that exit with angle θ will be attenuated if the optical depth traversed during the exiting leg, τ , exceeds unity. Let τ ⊥ be the perpendicular optical depth at the point of last scatter for a photon that would exit with an angle θ in the absence of absorption. The condition τ 1 implies a maximum perpendicular depth of such a photon is τ ⊥ max = cos θ (e.g., only shallow surface layers contribute photons escaping nearly parallel to the surface). Near the surface, the mean (Rybicki & Lightman 1979), pg. 320), and from our simulations we find that the prefactor is 2 for 1-D scattering. intensity is approximately constant for τ ⊥ < ∼ 1. This implies that the number of photons available to escape at an angle θ scales as ∝ τ ⊥ max , which implies Dss(θ) ∝ cos θ. Including both the effect of randomization and attenuation gives a distribution Dss(θ) ∝ sin θ cos θ, which, when normalized over θ ∈ [0, π/2], gives eqn. (28). In the case of dipole scattering there are deviations from this fit for grazing incident angles, but overall this fit is generic for surface scattering when the single particle scattering distribution is front-back symmetric. We also find that there is very little dependence on the absorption albedo ǫ. When the distribution Dss(θ) holds, the distribution of azimuthal angles is uniform over the interval [0, 2π). Finally, we note that if the slab is viewed at an angle θ obs (from the outward surface normal), then the observed intensity is ∝ Dss(θ obs ) cos θ obs . The extra cos θ obs factor is due to the dependence of the projected surface area on the viewing angle.

Number of Scatters for Escape
In §4.1 we present a general derivation of the average number of scatters, N , as a function of the escape fraction fe and the absorption albedo ǫ, given by eqn. (55). Application of this formula to surface scattering, where the escape fraction is fe = 1 − ǫc(ǫ), gives which is shown by the solid line in Figure 4. As the absorption albedo goes to zero, the average number of scatters diverges, although the median number of scatters does not seem to increase beyond ∼ 5. Clearly, the average is dominated by the rare photons that wander deep into the surface.

Surface Lyα Transfer
As in the case of coherent scattering considered above, we begin by calculating how thick a finite slab of gas must be before the surface scattering approximations apply. We then describe fits for the net absorption ǫc, scattering angle dis-tribution D(θ), the typical number of scatters, and the Lyα frequency redistribution R(xi, x) as a function of the incident frequency xi.

The Lyα Surface Approximation
The criterion for the surface scattering approximation to apply for incident Lyα photons is similar to that defined for coherent scattering. Consider Lyα photons with frequency xi incident on a finite slab with optical thickness τi = τ0Φ(xi), evaluated at the incident frequency. When xi is in the Lorentzian wing and if the slab is pure HI, Neufeld (1990) analytically derived that the fraction transmitted is: when τi ≫ 1. Our simulations confirm this result when there is very little dust, but find that f T can be substantially less than f T HI when a small amount of absorbing dust is present. To derive a fitting formula for f T we ran simulations for various incident frequencies and dust absorption cross-sections. To a large degree, f T depends only upon the incident single scattering albedo ǫi and the incident slab optical depth τi; almost all the frequency dependence is captured by these two parameters. This is extremely convenient: the transmitted fraction (and as we shall subsequently see, the reflected fraction) depends in a fairly simple way on the properties of the slab and the incident frequency. One might worry that due to frequency redistribution (which can be substantial; see §3.3.5), the frequency dependence becomes extremely complicated, but that appears not to be the case. In particular, the typical incident photon never "loses memory" of its initial frequency. Most photons that escape do so after a handful of scatters, ∼ 10. When xi > ∼ 4, the majority of photons do not wander into the core before escaping, and so most photons retain some memory of the incident frequency. As shown in Figure 5, a reasonable fit for photons initially in the line wing is: When ǫi → 0, the pure HI formula f T HI , eqn. (25), is recovered when τi ≫ 1. As can be seen in the figure, for photons initially in the Doppler core the transmitted fraction is slightly larger than this. The transmission will be negligible if either τi √ ǫi > ∼ 3 or τi > ∼ 12. This corresponds to a penetration column density for any incident frequency |xi| 3, with N pt 21 substantially smaller for photons in the Doppler core.

Surface Lyα Absorption
To derive a fit for the surface albedo ǫc, we ran again simulations for a wide range of incident frequencies and dust cross-sections. As in the coherent scattering case, we average ǫc over an isotropic incident direction. We find that ǫc mainly depends upon the single scattering albedo at the incident frequency, ǫi. As shown by the solid line in Figure 6, ǫc is well fit by which has the same form as for coherently scattered photons, eqn. (27). A slightly better fit is shown by the dashed line in Figure 6, Since eqn. (34) gives a slightly better fit only when ǫc is negligibly small, in practice we always use eqn. (33); this will simplify the subsequent multiphase analysis somewhat, for only a small loss in accuracy. Either of these formulas for ǫc applies even when xi is in the Doppler core (although eqn. (34) gives the better fit in this case). As long as the Lyα photon is in the line wing, the absorption probability is independent of the gas temperature. To show this, in Figure 7 we calculate ǫc using eqn. (33), as a function of the incident velocity shift ∆V in physical units, rather than Doppler units x. The figure shows the calculation for gas at temperatures T = 100 K and T = 10 4 K and for varying dust content σ a −21 . For photons in the line wing, the temperature dependence of ǫ drops out:  (1) The surface albedo ǫc is plotted against the single scattering albedo evaluated at the incident frequency, ǫ i , for gas at T = 10 4 K with dust parameters g d = 0.5 and ǫ d = 0.5. The inset shows an enlargement of the ǫ = 0.01 − 1 region. Simulations were run for four different dust absorption cross-sections and seven different incident frequencies. Each shade corresponds to a different σ a : from lightest to darkest σ a −21 = 0.01, 0.1, 1, 10. For each value of σ a , seven different incident frequencies x i were run: from right to left, x i = 20, 10, 5, 4, 2, 1, 0. Note that for x i = 0, the value of ǫ i for σ a −21 = 0.01 and σ a −21 = 0.1 is less than 10 −9 , and hence lies off the plot. The absorption probability ǫc depends only on ǫ i ; at fixed ǫ i , there is no independent variation with dust content σ a Focusing just on the temperature, a ∝ 1/ √ Therefore the combination x 2 /aσ0 has no dependence on the gas temperature. Physically, scattering in the Lorenztian wing is dominated by the quantum broadening of the cross-section, which does not depend upon the thermal motion of the scattering atoms. Since we have shown that ǫc mainly depends upon ǫi, it follows that ǫc is also temperature independent. Thus, the surface absorption probability does not depend upon the gas temperature

Escape Angles
As argued in §3.2.3 the scattering distribution Dss(θ), eqn. (28) is fairly generic when the single particle scattering is front-back symmetric (g = 0). Since Lyα scattering-which is either isotropic or dipole-is always front-back symmetric, the Lyα surface scattering angle distribution is also given by Dss(θ), independent of the incident angle θi.

Number of Scatters for Escape
The average number of scatters for Lyα photons to escape is more complicated than for continuum photons mainly because Lyα can be trapped in the Doppler core. As discussed in §3.1.1, Lyα photons in the Doppler core must scatter  (2) Simulation results for the net probability of surface absorption plotted as a function of the incident Lyα frequency shift off line center, in velocity units. The solid lines are for a gas temperature T = 10 4 K, while the dashed lines are for T = 10 2 K. We also show the effect of varying the dust content: from the thickest line to the thinnest, σ a −21 = 100, 10, 1, 0.1, 0.01, 0.001. The cloud albedo ǫc is independent of gas temperature if the photon is in the line wing, though the core-wing boundary is obviously temperaturedependent: Vc ≈ 3V dop ∝ T 1/2 (the vertical dashed lines show the core-wing boundary at T = 100 K, 10 4 K). The lines flatten out once the incident frequency shift lies within the Doppler core, because the Lyα scattering cross section approaches its line center value σ 0 . This causes the single scattering albedo to asymptote as V → 0. Since ǫc ∝ √ ǫ approximately holds even for incident photons in the Doppler core, the surface albedo will also asymptote as V → 0.
∼ 3 × 10 4 times before a rare scattering event brings the frequency into the line wing. In Figure 4 we show exact Monte-Carlo simulations of Lyα photons, where the core scatters are directly calculated, i.e. the Monte-Carlo acceleration scheme described in §3.1.1 is not used. The huge increase in scatterings over the continuum scattering result is obviously due to scattering in the Doppler core. The average number of scatterings in the line wing is (analogous to the continuum scattering result) while the average number of scatterings required to reach the Doppler core is ∼ 25[xi/5] 2 (see the discussion in §2.2). Thus, although the typical photon does not reach the core, the probability of reaching the core is large enough that the the average N is dominated by core scattering.

Surface Lyα Frequency Redistribution
In this section, we find a formula for the reflected frequency distribution R(xi, x) as a function of the incident frequency xi, for the pure dust-free HI case. When dust is added, we find that the distribution R dust (x) adheres closely to R(xi, x) as long as σ a −21 xHI < ∼ 10. Dust will have little ef- The same as the middle panel, but shown on a logarithmic scale and over a larger frequency range. In all three plots, the sharp peaks at x ≈ ±3 are an artifact of the acceleration scheme that we use; we do not calculate any core scatters, and photons which escape the core are all placed at x = ±3.3. Exact simulations indicate that there is indeed a "pile up" of photons just outside the core, but the peaks are not as sharp. In any event, the number of photons in the peaks are relatively insignificant.
fect on frequency redistribution, except for extremely dusty or highly ionized clouds. For photons incident at frequency xi on an optically thick slab, aτ0 > ∼ 10 3 , an analytic solution for the transmitted and reflected emission profile has been derived by Neufeld (1990), extending the earlier work of Harrington (1973) who obtained these results for the case xi = 0. By taking the τ0 → ∞ limit of eqn. (2.33) in Neufeld (1990), we derive the analytic result for the reflected spectrum (normalized to unity): When compared to simulations, R anlyt (x) is inaccurate in two respects. First, the actual peaks are shifted by an amount −2/xi (towards line center) compared to R anlyt . This can be compensated for by using a shifted incident frequencyxi in place of xi, wherẽ The same as the top panel but on a logarithmic scale and over a larger frequency range. Over the frequency range 0 < x < 15, the frequency redistribution function R(x, x i ) differs by < 20% between 0 < σ a −21 < 1; hence, to lowest order it is independent of dust content for the cases we are interested in.
Second, the peaks are significantly flatter than those given by R anlyt . It is not surprising that R anlyt is inaccurate, since the analytic results of Neufeld (1990) only hold when the photons traverse a large line center optical depth, aτ0 > ∼ 10 3 . The bulk of photons that reflect off a surface do not traverse such a large optical distance, and so the analytic result may not be accurate, which seems to be the case. However, one can obtain an excellent approximation to the redistribution function by modifying eqn. (37) to include an additional fitting parameter α (see Appendix B for details): Note that R(xi, x; α) is guaranteed to be normalized for all α. For example, R anlyt is recovered with α = 6 andxi = xi. We find that our simulations for pure HI are well fit by α = 45/2, withxi given by eqn. (38). To generate random frequencies x that obey the probability distribution R(xi, x; α), one draws a random univariate u ∈ [0, 1], and sets The frequency x is then given by functional inversion, x = F −1 (u). Carrying out these steps on eqn. (39), done in Appendix B, we find that the exiting frequencies x which obey the probability distribution R(x; α) can be generated by the equation where u ∈ [0, 1] is a random univariate.
When dust is included, the profile peaks become slightly sharper and the tails fall off slightly faster. However, as shown in Figure 9, the pure HI distribution closely matches the dusty distribution when σ a −21 < 10. In practice, we shall therefore always adopt eqn. (39) with α = 45/2 for frequency redistribution in the line wing, since it is accurate except for galaxies with highly supersolar metallicities, which is unlikely at high redshift.
When the incident frequency is in the Doppler core, the analytic fit, eqn. (39), breaks down, and the emission profile takes on an entirely different form. The emission profile roughly breaks down into two principal components: photons that escape after only a few scatters, and photons that scatter enough times that they reach line center before escape. The former photons retain some "memory" of their incident frequency xi, and produce emission peaks at x = xi and x = −xi. (the peak at −xi is from photons that have undergone a single particle back-scattering off atoms with velocity V = xiV dop along the photon propagation direction). The latter photons lose all memory of their initial frequency, and produce a broader emission peak centered on x = 0. Accordingly, we fit the core redistribution function R core (xi, x) with three gaussians, centered on −xi, 0, and xi, respectively. Let us define the gaussian distribution with r.m.s. frequency σ: By comparing to exact simulations, shown in Figure 10, we find that a decent fit is given by where The above fit works well when dust is included, since the effect of dust on photons in the Doppler core is small in absolute terms (see §3.1.1).
In summary, a simple analytic prescription for the surface frequency redistribution function R(xi, x) is to use equation (39) with α = 22.5 andxi = xi − 2/xi when the incident photon is in the line wing |xi| 3, while equation (43) can be used when the incident photon is in the Doppler core |xi| < 3.

Surface Kinematics
In this section we discuss two kinematic effects of surface scattering. First, we calculate the average scattering angle cosine g = cos θscat , where θscat is the angle between the incident and exiting direction. Second, we calculate the net frequency shift ∆V due to the Doppler shift induced by scattering off a moving surface. In each case, we consider isotropic and perpendicularly incident photons, and use the surface scattering angular distribution Dss(θ), eqn. (28). The case of perpendicularly incident photons is of interest because it applies to photons bouncing around inside a spherical shell; most photons that strike a point on the surface last reflected off the far side of the shell, and hence incident photons have a strong bias to lay along the surface's perpendicular. We use the following conventions: the outward normal to the surfacens defines theẑ direction, an incident photon has directionn1 such thatns ·n1 < 0 with polar angle θ1 and azimuthal angle φ1, an exiting photon has direction n2 such thatns ·n2 > 0 with polar and azimuthal angles θ2 and φ2. For example, in Cartesian co-ordinates,n1 = (cos φ1 sin θ1, sin φ1 sin θ1, cos θ1). The surface has a bulk velocity Vs, with a perpendicular component V ⊥ ≡ns · Vs. For an isotropic angular distribution of incident photons, the polar angle distribution of incident photons is Diso(θ1) = sin θ1, which is normalized to unity over θ1 ∈ [π/2, π], and φ2 is uniformly distributed over 2π. For exiting photons, eqn. (28), the polar angle distribution is given by Dss(θ2) = 2 sin θ2 cos θ2, which is normalized to unity over θ2 ∈ [0, π/2], and φ2 is uniformly distributed over 2π.

The Average Scattering Angle
The scattering angle cosine, also called the scattering asymmetry parameter, is defined by g ≡ cos θscat = n1 ·n2 . For an isotropic incident angle distribution, the average over the angles is Azimuthal symmetry eliminates all the terms inn1 ·n2 that have anx andŷ contribution, leaving just the term cos θ1 cos θ2. The integral over θ1 and θ2 separate, giving The same steps can be carried out for perpendicularly incident photons,n1 = −ẑ, resulting in Surface scattering is characterized by a net average backscatter, g < 0.

Frequency Shift Due to a Bulk Surface Velocity
If the surface has a bulk velocity, then the frequency of the scattered photons suffer a net Doppler shift, due purely to the surface motion. Consider a photon with frequency V striking a moving surface. In the surface rest frame, the photon has incident frequency V ′ = V −n1 · Vs. An exiting photon with frequency V ′′ in the surface rest frame has an exiting frequency V ′′′ = V ′′ +n2 · Vs in the original ("lab") frame. Therefore the surface motion induces a net frequency shift which is in addition to any frequency shift from scattering within the surface. Averaging over an isotropic incident angle gives For perpendicularly incident photons, we simply have n1 = −ẑ. For the exiting direction, averaging over Dss(θ2) gives n2 = 2 3ẑ .
Thus, the average frequency shift for isotropic and perpendicularly incident photons are, respectively, where, as stated above, we define V ⊥ ≡ns · Vs =ẑ · Vs. If the surface is moving away from (towards) the incident photons, V ⊥ < 0, then surface scattering causes a net redshift (blueshift) of the photon frequency.

ANALYTIC MULTI-PHASE TRANSFER
In this section, we build an analytic model for estimating the escape fraction and line width for Lyα escaping from a multi-phase region composed of dusty, optically thick clumps. Compared against Monte-Carlo simulations, these analytic estimates give remarkably accurate results. We consider both stationary clumps and clumps with a Maxwellian bulk velocity distribution, but postpone the discussion of bulk gas outflow/inflow until §5.2. We derive a general photon escape fraction formula in terms of two parameters: the mean number of surface scatters in the absence of absorption N0, and the cloud albedo ǫc, for the case of coherent scattering. The parameter N0 is purely geometric and independent of photon frequency, as long as clouds are very optically thick. We derive fits for N0 for a variety of multiphase geometries of astrophysical interest, and compare the resulting analytic escape fractions to Monte-Carlo simulations of coherent scattering. Applying the coherent scattering formulas to Lyα radiative transfer requires some care, since in this case ǫc is frequency dependent through the frequency dependence of the Lyα scattering cross-section. The cloud albedo, therefore, changes as the photon executes a random walk in frequency. Nonetheless, we find that the analytic formula for coherent scattering can be extended to Lyα scattering as long as ǫc is evaluated at the characteristic escape frequency. In §4.3, we study how the characteristic escape frequency depends upon the frequency redistribution per surface scatter and the random bulk motion of the clumps. We derive estimates of the escaping line profile r.m.s. width and FWHM, and compare these to simulations of repeated surface scatterings.

Escape Fractions
A more detailed prediction for the escape fraction can be given than the simple scaling laws in §2. We derive a generic escape fraction formula in terms of the probability of absorption per scattering, ǫ, and the average number of scatters for escape in the absence of absorption, N0. In the context of multi-phase transfer, each "scatter" refers to an entire surface scattering, in which case ǫ is the cloud albedo: ǫ ≡ ǫc.
Central to the analysis is the average number of interactions for escaping photons, N . Let ǫ be the average absorption probability per interaction, and define D(n) to be the probability distribution for photons to interact n times before escaping, when ǫ = 0 (i.e., in the absence of absorption). For a constant ǫ, the escape fraction can be written The average number of interaction for escaping photons can likewise be written as From these two expressions, we derive This can be inverted to express fe in terms of N : From eqn. (55), N0 ≡ N (0) is given by To derive escape fractions in an arbitrary geometry, let us consider two limits: narrowly beamed forward scattering, and front-back symmetric scattering. The latter refers to the g ≡ cos θsct = 0 case, where there is equal probability to scatter in the forward and backward directions; both isotropic and dipole scattering satisfy this, though dust grains (g d = 0.5) and cloud surfaces (giso = −1/3, from eqn. (46)), do not. As we shall show, however, the escape fraction for scattering for any value of g < 1 can be based upon the g = 0 case as long as there are sufficient scatters that the photon's trajectory is randomized.
The case of purely forward scattering is trivial: since the photon does not change direction, we have where we have used eqn. (57) in the final step. Note that the escape fraction is same as if there were no scatterers, since scattering does not alter the photon trajectory. For frontback symmetric scattering (g = 0), the escape fraction is (see Figure C1 in Appendix C): where we applied eqn. (57) to obtain N0 = 1 2 (τ 2 + 2τ ) in the second step 7 . Although this formula for the escape fraction is derived assuming g = 0, we show in Appendix C that this formula is valid for any g as long as there are enough scatters to randomize the photon's direction. Specifically, if N0 > n * (g) (where n * (g) is given by eqn. (C3)), then eqn. (59) can still be used. Otherwise, eqn. (58) is a better approximation to the escape fraction-since a situation with such few scatters and/or a strongly peaked forward scattering profile converges to the straight-line trajectory case. Note that in both these limits, the escape fraction only depends upon a single parameter, ǫN0. The average number of interactions required for a photon to be absorbed is Na = 1/ǫ, so the controlling parameter is equivalent to ǫN0 = N0/Na.
The average number of interactions for an escaping photon, N , can be derived by applying eqn. (55) to the appropriate escape fraction formulae, eqn. (58) or eqn. (59). For straight-line trajectories, the result is while for the random walk trajectories 7 Within the Eddington approximation, the escape fraction from a source in the mid-plane of a slab is fe = 1/ cosh √ 3 ǫ 1/2 τ , as derived, for example, in Neufeld (1991). Instead, we propose eqn. (59), which from our simulations is more accurate for 1-D scattering (see Appendix C).
In the limit √ ǫN0 ≪ 1, we find N → N0 − ǫN0( 2 3 N0 + 1), and in the limit √ ǫN0 ≫ 1, we have N → (1−ǫ) 2ǫ √ 2ǫN0. The entire effect of the cloud geometry is characterized by a single number, N0. Directly calculating N0 from a given geometry is not practical except for the simplest geometries. In general, given a clump geometry, N0 must be computed via a simulation, where one can use the exiting angle distribution, eqn. (28), for each surface reflection. In practice, we find that for many generic geometries expected to crop up in astrophysical applications, an appropriate "line of sight"averaged N0 can be accurately determined from simple geometric parameters (such as the cloud covering factor fC ). We now proceed to do so.

Example Geometries
We now test the accuracy of the analytic formula eqn. (59), and investigate its application in a variety of multi-phase geometries. For each geometry, we fit N0 as a function of the appropriate natural geometric parameter (such as cloud covering factor fC ) using Monte Carlo simulations. We then calculate the escape fraction using N0 and eqn. (59), and compare this to simulations. In several cases, we find that we obtain better fits by rescaling with an order unity fitting parameter κ, where we use κǫcN0 in place of ǫcN0 in the escape fraction formulae eqn. (59). In general, even with no correction factor, eqn. (59) is accurate to < ∼ 10% when fe 10%, and is generally correct to within a factor < ∼ 2 for fe < 10%. When the escape fraction is very small, the photons that do escape comprise the rare trajectories. As their transfer behavior can depend sensitively on the specifics of the geometry, it is not surprising that eqn. (59) breaks down when fe is very small. In any case, Lyα emission is undetectable for such small fe, so these cases are of little observational consequence.
Some notes about our Monte-Carlo simulations: for simplicity, we assume that the region in between the scattering surfaces is empty (we justify this assumption in §4.5). We also assume that scattering surfaces are extremely optically thick, so that the surface scattering approximations of §3 apply. In particular, when a photon hits a gas surface, it has a net probability ǫc of being absorbed; if it survives, then its exiting angle relative to the surface normal follows the distribution Dss(θ) = sin 2θ, and the exit location is the same as the point of incidence.

Spherical Clumps
The canonical example of a multi-phase geometry is a spherical region populated by randomly placed, optically thick spherical clumps. Such clumps tend to be cool-phase gas such as molecular clouds, which arose via thermal instability. The natural geometric parameter is the mean number of clouds intersected along a random line of sight, called the cloud covering factor fC . The covering factor is analogous to the interaction optical depth τ for homogenous media. For a central source, fC is measured from the region center to the edge. We computed N0 for various covering factors and clump radii distributions. As shown by the solid line in the top panel of Figure 11, we find that a general fit for any clump radii distribution is This highlights the fact that when surface scattering applies, the spherical clump geometry does not depend upon the size distribution of the clumps nor their volume filling fraction; the entire geometry is characterized by a single parameter, fC . This was postulated by Neufeld (1991); we have confirmed this insight in our simulations. The scalings in eqn.
While the covering fraction is an unknown free parameter, it is easy to see why fC ∼ 1 is reasonable. Suppose the cold clouds constitute a mass fraction fM of the galaxy, and are overdense by a factor δ relative to the intercloud medium. Assuming pressure balance, δ ∼ TICM/Tc ∼ 100, for Tc ∼ 10 4 K and TICM ∼ 10 6 K for the cloud and intercloud medium temperatures respectively. The volume filling factor of clouds is then fV = fM /δ. The number density of clouds is nc ∼ fV /Vc, where Vc ∼ r 3 c is the volume of a typical cloud. The cloud covering factor is:

Random Surfaces
An abstraction of the spherical clump geometry is to have the photons strike a surface after traveling a distance ℓ, where ℓ is drawn from a given probability distribution, and where each surface is randomly oriented with respect to the photon direction. We have investigated the exponential distribution exp(−fCℓ/R), where R is the radius of the region and fC is the covering factor. The photon escapes once it leaves the spherical region of radius R. When a photon traveling in the directionnp strikes a surface, the outward normal of the surfacens is drawn from an isotropic distribution such thatnp ·ns < 0. As shown by the top panel in Figure  12, a good fit for the scattering number is The random surfaces model is faster to simulate, and any results reliably apply to the spherical clump model when expressed in terms of N0. The bottom panel of Figure 12 compares the escape fraction formula eqn. (59) to simulations. This procedure for generating random surfaces can in fact be used to quickly simulate any arbitrary geometry, with a suitable characterization of the probability distribution of path lengths ℓ and surface orientationsns, though, of course, the fitting formulae for N0 will depend on these quantities.

Shell with Holes
Another geometry that frequently arises in astrophysics is that of a shell of material surrounding a photon sourcefor instance, in stellar and galactic outflows. Much work has been done on Lyα radiative transfer through opaque shells (Tenorio-Tagle et al 1999;Ahn, Lee, & Lee 2003;Ahn 2004), but the effect of gaps in the shell has not been investigated. Since a completely homogeneous shell of gas is rarely, if ever, observed, a shell with holes is an interesting geometry. The natural geometric parameter here is the fraction of the solid angle covered by the shell, fcov (i.e., the gaps comprise a total solid angle 4π(1 − fcov)). To estimate N0, assume that during each bounce the photon has a probability of (1 − fcov) to escape through a gap. This leads to the expression which is shown by the solid line in the top panel of Figure  13. We computed N0 for shells with a random placement of non-overlapping circular gaps, and found that many small gaps were equivalent to fewer large gaps with the same fcov, as to be expected. As shown in the bottom panel of Figure  13, the analytic form for fe does well where expected, and only breaks down when fcov is near unity and ǫc > ∼ 0.3. This agreement is quite remarkable, given that the escape fraction formula (59) was derived assuming a very different scattering particle geometry (homogeneous media). This gives us confidence that the impact of geometry on the escape fraction can indeed be encapsulated by a single parameter ǫcN0.

Open Ended Tube
To illustrate how well photons can escape through cracks and gaps in optically thick material, we consider photons escaping from the middle of an open-ended tube. This geometry may apply to star forming regions or AGNs where outflows have punched (bipolar) cavities through a surrounding gas cloud, allowing the photons to escape through the cavities (e.g. Shopbell & Bland-Hawthorn (1998)), although we do not investigate the dependence on the opening angle. The natural geometric parameter is the ratio L/R, where L is the total length of the tube and R is the tube's radius. The average length traveled along the tube per scatter is ℓ ∼ R. Therefore N0 is estimated from L/2 ≈ √ N0ℓ, which implies N0 ∼ (L/R) 2 . As shown by the solid line in the top panel of Figure 14, a decent fit for the number of surface scatters is N0 = 1 10 (L/R) 2 + 13 20 (L/R) As shown in the bottom panel of Figure 14, the escape fractions fit the data very well, except as ǫc → 1. As previously discussed, for such low escape fractions, rare escaping photons follow unusual trajectories not well captured by our formalism. In any case, for such low escape fractions, Lyα emission cannot be observed, and the results have no observational relevance.

Line Widths
In this section we shall consider two separate effects on the Lyα frequency as it escapes: the frequency redistribution due to thermal motion of atoms, as well as random bulk motion of the scattering surfaces. The effect of global outflow/inflow on the line profile is discussed separately, in §5.2. We consider two simple diagnostics of the profile: the r.m.s. frequency shift Σ and the FWHM Γ. We ran simulations of dust-free surface scattering to determine Σ and Γ as a function of the number of surface scatterings, nss. Although escaping photons vary in the number of scattering surfaces they encounter before escape, in practice we find that line profiles can be accurately characterized by the average num- ber of scattering surfaces encountered. We now derive formulas for Σ(nss) and Γ(nss) for resonant scattering frequency redistribution and random bulk surface velocities.

Resonant Scattering
By using the analytic approximation to the surface frequency redistribution function R(xi, x) derived in §3.3.5, we simulated repeated surface scattering for an ensemble of photons. Carrying this out using the exact Monte-Carlo simulation is not computationally feasible. For a few cases we compared the results of the analytic approximation with an exact Monte-Carlo simulation and a Monte-Carlo simulation that approximates the core scatterings ( §3.1.1), and all three are in good accord. This gives us some confidence that the inaccuracies in the analytic approximation do not compound significantly for the number of scattering surfaces we investigated. As shown in Figure 15, decent fits for the resonant scattering Σ and Γ are Γrs(nss) = 3 nss 1 + nss V dop .
The behavior can be easily understood. If atomic thermal motions are responsible for frequency redistribution, then the line profile quickly relaxes to a Gaussian whose FWHM is determined by the width of the Doppler core. However, the profile also acquires non-Gaussian tails from rare scattering events; these tails dominate the r.m.s. line width Σ, hence Σ ∝ nss. For resonant scattering only, the FWHM Γ is a more accurate measure of the line profile than Σ. Examples of the spectra after repeated surface scatters is shown in Figure 16.

Surface Motion
As shown in §3.4, when photons scatter off a moving surface there is a net frequency shift per surface scatter of ∆V ∼ V ⊥ , where V ⊥ is the velocity along the outward normal. If the moving surfaces have a Maxwellian velocity distribution with r.m.s. speed V c , then we expect that the induced r.m.s. frequency shift after nss scattering surfaces will scale as Σsm(nss) ∼ √ nss V c . For a Gaussian distribution, the FWHM and the standard deviation are related by Γ ≈ 2.2Σ, and so we expect Γsm(nss) ≈ 2.2 √ nss V c . To check this, we simulated repeated surface scattering assuming an isotropic incident angle distribution and the surface scattering exiting angle distribution, eqn. (28). We indeed find that which is shown in the middle panel of Figure 15.

Combined Effect
In the two bottom panels of Figure 15, we show the combined effects of resonant line broadening and surface motion. For most multi-phase geometries, V c ≫ V dop , so line broadening is dominated by cloud motions. In this regime, the line width can be accurately estimated by eqn. (68). Only if cloud motions are small and comparable to atomic thermal motions, V c < ∼ V dop , does the behavior change. In this case, the FWHM does not increase with the number of scatters; instead it "thermalizes" to the characteristic Doppler width. The line profile is more accurately described by eqn. (67). When V c > ∼ V dop , the total r.m.s. width is (note that the dispersions add linearly rather than in quadrature) Σ(nss) = 2Σrs(nss) + Σsm(nss) .
The FWHM has a more complex behavior because the Doppler core tends to retard the increase in Γ beyond the size of the Doppler core.

Lyα Escape Fractions
In this section we show how Lyα multiphase transfer can be handled analytically based on the results in the previous subsections. We shall first estimate the typical escape frequency of Lyα photons xe, which provides the typical surface absorption albedo ǫc(xe). With this estimate of ǫc and N0, we can calculate the typical Lyα escape fraction fe using eqn. (59). We find that an adequate approximation for the typical escape frequency is given by a slight modification of the r.m.s. line width Σ(nss), eqn. (69), where the surface motion can be either a random Maxwellian cloud velocity with r.m.s. speed V c or a bulk outflow with typical outflow speed V c (see §5.2 below). It is not correct to approximate the typical number of scattering surfaces nss with the average number of scattering surfaces in the absence of absorption, N0, since in general nss ≪ N0 once absorption is taken into account. A more appropriate measure of nss is given by N , eqn. (61), the average number of scattering surfaces encountered before escape for a fixed cloud albedo ǫc. However, since ǫc is frequency-dependent, we still need to estimate the escape frequency xe. We do so iteratively. We first estimate xe assuming that absorption is unimportant, which obviously overestimates xe. We can then estimate N as: where with Σ(nss) given by eqn. (69).
With the above estimate of nss = N , a better approximation to xe is given by the next iteration, which gives a better approximation for the typical cloud albedo At this point one could iterate again-using ǫ ′ c to find a better approximation to nss-but we find that stopping after one iteration provides escape fractions in good accord with simulations. From eqn. (59), the escape fraction given by In Figure 17 we compare this analytic escape fraction to simulations of radiative transfer through the random surfaces geometry (see §4.2.2) for both random cloud motions and a bulk cloud outflow, and for several amounts of dust. The bulk cloud outflow is purely radial, with the same speed at all radii, which approximates galactic winds outside the initial acceleration zone. As can be seen, the analytic approximation captures the simulated escaped fraction to ∼ 20% when cloud motion dominates atomic thermal motion, V c > ∼ 3V dop ∼ 40 km/s. As explained in §4.3, once the effects of the Doppler core become important, the r.m.s. frequency shift Σ is no longer a good measure of the 2) with f C = 3 were run for two types of bulk cloud motion: a Maxwellian velocity distribution (open symbols), and a purely radial outflow (filled symbols). In the former case the x-axis corresponds to the r.m.s. speed of the clouds, while in the latter case the x-axis corresponds to the outflow speed, which is the same at all radii. Two different dust contents were simulated, both with ǫ d = 0.5: σ d −21 = 2 (circles) and σ d −21 = 0.2 (squares). In all cases the gas temperature is 10 4 K. The analytic approximations (lines) are generated by the steps outlined in text. typical frequency: Σ will overestimate the typical escape frequency, and hence lead to an overestimate of the absorption. This is seen in the figure for V c < ∼ 40 km/s. However, since the amount of absorption is typically not significant when the escape frequency is < ∼ 40 km/s, for most purposes using Σ allows one to accurately estimate fe. Note that for bulk outflows with the same characteristic speed V c , the escape fraction is smaller than for random cloud motion. The line profile for random cloud motion is approximately Gaussian centered on the line center, with a standard deviation of V c . In contrast, the line profile for a bulk outflow is has a mean at V c . For a given V c , a bulk outflow produces more photons far from line center than random cloud motion, and hence the bulk outflow causes more Lyα absorption.
In summary, the escape fraction depends upon five parameters: the typical cloud speed V c , the number of surface scatters in the absence of absorption N0, the gas temperature T , and the dust parameters σ d and ǫ d .

Dust and Gas Between the Clumps
Thus far, we have only considered radiative transfer off opaque surfaces, and ignored absorption and/or scattering in the optically thin hot intercloud medium (ICM). We treat the hot intercloud medium as having a very low neutral HI fraction (for gas in coronal equilibrium at T ∼ 10 6 K, xHI ∼ 10 −5 ), as well as being dust-depleted, due to sputtering and other dust-destruction processes. ICM resonant scat-tering and absorption is easily incorporated in our Monte-Carlo simulations at little computational cost, and we have experimented with various prescriptions for the ICM in our runs. In most cases, we find that radiative transfer within the ICM can be neglected, and here we show some simple estimates demonstrating why this is the case.
Let us first consider dust absorption in the ICM, assuming that resonant scattering in the ICM is negligible. It is easy to show that for a multiphase medium, the relative fractional column densities of a species i in phases X, Y is simply given by the relative mass fraction of species i in that phase: where x i X is the mass fraction of species i in phase X. Thus, for instance, N c HI /N ICM HI ≈ fM,c[x ICM HI ] −1 ≫ 1 for x ICM HI ≪ 1, and the observed HI column density N obs HI,21 ∼ 1 is strongly dominated by gas in the cold phase.
We can use this to estimate scattering in the ICM. Due to its random walk as it scatters off optically thick surfaces, a Lyα photon traverses an ICM column density ∼ √ N times larger than for a straight line path. At a characteristic escape frequency V e 2 ∼ 1, it therefore encounters an ICM HI optical depth: Hence, resonant scattering in the ICM is negligible. It only becomes important when the photon is within the Doppler core V2 ∼ 0.26 for the parameters above). However, even in this rare case τ ICM HI ∼few, comparable to the number of surface scatterings N .
What about dust absorption in the ICM? Let us suppose that due to dust depletion, only f ICM d ∼ 0.05 of all the dust in the galaxy is in the ICM. If the total dust optical depth through the cold phase is τ c d ∼ 1 N tot HI,21 σ a −21 , then the total dust optical depth through the ICM is τ ICM The mean free path of a Lyα photon is l ∼ r gal / √ N . Hence, between each bounce, a Lyα photon traverses an optical depth τ bounce N , and has a probability of absorption in the ICM of P ICM absorb ≈ 1 − e −τ bounce a ≈ τ bounce a , or: By contrast, during each surface scatter, the photon has an absorption probability of: Thus, photons are more likely to be absorbed on cloud surfaces, rather than in the ICM, justifying our neglect of ICM dust absorption. Obviously, all of these statements are parameter dependent and there are cases when scattering and absorption in the ICM cannot be neglected (if the IGM dust or HI content is high). For instance, if the HI fraction in the ICM is high, then Lyα can be strongly quenched. This may be partly responsible for the variation in the Lyα equivalent widths amongst different galaxies (see discussion in §5.1).

Accelerated Radiative Transfer on a Grid
The approach of this paper is to identify all optically thick surfaces in a Monte-Carlo simulation and apply the scattering properties identified in §3 to them, thus affording vast computational speed-ups. This should be readily amenable to radiative transfer in numerical simulations. For this scheme to be accurate: (i) the surfaces must be sufficiently optically thick that transmission is negligible; i.e., they should satisfy equation (32). (ii) The approximation that the photon is either reflected or absorbed on the spot without significant spatial diffusion must hold. We now discuss this second requirement.
For the "on the spot" approximation to hold, the photon's mean free path ℓi should be significantly smaller than the grid cell size L cell . The photon typically moves a distance ∼ √ N surf ℓi ∼ 5ℓi (see §3.3.4 for discussion of N surf ) whilst scattering within the surface. Thus, we require ℓi < ∼ α L cell where α ∼ 1/10 is a constant that designates the desired level of accuracy for the approximation. The "on the spot approximation" is accurate if the HI density is larger than a critical density: n HI cell n * = 13 (α/0.1) where V i is the incident Lyα frequency shift off line center in the rest frame of the HI in the cell (in velocity units).
Alternatively, if one is willing to sacrifice some spatial resolution, one can consider a group of neighboring cells which are opaque across their total length (and thus satisfy equation (32)) even though an individual cell may not be opaque. For cubic blocks of N blk cells per side, the entire block surface will act like an absorbing mirror if n HI blk n * = 13 (α/0.1) where n blk is the average HI density within the block. The surface approximations break down if the block is strongly inhomogeneous, i.e., the mean free path changes over a length scale that is much shorter than the block length. This is equivalent to |n cell /∇(n cell )| βN blk L cell , where β ∼ 1/10 is a constant that designates the desired level of accuracy for the approximation. In terms of cells on a grid, let n(2) cell and n(1) cell be any two neighboring cells in the block. The second condition on the absorbing mirror approximation takes the form We look forward to implementing these ideas in numerical simulations in the future.

APPLICATIONS
We briefly discuss two examples of applications of our radiative transfer framework. Many more extensive studies are possible.

Lyα Equivalent Widths
Most Lyα photons are produced in the HII regions surrounding sources of ionizing radiation, where roughly 2/3 of the ionizing photons are converted into Lyα photons (under case B recombination). Hence, in the absence of radiative transfer effects, the equivalent width measures the number of ionizing photons emitted relative to the UV continuum near 1200Å. There are numerous examples of high-redshift sources which have equivalent widths which are too large to be produced by conventional stellar populations. About ∼ 2/3 of the SCUBA submm galaxies with accurate positions from radio detections have Lyα in emission, many with equivalent widths too great for stellar sources (Smail et al. 2004). The mysterious Lyα emitters at z ∼ 3.1 observed by Steidel et al (2000) have enormous Lyα fluxes, but no observed continuum. Finally, the LALA survey detects many high redshift z = 4.5, 5.7 sources with equivalent widths EW 150Å significantly in excess of any known nearby stellar population (Rhoads et al. 2003). An AGN origin is unlikely because follow-up observations show no signs of the X-rays and high-ionization lines expected for a type II quasar source (Wang et al. 2004;Dawson et al. 2004). Another possibility is that the Lyα emission is due to an extremely top-heavy population of massive PopIII stars. However, there are no signs of the strong HeII emission at 1640Å expected from metal-free stars (Dawson et al. 2004).
Another possibility for high Lyα EWs, originally suggested by Neufeld (1991), is radiative transfer effects. If the continuum is more absorbed than Lyα photons during the escape from the host galaxy, then the equivalent width of the transmitted spectra is larger than the equivalent width of the source. The initial and transmitted equivalent widths are basically related by the ratio of Lyα to continuum escape fractions, where EWsrc is the source equivalent width and EWout is the equivalent width for the escaping photons. In order for a "normal" starburst IMF with an intrinsic equivalent width of EWsrc ∼ 150Å to produce an observed equivalent width of EWout > ∼ 300Å, then radiative transfer must account for a "boost" by a factor of at least 2-3. For sources where no continuum is observed, the continuum must be preferentially extinguished by an even larger factor.
Let us now estimate equivalent width boosts in our multi-phase model to see if this is possible. For any multiphase medium where the gas resides in clumps that are very opaque to Lyα , the surface scattering approximations apply, and so the Lyα escape fraction can be analytically estimated as in §4.4. What about the continuum escape fraction? For simplicity, assume that each gas clump is not opaque to dust extinction: τ d c < ∼ 1, where τ d c is the dust extinction (scatter-ing+absorption) optical depth across a clump diameter 8 . Since the self-shielding effect of clumpy gas is therefore small for the continuum photons, the effective dust distribution is approximately homogenous for the continuum radiative transfer. The escape fraction for a photons injected in the middle of a homogenous medium, with an absorption albedo ǫ d ≈ 1/2, is approximately that given by eqn. (59), where τ d ≡ N σ d is the average dust extinction optical depth through a region with average HI column density N . Figure  18 shows how the ratio f Lyα e /f ctm e varies as a function of N0 for a fiducial set of multiphase gas parameters. A substantial "boost" in the transmitted equivalent width due to selective absorption of the continuum is quite reasonable, as long as two basic conditions are met: 1) there must be enough dust present to absorb a substantial fraction of the continuum. 2) this dust must be pre-dominantly located in dense neutral gas, so that the Lyα photons are shielded from absorption. The latter condition is discussed in §4.5 above, and so we turn to discussing the first condition.
The Lyα escape fraction depends only weakly on the overall dust content of the galaxy. In Fig 18, f Lyα e ∼ 0.2−0.9 over a wide range of parameters, with the escape fraction decreasing with the number of surface scatters N0 and bulk gas motion Vc. On the other hand, because continuum photons are not shielded from dust by resonant scattering, they see the full optical depth of dust absorption, and very approximately, f ctm e ≈ exp(−τa). A significant boost in the equivalent width therefore requires that τ d 1. If the average HI column density across the region is N , then τ d ∼ 1 requires σ d −21 ∼ 1/ N21 . A damped Lyα type system with an average column density N ∼ 10 22 cm −2 would require a dust extinction cross-section per hydrogen of σ d ∼ 10 −22 cm 2 /H, which roughly corresponds to a metalicity of ∼ 1/10 solar. For sources in which these differential radiative transfer effects are taking place, the equivalent width should statistically have a positive correlation with the FIR flux. This correlation could potentially break down in more developed galaxies at lower redshifts, where the Lyα shielding effect of the HI can be broken by higher gas speeds in deeper potential wells (rendering clumps optically thinner in Lyα ), as well as the build-up/survival of dust in low density interclump gas.

Multi-phase Outflows
We turn to discussing the effect of a multi-phase gas outflow (or inflow) on the Lyα emission line profile. The effects of an outflowing shell (e.g. Tenorio-Tagle et al (1999); Ahn, Lee, & Lee (2003); Ahn (2004)) and a Hubble like expansion of a uniform gas sphere (e.g. Loeb & Rybicki (1999); Zheng & Miralda-Escudé (2002)) on the Lyα emission line is a well studied problem. In both the expanding shell and expanding sphere scenarios, the generic effect is that a characteristic outflow speed V f produces a redshifted emission peak at V peak ∼ −V f , with an asymmetric shape that has a longer tail on the red side of the peak. The peak comprises photons that reflect off the far side of the expanding gas, which Doppler shifts the frequency by ∼ −V f (see §3.4). However, in order for these singly backscattered photons to escape, the intervening gas column density must be small enough for the photons to be transmitted, rather than be reflected a second time. In the rest frame of the near-side shell, the singly back-scattered photons have a frequency shift V ′peak ∼ −2V f . In §3.3.1, we showed that a non-negligible amount of Lyα photons will be transmitted through a slab if N21 < N pt 21 ≈ 0.05[V2] 2 [σ a −21 ] −1 . Setting V ∼ −2V f , we see that an outflow with speed of 200 km/s will only allow a non-negligible amount of singly back-scattered photons to be transmitted if the intervening column density is N 2 × 10 20 cm −2 . Observational estimates of the column density in galactic winds often exceed this, yet Lyα is still often seen.
The main distinguishing feature of a multiphase outflow is that it allows photons of any frequency to escape even when the intervening gas column depth is very large, N N pt . As in the homogeneous gas outflow models with  Figure 19. Outflowing Shell with Gaps The normalized emission profile as a function of the frequency shift off line center, for an outflowing shell with speed V f = 200 km/s. Two dust amounts were simulated: σ a = 0 (thin line, filled in) and σ a −21 = 1 (thick line). In all simulations, we used a gas temperature T 4 = 1. For high redshift galaxies, the blue side of the profile would be quenched by IGM absorption. The spike of photons that escape at line center is easily scattered out of the line of sight by a small (N 5 × 10 13 cm −2 ) intervening column density of HI, and thus not likely to be observed. smaller column densities N < N pt , we find that for multiphase outflows, the emission peak is redshifted by ∼ V f . However, emission is still detectible even when N ≫ N pt , as expected.
In particular, we investigated the emission profile for two basic types of multiphase outflow geometries: an outflowing shell with holes ( §4.2.3) and an outflowing ensemble of gas clumps modeled with the Random Surfaces geometry ( §4.2.2). In both cases, all surfaces were given a radially velocity with constant speed. This choice is meant to reflect galactic winds, where the gas reaches the asymptotic wind speed quickly. We placed a source of line center photons in the center of the region. Since the regime of optically thick gas has been given the least attention, we assume that the extreme case holds, where none of the photons penetrate through the gas. In this limit the surface scattering approximations of §3.3 apply in the rest frame of the scattering surface. The kinematics of Doppler shifts in and out of a moving surface can be found in §3.4. In order to distinguish the effects of outflow from the effects of random bulk gas motion, we assume that there is no random bulk motion, so that each gas surface has an exactly radial velocity, V s = V fr .

Outflowing Shell with Holes
For photons perpendicularly incident on the inner side of the shell, the net frequency shift induced by a moving scattering surface is ∆V = − 5 3 V f (see §3.4.2). Since most photons escape after encountering either zero or one scattering surface, the two dominant peaks in the emission profile will be at V (0) = 0 and V (1) ≈ − 5 3 V f . The maximum frequency shift after encountering one scattering surface is ∆Vmax = −2V f , and so the emission peak at V (1) should have a sharp cut-off at V = −2V f . These emission peak features are verified in Figure 19. The series of secondary peaks are at integer multiples of V (1) , but tend to become blended together to form an long red-side tail to the profile. A possible exception is the third peak at V (2) = 2V (1) = − 10 3 V f (composed mostly of photons that encounter exactly two scattering surfaces before escaping), which can also be prominent in the profile.

Outflowing Clumps
For a simple model of outflowing gas clumps we used the Random Surfaces geometry, which is described in detail in §4.2.2. In Figure 20, we show how the profile varies with the covering factor, for dust free gas σ a = 0 and for a Milky Way type dust content σ a −21 = 1. As for the case of an outflowing shell with holes, the inclusion of dust suppresses photons which have redshifted far from line center (so the HI no longer shields them from the dust), and sharpens the line profile. To provide some insight into how a clumpy outflow causes a redshift in the emission line, in Figure 21 we show how the photon radius, frequency, and cloud absorption albedo vary as a function of nss.

CONCLUSIONS
Our main technical, radiative transfer results are: • With the aid of Monte-Carlo simulations, we study the scattering properties of Lyα photons incident on an opaque, dusty, moving cloud. We derive fitting formulas for the absorption probability, frequency and angular redistribution functions of incident photons.
• These formulas can be incorporated into radiative transfer codes, affording a vast computational speed up, and making feasible otherwise intractable calculations.
• Analytically, a multiphase gas geometry can be accurately characterized by a single number, N0, the number of surface scatters in the absence of absorption. Other factorssuch as the cloud radii distribution for fixed N0-are generally unimportant.
• We derive analytic formulas for the Lyα escape fraction and line widths.
• Several archetypal geometries are explored: randomly placed spherical clouds, randomly placed surfaces (an abstraction of the prior geometry), a shell with holes, and an open-ended, cylindrical cavity.
• Constant speed, radial, outflows are analyzed for broken shells and random surfaces. The red-shifted peaks and widths are connected to the geometry and outflow speed.  Figure 19. From top to bottom, each panel is a different covering factor: f C = 1, 3, 5. The escape fractions for the σ a −21 = 1 simulations are indicated. The delta function emission spike at V e = 0 is composed of all the photons that escape freely without striking a clump. As noted in Figure 19, exact line center photons are likely to scattered out of the line of sight before being observed.
Our main results of direct observational relevance are: • Lyα can escape from multi-phase dusty galaxies for HI column densities where it would be strongly quenched in a single-phase medium.
• If most of the dust resides in a neutral phase which is optically thick to Lyα, the Lyα equivalent width can be strongly enhanced: while Lyα photons typically scatter off such surfaces (which shield the dust), continuum photons penetrate inwards and are preferentially absorbed.
• When the characteristic bulk gas speed exceeds ∼ 100 km/s, the Lyα line width is dominated by the gas motion, and resonant scattering frequency redistribution is subdominant effect.
• Multiphase outflows generically produce Lyα line profiles that have the characteristic asymmetric shape seen in many starburst galaxies and Lyα emitters.
• Multiphase outflows can produce line widths several times larger than the actual outflow speed.
The ISM of galaxies at both low and high redshift is almost certain to be both dusty and multi-phase: metal and dust production begin very early, given the short lifetime of massive stars, and thermal instability is almost inevitable under galactic conditions. Nonetheless, despite an extensive literature, to the best of our knowledge this is the first de- tailed numerical study of resonance-line radiative transfer in a multi-phase dusty medium. The ground is surprisingly rich, and many future applications are envisaged!

APPENDIX A: TESTS OF THE CODE
In this Appendix we show various tests of our Monte Carlo code. First, we test the code against known analytic solutions for optically thick slabs. Second, we compare the acceleration scheme described in §3.1.1 to exact simulations.

A1 Comparison to Analytic Solutions
We tested the exact Monte Carlo code against known analytic solutions for very opaque slabs with a source at the mid-plane. First, we compared the emission frequency profile when the slab is pure HI, so there is no absorption. Second, we compared the escape fractions when the slab contains a small amount of absorbing dust, where the absorption crosssection of the dust is frequency independent. For an optically thick (aτ h0 > 10 3 ), uniform slab without dust, where τ h0 is the center-to-surface hydrogen optical depth at line center 9 , Harrington (1973) derived the mean 9 In many papers on this subject, e.g. Harrington (1973) and Neufeld (1990), "τ 0 " refers to the mean optical depth, while we use the line center optical depth. In our notation, φ(x)τ 0 = Φ(x)τ h0 where φ(x) is the normalized Voigt profile, which means τ 0 = √ πτ h0 when a ≪ 1. The surface mean intensity J(x) as a function of frequency, for a central source of Lyα photons in an optically thick slab. For three different values of aτ h0 (labeled in the plot), the analytic surface intensity (grey line) is compared to exact Monte-Carlo simulations (circles). The optical depth is fixed at τ h0 = 10 6 . Three different gas temperatures are used: T = 10 K, 10 2 K, 10 4 K, which corresponds to aτ h0 = 10 4.17 , 10 3.67 , 10 2.67 , respectively. The relative error for each frequency bin is approximately 1/ √ N bin , where N bin is the number of photons in the bin.
intensity emission spectrum, for line center photons which are injected at the slab center (see Neufeld (1990) for various extensions) J(±τ h0 , x) = 6 1/2 24π 1/2 x 2 aτ h0 cosh π 3/2 54 1/2 |x| 3 aτ h0 where J is the mean intensity and τ0 is the line center hydrogen optical depth from center-to-surface. Figure A1 compares the simulation to the formula for slabs at three different temperatures.
For an optically thick slab (aτ h0 > 10 3 ) with a centerto-surface absorption optical depth τa, Neufeld (1990) derived the exact escape fraction for a central source, in the limit (aτ h0 ) 1/3 ≫ τa, as well as several approximation formulae. In particular, the escape fraction is well approximated by fe ≈ 1/ cosh 3 1/2 π 5/12 ζ (aτ h0 ) 1/3 τa , where ζ is an order unity fitting parameter. Neufeld (1990) found that the choice ζ = 0.525 gives a good approximation to the exact analytic escape fraction. For comparison, we have included in Figure A2 the escape fractions found by Hummer & Kunasz (1980) using numerical integration techniques and the escape fractions found by Ahn, Lee, & Lee (2000) using a Monte Carlo simulation similar to ours. We find that the choice of fitting parameter ζ = 0.5 gives a good approximation to our Monte Carlo results. Both Ahn, Lee, & Lee (2000) and our Monte Carlo simulations show slightly more absorption than the analytic formula from Neufeld (1990). The analytic treatments assume the Lorentzian wing profile all the way down to x = 0, neglecting the Gaussian core. This will underestimate the number Us. MC, with dust scattering. aτ 0 =10 3.67 Figure A2. Dusty Slab Escape Fraction, Central Source. The escape fraction from the middle of an optically thick, dusty slab is compared for several analytic and numerical methods. As has been shown analytically, the escape fraction depends mainly upon the combination (aτ h0 ) 1/3 τa. Including dust scattering, using ǫ d = 0.5 and g d = 0.6, did not significantly affect the escape fraction.
of scatters spent in the core. Although the absorption probability per interaction is small in the core, neglecting core bounces will cause a slight underestimate of the overall absorption probability.

A2 Testing the Acceleration Scheme
We tested the acceleration scheme, described in §3.1.1, by comparing the surface absorption probability against exact simulations. Exact simulations are computationally expensive, and so we could only test the acceleration scheme against a handful of exact cases. As shown in Figure A3, the acceleration is quite accurate, even for initial frequencies in the line core.

APPENDIX B: SURFACE Lyα FREQUENCY REDISTRIBUTION FORMULA
In this appendix we first show that R(xi, x; α) given in eqn.
(39) has a unit norm over x, as claimed. Second, we outline the steps used to derive eqn. (41), the generating function for R(xi, x; α).
To integrate R(x, x; α) over x ∈ (−∞, ∞), first change variables to u ≡ x 3 −x 3 i . Then the integral becomes The integral over the functional form (A + u 2 ) −1 is a standard integral: Use of this formula in eqn. (B1) shows that the integral over x equals one, and so R(xi, x; α) is correctly normalized over x as advertised. To randomly generate exiting frequencies x that obey the probability distribution R(xi, x; α), we use the transformation method (Press et al 1992). First, select a random univariate u ∈ [0, 1] and set The frequency x is then given by functional inversion x = F −1 (u). As above, the integral is best carried out by changing variables to u ≡ x 3 −x 3 i . Then F (x) is given by eqn. (B1), except that the upper limit is "x 3 −x 3 i " rather than "∞". By using equation (B2) to complete the integration, we find that The functional inversion gives the randomly drawn exiting frequency x: x = x 3 i −x 2 i tan (πu) 1/3 .
APPENDIX C: 1-D TRANSFER FOR AN ARBITRARY SCATTERING ASYMMETRY PARAMETER The escape fraction for arbitrary g can be approximated by the g = 0 formula, eqn. (59), as we demonstrate next. Define n * to be the average number of interactions required for a 50% chance of back-scattering. For interactions with scattering parameter g, the probability of a forward scatter is (1 + g)/2 and that of a back-scatter is (1 − g)/2. Therefore n * is defined to satisfy which leads to or equivalently Every n * interactions acts like a single g = 0 interaction. The probability of absorption after n * interactions is 1−(1−ǫ) n * . Consequently, the escape fraction is approximately the same as the g = 0 formula, eqn. (59), with a re-scaled N0 and absorption albedo, N * 0 and ǫ * , given by The approximate escape fraction for arbitrary g is where N0 = N0(g) is calculated for the given value of g and n * is a function of g through eqn. (C3). Note that for 1-D transfer when g = 0, then N0(g = 0) = 1 2 (τ 2 + 2τ ), which we have verified with simulations. If n * < N0 then there is enough scatterings for this approximation to hold. If ǫn * ≪ 1 also holds, then ǫ * ≈ n * ǫ and the rescaling leaves the product ǫN0 unchanged. In this limit, the escape fraction given in eqn. (59) is valid for any type of scattering, and represents, therefore, the generic escape fraction form for any type of "random walk" photon transfer. On the other hand, if n * N0 then this approximation breaks down, and the trajectory of the photon is more accurately characterized by straight line motion, with negligible back-scattering, eqn. (58). As shown in Figure C1, the approximation that the escape fraction is given by eqn. (59) works well when fe 1%, and gives a decent order of magnitude estimate when fe is lower (such cases are observationally unimportant). In all simulations, the center to edge extinction (scattering + absorption) optical depth is constant, τ = 10. Three different values of g were simulated (circles); from darkest to lightest g = 0, 1/2, −1/3. We computed N 0 (g) for each value of g: