CONSTRAINTS ON PERTURBATIONS TO THE RECOMBINATION HISTORY FROM MEASUREMENTS OF THE COSMIC MICROWAVE BACKGROUND DAMPING TAIL

A primary goal of cosmic microwave background (CMB) experiments after the discovery of the acoustic features in the angular power spectrum was to measure the damping tail. This was achieved by the CBI (Sievers et al. 2003, 2007) and ACBAR (Goldstein et al. 2003; Reichardt et al. 2009) experiments, and spectacularly verified by the Atacama Cosmology Telescope (ACT; e.g., Dunkley et al. 2011) and the South Pole Telescope (SPT; e.g., Keisler et al. 2011). There is even higher precision data on its way from ACT, SPT, and Planck. Initially, the goal was simply to confirm that the damping tail agreed with theoretical predictions based on cosmological parameters determined from the first set of peaks and troughs. Contemporary high-precision data open the opportunity for new constraints of cosmological parameters that specifically influence the damping tail and have little influence on larger scales. The parameters wagged the tail, and now the tail can wag the parameters. Examples of these parameters are the primordial helium abundance Y_p and the effective number of relativistic species, N_eff. Both experiments have hinted at a deficit of power in the damping tail through the Y_p and N_eff probes (Keisler et al. 2011; Dunkley et al. 2011; Hou et al. 2011). Other physical possibilities have also been investigated to explain this damping tail tension; these include dark radiation (Smith et al. 2012; Eggers Bjaelde et al. 2012), the annihilation or decay of dark matter particles (Giesen et al. 2012), cosmic strings (Lizarraga et al. 2012), or a high-frequency cosmic gravitational wave background (Sendra & Smith 2012).

CMB anisotropies are suppressed on small scales by shear viscosity and thermal diffusion when the photons and electrons are tightly coupled, and by higher order transport effects as the photons break out from their random-walk Thomson scatterings by electrons to approach free-streaming. The physics of Silk damping (Silk 1968) has been heavily explored, as reviewed in Bond (1996), Hu & White (1997), and Sievers et al. (2003), and is very accurately computed numerically in CMB transport codes (Seljak & Zaldarriaga 1996; Lewis et al. 2000). One can associate a characteristic damping wavenumber k_D(z) with the processes, determined by the steep variation in the mean free path (n_e(z)σ_T)⁻¹ with redshift, where n_e(z) is the free-electron number density at redshift z and σ_T is the Thomson scattering cross section. The overall effect on the CMB power spectrum C_ℓ can be approximately characterized by an exponentially damped envelope suppressing the baryon acoustic oscillations, encoded in an associated angular multipole damping scale ℓ_D ∼ k_D(z_dec)χ(z_dec) ∼ 1350, where χ(z_dec) is the comoving distance to the last scattering redshift z_dec ∼ 1088.

Our focus here is how the damping tail is impacted by perturbations of the free-electron fraction, x_e(z) = n_e/(n_p + n_H i), where n_p and n_H i are the number density of ionized and neutral hydrogen atoms, during cosmological recombination in a model-independent way. This complements model-independent studies of H(z) variations (Samsing et al. 2012). Opening the constraint to the functional form of x_e(z) broadens the range of non-standard physics that can be tested, such as energy injection from decaying particles or dark matter annihilation (Chen & Kamionkowski 2004; Padmanabhan & Finkbeiner 2005), which cause a delay of recombination (Peebles et al. 2000).

The present data are only suggestive of a damping deficit, so discrimination between physical causes is not yet possible, but we develop a framework for considering general modifications and the impact of physical requirements on the posterior distributions. The recombination history is a central quantity in the interpretation of CMB data, and it has evolved significantly over time as more and more effects have been included (Zeldovich et al. 1968; Peebles 1968; Seager et al. 1999), and, with extensive recent work, it may now have nearly converged (Chluba & Thomas 2011; Ali-Haïmoud & Hirata 2011). It is important, however, to test directly whether tensions may exist which suggest that there are missing elements to the recombination story. For example, residual uncertainties in x_e(z) directly affect our ability to distinguish between different inflationary scenarios (Hu et al. 1995; Lewis et al. 2006; Rubiño-Martín et al. 2010; Shaw & Chluba 2011).

In this paper we follow the approach of Farhang et al. (2012) to generate eigenmodes for perturbations to the high-redshift ionization history focusing on the epoch of recombination. The method is applied to the seven-year Wilkinson Microwave Anisotropy Probe (WMAP) data combined with the 2011 SPT and ACT data and their recent upgrades, using a Fisher matrix eigenanalysis (Section 2) of finely binned redshift-localized x_e-modes. In Farhang et al. (2012) the eigenanalysis was performed for simulated data sets consisting of CMB maps (or their spherical harmonic amplitudes a_ℓm). Here, the modes are defined from the measured CMB bandpower errors, which leads to a different structure for the Fisher matrix. Further extending Farhang et al. (2012), we jointly treat the foreground nuisance parameters of the ACT, SPT, and WMAP7 experimental data. We use these data to constrain the amplitude of deviations from the standard recombination scenario in Section 3. In Section 4, we investigate the impact of imposing electron number conservation (through x_e(z)) on the initial and final probability distributions of the parameters and quantify its effect relative to the information delivered by the data. We conclude with a brief discussion in Section 5.

