CMB Spectral Distortions from an Axion-Dark Photon-Photon Interaction

The presence of a plethora of light spin 0 and spin 1 fields is motivated in a number of BSM scenarios, such as the axiverse. The study of the interactions of such light bosonic fields with the Standard Model has focused mostly on interactions involving only one such field, such as the axion ($\phi$) coupling to photons, $\phi F \tilde F$, or the kinetic mixing between photon and the dark photon, $ F F_D$. In this work, we continue the exploration of interactions involving two light BSM fields and the standard model, focusing on the mixed axion-photon-dark-photon interaction $\phi F \tilde F_D$. If either the axion or dark photon are dark matter, we show that this interaction leads to conversion of the CMB photons into a dark sector particle, leading to a distortion in the CMB spectrum. We present the details of these unique distortion signatures and the resulting constraints on the $\phi F \tilde F_D$ coupling. In particular, we find that for a wide range of masses, the constraints from these effect are stronger than on the more widely studied axion-photon coupling.


Introduction
In the past several decades, overwhelming gravitational evidence for the existence of dark matter (DM) has been collected [1,2].However, we have yet to observe non-gravitational dark matter interactions with standard model particles.This has led to a wide range of models for particles that could describe dark matter.Of those, a class of very motivated models are ultralight bosonic dark matter models, where dark matter is an ultralight (m ≪ eV) scalar or vector field.
Ultralight bosonic particles, such as axions and dark photons, are well motivated from many points of view.Extra dimensional theories, such as string theory, typically predict a plethora of light scalars and vectors [3][4][5][6].In fact, already the pioneering work on extra dimensional models by Kaluza and Klein proposed that the electromagnetic gauge symmetry, and thus the photon, could be a consequence of extra-dimensions [7,8].Even without such motivations coming from extra dimensional constructions, both the axion and dark photon are highly motivated candidates for ultralight bosonic dark matter [9][10][11][12][13][14], and, in addition, the axion provides an elegant solution to the strong CP problem [15][16][17][18].There are many instances of dark matter consisting of ultralight dark photons produced non-thermally [19][20][21][22][23][24].One such method for non-thermal production of dark photons is from axions via the coupling in Eq. 1.1 [25], providing extra motivation to carefully consider this interaction.
Given the expectation that there might be many light bosons, one of which can play the role of dark matter, one expects interactions between these light particles themselves, as well as interactions between them and the Standard Model.In this paper, we consider models with two new particles; an axion and a dark photon with a coupling to the standard model of the form Where ϕ is the axion and F (F D ) is the photon (dark photon) field strength tensor.Generically, in models with axions and dark photons, interactions of the form ϕF F and ϕF FD , as well as kinetic mixing between the photon and dark photon (F F D ) would be present, and could be more relevant for detecting the light bosons.However, if one demands that there is a dark charge conjugation symmetry C D under which the axion and dark photon are odd, these other couplings are absent, or highly suppressed, if there is a small breaking of the symmetry.In Appx.A we present a simple model that exhibits such a symmetry and show that the coupling given in Eq. 1.1 is the leading one.This coupling has been studied in a variety of scenarios [26][27][28][29][30][31][32][33][34][35][36][37].
Even if one considers models with the ϕF F interaction, we show the bounds placed from the ϕF FD interaction can be stronger.To see this, note bounds placed on a ϕF F interaction from the Cosmic Microwave Background (CMB) in the presence of an axion dark matter background are placed from polarization measurements [38].Conversely, as we will see in this work, bounds on the interaction in Eq 1.1 are placed from measurements of the CMB frequency spectrum.The CMB spectrum has been measured more precisely than the CMB polarization, implying that the constraints on ϕF FD will be stronger than those on ϕF F .Thus, if one considers both interactions to have similar strength, the bounds placed on ϕF FD will be stronger.
The goal of this paper is to investigate how this coupling affects the CMB when either the axion or the dark photon is dark matter.In the early universe, before redshift of z ∼ 1100, the universe was hot and dense enough that photons, electrons and protons were all in thermal equilibrium with one another.Once the universe cooled to redshift z = 1100, effectively all electrons were bound to nuclei forming neutral atoms, making the universe transparent to photons, in an era called recombination.As the number of free electrons decreased with the lowering temperature, the mean free path of photons λ γ increased.Around the time of recombination (z ≈ 1100), λ γ ∼ H and the photons transitioned from being trapped in the electron-baryon plasma to being free streaming.Afterwards, these photons could propagate freely until being detected by CMB experiments.Due to the early thermal equilibrium of these photons, their power spectrum follows that of a blackbody.In the early 1990's, the Cosmic Background Explorer (COBE) [39] satellite equipped with Far Infrared Absolute Spectrophotometer (FIRAS) performed the most accurate measurement of the CMB monopole power spectrum.They found it matched a blackbody with temperature 2.7255 K to agree to about 1 part in 1,000 to 10,000 [40], making the CMB monopole power spectrum one of the most precisely measured cosmological observables.Any phenomena that would distort this spectrum is then highly constrained by COBE-FIRAS.Constraints on kinetic mixing [41] and dark matter interactions [42][43][44][45][46][47] from the COBE-FIRAS data have been placed using these spectral distortions 1 .
The interaction given in Eq 1.1 can cause CMB spectral distortions.In the presence of an axion dark matter background, this interaction allows photons to be converted into dark photons.Likewise, in a dark photon dark matter background, it allows photons to be converted to axions (see Ref. [35] for an early study of this effect in the resonant regime).These dark sector particles are invisible to us and thus the effect of Eq. 1.1 in both cases is to remove photons from the CMB spectrum.This removal of photons naturally distorts the observed CMB spectrum.In Sec. 2, we will show how the time at which these photons are removed gives rise to different types of distortions and argue that the size of these distortions depends on the probability of removing a photon from the spectrum.In Sec. 3, we will compute this probability from the interaction in Eq. 1.1 and in Sec. 4 use it to compute the various types of distortions.In Sec. 5 we place constraints on the coupling 1/f a by comparing these distortions to the COBE-FIRAS data and briefly comment on the possible shapes of the distortions.We conclude in Sec 6.
In this section, we describe how our interaction gives rise to CMB spectral distortions.The effect of our interaction, in the presence of dark matter, is to convert photons into a dark sector particle, X.For example, if the axion constitutes the dark matter, photons will interact with the axions and convert into dark photons so that X is the dark photon.Conversely if the dark photon is dark matter, the photon will convert into axions and, in this case, X is the axion.In order to speak generally about either axion or dark photon dark matter, we will refer to the particle the photon converts to as X throughout this paper.The implications for the CMB spectral distortion are the same in either scenario since the important effect is that the photons convert to an invisible dark sector particle X and are removed from the photon spectrum.
The removal of any photons from the bath can lead to a deviation from the blackbody spectrum.We can quantify that change by a frequency dependent distortion δ(ω) defined in Eq. 2.1.
The exact frequency dependence of the distortion will depend on when in cosmic history the photons were removed from the photon spectrum.As shown in Figure 1, we can define 5 different eras, the T era, the µ era, the µ − y transition era, the y era, and the free streaming era, in which the injection or removal of photons gives rise to different distortions.In the remaining of this section, we briefly review these different eras and discuss the characteristic effect of photon removal in each.As we will discuss, the final distortion to the blackbody spectrum can be parameterized by the impacts coming from different eras as δ Tot (ω) = (μ + μt )M (ω/T ) + (y + y t )Y (ω/T ) + δ Doppler (ω) + δ Free (ω). (2.2) We see there are 4 different types of distortions.The µ and y distortions have distinct shapes and are insensitive to the details of the model generating the distortion, while δ Doppler and δ f ree have a model dependent shape.The contribution of all the pre-recombination distortions are computed using the Green's function method described in Ref. [49], using the rate at which photons are converted into particle X, Γ γ→X (ω), as discussed in Sec. 4 and Appx.D.2.In computing this rate, we will need to calculate is the conversion probability P γ→X (ω).We also show that the distortion due to the post-recombination free streaming era distortion is directly related the the conversation probability P γ→X (ω).Thus, the central quantity we will need to compute all spectral distortions is the conversion probability.

