Effective conductance method for the primordial recombination spectrum

# Effective conductance method for the primordial recombination spectrum

Yacine Ali-Haïmoud Institute for Advanced Study, Einstein Drive, Princeton, New Jersey 08540
July 7, 2019
###### Abstract

As atoms formed for the first time during primordial recombination, they emitted bound-bound and free-bound radiation leading to spectral distortions to the cosmic microwave background. These distortions might become observable in the future with high-sensitivity spectrometers, and provide a new window into physical conditions in the early universe. The standard multilevel atom method habitually used to compute the recombination spectrum is computationally expensive, impeding a detailed quantitative exploration of the information contained in spectral distortions thus far. In this work it is shown that the emissivity in optically thin allowed transitions can be factored into a computationally expensive but cosmology-independent part and a computationally cheap, cosmology-dependent part. The slow part of the computation consists in pre-computing temperature-dependent effective “conductances”, linearly relating line or continuum intensity to departures from Saha equilibrium of the lowest-order excited states ( and ), that can be seen as “voltages”. The computation of these departures from equilibrium as a function of redshift is itself very fast, thanks to the effective multilevel atom method introduced in an earlier work. With this factorization, the recurring cost of a single computation of the recombination spectrum is only a fraction of a second on a standard laptop, more than four orders of magnitude shorter than standard computations. The spectrum from helium recombination can be efficiently computed in an identical way, and a fast code computing the full primordial recombination spectrum with this method will be made publicly available soon.

## I Introduction

A significant part of our knowledge about the universe at early times and on large distance scales is derived from the observation and analysis of its spatial inhomogeneities. In particular, observations of the temperature and polarization anisotropies of the cosmic microwave background (CMB) have allowed cosmologists to determine the geometry, contents, and initial conditions of the universe to an exquisite level of precision (see for example Ref. Komatsu et al. (2011)).

Additional, and perhaps complementary information, is hidden in the tiny but unavoidable spectral distortions to the nearly perfectly thermal radiation background Sunyaev and Zeldovich (1970); Sunyaev and Chluba (2009). On the one hand, broad spectral distortions can be generated by energy injection in the early universe, taking the form of chemical potential (type) or Compton (type) distortions (and more generally continuously interpolating between these two analytic cases Khatri and Sunyaev (2012)), depending on the redshift of energy injection. Physical mechanisms causing such energy injections include the dissipation of small-scale acoustic waves, which occurs in the standard cosmological picture (see e.g. Hu et al. (1994)), and possibly the decay or annihilation of dark matter into standard model particles that then deposit their energy into the primeval plasma Hu and Silk (1993). A vast literature exists on the subject, and we refer the interested reader to the recent works Khatri and Sunyaev (2012); Chluba et al. (2012); Pajer and Zaldarriaga (2012) and references therein.

On the other hand, the primordial recombinations of helium and hydrogen lead to a few distortion photons per atom, in the form of free-bound and bound-bound photons. The seminal papers on primordial recombination by Peebles and Zeldovich et al. Peebles (1968); Zel’dovich et al. (1969) already evaluated the distortion due to Ly- transitions and two-photon decays. Even though this represents a large distortion to the Wien tail of the CMB blackbody spectrum, it lies many orders of magnitude below the cosmic infrared background (CIB, Fixsen et al. (1998)), which renders its detection very challenging, if not hopeless.

In addition, exactly one free-bound photon per hydrogen atom and two per helium atom were emitted, as well as a few bound-bound photons per atom from transitions between highly excited states111Throughout this paper we refer to bound-bound transitions between “highly excited” states as those for which the lower state itself is excited, . Chluba and Sunyaev (2006a). It was first pointed out by Dubrovich Dubrovich (1975) that these primordial recombination lines could be observed today as broadened features in the centimeter to decimeter wavelength range. A few authors have since then tackled the computation of the recombination spectrum, with various degrees of approximation or numerical convergence Bernshtein et al. (1977); Lyubarsky and Sunyaev (1983); Rybicki and dell’Antonio (1993); Burgin (2003), and it is only quite recently that highly-accurate computations were carried out by Rubiño-Martín, Chluba and Sunyaev Rubiño-Martín et al. (2006); Chluba and Sunyaev (2006a); Chluba et al. (2007); Rubiño-Martín et al. (2008).

