Hints of dark photon dark matter from observations and hydrodynamical simulations of the low-redshift Lyman- α forest

Recent work has suggested that an additional < ∼ 6 . 9 eV per baryon of heating in the intergalactic medium is needed to reconcile hydrodynamical simulations with Lyman- α forest absorption line widths at redshift z (cid:39) 0 . 1. Resonant conversion of dark photon dark matter into low frequency photons is a viable source of such heating. We perform the ﬁrst hydrodynamical simulations including dark photon heating and show that dark photons with mass m A (cid:48) ∼ 8 × 10 − 14 eV c − 2 and kinetic mixing (cid:15) ∼ 5 × 10 − 15 can alleviate the heating excess. A prediction of this model is a non-standard thermal history for underdense gas at z > ∼ 3. In this Supplemental Material, we show results from a wider range of our hydrodynamical simulations that use diﬀerent values of the A (cid:48) mass and kinetic mixing. In contrast to Fig. 2 in the main text, rather than show the best ﬁt model, we show results from the individual simulations that span a range of A (cid:48) masses and kinetic mixings. As already discussed in the main text, this demonstrates the eﬀectiveness of the low-redshift Ly- α forest data for obtaining best ﬁt parameters for A (cid:48) , rather than just individuating a preferred region of the parameter space.

