Dark Photon Oscillations in Our Inhomogeneous Universe

A dark photon may kinetically mix with the ordinary photon, inducing oscillations with observable imprints on cosmology. Oscillations are resonantly enhanced if the dark photon mass equals the ordinary photon plasma mass, which tracks the free electron number density. Previous studies have assumed a homogeneous Universe; in this Letter, we introduce for the first time an analytic formalism for treating resonant oscillations in the presence of inhomogeneities of the photon plasma mass. We apply our formalism to determine constraints from Cosmic Microwave Background photons oscillating into dark photons, and from heating of the primordial plasma due to dark photon dark matter converting into low-energy photons. Including the effect of inhomogeneities demonstrates that prior homogeneous constraints are not conservative, and simultaneously extends current experimental limits into a vast new parameter space.

A dark photon may kinetically mix with the ordinary photon, inducing oscillations with observable imprints on cosmology. Oscillations are resonantly enhanced if the dark photon mass equals the ordinary photon plasma mass, which tracks the free electron number density. Previous studies have assumed a homogeneous Universe; in this Letter, we introduce for the first time an analytic formalism for treating resonant oscillations in the presence of inhomogeneities of the photon plasma mass. We apply our formalism to determine constraints from Cosmic Microwave Background photons oscillating into dark photons, and from heating of the primordial plasma due to dark photon dark matter converting into low-energy photons. Including the effect of inhomogeneities demonstrates that prior homogeneous constraints are not conservative, and simultaneously extends current experimental limits into a vast new parameter space.
Introduction.-A dark photon, A , may kinetically mix [1] with the ordinary photon, γ. Kinetic mixing is one of a few portals allowing new physics to couple to the Standard Model (SM) through a dimensionless interaction that can be manifest at low energies. A further motivation for dark photons is that they may constitute dark matter [2][3][4][5][6][7][8][9][10][11][12]. Very light dark photons decouple from experiments as a positive power of m A /E exp , where m A is the mass of the dark photon and E exp is the experimental energy scale. Because of this decoupling behavior, light dark photons with sizable interactions are consistent with current experimental constraints.
Kinetic mixing induces oscillations of photons into dark photons, γ → A , as well as the reverse process, A → γ. In the early Universe the photon has a plasma mass, m γ (z, x), that tracks the free electron number density, n e (z, x). The oscillation probability is resonantly enhanced if the plasma mass equals the mass of the dark photon, m γ (z, x) ≈ m A . For dark photon masses in the interval 10 −14 m A 10 −9 eV, the homogeneous (spatially averaged) value of the plasma mass, m γ (z), crosses the dark photon mass after recombination. In this regime there are powerful constraints [13,14] from γ → A oscillations distorting the shape of the Cosmic Microwave Background (CMB) energy spectrum measured by FI-RAS [15]. Dark matter composed of dark photons is also constrained by heating of the primordial plasma from resonant A → γ oscillations [16], producing low-energy photons that are efficiently absorbed by baryons. Additional constraints on dark photon dark matter are considered by Refs. [4,[17][18][19].
Previous studies of cosmological dark photon oscillations have assumed a homogeneous plasma mass, i.e., m γ (z, x) ≈ m γ (z). However, the plasma mass has perturbations that track inhomogeneities in the electron number density, which are a predicted consequence of the growth of structure in the early Universe. Consider a photon that propagates along a worldline through the primordial plasma. In the homogeneous limit, this photon may experience a level crossing at a specific redshift, z res , when m γ (z res ) ≈ m A . In reality, a photon's path traverses regions with overdensities and underdensities, and may pass through many different level crossings at redshifts that differ from z res . Fig. 1 shows a simulation of this process for a dark photon mass with z res ≈ 100; perturbations in the plasma mass induce resonant conversions over a wide range of redshifts, 90 z 110. The effect of inhomogeneities is especially dramatic for dark photons with masses m A 10 −14 eV, which experience no level crossings in the homogeneous limit, but in reality can experience level crossings in regions with lower-than-average electron number density.
In this work, we initiate the study of resonant oscillations between photons and dark photons in the presence of inhomogeneities in the photon plasma mass. We introduce an analytic formalism for calculating the probability that photons or dark photons oscillate as they travel through the inhomogeneous plasma. As applications of our formalism, we revisit bounds from the CMB energy spectrum on photons oscillating to dark photons, γ → A , and bounds from energy injection due to dark photon dark matter oscillating into ordinary photons, A → γ. We find that these bounds require significant revision: compared to the homogeneous limit, perturbations both induce new resonances in underdensities and overdensities, extending these bounds into a vast new parameter space, and can also wash out resonances, making the sensitivity derived in the homogeneous approximation an overestimate for certain dark photon masses. The homogeneous limit is therefore not a conservative approximation of our inhomogeneous Universe.
Dark photons with masses 10 −15 m A 10 −9 eV  [20,21], while Dark SRF [22,23] aims to produce and detect dark photons without assuming a cosmic abundance [24]. Our cosmological constraints inform the model space targeted by these experiments.
The remainder of this Letter is organized as follows. We begin by reviewing γ ↔ A oscillations. We then introduce our analytic formalism for treating these oscillations in the presence of perturbations of the photon plasma mass, leaving the details of our formalism to a companion paper [25]. Next, we apply our formalism to determine the constraints on γ → A oscillations from FIRAS data. We then show how inhomogeneities extend constraints on energy injection from dark photon dark matter to new dark photon masses. Our conclusions highlight additional possible applications and extensions of our formalism. Finally, our appendices elaborate on our analysis of FIRAS data and explore systematic uncertainties associated with our constraints. Throughout this work, we use units with = c = k B = 1, and the Planck 2018 cosmology [26]. For reproducibility, we provide links in the figure captions (ƭ) pointing to the code used to generate them. Resonant photon-dark photon oscillations.-We consider the following photon-dark photon Lagrangian, where is a dimensionless measure of kinetic mixing. The propagation of CMB photons in the primordial plasma leads to birefringent effects that are described by a mass term, m γ , in the photon dispersion relation.
There are positive and negative contributions to m 2 γ from scattering off free electrons and neutral atoms [13,14]: where ω(z) is the photon energy, and n e (z, x) and n HI (z, x) represent the local free electron and neutral hydrogen densities, respectively. We model the evolution of cosmological quantities using CLASS [27] interfaced with HyRec [28]. For 1, γ → A conversion is a resonant process that is efficient only when the dark photon mass is equal to the plasma mass; in this limit, we can apply the Landau-Zener approximation for non-adiabatic transitions [13], where i indexes times t i when m 2 γ (t i ) = m 2 A and the resonance condition is met. Eq. (3) assumes P γ→A 1, which applies throughout this work.
Eq. (3) describes the probability that a photon will convert along its path, which depends on m γ (t) along this path. In the homogeneous limit, this average is given simply by substituting the average photon plasma mass m γ (z) into Eq. (3). The photon plasma mass for a range of frequencies spanning the FIRAS measurements is shown in the top panel of Fig. 2. Although several level crossings are possible in general, the dependence of the transition probability on ω and m 2 γ means that the latest resonance typically dominates the total conversion probability. The effect of inhomogeneities.-Inhomogeneities in the photon plasma mass substantially affect the conversion probability of photons into dark photons and vice versa, allowing for efficient oscillations over a range of cosmic times rather than at a single epoch.
In the presence of plasma mass inhomogeneities, we need to take the average of Eq. (3) over different photon paths to account for transitions in locally overdense and underdense regions. This problem reduces to integrating over m 2 γ at each point in time, weighted by the probability density function of finding a region with plasma mass m 2 γ . Our formalism draws from Rice's formula for the average number of level crossings of a random field [29]. As we show in Ref. [25], the differential conversion probability reads where f (m 2 γ ; t) is a probability density function (PDF) of m 2 γ at time t, and δ D is the Dirac delta function. Neglecting perturbations in the free electron fraction x e (see Ref. [25] for a discussion on why this assumption is valid within the range of redshifts considered here), Eq. (2) shows that m 2 γ (z, x) ∝ n b (z, x), the baryon number density; this implies that where P(δ b ; t) is the one-point PDF of baryon density fluctuations δ b ≡ (n b − n b )/n b . The proportionality m 2 γ ∝ n b together with the definition of δ b implies that that 1 + δ b = m 2 γ /m 2 γ . Eq. (5) therefore ties the physics of γ ↔ A directly to a cosmological observable.
The concise but powerful formula shown in Eq. (4) is one of our main results, and can be applied to a variety of scenarios related to resonant conversions in inhomogeneous media. In particular, it allows us to study the problem of oscillations even at low redshifts, where nonlinear effects are important and the PDFs of the matter density fields are far from Gaussian.
We consider a few different possibilities for the onepoint PDF P(δ b ; t) in order to estimate the size of the theoretical uncertainty associated with the nonlinear distribution of matter at low redshifts z 6. First, we consider a log-normal distribution, which has long been used as a simple model for the distribution of the lowredshift matter density [34][35][36][37]. To inform the spectrum of fluctuations for this distribution, we use the baryonic power spectra P b (k) derived from hydrodynamic simulations [38][39][40][41] and extracted in Refs. [42,43]. Second, we adopt an analytical prescription [44], which extends the spherical collapse model [45,46] to perform a firstprinciples computation of the nonlinear matter PDF. We note that the former prescription models the power spectrum of electron number density fluctuations while using a simpler log-normal ansatz for the underlying PDF, while the latter provides a more sophisticated description of the matter PDF while not directly accounting for the bias between the distribution of matter and that of free electrons.
In the literature, P(δ b ; t) is typically defined as a function of a smoothing scale R over which densities are averaged in order to match observations and simulation results; furthermore, the width of the distribution can exhibit a log-divergence in R if P b (k) ∝ k −3 at large k. In our work, we assume that baryonic structures are suppressed on scales larger than the baryon Jeans scale R J ∼ 10 kpc; in practice, the log-normal PDF is computed with P b (k) which has this cut off at k ∼ 1/R J , while our analytic PDF is obtained with a smoothing scale R = R J . A complete description of our PDFs appears in our companion paper [25].
The differential transition probability (normalized to unity) for a few benchmark dark photon mass points m A = 2 × 10 −15 , 10 −13 , and 10 −12 eV is shown in the bottom panel of Fig The 95% confidence level constraints on the kinetic mixing parameter as a function of dark photon mass m A , assuming log-normal (red) or analytic (blue) PDFs. We also show the reach of the proposed PIXIE satellite [30] (dotdashed red) assuming a log-normal PDF. For comparison we show the previous limit assuming a homogeneous plasma (dotted gray), a constraint from the magnetic field of Jupiter [31,32] (shaded brown), and the projected reach of the Dark SRF experiment [22,23] (dot-dashed orange). ƭ (Right) Constraints on dark photon dark matter from anomalous heating of the IGM during the epoch of HeII reionization, for the same PDFs. Prior constraints (shaded brown) come from non-resonant heating of the IGM [16] and heating of the gas in the dwarf galaxy Leo T [19]. We also show the projected reach of DM Radio Stage 3 [20,21,33] (dot-dashed orange). Limits from changes to the dark matter density and from IGM heating during the dark ages assuming a homogeneous plasma mass have been derived in Ref. [16] (dotted orange); we expect that these bounds will receive large corrections after including perturbations. ƭ additional broad resonance at z ∼ 6 is also present, corresponding to conversions in overdensities in the plasma mass post-reionization. For m A = 2 × 10 −15 eV, no resonance exists in the homogeneous limit; remarkably, however, fluctuations in the plasma mass result in resonant transitions over a broad range of redshifts at z 20 due to underdensities in the plasma mass. This opens up the possibility of probing dark photon masses m A 10 −14 eV through previously-neglected cosmological conversions.
Dark photon oscillations in the CMB spectrum.-We first apply our formalism to analyze the intensity of the CMB as measured by the FIRAS instrument aboard COBE [15] for evidence of deviations from a blackbody spectrum due to γ → A oscillations. Notably in this case the dark photon does not need to be the dark matter. The spectrum of the FIRAS data is fit by the nearly perfectly Planckian spectrum B ω with temperature T CMB = 2.725 K. For a given dark photon model specified by its mass m A and mixing parameter , the spectral distortion to the CMB spectrum will be given by where P γ→A is the conversion probability for the given model corresponding to the present-day frequency ω 0 , obtained by integrating Eq. (4). Details of the data analysis are presented in App. A.
Erring on the conservative side, we do not consider fluctuations outside of the range 10 −2 < 1 + δ b < 10 2 , and as such our results do not rely on conversions in the tails of the PDF where uncertainties are large. Additionally, for all cases considered here conversions in the redshift range 6 < z < 20 have been excised, providing a conservative result while being agnostic to the uncertainties arising from the complex physics of reionization in this epoch. We explore the effects of these choices in App. B.
We observe no significant evidence for a signal. In the left panel of Fig. 3 we show our fiducial constraints at the 95% confidence level on the dark photon mixing parameter for a range of dark photon masses m A . We show constraints using both the log-normal and analytic description of the PDF. We also show the projected limits for a future measurement of the CMB energy spectrum such as the proposed PIXIE satellite [30] using the putative specifications from Ref. [14]. The traditional constraint assuming a homogeneous plasma mass as a function of redshift is also shown for comparison, together with constraints projected by the Dark SRF experiment [22,23] and existing constraints from an analysis of the magnetic field of Jupiter [31]. There are bounds from black hole superradiance for values of m A that overlap our bounds assuming = 0 [47][48][49], but it is unknown if these apply when > 0 implying interactions of A with plasma around the black hole.
Dark photon dark matter.-Additional constraints apply when dark photons constitute the dark matter component of the Universe [2,4,8,10,12]. In particular, Ref. [16] proposed using measurements of the temperature of the intergalactic medium (IGM) around the epoch of HeII reionization (2 z 6) [50][51][52][53][54][55] to constrain the dark photon dark matter scenario. These measurements show that during HeII reionization, the total heat input per baryon is on the order of 1 eV. A → γ conversion for light A s produce soft photons that are absorbed efficiently through free-free absorption, leading to an anomalous heating of the IGM. The derived bound in the homogeneous limit extends over a limited mass range, precisely where the dark photon mass matches the homogeneous plasma mass in the narrow redshift range 2 z 6. Our formalism accounting for inhomogeneities extends this treatment to a wider range of dark photon masses.
The total energy injected per unit baryon E A →γ from the dark matter can be computed as whereρ A /n b is the ratio of the homogeneous dark matter energy density to baryon number density, which is a time-independent quantity. The total energy injected is then obtained by performing an integral over 2 < z < 6.
Considering the same PDFs discussed in the previous section and imposing E A →γ < 1 eV, we derive the constraints shown in the right panel of Fig. 3, with the homogeneous limit shown for comparison. Also shown is the parameter space covered by existing constraints [16,19], as well as the projected constraints from DM Radio Stage 3 [20,21,33]. Finally, our limits can be rescaled as a function of the maximum E A →γ allowed and ρ A by noting that E A →γ ∝ 2 ρ A .
Conclusions.-We have introduced a framework for treating oscillations between dark photons and ordinary photons as they traverse the inhomogeneous plasma of our Universe. Our main results are Eqs. (4) and (6). Further details of their derivation will appear in our companion paper [25]. We have applied this framework to determine constraints from CMB photons oscillating into dark photons (Fig. 3, left panel) and from energy injection from dark photon dark matter (Fig. 3, right panel).
Prior studies have assumed the homogeneous limit and require significant revision because inhomogeneities both extend the mass reach, and either strengthen or weaken the sensitivity for masses constrained in the homogeneous limit.
An accurate modeling of perturbations in the plasma mass requires modeling the non-Gaussian PDF of baryons down to very small scales, k ∼ 1 kpc, and in the regime of very large overdensities, 1 + δ 1, or underdensities, 1 + δ 1. We anticipate that more accurate constraints on dark photon oscillations can be set by using high resolution hydrodynamic simulations to constrain the PDF, or by extending analytical calculations to relevant regimes. We leave further refinement of these important systematics for future work.
We anticipate broader applications of our framework. Perturbations in the photon plasma mass will modify resonant oscillations of photons into axion-like-particles, which can occur in the presence of primordial magnetic fields [56] or dark magnetic fields [57]. Here we have considered oscillations of non-relativistic dark photon dark matter, but relativistic dark photons (or axion-likeparticles) can also resonantly inject photons, which can be tested by 21 cm observations [57][58][59]. We have here considered global (sky-averaged) effects, but photonto-dark photon oscillations in an inhomogeneous background will imprint anisotropies in the CMB that may be testable by Planck [26] and/or next-generation probes of CMB anisotropies [60,61].
Additional details of the data analysis performed and a discussion of systematic effects is presented in the Appendix. The code used to obtain the results in this paper as well as digitized constraints are available at https:// github.com/smsharma/dark-photons-perturbations. We are especially grateful to Misha Ivanov for many enlightening discussions regarding the analytic PDF of density fluctuations utilized in this work. AC acknowledges support from the "Generalitat Valencian" (Spain) through the "plan GenT" program (CIDEGENT/2018/019), as well as national grants FPA2014-57816-P, FPA2017-85985-P, and the European projects H2020-MSCA-ITN-2015//674896-ELUSIVES. HL is supported by the DOE under contract DESC0007968. SM and JTR are supported by the NSF CAREER grant PHY-1554858 and NSF grant PHY-1915409. SM is additionally supported by NSF grant PHY-1620727 and the Simons Foundation. JTR acknowledges hospitality from the Aspen Center for Physics, which is supported by the NSF grant PHY-1607611. This work made use of the NYU IT High Performance Computing resources, services, and staff expertise. The authors are pleased to acknowledge that the work reported on in this paper was substantially performed using the Princeton Research Computing resources at Princeton University which is consortium of groups including the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology's Research Computing department. This research has made use of NASA's Astrophysics Data System. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. This research made use of the astropy [62,63], CAMB [64,65], CLASS [27], HyRec [28], IPython [66] Jupyter [67], matplotlib [68], nbodykit [69], NumPy [70], seaborn [71], pandas [72], SciPy [73], and tqdm [74] software packages.

Appendix A: COBE/FIRAS data analysis
We utilize the low frequency FIRAS monopole data, 1 consisting of 43 linearly spaced data points spanning the frequency range ν 0 = 2.27-21.33 cm −1 and construct data covariance matrices following Ref. [15]. The spectrum of the FIRAS data is illustrated in Fig. 4, fit by the nearly perfectly Planckian spectrum with temperature T CMB 2.725 K, and ω 0 = 2πν 0 the angular frequency.
Residuals between the FI-RAS data and this spectrum are shown in the bottom panel. For a given dark photon model specified by its mass m A and mixing parameter , the spectral distortion to the CMB spectrum will be given by I ω0 (m A , ; T CMB ) = B ω0 (1 − P γ→A ), where P γ→A is the conversion probability for the given model corresponding to the present-day frequency ω 0 , obtained by integrating Eq. (4). An illustration of the distortion to the blackbody spectrum induced by dark photons of mass m A = 6 × 10 −15 eV with mixing = 4 × 10 −6 is shown in dashed blue in the bottom panel of Fig. 4.
We construct a Gaussian log-likelihood as where ∆ I = I (m A , ; T CMB ) − I d is the residual between the distorted CMB spectrum I (m A , ; T CMB ) = {I ω1 , I ω2 , . . .} and the FIRAS data vector I d , and C I d is the data covariance matrix. We treat the CMB temperature as a nuisance parameter and profile over it by maximizing the log-likelihood for T CMB at each {m A , } point. We define our test-statistic as whereˆ is the value of that maximizes the log-likelihood for a given m A , and obtain our limit by finding the value of at which TS = −2.71 corresponding to 95% containment for the one-sided χ 2 distribution.
1 Available at https://lambda.gsfc.nasa.gov/product/cobe/ firas_monopole_get.cfm. Although we consider the incoming photons to have energies corresponding to the specified FIRAS frequency bands, in reality FIRAS has a finite spectral resolution resulting in a spread in energies over a finite range of the order of the bin size. We check the impact of this binning on the constraints presented by modeling the response of each frequency bin as a Gaussian centered on the central frequency with standard deviation corresponding to the bin size [75]. We find that this choice has a percent-level impact on the computed inhomogeneous oscillation probabilities in the lowest frequency bin, with smaller errors for larger frequencies. Since the mixing angle constraint scales as the square root of the oscillation probability, our constraints are not qualitatively impacted by finite binning effects.
Nevertheless, resonances relying on particular values of ω 0 can cause local enhancements in the homogeneous constraints at masses m A 4 × 10 −13 eV due to individual crossings with small characteristic derivatives d ln m 2 γ (t)/dt in Eq. (3). These sharp enhancements are likely artifacts of finite spectral binning, and we thus smooth the homogeneous constraints with a Savitzky-Golay filter above this mass.

Appendix B: Systematics and Additional Results
The behavior of γ ↔ A conversions in the inhomogeneous Universe depends critically on the distribution of density perturbations as a function of redshift. While significant uncertainty exists for this distribution, we have already shown in the Letter that using two radically different approaches to computing the probability density function (PDF) of the photon plasma mass f (m 2 γ ; t) does not lead to qualitative differences in our results. In this section, we discuss several other possible sources of uncertainty, more consistency checks of our fiducial limits, and additional results that are more optimistic or are of pedagogical interest.

Probability Density Functions
The two different prescriptions for the one-point PDF used to construct f (m 2 γ ; t) are the log-normal distribution and an analytic distribution based on ideas presented in Refs. [44,45,[76][77][78][79]. The log-normal distribution is a phenomenological PDF that can take the nonlinear baryon power spectra from simulations into account, while the analytic distribution has a theoretical basis in structure formation theory, but only models the matter distribution without baryonic effects. Fig. 5 demonstrates that these two PDFs have a similar behavior in the range 10 −2 < 1 + δ < 10 2 despite being very different approaches, giving us confidence that considering fluctuations only in this range is a reasonable choice, and explaining why the limits derived from our two fiducial prescriptions are similar.
We have also considered three further models for the PDFs, which we believe are useful checks for our results, but are unlikely to be more accurate than our two fiducial approaches. In this section, we provide a brief description of these PDFs; for more details, we refer the reader to Ref. [25]. Log-normal with bias.-For our fiducial log-normal PDF, we use the nonlinear baryonic power spectrum as an input to determine f (m 2 γ ; t). Another approach in the literature is to add a bias parameter b, a constant ratio between baryonic fluctuations and matter fluctuations [81], to the log-normal distribution. This has been shown to be effective in modeling the one-point PDF extracted from galaxy count surveys [37,82].
As an independent check of the nonlinear baryonic power spectrum that we obtained from simulations, we use the log-normal with bias PDF together with the Halofit nonlinear matter power spectrum provided by CLASS [83] with a bias parameter b within the range of fit values obtained in Ref. [82]. This serves as an independent way of modeling the PDF without relying on the simulation results that we use for the fiducial lognormal distribution. We stress however that this model is unphysical, since negative fluctuations in the baryon density are mathematically allowed, calling into question the accuracy of the PDF for underdensities. Voids.-The simulation and characterization of voids in our Universe has been the subject of ongoing interest [84][85][86][87][88][89][90], and can inform the PDF f (m 2 γ ; t) for values of m 2 γ significantly below m 2 γ . To construct f (m 2 γ ; t) from these studies, we rely on the simulation results in We show the fiducial log-normal (red) and analytic (green) PDFs, together with several other PDFs considered in the appendix, including a log-normal PDF with bias b = 1.5 (blue), a PDF constructed from a model of voids (purple) [80] and a Gaussian PDF (orange). Also shown are the fiducial 10 −2 < 1+δ < 10 2 boundaries (dashed gray). At z = 0, the homogeneous plasma mass is m 2 γ = 1.9 × 10 −14 eV. ƭ Ref. [80], which provides a PDF for the mean density of voids. Together with the PDF for the volume of voids and the number of voids in their simulation, we find that voids typically occupy f void ∼ 10% of their simulation volume, and we rescale the void density PDF by this number to obtain a model of the PDF of finding a void of a certain mean density in our Universe and hence f (m 2 γ ; t). This approach gives an estimate of f (m 2 γ ; t) only for values of m 2 γ below the homogeneous value and should not be used for overdensities. Even so, it is highly conservative for two reasons: first, not all underdensities are local minima in position space, which is the working definition of a void, and second, we do not account for the density profile of the void, which neglects the fact that the centers of voids are likely to be significantly less dense than the mean density. Nevertheless, comparing this PDF with our fiducial choices can give us confidence in our modeling of underdensities. Gaussian.-Fluctuations in densities deep in the linear regime (z 200) are well-described by a Gaussian random field. In this limit, f (m 2 γ ; t) takes on a particularly simple form, making it useful for an intuitive understanding of our results. With a Gaussian PDF, the differential conversion probability in Eq. (4) is given by [25] where σ is the variance of plasma mass fluctuations. At late times, σ becomes of the same order as and exceeds Jupiter γ → A PDF Tails 10 −1 < 1 + δ < 10 1 10 −2 < 1 + δ < 10 2 10 −3 < 1 + δ < 10 3 10 −4 < 1 + δ < 10 4

Homogeneous
Log-normal PDF Analytic PDF . Progressively darker lines corresponding to the inclusion of larger underdensities and overdensities, from 10 times larger and smaller than the mean plasma mass to 10 4 times larger and smaller than the plasma mass, respectively. (Top right) Sensitivity of the constraints to the choice of the PDF parameterization. Shown is our fiducial constraint with the log-normal PDF (solid red) and the analytic PDF (solid green). Constraints with the inclusion of a linear bias b = 1.5 between the dark matter and baryons in the log-normal prescription (solid blue), the inferred PDF of underdensities in voids from Ref. [80] (solid purple), and with a Gaussian PDF (dotted orange) are also displayed. (Bottom) Effect of excluding conversions over different redshift ranges on the constraints are presented. Excising a larger redshift range 6 < z < 30 (dashed orange) has no effect on the fiducial limits, obtained by excising 6 < z < 20 (solid red). Constraints relying solely on conversions at redshifts above z > 0.1(1) are shown in solid green(purple), and those relying on conversions at linear cosmological epochs (z 20) are shown for the log-normal and analytic PDFs in solid blue and dashed red, respectively. ƭ m 2 γ , defining the nonlinear regime. Hence, it is not suitable for describing fluctuations at low redshifts on which our results critically depend; results using the Gaussian PDF should be taken as pedagogically interesting, but incorrect. Fig. 5 shows a plot of the different PDFs discussed here. Within the range of 10 −2 < 1 + δ < 10 2 , we can see that the two fiducial PDFs agree very well; outside this range, however, significant differences develop across all PDFs. Despite being highly conservative, the PDF constructed from voids generally agrees well with both the analytic and log-normal PDFs in the range 10 −2 < 1 + δ < 10 −1 , while the log-normal with bias PDF shows good agreement with the fiducial log-normal PDF for 1 + δ 1, even though they use different power spectra as inputs. We stress that we expect neither the void PDF nor the log-normal PDF with bias to agree with our fiducial results outside of the respective ranges specified here. Existing Constraints Ω A = Ω DM 10 −2 < 1 + δ < 10 2 DM A → γ PDF Tails 10 −1 < 1 + δ < 10 10 −2 < 1 + δ < 10 2 10 −3 < 1 + δ < 10 3 10 −4 < 1 + δ < The effect of truncating the 1-point PDF at different cutoffs in 1 + δ, shown for the log-normal (solid blue) and analytic (dashed red) PDFs. Progressively darker lines corresponding to the inclusion of larger underdensities and overdensities, from 10 times larger and smaller than the mean plasma mass to 10 4 times larger and smaller than the plasma mass, respectively. (Right) Sensitivity of the constraints to the choice of the PDF parameterization. Shown is our fiducial constraint with the log-normal PDF (solid red), as well as constraints with the inclusion of a linear bias b = 1.5 between the dark matter and baryons in the log-normal prescription (solid blue), the analytic PDF (solid green), the inferred PDF of underdensities in voids from Ref. [80] (solid purple), and a Gaussian PDF (dotted orange). ƭ Fig. 7 show the limits on obtained for γ → A oscillations and dark matter A → γ oscillations, respectively. The limits are qualitatively similar between the different PDFs in the relevant ranges of 1 + δ, providing some reassurance that our fiducial choice is reasonably conservative.

Baryonic Power Spectrum
The log-normal distribution for f (m 2 γ ; t) is fully characterized by two statistics: the mean, given by the homogeneous value m 2 γ , and a variance in log-space σ 2 LN ; the latter is defined through the usual variance of the baryon density fluctuations σ 2 as σ 2 LN ≡ log(1 + σ 2 ). Since σ 2 is directly proportional to the integral over the baryonic power spectrum (which can safely be taken to describe fluctuation in the free electron number density [25]), uncertainties in the baryonic power spectrum translate into uncertainties in f (m 2 γ ; t). To assess how significant these uncertainties are, we adopt two extremal prescriptions for the baryonic power spectrum, PS min and PS max , corresponding to reasonable lower and upper bounds; PS min leads to a narrower distribution f (m 2 γ ; t), while PS max leads to a broader one. These are described in detail in Ref. [25], and take into account typical uncertainties between different hydrodynamic simulations [38][39][40][41]. In the main letter, we always show the more conservative bound of the ones obtained with the two prescriptions. A maximum difference of O(15%) in the mixing parameter constraint is obtained between the two different power spectrum prescriptions, which is expected since σ 2 LN only has a log-dependence on the integral of the power spectrum.

PDF Tails
If the log-normal or the analytic PDF continues to be a good description out to larger upward or downward fluctuations in m 2 γ , we can extend the acceptable range of 1 + δ for these PDFs. In the top left panel of Fig. 6 and the left panel of Fig. 7 we show the effect of truncating the one-point PDF at different cutoffs in 1 + δ for the γ → A dark photon and A → γ dark photon dark matter cases, respectively. These are shown for the log-normal PDF (solid blue) and analytic PDF (dashed red). Progressively darker lines corresponding to the inclusion of larger underdensities and overdensities, from 10 times larger and smaller than the mean plasma mass to 10 4 times larger and smaller than the plasma mass, respectively.
Since the analytic PDF shows a strong cut-off in the probability of fluctuations below 10 −2 , extending the range of 1 + δ does not significantly improve our lower mass reach. For the log-normal distribution, however, an order of magnitude improvement in mass reach can be obtained with 10 −4 < 1 + δ < 10 4 as compared to our fiducial results. Both PDFs promise significant improvements at m A > 10 −13 eV. This underscores the fact that pinning down the PDF of plasma fluctuations at the tails can significantly improve on the fiducial constraints pre-

Additional Redshift Variations
We show in the bottom panel of Fig. 6 the effect of excluding conversions over redshift ranges different from the ones considered in the main Letter. Excising a larger redshift range of 6 < z < 30 (dashed orange) has no effect on the fiducial limits, obtained after excising 6 < z < 20. This shows that our results are robust to the details of reionization; in particular, an earlier onset to reionization as allowed by Planck 2018 data with a FlexKnot reionization parameterization [91] does not lead to any change in our limits.
Constraints relying solely on conversions before the deeply nonlinear cosmological epoch (z 20) are also shown in Fig. 6 for both the log-normal (solid blue) and analytic (dashed red) PDFs. Similar results are obtained with either prescription, as expected-the spectrum of fluctuations in the linear regime become increasingly Gaussian and are not subject to large systematics.
We further show constraints relying on conversions above z > 0.1 (solid green) and z > 1 (solid purple). The former has minimal impact on our fiducial constraints while the latter, in coördination with our cut on the allowed fluctuation size 10 −2 < 1 + δ < 10 2 , restricts conversions for the lowest masses accessible to our fiducial analysis.

Dependence of Limits on Smallest Scale
An understanding of the various scales at which fluctuations influence constraints from conversions in inhomogeneities is crucial. We show in Fig. 8 the constraints in the fiducial log-normal prescription as a function of maximum cutoff on the scale of perturbations, k max , for a few benchmark masses. In all cases, constraints stabilize around the baryon Jeans scale, k J ∼ 300 h Mpc −1 .
In the left panel, we show results for dark photon constraints from γ → A , shown for benchmark masses m A = 2 × 10 −15 , 10 −13 , and 10 −12 eV, respectively. The constraint in the homogeneous limit is shown for the latter two benchmark mass points, while for m A = 2 × 10 −15 eV no conversions are accessible in the homogeneous limit. In the right panel we show dark photon dark matter constraints from A → γ, shown for benchmark masses m A = 2 × 10 −14 , 10 −13 , and 10 −12 eV in solid red, blue, and green lines, respectively. The constraint in the homogeneous limit is shown for m A = 10 −13 eV, while no conversions are accessible in the homogeneous limit for the other two benchmark mass points shown.