Thus far the study of the primordial recombination spectrum has remained the niche domain of a few aficionados. First and foremost, the observational prospects may seem meager, as the predicted signal in the cm-wavelengths lies a billion times below the undistorted CMB spectrum, far below the current best upper bounds on spectral distortions from FIRAS Fixsen et al. (1996). In addition, galactic and extragalactic foregrounds would need to be understood and subtracted very precisely. Finally, the machinery required to compute a recombination spectrum with the standard multilevel method is, although conceptually not very difficult, somewhat cumbersome to implement and too computationally expensive for a systematic analysis of its information content.

The reward in finding such a needle-in-a-haystack signal is, however, potentially significant. Ref. Chluba and Sunyaev (2008a) showed that the recombination spectrum is a sensitive thermometer and baryometer. It could also provide a clean measurement of the primordial helium abundance, before the formation of the first stars Rubiño-Martín et al. (2008); Sunyaev and Chluba (2009). Finally, pre-existing spectral distortions could lead to a significant increase of the recombination radiation even if the initial distortions are small in absolute value Lyubarsky and Sunyaev (1983); Chluba and Sunyaev (2009a). The recombination spectrum could therefore be a probe of non-standard physics such as dark-matter annihilations Chluba (2010).

Let us point out, in addition, that technological advances should make it possible to reach sensitivities three orders of magnitude below that of FIRAS, corresponding to distortions at the level of (see the proposal for the PIXIE instrument Kogut et al. (2011)). Spectral distortions from recombination are only one order of magnitude weaker than this sensitivity (and even get to the level around 10 GHz Chluba and Sunyaev (2006a)), and it is not unlikely that they will be within reach of the next generation instruments. The feasibility of foreground subtraction at the required level has yet to be demonstrated, but one may hope that the spatial isotropy of the signal and its very specific spectral features should allow us to disentangle it from foreground emission.

Before any of the truly challenging issues of instrumental sensitivity and foreground subtraction are addressed, it seems that the first task is to undertake a detailed quantitative study of the information content of the recombination spectrum. In order to do so, a fast and accurate computational method is required, so that the large space of cosmological parameters can be efficiently explored222It was brought to my attention by J. Chluba that Fendt Fendt (2009) conducted a preliminary study of cosmological parameter estimation from spectral distortions, using the fast interpolation algorithm Pico Fendt and Wandelt (2007).. Introducing such a fast method is the purpose of the present work.

The main task in obtaining the spectrum resides in computing the populations of the excited states, or, more precisely, their small departures from equilibrium with one another, since the line and continuum emissions are proportional to the latter. High-precision spectra require accounting for excited states up to principal quantum number of a few hundred, resolving the angular momentum substates. With the standard multilevel atom method, one has to invert a large matrix at each timestep in order to compute the populations of the excited states (although this matrix is sparse due to selection rules, so only elements are nonzero). Currently the fastest code using this method takes about one hour for and one full day for on a standard laptop, with the computation time scaling as Chluba et al. (2010).

The method that we introduce here allows us to factorize the problem into a computationally expensive part (for which large linear systems need to be solved) that is cosmology-independent and can be pre-computed once and for all, and a computationally cheap part that does depend on cosmology. This method builds on and extends the effective multilevel atom (EMLA) method introduced in an earlier work Ali-Haïmoud and Hirata (2010) (hereafter AH10; see also Refs. Burgin (2009, 2010)); it is, however, not a trivial extension, since the EMLA method is designed for computing the free electron fraction and essentially collapses all the transitions between excited states into effective transition rates into and out of and . In this process, information not directly necessary to the evolution of the free electron fraction is lost, whereas our present goal is to go beyond and compute the full recombination spectrum.

We shall lay down our method in detail in the remainder of this paper, but the main idea can be intuitively understood if one pictures the system of radiatively connected excited levels as a circuit, where the currents are the line intensities, the voltages are the departures of the excited states populations from Saha equilibrium, and the conductances are the transition rates (this insightful analogy is due to Chris Hirata Hirata (2010)). The linearity of Kirchhoff’s laws (the steady-state rate equations for the excited state) ensures that the “current” in any transition is proportional to the outer “voltages”, i.e. the departures from Saha equilibrium of the and states. The proportionality coefficients, which we shall call effective conductances, moreover only depend on the temperature of the ambient blackbody radiation that is nearly undistorted at the relevant frequencies (this can in principle be generalized to include simple parameterizations of the ambient spectrum, as well as collisional transitions). Once the effective conductances are pre-computed as a function of transition energy and temperature, the recombination spectrum can be computed very efficiently for any cosmology, by using the EMLA method to evaluate the recombination history and the outer “voltages” as a function of redshift. This method is very similar in spirit to the widely used line-of-sight integration method for CMB anisotropies Seljak and Zaldarriaga (1996), which factorizes the computation of the CMB power-spectrum into a geometric, cosmology-independent part and a cosmology-dependent but multipole-independent source term. We illustrate our method graphically in Fig. 1.

