Radiative Transfer Effect on Ultraviolet Pumping of the 21cm Line in the High Redshift Universe
Abstract
During the epoch of reionization the 21cm signal is sensitive to the scattering rate of the ultraviolet photons, redshifting across the Ly resonance. Here we calculate the photon scattering rate profile for a single ultraviolet source. After taking into account previously neglected natural broadening of the resonance line, we find that photons approach the resonance frequency and experience most scatterings at a significantly smaller distance from the source than naively expected , where is the initial frequency offset, and the discrepancy increases as the initial frequency offset decreases. As a consequence, the scattering rate drops much faster with increasing distance than the previously assumed profile. Near the source, ( comoving Mpc), the scattering rate of photons that redshift into the Ly resonance converges to . The scattering rate of Ly photons produced by splitting of photons that redshift into a higher resonance (Ly, Ly, etc.) is only weakly affected by the radiative transfer, while the sum of scattering rates of Ly photons produced from all higher resonances also converges to near the source. At , on scales of 0.01–20Mpc (comoving), the total scattering rate of Ly photons from all Lyman resonances is found to be higher by a factor of than obtained without full radiative transfer. Consequently, during the early stage of reionization, the differential brightness of 21cm signal against the cosmic microwave background is also boosted by a similar factor.
Subject headings:
cosmic microwave background – cosmology: theory – diffuse radiation – intergalactic medium – radio lines: general1. Introduction
The next generation of radio telescopes (e.g., LOFAR, MWA, SKA, and 21CMA)^{1}^{1}1Information on these telescopes can be found at http://www.lofar.org/, http://www.haystack.mit.edu/ast/arrays/mwa/, http://www.skatelescope.org/, and http://cosmo.bao.ac.cn/project.html/, respectively. promises to open a first observational window into the epoch preceding the end of reionization at . By measuring the redshifted 21cm signal from neutral hydrogen, the new telescopes can provide us with information on the history of reionization, the nature of the first radiation source, the spectrum of the primordial density perturbation field, the cosmological parameters, the physical properties of dark matter particles, etc. (e.g., Madau et al., 1997; Ciardi & Madau, 2003; Loeb & Zaldarriaga, 2004; Nusser, 2005; McQuinn et al., 2006; Barkana & Loeb, 2005; Chuzhoy et al., 2006; Shchekinov & Vasiliev, 2006).
The 21cm signal can be observed either in absorption or emission against the cosmic microwave background (CMB), when the hydrogen spin temperature, , is different from the CMB temperature, . The former is defined by the relative populations of the hyperfine states of hydrogen atoms
(1) 
where GHz is the frequency of hydrogen hyperfine transition. In the highdensity gas clouds, which so far have been the only detected sources of the 21cm signal, the decoupling of from is done by collisions between atoms, which induce direct transitions between the hyperfine states and couple the spin temperature to the hydrogen kinetic temperature, . In the intergalactic medium (IGM), on the other hand, after the density becomes too low for collisional coupling to be effective and the decoupling of from can be effectively done only by scatterings of Ly photons, which likewise couple to (Wouthuysen, 1952). The spin temperature is therefore a weighted function of and (Field, 1958)
(2) 
where is the collisional coupling constant (which we neglect throughout the paper). The radiative coupling constant, , is
(3) 
where is the oscillator strength of the Ly transition, K, is the spontaneous emission coefficient of the hyperfine transition and is the intensity at Ly resonance, when the backreaction on the incident photons caused by resonant scattering is neglected. For the unperturbed IGM the backreaction correction, , is (Chuzhoy & Shapiro, 2006)
(4) 
The Ly photons can be produced in several ways, including recombinations, line cooling, and collisional excitation of atoms by nonthermal electrons produced by Xrays (Chuzhoy et al., 2006; Chen & MiraldaEscudé, 2006). However, in case the reionization epoch was dominated by stellar ultraviolet (UV) sources (which at present is the most popular theory) most Ly photons in the neutral IGM originate as photons between Ly and Lylimit frequencies. Due to expansion of the Universe, a fraction of these photons (those emitted between Ly and Ly) would gradually redshift into the Ly resonance. The rest would redshift into one of the higher resonances (generally the one which is just below their initial frequency). The scatterings of high resonance photons produce electron excitations to the () states, which are shortly followed by deexcitations and emission either of the original photon (in case the cascade goes directly to the ground state, ) or of several lower energy photons (in case the cascade goes via some intermediate state). Since between and of cascades follow the latter route, after just a few scatterings most of the high resonance photons will be split. Depending on whether the last cascade goes via the or the resulting photons will or will not include a Ly photon.
In this paper, we calculate the Ly scattering rate as a function of distance from the UV radiation source. Contrary to the naive expectation, the scattering rate does not evolve as . While is a good approximation at low redshifts when hydrogen is mostly ionized, at high redshifts where the GunnPeterson optical depth to the Ly photons is extremely large, we find that the scattering rate profile becomes much steeper. Therefore, until the radiation intensity reaches the saturation levels (i.e., ), the fluctuations of the 21cm signal would be significantly stronger than previously estimated.
The paper is organized as following. In § 2, we describe our calculation of radiative transfer for “continuum” Ly photons (i.e., photons that gradually redshift into the Ly resonance). In § 3, we describe the radiative transfer of “injected” Ly photons (i.e., photons produced by splitting of high resonance photons). In § 4, we study the total scattering rate from “continuum” and “injected” Ly photons. In § 5, we summarize and discuss our results. Throughout the paper, we adopt a spatiallyflat CDM cosmological model with matter density and baryon density in units of the critical density and a Hubble constant , which is consistent with the constraints from the WMAP data (Spergel et al., 2007).
2. Radiative transfer – “Continuum” photons
The fate of the UV photon depends on the frequency at which it was emitted. Most photons originally emitted between Ly and Ly frequencies will travel a large distance, up to Mpc (comoving), before being scattered by one of the hydrogen atoms. As photon redshifts closer to the Ly resonance its scatterings become more frequent and the mean distance between subsequent scatterings rapidly drops (see Figure 1). Therefore almost all scatterings occur within a very small region^{2}^{2}2We find that at above 99% of the scatterings occur within kpc comoving.. Thus to make an accurate estimate of the scattering rate, it is in practice sufficient to count only scatterings occurring within this region, where photons are very close to the resonance and the number of scatterings is of the order of GunnPeterson optical depth. However, to derive the scattering rate profile, , one also needs to calculate the distance of this region from the radiation source, which by contrast is determined mainly by the first few scatterings before photons redshift into the Ly resonance.
In general, photon scattering crosssection is the function of both its frequency and the gas temperature. Neglecting hyperfine splitting, the crosssection without thermal broadening can be written as
(5) 
with the line profile being a Lorentz function
(6) 
where is the linecenter frequency, is the oscillator strength and is the spontaneous decay rate ( for Ly ). Thermal broadening introduces a core in the crosssection around the line center with a width of the order of thermal velocity in units of light speed , which is about for the unheated IGM in the redshift range we are interested. In most calculations of the Ly scattering profile around first sources, the line profile is simplified to a Dirac function and scatterings occur as a photon redshifts into linecenter frequency. We show in this paper that the frequency dependence of the crosssection affects the scattering rate profile.
Since in the high optical depth regime, the first scatterings happen when photon frequency is still significantly above resonance (i.e., in the Lorentz wing), the distance between the source and the point where photon enters into the resonance is nearly independent of gas temperature. Hence the normalized scattering rate profile is also independent of temperature [though not the total number of scatterings, see eq. (4)]. Since the line frequency offset relevant for our study is much larger than the quantum width and the thermal core width, the crosssection [eq. (5)], being in the wing regime, can be well approximated as
(7) 
where is formally the scattering crosssection at twice the linecenter frequency. We note that in this study, the initial frequency of a photon is blueward of line center, i.e., . This is different from that in Loeb & Rybicki (1999). They investigate the brightness and spectral distributions of escaping Ly radiation around sources before reionization, and in their case photons start at a frequency slightly redward of the resonance, i.e., .
As a photon travels in the neutral medium with Hubble expansion, its frequency redshifts. For a photon with initial frequency offset , the scattering optical depth to a distance is^{3}^{3}3We reduce the problem to a calculation for photons traveling in a medium with Hubble expansion velocity field with fixed Hubble constant at a given redshift. Strictly speaking, this approximation is only accurate if the distance traveled is much less than the Hubble radius. Such a condition is slightly broken for the largest scales in the calculation, of the Hubble radius (see Fig. 2), which would introduce slight distortions in the scattering rate profiles on these scales. However, for our main purpose, the scales in consideration are generally much smaller and the above approximation is always valid.
(8) 
where and are the neutral hydrogen number density and Hubble constant at the redshift in consideration, and . In what follows, we mainly adopt as the distance variable, which proves to be convenient. At high redshifts, the conversion from to comoving distance is simply
(9) 
The parameter can be regarded as a redshift variable for a given cosmology. Formally, it is the optical depth to the Hubble radius seen by a photon with frequency twice the linecenter frequency. For Ly photons, we have
(10) 
where the mass fraction of helium is taken to be one quarter.
Equation (8) can be inverted to find the distance traveled by a photon of initial frequency offset for a given optical depth ,
(11) 
For the solution to approach , the one corresponding to a Dirac crosssection, the frequency offset needs to be large and/or the parameter needs to be small. In the regime that , is proportional to , which implies that on sufficiently small scales the distance where most scatterings occur deviates from the expectation with crosssection. The dependence on in Equation (11) also means that, for a given initial frequency offset, the place where most scatterings occur is no longer at a single distance and instead it has a distribution.
The spherical symmetry and Hubble velocity field allow a simple Monte Carlo calculation of the distribution of Ly photons from continuum between Ly and Ly as a function of distance from the central source. The procedure is as follows:

