Constraining axion inflation with gravitational waves from preheating
Peter Adshead John T. Giblin Jr Mauro Pieroni Zachary J. Weiner





adshead@illinois.edu
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA


giblinj@kenyon.edu
Department of Physics, Kenyon College, Gambier, Ohio 43022, USA
CERCA/ISO, Department of Physics, Case Western Reserve University, Cleveland, Ohio 44106, USA


mauro.pieroni@uam.es
Instituto de Física Teórica UAM/CSIC, Calle Nicolás Cabrera 13-15, Cantoblanco E-28049 Madrid, Spain
Departamento de Física Teórica, Universidad Autónoma de Madrid (UAM) Campus de Cantoblanco, 28049 Madrid, Spain
Theoretical Physics, Blackett Laboratory, Imperial College, London, SW7 2AZ, United Kingdom


zweiner2@illinois.edu
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA


    We study gravitational wave production from gauge preheating in a variety of inflationary models, detailing its dependence on both the energy scale and the shape of the potential.
    We show that preheating into Abelian gauge fields generically leads to a large gravitational wave background that contributes significantly to the effective number of relativistic degrees of freedom in the early universe, N_eff.
    We demonstrate that the efficiency of gravitational wave production is correlated with the tensor-to-scalar ratio, r.
    In particular, we show that efficient gauge preheating in models whose tensor-to-scalar ratio would be detected by next-generation cosmic microwave background experiments (r ≳ 10^-3) will be either detected through its contribution to N_eff or ruled out.
    Furthermore, we show that bounds on N_eff provide the most sensitive probe of the possible axial coupling of the inflaton to gauge fields regardless of the potential.


Constraining axion inflation with gravitational waves from preheating
    Zachary J. Weiner
    December 30, 2023
=====================================================================







§ INTRODUCTION


While the qualitative predictions of inflation are well motivated and well understood <cit.>, 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 <cit.>, 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. <cit.> for reviews).
The rapid production of inhomogeneities during this phase generically sources a significant gravitational wave background <cit.>.
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 <cit.>, with recent studies exploring nonminimally coupled scalar fields <cit.>.
Alternatively, the inflaton could couple directly to gauge fields <cit.>.
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 <cit.>.
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 <cit.>, observable gravitational waves <cit.>, primordial black holes <cit.>, μ-distortions <cit.>, primordial magnetic fields <cit.>, and the generation of the baryon asymmetry <cit.>.

Preheating into gauge fields via a Chern-Simons coupling  was first studied within the context of chaotic inflation in Ref. <cit.>, and the model's viability for magnetogenesis was explored in Ref. <cit.>.
Subsequently, we extended this work in Ref. <cit.> 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. <cit.> (and <ref> 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. <cit.> 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. <cit.> 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 <cit.>.
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 <ref> 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 <ref>, we outline some analytic estimates in <ref>.
We investigate the effect of the energy scale of inflation and the shape of the inflationary potential on preheating and gravitational wave production in <ref>, respectively.
In <ref> we draw our conclusions.
In <ref> we detail the equations of motion of the system and summarize a linear analysis thereof.
<Ref> details our procedure for setting accurate initial conditions.
Finally, in <ref> 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).






§ GAUGE FIELDS DURING AND AFTER AXION INFLATION


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 <ref>.






 §.§ Models of inflation


We consider the action for a pseudoscalar inflaton ϕ minimally coupled to gravity and to an Abelian gauge field A_μ,

    S
               = ∫^4 x √(-g)[
                ^2/2 R - 1/2∂_μϕ∂^μϕ - V(ϕ) 
       
                - 1/4 F_μν F^μν - α/4 fϕ F_μνF̃^μν],

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.[The extension to different Abelian gauge groups [such as 𝒩 copies of a U(1) gauge group] is straightforward.
Conversely, the extension to non-Abelian gauge groups is nontrivial.
See, for example, Refs. <cit.>.]
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

    F̃^μν   = 1/2ϵ^μναβ F_αβ.

where ϵ^μνρσ is the Levi-Civita symbol with convention ϵ^0123 = 1 / √(-g).
We set c = ħ = k_B = 1 and denote by = 1 / √(8 π G_N) the reduced Planck mass.
The background spacetime is the mostly plus, conformal Friedmann-Lemaître-Robertson-Walker (FLRW) metric,

    s^2 = a(τ)^2 ( - τ^2 + 𝐱^2 ),

