The Wouthuysen-Field effect in a clumpy intergalactic medium

We show that, due to the high optical depth of the intergalactic medium to Lyman-alpha photons before the Epoch of Reionization, the Lyman-alpha scattering rate responsible for the Wouthuysen-Field effect from an isolated source will be negligible unless (1) there is sufficient time for the scattering photons to establish a steady state, or (2) the scattering gas is undergoing internal expansion or has a peculiar motion of tens to hundreds of km/s away from the source. We present steady-state solutions in the radiative diffusion approximation for the radiation field trapped in a clump of gas and show that this may result in an enhancement, by a factor of up to 10^6, of the strength of the Wouthuysen-Field effect over that obtained from the free-streaming limit. Solutions to the time-dependent diffusion equation, however, suggest that the timescales required to reach such a steady state will generally exceed the source lifetimes. In the presence of internal expansion, a steady state may be established as photons are redshifted into the red wing, and significant enhancement in the scattering rate may again be produced. Alternatively, a substantial scattering rate may arise in systems with a peculiar motion away from the source that redshifts the received radiation into the resonance line centre. As a consequence, at epochs z<30, when collisional decoupling is small except in dense regions, and prior to the establishment of any large-scale diffuse radiation field of resonance line photons, the 21cm signature from the Intergalactic Medium produced by the Wouthuysen-Field effect will in general trace the peculiar velocity field of the gas in addition to its density structure.


INTRODUCTION
Following the Recombination Era at z ≃ 1100, the baryons produced in the Big Bang were largely neutral. By z ∼ > 6, the spectra of high redshift Quasi-Stellar Objects (QSOs) show that the hydrogen had become ionized (Becker & et al. 2001). The history of structure formation in the wide redshift expanse between these epochs, when there were few if any sources of radiation, is largely unknown. In principle, 21cm emission from intergalactic hydrogen could reveal the evolution of structure in the baryons during these cos-⋆ Scottish Universities Physics Alliance mic "Dark Ages." Because the baryons are cold, they would closely trace the evolution of the dark matter, so that 21cm imaging could trace the growth of structure following the Recombination Era (Hogan & Rees 1979;Scott & Rees 1990). Because of strong coupling between the hyperfine structure of the hydrogen and the Cosmic Microwave Background (CMB), however, a mechanism that decouples the spin temperature of the hydrogen hyperfine structure from the CMB temperature must be active, otherwise the hydrogen is rendered invisible: it absorbs and re-emits the CMB radiation at the same rate, leaving it indistinguishable from the CMB. At redshifts z > 30, collisions between hydrogen and electrons and other hydrogen atoms are adequate to begin de-coupling the spin temperature from the CMB (Scott & Rees 1990;Madau et al. 1997). Except in dense regions, however, at later times collisional decoupling is inadequate.
Once the first sources of radiation begin to turn on, the scattering of Lyα photons by the neutral hydrogen offers an alternative means of decoupling, through the Wouthuysen-Field effect (WFE) (Wouthuysen 1952;Field 1958). As sources turn on in sufficient number to begin reionizing the Intergalactic Medium (IGM), they will produce a combined intensity of Lyα photons sufficient for coupling the spin temperature to the light temperature of the Lyα photons rather than the CMB, rendering the intergalactic hydrogen visible against the CMB in either emission or absorption (Madau et al. 1997).
The discovery of the End of the Dark Ages (EDA) and the onset of the Epoch of Reionization (EoR) has become one of the paramount goals of a new generation of radio telescopes, such as the LOw Frequency Array (LOFAR) 1 , the Murchison Widefield Array (MWA) 2 , the Primeval Structure Telescope/21 Centimeter Array (PaST/21CMA) 3 , the Precision Array to Probe EoR (PAPER) 4 , and a possible Square Kilometre Array (SKA) 5 . Reviews of this rapidly growing area are provided by Loeb & Barkana (2001), Fan et al. (2006) and Furlanetto et al. (2006). Madau et al. (1997) estimated the Lyα collision rate Pα that drives the WFE as the integrated intensity from cosmologically distributed sources, assuming that photons blueward of Lyα emitted from a source will contribute their full amount to Pα once they redshift into the local Lyα resonance. In fact, this provides a lower limit (assuming photons are not destroyed by, eg, dust absorption). The multiple scattering of resonance line photons will in general enhance the radiation field. For pure Doppler scattering, Field (1959a) argued the rate is enhanced by the number of scatters a photon undergoes before being randomly scattered sufficiently redward of line centre to escape. The estimate neglected scattering in the Lorentz wings, however. It also assumed that the IGM takes part in the homogeneous and isotropic expansion of the Universe. In fact the IGM is clumpy, with structures breaking away from the cosmological expansion. Efforts are underway to estimate the radiation field and scattering rate in a cosmological context through Monte Carlo simulations of the scattering of resonance line photons in an inhomogeneous medium (Chuzhoy & Zheng 2007;Pierleoni et al. 2007;Semelin et al. 2007). In this paper, we use analytical means to explore some of the consequences of intergalactic structure formation for the WFE as a means of generating an intergalactic 21cm signature.
In the next section, we discuss the optical depth of the IGM to resonance line photons and the implications for the scattering rate. Approximate steady-state solutions to the radiative transfer equation are derived in § 3, and timedependent solutions in § 4. A discussion of the solutions and applications are provided in § 5, along with a summary of our conclusions. 1 www.lofar.org 2 www.haystack.mit.edu/ast/arrays/mwa 3 web.phys.cmu.edu/∼past 4 astro.berkeley.edu/∼dbacker/eor 5 www.skatelescope.org The curves are for T = 10 K (solid), 100 K (dashed) and 1000 K (dot-dashed). Also shown is the approximation τν ≃ x 1 /x for the T = 10 K model (dotted line, nearly coinciding with the solid line). (Lower panel) The value of x at which τν = 1 for z S − 1 < z < z S . The curves are labelled as in the upper panel.