This paper is organized as follows. In Sec. II, we write down the general equations that need to be solved for the computation of the recombination spectrum. The effective conductance method is described in Sec. III. We discuss our numerical implementation and results in Sec. IV and give our conclusions and future research directions in Sec. V.

## Ii General equations

### ii.1 Notation

We denote by the total number density of hydrogen in all its forms (ionized and neutral), the ratio of the free-electron abundance to the total hydrogen abundance, the fraction of ionized hydrogen, and (or in some cases ) the fractional abundance of neutral hydrogen in the excited state of principal quantum number and angular momentum number . The matter and radiation temperatures are denoted by and , respectively. We denote emissivities by (with units of energy per unit time per frequency interval per unit volume per unit solid angle) and specific intensities by (with units of energy per unit time per frequency interval per unit area per unit solid angle). All our derivations are for hydrogen atoms but the generalization to helium is straightforward. We only consider radiative transitions here and neglect the effect of collisions.

### ii.2 Bound-bound emission from transitions between excited states

The emissivity due to bound-bound transitions between excited states is given by

 jbb(ν)=nHhν4π× ∑2≤n

where is the frequency of the transitions and represents the radiative transition rate from to . The recombination process adds at most a few photons per atom, and the transitions between excited states are mostly below the peak of the blackbody spectrum (except for Balmer transitions, but their energy is only a few times above the blackbody peak, where there is still a very large number of thermal photons per hydrogen atom). As a consequence, the radiation field mediating the transitions is, to accuracy, a blackbody at temperature . One can check the validity of this assumption a posteriori once the distortions are computed (see for example Fig. 2 of Ref. Rubiño-Martín et al. (2006)). Note, however, that small -distortions to the blackbody spectrum (at the level of ) may significantly enhance the hydrogen and helium line emission Lyubarsky and Sunyaev (1983); Chluba and Sunyaev (2009a). We defer the study of the effect of such pre-existing distortions on the recombination spectrum to future work, and shall here assume that transitions between excited states are mediated by thermal photons only.

With this assumption, the radiative transition rates satisfy the detailed balance relations,

 Rn′l′→nl=qnlqn′l′Rnl→n′l′, (2)

where we have defined

 qnl≡(2l+1)e−En/kTr, (3)

where is the (negative) energy of the bound states with principal quantum number .

We see that if the excited states were in Boltzmann equilibrium, so that , the net emissivity would vanish. The emissivity therefore scales linearly with the small departures from equilibrium. To make this apparent, let us define the departures from Saha equilibrium at temperature with the free electron and protons,

 Δxnl≡xnl−qnlqenHxexp, (4)

where we have defined

 qe≡(2πmekTrh2)3/2, (5)

where is the reduced mass of the electron-proton system. We can now rewrite the bound-bound emissivity in the following form:

 jbb(ν)=nHhν4π× ∑2≤n

where the Saha-equilibrium pieces have cancelled out.

### ii.3 Free-bound emission

The emissivity due to free-bound transitions to excited states is given by

 jfb(ν)=nHhν4π∑n≥2∑l

where is the differential recombination coefficient per frequency interval of the emitted photon and is the differential photoionization rate per frequency interval of the ionizing photon. Recombinations of the thermal electrons and protons to the excited states are mediated by blackbody photons; as a consequence, and are related through the relation (see for example Eq. (2) of AH10):

 dαnldν = (TrTm)3/2exp[(hν+En)(1kTr−1kTm)] (8) × qnlqedβnldν.

We have purposefully made the ratio and difference of the matter and radiation temperatures appear in this expression. At high redshifts when the recombination spectrum is emitted, the matter temperature is locked to the radiation temperature by Compton heating (see for example Ref. Seager et al. (2000)). Computing the coupled evolution of the ionization history and matter temperature with HyRec333http://www.sns.ias.edu/yacine/hyrec/hyrec.html Ali-Haïmoud and Hirata (2011), we find that the fractional difference between the two temperatures is less than for and less than for . We can therefore Taylor-expand the above expression in the small parameter (note that with this convention for ) and obtain

 dαnldν ≈ qnlqedβnldν[1+(32−h(ν−νcn)kTr)ΔTT] (9) ≡ qnlqedβnldν[1+γn(ν)ΔTT],