T Era
In the very early universe, at temperatures T ≳ 0.5 keV (redshifts z ≳ 2 • 10 6 ), a number of scattering processes involving photons are very efficient at driving the photon distribution towards an equilibrium distribution with zero chemical potential.The main processes, shown in Fig. 2, are Compton scattering, which can quickly redistribute the photon energy and equilibrate the electron and photon temperatures, as well as number changing processes, such as double Compton and bremsstrahlung.Due to these processes, in this era, even if photons are lost due to conversion into X, the distribution would be quickly driven back to that of a blackbody and the only effect would be a small change in the blackbody temperature.Since we don't know a priori the temperature of the CMB, injections in this era would not lead to a bound from the CMB spectrum (there could be bounds by studying the anisotropies or comparing to Big Bang Nucleosynthesis predictions which we will not explore in this work).

µ Era
Once the universe has cooled past T ≈ 0.5 keV (z ≈ 2 • 10 6 ), the higher order processes like double Compton scattering and bremmstrahlung are no longer efficient at setting the chemical potential to zero (although they can still be efficient for absorbing/emitting very low frequency photons).In this era, Compton scattering is still efficient at redistributing the energy, which drives the spectrum towards an equilibrium distribution.Because photon number is now conserved, any removal of energy will result in a small effective chemical potential term μ This distortion has a fixed shape, M (ω/T ), given in Appx.D.2.The size of this distortion is captured by the effective chemical potential, μ, which can be calculated from Γ γ→X , the rate at which photons are being converted to dark sector particles X as shown in Appx.D.2.COBE-FIRAS placed a bound of |μ| < 9 • 10 −5 [40].

y Era
At temperatures lower than 3 eV (z ≈ 10 4 ), Compton scattering is still efficient enough to trap photons, but is now inefficient at changing photon energy, and transferring energy between photons and electrons.This leads to two effects.Firstly, some small amount of energy can be still exchanged with the electrons, leading to a difference in the photon and electron temperatures.Subsequent scatterings of photons with an electron fluid at a different temperature, lead to a y-distortion via the Sunyaev-Zeldovich (SZ) effect [50].Secondly, energy injections/removals in a given frequency, can still be smeared due to Doppler broadening via Compton scattering, even if the process is no longer efficient at thermalizing the spectrum.This leads to two separate distortions: a y-distortion and a Doppler smeared distortion, The y-distortion has a fixed shape, Y (ω/T ), given in Appx.D.2 and a size determined by the small parameter y which can be computed from the photon loss rate Γ γ→X as described in Appx.D.1.COBE-FIRAS placed a constraint |y| < 1.5 • 10 −5 [40].On the other hand, the shape of the Doppler smeared distortion is model dependent, so we instead place a bound by comparing directly to the COBE-FIRAS data.

µ − y Transition Era
Once the temperature decreases below T ≈ 70 eV (z ≈ 3 • 10 5 ), Compton scattering, while still efficient at trapping photons, begins to become inefficient at redistributing energy for certain frequency modes of the photon spectrum.This signals the end of the µ era and the start of the µ − y transition era which lasts until T ≈ 3 eV (z ≈ 10 4 ).In this transition era, higher energy modes still redistribute energy efficiently through Compton scattering, while energy redistribution is inefficient for lower energy modes.At intermediate modes, energy redistribution is not efficient but is non-negligible.In order to exactly treat this very frequency dependent behavior, one would need to simulate the distortion numerically [49].However, as noted in Ref. [49], for the range of photon frequencies we are interested in, the distortion can be modeled to good accuracy as a pure energy injection as described in Ref. [51].The resulting spectral distortion is a combination of a µ distortion M (x) and a y distortion Y (x) the shapes of which are given in Appx.D.2.
The subscripts t on the coefficients μt and y t denote that these coefficients are calculated differently from μ in Eq. 2.3 and y in Eq. 2.4.They still however are calculated from the photon conversion rate Γ γ→X as shown in Appx.D.2.

