Patchy screening of the CMB from dark photons

We study anisotropic (patchy) screening induced by the resonant conversion of cosmic microwave background (CMB) photons into dark-sector massive vector bosons (dark photons) as they cross non-linear large scale structure (LSS). Resonant conversion takes place through the kinetic mixing of the photon with the dark photon, one of the simplest low energy extensions to the Standard Model. In the early Universe, resonant conversion can occur when the photon plasma mass, obtained as the photon propagates through the ionized interstellar and intergalactic media, matches the dark photon mass. After the epoch of reionization, resonant conversion occurs mainly in the ionized gas that occupies virialized dark matter halos, for a range of dark photon masses between 10-13 eV ≲ m A' ≲ 10-11 eV. This leads to new CMB anisotropies that are correlated with LSS, which we refer to as patchy dark screening, in analogy with anisotropies from Thomson screening. Its unique frequency dependence allows it to be distinguished from the blackbody CMB. In this paper, we use a halo model approach to predict the imprint of dark screening on the CMB temperature and polarization anisotropies, as well as their correlation with LSS. We then examine the two- and three-point correlation functions of the dark-screened CMB, as well as correlation functions between CMB and LSS observables, to project the sensitivity of future measurements to the kinetic mixing parameter and dark photon mass. We demonstrate that an analysis with existing CMB data can improve upon current constraints on the kinetic mixing parameter by two orders of magnitude with the two-point correlation functions, while data from upcoming CMB experiments and LSS surveys can further improve the reach by another order of magnitude with two- and three-point correlation functions.