where the conformal time coordinate τ is related to the cosmic time coordinate via τ =  t / a.
Throughout, primes denote derivatives with respect to conformal time, e.g., a' ≡∂ a / ∂τ, and dots denote derivatives with respect to cosmic time, e.g., ȧ≡∂ 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 <cit.>), 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 <ref> the inflaton is pseudoscalar.
[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 <cit.>.
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 <cit.>.]

Specifically, we explore the following classes of models:


  * Chaoticlike monomial models, which have the potential
	
    V(ϕ) = m^4-n|ϕ|^n.

    Notice that for n = 2 this simply reduces to chaotic inflation <cit.>, the case we consider here.
    For n = 2 the Planck normalization of the scalar power spectrum <cit.> sets m^2 = 1.9 × 10^-11 ^2, 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 ) ) <cit.>.
    In this case, with v = √(8 π) the Planck normalization sets V_0 = 5.9 × 10^-10 ^4, 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.


  * Starobinsky-like models, which have potentials given by <cit.>[For other models of this class see also Refs. <cit.>.]
	
    V(ϕ) = V_0 ( 1 - exp( |ϕ|/v) )^2.

    For v = 10  / 3 the Planck normalization sets V_0 = 6.2 × 10^-10 ^4, in which case n_s = 0.969 and r = 0.016.
    The class of α-attractor models, motivated by supergravity <cit.> (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.


  * Monodromy inflation, corresponding to the potential
    
    V(ϕ) = μ^3 ( √(ϕ^2 + ϕ_c^2) - ϕ_c ).

    This is part of a broad class of string-theory–motivated models in which large (super-Planckian) field displacement is obtained by wrapping the inflaton trajectory into a series of (sub-Planckian) fundamental circuits <cit.>.
    Note that deep in the inflationary phase (ϕ≫ϕ_c) the potential is well approximated by a linear potential (n=1) of the class in <ref>.
    For ϕ_c =  / 10 the Planck normalization sets μ = 6.0 × 10^-4, for which n_s = 0.975 and r = 0.067.


  * Hilltoplike models, with potentials given by <cit.>
    
    V(ϕ) = V_0 ( 1 - ( |ϕ|/v)^p )^2,

	where p ≥ 2.
    We consider a variety of p and v and list the corresponding model fits and predictions in <ref>.

    


  * D-brane models, with potentials of the form <cit.>
	
    V(ϕ) = V_0 ( 1 - ( v/|ϕ|)^p )^2.

    Fixing p = 2 and v = /2, the Planck normalization sets V_0 = 7.5 × 10^-11 ^4, in which case n_s = 0.975 and r = 2.2 × 10^-3.
 





 §.§ Gauge fields on the axion background


The equations of motion that result from the variation of <ref> 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 <ref>.

In a homogeneous axion background, the helical polarizations of the gauge field obey <cit.>

    A_±”(𝐤) + k ( k ∓ 2 ℋξ) A_±(𝐤) = 0,

where we define the instability parameter

    ξ≡α/2 fϕ'/ℋ

and the conformal Hubble parameter, ℋ = a' / 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 ω∼√(k(k - 2 ℋξ)) (assuming ξ > 0) for modes k < 2ℋξ.
The polarizations are amplified (relative to the conformally invariant radiation solution A_rad,±) by a factor of

    |A_±/A_rad, ±|∼ e^π/2ξ±π/2ξ.

The amplification is controlled by the parameter ξ, which depends on the inflaton velocity, ϕ' / ℋ.
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/ℋ≲ 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 <ref>).
We use numerical simulations to capture these nonlinear dynamics in later sections.





 §.§ Gravitational waves and delta N effective


The tachyonic production of gauge-field modes leads to large anisotropic stresses which results in the copious production of gravitational waves <cit.>.
These gravitational waves form a stochastic background with frequencies near the Hubble scale during reheating, corresponding to f ∼ 10^9 Hz today[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.] 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. <cit.>, may be constrained by CMB and Big Bang Nucleosynthesis 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

    Ω_gw,0 h^2
    	= ∫ln k  1/ρ_0ρ_gw,0/ln k.

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 <cit.>

    Ω_gw,0 h^2/Ω_γ,0 h^2   = 7/8( 4/11)^4/3Δ N_eff,

