Stability of smallscale baryon perturbations during cosmological recombination
Abstract
In this paper, we study smallscale fluctuations (baryon pressure sound waves) in the baryon fluid during recombination. In particular, we look at their evolution in the presence of relative velocities between baryons and photons on large scales (), which are naturally present during the era of decoupling. Previous work concluded that the fluctuations grow due to an instability of sound waves in a recombining plasma, but that the growth factor is small for typical cosmological models. These analyses model recombination in an inhomogenous universe as a perturbation to the parameters of the homogenous solution. We show that for relevant wavenumbers the dynamics are significantly altered by the transport of both ionizing continuum ( eV) and Lyman photons between crests and troughs of the density perturbations. We solve the radiative transfer of photons in both these frequency ranges and incorporate the results in a perturbed threelevel atom model. We conclude that the instability persists at intermediate scales. We use the results to estimate a distribution of growth rates in random realizations of largescale relative velocities. Our results indicate that there is no appreciable growth; out of these realizations, the maximum growth factor we find is less than at wavenumbers of . The instability’s low growth factors are due to the relatively short duration of the recombination epoch during which the electrons and photons are coupled.
I Introduction
The early universe is largely composed of atomic matter, or baryons, radiation and cold dark matter. The main resources available to study this era are the Cosmic Microwave Background (CMB) and Large Scale Structure (LSS). Primary anisotropies of the CMB are a result of the imprint of primordial fluctuations on radiation at early times Planck Collaboration et al. (2013a); Larson et al. (2011), while secondary anisotropies probe the matter distribution at late times Planck Collaboration et al. (2013b). LSS surveys are a complementary probe of the clustering of matter at late times Dawson et al. (2013).
The distribution and dynamics of baryons during early epochs of the Universe is poorly constrained by this data. The angular distribution of power in the CMB constrains them on large scales, through their coupling with the radiation and its effect on the Baryon Acoustic Oscillations. The CMB is welldescribed by a spectrum of adiabatic fluctuations at these scales – these are motions of both the baryon and radiation fluids. Tight bounds exist on the primordial fluctuations of solely the baryon fluid at these scales – the socalled isocurvature modes Planck Collaboration et al. (2013a).
This paper deals with the complementary limit of fluctuations in the baryon field on very small scales. In the CMB, this information is lost due to diffusion damping of the anisotropies. The spectral distortion associated with diffusion damping has been suggested as a probe of modes on these small scales Chluba et al. (2012). The proposed PRISM mission aims to study CMB spectral distortions André et al. (2014).
In the rest of this paper, we use the term “matter” to refer to baryons, for reasons of readability; we are not concerned with the dynamics of cold dark matter. We concentrate on smallscale fluctuations of the matter field, and their evolution through the epoch of recombination. In particular, we undertake a detailed study of an instability which can amplify subJeans length fluctuations at recombination suggested by Shaviv Shaviv (1998). The mechanism of interest is potentially applicable to wave numbers in the range Mpc comoving. This is at much smaller scales than the standing acoustic waves responsible for peaks in the CMB power spectrum and baryon acoustic oscillations in the matter power spectrum (e.g. Sakharov (1965); Eisenstein and Hu (1998)), which are damped below the Silk scale Silk (1968) Mpc. We expect the prerecombination amplitudes of modes at to be extremely small, but if an instability is present then a “seed” amplitude could be generated by nonlinear generation of smallscale isocurvature modes Press and Vishniac (1979), or even thermal fluctuations if the growth rate is fast enough.
Shaviv’s instability acts on sound waves propagating in a partially ionized gas, in the presence of a background flux of radiation. The scenario is illustrated in Fig. 1. The key observation is that the fraction of ionized atoms is different in overdense and underdense regions; the ionization fraction, , is lower in overdense regions where recombination proceeds faster due to the increased flux of free electrons seen by the ionized atoms.
Sound waves are propagating longitudinal waves in the matter fluid – if we orient ourselves along the wavevector, , the local velocity at a compression is in the forward direction, while the opposite is true for rarefactions. Thus, the earlier observation leads to a negative correlation between the ionization fraction and the local velocity in the region of propagation.
In the presence of a background flux of radiation in the matter’s bulk restframe, the radiative force acting on a mass element is related to the radiation flux, or alternatively its velocity , by the opacity, which is proportional in turn to the ionization fraction, . Over a timeperiod of the sound wave, the resulting force per unit mass performs an amount of work given by
(1) 
where in the second equation, the multiplicative factor involving the energy density of the radiation (), its interaction cross section with matter (), the particle mass (the hydrogen mass ) and the speed of light relates the force per unit mass to the ionization fraction. The net work done over a time period is nonzero due to the difference in ionization fractions during the forward and backward motion. From consideration of Fig. 1, the work integral of Eq. (1) is positive if the flux, , is directed opposite to the wavevector, .
The first estimate of the growth rates due to this mechanism, due to Shaviv Shaviv (1998), used the assumption of local thermal equilibrium (LTE) to derive the variations in the ionization fraction. Recombination in the real universe proceeds out of LTE, and most of the hydrogen first recombines to excited states before reaching the ground state Peebles (1968); Zel’dovich et al. (1969); Rybicki and dell’Antonio (1994); Hu et al. (1995); Dubrovich and Stolyarov (1995); Seager et al. (2000); AliHaïmoud and Hirata (2011); Chluba and Thomas (2011). Subsequent work Liu et al. (2001); Singh and Ma (2002) used the three level approximation to model nonLTE recombination, and incorporated the diffusion of microwave background photons, following which the expected growth rates were revised downward.
The standard treatment of recombination assumes that the ionization state is set by the local radiation field. This is valid in the homogenous case, since the transport of photons out of the region of interest is perfectly balanced by the influx from other regions. This is no longer true in the inhomogenous case, and these two components (the influx and outflux) do not balance each other. In particular, direct recombinations to the ground state, which did not affect the homogenous ionization fraction, , are important in determining its fluctuation, .
In this paper, we incorporate the transport of both continuum and Lyman () photons. We find simple analytical expressions for this “nonlocal” contribution to the evolution of the ionization fraction, and provide revised estimates for the growth rates of the smallscale soundwaves.
The paper is organized as follows: In Section II, we expand upon the simple estimate given above for the work done on the fluctuations, and estimate the associated growth rates. In Section III, we list the relevant background variables, and the various factors which determine their size during the epochs of interest.
In the subsequent sections, we write down equations of motion for the smallscale fluctuations. We start with the standard Newtonian equations for the density and velocity in Section IV. We estimate growth rates using a simple scaling relation for the ionization fraction fluctuation in Section V. We then move beyond this simple treatment, and study in detail the radiative transport of photons between different parts of the fluctuations – Sections VI and VII deal with the transport of continuum and Lyman photons respectively.
Finally, we bring all the pieces together and estimate the growth rates of the smallscale fluctuations in Section VIII, and find their distribution due to a stochastic background of largescale relative velocities in Section IX. We finish with a short discussion of our results and their implications in Section X. We collect some details which lie outside the main line of analysis, but provide some physical intuition, into the appendices.
Ii Motivation and simple estimate
This section closely follows the analysis of Shaviv (1998).
We use the two fluid approximation, where matter and radiation fluids are coupled by Thomson scattering of photons off free electrons. The characteristic response time, , is inversely related to the matter’s opacity per unit mass, . For a given relative velocity between the two fluids, , the force per unit mass is expressed in terms of the response time as
(2) 
where is the photon flux seen in the matter’s rest frame. This force, and the related response time, are most easily obtained by considering the Doppler shifted background radiation field in the matter’s rest frame. The result is Weymann (1965)
(3) 
where is the hydrogen ionization fraction, is the Thomson scattering crosssection and is the radiation energy density constant. The matter temperature, closely follows the radiation temperature, , at these times. With this understanding, we omit the subscript on the temperature in subsequent equations.
Primordial adiabatic fluctuations entering the horizon lead to largescale motions of the matter and radiation fluids. Their physical size, is at recombination. Due to the small but finite response time, , during this epoch, the matter velocity does not perfectly follow the local radiation velocity; this leads to a spectrum of relative velocities that can be estimated from the background cosmology Tseliakhovich and Hirata (2010).
We consider motions of the matter fluid alone, as contrasted with the largescale adiabatic modes involving both matter and radiation. In particular, we concentrate the evolution of very small wavelength modes though the epoch of recombination out to late redshifts of . We consider modes that are isothermal in nature, i.e., have a uniform matter temperature. As noted in the discussion (Section X), this condition restricts our analysis to modes with wavenumbers smaller than . The large scale adiabatic modes are effectively fixed on the timescales relevant to these smallscale modes, and provide a background radiation flux due to their associated relative velocity. The radiative force due to this flux is given by Eq. (2).
The ionization fraction and opacity vary with the local density during recombination. Thus smallscale fluctuations of the matter density are associated with a modulation of the of the local force, denoted by . The inphase component of feeds power from the largescale relative motions into smallscales.
The rest of this section estimates the size of this effect in a simplified scenario with direct recombination to the ground state of neutral hydrogen. With this assumption, the ionization fraction is given by the Saha equilibrium value, which we denote by . This is set by the balance between the recombination of free electrons to the ground state, and photoionization by microwave background photons.
(4) 
where is the ionization energy of a hydrogen atom in the ground state, and is the hydrogen number density. We take the logarithm of both sides of Eq. (4), and perturb it to estimate the powerlaw exponent relating the perturbed free electron fraction and hydrogen density as follows
(5) 
where we have used the assumption that the smallscale fluctuations do not perturb the temperature, . The Saha electron fraction is approximately at recombination, so the exponent .
Consider a region with a background relative velocity between matter and radiation, . The associated force per unit mass, , is related to the relative velocity by the response time , according to Eq. (2). The local matter density, velocity and force per unit mass are perturbed due to the smallscale fluctuation. For a sound wave, these perturbations are of the form
(6a)  
(6b)  
(6c) 
In the above relations, denotes the isothermal sound speed. It is determined by the matter temperature according to . Averaged over the phase of the wave, the power input into the fluctuation by the extra force, of Eq. (6c), is
(7)  
(8) 
The first line uses Eq, (6b) and (6c) for the velocity and force respectively, while the second uses Eq. (2) for the background force and the definition of the response time in Eq. (3). The energy per unit mass in the fluctuation is . Hence the growth rate for the amplitude, , can be estimated from the input power of Eq. (8) as
(9) 
The growth of the instability is maximal during the epoch with large relative velocities and moderate response times. Relative velocities of the order of the isothermal sound speed are needed to produce an appreciable growth rate. The last part of Section III deals with the distribution of large scale relative velocities in detail. In particular, Figure 3 shows the mean relative speed, and the isothermal sound speed, as a function of redshift, . We see that large relative velocities are much more probable in the postrecombination era; however, this effect is mitigated by the growing response time. Ultimately, the instability is limited by the relatively narrow duration of cosmic recombination.
Iii Background parameters
This section describes the relevant properties of the background on which the small fluctuations of interest live.
We assume a standard spatially flat cold dark matter cosmology with the Planck cosmological parameters Planck Collaboration et al. (2013c). The derived quantities of interest to us are the hydrogen number density and ionization fraction, and the relative velocities between matter and radiation on large scales due to adiabatic fluctuations of primordial origin.
The simplest of these to obtain is the hydrogen number density, , which is given by
(10) 
where is the Baryon fraction and is the Helium mass fraction.
It is considerably harder to estimate the hydrogen ionization fraction, , as a function of redshift. It is especially challenging to follow it through the epoch of recombination, when the universe transitions from a plasma of free electrons and hydrogen nuclei to a largely neutral phase with traces of free electrons that are strongly coupled to the cosmic microwave background (CMB) radiation.
This difficulty arises from the fact that direct transitions to the ground state of hydrogen contribute very little to recombination, since they produce ionizing photons themselves. Instead, recombination mainly proceeds through excited states of neutral hydrogen. In order to derive the evolution of the ionization fraction to subpercent level accuracy, we should follow the populations of a large number of excited states of the hydrogen atom AliHaïmoud and Hirata (2011); Chluba and Thomas (2011).
We eschew this sophisticated analysis for a conceptually simpler, and less accurate, model of recombination originally proposed in Refs. Peebles (1968); Zel’dovich et al. (1969). This is adequate for the purposes of this paper, since we follow fluctuations in the ionization fraction. The errors introduced in the fluctuations by using the approximate model should be at the fewpercent level.
This model approximates the hydrogen atom as a three level system; it assumes that the excited states of the true hydrogen atom are in thermal equilibrium with each other, and cascade down to the level through fast radiative decays. Atoms in the state reach the ground state when photons redshift through the line due to cosmological expansion, while those in the level deexcite through a twophoton process. Direct recombination via the redshift of continuum photons is much slower (by a factor of ) than through the channel Chluba and Sunyaev (2007). Hence we set the direct recombination’s contribution to zero in the background case. As Section VI demonstrates, this assumption is no longer valid in the perturbed case.
We add the recombination coefficients to the excited states to obtain an effective, or case B recombination coefficient, . We also have an effective rate of photoionization from this state, . With these definitions, the ionization fraction evolves according to
(11) 
where is the Peebles factor, which is the probability that an atom in the state reaches the ground state Peebles (1968). It is defined in terms of the escape rate, the – twophoton transition rate, and the rate of photoionization from the state. We derive explicit expressions for and the population of the state, , in Section VII.1.1.
The case B recombination coefficient and the effective photoionization rate are related by the principle of detailed balance Peebles (1968); Zel’dovich et al. (1969); Seager et al. (2000):
(12) 
We assume that the four sublevels of the level are equally occupied. Thus their occupation fractions are related by . This is justified by the high effective – transition rate at these redshifts ( AliHaïmoud and Hirata (2010, 2011)). This is much faster than both the case B recombination rate per hydrogen atom, and the photoionization rate, which are , and the twophoton decay rate, Goldman (1989).
In deriving Eq. (11), we assume that the population of the level is in steady state, i.e., we balance the net rate of recombination and photoionization against the escape of photons and twophoton decays. This is valid if the abundance of intermediate states is very small; in this case, can be estimated from the recombination codes themselves (e.g. Ref. AliHaïmoud and Hirata (2011)), and is typically of order .
The final piece needed is the spectrum of relative velocities between matter and radiation on large scales. We assume that velocities are irrotational, i.e., they are aligned with their wave vectors, . The velocity at any point in space is a Gaussian random variable, whose twopoint correlation function is
(13) 
where is the dimensionless power per log wavenumber of the component along the wavevector. This power spectrum is given by Ma and Bertschinger (1995)
(14) 
where is the primordial curvature perturbation, and and are the transfer functions for the matter and radiation velocity divergence respectively. We use the publicly available CLASS code to obtain these transfer functions Blas et al. (2011). Figure 2 shows the resulting power spectrum for the relative velocity. We observe that most of the power is in scales near .
We estimate the typical velocities from the distribution of Eq. (13). Figure 3 shows the both the average speed of the matter relative to the radiation, and the isothermal sound speed, as a function of redshift. We observe that these velocities are very small during the prerecombination era: the matterradiation response time, , is much smaller than the expansion age due to rapid scattering, which suppresses the relative velocities. During recombination the free electron fraction drops, and the response time becomes comparable to the expansion age, i.e., recombination leads to decoupling.
Iv Linear analysis of density and velocity fluctuations
Smallscale fluctuations of the matter field perturb the density, velocity and the ionization fraction. We denote the fractional matter overdensity by , the velocity by and the ionization fraction and its fluctuation by and respectively. In addition to these, we denote the perturbed gravitational potential by . In this section, we derive the evolution equations for the density and velocity. In what follows, is the position on a comoving grid, while a dot represents a derivative with respect to coordinate time.
There is a small amount of helium present in the early Universe: the He:H ratio by number, , is given in terms of the Helium mass fraction, , by . We consider late times, , where the helium is fully neutral, so that it does not contribute to the ionization fraction. The hydrogen mass fraction is also used in the equations below.
The matter density, velocity, and gravitational potential on subhorizon scales are governed by the Newtonian equations of motion – the equation of continuity, the NavierStokes equation and Poisson’s equation written in the comoving frame (as in Kolb and Turner (1990)). The linearized forms of these equations are:
(15a)  
(15b)  
(15c) 
The quantity is the Hubble rate of expansion, . The relative velocity force term, , depends on the flux of background radiation in the local matter rest frame. We use Eqs. (2) and (3) to write the force as , where is the inverse of the response time in the case where the hydrogen is completely ionized and the helium mass is neglected. Typical largescale relative velocities, , on comoving scales , appear nearly uniform to the smallscale matter fluctuations. By definition, the latter do not perturb the radiation field. Hence the force associated with the relative velocity is
(16) 
We decompose the velocity into scalar (curlfree) and vector (divergencefree) parts:
(17) 
Under the equation of motion, (15b), the vector part’s evolution depends on the scalar part through the latter’s modulation of the free electron fraction in the force term, but the reverse is not true. We focus on the scalar part in the rest of this paper.
We expand the restoring force due to the pressure up to first order in the fluctuation as follows
(18) 
We substitute the pressure and relative velocity force terms [Eqs. (18) and (16)] in the Newtonian equations [Eq. (15)], and eliminate the gravitational potential, . Assuming planewave forms for the perturbed quantities, , the final forms of the evolution equations for the matter density and velocity are:
(19a)  
(19b) 
V Ionization fraction fluctuation: Saha equilibrium scaling
In order to get a complete picture of the ionization fraction’s evolution, we need to study the transport of photons between different parts of the fluctuations. Before we deal with this problem in Sections VI and VII, we make a simple first estimate following Ref. Shaviv (1998).
The simplifying assumption in this section is that the ionization fraction scales with matter density in the same manner as the value calculated using local thermodynamic equilibrium (LTE, or Saha equilibrium). In subsequent sections, we consider nonequilibrium ionization. We note that perturbed nonequilibrium ionization in cosmology is one of the contributions to the CMB bispectrum and hence has been investigated as a potential contaminant to primordial nonGaussianity studies Novosyadlyj (2006); Khatri and Wandelt (2009); Senatore et al. (2009); Pitrou et al. (2010) and probe of new physics Dvorkin et al. (2013), however these studies did not consider the very high of interest in this paper and hence did not have to solve the nonlocal radiative transfer problem considered in Sections VI and VII.
Using the scaling of Eq. (5) for the ionization fraction fluctuation in Eq. (19b), we reduce the Newtonian evolution equations to
(20a)  
(20b) 
The instantaneous growth rate, , is the largest eigenvalue of the system of Eq. (20). Figure 4 plots this growth rate (normalized to a net elapsed coordinate time, , at the redshift of recombination, ) for various values of the largescale relative velocity, with the wave vector oriented along its direction.
Modes with comoving wavenumbers satisfying (or physical wavelength smaller than ) at recombination are unstable. The growth rate increases with wavenumber, until it saturates on very large wavenumbers: , or physical wavelength , or . The modes at the saturation scale grow by a factor of a few hundred. Since there is a large number of smallscale modes, it is worth considering mechanisms that can cut off the growth on these scales.
Photons in the continuum and line interact strongly with matter during this epoch. We have briefly considered the aspects of this interaction relevant to background recombination in Section III. Continuum photons produced in direct recombinations to the ground state are completely unimportant for the background at the level of accuracy of Section III. Their interaction cross section with neutral hydrogen atoms is so large that they are promptly reabsorbed. However, we should keep track of them in the inhomogenous case, since they can stream from one part of the fluctuation to another.
Figure 5 is a schematic diagram of the radiative transport processes relevant to perturbed recombination. Before we study the various processes in detail in subsequent sections, we clarify a few general points.
Under the assumptions of the three level model of the hydrogen atom, we only need to consider a single spectral line (). This greatly simplifies our analysis. The photons can be decoupled from the continuum due to their wide separation in frequency. In the rest of this paper, we neglect the homogenous population of the first excited state, (except in equations which compute transitions from the level), and assume . As discussed in Section III, it is completely negligible compared to the other populations.
A first step towards judging the relative importance of various arms of Fig. 5 is to look at the mean free paths (MFPs) of the photons at this redshift. If we use numbers for photons at the line center, the comoving wavenumbers corresponding to the MFPs are
(21) 
and
(22) 
Here is the photoionization cross section for a ground state hydrogen atom at the threshold frequency, while and are the Sobolev optical depth and the dimensionless Doppler width of the line at the redshift of recombination.
The MFP for continuum photons is very close to the saturation scale in Fig. 4. Moreover, as we show in Appendix A.1, the length scale for the diffusion of photons is much larger than this naive estimate. In fact, we will see in Section VII that transport is important for wavenumbers satisfying . We begin by studying the outer arm of Fig. 5 in the next section.
Vi Radiative transfer in the continuum
We study the transport of continuum photons in two stages – we first determine their perturbed phase space density, and then calculate its effect on the recombination rate. We approach the problem using the Fourierspace Boltzmann equation (as used in previous sections and in modern CMB codes Ma and Bertschinger (1995); Seljak and Zaldarriaga (1996); Lewis et al. (2000); Lesgourgues (2011)). We note that the similar problem of ultraviolet and Xray radiative transfer in the literature on highredshift 21 cm radiation is usually addressed by a Green’s function approach, i.e. by summing the contributions from individual point sources either analytically or numerically Barkana and Loeb (2005); Pritchard and Furlanetto (2007); Ahn et al. (2009); Holzbauer and Furlanetto (2012).
Let the phase space density (henceforth, the PSD) of continuum photons be . It evolves via the Boltzmann equation
(23) 
The second and third terms on the lefthand side account for the redshift of photons, and their advection respectively. Both the background cosmological expansion and the peculiar matter velocity contribute to the redshift term.
We assume that the PSD is not a dynamical variable and drop the explicit timedependence. This is valid both in the unperturbed and perturbed cases: in the former, because photons redshift through the frequency range much faster than a Hubble time; and in the latter, because the advection term dominates below the Jeans scale.
We neglect the redshift term in Eq. (23). This is equivalent to neglecting the background rate of recombination through the continuum channel. We consider the contributions of the absorption and emission of continuum photons to the right hand side of Eq. (23), and neglect the redistribution of photons within the frequency range due to resonant scattering – this is important within the Lyman lines.
Let , and denote the continuum photon absorption crosssection, the direct recombination coefficient, and the probability distribution for the emitted photons’ frequency respectively. These quantities are functions of radiation (absorption) and matter (recombination) temperature. The integrated or total recombination coefficient to the ground state is defined by
(24) 
The rates of absorption and emission of continuum photons are
(25)  
(26) 
where we have used the fact that every direct recombination is accompanied by the emission of a continuum photon, and multiplied by a factor of to convert the contributions per hydrogen atom to those for the PSD. Substitution in the Boltzmann equation yields
(27) 
In the homogenous case, with just the background parameters, this reduces to the balance between absorption and recombination contributions.
(28)  
(29) 
In the presence of smallscale fluctuations, we linearize the Boltzmann equation and simplify using the unperturbed solution, Eq. (29).
(30) 
Let the total number flux of continuum photons in a direction be . In terms of the PSD, it is given by
(31) 
The photoionization crosssection, , is discontinuous across the threshold frequency. It falls off with increasing frequency in a powerlaw fashion Berestetskii et al. (1971), while the PSD falls in an exponential manner in the UV part of the spectrum. Hence we neglect the frequency dependence of in all integrals. Using Eq. (30) and the definition (31), we get the equation for the transport of the number flux
(32) 
where the coefficients are
(33a)  
(33b)  
(33c) 
Note that the coefficient is the inverse of the mean free path for continuum photons at the threshold for photoionization.
We assume a planewave dependence for the fluctuation, following which the solution to Eq. (32) is
(34) 
The photoionization from and recombinations to the ground state together cause the free electron fraction to evolve as
(35) 
In the homogenous case, we approximate the small contribution to be zero, which gives us a relation between the absorption crosssection and the recombination coefficient.
We can obtain this relation by considering the alternative scenario of local thermal equilibrium (LTE) between a population of ionized and hydrogens, free electrons, and a blackbody distribution of photons above the threshold frequency. The free electron fraction then equals the Saha equilibrium value of Eq. (4). As earlier, we neglect powerlaw frequency dependence of prefactors in the integrals over frequency and obtain the relation
(36) 
In the inhomogenous case, we perturb Eq. (35) and retain terms up to the first order.
(37) 
We use detailed balance in the homogenous case, and the definition of the total flux in Eq. (31) to simplify this contribution to
(38)  
(39) 
To get to the second line, we used equation (32) for the the flux.
We use the solution (34) and evaluate the angular integral to obtain the final equation for the effect of continuum photon transport on the ionization fraction for a planewave fluctuation.
(40) 
where the coefficients and are given in Eq. (33). The MFP of continuum photons is ; as expected the continuum photons’ contribution goes to zero when the wavelength becomes much larger than this.
Vii Radiative transfer in Lyman
This section works out the radiative transfer of photons in an inhomogenous universe. The subject and details of this calculation are selfcontained, but impact the rest of the paper through the resulting perturbed recombination rates. This sections’ results are applicable over a wide range of length scales; we show that they reduce to expected values in the large and smallscale limits in Appendices A and B.
The PSD of photons evolves via the Boltzmann equation of Eq. (23). It is simplest to work in the matter’s rest frame, since the source terms on the righthand side take on simple forms. Absorption, emission and resonant scattering contribute to this source term; we describe each of these processes in detail below.
The scattering of photons off a hydrogen atom in the ground state is a two step process, involving an excitation to a virtual excited state through the absorption of the incident photon, and subsequent decay through the emission of the outgoing one. When the first photon is of very low frequency, this corresponds to classical Rayleigh scattering. When its frequency approaches the frequency (henceforth ), the intermediate state is long lived and other processes which deplete it become important.
In particular, the excitation of the state to higher bound states and its photoionization compete with spontaneous emission. We count the former as true absorptions, and the latter as coherent scattering events. The net photon number is unaffected by coherent scattering, but the frequency of the outgoing photon is related to that of the incident one.
The branching ratio for coherent scattering is set by the rate of spontaneous emission from the state
(41) 
where is the width due to all processes, and is the complementary branching ratio for absorption via twophoton processes. Coherent scattering is the dominant process, and the scattering probability is close to unity.
A useful definition is the Sobolev optical depth of the line. It is the net optical depth for the absorption of a photon over its path as it redshifts through the line due to cosmological expansion.
(42) 
The line is optically thick at the redshift of recombination, i.e. . We divide this optical depth into true absorption and scattering contributions as
(43) 
The rate of removal of photons per unit volume of phase space due to coherent scattering is
(44) 
In the above expression, is broadened from a delta function at the frequency, , due to the thermal motions of the absorbing atoms and the finite lifetime of the excited state. The resulting profile is a Voigt function, which is most easily expressed in terms of the deviation from the central frequency in Doppler widths Hummer (1962):
(45)  
(46) 
The Voigt parameter, , quantifies the relative strength of the radiative and Doppler broadening mechanisms, and is given by
(47) 
The outgoing photon follows a redistribution function, . This is defined as the probability of an outgoing photon conditioned on the incoming photon Hummer (1962). It is normalized as
(48) 
The rate of injection per unit volume of phase space due to coherent scattering is
(49) 
True absorptions are twophoton transitions to higher states, through an intermediate ‘virtual’ state. Direct photoionization from the state is formally included by letting the summation over the higher states run over the continuum states. However, the dominant transitions from are to the and levels. To the first approximation, the resultant absorption probability is
(50) 
In the first line, we have neglected the absorption contribution in the denominator, and assumed that the PSD for the second photon of lower energy is that of a blackbody at the radiation temperature. The rate of removal of photons due to true absorption is
(51) 
In a similar manner, true emission of photons is a twophoton process, in which the first photon is emitted in a transition from one of the higher levels (as earlier, largely from and ) to a ‘virtual’ level, and the second one during a subsequent decay to the ground state. We neglect the stimulated component of both transitions since the PSDs involved are much smaller than unity. The rate of injection due to true emission is
(52) 
In principle, twophoton transitions to and from the state can also inject or remove photons within the line. Depending on the frequency of the more energetic photon involved, these are Raman scattering or twophoton transitions between and the ground state. However, these transitions are much slower than those involving the state; in particular, their crosssection goes to zero at the central frequency, since there is no phase space available for the second photon (see Fig. 6). This statement is no longer true if we include stimulated emission, but the full transition rates are still much smaller than the ones to within the line Hirata (2008). Thus, the majority of photons produced in this manner are on the far red side of the line. We can safely neglect this channel while calculating the spectral distortion within a few hundred Doppler widths of the line center.
Fig. 6 shows the rates of the radiative processes described above which add or remove photons from the frequency range of interest.
vii.1 Solution of the Boltzmann equation
We solve the Boltzmann equation under a number of simplifying assumptions.