Introduction
Cosmic microwave background (CMB) and large scale structure (LSS) surveys have provided some of the most important evidence for physics beyond the Standard Model (BSM) of particle physicsdark matter and dark energy, as well as early-Universe cosmological models such as inflation.These observations also provide powerful constraints on BSM physics, such as the number of relativistic degrees of freedom, neutrino masses, and the strength of interactions during an inflationary era.The next generation of CMB experiments such as the upcoming Simons Observatory [1], as well as CMB-S4 [2] and CMB-HD [3], and LSS surveys like DESI [4], Euclid [5], and LSST [6] promise to further expand our ability to detect and characterize BSM physics.These new surveys motivate the exploration of new observables that can exploit their full potential.In particular, a promising avenue of current and future effort is to use the CMB as a back-light through which to study the intervening LSS (see e.g.[7,8]).
Within the Standard Model of cosmology and particle physics, CMB photons can interact with LSS via gravity (e.g.weak lensing or the integrated Sachs-Wolfe effect) and electromagnetism (e.g.Sunyaev Zel'dovich effects).These effects lead to additional temperature and polarization anisotropies in the CMB, the so-called secondary CMB, as well as new correlations with tracers of LSS such as galaxy surveys.The slew of associated new observables, and in particular cross-correlations between CMB and LSS, can be used to extract valuable information about the initial conditions and the formation and growth of structure in the early Universe.Several examples include: lensing reconstruction (see e.g.[9] for a review), kinetic Sunyaev Zel'dovich velocity reconstruction [10][11][12][13][14], moving-lens velocity reconstruction [15], and patchy reionization optical depth reconstruction [16].
In this paper, we discuss the possibility of using the CMB and its cross-correlation with tracers of LSS to extract information about BSM physics that manifest in the low-redshift Universe.In particular, we will discuss how a new type of CMB secondary anisotropy and its correlation with LSS can be used to extend the reach in the parameter space of kinetically mixed dark photons [17,18] by orders of magnitude.
The dark photon is a hypothetical vector boson that arises in various extensions of the Standard Model [19,20].An ultra-light dark photon is an essential ingredient in dark matter models, either as a light bosonic dark matter candidate [21][22][23], or as a mediator to a sector of dark matter particles (see [24] and reference within).Despite recent evidence of various collective effects which cast doubt on the validity of some of the models (vortex production during production of dark photon dark matter [25] and the two stream instability in the case of freeze-in dark matter models [26]), it is still of great interest to probe the existence of ultra-light dark photons regardless of their cosmological abundance through superradiance [25,[27][28][29][30], cosmology [31][32][33][34], stellar objects [35][36][37] and lab searches [38][39][40][41][42].
A dark photon and its coupling to the Standard Model can be described by the Lagrangian where A µ and A ′ µ are the photon and dark photon fields respectively, with F µν and F ′µν their field strengths, and J µ is the Standard Model electromagnetic current.The dark photon has a mass m A ′ and couples to the Standard Model photon through a kinetic mixing parameter ε.This simple coupling leads to a plethora of observable consequences (see [43] and references within).Most of these are based on the conversion, in particular resonant conversion, between the photon and the dark photon in a medium.In the lab, resonant conversion is facilitated with carefully prepared small scale experiments.In the early Universe, resonant conversion happens in different astrophysical and cosmological environments, as the dispersion relation of the photon (plasma frequency) is naturally scanned.When the plasma frequency of the photon m 2 γ over its trajectory ⃗ x matches the mass of the dark photon CMB photons resonantly convert into dark photons.
In [31][32][33][34]44], the resonant conversion between CMB photons and dark photons was studied both in the homogeneous early-time and inhomogeneous late-time limits.Cosmic expansion and the inhomogeneous distribution of ionized gas were identified as important scanners of the plasma frequency.In this paper, we examine the conversion from CMB photons into dark photons inside non-linear structure, after the epoch of reionization, where the amplitude of the density profile of ionized gas within dark matter halos provides the primary scanner.
Resonant conversion leads to a frequency-dependent disappearance of CMB photons that traces the distribution of matter in the Universe.Hence it can be treated as a frequency (ω) and angle (n) dependent optical depth τ (ε, ω, n), after integrating along the line of sight.The optical depth from resonant conversion can be extracted or constrained from cosmological data, and we present five methods to search for dark photons, along with the projected sensitivity on ε: • Spectral distortions of the CMB: The spatially averaged ⟨τ (ε, ω)⟩ manifests as a distortion of the blackbody spectrum of the CMB, and as a result is constrained by COBE/FIRAS [45].This effect, the late-time component of the effect studied in [31][32][33][34], scales as ε 2 .
• CMB temperature and polarization anisotropies: The optical depth correlation function ⟨τ (ε, ω, n) τ (ε, ω, n′ )⟩ can be obtained from the measured CMB through the screening of the temperature and polarization anisotropies by resonant conversion.Data from CMB experiments can be used to extract the amplitude of dark screening, which scales like ε 4 .We show how the large signal-to-noise ratio of CMB surveys, along with the characteristic frequency dependence of the screening signal, implies that this method outperforms the COBE/FIRAS constraint.
• Correlating CMB anisotropies with templates from LSS: As described in greater detail below, the morphology of the dark photon optical depth anisotropy depends on the distribution of ionized gas in halos.With assumptions about the galaxy-gas connection, a galaxy survey can be used to create a template for the dark photon optical depth field τ (ω, n).Cross-correlating the CMB measurement with this template ⟨τ (ε, ω, n) τ (ω, n′ )⟩ scales as ε 2 , improving greatly on the CMB-only reach.
• Correlation with Thomson screening: The standard optical depth due to Thomson scattering by free electrons (we will denote by τ Th ) is also present as a source of screening in the measured CMB.Since this anisotropic signal traces the same distribution of ionized matter as the dark screening component, the two signals will be correlated yet distinguishable due to the latter being frequency dependent.The cross-correlation ⟨τ (ε, ω, n) τ Th (n ′ )⟩ also scales as ε 2 .
• The CMB bispectrum and optical depth reconstruction: Both conversion to dark photons and Thomson screening induce non-Gaussian statistics in the CMB anisotropies.The combined effect can be modeled via three-point correlation functions (bispectra) that also scale as ε 2 .These bispectra hold additional information compared to two-point functions since there are more modes; the associated statistical anisotropy can additionally be used to reconstruct the dark screening optical depth, allowing for its study at the field level.
In this paper, we demonstrate the possibility of using the aforementioned methods to improve the reach on kinetically mixed dark photon in the mass range (10 −13 eV ≲ m A ′ ≲ 10 −11 eV).The paper is organized as follows.We first review resonant photon to dark photon conversion and compute the properties of conversion inside individual halos in Section 2, before summing over halos to obtain a frequency dependent dark screening optical depth in Section 3. In Section 4, we discuss the anisotropies of this dark screening optical depth, correlation functions, and CMB observables.In Section 5 we study the cross-correlation between this dark screening optical depth and the LSS of our Universe, and construct two-point cross-correlation functions between the CMB and LSS, as well as three-point cross-correlation functions of CMB observables.In Section 6, we present a forecast of the sensitivity of existing and future CMB data-sets to the various correlation functions studied in this paper.The result of these forecasts are shown in Section 7, along with a discussion of the prospect for constructing similar correlation functions in other new physics scenarios.In the appendix, we present details about the modeling of dark matter and gas halos (Appendix A) and the computation of correlation functions of dark screening (Appendix B), as well as a list of useful two-point correlation functions and quadratic estimators for the optical depth and other quantities (Appendix C).

Photon to dark photon conversion
In this section, we discuss how photons resonantly convert into dark photons within non-linear structure in the context of the halo model of LSS, where matter is organized into virialized dark matter halos populated by gas (see e.g.[46,47] for a review).Resonant conversion can be modeled via the same formalism that describes neutrino oscillation in medium, that is, the Mikheyev-Smirnov-Wolfenstein (MSW) effect [48,49].In the following, we review the resonant conversion of photons into dark photons and present our prescription for modeling the conversion probability as a sum over halos.

Resonant conversion probability
In an ionized medium, the Lagrangian in Eq. (1.1) leads to resonant conversion of photons to dark photons.This can be described by the Schrödinger equation [31]: where γ is an incident photon with frequency ω(t) that follows a trajectory ⃗ x parameterized by time t.The photon acquires an effective mass m 2 γ (⃗ x(t)) (plasma frequency) as it crosses an ionized medium due to its interaction with the collective oscillations in the free electron density.Hence, to first order, the mass depends on the number density of electrons n e (⃗ x(t)) along its trajectory: Here we assume all baryonic matter is ionized and therefore ignore an additional negative contribution due to interactions with neutral atoms [31].In the small-ε limit where conversion from dark photons back to photons can be safely neglected, the conversion probability is given by where t res are the times when the resonance condition m 2 γ (t res ) = m 2 A ′ is met along the path ⃗ x.This expression for the total probability is a combination of the conversion rate Γ res = πε m 2 A ′ /ω(t res ) and the resonance time scale t=tres .

Photon to dark photon conversion in non-linear structure
In the homogeneous and weakly inhomogeneous early Universe (z ≫ 10), the slowly diluting charged particle density caused by cosmic expansion provides a natural scanner of the dark photon mass, and ensures efficient conversion for a wide range of dark photon masses.At low redshift (z ≲ 10), the scanner is mainly provided by the spatially varying electron density inside non-linear structure, e.g.halos, which also ensures efficient conversion over a range of dark photon masses due to the large density contrast.Expression (2.3) is the integrated probability to convert along the line of sight.Working in the halo model for LSS [46], this expression becomes a sum over halos, where each term represents the probability that a photon converts within each.We re-write the probability per halo in terms of the mass and redshift: This expression holds for any type of photon, but in what follows we will focus on the conversion of CMB photons along their path from the surface of last scattering to the Earth.In the remainder of this sub-section we explain how to simplify the term in the modulus to account for a photon's path across each halo.The effective photon mass m 2 γ (⃗ x) depends on the baryon number density as well as the ionization fraction.In a galactic halo, baryonic matter represents a fraction Ω b /Ω c ∼ 0.19 of the total halo mass m.Since baryonic matter is predominantly protons by mass, and the Universe is electrically neutral, we approximate the number density of the electrons to be the same as the number density of baryons.Furthermore, we are interested in the period after reionization (z ≲ 6 − 10) therefore we assume throughout that the ionization fraction is unity everywhere.Where relevant, we treat reionization as instantaneous at a redshift in the range 6 < z < 10 to encapsulate uncertainties about the history of reionization.Finally, we neglect the impact of He reionization.Further details of our modeling of reionization and other assumptions can be found in Appendix A.
For the density profile of baryons, we use the Battaglia et al. 'AGN Feedback' gas density profiles introduced in [50], which are based on hydrodynamic cosmological simulations.We use a version of the profile where the fit parameters are based on simulations that include a sub-grid model for active galactic nuclei (AGN) feedback.The profile is given by an expression that parametrically resembles the NFW profile: (2.5) The quantity ρ c (z) is the critical density for a flat FRW Universe and R 200 is the radius where the gas density reaches 200ρ c .The exponents α, β, γ fix the slope in the regimes where x ∼ 1, x ≫ 1 and x ≪ 1, respectively.There are two fixed quantities γ = −0.2 and the core scale x c = 0.5 that control the central region in each halo.The remaining functions ρ 0 (z, m), α(z, m) and β(z, m) are fit with power laws: where the best-fit parameters {a, b, c} are in each case as follows: for ρ 0 {4000, 0.29, −0.66}, for α {0.88, −0.03, 0.19}, and finally for β {3.83, 0.04, −0.025}.These values are taken from Table 2 of [50]. 1s discussed in detail in later sections, the assumptions about how gas inhabits dark matter halos has a significant effect on the signal and the resulting sensitivity.An extreme case to contrast with is to assume that baryons track the dark matter density everywhere.The Standard parametric expression for the dark matter density in halos is the Navarro-Frenk-White (NFW) profile [51]: Each halo has a physical scale radius and density that depend on its mass and redshift, i.e. r s (z i , m i ) and ρ s (z i , m i ).The assumption that the baryons follow dark matter is expected to be reasonable in the outer regions of halos (e.g.beyond the scale radius).However it is no longer valid in the inner regions where baryonic feedback processes are non-negligible.We present a discussion on the model uncertainty in our sensitivity due to varying assumptions about the electron profile in Appendix A, using the NFW and AGN Feedback models as a proxy for the span of models.
The effective mass of a photon crossing a halo i centered at redshift z i with mass m i depends on the number density of electrons n e along its path as defined in Eq. (2.2).Assuming all baryonic mass inside halos is contained within protons, that are as numerous as electrons, we can write this in terms of the halo gas density profile: where the term in brackets is dimensionless and we have defined κ = 5.7 × 10 −38 eV 2 .For resonant conversion that occurs at a radius r res , the resonance time scale can be broken into a radial and an angular part: where the velocity term |dr(t)/dt| in Eq. (2.9) depends on the precise photon trajectory through the halo and encapsulates the angular dependence of the probability.As shown in Fig. 1, the direction to the center of halo i is ni and the direction of the test photon is n.We assume each halo is located at a single redshift z i throughout the photon's crossing time (e.g. the halo size does not encompass cosmologically relevant distances), and take advantage of azimuthal symmetry around the ni direction.

Figure 1:
The trajectory of a photon ⃗ x(t) through a dark matter halo centered at χ i , ni on an observer's sky.In a coordinate system whose origin is at the halo center, the photon trajectory follows r(⃗ x(t)) ≡ r(t).Resonant conversion to a dark photon occurs when the photon crosses through a spherical shell at r res over a timescale ∆t res ; each trajectory has two crossings.We define the boundary of each halo by the virial radius r vir .
In the small angle approximation where n||n i , the minimum comoving distance between the halo center and the photon's path is (1 + z i ) r min = χ i θ, where χ i is the radial comoving distance to the halo center and θ is a small angle.The velocity |dr(t)/dt| = cos φ(t), where sin φ(t) = r min /r(t).
Using this geometry, the conversion probability due to halo i is separable into a radial part and an angular part.The latter one is a measure of how long the photon spends inside the resonance region ∼ ∆t res .We have: where we define (2.11) The Heaviside step function Θ(r res − r vir ) arises since we consider only photon to dark photon conversion happening inside the boundary of each halo, which we take here to be the virial radius r vir .The step function is normalized so that Θ(r res − r vir ) = 2 for r res < r vir and Θ(r res − r vir ) = 1 for r res = r vir to account for the fact that a photon crosses a resonance twice, going in and then out of a halo, except when it exactly grazes the edge of the virial radius.The effect of the sharp truncation at the boundary of the halo is significant for masses m γ ∼ 10 −13 eV that probe low densities.Note that the function u(n − ni |χ i , m i ) has the apparent singularity when θ = (1 + z i )r res /χ i , where the conversion probability blows up.However, the integral over the profile is finite and equal to (2.12) We close this section by noting that within the halo model, we make the rather drastic assumptions that each halo (dark matter and gas) is spherically symmetric, identical at each mass and redshift, and has properties independent of their formation history and local environment.These assumptions will fail for individual halos.However, we expect that quantities dependent on the statistical properties of the full distribution of halos, such as power spectra, will be well-approximated (see e.g.recent analyses such as [52]).We are primarily concerened with such statistical quantities in the following.

The photon to dark photon conversion monopole
The total probability for a photon to convert over its trajectory is the sum of contributions from each individual halo.In this section we focus on the sky-averaged probability and explain why this effect induces a new type of optical depth.In the next section we introduce the two-point function of optical depth fluctuations owing to the halos' shapes on the sky.
The overall conversion probability is where d 3 χ = χ 2 dχd 2 n and we have separated the resonant conversion probability P γ→A ′ into its radial and angular components that depend on each halo's characteristics.To obtain the expression above we also identified the number density of halos of mass m at redshift χ as where ⟨. ..⟩ denotes a sky-wide ensemble average and the delta functions are evaluated at each halo position and mass.The halo number density per volume per halo mass n(χ, m) is the halo mass function.We use the Tinker mass function throughout [53].
Integrating the angular profile u(n|χ, m) over the sky we implicitly weight the probability by the effective projected area of each halo.Evaluating the angular integral and simplifying we obtain: where u 00 = d 2 n u(n)/4π is the monopole of u(θ), and we changed the integration variable to redshift z from radial comoving distance χ.Within our assumption of instantaneous reionization, we impose a sharp upper limit on the integral over redshift at z reio .Photon to dark photon conversion manifests itself as a frequency dependent optical depth, encoding the removal of CMB photons along their path from recombination to CMB telescopes on Earth.The sky-averaged magnitude of this optical depth is the integrated probability in Eq. (3.3): where from here onward the bar notation stands for the projected, i.e. integrated over redshift, sky average.We introduce a new notation for the dimensionful average optical depth such that η depends only on the dark photon mass m 2 A ′ and the cosmology.It is useful to make this distinction because both ε and ω are parameters we vary later.
In Fig. 2, we show several examples of the differential optical depth along the line of sight for a range of dark photon masses.Notice in this plot that the light dark photons are only produced at low redshift.Light dark photons probe the outer-most regions of halos and the abrupt fall-off is due to the truncation of halos at the virial radius.Physically, this falloff would be broadened by the softer boundaries between halos and the intergalactic medium.In contrast, the heaviest dark photons probe regions near the core of halos where the gas density is highest.Since the gas profile is nearly flat near  The average dimensionful optical depth defined in Eq. (3.5) as a function of redshift, for a range of dark photon masses.To obtain the full optical depth, one needs to integrate over redshift and multiply by the unknown mixing parameter divided by the frequency of the photon, ε 2 /ω.Notice that the lighter masses produce the strongest signal at low redshift.Meanwhile, heavier dark photons probe redshifts all the way to reionization but require crossing large densities in order for the resonant conversion to take place.These are found within more massive halos whose number density is suppressed.
the core of halos, only the heaviest and rarest halos contribute to the optical depth of photons with m A ′ ≥ 10 −12 eV.In Fig. 3 we plot the total (e.g.integrated over the line of sight) dimensionful optical depth assuming both z reio = 6 and z reio = 10.The change in magnitude is less than the thickness of the blue curve, meaning that our model is insensitive to uncertainties related to when the end of reionization takes place.From now on, unless otherwise stated, we assume that reionization takes place instantaneously at z reio = 6.Notice that, unless ε ≪ 1 photon to dark photon conversion is in the optically thick regime.For example, τ = η ε 2 /ω = 1 for ∼ 100 GHz photons (near the peak of the CMB blackbody) at ε ∼ 10 −4 .This is a preliminary indicator that photon to dark photon conversion at low redshift can be a sensitive probe.The dark photon mass range over which there is a significant effect spans roughly one order of magnitude, peaking at m A ′ ∼ 6 × 10 −13 eV.
Baryonic feedback can affect the accessible dark photon mass range and strength of the photon to dark photon conversion monopole, in principle contributing a source of modeling uncertainty in the expected signal.We explore this modeling uncertainty in more detail in Appendix A by comparing the AGN gas model used here [50] with the results obtained by assuming that baryons trace the dark matter NFW density profile [51].We show that the two models give nearly identical results for low masses m A ′ where the resonant conversion condition is met near the outermost region of a halo.This can be understood as a consequence of the fact that the AGN gas profile and the NFW profile are only different near the core of halos where the NFW density increases without bound, while the gas density achieves a maximum; beyond the scale radius, the two density profiles are nearly the same.Also for this reason, the predictions for the two models differ most at large dark photon mass.
Here, details about the gas population in each halo can have a significant effect on the projected signal and resulting sensitivity, and it is important to incorporate this modeling uncertainty in the interpretation of our results below.As we will demonstrate in more detail below, the monopole signal is more strongly influenced by this modeling uncertainty than correlation functions.).Under our assumptions, resonant conversion occurs once reionization is completed and here we plot the case where z reio = 6 and where z reio = 10.The maximum change in the amplitude of η is of order 10 2 eV, less than the thickness of the blue line, and we conclude that the model is insensitive to the details of reionization over the range of masses we are probing.We therefore assume z reio = 6 in the remainder of this work.At low mass, η → 0 due to the constraint that conversion only happens within the virial radius in each halo.The upper bound is set by the shape of the gas profile and details about the halo model.The contribution due to gas in the Milky Way η MW (black) is defined in Eq. (3.6).