where the present energy density in photons is Ω_γ,0h^2 = 2.47 × 10^-5.
From the Planck limit, |Δ N_eff|≲ 0.33 <cit.>, 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σ <cit.>, improving the upper limit by an order of magnitude to

    Ω_gw,0 h^2 ≲ 1.68 - 3.36 × 10^-7.

A separate analysis <cit.> 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 <cit.> and Euclid <cit.> 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 <cit.> prompts consideration of direct gravitational wave detectors as alternative probes of stochastic gravitational wave backgrounds.
In particular, the presently operating LIGO <cit.> and VIRGO <cit.> detectors (with KAGRA <cit.> expected to join this network of detectors in the near future), as well as future missions like LISA <cit.> and the Einstein Telescope <cit.>, 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 signal-to-noise ratio (SNR), given by <cit.>

    SNR^2
           = T ∫_f_min^f_max f  (Ω_gw/Ω_s)^2
       ≡ T ∫_f_min^f_max f  [ Ω_gw/4 π^2 S_n f^3 / (3 H_0^2)]^2,

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 <cit.>.
At the same time, current methods (see, for example, <cit.>) 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.






§ NUMERICAL METHODS


Our numerical approach is very similar to that employed in Ref. <cit.>.
We numerically integrate the classical equations of motion of the inflaton, <ref>, and the gauge fields, <ref>, in an FLRW background governed by <ref>.
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 fourth-order Runge-Kutta method.
Spacetime expansion is implemented self-consistently by computing the current energy density and pressure, <ref>, 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, <ref>, 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 <cit.>, 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 <cit.>.
Doing so requires evolving additional, redundant degrees of freedom, which in our case means evolving all four components of the gauge potential,

    A_ν”   = ∇^2 A_ν
        + η_βνα/f∂_αϕ(1/2ε^αβρσ F_ρσ),

obtained by applying the Lorenz gauge choice, ∂_μ A^μ = 0, to <ref>.
The results of <cit.> 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 <cit.>, the results of which reproduce those of Ref. <cit.> (and those presented in <ref> below).
In addition, in Ref. <cit.> 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. <cit.>.
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. <cit.>, 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 <ref> (that |ξ|≳ 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 <ref> were developed <cit.>), 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  ^-1 for chaotic inflation), these errors artificially offset the end of inflation by ∼ 0.5 e-folds.
As a result, the use of <ref> 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 <ref> via quadrature of the numerically integrated gauge-field fluctuations.
See <ref> for full details of this procedure.

Other than improving the gauge-field power spectra, our procedure for generating initial conditions is unchanged from Ref. <cit.>.
We initialize the lattice simulations one to two e-folds (depending on the model; see <ref>) 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, Gaussian-random 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 and projecting the polarization fields, A^±(𝐤), onto their vector components A_i(𝐤), which automatically yields ∂_i A_i = 0.

The second difference in our numerical implementation is that, in contrast to the pseudospectral solver used in Ref. <cit.>, we represent spatial derivatives with finite differencing, using fourth-order–accurate centered difference stencils.
In particular, whereas in Ref. <cit.> we evolved the equations of motion for the tensor metric fluctuations h_ij in Fourier space [computing and Fourier transforming the source term <ref> and projecting with <ref> at each time step], we here evolve the equation of motion in position space without applying the transverse-traceless projection <cit.>,

    u_ij” + 2 ℋ u_ij' - ∂_k ∂_k u_ij   = 2/^2 T_ij.

Instead, we apply the transverse-traceless projection to u_ij to obtain

    h_i j   = ( P_i l P_j m - 1/2 P_i j P_l m) u_l m,

