Angular Correlations of Cosmic Microwave Background Spectrum Distortions from Photon Diffusion

During cosmic recombination, charged particles bind into neutral atoms and the mean free path of photons rapidly increases, resulting in the familiar diffusion damping of primordial radiation temperature variations. An additional effect is a small photon spectrum distortion, because photons arriving from a particular sky direction were originally in thermal equilibrium at various spatial locations with different temperatures; the combination of these different blackbody temperature distributions results in a spectrum with a Compton $y$-distortion. Using the approximation that photons had zero mean free path prior to their second-to-last scattering, we derive an expression for the resulting $y$-distortion, and compute the angular correlation function of the diffusion $y$-distortion and its cross-correlation with the square of the photon temperature fluctuation. Detection of the cross-correlation is within reach of existing arcminute-resolution microwave background experiments such as the Atacama Cosmology Telescope and the South Pole Telescope.


INTRODUCTION
The average of two or more blackbody spectra at different temperatures is not a perfect blackbody.The Cosmic Microwave Background (CMB) photons arriving from a particular line of sight (LOS) scattered multiple times during recombination before their last-scattering into that LOS.Because CMB photons at the epoch of last scattering have energies that are far smaller than the electron rest mass, the fractional energy change of a scattered photon is small (of order   /   2 by the Klein-Nishina formula (Klein & Nishina 1929, as translated in Klein & Nishina (1994)).Therefore, to a good approximation, photon last-scattering does not change the photon's energy but only its propagation direction.
If the photons of the CMB came to us along each LOS unscattered from a "surface of last thermal emission," then they would be drawn from the thermal-equilibrium blackbody distribution in their local neighborhood of emission.Instead, the CMB photons are scattered with nearly no change in energy ★ E-mail:n.starkman@mail.utoronto.cainto the LOS at their point of last scattering.The photon energies are therefore drawn, approximately, from the blackbody distributions of the neighborhoods from which the photons originated.The energy distribution of the photons arriving along each LOS is approximately the average of many blackbodies, each with the temperature of a different neighborhood of photon emission.
A simple approximation to this averaging effect is to assume that each CMB photon scatters exactly once after its emission from a blackbody, or equivalently that the photon mean-free-path is negligible prior to its second-to-last scattering (2LS).We therefore take the point of emission, whose temperature distribution the photon samples, to be the point of second-last scattering.In this approximation, we provide an analytic expression for the contribution of any given Fourier mode of adiabatic perturbation to the photon spectrum from a given sky direction.The net effect of all such perturbation modes is an integral over the contribution of each mode, which can be evaluated numerically.The firstorder effect (in the difference between the temperatures in the second-last-scattering neighborhood and the global average temperature at that cosmic time) is a blackbody photon distribution with a temperature averaged over the region of second-to-last scatterings; this averaging approximates the mechanism behind the familiar diffusion damping of tem-perature anisotropies on small angular scales.The secondorder effect is a -distortion of the blackbody distribution first described (but not quantified) in (Zel'dovich et al. 1972), and recalled in (Chluba & Sunyaev 2004).More details were provided in (Khatri et al. 2012), but again without specific predictions.A calculation of the mean expected distortion may be among the effects included in (Chluba et al. 2012).Similarly, it is described in (Sunyaev & Khatri 2013;Chluba 2016).
We calculate this effect, showing that the -distortion angular correlation function is at a level that is potentially measurable in future experiments, but confirm that it is small compared to other expected signals (Chluba & Sunyaev 2004).We also calculate the angular cross-correlation function of this y-distortion and the square of the temperature fluctuations.The cross-correlation is likely detectable in current experiments, including the Atacama Cosmology Telescope (ACT) (Coulton et al. 2023) and the South Pole Telescope (SPT) (Bleem et al. 2022), and especially in anticipated experiments like the Simons Observatory (Galitzki et al. 2018) and CMB-S4 (Abitbol et al. 2017b).The authors of (Chluba et al. 2022) and (Kite et al. 2022) address related questions, including cross-correlations, although not specifically this effect.
The mixing of blackbody photon distributions is distinct from the averaging of temperatures over different lines of sight, whether through the finite width of telescope beams (e.g.Chluba & Sunyaev 2004) or through the angular integration inherent in the calculation of the coefficients of spherical harmonics (e.g.Lucca et al. 2020).Each of these effects also creates small -distortions in measured signals through temperature averaging, but the effect discussed here is due to physical processes during the epoch of last scattering and thus carries information about the universe rather than about a given measurement.
In § 2, we calculate the distribution of second-lastscattering locations.In § 3, we calculate the expected distortion Y of the CMB spectrum due to this diffusive averaging.In § 4, we calculate the expected angular correlation function of Y, as well as the expected cross-correlation between Y and the square of the observed temperature fluctuation (Δ) 2 .We discuss the detectability of this signal in § 5, and conclude in § 6.

SECOND-TO-LAST SCATTERING DISTRIBUTION
Photons arriving at a common point of last scattering x 1 come from a surrounding neighborhood of prior locations x 2 , where we use the subscript to label the scattering number, counting back in time from observation at scattering "0".Arriving at x 1 , the photons scatter into the LOS and their combined distribution is no longer that of a perfect blackbody; it is an average of blackbodies.The resulting distortion of the blackbody spectrum retains information about the distribution of temperatures in the second-last-scattering neighborhood that goes beyond the intensity of the photon emissions from that LOS.
To calculate the observational effects of this blackbody averaging, we must first calculate the probability density P (x 2 , x 1 ) that a photon arriving at the observer from a last- (2  P ).P peaks at  || ≃ 2.76 -just slightly higher than  recombination = 2.73 (dashed line) -and at  ⊥ ≲ 10 −3 (below which it is approximately constant).When sampling from P each decade in  ⊥ contributes approximately one-tenth the samples as the previous decade; thus we display only  ⊥ ≥ 10 −7 .The other boundaries of the domain of P are taken at sufficient distance from the higher probability regions of P that enlarging the domain does not impact the results.
scattering location x 1 had its second-last scattering at x 2 .While we refer to our location as "the observer," in order to remove the physically distinct effects of propagation through the low-redshift universe, we consider a reference point x 0 at co-moving distance  0 from us along the LOS, which is the line between the observer and x 1 .We take x 0 to be at a redshift  0 = 100, well after recombination and last scattering (at  ≃ 10 3 ) but well before the onset of cosmic acceleration.The comoving displacement of the j-th-last scatter from the i-th-last-scattering location is x ij ≡ x j − x i .We also write    for |x ij |.
To leading order (in fluctuations around the homogeneous background metric and stress-energy tensor), the probability distribution P (x 2 , x 1 ) is rotationally symmetric about the line-of-sight axis.P (x 2 , x 1 ) can thus be written in terms of: (i)  1 ( 01 ) d 01 , the probability that last scattering occurs at a comoving distance from x 0 between  01 and  01 + d 01 along the LOS where   (  ,   ) is the visibility function for Compton scattering between comoving radii   and   ≥   (see Appendix § B); (ii)   () d, the probability that the last scattering is through an angle between  and  + d (iii)   (), the probability distribution of the last scattering azimuth (4) The overbar on   in eq. 1 and 4 indicates that we are considering   to be a function only of redshift, i.e. a homogeneous function of location on a given redshift slice.The arguments of   are cosmological epochs, not positions; thus in eq. 4,  01 +  12 ≠  02 .From eq. 1, it is clear that   (  ,   ) is normalized such that ∫ d    (  ,   ) = 1.We note that  2 is independent of scattering angle  (to leading order).
Multiplying these four probabilities - 1 ,   ,   , and  2 -gives the probability density P for the last scattering to take place at comoving distance  01 from the observer and the second-last scattering to take place at comoving distance  12 from the last scattering, with scattering angle  and azimuthal angle .
We prefer to express P in terms of new dimensionless comoving coordinates  | | and  ⊥ measured respectively along the LOS (in the n direction) and perpendicular to it, instead of  12 and .Thus where  eq is defined in Appendix § A. We define the vector s, with origin at r 0 , i.e. along the LOS at redshift  0 .s is directed away from the reference point and away from the observer; it can be decomposed as s ⊥ is a two-vector of magnitude  ⊥ in the plane perpendicular to n.The direction of s ⊥ is specified by .
eq cos  and  ⊥ =  12  eq sin . (7) We can take the new set of variables to be  01 ,  ⊥ ,  | | , and , leading to a Jacobian factor In terms of these new variables The probability density The numerical evaluation of P ( | | ,  ⊥ , ) is described in Appendix § C, and the result is presented in Fig. 1.We observe that P is sharply peaked close to the position of recombination, i.e at  | | ≃  recombination , and with  ⊥ ≲ 10 −3 .

CALCULATING THE SIGNAL
Consider a sum of blackbodies, each with photon occupation number   ( ) = (  − 1) −1 for  ≡ ℎ/  .To third order in the temperature fluctuations   ≡   − T T : We can choose to normalize   so that    = 1.Taking T ≡      as one would expect, then      = 0.
It is conventional and convenient to reorganize eq. 12, writing: where We are free to choose T and do so such that This choice of T means eq. 13 is a pure -distortion: Clearly T ≠ T, but how different are they, or more importantly how different are δ2 ≡    2  and ε2 ≡     2  ?We can solve the quadratic eq. 15 for T and expand T around T to find It is straightforward to show that and thus For small fluctuations, we can therefore replace T by T, and   by   in the calculation of the signal, which proves considerably simpler.
In the CMB, photons were scattered into the LOS at lastscattering from different locations of second-last scattering.We take the photon distribution originating from each 2LS point to be a blackbody of the temperature characteristic of the plasma there.In this approximation, the resulting observable photon distribution is the weighted average of the 2LS blackbodies at all the accessible 2LS points.We must therefore account for the temperature at each location x(s, n) that can scatter into the LOS, and we must weight the sum by the probability P (s) that the photons we see come from that location. 2 The temperature at x(s, n) is given by where  0 is the mean CMB temperature. T (k) includes the effect of the transfer function T (|k|) on the primordial amplitude (k) of the Fourier mode of the primordial curvature fluctuation with wave vector k: The appropriate weighted mean temperature of the photons scattered into the LOS in direction n is The -distortion signal is where again ds Expectation values pass through integrals to apply only to the factors of Δ, so that ⟨Y (n)⟩ can be expressed in terms of the correlation function  (|x 2 − x 1 |) of the Gaussian Δ field: where we have used the fact that Unsurprisingly, ⟨Y⟩ is independent of direction, since we have assumed here that  is statistically isotropic, at least on the small scales over which diffusion takes place during recombination: The TT power spectrum   is the product of the initial power spectrum of the gauge-invariant curvature perturbations times the "early" temperature transfer function squared Here  ★ = 0.05ℎ −1 Mpc is an arbitrary, but conventional, pivot scale.The high-k cutoff  2LS , accounting for the damping up to second-last-scattering, can be written in terms of the conformal-time rate of change of the opacity (Jungman et al. 1996), or more simply read directly from Baumann (2022, eq.7.139) We take the transfer function T  () to be approximately the pure Sachs-Wolfe power spectrum, and, assuming that it changes slowly with time, evaluate it at the time of last scattering, rather than the time of second-last scattering (we use Baumann 2022, 7.112): The label  stands for early, i.e. during recombination, as opposed to the usual transfer function T  (), evaluated late, i.e. at redshift  0 = 100.Here and the sound horizon at emission (see Baumann 2022, Appendix C) is ≃ 144 Mpc . Thus, The result is independent of the low- cutoff for Since the variance in the (dipole-subtracted) CMB temperature anisotropies, var() we might have thought that the maximum signal would be However, the observed CMB fluctuations themselves are damped compared to their amplitude at second-last scattering.Without that damping (i.e.eliminating the exponential term in eq.41), the signal would have been Damping of short range fluctuations is thus responsible for a factor-of-31 suppression in the signal from its maximum possible value.

CORRELATION FUNCTION OF THE SIGNAL AND TEMPERATURE FLUCTUATIONS
Unlike the -distortions that arise from the time-evolution of the background physics during recombination, the distortions due to local blackbody averaging are anisotropic.We are therefore interested in calculating the angular correlation function of Y (n) with the temperature fluctuations themselves Δ (n ′ ).The particular form of this crosscorrelation is distinctive to the diffusion distortion signal, and can be used to disentangle diffusion distortion from foregrounds.
Since the fluctuations are nearly Gaussian, the expected correlation of Y (n) with Δ (n ′ ) is a three-point function, and nearly vanishes.However, the correlation of Y (n) with (Δ (n ′ )) 2 does not: (See Appendix § D for a full derivation.)Given a statistically isotropic universe,  Y (Δ ) 2 (n 1 , n2 ) depends only on n1 • n2 ≡ cos .Here with T  () obtained by evaluating T  () from ( 35) at  0 .

DETECTABILITY
The Y power spectrum is well below the detection capability of envisioned experiments, and is small compared to Y from low-redshift contributions.But the Y (Δ) 2 cross power spectrum is more promising, because it has both a significantly larger amplitude and a distinctive shape, including slight imprints of acoustic oscillations, which can distinguish it from other contributions.
We estimate the detectability of  Y ℓ (Δ) 2 by assuming that both Y and (Δ) 2 are Gaussian random fields on the sky.While this is not precisely the case, it is a large simplification and should be sufficient for a reasonable signal-to-noise estimate.If the CMB blackbody temperature in a given sky pixel can be measured with an uncertainty of   , then the uncertainty in (Δ) 2 is 2Δ   +  2  .Any map with a chance to detect the signal described here will have   ≪ Δ so the second term can be dropped.Since Δ is itself normally distributed, it can be replaced by its rms value, so the uncertainty in (Δ) 2 is   2 ≃ 2Δ rms   .For a map of Y, each pixel corresponds to a particular intensity distortion at each frequency, given by the usual thermal Sunyaev-Zeldovich formula.Taking 90 GHz as a typical frequency for ground-based experiments, an uncertainty in temperature   gives   = 1.6  / 0 as the uncertainty in Y.
Given these pixel errors in two (approximately) Gaussian random fields, the uncertainty in measuring each ℓ mode of the angular cross-power spectrum  Y (Δ ) 2 ℓ is analogous to determining    ℓ from maps of microwave background temperature and E-mode polarization.The variance, including both pixel noise and cosmic variance, is given by Kamionkowski et al. 1997, Eq. (3.26), as Here  −1  ≡ 4( 2  )/ pix is the inverse statistical weight per unit solid angle on the sky for a map of some quantity with  pix pixels and a pixel variance of  2 , and   ℓ = exp(−ℓ 2  2  /2) is the experiment beam profile in harmonic space, which is generally well approximated as a Gaussian with beam with   =  fwhm / √ 8 ln 2 = 1.2 × 10 −4 ( fwhm /1 ′ ).The signal-to-noise with which a given ℓ mode can be measured is then just where a factor of  sky , giving the sky fraction covered by a map, has been included.To the extent the fields Y and (Δ) 2 are Gaussian, each mode is statistically independent, and the total signal-to-noise can be obtained by summing over all ℓ values probed by a given map.
For these parameters, the terms in eq.54 containing and  Y Y ℓ can be neglected, along with the 1, so the expression simplifies to Summing from ℓ = 500 to ℓ = 1500 gives an expected total S/N of 12 for ACT.So the diffusion spectral distortion should be detectable at current experiment sensitivities.SPT has released (Bleem et al. 2022)  , and might therefore allow Y to be used as a cosmological probe analogous to polarization.
In forecasting the S/N, we have made several assumptions.Primarily, we have assumed that Y can be detected with uncertainty comparable to that in Δ.That requires robust foreground subtraction, since Y is large in the location of galaxy clusters, and the Y map will be dominated by distortions from foreground haloes.Since we are correlating with (Δ) 2 , modeling this signal due to foreground Y contributions will likely be the primary challenge in extracting the recombination-era signal.We have also worked in the approximation that Y and (Δ) 2 are Gaussian, which they are not, and a more careful statistical treatment is merited.Finally, we have included only the Sachs-Wolfe contribution to the transfer function, and used an analytic approximation.

DISCUSSION AND CONCLUSIONS
In this paper we provide a real-space description of the -type spectral distortion of the CMB that arises from the scattering of the photons into our line of sight (Zel'dovich et al. 1972).Each photon's scattering history is different, sampling both radially along and transverse to the line of sight through the recombination era.Because the photons' final scatterings make nearly no change to the photon energy, the resulting distribution of photon energies that we observe is a mixture of blackbody distributions of different temperatures, representing the inhomogeneity of the temperature in the region from which the photons originated.
This "diffusion spectral distortion" is, like the CMB intensity (temperature) and polarization, a probe of the acoustic modes responsible for the inhomogeneities in the universe during the epoch of recombination.Like the E-mode polarization and the temperature, in standard Λ Cold Dark Matter cosmology, the diffusion spectral distortion signal is partly, but not perfectly, correlated with these other signals.
For example, while the E-mode fluctuations probe the local quadrupole of the temperature distribution, the diffusion spectral distortion is also sensitive to its dipole.
Diffusion spectral distortion offers another independent probe of the physics at the end of recombination.Inherently, this implies the opportunity to reduce cosmic variance on existing measurements of quantities probed by the temperature and polarization.Additionally, because it is sensitive to the variation in the photon temperature along the line-of-sight, diffusion spectral distortion could potentially be the basis of a new Alcock-Paczynski test (Alcock & Paczynski 1979) at the epoch of recombination.It could also allow us to test statistical isotropy at the epoch of recombination, complementing the large scale tests traditionally done with temperature that have yielded anomalous results (Schwarz et al. 2016;Planck Collaboration et al. 2020b;Abdalla et al. 2022).Diffusion spectral distortion is therefore both a consistency check for the standard cosmological model and sensitive to new physics in a way that is complementary to other signals.
Polarization is also generated during recombination, and usually described as a measure of the local quadrupole at the last scattering.The polarization spectrum and its distortion away from blackbody are therefore due to a different weighted sum over blackbodies, with the potential for yet another signal complementary to the one we have described.This is another step towards a possible tomographic probe (Yadav & Wandelt 2005) of the universe through recombination.
Current and upcoming experiments will not be sufficiently sensitive to the -type distortion to measure the diffusion spectral distortion signal directly.In part, the challenge is insufficient signal-to-noise.However, a greater challenge is likely to be separation of the component of the -distortion due to recombination-era diffusion from other causes of spectral distortion in the presence of foregrounds with spectra that are imperfectly known and that are also anisotropic (Abitbol et al. 2017b;Hart et al. 2020).These foregrounds include our Milky Way, and both galaxies and galaxy clusters at all redshifts.The Y auto-correlation function, like the -type distortion, will not be detectable by current or upcoming experiments.The amplitude of the correlation signal is far smaller than for other sources of  auto-correlation at low redshifts, and so will be difficult to separate from foregrounds.
More promising than the Y signal itself or its autocorrelation function is the cross-correlation  Y (Δ ) 2 ℓ between the -distortion and the square of the temperature fluctuations.(The Y- correlation will be zero if the primordial photon perturbations are a Gaussian field.)This is calculated in § 4. Like the temperature fluctuations, the diffusion spectral distortion, and hence their cross-correlation function, contains acoustic features.The combination of the specific spectral shape and the correlation with the primordial temperature fluctuations should facilitate the separation of this correlation function from foreground signals.Galactic foreground confusion is uncorrelated with primordial temperature anisotropies (though of course it would be correlated with any residual unsubtracted Galactic temperature foreground).Extragalactic foregrounds are concentrated at galaxy clusters (which are detected at high significance) and at galaxies.Fortunately, the clusters can be masked and are relatively sparse on the sky.Galaxies are far more numerous and non-sparsely distributed; however, the amplitude of the galaxy confusion limit -distortion is also small and would have a different correlation function with the temperature fluctuations.
The three-point correlation function between Y (n 1 ), Δ (n 2 ), and Δ (n 3 ) at three different points may also be detectably large.Considering different configurations of this three-point function will give further handles on separating the diffusion distortion from various foreground contributions (Coulton et al. 2018).Calculation of this signal is more complicated than the two-point function and will be considered elsewhere.
In § 5 we determined that multi-frequency maps over a third of the sky with the 10K-arcmin sensitivity attained by the final ACT dataset (Coulton et al. 2023) have sufficient sensitivity to detect this signal, and SPT has achieved similar reach (Bleem et al. 2022).Sky maps with sensitivities approaching 1 K-arcmin are anticipated in the next decade (Abitbol et al. 2017a).Sufficient frequency coverage is required to separate the -distortion from the primary blackbody and other spectrum components.
Our signal and sensitivity calculations involved several simplifying assumptions.More rigorous calculations of the distortion and its correlation with the temperature and polarization anisotropies during diffusion damping are warranted.Ultimately, a complete calculation of the statistics of spectral distortions arising from physical processes around last scattering may reveal additional probes of the cosmological model, providing substantial consistency checks or additional handles on non-standard physics.The data availability statement is modified from one provided to showyourwork by Mathieu Renzo.

DATA AVAILABILITY
This study was carried out using the reproducibility software show your work!(Luger et al. 2021), which uses contin- uous integration to programmatically download the data, perform the analyses, create the figures, and compile the manuscript.Each figure caption contains two links: one to the dataset used in the corresponding figure, and the other to the script used to make the figure.The datasets are stored at https://zenodo.org/record/8400583.The git repository associated with this study is publicly available at nstarman/Temperature-Diffusion-Spectral-Distortion-Paper.
The overbar in μ,   , and   indicates that they are functions only of  and not position, as they refer to the unperturbed "background cosmology".
For numerical and analytic purposes, we want to "split" our background quantities, i.e. we rewrite   ( 1 ,  2 ) (where in terms of CLASS's     () which is anchored at the observer   : For this we write Similarly, for   : and so APPENDIX C: CALCULATING AND SAMPLING P ( | | ,  ⊥ , ) For calculating and sampling P, it is convenient to perform a change of coordinates.In eq.7 we define  | | ,  ⊥ , and in eq.A3 we define .We define normalized such that  eq ∫ d 2   ( 2 ,  1 ) = 1.Rewriting eq. 10 in terms of   and the CLASS-defined functions, which we evaluate as follows.For a given  | | ,  ⊥ we perform a cubic spline in  1 of The integrals with different   1 can all be done analytically: We take particular care to ensure that  1 =   −  | | is not one of the knots of the spline, because the integral is numerically unstable (though not analytically problematic) if an integral bound (knot point) is too close to this value.
To do the full integral we need to know the range of  | | and  ⊥ over which we should evaluate P. We therefore identify a range  min <  <  max over which   is above its value at  0 .Requiring that Figure 1.P ( || ,  ⊥ ,  = 0) over a range of  || ,  ⊥ , colored by log 10(2  P ).P peaks at  || ≃ 2.76 -just slightly higher than  recombination = 2.73 (dashed line) -and at  ⊥ ≲ 10 −3 (below which it is approximately constant).When sampling from P each decade in  ⊥ contributes approximately one-tenth the samples as the previous decade; thus we display only  ⊥ ≥ 10 −7 .The other boundaries of the domain of P are taken at sufficient distance from the higher probability regions of P that enlarging the domain does not impact the results.

Figure 3 .
Figure 3. Upper panel: Legendre expansion coefficients (cf.(52)) of the correlation functions  YY (  )/ 2 ACKNOWLEDGEMENTS N.S. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) -Canadian Graduate Scholarships Doctorate Program [funding reference number 547219 -2020].N.S. received partial support from NSERC (funding reference number RGPIN-2020-04712) and from an Ontario Early Researcher Award (ER16-12-061; PI Bovy).N.S. would also like thank Prof. Jeremy Webb for providing computational resources.G.D.S. was partially supported by DOE grant DESC0009946.G.D.S. thanks Imperial College London for hospitality while some of this work was completed.