where is the photoionization threshold from the -th shell and the second expression defines the dimensionless parameter . Again, we see that the free-bound emissivity would vanish if the excited states were in Saha equilibrium with the ionized plasma and if the matter and radiation temperature were identical. We can rewrite the free-bound emissivity in terms of small departures from equilibrium as follows:

 jfb(ν)=nHhν4π× ∑n≥2∑l

Strictly speaking, the computation of the specific intensity requires the knowledge of not only the emissivity but also the absorption coefficient Rybicki and Lightman (1986). However, the Sobolev optical depths of most bound-bound transitions are much lower than unity Hirata (2008); Grin and Hirata (2010), perhaps with the exception of the very high- transitions (), as can be seen from extrapolating Fig. 11 of Ref. Chluba et al. (2007). It is possible that a proper accounting of the nonzero optical depth in very high- transitions could lead to small modifications in the low-frequency part of the spectrum ( GHz); however, in that frequency range other effects that we are neglecting significantly affect the spectrum too, such as free-free absorption Chluba et al. (2007) and collisional transitions Chluba et al. (2010). We shall therefore assume the optically-thin limit for all bound-bound transitions between excited states, as well as free-bound transitions.

The radiative transfer equation in the optically thin regime in an expanding homogeneous universe takes the simple form

 ddt(Iνν3)≡[∂∂t(Iνν3)−Hν∂∂ν(Iνν3)]=cjνν3. (11)

In between resonances (where ), the quantity is conserved along a photon trajectory, so that

 Iν(z)=(νν′)3Iν′(z′), (12)

where

 1+z′=ν′ν(1+z). (13)

In the vicinity of a resonance line, , the solution to the radiative transfer equation is

 I−ν0=I+ν0+cJ0Hν0, (14)

where and are the specific intensities at the blue and red sides of the resonant line, respectively.

The general solution of the radiative transfer equation can be written in the following integral form

 Iν(z) = c∫∞νdν′(1+z1+z′)3jν′Hν′∣∣z′ (15) = c nH(z)∫∞νdν′Hν′jν′nH∣∣z′,

where the redshift and frequency are related through Eq. (13) and we used the fact that .

## Iii The effective conductance method

### iii.1 Populations of the excited states

In order to obtain the bound-bound and free-bound emissivities, we see that we need to evaluate the populations of the excited states, or, more precisely, their small departures from Saha equilibrium with the plasma.

The populations of excited states can be obtained to very high accuracy in the steady-state approximation, because of the large ratio of internal transition rates to the overall recombination rate. This assumption was checked explicitly in Ref. Chluba et al. (2010) and found to be extremely accurate, for the computation of both the recombination history and the recombination spectrum.

Following AH10, we separate the excited states in “interior” states, only connected radiatively to other excited states and the continuum, and “interface” states, essentially and (and potentially any additional “weak interface” states, such as …), which are radiatively connected to the ground state. We denote the populations of the former by a capital , where stands for both quantum numbers of the state, and those of the latter by , where (and … if needed).

The steady-state rate equation for the population of the interior state is

 0≈˙XK = nHxexpαK+∑ixiRi→K (16) + ∑L≠KXLRL→K−XKΓK,

where is the total recombination coefficient to the state (accounting for stimulated recombinations) and

 ΓK≡βK+∑L≠KRK→L+∑iRK→i (17)

is the total rate of transitions out of the state , where is the total photoionization rate from . Note that here again we have assumed that the Sobolev escape probability is unity in all transitions (or that the Sobolev optical depth is zero). If this were not the case the net transition rates would depend non-linearly on the state populations, which would significantly complicate matters.

Provided the transitions between excited states are mediated by blackbody photons, the transition rates satisfy the detailed balance relation Eq. (2). The recombination coefficients and photoionization rates are related through

 αK(Tm=Tr)=qKqeβK(Tr). (18)

We can now rewrite the system in terms of the small departures from Saha equilibrium with the continuum, and to linear order in :

 0 = −nHxexpΔT∂αK∂Tm∣∣Tm=Tr+∑iΔxiRi→K (19) + ∑L≠KΔXLRL→K−ΔXKΓK.

In the standard multilevel atom method, the large system (19) is solved at every timestep, which makes the computation very slow.

Following AH10, we define the matrix of elements

 MKL≡ΓKδKL−(1−δKL)RK→L. (20)

