Constraining early dark energy with gravitational waves before recombination
CConstraining early dark energy with gravitational waves before recombination
Zachary J. Weiner, ∗ Peter Adshead, † and John T. Giblin, Jr.
2, 3, ‡ Illinois Center for Advanced Studies of the Universe & Department of Physics,University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA Department of Physics, Kenyon College, Gambier, Ohio 43022, USA CERCA/ISO, Department of Physics, Case Western Reserve University, Cleveland, Ohio 44106, USA
We show that the nonperturbative decay of ultralight scalars into Abelian gauge bosons, re-cently proposed as a possible solution to the Hubble tension, produces a stochastic background ofgravitational waves which is constrained by the cosmic microwave background. We simulate thefull nonlinear dynamics of resonant dark photon production and the associated gravitational waveproduction, finding the signals to exceed constraints for the entire parameter space we consider. Ourfindings suggest that gravitational wave production from the decay of early dark energy may providea unique probe of these models.
Measurements of the current expansion rate H asinferred from the acoustic peaks in the cosmic microwavebackground radiation (CMB) are in tension with the valueobtained from local measurements [1–3]. Rather thanalter the expansion history at intermediate redshifts, newphysics introduced to resolve the tension must alter theabsolute scales of the low- or high-redshift anchors of thecosmic distance scale [4, 5]. Early-Universe resolutionsthus focus on altering the high-redshift anchor, the CMBsound horizon at recombination. For a recent review, seeRef. [6].In particular, an increased expansion rate before recom-bination decreases the sound horizon. However, simplyadding more radiation [7, 8] also changes the dampingscale in a way that is increasingly disfavored by high-precision measurements of the high-multipole dampingtail [2, 9, 10]. To avoid this, one class of proposed solu-tions, so-called early dark energy (EDE) models, postu-lates an additional energy component that is only tran-siently important near recombination [11–18]. These pro-posed solutions supposedly relieve the tension betweenthe early and late datasets; however, see Ref. [19].In the simplest EDE implementations, a scalar field isinitially frozen up its potential in a homogeneous configu-ration. The field’s mass is tuned such that it begins toevolve near matter-radiation equality, oscillating aboutthe minimum of its potential. In order to redshift awayfast enough (at least as fast as radiation), the potentialmust have no quadratic term about its minima. Froma particle physics perspective, this requires an explana-tion; to avoid such extreme fine-tuning, Ref. [18] insteadproposes a model of a decaying ultralight scalar (dULS).Instead of oscillating about a peculiar potential and red-shifting away, the EDE scalar field decays resonantlyto dark radiation during oscillations about a quadraticminimum.Nonperturbative or resonant particle production is a ∗ [email protected] † [email protected] ‡ [email protected] common feature of early-Universe preheating after in-flation (for reviews, see [20–22]). Substantial study hasestablished that these violent processes can also lead to thecopious production of gravitational waves [23–38]. Thewavelength of the produced gravitational waves (GWs) isdetermined by the characteristic scale of particle produc-tion, which must be smaller than the horizon size. Duringpreheating after inflation this restricts the production tofrequencies from MHz to GHz, well beyond the reach ofcurrent or planned detectors [39–45]. By contrast, in order to resolve the Hubble tension,particle production due to the decay of the ultralightscalar must occur near the time of matter-radiation equal-ity. Since resonant particle production occurs at scalesnear the horizon scale at that time, gravitational waveemission occurs at current-day frequencies near 10 − Hz.CMB anisotropies constrain stochastic gravitational waveswith peak sensitivity at present-day frequencies near10 − Hz [48–51]. In this paper, we confront models ofultralight scalar decay into dark photon with these con-straints.
Background. —Following Ref. [18], we study the reso-nant decay of an ultralight axion φ into a (dark) Abeliangauge field A µ described by the action S = (cid:90) d x √− g (cid:34) M R − ∂ µ φ∂ µ φ − V ( φ ) − F µν F µν − α f φF µν ˜ F µν (cid:35) . (1)Here f is the axion decay constant and α is a dimension-less coupling that parameterizes the rate and efficiencyof energy transfer. The field strength tensor of the darkphoton is F µν ≡ ∂ µ A ν − ∂ ν A µ , while its dual is ˜ F µν = (cid:15) µναβ F αβ /
2, where (cid:15) µνρσ is the Levi-Civit`a symbol with However, high-frequency gravitational waves contribute to theenergy budget of the Universe as radiation. Constraints on thegravitational-wave energy density from N eff [46, 47] can be usedto indirectly restrict preheating scenarios [36–38]. a r X i v : . [ a s t r o - ph . C O ] A ug convention (cid:15) = 1 / √− g . Following Ref. [18], we takethe standard axion potential V ( φ ) = m φ f (1 − cos φ/f ).We set c = (cid:126) = k B = 1 and use M pl = 1 / √ πG N to de-note the reduced Planck mass. We work with the “mostly-plus,” conformal Friedmann-Lemaˆıtre-Robertson-Walker(FLRW) metric, g µν = a η µν = a diag ( − , , , τ .The classical equations of motion, A (cid:48)(cid:48) ν = ∂ i ∂ i A ν + η βν αf ∂ α φ (cid:18) √− g(cid:15) αβρσ F ρσ (cid:19) (2) φ (cid:48)(cid:48) = ∂ i ∂ i φ − H φ (cid:48) − a d V d φ − a α f F µν ˜ F µν , (3)permit solutions in which fluctuations of the gauge fieldsare exponentially enhanced via a tachyonic instabilitysourced by a homogeneous, rolling axion. Initially, a cos-mological axion has some static homogeneous component (cid:104) φ (cid:105) = θf , expressed in terms of the initial misalignmentangle θ . As a result, the axion’s energy is dominatedby its potential, acting as a source of early dark energywhich could alleviate the Hubble tension [13]. On theother hand, to linear order the helical polarizations of A i obey [52] A ±(cid:48)(cid:48) ( k ) + k (cid:18) k ∓ αf (cid:104) φ (cid:48) (cid:105) (cid:19) A ± ( k ) = 0 . (4)Thus, once the axion begins to oscillate (when the Hubbleparameter drops below ∼ m φ / k < α/f × (cid:104) φ (cid:48) (cid:105) .As the axion crosses zero, (cid:104) φ (cid:48) (cid:105) changes sign and so am-plifies the other polarization. Eventually, the gauge fieldfluctuations become so large that nonlinear effects beginto fragment the axion background, ending the phase oftachyonic resonance. Both the initial exponential gaugefield production and subsequent nonlinear dynamics cansource a significant gravitational wave background. Grav-itational waves correspond to the tensor part of perturba-tions to the spatial part of the spacetime metric, h (cid:48)(cid:48) ij − ∂ k ∂ k h ij + 2 H h (cid:48) ij = 2 M T TT ij , (5)where T TT ij is the transverse and traceless part of thestress-energy tensor T ij .We employ numerical simulations in order to fully cap-ture resonance, the nonlinear dynamics which terminateenergy transfer, and the resulting production of gravita-tional waves. We solve the classical equations of motionEqs. (2), (3), and (5) in a homogeneous ΛCDM cosmology,self-consistently including the contribution of the dULSsector to the expansion rate. We discretize these equationsonto a three-dimensional, periodic, regularly-spaced grid,computing spatial derivatives via fourth-order centereddifferencing and utilizing a fourth-order Runge-Kutta − − − − − ρ A / m φ f α = 50 α = 55 α = 60 α = 65 α = 70 z c z eq z − Ω d U L S FIG. 1. The total energy in the gauge fields (top) and thefractional energy in the dULS sector (bottom) as a functionof redshift for various couplings α indicated in the legend. Allsimulations fix m φ = 10 − eV, f = 1 . × GeV, and θ = 2.Dashed and solid vertical lines respectively indicate matter-radiation equality and the favored redshift of the transitionfrom dark energy to matter from the analysis of Ref. [18]. method for time integration. All results presented usegrids with N = 768 points, side-length L = 10 /m φ ,and a timestep ∆ τ = ∆ x/
10 = L/ N . We implementsimulations using pystella [37, 38] and provide details onour algorithm, initial conditions, and convergence testsin the Appendices. Results .—In our simulations of the decaying ultralightscalar model, we consider benchmark scenarios fromRef. [18], taking m φ = 10 − − − eV so that thedULS sector transitions from dark energy to matter-likebehavior around the favored redshift z c ≈ f = 1 . × GeV so that the dULS sector makes up ∼ −
4% of the Universe’s energy budget at its peak. Wefix an initial misalignment angle θ = 2 for convenience.We comment later on the dependence of the gravitationalwave signal to slight changes in these choices.We first verify that the dynamics of the dULS sectorqualitatively reproduce the effective fluid description em-ployed in Ref. [18] (which we evaluate in more detail inthe Appendices). In Fig. 1 we display the energy in thegauge fields ρ A and the fractional energy in the dULSsector, Ω dULS = ( ρ A + ρ φ ) /ρ , as a function of redshift. Wevary the coupling α from 50 to 70, spanning values large − − − Ω g w , h α − − − f (Hz)10 − − − Ω g w , h m φ × − eV4 × − eV8 × − eV FIG. 2. The present-day gravitational wave spectrum emittedby recombination (i.e., evaluated at z ≈ m φ =10 − eV and varying α (top) and fixing α = 70 and varying m φ (bottom). All simulations set f = 1 . × GeV and θ = 2. In black is the upper bound as constrained by theCMB, computed by Ref. [50]. enough for resonance and GW production to terminatebefore recombination while small enough to be reliablyresolved by our grid.Figure 2 depicts the gravitational wave signal resultingfrom the resonant production of dark photons, evaluatedat z = 1100. We compare the signal to the constraintsfrom a recent analysis of current data in Ref. [50]. Whilethe peak of the fractional gravitational wave energy spec-trum, reaching above 10 − , resides at larger frequencies ∼ − Hz, the infrared tail of the spectrum exceedsconstraints at frequencies (cid:46) − Hz by roughly an orderof magnitude. We note that, though we are unable tosimulate larger volumes to capture even lower frequencies,we expect the signal should continue roughly as a powerlaw further into the infrared, Ω gw ∝ k/k (cid:63) , where k (cid:63) is Note that the constraints of Ref. [50] are not directly applicableto this scenario, since they are computed from adiabatic initialconditions—constant initial GW amplitude on superhorizon scales.By contrast, the GWs here are actively sourced, analogous tothose in defect scenarios, e.g., Ref. [54]. However, we expectconstraints on active modes to be competitive, if not more severethan those on adiabatic modes. the peak wavenumber, as suggested by recent analyticestimates of GW production in similar scenarios [55].The analysis of Ref. [18] determined that the dULSsector should begin to decay like radiation by a redshiftbetween ∼ z ≈ − m φ is relatively simple. Fromthe transfer function Eq. (A10), the present-day frequencyscales as f ∼ k/ (cid:112) HM pl ∼ (cid:112) m φ /M pl . While the axionmass scales out of the dynamics, it has an effect on theinitial amplitude of gauge field vacuum fluctuations. Alower mass sets initial fluctuations with lower amplitude,requiring a longer period of resonance to fully deplete theaxion’s energy. However, this effect is relatively unimpor-tant and is easily compensated for by a slight increasein the coupling α . At smaller masses, resonance beginslater; therefore, marginally larger couplings are requiredin order for the process to complete before recombination(so that the signals are detectable). As can be seen inthe bottom panel of Fig. 2, this condition is easily metwith α = 70 and a variety of masses between 10 − and10 − eV. Furthermore, the shape and amplitude of theGW signal itself is qualitatively independent of the axionmass m φ .The value of the axion decay constant f is set bythe requirement that the dULS sector comprises a frac-tion of the Universe’s energy between 3 − ρ φ /ρ ∼ m φ ( θf ) /H ∼ ( θf ) . In turn, the amplitudeof the resulting gravitational wave signal is directly pro-portional to the square of the fraction of the Universe’senergy residing in its source [56]. Since the signals we findhere exceed current constraints by an order of magnitude,we do not expect the constraining power of the GW signalto be sensitive to any uncertainty in the best-fit ρ φ /ρ .Since θ sets the amplitude of axion oscillations (andso (cid:104) φ (cid:48) (cid:105) ), its effect on the dynamics is, to linear order,degenerate with the coupling. However, nonlinear effects(e.g., rescattering of power to higher momenta) becomemore important with larger couplings α , and so our choiceof θ = 2 allows reliable simulations with smaller couplingsthat still transition the dULS sector to a radiation-likestate on the required timescales.As a final investigation, we study whether the gravi-tational wave signal is significantly polarized. The sameaxial coupling of gauge fields to the inflaton generates ahelical gravitational wave background during inflation [57–61], and can also imprint on the spectrum of gravitationalwaves produced during preheating [36]. However, in the f (Hz)10 − − − Ω g w , h − − − f (Hz)10 − − − Ω g w , h FIG. 3. The polarization components of the present-daygravitational wave spectrum emitted by recombination (i.e.,evaluated at z ≈ α = 50 (top) and 70 (bottom). Theplus and minus polarizations are in blue and red, respectively,while the total signal is portrayed in dashed black. Bothpanels fix m φ = 10 − eV, f = 1 . × GeV, and θ = 2. Insolid black is the upper bound as constrained by the CMB,computed by Ref. [50]. former case, the sign of the axion’s velocity is fixed dur-ing inflation. Preheating via the axial coupling can alsocomplete within one (or even half an) oscillation of theinflaton [37, 38, 52], but nonlinear effects can result in agravitational wave signal dominated by different helicitiesat different scales. The results presented in Fig. 3 followin spirit. For the lowest coupling we consider, α = 50, theaxion oscillates numerous times before gauge field pro-duction terminates, emitting an essentially unpolarizedgravitational wave background. For the largest coupling, α = 70, the spectrum is moderately polarized at largescales, consistent with a more substantial enhancement ofone polarization before the axion first crosses zero. Whilethe signal at lower frequencies arises predominantly fromthe initial phase of helical tachyonic resonance, higherfrequencies are sourced by nonlinear mode interactionswhich do not retain the same polarization. In summary,it is difficult to evaluate whether the gravitational wavebackground is sufficiently polarized on CMB scales toprovide a unique signature of this model; doing so likelyrequires thorough study via numerical simulations in amodel-dependent way. As an aside, we note that these findings are applicable to models of dark matter as massivedark photons, which are produced via the same resonantinstability considered here [62–64]. Conclusions. —The rapid production of inhomogeneitiesfrom resonant particle production can induce a significantGW background. The amplitude of the GW signal islargest (and so offers the most constraining potential)when the GW source comprises a significant fraction ofthe Universe’s energy budget and occurs close to thehorizon scale at the time of emission [56]. The modelconsidered here, which exhibits a tachyonic instability viaan axion coupled to dark photons, is especially efficient,as has been shown in the context of preheating afterinflation [37, 38] in which case up to the entire energybudget of the Universe may source gravitational waves.New physics which relies on the same mechanism laterin cosmological history is subject to constraints fromdirect probes of stochastic backgrounds of gravitationalwaves [48, 49]. Models of early dark energy proposed toalleviate the Hubble tension are a prime example, as theirsuccess hinges on the new sector making up a substantial( O (1%)) fraction of the Universe’s energy. Furthermore,the source of early dark energy must soon decay one wayor another before recombination, pinning the relevantlength scales to those probed by the cosmic microwavebackground.In this work, we demonstrate that the decaying ultra-light scalar model, as motivated in Ref. [18], produces abackground of gravitational waves with a peak spectralenergy fraction exceeding O (10 − ) at its peak, with apower-law tail extending into the region that is alreadyconstrained by the CMB [50]. While we have not studiedthe entire available parameter space, we show that therequirements for the model to successfully alleviate theHubble tension generally coincide with those for its gravi-tational wave signature to be constrained by the CMB.For the parameter space we considered, we find the signalexceeds constraints (on adiabatic modes) by an order ofmagnitude. Because GWs are actively sourced in thisscenario, these constraints are not directly applicable, andwe leave a detailed computation of the CMB signaturesof these models to future work. We also point out thatresonant particle production is not a unique feature tothis model. The original single-field models of EDE ex-hibit similar parametric instabilities which may also emitsignificant GW backgrounds [15]. More broadly, our find-ings suggest that stochastic backgrounds of gravitationalwaves could provide an orthogonal probe with which toconstrain models of early dark energy.We gratefully acknowledge Tom Clarke, Ed Copeland,and Adam Moss for sharing the data from Ref. [50], andwe thank Tristan Smith for useful discussions. The work ofP.A. was supported in part by NASA Astrophysics TheoryGrant NNX17AG48G. J.T.G. is supported by the NationalScience Foundation Grant No. PHY-1719652. Z.J.W. issupported in part by the United States Department ofEnergy Computational Science Graduate Fellowship, pro-vided under Grant No. DE-FG02-97ER25308. This workused the Extreme Science and Engineering Discovery En-vironment (XSEDE) [65], which is supported by NationalScience Foundation grant number ACI-1548562; simula-tions were run on the Comet cluster at the San Diego Su-percomputer Center through allocation TG-PHY180049.This work made use of the Illinois Campus Cluster, a com-puting resource that is operated by the Illinois CampusCluster Program (ICCP) in conjunction with the Na- tional Center for Supercomputing Applications (NCSA)and which is supported by funds from the University ofIllinois at Urbana-Champaign.Simulations in this work were implemented with pystella ,which is available at github.com/zachjweiner/pystellaand makes use of the Python packages PyOpenCL [66],
Loo.py [67], mpi4py [68, 69], mpi4py-fft [70],
NumPy [71],and
SciPy [72]. [1] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, andD. Scolnic, Astrophys. J. , 85 (2019), arXiv:1903.07603[astro-ph.CO].[2] N. Aghanim et al. (Planck), (2018), arXiv:1807.06209[astro-ph.CO].[3] W. L. Freedman, Nature Astron. , 0121 (2017),arXiv:1706.02739 [astro-ph.CO].[4] J. L. Bernal, L. Verde, and A. G. Riess, JCAP , 019(2016), arXiv:1607.05617 [astro-ph.CO].[5] K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghu-nathan, and W. K. Wu, Astrophys. J. , 4 (2019),arXiv:1811.00537 [astro-ph.CO].[6] L. Knox and M. Millea, Phys. Rev. D , 043533 (2020),arXiv:1908.03663 [astro-ph.CO].[7] M. Wyman, D. H. Rudd, R. Vanderveld, and W. Hu,Phys. Rev. Lett. , 051302 (2014), arXiv:1307.7715[astro-ph.CO].[8] C. Dvorkin, M. Wyman, D. H. Rudd, and W. Hu, Phys.Rev. D , 083503 (2014), arXiv:1403.8049 [astro-ph.CO].[9] W. Hu and M. J. White, Astrophys. J. , 30 (1996),arXiv:astro-ph/9602019.[10] M. Raveri, W. Hu, T. Hoffman, and L.-T. Wang,Phys. Rev. D , 103501 (2017), arXiv:1709.04877 [astro-ph.CO].[11] T. Karwal and M. Kamionkowski, Phys. Rev. D ,103523 (2016), arXiv:1608.01309 [astro-ph.CO].[12] V. Poulin, T. L. Smith, D. Grin, T. Karwal, andM. Kamionkowski, Phys. Rev. D , 083525 (2018),arXiv:1806.10608 [astro-ph.CO].[13] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski,Phys. Rev. Lett. , 221301 (2019), arXiv:1811.04083[astro-ph.CO].[14] M.-X. Lin, G. Benevento, W. Hu, and M. Raveri,Phys. Rev. D100 , 063542 (2019), arXiv:1905.12618 [astro-ph.CO].[15] T. L. Smith, V. Poulin, and M. A. Amin, Phys. Rev.
D101 , 063523 (2020), arXiv:1908.06995 [astro-ph.CO].[16] K. V. Berghaus and T. Karwal, Phys. Rev. D , 083537(2020), arXiv:1911.06281 [astro-ph.CO].[17] M. Braglia, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu,and K. Koyama, (2020), arXiv:2005.14053 [astro-ph.CO].[18] M. Gonzalez, M. P. Hertzberg, and F. Rompineve, (2020),arXiv:2006.13959 [astro-ph.CO].[19] J. C. Hill, E. McDonough, M. W. Toomey, and S. Alexan-der, (2020), arXiv:2003.07355 [astro-ph.CO].[20] R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine, andA. Mazumdar, Ann. Rev. Nucl. Part. Sci. , 27 (2010),arXiv:1001.2600 [hep-th].[21] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, andJ. Karouby, Int. J. Mod. Phys. D24 , 1530003 (2014),arXiv:1410.3808 [hep-ph]. [22] K. D. Lozanov, (2019), arXiv:1907.04402 [astro-ph.CO].[23] S. Khlebnikov and I. Tkachev, Phys. Rev. D , 653(1997), arXiv:hep-ph/9701423.[24] R. Easther and E. A. Lim, JCAP , 010 (2006),arXiv:astro-ph/0601617.[25] R. Easther, J. Giblin, John T., and E. A. Lim, Phys.Rev. Lett. , 221301 (2007), arXiv:astro-ph/0612294.[26] R. Easther, J. T. Giblin, and E. A. Lim, Phys. Rev. D , 103519 (2008), arXiv:0712.2991 [astro-ph].[27] J. Garcia-Bellido and D. G. Figueroa, Phys. Rev. Lett. , 061302 (2007), arXiv:astro-ph/0701014.[28] J. F. Dufaux, A. Bergman, G. N. Felder, L. Kof-man, and J.-P. Uzan, Phys. Rev. D , 123517 (2007),arXiv:0707.0875 [astro-ph].[29] J.-F. Dufaux, D. G. Figueroa, and J. Garcia-Bellido, Phys.Rev. D , 083518 (2010), arXiv:1006.0217 [astro-ph.CO].[30] L. Bethke, D. G. Figueroa, and A. Rajantie, Phys. Rev.Lett. , 011301 (2013), arXiv:1304.2657 [astro-ph.CO].[31] D. G. Figueroa and T. Meriniemi, JHEP , 101 (2013),arXiv:1306.6911 [astro-ph.CO].[32] L. Bethke, D. G. Figueroa, and A. Rajantie, JCAP ,047 (2014), arXiv:1309.1148 [astro-ph.CO].[33] D. G. Figueroa, J. Garca-Bellido, and F. Torrent,Phys. Rev. D , 103521 (2016), arXiv:1602.03085 [astro-ph.CO].[34] D. G. Figueroa and F. Torrenti, JCAP , 057 (2017),arXiv:1707.04533 [astro-ph.CO].[35] M. A. Amin, J. Braden, E. J. Copeland, J. T. Giblin,C. Solorio, Z. J. Weiner, and S.-Y. Zhou, Phys. Rev. D , 024040 (2018), arXiv:1803.08047 [astro-ph.CO].[36] P. Adshead, J. T. Giblin, and Z. J. Weiner, Phys. Rev.D , 043525 (2018), arXiv:1805.04550 [astro-ph.CO].[37] P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner,Phys. Rev. D , 083534 (2020), arXiv:1909.12842 [astro-ph.CO].[38] P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner,Phys. Rev. Lett. , 171301 (2020), arXiv:1909.12843[astro-ph.CO].[39] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. ,074001 (2015), arXiv:1411.4547 [gr-qc].[40] F. Acernese et al. (VIRGO), Class. Quant. Grav. ,024001 (2015), arXiv:1408.3978 [gr-qc].[41] K. Somiya (KAGRA), (2011), 10.1088/0264-9381/29/12/124007, arXiv:1111.7185 [gr-qc].[42] P. Amaro-Seoane et al. (LISA), (2017), arXiv:1702.00786[astro-ph.IM].[43] M. Punturo et al. , (2010), 10.1088/0264-9381/27/19/194002.[44] N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev.Lett. , 221103 (2001), arXiv:astro-ph/0108011.[45] V. Corbin and N. J. Cornish, Class. Quant. Grav. , , 283 (2000), arXiv:gr-qc/9909001.[47] L. Pagano, L. Salvati, and A. Melchiorri, Phys. Lett. B , 823 (2016), arXiv:1508.02393 [astro-ph.CO].[48] P. D. Lasky et al. , Phys. Rev. X , 011035 (2016),arXiv:1511.05994 [astro-ph.CO].[49] C. Caprini and D. G. Figueroa, Class. Quant. Grav. ,163001 (2018), arXiv:1801.04268 [astro-ph.CO].[50] T. J. Clarke, E. J. Copeland, and A. Moss, (2020),arXiv:2004.11396 [astro-ph.CO].[51] T. Namikawa, S. Saga, D. Yamauchi, and A. Taruya,Phys. Rev. D , 021303 (2019), arXiv:1904.02115 [astro-ph.CO].[52] P. Adshead, J. T. Giblin, T. R. Scully, and E. I.Sfakianakis, JCAP , 034 (2015), arXiv:1502.06506[astro-ph.CO].[53] D. J. E. Marsh, Phys. Rept. , 1 (2016),arXiv:1510.07633 [astro-ph.CO].[54] U. Seljak, U.-L. Pen, and N. Turok, Phys. Rev. Lett. ,1615 (1997), arXiv:astro-ph/9704231.[55] B. Salehian, M. A. Gorji, S. Mukohyama, and H. Firouz-jahi, (2020), arXiv:2007.08148 [hep-ph].[56] J. T. Giblin and E. Thrane, Phys. Rev. D , 107502(2014), arXiv:1410.4779 [gr-qc].[57] J. L. Cook and L. Sorbo, Phys. Rev. D , 023534(2012), [Erratum: Phys.Rev.D 86, 069901 (2012)],arXiv:1109.0022 [astro-ph.CO].[58] N. Barnaby, E. Pajer, and M. Peloso, Phys. Rev. D ,023525 (2012), arXiv:1110.3327 [astro-ph.CO].[59] M. M. Anber and L. Sorbo, Phys. Rev. D , 123537(2012), arXiv:1203.5849 [astro-ph.CO].[60] V. Domcke, M. Pieroni, and P. Bintruy, JCAP , 031(2016), arXiv:1603.01287 [astro-ph.CO].[61] N. Bartolo et al. , JCAP , 026 (2016), arXiv:1610.06481[astro-ph.CO].[62] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi,and F. Takahashi, Phys. Lett. B , 135136 (2020),arXiv:1810.07188 [hep-ph].[63] C. S. Machado, W. Ratzinger, P. Schwaller, and B. A.Stefanek, JHEP , 053 (2019), arXiv:1811.01950 [hep-ph].[64] C. S. Machado, W. Ratzinger, P. Schwaller, and B. A.Stefanek, (2019), arXiv:1912.01007 [hep-ph].[65] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither,A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D.Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr,Computing in Science & Engineering , 62 (2014).[66] A. Kl¨ockner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov,and A. Fasih, Parallel Computing , 157 (2012).[67] A. Kl¨ockner, in Proceedings of ARRAY ‘14: ACM SIG-PLAN Workshop on Libraries, Languages, and Compilersfor Array Programming (Association for Computing Ma-chinery, Edinburgh, Scotland., 2014).[68] L. Dalcn, R. Paz, M. Storti, and J. DEla, Journal ofParallel and Distributed Computing , 655 (2008).[69] L. Dalcn, R. Paz, and M. Storti, Journal of Parallel andDistributed Computing , 1108 (2005).[70] Dalcin, Lisandro and Mortensen, Mikael and Keyes, DavidE, Journal of Parallel and Distributed Computing (2019),10.1016/j.jpdc.2019.02.006.[71] T. E. Oliphant, A guide to NumPy , Vol. 1 (Trelgol Pub-lishing USA, 2006).[72] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haber- land, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson,W. Weckesser, J. Bright, S. J. van der Walt, M. Brett,J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J.Nelson, E. Jones, R. Kern, E. Larson, C. Carey, ˙I. Po-lat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Laxalde,J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero,C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pe-dregosa, P. van Mulbregt, and S. . . Contributors, NatureMethods , 261 (2020).[73] L. Husdal, Galaxies , 78 (2016), arXiv:1609.04979 [astro-ph.CO]. Appendix A: Equations and numericalimplementation
The dynamical system, our conventions, and our nu-merical implementation match those of Refs. [37, 38] (seethe appendices of Ref. [37] for full details), with relativelyminimal differences we describe here.The energy density and pressure which source FLRWexpansion comprise contributions from the standard back-ground ΛCDM model, the axion, and the gauge fields.The spatially-averaged energy density and pressure are ρ φ ( τ ) ≡ (cid:42) φ (cid:48) a + ( ∂ i φ ) a + V ( φ ) (cid:43) , (A1) p φ ( τ ) ≡ (cid:42) φ (cid:48) a − ( ∂ i φ ) a − V ( φ ) (cid:43) (A2)for the axion and ρ A ( τ ) ≡ (cid:10) E + B (cid:11) , (A3) p A ( τ ) ≡ (cid:10) E + B (cid:11) (A4)for the gauge fields. Above we defined the electric andmagnetic fields E i = 1 a (cid:0) A i (cid:48) − ∂ i A (cid:1) , B i = 1 a (cid:15) ijk ∂ j A k . (A5)In our simulations, we consistently include the contribu-tions of Eqs. (A1) to (A4) to the expansion rate alongsidethe background ΛCDM model.The axion–dark-photon coupling violates parity and sothe resonantly produced gauge bosons are helical (at leastinitially), which may in principle imprint on the resultinggravitational waves. For details on the decompositionof gravitational waves into the polarization basis, seeRef. [36].The present-day spectrum Ω gw ( a ) is related to thatat the time of emission, Ω gw ( a e ), viaΩ gw ( a )Ω gw ( a e ) = (cid:18) a e a (cid:19) ρ ( a e ) ρ ( a ) , (A6)since gravitational waves redshift like radiation. In turn,the present-day frequency of observation is determinedfrom the physical momentum k phys = k/a via f = k phys π = 12 π ka e a e a . (A7)The entropy of the Standard Model (SM) provides aconserved quantity by which we may compute a e /a atany time after the SM thermalizes. The entropy density is s = 2 π g (cid:63)S ( T ) T /
45, where T is the SM temperature and g (cid:63)S the number of effective degrees of freedom in entropy.Combining with ρ rad = π g (cid:63) ( T ) T /
30 (where g (cid:63) is theeffective number of degrees of freedom in energy density) allows one to express a e /a in terms of the energy densityin radiation, g (cid:63) , and g (cid:63)S . Approximating g (cid:63)S = g (cid:63) allowsus to evaluate Eq. (A6) at any point after thermalizationas Ω gw ( a ) h = Ω rad ( a ) h Ω rad ( a e ) (cid:18) g (cid:63) ( a ) g (cid:63) ( a e ) (cid:19) / Ω gw ( a e ) . (A8)Likewise, Eq. (A7) may be expressed as f = k/ πa e (cid:112) H ( a e ) M pl (cid:18) Ω rad ( a )Ω rad ( a e ) H M (cid:19) / (cid:18) g (cid:63) ( a ) g (cid:63) ( a e ) (cid:19) / . (A9)While in general g (cid:63) may be evaluated using, e.g., tab-ulated values from Ref. [73], the effective number of rel-ativistic degrees of freedom is fixed to 3 .
36 since wellbefore the time of our simulations. The present radia-tion fraction Ω rad ( a ) h ≈ . × − , and the radiationfraction at the time of emission is calculated during thesimulation. Plugging in H = h · . × − Hz and M pl = 3 . × Hz, f = k/a e (cid:112) H ( a e ) M pl Ω rad ( a e ) / × . × Hz . (A10)To set initial conditions, we begin by solving for thedynamics of the linearized system. We evolve the Fried-mann equations, the homogeneous and linear-order partsof the axion’s equation of motion, (cid:104) φ (cid:48)(cid:48) (cid:105) = − H (cid:104) φ (cid:48) (cid:105) − a d V d (cid:104) φ (cid:105) − a α f (cid:68) F µν ˜ F µν (cid:69) , (A11) δφ (cid:48)(cid:48) ( k ) = − H δφ (cid:48) ( k ) − (cid:32) k + a d V d (cid:104) φ (cid:105) (cid:33) δφ (cid:48) ( k ) , (A12)and the linearized equation of motion for the gauge fieldpolarizations, A ±(cid:48)(cid:48) ( k ) = − k (cid:18) k ∓ αf (cid:104) φ (cid:48) (cid:105) (cid:19) A ± ( k ) . (A13)While solving the linearized system, we compute the gaugefield contributions to the background dynamics by inte-grating over the evolved gauge field modes via numericalquadrature (as employed in Ref. [37]). We neglect thebackreaction of the gauge fields onto the axion fluctua-tions in Eq. (A12), which is negligible even well after weinitialize the full lattice simulation.We begin the linear evolution when H (cid:29) m φ (so thatall modes are far outside the horizon) and set (cid:104) φ (cid:105) = θf and (cid:104) φ (cid:48) (cid:105) = 0. The gauge fields are initialized in theBunch-Davies vacuum, (cid:68) | A ± ( k ) | (cid:69) = 12 √ k , (A14) (cid:68)(cid:12)(cid:12) A (cid:48)± ( k ) (cid:12)(cid:12) (cid:69) = √ k . (A15)The axion acquires a nearly scale-invariant spectrum offluctuations during inflation, (cid:68) | φ ( k ) | (cid:69) = 1 k H π , (A16)and δφ (cid:48) = 0. As a benchmark value, we take the inflation-ary Hubble scale H inf = 2 × GeV, but our main resultsare insensitive to whether the axion is even initializedwith any fluctuations.We evolve the linear system until H = m φ , at whichpoint we initialize the lattice simulation. The subsequentdynamics are identical to when the lattice simulation isinitialized when H = 10 m φ , and are relatively unchangedeven if initialized at H = m φ /
10. Using the backgroundvalues and power spectra obtained at this point ( H = m φ ),we initialize the lattice simulation with the proceduredescribed in Ref. [37]. Our simulations use a grid with N = 768 points and conformal box length L = 10 /m φ ,with a timestep ∆ τ = ∆ x/
10 = L/ N . Appendix B: Convergence tests
The simulations we present in the main text requirefairly large volumes to capture the IR tail of the gravita-tional wave spectrum which resides at CMB frequencies,while also sufficient resolution to capture both the initialtachyonic resonance band and the subsequent power trans-fer to higher momenta via nonlinear effects. As N = 768 represents the largest grid size possible with our currentresources, in lieu of an ideal convergence test (fixing thebox length L while increasing N ) we compare results fortwo box lengths, both with N = 768. We present thiscomparison for L = 10 /m φ (as used for main results) and L = 5 /m φ in Fig. 4, for α = 70 (the largest coupling forwhich we present results, which should require the mostresolution). Specifically, Fig. 4 depicts for each field thefinal dimensionless power spectra, i.e., for a field f ,∆ f ( k ) = 12 π V (cid:90) dΩ4 π k | f ( k ) | . (B1)The results show excellent consistency, with expecteddiscrepancy at the largest scales in the simulation (dueto statistical variance) and negligible differences at theNyquist frequency. Appendix C: Comparing to the two-fluid description