Free Streaming Era
Around T = 0.25 eV (z = 1100), most electrons have been captured to form neutral hydrogen, and the universe becomes transparent to photons.From this point on, the photons become free streaming and can travel unimpeded across the universe, giving rise to the CMB we observe today.However, the presence of our interaction leads to a probability P γ→X (ω) that a CMB photon with frequency ω will convert to an invisible dark sector particle X before reaching us.Because the photons are free streaming, there is no thermalization, or redistribution of energy.So, the resulting spectrum is the original spectrum multiplied by the survival probability of a photon to reach us without converting to X, We can see that the distortion is simply the conversion probability P γ→X (ω).The frequency dependence of this distortion is model dependent and as such we will have compute it and then constrain it directly with the COBE-FIRAS data to obtain a bound.

Transition Probability
As described in the previous section, to compute the distortions we will need to compute the probability, P (ω, t, t 0 ), of converting photons of frequency ω produced at time t 0 into dark sector particles X at some later time t.Because our dark matter is made of bosons of mass m ≲ meV, the number density is large enough to treat dark matter as a classical background field.Therefore, these probabilities can be computed using Feynman diagrams like the one shown in Fig. 3.In many instances the time interval t − t 0 will be of cosmic scale.For example, when applied to the free streaming distortion, t − t 0 will be the time between recombination and the present.As a consequence, we will need to compute these Feynman diagrams in a curved FRW background.We will work in conformal coordinates, described by the metric In a slowly expanding universe (m DM , T CM B ≫ H), we can easily expand any general scalar field ϕ and vector field A µ in terms of ladder operators by solving their equations of motion using the WKB approximation.The details of this process are given in Appx.B.1 and the result is where a k (a λ k ) are the ladder operators of the scalar(vector) field, ω c (η, k) ≡ |k| 2 + a 2 (η)m 2 is the conformal energy, and the polarizations ϵ λ µ (k) are given in Eq.B.6.We will also need to include plasma effects on the photon due to its impact on photon propagation and mixing.At all times of relevance, electrons are non-relativistic and lead to a plasma frequency, ω p (η), given by where α e is the fine structure constant and n e (η) is the number density of electrons which changes as the universe expands.We are working in the limit where ω p (η) is much smaller than the frequency of the CMB photons ω ∼ T CM B .In this limit, the effects of the plasma can be reduced to the photons acquiring a small mass m γ (η) = ω p (η) ≪ T CM B .Despite the plasma frequency giving rise to an effective mass for transverse modes, at such large frequencies there are no longitudinal modes of the photon (plasmons) [52].We use the redshift dependent plasma frequency from [41].
Due to the non-trivial time dependence of the fields given in Eq. 3.2 and Eq.3.3, we will only Fourier transform the diagrams in space.In this regime, our diagrams are transition amplitudes with a time dependent interaction and so we can expand our amplitudes to leading order using a Dyson series, V I (η ′ ) is our interaction potential given by Note the momentum eigenstates in Eq. 3.5 are normalized such that ⟨k ′ |k⟩ = (2π) 3 δ(k − k ′ ) which differs from the usual Lorentz invariant normalization by a factor of 2ω.We can simplify V I (η ′ ) by using the fact that our dark matter field is nonrelativistic to ignore gradients in favor of time derivatives which simplifies the interaction to Axion DM: Next we insert the expansion of the fields for the photon's magnetic field B (Eq. 3.3) and the outgoing particle field X (either Eq. 3.3 for an outgoing dark photon or Eq.3.2 for an outgoing axion) in terms of creation and annihilation operators.As discussed earlier, we will treat the dark matter as a classical background field.These potentials can be inserted in Eq. 3.5 to compute the transition probability.After some simplifications described in Appx.B.1, this probability takes the form, with where the first (second) equation corresponds to the axion (dark photon) dark matter scenario.
In the above equations, x(t ′ ) is the position of the photon at time t ′ , the dots represent time derivatives with respect to comoving time t, and v(k, t) represents the velocity of a given particle.We have approximated v γ = 1.The ω's are now the physical energies defined as: Finally, we must determine what form our dark matter background takes.By solving the equations of motion for the dark matter fields (Eq.B.3 and Eq.B.4) in the non-relativistic limit, keeping terms up to O(v DM ), and demanding that the energy density is ρ 0 DM /a 3 (t), where ρ 0 DM is the energy density of dark matter of the universe at the present time, we find Both fields get a spatially dependent phase, β( x), while the vector dark matter field gets an additional spatially dependent polarization unit vector ϵ( x).Both of these quantities vary spatially on the scale of the dark matters de Broglie wavelength (m DM v DM ) −1 with v DM ≪ 1.Additionally, they vary in time, on timescales (m DM v 2 DM ) −1 .Since the time dependence is suppressed by a factor v DM relative to the spatial dependence it will be ignored.In Appx.B.3, we show that for all distortions, we average x over many de Broglie wavelengths of the dark matter field, which means that we can average all of these spatially dependent quantities.We will leave the averaging over the phase for later in the computation, but in Appx.B.3 we show that we can effectively replace Physically, this factor 1 √ 3 is reflecting the fact that the interaction ϵ λ (k) • ȦD in Eq. 3.10 is picking out one particular polarization of the vector dark matter.After averaging over x, this particular polarization must make up 1/3 of the total dark matter by isotropy, effectively sending ρ 0 DM → ρ 0 DM /3.Using Eqs.3.12-3.14to simplify Eqs.3.8-3.10,we can write the transition probabilities as where L is a length scale defined as ) , where ∆ω γ→X is the change in energy from a photon converting into particle X at momentum k, The ⟨⟩ β indicates the remaining average over the dark matter phase β(t ′ ) = β( x(t ′ )) which is handled for pre-recombination and free streaming distortions separately in Appx.D.1 and C respectively.From Eq. 3.15, we can see that the only difference between scalar dark matter and vector dark matter is the overall factor of 1/3 in the conversion probability from the effect described above.This means that the coupling to dark photon matter is effectively 1/ √ 3 that of the coupling to axion dark matter and so the bounds placed on the coupling in the dark photon dark matter will be weaker than the bounds for the axion dark matter by a factor √ 3 .For simplicity, we will only consider axion dark matter going forward, knowing that we can translate any result to dark photon dark matter by multiplying by 1/ √ 3.