only when computing the gravitational wave spectrum Ω_gw(k) via <ref>.
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 ,[https://github.com/zachjweiner/pystellagithub.com/zachjweiner/pystella] a Python-based, MPI-parallel and GPU-accelerated code making use of  <cit.> and  <cit.> for the generation of OpenCL code to run on GPUs.
Because  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. <cit.>; we also checked a representative sample of our results against simulations with 512^3 points to show convergence.
The Friedmann constraint, <ref>, is satisfied to a precision of 𝒪(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 <ref>.
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^2 = 1.






§ ANALYTIC ESTIMATES


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 <ref> 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 <cit.> to make contact between the characteristics of the inflationary potential and the efficiency of gravitational wave production from preheating.





 §.§ Efficiency of preheating


From the Friedmann equation, <ref>, by making use of the approximation for the energy density in the gauge fields, <ref>, we can derive the approximation

    H^2/^2    = 1/( 3 - ϵ_ϕ)V/^4 [ 1 + V/^4 1.4 × 10^-4  e^2 πξ/ξ^3 ( 3 - ϵ_ϕ)^2 ] 
       
            + 𝒪(V/^4)^3,

where we introduce the slow-roll parameters

    ϵ_ϕ≡ϕ̇^2/2 ^2 H^2 ,   ϵ_H≡ -Ḣ/H^2.

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.
<Ref> 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 ϕ'/ℋ, 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 <ref>, the gauge fields induce a new friction term in the equation of motion of the inflaton.
Substituting the spatially averaged interaction, <ref>, into the equation of motion for the inflaton, <ref>, while making use of the lowest order approximation of <ref> and assuming ϕ”/ℋ^2 ≪, we can derive the usual slow-roll condition, now including the effect of the gauge-field backreaction,

    ^2 V_, ϕ/V   ≃ -ϕ̇/H[ 1 + V/^4 1.4 × 10^-4  e^2 πξ/ξ^3 ( 3 - ϵ_ϕ)^2 ] 
       
        + 2.4 × 10^-4α/fV/^2 e^2 πξ/ξ^4 ( 3 - ϵ_ϕ)^2.

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 <ref>.
In this scenario, gauge-field–induced friction provides an additional mechanism to enforce the flatness of the potential required for slow-roll inflation <cit.>.

There are two effects at play in <ref>.
The first is that the backreaction changes the Hubble rate, increasing Hubble friction (on the inflaton), while the second is the backreaction of the 𝐄·𝐁 term.
The impact of these two effects can be understood by comparing <ref> with <ref>.
As the same (V-dependent) term is compared with unity (for the backreaction on H) and with √(2 ϵ_ϕ) (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 / ^4) 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 <ref> 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 ξ.
However, if V 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 <ref> (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,

    ρ_gauge/ρ   ≈1.4 × 10^-4/3(H / )^2/ξ^3 e^2 πξ.

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 / 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

    ρ_gauge, 1/ρ_gauge, 2
        = ( H_1/H_2)^2 ( ξ_2/ξ_1)^3 e^2 π (ξ_1 - ξ_2) = 1

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

    ( H_1/H_2)^2 ( ξ_2/ξ_1)^3
        ≈ 1,

yielding ξ_2 = ξ_1 (H_2 / H_1)^2/3.
Alternatively, if ξ_1 and ξ_2 are both large,

    ( H_1/H_2)^2 e^2 π (ξ_1 - ξ_2)≈ 1,

in which case ξ_2 = ξ_1 + ln (H_1 / H_2) / π.
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 <cit.> 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 9  ^-1.
With this as a baseline, we apply <ref> 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,

    V(ϕ) = 1/2 m_ϕ^2 ϕ^2.

In this case, the Hubble parameter scales as H ∼ m_ϕ, since at the end of inflation

    H^2 = ρ/3 ^2≈V(ϕ)/3 ^2∼ m_ϕ^2.

In terms of the coupling, <ref> tells us that, for some other mass m_ϕ, 2, we require

    α_2/f = α_1/f + 1/π∂ϕ / ∂ Nln( m_ϕ, 1/m_ϕ, 2)

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.4, in which case <ref> reduces to

    α_2/f = α_1/f + 1.1/log_10( m_ϕ, 1/m_ϕ, 2).

Thus, for comparable preheating efficiency we expect to need to increase the axion-gauge coupling α / f by ∼ ^-1 for each order of magnitude we reduce m_ϕ.





 §.§ Gravitational wave production and the “rule of thumb”


While regimes of highly efficient preheating exist for all models (as we show in <ref>), 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. <cit.> 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. <cit.> estimates the peak amplitude as

    Ω_gw, 0≈ 2.3 × 10^-4α^2 β w^2 k_∗/σ(H_∗/k_∗)^2,

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 <ref>: 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, k_∗ / H_∗∼ m_ϕ / H_∗∼√(V” / V)—i.e., 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 <ref>.






§ NUMERICAL RESULTS







 §.§ Dependence on inflationary scale


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 <ref>.
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 <ref>.
In <ref> 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 we use a box length L = 7.5   m_ϕ^-1 in order to ensure sufficient short-wavelength resolution in all cases.
We immediately observe that <ref> 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 <ref> 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, as predicted by the analytic estimate <ref>.

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, corresponds to 9  ^-1≲α / f ≲ 10.4  ^-1), ρ_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 <ref> 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,

    ρ_δϕ   ≡ρ_ϕ - ρ_⟨ϕ⟩
       = ⟨1/2 a^2ϕ'^2 + 1/2 a^2∂_i ϕ∂_i ϕ + V(ϕ) ⟩
       
            - ( 1/2 a^2⟨ϕ⟩'^2 + V(⟨ϕ⟩) ),

to the energy in the gauge fields, ρ_gauge.
In <ref> 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  ^-1), 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  ^-1≲α / f ≲ 10.2  ^-1, the inflaton becomes totally fragmented and the energy in 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  ^-1 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 <ref>.

Increasing α / f past 11.6  ^-1, 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 to 100%, as shown in <ref>, 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. <cit.>) are very similar to those presented more recently in Ref. <cit.> 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  ^-1).
As observed in Ref. <cit.>, the inflaton is momentarily “trapped,” as depicted in <ref>.

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  ^-1), 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 <cit.>.

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 <ref> 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 <ref>.
For the strongest couplings we simulate, the total integrated Ω_gw, 0 h^2 exceeds 10^-6, which is already ruled out by Planck data <cit.>.
To quantify this claim, in <ref> 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 <ref> we also observe that gravitational wave production is moderately less effective at lower inflationary scales.
Returning to the arguments of <ref>, 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 / a H < ξ∼α / f, from the rule of thumb, <ref>, we expect Ω_gw to decrease with (k_∗ / H)^2 ∼ (α / f)^2.
To verify this scaling, in the lower panel of <ref> we scale Ω_gw, 0 h^2 by the squared ratio of each simulation's coupling α / f to its “adjusted” value from <ref>.
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 <ref> demonstrates this observation, and also depicts the √(m_ϕ) scaling of the characteristic present-day frequencies at which these signals would be observed.[Note that the ultraviolet parts of the spectra in <ref>, which grow as k^4, are not physical signals but the product of vacuum modes in the simulation.]