2.1. Eigenmode Construction

The functional form of deviations from the standard recombination history can be decomposed into a set of uncorrelated functions ranked by their significance in the data, here the increasing order of their forecasted errors. We construct these functions from the inverse of the Fisher information matrix F for the given combination of data used,

$$\begin{equation} {\bf F}_{ij} =\sum _{b, b^{\prime }}\frac{\partial C_{b}}{\partial q_i} {\rm Cov}_{bb^{\prime }}^{-1}\frac{\partial C_{b^{\prime }}}{\partial q_j}, \end{equation} \tag{ 1 }$$

where the C_b's are the simulated bandpower measurements for a fiducial set of cosmic parameters and for a given experiment, and ${\rm Cov}_{bb^{\prime }}$ is the bandpower covariance matrix for that experiment. The parameters q_i include recombination perturbation parameters as well as the standard cosmological parameters and the various nuisance parameters for the experiments and for secondary effects. This form for F_ij assumes uniform prior distributions for the parameters, Gaussian likelihoods for the bandpower data points, and that information on the parameters mainly comes from the mean of the bandpowers, C_b, rather than their covariance, Cov. The latter can be verified for the high multipoles of interest by comparing the contributions to the Fisher matrix from the mean and the covariance of data (see, e.g., Tegmark et al. 1997). Gaussianity is a reasonable approximation again for the high multipoles which are of interest in this work (see, e.g., Verde et al. 2003). Further, it is the assumption adopted in the released likelihood codes for SPT and ACT.

We use δln(x_e(z)) as the perturbation parameter. x_e is the parameter of direct relevance to probe the atomic physics involved at recombination. It is also straightforward to limit its variations within physical ranges. However, it does not affect the CMB anisotropies as closely as, e.g., n_e and τ do. Using n_e would also decrease possible degeneracies with baryonic matter density. Using the log in the expansion balances the low- and high-z regimes, but the data define the region of dominant impact, namely, around decoupling, rather solidly in the hydrogen recombination regime. These perturbations can be represented interchangeably by any basis that represents their degrees of freedom and does not produce numerical errors (see Farhang et al. 2012 for a discussion on other parameterizations as well as various extended and localized basis functions). Here we have used the cubic M₄ spline (see, e.g., Monaghan 2005) representation, and confirmed convergence against an increasing number of basis functions, for the explicit Fisher form, Equation (1), for these experiments.

As mentioned above, the Fisher matrix considered here includes the x_e(z) amplitudes, the instrument-dependent nuisance parameters, and the six primary cosmological parameters $(\Omega _{\rm b}h^2,\Omega _{\rm dm}h^2,\theta _{\rm s},\tau , n_{\rm s},\Delta ^2_{\mathcal{R}})$, respectively, describing the physical baryon density, physical dark matter density, the angular size of sound horizon at the last scattering surface, the reionization optical depth, the scalar spectral index, and the curvature fluctuation amplitude. In taking the full Fisher inverse and focusing on sub-blocks of it, we are in effect marginalizing over the parameter directions not in the block in the approximation of a Gaussian posterior.

To analyze multiple experiments with different sets of nuisance parameters, their individual Fisher matrices should be constructed and marginalized to give effective matrices containing only standard cosmological and x_e-perturbation parameters that are common among all experiments. The effective Fisher matrices, now with the same dimension, are added to get the total Fisher matrix $\tilde{\bf F}$. The x_e-perturbation eigenmodes are the eigenvectors of the perturbation block of ${\tilde{\bf F}}^{-1}$ with their uncertainties estimated from the roots of the corresponding eigenvalues, assuming a Gaussian distribution for the mode amplitudes in the vicinity of their maximum likelihoods. Henceforth, the x_e-eigenmodes are unambiguously referred to as the eigenmodes or simply the modes. Details on the Fisher eigenanalysis for multiple experiments with different nuisance parameters are discussed in the Appendix. To generate the bandpowers required in the Fisher matrix construction we modified the publicly available code Camb⁴ (Lewis et al. 2000) to include the more general ionization histories required for this work.

2.2. Data Sets and Their Eigenmodes

Throughout this work, the SPT (Keisler et al. 2011) and ACT (Dunkley et al. 2011) measurements of the CMB temperature and the seven-year WMAP measurement of CMB temperature and polarization (Larson et al. 2011) are used to study perturbations around the ionization history of the universe during recombination.⁵ The SPT and ACT data consist of observation of 790 and 296 deg² of the sky, at 150 and 148 GHz, during 2008–2009 and 2008 seasons, respectively. For simplicity, we neglect non-CMB cosmological constraints.

The left plot in Figure 1 shows the first two ionization eigenmodes generated for the combined SPT+WMAP7 data set. The right and middle plots show the impact of these modes (with an amplitude equal to their 1σ error bars) on the CMB temperature power spectrum and on the Thomson differential visibility g(z) = n_eσ_Tc(1 + z)⁻¹e^−τ, where the Thomson depth to redshift z is τ(z) = ∫n_eσ_T/Hdln(1 + z). The first two ACT+WMAP7 modes look very similar, so for definiteness, unless stated otherwise, modes refer to these and lower significance SPT+WMAP7 eigenmodes in this paper. Indeed, it turns out that the first two x_e-perturbation modes for Planck+ACTPol data and for a high-resolution cosmic variance limited experiment also have similar shapes (Farhang et al. 2012), with their order reversed in some cases. We now show that the physical significance of the dominant first mode is, not surprisingly, intimately related to basic perturbative features in the differential visibility.

