Gravitational Wave Symphony from Oscillating Spectator Scalar Fields

We investigate a generic source of stochastic gravitational wave background due to the parametric resonance of oscillating scalar fields in the early Universe. By systematically analyzing benchmark models through lattice simulations and considering a wide range of parameters, we demonstrate that such a scenario can lead to detectable signals in gravitational wave detectors over a broad frequency range and potentially address the recent findings by pulsar timing array experiments. Furthermore, these models naturally yield ultralight dark matter candidates or dark radiation detectable by cosmic microwave background observatories.

Introduction.The discovery of gravitational waves (GWs) by the LIGO-Virgo collaboration in 2015 [1] heralded the beginning of a new era of observational astronomy, as well as unprecedented opportunities for new discovery in particle physics and cosmology.The very recent compelling evidence for a stochastic GW background (SGWB) reported by the pulsar timing array (PTA) experiments such as NANOGrav [2], EPTA+InPTA [3,4], PPTA [5], and CPTA [6] with plausible cosmological origins [7], further propels the emerging research field of using SGWB to probe the pre-Big Bang Nucleosynthesis "primordial dark age" [8][9][10][11][12] and the associated fundamental physics.Among the well-motivated sources of cosmogenic SGWBs, a few classes have been widely and systematically studied, such as those originating from inflation, cosmic strings, and phase transitions [14,18].Another generic source for SG-WBs is parametric resonance excited during the oscillation of a scalar field, leading to rapid, nonperturbative particle production and thus large time-dependent field inhomogeneities [15][16][17][18].Existing studies in this direction mostly focus on specific scenarios, such as those related to inflationary preheating or Affleck-Dine baryogenesis [5-8, 21, 22, 25, 26].Despite being highly motivated, the additional constraints associated with these models (e.g. on inflation or low-energy supersymmetry breaking scale) limit the detectability of their GW signals: either the frequency is too high, or the signal strength is too small, for observations in foreseeable future (for exceptions, e.g., by utilizing axion-gauge or axion-gravity coupling, see [27][28][29][30]).
In this Letter, we systematically expand the studies on GWs sourced through parametric resonance by considering a generic spectator scalar field displaced from its true vacuum during inflation.After inflation, it rolls toward its potential minimum and oscillates coherently until (some of) its energy dissipates into particle excitations and GWs.We investigate simple, minimal benchmark models with renormalizable scalar field potentials, allowing model parameters (masses and couplings) to vary over a wide range.We found detectable SGWB signals with peak frequency lying in a wide range from O(0.01) nHz to GHz, and peak signal strength up to Ω GW ∼ 10 −9 .Notably, when the (effective) mass of the scalar field is O(10 −13 ) eV, this class of models may provide another plausible cosmological explanation for the recent observation by pulsar timing arrays, complementing the existing literature on this timely topic.Interestingly, in such PTA-relevant models, the massive scalar field also serves as an ultralight dark matter (DM) candidate.Furthermore, models involving massless scalar field(s), a form of dark radiation, potentially lead to ∆N eff observable with upcoming cosmic microwave background (CMB) and largescale structure experiments.Therefore, our study demonstrates an intriguing complementarity between GW signal and dark matter or CMB physics that generically emerges from these types of models.
The rest of the Letter is organized as follows.In the next section, we present benchmark models to be examined, along with analytic expectations for GW production.Following that, we show the numerical results for GW signals and the complementary observables (DM and/or N eff ).We conclude with an outlook.
Benchmark Models and Estimates for GW Signals.The equation of motion for a scalar field oscillating in an Friedman-Robertson-Walker background is given by For a massive scalar field, the last term is ∂V /∂ϕ ⊃ m 2 ϕ ϕ.For the remainder of this Letter, we assume that the Universe is radiation dominated throughout the time range of interest, i.e. ρ ϕ is a subleading component, Pl , and H = 1/2t (M Pl is the Plank scale).When H ≫ m ϕ , the field is frozen and only starts rolling or oscillating significantly after H osc ∼ m ϕ .For a massless scalar field with V ⊃ λϕ 4 , the field is again frozen at early times, when H dominates over V ,ϕϕ , and starts rolling appreciably when H osc ≃ |V ,ϕ /ϕ| = λ ϕ ϕ in , where ϕ in is the initial field amplitude (see, e.g., [31]).We take the initial field displacement, ϕ in , as set during inflation.It typically needs to be up to ∼ M Pl for sufficient GW signal, as we found, which can be realized in various ways.First, it could be acquired due to quantum fluctuations during inflation, where ϕ in ∼ M Pl requires many e-folds of inflation (N ) and a small quartic coupling λ in the case of a self-interacting ϕ (see Supplemental Material for a detailed derivation [32]).Second, such a large displacement can be realized by introducing nonrenormalizable terms in the potential, which shift the vacuum expectation value of ϕ during inflation [33].Finally, a large ϕ in could be simply due to an ad hoc initial condition, similar to the case of the inflaton.Our result is independent of the specifics of how ϕ in is realized.
We consider the possibilities of both a massive and massless scalar field ϕ and choose two different models in each case, including only renormalizable terms (up to quartic).These can be summarized as follows In models A and B, a free massive scalar ϕ couples to a massless scalar χ through a quartic (A) or trilinear (B) term.The second scalar field χ is necessary in this case, as a free massive scalar does not exhibit parametric self-resonance.In contrast, for a self-interacting scalar, such a coupling is not necessary, and the quartic scalar potential is sufficient to generate parametric resonance and GWs (C).Finally we consider a massless quartic field ϕ coupled to a massless field χ through a quartic interaction (D).In light of the constraints related to BBN, we require that H osc > H BBN , leading to m ϕ ≳ 10 −16 eV and λ ϕ ϕ in ≳ 10 −16 eV for the massive and massless cases, respectively.The lower bound of m ϕ ∼ 10 −16 eV corresponds to a GW frequency of O(0.01) nHz.
In order to maximize the efficiency of parametric resonance, we keep the χ field VEV at zero.During inflation, the quantum fluctuations of χ can be suppressed by the coupling to the inflaton, stabilizing it at the origin with a Hubble-scale mass [33].Even if χ acquires a VEV during inflation, a coupling to the thermal plasma after reheating can lead to a large thermal mass that dynamically drives the VEV to zero.This provides interesting avenues for complementarity studies by considering various couplings of the χ to the standard model (SM), which we leave for future work.
For a simpler classification of the models into massive and massless and to minimize the number of model parameters, we did not include combinations such as λ ϕ ϕ 4 /4 + m 2 ϕ ϕ 2 /2.This does not reduce the general applicability of our models for the following reasons.If, during the rolling of the ϕ field, the motion is dominated by the quadratic term, the system is well-described by models A and B. If the two terms are comparable at first, the first few oscillations will be anharmonic, but the quartic term would quickly become subdominant to the quadratic.If, alternatively, the quadratic term is negligible relative to the quartic, we can follow the evolution of models C and D until the mass becomes important and the ϕ particles become nonrelativistic.It is straightforward to adapt our analysis to include these variations.
We separate the field ϕ as ϕ(x, t) = ϕ bg (t)+δϕ(x, t), where ϕ bg (t) is the background value, which follows Eq. (S6) with ∇ 2 ϕ bg = 0.The fluctuation component can be expanded in Fourier modes as The χ field background value is assumed to be 0, thus χ(x, t) = δχ(x, t), and can be similarly decomposed into modes δχ k (t).
As ϕ bg oscillates around the minimum of its potential, the fluctuations δϕ and δχ exhibit parametric resonance [34], similar to preheating following inflation [35,36].At the linear level, one needs to solve Eq. (S6) for ϕ bg (with ∇ 2 ϕ bg = 0), and the equations of motion of the fluctuation are Parametric resonance amplifies a certain range of wave numbers, typically up to O(m ϕ ).We can thus approximate the spectrum of scalar fluctuations as peaking around a characteristic scale k c , which depends on model parameters and is proportional to H osc .GW production from peaked sources can be estimated by the rule of thumb as given in [37] In these formulas, α is the fraction of the energy in the GW source (the amplified scalar fluctuations in our case) relative to the total energy density of the Universe, and β encodes the anisotropy of the source, which can be extracted from simulations, with β = O(0.01− 0.1) as a reasonable range [38].
The equation of state of the Universe at the time of GW emission is parametrized by w = 1/3 for production during the radiation-dominated era.σ is the width of the source in momentum space, which is estimated to be σ ∼ k c for peaked sources.For a massive field, we typically get k c ∼ m ϕ ≃ H osc .However, the Hubble scale during GW emission H c (close to the end of parametric resonance) is smaller than H osc due to the expansion of the Universe during the span of the parametric resonance era.The exact ratio H osc /H c is model-dependent but can be as large as 10 4 as per our numerical finding.Based on our simulations, we found ν peak GW ∼ O(10 12 ) H osc /M Pl Hz to be an accurate estimate.For a massless field, our simulation results agree with Eq. ( 9) taking H c = H osc .Note that ν peak GW includes the redshift after production, which depends on thermal history.We assume the standard radiation-dominated for concreteness, while a prolonged matter-dominated phase would modify our results.
The transverse-traceless tensor perturbations h ij follows the linearized equation where Π T T ij represents the transverse-traceless part of the effective anisotropic stress tensor Π ij = ∂ i ϕ∂ j ϕ + ∂ i χ∂ j χ.The quantity of interest is the energy density power spectrum of the GWs normalized by the critical energy density of the Universe today, ρ c .
where the differential dΩ k is the solid angle element in Fourier space and the integral is over the Fourier modes h ij (k, t).The simulations were performed using the publicly available code CosmoLattice [39,40] (the numerical setup and convergence tests are presented in Sec.B.1, B.4 of the Supplemental Material [32]) Numerical results of the GW signals and complementary phenomenology: Dark matter and dark radiation.In Fig. 1 and Table I, we illustrate the numerical results from our simulation for a set of benchmark models, which demonstrate detectable GW signals over a wide frequency range.In the following, we elaborate on the relevant dynamics to provide physics insights for interpreting these results.
For the massive models A and B, ϕ bg starts rolling and oscillating when H osc ≃ m ϕ , and the ratio of the energy density in the ϕ field over the total energy density [roughly equals the parameter α in Eq. ( 10 Pl ), leading to, e.g., ϕ in /M Pl = 0.77, 0.25 for 10%, 1% ratio, respectively.
In order to get sufficient GW signals, the transfer of a significant fraction of ρ ϕ into field fluctuations is required.The modes δχ k undergo parametric resonance since the interaction with the inflaton or thermal bath that stabilized χ at zero becomes inactive at the late times when ϕ oscillation starts.Because of rescattering effects, δϕ k modes will also be generated; however, they can be at most comparable to δχ k in power.Thus, we can safely neglect the δϕ k modes in our analytic estimates (but included in the numerical simulation).By measuring t, k in units of m ϕ , the last term of Eq. ( 8) becomes q A,B f (t)δχ k , where q A ≡ gϕ 2 in /m 2 ϕ and q B = σϕ in /m 2 ϕ are dimensionless quantities controlling the strength of parametric resonance in models A and B, respectively.f (t) is a decaying periodic function, bounded by ±1 and encodes the oscillating background field ϕ bg .
The dimensionless form of Eq. ( 8) indicates that during parametric resonance the amplification of δχ k is independent of m ϕ and evolves as δχ k ∝ e µ k t , where µ k is called the Floquet exponent.After a time t, the ratio of the energy density of δχ modes to the ϕ energy density scales as Pl ∼ (m ϕ e µt /M Pl ) 2 , since we take ρ ϕ to be a fraction of the total energy of the Universe and H ∼ m ϕ during the (early stages of) parametric resonance.We thus see that the amplification required to achieve a similar energy transfer efficiency must be larger for smaller values of m ϕ .By measuring both µ and time in units of m ϕ , we deduce that for smaller values of m ϕ , parametric resonance needs to be active for more oscillation periods of the background field ϕ bg , which agrees with our numerical simulations.Consequently, the endvalue of ρ δχ /ρ tot is largely independent of m ϕ ; thus we expect Ω GW also to be m ϕ -independent while ν GW can vary without affecting Ω GW [41].In reality, there is a mild dependence on m ϕ , due to the different redshift behaviors of matter and radiation during the prolonged parametric resonance, which manifests itself as a small suppression of Ω GW for nHz frequencies (see Sec. B.2, B.3 of the Supplemental Material for more details [32]).
Models C and D exhibit similar dynamics as models A and B. The ϕ field starts rolling or oscillating at H osc ≃ λ ϕ ϕ in , whereas the initial field amplitude is given by the energy ratio Pl ).The initial field amplitude is ϕ in /M Pl = 1.17, 0.35 for 10%, 1% energy ratio.In model C, parametric resonance amplifies δϕ fluctuations without the need for a second field χ.The resulting GW spectrum exhibits a multipeak structure, which is absent in all other models and provides a distinct feature.The origin of the multipeak structure in model C can be best understood by comparing it to the similar model D. We see that these features are accompanied by an overall reduction in the peak GW amplitude.Self-resonance in model C is weaker than the parametric resonance of χ in model D (and similarly in models A and B).A weak resonance amplifies specific waven umbers instead of a broad range of modes, leading to a multipeak structure, which would be washed out for more efficient parametric resonance.A similar effect appears in oscillon preheating, where strong resonance washes out multi-peak features in the GW signal [17].
We find that for all the models we consider, in order to produce sufficient GW signals for detection, the initial displacement ϕ in is typically driven to be close to M Pl (see ϕ in dependence of ρ ϕ /ρ tot ).This requires a small self-coupling for models C, D in order to yield ν GW in the observable range.Similarly, for models A and B, the dimensionless couplings between the fields, g or σ/m ϕ , also must be small.Since the resonance parameters must be q A , q B ∼ O(10 − 10 3 ) for efficient parametric resonance, g ∼ m 2 ϕ /M 2 Pl ≪ 1 and σ/m ϕ ∼ m ϕ /M Pl ≪ 1.These small dimensionless couplings, as preferred, are nevertheless technically natural and thus innocuous theoretically.Furthermore, certain UV completions may provide even more satisfying ways to address the naturalness of such small couplings.For example, the quartic self-coupling of ϕ may arise from a loop diagram involving Yukawa coupling of ϕ to heavy intermediate sterile fermions.A moderately small Yukawa coupling of values as present in the SM, along with the loop factors and mass suppression from the heavy fermion, can readily yield the very small ϕ couplings of our interest.Alternatively, the preferred small quartic couplings may naturally arise from exponential factors (with originally mild hierarchy in parameters) in the framework of a deformed CFT with a dual AdS interpretation, where the scalar is a pseudo-Goldstone boson (dilaton or radion), as suggested in [42].Note that by identifying ϕ as a Goldstone boson or embedding it in a supersymmetric framework, the concern about radiative stability of m ϕ , which is generically present for scalar fields including the SM Higgs, can also be addressed.Nevertheless, the detailed embedding into UV realizations to address naturalness considerations is beyond the scope of this Letter, which focuses on the intriguing phenomenology arising from generic effective field theory models.Next, we discuss other phenomenological aspects that naturally complement the GW signals originating from these models.A particularly intriguing possibility arises from models A and B involving a massive scalar field: the scalar field itself can be identified as dark matter.Let us consider the state of the scalar field at the time t end when parametric resonance production of χ particles ends, and t end /t osc ∼ 10 3 − 10 5 , depending on m ϕ .The energy density of the ϕ condensate at t end is ρ ϕ,end = 1 2 m 2 ϕ ϕ 2 end and afterward redshifts as matter, i.e., ρ ϕ = ρ ϕ,end (a end /a) 3 (a: scale factor).Numerically we find e.q.ϕ end ∼ O(10 −5 )ϕ in for efficient parametric resonance in model A. By entropy conservation, g * a 3 T 3 is constant (T : cosmic temperature).The relic abundance of ϕ particles today can then be estimated as where the subscripts 0 denote current-day values, we have assumed that the ϕ field is stable, in particular, unable to decay to SM particles.This stability assumption is consistent with the models presented and can be reinforced with symmetry protection.For model A, Ω ϕ,0 monotonically increases with m ϕ , and for m ϕ ≳ 10 −13 eV it is larger than the observed Ω DM ≈ 0.2, causing overclosure problem.Nevertheless, an interesting, viable DM scenario emerges for ultralight ϕ with m ϕ ∼ 10 −13 eV, which also produces a signal in the nHz band in addition to Ω ϕ,0 ∼ Ω DM .Meanwhile, the larger m ϕ range for model A, e.g., those relevant for LISA and LIGO (m ϕ ∼ 10 −2 eV for ν GW ∼ 10 −3 Hz and m ϕ ∼ 10 6 eV for ν GW ∼ 10 Hz), can be viable if ϕ is not DM, e.g., by introducing a coupling of ϕ to sterile neutrinos as a minimal extension of the model, thus allowing ϕ to decay into radiation modes (possibly SM neutrinos at the final stage).
Unlike model A, model B allows for ϕ decay into pairs of χ particles with a decay rate Γ ϕ→2χ = σ 2 /(8πm ϕ ).For Γ ϕ→2χ ≲ H 0 , model B suffers from the same overclosure problem for m ϕ ≳ 10 −13 eV.If Γ ϕ→2χ ≃ H BBN , ϕ would be subject to BBN constraint, which imposes a condition on model parameters.For m ϕ ≳ σ ≳ 5 × 10 −8 m ϕ × eV the model satisfies both the perturbative unitarity and the condition that ϕ decays before BBN, which we elaborate in the Supplemental Material [32].However, increasing σ beyond a certain value can lead to suppression of parametric resonance, making it challenging to find a viable parameter set.We leave such a dedicated parameter scan for future work.
All four models we consider yield new relativistic degrees of freedom beyond the SM, i.e., dark radiation, including the GWs and/or massless scalar(s), potentially contributing to ∆N eff relevant to BBN prediction and CMB observations.The current limit of |∆N eff | ≲ 0.29 at 95% C.L. [13,14] is set by the Planck satellite.Next-generation CMB and large-scale structure experiments (CMB-S4 [15], COrE [16], Euclid [17]) will be able to probe ∆N eff = O(0.01).Our models can yield ∆N eff up to O(0.1), while being consistent with the current bound, making these simple models relevant for both the upcoming GW and CMB experiments, offering complementary signals (see Supplemental Material for details on estimating ∆N eff [32]).Discussion and Conclusion.The presence of oscillating scalar fields in the early Universe naturally arises in many wellmotivated particle physics theories beyond the SM, with the potential of sourcing SGWBs through parametric resonance effects.In this Letter, we conducted a general study of this scenario with a focus on simple representative models and surveyed a wide range of mass and coupling parameters using lattice simulations.Depending on model specifics, we found that such models can give rise to GW signal with amplitude up to Ω GW ∼ 10 −9 , with frequency ranging from the band of PTA experiments to that of LISA and LIGO/ET/CT [48][49][50][51] and beyond [2,[52][53][54][55][56], and with a distinct, multipeak spectrum in some cases.Intriguing complementary phenomenology  The energy density of the scalar field is assumed to be 10% of the total energy density of the Universe at the start of the simulation (solid lines in Fig. 1).ΩGW quadratically dependents on the ratio ρ ϕ /ρtot.Finally, the points denoted with a * lead to DM overproduction for stable ϕ, but can be viable with additional couplings that allow ϕ to decay (see main text for details).
is also identified: in the models with a massive scalar ϕ, the residual relic ϕ field can naturally serve as DM candidate, and in particular m ϕ ∼ 10 −13 eV marks a viable DM mass range that correlates with GW signal covered by PTA experiments; in the models with massless scalars, the relic scalar radiation contributes to CMB N eff observable (in addition to the generally subleading contribution from GWs).Further explorations of these scenarios may lead to even richer findings.For instance, while we focused on considering the standard cosmic history of radiation domination during the ϕ oscillation period, early matter domination driven by the ϕ condensate may amplify the GW signal and affect DM structure formation and DM detection today [57,58].In addition, it is worth further investigating the natural connection with ultralight DM as emerging in some benchmarks and its ramifications, as well as the potential complementarity with CMB observations [59] We thank Anish Ghoshal for collaboration during the early stage of this project.We also thank Soubhik Kumar and Raman Sundrum for helpful discussions.YC is supported by the US Department of Energy under award number DE-SC0008541.EIS is supported by a fellowship from "la Caixa" Foundation (ID 100010434) and from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 847648, the fellowship code is LCF/BQ/PI20/11760021. PS is supported by Grant-in-Aid for Scientific Research (B) under Contract No. 23H01177 and was supported by Seoul National University of Science and Technology during the earlier stage of this work.Numerical computation in this work was carried out at the Yukawa Institute Computer Facility.Some earlier simulations were performed at the Dark cluster at Seoultech.

A. Inflationary fluctuations and Initial conditions
We consider two scalar fields ϕ and χ coupled to each other during inflation but not to the inflaton.As light spectator fields, they acquire a vacuum expectation value through de Sitter fluctuations, as discussed in Ref. [1] for a single self-interacting scalar field.Consider the most general potential that encompasses all models A, B, C, and D: The equations of motion are We then multiply each equation by ϕ and χ, respectively, and take the average values.Following Ref. [1] we get where we considered ⟨ϕ⟩ = 0 and thus the trilinear coupling σϕχ 2 does not appear in the equation for the variances.Furthermore, we wrote the equations in terms of the number of e-folds N = N t.Both fields start at ⟨ϕ 2 ⟩ = 0 = ⟨χ 2 ⟩ and initially start growing as (H 2 /4π 2 )N for m ϕ ≪ H.We can solve for the late-time (asymptotic) values by setting ∂/∂N = 0. We assume g > 0 for simplicity, and obtain If the two fields are decoupled (g = 0), we recover the familiar single-field result In our simulations, we require only one of the fields, ϕ, to have a non-zero VEV, while the other field, χ, which is present in models A, B, and D, to be at the origin.If both fields are light spectators during inflation, this will not be true.We examined some numerical examples with a non-zero initial VEV for χ, which are found to lead to a suppression of parametric resonance.Therefore, ⟨ϕ 2 ⟩ ≫ ⟨χ 2 ⟩ is required to achieve efficient parametric resonance.As discussed in the main text, this condition may not be ad hoc but rather can be realized naturally realized by coupling the χ field to inflaton or particles in the thermal bath (including SM states).Such couplings could introduce a large mass term for χ, stabilizing it at the origin during or after inflation but before the ϕ field starts rolling and oscillating.Now, let us consider the dynamics of the ϕ field alone.As the ϕ mass can be neglected during inflation, the field variance ⟨ϕ 2 ⟩ grows linearly with the number of e-folds of inflation.A requirement for reaching a VEV of the order of the Planck mass is N ∼ 100(M Pl /H I ) 2 ≳ 10 12 , where we took H I ≲ 10 −5 M Pl as constrained by the CMB data.We thus see that a very large number of e-folds are required.If a quartic self-coupling is present, the distribution of ϕ reaches an equilibrium of Eq. (S5), whereby for a Planckian field VEV, the self-coupling must be λ ϕ ≲ 10 −20 , depending on the ratio H I /M Pl .In the case of a massive light free scalar field during inflation, the equilibrium distribution leads to ⟨ϕ 2 ⟩ ∼ H 2 /m ϕ , which can also reach Planckian values for m ϕ ≪ H.As an example, we can consider one of the benchmark points for model B as shown in the main text, with m ϕ = 10 −13 eV.Assuming that ϕ in arises solely due to de-Sitter fluctuation, we would be lead to a definite prediction for the Hubble scale of inflation.In particular, in order to produce ⟨ϕ 2 ⟩ ∼ M 2 Pl as we found for producing sufficient Ω GW , we found that H I ∼ 10 −2 GeV is required, along with the necessity of very many e-folds of inflation.Nevertheless, as noted, ϕ in may be dominantly determined by an initial condition, without imposing specific requirements on inflation.Finally, we must note that the quantum fluctuations that generate ϕ in are uncorrelated with the density perturbations sourced by the inflaton, leading to isocurvature perturbations if ϕ is stable.For an initial displacement ϕ in , the isocurvature fluctuations of the ϕ field scale as P S ∼ (δϕ/ϕ in ) 2 ∼ (H I /ϕ in ) 2 (see e.g.Ref. [4]).For the desirable value of ϕ in ∼ M Pl we get P S ∼ H 2 /M 2 Pl .CMB data constrains the amount of isocurvature perturbations compared to the adiabatic density perturbations as P S ≲ 0.04P ζ /f , where P ζ = 2.2 × 10 −9 , f ≤ 1 is the fraction of ϕ energy in the total dark matter (model A and B) or total radiation energy (model C and D).Putting everything together, we found the resultant constraint is H I ≲ 10 −5 M Pl /f , which imposes a constraints on H I that is comparable to or weaker than what is already required by tensor to scalar ratio in CMB data.Our benchmark model parameters are fully consistent with this bound.

B. Lattice simulation details 1. Lattice simulations set-up
We solve the evolution equations for the two fields tat each lattice point, with the scale factor for the fixed background expansion where a(t * ) and H(t * ) are, respectively, the scale factor and Hubble parameter at the beginning of the simulation at t = t * .For the case of expansion in the radiation-dominated era, we consider the constant equation of state parameter (w) for the external fluid sourcing this expansion to be w = 1/3.This expression can be readily used to get the Hubble expansion rate for any given background of cosmic evolution, such as radiation or matter domination.The GWs are computed using the linearized equation of motion of the transverse-traceless tensor perturbations h ij : where Π T T ij represents the transverse-traceless part of the effective anisotropic stress tensor Π ij = a ∂ i ϕ a ∂ j ϕ a for the system of scalar fields ϕ a , {a = 1, 2, • • • }.The GW energy density is given by where ⟨• • • ⟩ V denotes a spatial average over V. Since we are in a radiation-dominated Universe, assuming that the expansion of the Universe obeys entropy conservation from the end of simulation (subscript "e") to the present time (subscript "0"), the spectrum of the energy density of GWs (per logarithmic momentum interval) observable today is [5][6][7][8]: where the critical density of the Universe is ρ crit = 3M 2 Pl H 2 and Ω rad,0 h 2 = h 2 ρ rad,0 /ρ crit = 4.3 × 10 −5 .Solving these equations numerically on a lattice requires appropriately rescaling the variables so that the equations are reformulated in terms of dimensionless quantities for numerical stability.We worked with the following rescaling for the field and spacetime variables: where the subscript "pr" refers to "program" variables.The parameter α is chosen from the exponent of the dominating (classical) scalar component, i.e., for a scalar potential whose leading term around the potential minima is V (ϕ) ∝ ϕ n , we have α = 3(n − 2)/(n + 2).The effective mass is defined as: The tensor perturbations are rescaled in a similar way.

Dimensionless equations
We can understand the dynamics of parametric resonance and the scaling between different parameter choices (leading to different GW frequencies) by using properly rescaled equations.We consider Model A for concreteness, where the equations for the background ϕ field and δχ k modes are where we divided the equation of motion for ϕ by m 2 ϕ ϕ in and the equation of motion for δχ k by m 2 ϕ .We first see how the resonance parameter q arises in this model.Furthermore, by defining t = m ϕ t, H = H/m ϕ , k = k/m ϕ and defining a function f ( t) which contains the time-dependence of the background, such that gϕ 2 /m 2 ϕ ≡ q f ( t), we see that the equation of motion for the χ fluctuations is written as The equation has "lost" all information about specific scales thus, the growing modes will behave as δχ k ∝ e μk t, where μ is the Floquet exponent of Eq. (S17) and is dimensionless.In physical units μt = (μm ϕ )t, where the dimensionful physical Floquet exponent is µ = μ • m ϕ .Thus, we can say that we measure time, Floquet exponents and wavenumbers in units of m ϕ .The only "memory" of the mass-scale m ϕ comes through the initial conditions of the χ fluctuations, which follow the Bunch-Davies vacuum and thus without amplification from parametric resonance ρ δχ ∼ O(H) 0 ϕ .The amplification will effectively multiply this result by e 2μ t and parametric resonance completes when ρ δχ ∼ ρ ϕ ∼ m 2 ϕ M 2 Pl or μt ∼ log(M Pl /m ϕ ).Because of this fact, parametric resonance will take longer to transfer the energy from the ϕ field to χ fluctuations for smaller values of m ϕ , corresponding to smaller GW frequencies.

Parameter choice
We can demonstrate the various parameter choices with a specific example.Let us fix the dynamics to take place around H = 1 GeV.In a massive scalar field model (e.g., model A), we choose m ϕ = 1 GeV.Choosing when the field starts rolling largely determines the frequency of GWs.Since we want the scalar field to have subdominant but significant fraction of the energy  [9][10][11][12], the UV parts of the spectra for some choices of kIR showing a growth as k 4 are numerical artifacts and do not correspond to any physical result.With an optimal choice of lattice parameters (for instance, the blue dash-dotted line in this figure), such UV "peaks" can be avoided.than the SM due to asymmetric reheating and never thermalized with the SM due to feeble interactions [20].Therefore we choose to focus on the latter, non-thermal component.In order to estimate this contribution, we make a few simple, reasonable assumptions.First, after production, the dark radiation does not thermalize with the SM, which is consistent with our minimal models where non-gravitational coupling to the SM is absent.Second, we assume that particle production from the oscillation of the ϕ condense completes quickly (i.e., within O(1) Hubble time) after the onset of oscillation at T osc (in analogy to "instant reheating"), and the fraction of energy released in the form of dark radiation is simply parametrized by a constant ξ.This assumption is generally consistent with the efficient parametric resonance effect we consider here for detectable GW signal and enables a simple yet reasonable estimate of the effect.For models C and D, ξ = 100% is easily justified, while for models A and B, ξ ≲ 1 due to the possible presence of residual massive ϕ, in the form of a condensate or fluctuations.Nevertheless, for models A and B, ξ = 100% can be taken for a conservative estimate for the maximal ∆N eff .A more accurate determination of ∆N eff from massless χ (ϕ) from ϕ condensate requires a dedicated numerical tracking of the evolution of the particles following production, which we will defer for further investigation.The estimate given here is sufficient to provide the rough estimate/upper bounds necessary for our current analysis.We may now proceed to estimate ∆N eff following the arguments of [21,22].Note that most existing literature considers the scenario where the dark radiation sector was once in thermal equilibrium with the SM, then thermal decoupling occurs at a later time, after which comoving entropy is conserved separately in each sector.The scenario we consider here is very different in that we consider non-thermally produced dark radiation that never equilibrates with the SM and may or may not self-thermalize following the production (depending on the model couplings); thus, entropy conservation may or may not apply well for the dark sector.Instead of tracking the energy or entropy evolution in the two sectors starting from the time of thermal decoupling, here, the proper starting point is the onset of ϕ oscillations at T osc .For a straightforward application of the procedure of ∆N eff as given in [22], we assume that self-coupling is effective enough to thermalize the dark radiation species among themselves, thus entropy conservation can apply to the dark sector.By applying entropy conservation in the SM and dark sector starting from T osc , we find: where the quantities withˆ(all in the numerator) are of the dark sector, while those withoutˆ(in the denominator) are of the SM; g * , T are evaluated at any time after T osc .Strictly speaking, the effective number of relativistic degrees of freedom g * in the above equation should be g * s ; however, for the estimate we do here, it is the same as g * .Next, we apply the already introduced input parameter α, which is the ratio of total energy density stored in ϕ condensate over the total cosmic energy density (assumed to be radiation domination) at T osc , as well as the ξ parameter introduced earlier as the fraction of ϕ energy released into dark radiation, and find:  10)% yet much lower T osc , ∆N eff could be in conflict with the Planck bound.Nevertheless, keep in mind that for models A and B, ∆N eff as estimated here is a conservative upper bound.For α ∼ O(1)%, ∆N eff ∼ O(0.01), thus is generally safe from existing bound while being detectable with future observations.In the main text, we have given the estimates of H osc for the models we consider, which relates to T osc according to H 2 osc = (π 2 g/90)T 4 osc /M 2 Pl .Therefore, the result given here can be readily converted to the constraint/region of interest for model parameters such as m ϕ and the couplings that determine H osc .In the following, we will give a concrete example based on a simplified numerical analysis, which is consistent with the estimates given above.Consider model B with m ϕ ∼ 10 −13 eV, leading to a GW frequency in the nHz range.Here χ represents the dark radiation contributing to N eff .By numerically tracking our simulation, we find that approximately 60% of the total ϕ energy is transferred to χ at the end of the simulation for typical values of coupling, as shown in Table 1.Hence, we take ξ ≃ 0.6 and neglect the variation of the effective number of relativistic d.o.f's, and find ∆N eff = 0.016 assuming α = 0.01 initially.