Contribution from the Milky Way
The gas halo surrounding the Milky Way is also a source of resonant photon to dark photon conversion.
In this section we model its contribution to the overall optical depth.The Earth's distance from the galactic center is well below the scale radius of the Milky Way halo, therefore in our model we assume that there is only an appreciable contribution to the optical depth monopole and not to optical depth anisotropies.We further assume that the Milky Way is an average spherical halo, with average AGN feedback and model it with the same gas profile of Eq. (2.5) computed for z MW = 0 and m MW from [54].In reality, the details of the gas distribution can affect both the monopole and anisotropies in the optical depth.Nevertheless, we can estimate the relative magnitude of the galactic versus extra-galactic conversion under these assumptions.
The expression for the optical depth is simply the contribution from a single halo with appropriate properties: where we discard the angular component by assuming the Earth is near the center of the halo, and ω is the frequency of the CMB photons today.Assuming a virial radius and virial mass for the Milky Way as found in [54], as well as the concentration-mass relation at redshift zero from [55], we calculate the conversion probability for the relevant range of dark photon masses.The result is plotted in Fig. 3. Given the assumptions made in this paper, the extra-galactic optical depth dominates the contribution due to conversion in the Milky Way for all dark photon masses considered.Since the Milky Way does not host an AGN, the gas profile assumed here may be too diffuse.In the extreme scenario where gas traces dark matter, the magnitude of the Milky Way contribution increases monotonically for higher dark photon masses.The increase is cutoff at the core when gas no longer traces dark matter.The high dark photon mass regime is where the extra-galactic contribution to the optical depth becomes small, and the Milky Way therefore introduces additional modeling uncertainties at the upper end of dark photon masses we consider.Note that such effects will only increase the reach in sensitivity to conversion, making the neglect of contributions from the Milky Way a conservative assumption.

Patchy dark screening
When photon to dark photon conversion occurs in non-linear structure, the associated optical depth for conversion is strongly anisotropic on the sky.These anisotropies in the optical depth serve as a screen of varying opacity through which the CMB must propagate on the way from decoupling to our telescopes here on Earth.In addition to the sky-averaged suppression in the intensity of CMB photons, new temperature and polarization anisotropies are introduced due to the different conversion probability of CMB photons to dark photons across different lines of sight.
The analogous effect in the standard cosmological model is the 'screening' of CMB anisotropies due to the Thomson scattering of CMB photons by free electrons.For Thomson scattering, anisotropies in the optical depth couple to anisotropies in the CMB temperature and polarization.This is known as 'patchy screening' of the CMB.Notably, there is no coupling of optical depth anisotropies to the CMB temperature monopole -Thomson scattering doesn't change the energy of photons and for every photon scattered out of the line of sight, another is scattered into the line of sight.Therefore, patchy screening of the CMB is always a small effect, i.e. second order in perturbations of the CMB and optical depth anisotropies.Although small, the detection of patchy screening during [16] and after [56] reionization is within the reach of future CMB experiments.
Photon to dark photon conversion produces patchy dark screening, which has two crucial differences to Thomson screening: conversion does not preserve the blackbody spectrum and conversion only removes photons from the line of sight.Therefore, patchy dark screening couples the CMB monopole to fluctuations in the optical depth, and is a 1st order effect (in anisotropies).Since the CMB monopole is ∼ 10 4 times larger than the temperature anisotropies, patchy dark screening is far stronger at fixed optical depth than Thomson screening.Furthermore, the characteristic frequency-dependence of patchy dark screening can be used to separate it from the primary CMB and astrophysical foregrounds.In the remainder of this section and the next section, we derive various correlation functions that will be used to forecast the sensitivity of CMB experiments to the kinetic mixing parameter and dark photon mass, assuming the frequency dependent dark screening anisotropies can be separated from the primary CMB anisotropies.A detailed discussion of how well this separation can be performed is presented in Section 6.2.

Anisotropic screening
Before computing CMB correlators, we must first describe anisotropies in the photon to dark photon conversion optical depth, the anisotropies in the Thomson optical depth, and the cross-correlation between these two fields.The halo model for large-scale structure is a useful tool for these computations, since it is straightforward to populate dark matter halos with the electron density (for Thomson screening) and the conversion probability (for dark screening).
The standard optical depth to reionization is the integrated electron density along the line of sight: where σ T is the Thomson cross-section and a is the local scale factor.The inhomogeneous matter distribution introduces spatial fluctuations in the electron number density n e (χ, n).This can be measured via the directional dependence they introduce on the optical depth field.The local perturbations in the electron density induce small fluctuations δτ Th (n) in the optical depth profile where τ Th is the standard cosmological parameter.Once again we work under the assumptions of the halo model where the integral above can be written as an average over halos whose number density depends on mass and redshift, where the electron number density fluctuations are related to the gas density profile.In the Limber approximation, the multipole expansion of the electron density from one halo centred on the north celestial pole (such that the azimuthal angular momentum m = 0) is [57]: where the ρ gas is the Battaglia et al.AGN gas density profile [50] and ρ e = n e m p .The Limber approximation works best for ℓ ≫ 1, which is the regime we are in.We now introduce the optical depth for photon to dark photon conversion.Consider small anisotropic perturbations to the average extra-galactic optical depth τ : Note that with this convention δτ (n) is independent of ε and ω.The next step is to compute the two-point angular power spectra for these fields.Projecting onto spherical harmonics, we define the power spectrum In the context of the halo model, the power spectrum can be described by a sum of an intra-halo (1-halo) and an inter-halo (2-halo) contributions: The quantity C τ τ ℓ is the angular power spectrum of the optical depth due to photon to dark photon conversion (which depends on ε and ω).The terms on the right hand side of Eq. (4.6) are computed in detail in Appendix B. The final expressions are: (4.7)The quantity j ℓ (kχ(z)) is the spherical Bessel function and b(z, m) is the linear halo bias.P lin (k) is the linear matter power spectrum calculated using CAMB [58].All quantities are computed for fixed cosmology.Throughout this work we used the following set of parameters: dark matter density Ω c h 2 = 0.12, baryon density Ω b h 2 = 0.022, Hubble constant H 0 = 67.3km s −1 Mpc −1 , scalar spectral index n s = 0.96, scalar amplitude A s = 2.2 × 10 −9 , and optical depth to reionization τ Th = 0.06.
In Fig. 4 we plot the optical depth power spectrum of δτ fluctuations for a range of dark photon masses, assuming z reio = 6.In the left panel, we show the relative importance of the 1-halo term compared to the total.The 1-halo term dominates at high ℓ, while the 2-halo term dictates the shape and amplitude on large scales ℓ ≲ 1000.The 1-halo term is larger at the upper end of the dark photon mass range we consider.This is consistent with the fact that conversion happens near the halo core in this regime.In the right panel we show the total power spectrum C δτ δτ ℓ for 4 choices of m A ′ that span the parameter space we probe for.Any changes can be attributed to the radius of the resonance surface at that mass given the gas density profile ρ gas as well as details about the halo model, for e.g. the population of halos n(z, m), which also dictates the bias function b(z, m).The total magnitude of the power spectrum, C τ τ ℓ , depends strongly on the dark photon mass through η.Hence, it will be maximized by masses near the peak in the monopole τ (m A ′ ), which was depicted in Fig. 3.
Repeating the computation presented in Appendix B for the dark screening case, we find the 1-halo and 2-halo contributions to the Thomson screening auto-power spectrum.The full expression  , which is independent of ε and ω, for a range of dark photon masses that span the parameter space we are probing.In the left panel, we compare the 1-halo term with the sum in Eq. (4.6) to show that most power on large scales is due to the 2-halo term, while on small scales the 1-halo term dominates.At large mass, when the transition from photon to dark photon occurs nearer to the halo core where gas densities are largest, the two terms become comparable and the spectrum is scale invariant.This feature also affects the signal amplitude hierarchy: although at large m A ′ > 10 −12 eV the monopole is subdominant (e.g.Fig. 3), the additional structure on small scales causes C τ τ ℓ in this regime to be equivalent in terms of constraining power for ε to spectra for m A ′ < 10 −12 eV.This will be relevant in Section 6 where we present our forecasts.is The quantities have the same meaning as in expressions (4.7) above.Anisotropies in the Thomson screening and dark screening optical depth fields trace similar matter density profiles over the sky and across the line of sight.It is therefore expected that the fluctuations in either field are correlated.The two-point function τ τ Th ∼ ε 2 has the following angular power spectrum: where all quantities are real.The Thomson auto-power spectrum is shown in Fig. 5, alongside examples for the power spectrum of the dark screening optical depth and their cross-correlation for fixed m A ′ and ε ( C τ Th τ Th ℓ is independent of ε).The power spectra involving dark screening scale differently with ε, so that for ε around 7 ℓ .The magnitude of C τ τ ℓ varies over the range of dark photon masses we study such that it peaks around m A ′ ≃ 10 −12 eV and falls abruptly for mass values towards both ends of the interval./C τ τ ℓ (pink ), which will be useful for the sensitivity forecast in Section 6.

Patchy dark screening of the CMB
We now describe in detail how the optical depth anisotropies of the previous section manifest as anisotropic spectral distortions in the CMB temperature and polarization.The largest effect arises from patchy dark screening of the CMB temperature monopole.In addition to this, there are new spectral anisotropies from the screening of temperature and polarization anisotropies.An additional novel signature of patchy dark screening is the production of B-modes (curl) from pure E-mode (divergence) polarization anisotropies, as occurs in any scenario with screening [59].This leads to statistically anisotropic correlations between the temperature, E modes, and B modes dependent on the specific realization of the anisotropic optical depth and CMB anisotropies in our Universe.
In analogy with patchy screening from Thomson scattering, the photon to dark photon optical depth fluctuations suppress the CMB temperature fluctuations and polarization Stokes parameters.The combined effect of dark and Thomson screening on the observed temperature and polarization are where (4.12) The label obs stands for observed anisotropies, Sc for screening from Thomson scattering, and dSc for dark screening.We have explicitly isolated the dependence on the CMB blackbody temperature monopole T = 2.725K, where in our notation the sky average of T (n) is zero while the sky average of T (obs) is not.As our first approximation, we assume that T (n) and (Q ± iU )(n) are the lensed CMB temperature and polarization anisotropies.This neglects the lensing of the screened CMB along the line of sight, which is a small higher order effect.We neglect extra-galactic foregrounds, such as the cosmic infrared background, point sources, and Sunyaev Zel'dovich effects.We assume that no significant kinetic mixing happened between recombination and the end of reionization.Finally, we work in the limit where τ (n, ω), τ Th (n) ≪ 1.
The total spectral distortion to the blackbody spectrum due to photons' conversion into dark photons is given by where P γ→A ′ is the overall conversion probability and B 0 (ω, T ) represents the intensity of the theoretical Planckian spectrum.In natural units this is: The dependence on frequency ω in Eq. (4.12) is relative to the blackbody spectrum of the CMB in units of temperature, not intensity (see, e.g. also [60]).To convert, we take the leading order expansion in temperature fluctuations δT (n) of Eq. (4.13): Therefore, in CMB temperature units, the temperature fluctuation function of frequency is where x = ω/ T in natural units.A factor of appears when converting the frequency dependence in the absorption optical depth to the frequency dependence in the standard temperature units, and the full frequency dependence of the dark screening signal in units of temperature is ζ(ω)/ω.Next, we decompose the lensed CMB temperature and optical depths into spherical harmonics and the Stokes parameters into spin-2 spherical harmonics: where E ℓm and B ℓm are the moments of the E-and B-mode polarization anisotropies.Note that the latter is induced only by lensing in the absence of primordial gravitational waves and is therefore substantially smaller than E-mode polarization.We neglect the dark screening of lensing B-modes below.
Given multi-frequency observations of the CMB, it is possible to separate the blackbody and dark screening components of the temperature and polarization anisotropies.Assuming perfect separation of the blackbody and patchy dark screening components (we discuss the scenario where this separation is imperfect in Section 6.2), we construct correlation functions between Thomson-screened and darkscreened temperature and polarization anisotropies, that is, T Sc and T dSc , respectively.Note that there will be both statistically isotropic and statistically anisotropic components of the correlation functions.