Zoom In Zoom Out Reset image size

2.3. The Damping Physics of the Low-order Recombination Modes

Since the rank-ordered eigenmodes are direct probes of the map from x_e-trajectories to C_ℓ, the data-sensitive top-ranked modes should reflect the most basic C_ℓ-sensitive recombination effects, namely, through the damping tail, which is intimately related to the sharply peaked differential visibility. We find that the first two SPT+WMAP7 modes confirm this: they largely describe shifts in the decoupling redshift (defined as the peak of g(z), z_dec ≈ 1088) and shifts in the decoupling width, σ_z,dec, the "1σ" spread in g. We find that a +1σ amplitude for the first mode changes the visibility by ∼ − 1.4% in the width and by ∼0.4% in the position of the peak; for the second mode, the width increases by ∼6.2% and the peak by ∼0.8%.

The physical processes that define the structure of the damping tail have been well understood for a long time (for a review, see Bond 1996; Hu & White 1997 and references therein), and were discussed in relation to the experimental unveiling of the damping tail, first by CBI (Sievers et al. 2003) and then by ACBAR (Goldstein et al. 2003). Not surprisingly, the tail is controlled by the Compton scattering rate, n_eσ_Tc, and the way it runs as the baryon density n_b drops, characterized by the local power-law index, p_e = 3dln(n_e/n_b)/dln n_b =dln x_e/dln(1 + z) (Bond 1996). The basic recombination quantities can be related to p_e, which is zero at low and high z, has a maximum of about 12 and is about 9 at z_dec for ΛCDM. The peak of the differential visibility g(z) occurs when n_eσ_T/H(z) = p_e + 2, and the "Gaussian" width of decoupling in ln(1 + z)/ln(1 + z_dec) is σ_z,dec ∼ (p_e + 2)⁻¹. Thus for ΛCDM, the Compton time is about one-twentieth of the horizon size at z_dec, about one-fifth of the sound crossing time, and the relative width is about 0.09.

Earlier than decoupling, the photons and baryons are so tightly coupled by Thomson scattering that they can be treated as a single fluid with sound speed $c_{\rm s}= (1+R)^{-1/2} c/\sqrt{3}$, lowered by the extra inertia of the baryons, $R \equiv ({3 \bar{\rho }_{\rm b} / 4\bar{\rho }_\gamma })$, a photon+baryon kinematic shear viscosity (4/5)c²_s(n_eσ_Tc)⁻¹ (in a full treatment of Thomson scattering including angular anisotropy and polarization effects), zero bulk viscosity, and thermal conductivity κ_γ = n_bs_γ(n_eσ_Tc)⁻¹, where s_γ ∼ 10^9.8 is the photon entropy per baryon. In this tightly coupled regime, a WKB treatment of the perturbed photon density shows the baryon acoustic oscillations are exponentially damped, ∝exp (− ∫Γ/Hdln a), where the sound wave damping rate relative to the Hubble expansion rate is (see Section C.3.1 in Bond 1996; Kaiser 1983)

$$\begin{eqnarray} &&\Gamma /H = {1\over 2} (kc_{\rm s} /Ha)^2 {H \over n_{\rm e}\sigma _{\rm T} c} {16\over 15} \left[1+{\Gamma _{\rm diff}\over \Gamma _{\rm visc}} \right], \nonumber \\ && {\Gamma _{\rm diff}\over \Gamma _{\rm visc}} = {15R^2\over 16 (1+R)} \,. \end{eqnarray} \tag{ 2 }$$

Here a = (1 + z)⁻¹. kc_s/Ha multiplies the comoving wavenumber k of the acoustic oscillations by (approximately) the comoving distance sound travels in a Hubble time, c_s(Ha)⁻¹. The contribution of thermal diffusion relative to that of shear viscosity is Γ_diff/Γ_visc ≈ 0.22.

To relate this to a WKB damping envelope acting on C_ℓ, we replace k by ℓ/χ_dec, where χ_dec is the comoving distance from us to decoupling, and integrate up to ln a_dec. The damping scale obtained is

$$\begin{eqnarray} \ell _{\rm D} &\approx & 1.7 (p_{\rm e}+2) (1+z_{{\rm dec}})^{1/2} (c/{c}_{{\rm s, dec}}) \nonumber\\ && \times [1+\Gamma _{\rm diff}/\Gamma _{\rm visc}]^{-1/2} (1+a_{\rm eq}/a_{{\rm dec}})^{1/2} \nonumber \\ && \times ({c}_{\rm s,dec}/\bar{c}_{\rm s,dec})\sqrt{1+{1/2 \over p_{\rm e}+2}}\, . \end{eqnarray} \tag{ 3 }$$