The – transition rate is high enough so that all their sublevels are equally occupied. Consequently we neglect the fast transitions between these sublevels.

The line profile, , dominates the frequency dependence of the absorption and emission terms. Thus we replace all factors of multiplying the profile with the central frequency, .

The rates of radiative processes are large compared to the Hubble rate, so the PSD and excited level populations are effectively in steady state. This is valid within the line profile due to the high scattering rate.

The absorption and emission profiles are identical. Under this approximation, factors of are approximately equal to unity. This is valid if we restrict ourselves to frequencies which satisfy
(53) This is satisfied within the frequency range of interest, since the wings are optically thick to true absorption only up to Doppler widths at this redshift AliHaïmoud et al. (2010).

On the far blue side of the line, we take the PSD to equal that of a blackbody at the radiation temperature.

The redistribution function, , is isotropic. We condense it to the the notation .
We use the steady state approximation to balance the rate of processes which populate the level – downward transitions from higher levels and upward transitions from the level – with its net rate of depletion.
(54) 
where is the average of the phasespace density over the line profile, . We use this along with the definition of the scattering probability in Eq. (41) to rewrite the emission term of Eq. (52) as
(55) 
where we have introduced the equilibrium PSD, , which is defined as
(56) 
vii.1.1 Homogenous case
If the background ionization state and density are homogenous, the PSD is independent of direction and position. Under the assumptions listed above, the Boltzmann equation of Eq. (23) reduces to
(57) 
This is easily solved if the redistribution due to coherent scattering is unimportant, i.e., , or independent of the incoming frequency, i.e., . The PSD is then given by the Sobolev solution. Complete redistribution is a good approximation within the Doppler core (up to Doppler widths away from at AliHaïmoud et al. (2010)).
However, redistribution due to coherent scattering is nontrivial in the wings, since the average change in frequency between the incident and outgoing photons is only a few Doppler widths. We implement the resulting diffusion in frequency using a secondorder differential operator. This is commonly known as the FokkerPlanck approximation Unno (1955); Rybicki and dell’Antonio (1994); AliHaïmoud et al. (2010). It is well suited for describing the partial redistribution in the wings. Due to the high scattering rates near the line center, the PSD sets itself to the equilibrium value, and the particular prescription used becomes unimportant, as long as it yields a small result. Under this approximation, the rates of injection and removal due to scattering are
(58) 
The operator above does not account for the effect of atomic recoil; this is consistent with the approximation of equal absorption and emission profiles (assumption 4). Using this in Eq. (57), we get a second–order ordinary differential equation (ODE) for the phasespace density
(59) 
We numerically solve this differential equation in a frequency range extending out to Doppler widths on either side of , with bins per Doppler width. We set the PSD to a blackbody on the far blue side, and use a Neumann boundary condition on the far red side, where we set the derivative to zero. The latter is designed to kill an unphysical solution where the PSD grows catastrophically as we approach the red side of the line.
Technically, this region is larger than the domain of validity for some of our approximations, but we formally extend the equation out to this region in order to reduce boundary effects. We evaluate the Voigt profile using Gubner’s series in the core, and a fourth order asymptotic expansion in the wings Gubner (1994).
In order to evaluate the equilibrium PSD, , we need the occupancies of the ground () and excited () states. The rates of their depletion and population depend on the PSD itself, so to be completely selfconsistent, we need to solve for the level populations together with the PSD. Instead, we use the three level model of recombination of Section III, which assumes the Sobolev solution. The error introduced by doing so is small, because the most significant effect of the redistribution is to broaden the jump in the PSD, rather than change its amplitude.
Figure 7 shows the resulting spectral distortion, which is defined via the PSD as the number of excess photons over a blackbody distribution per hydrogen atom per logarithmic frequency interval. Also shown are the true distortion (as calculated by the publicly available HyRec code AliHaïmoud and Hirata (2011)), and the Sobolev approximation to it, which neglects redistribution due to coherent scattering. HyRec’s treatment of recombination and radiative processes is significantly more sophisticated than ours – it does not assume a steady state or equal emission and absorption profiles, follows the population of the higher levels, and accounts for twophoton and Raman transitions which are nonresonant with the transition.
The rate of recombination through the channel is the difference between the downward and upward transition rates
(60) 
We get an expression for the average monopole, , and hence the recombination rate through the Ly channel by integrating Eq. (57) over frequency, and using the normalization of the redistribution probability.
(61) 
where the notation respresents the jump in a quantity across the line, . Using this in Eq. (60), we recover the background recombination rate in the Sobolev approximation with large optical depth
(62) 
Typically the PSD on the red side, , sets itself to the equilibrium value, , due to the high optical depth. On the far blue side, we take to equal the blackbody value to maintain consistency with assumption 5 and the numerical solution.
A significant fraction of atoms reach the ground state via twophoton decays from the level. From Fig. 6, we see that the more energetic of the emitted photons is largely on the far red side of the line. The effect of absorption of the background spectral distortion in this region is largely canceled by that of the stimulated emission of the low energy photon Hirata (2008). Thus, we compute the twophoton decay rate using the blackbody PSD.