In 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 <ref>) is 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  ^-1 for m_ϕ = 6.16 × 10^-6).
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 <ref>.
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 <ref> accounts for (which assumes the Universe was radiation dominated from the time of emission until matter-radiation equality).
In <ref> 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.





 §.§ Dependence on the shape of the potential


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 <ref> and discuss the extent to which the results of <ref> are modified.

In <ref>, 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 <ref>.
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 <ref>, 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.

In <ref> 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 <ref>.
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 <ref>, low-scale models require larger coupling for efficient preheating and present correspondingly weaker detection prospects.
This observation, together with <ref>, leads us to make the broad claim that for models of inflation with tensor-to-scalar 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 <ref> 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 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 <ref> 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 <ref>, increasing v decreases k_∗ / H_∗, which (referring to <ref>) correlates to larger tensor-to-scalar ratios, r.
From <ref> 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 <cit.>) still provide significantly weaker constraints on α / f for small-r (and so small-ϵ_H) models, as ∂ϕ / ∂ N ∼√(ϵ_H).
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.






§ CONCLUSIONS


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. <cit.>, 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 gauge-field 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  ^-1 and 19.6  ^-1 for the chaotic and monodromy potentials, respectively, which CMB-S4 will improve to α / f ≲ 9  ^-1 and 13  ^-1.
[These constraints are often phrased in terms of the instability parameter ξ, defined in <ref>, 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.
]
These results improve upon bounds from larger scales: non-Gaussianity limits α / f ≲ 32.3 ^-1 and α / f ≲ 46.5  ^-1, respectively <cit.>, while for primordial black hole (over) production these constraints are α / f ≲ 21.9 ^-1 - 24.9  ^-1 and α / f ≲ 35.9  ^-1 <cit.>.
Our results place similarly tight constraints for Starobinsky inflation (evident in <ref>) and natural inflation (similar to those for chaotic inflation).

The couplings we constrain here do not exclude regimes that could be relevant for magnetogenesis <cit.> or gravitational leptogenesis <cit.>.
Further, non-Gaussianities can be generated during preheating <cit.>; 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. <cit.>.
Additionally, we could investigate the extension of our results to more complicated gauge theories <cit.>.
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 <cit.>.
We leave a detailed investigation of these questions for future work.