The dominant first line gives the main dependences, (1 + z_dec)^1/2(p_e + 2)(1 + R)^1/2. With p_e,dec ≈ 9 and z_dec ≈ 1088, $c_{{\rm s,dec}}\approx 0.79 c/\sqrt{3}$ at decoupling, the first line gives ℓ_D ∼ 1360. The terms in the second and third line are subdominant. The first adds the thermal diffusion contribution to the viscous one, giving a ∼10% decrease; the second from the relativistic matter contribution to H gives a ∼14% rise. With these two, ℓ_D ∼ 1410. The fourth correction accounts for the decoupling sound speed being about 10% lower than the speed averaged over the z > z_dec range; and the last 2.4% correction occurs if we use a sharp integration down to z_dec, then stop. The third line terms change ℓ_D to 1290, but are not there if we just replace ∫Γ/Hdln a by Γ/Hσ_z,dec.

Of course a full transport treatment taking into account multipoles beyond the three that enter tight coupling (density, velocity, and anisotropic stress) is required to get an accurate damping rate. The phenomenology adopted by Hu & White (1997) estimated damping envelope functions multiplying "undamped" acoustic C_ℓ's from numerical C_ℓ-results, fitting them to a form $\exp [-(\ell /\ell _{\rm D})^{m_{\rm D}}]$, similar to the WKB approximation but with a floating slope to allow for a slower falloff reflecting complexities beyond WKB physics (such as the less severe damping associated with fuzziness of last scattering reflected in the g(z) structure cf. the stronger viscous damping; the breakout into higher temperature multipoles in the Thomson-thick to Thomson-thin transition). For ΛCDM parameters, we get ℓ_D ≈ 1345 and m_D ≈ 1.26, in better-than-expected accord with the WKB estimate.

Apart from the residual memory of the acoustic oscillations, the rise in δln C_ℓ of the first mode seen in Figure 1 (right) conforms to the $(\ell /\ell _{\rm D})^{m_{\rm D}}$ form. The fluctuations in C_ℓ are dominated by those in ℓ_D, with less sensitivity to m_D. These are related to the fluctuations in the peak position and width (which is in turn related to p_e,dec) by Equation (3),

$$\begin{eqnarray*} && \delta {C}_\ell /C_{{\rm f}\ell } \sim (\ell /\ell _{\rm fD})^{m_{\rm D}} m_{\rm D} \delta \ell _{\rm D} /\ell _{\rm fD} \\ && \delta \ell _{\rm D} /\ell _{\rm fD} \sim -{\delta \sigma _{z,{\rm dec}}\over \sigma _{{\rm f}z,{\rm dec}}} +{1\over 2} {\delta z_{{\rm dec}} \over (1+z_{\rm f,dec})} \,, \end{eqnarray*}$$

with respect to the fiducial values with subscript f. So we can interpret the first, most significant, mode as primarily due to ℓ_D variations. Similarly we can understand the sign change in Δg for the first mode as being in response to δz_dec.

If we use the same approach for the influence of helium abundance fluctuations, the effect would be the smaller number of hydrogen nuclei near decoupling, suggesting a $\delta C_\ell /C_{{\rm f}\ell } \sim - (\ell /\ell _{\rm fD})^{m_{\rm D}} m_{\rm D} \delta {Y_{\rm p}}/(1-{Y_{\rm p}})$ form, in accord with what we see in Figure 4.

To search for perturbations in the standard recombination history (Recfast; Seager et al. 1999) with recent recombination corrections included (Chluba & Thomas 2011; Ali-Haïmoud & Hirata 2011), we use the amplitude of the modes introduced in Section 2.2 as a set of new parameters, and estimate them jointly with the six primary and three nuisance parameters. The nuisance parameters follow Keisler et al. (2011) and Dunkley et al. (2011) and are the amplitudes of the Poisson and clustered power from point sources, and a template for the total thermal and kinetic Sunyaev–Zel'dovich power. The shapes of the associated C_ℓ templates do not look like our x_e-modes, and the data can differentiate what is nuisance from what may be standard recombination deviation, albeit with correlations that are fully taken into account in the statistics. This statement holds even though we restricted ourselves to single effective frequency analyses for SPT and ACT. The priors used for the nuisance parameters are taken from multi-band data particular to the flux cut for point-source removal. Unless stated otherwise, throughout the analysis, Y_p and N_eff are fixed to Y_p = 0.2478 (from SPT data for the ΛCDM model and standard BBN; Keisler et al. 2011) and N_eff = 3.046 (from the standard model of particle physics; Beringer et al. 2012).

For parameter estimation we use the publicly available code, CosmoMC,⁶ and modify it to include estimation of the mode amplitudes. We use the versions of CosmoMC adapted for SPT⁷ and ACT⁸ data set likelihood functions. We checked that when the mode amplitudes are fixed to zero, the modified CosmoMC recovers the reported SPT and ACT parameter measurements (Keisler et al. 2011; Dunkley et al. 2011). Lensing of the C_ℓ's has been included throughout this work.

Table 1 summarizes the results when one and two modes are used in the analysis, for SPT+WMAP7 and ACT+WMAP7, respectively, and compares them to the standard six-parameter model. Both experiments find a non-zero (but different) amplitude for the first mode, but they are only of 1.7σ significance and so do not correspond to a detection; they have Δχ² ∼ 2. ACT+WMAP7 also has a non-zero mean for the second mode, but only at 1.3σ. When the second mode is added to the ACT+WMAP7 analysis, the shift and uncertainty of the first mode change slightly.