Introduction.-The Lyman-α forest, a series of absorption features that arise from the distribution of intergalactic gas, is a powerful tool for investigating the properties of dark matter (DM). The absorption is typically observed in the spectra of distant, luminous quasars, where resonant scattering along the line of sight occurs as photons redshift into the rest frame Ly-α (n = 1 → 2) transition of intervening neutral hydrogen [1]. The 1D power spectrum of the Ly-α forest transmitted flux is an excellent tracer of the underlying DM distribution on scales ∼ 0.5-50 comoving Mpc, and has been routinely used to place tight constraints on warm dark matter [2][3][4][5][6], fuzzy DM [7][8][9][10], as well as primordial black hole (PBH) DM [11].
The physics behind Ly-α forest temperature measurements is straightforward. The hydrogen atoms in the IGM will in general not be at rest, but will undergo thermal motion described by a Maxwellian velocity distribution, leading to a line width ∆ν = ν α (b 2 th + b 2 nth ) 1/2 /c, where ν α is the resonance line frequency, b th = (2k B T /m H ) 1/2 is the Doppler parameter due to thermal motion, m H is the hydrogen atom mass and k B is Boltzmann's constant. Here, b nth accounts for any additional, non-thermal line broadening, including, smoothing of the absorbing structures by gas pressure, smallscale turbulence, peculiar motion or expansion with the Hubble flow. b nth is in general non-zero for all but the narrowest Ly-α lines. Hence, given a forward model for b nth , typically provided by cosmological hydrodynamical simulations, a determination of the Ly-α spectral widtheither by directly measuring Doppler parameters or using another statistic that is sensitive to the thermal broadening kernel -allows a measurement of the IGM temperature [33][34][35][36][37][38][39][40][41].
Thus far, all the calorimetric constraints on new physics from the Ly-α forest have relied on observations at redshifts z > ∼ 2. In this Letter, we pioneer the use of low-redshift (z 0.1) Ly-α forest observations for this purpose. Three independent studies [42][43][44] have now highlighted a discrepancy between the widths of Lyα forest absorption lines measured from Hubble Space Telescope/Cosmic Origins Spectrograph (COS) data at z 0.1 [45,46] and the predictions from detailed cosmological hydrodynamical simulations. The simulated line widths are always too narrow compared to the observations. Appealing to enhanced photoelectric heating alone is not a viable solution, as the integrated UV background would require an unphysically hard spectrum [47]. This implies there is a non-canonical heating process in the IGM neglected in the simulations, such that an additional < ∼ 6.9 eV per baryon is deposited into typical Lyα forest absorbers by z = 0.1, and/or the non-thermal line broadening has been underestimated [47]. Additional turbulence below the simulation resolution limit is possible [44,48,49], although current models have line widths that become narrower, rather than broader, with increasing resolution. Here, we instead explore the role of additional heating, and suggest that ultralight dark photon DM with a small mixing with the Standard Model (SM) photon provides an intriguing solution to the line width discrepancy. Dark photons can undergo resonant conversions into SM photons, which are subsequently absorbed by the IGM. The condition for resonant conversion is set by the dark photon mass and the local electron density, allowing dark photons to naturally explain the low-redshift Ly-α forest data without altering the established agreement with IGM temperature measurements at z > 2.
Dark photon dark matter.-The dark photon A is a minimal and well-motivated extension of the Standard Model (SM) which kinetically mixes with the ordinary photon, γ [50][51][52][53][54][55][56][57]. Ultralight dark photons are also an attractive cold DM candidate, with several earlyuniverse production mechanisms that have been studied extensively in the literature that are capable of producing A non-relativistically [58][59][60][61][62][63][64][65][66][67][68][69]; its perturbations are therefore expected to be well-described by the standard ΛCDM matter power spectrum. The photon-dark photon Lagrangian reads where is the dimensionless kinetic mixing parameter, and m A is the A mass. F and F are the field strength tensors for γ and A respectively. The presence of kinetic mixing, and the resulting mismatch between the interaction and propagation eigenstates, induce oscillations of dark photons into photons, A → γ, and vice-versa. In the presence of a plasma, ordinary photons acquire an effective mass, m γ , given primarily by the plasma frequency of the medium [31,70,71]. At every point in space x and redshift z, the effective plasma mass is given by where n e is the free-electron number density. At points in space where m 2 γ (z, x) = m 2 A , the probability of conversion is resonantly enhanced [30,31,70,71]. This process can be understood as a two-level quantum system with an energy difference that is initially well-separated, but with one state having its energy evolve with time. Whenever the energy difference passes through zero, a nonadiabatic transition between the two states occurs, with transition probability described by the Landau-Zener formula [31,70,72,73].
If A constitutes the DM, the probability of conversion of A into photons is [29][30][31] where t res is the time at which the resonant condition is met. Here, h P ≡ 2πh is Planck's constant. For m A between 10 −15 -10 −12 eV c −2 , A DM converts into low-frequency photons with frequency ν = m A c 2 /h P , which rapidly undergo free-free absorption in the ionized IGM. The mean free path, λ ff , is given approximately by [74,75] where α 1/137 is the fine-structure constant, σ T is the Thomson cross section, T is the temperature of the IGM, and g ff (ν, T ) is the Gaunt factor, which only has a mild dependence on ν and T . Numerically, taking g ff = 15.7 [76], we find the local overdensity of baryons at which A conversion occurs. We assume that overdensities in both free electrons and baryons are equal, which is true for the IGM after reionization is complete. Since λ ff is much smaller than the typical size of a Ly-α forest absorber (∼ 100 kpc), we can safely assume that the absorption of these photons occurs locally in our simulations.
Since A → γ conversions only occur when the resonance condition is met, at each point in time, heating only takes place in regions of specific ∆ b , such that m 2 γ (z, x) = m 2 A is satisfied. Moreover, the probability of conversion is governed by the rate at which ∆ b is evolving in each of those regions. In regions where conversions are happening, the energy deposited per baryon E A →γ due to such a resonance is simply where n b is the local number density of baryons, and ρ A is the mass density of A DM.
We can estimate E A →γ by assuming that the number density of baryons everywhere evolves only through adiabatic expansion, i.e. n e ∝ (1 + z) 3 is the Hubble parameter, and the probability of conversion simplifies to P A →γ = π 2 m A c 2 /(3H(z res )h), with z res indicating the redshift at which the resonance condition is met. For m A ∼ 8 × 10 −14 eV c −2 , we expect baryons at the mean cosmological density to experience resonant conversion at z ∼ 2. Further assuming matter domination and ∆ dm = ∆ b , we obtain approximately where −14 ≡ /10 −14 , and m −13 ≡ m A /(10 −13 eV c −2 ). We see that even a mixing as small as −14 ∼ 1 leads to the absorption of enough photons to heat the IGM by several eV per baryon [29][30][31][32].
The resonant nature of A heating makes it an attractive model for reconciling low and high redshift Ly-α forest data. The smaller the value of m A , the later the resonance condition is met at fixed ∆ b . For a sufficiently light A , most of the resonant conversions occur at z < ∼ 2, broadening the absorption line widths at low redshifts, without depositing a significant amount of heat at z > 2. Furthermore, the resonance condition is set by during matter domination, giving approximately the density dependence required by observations [47].
Simulations.-To fully test the viability of heating from A dark matter, we now turn to implementing this model in Ly-α forest simulations for the first time. Cosmological hydrodynamical simulations of the Ly-α forest were performed with a version of P-Gadget-3 [77], modified for the Sherwood simulation project [14]. Following [47], we use a simulation box size of L = 60h −1 cMpc with 2 × 768 3 gas and dark matter particles, giving a gas (dark matter) particle mass of M gas = 6.38 The simulations were started at z = 99, with initial conditions generated on a regular grid using a ΛCDM transfer function generated by CAMB [78]. The cosmological parameters we use are Ω m = 0.308, Ω Λ = 0.692, h = 0.678, Ω b = 0.0482, σ 8 = 0.829 and n = 0.961 [79], with a primordial helium fraction by mass of Y p = 0.24. All gas particles with overdensity ∆ b > 10 3 and temperature T < 10 5 K are converted into collisionless star particles, and photoionization and heating by a spatially uniform UV background is included [13]. Mock Ly-α forest spectra were extracted from the simulations and processed to resemble COS observational data following the approach described by [47], where further details and tests of our numerical methodology can be found.
In contrast to [47], however, in this work we also implement dark photon heating in our hydrodynamical simulations using Eq. (3). It will be numerically convenient to assume the baryons closely trace the dark matter in the IGM, and indeed, on scales exceeding ∼ 100 kpc (the pressure smoothing scale in the IGM), this is a good approximation. Defining the overdensity of a given cosmological species as ∆ i = ρ i / ρ i , we thus assume ∆ dm = z) is determined for each gas particle at position x and redshift z in our simulations. For each gas particle, we then set For a fixed value of m A , we may therefore determine where and when a resonant conversion happens for each gas particle (i.e. when the condition m 2 γ ( x, z res ( x)) = m 2 is met). The converted energy per baryon E A →γ is then calculated using Eq. (3), and directly injected into the gas particles.
Results.-We first demonstrate the effect that dark photons have on the IGM temperature by considering the redshift evolution of gas parcels at fixed overdensity, ∆ b , heated by both the UV background and E A →γ . We adapt the non-equilibrium ionization and heating calculations performed by [47] for this purpose. In the upper panels of Fig. 1 we show the thermal history of gas at the mean density, ∆ b = 1. The solid gray curves correspond to UV heating only (i.e. no A ), following the synthesis model presented in [13]. The data points are IGM temperature measurements at the mean density derived from the Ly-α forest [38,41]. All other curves in Fig. 1 include A heating. These exhibit a sharp rise in the gas temperature when the resonance condition is met. In the top left panel, we have fixed the A mass to m −13 = 0.8 and varied the kinetic mixing parameter, . In this case, the timing of the energy injection does not change, but the amplitude of the temperature peak increases with .
The top right panel instead shows the results for a fixed Contours show ∆χ 2 = χ 2 − χ 2 min = 1 and 4, corresponding to the projection of the 68% and 95% intervals for the individual parameters m A and . The dashed red curves show the results for a tight prior on the z = 2 IGM temperature at mean density from [41]. The solid blue curves show the effect of a weaker prior, where the 1σ uncertainty from [41] has been increased by a factor of four for consistency with the independent temperature measurement from [38]. Lower panels: The corresponding best-fit models compared to the COS observational data. The solid gray curve shows the UV heating model with no dark photon heating [13]. kinetic mixing, −14 = 0.5, for different A masses. In this case, smaller A masses result in later injection of heat into the gas parcel. From this, we may already conclude that A masses with m −13 > ∼ 0.9 and −14 = 0.5 are excluded by the data from [38,41] at 2 < z < 4. This is consistent with the bounds derived analytically in [30] using earlier measurements of the mean density IGM temperature (see their Fig. 9, as well as [29,32]).
The lower panel of Fig. 1 shows the gas thermal evolution for a fixed pair of parameters [m −13 , −14 ] = [0.8, 0.5], but now for varying gas overdensities. For the adopted value of the A mass, m −13 = 0.8, the resonance at the mean background density occurs at z res = 1.8 (black solid curve). Later energy injection occurs for overdensities (log 10 ∆ b > 0), while earlier energy injection occurs for underdensities (log 10 ∆ b < 0). This density dependence allows low-redshift Ly-α forest observa-tions to place a bound on A heating that is complementary to constraints obtained at z > 2. The Ly-α forest at z = 0.1 is sensitive to gas at overdensities of ∆ b 10 [47], whereas at z > 2, it probes gas close to the mean density [36]. Hence, if appealing to A heating as a possible resolution to the COS line width discrepancy, we require m −13 > 0.6 for resonant conversion to occur in gas with ∆ b = 10 by z = 0.1; smaller masses would inject the energy too late. Coupled with the upper bound from data at 2 < z < 4, this implies 0.6 < ∼ m −13 < ∼ 0.9 for all kinetic mixing parameters that heat the IGM by more than a few thousand degrees. Note also that heating of the underdense IGM is expected at z > 2, even if the mean density IGM temperature constraints at 2 < z < 4 are fully satisfied. Intriguingly, there are indeed hints from the distribution of the Ly-α forest transmitted flux that underdense gas in the IGM at z = 3 is hotter than expected in canonical UV photo-heating models [80,81]. We intend to explore this further in future work.
We now turn to obtaining a best-fit dark photon dark matter model from the low-redshift Ly-α forest assuming maximal A heating, i.e. with no other sources of noncanonical heating. Using the recipe outlined in the pre- Voigt profiles were fit to the mock Ly-α spectra, giving the Ly-α line column densities, N HI , and Doppler parameters, b. The b-distribution and column density distribution function (CDDF) were then constructed for each simulation following [47], and the model grid was used to perform a χ 2 minimisation on the covariance matrix derived from the COS observations.
The resulting best fit A parameters are presented in Fig. 2. The ∆χ 2 = χ 2 − χ 2 min contours corresponding to the 68% and 95% intervals for individual parameters are displayed in the upper panel, while the lower panels show the best-fit models (corresponding to the crosses in the upper panel). We assume two different priors on the temperature of the IGM at z = 2: a tight prior of T 0,z=2 = 9500 ± 1393 K based on [41] (red dashed contours), and a weaker prior where the 1σ uncertainty on T 0,z=2 has been increased by a factor of four for consistency with the independent measurement of T 0,z=2 = 13721 +1694 −2152 K from [38] (blue solid contours). As expected from Fig. 1, a weaker prior has the effect of increasing the best-fit A mass. For the weak (tight) z = 2 prior, the best-fit model has m −13 = 0.84±0.06 (m −13 = 0.80 ± 0.04) and −14 = 0.46 +0.05 −0.04 ( −14 = 0.46 +0.06 −0.05 ) (1σ) for χ 2 min /ν = 1.13 (χ 2 min /ν = 1.37) and ν = 16 degrees of freedom, with a p-value of p = 0.32 (p = 0.14). For comparison, a model with UV photon heating only [13], shown by the gray curves in the lower panels of Fig. 2, has χ 2 min /ν = 3.50 and p = 2.5 × 10 −6 assuming the weak z = 2 thermal prior. The addition of A heating considerably improves the fit and leads to very good agreement with the COS data. For the weak z = 2 prior, the bestfit A parameters deposit an extra 5.3 eV per baryon into gas with ∆ b = 10 by z = 0.1, consistent with the limit of < ∼ 6.9 eV per baryon obtained by [47].
Conclusions.-In this Letter we have pioneered the use of the Ly-α forest at redshift z 0.1 as a calorimeter for investigating properties of the dark sector. Specifically, we studied a model of ultralight dark photons, A , that can naturally alleviate the tension [42][43][44]47] between the (too narrow) Ly-α absorption line widths predicted in hydrodynamical simulations compared to observational data at z = 0.1. Assuming a maximal contribution from A heating and a thermal prior of T 0 = 9500 ± 5572 K at z = 2 [38,41], our best-fit model has A mass m A = 8.4 ± 0.6 × 10 −14 eV c −2 and kinetic mixing parameter = 4.6 +0.5 −0.4 × 10 −15 (1σ). Although astrophysical sources, such as turbulent broadening or other non-canonical heating processes may also explain the line width discrepancy, our study is a first clear indication that DM energy injection can be a compelling alternative. We also highlight that our best-fit A parameters will have testable consequences for the temperature of the underdense IGM at z = 3, where there are already hints of missing heating in Ly-α forest simulations [80,81].
Finally, we note that dark photons in our mass range of interest can be produced around spinning black holes (BHs) through the superradiance instability. The superradiant cloud can affect the spin-mass distribution of BHs or produce gravitational waves, giving a promising way to look for A with m A ∼ 8 × 10 −14 eV c −2 . At the moment, BH measurements appear to be in tension with this A mass [82][83][84], but are currently subject to significant uncertainties [84,85] and model dependence [82,[86][87][88]. Sharpening our understanding of the superradiance phenomenon, as well as improving the experimental searches, will be pivotal in testing our model.
Acknowledgments In this Supplemental Material, we show results from a wider range of our hydrodynamical simulations that use different values of the A mass and kinetic mixing. In contrast to Fig. 2 in the main text, rather than show the best fit model, we show results from the individual simulations that span a range of A masses and kinetic mixings. As already discussed in the main text, this demonstrates the effectiveness of the low-redshift Ly-α forest data for obtaining best fit parameters for A , rather than just individuating a preferred region of the parameter space.  [42]. Solid gray curves show a simulation with UV photon heating only [13]. Note that the dark photon heating impacts both the Ly-α line widths and the shape of the CDDF, making a joint fit particularly powerful. Fig. S1 shows the effects of dark photon energy injection on the z = 0.1 Ly-α forest b-parameter distribution (left panel) and column density distribution function (CDDF, right panel) predicted by the cosmological hydrodynamical simulations. Note that, following [42,47], the neutral hydrogen densities in the models have been linearly rescaled to approximately match the amplitude of the CDDF at log 10 (N HI /cm −2 ) ∼ 13.5 (corresponding to a baryon overdensity ∆ b ∼ 10). The only difference here is we consider Doppler parameters in the range 20 km s −1 ≤ b ≤ 80 km s −1 (cf. 20 km s −1 ≤ b ≤ 90 km s −1 in [47]); the largest Doppler parameters are highly suprathermal and are insensitive to A heating.
From the b-parameter distribution in the left panel of Fig. S1, it is clear that A with mass m −13 = 0.6 does not give substantially better agreement with the data compared to the case without A heating (gray line); both models overshoot the data for b < 25 km s −1 . In this particular case, the A heating occurs too late to affect the typical gas densities probed by the Ly-α forest at z = 0.1 (for m −13 = 0.6 the resonance redshift is z res = 0.06 for ∆ b = 10). Notice also that some of the A parameters excluded by the b-parameter distribution provide a reasonable fit to CDDF, while, on the other hand, some parameters that yield a reasonable fit to b-parameter distribution data are inconsistent with the CDDF data (see for example the orange dot-dashed curve in Fig. S1). In the particular case of m −13 = 0.8 (corresponding to z res = 0.28 for ∆ b = 10), gas with log 10 (N HI /cm −2 ) ∼ 13.5 (∆ b ∼ 10) is hotter in the models with larger kinetic mixing. However, gas at ∆ b > ∼ 16 (i.e. densities where z res < 0.1 for m −13 = 0.8) corresponding to the higher column density Ly-α absorbers) is not heated because the resonance threshold has not been crossed. This has the effect of flattening the gradient of the CDDF; if the lower density gas is hotter, it is also more ionized (the case-A recombination rate for hydrogen, α A ∝ T −0.72 , where T is the gas temperature). This results in relatively more of the stronger lines with log 10 (N HI /cm −2 ) > 13.8 as the kinetic mixing is increased. This demonstrates the interplay between these two observables which, when combined, provide a more powerful test of A heating than either measurement alone.