Constraining reionization with the first measurement of the CMB optical depth fluctuation - Compton-y cross-correlation
Toshiya Namikawa, Anirban Roy, Blake D. Sherwin, Nicholas Battaglia, David N. Spergel
CConstraining reionization with the first measurementof the CMB optical depth fluctuation – Compton-y cross-correlation
Toshiya Namikawa, Anirban Roy, Blake D. Sherwin,
1, 3
Nicholas Battaglia, and David N. Spergel
4, 5 Department of Applied Mathematics and Theoretical Physics,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, Unite Kingdom Department of Astronomy, Cornell University,Ithaca, NY 14853, USA Kavli Institute for Cosmology, University of Cambridge,Madingley Road, Cambridge CB3 OHA, United Kingdom Center for Computational Astrophysics, Flatiron Institute,162 Fifth Avenue, New York NY 10010, USA Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton NJ 08544, USA (Dated: February 2, 2021)We propose a new reionization probe that uses cosmic microwave background (CMB) observations:the cross-correlation between fluctuations in the CMB optical depth which probes the integratedelectron density, δτ , and the Compton y -map which probes the integrated electron pressure. Thiscross-correlation is much less contaminated than the y -map power spectrum by late-time clustercontributions; in addition, this cross-correlation can constrain the temperature of ionized bubbleswhile the optical-depth fluctuations and kinetic SZ effect can not. We measure this new observableusing a Planck y -map as well as a map of optical-depth fluctuations that we reconstruct fromPlanck CMB temperature data. We use our measurements to derive a first CMB-only upper limiton the temperature inside ionized bubbles, T b (cid:46) . × K (2 σ ). We also present future forecasts,assuming a fiducial model with characteristic reionization bubble size R b = 5 Mpc and T b = 5 × K. The signal-to-noise ratio of the fiducial cross-correlation using a signal dominated PICO-like y -map becomes (cid:39) δτ and (cid:39)
13 with CMB-HD δτ . For the fiducial model, wepredict that the CMB-HD − PICO cross-correlation should achieve an accurate measurement of thereionization parameters: T b (cid:39) +4500 − K and R b (cid:39) . +0 . − . Mpc. Since the power spectrum ofthe electron density fluctuations is constrained by the δτ auto spectrum, the temperature constraintsshould be only weakly model-dependent on the details of the electron distributions and should bestatistically representative of the temperature in ionized bubbles during reionization. This cross-correlation could, therefore, become an important observable for future CMB experiments. I. INTRODUCTION
Several hundred million years after the hot big bang,the universe was reionized. This process is poorly un-derstood and there are many open questions about thereionization epoch: when did reionization start? howlong did it last? what objects reionized the universe?was this a smooth or highly inhomogenuous process?There are many different observational constraints onthe timing and duration of the reionization epoch. Ob-servations of the Ly α absorption in high-redshift quasarspectra indicate that reionization had completed between z = 5 and 6 [1–3]. The CMB optical depth, which isinduced by the free electrons along the line of sight, hasbeen constrained by the angular power spectrum of CMBanisotropies [4, 5]. Recent measurements of the CMBtemperature and polarization spectra [4] suggest that theionization fraction at z (cid:38)
10 is less than 10%, so that anearly onset of reionization above is highly disfavored.In contrast, constraints on the properties of ionizedbubbles such as their temperature and size are still verypoor. During the reionization epoch, the ionized bub-bles in the intergalactic medium were formed by radia-tion from the first star-forming galaxies and their sizesgrew toward the end of reionization. The ionized bubbleswere heated significantly, with their temperature eventu- ally rising to ∼ K [6]. Constraints on the propertiesof the ionized bubbles can therefore provide insights intothe mechanism of cosmic reionization.An inhomogeneous reionization process produces sig-nificant spatial fluctuations of the electron density at thereionization epoch which lead to anisotropies in the op-tical depth to the CMB. The anisotropy depends on thereionization history, the characteristic radius of ionizedbubbles and their distribution [7, 8]. The anisotropies ofthe CMB optical depth have been constrained by severalanalyses [9–11]. However, the measurements so far havebeen limited by significant reconstruction noise. Assum-ing the analytic model of [7], measurements of the an-gular power spectrum of the optical depth have put anupper bound on the characteristic bubble size as R b (cid:46)
10 Mpc at 2 σ [10], while theoretical and simulation-basedstudies [12–14] suggest R b ∼ O (1 −
10) Mpc. Future ex-periments such as CMB-S4 as well as future 21-cm obser-vations such as the Square Kilometre Array (SKA) willbe able to constrain the bubble size down to ∼ −
10 [15, 16]. Progress in such ob-servations of inhomogeneous reionization will allow us to We note that the recent measurement of the global 21cm sig-nal from reionization by EDGES could, if confirmed, indicate a r X i v : . [ a s t r o - ph . C O ] F e b test our current knowledge of reionization and gain newinsights into the reionization epoch.The Sunyaev Zel’dvich (SZ) effect generated by hotelectrons during the reionization epoch can be also usedto explore the reionization epoch. In particular, the ther-mal SZ (tSZ) effect is sensitive to the temperature of theionized bubbles, which makes the tSZ signal a uniqueprobe of the gas temperature in the reionization era.However, a challenge is that the tSZ effect is dominatedby the late-time contributions from clusters [18]. Ref. [19]proposes to use the cross-correlation between tSZ andhigh-redshift galaxies for another reionization probe. Al-ternatively, the kinetic SZ (kSZ) effect generated by theradial motion of ionized bubbles during reionization couldbe detected by future CMB experiments, which will po-tentially constrain the redshift and duration of reioniza-tion [20–23]. However, the kSZ effect does not directlyconstrain the temperature of the ionized bubbles.In this work, we propose a new approach to probereionization via the cross-correlation of the anisotropiesof the CMB optical depth and the tSZ effect. The anal-ysis of this cross spectrum has the following advantages:the cross correlation can constrain the temperature of theionized bubbles, but the late-time cluster tSZ contribu-tion to the cross correlation is less significant than for atSZ measurement alone.We apply this cross correlation method to Planck dataand derive constraints on reionization from our measure-ment; we also perform forecasts for future experimentswith this method.Throughout this work we assume a flatΛCDM universe with the cosmological parameters de-fined by Planck
TT, TE, EE+lowE+lensing [5].
II. THEORETICAL MODEL
The change in temperature at each frequency ν due tothe tSZ effect is expressed as:∆ T ν T CMB = g ( ν ) y , (1)where g ( ν ) = x coth( x/ − x = hν/ ( k B T CMB ) ig-noring the relativistic corrections, T CMB is the CMB tem-perature, k B is the Boltzmann constant, h is the Planckconstant, and y is the Compton y parameter.Temperature anisotropies due to the tSZ effect havecontributions from both reionization and clusters. Hotelectron gas inside the ionized bubbles is parameterizedeither by the electron temperature, T e and the electronnumber density n e separately or by the pressure, P e = k B n e T e . The Compton y parameter is given as: y ( ˆ n ) = cσ T m e c (cid:90) d χ ak B T e ( ˆ n , χ ) n e ( ˆ n , χ ) , (2) that the true reionization process significantly deviates from thecurrent standard scenario of reionization [17]. where ˆ n is the unit vector along the line of sight, m e isthe electron mass, σ T is the Thomson scattering cross-section, c is the speed of the light in free space, a isthe scale factor and χ is the comoving distance. Theionization fraction, x e , traces the phase transition fromneutral to ionized states of the Universe, and is definedas the ratio of the free electron density to the neutralhydrogen density. In particular, if n p0 is the density ofprotons at present, we have n e = n p0 x e /a . Thus, the y parameter depends on both the electron temperature, T e , and ionization fraction, x e .The optical depth measures the integrated electrondensity along the line of sight: τ ( ˆ n ) = cσ T (cid:90) d χ an e ( ˆ n , z ) . (3)If the temperature of ionized regions in any specific di-rection ˆ n is constant during the epoch of reionization(EoR), y could be determined through the optical depthas y = ( k B T e /m e c ) × τ (cid:39) . × − ( τ / . T e / K)[24].Denoting the fluctuation part of τ as δτ , the angularcross-spectrum between δτ and y due to inhomogeneousreionization, C τyL , is given under the Limber approxima-tion as: C τyL = k B σ n m e c (cid:90) d χa χ T e ( χ ) P x e x e (cid:18) k = L + 1 / χ , χ (cid:19) , (4)and the angular power spectra of δτ can be expressed as: C ττL = σ n (cid:90) d χa χ P x e x e (cid:18) k = L + 1 / χ , χ (cid:19) . (5)In the above equations, P x e x e is the three-dimensionalpower spectrum of density-weighted ionization fractionfluctuations at comoving distance χ . This quantity cap-tures the physics and morphology of reionization, and inparticular, the size and distribution of ionized bubbles.As reionization is a complex and poorly constrained pro-cess, it is a difficult task to model the distribution of freeelectrons in different redshift bins during EoR. Therefore,we use the halo model approach to model P x e x e given by[7] in which the ionized bubbles are associated with darkmatter halos. We refer the reader to [7, 8, 25] for furtherdetails on the modelling of P x e x e and [26] for comparisonwith simulations. The model requires a size distributionof the ionized bubbles. We assume the following log-normal distribution for the bubble size distribution [8]: P ( R ) = 1 R (cid:112) πσ exp (cid:26) − [ln ( R/R b )] σ (cid:27) . (6)Here, R b is the characteristic radius of ionized bubblesand σ lnr is the width of the log-normal distribution.Note that we ignore the perturbations to the electrontemperature. The fluctuations of T e could be generatedif different bubbles/patches are reionized at slightly dif-ferent times and different bubbles have slightly differenttemperature. In this case, the two halo term is modifiedby the correlation between the temperature fluctuations.However, the one halo term dominates the angular powerspectra at all scales and such an effect should be negligi-ble. III. METHOD
The y -map can be extracted from multi-frequency tem-perature maps via the frequency dependence of the tSZeffect. We use the y -maps provided by Planck, whichare based on internal linear combination techniques suchas the Needlet Internal Linear Combination ( NILC ; [27])and the Modified Internal Linear Combination Algorithm(
MILCA ; [28]).We reconstruct the optical depth anisotropies as fol-lows. The non-zero optical depth leads to a suppressionof small-scale temperature fluctuations by e − τ . Denot-ing the temperature fluctuations in the absence of δτ asΘ, the observed temperature anisotropies become: (cid:98) Θ( ˆ n ) = e − δτ (ˆ n ) Θ( ˆ n ) = Θ( ˆ n ) − δτ ( ˆ n )Θ( ˆ n ) + O ( δτ ) . (7)The second term leads to mode coupling between thetemperature anisotropies ( (cid:96) (cid:54) = (cid:96) (cid:48) , m (cid:54) = m (cid:48) ) [8]: (cid:104) (cid:98) Θ (cid:96)m (cid:98) Θ (cid:96) (cid:48) m (cid:48) (cid:105) CMB = (cid:88) LM (cid:18) (cid:96) (cid:96) (cid:48) Lm m (cid:48) M (cid:19) f τ(cid:96)L(cid:96) (cid:48) δτ ∗ LM , (8)where (cid:104)· · ·(cid:105) CMB denotes an ensemble average where therealization of the optical depth anisotropies is held fixedand f δτ(cid:96)L(cid:96) (cid:48) is defined as: f δτ(cid:96)L(cid:96) (cid:48) = − ( C ΘΘ (cid:96) + C ΘΘ (cid:96) (cid:48) ) γ (cid:96)L(cid:96) (cid:48) (cid:18) (cid:96) L (cid:96) (cid:48) (cid:19) , (9)with C ΘΘ (cid:96) being the temperature power spectrum and γ (cid:96)L(cid:96) (cid:48) = [(2 (cid:96) + 1)(2 L + 1)(2 (cid:96) (cid:48) + 1) / π ] / . From Eq. (8),we can derive a quadratic estimator to reconstruct theoptical depth anisotropies expressed as a convolution oftwo temperature fluctuations with appropriate weights;this process is very similar to CMB lensing reconstruction[8].We first compute the following unnormalized estima-tor: δτ ∗ LM = 12 (cid:88) (cid:96)m(cid:96) (cid:48) m (cid:48) (cid:18) (cid:96) (cid:96) (cid:48) Lm m (cid:48) M (cid:19) f δτ(cid:96)L(cid:96) (cid:48) Θ (cid:96)m Θ (cid:96) (cid:48) m (cid:48) . (10)The filtered multipoles are given by Θ (cid:96)m = Q (cid:96) { C − (cid:98) Θ } (cid:96)m where C is the covariance of the observed temperatureanisotropies and Q (cid:96) is the quality factor. We follow [29]to compute Θ (cid:48) (cid:96)m ≡ { C − (cid:98) Θ } (cid:96)m in which the signal covari-ance is diagonal in harmonic space and the noise covari-ance is diagonal in pixel space with a white noise level of 30 µ K-arcmin for unmasked pixels and infinite noise formasked pixels. Denoting N (cid:96) as a noise power spectrumof the 30 µ K-arcmin white noise divided by the squareof the beam function, we choose the quality factor as Q (cid:96) = 1 / ( C ΘΘ (cid:96) + N (cid:96) ) /C Θ (cid:48) Θ (cid:48) (cid:96) [29]. . The temperaturefluctuations at 100 ≤ (cid:96) ≤ (cid:104) δτ LM (cid:105) , is nonzero dueto e.g. Galactic and point source masks, and inhomoge-neous instrumental noise. We compute it by averagingover simulation realizations. The unbiased estimator isthen given by: (cid:99) δτ LM = A τL ( δτ LM − (cid:104) δτ LM (cid:105) ) . (11)Here, the estimator normalization A τL is defined as: A δτL = (cid:34) L + 1 (cid:88) (cid:96)(cid:96) (cid:48) F (cid:96) F (cid:96) (cid:48) ( f δτ(cid:96)L(cid:96) (cid:48) ) (cid:35) − , (12)with F (cid:96) ≡ Q (cid:96) / ( C ΘΘ (cid:96) + N (cid:96) ). Using a non-zero δτ LM sim-ulation, we checked that the cross-spectrum between re-constructed and input δτ LM agrees with the input δτ power spectrum within (cid:46) (cid:99) δτ BH LM [30]. The estimator is ex-pressed as a linear combination of three estimators; theoptical depth, (cid:99) δτ LM , lensing potential, (cid:98) φ LM , and squareof point-source fluctuations, (cid:98) σ LM , so that (cid:99) δτ BH LM is im-mune to the bias from lensing and point-source fields(see [30–33] for the details of (cid:98) φ LM and (cid:98) σ LM ). The bias-hardened estimator is given by [34]: (cid:99) δτ BH LM = (cid:88) x = δτ,φ,σ { R − L } δτ,x (cid:98) x LM . (13)The elements of the above matrix are defined as: { R L } a,b = A aL L + 1 (cid:88) (cid:96)(cid:96) (cid:48) F (cid:96) F (cid:96) (cid:48) f a(cid:96)L(cid:96) (cid:48) f b(cid:96)L(cid:96) (cid:48) , (14)where a, b are either δτ , φ or σ , and f φ(cid:96)L(cid:96) (cid:48) and f σ(cid:96)L(cid:96) (cid:48) are thefunctions characterising the mode mixing in temperatureanisotropies by lensing and point-source fields, respec-tively [32, 33]. The estimator normalization, A φL and A σL ,are given by Eq. (12) but replacing f δτ(cid:96)L(cid:96) (cid:48) with f φ(cid:96)L(cid:96) (cid:48) and f σ(cid:96)L(cid:96) (cid:48) , respectively. A bias-hardened estimator which onlymitigates the lensing contribution can be obtained by set-ting { R L } δτ,σ = { R L } φ,σ = { R L } σ,δτ = { R L } σ,φ = 0 For details, see Appendix A.1 of [29]. and (cid:98) σ LM = 0 in Eq. (13). The bias-hardened estima-tor minimizes any contamination by lensing and pointsources; minimizing lensing contamination is especiallycritical because the correlation between the y -map andlensing is much larger than C τyL . Note that at the resolu-tion of Planck, the pressure profile of all but the largestand most nearby clusters is not resolved, so that theyappear as point sources and can be accounted for by thisprocedure. The bias-hardened estimator has a larger re-construction noise than the non bias-hardened estimator.In our case, using the bias-hardened estimator increasesthe reconstruction noise by ∼ −
40% depending on scale.
IV. DATA AND SIMULATION
We use CMB data from the Planck 2015 release in ouranalysis . The procedure for measuring δτ LM is simi-lar to that described in [10] and the Planck 2015 lensingpaper [29]. We use temperature anisotropies from thePlanck component separated maps ( SMICA or NILC ) inorder to reconstruct δτ LM . Following [29], we use tem-perature multipoles between (cid:96) = 100 − W CMB ( ˆ n ), used in thelensing trispectrum analysis to avoid Galactic foregroundand point sources, and make use of the simulations ofCMB signal and noise realization publicly available on NERSC .In addition, we use the Planck NILC or MILCA y -mapfor our cross-correlation measurement. In our baselineanalysis, we use the full NILC y -map and apply the maskobtained by combining the Galactic mask that retains60% of the sky and point-source mask provided by Planck[35]. To verify the stability of our results, we test thedependence of our analysis on the component separationmethod ( NILC / MILCA ) and the masking of the y -map. Toevaluate the errors of the measured cross spectrum, wegenerate Gaussian random field simulations of the y -mapwhose power spectrum is consistent with the measuredspectrum. The expected correlation with the δτ mapis sufficiently small that we may safely neglect it in oursimulations.In our cross-spectrum measurement, we multiply the y -map by an apodized mask, W y ( ˆ n ), calculate the spheri-cal harmonic coefficients, y LM , and deconvolve them witha 10 arcmin Gaussian beam. We finally cross-correlate (cid:99) δτ BH LM and y LM to obtain C τyL . The cross-spectrum is di-vided by (cid:82) d ˆ n [ W CMB ( ˆ n )] W y ( ˆ n ) to correct the normal-ization for masking effects. We also use the auto powerspectrum of the optical depth anisotropies to constrainreionization. The methodology for measuring C ττL from http://pla.esac.esa.int/pla/ https://crd.lbl.gov/departments/computational-science/c3/c3-research/cosmic-microwave-background/cmb-data-at-nersc/ the reconstructed optical depth anisotropies is the sameas that of [10]. V. RESULTS
Fig. 1 contains one of the main results of our analy-sis: the measurement of C τyL . We show three differentcross-spectra in this plot, which differ in how they recon-struct δτ : δτ is reconstructed either with the standardquadratic estimator, with a bias-hardened estimator tomitigate lensing, or with a bias-hardened estimator thatmitigates both lensing and point-source contributions.It can be seen that the lensing contamination of (cid:99) δτ LM produces a significant bias in the cross spectrum mea-surement since the lensing and y -map are well correlated[36]. There is almost no difference in the cross spectrumfrom the lensing-hardened and lensing-source-hardenedestimators, indicating that the point-source-like extra-galactic foreground contribution (e.g., from tSZ and CIBleaking through the δτ estimator and correlating withthe y -map) is negligible in the cross-spectrum measure-ment. We compute the chi-squared probability-to-exceed( χ -PTE) for the baseline spectrum with respect to null,finding that the value is 0 .
18 and the spectrum is consis-tent with null.Fig. 2 shows the measured cross spectra for differentanalysis choices. Using the different component sepa-ration methods for the y -map and CMB temperatureanisotropies, we can test for contamination by extra-galactic foregrounds such as tSZ, CIB and point sources.Using the different Galactic masks, we can also test forGalactic foreground contaminants. Overall, we only findsmall shifts to the bandpowers obtained when we varythe analysis choices. We compute the χ -PTE for eachcase, finding that the values are in the range between0 .
21 and 0 .
74 and that all these measured cross spectraare consistent with null within 2 σ . VI. DATA INTERPRETATION
We assume that both reionization at high redshift andclusters at low redshift contribute to the C τyL signal. Thetotal signal can be expressed as: C τyL = C τy, re L + α × C τy, lowz L , (15)where the first term is equivalent to Eq. (4) and C τy, lowz L is the contribution from clusters calculated from cosmo-logical simulations [37] and diffuse sources and α is a freeparameter which takes into account the uncertainties inthe estimation.We perform Markov chain Monte Carlo (MCMC) anal-ysis using emcee software package [38] to constrain thereionization parameters such as the average temperatureof ionized bubbles T b , characteristic radius of ionizedbubbles R b , and their distribution σ lnr [7]. We also
250 500 750 1000 1250 1500 1750 2000 L − . − . . . . × L C τ y L baseline (lens-hardening)no bias-hardeninglens & source-hardening FIG. 1. Optical depth – Compton-y cross-spectrum mea-surements with different estimators of the CMB optical depthfluctuations δτ ; the bias-hardened estimator to mitigate lens-ing bias (red, baseline) and the bias-hardened estimator tomitigate both lensing and point source bias (green, lens &source hardening). A signal is not significantly detected butthis measurement can be used to place constraints on reion-ization parameters. Note that, to illustrate the importanceof bias-hardening, we also plot the result of naively using asimple quadratic estimator for the optical depth (grey, nobias-hardening).
250 500 750 1000 1250 1500 1750 2000 L − . − . . . . × L C τ y L BaselineMILCA y-mapG50 mask for y-mapNILC CMB mapNo SZ
FIG. 2. Cross-spectrum measurements as in Fig. 1 butwith variations of the standard analysis choices, i.e., the useof a MILCA y-map, an aggressive Galactic mask (G50 mask)for the y-map, a NILC component separated CMB map, or anSZ-deprojected (No SZ) CMB map. The stability of the resultto variations in these analysis choices suggests that any con-tamination of the measurement by galactic and extragalacticforegrounds is small. marginalize over the normalization factor of the low- z contribution, α . Our model of C τyL depends on theseparameters as well as on the reionization history andthe bias of ionizing sources, b . For simplicity, we keep b = 6 fixed and check the robustness of our results against b below. We place the following flat logarithmic pri-ors on 0 .
01 Mpc < R b <
10 Mpc, 0 . < σ lnr < . < log( T b ) <
7, and 0 . < α < . . . . . . z . . . . . . . x e Kulkarni19 z re = 6 . z re = 7 . z re = 7 . FIG. 3. Evolution of the ionization fraction of the Uni-verse inferred from the dynamic range radiative transfer sim-ulation [3]. For comparison, we show the tanh reionizationmodel varying the redshift of reionization, z re . Error barsare shown from the study of the spectral shape of the twohighest-redshift Quasars [39, 40] and Ly α and Ly β forests [1]. ionization fraction. In Fig. 3, we show the reionizationhistory from the Sherwood simulation [3] which givesrise to ¯ τ = 0 . tanh models of ionization histories parameterized by the red-shift of reionization z re and the duration of reionization∆ z as x e ( z ) = ( f / y re − y ) / ∆ y re )], where y ( z ) = (1 + z ) . , y re = y ( z re ) and ∆ y re = √ z re ∆ z and f = 1 .
08 for the first ionization of helium atoms.For a fixed ∆ z = 1 .
5, we varied the redshift of reion-ization z re to 6 .
6, 7 . .
9; the isotropic CMB op-tical depths that correspond to these reionization his-tories are ¯ τ = 0 . , .
054 and 0 .
061 respectively. Un-less mentioned specifically, we have used the Kulkarni19[3] reionization history for C τyL calculation from semi-analytic model.In Fig. 4, we compare the C τyL signal from the reioniza-tion and cluster contributions. We use R b = 10 Mpc and σ lnr = ln(2) for this C τyL calculation. The cluster contri-bution in C τyL is estimated using cosmological simulations[37]. Here, δτ and y maps were made at correspond-ing redshift slices and cross-correlated. The cross powerspectra were then combined following the procedure in[42]. The reionization contribution becomes larger atlarge scales than that of cluster whereas the cluster con-tribution becomes dominant on small scales. For fixed R b and σ lnr , the ionized bubble temperature scales theamplitude of the signal linearly keeping the shape of thespectra unchanged. α = 0 . C τyL is ten times smaller than thatestimated from these simulations.In Fig. 5, we show the constraints on R b and σ lnr through a triangle plot showing 1 σ and 2 σ contours. Weperformed the analysis by fitting our reionization modelwith C τyL alone and also for a joint analysis of C τyL and C ττL . We quote the 2 σ upper limit on these parameters inTable I. We impose a hard prior 0 .
01 Mpc < R b <
10 Mpc
250 500 750 1000 1250 1500 1750 2000 L L ( L + ) C y L / Reionization, T b =10 K Reionization, T b =10 K Low z, =1Low z, =0.1
FIG. 4. The cross-correlation signal between δτ and y arisingfrom reionization (blue) and clusters (red). We assume a char-acteristic ionized bubble radius R b = 10 Mpc and distributionwidth σ lnr = ln(2), and we vary the ionization temperature T b and α (a parameter rescaling the low-redshift cluster con-tribution to the cross-spectrum). The cluster contribution iscomputed from [41]. since the simulation studies suggest R b (cid:46)
10 Mpc. Al-though the 2 σ upper bound on R b does not significantlychange for the joint analysis with this prior, the tem-perature bound decreases by 45% for the joint anal-ysis relative to the C τyL only analysis. To check thevalidity of our results we also assume a flat prior on0 .
01 Mpc < R b <
50 Mpc. We show dependencies of ourresult on the choice of priors on R b in Fig. 6. The con-straint on T b becomes ∼
10% tighter if we assume a prior0 .
01 Mpc < R b <
50 Mpc. However, the upper limit onthe temperature of the ionized bubbles is more than oneorder of magnitude larger than the measurement fromthe Ly α forest [43].We quote the 2 σ upper limit on these parameters in Ta-ble I. The upper limits on R b and T b increase by ∼
30 %and ∼
10 %, respectively, if z re varies from 7 . . R b - σ lnr relation.For different choices of Gaussian prior on b we did notfind any significant changes on the 2 σ upper limit on R b and T b . Similarly, we use different reionization historiesas shown in Fig. 3 and find that the the 2 σ bound on R b and T b changes by, at most, 10% and 2% respectively.In our analysis, we consider only the temperature in-crease of the IGM due to the ionization of hydrogen andsingly ionized helium. It should be noted that the IGMtemperature is also heated up by ∼ − z > . − R b l o g ( T b ) l n r lnr log( T b ) yy + FIG. 5. 1 σ and 2 σ constraints on reionization parametersarising from our Planck measurements of C τyL shown in Fig. 1(as well as C ττL ). We show constraints on the temperature T b and radius R b of ionized bubbles and their distribution alongwith the normalization factor of the cluster contribution α .Priors on R b [0 . ,
10] [0 . , C τyL C τyL + C ττL C τyL R b . . . σ lnr .
81 0 .
83 0 . T b ) 6 .
11 5 .
85 6 . α . . . σ upper bound on parameters de-scribing reionization derived (in part) from our baseline C τyL measurement from Planck data shown in Fig. 1. this effect would further tighten our constraints on theIGM temperature. We also do not include the contribu-tion of the cooling of different metal lines present in theIGM, and the cooling due to the adiabatic expansion andinverse Compton scattering off of the CMB. The temper-ature of the IGM depends on the shape of the radiationspectrum of ionizing sources; one part of it ionizes theneutral hydrogen atoms and rest of it increases the ki-netic energy of the electron gas. In reality, this process isvery complex as the shape of the spectral energy distri-bution varies from one source to another source. It alsoinvolves the modelling of adiabatic cooling due to expan-sion, the temperature fluctuation in the IGM due to theinhomogeneous nature of reionization, and the influenceof low mass galaxies in the reionization process.
10 20 30 40 R b l o g ( T b ) l n r lnr log( T b ) R b <50Mpc0.1Mpc< R b <10Mpc FIG. 6. Constraints as in Fig. 5 ( C τyL -derived constraints),but we explore the sensitivity to priors on the characteristicbubble radius R b by imposing two different flat logarithmicpriors on 0 .
01 Mpc < R b <
10 Mpc and 0 .
01 Mpc < R b <
50 Mpc. Most parameters appear insensitive to the details ofthis prior choice.CMB experiment f sky θ [arcmin] σ P [ µ K-arcmin] α del S4-Wide 0 . . . . .
04 2 . .
42 0 . . . . . δτ . f sky , θ , σ P and α del denote the observed skyfraction, angular resolution, white noise level in the polariza-tion map, and residual fraction of lensing B -mode spectrumafter delensing, respectively. VII. FORECASTS
In this section we investigate how the constraintson reionization parameters, such as the temperatureand ionized bubble size, will improve for future high-resolution and low-noise CMB experiments.In our forecasts, we define a fiducial model for C τyL and C ττL with the parameters R b = 5 Mpc, σ lnr = ln(2), T b = 5 × K, and α = 1. For the y -map, we considerthe Planck and PICO [45] cases. We use the observed C yyL for the Planck case. In the PICO case, we assume anoiseless y -map spectrum since the PICO y -map is signaldominated at the scales of our interest [45]. For recon-structing δτ , we consider the CMB experiments summa-rized in Table II; “S4-Wide” corresponds to the CMB-S4wide area legacy survey. “S4-Deep” denotes the CMB-S4delensing survey, i.e. a sky fraction of few percent butwith very deep observations to effectively perform delens-ing. The “HD” survey is similar to CMB-HD [46]. We do not consider the PICO δτ since the expected noise level ishigher than S4 due to the limited angular resolution. Forsimplicity, we adopt approximate noise levels and angularresolutions which mimic the above experiments. For thepurpose of forecasts, we simplify the covariance of the C ττL and C τyL spectra as being described by a diagonalmatrix, i.e., σ ( C τyb , C τyb (cid:48) ) = δ bb (cid:48) (cid:34)(cid:88) L ∈ b f sky (2 L + 1) A δτL C yy, obs L (cid:35) − , (16) σ ( C ττb , C ττb (cid:48) ) = δ bb (cid:48) (cid:34)(cid:88) L ∈ b f sky (2 L + 1)2( A δτL ) (cid:35) − . (17)We ignore C τyL and C ττL in σ ( C τyb , C τyb (cid:48) ) because C τyL (cid:46) . × ( A δτL C yyL ) / and C ττL (cid:46) . × A δτL , respectively,for our fiducial model. The cross-covariance between C ττL and C τyL is set to zero. The reconstruction noise is com-puted from the EB -estimator with CMB multipoles at100 ≤ (cid:96) ≤ B -modes (delensing), the lensing B -modespectrum is scaled by α del in Table II whose values aredetermined by performing the iterative delensing proce-dure proposed by [47]. In this calculation, we ignore thefact that in practice we need a joint iterative estimate of δτ and lensing from a given CMB map [48]. Althoughthe joint estimate could enhance the noise level of δτ , theimpact is expected to be small; the two effects, optical-depth anisotropies and lensing, on CMB anisotropies arevery different and the degradation of the noise level whenbias-hardening against lensing is typically 10% in the caseof the quadratic estimator [30]. A thorough analysis ofthis joint estimate using optimal reconstruction methodsis beyond the scope of our paper and will be addressedelsewhere.Fig. 7 shows the cross-power spectrum, C τyL , with theexpected error bars computed from Eq. (16). Futureexperiments have much better sensitivity to detect thecross spectrum due to both improving the reconstructionnoise level of δτ with lower polarization map noise andthe y -map with better component separation from mul-tiple frequency data. Combining δτ from CMB-S4 anda y -map from PICO, and assuming our fiducial model,the signal-to-noise of C τyL becomes (cid:39)
7. Note that thesignal-to-noise linearly scales with T b . Using δτ fromCMB-HD and y -map from PICO further improves thesignal-to-noise of C τyL to be (cid:39)
13. In Fig. 8 we forecastthe constraints on the reionization model parameters forfuture CMB experiments. In this analysis, we assume aGaussian prior on α = 1 ± .
1. Current measurementsof density and pressure profiles from kSZ and tSZ cross-correlations, respectively, with luminous red galaxies arealready observed at ∼ σ or better [49, 50], so it reason-able to assume that future measurements will improveour understanding of the low-z contribution to C τyL at the10% level. In Table III we quote the reionization parame-ter constraint forecasts. We find that for current y mapsfrom Planck, the constraints on reionization parameters
250 500 750 1000 1250 1500 1750 2000 L − − L ( L + ) C τ y L / π PICO x S4-WidePICO x S4-DeepPICO x HD
FIG. 7. Forecast cross-power spectrum between optical-depthfluctuations and a y -map, C τyL , with expected error barsshown for some of the experimental cases we consider. We as-sume R b = 5 Mpc, σ lnr = ln(2), and T b = 5 × K. It can beseen that for future CMB experiments, our forecasts indicatethat significant detections of C τyL appear possible (assumingour fiducial reionization model). are fairly weak, unless a futuristic δτ measurement fromCMB-HD is provided. However, if a signal-dominated y -map, e.g. from PICO, is available, combining the opticaldepth anisotropies measured from CMB-S4 could tightlyconstrain reionization parameters, allowing a detectionof the temperature of 5 × K at more than 3 σ .Note that we have ignored any possible contamina-tion in the cross-spectrum measurements in the fore-cast. For example, as we discussed above, leakage fromCIB and kSZ propagating through the quadratic δτ es-timator could bias the cross-spectrum. Although thisbias should be negligible in our measurement with cur-rent data, in future CMB experiments, we will use CMBanisotropies up to very small scales, which could lead toa non-negligible bias in the measurement. Point-sourcehardening can capture some amount of this bias, but itis not clear how well point-source hardening will workfor higher-precision experiments. However, even if sim-ple point-source hardening does not completely subtractthe bias from CIB and kSZ, we can construct a bias-hardened estimator with source profiles as in [51], or em-ploy more sophisticated multifrequency mitigation tech-niques. Further, in future experiments, the signal-to-noise of δτ is almost determined by precision of polariza-tion data. The polarization-based reconstruction of δτ should be much less affected by the above contaminants.Therefore, astrophysical contaminants are not expectedto significantly bias C τyL measurements in future experi-ments. VIII. DISCUSSION
The temperature measurement of the IGM carries im-portant information about the source of the reionization.
Experiments R b [Mpc] σ lnr T b [K]S4-Deep × Planck 9 . +2 . − . . +0 . − . +23000 − S4-Wide × Planck 9 . +2 . − . . +0 . − . +13000 − HD × Planck 5 . +0 . − . . +0 . − . +14000 − S4-Deep × PICO 7 . +1 . − . . +0 . − . +12000 − S4-Wide × PICO 7 . +1 . − . . +0 . − . +9500 − HD × PICO 5 . +0 . − . . . − . +4500 − TABLE III. 1 σ uncertainties on the reionization parametersforecast for future CMB experiments. We here consider re-constructing the fluctuations in optical depth for configura-tion of the first experiment and then cross-correlating withthe Compton y parameter for the second experiment. Various ionizing sources emit photons with different ener-gies which later heat up the surrounding neutral mediumby the propagation of ionization front. The temperatureof the IGM depends on the velocity of the ionization frontand the energy of the ionizing photons that determinesthe volume average ionization fraction at a particularredshift [52]. Radiative transfer simulations show thattemperature of IGM can reach 25000 − km/s [53–55].Observation of the Ly α forest around Quasar prox-imity zones is an alternative and powerful probe of theIGM temperature. However, it does not necessarily rep-resent the global properties of the IGM temperature athigh redshift as there are very few Quasars at z > α forest power spectrum, which means thatthe temperature is, to some extent, model dependent.Recent measurements of temperature from the high res-olution Ly α forest spectra show that the temperature is ∼ ± z ∼ .
8. The first constraint ofthe IGM temperature at z = 6 is 23600 +9200 − K [57].In the future, the precise measurement of the tempera-ture fluctuations of individual HII bubbles will enrich ourunderstanding about the metallicity of the ionizing stars.The temperature of HII regions ionized by the interme-diate metallicity stars is 2-3 times higher than the HIIbubbles ionized by the low metallcity stars [58]. On theother hand, to keep the gas ionized during the EoR, thevirial temperature of the halos for star formation shouldbe higher than 10 K. At z > M (cid:12) are ∼ C τyL .
10 20 30 40 R T b / l n r lnr T b /10 CMB-S4D × PlanckCMB-S4W × PlanckCMBHD × PlanckCMB-S4D × PICOCMB-S4W × PICOCMBHD × PICO
FIG. 8. Forecast of reionization parameter constraints ex-pected from future experiments (the experiment named first isthe experiment providing the δτ map, the experiment namedsecond provides the y -map). We show 68% and 95% contoursderived from a joint analysis of C τyL and C ττL . It can be seenthat precise constraints appear possible, especially when asignal-dominated PICO-like y map is available. IX. CONCLUSIONS
In this paper, we present the first constraint on theoptical depth – Compton-y cross-correlation C τyL fromPlanck CMB data and use it to place bounds on severalimportant properties of reionization: the temperature ofionized bubbles T b as well as the bubbles’ characteristicradius R b and size distribution.We also present forecasts of how well this new observ-able can constrain reionization parameters with futureexperiments. Our forecasts show that future experimentswill be able to measure C τyL and C ττL signals and thusreionization parameters such as T b and R b with muchbetter sensitivity. In particular, for a fiducial reioniza-tion model, we predict that the CMB-HD – PICO cross-correlation should enable a precise measurement of thereionization parameters: T b (cid:39) +4500 − K and R b (cid:39) . +0 . − . Mpc; with a CMB-S4 – PICO cross-correlation,informative constraints on reionization paramters (thatare only ∼ z <
2; we defer a detailed discussion of this possibilityto future work.The unique aspect of the C τyL signal is that it can beused as a direct probe of the reionization temperature.In a joint anaysis with C ττL , the details of the reionization morphology do not affect the temperature constraint sub-stantially since the model of the reionization morphologyis constrained significantly by C ττL .Future measurements of this C τyL signal will hence com-plement other CMB constraints on reionization. Thusfar, CMB data have constrained the duration and redshiftof reionization via the sky-averaged optical depth [60] andthe upper limits on the kSZ power spectrum [4, 61]. Inthe future, high resolution, low noise observations fromground-based, CMB experiments such as AdvACT [62],SPT-3G [63], SO [64], and CMB-S4 [65] will provide newinsights into reionization. Recent studies by [11, 66] pro-pose using cross-correlation between δτ and large-scalestructure tracers, showing that polarization data fromfuturistic experiments will be useful to probe reioniza-tion properties. In addition, LiteBIRD [67] will achievea cosmic-variance limited measurement of large-scale E modes, providing a constraint on the isotropic part of theoptical depth, ¯ τ , with σ (¯ τ ) (cid:39) × − . Looking ahead,combined analyses of C τyL , the kSZ power-spectrum andtrispectum [23], and the large-scale polarization can beexpected to significantly improve our knowledge of reion-ization from CMB observations alone. These CMB ob-servations will be independent of and complementary tonew insights gained from 21cm experiments. ACKNOWLEDGMENTS
We thank Simone Ferraro, Colin Hill, Vid Irˇsiˇc, GirishKulkarni, and Daan Meerburg for providing helpful andconstructive comments on the draft. TN and BDS ac-knowledge support from an Isaac Newton Trust EarlyCareer Grant and from the European Research Coun-cil (ERC) under the European Unions Horizon 2020 re-search and innovation programme (Grant agreement No.851274). BDS further acknowledges support from anSTFC Ernest Rutherford Fellowship. NB acknowledgessupport from NSF grant AST-1910021. For numericalcalculations, we used resources of the National EnergyResearch Scientific Computing Center (NERSC), a U.S.Department of Energy Office of Science User Facility op-erated under Contract No. DE-AC02-05CH11231. TheFlatiron Institute is supported by the Simons Founda-tion.
Appendix A: Bias-hardened estimator for CMBoptical depth
Here, we derive the bias-hardened estimator which isused for the mitigation strategy of the lensing and point-source contributions in estimating δτ .We first note that the lensing produces similar modecouplings in the temperature anisotropies. The quadraticestimator for the optical depth anisotropies given aboveis then the same as that for the CMB lensing potential, φ , but with a different weight function. In the case of0CMB lensing, the mode coupling of Eq. (8) is given by[32]: (cid:104) (cid:98) Θ (cid:96)m (cid:98) Θ (cid:96) (cid:48) m (cid:48) (cid:105) CMB = (cid:88) LM (cid:32) (cid:96) (cid:96) (cid:48) Lm m (cid:48) M (cid:33) f φ(cid:96)L(cid:96) (cid:48) φ ∗ LM , (A1)where f φ(cid:96)L(cid:96) (cid:48) is defined as f φ(cid:96)L(cid:96) (cid:48) = − γ (cid:96)L(cid:96) (cid:48) (cid:112) L ( L + 1) (cid:20)(cid:112) (cid:96) (cid:48) ( (cid:96) (cid:48) + 1) (cid:32) (cid:96) L (cid:96) (cid:48) − (cid:33) C ΘΘ (cid:96) (cid:48) + (cid:112) (cid:96) ( (cid:96) + 1) (cid:32) (cid:96) (cid:48) L (cid:96) − (cid:33) C ΘΘ (cid:96) (cid:21) . (A2)The estimator for the CMB lensing potential is then ob-tained in the same form as Eq. (10) but replacing f δτ(cid:96)L(cid:96) (cid:48) with f φ(cid:96)L(cid:96) (cid:48) . The estimator normalization is also obtainedin the same way.The mode mixing is also induced by point sourcesand we can estimate the point-source fields using thequadratic estimator. The temperature fluctuations havethe following additive term in the presence of the pointsource, and more generally any additive fields, s ( ˆ n ), in-cluding inhomogeneous noise: (cid:98) Θ( ˆ n ) = Θ( ˆ n ) + s ( ˆ n ) . (A3)Assuming (cid:104) s ( ˆ n ) s ( ˆ n (cid:48) ) (cid:105) ≡ δ ( ˆ n − ˆ n (cid:48) ) σ ( ˆ n ) [30], the modecoupling between different temperature multipole has thefollowing form: (cid:104) (cid:98) Θ (cid:96)m (cid:98) Θ (cid:96) (cid:48) m (cid:48) (cid:105) CMB = (cid:88) LM (cid:32) (cid:96) (cid:96) (cid:48) Lm m (cid:48) M (cid:33) f σ(cid:96)L(cid:96) (cid:48) σ ∗ LM , (A4)where σ LM is the harmonic coefficients of σ ( ˆ n ) and theweight function is defined as [29, 33]: f σ(cid:96)L(cid:96) (cid:48) = γ (cid:96)L(cid:96) (cid:48) (cid:32) (cid:96) L (cid:96) (cid:48) (cid:33) . (A5)Eqs. (A1) and (A4) mean that the lensing and pointsource fields also contribute to the estimator of the opti-cal depth anisotropies in Eq. (10). Using Eqs. (8), (A1)and (A4), the estimator of the optical depth anisotropieshas the following biases from the lensing and point sourcefields: (cid:104) (cid:99) δτ LM (cid:105) CMB = δτ LM + R δτ,φL φ LM + R δτ,σL σ LM . (A6)Here, the response function is given by Eq. (14).To mitigate the biases in the estimator of the opticaldepth anisotropies, let us consider biases to the lensingand point-source estimators, (cid:98) φ LM and (cid:98) σ LM . The lens-ing and point source estimators are obtained by simplyreplacing the weight function f δτ(cid:96)L(cid:96) (cid:48) with f φ(cid:96)L(cid:96) (cid:48) or f σ(cid:96)L(cid:96) (cid:48) inEq. (10). Using Eqs. (8), (A1) and (A4), the lensing andpoint-source estimators are biased by: (cid:104) (cid:98) φ LM (cid:105) CMB = φ LM + R φ,δτL δτ LM + R φ,σL σ LM , (A7) (cid:104) (cid:98) σ LM (cid:105) CMB = σ LM + R σ,φL φ LM + R σδτL δτ LM . (A8) TABLE IV. The PTE values for our measured cross-spectrum, C δτyL , with variation of analysis choices.PTEBaseline 0 . . y -map 0 . y -map 0 . . Eqs. (A6) to (A8) are summarized as the following matrixform: (cid:104) (cid:99) δτ LM (cid:105) CMB (cid:104) (cid:98) φ LM (cid:105) CMB (cid:104) (cid:98) σ LM (cid:105) CMB = R δτ,φL R δτ,σL R φ,δτL R φ,σL R σ,δτL R σ,φL δτ LM φ LM σ LM ≡ R L δτ LM φ LM σ LM . (A9)Motivated by the above equation, we can construct a setof “bias-hardened” estimators for reconstructing δτ , φ and σ : (cid:99) δτ BH LM (cid:98) φ BH LM (cid:98) σ BH LM = R − L (cid:99) δτ LM (cid:98) φ LM (cid:98) σ LM . (A10)The above estimators are unbiased by construction. If weonly consider mitigation of the lensing contribution in es-timating τ , the “lensing-hardened” τ estimator becomes[30]: (cid:99) δτ BH , lens LM = (cid:99) δτ LM − R δτ,φL (cid:98) φ LM − R δτ,φL R φ,δτL . (A11)The above estimators simultaneously estimate τ and lens-ing (and point source), and looses small amount of signal-to-noise compared to the standard quadratic estimator.The reconstruction noise level of (cid:99) δτ BH , lens LM is increased by ∼ −
40% depending on scales compared to (cid:99) δτ LM . Thenoise levels of (cid:99) δτ BH LM and (cid:99) δτ BH , lens LM are almost the same. Appendix B: Chi-squared PTE of measured spectra
Here we show the χ -PTE for the measured spectrum.The χ for the null spectrum is defined by (e.g. [68]): χ = (cid:88) bb (cid:48) (cid:98) C τyb { Cov − } bb (cid:48) (cid:98) C τyb (cid:48) , (B1)where b, b (cid:48) are the index of the multipole bin, (cid:98) C τyb is themeasured cross spectrum, and Cov is the band-power1covariance of the cross spectrum. The covariance is eval-uated from the simulation. We evaluate the PTE of theobserved χ value against the distribution of the simula-tion. Table IV summarizes the χ -PTE values for the mea-sured spectra shown in Fig. 2. From these values, all ofthe cross-spectra are consistent with null within 2 σ . [1] I. McGreer, A. Mesinger, and V. D’Odorico, Mon. Not.Roy. Astron. Soc. , 499 (2015), arXiv:1411.5375[astro-ph.CO].[2] J. E. Gunn and B. A. Peterson, Astrophys. J. , 1633(1965).[3] G. Kulkarni, L. C. Keating, M. G. Haehnelt, S. E.Bosman, E. Puchwein, J. Chardin, and D. Aubert,Mon. Not. Roy. Astron. Soc. , L24 (2019),arXiv:1809.06374 [astro-ph.CO].[4] Planck
Collaboration, Astron. Astrophys. , A108(2016), 1605.03507.[5]
Planck
Collaboration, Astron. Astrophys. (2018),1807.06209.[6] M. McQuinn, Ann. Rev. Astron. Astrophys. , 313(2016), arXiv:1512.00086 [astro-ph.CO].[7] X. Wang and W. Hu, Astrophys. J. , 585 (2006),arXiv:astro-ph/0511141 [astro-ph].[8] C. Dvorkin and K. M. Smith, Phys. Rev. D , 043003(2009), 0812.1566.[9] V. Gluscevic, M. Kamionkowski, and D. Hanson, Phys.Rev. D , 047303 (2013), 1210.5507.[10] T. Namikawa, Phys. Rev. D , 063505 (2018),1711.00058.[11] C. Feng and G. Holder, Phys. Rev. D , 123502 (2019),1808.01592.[12] S. Furlanetto, M. McQuinn, and L. Hernquist, Mon. Not.R. Astron. Soc. , 115 (2006), astro-ph/0507524.[13] O. Zahn et al. , Astrophys. J. , 12 (2006), astro-ph/0604177.[14] N. Battaglia, H. Trac, R. Cen, and A. Loeb, Astrophys.J. , 81 (2013), arXiv:1211.2821 [astro-ph.CO].[15] J. B. Wyithe and A. Loeb, Nature , 194 (2004),arXiv:astro-ph/0409412.[16] G. Mellema, L. Koopmans, H. Shukla, K. K. Datta,A. Mesinger, and S. Majumdar, in Advancing Astro-physics with the Square Kilometre Array (AASKA14) (2015) p. 10, arXiv:1501.04203 [astro-ph.CO].[17] J. D. Bowman, A. E. E. Rogers, R. A. Monsalve,T. J. Mozdzen, and N. Mahesh, Nature , 67 (2018),arXiv:1810.05912 [astro-ph.CO].[18] J. C. Hill, N. Battaglia, J. Chluba, S. Ferraro, E. Schaan,and D. N. Spergel, Phys. Rev. Lett. , 261301 (2015),arXiv:1507.01583 [astro-ph.CO].[19] E. J. Baxter, L. Weinberger, M. Haehnelt, V. Irˇsiˇc,G. Kulkarni, S. Pandey, and A. Roy, Mon. Not. R. As-tron. Soc. 10.1093/mnras/stab016 (2021), 2006.09742.[20] A. Mesinger, M. McQuinn, and D. N. Spergel, Mon.Not. R. Astron. Soc. , 1403 (2012), arXiv:1112.1820[astro-ph.CO].[21] N. Battaglia, A. Natarajan, H. Trac, R. Cen, andA. Loeb, Astrophys. J. , 83 (2013), arXiv:1211.2832[astro-ph.CO].[22] K. M. Smith and S. Ferraro, Phys. Rev. Lett. ,021301 (2017), 1607.01769.[23] M. A. Alvarez, S. Ferraro, J. C. Hill, R. Hloˇ zek, and M. Ikape, arXiv e-prints (2020), arXiv:2006.06594 [astro-ph.CO].[24] G. De Zotti, M. Negrello, G. Castex, A. Lapi, andM. Bonato, JCAP (3), 047, arXiv:1512.04816 [astro-ph.CO].[25] A. Roy, A. Lapi, D. Spergel, and C. Baccigalupi, J. Cos-mol. Astropart. Phys. , 014 (2018), arXiv:1801.02393[astro-ph.CO].[26] A. Roy, G. Kulkarni, P. D. Meerburg, A. Challinor,C. Baccigalupi, A. Lapi, and M. G. Haehnelt, arXiv e-prints (2020), arXiv:2004.02927 [astro-ph.CO].[27] M. Remazeilles, J. Delabrouille, and J.-F. Cardoso, Mon.Not. R. Astron. Soc. , 2481 (2011), 1006.5599.[28] G. Hurier, J. F. Mac´ıas-P´erez, and S. Hildebrandt, As-tron. Astrophys. , A118 (2013), 1007.1149.[29] Planck
Collaboration, Astron. Astrophys. , A15(2015), 1502.01591.[30] T. Namikawa, D. Hanson, and R. Takahashi, Mon. Not.R. Astron. Soc. , 609 (2013), 1209.0091.[31] W. Hu and T. Okamoto, Astrophys. J. , 566 (2002),astro-ph/0111606.[32] T. Okamoto and W. Hu, Phys. Rev. D , 083002 (2003),astro-ph/0301031.[33] S. J. Osborne, D. Hanson, and O. Dore, J. Cosmol. As-tropart. Phys. , 024 (2014), 1310.7547.[34] T. Namikawa and R. Takahashi, Mon. Not. R. Astron.Soc. , 1507 (2014), 1310.2372.[35] Planck
Collaboration, Astron. Astrophys. , A22(2016), 1502.01596.[36] J. C. Hill and D. N. Spergel, J. Cosmol. Astropart. Phys. , 030 (2014), 1312.4525.[37] N. Battaglia, J. R. Bond, C. Pfrommer, J. L. Sievers, andD. Sijacki, Astrophys. J. , 91 (2010), arXiv:1003.4256[astro-ph.CO].[38] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Good-man, Publications of the Astronomical Society of the Pa-cific , 306 (2013), arXiv:1202.3665 [astro-ph.IM].[39] F. B. Davies et al. , Astrophys. J. , 142 (2018),arXiv:1802.06066 [astro-ph.CO].[40] B. Greig, A. Mesinger, Z. Haiman, and R. A. Sim-coe, Mon. Not. Roy. Astron. Soc. , 4239 (2017),arXiv:1606.00441 [astro-ph.CO].[41] N. Battaglia, JCAP , 058, arXiv:1607.02442 [astro-ph.CO].[42] N. Battaglia, J. R. Bond, C. Pfrommer, and J. L. Sievers,Astrophys. J. , 75 (2012), arXiv:1109.3711 [astro-ph.CO].[43] J. S. Bolton, G. D. Becker, J. S. B. Wyithe, M. G.Haehnelt, and W. L. W. Sargent, Mon. Not. R. Astron.Soc. , 612 (2010), arXiv:1001.3415 [astro-ph.CO].[44] M. McQuinn, A. Lidz, M. Zaldarriaga, L. Hernquist, P. F.Hopkins, S. Dutta, and C. A. Faucher-Giguere, Astro-phys. J. , 842 (2009), arXiv:0807.2799 [astro-ph].[45] S. Hanany et al. (NASA PICO), (2019), 1902.10541.[46] N. Sehgal et al. , (2019), arXiv:1906.10134 [astro-ph.CO]. [47] K. M. Smith, D. Hanson, M. LoVerde, C. M. Hirata, andO. Zahn, J. Cosmol. Astropart. Phys. , 014 (2012),1010.0048.[48] E. Guzman and J. Meyers, (2021), 2101.01214.[49] E. Schaan, S. Ferraro, S. Amodeo, N. Battaglia, S. Aiola,J. E. Austermann, J. A. Beall, R. Bean, D. T. Becker,R. J. Bond, E. Calabrese, V. Calafut, S. K. Choi, E. V.Denison, M. J. Devlin, S. M. Duff, A. J. Duivenvoorden,J. Dunkley, R. D¨unner, P. A. Gallardo, Y. Guan, D. Han,J. C. Hill, G. C. Hilton, M. Hilton, R. Hloˇzek, J. Hub-mayr, K. M. Huffenberger, J. P. Hughes, B. J. Koop-man, A. MacInnis, J. McMahon, M. S. Madhavacheril,K. Moodley, T. Mroczkowski, S. Naess, F. Nati, L. B.Newburgh, M. D. Niemack, L. A. Page, B. Partridge,M. Salatino, N. Sehgal, A. Schillaci, C. Sif´on, K. M.Smith, D. N. Spergel, S. Staggs, E. R. Storer, H. Trac,J. N. Ullom, J. Van Lanen, L. R. Vale, A. van Enge-len, M. Vargas Maga˜na, E. M. Vavagiakis, E. J. Wol-lack, and Z. Xu, arXiv e-prints , arXiv:2009.05557 (2020),arXiv:2009.05557 [astro-ph.CO].[50] S. Amodeo, N. Battaglia, E. Schaan, S. Ferraro,E. Moser, S. Aiola, J. E. Austermann, J. A. Beall,R. Bean, D. T. Becker, R. J. Bond, E. Calabrese,V. Calafut, S. K. Choi, E. V. Denison, M. Devlin,S. M. Duff, A. J. Duivenvoorden, J. Dunkley, R. D¨unner,P. A. Gallardo, K. R. Hall, D. Han, J. C. Hill, G. C.Hilton, M. Hilton, R. Hloˇzek, J. Hubmayr, K. M. Huf-fenberger, J. P. Hughes, B. J. Koopman, A. MacIn-nis, J. McMahon, M. S. Madhavacheril, K. Moodley,T. Mroczkowski, S. Naess, F. Nati, L. B. Newburgh,M. D. Niemack, L. A. Page, B. Partridge, A. Schillaci,N. Sehgal, C. Sif´on, D. N. Spergel, S. Staggs, E. R. Storer,J. N. Ullom, L. R. Vale, A. van Engelen, J. Van La-nen, E. M. Vavagiakis, E. J. Wollack, and Z. Xu, arXive-prints , arXiv:2009.05558 (2020), arXiv:2009.05558[astro-ph.CO].[51] N. Sailer, E. Schaan, and S. Ferraro, Phys. Rev. D ,063517 (2020), arXiv:2007.04325 [astro-ph.CO].[52] A. Lidz and M. Malloy, Astrophys. J. , 175 (2014), arXiv:1403.6350 [astro-ph.CO].[53] A. D’Aloisio, M. McQuinn, O. Maupin, F. B. Davies,H. Trac, S. Fuller, and P. R. Upton Sanderbeck, As-trophys. J. , 154 (2019), arXiv:1807.09282 [astro-ph.CO].[54] M. McQuinn, Mon. Not. R. Astron. Soc. , 1349(2012), arXiv:1206.1335 [astro-ph.CO].[55] K. Finlator, L. Keating, B. D. Oppenheimer, R. Dav´e,and E. Zackrisson, Mon. Not. Roy. Astron. Soc. , 2628(2018), arXiv:1805.00099 [astro-ph.CO].[56] P. Gaikwad et al. , Mon. Not. R. Astron. Soc. , 5091(2020), arXiv:2001.10018 [astro-ph.CO].[57] J. S. Bolton, G. D. Becker, J. B. Wyithe, M. G. Haehnelt,and W. L. Sargent, Mon. Not. R. Astron. Soc. , 612(2010), arXiv:1001.3415 [astro-ph.CO].[58] H. Katz, T. Kimm, M. G. Haehnelt, D. Sijacki, J. Ros-dahl, and J. Blaizot, Mon. Not. R. Astron. Soc. ,1029 (2019), arXiv:1806.07474 [astro-ph.CO].[59] J. Tumlinson and J. Shull, Astrophys. J. Lett. , L65(2000), arXiv:astro-ph/9911339.[60] L. Pagano, J.-M. Delouis, S. Mottet, J.-L. Puget,and L. Vibert, Astron. Astrophys. , A99 (2020),arXiv:1908.09856 [astro-ph.CO].[61] C. Reichardt et al. (SPT), arXiv e-prints (2020),arXiv:2002.06197 [astro-ph.CO].[62] S. Henderson et al. , J. Low Temp. Phys. , 772 (2016),arXiv:1510.02809 [astro-ph.IM].[63] B. Benson et al. (SPT-3G), Proc. SPIE Int. Soc.Opt. Eng. , 91531P (2014), arXiv:1407.2973 [astro-ph.IM].[64] The Simons Observatory: Astro2020 Decadal ProjectWhitepaper , Vol. 51 (2019) 1907.08284.[65]
CMB-S4 Science Case, Reference Design, and ProjectPlan (2019) 1907.04473.[66] C. Feng and G. Holder, Phys. Rev. D , 123523 (2018),1801.05396.[67] M. Hazumi et al. , J. Low. Temp. Phys. , 443 (2019).[68] Bicep2
Collaboration, Phys. Rev. Lett.112