Table 1. The Constraints on Cosmological Parameters with Different Sets of Parameters Used, as Measured by SPT+WMAP7 and ACT+WMAP7

Parameters	SPT+WMAP7			ACT+WMAP7
	6s	+ mode 1	+ mode 2	6s	+ mode 1	+ mode 2
100Ωbh2	2.221 ± 0.042	2.253 ± 0.046	2.249 ± 0.047	2.219 ± 0.051	2.240 ± 0.050	2.236 ± 0.053
Ωch2	0.1110 ± 0.0048	0.1123 ± 0.0049	0.1118 ± 0.0052	0.1121 ± 0.0052	0.1155 ± 0.0056	0.1121 ± 0.0061
100θs	1.041 ± 0.002	1.041 ± 0.002	1.040 ± 0.003	1.039 ± 0.002	1.039 ± 0.002	1.035 ± 0.004
τ	0.086 ± 0.015	0.089 ± 0.015	0.089 ± 0.015	0.086 ± 0.015	0.089 ± 0.015	0.0875 ± 0.015
ns	0.964 ± 0.011	0.977 ± 0.013	0.975 ± 0.016	0.963 ± 0.013	0.976 ± 0.015	0.960 ± 0.019
$10^9\Delta ^2_{\mathcal{R}}$	2.43 ± 0.10	2.40 ± 0.10	2.40 ± 0.10	2.45 ± 0.11	2.43 ± 0.11	2.45 ± 0.11
μ1a	(0)	−0.77 ± 0.46	−0.76 ± 0.47	(0)	−1.27 ± 0.74	−1.67 ± 0.86
μ2	(0)	(0)	−0.39 ± 1.09	(0)	(0)	−3.5 ± 2.7
σ8(derived)	0.807 ± 0.024	0.825 ± 0.027	0.818 ± 0.032	0.814 ± 0.028	0.841 ± 0.031	0.802 ± 0.040
δzdec/zdecb	⋅⋅⋅	−0.6%	−0.7%	⋅⋅⋅	−1.0%	−1.7%
δσz,dec/σz,decc	⋅⋅⋅	1.5%	−0.5%	⋅⋅⋅	2.6%	−14.0%
(|δxe|/xe)maxd	⋅⋅⋅	5%(z ∼ 1196)	5%(z ∼ 1039)	⋅⋅⋅	8%(z ∼ 1006)	31%(z ∼ 1076)e
Δχ2	⋅⋅⋅	2.5	2.5	⋅⋅⋅	2.1	2.5

Notes. μ₁ and μ₂ refer to the amplitudes of the first and second modes. ^aThe mode amplitudes and errors in this table (and throughout the paper) should be interpreted with respect to the normalized version of the modes as plotted in Figure 1. So, e.g., perturbations with μ₁ = 1 correspond to x_e changes in the form of mode 1 and with an amplitude exactly as plotted in Figure 1. ^bRelative change in the redshift of maximum visibility, where z_dec = 1088 is the fiducial maximum visibility point. ^cRelative change in the width of the visibility function. ^dMaximum relative change in the ionization fraction. The redshift corresponding to this maximum change is also included. ^eThis large deviation, though looking curious, is not statistically significant. This point is understandable given the relatively large estimated values for μ₁ and μ₂ and their uncertainties.

Download table as:  ASCIITypeset image

While this paper was in press, improved measurements of the damping tail were released by SPT (referred to as SPT12; Hou et al. 2012) and ACT (referred to as ACT13; Dunkley et al. 2013; Sievers et al. 2013). We present our new results using these recent datasets for the SPT11+WMAP7 modes in Table 2. We verified our SPT12 and ACT12 implementations in COSMOMC by checking that we reproduce the Hou et al. (2012) and Sievers et al. (2013) determinations for the standard six cosmological cases and also for the cases including Y_p or the Neff variation. The table shows that the mode amplitudes for SPT12 are highly consistent with those of SPT11, with a slight error decrease, with a mean for the first mode amplitude now 1.9σ away from zero. As our discussion of the similarity of the first mode δC_l shape and the perturbed Helium abundance shape illustrates, this is as expected from the new Hou et al. (2012) SPT measurement of Y_p = 0.300 ± 0.025 (and N_eff = 3.62 ± 0.48). ACT12 reported values of these parameters that were consistent with the unperturbed values of the basic-six parameter case, and as expected we see consistency with zero for the ACT12 mode amplitudes.

Our perturbed ionization history is required to satisfy electron conservation through

$$\begin{equation*} 0 \le {x_{\rm e}}(z) \le x_{\rm e}^{\rm max}, \end{equation*}$$