Computing the Distortions
In this section, we will use Eq.3.15 to determine the strength of the distortions arising from the various eras.This will be very different for distortion generated pre-recombination versus in the free streaming regime, so we consider them separately.

Free Streaming Distortion
In Sec.2.5, we showed that the free streaming distortion, δ f ree , is equal to the probability of converting the photon to dark sector particle X between recombination, t 0 and today, t.Thus, we need to compute From Eq. 3.16, it is clear that we can write , where The remaining integral is an oscillatory integral with frequency Notice that all of the time dependent quantities change on the Hubble scale H(t) due to the expansion of the universe2 There are two limits in which this integral can be computed.The first is the fast oscillation limit where the oscillation frequency Ω ± (t) is approximately constant over many oscillations.This is the limit where Ω ± (t) ≫ H(t). (4.4) In this limit, all of the time dependent quantities in Eq 4.2 become approximately constant up to corrections of order H/Ω ± and the integral can be computed analytically.The second limit is the resonant limit.In this limit, there is a time (or possibly multiple times), t r , where there is a stationary phase in the exponential (Ω ± (t r ) = 0).Since β(t) ≪ m DM to a good approximation these resonant times can be found by solving the equation Physically, this corresponds to times in which the dark matter particle being absorbed/emitted by the photon is on shell, leading to an enhancement in the conversion probability.The stationary phase approximation is used to compute the integral in this limit.Appx.C contains the details of computing the distortions in both of these limits.In the end, we find, Here, a ± are all solutions to Eq. 4.5, and a * = (1090) −1 is the redshift at recombination.Notice that the resonant distortion is enhanced by a factor of |∆ω±m DM | 2 ∂t∆ω ∼ m DM H with respect to the non-resonant distortion.Even for the smallest possible dark matter masses, m DM ∼ 10 −20 , this is an enhancement by a factor of O(10 8 ).Thus, we expect our bounds on the coupling 1/f a to be enhanced by orders of magnitude in regions of parameter space where these resonances happen.

Pre-Recombination Distortions
The application of Eq. 3.15 to pre-recombination distortions is not as straightforward as for the free streaming distortions.The distortions can be computed via the Green's function method outlined in Ref. [49] where the distortion, δ(x) is given by where x is the dimensionless frequency x = ω/T (a) and Γ γ→X (x ′ , a) is the rate at which photons of frequency x ′ T are converted to dark sector particle X.The Green's function G(x, x ′ , a) describes how photons injected into mode x ′ at time when the scale factor is a are redistributed to mode x.G(x, x ′ , a) is given for the various eras in Appx.D.2.In this section we will describe how to compute the rate Γ γ→X (x ′ , a) appearing in Eq. 4.7.
We can compute Γ(a, k) from P (k, t 0 , t) as follows.Consider a photon that scatters off of an electron at time t 0 and travels some time τ before scattering off of another electron.P (k, t 0 , t 0 +τ ) is then the probability of the photon converting to X between these scatterings.For an ensamble of photons scattering with time τ between scatterings, the rate at which those photons are converted to X is then Since the photon is traveling through a very dense medium of electrons, it has a certain probability p(τ ) of traveling a distance τ characterized by its mean free path λ γ , given by In order to find the average rate for all photons, we average Γ(t 0 , τ, k) over this path length distribution Because this era is before recombination, we have Hλ γ ≪ 1.In this limit, we can treat space as static and all parameters that change due to the expansion of the universe as constant in the integral, and Eq.4.10 can be computed analytically.The details are given in Appx.D.1, and the end result is where 2  .
This rate can be plugged into Eq.4.7 and integrated numerically using the Green's functions given in App D.2 find the distortion from the various pre-recombination eras.