THE WFE IN THE INTERGALACTIC MEDIUM
The optical depth through a homogeneous and isotropic expanding IGM of a photon emitted by a source at redshift zS and received at redshift z at frequency ν > ν0, where ν0 is the resonance line frequency, is (Field 1959a;Gunn & Peterson 1965) where n l (z ′ ) is the number density of scattering atoms in the lower level at epoch z ′ , σ = πe 2 f lu /(mec) ≃ 0.0110 cm 2 Hz is the total resonance line cross section, where f lu ≃ 0.4162 is the upward oscillator strength for hydrogen Lyα, ϕV (a, ν) is the Voigt line profile normalized to R dνϕV (a, ν) = 1, a ≃ 0.0472T −1/2 is the ratio of the decay rate to the Doppler width ∆νD = ν0b/c, where b = (2kBT /mH) 1/2 is the Doppler parameter for hydrogen gas at temperature T and c is the speed of light, and lp is the proper path length. In the Lorentz wing, expressed as a function of x = (ν−ν0)/∆νD, the dimensionless Voigt profile φV (a, x) = (∆νD)ϕV (a, ν) is well-approximated as φV (a, x) ≃ a/(πx 2 ). The differential proper line element evolves according to dlp/dz ≃ (c/H0)Ω −1/2 m (1 + z) −5/2 in a flat universe at redshifts for which Ωm(1 + z) 3 dominates the contribution from the vacuum energy, where Ωm is the ratio of the total mass density to the critical Einstein-deSitter density, and where H0 = 100h km s −1 Mpc −1 is the Hubble constant today.
Photons emitted by the source at a frequency νe > ν0(1 + zS)/(1 + z) will scatter in the blue Lorentz wing (except for those for which the inequality is nearly an equality, in which case they will scatter through the Doppler core). In the blue wing the optical depth due to Lyα scattering by neutral hydrogen is τν ≃ 1.75 × 10 5 cn l (0) H0Ω where y = ν/ν0 and u = y(1 + zS)/(1 + z), and all quantities are assumed to be in cgs units. For 0 < xb/c << (zS−z)/(1+ z) << 1, the expression simplifies to where x1 = a π σc ν0 n l (0) for a mean line-of-sight IGM temperature TIGM for a universe with Ωm = 0.3 and baryon density Ω b h 2 = 0.022 (O'Meara et al. 2006;Spergel & et al. 2007). As an illustration, τν(x) for a source at zS = 8 is shown at z = 7 in Fig. 1 (upper panel). The optical depth is extremely high until ν is well displaced from the line centre. The values of x1, at which τν = 1, are shown in the lower panel of Fig. 1 for a range of redshifts. Photons emitted by the source with frequencies between ν0 < νe < ν0(1 + zS)/(1 + z) will pass through the resonance line frequency en route to the gas at z. Upon passing through the resonance line, the photons will essentially be completely scattered out of the line of sight. As a consequence, essentially no photons will be received in the frequency range ν0(1 + z)/(1 + zS) < ν < ν0 (the Gunn-Peterson effect). A general expression for the optical depth through the Doppler core is provided in Appendix A.
Photons emitted sufficiently redward of ν0 to avoid the Doppler core may still be scattered in the Lorentz red wing of the surrounding hydrogen. This will produce a halo of scattered Lyα photons around a central source like a QSO (Loeb & Rybicki 1999). For frequencies x < 0 near line centre, Loeb & Rybicki (1999) show that the photons will build up an energy density as they diffuse through the surrounding IGM of the form uν ∼ (−x) −9/2 e −[r/rc(x)] 2 , where rc(x) = (−x) 1/2 (2/3)[(b/H)λ mfp (x)] 1/2 scales like the harmonic mean between the distance over which Hubble expansion produces a velocity difference matching the Doppler parameter b and the scattering mean free path of a resonance line photon in the Lorentz wing where n l is the number density of atoms in the lower level. The energy density is exponentially suppressed for x → 0 − at a fixed position, resulting in a negligible Lyα scattering rate except very near the QSO. Numerically, so that a substantial Lyα scattering rate is confined to a small region around the source. The region is so small in fact, that the assumption that the surrounding IGM takes part in the Hubble flow will not be valid. More likely the gas will have been ionized and violently disturbed by the QSO itself. Finally, photons emitted sufficiently redward of ν0 at the source to escape scattering in the Lorentz wing will pass freely through the IGM, where they will be received at redshift z at frequencies The vast majority of photons from the source emitted blueward of, but near, the resonance line frequency never arrive near the line centre, having been scattered out of the path along the way. The resonance line radiation impinging on the absorption site 6 from the source scatters in the Lyα resonance line at the rate per neutral atom in the lower state where hP is Planck's constant and uν is the local specific energy density of photons received at z from the source at zS with specific luminosity Lν where rL is the luminosity distance between the source and the gas at redshift z. The scattering rate assuming none of the photons have been scattered out of the line of sight en route to the absorption site (the free-streaming limit) is The ratio of the scattering rate for τν ≃ x1/x to the freestreaming value is then given by where Ta is the temperature of the absorbing gas. The low value poses a fundamental problem to the effectiveness of the WFE as a means of decoupling the spin state of the gas from statistical equilibrium with the CMB. Only if photons have sufficient time to diffuse across the line centre will the scattering of Lyα photons be an effective means of decoupling the spin temperature from the CMB in comoving objects.

Steady-state solutions without atomic recoil
The assumption of a homogeneous and isotropic universe, while a fair approximation over large scales, breaks down on small scales: the IGM is clumpy at the redshifts from which the 21cm signal is expected (Tozzi et al. 2000;Gnedin & Shaver 2004). For a sufficiently short mean free path, photons will be trapped within a discrete system and diffuse in frequency and space as they scatter within the system. In this section, we consider the energy distribution the radiation will reach given adequate time to relax to a steady state. For a Lyα photon to become trapped within an absorber of lengthscale La, the mean free path must minimally satisfy λ mfp < La, for which τν = La/λ mfp > 1. As the optical depth in the wing may be expressed as τν = π −1/2 (aτ0)x −2 , where τ0 = n l Laπ −1/2 σ(∆νD) −1 is the optical depth at the line centre, this corresponds to photons escaping with x > π −1/4 (aτ0) 1/2 (Osterbrock 1962). Radiative transfer solutions for sources embedded within a slab suggest this may be somewhat too restrictive for systems with very high line centre optical depths, in which case the photons escape through spatial diffusion rather than frequency diffusion. Steady-state solutions in the diffusion approximation for scattering in the Lorentz wings in very high optical depth systems show that photons escape the slabs with frequencies typically on the order of (aτ0) 1/3 (Adams 1972;Harrington 1973;Neufeld 1990). An escape frequency of xesc = x * (aτ0) 1/3 corresponds to a mean free path smaller than the thickness of the slab by the factor , where Ta is the temperature of the absorber and N19 is the column density in units of 10 19 cm −2 . (This corresponds to a region 24 kpc across for the mean intergalactic hydrogen density at z = 8, and smaller for overdense systems.) As f depends only weakly on Ta and NHI, while x * depends on aτ0 and the radiation source parameters of any particular problem, we express the escape frequency simply as noting that f will in general vary with the depth within the absorber, and must be solved for according to each particular configuration, but will be of order 0.01 − 1 for applications we consider. Comparison of Eq. (12) with Eq. (4) for x1 shows that photons with τν < 1 in the IGM will generally not be trapped within a structure comoving with the expansion of the Universe unless the structure has a neutral hydrogen column density N19 > N19,crit, where If x1 < xesc, the absorber will behave as a photon bucket, allowing the energy density of the radiation field at frequencies x < xesc to build up with time as radiation from the source becomes trapped. As the photons scatter within the absorber, they will diffuse in frequency and space. Because photons with |x| > xesc will escape, the energy density of the radiation trapped within the absorber will reach a steady state when the rate of incoming photons balances the rate of escape. In principle, the energy density that develops could produce a substantial scattering rate, even exceeding the rate estimated assuming no photon losses from the source due to scattering by the intervening IGM. The depth to which photons at frequency x1 penetrate an absorber depends on the optical depth of the absorber. If it is optically thin at x1, the photons will stream through.
If optically thick, the photons will scatter near the surface of the absorber. Photons blueward of x1, however, will penetrate more deeply, as the cross section diminishes like 1/x 2 . In general, photons will be injected at varying layers within the absorber, with bluer photons injected at increasing depths.
We would like to estimate the evolution of the radiation field within an absorber. No solution to this problem exists in the literature. The full radiative transfer problem through the slab is an involved one, as the radiation field will vary with both the optical depth through the slab and with frequency. It is particularly important in our application to obtain the solution across the Doppler core to estimate the Lyα scattering rate. Instead of solving the full problem, we seek approximate solutions in which we treat the scattered radiation field as locally isotropic, neglecting the spatial diffusion of the radiation through the absorber. We show in Appendix B that the resulting radiative transfer equation in the diffusion approximation is formally identical to the equation of Harrington (1973) for the problem of a slab with a uniform distribution of spectrally flat sources. The solution we obtain neglecting the spatial diffusion of the radiation agrees well with Harrington's solution which includes both spatial and frequency diffusion. A detailed quantitative description will require a more exact solution to the problem, but our approximate approach allows us to explore qualitatively several effects of interest. We would like to estimate the energy density that may be achieved as a function of xesc. We also wish to examine the effect of internal motions on the radiation energy density across the line centre and the possible role of atomic recoil. Lastly, we would like to estimate the time it takes the radiation field to establish a steady state. All of these effects extend well beyond those explored by existing slab solutions.
We estimate the radiation density that builds up within an absorber by assuming a steady state between the rate of photons injected into the absorber and the rate of diffusion of the photons across the escape frequency, imposing the boundary condition u(x) = 0 for |x| > xesc, where u(x) = ∆νDuν . We assume isotropic scattering within a homogeneous absorber, and use the diffusion approximation to describe the resonance line scattering of the photons. In equilibrium, where n(x) = ∆νDnν , nν = uν/hPν,S(x) = (∆νD) 2 S(ν)/(n l cσ), where S(ν) is a source function, and D(a, x) is the diffusion coefficient. The second order moment of the frequency redistribution function gives for the diffusion coefficient D(a, x) = φV (a, x) + (1/3)d 2 φV (a, x)/dx 2 (Rybicki & dell'Antonio 1994). The derivation of the diffusion equation, however, does not conserve photon number to the order of the approximation. As a consequence, the diffusion coefficient is left ambiguous. Rybicki & dell'Antonio (1994) advocate adopting D(a, x) = φV (a, x) for its simplicity, and we shall do so here. The source S(ν) represents the rate at which photons received from the source are incident on the absorber at x < xesc, where they are trapped and diffuse in frequency. It is given by the scattering emissivity allowing photons with |x| > xesc = 200 to escape, as given by the diffusion approximation, Eq. (14), with D(a, x) = φ V (a, x). An absorber temperature Ta = 100 K is assumed.
where R(ν ′ , ν) is the photon redistribution function for resonance line scattering (Mihalas 1978). Since scattering in the wings is coherent, The dependence on x results in a sharp peak at x = x1/2, so that the sourceS(x) is well approximated bỹ where the normalization A = 2.5 × 10 −6 h(1 + z) −3/2 T 1/2 IGM [(Lν /hPν)/4πr 2 L ] (with all quantities expressed in cgs units), matches the source rate of Eq. (16) integrated over all frequencies. Note that the dependence on the absorber temperature cancels. Here, the frequency at which photons are injected into the absorber is represented generally as xinj to allow for any peculiar motion of the absorber, which could shift the peak of the received radiation either redward or blueward, and to allow for scattering within the absorber which will shift the peak of the radiation field penetrating to deeper layers further toward the blue. For a comoving absorber, xinj = x1/2 near the surface. The solution to Eq. (14) using the wing approximation for φV (a, x) is The solution is shown in Fig. 2. We find that adopting D(a, x) = φV (a, x) + (1/3)d 2 φV (a, x)/dx 2 produces a nearly identical solution, agreeing to within a fraction of a percent for xinj well out of the core. Numerically integrating the diffusion equation using the full Voigt profile produces a solution that also differs negligibly from Eq. (18). The ratio of the collision rate P l , given by using Eq. (18) in Eq. (8), to the free-streaming collision rate P l (0) of Eq. (10) is then where we have used the approximation Eq. (A1) below for the Voigt profile with xesc > xinj ≫ xm assumed, and neglected the terms following the leading. The enhancement may be substantial. For N19 = N19,crit, P l,ss /P l (0) ≃ 1.5 × 10 4 h −2 (1 + z) 3 T −1 IGM (1 − x 3 inj /x 3 esc ), so that even for a warm IGM temperature of TIGM = 1000 K a substantial boost in the scattering rate will result.
The solution Eq. (18) assumes velocities internal to the cloud are negligible. To estimate the effect internal motion may have on the frequency distribution of the photons, we approximate the motion as an isotropic expansion (or contraction) within the cloud characterised by a uniform velocity gradient Ha (which in general will not be the same as the Hubble parameter). In the diffusion approximation, Eq. (14) is modified to (Rybicki & dell'Antonio 1994)  tures are typically contracting, but the internal motions along a filament may expand along the filament in regions (Zhang et al. 1998). For an overdensity ten times the cosmic mean at z = 8, an internal velocity gradient of 100 km s −1 over 100 kpc corresponds to γ ≃ 1.7 × 10 −7 . Contraction by the same magnitude corresponds to γ ≃ −1.7 × 10 −7 . The solution for the sourceS(x) = AδD(x − xinj) is where gγ (x, y) = (A/γ)fγ (xesc, xinj)fγ (x, y)/fγ (xesc, y) and fγ (x, y) = 1 − exp[− 2πγ 3a (x 3 − y 3 )]. The solution is extremely sensitive to γx 3 esc /a. For γ ≫ γcrit ≡ 3a/(2πx 3 esc ), a broad wing of amplitude A/γ redward of xinj will result as photons are redshifted across the line centre. The distribution will cut off sharply at x > xinj. For an absorber with Ta = 100 K (a ≃ 0.00472), and xesc = 200, γcrit ≃ 2.8 × 10 −10 . The solutions for γ = 1.7 × 10 −7 with (xinj, xesc) = (100, 200) and (25, 50) are shown in Fig. 3. A substantial decrease in n(x) is produced compared with the γ = 0 case. A significant enhancement of the Lyα collision rate over the free-streaming limit, however, may still result, with P l,ss /P l (0) ≃ a/(πγx1). For γ < 0, a broad wing blueward of xinj will form as photons are blueshifted in the contracting gas. A complex profile, however, may form redward of xinj as redward diffusion is resisted by blueshifting. Such a case is illustrated in Fig. 3 for γ = −1.7 × 10 −7 with xinj = 25 and xesc = 50.

Steady-state solution with atomic recoil
The above solutions neglect the effect of atomic recoil on the evolution of the radiation field. Since the photon momentum (∼ hPν0/c) is small compared with the momentum of the hydrogen atoms (∼ mHb), the effect should be small, but may result in a significant heating rate when the gas is still cold (Madau et al. 1997 where the recoil parameter ǫ = hPν0/(2kBT mHc 2 ) 1/2 ≃ 0.025T −1/2 is the ratio of the momentum of a resonance line photon to the momentum of an atom moving at the thermal velocity (Rybicki & dell'Antonio 1994;Meiksin 2006). Eq. (22) expresses only the first order correction in ǫ. The solution for a sourceS(x) = AδD(x − xinj) is where gǫ(x, y) and fǫ(x) = (1 − 2ǫx + 2ǫ 2 x 2 ) exp(2ǫx). Here, the wing approximation for φV (a, x) is adopted, but we find a nearly identical result allowing for the Doppler core as well.
The solution for xinj = 100, xesc = 200 and ǫ = 0.0079, corresponding to T ≃ 10 K, is shown in Fig. 4. The result differs substantially from the solution for ǫ = 0. In fact the solution cannot be trusted, as it is not consistent with the original order of Eq. (22) in ǫ. This may be demonstrated by adding to 2ǫφV (a, x)n(x) the quantity ǫ 2 y(x), where y(x) is an unspecified function representing the higher order corrections to the recoil, and which in general may also depend on n(x). Inserting the series solution n(x) = Σ ∞ i=0 ǫ i ni(x) into Eq. (22) and equating equal orders in ǫ recovers n0(x) as given by Eq. (18), and produces the first order equation (d/dx)[2φV (a, x)n0(x) + D(a, x)dn1(x)/dx] = 0 and the second order equation (d/dx)[2φV (a, x)n1(x) + y(x) + D(a, x)dn2(x)/dx] = 0. The first order equation gives n1(x) in terms of n0(x), subject to the boundary conditions n1(−xesc) = n1(xesc) = 0. The second order equation shows that n2(x) depends on the unspecified (neglected) function y(x). Thus the solution is valid only up to n(x) = n0(x) + ǫn1(x). If higher order terms are substantial, then the solution for n(x) cannot be trusted. The solution n0(x) + ǫn1(x) for the example above differs substantially from the solution n(x), and indeed is not even positive-definite, as shown in Fig. 4, demonstrating the solution is not valid.
Including higher order terms involves a non-trivial expression for the redistribution function including non-linear terms in the recoil parameter (Basko 1981;Meiksin 2006). (Additional higher order terms in the recoil parameter arise from relativistic corrections, but these are reduced by factors of b/c.) For temperatures T > 100 K, the discrepancy between the full solution and the first order solution is smaller than about 20% for xinj = 100 and xesc = 200. In general, the correction will be small for ǫxesc ≪ 1, and we caution that if this is violated the steady-state solutions may not be valid. We do not pursue further consequences of the atomic recoil solution here.

The role of bulk peculiar velocities
Absorption systems with a bulk peculiar velocity component away from the source will see the photons redshifted into line centre. In this case, the scattering rate of resonance line photons may be appreciable even in an optically thin absorber (N19 < N19,crit), so that the build-up of a strong radiation field at line centre through photon capture and Figure 5. Surface of constant Lyα scattering rate P l produced by a continuum source to the left within a collapsing halo of neutral hydrogen. The inner surface terminates at the accretion shock. The circle represents the turnaround radius rta of the halo, with the solid line representing the hemisphere of gas with peculiar velocity redshifted relative to the source and the dashed line the hemisphere with blueshifted peculiar velocity. The surface shown corresponds to a maximum radius of 0.59rta. The figure is tilted by 5 • to reveal the 3D structure of the surface. diffusion is not necessary for significant WFE decoupling of the spin state from the CMB. The peculiar velocity required to bring the absorber into the peak of the source function S(ν) in Eq. (16) is given by vpec ≃ 1 2 x1ba, where x1 is given by Eq. (4). For pure thermal broadening (assuming negligible internal velocity structure within the absorber), this corresponds to vpec ≃ 14h −1 (1 + z) 3/2 " Ta TIGM This is comparable to the typical peculiar velocities of nonlinear cosmological structures at the epochs of interest. As a result, the WFE will produce a 21cm signature that maps the peculiar velocity field of the structures within the IGM.
As an illustration, we consider the signal produced by a collapsing halo, such as around a forming galaxy or galaxy cluster, illuminated by an external radiation source, as depicted in Fig. 5.
The collision rate from a source of specific luminosity Lν on an absorber with Doppler parameter ba at luminosity distance rL moving away from the source with peculiar velocity v is using Eq. (3), approximating φV (a, x) ≃ π −1/2 exp(−x 2 ), noting that the resulting exponential in the integrand peaks sharply near x = v/ba assuming v0/ba ≃ x1 or larger, Taylor expanding the argument of the exponential to second order in x − v/ba, and expressing the peculiar velocity as v = v0 cos(θ), the projected accretion velocity along the line of sight to the source. In general, v0 will vary with distance from the centre of the collapsing halo. As a definite case, we approximate the peculiar velocity of the accreting gas using the self-similar solution of Bertschinger (1985) for adiabatic accretion of a γ = 5/3 gas of negligible density onto a dark matter halo with a comoving centre in an Einstein-de Sitter cosmology. The velocity field in the region between the accretion shock and the turn-around radius at cosmological time t is wellapproximated by v0(r) ≃ rta t Vmin log 10 λs log 10 λ, where rta is the turn-around radius, λ = r/rta, λs = rs/rta ≃ 0.3472, where rs is the radius of the accretion shock, and Vmin ≃ −1.433 characterises the inflow velocity of the gas just before passing through the accretion shock. Contours of constant values of P l correspond to values of constant projected accretion velocity. If at a radius r towards the source (θ = 0), the accretion velocity is a fraction f of its maximum infall velocity, then contours of fixed P l will correspond to the surfaces extending over the angular range 0 ≤ θ ≤ arccos(f ).
The surface corresponding to f = 0.5 is shown in Fig. 5, with a maximum angle θmax = 60 • and maximum radius λ0 = λ f s = 0.59. The loci of constant P l correspond to arcing sheets. The exponential sensitivity of P l to v for v0/ba ≃ x1 (see Eq. [25]) ensures a surface of narrow width will dominate the 21cm signal produced. Such a surface could be mistaken for the region of neutral hydrogen just outside an H II region surrouding a source of ionizing radiation. A varying gas temperature, and therefore varying ba, in the accreting region will further complicate the signal. If v0/ba >> x1, then the entire solid hemisphere will be illuminated by the source (except for a thin circular wafer at θ < ∼ 90 • ). A similar effect may apply to an expanding void. Since voids expand faster than the Hubble expansion, the peculiar velocity on the far side of a void from a distant source could redshift the received radiation into line centre. A shell of direct scattering may then result even without including the reddening effects of radiative transfer within the void.

THE EVOLUTION OF THE WFE IN DISCRETE OBJECTS
It was assumed for the solutions in § 3 that the photons are able to achieve a steady state as they scatter within the absorbing system. In this section, we compute the timescales required to establish a steady state by solving the timedependent diffusion equation where τ = t/ts, and ts = ∆νD/(n l σc) ≃ 3.2n −1 l T 1/2 a s is the scattering time of the resonance line photons at line centre.
We solve Eq. (28) using a Crank-Nicholson scheme. We tested the scheme against the time-dependent solutions for a  constant δ-function source at line centre in an infinite homogeneous medium and a flat continuum source in the diffusion approximation (Rybicki & dell'Antonio 1994). The results applied to the sourceS(x) = δD(x−xinj) with xinj = 100 and boundary condition n(x) = 0 for |x| > xesc = 200 are shown in the upper panel of Fig. 6. The photon frequency distribution takes on the approximate shape of the final steady-state solution Eq. (18) as it converges towards it. The evolution of the photon frequency distribution at line centre (x = 0) is shown in the lower panel of Fig. 6. For sources that have been active only a small fraction of the time required to achieve full convergence, the scattering rate may be negligibly small.

DISCUSSION AND CONCLUSIONS
The role played by the Wouthuysen-Field effect in decoupling the spin temperature of the neutral hydrogen from the CMB temperature depends most crucially on the scattering rate of Lyα photons. The spin temperature TS is generally a weighted mean of the colour temperature Tα, kinetic temperature TK of the gas, and brightness temperature TR of any incident radiation (such as the CMB) at the 21cm frequency, where yα ≡ P10 A10 T * Tα and yc ≡ C10 A10 T * TK are the Lyα and collisional pumping efficiencies, respectively (Field 1958;Madau et al. 1997), A10 = 2.85 × 10 −15 s −1 is the spontaneous decay rate of the 21cm transition, and T * ≡ hPν10/kB, where ν10 is the frequency of the 21cm transition. Here P10 = 4P l /27 is the indirect de-excitation rate of the triplet hyperfine state induced by Lyα photon scattering. The colour temperature of the radiation field is the harmonic mean temperature Tα = 1/ T −1 u (ν) , where Tu(ν) = −(hP/kB)(d ln uν /dν) −1 , weighted by uν ϕV (a, ν) for an energy density uν of resonance line photons (Meiksin 2006). The coefficient C10 is the de-excitation rate of the triplet hyperfine state induced by atomic collisions.
For a Lyα scattering rate P l exceeding the critical thermalization rate where TCMB ≃ 2.73(1 + z) (Mather & et al. 1994) is the CMB temperature, the spin temperature will be driven to the colour temperature in the absence of strong collisional de-excitation (Madau et al. 1997). The radiation field will typically rapidly thermalize with the neutral hydrogen, resulting in Tα → TK (Field 1959b;Meiksin 2006). The neutral hydrogen will then appear in either absorption or emission against the CMB, depending on whether TK is below or above TCMB, respectively.
Before the Epoch of Reionization, when the first radiation sources turn on at the End of the Dark Ages, regions around the sources may become visible through their 21cm signature. Because the photons produced by a source, such as a QSO or star-forming galaxy, near the Lyα resonance will be scattered out of the line of sight over short distances, the scattering rate of Lyα photons will be much smaller than P th except extremely near the source. Only photons emitted sufficiently blueward of the Lyα resonance frequency will survive to any distant structures, where they will be received well in the blue Lorentz wing. If a structure has a sufficient optical depth in the wing to the photons received, it will trap them and rescatterings within the structure will allow the photons to diffuse in frequency across the resonance line. The radiation field may be built up in this way until it reaches a steady state in which the rate of incoming photons is balanced by the rate at which photons diffuse sufficiently far from the resonance line frequency to escape in the wings where the structure is optically thin. The resulting photon scattering rate in the diffusion approximation may greatly exceed the estimate P l (0) in the free-streaming limit, for which it is assumed no resonance line photons emitted by the source have been scattered out of the line of sight before reaching the absorber.
We identify several factors, however, that may reduce the energy density of the resonance line photons from the steady-state diffusion approximation estimate for static structures: (1) Internal expansion within the structure, characterised by the dimensionless expansion parameter γ, will redshift the photons into a red wing, where they will escape, on a timescale much faster than the diffusion time, resulting in a much reduced energy density. In the diffusion approximation, the scattering rate will be reduced compared with the γ = 0 case by the factor (3a/πγ)/(x 3 esc − x 3 inj ) for γ ≫ γcrit = 3a/(2πx 3 esc ), where xesc is the dimensionless escape frequency at which photons escape the structure, and xinj is the frequency at which photons are predominantly received from the source. Contraction within the structure may altogether prevent photons from diffusing across the resonance line frequency.
(2) The scattering rate may be reduced by atomic recoil. Estimating the magnitude of this effect, however, generally requires including non-trivial higher order terms in the recoil parameter ǫ than first order in the steady-state equation. We caution against introducing atomic recoils into Monte Carlo calculations and other steady-state solutions, as the results obtained may not be physically self-consistent if the solutions differ substantially from those obtained for ǫ = 0.
(3) The time scale to establish a steady state will typically exceed the duration of high star formation rates in galaxies or the lifetime of QSOs. As a consequence, the collision rate will be reduced by factors of tens to thousands compared with the steady-state value.
As an illustration, consider the region around a bright QSO with specific luminosity at the Lyman edge of 10 31 ergs s −1 Hz −1 pre-heated by soft x-rays ahead of the ionization front to a mean temperature of TIGM ≃ 1800 K at a distance of 10 Mpc (Madau et al. 1997). According to Eq. (4), x1 ≃ 200 at z = 8 (assuming h = 0.7). A structure 250 kpc across with a moderate overdensity of 3 and internal temperature of Ta ≃ 120 K will have a neutral hydrogen column density of NHI ≃ 3 × 10 20 cm −2 and escape frequency xesc ≃ 200, from Eq. (12). For a power law spectrum ν −1.5 , Eq. (10) gives for the Lyα collision rate in the free-streaming limit P l (0) ≃ 8.7 × 10 −13 s −1 , about an order of magnitude smaller than the thermalization rate P th . The steady-state value from Eq. (19) then gives P l,ss ≃ 10 −8 s −1 , well in excess of P th . The time to establish a steady state, however, is long. The scattering time of Lyα photons in the structure is about ts ≃ 8.5 × 10 4 s. According to Fig. 7, the time for the radiation field to achieve a steady state is then about 10 17 s, or 3 Gyr, greatly exceeding the expected lifetime of the QSO. If the QSO had illuminated the structure for as long as 0.1 Gyr, according to Fig. 6 the steady-state value must be reduced by a factor of 0.05. The resulting scattering rate will still exceed P th . If the system has a peculiar velocity away from the QSO, the rate could be quite different. Eq. (24) shows that a peculiar velocity of 140 km s −1 , a plausible value, would shift the peak of the radiation from the source received by the structure to the resonance line frequency. The radiation field would then rapidly reach its steady-state value, and P l,ss would greatly exceed P th .
A much more overdense structure at the same distance from the QSO would be able to achieve the steady-state scattering rate. For example, a virialized structure 35 kpc across with an overdensity of 200 and temperature Ta ≃ 1000 K would have a neutral hydrogen column density of 3 × 10 21 cm −2 and xesc ≃ 220. The time for the radiation field to reach a steady state would then be about 0.2 Gyr, a plausible lifetime for the QSO.
By contrast, it may be difficult to achieve a steady state in a void. For an IGM temperature TIGM ≃ 100 K, x1 = 825 at z = 8. A region 3 Mpc across with a temperature of Ta = 10 K, and underdense by a factor of 3 will have a neutral hydrogen column density of 4 × 10 20 cm −2 , xesc ≃ 800 and a scattering time at line centre of ts ≃ 2.2 × 10 5 s. The time to reach a steady state is then 2100 Gyr. Even if the void is illuminated for as long as the age of the Universe, the radiation field will be negligibly small compared with its steady-state value. Since the void will expand, however, the incident radiation will be redshifted across the resonance line. Rybicki & dell'Antonio (1994) estimate a characteristic time for redshifting to dominate diffusion in an infinite homogeneous medium to be τγ ≃ (a/γ 4 ) 1/3 . Voids expand somewhat faster than the mean Hubble expansion. Adopting the Hubble constant at z = 8 for h = 0.7, γ ≃ 5 × 10 −6 . A steady state should be established for t ≫ τγ ts ≃ 2 × 10 4 yr, and the steady-state scattering rate will then match the freestreaming value, with P l,ss /P l (0) ≃ a/(πγx1) ≃ 1, which may be adequate for producing a 21cm signature in the presence of sufficient sources.
The examples above show that the effect of peculiar motions may lead to an enhancement in the Lyα scattering rate over the free-streaming limit and produce the required decoupling of TS from TCMB. Even if the scattering rate is inadequate for inducing decoupling throughout most of the structure, it may be adequate along narrow loci within which the velocity field either permits the resonance line photons to accumulate near the resonance line frequency or shifts the resonance line frequency of the absorber to matching the incident radiation field from the source as filtered through the IGM. We presented an example of an accreting halo which produces a 21cm signature along a curved surface for which the impinging photons are received at the resonance line frequency due to the infall velocity of the accreting gas. Such a curved structure could be mistaken for a rim of neutral gas outside an H II region in radio maps of the IGM. Similar curved surfaces may form on the far side of voids from a distant source.
As a consequence of the above effects, the 21cm signature of the IGM will not simply trace the density structure of the IGM, but will trace the peculiar velocity structure as well. We have not considered multiple sources. In the presence of a large number of nearby sources, a greater fraction of the gas will be able to match the incident radiation field through their peculiar motions. In addition, the escaping radiation from the absorbers will contribute to the ambient radiation field. In due course, a metagalactic diffuse radiation field may grow sufficiently strong to ensure the WFE is capable of decoupling the spin temperature from the CMB temperature throughout the IGM, and so ensure the production of a 21cm signature everywhere. The effects described here suggest that in the early stages, before any such metagalactic field has been established (and it has yet to be demonstrated that such a field will be established), the 21cm signatures produced at the End of the Dark Ages, and possibly extending into the Epoch of Reionization, may only be interpreted with a thorough understanding of the evolution of the energy density of resonance line photons within a clumpy medium, including the effects of peculiar motions.