where x^max_e = 1 + (2Y_p/3.97(1 − Y_p)) is the maximum total electron fraction, using m_He/m_H ≈ 3.97, obtained when helium and hydrogen are fully ionized. When the mode amplitudes are poorly determined by data, the reconstructed x_e(z) could break through this bound, which of course we do not allow. Thus, although our starting prior may have been uniform with a wide possible spread in the amplitudes, the true prior distribution can only be determined with a full suite of Monte Carlo calculations restricting the allowed range. The Fisher analysis does not catch this because the amplitudes are supposed to be tiny. They are not in our case for which the data do not have strong discriminatory power so allowed variations in the ionization history can be very broad. Intuitively, if the volume spanned by the prior space is comparable to or smaller than the volume of the likelihood space (at a given significance level), the posterior will be influenced by the physical constraint, and the Fisher matrix analysis will be a poor approximation to the full analysis. Experiments with higher sensitivities will provide information about a larger number of modes before running into this condition (see simulations for Planck+ACTPol-like observations in Farhang et al. 2012).

Figure 2 shows the marginalized one-dimensional (1D) distributions of the amplitudes of the first three modes, μ₁ to μ₃, in an analysis with three modes included, and for two experimental setups: the posterior distributions of the SPT+WMAP7 case (solid black lines) and the prior-only simulations (dotted black lines) by ignoring the likelihood, i.e., assuming infinite errors in the data. The overplotted blue lines (in the left plot only) correspond to a case with only one mode being included in the analysis. Note that the prior distributions of the first and second modes are skewed toward negative values for the case with three modes in the analysis (black lines). However, the prior distribution for the first-mode-only case (the quite symmetric blue dashed curve) shows that the μ₁ measurement is not prior-driven.

Zoom In Zoom Out Reset image size

The very narrow posterior distributions of μ₁ and μ₂ relative to their priors illustrate the power of the ACT/SPT data in constraining these parameters, although they are not found to be significantly different from zero. For μ₃, on the other hand, the comparable widths of the prior and posterior distributions imply that the data set under consideration hardly contains more information about this parameter than the limits set by electron conservation. The measured errors on the fourth and higher modes differ significantly from their Fisher forecasts, not even keeping their ranking. The insensitivity of the data to higher modes explains why we have limited our study to the first two modes.

We can quantify the impact of the prior by measuring the Shannon entropy decrease in the measurement of n-parameters, q = {q_i}, associated with the transition from the prior distribution p_i to the posterior distribution p_f when data are added,

$$\begin{equation} R(n) \equiv S_{\rm i}(n)-S_{\rm f}(n) \equiv \langle {\ln p}_{\rm f} \rangle _{\rm f} - \langle {\ln p}_{\rm i} \rangle _{\rm i} . \end{equation} \tag{ 4 }$$

For us the relevant q_i's are the amplitudes of the first n modes.

Although the full calculation can be made, we have found that a Gaussian approximation works reasonably well, and does not have the numerical challenges associated with an accurate full calculation. The posterior p_f is closer to Gaussian than is the prior p_i. With an nD-Gaussian distribution with zero mean and covariance matrices C, S = 1/2ln det(C) + nln(2π)/2 + 1/2〈q^TC⁻¹q〉. The last term is one-half of the mean χ² associated with the measurement, hence is n/2 since 〈q^TC⁻¹q〉 = n. The entropy difference is then the ratio of mean log-volumes of the parameter space in question, namely,

$$\begin{equation} R(n)=\frac{1}{2}\ln ({\rm det}({\bf C}_{{\rm i}})/ {\rm det}({\bf C}_{{\rm f}})), \end{equation} \tag{ 5 }$$

where C_i,f are the prior and posterior parameter covariance matrices. We checked the Gaussian approximation by comparing Equation (5) with estimates of the integral form Equation (4) of the entropies. The integral was calculated from the nearest neighbor entropy estimate (Singh et al. 2003). This non-parametric entropy estimation method is based on the distribution of the nearest neighbor distance of the samples, here the MCMC chain outputs, and is used for parameter spaces with more than one dimension. We found that the results from Gaussian approximation agree well with those from full integration, and are less noisy as the dimensionality of the parameter space increases. Using the determinant ratio to measure the level of improvement with improved data is familiar as a figure of merit (see, e.g., Mortonson et al. 2010). Although we have found for our application for deciding which modes to include that Equation (5) is adequate, Equation (4) is the better expression for a more accurate figure of merit (Farhang et al. 2011).

As n increases, the data add less information about the parameters relative to the prior. Thus, the difference between successive R's gradually decreases. By adding new parameters, the volumes of the posterior and prior spaces change by a similar prior-dominated factor. In the limit of very large n, R(n) → constant. This is shown in Figure 3 where we compare R(n) for different data sets and various numbers of modes included, n. For this plot, the modes of each curve are the eigenmodes constructed for the corresponding experiment. For ACT+WMAP7 and SPT+WMAP7, the difference between one and two parameters is greater than the difference between other subsequent modes. This shows that these data sets are much more informative about the second mode than the higher modes, which are entering the prior-dominated regime. This difference between the first two modes and the higher ones is also evident from Figure 2. With higher precision data sets, we expect the transition from likelihood to prior-dominance to happen at a higher mode number. This prior-likelihood dominance transition hints to a natural criterion for mode-hierarchy truncation. However, one should note that choosing a quantitative mode selection criterion can be rather subjective and not necessarily applicable to all data sets.

Zoom In Zoom Out Reset image size