Results
The total distortion for a given set of masses m X and m DM and coupling 1/f a is given by Eq. 2.2.As described in Sec. 2, temperature shift distortions are undetectable by COBE-FIRAS.Therefore we should add an arbitrary temperature shift, T (x) (defined in App.D.2), to the distortion and do a best fit to COBE-FIRAS data [40] with both the coupling, 1/f a and the size of the temperature shift, α, as free parameters.
δ Tot (ω) = (μ + μt )M (ω/T ) + (y + y t )Y (ω/T ) + δ Doppler (ω) + δ Free (ω) − αT (ω/T ).( 5.1) However we can simplify this by demanding that the number density of CMB photons be the the same as that of a perfect blackbody at the measured temperature T CM B = 2.35 • 10 −4 eV.This is exactly the procedure that is commonly done for µ and y distortions as described in [53].This constraint fixes the size of the temperature distortion α.For fixed m D and m ϕ , we do a χ 2 fit of our distortion to the COBE-FIRAS data with a single free parameter, 1/f a .By demanding that the distorted spectrum matches the measured spectrum to within 2σ, we obtain bounds for the coupling as a function of the dark photon and axion masses m D and m ϕ .Fig. 4 shows a contour plot of these bounds as a function of the axion and dark photon masses.We show contours for both axion dark matter and dark photon dark matter.These bounds are plotted against the leading best bound on this coupling from red giant cooling constraints.Ref. [37] shows that this coupling leads to a novel cooling mechanism in red giants due to plasmon decay and deduced that the cooling from this coupling is equivalent to the cooling from a neutrino magnetic dipole moment µ ν = 1 2fa .Then using the bound placed on neutrino magnetic dipole moments found in [54] from red giant cooling they were able to place a bound 1/f a < 7.1 • 10 −10 GeV on the axion-photon-dark photon coupling.As seen in Fig. 5, in a large portion of parameter space, roughly m X < 10 −2 eV and m DM < 10 −8 eV, our bounds beat this red giant bound by several orders of magnitude.It is worth noting that while the red giant bound is the most stringent bound in this region of parameter space (aside from the bounds placed in this work), there have been numerous other constraints placed on the coupling in this region of parameter space.The bounds placed on this coupling from stellar evolution [29], Horizontal Branch stars [35], and white dwarfs [36] are all with in an order of magnitude of those from red giants.For simplicity, we only include the red giant bound in our plots.
Fig. 5 shows our bounds as a function of the axion dark matter mass (m DM ) for selected values of the dark photon mass (m X ) and shows the contribution to these bounds from each distortion era.While the bounds in Fig. 5 are shown for axion dark matter, the equivalent bounds for dark photon dark matter can be found by scaling the bounds up by a factor of √ 3 as discussed at the end of section 3.As can be seen from the colored dashed lines in Fig. 5 each constraint from each distortion has roughly the same behavior: constant for small m DM and increasing linearly in m DM for large m DM with an enhanced region in between.We can understand why the bounds have this behavior.Firstly, the enhanced region is the region of parameter space where photons during that particular era are able to resonantly convert.The other two limits can be understood by considering the integral L 2 defined in Eq. 3.16 in the non-resonant regime.Here L 2 , the effective oscillation length, is the square of an oscillatory integral and thus should scale as L 2 ∼ Ω −2 where Ω is the fastest oscillation frequency in the integral.For sufficiently small m DM , the coherent oscillation of dark matter is unimportant.As such, the oscillation length is the standard ∆ω γ→X present for well known systems such as neutrino oscillations.For larger m DM , the fast oscillation of dark matter dominates and the oscillation frequency is m DM .This combined with Eq. 3.15, shows that the conversion probability P γ→X ∝ ρ DM f 2 a ∆ω 2 for small m DM and

DM
for large m DM .Since the distortions all scale with the conversion probability, its clear that the bounds at low m DM are independent of m DM and linearly proportional to m DM for large m DM .
Given that the COBE-FIRAS data was collected over 30 years ago, current technology could measure the CMB spectrum to higher precision.In fact, there are proposals for experiments, like PIXIE [55], that aim to measure the spectrum to within a factor of 10 −8 − 10 −9 , an improvement of around a 3 to 4 orders of magnitude from COBE-FIRAS.These future experiments could potentially measure a distortion in the CMB frequency spectrum and thus it is interesting to ask what such a measured distortion could tell us about our dark matter models.
Specifically, we will discuss qualitatively whether a distortion produced from our dark matter model(s) could potentially be distinguished from other distortion sources.Energy injection or removal into or from the background electron plasma before recombination produces primarily a µ and/or y distortion and thus these types of distortions which are essentially model independent.However, because of their non-thermal origin, the Doppler distortion arising pre-recombination, and free streaming distortion post-recombination have model dependent spectral shapes and do provide a distinctive signature.For simplicity, we can focus .Bounds on the coupling 1/f a as a function of the axion dark matter mass m DM plotted for various values of the dark photon mass, m X .The purple, blue, green and orange dashed lines are the individual µ, y, Doppler and free streaming bounds respectively while the solid black line is the total bound from all distortions.The grey line represents the current best bound on the coupling derived from red giant cooling [37].The bounds shown here are for the axion dark matter case, but the bounds for the dark photon dark matter case can be obtained by scaling the bounds up by a factor of √ 3. The bounds for m X < 10 −13 eV are exactly those given in the upper-left plot.
on the free streaming distortion to get a sense of the various types of shapes this distortion can take.To start, one can take the large and small m DM limits (m DM ≫ ∆ω and m DM ≪ ∆ω respectively) of Eq. 4.6 and see that Small m DM limit: δ ∼ ω 2 Large m DM limit: δ ∼ Constant. (5.2) Figure 6 shows these quadratic and constant distortions plotted against the µ and y distortions.The amplitudes of these distortions in Fig. 6 are chosen so that each distortion disagrees with COBE-FIRAS at 2σ.In this sense we can think of these distortions as being of equal strength.We can see that the large m DM matches very well with the µ-distortion and the small m DM distortion, matches well with a y distortion making them difficult to distinguish from the generic µ and y distortions respectively. .Each distortion is plotted with an amplitude such that it disagrees with the COBE-FIRAS data at 2σ.Thus all the distortions are effectively the same strength.The purple and blue lines are µ and y distortions respectively while the dashed green and orange lines depict Doppler and free-streaming distortion shapes for different choices of masses m DM and m X .On the left, we show the distortions in the non-resonant limits given in Eq. 5.2.On the right, we show two of the many possible shapes the distortions can take when there is a resonance in either the Doppler or free distortions.These resonant distortions have shapes distinct from the µ and y distortions.It is also worth noting that the differences in shape between the Doppler and free distortion in the right-hand plot are due to the difference in choice of parameters, m DM and m X , rather than a difference in distortion type.
The shapes of the free streaming distortion get more interesting if we consider the resonant region of parameter space (m DM ∼ ∆ω).Here the distortion depends on the time of the resonance a ± which is found by solving Eq. 4.5.Because ∆ω depends of the frequency of the photon, so does the resonant time through Eq. 4.5.Thus different frequency modes can have different resonance times which can lead to very distinctive frequency dependencies in the distortion easily distinguishable from the standard µ/y-distortion.In particular, it is possible that some frequency modes undergo resonance, while other modes don't.This leads to especially unique distortions, where some frequency modes are distorted while others, effectively, are not.Such extreme distortions are shown in the orange and green dashed lines in figure Fig. 6.As can be seen, lower frequency photons never resonate, and thus are effectively undistorted while higher frequencies do experience a distortion due to resonance.Additionally, the frequency dependence of the resonant piece of the distortion is not a simple power law of ω due to the dependence of the resonant time on the frequency and in turn the non-trivial dependence of the distortion on the resonant time.It is also important to note that these these qualitative features can also arise from the Doppler distortion.From Eq. 4.11 and 4.12 one can derive the same small and large m DM behavior and show similar types of resonant behavior are possible.The green dashed line in Fig. 6 shows one such distinctive resonant shape for the Doppler distortion.This makes the prospect of observing these unique distortions even more likely since, as shown in Fig. 5, there are significant regions of parameter space where the Doppler distortion leads to the strongest bound, which shows that in such regions of parameter space it is the most observable effect.It is worth pointing out that difference in shape and severity of the jumps of the Doppler and free distortions in Fig. 6 is not due to an inherent difference between the free and Doppler distortion, but rather a difference in the parameters m DM and m X at which these distortions are evaluated.These set of parameters were chosen to highlight the difference in distortion shapes achievable by either the free or Doppler distortions rather than an inherent difference between them.

