While the qualitative predictions of inflation are well motivated and well understood [1][2][3][4][5], there is as yet no unique, complete model which connects inflation to the standard model (SM) of particle physics. A crucial component of such a model is the subsequent reheating process [6][7][8][9], which must realize a phase transition from the cold postinflationary state to the hot big bang. In this phase, the inflaton decays into other species to repopulate the Universe and begin the radiation-dominated era.
During an initial stage of preheating, the coherent oscillation of the inflaton field induces explosive production of bosons via parametric or tachyonic resonance (see Refs. [10,11] for reviews). The rapid production of inhomogeneities during this phase generically sources a significant gravitational wave background [12][13][14][15][16][17][18][19][20][21][22][23]. Furthermore, preheating can have important consequences for observable predictions of inflationary models in the cosmic microwave background (CMB): the evolution of the equation of state after inflation affects the mapping of fluctuations' present-day length scales to the times those scales first left the horizon during inflation. In particular, if preheating is sufficiently efficient, then the onset of radiation domination occurs nearly instantaneously after inflation ends.
Most of the early work on preheating focused on models which couple the inflaton to another canonical scalar field [12,[24][25][26][27][28][29][30][31][32][33], with recent studies exploring nonminimally coupled scalar fields [34][35][36][37][38]. Alternatively, the inflaton could couple directly to gauge fields [39][40][41][42][43][44][45][46][47]. One such coupling of particular theoretical and phenomenological interest is that of a (pseudoscalar) inflaton ϕ coupled to the Chern-Simons density F F of a gauge field [48][49][50]. From the theoretical point of view, the (approximate) shift symmetry of a pseudoscalar inflaton (axion) inherently protects its potential from large corrections, ensuring the flatness required for a successful inflationary phase. The phenomenology of these models is extremely rich, with possible observable signatures including the production of sizable non-Gaussianities [51][52][53][54][55], observable gravitational waves [52,54,[56][57][58][59], primordial black holes [42,55,[60][61][62][63][64], μ-distortions [57,65], primordial magnetic fields [50,[66][67][68][69][70], and the generation of the baryon asymmetry [59,[71][72][73][74].
Preheating into gauge fields via a Chern-Simons coupling was first studied within the context of chaotic inflation in Ref. [40], and the model's viability for magnetogenesis was explored in Ref. [41]. Subsequently, we extended this work in Ref. [47] to compute the gravitational wave spectrum produced by the dynamics of gauge preheating. The surprising result was that the (over) production of gravitational radiation provides the strongest probe of (or constraints on) the coupling scale between the axion and gauge fields. As detailed in Ref. [41] (and Sec. II C below), next-generation CMB experiments will probe the radiation content of the Universe to a precision sufficient to rule out most of the interesting region of parameter space in these models. Here and in Ref. [75] we use this result to put important constraints on the axial coupling between the inflaton and gauge fields in a variety of well-motivated inflationary models.
While the chaotic inflation scenario considered in Ref. [47] models the inflaton's coherent oscillations about the minimum of its potential to leading order, this model is disfavored at the 95% confidence level, primarily by constraints on the tensor-to-scalar ratio [5]. To more completely understand the role gauge preheating could play in constraining axion inflation, we investigate the dependence of preheating dynamics and gravitational wave production on the details of the inflationary potential. While the efficiency of preheating and the production of gravitational waves is qualitatively generic, both the energy scale of inflation and the shape of the potential significantly alter the quantitative results as we demonstrate below.
This paper is organized as follows. In Sec. II we review models in which an axion or pseudoscalar inflaton is coupled to an Abelian gauge field and introduce the constraints on stochastic backgrounds of gravitational waves from the CMB. After describing the numerical prescription employed for simulations in Sec. III, we outline some analytic estimates in Sec. IV. We investigate the effect of the energy scale of inflation and the shape of the inflationary potential on preheating and gravitational wave production in Secs. VA and VB, respectively. In Sec. VI we draw our conclusions. In Appendix A we detail the equations of motion of the system and summarize a linear analysis thereof. Appendix B details our procedure for setting accurate initial conditions. Finally, in Appendix C we verify that the Universe remains radiation dominated after preheating to establish the robustness of the gravitational-wave transfer function (for the models we consider).
In this section we introduce the models of interest and briefly review the dynamics of axially coupled gauge fields during and after axion inflation, relegating many details to Appendix A.
We consider the action for a pseudoscalar inflaton ϕ m i n i m a l l yc o u p l e dt og r a v i t ya n dt oa nA b e l i a ng a u g e field A μ ,
where f is a mass scale (the axion decay constant), α is a dimensionless coupling constant, and F μν ≡∂ μ A ν -∂ ν A μ is the field strength tensor of a U(1) gauge theory. 1 We do not identify the gauge field as, e.g., that of standard model hypercharge, nor do we include any charged fields in our model. The dual tensor is defined by
where ϵ μνρσ is the Levi-Civita symbol with convention ϵ 0123 ¼ 1= ffiffiffiffiffiffi -g p .W es e tc ¼ ℏ ¼ k B ¼ 1 and denote by
the reduced Planck mass. The background spacetime is the mostly plus, conformal Friedmann-Lemaître-Robertson-Walker (FLRW) metric,
where the conformal time coordinate τ is related to the cosmic time coordinate via dτ ¼ dt=a. Throughout, primes denote derivatives with respect to conformal time, e.g., a 0 ≡∂a=∂τ, and dots denote derivatives with respect to cosmic time, e.g., _ a ≡∂a=∂t. We use repeated lower indices to indicate contractions with the Kronecker delta.
In this work, we are interested in the effects of varying the shape of the potential VðϕÞ on the efficiency of preheating and the subsequent gravitational wave production. For this purpose we consider various classes of inflationary potentials (the most representative among the ones typically considered by Planck [5]), reporting our parameter choice(s) and the corresponding scalar tilt n s and tensor-to-scalar ratio r (evaluated at a pivot scale which exited the horizon 60 e-folds before the end of inflation). It is worth stressing that in the context of standard inflationary model building most of these models assume the inflaton to be a scalar field. Conversely, for the models described by Eq. ( 1) the inflaton is pseudoscalar. 2 Specifically, we explore the following classes of models:
(1) Chaoticlike monomial models, which have the potential
The extension to different Abelian gauge groups [such as N copies of a U(1) gauge group] is straightforward. Conversely, the extension to non-Abelian gauge groups is nontrivial. See, for example, [76][77][78][79][80][81][82][83][84].
In the context of supergravity this can be achieved by, e.g., using a Kähler potential which is shift symmetric under the imaginary part of a chiral superfield [85]. In order to generate a potential for the inflaton the shift symmetry must be broken, which can be achieved, for example, by coupling the inflaton to a 3-form [86][87][88].
Notice that for n ¼ 2 this simply reduces to chaotic inflation [4], the case we consider here. For n ¼ 2 the Planck normalization of the scalar power spectrum [5] sets m 2 ¼ 1.9 × 10 -11 M 2 pl , with n s ¼ 0.966 and r ¼ 0.13. Another set of models, natural inflation, have potentials of the form VðϕÞ¼V 0 ð1þ cosðϕ=vÞÞ [89]. In this case, with v ¼ ffiffiffiffiffi ffi 8π p M pl the Planck normalization sets V 0 ¼ 5.9 × 10 -10 M 4 pl , for which n s ¼ 0.952 and r ¼ 0.033. However, we find that the results in this case are virtually identical to those for chaotic inflation, so we omit them below.
(2) Starobinsky-like models, which have potentials
pl , in which case n s ¼ 0.969 and r ¼ 0.016. The class of α-attractor models, motivated by supergravity [94,95] (typically divided into E-models and T-models) have exponentially flat potentials whose steepness is controlled by the parameter α, related to the curvature of the Kähler manifold. Since by varying this parameter it is possible to interpolate between chaoticlike models and Starobinsky-like models, we omit the results for this class.
(3) Monodromy inflation, corresponding to the potential
This is part of a broad class of string-theorymotivated models in which large (super-Planckian) field displacement is obtained by wrapping the inflaton trajectory into a series of (sub-Planckian) fundamental circuits [96][97][98]. Note that deep in the inflationary phase (ϕ ≫ ϕ c ) the potential is well approximated by a linear potential (n ¼ 1) of the class in Eq. ( 4).F o rϕ c ¼ M pl =10 the Planck normalization sets μ ¼ 6.0 × 10 -4 M pl , for which n s ¼ 0.975 and r ¼ 0.067. (4) Hilltoplike models, with potentials given by [99] VðϕÞ¼V
where p ≥ 2. We consider a variety of p and v and list the corresponding model fits and predictions in Table I. (5) D-brane models, with potentials of the form [100][101][102][103] VðϕÞ¼V
The equations of motion that result from the variation of Eq. ( 1) are difficult to solve analytically due to the nonlinear interactions of the axion and gauge field. Furthermore, the homogeneous axion background sources a tachyonic instability in the gauge field which quickly renders a linear analysis invalid. In this subsection, we briefly outline the well-known linear treatment of the gauge-field fluctuations in the inflationary phase and during early stages of resonance. Full details may be found in Appendix A.
In a homogeneous axion background, the helical polarizations of the gauge field obey [40] A 00 AE ðkÞþkðk ∓ 2HξÞA AE ðkÞ¼0; ð9Þ
where we define the instability parameter,
and the conformal Hubble parameter, H ¼ a 0 =a. The interaction between the gauge field and the homogeneous, rolling axion leads to the exponential growth of one of the polarizations of the gauge field, understood as an imaginary ("tachyonic") effective frequency ω ∼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi kðk -2HξÞ p (assuming ξ > 0) for modes k<2Hξ. The polarizations are amplified (relative to the conformally invariant radiation solution A rad;AE ) by a factor of
The amplification is controlled by the parameter ξ, which depends on the inflaton velocity, ϕ 0 =H. The largest effects therefore occur near the end of inflation and during reheating when the inflaton velocity is largest. Modes with wave numbers in the band 1 ≲ k=H ≲ 2ξ are significantly enhanced by the axion. Shortly after their production during reheating these modes can rescatter, generating the subdominant helicity as well as inflaton quanta (see Fig. 2). We use numerical simulations to capture these nonlinear dynamics in later sections.
The tachyonic production of gauge-field modes leads to large anisotropic stresses which results in the copious production of gravitational waves [47]. These gravitational waves form a stochastic background with frequencies near the Hubble scale during reheating, corresponding to f ∼ 10 9 Hz today 4 for inflation occurring near the scale of grand unification, Λ ∼ 10 16 GeV.
Gravitational waves with wavelengths that are much shorter than the Hubble scale contribute energy density that gravitates like radiation. During preheating, significant energy density can be deposited into subhorizon gravitational waves which, as in Ref. [47], may be constrained by CMB and big bang nucelosynthesis measurements of the total radiation density in species beyond the standard model, ΔN eff ¼ N eff -3.046.
The net energy density in gravitational waves is given by
By conservatively assuming that all extra radiation density present during the formation of the CMB (beyond the standard model) is composed of gravitational waves, a bound on N eff directly constrains Ω gw;0 h 2 via [104]
where the present energy density in photons is Ω γ;0 h 2 ¼ 2.47 × 10 -5 . From the Planck limit, jΔN eff j ≲ 0.33 [105,106], we obtain a bound Ω gw;0 h 2 ≲ 1.85 × 10 -6 . Next-generation CMB experiments (e.g., CMB-S4) will probe ΔN eff ≤ 0.03 at 1σ and ΔN eff ≤ 0.06 at 2σ [107], improving the upper limit by an order of magnitude to
A separate analysis [108] using Planck data provides a stronger constraint, Ω gw;0 h 2 < 1.2 × 10 -6 at 95% confidence. The same study projected that next-generation satellite missions COrE [109] and Euclid [110] will place 2σ bounds of ΔN eff < 0.013 or Ω gw;0 h 2 < 7.6 × 10 -8 . In the results below, we explore the relationship between gravitational wave production at these levels and the preheating process itself.
The first direct detection of gravitational waves by LIGO [111] prompts consideration of direct gravitational wave detectors as alternative probes of stochastic gravitational wave backgrounds. In particular, the presently operating LIGO [112] and VIRGO [113] detectors (with KAGRA [114] expected to join this network of detectors in the near future), as well as future missions like LISA [115] and the Einstein Telescope [116], could detect stochastic gravitational wave backgrounds generated in the very early universe. In order to detect a stochastic gravitational wave background with a direct gravitational wave detector, the corresponding signal must have a sufficiently large signalto-noise ratio (SNR), given by [104,117]
where Ω s is used to denote the detector sensitivity (expressed as a cosmological parameter), S n is the strain sensitivity of the instrument, T denotes the observation time, and f min and f max are, respectively, the minimal and maximal frequencies to which the instrument is sensitive. The fact that Ω s is proportional to f 3 (times the strain) makes direct detection of high frequency gravitational waves difficult; however, if inflation occurred at a lower scale, therefore producing a lower-frequency signal, direct detection is possible [14,15]. At the same time, current methods (see, for example, [118][119][120]) for the direct detection of high-frequency gravitational waves are unfortunately not expected to be sensitive enough to place interesting bounds on stochastic backgrounds of cosmological gravitational waves.
Our numerical approach is very similar to that employed in Ref. [47]. We numerically integrate the classical equations of motion of the inflaton, Eq. (A22), and the gauge fields, Eq. (A19), in an FLRW background governed by Eqs. (A1) and (A2). Specifically, we discretize the evolution equations onto a three-dimensional (3D), regularly spaced grid with periodic boundary conditions and step the coupled system of equations through time using the fourthorder Runge-Kutta method. Spacetime expansion is implemented self-consistently by computing the current energy density and pressure, Eqs. (A3) and (A4), averaged over the simulation volume.
Numerical treatments of gauge fields, compared to those for scalar fields, require particular care. Specifically, numerical evolutions of gauge fields must be stable with respect to their constraints. These constraints are Gauss's law, Eq. (A20), and the prescribed gauge condition (in our case Lorenz gauge, ∂ μ A μ ¼ 0).
One numerical method for evolving gauge theories is lattice gauge theory, which recasts the gauge fields as a system of "link variables" representing the discrete connection between adjacent lattice sites. The advantage of this method is that discrete gauge invariance is an exact symmetry of the discrete system. The system of dynamical equations of motion that results from the direct variation of the discrete lattice action yields an evolution scheme which preserves the gauge condition and (degree of satisfaction of) Gauss's law. However, as noted by [46,121], an appropriately gauge-invariant representation of the axial coupling term Fμν F μν renders the equations of motion implicit, requiring a computationally expensive iterative solution technique.
We take an alternative approach and evolve the Euler-Lagrange equations of the continuum theory for all of the components of the gauge field A μ in Lorenz gauge. That is, rather than evolve the evolution equations of a discretized theory, we numerically integrate the equations of motion of the continuum theory by discretizing the dynamical equations themselves (as often employed for scalar fields). While gauge invariance is not exact in this case, the critical property required for robust results is stability. One can recast the equations of motion such that constraint violations, while dynamical, remain bounded in time [122,123]. Doing so requires evolving additional, redundant degrees of freedom, which in our case means evolving all four components of the gauge potential,
obtained by applying the Lorenz gauge choice, ∂ μ A μ ¼ 0, to Eq. (A18). The results of [39,40] demonstrate the stability of the satisfaction of the gauge constraint under both the axial coupling considered here and a dilatonic coupling to the gauge-field kinetic term. Recent work used lattice-gauge-theory-based simulations employing an iterative scheme [46], the results of which reproduce those of Ref. [40] (and those presented in Sec. VA below). In addition, in Ref. [45] we compared simulations of preheating into three U(1) gauge fields using the approach detailed above as well as preheating into a set of SU(2) gauge fields (with the gauge self-coupling tuned so that internal interactions are negligible) using lattice gauge theory, finding near-perfect agreement between the two methods.
There are two differences in our numerical implementation relative to that in Ref. [47]. First, we have improved our procedure for obtaining the power spectra of the gauge fields at the end of inflation. To capture the tachyonic enhancement of (one polarization of) the gauge fields during inflation, we numerically evolve the background spacetime, the homogeneous mode of the inflaton, and the linearized equations of motion of the gauge-field polarizations. In contrast to Refs. [40,47], here we more carefully account for the backreaction of the gauge fields onto the background quantities. At the range of couplings we explore, the approximations used to derive the analytic results Eqs. (A35) and (A36) (that jξj ≳ 4 and is constant) are inaccurate during the final few e-folds of inflation. While these approximations are valid earlier during inflation and at larger couplings [the regimes for which Eqs. (A35) and (A36) were developed [59,124] ], they overestimate the backreaction of the gauge fields onto the background at the end of inflation. At the largest couplings we consider here (e.g., α=f ∼ 15 M -1 pl for chaotic inflation), these errors artificially offset the end of inflation by ∼0.5 e-folds. As a result, the use of Eqs. (A35) and (A36) to compute the duration of inflation contributes errors to the relationship between present-day length scales and the number of e-folds before the end of inflation when those scales exited the horizon. To mitigate these errors, we perform the integrals Eqs. (A32) to (A34) via quadrature of the numerically integrated gauge-field fluctuations. See Appendix B for full details of this procedure.
Other than improving the gauge-field power spectra, our procedure for generating initial conditions is unchanged from Ref. [47]. We initialize the lattice simulations one to two e-folds (depending on the model; see Table II) before the end of inflation, where the inflaton is initialized with a mean value and velocity given by the background evolution. We seed the fluctuations in each field and its (conformal-time) velocity as independent, Gaussianrandom fields by drawing, for each Fourier mode on the lattice, an amplitude from a Rayleigh distribution and a uniform random phase. The variance of the Rayleigh distribution for each mode is set according to the power spectrum obtained from the numerical solution of the linearized equation of motion for the field fluctuations evolved through inflation. Finally, the chosen Lorenz gauge condition is satisfied on the initial slice by setting A 0 0 ¼ 0 and projecting the polarization fields, A AE ðkÞ, onto their vector components A i ðkÞ, which automatically yields
The second difference in our numerical implementation is that, in contrast to the pseudospectral solver used in Ref. [47], we represent spatial derivatives with finite differencing, using fourth-order-accurate centered difference stencils. In particular, whereas in Ref. [47] we evolved the equations of motion for the tensor metric fluctuations h ij in Fourier space [computing and Fourier transforming the source term Eq. (A12) and projecting with Eq. (A10) at each time step], we here evolve the equation of motion in position space without applying the transverse-traceless projection [125],
Instead, we apply the transverse-traceless projection to u ij to obtain
only when computing the gravitational wave spectrum Ω gw ðkÞ via Eq. (A14). The advantage of this procedure is that it requires no Fast Fourier Transforms at each time step, which scale poorly to distributed-memory systems.
We have verified that the axion-gauge-field dynamics and the gravitational wave spectra Ω gw ðkÞ reproduce the results of the pseudospectral method extremely well. The software we developed for this purpose is PYSTELLA,foot_2 a PYTHON-based, MPI-parallel, and GPUaccelerated code making use of PyOpenCL [126] and Loo.py [127] for the generation of OpenCL code to run on GPUs. Because PYSTELLA may thus run on multiple GPUs (or any architecture with OpenCL support), we evolve lattices with N 3 ¼ 384 3 points, enabling reliable simulations of larger couplings α=f than in Ref. [47]; we also checked a representative sample of our results against simulations with 512 3 points to show convergence. The Friedmann constraint, Eq. (A1), is satisfied to a precision of Oð10 -3 Þ or better for all simulations we present.
The inflationary models we consider exhibit different postinflationary Hubble scales relative to the oscillation timescale about the minima of their respective potentials, m ϕ (defined by m 2 ϕ ¼ ∂ 2 V=∂ϕ 2 evaluated at the minimum of the potential). For this reason, we choose different comoving box lengths L for different models. In all cases, we have checked our results are insensitive to the precise choice of L, i.e., that the simulations have converged. We also tune the initialization time (relative to the end of inflation) based on the model, starting sufficiently early to capture any nonlinear effects at the end of inflation (but no earlier than necessary to save computational expense). We tabulate these choices, as well as other relevant parameters for each model, in Table II. We use a time step Δτ ¼ Δx=10 ¼ L=N=10 in all cases and set the scale factor a ¼ 1 at the end of inflation, when ϵ H ≡ -_ H=H 2 ¼ 1.
To establish expectations for the effects of the inflationary model on preheating and gravitational wave production, we begin with some analytical estimates. While the linear theory reviewed in Appendix A is not valid during preheating, it can be used to gain intuition for the scaling of the backreaction on the expansion rate and on the motion of the inflaton during the initial phase of preheating which is typically the most violent. We additionally use a "rule of thumb" for stochastic gravitational wave production from cosmological processes [128] to make contact between the characteristics of the inflationary potential and the efficiency of gravitational wave production from preheating.
From the Friedmann equation, Eq. (A1), by making use of the approximation for the energy density in the gauge fields, Eq. (A35), we can derive the approximation TABLE II. The specific parameters chosen for each inflationary model under consideration: the effective inflaton mass, the simulation box length, the simulation start time in terms of the number of e-folds relative to the end of inflation, the Hubble parameter at the end of inflation H e , the ratio of the lattice's infrared cutoff to the comoving Hubble scale at the end of inflation, equal to ð2π=LÞ=H e , and the energy scale at the end of inflation.
where we introduce the slow-roll parameters
In standard single-field slow-roll inflation, the Einstein equations imply that ϵ ϕ ¼ ϵ H . However, in the regime of strong gauge-field backreaction, this relation does not hold. Equation ( 20) quantifies the contribution of the gauge fields to the Hubble rate relative to the contribution from the inflaton. As such, it provides a rough estimate of the value of ξ (and therefore ϕ 0 =H, once α=f is fixed) required in order for the gauge-field contribution to the energy density to be comparable with the inflaton's.
Beyond modifying the value of the Hubble parameter as in Eq. ( 20), the gauge fields induce a new friction term in the equation of motion of the inflaton. Substituting the spatially averaged interaction, Eq. (A35), into the equation of motion for the inflaton, Eq. (A22), while making use of the lowest order approximation of Eq. ( 20) and assuming ϕ 00 =H 2 ≪ M pl , we can derive the usual slow-roll condition, now including the effect of the gauge-field backreaction,
Note that the gauge-field backreaction changes the relationship between the slope of the potential relative to its magnitude and so also the slow-roll parameter in Eq. (21). In this scenario, gauge-field-induced friction provides an additional mechanism to enforce the flatness of the potential required for slow-roll inflation [52,57,124].
There are two effects at play in Eq. ( 22). The first is that the backreaction changes the Hubble rate, increasing Hubble friction (on the inflaton), while the second is the backreaction of the E • B term. The impact of these two effects can be understood by comparing Eq. ( 22) with Eq. (20). As the same (V-dependent) term is compared with unity (for the backreaction on H) and with ffiffiffiffiffiffiffi 2ϵ ϕ p (for the backreaction on the equation of motion of the inflaton), respectively, the gauge-field-induced modification of the Hubble parameter is higher order (in V=M 4 pl ) than the new friction term in the equation of motion of the inflaton. In addition, while both the gradient of the potential and the (lowest order expression of the) Hubble friction term in Eq. ( 22) do not depend explicitly on the scale of V, the gauge-field friction depends linearly on the inflationary energy scale. As a consequence, low-scale models are expected to require a larger value of ξ in order for the gauge fields to become important.
Finally, while V is nearly constant (which is the case during inflation), the gauge-field-induced backreaction grows exponentially with ξ.H o w e v e r ,i fV is not constant (as is the case during preheating where it decreases as fast as, or faster than, e -2πξ ) then the gauge-field-induced backreaction may either never become relevant or be shut off. Such a possibility is realized in potentials without minima, such as a tanhðϕ=MÞ shape. We leave investigations of this class of potentials to future work.
To set expectations for the range of couplings α=f relevant for preheating in a given model, we make use of Eq. ( 20) (to leading order in slow roll) to obtain scaling relations between α=f and the Hubble scale at the onset of preheating (and so the parameters of the model). We quantify efficiency in terms of the fraction of the Universe's energy residing in the gauge fields,
In the regime of this approximation's validity, we may ask how much one must tune ξ (or α=f) to compensate for a reduction in H=M pl in order to keep ρ gauge =ρ fixed.
Considering two sets of parameters ðH 1 ; ξ 1 Þ and ðH 2 ; ξ 2 Þ, this amounts to solving the nonlinear equation
for ξ 2 in terms of H 1 =H 2 and ξ 1 . When ξ 1 and ξ 2 are both small (and so ξ 1 -ξ 2 is as well), we solve
From this we read that the coupling must increase by an additive amount proportional to the logarithm of the ratio of Hubble scales.
In our analysis below we demonstrate that this scaling relation is sufficiently accurate for our estimates.
References [40,41,47] demonstrate that in chaotic inflation models with m ϕ ≈ 10 -6 , the threshold value of α=f for which preheating is complete (i.e., max ρ gauge =ρ ∼ 80%)is 9M -1 pl . With this as a baseline, we apply Eq. ( 24) to estimate how this threshold scales as we change the energy scale of inflation (and, later, the inflaton model itself). Consider the simplest case of a quadratic potential, an approximation useful for all of our models during the preheating phase,
In this case, the Hubble parameter scales as H ∼ m ϕ , since at the end of inflation
In terms of the coupling, Eq. ( 26) tells us that, for some other mass m ϕ;2 , we require
for comparable preheating efficiency. Since we are interested in the efficiency of preheating (rather than the inflationary production of gauge bosons), in this expression we evaluate ∂ϕ=∂N at the end of inflation, N ¼ 0. At the end of chaotic inflation, ∂ϕ=∂N ≈ 1.4M pl , in which case Eq. ( 29) reduces to
Thus, for comparable preheating efficiency we expect to need to increase the axion-gauge coupling α=f by ∼M -1 pl for each order of magnitude we reduce m ϕ .
B. Gravitational wave production and the "rule of thumb"
While regimes of highly efficient preheating exist for all models (as we show in Sec. VB), the structure of resonance varies from model to model. To understand the relationship between the scales at which gravitational waves are produced and their resulting amplitude, we use the "rule of thumb" developed by the authors of Ref. [128] to estimate the stochastic gravitational wave production from cosmological processes. By approximating the source as a Gaussian of width σ peaked at wave number k à , Ref. [128] estimates the peak amplitude as
where α is the fraction of the energy in the gravitational wave source relative to the Universe's total energy density at that time, β encodes the anisotropy of the source, and w is the equation of state of the Universe at that time.
Observe that Ω gw;0 h 2 decreases quadratically with the ratio of the peak wave number k à to the Hubble parameter at that time, H à . The oscillation frequency of the inflaton background (i.e., its effective mass) sets the scales of efficient resonance, while ξ ∼ α=f controls how far inside the horizon tachyonic resonance occurs. The former depends on the shape of the potential, while the latter can depend on the inflationary scale via the arguments of Sec. IVA: lower inflationary scales require larger couplings for complete preheating. To quantify the effect of the potential shape, note that the effective inflaton mass, defined by m 2 ϕ ≡∂ 2 V=∂ϕ 2 evaluated at the minimum of the potential, sets the scales of interest during preheating. As such,
, the shape of the potential (about its minimum) determines the scales at which preheating occurs, which in turn affects the size of the gravitational wave signal according to Eq. ( 31).
In this section, we present simulations detailing the effect of the energy scale of inflation on the efficiency of preheating and the subsequent generation of gravitational waves. For this analysis we fix a chaotic inflationary potential while tuning m ϕ to study different energy scales of inflation (i.e., for the time being, we ignore the fact that m ϕ should be chosen to fit the normalization of the scalar power spectrum during inflation). This simplification allows us to separate the effect of the scale of the potential from that of its shape, which we consider in Sec. VB. Further, preheating studies are often restricted to chaotic inflationary models, a choice justified because preheating probes the inflaton's oscillation about the minimum of its potential, which is quadratic to leading order.
We first seek to determine how the axion-gauge coupling α=f must be tuned in order to achieve complete preheating as we lower the inflationary scale and to evaluate the accuracy of the analytic estimates made in Sec. IVA.I n Fig. 1 we depict the relationship between the efficiency of preheating and the net gravitational wave production across five decades of m ϕ . Because larger couplings result in power transfer to modes with larger momentum, for values of m ϕ < 10 -6 M pl we use a box length L ¼ 7.5m -1 ϕ in order to ensure sufficient short-wavelength resolution in all cases. We immediately observe that Eq. ( 30) does indeed accurately predict the range of couplings for which preheating becomes efficient and even complete (ρ gauge =ρ ≳ 80%). It is reassuring that this scaling argument-derived from linearized approximations-is applicable and that it breaks down in the large-coupling limit when the approximations are least valid. To visualize this scaling, in the right panel of Fig. 1 we plot the preheating efficiency and net gravitational wave production as functions of an "adjusted" coupling, i.e., the value of α=f that would be required to obtain the same efficiency were m ϕ ¼ 6.16 × 10 -6 M pl , as predicted by the analytic estimate Eq. ( 29).
The efficiency of gauge-field production becomes more complicated at larger couplings where nonlinear processes become important. After a regime of α=f in which preheating remains comparably efficient (which, for m ϕ ¼ 6.16 × 10 -6 M pl , corresponds to 9M -1 pl ≲ α=f ≲ 10.4 M -1 pl ), ρ gauge =ρ begins to decrease. At even stronger coupling, energy transfer to the gauge fields continues, gradually increasing ρ gauge =ρ to 100%. Regardless, the inflaton condensate is totally depleted, as shown by the dotted lines in Fig. 1 which indicate the (maximum) fraction of energy in either gauge fields or fluctuations of the axion. To examine the physics that realizes this trend, consider comparing the energy in inflaton fluctuations,
to the energy in the gauge fields, ρ gauge . In Fig. 2 we plot the ratio of these two quantities, ρ gauge =ρ δϕ , as well as ρ δϕ =ρ ϕ which measures the degree to which the axion has fragmented. In the top row, we observe that as preheating approaches ∼80% efficiency (at α=f ∼ 9 M -1 pl ), the inflaton also becomes more fragmented. Since the axion's equation of motion is linear for a quadratic potential, any inflaton particle production must occur through backscattering from produced gauge bosons, an inherently nonlinear process. As the coupling increases to 9 M -1 pl ≲ α=f ≲ 10.2 M -1 pl , the inflaton becomes totally fragmented and the energy in FIG. 1. Left panels: Preheating efficiency, quantified by the maximum ρ gauge =ρ over the simulation (top panels), and the total fractional energy in gravitational waves today, Ω gw;0 h 2 (bottom panel), as functions of axion-gauge coupling α=f. Each curve fixes a quadratic potential and corresponds to various values of m ϕ , decreasing by factors of 10 from m ϕ ¼ 6.16 × 10 -6 M pl (red curves) to m ϕ ¼ 6.16 × 10 -10 M pl (purple curves). The dashed lines (which follow the same color scheme) correspond to the maximum total energy in fluctuations, i.e., the maximum amount of energy in the simulation outside of the inflaton condensate. Lines indicating the current ΔN eff bounds on Ω gw;0 h 2 from Planck and CMB-S4 from Ref. [108] are plotted in solid and dashed black curves, respectively, while the region between CMB-S4's 1σ and 2σ projections [107] is shaded grey. Right panels: Same, but with the horizontal axis adjusted according to Eq. (29). That is, the horizontal axis corresponds to the actual coupling α=f shifted to the coupling that would be required for the same preheating efficiency for m ϕ ¼ 6.16 × 10 -6 M pl (the value which fits the amplitude of the scalar power spectrum), as predicted by the analytic estimate Eq. (29).
gauge fields relative to inflaton fluctuations does not increase. Preheating thus remains comparably efficient in this coupling regime. Evident in the middle row is that, as α=f is increased past 10.2 M -1 pl and inflaton fragmentation occurs more rapidly, the energy in fluctuations of the axion approach roughly half that of the gauge fields. The exponential amplification of the gauge fields likewise ends earlier and earlier, as the tachyonic resonance is driven by the motion of the (now-depleted) inflaton condensate. Thus, preheating becomes slightly less efficient, as apparent in Fig. 1.
Increasing α=f past 11.6 M -1 pl , we see that while the initial phase of resonance continues to end earlier and earlier, the inflaton also becomes less fragmented during this phase. The backreaction of the gauge fields onto the axion's background dynamics now exerts a dramatic amount of friction before the axion even first crosses through the minimum of its potential. Tachyonic amplification of the gauge fields recommences and proceeds at a slower and slower rate as the coupling increases (as this increases the gauge-field friction). This process terminates later and later at values of ρ gauge =ρ increasingly close FIG. 2. Left: Ratio of the energy in gauge fields, ρ gauge , to the energy in inflaton fluctuations, ρ δϕ [defined in Eq. ( 32)], as a function of e-folds after inflation, N ¼ ln a. Right: Fragmentation of the inflaton, measured by ρ δϕ =ρ ϕ . Each row depicts a different coupling regime, with couplings α=f labeled in the legend of the left column (in units of M -1 pl ). The displayed simulations fix a chaotic inflationary potential with m ϕ ¼ 6.16 × 10 -6 M pl to fit the amplitude of the scalar power spectrum.
to 100%, as shown in Fig. 1, at which point the inflaton condensate is totally depleted. Indeed, preheating is 100% efficient in these cases because the gauge-field friction is strong enough to ensure that tachyonic amplification continues until preheating is complete (in contrast to the axion's oscillations allowing for backscattering effects to become important). Note that these results (which reproduce those originally obtained in Ref. [40]) are very similar to those presented more recently in Ref. [46] for the case of chaotic inflation, and our analyses of these regimes are likewise similar.
The inflaton halts higher and higher up its potential as the gauge-field friction becomes more important, and eventually accelerated expansion (w<-1=3) recommences (for α=f ≳ 13.2 M -1 pl ). As observed in Ref. [41], the inflaton is momentarily "trapped," as depicted in Fig. 3. Its vanishing velocity momentarily shuts off the tachyonic instability, so that the gauge fields redshift and the inflaton starts to roll again (restarting the resonant enhancement of the gauge fields). At the highest couplings we simulate here (α=f ∼ 14.8 M -1 pl ), the inflaton is trapped before the background stops accelerating (i.e., before w ≥-1=3), resembling models of inflation which achieve slow roll via the same axion-gauge-field coupling [124].
Regardless of the value of m ϕ , each set of simulations exhibits the same qualitative features as the coupling increases. However, as m ϕ is decreased, the range of couplings α=f spanned by these features broadens in a manner that cannot be described by the analytic estimates above. The backscattering effects we posit as responsible for the dip in efficiency depend on the amplitude of axion fluctuations, which scales linearly with m ϕ . Thus, in order for backscattering to be important to a comparable degree, the coupling has to be increased by a greater relative amount to compensate for the smaller amplitude of vacuum fluctuations at lower values of m ϕ . From the analytical estimates presented at the beginning of this section, we expect the backreaction of the gauge fields onto the axion condensate to scale with the energy scale of the inflationary potential. Thus, models with lower m ϕ require a larger relative increase in α=f to enter the regime of slow tachyonic resonance due to the gauge fields' friction on the axion.
The approximate universality we observe also extends to the gravitational waves produced by the dynamics of preheating. Of note in Fig. 1 is that any scenario which exhibits (near-)complete preheating (namely, ρ gauge =ρ ∼ 80%) results in a net gravitational wave production that could be detected (or ruled out) by CMB-S4, as discussed in Sec. II C. For the strongest couplings we simulate, the total integrated Ω gw;0 h 2 exceeds 10 -6 , which is already ruled out by Planck data [108]. To quantify this claim, in Fig. 4 we scatterplot the fractional energy density in gravitational waves today versus the efficiency of preheating. The cluster of data points around and above max ρ gauge =ρ ∼ 0.8 mostly lies above the CMB-S4 bound, indicating that next-generation experiments will place bounds on the axion-gauge-field coupling α=f, regardless of m ϕ .
In the upper panel of Fig. 4 we also observe that gravitational wave production is moderately less effective at lower inflationary scales. Returning to the arguments of Sec. IV B, the larger couplings required for comparable preheating at lower m ϕ push the resonance further inside the horizon. As the momenta undergoing the tachyonic instability are k=aH < ξ ∼ α=f, from the rule of thumb, Eq. ( 31), we expect Ω gw to decrease with ðk à =HÞ 2 ∼ ðα=fÞ 2 . To verify this scaling, in the lower panel of Fig. 4 we scale Ω gw;0 h 2 by the squared ratio of each simulation's coupling α=f to its "adjusted" value from Eq. (29). With this multiplicative factor, the trends of Ω gw;0 h 2 versus max ρ gauge =ρ line up closely, confirming our hypothesis.
The shape of the gravitational wave spectra (which is fairly generic and qualitatively comparable to that produced by tachyonic resonances) is similar for all values of m ϕ (when comparing couplings which yield comparable preheating efficiency). Plotting the signals that would be observed today in Fig. 5 present-day frequencies at which these signals would be observed. 6In general, Ω gw;0 h 2 increases exponentially with α=f, which is to be expected as the tachyonic resonance is exponential in ξ ∝ α=f. The actual rate at which Ω gw;0 h 2 increases as a function of expðα=fÞ falls into two regimes: the range of couplings where the actual efficiency of preheating increases exponentially with α=f, and those for which preheating is always complete (the latter of which is slower than the former). In the first regime the gravitational wave source [parametrized by α in Eq. ( 31)]i s growing exponentially. In the second, while the simulations all completely transition to radiation domination, gravitational waves are continually sourced by the second phase of slow tachyonic resonance at an efficiency which still increases with the coupling strength.
Last, we note that, despite the complicated trend of ρ gauge =ρ as α=f increases, the inflaton condensate always ends up depleted as the coupling increases past the critical value where efficient preheating is first achieved (e.g., α=f ∼ 9.6 M -1 pl for m ϕ ¼ 6.16 × 10 -6 M pl ). Thus, in these cases the end state of the simulations is always radiation domination. However, the proportion of that radiation composed of axion fluctuations varies with coupling, as depicted in Fig. 2. Ignoring any decay channels, the axion's fluctuations redshift after preheating until their physical momenta drop below their mass, at which point they become nonrelativistic. From this point on the axion energy density redshifts as matter. If the axion's lifetime is sufficiently long, its energy eventually dominates over the gauge fields (which, being radiationlike, decay faster than matter). Any deviation from an equation of state of radiation, w ≡ p=ρ ¼ 1=3, suppresses the gravitational wave density observed today, Ω gw;0 h 2 , relative to what the transfer function Eq. (A16) accounts for (which assumes the Universe was radiation dominated from the time of emission until matter-radiation equality). In Appendix C we demonstrate that Bose enhancement resulting from the larger occupation numbers from preheating ensures that perturbative decays happen quickly enough that the Universe remains radiation dominated. As such, we expect little to no suppression of Ω gw;0 h 2 relative to the values we report.
We now explore the dependence of our results-in particular the amplitude of the resulting gravitational wave spectrum-on the shape of the potential during the reheating phase. We simulate preheating in the inflationary models detailed in Sec. II A and discuss the extent to which the results of Sec. VA are modified.
In Fig. 6, we plot the efficiency of preheating and corresponding gravitational wave production over a range of axion-gauge-field couplings α=f. We first observe that the relationship between max ρ gauge =ρ and α=f follows the same general trend as presented in Sec. VA. That is, once α=f is large enough for preheating to be efficient, max ρ gauge =ρ remains roughly 90% until backscattering effects become important, at which point we observe a dip in efficiency. Finally, even larger couplings lead again to a regime of strong backreaction leading to slow tachyonic resonance, resulting in near-completely efficient preheating. Turning to the lower panel of Fig. 6, we see that (i) depending on the coupling strength, preheating in all models can yield a net production of gravitational waves that would be probed by CMB-S4 measurements of N eff , and (ii) models with r ≳ 10 -2 are already constrained by Planck data. While all models we study here reach a regime of gravitational wave production detectable by future experiments such as CMB-S4, for those with r ≳ 10 -2 the entire regime of efficient preheating could be ruled out by a null detection of ΔN eff . FIG. 4. The total gravitational wave energy today versus the final preheating efficiency ρ gauge =ρ, plotted for a chaotic potential with varying m ϕ . The bottom panel depicts the same data, but scaled by the square of the ratio of α=f (for a given simulation) to the coupling required for equal preheating efficiency for m ϕ ¼ 6.16 × 10 -6 M pl , as estimated by Eq. ( 29).
In Fig. 7 we display the energy density in gravitational waves today, Ω gw;0 h 2 , as a function of the efficiency of gauge preheating, which, as above, we quantify by the maximum fraction of energy in gauge-field fluctuations during the simulation. These results exhibit the sensitivity of Ω gw;0 h 2 to the details of the potential. That is, while the general relationship between gravitational wave production and preheating efficiency follows a similar trend, the overall scaling of Ω gw;0 h 2 differs from model to model as observed in Sec. VA. Again, in general the larger coupling α=f required for complete preheating for a particular model, the less efficient the subsequent gravitational wave production. Referring to Table II, low-scale models require larger coupling for efficient preheating and present correspondingly weaker detection prospects. This observation, together with Figs. 6 and7, leads us to make the broad claim that for models of inflation with tensor-toscalar ratios observable by CMB-S4 experiments, preheating into gauge fields could be simultaneously probed via the contribution of gravitational waves to ΔN eff .
A positive detection of r ≳ 10 -3 together with ΔN eff by CMB-S4 experiments could provide evidence for a pseudoscalar inflaton (axion) reheating the Universe through preheating to gauge fields. Alternatively, detection of r with a measurement of ΔN eff consistent with zero would provide stringent bounds on the axion-gauge coupling α=f-in particular, the regime in which preheating is the sole mechanism by which the Universe was reheated would be ruled out. In the former case, a precise prediction of the end of inflation (relative to the pivot scale used to parametrize CMB observables) may require the methods of Appendix B to accurately model the backreaction of gauge fields onto the inflationary background toward the end of inflation.
To evaluate the claim that the efficiency of gravitational wave production is correlated with r, we vary the free FIG. 5. Gravitational wave spectrum observed today for various inflationary energy scales. Each curve corresponds to a different simulation of preheating in chaotic inflation, with coupling α=f and inflaton mass m ϕ labeled in the legend. In each case, the total amount of gravitational wave production Ω gw;0 ∼ 10 -7 . FIG. 6. Preheating efficiency, quantified by the maximum ρ gauge =ρ over the simulation (top panel), and the total fractional energy in gravitational waves today, Ω gw;0 h 2 (bottom panel), as functions of axion-gauge coupling α=f. Each color denotes a different inflationary potential, as labeled in the legend. The dashed lines (which follow the same color scheme) correspond to the maximum total energy in fluctuations, i.e., the maximum amount of energy in the simulation outside of the inflaton condensate. Lines indicating the current ΔN eff bounds on Ω gw;0 h 2 from Planck and CMB-S4 from Ref. [108] are plotted in solid and dashed black, respectively, while the region between CMB-S4's 1σ and 2σ projections [107] is shaded grey. Note that, to save clutter, we only plot a subset of the hilltop parameter points listed in Table I.
parameter v of the hilltop model, which allows us to tune the flatness of the potential during inflation, and so both the energy scale of inflation and r. In Fig. 8 we observe the exact trend we noted in our survey of different potentials: while the relationship between max ρ gauge =ρ and Ω gw;0 h 2 is clearly universal, as v increases, gravitational wave production is more efficient (as is preheating, which requires increasingly lower α=f to be efficient). Consulting Table II, increasing v decreases k à =H à , which (referring to Table I) correlates to larger tensor-to-scalar ratios, r. From Fig. 8 we extrapolate that, regardless of the choice of v, at sufficiently large coupling the hilltop models will exhibit complete preheating (entering the regime of slow tachyonic resonance sustained by gauge-field backreaction). However, the resulting gravitational wave production will be suppressed, reducing the constraining power of ΔN eff in this scenario. We note, however, that existing constraints on the instability parameter ξ during inflation (the strongest of which are due to primordial black hole production [55,[60][61][62][63][64]) still provide significantly weaker constraints on α=f for small-r (and so small-ϵ H ) models, as ∂ϕ=∂N ∼ ffiffiffiffiffi ffi ϵ H p . At the same time, for larger values of v there is a regime of detectable gravitational waves from preheating (while still r ≲ 10 -3 , below the target for CMB-S4). As such, a detection of nonzero ΔN eff but not nonzero r would be consistent with axion inflation in models with r ≲ 10 -3 . Finally, we note that we have repeated this analysis for the D-brane model, finding the same relationship as for the hilltop model.
In this paper, we have studied the effect of the inflationary potential on the dynamics of gauge preheating and gravitational wave production after axion inflation. Abelian gauge fields coupled via ϕF F are amplified via a tachyonic resonance from the postinflationary axion condensate, which rapidly and efficiently reheats the Universe and copiously produces a stochastic background of gravitational waves. For strong enough coupling, these signals, first computed in Ref. [47], contribute so greatly to the radiation content of the Universe that Planck data already places bounds on the axion-gauge coupling α=f, which CMB-S4 will improve upon significantly. Here we have extended this study to consider a variety of inflationary models, finding that such a level of gravitational wave production typically occurs in scenarios where preheating completely reheats the Universe.
While the results are qualitatively similar across models, the efficiency of preheating (at a given coupling α=f) varies depending on both the energy scale of inflation and the shape of the potential. The scale of the potential (as well as the coupling parameter α=f) controls the significance of backreaction of the gauge fields onto the axion's background dynamics. At strong enough couplings, this gaugefield friction supports an extended phase of slow tachyonic resonance which persists until 100% of the axion's energy has transferred to the gauge fields. As such, lower-scale inflation models require larger coupling to enter into this regime, which suppresses the resulting gravitational wave amplitude by moving the source further inside the horizon.
The tensor-to-scalar ratio is sensitive to the flatness of the potential at CMB, measured by ϵ H , and (once the amplitude of scalar fluctuations is set to match CMB data) to the energy scale of inflation. As a result, gravitational wave production is comparatively less efficient in models with smaller tensor-to-scalar ratios. However, when preheating completely transitions the Universe to radiation domination, generally the level of gravitational wave production contributes to the effective number of radiationlike degrees of freedom N eff at a level that could be detected (or ruled out) by next-generation CMB experiments. Should CMB-S4 (which targets r ≳ 10 -3 ) detect nonzero r, a simultaneous detection of nonzero ΔN eff could be an indication that preheating occurred via strong coupling to gauge fields. Conversely, the lack of observed gravitational waves via ΔN eff would place severe constraints on α=f, in particular ruling out most of the parameter space in which preheating in this model is solely responsible for the transition to radiation domination. Notably, for these models CMB-S4 could rule out the entire regime of couplings for which preheating alone reheated the Universe.
Planck data currently set α=f ≲ 14 M -1 pl and 19.6 M -1 pl for the chaotic and monodromy potentials, respectively, which CMB-S4 will improve to α=f ≲ 9 M -1 pl and 13 M -1 pl . 7 These results improve upon bounds from larger scales: non-Gaussianity limits α=f ≲ 32.3 M -1 pl and α=f ≲ 46.5 M -1 pl , respectively [51,53], while for primordial black hole (over) production these constraints are α=f ≲ 21.9 M -1 pl -24.9 M -1 pl and α=f ≲ 35.9 M -1 pl [55,60]. Our results place similarly tight constraints for Starobinsky inflation (evident in Fig. 6) and natural inflation (similar to those for chaotic inflation).
The couplings we constrain here do not exclude regimes that could be relevant for magnetogenesis [41] or gravitational leptogenesis [47]. Further, non-Gaussianities can be generated during preheating [129][130][131]; however, there is no indication the scenario studied here would have a significant effect on observable scales.
The significant gravitational response (of the tensor part of the metric) to the preheating dynamics studied here motivates a similar investigation into the metric's scalar degrees of freedom (i.e., curvature perturbations) and the associated potential for primordial black hole formation. Nonlinear gravitational effects could potentially alter our findings here, which may be studied in a framework similar to that employed by Ref. [132]. Additionally, we could investigate the extension of our results to more complicated gauge theories [45]. Finally, we have ignored the effects of the backreaction of any charged matter on the gauge preheating process. Charged particles are produced via the Schwinger effect, leading to a nonzero conductivity in the resulting plasma which could possibly damp gauge-field production [73,133]. We leave a detailed investigation of these questions for future work.
The evolution of the background FLRW spacetime [Eq. ( 3)] is governed by the Einstein equations
where τ is the conformal time and H ¼ a 0 =a is the conformal Hubble parameter related to the usual Hubble parameter via H ¼ aH. For the action in Eq. ( 1), the spatially averaged energy density and pressure are
7 These constraints are often phrased in terms of the instability parameter ξ, defined in Eq. (A30), evaluated when CMB modes left the horizon. The Planck limits using our results are ξ CMB ≲ 0.9 for both the chaotic and monodromy potentials; from CMB-S4 projections these are both ∼0.6.
where brackets hÁ Á Ái denote a spatial average and we have defined the electric and magnetic fields
Finally, our Fourier convention is set by
To study the production of gravitational waves, we consider an FLRW metric including tensor (but not scalar or vector) perturbations,
where ∂ i h ij ¼ h ii ¼ 0 is the transverse-traceless perturbation of the spatial metric. In this work we compute the dynamics of h ij passively-that is, we determine the amount of gravitational waves sourced by the axion and the gauge fields but we neglect their backreaction onto the system. From the linearized Einstein equations we obtain the equation of motion for h ij ,
In this expression, T TT ij is the transverse-traceless component of the stress-energy tensor,
where the transverse-traceless projector is
The background Einstein equations show that 2ð2H 0 þ H 2 Þ is proportional to the background pressure pðτÞ. Since usually all the modes in our simulations are subhorizon, this term in Eq. (A9) induces a negligible amount of dispersion, and so we neglect it. Under the approximations discussed after Eq. (A8), the stress-energy tensor is expressed as
where ḡμν denotes the unperturbed FLRW metric. Moreover, since the terms proportional to ḡij are pure trace, they do not contribute to Eq. (A10) (they are projected out by P ij ). As a consequence, only the first two terms in Eq. (A12) are responsible for gravitational wave generation. The stress tensor of gravitational waves is [135] T gw μν ¼
with the sum over i and j implied. The corresponding fractional energy density in gravitational waves as
The spectrum redshifted to today is related to the spectrum at emission by the transfer function [14] Ω gw;0 ðfÞh
where k phys is the physical wave number and H is the Hubble parameter evaluated at the time when the spectrum is being computed. Above, g 0 =g à is the ratio of thermal degrees of freedom today to matter-radiation equality. Here we assume that the Universe was radiation dominated from the moment of emission until matter-radiation equality.
The Euler-Lagrange equations for the action Eq. ( 1) yield the dynamics of the gauge fields,
In terms of the gauge potentials, the ν ¼ i equations are the dynamical equations of motion
while the ν ¼ 0 equation is the Gauss constraint
The four components of the gauge field are not physical: the theory is invariant under gauge transformations where
where β is an arbitrary function. This freedom allows us to eliminate one degree of freedom by fixing a gauge.
Together with Gauss's law, the gauge sector represents only two physical, dynamical degrees of freedom.
Likewise, the Euler-Lagrange equation for the inflaton provides its equation of motion,
The system described by Eqs. (A18) and (A22) is difficult to solve analytically due to the nonlinear interactions between the axion and gauge field. However, provided the coupling α=f is not too large, a linear treatment provides some insight into the dynamics of this system during inflation and the early phase of preheating.
We begin by linearizing the system of equations, expanding the axion about its homogeneous background ϕ ¼ φðτÞþδϕ and treating the gauge field as a first-order perturbation. We choose the temporal gauge, A 0 ¼ 0. Note that, at linear order this gauge is equivalent to Coulomb gauge ∂ i A i ¼ 0 via the gauge constraint; the gauge field may be taken to be purely transverse. Expanding into Fourier modes
where the polarization vectors ε AE ðkÞ form an orthogonal basis of polarizations transverse to the momentum k,
and satisfy the relations
At linear order, the equation of motion, Eq. (A19), reduces to
where we have defined the instability parameter
Equation (A29) shows that the interaction between the inflaton and the gauge fields induces a tachyonic instability (for modes with k=H < 2jξj) for one of the two helicity states (A þ if ϕ 0 > 0, or conversely A -if ϕ 0 < 0). The modes within this band experience (for values of ξ ≳ 2) exponential enhancement slightly before horizon crossing and then approach a constant value on superhorizon scales.
In the limit of a nearly constant ξ, the solutions to this equation are Whittaker functions (see, e.g., [40]). Assuming ϕ 0 > 0, near horizon crossing (k=aH ∼ 1) these solutions are well approximated by
These results evince the importance of the instability parameter ξ which controls the exponential enhancement of the gauge field.
Computing the backreaction of the gauge field onto the background dynamics requires the integrals
where A λ ðkÞ denotes the helicity modes of the gauge potentials. [Note that Eqs. (A32) to (A34) assume the temporal gauge.] For ξ ≳ 4 [59], fairly accurate approximations of these quantities are given by [124]
Substituting Eq. (A35) into Eq. (A22) demonstrates that gauge fields exert an additional friction term for the inflaton [see Eq. ( 22)]. Since ξ is proportional to the inflaton background's velocity, it is expected to grow during slow-roll inflation. As such, backreaction may significantly alter the evolution of the inflaton toward the end of inflation, even before preheating. In Appendix B we present our procedure to accurately account for these effects during inflation (as used to set initial conditions for the lattice simulations).
As implemented in Refs. [40,41,47], in order to capture the tachyonic enhancement of (one polarization of) the gauge fields during inflation, we numerically integrate the linearized equations of motion for the fluctuations of the gauge-field helicity modes, Eq. (A29), during inflation. These fluctuations are integrated alongside the background dynamics, i.e., Eq. (A2), and the homogeneous part of the inflaton's equation of motion, Eq. (A22). To set initial conditions for the subsequent lattice simulation, we evaluate the background quantities and obtain power spectra from the integrated gauge-field modes between one and two e-folds before ϵ H ¼ 1, marking the end of inflation, as described in Sec. III. By beginning the lattice simulation sufficiently early before the end of inflation we capture nonlinear effects that become important, and thus only use the solutions to the linearized equations while they remain valid.
To account for the gauge fields' backreaction onto the background dynamics during inflation (as described in Appendix A3), Refs. [40,41,47] use the approximations Eqs. (A35) and (A36). The accuracy of these expressions is sensitive to two main assumptions. First, the solution of Eq. (A29) is well approximated by the Whittaker function only in the limit that ξ is constant, and H ¼ -1=τ (in nearde Sitter space). Second, Eqs. (A35) and (A36) themselves are approximations to the integrals Eqs. (A32) to (A34) over Eq. (A31), accurate for ξ ≳ 4 [59]. When ξ is nearly constant but small, the approximate integrals Eqs. (A35) and (A36) are inaccurate, but the Whittaker solution is still valid. However, if ξ is not approximately constant, then the Whittaker solution itself is inaccurate, meaning the integrands of Eqs. (A32) to (A34) themselves are inaccurate, regardless of the size of ξ.
At low couplings (which roughly correspond to those which do not achieve complete preheating), the gauge fields, while amplified, have a negligible effect on the background evolution during inflation. Thus, while ξ is not large enough for Eqs. (A35) and (A36) to be valid, their effect is small enough that this error is unimportant. At larger couplings for which preheating is complete, the gauge fields have a more substantial impact on the background evolution during the final e-folds of inflation, but these effects can in principle be captured by linearized calculations (without further approximations). However, in this regime analytic approximations are in fact never simultaneously valid: during slow roll when the Whittaker solution is approximately valid, jξj < 4, meaning the approximation to the analytic integral is not trustworthy [59]. In fact, by the time ξ is large enough to validate the approximation to the integral, the inflaton is no longer slowly rolling, so the Whittaker solution is no longer valid. Further, the backreaction of the gauge fields onto the background dynamics is no longer negligible (during the final e-folds of inflation), meaning the background evolution using these approximations may be inaccurate. At even stronger couplings, nonlinear dynamics have an important impact on the background evolution during the final e-folds of inflation. The friction exerted by the gauge fields on the inflaton background (beyond what is captured by a linear treatment) can even postpone the end of inflation. In this regime we might not trust any calculation using only the linear equations of motion toward the end of inflation, so we leave a detailed study of such strong coupling to future work.
Here we improve upon this procedure by computing the gauge-field integrals Eqs. (A32) to (A34) via numerical quadrature of the numerically evolved gauge fields. As chosen for the analytic results Eqs. (A35) and (A36),w e take as a UV cutoff for these integrals k ¼ 2aHξ, the upper end of the tachyonic instability band; we find that, in the regime where the gauge fields' contribution to the background evolution is non-negligible, the integration is insensitive to slight variation about this cutoff. With this method, the evolution of the background in tandem with the linearized fluctuations agrees very well with the background dynamics of the subsequent lattice calculation (for the overlapping final one to two e-folds of inflation).
In practice, it is prohibitively expensive to include the full dynamic range of fluctuations which are amplifiedthe band of tachyonic instability k ≲ 2aHξ is proportional to the comoving horizon H ¼ aH, which increases exponentially during inflation. To address this, we only include modes which experience amplification during the last ∼20 e-folds of inflation, well before backreaction has any impact on the background dynamics of inflation (for the values of α=f we consider). We sample this range of modes (spanning 8-10 orders of magnitude) logarithmically, yielding good results from quadrature of Eqs. (A32) to (A34) over ln k. A final complication is that fluctuations deep inside the horizon oscillate extremely rapidly compared to the Hubble rate. However, this regime (k ≫ 2aHξ) is precisely the regime where the Wentzel-Kramers-Brillouin (WKB) approximation
with ω AE ðk; τÞ 2 ¼ k 2 ∓ 2aHξk is extremely accurate. Thus, we only numerically integrate a given mode k when k ≤ 50 × 2aHξ, using the WKB solution until this point.
Using an adaptive ordinary differential equation integrator (in practice, SciPy's DOP853 routine [136,137]) ensures that we take as large of time steps as possible while maintaining a prescribed relative accuracy (typically to one part in 10 11 ). A similar strategy was employed in Ref. [138] to study this model's dynamics during both inflation and preheating; however, the validity of such a method applied to study preheating is unclear. In Fig. 9 we investigate the difference between the analytic approximations [Eqs. (A35) and (A36)] and the numerical computations described above, fixing a chaotic inflaton potential. We consider couplings where the background evolutions agree well (α=f ¼ 8 M -1 pl ) and disagree (α=f ¼ 14 M -1 pl ). We additionally compare to the numerical quadrature of Eqs. (A32) to (A34) using the actual Whittaker-function solutions, which demonstrates agreement with the approximations Eqs. (A35) and (A36) once jξj ≳ 3-4. Earlier in inflation, when exact slow roll (and so the Whittaker-function solution) is a good approximation, we observe agreement between the exact integrals of the Whittaker function solutions and the integrals over numerically evolved gauge-field perturbations. However, in neither case do the analytic approximations agree with the results of numerical evolution and quadrature at any time.
While these errors have no significant effect on the inflationary background deep in inflation, Fig. 9 indicates that the analytic approximations greatly overestimate hE 2 þ B 2 i and hE • Bi during the final two e-folds of inflation, regardless of the coupling to the axion. In the strong-coupling case (where gauge-field backreaction is significant), the analytic result's overestimation of the friction exerted by the gauge fields onto the axion background leads to inflation ending artificially late, as indicated by the vertical lines (marking when ϵ H ¼ 1)inFig.9. In this case, with the analytic-based evolution we would end up initializing the lattice simulation with background values roughly half an e-fold later than intended (depending on the coupling). Note that α=f ¼ 8 M -1 pl is roughly the largest coupling for which the background evolutions visibly agree.
In this Appendix we study the dynamics after preheating to evaluate whether the Universe remains radiation FIG. 9. Comparison of the analytic approximations (solid blue lines) to hE 2 þ B 2 i (top panels) and hE • Bi (bottom panels) (i.e., Eqs. (A35) and (A36) to those obtained via direct numerical quadrature of gauge-field fluctuations evolved with their linearized equations of motion (dashed orange lines). Both background evolutions are initialized identically, with N ¼ 0 corresponding to the time inflation ends (i.e., ϵ H ¼ 1) in the evolution using the analytic approximations. The vertical lines indicate the point at which inflation ends according to both evolutions. We additionally plot the result of direct numerical quadrature of Eqs. (A32) to (A34) using the Whittaker-function solutions to evaluate the accuracy of Eqs. (A35) and (A36) as approximations to this exact integral. dominated. In particular, we seek to quantify the amount by which the fractional energy density in gravitational waves is suppressed due to any deviation from radiation domination between preheating and the end of reheating.
As discussed in Sec. VA, at the end of all of our simulations a small amount of the total energy remains in axion fluctuations [Oð1%Þ for the cases where preheating is most efficient]. While these are typically relativistic at this point, because the axion is rather massive these fluctuations become nonrelativistic within just a few e-folds and gravitate as pressureless matter. Depending on their lifetime, these massive axions can come to dominate the energy density of the Universe, leading to a period of matter domination. Even if the Universe does not become matter dominated, the Universe could still depart significantly from radiation domination. The resulting effect on the gravitational wave transfer function depends on (i) how much energy remains in the axion, Ω ϕ , and (ii) its decay rate into radiation, Γ.
Noting that ρ GW ∝ a -4 for the modes generated inside the horizon during preheating, the fractional density in gravitational waves scales as
where a subscript 0 denotes some reference time (which for our purposes is the end of preheating) and w ≡ p=ρ is the time-dependent equation of state. Thus, if the Universe remains radiation dominated (w ¼ 1=3), then Ω gw remains constant (as it redshifts at the same rate as the rest of the Universe). On the other hand, if the Universe is dominated by matter (w ¼ 0) or by a mixture of matter and radiation (0 <w<1=3), then Ω gw decreases with the scale factor. Thus, any suppression of gravitational waves after preheating depends on both the (evolution of the) equation of state w and the duration for which w<1=3.
In the standard reheating scenario (without a preheating phase), the equation of state of the Universe is that of matter due to the inflaton condensate's oscillation about the minimum of its potential. 8 Matter domination persists until the decay of the inflaton into relativistic species (in our case, bosons) becomes efficient, which occurs when the inflaton's decay rate Γ becomes comparable to present the Hubble scale H. Accounting for preheating affects this description in two significant ways: the equation of state is not (initially) w ¼ 0 and the coupled sector is already (highly) occupied. The former would delay the possible onset of matter domination before reheating completes, while the latter significantly alters the inflaton decay rate due to Bose enhancement. Accounting for the nontrivial phase space distribution of the gauge fields f A ðpÞ in our models enhances the decay rate by a factor ∼1 þ 2f A ðm ϕ =2Þ [139].
For simplicity, we study the dynamics after preheating with the standard Boltzmann equations describing threebody decay of inflaton particles at rest into two relativistic daughter particles [140],
Rather than fixing the decay rate Γ and the initial (i.e., postpreheating) fractional energy in the axion Ω ϕ according to our model (and the final state of our simulations), we instead study a wide range of parameter space and use order-of-magnitude estimates to determine whether this period significantly affects our results.
In Fig. 10 we numerically integrate Eqs. (C2) and (C3) over a number of decades in Γ and Ω ϕ , plotting the amount by which Ω gw would reduce due to w being less than 1=3 during this epoch [using Eq. (C1)]. In this figure, the vertical slice at Ω ϕ ¼ 1 corresponds to reheating without any preheating phase. In our model, the zero-temperature decay rate of the axion into gauge fields is [141]
To estimate f A ðm ϕ =2Þ we extract the occupation number of the gauge fields from the simulations, which varies between ∼10 8 and 10 12 , in line with results from tachyonic resonance in scalar-field preheating [142]. Thus, the relevant portion of gravitational waves would drop by no more than ∼10%, and likely by less than 1%. Naturally, the regime for which Ω gw;0 h 2 exceeds the projected CMB-S4 bounds on ΔN eff corresponds to both the smallest Ω ϕ (approaching 1%) and the largest occupation numbers. Therefore, our results are most robust for the coupling regimes which would be probed or excluded by CMB-S4.
In principle, the true dynamics depend on the time dependence of the entire phase-space distribution fðpÞ, as well as other possible thermal effects due to sectors coupled to the daughter particles. We have further approximated that the axion particles are at rest (while at the end of the simulations all would have nonzero momentum). However, Fig. 10 makes clear that we are far from the regime where Ω gw is suppressed by any appreciable amount: the effective decay rate could drop by ∼2-8 orders of magnitude before Ω gw would redshift to half of its original value. As such, we conclude that the Universe would remain radiation dominated after preheating, and so our reported Ω gw;0 h 2 are robust.
Note that, in an abuse of notation, we use f to denote both the axion decay constant as well as the frequency of the gravitational waves.
https://github.com/zachjweiner/pystella.
Note that the ultraviolet parts of the spectra in Fig.5, which grow as k 4 , are not physical signals but the product of vacuum modes in the simulation.
More precisely, the equation of state of a coherently oscillating scalar field time averages to zero.