We thank Mustafa Amin and Jessie Shelton for useful discussions and Valerie Domcke for comments on our draft.
Z.J.W. thanks Andreas Kloeckner for generous support and advice on the development of .
The work of P.A. was supported in part by NASA Astrophysics Theory Grant No. NNX17AG48G.
J.T.G. is supported by the National Science Foundation Grant No. PHY-1719652.
M.P. acknowledges the support of the Spanish MINECOs “Centro de Excelencia Severo Ochoa” Programme under Grant No. SEV-2016-059.
This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 713366.
Z.J.W. is supported in part by the United States Department of Energy Computational Science Graduate Fellowship, provided under Grant No. DE-FG02-97ER25308.
The development of  made use of the Extreme Science and Engineering Discovery Environment (XSEDE) <cit.> through allocation TG-PHY180049, which is supported by National Science Foundation Grant No. ACI-1548562, and also made use of hardware purchased by the National Science Foundation, Kenyon College, and the Kenyon College Department of Physics.
This work made use of the Illinois Campus Cluster, a computing resource that is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with the National Center for Supercomputing Applications (NCSA) and which is supported by funds from the University of Illinois at Urbana-Champaign.
P.A. acknowledges the hospitality of the Yukawa Institute for Theoretical Physics at Kyoto University, where some of this work was completed during the YITP-T-19-02 on “Resonant instabilities in cosmology.”








§ EQUATIONS OF MOTION AND LINEAR ANALYSIS


In this appendix, we write down the equations of motion for the system and collect some well-known results about the linear evolution of the system during inflation and the backreaction of the gauge field on the inflaton.





 §.§ Gravitational sector


The evolution of the background FLRW spacetime [<ref>] is governed by the Einstein equations

    ℋ(τ)^2
           = a(τ)^2/3 ^2ρ(τ) 
    ℋ'(τ)
           = - a(τ)^2 /2 ^2( 3 ρ(τ) + p(τ) ),

where τ is the conformal time and ℋ = a'/a is the conformal Hubble parameter related to the usual Hubble parameter via ℋ = a H.
For the action in <ref>, the spatially averaged energy density and pressure are

    ρ(τ)
           ≡⟨ϕ'^2/2 a^2 + (∂_i ϕ)^2/2 a^2+ V(ϕ) + 1/2( 𝐄^2 + 𝐁^2 ) ⟩, 
    
        p(τ)
           ≡⟨ϕ'^2/2 a^2 - (∂_i ϕ)^2/6 a^2 - V(ϕ) + 1/6( 𝐄^2 + 𝐁^2 ) ⟩,