Conclusion
In this paper, we studied the effects from an axion-photon-dark-photon coupling to the Cosmic Microwave Background if either the axion or the dark photon is dark matter.This interaction, in a dark matter background, induces mixing of the photon with a new light boson, and can remove photons from the universe either before or after recombination.Removing photons from the baryon-photon plasma before recombination produces the well known, model independent, µ distortion or y distortion, as well as a model dependent distortion due to Doppler broadening.The size of the distortion is determined by the rate at which photons are removed from the spectrum.Removing photons after recombination naturally changes the frequency spectrum of the CMB and thus produces a distortion with a new spectral shape.We computed these individual distortions and the corresponding total distortion that would be produced by an interaction of this type in the presence of either axion or dark photon dark matter.
The distortion produced in the presence of dark photon dark matter is smaller than that produced by axion dark matter by a factor of 1/3, but otherwise identical in terms of dependence on the model parameters.This is due to the energy of dark matter being spread over the 3 polarizations of dark photon dark matter, as opposed to the single polarization for axion dark matter.The interaction we study only couples a photon of given polarization to a single polarization of dark photon dark matter.Thus, a given photon effectively only couples to 1/3 of the total dark matter background.This effectively leads to an interaction strength for dark photon dark matter that is 1/3 of that for axion dark matter.
By comparing our computed distortions with the COBE-FIRAS data measuring the CMB frequency spectrum, we were able to place very restrictive bounds on our coupling 1/f a .These bounds are a significant improvement of several orders of magnitude over the previous best bound on this coupling from red giant cooling constraints [37] as shown in Fig. 5. Additionally one can compare these bounds to those placed on the ϕF F coupling in [38] and see that our bounds are several orders of magnitude stronger.
We also briefly considered the possibility for future measurements of the CMB frequency spectrum to detect distortions produced by this model.While the µ and y distortions are produced by any mechanism that adds or removes energy from the photon spectrum before recombination, the distortions generated in the free streaming era when the resonance condition is met lead to much more distinctive spectral features.Thus, such resonant distortions offer a promising avenue to single out the axion-photon-dark-photon interaction interaction if the next generation experiments measures distortions on the CMB spectrum.

A A Model
Here we summarize a simple model containing an axion, photon and dark photon where the leading interaction is that given in Eq. 1.1.This model contains a dark sector with a complex scalar Φ a dark sector U (1) gauge boson A D and two sets of two left handed Weyl fermions ξ, ξ c , χ and χ c .These fermions are charged under both electromagnetism and the dark U (1) gauge group with charges shown in Table 1.The scalar Φ is uncharged under both electromagnetism and the dark U (1).
Table 1.The charges for our dark sector fermions.Q EM is the particles electric charge in units of the fundamental electric charge and similarly, Q D is the particles charge in units of the fundamental dark charge.
In addition to these charges we demand that our model obey a dark charge conjugation symmetry C D defined by and a Z 8 symmetry defined by This Z 8 , which is a subgroup of the would be U (1) P Q associated with the axion, may seem troublesome since it corresponds to a chiral rotation of the χ and ξ by −π/4 and π/4 respectively and thus should generate an anomalous term through the chiral anomaly.However, with the charges defined in Table 1, one can show that this is equivalent to a rotation of a single Weyl fermion with charges Q EM = Q D = 1 by an angle of 2π.Our theory should be consistent with the addition of a Dirac fermion with fundamental electric and dark charge and so a chiral rotation of this fermion by 2π must leave the theory invariant.Thus the Z 8 symmetry is nonanomalous.Additionally, this Z 8 symmetry forbids any operators in the potential V (Φ) that could give the axion a mass up to dimension 8.With these symmetries in mind, we write all possible terms in our Lagrangian up to dimension 4, where V is the scalar potential, the last term is a Yukawa coupling for the fermions to the scalar and L kinetic contains all of the kinetic terms for the fermions, Φ, and A D with the gauge couplings packaged in covariant derivatives.Now let us suppose that Φ under goes symmetry breaking and obtains a VEV, Φ(x) = (v + h(x))e iϕ(x)/f .(A.4) Amongst other changes, the Yukawa piece becomes These phases can be eliminated by a chiral rotation, which naturally generates anomalous terms in the Lagrangian due to the chiral anomaly.It is a simple exercise to show that these terms are We get two terms, one for each rotation.Plugging in the charges given in Table 1 we can easily see that the F F and F D FD terms vanish while the F FD term remains (A.8) Figure 7.At 1-loop order, it is easy to see that the diagrams for kinetic mixing exactly cancel due to the χ and ξ particles having opposite charges.
A similar cancellation happens in the kinetic mixing term.At one loop, as shown in Fig. 7, there are two diagrams for the kinetic mixing term which exactly cancel due to the χ and ξ particle's opposite dark charges but identical masses from the Φ VEV.At the heart of this cancellation is the dark charge conjugation symmetry, C D .This symmetry can easily be seen to forbid the generation of a kinetic mixing terms, which shows the cancellation observed at 1-loop occurs to all orders in perturbation theory.