Dark screening auto-correlation functions
As a warm up, in the Standard Model (ε = 0 and hence T dSc = 0), the anisotropies in optical depth τ Th lead to temperature anisotropies, which are captured in the auto-correlation function for the Thomson-screened blackbody anisotropies: where we have defined The statistically isotropic component of the correlator is [59] is the lensed primary CMB temperature power spectrum, C τ Th τ Th ℓ is the Thomson optical depth power spectrum, and N T Sc T Sc ℓ encompasses all other contributions to the blackbody CMB such as instrumental noise, foregrounds, etc.The statistically anisotropic component of the correlator is induced by the particular realization of patchy Thomson screening in our Universe.We have explicitly kept this term, which would vanish for a full ensemble average over all fields.As we describe in detail below, this statistical anisotropy can be used to reconstruct the anisotropies of optical depth τ Th , and induces various three-point correlation functions.
Similarly, anisotropies in the dark screening optical depth τ (n, ω) lead to anisotropies in T dSc (ω), which is captured by the auto-correlation function of the dark screening component of the observed CMB temperature: ) where C τ τ ℓ (ω, ω ′ ) is the dark screening optical depth, for which we have retained the frequency dependence with respect to the blackbody.The statistically isotropic contribution to the correlator is The contribution from instrumental noise and foregrounds is N T dSc T dSc ℓ (ω, ω ′ ).Because the CMB anisotropies are so small in comparison to the monopole, the first term completely dominates the statistically isotropic contribution to the correlator.As a result, measuring C T dSc T dSc ℓ can be extremely sensitive to photon to dark photon conversion.However, such an auto-correlation function is proportional to the small kinetic mixing parameter ε 4 , which limits the reach of measurements of the dark-screened CMB power spectrum.We forecast the reach of such an analysis in Section 6.3.
Similar dark-screened auto-correlations can be computed for the T E, EE and BB CMB spectra.These are listed in Appendix C. In determining the sensitivity of CMB experiments to ε in Section 6, we use all isotropic components of the two-point auto-correlators.However, the temperature autocorrelation provides the best sensitivity because of the coupling of τ to T .
The statistically anisotropic component of the dark screening auto-correlation is proportional to the un-screened temperature anisotropies, and we explore in the next section how this can be used to search for photon to dark photon conversion.We can also construct cross-correlations functions in the form of ⟨T dSc ℓm (ω)T Sc ℓ ′ m ′ ⟩, the discussion of which we also leave to the next section.

Cross-correlating dark screening
In this section, we study the cross-correlation between observables that contain the dark screening optical depth τ (n, ω) and those that do not.The essential qualitative understanding that motivates the construction of these correlation functions is the following: Dark screening occurs in halos, and is therefore correlated with observables (within the Standard Model) that are sensitive to either the halos' locations or their electron density distributions.As we showed, both Thomson screening and photon to dark photon conversion depend on the electron density profile.Therefore, there is a nonzero cross-correlation between the Thomson-screened and dark-screened temperature and polarization anisotropies.In the following, we discuss two ways of combining these maps: we consider first crosscorrelation between the Thomson-screened and dark-screened CMB, and then their correlation with LSS tracers.We compute the relevant two-and three-point correlation functions for each method and identify the ones that are most sensitive to dark screening.