We also define the vector of elements and the source vector of elements

 ΔSK≡−nHxexpΔT∂αK∂Tm∣∣Tm=Tr+∑iΔxiRi→K. (21)

The system (19) can be rewritten in compact matrix form as

 MT(ΔX)=ΔS, (22)

where is the transpose of . We showed in AH10 that the matrix is nonsingular, and this system has the formal solution

 ΔX=(MT)−1(ΔS)=(M−1)T(ΔS), (23)

i.e., explicitly,

 ΔXK = ∑L(M−1)LKΔSL (24) = −nHxexpΔT∑L(M−1)LK∂αL∂Tm∣∣Tm=Tr + ∑iΔxi∑L(M−1)LKRi→L.

We showed in AH10 that detailed balance relations between radiative transition rates ensure that

 (M−1)LK=qKqL(M−1)KL. (25)

Using the detailed balance relation for , we may rewrite the last term of Eq. (24) as

 ∑L(M−1)LKRi→L = qKqi∑L(M−1)KLRL→i (26) ≡ qKqiPiK,

where is the probability that an excited atom initially in the interior state eventually reaches the interface state (after possibly many transitions in the interior), the complementary events being to eventually reach one of the other interface states or to be photoionized. The probabilities were defined in AH10, where they were an intermediate step to compute the effective recombination coefficients to and effective transition rates between the interface states. In the last line of Eq. (26), we have used the formal solution of the defining equation for the probabilities .

Let us now deal with the first term of Eq. (24). We define the dimensionless coefficients

 γK≡−∂logαK∂logTm∣∣Tm=Tr. (27)

Using again detailed balance relations, we obtain

 ∑L(M−1)LK∂αL∂Tm∣∣Tm=Tr = −1TrqKqe∑L(M−1)KLβLγL (28) ≡ −1TrqKqe~PeK,

where the last equality defines the dimensionless coefficient . If all the coefficients were equal to unity, then we would have , the probability that an excited atom initially in the interior state eventually gets photoionized before reaching an interface state. In general, however, (but is in general positive and of order unity), so the numbers do not have a clear physical significance but are numerically of the same order as the .

We therefore end up with the following compact expression for the departures from Saha equilibrium:

 ΔXK=qKqe~PeKnHxexpΔTT+∑iqKqiPiKΔxi. (29)

Looking more closely at equation (29), we see that the departure from Saha equilibrium of any excited state is proportional to the difference between matter and radiation temperature (so long as this difference is small) times , and to the departures from Saha equilibrium of the small set of interface states. The proportionality coefficients are functions of radiation temperature only. If we define and for interface states , we can write a general equation valid for any excited state (including the interface states) in the form

 Δxnl=qnlqe~PenlnHxexpΔTT+∑iqnlqiPinlΔxi. (30)

### iii.2 Effective conductances

We can now rewrite the net decay rate in the transitions (with ), per hydrogen atom, in the form

 ∑l,l′[qnlqn′l′Δxn′l′−Δxnl]Rnl→n′l′ =Gen′nnHxexpΔTT−∑iGin′nΔxi, (31)

where we have defined the coefficients

 Gen′n(Tr) ≡ ∑l,l′qnlqe[~Pen′l′−~Penl]Rnl→n′l′, (32) Gin′n(Tr) ≡ ∑l,l′qnlqi[Pinl−Pin′l′]Rnl→n′l′. (33)

Note the minus sign in Eq. (31): we have chosen this convention because the excited states are in general under-populated with respect to Saha equilibrium, and the coefficients defined in Eq. (33) are positive (in general the probability of reaching interface states decreases as increases).

If one thinks of the departures from Saha equilibrium as voltages and of the net decay rate in the transitions as a current, then the coefficients can be thought of as effective conductances linearly relating the two. A similar analogy can be made for the first term: there, the voltage is the fractional temperature difference and the conductance would be . Note that the have units of cm s whereas the have units of s.

Similarly, the net free-bound decay rate per unit frequency can be rewritten as

 ∑n≥2∑l

where we have defined the differential effective conductances

 dGefbdν ≡ ∑n≥2∑l

Here again we have defined the effective conductance with a positive sign, such that the current is positive when the interface states are underpopulated with respect to Saha equilibrium.

### iii.3 Bound-bound and free-bound emissivities

The expressions for the emissivities can now be rewritten in terms of the effective conductances:

 jbb(ν)=nHhν4π× ∑n