B Computation of the conversion probability
In this appendix we summarize the details of our conversion probability computation.

B.1 Quantized FRW Fields
Here we give a brief description of general massive scalar and vector fields in an expanding FRW background metric dτ 2 = a 2 (η)(dη 2 −dx 2 ).We can expand the fields in ladder operators with mode functions u(η, k) for the scalar and v µ λ (η, k) for the vector where λ is a polarization index.
The mode functions u(η, k)e ik•x and v λ µ (k)e ik•x satisfy the equations of motion for scalar and vector fields respectively.The equations of motion are, To simplify, we make use of the fact that rate at which the universe is expanding (H ≡ ∂ η ln(a(η))) is much slower than the rate at which our fields are oscillating, which is roughly the comoving CMB temperature T CM B .In this limit one can show that up to corrections of O(k/H), The mode functions are normalized by making sure they reduce to the familiar flat space mode functions in the flat space limit.Here the ϵ λ µ represent the 3 different polarizations for the vector which are the usual transverse and longitudinal polarizations with m → a(η)m.Inserting Eq.B.5 into Eq.B.1 gives Eq. 3.2 and 3.3.

B.2 The Interaction Potential
Here we simplify the interaction potentials given in Eq. 3.7.We begin as described, by inserting field operators for the incoming photon and outgoing X particle, and a classical field background for the dark matter field, leading to and recalling that V ϕ (V D ) is the interacting potential when the axion (dark photon) is dark matter.
The terms not explicitly written in Eqs.B.7 and B.8 represent different combinations of the ladder operators that will be irrelevant for us since we only wish to consider photons as the initial state and outgoing axions/dark photons in the final state.In the regions of parameter space we will be interested in, the dark matter mass will be smaller than the CMB temperature.Since dark matter is also non-relativistic, this means that the momentum transfer from the dark matter, q ≡ p − p ′ , must be small with respect to the photon momentum, so we can do an expansion in small |q|.To lowest order, this means setting p = p ′ everywhere.However, we must keep q to linear order in the x • (p − p ′ ) order since q • x is not necessarily small since we will be integrating over all x. x Now, the only q dependence is in the exponent and we can shift variables d3 p ′ → d 3 q.After integrating over d 3 q we get a delta function δ 3 (x − x) where x = pη. 3This can then be used to eliminate the d 3 x integral.
Now, in order to simplify the cross products, we work in the helicity basis for photon polarizations where the following identities hold.
where λ = ±1 is the helicity of the photon.This helps simplify the expressions to Next, we can look at V I 's matrix elements with momentum eigenstates.After some simplification, it is easy to see these matrix elements take the form, where, the (λ ′ ) and (δ λλ ′ ) are meant to be included if the final state is a dark photon and M λ (η ′ , k) are given by Eq. 3.9 and 3.10.After dropping the factor (2π) 3 δ 3 (k−k ′ ) due to state normalization, Eq.B.17 can then be inserted into Eq.3.5 and squared to yield Eq. 3.8-3.10for the conversion probability.

B.3 Spatial Averages
In this section, we argue that in all distortions we are averaging the interaction position x over many de Broglie wavelengths of the dark matter field.This can be easily justified given that we are interested in the monopole spectrum, and thus will average distortions over all directions.This effectively means we will be averaging over many de Broglie wavelengths of the dark matter field.For the pre-recombination distortions, we are interested in the averaged conversion rate as a function of redshift, which depend on the conversion probability between photon scatterings.Thus, we average over all possible photon trajectories, and thus x everywhere in space.
For the free streaming case, one can make a more general argument, which shows that even if one is interested in anisotropies, this averaging is justified.First note that the smallest dark matter mass we consider is 10 −22 eV, and so the largest de Broglie wavelengths we must consider are ≲ 10 kpc, which is much smaller than the horizon size today.Because the probability conversion depends on the dark matter density, it is dominated at larger redshifts.This means that for photon conversions happening at similar times by directions separated by δθ ∼ 1, the distance between the transition points is ∼ H −1 0 * a t ≫ 10 kpc, where a t is the redshift of the transition.
Finally, let us average the dark photon dark matters polarization ϵ(x) over x and derive the replacement given in Eq. 3.14.To start notice that when plugging Eq. 3.10 into Eq.3.8, There will be a factor that looks like If we call the k direction the z direction, then this is simply, Now we average over all possible z-compontents of the dark photons polarization which gives, Thus the effect of averaging is simply to send ϵ λ (k) • ϵ( x) → 1/ √ 3 as described in Eq. 3.14.

C Free Distortion Computation
In this section we will detail the computation of the integral L 2 defined in Eq. 3.16 which we showed can be written as The strategy to computing this integral is to break it into intervals in which we can use either the fast oscillation limit or the stationary phase approximation.By choosing the boundaries of this regions appropriately, we can piece these intervals together to get the full result.We first identify any resonant times t r by solving Eq. 4.5 for t r .We can then break up the time interval into sub intervals as shown in Fig. 8.We will choose the endpoints of these intervals t r ± ∆t so that the following 3 conditions are true.
1.The stationary phase approximation should be valid everywhere inside the interval (t r − ∆t, t r + ∆t) such that we can expand the phase of the exponential to second order In order for this approximation to be valid, we must be able to ignore the third order term.This means we need 2. In order to match regions where the stationary phase approximation is valid to regions where the fast oscillation condition is valid, we want both the fast oscillation condition, Eq. 4.4, and the stationary phase approximation to be valid at the endpoints t r ± ∆t.Given that the stationary phase approximation is valid, we can write Then since we are near resonance Ω(t r ) ∼ m DM H and we find that Eq. 4.4 requires 3. Finally for computational ease, in the stationary phase integrals, we want to be able to take the limit ∆t → ∞.More precisely this will require that ∆t is much larger than the spread of the Gaussian integrand, Ω± ∼ √ m DM H.We are then able to take the ∆t → ∞ limit as long as, It is easy to see that all 3 conditions are satisfied if Using m DM ≥ 10 −22 eV and H ≤ 10 −29 eV, this translates to 1 ≫ H(t r )∆t ≫ 10 −4 which is easily satisfied.This shows that we are able to choose endpoints that satisfy all of the 3 conditions.Then from the first and second condition, as shown in Fig. 8, the total L ± can be broken into a series of alternating resonant and fast contributions.However, as we will see, the resonant pieces will always dominate over the fast contributions.We will find that these terms take the form Where L f ast ± (t 0 ) and L res ± (t r ) do not depend on the dark matter phase β.This allows us to square and average over the phase which simply eliminates any terms that get a non-trivial phase.This eliminates not only cross terms between the L + and L − pieces, but any cross terms between different resonances.The end result is (C.9) Note that L + and L − have different sets of resonant times.Finally, we must compute the fast and resonant integrals to find |L f ast ± | and |L res ± |.

