Gaugepreheating and the end of axion inflation
Abstract
We study the onset of the reheating epoch at the end of axiondriven inflation where the axion is coupled to an Abelian, , gauge field via a ChernSimons interaction term. We focus primarily on inflation and explore the possibility that preheating can occur for a range of coupling values consistent with recent observations and bounds on the overproduction of primordial black holes. We find that for a wide range of parameters preheating is efficient. In certain cases the inflaton transfers all of its energy to the gauge fields within a few oscillations. In most cases, we find that the gauge fields on subhorizon scales end preheating in an unpolarized state due to the existence of strong rescattering between the inflaton and gaugefield modes. We also present a preliminary study of an axion monodromy model coupled to gauge fields, seeing a similarly efficient preheating behavior as well as indications that the coupling strength has an effect on the creation of oscillons.
a]Peter Adshead b,c]John T. Giblin, Jr. b]Timothy R. Scully a]Evangelos I. Sfakianakis Prepared for submission to JCAP
Gaugepreheating and the end of axion inflation

Department of Physics, University of Illinois at UrbanaChampaign, Urbana, Illinois 61801, U.S.A.

Department of Physics, Kenyon College, Gambier, Ohio 43022, U.S.A.

Department of Physics, Case Western Reserve University, Cleveland, Ohio 44106, U.S.A.
Contents
1 Introduction
Inflationary model building has entered a particularly exciting phase with the demonstration by the BICEP2 experiment [1] of the sensitivity to Bmode polarization of the cosmic microwave background (CMB) at a level where interesting constraints can be, or soon will be, placed on the inflationary model space. The recent Planck results on dust emission [2] combined with the joint Planck and BICEP2/KECK Array analysis [3] mean that dustforegrounds need to be accurately characterized in order to determine if any of the Bmode polarization observed by the BICEP2 experiment is due to primordial gravitational waves. If confirmed, the observation of these gravitational waves present a spectacular confirmation of one of the early observational predictions for slowroll inflation [4]. However, they also present challenges for inflationary model building. The large primordial gravitational wave amplitudes required to explain the BICEP2 signal generically require that the scalarinflaton field rolls over a distance in field space that is large compared to the fourdimensional Planck scale [5]. This makes the theory extremely sensitive to unknown physics at short wavelengths, e.g. in the ultraviolet (UV), at or near the Planck scale, and leads to a loss of predictive power, if not a loss of the inflationary mechanism itself. A possible way around this UV sensitivity problem is to find a symmetry powerful enough to forbid interactions between the sector driving the inflationary expansion and other unknown physics.
A promising candidate for the inflaton field is a pseudoscalar or axion. These fields enjoy shiftsymmetries that protect their role as inflatons from being spoiled by coupling to unknown UV physics. Shiftsymmetries require that the theory is invariant under a constant shift of the field value, and severely restrict the form of possible interactions with other fields. One of the earliest proposed axion inflation models was Natural inflation [6, 7]. In this scenario, a cosine potential for the axion is generated by the condensation of a nonAbelian gauge group. Slowroll inflation is achieved by the hierarchy between the height and width of the potential, the form of which is protected by the shiftsymmetry. In order to generate density fluctuations with a spectrum that reproduces the observed fluctuations in the CMB, the axion is required to have a periodicity, or associated mass scale, larger than the Planck scale. This superPlanckian periodicity makes it extremely difficult to embed Natural inflation in a fundamental theory such as string theory [8]. Model builders sought to circumvent this obstruction in a model dubbed flation [9, 10, 11] by using a large number of axions, each with small periodicities, to implement assisted inflation [12]. While each of the axions only rolls a small distance in field space, their collective motion is responsible for inflation and effectively traverses a large distance. Unfortunately the axions in flation have the effect of renormalizing the Planck scale to a lower value, and the theory ultimately suffers from similar pathologies as the original Natural inflation model it was supposed to fix. Other formulations making use of misaligned axions have also been proposed [13, 14, 15]. More recently, the observation that a single axion undergoing monodromy can have a small periodicity while ultimately traversing a large field range [16, 17] has seen resurgent interest in axion inflation [18, 19, 20, 21, 22, 23]. For a recent review of axion inflation see Ref. [24] and its realizations in string theory see Ref. [25].
At the end of the inflationary phase, the Universe must undergo a phase transition, from its supercooled state to a state filled with radiation and ultrarelativistic matter, to begin the hot big bang. The physics of this phase transition is thought to be highly nonlinear, and its details are unknown (see e.g. Ref. [26] for a recent review). The shiftsymmetry in axiondriven inflation that is so effective at protecting the form of the inflationary sector from unknown UV physics also severely constrains the form of its couplings to the visible sector. These couplings to the standard model of particle physics, either directly or indirectly via intermediaries, are required in order to transfer the inflaton energy into radiation and ultrarelativistic matter. The shift symmetry dictates that couplings to matter fields must be derivative interactions, therefore a class of allowed interactions are those in which the axion is coupled to a gauge field through a dimensionfive ChernSimons term or Pontryagin density. A coupling of this form is allowed by the symmetries and, from the viewpoint of an effective field theory, must be present. Furthermore, this coupling provides a perturbative decay channel for the axion into gauge bosons which guarantees that reheating eventually completes through perturbative decays alone, and thus provides a viable pathway through which the Universe can transition from inflation into the hot big bang.
The effects of the coupling of axions to gauge fields during inflation is, by now, a well studied field. It has been known for a long time that axiongauge field couplings lead to parityviolating gaugefields that are amplified during slowroll inflation [27, 28]. The behavior of these gauge fields during inflation and their influence on the inflationary dynamics has also been extensively studied [29, 30, 31, 32, 33, 34]. Further, the authors of Refs. [35, 36] studied metric fluctuations generated by a rolling auxiliary pseudoscalar during inflation. Some work has also been done on the reheating of axiondriven inflation. The authors of Refs. [37] considered stringy reheating of a monodromy scenario while Ref. [21] studied perturbative decay of an axion inflaton to photons/gauge bosons. Furthermore, previous studies have considered perturbative analyses of the Mathieu equation for gauge fields coupled to an oscillating scalar or collection of scalars [38, 39]. However, these studies were focussed primarily on ranges of parameters considered natural for scenarios of Natural inflation and flation and concluded that highly nonlinear effects and parametric resonance were unimportant in these models.
Recently, larger couplings between the axion and gauge sectors have been considered [30, 32], which can lead to observable effects during inflation due to the rescattering of the gauge fields off the axion condensate [40]. Effects such as nonGaussianity of the density fluctuations, chiral gravitational waves, and the production of primordial black holes [40, 31, 32, 41, 42] place upperbounds on the strength of the couplings, the most stringent being due to restrictions on the production of primordial black holes at the end of inflation.
Our current work addresses a specific problem; given that highscale inflation requires highly constrained couplings, can the Universe undergo a preheating phase (either tachyonic or resonant) following axion inflation? Numerical investigations of reheating in canonical [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55] or noncanonical [56] scalar field scenarios is becoming routine. Furthermore, couplings of scalar fields to Abelian [57, 58, 59] or nonAbelian gauge fields in cosmological settings have been studied in recent years [60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75]. The question we address here is new. We employ lattice simulations, using the same numerical technique of Ref. [59], to study the possibility of tachyonic or parametric amplification of gauge fields after inflation due to the presence of a ChernSimons term. As in most studies of preheating, we remain agnostic as to the specific nature of the coupled gauge field and do not necessarily identify it as a gauge field from the standard model.
Our results can be summarized as follows. We find that for reasonable ranges of the axiongauge field coupling, nonlinear effects can be very important at the end of inflation. In particular, at the middle to upper range of the couplings allowed by black hole abundance, we find that reheating is essentially instantaneous, proceeding via a phase of tachyonic resonance [76] and completing within a single oscillation of the axion. Despite the asymmetry in the equations of motion for the two polarizations of the gauge fields, on subhorizon scales, rescattering of the gauge bosons off the axion condensate is efficient at generating the second polarization. On scales larger than the horizon at the end of inflation, an asymmetry between the gauge field polarizations remains. The Universe that results in these cases is radiation dominated and is characterized by a very high reheating temperature. As the coupling is decreased, this tachyonic resonance is weakened and the axion oscillates multiple times before reheating completes. During these multiple oscillations equal levels of both polarizations of the gauge field are excited. Decreasing the coupling further yields a brief window where parametric resonance effects become important before preheating abruptly shuts off and nonlinear effects cease to be important. At these lower couplings, nonlinear effects are negligible and the Universe reheats via perturbative decay of the axion into gauge bosons.
We also investigate how these preheating effects might depend on the shape of the inflationary potential. As a second test, we subject the axion to a monodromytype potential [16], and show that the range of couplings for which efficient preheating can occur is slightly different, although the same orderofmagnitude. We conclude by presenting an intriguing set of data that suggest gauge fields might play a role in the creation and evaporation of oscillons in this scenario.
We work in natural units where , however, we retain the Planck mass, .
2 Background and conventions
We begin with the usual action for axiondriven inflation
(2.1) 
where is a pseudoscalar (axion) and is a potential that supports slowroll inflation. For definiteness, we consider the potentials for the simplest type of chaotic inflation [77],
(2.2) 
and the simplest type of axion monodromy inflation [16]
(2.3) 
which is well described by a linear function of for large field values. The amplitude of the scalar spectrum fixes the parameters and to be [78, 79, 80]
(2.4) 
and
(2.5) 
The parameter in the monodromy potential, Eq. 2.3, has a negligible effect on the spectrum of curvature fluctuations and gravitational waves on the scales that are observable in the CMB (provided, of course, that is much smaller than the field value where the CMB fluctuations are generated) and, consequently, is unconstrained by data. However, becomes important near the end of inflation and has a small effect on the value of when inflation ends. Further, for small field values, , the potential can be expanded as
(2.6) 
and the resulting dynamics of in this region depend strongly on the value of [81].
In addition to the axion, we consider a gauge field coupled to the axion
(2.7) 
where is a dimensionless coupling constant of order unity and is a mass scale associated with the pseudoscalar (axion). The field strength and its dual are given by the standard expressions, and , where is the completely antisymmetric tensor, and our convention is
(2.8) 
Greek letters here and throughout denote four dimensional spacetime indices and Roman letters from the middle of the alphabet are used to denote spatial indices. Repeated lower spatial indices are summed using the Kronecker delta. We work with the FriedmannLemaîtreRobertsonWalker (FLRW) metric in conformal time with mostlyplus conventions
(2.9) 
The equation of motion for the pseudoscalar field is the KleinGordon equation sourced by the ChernSimons density of the gauge field
(2.10) 
where and where here and below . The equations of motion for the gauge field are
(2.11) 
The equation is the Gauss’ law constraint
(2.12) 
while the equations are the field equations for the spatial components of the gauge field
(2.13) 
Finally, assuming the metric is unperturbed, the scale factor satisfies Einstein’s equations
(2.14) 
and
(2.15) 
The pressure, , and energy density, , are found from the stressenergy tensor
(2.16) 
which can be explicitly written as
(2.17) 
and
(2.18) 
Note that the axiongauge field coupling does not contribute directly to the stressenergy tensor.
3 Gaugefield production during inflation
During inflation, the coupling of the gauge field to the axion results in exponential production of one polarization of the gauge field over the other [27, 28, 29]. To see this, we first fix the gauge by choosing Coulomb, or transverse, gauge . The Gauss’ law constraint, Eq. 2.12, then implies that at linear order in fluctuations.^{1}^{1}1Note that at this order (linear) in fluctuations, Coulomb gauge and temporal gauge () are equivalent. This can be seen trivially from Eqn. (2.12), the Gauss’ law constraint. In the linear regime, this equation reads . In Coulomb gauge, assuming , this constraint reads . In temporal gauge, Gauss’s law reads , which for , implies .
At linear order in fluctuations, with this choice of gauge, the equation of motion for the gauge field becomes
(3.1) 
To study the fluctuations of the gauge field, it proves most convenient to work in Fourier space, where our convention is
(3.2) 
The Fourier components are then expanded in a basis of helicity states
(3.3) 
where the polarization vectors satisfy the orthogonality and normalization relations
(3.4) 
In this situation, conformal time is defined to be a negative, increasing quantity during inflation
(3.5) 
where the last approximation is exact in the de Sitter limit, , where the slowroll parameter, is defined as . Here and throughout, an overdot is used to denote a derivative with respect to cosmic time, .
We can now quantize the modes by introducing the creation and annihilation operators, and satisfying the canonical commutation relations
(3.6) 
which allows us to expand the modefunctions as
(3.7) 
With our conventions, the gauge field equation of motion Eq. 2.13 becomes a separate equation for each polarization, depending only on the magnitude of the momenta
(3.8) 
where we have also made use of the de Sitter approximation for the scale factor during inflation in the last term. During inflation, and, after changing variable to , the equation of motion is transformed to the Whittaker equation
(3.9) 
In our case, we have
(3.10) 
and thus and where we have defined
(3.11) 
The general solution of the Whittaker equation can be written in terms of the Whittaker Wfunction
(3.12) 
The constants of integration, and , are set by canonical quantization which amounts to normalizing the modes according to the Wronskian condition
(3.13) 
and demanding that the modefunction approaches the Minkowski vacuum in the limit,
(3.14) 
The properly normalized solutions are
(3.15) 
In the limit , the asymptotic form of the mode function is
(3.16) 
Compared to the conformally invariant radiation solution, the circularly polarized modes get amplified by a factor
(3.17) 
where we have used the Stirling formula
(3.18) 
and we have assumed that . Note that this means that when (), the mode () gets amplified by a factor while the other mode is unchanged. For this work, we focus on largefield inflationary models, and assume that so that the mode is amplified during inflation.
The exponentially enhanced gauge fields have important effects during inflation due to their rescattering off the inflaton condensate and their interactions with the metric. The former leads to the production of fluctuations of the inflaton which are statistically nonGaussian, while the latter leads to the production of gravitational radiation [32]. The Gaussianity of the observed density fluctuations by Planck [82] then implies that the quantity [24], where is the quantity in Eq. 3.11 evaluated during the time when the modes that form the CMB leave the horizon.
During inflation, and thus the ratio increases. This means that shorterwavelength modes that leave the horizon later in inflation are amplified more than their longerwavelength counterparts that leave the horizon earlier. The largest effects occur when is near unity near the end of inflation. In the limit that (which for models that satisfy only possibly occurs near the end of inflation) the energy density in the gauge fields becomes important and the gaugefield fluctuations begin to backreact on the homogeneous background equations of motion. In this limit, in the Hartree approximation, the Friedmann (Eq. 2.14) and KleinGordon equations (Eq. 2.10) become
(3.19) 
(3.20) 
and
(3.21) 
where the electric and magnetic fields^{2}^{2}2We refer to these fields as electric and magnetic, however, they need not be the electromagnetic fields of the standard model. associated with the gauge field are and . In this limit, up to an irrelevant constant phase, the gauge field mode that is amplified is approximated near horizon crossing by [30]
(3.22) 
while the other mode is unaffected and is negligible. The expectation values of the quantum fields are well approximated by [30]^{3}^{3}3Note that the axion velocity assumed here is opposite to that assumed by [30]. This means relative to this work, the other gauge mode is amplified and consequently the sign of is opposite.
(3.23) 
Toward the end of inflation, for large values of , the back reaction of the gauge fields on the rolling axion becomes important and inflation is prolonged [32]. During this phase, the primordial density fluctuation spectrum is expected to be dominated by rescattering and large, nonGaussian density fluctuations are predicted. These large amplitude density fluctuations can produce primordial blackholes [83, 84, 85] and ensuring that they are not overproduced requires that for [41, 42] and for monodromy [42], which is tighter than the current bounds from the Gaussianity of the CMB fluctuations. These limits are model dependent, and are somewhat sensitive to the form of the potential. In this work we consider only values of the coupling such that . This bound translates to roughly for the potential and for the simple monodromy potential of Eq. 2.3. In practice, we do not approach this threshold in our simulations, keeping the couplings used for our simulations lower by around a factor of two.
The bounds from primordial blackhole abundances rely on the above approximations of Eqs. 3.19  3.21 together with Eq. 3.23 being a good description of the system in the region of strong back reaction [41]. However, these approximations do not yield a self consistent set of equations because the resulting stressenergy tensor is not covariantly conserved. This means that in the region where back reaction becomes significant, the above equations are inaccurate. For small couplings, we expect that the above approximations are accurate enough to capture the onset of back reaction. In what follows we initialize our lattice simulations using the values of the fields found from the numerical evolutions of the Eqs. 3.19  3.21. We make use of the approximations of Eq. 3.23 for the expectation values of the energy density in gauge field fluctuations and the Pontryagin density respectively. The use of these approximations implies that the initialization of our simulations is less and less accurate for larger values of the coupling between the gauge field and the axion. To minimize this error, we initialize our simulations two efoldings before inflation ends. The modes that are important in the reheating era are well within the horizon at this time, and well described by the linear approximation at the start of our simulations (See Table 1 and Fig. 6).
4 Perturbative reheating
Before we move to the nonlinear regime, we briefly revisit the perturbative reheating case. The coupling of the axion to gauge fields provides a natural decay channel for the axion to produce gauge bosons. Even in the absence of nonlinear effects due to the homogenous motion of the inflaton condensate, the Universe will reheat to gauge bosons via perturbative decays of the inflaton into a pair of gauge bosons. When the Hubble rate becomes comparable to the perturbative decay rate, the decay of the inflaton into gauge bosons proceeds rapidly. The contribution of these newly created particles to the energy density quickly dominates and reheats the universe. This represents a lower bound on the temperature of the Universe at reheating.
The coupling allows the axion to decay to two gauge bosons with a rate that is well known (see e.g. [86])
(4.1) 
where is the mass of the axion about its vacuum. This rate sets a lower bound on the reheating temperature which is found by comparing the decay rate to the Hubble rate; perturbative reheating finishes when , which results in a reheating temperature
(4.2) 
where is the number of relativistic degrees of freedom. For the simplest version of chaotic inflation, we can write
(4.3) 
For small values of the coupling, this is the dominant effect. Reheating in this case proceeds by perturbative decay of the axion into gauge bosons.
5 Instabilities and resonance
We are interested in what happens immediately following the inflationary epoch. We assume that the pseudoscalar begins to oscillate about the minimum of its potential, which we initially take to be a quadratic function, Eq. 2.2. While this form is exact for the simplest model of inflation, there are anharmonic corrections to the potential for other important axion inflation scenarios, such as monodromy inflation.
It is instructive to first consider the behavior of the system by temporarily neglecting the expansion of space and considering only the linear theory to gain intuition for the system and to identify features present in the fully nonlinear treatment. In this regime the axion satisfies the equation
(5.1) 
where overdots denote derivatives with respect to cosmic time . This equation has the simple solution
(5.2) 
where is (approximately) the field value at which the slowroll conditions are violated and inflation ends. To estimate the effect of parametric resonance (while still neglecting the expansion of space) we write the equation of motion for the mode amplitudes, Eq. 3.8, as
(5.3) 
With the solution for from Eq. (5.2), and after redefining time , this equation can be recast as
(5.4) 
The fact that each helicity obeys a different equation is irrelevant here because the sign of the term proportional to in Eq. 5.4 can be reversed by a constant shift of its argument . In this approximation both modes are expected to grow equally after each inflaton oscillation. We show how this symmetry between the two helicities is broken once the expansion of the Universe is taken into account.
We can compare Eq. 5.4 to the normal Mathieu equation,
(5.5) 
to see that^{4}^{4}4This definition of the two parameters of the Mathieu equation as and is typically used in models where the inflaton decays to scalars. In these models does not depend on the wavenumber. In the case of gauge fields does depend on , but we refrain from using a subscript to be consistent with prior literature.
(5.6) 
From Eq. 5.5 we see two important thresholds that define the behavior of the solutions. These thresholds are the tachyonic resonance threshold, which is set by , and the broadtonarrow resonance threshold, which is determined by the size of . Broad resonance occurs for while narrow resonance occurs for . Fig. 1 shows the Mathieu instability diagram for our process, along with three curves for fixed coupling, , and varying wavenumber. The line is the diagonal for our choice of axes range. It can be clearly seen that in the regime where , the instability bands are much broader and the imaginary parts of the Mathieu exponents are larger. It is also interesting to note that this system can be cast in terms of only two (dimensionless) combinations: the ratio of wavenumber to mass scale, , and the product of the coupling strength and initial axion amplitude, . The curves in Fig. 1 are defined by a fixed coupling and initial field amplitude , and are parameterized by wavenumber.
When the combination , which occurs for , there is a tachyonic instability in the equation of motion [76]. This condition defines a set of modes
(5.7) 
whose masssquared is negative ( from Eq. 5.5). There is always a set of such modes, as long as the homogeneous mode of the axion is oscillating.
The details of this tachyonic regime are seen directly from the equations of motion for the two polarizations, Eq. 5.3. When the axion velocity changes sign, the combination
(5.8) 
is negative for only one of the two polarizations, depending on the sign of . When the axion velocity is positive (negative), the () mode is tachyonically amplified. This regime is obviously most important for large couplings and when the amplitude of the field oscillation is large. In these large coupling cases, preheating is extremely efficient and can complete after only one or two oscillations of the homogeneous inflation condensate. These tachyonic instabilities disappear when the homogeneous mode of the inflation breaks down, which occurs due to back scattering of the gauge modes, or self resonance of the axion itself (when the axion potential includes anharmonic terms). This should happen very quickly for large couplings – in some cases, only one polarization is amplified by this tachyonic resonance. However, as we show, in these cases rescattering or backscattering effects are extremely efficient and generate equal amounts of the nontachyonic polarization. As can be seen from Fig. 1, the wavenumbers corresponding to the tachyonic regime (Eq. 5.8) lead to much larger Mathieu exponents, so they dominate the behavior of the system, hence we concentrate on them.
Having discussed the behavior of Eq. 5.5 for , we move to the second threshold which is defined by the size of . If the curve (with constant ) intersects the Mathieu bands for values of , we are in the broad resonance regime [87, 76] characterized by the curve intersecting large instability bands. In the opposite regime, , the modes that are amplified have frequencies comparable to those of the oscillating inflaton. This is called the narrow resonance regime, since the instability bands of the Mathieu chart have a narrow width. The growth rate of gauge field modes in this regime can be fully analyzed using the methods of parametric resonance based on Floquet theory, as described for example in [87] and applied to narrow resonance for gauge fields in [38]. Since narrow parametric resonance is only present for small values of the coupling where the Universe does not completely preheat, we do not consider it for the remainder of this work.
Previous studies of reheating after axion inflation inflation through gaugefield production, such as [38], focused on the region where , and seem to have missed the efficient, tachyonic preheating phenomenon on which we focus here. As can be seen in Fig. 2, both the range of wavenumbers, as well as the size of the Floquet exponent are much larger for the range of couplings we study.
These tachyonic instabilities exist for long wavelength modes of the gauge field, but the assumption of a homogeneous inflaton, , breaks down quickly once backreaction occurs. To see whether rapid production of gauge field quanta occurs, and/or whether this final state is polarized, quickly becomes a numerical question. However, there is some progress to be made using semianalytic methods by reintroducing the expansion of the Universe into the equation of motion for the gauge fields, as explained in the next section.
5.1 Semianalytic treatment
A detailed description of tachyonic resonance in the staticuniverse approximation can be found in [76], where it is shown that the agreement between these analytic results and numerical simulations (neglecting rescattering and back reaction) is excellent (for and ). In our current study, we go beyond this approximation. The richest phenomenology comes from the period between the end of inflation and the first few oscillations of the inflaton. During this epoch, the scalefactor is evolving nontrivially in time and cannot be approximated by an exponential, as it could during slowroll inflation, or a powerlaw, as it is in a radiationdominated universe after reheating. In order to proceed with a semianalytical treatment, as an approximation, we model the evolution of the Universe during this time by the evolution of the classical axion in a background FRW spacetime, neglecting the effect of back reaction of the gauge fields on the expansion. In Fig. 3 we plot the evolution of the inflaton field, its velocity, and the Hubble parameter in units of the inflaton mass for simple chaotic inflation. Without loss of generality, consistent with Section 3, we chose the value of the axion to be positive during inflation, , and its time derivative to be negative, (note that here ). This then determines the signs of these two quantities immediately at the end of inflation. Note that during the first oscillation the amplitude falls by about a factor of two, and thus approximating the inflaton background by a sinusoidal function with a constant amplitude is insufficient.
We need to extend the method of [76] in order to treat the case of a timedependent but nonharmonic effective frequency. We are most interested in the cases where the growth of the gauge field modes is due to tachyonic effects near the end of inflation. In these cases, preheating is extremely efficient, and we focus our attention on the growth of the gauge field modes during the first oscillation of the axion. Due to the strong growth of the fluctuations, analytical methods quickly fail and we enter the regime of strong backreaction, dominated by nonlinear physics. In this regime we have to rely on full numerical simulations, which we present in the following section. Extending the (linear) analysis beyond the first axion oscillation is straightforward, and can be achieved by combining the procedure we describe here with the method of [76].
We start from Eq. 3.8 written in cosmic time
(5.9) 
and redefine the gauge field,
(5.10) 
which leads to the equation of motion for this rescaled gauge field
(5.11) 
We use the WentzelKramersBrillouin (WKB) approximation to describe the solution to Eq. 5.11 in the regions where the frequency
(5.12) 
is varying slowly, that is, . Following the general procedure of the WKB approximation, we distinguish three regimes based on the behavior of the effective frequency. We track a mode as it goes from having a positive frequencysquared, , (regime I) to a negative frequencysquared (regime II), and then back to a positive frequencysquared (regime III). At the interface of these regions, the frequency vanishes (), and the condition is maximally violated. In practice, this condition is violated for some time interval around the points where the frequencysquared changes sign.
Assuming that the frequency is slowly varying, we begin by writing the lowestorder WKB approximation to the solution of Eq. 5.11 in each region [88]
(5.13)  
(5.14)  
(5.15) 
where , and and the integrals in the expression for are performed for the region where . We are ultimately interested in the two coefficients and which describe the amplitude of the mode after its brief growth due to the tachyonic instability. In regime I we match the solution to the asymptotic past (BunchDavies vacuum), giving us and . After matching at the two interfaces between regimes I and II and regimes II and III (by Taylorexpanding near the turning points and using the asymptotic form of the Airy functions) the coefficients are
(5.16) 
where
(5.17) 
Before proceeding, we note that the WKB approximation is accurate only in the broad resonance regime, . This translates to the condition that
(5.18) 
For typical values of the coupling and the inflaton amplitude, this condition restricts the validity of the WKB approximation to . As shown in Fig. 5, the maximum amplification occurs well within the region of validity of our WKB calculation, especially for higher couplings.
We study each polarization individually, starting with the mode . In the absence of back reaction, the axion velocity, , takes its maximum positive value at time after inflation ends,^{5}^{5}5In this section we denote by a subscript the value of a quantity at the end of the inflationary phase where for the first time. which sets the interval during which is tachyonic for the first time. The evolution of the inflaton field is close to sinusoidal at this stage, meaning that we can use a modified staticuniverse approach. However, we employ our full expandinguniverse WKBmethod as a way of testing its accuracy and providing a unified treatment of the two gaugefield polarizations.
The evolution of the second polarization, , is more involved due to the fact that regime II, characterized by , starts while the Universe is still inflating and continues through the end of inflation into the first oscillation of the inflaton. This complication does not significantly alter our analysis. To proceed, we simply need to initialize the mode in regime I, sufficiently early during the inflationary stage, before its tachyonic transition and follow it, using the WKB approximation, through this transition and the end of inflation. The final formulas are exactly the same in this case.
We can evaluate the validity and accuracy of the semianalytical method described above by comparing the results directly with a full numerical solution of the linear equations of motion. Using Mathematica, we solve the homogeneous KleinGordon and Friedman equations for the evolution of the background axion and spacetime. On this background, we follow three gaugefield modes for whose wavelengths are equal to the horizon, and one and two efoldings smaller than the horizon at the end of inflation, i.e. respectively at the end of inflation. We also track three modes for , one that exits the horizon one efold before the end of inflation, a mode whose wavelength is equal to the horizon, and a third that is one efold smaller than the horizon at the end of inflation, i.e. respectively at the end of inflation. In order to facilitate the comparison, we take the mode amplitudes to have unit size at the start of our simulation. The WKB condition is violated around the points where for each mode and during the tachyonic regime. This limits the accuracy of the approximation. The approximation could be improved by making use of a transformation of variables similar to [89], however, given the good agreement with the numerical results shown in Fig. 4, and the fact that we are trying to understand the results of full lattice simulations rather than substitute them, we do not attempt to refine the WKB method used here.
In Fig. 4 we show the excellent overall agreement between the numerical results from direct numerical solution in Mathematica and those obtained from the WKB approximation. As expected the approximation diverges at the points where , and closely follows the curves as one moves away from these points. Furthermore, we can see that for both and the accuracy of the approximations decreases with decreasing wavenumber (leading to decreasing ), especially after the tachyonic transition. Although we are pushing the limits of the WKB approximation, the accuracy of the resulting growth rates is remarkably robust.
These results demonstrate that this semianalytical method can be used to accurately estimate the growth of the gauge fields during the first inflaton oscillation, if one neglects rescattering and backreaction effects. In Fig. 5, we plot the growth factor , which shows by how much the amplitude of each gauge mode has grown after the first tachyonic regime. Note that this WKB method breaks down at small wavenumbers.
For the mode we can make a comparison to the staticuniverse calculation by using rescaled parameters. The amplitude of the axion oscillations has decreased due to the expansion of the Universe after the end of inflation (as shown in Fig. 3). For the tachyonic halfperiod of interest, the behavior of the axion field is well approximated by
(5.19) 
where is a phase offset that allows the time of zerocrossing to correspond to our model. The wavenumbers of the modes under consideration are also redshifted. The modes used for our full WKB calculation were measured in terms of the scale factor at the end of inflation , while the Universe has grown by a factor of by the middle of the first tachyonic regime for . We thus rescale the wavenumber used in our staticuniverse calculation by . As we can see in Fig. 5 the results are very close, giving us a simple physical way to understand the result of using the WKB method in an expanding universe.
The polarization is more complicated due to the fact that it becomes tachyonic during the inflationary phase. From Fig. 5 note that the growth factors for a given mode are significantly larger, and a larger range of wavenumbers get amplified. We do not compare this case with a staticuniverse approximation, since is not well approximated by a sinusoidal function in the regime where is tachyonic, even after inflation has ended. This is due to the Hubble friction term, as shown in Fig. 3.
There is one further complication we need to keep in mind when comparing the semianalytic predictions of this section with the full lattice simulations presented in Section 6. Once the coupling to the gauge fields is turned on, these fields begin to backreact. This causes inflation to end at a slightly different value of for each coupling , as shown in Table 1. We also show the values of the field and field velocity two efoldings before inflation ends. Note that these values are all close to the zero coupling case, which indicates that backreaction has not yet become important yet.The WKB analysis of this section was performed using the evolution of the axion in the limit when we neglect back reaction from the gauge fields, equivalently in the limit of .
6 Nonlinear lattice simulations
The results of Section 5.1 imply that at the couplings of interest, the dynamics of the axiongauge field system very quickly become highly nonlinear and require numerical analysis. Numerical methods that evolve scalar fields in an expanding background have been around for a couple of decades. Many numerical methods have been developed [90, 91, 92, 93], and all have regimes in which they are most successful. Here we use GABE [94] (as in [59]), where the axion and the gauge field are defined on a discrete lattice (grid) with points. The software uses a secondorder RungeKutta integration method to solve Eqs. 2.10, 2.12 and 2.13 alongside the self consistent expansion of spacetime Eq. 2.14. We work in Lorenz gauge, , to evolve the gauge fields. In this gauge, the Gauss law constraint becomes a dynamical equation for which we solve in parallel with Eqn. 2.13. Note that, although we initialize our fields using solutions of the linear equations of motion in Coulomb gauge, these gauge fields are equivalent to the fields in Lorenz gauge. To see this, note that if one begins with fields in Coulomb gauge, and performs a gauge transformation to Lorenz gauge , the Lorenz gauge condition () becomes , which is simply the residual gauge symmetry of the Lorenz gauge. Unless otherwise noted, all simulations use a box size that is at the end of inflation and are run using the parallel processing standard OpenMP with 12 threads. The equations of motion for the gauge degrees of freedom, , as well as the axion, , are integrated using a secondorder (midpoint method) RungeKutta integration scheme. We normalize quantities so that wave numbers are expressed in terms of the size of the box at the end of inflation, , and (at the end of inflation).
We begin our simulations two efoldings before inflation ends, defining this point as . As described at the end of Section 3, we determine the value of the homogeneous field and its derivative by numerically evolving (using Mathematica) the system of Eqs. 3.19  3.21 together with the approximations of Eq. 3.23. At this point the boxsize is which is just larger than the Hubble Scale, , where the final approximation varies slightly from coupling to coupling. We initialize the power in the modes by numerically evolving Eq. 3.8 for each physical mode, tracking it from when it was well inside the horizon () until two efoldings before the end of inflation. Since we have no analytic form for the fluctuations of at this time, we initialize it in the BunchDavies vacuum, , where it is an excellent approximation due to the fact that almost all of our modes are subhorizon. Using this procedure, the modes that leave the horizon during the final two efoldings of inflation are done so selfconsistently and the spectrum of perturbations for is consistent with our equations of motion. Fig. 6 shows the spectra at the beginning of the simulation and the end of inflation showing the amplification of large wavelength modes during the final stages of inflation.
After we set the initial spectrum of in momentum space we project these onto the gauge fields
(6.1) 
and (inverse) Fourier transform them into configuration space using a set of projection operators, , that satisfy the relations
(6.2) 
These relations set only the spatial components of the gauge field, , on the initial surface. Since we are numerically tracking the values of the full fourpotential, , we must check to make sure the Lorenz gauge condition, , is obeyed in configuration space as required by our equations of motion. The definition of the polarizations, Eq. 6.2, requires that () on the initial slice. Therefore any choice of (with the choice ) obeys the gauge condition. We are, of course, neglecting any effect that the initial conditions of have on on the initial surface. However, we begin our simulations during inflation where all modes of interest are subhorizon. After the fields are initialized, there is an amplification of modes of that reaches equilibrium well before the first zerocrossing of the field. To make sure that we don’t depart from the gaugecondition surface, obeying Lorenz gauge, we track the size of
(6.3) 
as in [59]. Fig. 7 shows a plot of averaged over the box and the RMS averaged value, , for a simulation with . We compare this at two resolutions, and . The reader should understand this as a proxy for staying on the gauge surface—when the simulation diverges from satisfying the gauge condition, this measure quickly grows large. At early times, since we do not ‘cutoff’ the initial spectra, the power in the very highest frequency modes contribute to variations in the finitederivatives we use to calculate are artificially important, and contribute to the size of the RMS value, since the field values are so small. Once the modes of the gauge field are populated, at the end of inflation, this measure is increasingly good.
The goals of our simulations are to see if the tachyonic and/or parametric instabilities identified in Section 5 lead to the efficient generation of gauge modes, whether the energy deposited in those gauge modes is enough to preheat the Universe, and whether that final state has any anisotropy in the two polarization states of the gauge field, as a naive interpretation of Fig. 5 would suggest. To parameterize the success of the first and second of these questions, we calculate the total energy in the gauge field (as defined in Eq. 3.20)
(6.4) 
although we express this as a fraction of the total energy density of the Universe, , for various couplings.
As the simulations progress, we see that energy is, generically speaking, quickly and efficiently transferred into the gauge fields. Fig. 8 shows the evolution of the ratio of the energy in the gauge fields to the total energy as a function of time through the simulation for different values of the couplings.
In all cases where the Universe completely preheats and al