We can rewrite the total emissivity formally as

 jν≡jbb(ν)+jfb(ν)=nHhν4π× [dGedν(ν,Tr) nHxexpΔTT−∑idGidν(ν,Tr)Δxi], (39)

where

 dGe,idν≡∑n

Equation (39), along with the definitions of effective conductances given above, constitutes the fundamental result of this paper. What it means is that the emissivity can be factored into a cosmology-independent part, embodied in the effective conductances, that is computationally expensive but can be pretabulated as a function of temperature, and a simple cosmology-dependent part, entering through the departures from Saha equilibrium of and and temperature differences, which are straightforward to compute for each particular cosmology with the effective multilevel atom method developed in AH10.

Before proceeding further, let us point out that even though we have Taylor-expanded our expressions in the small difference between matter and radiation temperature, this is not required. One can very easily obtain more general expressions for arbitrary , in which the coefficient would become a function of radiation and matter temperatures. We leave it as an exercice for the interested reader to derive the exact expressions.

## Iv Implementation and results

### iv.1 Computation of effective conductances

The computation of the effective conductances requires solving large (formally infinite) linear systems. One must impose some truncation criterion to make the system finite and tractable. Following previous works, we simply ignore all excited levels above some cutoff value of the principal quantum number444Other truncation schemes could be imagined, such as assuming that the excited states above some threshold are in Saha equilibrium with the continuum; in practice, as long as a clear convergence is exhibited with increasing , we need not worry about the detailed truncation prescription. .

We tabulated the effective conductances on a grid of temperatures for several values of ranging from 60 to 500. Matrix elements for bound-bound transitions were computed as in Ref. Hey (2006) and those for free-bound transitions as in Ref. Burgess (1965). Equations (26) and (28) for the probabilities and the dimensionless numbers represent the time consuming part of the problem, as they require solving large (of order ) matrix equations. Due to selection rules for radiative transitions, these large matrices are very sparse, and only have of order nonvanishing elements. We can therefore use a sparse matrix technique identical to the one introduced in Ref. Grin and Hirata (2010).

We show some of the computed effective conductances in Fig. 2, for a temperature K (corresponding roughly to the peak of the emission).

### iv.2 Practical simplification

Equation (39) is a rather remarkable result in its raw form, but its implementation without further simplification could be somewhat cumbersome: if considering excited states up to principal quantum number , one would in principle need to tabulate of order functions as a function of temperature. For of a few hundred, required for a fully converged spectrum in the GHz region, one would need to store tens of thousands of coefficients on a fine grid of temperature values, interpolate them at each redshift, and compute the recombination spectrum on a grid fine enough to resolve all the resonances. In addition, whereas hydrogen and singly ionized helium benefit from the accidental energy degeneracy between angular momentum substates, this is not the case for neutral helium, which as a consequence has a much larger set of lines.

In order to save a significant amount of memory with negligible cost in accuracy (as we shall demonstrate in the next section), we group resonances into bins of finite width , so we use

 dGe,idν∣∣used=∑bGe,ibδ(ν−νb), (41)

where are the bin centers, and

 Ge,ib≡∑n

where is unity if falls inside the bin and zero elsewhere. The characteristic error resulting from this simplification should be of the order of a few times , the log-spacing between bins, since the recombination timescale is a few times shorter than the Hubble time555This can be understood from Eq. (13): a fractional error in the rest-frame frequency translates into the same fractional error in the redshift of emission, hence a fractional error on the emissivity, where few is the characteristic time of evolution of the populations.. In order to reduce the error induced by this simplification, we enforce that the bin centers coincide with the lowest order transition they may contain (and are logarithmically spaced otherwise). As can be seen from Fig. 2, effective conductances indeed decrease with increasing transition order (see also Table 1 of Ref. Rubiño-Martín et al. (2006)).

With this discretization, the emissivity for bound-bound and free-bound transitions between excited states becomes

 1nHjν∣∣high-nused=hν4π×∑b{Geb nHxexpΔTT −∑i=2s,2pGib Δxi}δ(ν−νb). (43)

### iv.3 Lyman-α and 2s−1s emission

The net emissivity in the Lyman- line is given by

 jν∣∣LyαnH = hν4π(x2p−3x1se−hνLyα/kTr) (44) × PescA2p,1sδ(ν−νLyα),

where is the escape probability for the optically thick Lyman- line. In the Sobolev approximation (see for example Ref. Ali-Haïmoud et al. (2010) for a detailed derivation), it is given by

 Pesc=8πHν33c3nHx1sA2p,1s. (45)