D. Trilinear coupling, backreaction, and DM
It is worth studying the evolution of the background field ϕ in more detail.The equation of motion for the ϕ field is When considering only the spatially averaged field value ϕ av , we drop the gradient term.The treatment of the term χ 2 requires more care.Initially, it is a subleading term when the field χ starts in its Bunch-Davies vacuum.As χ grows, at some point, the variance of χ becomes comparable to ϕ 2 av , and back-reaction becomes important.We can then use the Hartree approximation to write [23] φav + 3H φav + m 2 ϕ ϕ av + σ⟨χ 2 ⟩ = 0 (S24) Neglecting the second time derivative and the Hubble friction, we find thus ϕ av follows the evolution of ⟨χ 2 ⟩.Another interesting feature of model B is the existence of a perturbative decay channel of ϕ into pairs of χ particles with a rate where we used the definition of the resonance parameter q B ≡ σϕ in /m 2 ϕ .To achieve efficient parametric resonance, one must choose q B > 1.The simulations for the benchmarks presented in the main text were performed using q B = 100.By using ϕ 0 ∼ M Pl and σϕ 0 /m 2 ϕ = 100, we can eliminate σ from the expression of the decay rate, leading to Γ ϕ→2χ ≃ 10 3 m 3 ϕ /M 2 Pl .For m ϕ ≳ 10 3 GeV, ϕ decays before BBN, yielding a viable scenario.However, the GW frequency corresponding to such masses would be typically higher than the sensitive range of LIGO/CE/ET.For m ϕ ≲ 10 7 eV, the ϕ field has not decayed by today, causing concern of over-closure as for part of the parameter region in model A, which can be alleviated by invoking additional coupling of ϕ to sterile neutrinos.For much smaller masses of m ϕ ≲ 10 −13 eV, ϕ is potentially a DM candidate.