Two-point cross-correlation
Based on the discussions in the last section, an obvious candidate for the correlation function we want is: where C τ τ Th ℓ (ω) was defined in Eq. (4.9) and Note that for this correlator there are statistically anisotropic contributions proportional to both the dark screening optical depth and the un-screened temperature anisotropies.The cross-correlation between the Thomson-screened and dark-screened temperature scales as τ ∼ ε 2 .This is more favorable than the ε 4 scaling found in Eq. (4.22).However, unlike for Eq.(4.22),only the second statistically anisotropic term depends on the temperature monopole T (Thomson screening doesn't couple to the temperature monopole).The consequence is that this cross-correlation is less competitive than the dark-screened auto-correlation at fixed noise.Before moving on, let's turn to polarization.E and B modes are defined with the relation in Eq. (4.18).Combined with the assumption that the Thomson-screened temperature anisotropies can be separated from the dark-screened anisotropies through the frequency dependence, we can construct a variety of two-point correlation functions -12 in total!These can be found in Appendix C. Unlike for temperature auto-and cross-correlations, some correlators involving polarization vanish in the absence of screening, and only receive statistically anisotropic contributions.These include and where C EE ℓ is the lensed primary CMB E-mode power spectrum and These correlations are also potential sensitive probes of patchy dark screening, as demonstrated in greater detail below.

Correlating patchy dark screening with LSS
Expanding our focus beyond CMB observables, photon to dark photon conversion happens inside halos, and therefore the patchy dark screening signal is correlated with the LSS.It is most natural to look for correlations between patchy dark screening and tracers of LSS, such as galaxy redshift surveys.Given a tracer and various model assumptions, one can build a template for patchy dark screening, which we will denote by τ (ε 0 ).In this section we build the intuition of how to use such a template to improve our sensitivity to ε.As one example, the template can be built in the following way: where g ℓm are the moments of the galaxy overdensity field and the vector notation denotes the redshift information; C τ τ ℓ is the model photon to dark photon optical depth power spectrum, C gg ℓ is the redshift × redshift galaxy overdensity covariance matrix, and C gτ ℓ is a vector of the dark photon optical depth × galaxy overdensity cross-spectra at each redshift.We have explicitly indicated that the model power spectra involving the patchy dark screening optical depth depends on the fiducial choice ε 0 for the kinetic mixing parameter.Note that the template defined in Eq. (5.6) can be improved by going beyond this simple linear filter, for example using machine learning techniques as in Ref [61].
The largest contribution to the cross-correlation of the template with the patchy dark screening component of the CMB is statistically isotropic and given by where C τ τ g ℓ is the cross-power spectrum between the template and real photon to dark photon optical depth.Importantly, C τ τ g ℓ (ε, ε 0 , ω) ∝ ε 2 , and therefore this quantity scales more favorably with ε in the small-ε limit than the monopole contribution to the temperature auto-spectrum (which scales ∝ ε 4 ).In addition to the more favorable scaling with the mixing parameter, cross-correlation with a template can be beneficial for mitigating systematic effects and galactic or extra-galactic foregrounds in the observed CMB.We forecast the reach of such analysis in Section 6.4.

Reconstruction
Coming back to the CMB, the discussion in the previous section suggests searching for correlation functions that would allow us to reconstruct the map τ Th (n) from CMB observables.τ Th (n), depending on the same electron density distribution in the Universe as τ (n, ω), would be correlated with T dSc , just like τ g .The statistically anisotropic components of the two-point correlation functions in Eq. (4.19) can be used to construct quadratic estimators for patchy Thomson screening optical depth τ Th (n).Similarly, patchy dark screening, as well as the un-screened primary CMB temperature and polarization anisotropies can be reconstructed e.g. from equations (4.22), (5.1), (5.4).An exhaustive list of quadratic estimators is presented in Appendix C. Similar quadratic estimators are used to measure weak lensing of the CMB (see e.g.[9] for a review), and are employed in a wide variety of other contexts in CMB science, in particular for the reconstruction of patchy Thomson screening during reionization [16] and kinetic Sunyaev Zel'dovich velocity reconstruction [10][11][12][13][14].
Starting again with a Standard Model example, the reionization optical depth can be found from Eq. (4.19) as: where (5.9) The weights and pre-factor are chosen such that the estimator is unbiased if the input power spectra provide an accurate model, i.e.
as well as has minimum variance when all fields in the problem are Gaussian: (5.11) The pre-factor N is the noise on the reconstruction.Similarly, one can reconstruct τ Th from measurements of polarization, notably from E Sc B Sc .The relevant quadratic estimators are presented in Appendix C. Estimators for the Thomson optical depth found there are equivalent those presented in Ref. [16].The two-point function that would result from cross-correlating the dark-screened CMB temperature with the Thomson optical depth map, i.e. ⟨τ Th T dSc ⟩, is another means to construct the correlation function in Eq. (5.7).
More quadratic estimators can be built from the dark screening anisotropies T dSc .The leading term that contains the dark screening optical depth is: where We used the fact that the product of the screened auto-spectra is larger than the product of the screened cross-spectra in the second line.
It is important to note that in constructing the weights G τ ;T Sc T dSc ℓℓ ′ L and the pre-factor N τ ;T Sc T dSc L that the input power spectra come from a theoretical model for the un-screened CMB temperature power spectrum C T T ℓ as well as contributions from noise and foregrounds to both can be checked against the measured power spectra of these maps.However, since we cannot directly measure C T T ℓ , there is inevitably some residual model uncertainty.This manifests as a bias on the reconstruction: (5.15) When the model and actual power spectra are identical, the bias factor is unity.Importantly, it is possible to measure this bias by comparing e.g.τ qe LM to τ g LM .In principle, the bias can be measured without cosmic variance since the comparison is done at the level of the modes and not the power spectra; this is an example of 'sample variance cancellation' (e.g.[62]).This procedure is elucidated in greater detail in Section 6.
As a third example, a quadratic estimator for the un-screened temperature anisotropies is where (5.17) The weights and reconstruction noise in this case depend on a theoretical model for C τ τ Th ℓ , an object we have no prior knowledge of and hope to search for.This implies that the reconstruction of the un-screened temperature anisotropies will be significantly biased.However, as described above, one can measure the bias by comparing to a template for T LM , which at least on large angular scales, can be provided by T Sc .That is, a correlation function ⟨ T T Sc ⟩, can also be used to search for anisotropic dark screening.

Three-point correlation functions
The discussion in the previous section suggests that one should correlate maps reconstructed from CMB observables, such as τ Th , τ and T , with CMB observables (e.g.T Sc or T dSc ) or itself.In terms of CMB observables, these correlation functions will be three-point or four-point correlation functions.In other words, the statistically anisotropic contributions to the two-point correlation functions in the previous section imply that there are many non-vanishing three-point correlation functions even in the case where the temperature, polarization, and optical depth fields are Gaussian.
For example, the correlation functions both come from the three-point correlation function ⟨T dSc T Sc T Sc ⟩.Therefore, rather than working with the more intuitive two-point correlation functions involving reconstructed maps, we can forecast the sensitivity directly using three-point functions; we present this forecast in Section 6.5.Three-point functions are described by the angle-average bispectrum, defined as There are many bispectra to consider between temperature, polarization, and templates for the optical depth.Here we focus on the largest bispectra that scale like ε 2 , since these will be most sensitive to photon to dark photon conversion.The most relevant bispectra involving only CMB temperature and polarization are (see Eq. (5.1) and Eq.(5.3)) (5.21) Note that both are proportional to the CMB monopole.The most relevant bispectra that scale like ε 2 and involve both the CMB and an optical depth template are (5.23) These bispectra, similar to the cross-correlation with LSS in Eq. (5.7), are proportional to the CMB monopole T , while at the same time, scales as ε 2 .This, as we will demonstrate in more detail in Section 6, makes these bispectra almost as sensitive as the two-point functions ⟨T dSc T dSc ⟩ (scaling as ε 4 ), and more sensitive than ⟨T dSc T Sc ⟩ (proportional to C T T ℓ ).

Forecast
In this section we forecast the projected sensitivity of several CMB experiments to the mixing parameter ε over a range of fixed values for the dark photon mass m A ′ .We consider each of the five techniques mentioned in the introduction, and identify the most promising observables for each experimental configuration.These relevant correlation functions are summarized in Table 1.We first compute the constraint on ε from the monopole spectral distortion using COBE/FIRAS.We then describe our forecast assumptions, and compute the reach of existing and future CMB anisotropy experiments.
Table 1: The scaling of various correlation functions with the small parameters: kinetic mixing ε, primary CMB power C T T ℓ , primary E-mode polarization power C EE ℓ and optical depth of Thomson screening τ Th .The correlation functions ⟨T Sc T Sc ⟩ and ⟨T Sc T dSc ⟩ are also shown for comparison.

FIRAS bounds
First we look at the CMB monopole constraint given by the COBE satellite.This method has been used in the past to forecast the constraint on ε from different models for the distribution of ionized gas since recombination [31,34,44].
The CMB monopole was measured by the FIRAS instrument on the COBE satellite and was discovered to have a near perfect blackbody spectrum [45].The data consists of 43 measurements of the sky-averaged temperature over frequencies in the range ω = 68.05− 639.46 GHz2 .This gives a best fit blackbody temperature of T = 2.725 ± 0.002 K, with residuals of order roughly 10 −4 below the peak intensity and 1σ uncertainties of the same magnitude.This remarkable precision already gives a constraint on the amplitude of the conversion probability for CMB photons of order the uncertainty P γ→A ′ ≲ 10 −4 .In this section we improve this bound by considering the full available frequency spectrum.To constrain our model we use the method described in [31].
Assuming an isotropic conversion of photons into dark photons, the CMB blackbody spectrum is distorted according to where B 0 (ω, T ) is the theoretical blackbody spectrum defined in Eq. (4.14) and τtot , the total dark screening monopole including the galactic component, implicitly depends on ω and ε.The reduced χ 2 estimator is an average over all available frequency channels of the difference between the measured data and the expected signal B exp in each frequency bin: This estimator is minimized by some value of T at each point in the plane spanned by ε and m A ′ .In Fig. 6 we show the 95% and 99% confidence limit contours in this parameter space, for the χ 2 of a distribution with 42 degrees of freedom.The exclusion regions are similar in both cases, indicating that the regression method is robust, i.e. the χ 2 changes rapidly with temperature around the minimum.Over the accessible range of dark photon masses, the constraint on ε reaches up to ≃ 10 −6 , roughly two orders of magnitude better than the naive limit set by the error bar on the blackbody temperature.In Appendix A we explore the model uncertainty implicit in our constraint, which arises primarily from our lack of knowledge of the gas profile.Comparing the constraints obtained in our fiducial model for gas profiles with those obtained for a model in which gas perfectly tracks dark matter (e.g. a model without any baryonic feedback), we conclude that the constraints at low dark photon mass are robust.At high dark photon mass, the constraint is strengthened in a model where baryons perfectly track dark matter.We can extrapolate that if our fiducial gas model has too little baryonic feedback, the constraints could further weaken at high dark photon mass.Incorporating information from measurements of the Sunyaev Zel'dovich effect (e.g.[63]), the dispersion measure of fast radio bursts (e.g.[64,65])), 21 cm intensity mapping (e.g.[66]), and other tracers of baryons will be helpful in mitigating this modeling uncertainty and will be explored in future work.[45] data.We see that the uncertainty in the measurement of T gives a constraint on the magnitude of the dark-screened temperature monopole T dSc .Here we used the total dark screening optical depth due to both galactic and extra-galactic contributions τ +τ MW .This constraint provides an upper bound on ε of at most 10 −6 .

Forecast assumptions for CMB anisotropy experiments
In Section 4.2, we assumed that the frequency-dependent anisotropies due to resonant conversion (T dSc , E dSc , B dSc ) could be perfectly separated from the blackbody Thomson-screened anisotropies (T Sc , E Sc , B Sc ).Here, we explore the degree to which this separation can be made in the presence of instrumental noise and measurements in only a small number of frequency channels.We estimate the residual noise on the dark-screened and Thomson-screened maps achievable with existing and future CMB experiments, which is used in the following forecasts based on two-and three-point correlation functions.
We consider three different CMB experiments: the combination of the Low Frequency Instrument (LFI) [67] and High Frequency Instrument (HFI) [68] on the Planck satellite, CMB Stage-4 [2, 69] and CMB-HD [3,70].In the context of our analysis, a CMB experiment is characterized by the sensitivity and resolution at a set of measured frequencies.Throughout, we assume Gaussian beam and white uncorrelated noise for all instruments, as well as full-sky coverage and no foregrounds.
Before proceeding, it is important to comment on the potential impact of CMB foregrounds.There are a variety of galactic foregrounds (see e.g.[71] for an overview) that fall with frequency, with a similar power law in the power spectra to the patchy dark screening signal (falling as ω −2 − ω −3 ) including: synchrotron, free-free, and spinning dust.These are strongest in the galactic plane, and their influence can be mitigated by masking the most contaminated regions of the sky, or by incorporating information about the morphology of the signal.Nevertheless, these foregrounds can in principle add significant extra power at low-frequencies and on large angular scales.Extra-galactic radio point sources, which are dominated by synchrotron emission, are also a potentially important foreground to consider.Resolved point sources can again be dealt with by masking, however the unresolved point sources can add power at low frequencies and on small angular scales.This signal is also correlated with other tracers of LSS, limiting the power of cross-correlations to mitigate foregrounds.The extent to which galactic and extra-galactic foregrounds degrade the forecast we present below requires a detailed analysis, which we defer to future work.
The instrumental noise considered throughout our analysis is modeled as: Here, ∆ T [µK rad] represents the sensitivity in temperature, while θ FWHM [rad] is the full width at half maximum of our assumed Gaussian beam, which characterizes the resolution of the instrument.The sensitivity and resolution vary with frequency.Furthermore, ground-based measurements are subject to atmospheric contamination on large angular scales.To account for this effect, in the analysis for both CMB-S4 and CMB-HD we include the additional 'red noise' term in the second bracket with α knee = −3 and ℓ knee = 100 in all frequency bins.This contribution diverges at low-ℓ, and becomes increasingly irrelevant for ℓ > ℓ knee .The values we choose to represent each experiment are shown in Table 2.
Frequency Table 2: Noise parameters for Fisher forecasts.We model the noise covariance as in Eq. ( 6.3) where in the case of ground-based CMB-S4 and CMB-HD we include the red noise term with parameters α knee = −3 and ℓ knee = 100 in all frequency bins.
Later, when we compute the Fisher information, we will need to estimate the noise covariance for the dark screening two-point functions of Section 4.3.These power spectra have an intrinsic inverse frequency squared dependence C X dSc X dSc ℓ ∝ ε 4 /ω 2 .It is possible to disentangle this signal from the measured CMB by cross-correlating measurements across multiple frequency channels.We do this by applying a harmonic Internal Linear Combination (ILC) algorithm [72] which we now describe.
Recall that the expected isotropic measured CMB signal consists of the primary CMB which is screened by the inhomogeneous field τ plus instrumental noise.For example, in the temperature case we have from equations (4.11), (4.12), (4.23): The instrumental noise term is frequency dependent as displayed in Table 2, where the i and j labels above denote each available channel.Let us re-write the dark screening two-point function term in the frequency explicit form C τ τ ℓ (ω i = 1, ω j = 1)/ω i ω j .To build the covariance matrix for all temperature measurements we simply choose a frequency for the signal, ω 0 , and measure all entries in reference to it.Overall we find: The frequency-dependent matrix has entries ) and e = (1, 1, . . ., 1).The ILC method consists of weighting each matrix element appropriately in order to minimize a frequency-independent residual.This is constructed like where the weights w ℓ satisfy

.7)
In our analysis, we use a fixed baseline frequency ω 0 = 30 GHz.For the C XX ℓ terms of the primary CMB we compute the lensed temperature and polarization power spectra using CAMB [58], using the cosmological parameters listed in Section 4. Fig. 7 shows the weights computed for CMB-S4 specifications.On large scales, the 93 and 145 GHz maps, which are the lowest noise, are used to subtract off the blackbody CMB, with the other channels weighted inversely with frequency.
A similar computation can be used to find the weights and residual noise for polarization.The ILC-cleaned spectrum of the temperature map includes the frequency-dependent term, computed at the baseline frequency ω 0 , plus the noise residual: Henceforth, the .notation will refer to an ILC-cleaned map, or its associated residual defined as in Eq. (6.6).Note the factor of ζ(ω) defined in Eq. (4.17).The CMB peaks in intensity around 160 GHz.The factor of ζ(ω) approaches unity for frequencies below the CMB peak, and vanishes for frequencies ω ≫ 160 GHz.Since a ω −1 scaling weighs small frequencies more strongly, the additional ζ(ω) factor does not affect the analysis significantly.Not only does Ñ T dSc T dSc ℓ do a good job of removing the primary CMB signal, it also gives a lower noise amplitude at high ℓ compared to the noisiest frequency channels.This can be seen by comparing the solid purple line in Fig. 8 depicting Ñ T dSc T dSc ℓ with the dotted lines that show N T T ℓ (ω) for CMB-S4.In the limit of infinitely many frequency measurements, one could in principle fully isolate the frequency-dependent signal from the blackbody component.
Finally, to compute the ILC for the strictly blackbody signal of the Thomson-screened CMB of Eq. (4.19), we set Ω −1 ij = 1.The residual Ñ T Sc T Sc ℓ in this case is a linear combination of the instrumental noise spectra, once again weighted heavily by the middle frequencies with lowest magnitude, but this time all weights are positive.Again, if we had access to infinitely many frequency channels, the ILC would perfectly recover the Thomson-screened CMB. Figure 7: Illustration of the weight functions w ℓ as defined in Eq. (6.7) for the frequencies of CMB-S4 with noise parameters defined in Table 2. Notice that the dominant frequency is the ω = 93 GHz channel, but this changes at higher ℓ.This is because the ILC favors the N ℓ (ω) with the lowest magnitude at a given ℓ.  in green is comparable in magnitude to the orange curve around a value ε ∼ 8 × 10 −11 .In the case of EE, the switch happens around ε ∼ 7 × 10 −7 .The lensed primary CMB spectra (blue) and the ILC leftover noise (purple) computed at the baseline ω 0 = 30 GHz are shown for comparison.The dotted lines represent the instrumental noise N ℓ (ω) for each channel in CMB-S4, as defined in Table 2.The colors of the dotted lines are the same as their corresponding weights' in Fig. 7.