The two-photon emissivity is given by

 jν∣∣2γnH=hν4πdΛ2s,1sdν[x2s(1+fν′)(1+fν)−x1sfν′fν], (46)

where and are the values of the photon occupation number at . This expression properly accounts for stimulated decays and absorption of non-thermal photons (emitted in the Ly line for example). In order to be consistent with our simple “standard” treatment of two-photon decays (see next section), we shall neglect these two effects here and use the approximate expression

 jν∣∣2γnH≈hν4πdΛ2s,1sdν[x2s−x1se−hνLyα/kTr]. (47)

We approximated the differential two-photon decay rate with the fitting formula of Ref. Nussbaumer and Schmutz (1984).

Note that neither Eq. (44) nor Eq. (47) are accurate at the percent level. If needed, it would be relatively straightforward to include the appropriate corrections.

### iv.4 Fast part of the computation

In this section we briefly recall how the recombination history can be very efficiently computed with the EMLA method Ali-Haïmoud and Hirata (2010) and give explicit equations for the “voltages” that source the emissivities.

The EMLA formulation of the problem collapses all the fast and thermally-mediated interior transitions into the effective recombination coefficients , photoionization rates and transition rates for . All the complication of the recombination computation then resides in the slow transitions to the ground state, where a proper time-dependent radiative transfer calculation is required for percent-level precision. Here we are not worried about such subtleties, and use the following simple prescriptions for decay rates to the ground state. For the two-photon transitions from , we neglect stimulated decays and absorptions of non-thermal photons, so the net decay rate is

 ˙x1s∣∣2γ=−˙x2s∣∣2γ=Λ2s,1s[x2s−x1se−E21/kTr], (48)

where s is the spontaneous two-photon decay rate. For the net decay rate in the Lyman- line, we use the Sobolev approximation,

 ˙x1s∣∣Lyα=−˙x2p∣∣Lyα=RLyα[x2p−3x1se−E21/kTr], (49)

where is the net decay rate in the Ly- line, accounting for the small escape probability of resonant photon, given in Eq. (45).

 0≈˙x2s = nHxexpA2s+x1se−E21/kTrΛ2s,1s (50) + x2pR2p→2s−Γ2sx2s, 0≈˙x2s = nHxexpA2p+3x1se−E21/kTrRLyα (51) + x2sR2s→2p−Γ2px2p,

where the effective inverse lifetimes of the and states are given by

 Γ2s ≡ B2s+R2s→2p+Λ2s,1s, (52) Γ2p ≡ B2p+R2p→2s+RLyα. (53)

This simple 2 by 2 system (which is exact in the limit that and are the only interface states) can be solved analytically. Using detailed balance relations satisfied by the effective rates, we find that the departures from Saha equilibrium are given by

 Δx2s = s2s+s2pR2p→2sΓ2pΓ2s−R2s→2pR2p→2sΓ2p (54) Δx2p = s2p+s2sR2s→2pΓ2sΓ2p−R2p→2sR2s→2pΓ2s, (55)

where we have defined

 s2s ≡ nHxexp ΔA2s+Δx1s e−E21/kTrΛ2s,1s, (56) s2p ≡ nHxexp ΔA2p+3 Δx1s e−E21/kTrRLyα, (57)

where and is the departure of the ground state population from Saha equilibrium with the plasma, defined as in Eq. (4). Similar expressions can easily be obtained for the departures from Boltzmann equilibrium with the ground state, needed for the Lyman- and emissivities.

The rate of change of the free-electron fraction can be obtained, for example, from the departures from Saha equilibrium computed above:

 ˙xe=−∑i=2s,2p[nHxexpΔAi−ΔxiBi]. (58)

Finally, the matter temperature is obtained from the Compton-heating equation,

 ˙Tm=−2HTm+8xeσTarT4r3(1+xe+fHe)mec(Tr−Tm), (59)

where is the Thomson cross-section, is the radiation constant, is the electron mass and is the He:H abundance ratio. At high redshifts where is locked to , one can use the quasi-steady-state solution Hirata (2008)

 ΔTT≈3H(1+xe+fHe)mec8xeσTarT4r. (60)

As an aside, we point out that it is straightforward to include the effect of dark matter annihilations in this system of equations (see for example Refs. Chluba (2010); Giesen et al. (2012)). One simply has to add an additional photoionization term to the free-electron fraction evolution equation, a heating term to the matter-temperature evolution equation, and properly account for the additional excitations when solving for the populations of the excited states.

### iv.5 Numerical solution of the radiative transfer equation