where brackets ⟨⋯⟩ denote a spatial average and we have defined the electric and magnetic fields

    E_i   = 1/a^2( A_i' - ∂_i A_0 ),   
        B_i
        = 1/a^2ϵ_i j k∂_j A_k.

Finally, our Fourier convention is set by

    f(𝐤)
           = ∫^3 x   f(𝐱) e^- i 𝐤·𝐱, 
    
        f(𝐱)
           = ∫^3 k/(2π)^3  f(𝐤) e^i 𝐤·𝐱.






 §.§ Gravitational wave dynamics


To study the production of gravitational waves, we consider an FLRW metric including tensor (but not scalar nor vector) perturbations,

    s^2
           = a(τ)^2 ( - τ^2 + (δ_ij + h_ij)  x^i  x^j ),

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,

    h”_ij - ∂_k ∂_k h_ij + 2 ℋ h'_ij
        - 2 ( 2 ℋ' + ℋ^2 ) h_ij   = 2/^2 T_ij^TT.

In this expression, T_ij^TT is the transverse-traceless component of the stress-energy tensor,

    T_ij^TT = (P_il P_jm - 1/2 P_ij P_lm) T_lm,

where the transverse-traceless projector is

    P_ij   = δ_ij - k_i k_j/k^2.

The background Einstein equations show that 2 ( 2 ℋ' + ℋ^2 ) is proportional to the background pressure p(τ).
Since usually all the modes in our simulations are subhorizon, this term in <ref> induces a negligible amount of dispersion and so we neglect it.

Under the approximations discussed after <ref>, the stress-energy tensor is expressed as

    T_ij   = ∂_i ϕ∂_j ϕ + F_iα F_jβg̅^αβ
       
            - g̅_ij[ 1/2∂_μϕ∂^μϕ + V(ϕ) + 1/4 F_μν F^μν],

where g̅_μν denotes the unperturbed FLRW metric.
Moreover, since the terms proportional to g̅_ij are pure trace, they do not contribute to <ref> (they are projected out by P_ij).
As a consequence, only the first two terms in <ref> are responsible for gravitational wave generation.

The stress tensor of gravitational waves is <cit.>

    T_μν^gw = ^2/4⟨ h_ij,μ h_ ij,ν⟩,

with the sum over i and j implied.
The corresponding fractional energy density in gravitational waves as

    Ω_gw(k)
    	   ≡1/ρρ_gw/ln k
       = 1/24π^2 L^3k^3/ℋ^2∑_i, j| h_ij^'(k, τ)) |^2.

The spectrum redshifted to today is related to the spectrum at emission by the transfer function <cit.>

    Ω_gw,0(f) h^2
    	   = Ω_gw,e(f) (g_0/g_∗)^1/3Ω_r,0 h^2

at frequencies

    f ≈ 2.7 × 10^10k_phys/√( H) Hz,

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.





 §.§ Field equations


The Euler-Lagrange equations for the action <ref> yield the dynamics of the gauge fields,

    ∂_μ( √(-g) F^μν + √(-g)α/fF̃^μν) = 0.

In terms of the gauge potentials, the ν = i equations are the dynamical equations of motion

    0    =
            A_i” - ∂_i A_0' - ∂_j ∂_j A_i + ∂_i ∂_j A_j 
       
            - α/fϵ_i k lϕ' ∂_k A_l + α/fϵ_i k l∂_k ϕ( A_l' - ∂_l A_0 ),

while the ν = 0 equation is the Gauss constraint

    ∂_i A_i' - ∂_j ∂_j A_0
           = - α/fϵ_i j k∂_kϕ∂_i A_j.

The four components of the gauge field are not physical: the theory is invariant under gauge transformations where

    A_μ→Ã_μ = A_μ - ∂_μβ,

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,

    ϕ” - ∂_i ∂_i ϕ + 2 ℋϕ' + a^2  V/ϕ   = - a^2 α/4 f F_μνF̃^μν.

Note that F_μνF̃^μν / 4 = a^4 𝐄·𝐁.






 §.§ Linear theory


The system described by <ref> 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

    A_i(τ,𝐱)
           = ∑_λ = ±∫^3 k/(2π)^3
            A_k^λ(τ) ε^λ(𝐤) e^- i 𝐤·𝐱

where the polarization vectors ε^±(𝐤) form an orthogonal basis of polarizations transverse to the momentum 𝐤,

    ϵ_ijk k_j ε^±_k(𝐤) = ∓ i k ε^±_i(𝐤),

and satisfy the relations

    ε^λ_i(𝐤) ε^λ'_i(𝐤)^∗   = δ^λλ', 
    
    	k_i ε^±_i(𝐤)
           = 0, 
    ε^±_i(-𝐤)
           = ε^±_i(𝐤)^∗, 
    ε^±_i(𝐤)^∗   = ε^∓_i(𝐤).

At linear order, the equation of motion, <ref>, reduces to

    A_±”(𝐤) + k ( k ∓ 2 ℋξ) A_±(𝐤) = 0,

where we have defined the instability parameter

    ξ≡α/2 fϕ'/ℋ.

<Ref> shows that the interaction between the inflaton and the gauge fields induces a tachyonic instability (for modes with k/ℋ < 2 |ξ|) for one of the two helicity states (A_+ if ϕ' > 0, or conversely A_- if ϕ' < 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., <cit.>).
Assuming ϕ' > 0, near horizon crossing (k / a H ∼ 1) these solutions are well approximated by

    A_+ ≃1/√(2k)( k/2 ℋξ)^1/4 e^πξ - 2 √(2 ξ k/ℋ).

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

    1/2⟨𝐄^2⟩   =1/4 π^2 a^4∑_λ=±∫ k   k^2| A_λ^'(k)|^2
    1/2⟨𝐁^2⟩   =1/4 π^2 a^4∑_λ=±∫ k   k^4| A_λ(k)|^2, 
    ⟨𝐄·𝐁⟩   = - 1/4 π^2 a^4∑_λ=±λ∫ k   k^3∂_τ| A_λ(k) |^2,

where A_λ(k) denotes the helicity modes of the gauge potentials.
(Note that <ref> assume the temporal gauge.)
For ξ≳ 4 <cit.>, fairly accurate approximations of these quantities are given by <cit.>

    ⟨𝐄·𝐁⟩   ≃ 2.4 × 10^-4 ℋ^4/a^4 ξ^4 e^2 πξ, 
    1/2⟨𝐄^2 + 𝐁^2⟩   ≃ 1.4 × 10^-4 ℋ^4/a^4 ξ^3 e^2 πξ.

Substituting <ref> into <ref> demonstrates that gauge fields exert an additional friction term for the inflaton (see <ref>).
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 <ref> we present our procedure to accurately account for these effects during inflation (as used to set initial conditions for the lattice simulations).





§ INITIAL CONDITIONS


As implemented in Refs. <cit.>, 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, <ref>, during inflation.
These fluctuations are integrated alongside the background dynamics, i.e., <ref>, and the homogeneous part of the inflaton's equation of motion, <ref>.
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 <ref>.
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 <ref>), Refs. <cit.> use the approximations <ref>.
The accuracy of these expressions is sensitive to two main assumptions.
First, the solution of <ref> is well approximated by the Whittaker function only in the limit that ξ is constant, and ℋ = -1/τ (in near–de Sitter space).
Second, <ref> themselves are approximations to the integrals <ref> over <ref>, accurate for ξ≳ 4 <cit.>.
When ξ is nearly constant but small, the approximate integrals <ref> 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 <ref> 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 <ref> 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, ξ < 4, meaning the approximation to the analytic integral is not trustworthy <cit.>.
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 <ref> via numerical quadrature of the numerically evolved gauge fields.
As chosen for the analytic results <ref>, we take as a UV cutoff for these integrals k = 2 a H ξ, 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 amplified—the band of tachyonic instability k ≲ 2 a H ξ is proportional to the comoving horizon ℋ = a H, 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 to 10 orders of magnitude) logarithmically, yielding good results from quadrature of <ref> 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 ≫ 2 a H ξ) is precisely the regime where the Wentzel-Kramers-Brillouin (WKB) approximation

    A_±(k, τ) ≈1/√(2 ω_±(k, τ)) e^i ω_±(k, τ) τ