CMB auto-correlation
The first method we look at considers the cross-correlation of all dark-screened CMB auto-power spectra.The covariance matrix takes the following form: where each noise term contained in Cℓ is the ILC residual for the particular T , E or B measurement.The Fisher matrix is defined as: where (C ℓ ) −1 is the matrix inverse of C ℓ evaluated at a fiducial value for parameters i and j; ∂ i C ℓ is the derivative of the covariance matrix with respect to parameter i evaluated at the fiducial parameters i and j.The observed sky-fraction is assumed to be f sky = 0.7 for Planck and f sky = 0.5 for CMB-S4 and CMB-HD.We are interested in expanding the reach for the mixing parameter ε at a variety of values for m A ′ .Since in this case we only have access to dark screening power spectra, we are constraining ε 4 .The Fisher 'matrix' has a single entry given by F ε 4 ε 4 .To find the exclusion region, we compute the the 1-sigma constraint on ε 4 , given a fiducial value ε 4 = 0: The Fisher matrix simplifies in this case to: .11)and the 1-sigma sensitivity is roughly given by σ ε 4 ∼ 1/(S/N ) ε 4 =0 .
To obtain the variance on ε, we assume a Gaussian posterior for the probability distribution of the positive real-valued ε 4 .The sensitivity on ε is then well approximated by: . (6.12) Notice from the expression above that the value of σ ε improves as ∼ ℓ −1/4 with the number of modes with significant S/N.This feature is general, regardless of the observable we use in the forecast.In terms of signal, most of the sensitivity of the CMB auto-correlation is due to the C T dSc T dSc ℓ ∼ T 2 term.Finally, the shape of the sensitivity as a function of mass m A ′ will trace the optical depth monopole as η−1/2 .This sensitivity is shown in Fig. 9 for Planck and CMB S4.As explained above, the sensitivity is bounded at low dark photon mass from imposing a hard cutoff at the virial radius in each halo, while at high mass it is given by the fact that the halo mass function n(z, m) falls to zero.The boundaries of the contours are also sensitive to our assumptions about the lower bound on the mass of halos described by the AGN gas profiles.We comment on the impact of this modeling uncertainty in Appendix A. Overall, the sensitivity is superior to the FIRAS bound due to the dependence on the CMB monopole T .

CMB cross-correlation with a conversion template
Cross-correlating the measured CMB with other probes increases the forecasted sensitivity to ε.To investigate the degree of improvement that could eventually be possible, we assume in this section that a perfect template for patchy screening occurring at z < 2 is available.Such a template could be created from a massive galaxy survey and a detailed model for the relation between galaxy density and ionized gas density, as described in Section 5.2.The choice of z < 2 is motivated by the redshift range that will be covered by near-term surveys.When referring to the template, we will write τ and its amplitude will depend on a fiducial choice of coupling strength denoted by ε 2 0 .In the presence of a dark photon, the measured dark-screened CMB depends on an unknown 'true' τ ∝ ε 2 .One can write the first order in δτ contributions to the two-point function between the template τ and the dark-screened temperature T dSc from Eq. (4.12) as: Notice that another advantage of this method is that it allows us to be sensitive to ε 2 directly, which is the same power of the coupling that appears in the dark screening optical depth monopole (in Section 6.1).There is no statistically isotropic component of the cross-correlation between polarization and the template, and so we do not discuss them here.Note that for the template power spectrum, anisotropies are calculated up to z = 2 but the monopole τ corresponds to the full contribution up to reionization, i.e. τ = τ .Breaking up the two-point function into an integral over redshift, one can see that for a perfect template.For an imperfect template, this result would include a (scale-dependent) correlation coefficient describing the imperfect overlap of the template with the actual dark screening optical depth.
The covariance matrix assuming a perfect template is the following: With a fiducial ε → 0, the Fisher estimator in this case simplifies to: and the uncertainty on the coupling constant is This sensitivity is plotted in Fig. 9 for the various CMB experiments we consider.Notice here that the sensitivity contour corresponding to this estimator is improved by around one order of magnitude compared to the dark-screened CMB-only result.In short, this is due to the ε 2 scaling instead of ε 4 that brings about a more favorable scaling ∼ ℓ −1/2 with the number of modes that are measured with appreciable S/N.

Correlations with Thomson screening, the bispectrum and reconstruction
Recall from Section 4.3 that there are statistically anisotropic contributions to the two-point function between the CMB temperature and polarization anisotropies.These statistically anisotropic contributions form the basis of the reconstruction of the un-screened CMB temperature and polarization anisotropies as well as the Thomson and photon to dark photon optical depth introduced in Section 5.3 and enumerated in Appendix C.These statistical anisotropies also imply the existence of the three-point functions introduced in Section 5.4.In this section, we explore how this non-Gaussian information can be used to search for dark photons.We first consider the T dSc T Sc T Sc bispectrum, which is proportional to C τ τ Th ℓ (ω).Factoring out the dependence on ε: we see that forecasting the limits on ε 2 is a straightforward exercise in estimating the amplitude of this bispectrum, a problem whose optimal solution is already known from studies of primordial non-Gaussianity (for an overview, see e.g.[73]).The simplest bispectrum estimator is where the resulting constraint is related to the estimator variance by Note that just as for the quadratic estimators discussed above, the weights in the bispectrum estimator are constructed from models for C T T ℓ and C τ τ Th ℓ (ω, ε = 1).An alternative starting point would be to first use T dSc and T Sc to reconstruct T (see Eq. (5.16)) and then correlate this with T Sc .If we assume that the model for C τ τ Th ℓ is known up to the value of ε 2 then the reconstruction of T will have a bias of ε 2 such that: where we have neglected small contributions that appear at higher order in ε, higher order in τ Th or lower order in T .Computing the Fisher matrix we have: It is informative to compare this result to the variance on the bispectrum estimator in Eq. (6.19).Since ℓ 2 C T T ℓ falls with ℓ (at sufficiently high ℓ), the two results agree in the limit where the dominant contributions to the sum in Eq. (6.19) come from squeezed configurations with ℓ ≫ 1 -the triangle rule then implies that either the term proportional to C T T ℓ ′ or C T T ℓ ′′ dominates the bispectrum.Said differently, because the reconstruction of T improves by measuring many small-scale modes, we mainly capture information about squeezed configurations of the bispectrum where T Sc is evaluated at low-ℓ while the other power of T Sc and T dSc are evaluated at high-ℓ.We expect to be in this regime, since screening occurs mainly on small scales (e.g. it is associated with halos).Note that a completely analogous situation arises in kinetic Sunyaev Zel'dovich velocity reconstruction, as described in detail in Ref. [74].
Including polarization, the best sensitivity on ε can be obtained from the T dSc E Sc B Sc bispectrum, where the estimator variance is where was defined in Eq. (5.21).This bispectrum can yield a competitive sensitivity compared to the temperature-only bispectrum above since The sensitivities to ε using the T dSc T Sc T Sc and T dSc E Sc B Sc bispectra are shown in Fig. 9.The sensitivity of the bispectrum estimator for the experimental configurations studied here is slightly weaker than the result from the CMB auto-correlation.This is due to the smallness of τ Th , which, parametrically, suppresses the sensitivity compared to CMB auto-correlation as well as LSS crosscorrelation by an ε-independent factor of (C τ τ Th ℓ ) 2 /C τ τ ℓ ∼ 10 −11 − 10 −12 , depending on m A ′ (see Fig. 5).However, this estimator also scales as ε 2 , and brings about the most favorable scaling ∼ ℓ −3/4 with the number of modes measured at significant S/N.
To explore the improved sensitivity on ε from small-scale modes we evaluate he bispectrum estimators (6.19) and (6.22) for Planck and CMB-HD.For Planck we capture most of the SNR by summing from ℓ = 2 up to ℓ = 5000, while for CMB-HD there is significant SNR up to ℓ = 9000.From Fig. 4, at large dark photon mass, there is a lot of structure at high ℓ in the dark screening optical depth coming from the 1-halo term.Due to the favorable scaling with ℓ in the bispectrum constraints, the CMB-HD sensitivity is greatly enhanced by the additional support on small scales of the cross-power spectrum, despite the decreasing monopole τ and the smallness of τ Th .Note that if we extended the range of halo masses in our halo model below the conservative 10 11 M ⊙ , the sensitivity on ε would increase.However, this would also bring about larger uncertainties in what the appropriate gas profile to be used on those small scales might be.