The mode selection criterion introduced here is much stronger than the Occam's razor argument developed in Farhang et al. (2012), where the truncation of the mode hierarchy was based on the change in the information content as more ordered modes were added. Here, we have used the posterior information that the modes have in excess of the prior, a change of perspective motivated by the weak constraints from current data sets. The analysis enters a prior-dominated regime beyond just one or two modes and there is no need to consider the more sophisticated Farhang et al. (2012) criteria.

In this work we studied how allowing for some freedom in the recombination history gives a better fit to the damping tail of the CMB power spectrum as measured by SPT and ACT, compared to the primary six-parameter model. The red solid line in Figure 4 shows the relatively small shift between the best-fit SPT+WMAP7 C_ℓ's, one with the basic six parameters fixed at their best-fit values with the x_e-perturbations on, and the other with the basic six parameters fixed at their unperturbed values. That is because the non-zero μ₁ is accompanied by compensation in the values of the standard parameters, most significantly shifting n_s and Ω_bh² (see Figure 5) to give a small net ΔC_ℓ. (Although the modes are marginalized over standard parameters, they are generally correlated with them.) When the basic six parameters are set to their x_e-unperturbed values, and only μ₁ varies, the red dashed line is obtained. Apart from the oscillation difference, the shape and value look rather like that for Y_p variation for the reasons discussed in Section 2.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The goodness of fit with μ₁ added to the analysis is comparable to cases where the recombination history is assumed to be perfectly known (see Chluba & Thomas 2011; Ali-Haïmoud & Hirata 2011 for the most recent calculations) and instead either Y_p or N_eff are allowed to vary (Keisler et al. 2011; Dunkley et al. 2011). (The corresponding best-fit C_ℓ difference between these two cases with the standard six-parameter case is not shown here but is similar to the red solid line for the μ₁ case.) As mentioned, including Y_p has a similar impact on the standard parameter measurements to that of μ₁, as illustrated in Figure 5, where n_s and Ω_bh² have the most significant shifts. N_eff, on the other hand, behaves differently, with a different C_ℓ-shape, and it also has a large impact on Ω_dmh² and H₀ (see Keisler et al. 2011; Dunkley et al. 2011).

It is also noteworthy that the measurements of σ₈ (the amplitude of linear matter fluctuations at z = 0 on scales of 8h⁻¹Mpc) for the model with N_eff for ACT data (σ₈ = 0.906 ± 0.059; Dunkley et al. 2011) and for SPT data (σ₈ = 0.859 ± 0.043; Keisler et al. 2011) are currently slightly disfavored by σ₈ inferred from clusters (σ₈(Ω_m/0.25)^0.47 = 0.813 ± 0.013 ± 0.024; here the second set of errors is systematic and due to the uncertainty in cluster masses; Vikhlinin et al. 2009). Including μ₁ and especially μ₂ in the analysis, with fixed N_eff = 3.046, brings σ₈ toward lower values (see Table 1), consistent with external data sets.

The first principal component x_e-mode, Y_p, and N_eff are strongly correlated due to their similar effects on the damping tail of the power spectrum, as illustrated in Figure 4. The correlation coefficients of μ₁ with the other two are corr(μ₁, Y_p) = 0.95 and corr(μ₁, N_eff) = 0.72, based on results from CosmoMC for SPT+WMAP7 bandpower data. Due to this partial degeneracy, varying more than one of these parameters simultaneously would significantly increase their uncertainties. The improvement in errors forecasted for near-future data and the addition of polarization information will modestly reduce the degeneracy of μ₁ with Y_p and N_eff. For example, for a Planck+ACTPol-like scenario the correlation coefficients for an ℓ-by-ℓ analysis from CosmoMC are forecasted to be corr(μ₁, Y_p) = 0.77 and corr(μ₁, N_eff) = 0.60. (An ℓ-by-ℓ analysis is more sensitive to acoustic oscillation phase information than data with wide bandpowers.)

The perturbative x_e(z)-eigenmodes are, by definition, those that are most constrained by the CMB. The perturbations in x_e(z) induced by varying Y_p and N_eff, are shown in Figure 6. They look very different than our most significant mode. What is interesting is the response in differential visibility. The first mode and the Y_p-induced perturbations look somewhat similar, although the differences are important. In C_ℓ the dominant damping tail behavior of ΔC_ℓ is even closer, although the details of the peak-trough oscillations about it differ. The full Y_p-induced perturbation involves a coherent sum over many eigenmodes, but the data are mostly trying to constrain its component in the first mode, as Figure 4 shows. The low-order mode amplitudes are therefore a way to efficiently transfer information from the CMB into constraints on ionization perturbations. Alternative high-z ionization history models due to specific physical effects (and priors they may impose) could then be differentiated in, e.g., the μ₁ − μ₂ plane.

Zoom In Zoom Out Reset image size

There are several alternative effects that can cause modifications to the CMB power spectrum in the damping tail. These include possible modifications to the physics of recombination, dark radiation (e.g., Archidiacono et al. 2011 and references therein), changes in the fine structure constant α (Kaplinghat et al. 1999; Battye et al. 2001; Scóccola et al. 2008), high-frequency cosmic gravitational wave background (Smith et al. 2006), dark matter annihilation, and particle decay (Chen & Kamionkowski 2004; Padmanabhan & Finkbeiner 2005; Zhang et al. 2006, 2007; Hütsi et al. 2009, 2011; Galli et al. 2009, 2011; Giesen et al. 2012). Deviations from standard recombination may be differentiated or corroborated by non-CMB measurements. Apart from invoking additional physical processes to explain damping tail measurements, it is important to note that the tail is sensitive to experimental issues, such as the instrumental beam, point sources in the maps, and detector time constants.