with ω_±(k, τ)^2 = k^2 ∓ 2 a H ξ k is extremely accurate.
Thus, we only numerically integrate a given mode k when k ≤ 50 × 2 a H ξ, using the WKB solution until this point.
Using an adaptive ordinary differential equation integrator (in practice, 's DOP853 routine <cit.>) 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. <cit.> 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 <ref> we investigate the difference between the analytic approximations [<ref>] and the numerical computations described above, fixing a chaotic inflaton potential.

We consider couplings where the background evolutions agree well (α / f = 8  ^-1) and disagree (α / f = 14  ^-1).
We additionally compare to the numerical quadrature of <ref> using the actual Whittaker-function solutions, which demonstrates agreement with the approximations <ref> once ξ≳ 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, <ref> indicates that the analytic approximations greatly overestimate ⟨𝐄^2 + 𝐁^2 ⟩ and ⟨𝐄·𝐁⟩ 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) in <ref>.
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  ^-1 is roughly the largest coupling for which the background evolutions visibly agree.





§ POST-PREHEATING DYNAMICS


In this appendix we study the dynamics after preheating to evaluate whether the Universe remains radiation 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 <ref>, at the end of all of our simulations a small amount of the total energy remains in axion fluctuations [𝒪(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

    Ω_gw≡ρ_gw/ρ∝1/(a / a_0)^4 (H / H_0)^2∝( a/a_0)^3 w - 1,

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.[More precisely, the equation of state of a coherently oscillating scalar field time averages to zero.]
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 + 2 f_A(m_ϕ / 2) <cit.>.

For simplicity, we study the dynamics after preheating with the standard Boltzmann equations describing three-body decay of inflaton particles at rest into two relativistic daughter particles <cit.>,

    ρ_ϕ/ t + 3 H ρ_ϕ   = - Γρ_ϕ
    ρ_γ/ t + 4 H ρ_γ   = Γρ_ϕ.

Rather than fixing the decay rate Γ and the initial (i.e., post-preheating) 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 <ref> we numerically integrate <ref> 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 <ref>).

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 <cit.>

    Γ_0/H = 1/Hα^2 m_ϕ^3/64 π f^2∼ 10^-9.

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 <cit.>.
Thus, the relevant portion of <ref> is 10^-2≲Ω_ϕ≲ 10^-1 and Γ≳ 10^-1, for which the fraction energy density in 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, <ref> 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 to 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.