Conclusion and remarks
In this paper, we have studied the resonant conversion of CMB photons to a hypothetical dark photon inside LSS at low redshift.Conversion leads to a frequency dependent patchy 'dark' screening of the CMB.Observationally, dark screening manifests as an anisotropic and frequency dependent optical depth τ (ω, n), which can in principle be extracted from CMB data.The sensitivity to ε for a variety of two-and three-point correlation functions are shown in Fig. 9, which we summarize in the following (see also Table 1).The projected sensitivity of several estimators for Planck, as well as future surveys CMB-S4 and CMB-HD.The gray shaded region is excluded from analysis [31,34] with data from COBE/FIRAS [45].The blue shaded region shows the equivalent constraint using our model, as explained in Section 6.1.The solid lines show the projected sensitivity of the CMB auto-correlation functions in Section 4.3 and 6.3 with uncertainty given by Eq. (6.12).The dashed contours were computed using Eq.(6.16) and show the projected sensitivity from the cross-correlation between CMB and LSS in Section 5.2 and 6.4.The dotted and dot-dashed contours given by Eq. (6.19) and Eq.(6.22) respectively show the projected sensitivities of the bispectra presented in Section 5.4 and 6.5.The projected sensitivity from CMB-S4 and CMB-HD are similar for the CMB auto-correlation and cross-correlation between CMB and LSS, whereas bispectra sensitivity with CMB-HD is superior to CMB-S4.We used f sky = 0.7 for Planck and f sky = 0.5 for CMB-S4 and CMB-HD.
The patchy dark screening optical depth τ (ω, n) can be measured using CMB data alone.The global signal, i.e. the τ (ω) monopole, leads to a spectral distortion which is constrained by COBE/FIRAS [45] and can potentially be measured with future experiments targeting spectral dis-tortions.The constraint from conversion at low-redshift in non-linear structure that we obtain using existing FIRAS data [45] is consistent with previous limits obtained from conversion over a wider range of redshift and a different treatment of inhomogeneities.
Extending previous analyses, we have demonstrated that CMB and LSS correlation functions are in principle a far more powerful probe of photon to dark photon conversion.The two-point function of dark screening (⟨τ (ω, n)τ (ω, n′ )⟩) can be extracted from existing and future CMB data by measuring ⟨T dSc (ω, n)T dSc (ω, n′ )⟩.The patchy dark screening map T dSc (ω, n) can be separated from the blackbody CMB by taking advantage of multi-frequency observations of the CMB.As shown in Fig. 9, such a correlation function can be more sensitive than the existing FIRAS constraints by 2 − 3 orders of magnitude despite scaling with the small kinetic mixing parameter as ε 4 .
A better reach on the mixing parameter can be obtained by cross-correlating τ (ω, n) with other observables that are sensitive to the underlying distribution of electron density in the Universe.One such correlation function is ⟨T dSc (ω, n)τ (n ′ )⟩, where τ (n ′ ) is a template for the dark screening optical depth.This correlation function takes advantage of the fact that for a dark photon with mass ≲ 10 −12 eV, conversion mostly happens at late times, creating a strong correlation between tracers of LSS and patchy dark screening.As shown in Fig. 9, the reach obtained from this correlator can be about 1 − 2 orders of magnitude better than the CMB-only result.
Cross-correlation functions can also be constructed from the dark-screened and Thomson-screened CMB alone.Qualitatively, these correlation functions can be understood as a correlation between the dark screening τ (ω, n) and the Thomson screening optical depth τ Th from halos.Three-point correlation functions offer the best sensitivity to this correlation.In particular, the ⟨T dSc T Sc T Sc ⟩ and ⟨T dSc E Sc B Sc ⟩ bispectra (both shown in Fig. 9) offer the best sensitivity, and are comparable to the reach anticipated with CMB auto-correlation functions.Compared to the cross-correlation with a template, the reduction in sensitivity is mainly a result of the smallness of the Thomson screening optical depth τ Th , which suggests that these three-point correlation functions could be better probes of photon to dark photon conversion in the weakly inhomogeneous Universe around recombination, when τ Th was much larger.Compared to the CMB auto-correlation functions, the sensitivity of these three-point correlation function scales much more favorably with ℓ, and hence improves with increased sensitivity and high resolution -the regime targeted by future surveys.The similarity in the reach of the CMB auto-correlation functions and these three-point correlation functions for future CMB survey is a numerical coincidence, and the relative strength of these two methods in a real data analysis likely depend on systematics and foregrounds, an investigation that we postpone to future work.Finally, we note that in the event of a detection, a combination of the two-point and three-point functions can be used to break degeneracies between the dark photon mass, kinetic mixing parameter and electron density profile, which is essential for extracting detailed information about the dark screening optical depth and how it correlates with the distribution of ionized gas.The methodologies we developed in this paper can be used to search for conversions of photon to dark photon in various other environment in the early universe, the details of which we will work out in a few follow up studies.
The study presented in this paper is a novel example of using cross-correlations between an observable in the Standard Model (SM) of cosmology (ΛCDM) and signals of a model beyond the Standard Model (BSM) of particle physics.Experimentally, the measurement of these correlators is enabled by the rapid improvement of cosmological experiments.Theoretically, these ⟨SM × BSM⟩ correlators allow us to use the ultra-high precision cosmological data on the anisotropic Universe to study BSM signals at the same order of the small BSM parameter as the monopole signal ⟨BSM⟩.We expect similar ⟨SM × BSM⟩ correlators will allow us to better probe other BSM signals with the rapidly improving cosmological CMB and LSS datasets, and take advantage of the synergy between the upcoming CMB experiments like the Simons Observatory [1], as well as future experiments CMB-S4 [2] and CMB-HD [3], with upcoming LSS surveys like DESI [4], Euclid [5], and LSST [6].Constructing new observables of this kind can allow us to better search for new interactions between the Standard Model and dark sectors, including dark matter annihilation, decay, and mixing between the visible and dark sector particles.
A Modeling dark screening in a dark matter halo

A.1 Dark matter halo models
The halo model of large scale structure assumes that all matter in the Universe is stored in virialized halos whose physical properties are fully described by the mass contained within their boundary (defined e.g. by the virial radius).Galaxies occupy dark matter halos and in doing so they act as tracers for the underlying dark matter distribution.The halo model is a semi-analytical framework used for understanding the non-linear structure of the matter distribution.There are two principal quantities needed to make predictions: the halo mass-function and the halo density profile.The former describes the halo number density as a function of mass and redshift.The latter describes how mass is distributed within each halo.Unlike the mass-function, it is not universal, meaning that it depends on cosmology and astrophysics.To make matters more simple, these expressions are assumed to be a function of a few variables such as mass, redshift and halo radius, and the parameters that enter these expressions are generally obtained from a mix of analytic predictions, simulations and even data.Other quantities that need to be specified in a halo model are the halo bias function, which to first order is fully determined by the mass function, as well as a concentration-mass relation that gives a characteristic scale radius for the halo density profile.Useful reviews on halo models are e.g.[46,47].
To perform halo model computations throughout this paper we assume the mass-function of [53] that fixes the bias function [75], and the concentration-mass relation from [76], which fixes the free parameters in the halo density profile.We also work under the assumption that the halo boundary is its virial radius, so that the halo mass is defined in a sphere of radius r vir .Our computations use the code hmvec3 , which we used throughout this work.A detailed description of the assumptions that enter the code can be found in the Appendix B of [74].
For our modeling, we use 50 redshift bins of equal comoving radial width in the range z = [0.01,z max ], where the reionization redshift z max ∈ {2, 6, 10}.The first case was used with the purpose of obtaining a template angular power spectrum for the distribution of galaxies, as measured by a futuristic LSS survey.The latter two provide a conservative range for when reionization was completed.Additionally, the cosmology (to be precise, the linear matter power spectrum) is defined for 10 4 comoving wavenumber bins k logarithmically spaced in the 10 −4 − 10 3 Mpc −1 range.
Throughout the paper, unless m 200 is written explicitly (e.g. when we define the gas profile in Eq. (2.5) from [50]), the symbol m denotes the halo virial mass.In terms of where resonant dark photon conversion can take place, we considered 100 halo mass bins logarithmically spaced in the 10 11 − 10 17 M ⊙ interval.The lower bound is a conservative mass limit of halos with feedback processes significant enough to disrupt the gas profile.For the upper bound, the number density of halos in the halo model is exponentially suppressed with mass, and using the halo mass function in [53], one can expect less than one halo with a mass > 3 × 10 16 M ⊙ in a volume the size of the Hubble sphere.In other words, we consider all halos with mass ≥ 10 11 M ⊙ .

A.2 Charged particle density profiles
In our calculations we considered an idealized scenario where reionization takes place instantaneously: the ionization fraction goes from zero to unity at exactly z reio .We compared the case where z reio = 6 and z reio = 10 and found no significant difference e.g.see Fig. 3.However, if we relax the assumption about the lower mass bound of halos where conversion can happen, this no longer holds true.In fact, we found that there is an order 1 − 10% difference in the magnitude of the optical depth monopole τ for masses above m A ′ > 10 −12 eV when we consider halos with mass down to m = 10 9 M ⊙ .This would not affect the contours in Fig. 9 significantly, but represents an example of a source of error in the modeling.
A potentially more important assumption has to do with the precise choice of density profile for the charged electrons on small scales.We explore this in more detail in the present section.As an example, the density of dark matter is consistently greater towards the core of halos (in fact, it is unbounded at r → 0) than the gas density.Repeating all calculations under the assumption that electrons follow the NFW profile [51] in Eq. 2.7 instead of expression (2.5) for gas from [50] that we have been using in the main text, we obtain the results shown in Fig. 10 for the sky-averaged dark screening optical depth monopole, and in Fig. 11 for the differential (dimensionful) monopole.Notice that the greatest disagreement between the curves lies in the upper half of the dark photon masses we considered.This is related to the differences in the two density profiles on small scales which we elaborate on further.
Since the NFW profile is unbounded at r → 0, we imposed that no resonant conversion happens below the scale radius r s of any halo, effectively adding a factor of Θ(r s − r res ) in the expression for the radial probability in (2.11).The NFW profile is monotonically decreasing with radius, hence this approximation excludes halos where the resonance condition is met near their core.As is evident in Fig. 10, the constraint is most relevant for the production of heavy dark photons, which require the largest overdensities.The same quantitative difference between density profiles shows up in the shape of the angular two-point functions, which is depicted in Fig. 12.Here, the relevant quantity is the ratio between the 1-halo and 2-halo terms.The distribution of power on small angular scales is influenced by the shape of the density profiles.However, the overall magnitude is given by the corresponding monopole at every mass.
For consistency reasons, we also assume that the NFW model breaks down at the scale radius of the Milky Way, which leads to a hard upper boundary on the range of dark photon masses that can be considered in this case, given by m A ′ ∝ ρ NFW (r MW s ) ≈ 2.86 × 10 −12 eV.As mentioned in the main text, for the Milky Way we assuming a virial radius and virial mass from [54], and the concentration-mass relation at z = 0 from [55] which can be used to compute the scale radius via r MW vir = c MW (m MW )r MW s .Overall, it appears that the effect of changing the density profile modeling on the overall sensitivity on ε is minimal, as shown in Fig. 13.This suggests that our projection is relatively insensitive to the exact electron density model for masses below the mass limit set by the Milky Way.This is a indication that the assumptions going into our forecast analysis are reasonable and robust, at least on scales above the scale radii of halos.However, the precise modeling of the gas profile around the core regions of halos will ultimately dictate the sensitivity on ε at larger dark photon masses.On the other hand, for lower dark photon masses, an improved sensitivity can be reached with improving knowledge of the electron density profile in the region right outside the virial radius of a halo [63].Improving the reach by more than half an order of magnitude in mass towards lower masses is likely with new data from ACT DR6 and DESI.
To conclude, in recent years there has been significant progress in our understanding of the electron density distributions and fluctuations in halos, from new experimental measurements [63] and numerical simulation [77].Prospective studies with cosmological data and CHIME/FRB [65] will further shrink these uncertainties, and potentially allow us to have more robust projections on the sensitivity at higher dark photon masses.