FIG. 1 .
FIG. 1. Example resultsfor Models A, B, C, and D (red, blue, brown, and orange), corresponding to the parameter choices given in TableI.Solid and dashed lines indicate the results when the initial energy in the scalar is assumed to be 10% and 1% of the total energy density of the Universe, respectively.Note that model C has multipeaked oscillatory features.The constraints and sensitivities of various GW experiments are also shown.

3 N = 512 3 ,
FIG.S1.The GW spectrum for Model B for different box sizes, exhibiting convergence in the central region of wave-numbers, allowing us to resolve the k 3 scaling at low wave numbers.A noted earlier in the literature[9][10][11][12], the UV parts of the spectra for some choices of kIR showing a growth as k 4 are numerical artifacts and do not correspond to any physical result.With an optimal choice of lattice parameters (for instance, the blue dash-dotted line in this figure), such UV "peaks" can be avoided.
FIG. S2.The value of ϕ averaged over the simulation box as a function of time for model B. With the introduction of the trilinear term, the back-reaction from the χ quanta leads to a shift in the global minima of the potential, causing ϕ to acquire a non-zero vacuum expectation value (VEV), as evident from the simulation results shown in this figure.

TABLE I
. Example model parameters leading to GW production with frequencies in observationally relevant range, from NANOGrav to LISA to LIGO.The values of ΩGW indicate the peak amplitude.
S21, and evaluating just above SM neutrino coupling at T ∼ O(10) MeV when g * = 10.75, we derive ∆N eff (applicable at the CMB time) as: Although Eq.S22 was formally derived assuming entropy conservation in a thermalized dark sector, essentially the same result can be reached even if the dark sector is not thermalized: in the latter case we would instead directly apply the redshift of the energy density of dark radiation and take ĝ * = ĝ * ,osc = 1 or 2 depending on the model (ĝ and ĝ * ,osc factors then cancel in Eq.S22).Plugging in g * = 10.75 and the model-dependent T osc and g * ,osc , we can find the value of ∆N eff .For T osc ≳ T EW ∼ 100 GeV, g * ,osc = 106.75, with α ∼ 10%, ξ ∼ 1, we find ∆N eff ∼ O(0.1), readily compatible with existing constraint and potentially detectable with near-future experiments.With α ∼ O(