We first solve for the ionization history and matter temperature as described in Section IV.4, using the code HyRec in EMLA mode, which we have adapted to also extract the populations of the excited states (more precisely, their departures form equilibrium). We store these quantities on a fine redshift grid ranging from to for future interpolation.

We follow the spectral distortion from to on a constant energy range eV, where the upper limit corresponds to the Lyman- frequency. We assume a purely thermal spectrum blueward of Ly-. Below , we set the emissivity to zero and freely redshift the distortion down to .

We discretize the emissivity on the same grid as the high- transitions emissivity. Our total discretized emissivity therefore has the form

 jν∣∣used=∑bJbδ(ν−νb). (61)

At each time step, we first redshift the pre-existing spectral distortion from one bin to the next lower one, using Eq. (12). We then update it by adding the emission from the “lines” at frequencies , as in Eq. (14). This procedure was used (with additional complications due to large optical depths and frequency diffusion) in Refs. Hirata (2008); Ali-Haïmoud and Hirata (2011). This method forces the timestep to be no greater than the smallest bin separation .

In our fiducial computation, we have grouped effective conductances in bins of width (and we recall that the central frequency assigned to each bin is chosen to coincide with the frequency of the lowest-order line it contains), and used a time-step . We have checked that reducing the bin width and timestep by a factor of 10 leads to maximum changes of at most a percent over the whole range of frequencies considered, with a root-mean-square difference of the order of 0.15%. We also checked that the spectrum is converged with respect to the redshift range over which it is computed – this stems from the fact that the emissivities are relatively sharply peaked around , as we show in the next section.

### iv.6 Results

All quantities shown in this section are evaluated for a flat universe with fiducial cosmological parameters consistent with the latest WMAP results Komatsu et al. (2011): K, , , , , .

In Fig. 3, we show the number of photons emitted per hydrogen atom per logarithmic redshift interval or per logarithmic interval of observed frequency, for the first few transitions of the Balmer series, and for the -transitions of the first few series. We see that in general the number of emitted photons decreases rapidly with the order of the transition within a series, and decreases as well (but less rapidly) for higher series. Note that some transitions may show absorption, as the H transition Rubiño-Martín et al. (2006). We find that a total of 0.63, 0.019, 0.036, 0.32 and 0.13 photons are emitted per hydrogen atom in the H, H, H, P and Br transitions, respectively.

Figure 4 shows the convergence of the high- bound-bound and free-bound spectra with . We find that the fractional difference in the total spectrum between and is less than a percent for GHz. Nothing formally limits us from going beyond with our method, since the tabulation of effective conductances needs to be done only once. However, consistently computing the spectrum at low frequencies would also require accounting for free-free absorption and collisional transitions Chluba et al. (2007, 2010), which we do not include here. We therefore limit ourselves to for now, keeping in mind that the spectrum obtained is not accurate below a few tenths of GHz.

Figure 5 shows the total recombination spectrum today, as well as its subcomponents: bound-bound, free-bound, two-photon emission from decays, and Lyman- emission. Note that in the present work we have not accounted for any of the radiative transfer effects recently investigated with the purpose of obtaining high-accuracy recombination histories (see for example Refs. Chluba and Sunyaev (2006b); Kholupenko and Ivanchik (2006); Grachev and Dubrovich (2008); Chluba and Sunyaev (2008b); Hirata (2008); Chluba and Sunyaev (2009b); Hirata and Forbes (2009); Ali-Haïmoud et al. (2010) and many more references therein). It is to be expected that these corrections will lead to a few percent corrections to the recombination spectrum; nevertheless, we are far from even detecting the recombination spectrum, and such refinements are not yet needed for this purpose. The effective conductance method is, moreover, oblivious to all the complications that may occur between interface states and the ground state, and could be very easily adapted to include these effects666A subtlety might arise when dealing with two-photon decays from higher levels, as one must avoid double-counting of the low-frequency photons.. If needed for some reason, a high-accuracy Lyman-lines and spectrum can be extracted from the modern recombination codes HyRec Ali-Haïmoud and Hirata (2011) and CosmoRec777http://www.cita.utoronto.ca/jchluba Chluba and Thomas (2011), that do account for these radiative transfer effects.

Let us point out that even though we have performed all computations with the correct matter temperature, we found that setting (and doing so consistently, including when computing the ionization history and departures from Saha equilibrium), leads to an error of at most for GHz. It would therefore be sufficient, at the percent level of accuracy, to assume