C.1 Fast limit
In the limit of fast oscillations we assume Ω ± ≫ H for the entire integral.Then we can rewrite Eq.C.1 as This can then be integrated by parts Since the time derivative in the second term is hitting quantites that change on the Hubble scale, this second term represents an O(H/Ω ± ) correction and can be ignored.In the first term, the piece evaluated at t can be ignored due to the scale factor in the denominator.Finally, since β ≪ m DM we can ignore this term in the denominator, leaving We see we get the exact β dependence predicted in Eq.C.8.We find then

C.2 Resonant Limit
Now we look at the resonant integral.
Expanding the integrand to leading order about the resonance time t r gives This is simply a gaussian integral and can be easily computed in the limit ∆t → ±∞.

D Pre-recombination Distortion Computation
Here we give some of the details of the pre-recombination distortion.First, we describe how to use the probabilities given in Eq. 3.15 to compute the rate of photon conversion in Eq. 4.11.
Second we describe the Green's function method for using this rate to compute the distortions in different eras.

D.1 Photon Conversion rate
Here we will give the details of computing the conversion rate for photons into particle X given in Eq. 4.10.Because the photons are not free streaming, our integral is over a small time interval, τ , with respect to the expansion rate H. Therefore, quantities which depend on time through the expansion of the universe are approximately constant.We then can write, P (k, t 0 , t 0 + τ ) = ρ 0 DM v X (a(t 0 )) 8f This can be easily integrated.
ℓ ± = 2e iΩ ± (t 0 +τ /2) sin(Ω ± τ /2) Now we will square this and average over the dark matter phase β as discussed in Appendix B.3.From the Ω ± in the exponential we get a factor of e ∓iβ(t 0 +τ /2) .Everywhere else we can ignore β because it is sub-leading.Then just as for the free streaming distortion, the phase averaging eliminates the cross terms between ℓ + and ℓ − .This gives, Finally, one can plug this into Eq.4.10 and evaluate that integral analytically to get Eq.4.11.

D.2 Green's Function Method
Here we present the Green's function for the different eras used in Eq. 4.7 and explain how to compute the parameters μ, μt , y t and y in Eq 2.2.The Green's functions used are taken from Ref. [49] for the µ and y era and from Ref. [51] for the µ − y transition era and modified by absorbing and moving a few factors to fit with the definition in Eq 4.7.To start let us define the temperature shift function T (x), the µ distortion shape M (x), and the y distortion shape Y (x) x 3 e x −1 is the unit-normalized blackbody energy spectrum.J * (a) is called the visibility function and captures how inefficient bremsstrahlung and double Compton scattering are at changing the number density.J * (a) goes to 0 for a ≪ 10 −6 and quickly goes to zero at early times.Its analytic form can be found in Eq. 13 of Ref. [49].Note that the M (x) factors out completely and we can write the distortion as where the a integral runs from 0 to a = 3.3 • 10 −6 .

µ-y Transition Era
The Green's function for the transition era is similar to that of the µ era with the addition of a y-distortion piece.The additional factors J µ and J y smoothly transition the Green's function from having mostly µ distortion at early times in the era to mostly y-distortions late in the era.Their analytical form can be found in Eq. 5 of Ref. [51].Much like in the µ era we can write this distortion as Here the a integral runs from a = 2 • 10 −5 to the time of recombination at a * .

Figure 1 .
Figure 1.A timeline of the types of relevant distortion eras.The timeline is presented with time described by decreasing temperature T .

Figure 2 .
Figure 2. A sample diagram for each of the processes holding photons in equilibrium with the electrons.

Figure 3 .
Figure 3. Diagrams for the probability of a photon produced at time t 0 to have converted to a dark photon or axion by a time t.The vertex indicates the interaction with the background dark matter field.λ and λ ′ represent the polarizations of the relevant particles.

Figure 4 .
Figure 4. Bounds on the coupling 1/f a (GeV −1 ) plotted as a contour plot as a function of the dark photon mass (m D ) and the axion dark matter mass (m a ).The plot on the left shows the bounds for the axion dark matter case and the plot of the right shows the dark photon dark matter case.The grey region represents the region where the previous best bound on the coupling derived from red giant cooling [37] is stronger.

Figure 5
Figure 5. Bounds on the coupling 1/f a as a function of the axion dark matter mass m DM plotted for various values of the dark photon mass, m X .The purple, blue, green and orange dashed lines are the individual µ, y, Doppler and free streaming bounds respectively while the solid black line is the total bound from all distortions.The grey line represents the current best bound on the coupling derived from red giant cooling[37].The bounds shown here are for the axion dark matter case, but the bounds for the dark photon dark matter case can be obtained by scaling the bounds up by a factor of √ 3. The bounds for m X < 10 −13 eV are exactly those given in the upper-left plot.

Figure 8 .
Figure 8. Breaking the time interval between t and t 0 into subintervals that either contain or do not contain a resonance.Note that this is easily extended to the case of multiple resonances.The size of the intermediate integral is greatly exaggerated so that it is visible.