As shown in Table 2, we have now found that the damping tail tension in the SPT11 and ACT11 data of Table 2 persists in the full 2500 deg² SPT data reported in Story et al. (2012) and Hou et al. (2012), but not in the three-season ACT12 data. The Planck data, with its much larger sky coverage, should be able to address whether or not tension exists. Our forecasts for the Planck 2.5-year data show that the current SPT11 best-fit amplitude of the first mode, μ₁ = −0.77 (which is only 1.7σ now), could be detected at more than 10σ.

If such a deviation were found, it would be well beyond the levels of the standard recombination corrections which have been discussed extensively over the past several years (see Dubrovich & Grachev 2005; Chluba & Sunyaev 2006; Kholupenko & Ivanchik 2006; Switzer & Hirata 2008; Hirata 2008; Grin & Hirata 2010; Ali-Haïmoud & Hirata 2010; and additional references) and would have to indicate a new and possibly non-standard process at work. These physical mechanisms, with their high level of nuance and theoretical concordance, result in a modification to recombination which is 3.5 times smaller than the 1σ errors found here (Dunkley et al. 2011). In the far future, measurements of the cosmological recombination radiation from hydrogen and helium (see Sunyaev & Chluba 2009 for an overview) may provide another way to investigate this question and break some of the expected degeneracies.

If a detection of deviation is made with high significance, further explorations could be done with the same sort of analysis as that given here, but with modes weighted toward different redshift regimes. Indeed, although we have focused here on just the recombination epoch, viewing recombination, and reionization as a connected x_e-trajectory has appeal, since CMB data (though at large rather than small angular scales) inform the latter. In these extended studies one needs to explore other possible variables to linearly expand in rather than δln(x_e(z)), as discussed in Farhang et al. (2012). The merit of x_e(z) expansions is that one can weight them to concentrate on specific recombination regions, e.g., at higher z where helium recombines, or at lower z as x_e-freezeout is approached. Ultimately showing in a model-independent way that the allowed x_e-trajectories do not compromise our determination of cosmological parameters would further increase our confidence in conclusions drawn from CMB data sets.

This work was accomplished by support from NSERC and the Canadian Institute for Advanced Research.

The goal of this work is to search for deviations from the standard ionization scenario at high redshifts, around the epoch of recombination. For this purpose, we search for the perturbation patterns in x_e best constrained by data (see Farhang et al. 2012 for more details). The most constrained perturbation parameters are the eigenvectors of the x_e-perturbation block of the inverse of Fisher information matrix. The Fisher matrix for each data set under consideration is

$$\begin{equation*} {\bf F}_{ij}=\sum _{b,b^{\prime }}\frac{\partial C^T_{b}}{\partial q_i} {\rm Cov}_{bb^{\prime }}^{-1}\frac{\partial C_{b^{\prime }}}{\partial q_j}, \end{equation*}$$

where the q's represent any of the standard, the x_e-perturbation, the secondary, or experimental nuisance parameters. The bandpowers C_b are

$$\begin{equation*} C_b=\sum W_{b\ell } \tilde{C}_{\ell } , \tilde{C}_{\ell } = \frac{\ell (\ell +1)}{2 \pi } C_\ell , \end{equation*}$$

where the window functions, W_bℓ's, are specific to the experiment and ${\rm Cov}_{bb^{\prime }}=\langle \delta C_b \delta C_{b^{\prime }}\rangle$ is their covariance matrix. The derivatives are calculated at the fiducial q_i's. If multiple experiments are used, the total Fisher matrix is the sum of the individual Fisher matrices constructed for each experiment—if they are statistically independent—and marginalized over their nuisance, experiment-dependent parameters. To make this marginalized Fisher matrix, we divide F into blocks of nuisance parameters (represented by n) and the cosmologically interesting parameters, represented by y. We then have

$$\begin{eqnarray*} \bf {F}&= \left({\begin{array}{@{}cc@{}}\bf F_{yy} & {\bf F}_{yn}\\ \ms {\bf F}_{ny} & {\bf F}_{nn}\\ \end{array}} \right), {\bf F}|_{m}={\bf F}_{yy}-{\bf F}_{yn}{\bf F}^{-1}_{nn}{\bf F}_{ny}, \end{eqnarray*}$$

where F|_m is the Fisher matrix (including only cosmologically interesting parameters) marginalized over nuisance parameters. The individual F|_m's should be added to get the total Fisher matrix. The x_e-perturbation eigenmodes are the eigenvectors of the perturbation block of the total Fisher matrix after it has been marginalized over the standard parameters, similar to the above marginalization.

Another approach is to calculate Fisher matrices that include all nuisance parameters. Then, for those experiments that provide no constraint on a set of nuisance parameters, set those matrix entries to zero, sum the matrices over experiments and marginalize over nuisance and standard cosmic parameters.