B Correlation functions of dark screening
In this appendix, we derive in detail the angular power spectrum of the two-point auto-correlation of optical depth fluctuations.We use the notation conventions from the main text.First, we recall that Comparison between the average dimensionful optical depth for the NFW and gas profiles from [51] and [50], respectively.Here two cases are plotted: the top edge in each line corresponds to taking z reio = 10, while the bottom edge is computed for z reio = 6.In either case, assumptions about the end of reionization are of little importance.As described in the main text, the shape and magnitude of the monopole η plays an important role in determining the shape and reach in sensitivity of ε, for any experiment and forecast method.Therefore, the difference in magnitude between the two curves depicted here illustrates the importance of modeling the charged electron density in halos.Notice that the range of relevant redshift bins is sensitive to this choice, as well as the magnitude of the signal.For example, if we have an LSS template up to z = 2 as assumed in the main text, then we expect the cross-correlation with CMB dark screening to be less strong at m A ′ ≳ 10 −12 eV if gas traces dark matter in the absence of AGN feedback processes.The sharp falloff in the NFW case at the high z end is due to the assumption that no conversion happens below the scale radius in any halo.
the optical depth can be written as a sum of a homogeneous and an anisotropic part:  The left panel shows the dark screening power spectra for three choices of dark photon mass, given a conversion model that assumes electrons in halos follow ρ gas , as in the main text.On the right we show the equivalent spectra where ρ NFW is assumed.The ρ gas profile is flat near the core and drops slower with radius than ρ NFW , hence at low m A ′ we expect more power on larger scales on the left side.Meanwhile, ρ NFW favours larger densities near the halo cores, so at large mass m A ′ we expect more power on smaller scales.All spectra are multiplied by their respective monopole ηgas and ηNFW depicted in Fig. 10.13: Comparison between sensitivity given by two different profiles for the distribution of electrons inside halos, assuming CMB-S4-like noise.The green contours are the same as those shown in Fig. 9 for the 'AGN feedback' gas profile [50], while the black contours show the equivalent constraints starting from the NFW density profile [51].As previously, the FIRAS constraint is computed taking into account the total average optical depth, including the Milky Way contribution, which for NFW is monotonically increasing with mass and dominates over the extra-galactic component at the upper mass end.The solid black boundary at m A ′ ≈ 3 × 10 −12 eV is imposed by assuming the NFW model breaks down where m 2 A ′ meets the resonance condition on scales less than the Milky Way scale radius r MW s .For the contours that take into account the CMB anisotropies, we point out that the match between the contours shows that our model assumptions are robust on scales between the virial radius and the scale radius.However, the precise modeling of the gas profile around the core regions of halos will ultimately dictate the sensitivity on ε at larger dark photon masses.
To fix the notation, we define the volume elements d 3 k = k 2 dk d 2 nk and d 3 χ = χ 2 dχ d 2 n, where the solid angle d 2 n = sin θdθdϕ.The two-point function of the overall optical depth is defined in configuration space as follows: where we introduced the notation ξ δτ δτ = ⟨δτ * (χ 1 )δτ (χ 2 )⟩ to represent the two-point auto-correlation of anisotropies.The physical interpretation is the following: given that a photon following trajectory n1 on the observer's sky undergoes resonant conversion in a halo with mass m 1 at redshift z 1 , what is the chance that another photon traveling along n2 converts as well?In the halo model, this probability is a sum of two terms: The first describes the case where the two photon conversions happen within the same halo centered at χ(z 1 ), and the second term represents the contribution from two different halos at χ(z 1 ) and χ(z 2 ).
To derive the expressions for the sky-averaged power spectra, we will need to perform a sum over all halos in the halo model, project onto spherical harmonics, then integrate over the Hubble volume.
In the derivations below we sometimes leave the mass and redshift dependence implicit to simplify notation.We also use redshift z and comoving distance χ(z) interchangeably to denote the parametric dependence on time of various quantities.

B.1 One-halo term
We start by computing the 1-halo term.Labeling all halos within each redshift/comoving distance bin by i, we get: In the last line, we used the expression for the average halo number density from Eq. (3.2).Only the second integral is now angle-dependent, so we project it onto spherical harmonics: To get to the second equality we assumed that the halo profile u(n) ∈ R has azimuthal symmetry implying that m u ℓm Y ℓm = u ℓ0 Y ℓ0 .Next we used the definition of the Wigner D-matrices to rotate the spherical harmonics.Then, we enforced the orthonormality condition d In the last step we used the symmetries of the Wigner D-matrices and their relationship to spherical harmonics and Legendre polynomials P ℓ to simplify the equation.
We perform a multipole expansion on the left hand side of Eq. (B.4): We introduced the angular power spectrum of the 1-halo term C 1−halo ℓ , here as a function of redshift.To simplify the right-hand side we use the property that 4π 2l+1 m Y * ℓm Y ℓm = P ℓ .Next we match the terms in ℓ from both expressions above: To obtain the full sky-averaged power spectrum we need to integrate over the comoving radial coordinate χ.Changing also the integration variable to redshift, i.e. dχ = dz/H(z) where H(z) is the Hubble constant, we obtain the final expression for the 1-halo term:

B.2 Two-halo term
For the second term we follow similar reasoning and steps as in the previous case, except now we need to consider the conversion probability inside two different halos.Labeling these by i and j, located at arbitrary comoving distances χ 1 and χ 2 we find the average configuration space two-point function as: The term in the brackets is related to the number density of each halo's characteristic m and z, as well as to the correlation between positions of halos ξ hh in the following way: Known as the halo-halo auto-correlation function, ξ hh is proportional to the linear matter two-point function.To first order in linear theory, the following is true: where the bias function b(z, m) is a deterministic function of the halo mass and redshift.We can further relate this to the linear matter power spectrum by doing a Fourier expansion over comoving wavenumbers: On large scales, the linear matter power spectrum is well approximated by Consider next the multipole expansion of a plane wave: where j ℓ (kχ) ∈ R is the spherical Bessel function.Performing the multipole expansion on ξ hh from (B.11) and putting everything together, we obtain the halo-halo angular power spectrum to first order in linear theory: We come back to the expression for the 2-halo two-point function in Eq.B.9 and project the right hand side onto spherical harmonics.For this, we generalize a result we obtained in the previous section for the angle-dependent integrand: Multiplying by the halo-halo auto-correlation and using the properties of Y ℓm 's and D ℓ mm ′ 's to simplify the expression we get (B.17) The multipole expansion for the general 2-halo two-point function is

C Two-point correlators and quadratic estimators
In this appendix we enumerate the various two-point correlation functions among the temperature and polarization anisotropies.We then list all quadratic estimators for the un-screened CMB temperature and polarization anisotropies as well as the Thomson optical depth and photon to dark photon optical depth.
The temperature correlators are: The E-mode correlators are: .

(C.50)
There are five estimators for the photon to dark photon optical depth:

Figure 2 :
Figure2: The average dimensionful optical depth defined in Eq. (3.5) as a function of redshift, for a range of dark photon masses.To obtain the full optical depth, one needs to integrate over redshift and multiply by the unknown mixing parameter divided by the frequency of the photon, ε 2 /ω.Notice that the lighter masses produce the strongest signal at low redshift.Meanwhile, heavier dark photons probe redshifts all the way to reionization but require crossing large densities in order for the resonant conversion to take place.These are found within more massive halos whose number density is suppressed.

Figure 3 :
Figure 3: The dimensionful optical depth monopole as a function of dark photon mass m A ′ .The extra-galactic contribution η (blue) is defined by Eq. (3.5).Under our assumptions, resonant conversion occurs once reionization is completed and here we plot the case where z reio = 6 and where z reio = 10.The maximum change in the amplitude of η is of order 10 2 eV, less than the thickness of the blue line, and we conclude that the model is insensitive to the details of reionization over the range of masses we are probing.We therefore assume z reio = 6 in the remainder of this work.At low mass, η → 0 due to the constraint that conversion only happens within the virial radius in each halo.The upper bound is set by the shape of the gas profile and details about the halo model.The contribution due to gas in the Milky Way η MW (black) is defined in Eq. (3.6).

Figure 4 :
Figure 4: The angular power spectrum of dark screening fluctuations, C δτ δτ ℓ

Figure 5 :
Figure5: We show the angular power spectrum of the two-point function of Thomson screening (blue), dark screening (orange), as well as their cross-correlation (green).We chose the value of ε = 10 −6 and fixed m A ′ = 6 × 10 −13 eV and ω = 30 GHz.We also present the ε-independent combination C τ τ Th ℓ

Figure 6 :
Figure6: Exclusion contours on coupling constant ε function of dark photon mass m A ′ from COBE/FIRAS[45] data.We see that the uncertainty in the measurement of T gives a constraint on the magnitude of the dark-screened temperature monopole T dSc .Here we used the total dark screening optical depth due to both galactic and extra-galactic contributions τ +τ MW .This constraint provides an upper bound on ε of at most 10 −6 .

Figure 8 :
Figure 8: Scale comparison.We fixed m A ′ = 10 −12 eV and ε as shown in each panel.We show the isotropic component to the XX = T T, EE power spectra due to dark screening C X dSc X dSc ℓ ∝ ε 4 (orange), Thomson screening C X Sc X Sc ℓ

Figure 9 :
Figure9: The projected sensitivity of several estimators for Planck, as well as future surveys CMB-S4 and CMB-HD.The gray shaded region is excluded from analysis[31,34] with data from COBE/FIRAS[45].The blue shaded region shows the equivalent constraint using our model, as explained in Section 6.1.The solid lines show the projected sensitivity of the CMB auto-correlation functions in Section 4.3 and 6.3 with uncertainty given by Eq. (6.12).The dashed contours were computed using Eq.(6.16) and show the projected sensitivity from the cross-correlation between CMB and LSS in Section 5.2 and 6.4.The dotted and dot-dashed contours given by Eq. (6.19) and Eq.(6.22) respectively show the projected sensitivities of the bispectra presented in Section 5.4 and 6.5.The projected sensitivity from CMB-S4 and CMB-HD are similar for the CMB auto-correlation and cross-correlation between CMB and LSS, whereas bispectra sensitivity with CMB-HD is superior to CMB-S4.We used f sky = 0.7 for Planck and f sky = 0.5 for CMB-S4 and CMB-HD.
Figure10: Comparison between the average dimensionful optical depth for the NFW and gas profiles from[51] and[50], respectively.Here two cases are plotted: the top edge in each line corresponds to taking z reio = 10, while the bottom edge is computed for z reio = 6.In either case, assumptions about the end of reionization are of little importance.As described in the main text, the shape and magnitude of the monopole η plays an important role in determining the shape and reach in sensitivity of ε, for any experiment and forecast method.Therefore, the difference in magnitude between the two curves depicted here illustrates the importance of modeling the charged electron density in halos.

Figure 11 :
Figure 11:  Comparison between the differential average optical depth per redshift bin for two different choices of density profile.Notice that the range of relevant redshift bins is sensitive to this choice, as well as the magnitude of the signal.For example, if we have an LSS template up to z = 2 as assumed in the main text, then we expect the cross-correlation with CMB dark screening to be less strong at m A ′ ≳ 10 −12 eV if gas traces dark matter in the absence of AGN feedback processes.The sharp falloff in the NFW case at the high z end is due to the assumption that no conversion happens below the scale radius in any halo.

Figure 12 :
Figure12: The left panel shows the dark screening power spectra for three choices of dark photon mass, given a conversion model that assumes electrons in halos follow ρ gas , as in the main text.On the right we show the equivalent spectra where ρ NFW is assumed.The ρ gas profile is flat near the core and drops slower with radius than ρ NFW , hence at low m A ′ we expect more power on larger scales on the left side.Meanwhile, ρ NFW favours larger densities near the halo cores, so at large mass m A ′ we expect more power on smaller scales.All spectra are multiplied by their respective monopole ηgas and ηNFW depicted in Fig.10.

Figure
Figure 13: Comparison between sensitivity given by two different profiles for the distribution of electrons inside halos, assuming CMB-S4-like noise.The green contours are the same as those shown in Fig.9for the 'AGN feedback' gas profile[50], while the black contours show the equivalent constraints starting from the NFW density profile[51].As previously, the FIRAS constraint is computed taking into account the total average optical depth, including the Milky Way contribution, which for NFW is monotonically increasing with mass and dominates over the extra-galactic component at the upper mass end.The solid black boundary at m A ′ ≈ 3 × 10 −12 eV is imposed by assuming the NFW model breaks down where m 2 A ′ meets the resonance condition on scales less than the Milky Way scale radius r MW