Modeling Dark Photon Oscillations in Our Inhomogeneous Universe
Andrea Caputo, Hongwan Liu, Siddharth Mishra-Sharma, Joshua T. Ruderman
MModeling Dark Photon Oscillations in Our Inhomogeneous Universe
Andrea Caputo, ∗ Hongwan Liu,
2, 3, † Siddharth Mishra-Sharma, ‡ and Joshua T. Ruderman § Instituto de F´ısica Corpuscular, CSIC-Universitat de Valencia, Apartado de Correos 22085, E-46071, Spain Center for Cosmology and Particle Physics, Department of Physics,New York University, New York, NY 10003, USA Department of Physics, Princeton University, Princeton, NJ 08544, USA (Dated: April 28, 2020)A dark photon may kinetically mix with the Standard Model photon, leading to observable cos-mological signatures. The mixing is resonantly enhanced when the dark photon mass matchesthe primordial plasma frequency, which depends sensitively on the underlying spatial distribution ofelectrons. Crucially, inhomogeneities in this distribution can have a significant impact on the natureof resonant conversions. We develop and describe, for the first time, a general analytic formalismto treat resonant oscillations in the presence of inhomogeneities. Our formalism follows from thetheory of level crossings of random fields and only requires knowledge of the one-point probabilitydistribution function (PDF) of the underlying electron number density fluctuations. We validateour formalism using simulations and illustrate the photon-to-dark photon conversion probability forseveral different choices of PDFs that are used to characterize the low-redshift Universe. (cid:135)
I. INTRODUCTION
A dark photon A (cid:48) which kinetically mixes with theStandard Model (SM) photon, γ , is one of the sim-plest extensions of the SM [1]. The range of possible A (cid:48) masses m A (cid:48) spans many orders of magnitude, and anintense theoretical and experimental program is ongo-ing to constrain and test dark photon models. At lowmasses ( m A (cid:48) (cid:46) − eV), the Compton wavelength of m A (cid:48) starts to exceed the size of typical experiments, andterrestrial probes start to become increasingly insensitiveto the presence of A (cid:48) , motivating probes on larger lengthscales. Light dark photons in this mass range are alsoa well-motivated candidate for dark matter [2–12], whilerelativistic A (cid:48) produced by decaying dark matter whichthen resonantly convert into γ has also been proposed asa new-physics explanation [13, 14] and can be detectedby 21-cm observations. Probes of the dark photon overcosmological scales are therefore critical to constrainingits properties.Existing experimental measurements are sensitive tooscillations of cosmic microwave background (CMB) pho-tons into dark photons, γ → A (cid:48) , or to oscillations of darkphoton dark matter into low-energy photons, A (cid:48) → γ .The probability of these conversions at a particular red-shift z and position in space (cid:126)x depends on the photonplasma mass at that point, m γ ( z, (cid:126)x ), and becomes res-onantly enhanced whenever it becomes equal to m A (cid:48) . γ → A (cid:48) conversions can leave a distortion in the energyspectrum of the CMB due to a disappearance of photonsfrom the spectrum, while A (cid:48) → γ conversions for lowmass dark photons produce SM photons that are readily ∗ [email protected]; ORCID: 0000-0003-1122-6606 † [email protected]; ORCID: 0000-0003-2486-0681 ‡ [email protected]; ORCID: 0000-0001-9088-7845 § [email protected]; ORCID: 0000-0001-6051-9216 absorbed by baryons and electrons, resulting in an in-crease of the intergalactic medium (IGM) temperature.Under the assumption of a completely homogeneousUniverse, constraints on the kinetic mixing parameter (cid:15) for the case of γ → A (cid:48) were obtained using theCOBE/FIRAS [15] measurement of the CMB energyspectrum, which shows no significant evidence of distor-tion from a pure blackbody spectrum [16, 17]. More re-cently, Ref. [18] presented new homogeneous constraintsfor A (cid:48) → γ in the case of dark photon dark matter, find-ing strong limits on (cid:15) using IGM temperature measure-ments during HeII reionization, among other novel cos-mological constraints.This paper is part of a pair of companion papers withthe overarching goal of establishing a new formalism forunderstanding both the physics and the experimentalconsequences of γ → A (cid:48) and A (cid:48) → γ oscillations in ourinhomogeneous Universe. In Ref. [19], hereafter referredto as Paper I, we briefly introduce our formalism andpresent (i) the γ → A (cid:48) CMB spectral distortion and (ii) A (cid:48) → γ dark photon dark matter IGM tempera-ture constraints on the kinetic mixing parameter (cid:15) . Wefind that limits derived under the assumption of a homo-geneous photon plasma were not conservative, and thatincluding inhomogeneities allows for constraints to be setover a much broader mass range of A (cid:48) . In this paper, weprovide a detailed description of the formalism and itsmathematical derivation, as well as an elaboration on thecosmological inputs that go into the A (cid:48) limits obtainedin Paper I.This paper is organized as follows. We begin Sec. IIwith a quantum mechanical derivation of the oscillationprobability of γ ↔ A (cid:48) for both relativistic and nonrel-ativistic A (cid:48) with multiple resonance crossings. We thenintroduce our analytic formalism for computing the ex-pected probability of conversion for both γ → A (cid:48) and A (cid:48) → γ , taking as input the one-point probability densityfunction (PDF) of baryon inhomogeneities in our Uni-verse, described in Sec. III. In Sec. IV, we explore our a r X i v : . [ a s t r o - ph . C O ] A p r formalism in the regime where fluctuations are Gaussianto gain some analytic understanding. We then move onto describe the two main cosmological inputs that areneeded to apply our results to our Universe: the one-point probability density functions of baryon fluctuationsin Sec. V, and the variance of fluctuations (characterizedby power spectra for the number density fluctuation ofbaryons and free electrons), in Sec. VI. We validate ourformalism against several simulations of baryon fluctu-ations and γ ↔ A (cid:48) conversions, which we describe inSec. VII. Some results for the γ ↔ A (cid:48) conversion prob-ability obtained from our analytic formalism for variouscosmological inputs are presented in Sec. VIII. We fi-nally conclude in Sec. IX. In our appendices, we providea comparison between our work and several recent pa-pers treating inhomogeneities [20–22], along with otherdetails of the formalism.Throughout this work, we use natural units with (cid:126) = c = k B = 1, as well as the Planck ) pointing to the Jupyter notebooks used togenerate them. II. OSCILLATIONS γ ↔ A (cid:48) oscillations are described by the same formal-ism as neutrino flavor oscillations, which have been stud-ied extensively in the literature. In this section, we followthe neutrino discussion of Ref. [24] closely, first review-ing γ → A (cid:48) oscillations and highlighting any differencesbetween γ ↔ A (cid:48) and neutrino oscillations whenever theyarise. A (cid:48) → γ oscillations are similar, and are discussedat the end of this section.Consider a single photon passing through some world-line from the early Universe to us. Along this path,parametrized by t , there are variations in the numberdensities of free electrons and neutral atoms, leadingto variations in the plasma properties, giving rise to aplasma mass m γ ( t ) [16]: m γ ( t ) (cid:39) πα EM n e ( t ) m e − ω ( t ) ( n HI ( t ) − (cid:39) . × − eV (cid:18) n e ( t )cm − (cid:19) − . × − eV (cid:18) ω ( t )eV (cid:19) (cid:18) n HI ( t )cm − (cid:19) . (1)Here, α EM is the electromagnetic fine structure constant, m e is the electron mass, n HI is the refractive index ofmonatomic hydrogen [25], ω ( t ) is the photon energy, and n e ( t ) and n HI ( t ) are the local free electron and neutralhydrogen densities along the path. We similarly define m γ as the homogeneous value of m γ , evaluated with the Our expression clarifies the actual species densities that enter the mean cosmological values of n e and n HI . We neglect he-lium, which makes up only 8% by number density andhas a smaller index of refraction. If fluctuations in freeelectron density x e are small, i.e. , x e has essentially thesame value everywhere in space at each point in time,then fluctuations in n e and n HI track fluctuations in thenumber density of baryons, n b , m γ ( t ) m γ ( t ) = n b n b . (2)Further discussion of this proportionality and the effectof fluctuations in x e can be found in Sec. VI. Figure 1illustrates the variation of the photon plasma mass as afunction of redshift for several representative values ofthe present-day photon frequency ω .The kinetic mixing between γ and A (cid:48) induces oscilla-tions between these two interaction eigenstates, describedby the Schr¨odinger equation i dd t (cid:18) γA (cid:48) (cid:19) = H (cid:18) γA (cid:48) (cid:19) , (3)where H is the Hamiltonian (assuming all particles arerelativistic) [24] H = 14 ω ( t ) (cid:18) m γ ( t ) − m A (cid:48) (cid:15)m A (cid:48) (cid:15)m A (cid:48) − m γ ( t ) + m A (cid:48) (cid:19) . (4)The plasma mass is an in-medium effect similar to theMikheyev-Smirnov-Wolfenstein (MSW) effect in the caseof neutrino oscillations [26, 27], leading in our case tothe familiar correction of m γ / ω relative to the propa-gation phase of a massless particle [28, 29]. H can beconveniently written in terms of Pauli matrices, H = φ ( t ) σ + η ( t ) σ , (5)where φ ( t ) ≡ m γ ( t ) − m A (cid:48) ω ( t ) , η ( t ) = (cid:15)m A (cid:48) ω ( t ) . (6) φ ( t ) has the intuitive interpretation of being half the rel-ative phase between γ and A (cid:48) .Starting with an initial state of γ , the Schr¨odingerequation can be solved perturbatively in (cid:15) , with η ( t ) σ asan interaction Hamiltonian. To first order in (cid:15) , we obtain A (cid:48) ( t ) = − ie iα ( t ) (cid:90) t d ξ η ( ξ ) e − iα ( ξ ) + O ( (cid:15) ) , (7)where we have defined α ( s ) ≡ (cid:90) s d ξ φ ( ξ ) , (8) plasma mass expression in Ref. [16], and corrects earlier expres-sions for the photon mass, which mistakenly used the refractiveindex of diatomic hydrogen gas and not monatomic hydrogen. the accumulated phase between 0 and s . This leads tothe probability of disappearance at t , given by | A (cid:48) ( t ) | ,or explicitly, P γ → A (cid:48) ( t ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) t d ξ η ( ξ ) e − iα ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) + O ( (cid:15) ) . (9)Away from regions of space where m γ ( t ) ∼ m A (cid:48) , φ ( t )is given parametrically by φ ( t ) ∼
200 kpc − z ( t ) (cid:32) (cid:12)(cid:12) m γ ( t ) − m A (cid:48) (cid:12)(cid:12) − eV (cid:33) (cid:18) ω /T CMB , (cid:19) , (10)where T CMB , is the temperature of the CMB today, and ω ( t ) = ω (1 + z ( t )), with z ( t ) being the cosmological red-shift at t . The FIRAS experiment detects photons in therange 1 . (cid:46) ω /T CMB , (cid:46) . φ ( t ) therefore oscillates rapidlywith t , except when m γ ( t ) ∼ m A (cid:48) ; we can therefore eval-uate the integral in Eq. (9) over the entire worldline of thephoton using the stationary phase approximation, giving P γ → A (cid:48) = π (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) i η ( t i ) (cid:112) | φ (cid:48) ( t i ) | e − iα ( t i ) e iβ i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + O ( (cid:15) ) , (11)where i indexes positions t i where m γ ( t i ) = m A (cid:48) , φ (cid:48) isthe derivative of φ , and β i = ± π/
4, with the sign givenby the sign of φ (cid:48) ( t i ).If there is only one point t r where m γ ( ξ ) = m A (cid:48) , thenoscillations from γ to A (cid:48) occur resonantly at t r , giving P γ → A (cid:48) (cid:39) πη ( t r ) | φ (cid:48) ( t r ) | = π(cid:15) m A (cid:48) ω ( t r ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ ( t )d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − t = t r . (12)This is the same expression derived using the Landau-Zener approximation for non-adiabatic transitions of γ → A (cid:48) in Ref. [16]; it is also similar to expressions for the A (cid:48) production rate in stars under the narrow width approx-imation [30, 31]. The stationary phase approximationhas also been used to calculate appearance and disap-pearance probabilities with resonant oscillations in thecontext of neutrino oscillations [24] and in axion-photonconversions in magnetic fields [28, 29, 32].When multiple resonances exist, the probability be-comes P γ → A (cid:48) (cid:39) (cid:88) i πη ( t i ) | φ (cid:48) ( t i ) | + (cid:88) i 10 kpc. Second, the size of fluctuations in m γ and hence φ is determined by the standard deviation ofbaryon density fluctuations σ b (see Fig. 2 for some typicalvalues of these fluctuations); in other words, we expectthat between resonances when m γ = m A (cid:48) , fluctuationsin m γ can typically reach values of around (1 ± σ b ) m A (cid:48) .These two estimates and Eq. (10) show that the phasedifference between consecutive resonances is roughly θ ( t i , t i +1 ) (cid:38) × (1 + z h ) (cid:18) R J 10 kpc (cid:19) min [ σ b ( z h ) , × (cid:16) m A (cid:48) − eV (cid:17) (cid:18) ω /T CMB , (cid:19) , (14)where z h is the lowest redshift at which m γ = m A (cid:48) , sinceall of the resonances occur in a redshift window centeredat redshifts where m γ = m A (cid:48) , and transitions are moreadiabatic at lower redshifts. The relative phase between A (cid:48) s produced at any two resonance points is thus manytimes larger than 2 π throughout the history of the Uni-verse, so that cos θ ( t i , t j ) is expected to be uncorrelatedwith the location of the resonances. We will ultimatelybe interested in the mean value of P γ → A (cid:48) across all pos-sible photon worldlines, such that uncorrelated interfer-ence effects average out. We therefore do not expect thesecond summation in Eq. (13) to contribute to the overallprobability of conversion, obtained by averaging over allworldlines, each with a different distribution of resonancepoints. The total probability of oscillations is thus ob-tained by summing up the conversion probability of eachresonance, each given by the Landau-Zener expression: P γ → A (cid:48) (cid:39) (cid:88) i π(cid:15) m A (cid:48) ω ( t i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ ( t )d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − t = t i . (15)We note that the Landau-Zener approximation holds forany crossing encountered by the photon. For (cid:15) (cid:28) τ res ∼ (cid:15) | d ln m γ / d t | − , ismuch smaller than the timescale over which m γ changes, | d ln m γ / d t | − , at any crossing, allowing the use of theLandau-Zener approximation of taking the density pro-file over the resonance to be linear. The suitability ofthe Landau-Zener approximation in the context of neu-trino oscillations, starting from a similar Hamiltonian toEq. (4), is derived in Ref. [35].Following a similar derivation, we can show that rel-ativistic dark photons undergoing A (cid:48) → γ conversionswill also have a conversion probability that is identical − − z − − − − − − − − − − m γ [ e V ] Photon plasma mass ω = 10 − eV ω = 5 × − eV ω = 10 − eV ω = 5 × − eVGaussian σ m γ Log-normal 1- σ containment FIG. 1. The photon plasma mass as a function of redshift forseveral values of the present-day photon energy ω . The Gaus-sian standard deviation of plasma mass fluctuations σ m γ , in-formed by the linear baryon power spectrum for illustration,is shown as the blue band. The equivalent middle-68% con-tainment of fluctuations assuming a log-normal description ofthe PDF is shown as the red band. to Eq. (15), with ω now specifying the A (cid:48) energy. If A (cid:48) isthe dark matter, however, the assumption of relativisticparticles assumed in Eq. (4) breaks down. Nevertheless,there are several ways to see that the conversion prob-ability P A (cid:48) → γ is identical to P γ → A (cid:48) with ω ( t i ) → m A (cid:48) .First, it can be derived in thermal field theory [4] by ap-plying a narrow-width approximation (see App. C). Sec-ond, the probability of conversion P A (cid:48) → γ , shown on theright-hand side of Eq. (15), is Lorentz invariant, as alltransition probabilities should be. Evaluating the prob-ability in the rest frame of the dark matter A (cid:48) gives P A (cid:48) → γ (cid:39) (cid:88) i π(cid:15) m A (cid:48) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ ( t )d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − t = t i , (16)consistent with the result in Ref. [4]. Under standard cos-mology scenarios where the magnitude of δ b grows mono-tonically with redshift, each value of m A (cid:48) has at most oneresonance transition point; our formalism, however, doesnot rely on this assumption.The results in Eqs. (15) and (16) form the startingpoint for understanding γ ↔ A (cid:48) conversions along a sin-gle worldline, as well as for all of the results presentedin Paper I. III. FORMALISM In the presence of inhomogeneities, the resonance con-dition can be met many times along a path, even at timeswhen the homogeneous plasma mass m γ is far from m A (cid:48) and no resonance is present in the homogeneous limit. Each worldline passes through a different series of per-turbations, leading to conversions that vary significantlyin number and in distance from the observer. A. γ → A (cid:48) oscillations We will now discuss how to determine the expectedprobability of γ → A (cid:48) conversion, (cid:104) P γ → A (cid:48) (cid:105) . The deriva-tion of our results is closely related to the derivation ofthe mean number of times a stationary process crosses afixed level per unit time [36, 37].To average over all worldlines, we first begin by rewrit-ing the probability of conversion along a worldline asd P γ → A (cid:48) d t = πm A (cid:48) (cid:15) ω ( t ) δ D ( m γ ( t ) − m A (cid:48) ) m γ ( t ) , (17)where δ D is the Dirac delta function. We can check thatEq. (15) is recovered by performing the substitutiond t = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − d m γ m γ (18)and integrating the delta function over the entire world-line. The mean value of P γ → A (cid:48) is then obtained by inte-grating over all possible values of m γ at each point alongthe path, weighted by the probability density function(PDF) f ( m γ ; t ) of m γ :d (cid:104) P γ → A (cid:48) (cid:105) d z = πm A (cid:48) (cid:15) ω ( t ) (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) × (cid:90) d m γ f ( m γ ; t ) δ D ( m γ − m A (cid:48) ) m γ . (19)Note that the PDF evolves with time since m γ tracks thebaryon density (in the limit of small fluctuations in thefree electron fraction), as shown in Eq. (2). We can nowperform the integral to gived (cid:104) P γ → A (cid:48) (cid:105) d z = πm A (cid:48) (cid:15) ω ( t ) (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) f ( m γ = m A (cid:48) ; t ) . (20)The problem of determining the averaged probabilitytherefore reduces to finding the PDF of m γ , which we dis-cuss in detail in subsequent sections. Note that Eqs. (19)and (20) both apply equally to relativistic A (cid:48) → γ oscil-lations as well.As an example, let us consider the homogeneous limitwhere m γ = m γ everywhere; in this case, the PDF istrivially given by f h ( m γ ; t ) = δ D ( m γ − m γ ( t )) . (21)We therefore see that the mean homogeneous conversionprobability is (cid:104) P γ → A (cid:48) (cid:105) h = (cid:90) d t πm A (cid:48) (cid:15) ω ( t ) δ D ( m γ − m γ ( t ))= (cid:88) i πm A (cid:48) (cid:15) ω ( t i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ ( t )d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − t = t i , (22)where i indexes times t i when m γ ( t i ) = m A (cid:48) , and we haveagain made use of the substitution shown in Eq. (18).This recovers the homogeneous limit expressions foundin Eq. (15) and Ref. [16]. B. A (cid:48) → γ oscillations For A (cid:48) → γ conversions with A (cid:48) dark matter, in therange of m A (cid:48) of interest, the converted photons are ab-sorbed quickly by electrons in the IGM via free-free ab-sorption [18], heating the IGM. The quantity of in-terest is therefore the average energy injected into theplasma per baryon, (cid:104) E A (cid:48) → γ (cid:105) . The derivation of (cid:104) E A (cid:48) → γ (cid:105) proceeds in a similar fashion, except that the energyinjected per volume along the worldline is given by P A (cid:48) → γ ( t ) ρ A (cid:48) ( t ), where ρ A (cid:48) ( t ) is the mass density of A (cid:48) dark matter at the point on the worldline t . The rateof energy injected per baryon along the worldline of themassive dark photon is therefored E A (cid:48) → γ d t = πm A (cid:48) (cid:15) ρ A (cid:48) n b ρ A (cid:48) ( t ) ρ A (cid:48) ( t ) δ D ( m γ ( t ) − m A (cid:48) ) m γ ( t ) , (23)where n b is the homogeneous baryon number density,with ρ A (cid:48) /n b being a time-independent quantity.To obtain the mean value, we technically need to per-form an integral over the joint distribution of both m γ and ρ A (cid:48) . However, two points make this unnecessary.First, if fluctuations in the free electron fraction aresmall, then as we argued in Eq. (2), m γ ∝ n b . This as-sumption is true during the period of HeII reionization,the regime we study in Paper I to obtain limits on (cid:15) inthe case of A (cid:48) dark matter, since the Universe is almostcompletely ionized at this time except for HeII, whilefluctuations in baryon density are large compared to themean. Second, we adopt the standard assumption thatbaryon density fluctuations track matter density fluctu-ations ρ m with a bias b ∼ O (1). This means that ρ m ρ m ( t ) = 1 b n b n b = 1 b m γ m γ ( t ) , (24)where in the case of A (cid:48) dark matter, ρ m (cid:39) ρ A (cid:48) . Note thatin Paper I, we assumed b = 1 for simplicity, althoughincluding a small bias consistent with values reported inRef. [38] does not change the result significantly. Withthis relation, we findd (cid:104) E A (cid:48) → γ (cid:105) d z = πm A (cid:48) (cid:15) ρ A (cid:48) b n b (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) × (cid:90) d m γ m γ m γ ( t ) f ( m γ ; t ) δ D ( m γ − m A (cid:48) ) m γ , (25)and as before we can perform the integral to obtaind (cid:104) E A (cid:48) → γ (cid:105) d z = πm A (cid:48) (cid:15) m γ ( t ) ρ A (cid:48) b n b (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) f ( m γ = m A (cid:48) ; t ) . (26) This treatment implicitly assumes that the conversionprobability of A (cid:48) → γ is small, which is required if A (cid:48) is all of the dark matter. A more general treatment ispossible by allowing b ( z ) to vary as a function of thetotal conversion up to z .In deriving Eq. (25), we have assumed that the en-ergy deposited by the conversion is distributed uniformlyacross all baryons, enabling us to characterize the entireplasma with a single temperature. This is in contrast tothe assumption made in Ref. [22], where energy deposi-tion is local. The corresponding expression under thisassumption can be obtained by replacing n b → n b insidethe integral,d (cid:104) E A (cid:48) → γ (cid:105) local d z = πm A (cid:48) (cid:15) ρ A (cid:48) b n b (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) × (cid:90) d m γ f ( m γ ; t ) δ D ( m γ − m A (cid:48) ) m γ , (27)which we can integrate to obtaind (cid:104) E A (cid:48) → γ (cid:105) local d z = πm A (cid:48) (cid:15) ρ A (cid:48) b n b (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) f ( m γ = m A (cid:48) ; t ) . (28)These results agree with the analogous expression inRef. [22]. We leave a detailed comparison of our resultsto App. A.Eqs. (19) and (25) were presented in Paper I, and withseveral different choices of the PDF, f ( m γ ; t ), were usedto derive all of the relevant bounds on the existence on A (cid:48) . The rest of the paper will now focus on determiningthe analytic form of f ( m γ ; t ), and checking these resultswith simulation. IV. UNDERSTANDING THE FORMALISM We are now in a position to evaluate Eqs. (19) and (25)numerically. To gain some intuition regarding our for-malism and highlight some important physics, we beginour discussion assuming Gaussian fluctuations, a validassumption at redshifts z (cid:29) 20, where density perturba-tions are well described by linear perturbation theory. Inthis limit, (cid:104) P γ → A (cid:48) (cid:105) and (cid:104) E A (cid:48) → γ (cid:105) have analytic solutions,which serve as a useful pedagogical example for our fulltreatment. We will first discuss the various inputs thatdetermine f ( m γ ; t ), before discussing the analytics of theresult in the Gaussian regime. A. PDF, variance of fluctuations and powerspectrum We begin by taking the limit where we neglect fluctu-ations in the free electron fraction, as in Eq. (2). Thebaryon density fluctuation δ b ( (cid:126)x ) at each point in spaceis defined as δ b ( (cid:126)x ) ≡ ρ b ( (cid:126)x ) − ρ b ρ b , (29)where ρ b ( (cid:126)x ) is the baryon mass density at (cid:126)x and ρ b is themean, homogeneous baryon mass density. In the linearregime, the fluctuations follow a Gaussian distribution,given by the one-point PDF of baryon density fluctua-tions, P G ( δ b ; z ) = 1 (cid:112) πσ ( z ) exp (cid:18) − δ σ ( z ) (cid:19) , (30)with the variance of the distribution σ directly relatedto the baryon (auto) power spectrum, P bb ( k, z ) through σ ( z ) = (cid:90) d (cid:126)k (2 π ) P bb ( k, z ) . (31)In linear perturbation theory, P bb is the linear baryonpower spectrum, P bb,L ( k, z ). Fig. 2 shows σ b ( z ),computed using the value of P bb,L ( k, z ) produced byCLASS [39]. With this function, we have fully specifiedthe one-point PDF: f ( m γ ; t ) = d δ b d m γ P ( δ b ; t ) = P ( δ b ( m γ ); t ) m γ ( t ) , (32)directly relating the PDF for m γ to a cosmological ob-servable. We discuss the issue of perturbations in x e inSec. VI. The blue band in Fig. 1 shows the standard de-viation of plasma mass fluctuations induced by baryonGaussian fluctuations, for illustration. B. Jeans scale and sensitivity to small scales In linear perturbation theory, the linear matter powerspectrum P mm,L ( k, z ) scales as k − at large k , so that thevariance in matter fluctuations, calculated using Eq. (31)with P mm,L ( k, z ), theoretically exhibits a log k ultravioletdivergence. This divergence is regulated by the fact thatmeasurements and simulations of matter density are al-ways averaged over some smoothing scale R ; P mm,L ( k, z )needs to be convolved with a windowing function ( e.g. ,a top-hat function) with characteristic size R , giving avariance as a function of R . For baryons in the linearregime, baryonic structures have the Jeans length as aphysical cut-off scale: the formation of structures withcomoving size less than R J is suppressed due to gas pres-sure counteracting the gravitational collapse, defined by R J ( z ) = 2 √ π √ z ) H ( z ) (cid:115) γT b ( z ) µm p , (33)where γ = 5 / µ = 1 . 22 is the mean molecular weight of − z − − − − σ b Variance of linear fluctuations FIG. 2. Standard deviation of baryon fluctuations σ b in linearperturbation theory (red). The dashed line indicates wherethe typical size of fluctuations becomes comparable to themean density. the neutral IGM, m p is the proton mass, T b is the baryontemperature, c s ( z ) is the baryon sound speed, and H ( z )is the Hubble parameter. Numerically, this is R J ( z ) ∼ . (cid:18) . 01 + z (cid:19) / (cid:18) T b K (cid:19) / , (34)with a minimum value of R J,min ∼ − Mpc at z ∼ 20 with T b ∼ 10 K, before reionization heats baryonssignificantly. In terms of wavenumber, the Jeans lengthensures that P bb,L ( k, z ) is suppressed above k J ∼ π/R J ,which lies between 10 and 10 Mpc − for z (cid:38) k J de-creases rapidly due to the increase in baryon temper-ature. Fluctuations also become increasingly nonlinearduring this epoch. On the other hand, Boltzmann codeslike CLASS [39] and CAMB [40] compute the linearbaryon power spectrum P bb,L ( k, z ) with a suppression at k J without reionization sources included when computing T b , leading to a suppression scale of k J ∼ h Mpc − ,instead of k J ∼ h Mpc − as estimated from Eq. (34).However, power above k J ∼ h Mpc − is actually un-suppressed due to the increasingly nonlinear behavior ofbaryons at late times; this lack of suppression is con-firmed by baryon power spectra extracted from high-resolution hydrodynamic N -body simulations with bary-onic physics included [41]. In light of this, we continue toadopt the linear power spectrum computed by CLASS for P bb,L with power suppressed above roughly 700 h Mpc − ,and defer a complete discussion of this to Sec. VI. Wewill also refer to the Jeans scale and corresponding Jeanslength as the value of k at which the linear power spec-trum of CLASS shows a suppression of power relative tothe matter power spectrum, instead of Eq. (33).Since the baryon power spectrum P bb,L like P mm,L alsoscales as approximately k − at large k up to k J , and non-linear effects usually lead to the baryon power spectrum P bb exceeding P bb,L at large k , the value of σ and hencethe probability of conversion is sensitive to the small-est unsuppressed length scales in P bb ( z ). This exhibitsone of the key peculiarities of dark photon oscillations inthe presence of inhomogeneities: the resulting physics issensitive to small-scale perturbations, depending on thedetails of the baryon power spectrum at scales as smallas Mpc − , providing a rare example of a cosmologi-cal phenomenon that is ultraviolet-sensitive to perturba-tions. We will discuss our treatment of the baryon powerspectrum beyond the linear regime in significant detail inSec. VI.Finally, although the Gaussian distribution is well-motivated at high redshifts when fluctuations are small,the Gaussian PDF shown in Eq. (30) breaks down once σ b ∼ 1, since large negative fluctuations which lead to anoverall negative density is assigned a sizable probability.Fig. 2 shows that the applicability of the Gaussian PDFstarts becoming questionable once z (cid:46) C. Analytics Substituting the expression for f ( m γ ; t ) in Eq. (32)into Eq. (16) givesd (cid:104) P γ → A (cid:48) (cid:105) G d z = πm A (cid:48) (cid:15) m γ ( z ) ω ( z ) (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) × (cid:112) πσ ( z ) exp (cid:34) − ( m A (cid:48) /m γ ( z ) − σ ( z ) (cid:35) , (35)where the subscript ‘G’ stands for Gaussian. The corre-sponding energy deposited per baryon isd (cid:104) E A (cid:48) → γ (cid:105) G d z = πm A (cid:48) (cid:15) m γ ( z ) m A (cid:48) m γ ( z ) ρ A (cid:48) bn b (cid:12)(cid:12)(cid:12)(cid:12) d t d z (cid:12)(cid:12)(cid:12)(cid:12) × (cid:112) πσ ( z ) exp (cid:34) − ( m A (cid:48) /m γ ( z ) − σ ( z ) (cid:35) . (36)Given σ b ( z ) and m γ ( z ) from Eq. (1), these compactresults can now be integrated numerically to obtain (cid:104) P γ → A (cid:48) (cid:105) .In the σ → σ , the characteristicredshift width ∆ z over which transitions occur is givenby ∆ z ∼ σ b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ d t d t d z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − , (37)which during periods when x e does not change signifi-cantly ( e.g. , before recombination, during the dark ages and after reionization is complete) is approximately∆ z ∼ . (cid:18) z h (cid:19) (cid:18) σ b ( z h )0 . (cid:19) , (38)where z h is the redshift at which m γ = m A (cid:48) . In thelinear regime, fluctuations grow linearly with the scalefactor, and thus σ b ∝ / (1+ z ); this implies that ∆ z staysrelatively constant throughout the dark ages. We can seethat the range of redshifts over which conversions canhappen can be very large, with ∆ z (cid:38) z at low redshifts.Similarly, a range of m A (cid:48) can now convert with signif-icant probability at any given redshift z . At a particularvalue of z h , this range is roughly∆ m A (cid:48) ∼ ± σ b m γ ( z h ) . (39)We note that when σ b exceeds one at z (cid:46) 20, this rangeof m A (cid:48) includes negative values, highlighting the fact thatthe Gaussian PDF becomes unphysical in this range, aswe discussed above. However, the lesson here is clear: thepresence of under- and overdensities allows conversionswell above and below the homogeneous value m γ ( z h ), al-lowing (i) conversions with m A (cid:48) (cid:46) − eV, i.e. , belowthe homogeneous plasma mass at any point in the his-tory of the Universe, and (ii) lower redshift conversionsfor 10 − eV (cid:46) m A (cid:48) (cid:46) − eV, which have a higherprobability of conversion.In the Gaussian limit, we can derive the ratio of theprobability calculated under the homogeneous assump-tion to the probability given a Gaussian PDF analyti-cally. We begin by defining the variable ∆ ≡ m A (cid:48) /m γ − (cid:104) P γ → A (cid:48) (cid:105) G = (cid:90) ∆ − d∆ g (∆) (cid:112) πσ exp (cid:18) − ∆ σ (cid:19) , (40)where we have defined g (∆) ≡ πm A (cid:48) (∆ + 1) (cid:15) ω (∆) d t d∆ , (41)and ∆ = m A (cid:48) /m γ ( z = 0) − 1, with ω now being afunction of ∆. Observe that g (0) = (cid:104) P γ → A (cid:48) (cid:105) h provided∆ ≥ 0, where (cid:104) P γ → A (cid:48) (cid:105) h is the homogeneous conversionprobability. Since the contribution to the integral is cen-tered at ∆ = 0, we can set g (∆) ≈ g (0) + g (cid:48) (0)∆ andtake σ b to be constant, giving (cid:104) P γ → A (cid:48) (cid:105) G (cid:104) P γ → A (cid:48) (cid:105) h (cid:39) (cid:34) erf (cid:32) (cid:112) σ (cid:33) + erf (cid:32) ∆ (cid:112) σ (cid:33)(cid:35) + g (cid:48) (0) g (0) σ b √ π (cid:20) exp (cid:18) − σ (cid:19) − exp (cid:18) − ∆ σ (cid:19)(cid:21) (42)for ∆ > 0. In the limit of constant x e and a mat-ter dominated Universe, g (cid:48) (0) /g (0) = 5 / 6. The ratioof probabilities would be greater than one if the homo-geneous assumption is conservative with respect to the − − − − m A [eV] − h P γ → A i G / h P γ → A i h Homogeneous vs. Gaussian Analytic estimateCalculated FIG. 3. An analytic estimate for the ratio of the probabilityof conversion with the Gaussian PDF v.s. that of the homoge-neous assumption (red) for conversions that happen at z < m A (cid:48) (cid:46) × − eV. Gaussian result. Moreover, in the limit when σ b → m A (cid:48) . The analytic estimate in Eq. (42) isevaluated with σ b at the homogeneous resonance redshift z h , and is shown for homogeneous conversions that occurat z < 6. We also include the exact probability ratio com-puted numerically. Large values of the ratio of Gaussianto homogeneous conversion probabilities occur for valuesof m A (cid:48) where the homogeneous limit resonance is deep inthe dark ages, but overdensities allow for significant con-versions at z ∼ m A (cid:48) , the Gaussian and homogeneousconversion probabilities rapidly converge as the varianceof fluctuations decreases.This ratio is significantly less than one for later conver-sions, i.e. , lighter m A (cid:48) . Qualitatively, the Gaussian PDFspreads out the probability of conversion over a range∆ z given in Eq. (38) compared to the homogeneous as-sumption; for small z h , this can mean that most of theprobability of conversion lies in the future, even though z h > 0. For sufficiently large z h , however, the probabil-ity of conversion in the future is negligible while the to-tal conversion probability in the Gaussian limit is larger,since conversions happening below z h have higher valuesof d P/ d z , increasing the overall integrated probability. D. Main takeaways Having gone through the example of a Gaussian PDF,we are now ready to understand how to arrive at a numer- ical result for (cid:104) P γ → A (cid:48) (cid:105) and (cid:104) E A (cid:48) → γ (cid:105) in general. We needtwo inputs, both of which need to be evaluated correctlyin the nonlinear regime:1. Functional form for baryon one-point PDF .In the linear regime, the PDF has a Gaussian form,but outside of the linear regime ( z (cid:46) The variance of baryon fluctuations . Whilethe mean of the PDF is fixed to be zero by the factthat the average baryon density must be the ho-mogeneous baryon density, the variance is not de-termined. The variance of baryon fluctuations willultimately be determined by the power spectrum ofmatter or baryons as a function of redshift. Out-side of the linear regime, one can no longer rely onBoltzmann codes to calculate these power spectra,and must instead make use of results informed by N -body simulations to obtain this information. Inall cases, the variance is ultraviolet-sensitive to thepower spectrum at small scales, but this UV sen-sitivity is cut off by the Jeans scale, k J , which weobtain from the CLASS linear baryon power spec-trum.We will devote Sec. V to examining more realistic alter-native one-point PDFs to the Gaussian, and Sec. VI to adiscussion of how to obtain the variance of the PDF deepin the nonlinear regime. V. ONE-POINT PROBABILITY DENSITYFUNCTIONS Table I shows a summary of all of the baryon one-pointPDFs considered in this paper and in Paper I, and Fig. 4shows a plot of these PDFs at a range of redshifts. Be-yond the linear regime, the log-normal PDF has beenproposed as a phenomenological fit to the total matterdistribution [42, 43] for both observations [38, 44–46] and N -body simulations [47–49]. The introduction of a biasparameter to the log-normal distribution has also beenshown to produce good fits phenomenologically [38, 46].There has also been a significant effort to calculate thematter one-point PDF from first principles [50, 51] withthe linear regime as a starting point, especially using apath-integral approach [52–55]. Finally, the study of cos-mic voids has shed some light on the underdense tail ofthe one-point PDF [56–61], and simulation results can beturned into a reasonable PDF at low densities. In this paper, we use only PDFs with functional forms that arefully defined by the mean and variance. Higher order statisticscould play an important role in a full characterization of baryonfluctuations. PDF Equation Power spectrum Remarks Log-normal(fiducial) P LN ( δ b ; z ) (1+ δ b ) − √ π Σ ( z ) exp (cid:104) − [ln(1+ δ b )+Σ ( z ) / ( z ) (cid:105) Nonlinear baryon Σ ( z ) = ln[1 + σ ( z )].Analytic P an ( δ b ; z ) ˆ C ( δ b ) (cid:113) πσ R J ( z ) exp (cid:20) − F ( δ b )2 σ R J ( z ) (cid:21) Linear matter, smoothed overbaryon Jeans length R J ˆ C ( δ b ) and F ( δ b ) defined in App. B.Log-normalwith bias P b LN ( δ b ; z ) b P LN (cid:16) δ b b ; z (cid:17) Nonlinear matter, withbaryon Jeans scale cut-off We adopt b = 1 . P voids ( δ b ; z ) φ voids ( z ) g voids (1 + δ b ; z ) – φ voids is the fractional volume of thesimulation in a void, g voids is the PDFof the mean 1 + δ b in voids [61, 62].Only used for underdensities.TABLE I: Summary of the baryon one-point PDFs used in this paper and in Paper I, with their defining equations and inputpower spectra used to determine the variance of these PDFs, where applicable. To understand γ ↔ A (cid:48) oscillations, we need a PDFthat is able to: (i) capture baryonic effects, and notjust the overall matter distribution; (ii) capture the dis-tribution of large overdensities and underdensities cor-rectly, and (iii) capture the behavior of baryonic fluctu-ations down to the Jeans scale of k ∼ – 10 Mpc − .Existing studies of the one-point PDF cannot meet allthree of these criteria simultaneously: first-principle, an-alytic results only apply to cold dark matter and donot account for baryonic effects, while the log-normalphenomenological fits have only been applied to simu-lations or data that have an effective smoothing scalemuch larger than the Jeans scale. Almost all results arevalidated with observations in a small range of densityfluctuations (10 − (cid:46) δ b (cid:46) e.g. , voids). These uncertainties surrounding thedistribution of baryonic fluctuations make it a challengeto arrive at a rigorous conclusion regarding constraintson γ ↔ A (cid:48) oscillations.Our approach is to adopt several independent modelsof the baryonic one-point PDF, in an attempt to capturethe systematic uncertainties discussed here. In our fidu-cial approach, we adopt a log-normal functional form forthe one-point PDF, with the variance of this distributiondetermined by baryonic power spectra obtained from acombination of different hydrodynamic N -body simula-tion results, which we detail in Section VI B. We truncatethe PDF to the range 10 − ≤ δ b ≤ to avoid thelarge uncertainties in the tails of the PDF. Our secondapproach relies on analytic results described in Ref. [55],which takes as input the linear matter power spectrumand computes the one-point PDF for matter fluctuationsas a function of redshift due to spherical collapse, whichwe then take to be equal to the baryon one-point PDF.We find that at low redshifts, these two approaches leadto similar PDFs in the range 10 − (cid:46) δ b (cid:46) at z = 0, as shown in Fig. 4; this range decreases to 10 − (cid:46) δ b (cid:46) 10 at z = 6. Restricting the PDFs tothe range 10 − (cid:46) δ b (cid:46) , the constraints on (cid:15) de-rived from γ ↔ A (cid:48) in Paper I differ by at most a factor ofapproximately three at m A (cid:48) ∼ − eV between our twoprescriptions, suggesting that we have reasonable controlover the uncertainties on the baryon PDF.In addition to the log-normal PDF and the analyticallyderived PDF, we also use two other PDFs as cross checksto our results. First, we use a log-normal distributionwith a bias parameter b , with the variance of the distri-bution given by the nonlinear matter power spectrum.This approach models the baryonic fluctuations as sim-ply a factor b times the overall matter fluctuations, givingus an estimate of how reliant we are on baryonic physicsmodeled by the simulations we used to obtain the bary-onic power spectrum for our fiducial log-normal PDF.Second, we use results from Ref. [61] for the probabilitydistribution of finding voids of a certain volume with acertain underdensity in their simulations, and constructa PDF of underdensities to test the underdense tails ofour PDFs. Both of these cross checks show that the con-straints we derive in Paper I are likely to be robust todifferences in systematics in the PDFs, and may improveif we can trust these PDF distributions to much largerunderdense and overdense fluctuations. A. Log-normal PDF Our fiducial choice for the PDF in this paper is thelog-normal PDF P LN ( δ b ; z ), given by P LN ( δ b ; z ) = (1 + δ b ) − (cid:112) π Σ ( z ) × exp (cid:18) − [ln(1 + δ b ) + Σ ( z ) / ( z ) (cid:19) , (43)0 − − − P ( δ b ) = m γ f ( m γ ) z = 0 . z = 0 . z = 0 . z = 0 . z = 0 . z = 3 . z = 3 . z = 3 . z = 3 . z = 3 . Probability density functions z = 6 . z = 6 . z = 6 . z = 6 . z = 6 . − − − − − − P ( δ b ) = m γ f ( m γ ) z = 20 z = 20 z = 20 z = 20 z = 20 − − − δ b = m γ /m γ z = 50 z = 50 z = 50 z = 50 z = 50 − − − z = 150 z = 150 z = 150 z = 150 z = 150 Log-normal(fiducial)LN b = 1 . FIG. 4. One point PDFs P ( δ b ; z ) at six different redshifts. The fiducial log-normal P LN (red), analytic P an (green) PDFs,the log-normal PDF with bias b = 1 . P . (blue), the PDF constructed from a model of voids P void (purple) [61], and theGaussian PDF P G (orange). Also shown are the fiducial 10 − < δ < boundaries (dashed gray). with Σ ( z ) = ln[1 + σ ( z )] as defined in Eq. (31). Thevariable ln(1+ δ b ) has a Gaussian distribution with mean − Σ / . As an immedi-ate consequence, unphysical fluctuations of δ b < − δ b . With thischoice of Σ, P LN satisfies (cid:90) ∞− d δ b P LN ( δ b ; z ) = 1 , (44) (cid:90) ∞− d δ b δ b P LN ( δ b ; z ) = 0 , (45) (cid:90) ∞− d δ b δ P LN ( δ b ; z ) = σ ( z ) , (46) i.e. , P LN is correctly normalized, with (cid:104) δ b (cid:105) = 0 and (cid:104) δ (cid:105) = σ , as required. These normalization conditionsmean that as a function of ln(1 + δ b ), the log-distributionis symmetric about ln(1 + δ b ) = − Σ / σ (cid:28) δ b (cid:28) 1, the log-normalPDF in Eq. (43) reduces to the Gaussian PDF to O ( δ b )and O ( σ ); in the linear regime, with σ (cid:28) δ b having an extremely low probability of approaching one,the fluctuations drawn from both the Gaussian and log-normal PDFs are virtually identical. The red band inFig. 1 illustrates the middle-68% containment of the in-homogeneous photon plasma mass assuming a log-normalPDF for the perturbations. Unlike in the case of a Gaus-sian PDF description (illustrated by the blue band), un- physically negative fluctuations are forbidden in this case.For our fiducial PDF, we limit the range of the PDF to10 − ≤ δ b ≤ , removing the highly uncertain PDFtails. B. Analytic PDF Computing the PDF of matter fluctuations from firstprinciples has been effectively studied in the language ofpath integrals, giving expressions that have been shownto be reliable in the nonlinear regime, even at large over-densities [52–55]. Here, we provide only a brief outlineof the derivation of such an analytic PDF, and refer thereader to Ref. [55] for the details of the calculation.Consider a spherical volume of radius r ∗ at some red-shift z containing some density fluctuation δ ∗ obtainedby integrating the spherical volume over a top-hat func-tion. This fluctuation was formed from some field con-figuration δ i ( (cid:126)x ) deep in the linear regime undergoinggravitational collapse, where δ i ( (cid:126)x ) can be described asa Gaussian random field. If the evolution of fluctuationsis purely linear, then the size of linear fluctuations at We will only consider an averaging procedure using a top-hatwindowing function, although more general arguments can bemade for any arbitrary windowing function [55]. z is δ L = (1 + z i ) δ i / (1 + z ), since lin-ear fluctuations grow in proportion to the scale factor ofthe Universe during matter domination. The statisticalproperties of a Gaussian random field are governed en-tirely by the two-point correlation function ξ ( (cid:126)x − (cid:126)y ) ≡(cid:104) δ L ( (cid:126)x ) δ L ( (cid:126)y ) (cid:105) , which is related by the Fourier transformto the linear matter power spectrum P mm,L ( k ). If themapping between overdensities δ ∗ in a cell of size r ∗ and field configurations in the linear regime δ L is well-understood, then the PDF of finding such an overden-sity can be mapped onto the statistical properties of theGaussian random field.Concretely, let us define the functional δ W [ δ L ] whichtakes a given Gaussian field configuration δ L expectedby linear evolution to redshift z of an initial (Gaussian)field configuration δ i , and maps it to the actual densitycontrast δ ∗ averaged over some spherical volume of ra-dius r ∗ , produced by the actual gravitational evolutionof δ i . Then the PDF of δ ∗ is given by a path integralover all Gaussian field configurations δ L with a Gaussianweight [52]: P ( δ ∗ ) = N − (cid:90) D δ L e − S G [ δ L ] δ D ( δ ∗ − δ W [ δ L ]) , (47)where S G [ δ L ] ≡ (cid:90) d (cid:126)x (cid:90) d (cid:126)y δ L ( (cid:126)x ) ξ − ( (cid:126)x − (cid:126)y ) δ L ( (cid:126)y ) , (48)with ξ − defined as the functional inverse of ξ , (cid:90) d (cid:126)z ξ − ( (cid:126)x − (cid:126)z ) ξ ( (cid:126)z − (cid:126)y ) = δ (3)D ( (cid:126)x − (cid:126)y ) . (49)The overall normalization factor is simply N = (cid:90) D δ L e − S G [ δ L ] . (50)Taking the Fourier transform of the integrand in Eq. (48)gives [52] S G [ δ L ] = 12 (cid:90) d (cid:126)k (2 π ) | ˜ δ L ( (cid:126)k ) | P mm,L ( k, z ) , (51)where ˜ δ L ( (cid:126)k ) is the Fourier transform of the field configu-ration δ L .Ref. [55] showed that Eq. (47) can be integrated us-ing the saddle point approximation, by showing that thesaddle point configuration is spherically symmetric, andby making use of the fact that the spherical collapsemodel provides a mapping F between δ ∗ and δ L ( R ∗ ),where δ L ( R ∗ ) is the mean density of the configuration δ L smoothed over a radius R ∗ ≡ r ∗ (1 + δ ∗ ) / , with F ( δ ∗ ) ≡ δ L ( R ∗ ) . (52) Translational and rotational invariance means that ξ ultimatelyonly depends on the magnitude | (cid:126)x − (cid:126)y | . With this, they were able to show that taking into ac-count only spherically-symmetric fluctuations, the prob-ability distribution function is P ( δ ∗ ; z ) = ˆ C ( δ ∗ ) (cid:113) πσ R ∗ ( z ) exp (cid:18) − F ( δ ∗ )2 σ R ∗ ( z ) (cid:19) , (53)where σ R ∗ is the variance of linear matter fluctuationssmoothed with a top-hat of radius R ∗ , σ R ∗ ( z ) = (cid:90) d (cid:126)k (2 π ) P mm,L ( k, z ) | W th ( kR ∗ ) | , (54)with W th being the Fourier transform of the top-hat, W th ( x ) ≡ j ( x ) /x .The intuition behind this result is clear: a density fluc-tuation δ ∗ within a sphere of radius r ∗ at redshift z isformed through spherical collapse of some initial linear density fluctuation, which under linear evolution corre-sponds to a linear density fluctuation of size F ( δ ∗ ) in asphere of radius R ∗ at the same redshift z . Since thelinear density fluctuations follow a Gaussian distributionwith variance σ R ∗ ( z ), P ( δ ∗ ; z ) is also Gaussian with re-spect to F ( δ ∗ ).Several further comments are in order before we areready to use this PDF in our analysis:1. Although Ref. [55] introduces an O (1) asphericalfactor that includes the effects of aspherical fluctu-ations, this factor was not computed for the smallscales of interest to this work. Since we are mostlyinterested in understanding the systematics associ-ated with the use of different PDFs, for simplicity,we neglect this aspherical factor throughout. Inprinciple this prefactor can be computed from the-ory, allowing an improvement to the PDF. Never-theless, this will be a small correction compared tothe baryonic bias with respect to the matter fluc-tuations, which is not included in the analytic cal-culation at the moment. We neglect all other bary-onic effects that may cause a difference between P bb,L ( k, z ) and P mm,L ( k, z ), and take δ ∗ = δ b .2. The PDF as defined in Eq. (53) for δ ∗ is definedwith respect to a sphere of size r ∗ . This is crit-ical in light of the UV divergence exhibited by P mm,L ( k, z ), as discussed in Sec. IV, which leadsto a divergence in σ R ∗ as R ∗ → 0. As we arguedin Sec. IV, baryons naturally have a cut-off lengthscale given by the Jeans length R J , below which thepower spectrum is suppressed. We therefore set thesmoothing scale R ∗ = R J to approximately repro-duce this suppression of power, and take the resultto be the PDF for baryon density fluctuations.In summary, the analytic PDF for baryon fluctuationsthat we adopt in this paper is P an ( δ b ; z ) ≡ ˆ C ( δ b ) (cid:113) πσ R J ( z ) exp (cid:20) − F ( δ b )2 σ R J ( z ) (cid:21) . (55)2We show the full expression for the terms ˆ C and F inApp. B. C. Log-normal PDF with bias The log-normal PDF can be generalized to include anadditional parameter b , known as the bias [63]. Thisdistribution is given by P b LN ( δ b ; z ) ≡ b P LN (cid:18) δ b b ; z (cid:19) , (56)where the choice of b = 1 gives us the log-normal PDFdiscussed in Sec. V A. For this distribution, however, wechoose Σ = ln[1 + σ ( z )] where σ ( z ) = (cid:90) d (cid:126)k (2 π ) P mm ( k, z ) (57)is the variance of the matter power spectrum. The biasparameter is a constant factor relating matter densityfluctuations δ m to baryonic density fluctuations δ b , i.e. , δ b = bδ m . With this in mind, the normalization condi-tions are now (cid:90) ∞− b d δ b P b LN ( δ b ; z ) = 1 , (58) (cid:90) ∞− b d δ b δ b P b LN ( δ b ; z ) = 0 , (59) (cid:90) ∞− b d δ δ P b LN ( δ b ; z ) = b σ . (60)These normalization conditions follow naturally fromhaving matter fluctuations − ≤ δ m < ∞ , and the factthat δ b = bδ m implies σ = b σ . For b > δ b can havedownward fluctuations of up to − b , which are clearly un-physical; however, P b LN has been shown to be a reason-able fit to data [38, 46], and we are once again using thePDF only as a way of capturing systematic uncertainties.In particular, P b LN relies on the distribution of matter and not baryons, allowing us to arrive at a log-normal-like PDF without relying on N -body simulations withbaryonic feedback included, using instead P mm from N -body simulations with cold dark matter only. We againuse the Jeans scale as a UV cut-off for P mm to regulatethe power spectrum. We will adopt the value of b = 1 . D. PDF from voids In Refs. [61, 62], a ΛCDM N -body simulation was per-formed in a box of volume V sim = 500 h − Mpc overthe redshift range 0 ≤ z ≤ 12. The number of voids N voids ( z ), the PDF f voids ( V ; z ) of the volume V of voids,and the PDF g voids ( ρ/ρ ; z ) of the ratio of the mean mat-ter density in voids to the mean cosmological matter den-sity ρ/ρ are all reported. We can now construct a PDF for baryonic fluctuations by making the following sim-plifying assumptions: (i) all underdensities are found invoids that are successfully detected by the simulation; (ii) the density in the void is constant, and is given bythe mean matter density in the void, and iii) no conver-sions happen outside of voids. First, we can work outthe fractional volume of the simulation that is in a void,given by φ voids ( z ) = N voids ( z ) V sim (cid:90) d V V f voids ( V ; z ) . (61) φ voids ∼ . P voids ( δ b ; z ) ≡ φ voids ( z ) g voids (1 + δ b ; z ) . (62)The normalization of P voids is φ voids < 1; in obtainingthe ensemble average in Eqs. (19) and (25), this is equiv-alent to discarding all worldlines at redshift z that are notin voids. This PDF therefore is, by construction, aimedat modeling only underdensities. The assumptions madehere can certainly be improved: not all underdensitiesare found in voids, which necessarily must have a localminimum in density in 3D space, and the void densityprofile should also be taken into account. However, themain purpose of constructing this PDF is less about get-ting an accurate model for the density fluctuations andmore to provide a sanity check on our modeling of un-derdensities using the log-normal or analytic PDFs. VI. VARIANCE OF FLUCTUATIONS A key input to calculating the photon-to-dark photonoscillation probability in the presence of inhomogeneitiesis a description of the spectrum of fluctuations of thephoton plasma. A particular challenge at late times isposed by nonlinear effects, which can be quantified usinginput from N -body simulations. At early times post-recombination on the other hand, spatial fluctuations inthe fraction of free electrons come into play and have tobe accounted for. We now describe in turn the calculationof the variance of fluctuations and relevant inputs in eachregime. A. Free electron fraction perturbations Eq. (1) shows that there are two sources of fluctuationsfor m γ ( t ): fluctuations in the baryon density, as well asfluctuations in the free electron fraction, which we defineas x e ≡ n e /n H , where n H is the number density of bothneutral and ionized hydrogen atoms. So far, we haveneglected fluctuations in x e ; we will now show how fluc-tuations in m γ are related to fluctuations in both baryondensity and x e , and discuss the conditions under which x e can be neglected.3 − − − − k [ h Mpc − ] − − − − − − P i ( k ) [ M p c h − ] P ee = P χχ + P bb + 2 P χ b Electron/baryon power spectra, z = 200 P bb P χχ P ee − P χ b FIG. 5. Baryon (red), electron (green), and free electron frac-tion (blue) power spectra; and negative of the baryon-freeelectron fraction cross-power spectrum (purple) at z = 200.The electron fluctuations are reduced compared to the baryonones due to the baryon and free electron fraction densities be-ing anti-correlated. Consider a point t along a worldline of a photon withsome HI density n HI ( t ) and free electron density n e ( t ),each with a fluctuation from the mean values n HI and n e given by δ HI and δ e respectively, so that n HI = (1 + δ HI ) n HI , n e = (1 + δ e ) n e . (63)We can further rewrite δ e in terms of baryon density fluc-tuations δ b and free electron density fluctuations δ χ ≡ x e x e − . (64)Writing n e (1 + δ e ) = x e (1 + δ χ ) n H (1 + δ b ), δ e = δ b + δ χ + δ χ δ b . (65)We can see that as long as δ χ (cid:28) δ b and δ χ (cid:28) 1, we have δ e = δ b to leading order, i.e. , perturbations in the freeelectron density are given entirely by fluctuations in thebaryon density when free electron fraction perturbationsare small, even in the nonlinear regime. On the otherhand, if δ χ ∼ δ b (cid:28) 1, then δ e = δ χ + δ b . (66)With this new notation, we can rewrite the plasmamass fluctuation δ m γ as δ m γ m γ ≡ m γ − m γ = Aδ e n e − Bω δ HI n HI , (67)where we have defined for convenience the constants A ≡ . × − eV cm , B ≡ . × − cm . (68) In the linear regime, with δ e and δ HI being small andGaussian, m γ is also Gaussian: f ( m γ ; z ) = 1 (cid:113) πσ m γ exp (cid:34) − (1 − m γ /m γ ) σ m γ (cid:35) , (69)where σ m γ ≡ (cid:104) δ m γ δ m γ (cid:105) . (70)We can now make use of Eq. (67) to obtain an expressionfor this variance. For simplicity, we consider the redshiftrange 20 (cid:46) z (cid:46) n HI = (1 − x e ) n H . We find m γ σ m γ = ( A + Bω ) n (cid:104) δ e δ e (cid:105) + B ω n (cid:104) δ b δ b (cid:105)− A + Bω ) Bω n e n H (cid:104) δ e δ b (cid:105) , (71)where (cid:104) δ i δ j (cid:105) = (cid:90) d (cid:126)k (2 π ) P ij, L ( k ) , (72)where P ij, L is the linear (auto) power spectrum of i for i = j , and the cross power spectrum for i and j for i (cid:54) = j ,with i, j = b or e.A more mathematically transparent form of Eq. (71)is obtained by rewriting δ e in terms of δ χ and δ b , whichin the linear regime is simply given by Eq. (66). Thisimmediately leads to the following relation between auto-and cross-power spectra: P ee = P χχ + P bb + 2 P χ b , (73) P eb = P χ b + P bb . (74)Putting together these results, we find m γ σ m γ = m γ (cid:104) δ b δ b (cid:105) + (cid:0) A + Bω (cid:1) n (cid:104) δ χ δ χ (cid:105) + 2 n e m γ (cid:0) A + Bω (cid:1) (cid:104) δ χ δ b (cid:105) . (75)The power spectra that enter into Eq. (75) are all cal-culable in the linear regime after photons decouple frombaryons at z ∼ x e leads to the previous result, σ m γ = σ .The coefficients for the terms on the right-hand side ofEq. (75), however, are of comparable size, and hence the Outside of this range, one must take into account that x e canexceed one, which would require a simple modification to theresults shown here; we omit these modifications since fluctuationsin x e are not important outside the specified redshift range. δ χ → δ χ (cid:28) δ b . To get a sense of how important theseterms are, we plot the power spectra required to computethe two-point correlations shown in Eq. (75) in Fig. 5 at z = 200. Since the baryon δ b and free electron δ χ fluctu-ations are anti-correlated, the presence of free-electronfluctuations causes a reduced variance in electron fluctu-ations (cid:104) δ e δ e (cid:105) at higher redshifts. We see that at z ∼ P χχ < | P χ b | < P bb , with the spectra becomingmore comparable in magnitude for z > z < to extract the transfer functions associated with pertur-bations in the free electron fraction.With this, we can now discuss the importance of δ χ onour results at the following redshifts:1. z (cid:38) . The Universe is completely ionized priorto recombination, and there are no significant per-turbations in x e . We may neglect δ χ ;2. (cid:46) z (cid:46) . At this time, δ χ ∼ δ b , both per-turbations are small, and aside from differences inthe functional form of d (cid:104) P γ → A (cid:48) (cid:105) / d z , this redshiftrange is well approximated by the homogeneouslimit;3. (cid:46) z (cid:46) . During this period, δ χ (cid:28) δ b , andwe may once again neglect δ χ to a good approxi-mation;4. (cid:46) z (cid:46) . This is the period of reionization, anincreasingly nonlinear regime where the behaviorof δ χ depends on the details of reionization, andcan have potentially large effects on the PDF ofplasma mass fluctuations. In principle, δ χ canbe calculated from reionization codes like 21cm-FAST [64, 65], but to avoid this complication,we neglect any γ ↔ A (cid:48) transitions in this epochthroughout our work; and5. z (cid:46) . Reionization is complete, and once againthere are no significant perturbations in x e . Wemay once again neglect δ χ , even though baryondensity fluctuations are highly nonlinear.In summary, we avoid the redshift regime during whichreliably predicting the effect of x e perturbations is non-trivial, staying in regimes where the effect is either ab-sent, or has a minimal and calculable effect on the to-tal conversion probability. This latter regime, 200 (cid:46) z (cid:46) The anticorrelation is due to the fact that recombination is moreefficient when there are more hydrogen atoms present [40]. Available at https://github.com/smsharma/class_public . z − − − k [ h M p c − ] CLASS P bb , L ( k ) 2D interpolation ExtrapolationExpanded envelope2D interpolation × envelope Extrapolation × envelope S i m u l a t i o nd a t a Extrapolations Nonlinear baryon power spectra FIG. 6. Illustration of the scheme used to construct an enve-lope of the nonlinear baryon power spectra at low redshifts,0 (cid:46) z (cid:46) 6, in different redshift z and scale k regimes. Weuse as input the CLASS linear baryon power spectrum P bb,L as well as the envelope of simulation data from Refs. [41, 66],and linearly extrapolate the bias P bb /P mm into regions with-out data (red arrows). For k ≤ . h Mpc − , we use theCLASS linear baryon power spectrum (red). In the range0 . h Mpc − < k < h Mpc − and 0 ≤ z ≤ 3, a 2D inter-polation over available data is performed (blue). We thenextrapolate into the region 3 < z ≤ 6, multiplying the result-ing envelope by a factor of 3 (green). For k > h Mpc − , weextrapolate the power spectra using the CLASS linear baryonpower spectrum as a guide. We then perform a 2D interpola-tion in the range 0 ≤ z ≤ 3, taking as an envelope a factor of3 above and below the central value of the interpolated bias(purple), and then extrapolate this into 3 < z ≤ B. Low-redshift power spectra As described in the last section, at late times z (cid:46) matter fluctuations havebeen extensively studied in the literature, the distinctionbetween baryonic and total matter fluctuations must betaken into account as the two components (baryons anddark matter) evolve separately and baryonic effects be-come increasingly important at late times, especially atthe smaller scales of interest here. In this subsection,we describe our approach for constructing the nonlinearbaryonic power spectra at low redshifts z < − − − − − − − − P i ( k ) [ M p c h − ] z = 0 Linear matterLinear baryonNonlinear baryon (simulation-informed) − − − − − − − − z = 1 − − − k [ h Mpc − ] − − − − − P i ( k ) [ M p c h − ] z = 3 − − − k [ h Mpc − ] − − − − − z = 50 Simulation-informed baryon power spectra FIG. 7. Simulation-informed baryon power spectra at low redshifts, bracketed with the green band and obtained using themethod outlined in Sec. VI B, shown at redshifts z = 0 , 1, and 3. Solid green lines correspond to baryon power spectra fromindividual hydrodynamic simulations as obtained in Ref. [66]. Also shown for comparison are the linear matter and baryonpower spectra as the solid red and blue lines, respectively, also at z = 50. Suppression due to the baryonic Jeans scale canclearly be seen. IllustrisTNG [67], Illustris [68], EAGLE [69], and BA-HAMAS [70] up to k ∼ h Mpc − at the discrete red-shifts z = 0 , , 2, and 3, with Ref. [41] further providingbaryonic spectra from the BAHAMAS simulation at red-shift z = 0 up to k = 500 h Mpc − . We use the follow-ing algorithmic procedure for constructing the nonlinearbaryonic power spectra from these. We first constructlower and upper envelopes encoding the uncertainty onthe power spectra extracted from simulations. Wherefewer than three simulations are available, we obtain themedian spectra over the available simulations and mul-tiply and divide these by a factor of 3 to obtain upperand lower uncertainty envelopes, respectively, motivatedby the magnitude of the typical spread in the regimewhere the full suite of simulations is available. Wherethree or more simulations are available, we use the ex-tremal values over those simulations to construct the en-velopes. At large scales (cid:46) . h Mpc − where simula-tions are not available, we use the well-constrained linearpower spectrum from CLASS. At smaller scales and red-shifts 0 < z < z > 3, we linearly extrapolate thenonlinear baryonic bias, multiplying and dividing the re-sulting power spectra by a factor of 3 to obtain the uncer-tainty envelope. In the regime above z > 20, we simplyuse the linear baryonic power spectrum from CLASS.An illustration of this algorithmic procedure is pro-vided in Fig. 6, showing how the nonlinear baryon powerspectra are estimated at different redshifts z and scales k . The resulting baryon power spectra at several differ-ent redshifts obtained using this procedure are shown inFig. 7 (green envelopes), with the power spectra fromindividual simulations shown as green lines for reference.The inferred variance of fluctuations as a function ofredshift is shown in Fig. 8. At late times z < 6, the vari-ance is informed by the nonlinear baryon power spectrumextracted from hydrodynamic simulations and is shownbracketed by the green band. The variance from the lin-ear baryon power spectrum in this regime is shown as theblue line for comparison. Pre-reionization, the varianceof photon plasma mass fluctuations is given by Eq. (75)and involves the (linear) baryon and free electron pertur-bations, shown as the red line. The variance due to justbaryon perturbations, ignoring the effects of free electronperturbations, is shown as the dashed blue line for com-parison.6 − − z − − − − σ i Variance of fluctuations Photon plasmaLinear baryon R e i o n i z a t i o n ( e x c i s e d ) FIG. 8. Variance of fluctuations as a function of redshift forthe various power spectra configurations considered in thiswork. The photon plasma mass variance is informed by thenonlinear baryon power spectrum from simulations at latetimes z < z > 20 it is informed by the linear baryon andfree electron fraction perturbation spectra. The variance oflinear baryon fluctuations is shown as the dashed blue line,for comparison. VII. SIMULATION STUDIES We use Gaussian and log-normal simulations, whichare relatively cheap to generate, to validate key aspectsof the analytic approach presented in this paper. In par-ticular, we verify that:1. The width of the oscillation probability is describedby Eq. (20), even when the fluctuations in theplasma are non-Gaussian, and2. Averaged over a large number of photon paths, thedifferential transition probability depends only onthe one-point PDF of the underlying plasma den-sity field, and not higher-order moments such astwo-point correlations.We note that the simulations we generate in this sec-tion are fundamentally different from the N -body simula-tions used to inform the baryon power spectra in the pre-vious section—these simulations simply produce a Gaus-sian or log-normal random field with statistics consistentwith a given input power spectrum.We create realizations of the perturbed plasma mass bycreating instances of baryon density fluctuations 1 + δ b ,described either as a Gaussian or log-normal field, andthen obtaining the perturbed plasma mass as m γ = m γ (1 + δ b ), where m γ is the homogeneous plasma mass.Gaussian random fields consistent with the baryon powerspectrum described in the last section are generated us-ing nbodykit [71], and log-normal fields as described in Sec. V A are generated by rescaling these asln(1 + δ LNb ) = − Σ δ b σ × Σ (76)where δ b are the Gaussian overdensities and δ LNb the cor-responding log-normal overdensities. This transforma-tion ensures that the resulting log-normal field has thesame mean and variance as the initial Gaussian field, asin Eqs. (58) and (60).We choose a benchmark dark photon mass of m A (cid:48) =10 − eV, which would correspond to a broad resonancearound z ∼ < z < 6, going up toscales of k max = 20 h Mpc − and up to n points = 100points in each of the simulated boxes. Several boxes arecreated within the specified redshift range for computa-tional efficiency and also to capture the redshift depen-dence of the power spectrum of fluctuations. While thisdoes not capture the full spectrum of fluctuations rele-vant to oscillations (since k max < k J ), the realized fieldshave large enough fluctuations ( δ < − 1) so as to notbe physically describable as Gaussian. We additionallyimpose a top-hat filter of 4 times the grid size in orderto mitigate against the effects of finite gridding at thesmallest simulated scales.An example 2D section through a Gaussian randomfield box generated with this procedure is shown in themiddle panel of Fig. 9, with the corresponding sectionthrough a log-normally-transformed field in the rightpanel. Blue and red patches correspond to positive andnegative (unphysical) values of the resulting field. Theleft panel shows the PDF of fluctuations in both boxes.The Gaussian random field description leads to frequentunphysical, negative fluctuations in this case.The perturbed squared plasma mass over the consid-ered redshift range for one particular sequence of boxesis shown in the left panel of Fig. 10, for the Gaus-sian (blue) and log-normal (red) descriptions. The ho-mogeneous plasma mass is shown as the dashed blackline. Again, frequent unphysically negative values of thesquared plasma mass can be seen in the Gaussian de-scription. We obtain the averaged conversion probabilityby creating a large number of such simulations, draw-ing photon paths separated by at least twice the sizeof the top-hat filter (to ensure they are sufficiently un-correlated) through them, and numerically calculatingtransition probabilities at each crossing using Eq. (15).Probabilities over a large number of photon paths arethen histogrammed to obtain the numerical estimates ford (cid:104) P γ → A (cid:48) (cid:105) / d z , shown in the right panel of Fig. 10 as thedashed red line for the log-normal case. The analytically-computed differential conversion probability for this con-figuration is shown in solid red, and provides a goodmatch to the numerical results. The analytic Gaussiandescription, shown in blue, does not accurately describedthe conversion probability in this regime.At higher redshifts and in the linear regime, on the7 − δ z = 5 Simulation PDF histogram Log-normalGaussian D c [Mpc] D c [ M p c ] Gaussian simulation D c [Mpc] Log-normal simulation − . − . − . − . . . . . . + δ FIG. 9. 2D section through a Gaussian random field simulated at z ∼ other hand, the Gaussian PDF is an excellent descriptionof the plasma mass fluctuations. Fig. 11 shows a compar-ison of the analytically-computed differential conversionand the probability derived by considering photon pathsthrough Gaussian random field-simulations of the plasmamass, showing once again good agreement between thetwo. VIII. SYSTEMATICS OF CONVERSIONPROBABILITY AND ENERGY INJECTION Given a PDF of density fluctuations and a descriptionof the fluctuations through the power spectra, the differ-ential conversion probability d (cid:104) P γ → A (cid:48) (cid:105) / d ln z at a givenredshift, for a given dark photon mass, can be computed.This is the main deliverable of this paper, and is plottedin the top rows of Fig. 12 and Fig. 13 for various PDFdescriptions and benchmark masses m A (cid:48) = 4 × − eV(red), 10 − eV (blue), and 10 − eV (green). The cu-mulative probabilities above (below) a given redshift areplotted in the middle(bottom) panels of these figures.Fig. 12 shows various log-normal PDFs—including alloverdensities and underdensities (dashed lines), impos-ing 10 − (cid:46) δ b (cid:46) (solid lines), and additionallywith a bias b = 1 . − (cid:46) δ b (cid:46) (solid lines),and the voids PDF (dotted lines). For ease of compari-son, these are normalized such that the cumulative prob-abilities for the fiducial 10 − (cid:46) δ b (cid:46) -boundedlog-normal PDF cases are unity. The primary focus hereis on dark photons of masses (cid:46) − eV, where the con-version probability is dominated by a broad efficiency ofconversions at late times z (cid:46) 6. The lower uncertaintyenvelope of the simulation-informed power spectrum de-scribed in Sec. VI B was used to inform the variance forthe PDFs in these plots; using the power spectrum cor- responding to the upper uncertainty envelope producesqualitatively similar results.In order to illustrate how the total γ → A (cid:48) conver-sion probability is affected by various PDFs for differ-ent dark photon masses, the total conversion probabil-ity per squared kinetic mixing parameter (cid:15) is shown inthe left panel of Fig. 14 for the different PDFs we haveconsidered. Log-normal (dashed red), log-normal impos-ing 10 − (cid:46) δ b (cid:46) (solid red), log-normal withbias b = 1 . e.g. , inthe case of dark photon dark matter) is shown in theright panel of Fig. 14. In each case, the correspondingquantities under the assumption of a homogeneous pho-ton plasma are shown in dotted gray. It can be seen thatinhomogeneities have a significant effect on the nature ofphoton-to-dark photon oscillations, either underestimat-ing or overestimating the total conversion probability andenergy deposition depending on the dark photon masspoint considered. Variation is also observed across thedifferent PDFs considered; however, after restricting tofluctuations of size 10 − (cid:46) δ b (cid:46) , the log-normaland analytic PDFs show quantitatively similar behavior,with the log-normal PDF being somewhat more conser-vative. For this reason, henceforth in this paper andin Paper I, we use the log-normal PDF with variance in-formed by hydrodynamic simulations as the benchmarkfor computing the effects of γ ↔ A (cid:48) conversions. In theabsence of dedicated PDFs capturing baryonic effects andtheir uncertainties to the smallest relevant scales, we ad-vocate for its use in applications beyond those consid-ered in these papers where the effects of inhomogeneitiesin the nonlinear regime on γ ↔ A (cid:48) conversions may beimportant.Conversions at earlier times z (cid:38) 100 can be well-described by a Gaussian in redshift with a weakly8 . . . . . z m γ [ × − e V ] Gaussian vs. log-normal simulation Log-normal simulationGaussian simulationHomogeneous . . . . . z . . . . . . . d h P γ → A i / d z k max = 20 h Mpc − , r filt = 2 . h − Simulation vs. analytic probability m A = 10 − eV Analytic, Log-normalSimulations, Log-normalAnalytic, Gaussian FIG. 10. (Left) 1D sections through realizations of Gaussian (blue) and log-normal (red) perturbations in the squared plasmamass. The homogeneous plasma mass is shown in dashed gray. (Right) The log-normal differential oscillation probabilityaveraged over a large number of photon paths drawn through simulations (dashed red) and derived analytically (solid red),with good agreement between the two. The analytic Gaussian description in shown in solid blue. 96 98 100 102 104 z . . . . . . . d h P γ → A i / d z k max = 20 h Mpc − r filt = 2 . h − Simulation vs. analytic probability AnalyticGaussian simulation FIG. 11. Differential conversion probability obtained bydrawing photon paths through Gaussian random field sim-ulations (dashed red) and computed analytically (solid red),for a resonance around z = 100. Good agreement betweensimulations and the analytic description can be seen. redshift-dependent variance, described in Eq. (35). Ex-ample differential conversion probabilities are shown inthe left panel of Fig. 15 for resonance redshifts spanning100 ≤ z res ≤ z between redshifts where thesquared plasma mass is ± σ m γ / σ containment interval. The presence of spatial perturbations in the freeelectron fraction becomes increasingly important closerto the redshift of recombination, although the relativewidth of the conversion feature is already less than onepart in 10 − by z res = 600.Due to the sensitive dependence of the conversion prob-ability on small-scale physics as discussed in Sec. IV,it is illustrative to see how the total conversion proba-bility depends on the maximum scale k max considered.This is illustrated in Fig. 16 for our benchmark masses,shown as the ratio of the total probability consideringscales up to k max to the asymptotic probability. Wesee that the total probability approaches the asymptoticvalue around the characteristic baryon Jeans scale at latetimes, k J ∼ h Mpc − . Note that neglecting the effectof small scales is not necessarily conservative and maysignificantly underestimate or overestimate the conver-sion probability.Finally, although we advocate restricting to fluctua-tions in the range 10 − (cid:46) δ b (cid:46) where the differ-ent PDF descriptions considered show qualitative agree-ment, it is instructive to ask how expanding this rangeand including larger underdensities and overdensities inthe tails of the PDFs can affect the oscillation physics.In Fig. 17, we show the total conversion probability as afunction of dark photon mass varying the range of fluc-tuations from 10 − (cid:46) δ b (cid:46) 10 to 10 − (cid:46) δ b (cid:46) for the log-normal (solid red lines) and analytic (dashedblue lines) PDFs. Although the two descriptions disagreefor fluctuations beyond 10 − (cid:46) δ b (cid:46) , in eithercase larger conversion probabilities over a much widerrange of dark photon masses can be seen when includingconversions from fluctuations deeper in the tails of thePDFs. This motivates the need for a better understand-ing of the nonlinear baryon PDF at late times. A similar9 − − − d h P γ → A i / d l n z Log-normal PDFs R e i o n i z a t i o n ( e x c i s e d ) − − − h P γ → A i ( > z ) Log-normal PDF+ 10 − < δ < (fiducial)+ b = 1 . R e i o n i z a t i o n ( e x c i s e d ) − z − − h P γ → A i ( < z ) m A = 4 × − eV m A = 10 − eV m A = 10 − eV R e i o n i z a t i o n ( e x c i s e d ) FIG. 12. The differential conversion probability d (cid:104) P γ → A (cid:48) (cid:105) / d ln z (top row), cumulative conversion probability above a givenredshift z (middle row), and cumulative conversion probability below a given redshift z , shown for a log-normal PDF (dashedlines), our fiducial log-normal PDF with 10 − (cid:46) δ b (cid:46) (solid lines), and additionally with a bias b = 1 . m A (cid:48) = 4 × − eV (red), 10 − eV (blue), and 10 − eV (green) are shown. Lines are normalized such that thecumulative probabilities for the 10 − (cid:46) δ b (cid:46) -bounded log-normal PDF cases are unity. conclusion can be drawn for A (cid:48) → γ dark-photon darkmatter conversions, also shown in Fig. 17. IX. CONCLUSIONS In this paper, we have studied photon-dark photonoscillations in the early Universe, deriving a formalismfor computing the averaged probability of conversions inboth directions, taking into account the effect of inho-mogeneities in the photon plasma. We found that theaverage probability of γ ↔ A (cid:48) and the average energy in-jected per baryon for A (cid:48) → γ for dark photon dark matterare completely specified given the standard ΛCDM pa-rameters as well as three inputs: (i) a description of theone-point PDF of baryon fluctuations, (ii) the baryonpower spectrum which, to a good approximation in thelow-redshift Universe, provides the variance of plasmamass fluctuations, and (iii) fluctuations in the free elec-tron fraction, which contributes to the variance of plasmamass fluctuations at high redshift. To understand thesystematic uncertainties associated with the PDF and thevariance of fluctuations, we studied several independentchoices of the one-point PDF. We also constructed a non-linear baryon power spectrum that is informed by high-resolution hydrodynamic N -body simulations, allowing us to characterize the behavior of baryons at small scales.Finally, we also performed a series of Gaussian and log-normal random field simulations in order to validate ouranalytic results, finding agreement between theory andsimulations.In our companion work Paper I, we have applied thisformalism in order to derive constraints on the dark pho-ton kinetic mixing parameter (cid:15) by through the effect of γ → A (cid:48) conversions on the CMB spectrum as measuredby COBE/FIRAS in the general case, as well as dedicatedconstraints for the case of dark photon dark matter ob-tained by computing the amount of IGM heating due to A (cid:48) → γ conversions. We found that previous constraintsassuming a homogeneous plasma were not conservative,and were able to expand the mass range over which reso-nant oscillations are possible due to conversions in plasmaunderdensities and overdensities. We also found goodagreement between constraints obtained using differentPDFs and power spectra, showing that we have a suffi-ciently good understanding of baryon fluctuations to setreliable constraints.The formalism that we have developed across bothpapers has additional applications. For example, per-turbations in the photon plasma mass will also modifyresonant oscillations of photons into axion-like-particles,which can occur in the presence of primordial magnetic0 − − − d h P γ → A i / d l n z Analytic and voids PDFs R e i o n i z a t i o n ( e x c i s e d ) − − − h P γ → A i ( > z ) Analytic PDF+ 10 − < δ < Voids PDF R e i o n i z a t i o n ( e x c i s e d ) − z − − h P γ → A i ( < z ) m A = 4 × − eV m A = 10 − eV m A = 10 − eV R e i o n i z a t i o n ( e x c i s e d ) FIG. 13. The same as Fig. 12, shown for the analytic PDF (dashed lines), additionally imposing 10 − (cid:46) δ b (cid:46) (solidlines), and the voids PDF (dotted lines). − − − − − − m A [eV] h P γ → A i / (cid:15) γ → A total conversion probability HomogeneousLog-normal+ 10 − < δ < (fiducial)+ b = 1 . − − − − − m A [eV] h E A → γ i / (cid:15) × Ω D M / Ω A [ e V ] DM A → γ deposited energy per baryon HomogeneousLog-normal+ 10 − < δ < (fiducial)+ b = 1 . FIG. 14. (Left) The total γ ↔ A (cid:48) conversion probability as a function of dark photon mass, and (Right) The A (cid:48) → γ darkphoton dark matter energy deposited per baryon as a function of dark photon mass, shown for different choices of PDFsexplored in this work: log-normal (red dashed), the fiducial log-normal with 10 − (cid:46) δ b (cid:46) (red solid), log-normal witha bias b = 1 . fields [72]. Moreover, relativistic dark photons can alsoresonantly inject photons, which can be tested by 21-cmobservations [13, 14, 73]. Photon-to-dark photon oscil-lations in an inhomogeneous background will also im-print anisotropies in the CMB that may be testable by Planck [23] or future CMB probes [74], as also explored in Ref. [21].A comparison with our results and methodology withthose presented in related recent studies, auxillary in-formation about the analytic PDF employed in thiswork, and a complementary derivation of the Landau-Zener formula for resonant conversions in the lan-1 − − z − z res . . . . . . . d h P γ → A i / d z [ n o r m a li ze d ] Differential conversion probability z res = 100 z res = 200 z res = 300 z res = 400 z res = 500 z res = 600 200 400 600 800 z res − − − − − ∆ z / z r e s Width of resonance FullWithout x e perturbationsApproximation FIG. 15. (Left) The differential conversion probability d (cid:104) P γ → A (cid:48) (cid:105) / d z for resonant conversion at higher redshifts z res = 100 to600, shown centered on the resonant redshift z res . (Right) Relative width of the resonance as a function of resonance redshift z res . Shown with (red) and without (blue) accounting for perturbations in the electron ionization fraction x e . The dottedgreen line shows the approximate width as given by Eq. (37), showing good agreement with the numerical estimate withoutaccounting for x e perturbations. − − k max [ h Mpc − ] h P γ → A i / h P γ → A i k m a x = ∞ Dependence on k max m A = 4 × − eV m A = 10 − eV m A = 10 − eV FIG. 16. For the fiducial log-normal PDF, the ratio of thetotal conversion probability using fluctuations only up to agiven scale k max and its asymptotic value, shown for masses m A (cid:48) = 4 × − eV (red), 10 − eV (blue), and 10 − eV(green). guage of thermal field theory is provided in the ap-pendices. The code used to obtain the results in bothpapers is available at https://github.com/smsharma/dark-photons-perturbations . ACKNOWLEDGMENTS We thank Yacine Ali-Ha¨ımoud, Masha Baryakhtar,Asher Berlin, Julien Lesgourgues, Sam McDermott, Alessandro Mirizzi, Julian Mu˜noz, Stephen Parke,Maxim Pospelov, Josef Pradler, Javier Redondo, Ro-man Scoccimarro, Anastasia Sokolenko, Alfredo Ur-bano, Edoardo Vitagliano, Sam Witte, and Chih-LiangWu for helpful conversations. We thank Marcel vanDaalen for providing baryonic power spectra from high-resolution BAHAMAS simulations. We are especiallygrateful to Misha Ivanov for many enlightening dis-cussions regarding the analytic PDF of density fluctu-ations utilized in this work. AC acknowledges sup-port from the “Generalitat Valencian” (Spain) throughthe “plan GenT” program (CIDEGENT/2018/019), aswell as national grants FPA2014-57816-P, FPA2017-85985-P, and the European projects H2020-MSCA-ITN-2015//674896-ELUSIVES. HL is supported by the DOEunder contract DESC0007968. SM and JTR are sup-ported by the NSF CAREER grant PHY-1554858 andNSF grant PHY-1915409. SM is additionally supportedby NSF grant PHY-1620727 and the Simons Founda-tion. JTR acknowledges hospitality from the Aspen Cen-ter for Physics, which is supported by the NSF grantPHY-1607611. This work made use of the NYU ITHigh Performance Computing resources, services, andstaff expertise. The authors are pleased to acknowledgethat the work reported on in this paper was substan-tially performed using the Princeton Research Comput-ing resources at Princeton University which is a con-sortium of groups including the Princeton Institute forComputational Science and Engineering and the Prince-ton University Office of Information Technology’s Re-search Computing department. This research has madeuse of NASA’s Astrophysics Data System. We acknowl-edge the use of the Legacy Archive for Microwave Back-ground Data Analysis (LAMBDA), part of the High En-2 − − − − − − − m A [eV] h P γ → A i / (cid:15) γ → A dependence on PDF tails − < δ < − < δ < − < δ < − < δ < HomogeneousLog-normalAnalytic − − − − − − − m A [eV] h E A → γ i / (cid:15) × Ω D M / Ω A [ e V ] DM A → γ dependence on PDF tails − < δ < − < δ < − < δ < − < δ < HomogeneousLog-normalAnalytic FIG. 17. Dependence of (Left) the total γ ↔ A (cid:48) conversion probability and (Right) the energy injected by A (cid:48) → γ dark-photondark matter on the tails of the plasma mass PDF, shown for the log-normal (blue) and analytic (dashed red) PDFs. The totalhomogeneous probability is shown as dotted grey, for comparison. ergy Astrophysics Science Archive Center (HEASARC).HEASARC/LAMBDA is a service of the AstrophysicsScience Division at the NASA Goddard Space FlightCenter. This research made use of the astropy [75, 76],CAMB [77, 78], CLASS [39], HyRec [79], IPython [80],Jupyter [81], matplotlib [82], nbodykit [71], NumPy [83], seaborn [84], pandas [85], SciPy [86], and tqdm [87] soft-ware packages. Appendix A: Comparison with other work In this section we present a comparison of the formal-ism and results described in this work and in Paper I withthose presented in several recent studies which also at-tempt to model inhomogeneous γ ↔ A (cid:48) oscillations andtheir observational consequences.In Refs. [20, 21], the conversion probability as photonspass through inhomogeneities was determined throughthe use of the EAGLE simulation [69] with baryons.Lines were drawn at random for each redshift snapshot inthe simulation, and one hundred continuous lines-of-sightin the range 0 < z < γ → A (cid:48) conversion with the inhomogeneities encountered in thesimulation, and used to set limits on the kinetic mixingparameter (cid:15) . Ref. [21] found good agreement betweentheir results and those presented in Paper I. They alsouse a similar approach to obtain CMB power spectrumconstraints by comparing the fluctuation in conversionprobability between line-of-sights, finding a weaker limitthan that obtained from the COBE/FIRAS energy spec-trum measurement.We note that while we also use input from the sameEAGLE simulation [69], we only rely on the baryon powerspectrum from this and other simulations, rather than the full spatial information. This significantly simplifiesthe process of understanding γ ↔ A (cid:48) oscillations, and al-lows us to do two things: (i) avoid the need to smooth thesimulation excessively, and (ii) capture the uncertaintyassociated with different choices of the one-point PDF.We will now discuss each point in turn:1. Smoothing. N -body simulations have a finite res-olution, and it is often the case that some smooth-ing of the data needs to be done prior to anal-ysis. Finite resolution effects and smoothing ul-timately introduce an effective cut-off k res in thepower spectrum of fluctuations. For values of k res (cid:46) kpc, Fig. 16 shows that the calculated con-version probability can deviate significantly fromthe asymptotic value we infer using the proceduredescribed in Sec. VI B. In Ref. [20], the lines-of-sight are smoothed over a comoving pixel size of20 kpc × 20 kpc × 250 kpc, while in Ref. [21], this isreduced to 20 kpc × 20 kpc × 25 kpc, with the au-thors of Ref. [21] finding no difference in the conver-sion probability between the two smoothing scales.We have checked that performing this anisotropicsmoothing over comoving 20 kpc × 20 kpc × 250 kpcpixels produces a variance of fluctuations σ b thatis similar to having k res ∼ h Mpc − in the red-shift range 0 < z < 6. This should therefore leadto similar results for the conversion probability, asshown in Fig. 16. This also explains why Ref. [21]observes no difference in results between the twopixel sizes. In general, however, smoothing mustbe used with caution due to the ultraviolet diver-gence of the variance of fluctuations, as describedin Secs. IV B and VIII. Too large of a smoothingscale, either due to the finite resolution of a sim-ulation or post-processing of the results, may leadto very different and incorrect (not necessarily con-3servative) outcomes. It is important to use highresolution results and smooth as little as possible.2. Capturing uncertainties. As we showed inSec. VII, the full simulation data is not neces-sary to determine the γ → A (cid:48) conversion probabil-ity in the presence of inhomogeneities; knowledgeof the one-point PDF alone is sufficient for that.Our work therefore represents a significant sim-plification compared to constructing lines-of-sightthrough simulation results. In particular, we do notneed to rely on the outcome of a single simulationto extract our results, as was done in Refs. [20, 21];we have shown how our results change dependingon our choice of one-point PDFs and baryon powerspectra, allowing us to study the uncertainty as-sociated with these inputs based a large array oftheoretical and simulation results. This is particu-larly important for conversions in large under- andoverdensities, where the PDFs are highly uncertain.The authors of Ref. [22] on the other hand reconsideredthe bounds on dark photon dark matter A (cid:48) → γ conver-sions, obtained from Ly- α observations of the IGM tem-perature, in the presence of inhomogeneities. Their over-all approach to the problem is similar to ours, althoughthey do not generalize their results to treat γ → A (cid:48) as wedo in our work, where a CMB photon passes through mul-tiple level crossings along its path at which m γ = m A (cid:48) .Our results, however, differ from Ref. [22] for the follow-ing reasons:1. Value of the Jeans scale. The authors ofRef. [22] adopt a value of the Jeans scale close to R J ∼ T b ∼ K. However, as we discussed inSec. IV B, a suppression at these scales is not seenin any of the N -body simulations (with baryonicphysics included) we used to infer our cut-off scale k J , which is then smaller by roughly two orders ofmagnitude. This is due to the increasingly nonlin-ear behavior of baryons at late times, which makesdifficult to analytically predict the scale at whichstructure formation is suppressed. Their choice ofthe Jeans scale is therefore an underestimate, lead-ing to overly narrow d P/ d z as a function of red-shift. This can have a large effect on the derivedconstraints, as we show in Fig. 16.2. Ly- α observations sensitivity. The authors ofRef. [22] note that IGM temperature measurementsfrom Ly- α observations are not sensitive to largeunder- or overdensities, which is of particular im-portance if the energy injection is deposited locally(see the following point). Too large values of δ b leadto a large optical depth of the IGM medium, lead-ing to near-total absorption of Ly- α photons, pre-venting us from learning anything about optically thick regions; on the other hand, too low δ b wouldmean no absorption lines at all, which is requiredto deduce the IGM temperature [88]. Ref. [22] pro-posed two heuristic ways of correcting for this; theirfiducial method, for example, rescales the energydeposited by a factor proportional to the deriva-tive of the Ly- α absorption probability, while theiralternative method simply assumes that no temper-ature measurements are possible outside of someoptical depth range. Both prescriptions adoptedin Ref. [22] are reasonable, but nevertheless onlyheuristic, and have many caveats. They depend,for example, on the IGM temperature-density rela-tion, assumed to be T ∝ (1 + δ b ) γ − where γ ∼ . Energy injection. We worked under the assump-tion that energy injection is a global phenomenon, i.e. , energy injected from A (cid:48) → γ conversions isshared evenly among all baryons. The authors ofRef. [22], on the other hand, assume local energyinjection, where the energy is deposited only intobaryons at the point where conversions occur. Forcompleteness, we have also derived the energy de-position per baryon under the local assumptions,shown in Eq. (27). This expressions agrees with theexpression derived in Ref. [22], although we showthat it reduces to a much simpler form shown inEq. (28) within our framework, as compared to theresults shown in Ref. [22]. The authors of Ref. [22]justify the local assumption by noting that the elec-trons that absorb this energy are nonrelativistic,and so the energy transport timescale has to bemuch longer than the age of the Universe.We expect the transport of energy from A (cid:48) → γ conversions to lie somewhere in between bothregimes. The argument in Ref. [22] about nonrela-tivistic electrons applied to reionization, for exam-ple, would seem to preclude the possibility of com-plete reionization across the entire Universe. In-stead, as in the process of reionization, we expectphotons with energy above the ionization thresholdof HI to play a large role in energy transport. Dur-ing HeII reionization, the epoch in which we deriveour constraints in Paper I, the IGM is already at T b ∼ K and will be heated beyond that due tothe A (cid:48) conversion. The blackbody spectrum of theIGM contains ionizing photons, which have a longinteraction path length, potentially comparable tothe size of the Universe at redshifts 2 (cid:46) z (cid:46) 6. Thismay allow for energy transport over large distances.Whether or not the energy injection is local is a non-trivial problem which requires a more involved treatmentof the complete transport equations describing the sys-tem under consideration; we defer such an effort to fu-ture work. To account for general uncertainties regardinglarge under- and overdensities, especially with regard to4uncertainties in the tails of the baryon one-point PDFs,we presented our limits on (cid:15) as a function of the ex-pected range in δ b in Paper I. In addition, we showthe A (cid:48) → γ dark photon dark matter constraints de-rived from Ly- α temperature measurements of the IGMduring HeII reionization in Paper I, neglecting densitieswhich lead to an optical depth for Ly- α photons that sat-isfy exp( − τ ) < . 05 or exp( − τ ) > . 95, the ‘alternate’method adopted by Ref. [22]. Our results broadly agreewith those obtained in Ref. [22]. Appendix B: Functions for the analytic PDF Following Ref. [55], the function F ( δ ∗ ) is defined as thecomposition of two functions F ≡ G ◦ F − , (B1)where G ( θ ) ≡ 320 (6[ θ − s ( θ )]) / , (B2)and F ( θ ) ≡ s ( θ ) − θ ] c ( θ ) − − , (B3)with s ( θ ) ≡ (cid:40) sin θ , δ ∗ > , sinh θ , δ ∗ ≤ , c ( θ ) ≡ (cid:40) cos θ , δ ∗ > , cosh θ , δ ∗ ≤ . (B4)ˆ C ( δ ∗ ) is then defined asˆ C ( δ ∗ ) ≡ F (cid:48) ( δ ∗ ) + F ( δ ∗ )1 + δ ∗ (cid:18) − ξ R ∗ σ R ∗ (cid:19) , (B5)and ξ R ∗ ≡ π (cid:90) d k k sin( kR ∗ ) kR ∗ W th ( kR ∗ ) P m,L ( k ) , (B6) where W th is the Fourier transform of the top-hat func-tion defined in Sec. V B, and P m,L is the linear matterpower spectrum. Appendix C: Thermal Field Theory derivation ofLandau-Zener probability Here we give a brief derivation of the Landau-Zenerformula using thermal field theory techniques. Indeed theconversion of CMB photons to dark photons can be seenas the production of dark photons from a thermal bathof photons following a blackbody spectrum. FollowingRefs. [4, 30, 31], we can write the production rate of darkphotons asΓ prod = (cid:18) e ω/T − (cid:19) (cid:15) m A (cid:48) Γ ω Γ + ( m γ − m A (cid:48) ) ≡ f γ ( ω, T ) (cid:15) m A (cid:48) Γ ω Γ + ( m γ − m A (cid:48) ) , (C1)where Γ is the damping rate of the plasmon quanta, m γ isthe plasma mass acquired by the photons in the plasma.The first factor f γ ( ω, T ) is the photon occupation num-ber, with T being the CMB temperature. The secondfactor is the probability of conversion per unit time. Inthe limit of the narrow width approximation, assumingthat the plasmons are weakly damped, the probability ofconversion reduces toΓ prod f γ ( ω, T ) → (cid:15) m A (cid:48) ω δ D (cid:16) m γ − m A (cid:48) ω (cid:17) , (C2)where we used the definition of the Dirac δ D -functionlim α → αα + x = δ D ( x ) . (C3)We can then integrate it over time along the photonpath to find P γ → A (cid:48) = (cid:90) d t (cid:15) m A (cid:48) ω δ D (cid:16) m γ − m A (cid:48) ω (cid:17) = (cid:88) i (cid:15) m A (cid:48) ω ( t i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ln m γ d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − t = t i , (C4)which is indeed in agreement with Eq. (15). [1] B. Holdom, Phys. Lett. , 196 (1986).[2] J. Redondo and M. Postma, JCAP , 005 (2009),arXiv:0811.0326 [hep-ph].[3] A. E. Nelson and J. Scholtz, Phys. Rev. D84 , 103501(2011), arXiv:1105.2812 [hep-ph].[4] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Re- dondo, and A. Ringwald, JCAP , 013 (2012),arXiv:1201.5902 [hep-ph].[5] A. Fradette, M. Pospelov, J. Pradler, and A. Ritz, Phys.Rev. D90 , 035022 (2014), arXiv:1407.0993 [hep-ph].[6] H. An, M. Pospelov, J. Pradler, and A. Ritz, Phys. Lett. B747 , 331 (2015), arXiv:1412.8378 [hep-ph]. [7] P. W. Graham, J. Mardon, and S. Rajendran, Phys.Rev. D93 , 103520 (2016), arXiv:1504.02102 [hep-ph].[8] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi,and F. Takahashi, Phys. Lett. B801 , 135136 (2020),arXiv:1810.07188 [hep-ph].[9] J. A. Dror, K. Harigaya, and V. Narayan, Phys. Rev. D99 , 035036 (2019), arXiv:1810.07195 [hep-ph].[10] R. T. Co, A. Pierce, Z. Zhang, and Y. Zhao, Phys. Rev. D99 , 075002 (2019), arXiv:1810.07196 [hep-ph].[11] M. Bastero-Gil, J. Santiago, L. Ubaldi, and R. Vega-Morales, JCAP , 015 (2019), arXiv:1810.07208 [hep-ph].[12] A. J. Long and L.-T. Wang, Phys. Rev. D99 , 063529(2019), arXiv:1901.03312 [hep-ph].[13] M. Pospelov, J. Pradler, J. T. Ruderman, andA. Urbano, Phys. Rev. Lett. , 031103 (2018),arXiv:1803.07048 [hep-ph].[14] K. Choi, H. Seong, and S. Yun, (2019), arXiv:1911.00532[hep-ph].[15] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather,R. A. Shafer, and E. L. Wright, Astrophys. J. , 576(1996), arXiv:astro-ph/9605054 [astro-ph].[16] A. Mirizzi, J. Redondo, and G. Sigl, JCAP , 026(2009), arXiv:0901.0014 [hep-ph].[17] K. E. Kunze and M. ´A. V´azquez-Mozo, JCAP , 028(2015), arXiv:1507.02614 [astro-ph.CO].[18] S. D. McDermott and S. J. Witte, (2019),arXiv:1911.05086 [hep-ph].[19] A. Caputo, H. Liu, S. Mishra-Sharma, and J. T. Ruder-man, (2020), arXiv:2002.05165 [astro-ph.CO].[20] K. Bondarenko, J. Pradler, and A. Sokolenko, (2020),arXiv:2002.08942 [astro-ph.CO].[21] A. A. Garcia, K. Bondarenko, S. Ploeckinger, J. Pradler,and A. Sokolenko, (2020), arXiv:2003.10465 [astro-ph.CO].[22] S. J. Witte, S. Rosauro-Alcaraz, S. D. McDermott, andV. Poulin, (2020), arXiv:2003.13698 [astro-ph.CO].[23] N. Aghanim et al. (Planck), (2019), arXiv:1907.12875[astro-ph.CO].[24] B. Dasgupta and A. Dighe, Phys. Rev. D75 , 093002(2007), arXiv:hep-ph/0510219 [hep-ph].[25] L. Pauling and E. Wilson, Introduction to Quantum Me-chanics: With Applications to Chemistry , Dover Bookson Physics (Dover Publications, 1985).[26] L. Wolfenstein, Phys. Rev. D17 , 2369 (1978),[,294(1977)].[27] S. P. Mikheyev and A. Yu. Smirnov, Sov. J. Nucl. Phys. , 913 (1985), [Yad. Fiz.42,1441(1985), 305(1986)].[28] A. Hook, Y. Kahn, B. R. Safdi, and Z. Sun, Phys. Rev.Lett. , 241102 (2018), arXiv:1804.03145 [hep-ph].[29] R. A. Battye, B. Garbrecht, J. I. McDonald, F. Pace, andS. Srinivasan, (2019), arXiv:1910.11907 [astro-ph.CO].[30] E. Hardy and R. Lasenby, JHEP , 033 (2017),arXiv:1611.05852 [hep-ph].[31] J. Redondo and G. Raffelt, JCAP , 034 (2013),arXiv:1305.2920 [hep-ph].[32] G. Raffelt and L. Stodolsky, Phys. Rev. D37 , 1237(1988).[33] G. L. Fogli, E. Lisi, A. Mirizzi, and D. Montanino, JCAP , 012 (2006), arXiv:hep-ph/0603033 [hep-ph].[34] A. Friedland and A. Gruzinov, (2006), arXiv:astro-ph/0607244 [astro-ph].[35] S. Petcov, , 427 (1987). [36] S. O. Rice, Bell System Technical Journal , 282 (1944).[37] G. Lindgren, Stationary stochastic processes : theory andapplications (CRC Press, Taylor & Francis Group, BocaRaton, FL, 2013).[38] L. Hurtado-Gil, V. J. Mart´ınez, P. Arnalte-Mur, M. J.Pons-Border´ıa, C. Pareja-Flores, and S. Paredes, Astron.Astrophys. , A40 (2017), arXiv:1703.01087 [astro-ph.CO].[39] D. Blas, J. Lesgourgues, and T. Tram, JCAP , 034(2011), arXiv:1104.2933 [astro-ph.CO].[40] A. Lewis, Phys. Rev. D76 , 063001 (2007),arXiv:0707.2727 [astro-ph].[41] M. P. van Daalen, I. G. McCarthy, and J. Schaye,Mon. Not. Roy. Astron. Soc. , 2424 (2020),arXiv:1906.00968 [astro-ph.CO].[42] E. Hubble, Astrophys. J. , 8 (1934).[43] P. Coles and B. Jones, Mon. Not. Roy. Astron. Soc. ,1 (1991).[44] L. Clerkin et al. (DES), Mon. Not. Roy. Astron. Soc. ,1444 (2017), arXiv:1605.02036 [astro-ph.CO].[45] D. Gruen et al. (DES), Phys. Rev. D98 , 023507 (2018),arXiv:1710.05045 [astro-ph.CO].[46] V. Wild et al. (2dFGRS), Mon. Not. Roy. Astron. Soc. , 247 (2005), arXiv:astro-ph/0404275 [astro-ph].[47] L. Kofman, E. Bertschinger, J. M. Gelb, A. Nusser,and A. Dekel, Astrophys. J. , 44 (1994), arXiv:astro-ph/9311028 [astro-ph].[48] I. Kayo, A. Taruya, and Y. Suto, Astrophys. J. , 22(2001), arXiv:astro-ph/0105218 [astro-ph].[49] A. Klypin, F. Prada, J. Betancort-Rijo, and F. D. Al-bareti, Mon. Not. Roy. Astron. Soc. , 4588 (2018),arXiv:1706.01909 [astro-ph.CO].[50] F. Bernardeau, Astrophys. J. , 1 (1992).[51] F. Bernardeau, S. Colombi, E. Gaztanaga, andR. Scoccimarro, Phys. Rept. , 1 (2002), arXiv:astro-ph/0112551 [astro-ph].[52] P. Valageas, Astron. Astrophys. , 412 (2002),arXiv:astro-ph/0107126 [astro-ph].[53] P. Valageas, Astron. Astrophys. , 477 (2002),arXiv:astro-ph/0109408 [astro-ph].[54] S. Matarrese, L. Verde, and R. Jimenez, Astrophys. J. , 10 (2000), arXiv:astro-ph/0001366 [astro-ph].[55] M. M. Ivanov, A. A. Kaurov, and S. Sibiryakov, JCAP , 009 (2019), arXiv:1811.07913 [astro-ph.CO].[56] Ya. B. Zeldovich, J. Einasto, and S. F. Shandarin, Na-ture , 407 (1982).[57] M. Plionis and S. Basilakos, Mon. Not. Roy. Astron. Soc. , 399 (2002), arXiv:astro-ph/0106491 [astro-ph].[58] J. Einasto et al. , Astron. Astrophys. , A128 (2011),arXiv:1105.2464 [astro-ph.CO].[59] E. Jennings, Y. Li, and W. Hu, (2013),10.1093/mnras/stt1169, [Mon. Not. Roy. Astron.Soc.434,2167(2013)], arXiv:1304.6087 [astro-ph.CO].[60] K. C. Chan, N. Hamaus, and V. Desjacques, Phys. Rev. D90 , 103521 (2014), arXiv:1409.3849 [astro-ph.CO].[61] E. Adermann, P. J. Elahi, G. F. Lewis, andC. Power, Mon. Not. Roy. Astron. Soc. , 4861 (2018),arXiv:1807.02938 [astro-ph.CO].[62] E. Adermann, P. J. Elahi, G. F. Lewis, andC. Power, Mon. Not. Roy. Astron. Soc. , 3381 (2017),arXiv:1703.04885 [astro-ph.CO].[63] A. Dekel and O. Lahav, Astrophys. J. , 24 (1999),arXiv:astro-ph/9806193 [astro-ph].[64] A. Mesinger, S. Furlanetto, and R. Cen, Mon. Not. Roy. Astron. Soc. , 955 (2011), arXiv:1003.3878 [astro-ph.CO].[65] J. B. Mu˜noz, Phys. Rev. D100 , 063538 (2019),arXiv:1904.07881 [astro-ph.CO].[66] S. Foreman, W. Coulton, F. Villaescusa-Navarro, andA. Barreira, (2019), arXiv:1910.03597 [astro-ph.CO].[67] D. Nelson et al. , (2018), arXiv:1812.05609 [astro-ph.GA].[68] S. Genel, M. Vogelsberger, V. Springel, D. Sijacki,D. Nelson, G. Snyder, V. Rodriguez-Gomez, P. Torrey,and L. Hernquist, Mon. Not. Roy. Astron. Soc. , 175(2014), arXiv:1405.3749 [astro-ph.CO].[69] S. McAlpine et al. , Astron. Comput. , 72 (2016),arXiv:1510.01320 [astro-ph.GA].[70] I. G. McCarthy, J. Schaye, S. Bird, and A. M. C.Le Brun, Mon. Not. Roy. Astron. Soc. , 2936 (2017),arXiv:1603.02702 [astro-ph.CO].[71] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Sel-jak, and Z. Slepian, Astron. J. , 160 (2018),arXiv:1712.05834 [astro-ph.IM].[72] A. Mirizzi, J. Redondo, and G. Sigl, JCAP , 001(2009), arXiv:0905.4865 [hep-ph].[73] T. Moroi, K. Nakayama, and Y. Tang, Phys. Lett. B783 ,301 (2018), arXiv:1804.10378 [hep-ph].[74] K. N. Abazajian et al. (CMB-S4), (2016),arXiv:1610.02743 [astro-ph.CO].[75] A. M. Price-Whelan et al. , Astron. J. , 123 (2018),arXiv:1801.02634.[76] T. P. Robitaille et al. (Astropy), Astron. Astrophys. ,A33 (2013), arXiv:1307.6212 [astro-ph.IM].[77] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. , 473 (2000), arXiv:astro-ph/9911177 [astro-ph].[78] A. Lewis and S. Bridle, Phys. Rev. D66 , 103511 (2002),arXiv:astro-ph/0205436 [astro-ph].[79] Y. Ali-Haimoud and C. M. Hirata, Phys. Rev. D83 ,043513 (2011), arXiv:1011.3758 [astro-ph.CO].[80] F. Perez and B. E. Granger, Computing in Science andEngineering , 21 (2007).[81] T. Kluyver et al. , in ELPUB (2016).[82] J. D. Hunter, Computing In Science & Engineering , 90(2007).[83] S. van der Walt, S. C. Colbert, and G. Varoquaux,Computing in Science and Engineering , 22 (2011),arXiv:1102.1523 [cs.MS].[84] M. Waskom et al. , “mwaskom/seaborn: v0.8.1 (septem-ber 2017),” (2017).[85] W. McKinney, in Proceedings of the 9th Python in Sci-ence Conference , edited by S. van der Walt and J. Mill-man (2010) pp. 51 – 56.[86] P. Virtanen et al. , Nature Methods (2020),https://doi.org/10.1038/s41592-019-0686-2.[87] C. O. da Costa-Luis, JOSS , 1277 (2019).[88] G. D. Becker, J. S. Bolton, M. G. Haehnelt, and W. L.Sargent, Mon. Not. Roy. Astron. Soc. , 1096 (2011),arXiv:1008.2622 [astro-ph.CO].[89] J. S. Bolton, M. Viel, T. S. Kim, M. G. Haehnelt, andR. F. Carswell, Mon. Not. Roy. Astron. Soc. , 1131(2008), arXiv:0711.2064 [astro-ph].[90] A. Rorai et al. , Mon. Not. Roy. Astron. Soc.466