Start a photon from the center () with frequency offset drawn according to the slope of the UV spectrum of the central source.

Draw a scattering optical depth according to the exponential distribution.

Find the distance of traveling before scattering according to equation (11).

Determine the position of scattering, . That is, , where is drawn uniformly between and 1.

(i.e., redshifted) and .

Repeat ii.–v. until approaches the width of the thermal core (a few times for situations here). Then start from i. again until the desired number of photons have been drawn.
The cosmology and redshift are fully encoded in a single parameter [eq. (10)]. In the following discussions, we will assume for Ly photons, which corresponds to . In step iv., we simply assume the direction after scattering is isotropic. We have tested that a more realistic angular distribution, such as a dipole distribution, has little effect on the resultant spatial distribution of Ly photons.
In Figure 1, we show the spatial distribution of Ly photons, , after scatterings for two initial frequency offsets , based on the above Monte Carlo procedure. The distribution has been normalized so that . In each panel, the dotted curve, dashed curve, filled circles, and solid circles are the distributions around the central source after the 1st, 2nd, 10th and 100th scattering, respectively. As photons redshift towards the line center, the scatterings become more frequent and the distribution approaches an asymptotic one, (thick solid curve). Hence, in practice, to derive the total scattering rate of Ly photons, it is sufficient to perform the Monte Carlo simulation for the first few scatterings.
If the Ly crosssection is assumed to be a Dirac function, for an initial photon, all the scatterings would happen at a distance . However, this is not the case if the frequency dependence of the crosssection is taken into account, as shown in Figure 1 (thick solid curves). The position where most of the scatterings happen is broadly distributed at distances smaller than . Furthermore, the peak of the distribution depends on the initial frequency, the smaller the initial frequency offset , the farther away of the peak from the position expected from a function crosssection. For example, as shown in Figure 1, for , the peak of the distribution is at a distance of that expected from a function crosssection, while for it is only . As already mentioned, this dependence on initial frequency can be understood by considering the first scattering with equation (11): for a given optical depth , becomes increasingly smaller than as becomes smaller. Equation (11) also shows that for smaller , the scattering position would become close to that expected from a function crosssection, which we discuss more about in the next section.
To get the total scattering rate of the “continuum” Ly photons for a general UV source, one needs to integrate over the frequency range between Ly and Ly
(12) 
where is the GunnPeterson optical depth to Ly photons, is the backreaction correction factor [eq. (4)], and is the luminosity of the source in terms of number of photons per unit frequency per unit time.
The above dependence of the scattering position on the initial frequency implies that the radial profile of the scattering rate would be steeper than the drop. With the Monte Carlo technique, we calculate the scattering rate profile from “continuum” photons, assuming a flat UV spectrum of the central source (“flat” here means equal number of photons per unit frequency per unit time). In Figure 2, the dashed blue curve is the result from the Monte Carlo simulation, while the dotted blue curve is the corresponding curve of the drop with the same normalization of the UV spectrum of the central source. The true profile is steeper than and the amplitude can be an order of magnitude higher at small distances.
Close to the source the scattering rate scales as . This behavior can be explained by the following argument. The mean free path of the photon emitted with frequency offset (i.e., the average distance between subsequent scatterings) scales as [eq.(11) in the limit of ]. The change of frequency offset between two subsequent scatterings is simply proportional to . Therefore the number of scatterings that would occur before the photon frequency moves closer to resonance and its mean free path drops significantly, scales as . Since the photon direction changes almost at random after each scattering, the total distance it travels until it redshifts into the resonance scales as . Conversely, we can say that photons reaching the resonance within distance from the source, are emitted with frequency within from Ly resonance. Since the total number of photons within this frequency range increases in direct proportion with , the photon scattering rate, which is proportional to the number of photons reaching the resonance within some region divided by its volume, scales as .
3. Radiative transfer – “Injected” photons
In addition to photons originally emitted between Ly and Ly frequencies, which redshift into the Ly resonance, Ly photons are also produced by splitting of photons originally emitted between Ly and Lymanlimit frequencies. These photons will be first scattered as they redshift into the closest resonance. In case, following their absorption by a hydrogen atom, the excited electron cascades directly to the ground state, the original photon would be reemitted. Alternatively, the electron can cascade via some intermediate state, in which case the original photon is split into several photons. Depending on the path of the cascade, the cascade products may include the Ly photon. The fraction of Ly photons made up by cascade of high resonances is typically less than of the total (Hirata, 2006; Pritchard & Furlanetto, 2006; Chuzhoy & Shapiro, 2007). However, as the distance traveled by photon before it redshifts into the closest resonance is of order , the “injected” Ly photons are produced within much smaller volume and within Mpc (comoving) from the source outnumber the “continuum” Ly photons (Chuzhoy et al., 2006).
Since the “injected” Ly photons are injected directly into the resonance (hence their name), to derive their scattering rate profile it is sufficient to follow the path of their high resonance progenitors. The Monte Carlo simulation is similar to what we perform for “continuum” photons, but we need to use crosssections for higher resonance lines [i.e., different and parameters in eq. (11) for different lines] and record the probability of Ly production at each scattering. In the calculation, we use the data in Table 1 of Pritchard & Furlanetto (2006) for the probability of decay from the state to the ground state and the probability of producing a Ly photon through cascades from the level (also see Hirata, 2006). The scattering crosssection for higher resonance Lyman lines are calculated from oscillator strength and Einstein coefficient for transition listed in Morton (2003) (and extrapolations are used for higher ).
The total scattering rate of Ly photons produced by cascade from state with initial frequency is
(13) 
where is the normalized spatial distribution for Ly photons that experience scatterings, is the probability that the original Ly photon is destroyed, is the chance that the destruction of the Ly photon results in production of the Ly photon (i.e., that the electron cascade goes via rather than state), and other symbols have the same meaning as in equation (12). Since at each scattering there is a significant probability, , that the original Ly photon is destroyed, at large values of , approaches zero. That is, the total distribution of “injected” Ly photons can converge by considering those photons that experience only a small number of scatterings.
Figure 3 illustrates the evolution of Ly photon distribution for higher resonance photons (Ly and Ly, with initial frequency offset ) that experience no more than scatterings (1, 2, 10, and 20) before destroyed. As expected, the distribution approaches quickly to the total distribution (solid curve), similar to the “continuum” case. At the same initial frequency offset, higher resonance photons lead to a peak position closer to that expected from function crosssection. The reason is that a higher resonance line has a smaller crosssection [thus a smaller parameter in eq. (11)].
In contrast to “continuum” photons that typically reach the resonance after multiple scatterings, a photon emitted with frequency blueward of a higher Lyman resonance can produce a resonant Ly photon just after a single scattering. If such a photon has mean free path as the “continuum” photon case (see § 2), the distance it travels before producing a Ly photon is just and one would expect the scattering rate profile . However, this can only formally happen extremely close to the central source [, see eq.(11)] because of the fast drop of the corresponding , and proximity effect would take over at such scales. Since both the oscillator strength and the spontaneous decay rate drops as at large , the characteristic scattering optical depth in equation (11) drops as . Even for , has dropped to a level. The low value of means that, at a given initial frequency offset, the scattering rate distribution for injected photons is almost the same as that expected from function crosssection with the peak position approaching . For scales of interest (i.e., in the regime of ), the mean free path instead of for higher resonance photons, thus each individual follows the profile.
The solid magenta curve in Figure 2 shows the Ly scattering rate profile from injected photons for the same flat UV spectrum of the central source as in the “continuum” case. Each of the steps in the curve reflects the place that a higher resonance line starts to contribute to produce injected Ly photons. At a given initial frequency offset, the scattering rate of injected Ly photons is dominated by those produced by higher level resonance lines. As mentioned above, the location of the injected Ly photons from these resonance lines is close to the expectation of form crosssection as the corresponding decreases fast. Therefore, the true scattering rate profile is quite close to that from form crosssection, which is a sum of a series of truncated functions (dotted magenta curve in Figure 2, almost on top of the solid magenta curve). At a distance , the true profile is only higher by 4% in this particular case.
Although each individual follows a truncated profile, the overall profile from summing over all “injected” photon contributions no longer follows a profile because the truncation place for the individual component depends on . It can be shown that at small scales, the overall profile (Chuzhoy et al., 2006), interestingly the same dependence as the “continuum” photon case (see § 2). To see this, we note that, with , the individual component is simply for a flat source spectrum (in photon number per unit frequency per unit time). The slope of the overall profile is then . Since , we have and in the large limit. Therefore, the overall profile , which is what is seen in Figure 2.
4. The Total Scattering Rate
The total scattering rate of Ly photons is obtained after summing over all Lyman series:
(14) 
The solid black curve in Figure 2 shows the profile of the total scattering rate. For the case shown here (), “continuum” (“injected”) photons dominate at distances to the center greater (less) than (Mpc comoving). At smaller distance, even though the amplitude from “continuum” photons can be an order of magnitude higher than the naive rate given by the form crosssection, the domination of “injected” photons reduces the difference in the total scattering rate. As for the slope, the overall scattering rate is slightly steeper than the sum of a series of truncated drops (dotted black curve). To the first order, the overall scattering rate is about 40% higher than the naive calculation in a wide range of distances (see the thick black curve in the lower panel).
The calculation of the scattering rate profile becomes much easier by assuming a form scattering crosssection. From the above example, we see that to the first order, the exact profile that takes into account the frequency dependence of the crosssection and the radiative transfer effect can be obtained by applying a correction to that from the calculation with crosssection. This correction is almost a constant, increasing the amplitude by a few tens of percent. Apparently, the correction factor depends on redshift, larger at higher redshift because of an increase in the parameter [eq. (11)]. We perform simulations at a series of redshifts and find that the fractional difference between the exact and the naive calculations of the scattering rate can be approximated as a constant at a given redshift,
(15) 
This approximation underestimates (overestimates) the correction at small (large) scales. For example, at , =44% (32%) at (). A more accurate fit that accounts for such a tilt is
(16) 
which works almost perfectly for and (roughly corresponding to 0.01–20 Mpc comoving). The first term on the right hand side always dominates and the second term accounts for the slight difference in the slope of the two profiles. With this fitting formula, the profile from the naive calculation can be corrected by multiplying the scattering rate by a factor of .
All the above examples assume a flat UV spectrum, with between Ly and Lyman limit, where is the luminosity in terms of number of photons per unit frequency per unit time. The spectrum slope depends on the nature of the central source. For example, for Population III stars, , while for Population II stars (Barkana & Loeb, 2005). Given the narrow frequency range between Ly and Lyman limit, we do not expect a strong dependence of the scattering rate profile on . In the case of adopting the crosssection, the scattering rate for an individual component is simply , where the relation is used. The departure from the drop with respect to the flat spectrum is introduced through and and both are small given the narrow frequency range of Lyman series and .
We perform simulations for two cases with quite different spectral slopes, and . In Figure 4, we plot ratios of scattering rates for the two cases. The central source luminosity is normalized so that the numbers of photons emitted between Ly and Lyman limit for the two cases are the same. That’s why the case has a higher rate from “injected” photons (thin solid curve) but lower rate from “continuum” photons (dashed curve). The ratios for the individual components are similar to the crosssection results (dotted curves). However, the steepening of the slope of the scattering rate for the “continuum” photons at small scales, which increases its contribution to the total scattering rate, reduces the difference in the total scattering rates for and (thick solid curve). Even with the large difference in the spectral slope, the difference in the overall scattering rates is only at a level less than 10% for . A comparison between cases using spectral slopes of Population II () and Population III () stars shows a pattern similar to what is seen in Figure 4 with the difference in the overall scattering rates being at a level of 5%.
5. Summary and Discussion
We investigate the effects of frequency dependence of scattering crosssection and radiative transfer on the scattering rate profile of Ly photons around a central UV source, which is relevant for the pumping of the 21cm line in the high redshift universe.
Because of the frequency dependence of the scattering crosssection, a photon between Ly and Ly frequency (“continuum” Ly photon) experiences a small number of wing scatterings before it redshifts to the core frequency and starts core scattering. Core scatterings, which determine the scattering rate, happen at distances much less than that expected from a scattering crosssection. For scatterings of higher resonance photons that produce “injected” Ly photons, the effect is weak owing to the fast drop in the crosssection. We have shown that for a single UV source the scattering rate profile of “continuum” Ly photons is significantly steeper than (previously expected from scattering crosssection), approaching near the source. The scattering rate profile of “injected” Ly photons from an individual high Lyman resonance closely follows a truncated profile, while the overall profile from all higher Lyman resonances coincidentally also approaches near the source. At , the total scattering rate from “continuum” and “injected” photons is 30–50% higher than the naive calculation that adopts crosssections, on scales of 0.01–20Mpc (comoving). We also find that, when radiative transfer effects are properly accounted for, the slope of the scattering rate profile does not have a strong dependence on the spectral slope of the central source.
In our calculations we have assumed that during the early stages of reionization the gas temperature changes adiabatically. If instead, the gas was significantly heated, the Doppler core for Ly photons would be increased and their scattering rate profile would be even steeper. However, since the temperature of neutral hydrogen does not exceed K, the scales above kpc (comoving) would be virtually unaffected.
As a consequence of the higher scattering rate with respect to the naive calculation, the spatial fluctuations in the pumping radiation field produced by multiple UV sources would be significantly higher as well. During the early stages of reionization, when the UV intensity is relatively low, the differential brightness of the 21cm signal scales almost linearly with the scattering rate of the Ly photons and the corrections to the scattering rate translate into similar corrections to the 21cm brightness. Such a correction would increase the size of the expected Ly spheres around first sources (e.g., Cen, 2006; Chen & MiraldaEscudé, 2006). Moreover, if the radiative transfer effect is overlooked, the fluctuation power spectrum estimated from 21cm observations would be skewed on scales up to a few tens of Mpc, which would lead to inaccurate inference of the matter fluctuation power spectrum (such as the amplitude and the running of the spectral index). Therefore, until the epoch when the intensity of pumping radiation reaches a saturation level (if such epoch in fact exists), the correct interpretation of the 21cm signal requires taking into account radiative transfer (mainly for “continuum” Ly photons), as described in this paper.
References
 Barkana & Loeb (2005) Barkana, R., & Loeb, A. 2005, 626, 1
 Cen (2006) Cen, R. 2006, ApJ, 648, 47
 Chen & MiraldaEscudé (2006) Chen, X., & MiraldaEscudé, J. 2006, ArXiv Astrophysics eprints, arXiv:astroph/0605439
 Chuzhoy et al. (2006) Chuzhoy, L., Alvarez, M., & Shapiro, P. R. 2006, ApJ, 648, L1
 Chuzhoy & Shapiro (2006) Chuzhoy, L. & Shapiro, P. R. 2006, ApJ, 651, 1
 Chuzhoy & Shapiro (2007) Chuzhoy, L. & Shapiro, P. R. 2007, ApJ, 655, 843
 Ciardi & Madau (2003) Ciardi, B., & Madau, P., 2003, ApJ, 596, 1
 Field (1958) Field, G. B. 1958, Proc. IRE, 46, 240
 Hirata (2006) Hirata, C.M. 2006, MNRAS, 367, 259
 Loeb & Rybicki (1999) Loeb, A., & Rybicki, G. B. 1999, ApJ, 524, 527
 Loeb & Zaldarriaga (2004) Loeb, A., & Zaldarriaga, M. 2004, Phys. Rev. Lett., 92, 211301
 Madau et al. (1997) Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429
 McQuinn et al. (2006) McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815
 Morton (2003) Morton, D. C. 2003, ApJS, 149, 205
 Nusser (2005) Nusser, A. 2005, MNRAS, 364, 743
 Pritchard & Furlanetto (2006) Pritchard, J. R., & Furlanetto, S. R. 2006, MNRAS, 367, 1057
 Shchekinov & Vasiliev (2006) Shchekinov, Yu. A., & Vasiliev, E. O. 2006, astroph/0604231
 Spergel et al. (2007) Spergel, D. N., et al. 2007, ApJS, 170, 377
 Wouthuysen (1952) Wouthuysen, S. A. 1952, AJ, 57, 31