Axion inflation in the strong-backreaction regime: decay of the Anber-Sorbo solution

Axion inflation coupled to Abelian gauge fields via a Chern-Simons-like term of the form $\phi F\tilde{F}$ represents an attractive inflationary model with a rich phenomenology, including the production of magnetic fields, black holes, gravitational waves, and the matter-antimatter asymmetry. In this work, we focus on a particular regime of axion inflation, the so-called Anber-Sorbo (AS) solution, in which the energy loss in the gauge-field production provides the dominant source of friction for the inflaton motion. We revisit the AS solution and confirm that it is unstable. Contrary to earlier numerical works that attempted to reach the AS solution starting from a regime of weak backreaction, we perform, for the first time, a numerical evolution starting directly from the regime of strong backreaction. Our analysis strongly suggests that, at least as long as one neglects spatial inhomogeneities in the inflaton field, the AS solution has no basin of attraction, not even a very small one that might have been missed in previous numerical studies. Our analysis employs an arsenal of analytical and numerical techniques, some established and some newly introduced, including (1) linear perturbation theory along the lines of arXiv:2209.08131, (2) the gradient expansion formalism (GEF) developed in arXiv:2109.01651, (3) a new linearized version of the GEF, and (4) the standard mode-by-mode approach in momentum space in combination with input from the GEF. All these methods yield consistent results confirming the instability of the AS solution, which renders the dynamics of axion inflation in the strong-backreaction regime even more interesting than previously believed.


Introduction
Primordial inflation [3][4][5][6] explains some of the puzzles of standard hot big-bang cosmology, and is well supported by observations [7,8].However, the specific particle physics realization of inflation remains unknown, and concrete implementations encounter theoretical challenges.For instance, explaining how the inflaton potential V (ϕ) can maintain flatness against radiative corrections poses a problem.Additionally, in large-field models of inflation, there is the issue of how the inflaton can span a trans-Planckian range.These challenges become even more significant within the swampland program, where bounds on the curvature of scalar potentials and on the excursion of (pseudo)scalar fields [9] have been conjectured.
Axion inflation offers a potential solution to the first problem by offering a mechanism that protects the inflaton potential from large radiative corrections [10].Moreover, the most natural coupling of the axion inflaton to gauge fields provides additional elements that might help address both problems [11].It is important to note that the simplest model of axion inflation, called natural inflation [10], has been ruled out by data [7,8,12].However, certain nonminimal realizations, such as the two-axion model of aligned inflation [13], or models of monodromy inflation with a flattened potential [14] have solutions compatible with current data [15,16], that will be probed by future CMB experiments, such as CMB-S4 [17].In these realizations, the inflaton excursion ∆ϕ is parametrically greater than the axion scale f that characterizes the potentials of the axion(s).Parametrically, this same scale governs the coupling of the inflaton to gauge fields, where F is the gauge field strength, F its dual, and α ϕ is a dimensionless coupling that, depending on the specific UV completion of the model, can be expected to be of order one.
The motion of the inflaton amplifies one gauge field helicity during inflation [11] with a magnitude that is exponentially sensitive to the parameter where the dot denotes the derivative with respect to cosmic time; H is the Hubble rate; ϵ V is a slow-roll parameter, ϵ V ≡ M 2 P (dV /dϕ) 2 / 2 V 2 ; and M P is the reduced Planck scale, related to Newton's constant by M P = (8πG N ) −1/2 .In the most straightforward case, ∆ϕ = O (f ), the axion scale needs to be close to Planckian, and ξ is suppressed by the smallness of ϵ V .On the contrary, models of aligned natural inflation or of monodromy inflation accommodate an axion scale parametrically smaller than M P , so the parameter ξ in (1.2) can be naturally of order one.
For ξ = O (1), the gauge-field amplification induced by its coupling (1.1) to the inflaton can result into a very interesting phenomenology, as the amplified gauge modes can scatter to produce primordial scalar perturbations and gravitational waves (GWs) [18,19].Interestingly, this stochastic GW background (SGWB) is circularly polarized, as the coupling (1.1) amplifies only one of the two gauge-field polarizations, and these modes produce one GW polarization much more significantly than the other one.Unfortunately, the strong limits on the gauge-field amplification enforced by the non-observation of the scalar perturbations that they produce prevent this GW signal from being observable at CMB scales. 1  The situation might be more interesting at smaller scales.As long as its evolution is described by Eq. (1.2), the parameter ξ grows during inflation.The exponential sensitivity of the gauge-field production to ξ then results in a GW signal that is naturally much greater at scales smaller than the CMB scale [27][28][29], so that the produced GWs might be observable by a variety of GW observatories [30].Also in this case, one must investigate whether this potential GW signature can take place without a simultaneous overproduction of scalar perturbations, which would lead to too many primordial black holes (PBHs) [31][32][33][34][35].For reasons that we discuss in the remainder of this Introduction, we still do not have a reliable answer to this question.
When the parameter ξ is sufficiently large, the amplified gauge fields significantly backreact on the background dynamics.The backreaction occurs via an α ϕ ⟨E • B⟩ /f term in the equation of motion (EOM) of the inflaton and via the gauge-field energy density ∝ ⟨E • E⟩ + ⟨B • B⟩ in the Friedmann equation for the Hubble rate. 2 The former effect is typically more relevant, as it is easier to impact the dynamics of the inflaton, which is slowroll suppressed, than that of the scale factor.Anber and Sorbo (AS) investigated this model for the case where the background evolution is always in a regime of strong backreaction, in which the dissipation due to the gauge-field production provides the dominant source of friction (over the standard 3H φ Hubble friction term) for the inflaton motion [11].Ref. [28] considered instead the case where the evolution is in a regime of weak backreaction (in which the gauge-field amplification negligibly impacts the background dynamics) at the time when the CMB modes exit the horizon, followed by a smooth transition to the regime of strong backreaction, causing the generation of a visible GW signal at smaller scales.In these and in several successive works, the strong-backreaction AS regime was studied under the assumption that the inflaton speed, and hence the parameter ξ (t), follow a slow and monotonic steady-state evolution, characterized by the friction due to gauge-field amplification perfectly balancing the gradient force from the inflaton potential at all times.This is the typical behavior that is expected in a realization of warm inflation [36].
In the past few years, this system has been studied with different numerical schemes of increasing precision and sophistication.Refs.[34,[37][38][39][40] numerically integrated the evolution for the case of a homogeneous inflaton coupled to a large set of gauge-field modes A λ (τ, k).The authors of Ref. [41] adopted a recursive integration approach to study the same system of equations.They initially integrated the equations for the gauge modes using an "external" inflaton and scale factor evolution, where the backreaction of the gauge modes is neglected.Next, they employed these gauge-mode solutions as "external functions" for the backreaction in the evolution equations for the inflaton and scale factor.In this way, they obtained improved solutions for these two quantities from which they could then obtain improved solutions for the gauge modes.By iteratively updating this procedure, they achieved convergence towards a consistent solution encompassing all these quantities.Refs.[2,42,43] took a different approach [44] by considering a set of equations for the two-point correlators of the "electric" and "magnetic" combinations.These equations involve other two-point correlators that include spatial derivatives, such as ⟨E • rot n B⟩.By constructing a hierarchy of equations involving correlators with an increasing number of spatial derivatives, one can numerically solve them after truncating the hierarchy at a certain order.This allows for the computation of the correlators and provides a systematic method for studying the dynamics of the system.This gradient expansion formalism (GEF) is extended to a computation of linear perturbations in this work.
While all these studies assumed a homogeneous inflaton, the equations for a gauge field coupled to a spacetime-dependent inflaton were studied in Refs.[45,46] on the lattice. 3The lattice simulations conducted in Ref. [45] yielded results for the inflaton power spectrum and bispectrum in the weak-backreaction regime, which exhibited excellent agreement with the analytical computations of Refs.[18,55].However, the findings from the numerical studies [2,34,[37][38][39][40][41][42][43][44][45][46] of the strong-backreaction regime contradicted the analytical expectations.It was discovered that the evolution of the inflaton does not occur in a steady-state regime; instead, the parameter ξ (t) undergoes large oscillations, with a period of approximately ∼ 5 e-folds, around the steady-state analytical solution.Ref. [41] showed that these oscillations can be attributed to a memory effect.Specifically, a gaugefield mode begins to undergo amplification when its reduced wavelength is approximately a factor of 2ξ smaller than the Hubble horizon, λ ≡ λ/ (2π) ≃ H −1 / (2ξ).The amplification ceases shortly after the mode crosses the horizon at λ = H −1 , and the energy of the mode is subsequently diluted by the expansion of the Universe.This introduces a sensitivity of the backreaction at a given time t to the evolution of the system during the previous few e-folds.This sensitivity gives rise to an oscillatory behavior in the derivative of the inflaton field, φ(t).These oscillations aim to "adjust" the cumulative effect of the gauge fields amplified during the preceding few e-folds, aligning it with the slope of the potential at the specific moment in question. 4 This interpretation was confirmed by the analytical study in Ref. [1], which solved the linearized set of equations for the homogeneous inflaton perturbation, δϕ, and gauge-mode perturbations about an AS solution with constant H and φ.The equation for the gaugemode perturbations can be formally solved in terms of a Green function (constructed from the gauge modes of the unperturbed AS solution) and the inflaton perturbation.These formal solutions were then substituted back into the equation for δϕ, resulting in an integrodifferential equation where the memory effects are encoded in the kernel of the integral.Through suitable simplifications, Ref. [1] reduced this equation to an algebraic equation, and the roots of this equation were then determined numerically.In addition to providing a relationship between the growth and period of the oscillations and the model parameters, this analysis differs from the previous ones, as it addresses the instability of the AS solution assuming it as an initial condition rather than considering an evolution that started in the weak-backreaction regime and that was expected to evolve into the AS solution.This explored the stability of the AS solution itself, excluding the possibility that it may only have a small basin of attraction that was not reached by the existing studies.
The findings of Ref. [1] confirmed the instability previously observed in numerical studies.These analytical results were, however, obtained with considerable simplifications (particularly, with respect to the form of the Green function), and it is therefore important to confirm them with a numerical study that, contrary to the existing ones, assumes the 3 While this discussion is focused on the inflationary evolution, lattice simulations of reheating at the end of inflation can be found in Refs.[47][48][49][50][51][52][53][54]. 4 To avoid this problem, and generate a stable steady-state dissipative regime, Ref. [56] recently provided a construction, based on scalar field interactions, that can generate an "instantaneous" sensitivity of the backreaction to the particle production.
AS solution as initial condition.This is the goal of the present work.We do this with two independent methods, both based on the GEF.Firstly, we consider the linearized system of perturbations considered in Ref. [1] and rewrite them in GEF language.We do not apply any analytical approximation to these equations, but we rather solve them numerically, using the same general ansatz for δϕ adopted in Ref. [1].Our solutions provide more accurate results compared to Ref. [1], albeit at the cost of a less direct connection with the model parameters.These improved results confirm the presence of unstable modes observed in Ref. [1], yielding more precise values for the growth and period of the oscillations.Secondly, we employ the GEF equations to conduct a full numerical study, starting from the AS solution as initial condition.This approach allows us to consistently incorporate variations in the Hubble rate and the inflaton speed, which occur in concrete inflationary models but were neglected in the analytical computations of Ref. [1] and our first method.Importantly, this represents the first numerical evidence of the instability of the AS solution in the strong-backreaction regime, assuming its validity at the beginning of the evolution.The work is organized as follows.In Sec. 2, we review the basics of the axion inflation model that we are going to be interested in, presenting all the equations and tools that are generally used to study the dynamics of the inflaton and gauge fields in this model.We also precisely define the AS solution and discuss in which part of parameter space we expect this solution to be realized.In Sec. 3, we consider the simplified case of a constant Hubble rate and a constant inflaton dragging force, such that the AS solution corresponds to a constant inflaton velocity.After discussing two methods that allow us to study the stability of the AS solution in the linear perturbation regime (the method in Ref. [1] as well as the linearized GEF), we present the numerical results showing the spectrum of Lyapunov exponents of the linear system as well as the survival time of the AS solution and the latetime behavior of the exact solution in the nonlinear regime.In Sec. 4, we consider the case of a realistic inflationary model where the Hubble rate is now consistently determined by the Friedmann equation and analyze the stability of the AS solution in this case.Sec. 5 is devoted to our conclusions.In Appendix A, we list some auxiliary formulas for the bilinear functions in the GEF.In Appendix B, we explain how we impose the initial conditions for the GEF system for the purposes of the analysis in Sec. 3, and we study the dependence of the survival time of the AS solution on the axion inflation model parameters and on the initial conditions.Finally, in Appendix C, we give more details on a novel self-correction procedure that we employ in our GEF computations and which allows us to extend the applicability of the GEF to later times.Throughout the work, we use natural units and set ℏ = c = 1; then, the reduced Planck mass reads M P = 2.43 × 10 18 GeV.We assume that the Universe is described by a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric with line element (in terms of cosmic time t and conformal time τ ) (1.3) 2 Model and tools

Axion inflation
In the axion inflation model, the pseudoscalar inflaton field ϕ couples to a U (1) gauge field A µ via a Chern-Simons-like term.The corresponding action has the following form: where V (ϕ) is the inflaton potential; is the associated dual tensor with totally antisymmetric Levi-Civita symbol, ε µναβ with ε 0123 = 1; and β ≡ α ϕ M P /f is a dimensionless axion-vector coupling constant.From Eq. (2.1), we can compute the inflaton and gauge-field EOMs, where the latter equation is supplemented by the Bianchi identity for the dual tensor, In temporal gauge, the gauge field is written as A µ = (0, A).Then, the three-vectors of the physical electric field E and magnetic field B in the comoving frame are defined as Correspondingly, the components of the gauge-field tensor and its dual are expressed in terms of the components of E and B in the following way: where ε ijk is the totally antisymmetric Levi-Civita symbol in three spatial dimensions.The energy-momentum tensor following from the action (2.1) reads which, assuming a spatially homogeneous inflaton field, yields an energy density where the angle brackets around T 0 0 and E 2 +B 2 denote the expectation value during inflation.The energy density determines the Hubble expansion rate H through the Friedmann equation, (2.9) Finally, we rewrite Eqs.(2.2)-(2.4) in three-vector form, div E = 0, div B = 0 . (2.13) The system of equations (2.9)-(2.13) is a complete set of equations governing the joint evolution of the homogeneous inflaton field, scale factor, and gauge fields in position space during axion inflation.
Let us now switch to momentum space and consider the quantized gauge field with mode functions A λ (t, k), polarization three-vectors ϵ λ (k), and annihilation (creation) operators âk,λ (â † k,λ ) for electromagnetic modes with momentum k and circular polarization λ = ±, and k = |k|.The polarization vectors satisfy the relations Note that the first relation is equivalent to choosing Coulomb gauge in addition to temporal gauge, which we can impose since div E = 0.The creation and annihilation operators obey the canonical commutation relations (2.16) For the vector potential in Eq. (2.14), the Maxwell equations (2.12)-(2.13)are identically satisfied, while Eq.(2.11) leads to an EOM for the mode functions, which takes a slightly simpler form when written in conformal time τ = t dt ′ /a(t ′ ), Deep inside the horizon, kτ → −∞, the first term inside the square brackets dominates, and Eq.(2.18) takes the form of a simple harmonic-oscillator equation.This allows us to impose the Bunch-Davies boundary condition [57] in the asymptotic past, which amounts to selecting the flat-space positive-frequency solutions for modes deep inside the horizon, Next, we use Eq.(2.14) to compute the vacuum expectation value of E • B on the right-hand side of Eq. (2.10).Then, in conformal time, the Klein-Gordon equation reads Finally, we express the gauge-field energy density in terms of the mode functions.Then, the Friedmann equation for the Hubble rate H = (da/dτ )/a 2 takes the form

Gradient expansion formalism
An alternative way to treat axion inflation in position space is the gradient expansion formalism (GEF) [2].Let us introduce the following set of bilinear gauge-field functions: and recast the Maxwell equations (2.11) and (2.12) as an infinite tower of equations, (2.28) Here, the source terms on the right-hand side of the equations denote boundary terms that account for the fact that the number of physically relevant modes changes in time during inflation, as a consequence of the fact that the cutoff scale in Eq. (2.22) changes with time.These boundary terms were derived in Ref. [2] and are given by the expressions where the auxiliary functions E λ , G λ and B λ are given in terms of Whittaker functions, and where the gauge-field production parameter ξ was already introduced in Eq. (1.2), The quantity k h (t) in Eqs.(2.29)-(2.31) is the wavenumber of the highest-momentum mode that has ever become tachyonically unstable during the entire evolution of the system leading up to the moment of time t, .36)Note that this definition slightly deviates from the momentum scale k h defined in Eq. (2.22); in fact, k h (t) in Eq. (2.36) represents an improved version of k h in the sense that it accounts for the fact that the evolution of the right hand side of Eq. (2.22) is not monotonic when the inflaton velocity oscillates.Employing k h in Eq. (2.22) as the upper integration boundary, therefore, leads to situations where modes that already experienced the tachyonic instability, and which should thus be regarded as physically excited, fall into the region above the UV cutoff and are hence removed from the momentum integrals.With the improved definition, a mode is accounted for at all times after it has become unstable for the first time.In the following, we will exclusively work with the definition in Eq. (2.36).Finally, in order to solve the system of equations numerically, we need to truncate it at some finite order n cut .This can be done by expressing the quantities of order (n cut + 1) through expressions of lower order.One of the simplest ways to truncate the system was discussed in Ref. [2] and it is based on the following relation: and analogously for G (ncut+1) and B (ncut+1) .This truncation condition can be justified by the following consideration.For sufficiently large order n, the spectrum of bilinear quantities E (n) , G (n) , and B (n) is always blue and, therefore, the dominant contribution stems from modes with momenta k ≃ k h .Then, the mean value theorem for the integral over momentum leads to the condition (2.37) [2].Knowing the behavior of the spectrum near the cutoff momentum k h , one could estimate the error of the truncation introduced by the truncation condition.In practice, however, the truncation order n cut is chosen in such a way that increasing n cut further does not have an impact on the solution.

Anber-Sorbo solution
If the gauge field is absent (or, at least, sufficiently weak such that its contributions to the Friedmann and Klein-Gordon equations can be neglected) and the potential V (ϕ) is sufficiently flat, the inflaton ϕ follows the slow-roll attractor solution where the Hubble friction term is almost exactly compensated by the gradient force from the inflaton potential, This describes the usual case of slow-roll inflation, which is a true attractor solution: any initial deviation quickly tends to zero, leading the system into the slow-roll regime.
The idea of Anber and Sorbo in Ref. [11] was to realize the inflationary stage even with a steep potential V (ϕ) (for which the standard slow-roll regime is impossible) due to the backreaction from the produced gauge fields.In this case, the potential-gradient force is compensated by the gauge-field friction term on the right-hand side of the Klein-Gordon equation, However, this solution is now known to be unstable, as discussed in detail in the Introduction.A first analytical computation demonstrating the instability of the AS solution was recently presented by two of us (M.P. and L. S.) in Ref. [1].In the present paper, we shall substantiate this computation by a comprehensive numerical analysis that will allow us to achieve two results: (i) precisely determine the growth rate of instabilities around the AS solution and hence (ii) measure the survival time of the AS solution, which we define as the time when the relative deviation from the AS solution becomes of O (1).First of all, let us give a precise definition of the AS solution that we will use throughout the paper.Let us assume for a moment that the AS solution is indeed an attractor solution, as initially proposed in Ref. [11] and subsequently assumed as a working hypothesis in many papers in the literature.Then, under this assumption, we expect the system to slowly evolve in time because of the strong friction coming from the gauge-field backreaction.In this regime, it is natural to assume that the value of the gauge field at a given moment of time is determined by the inflaton velocity at the same moment of time, i.e., there is no retardation between the inflaton evolution and the gauge-field response.In this case, the Klein-Gordon equation (2.10) represents a closed equation for the inflaton field, where we now emphasize that ⟨E • B⟩ is a functional of the inflaton velocity φ and where the Hubble parameter follows from the Friedmann equation in the usual way, In order to find the explicit form of the functional dependence of the energy densities and the Chern-Pontryagin density on φ, we assume again that the inflaton is rolling slowly and the Universe is expanding quasi-exponentially.Therefore, on the timescale of a few e-folds, one can assume that the Hubble parameter H and the inflaton velocity φ are approximately constant.This significantly simplifies the mode equation (2.18) in conformal time, where we used the relation a(τ ) = −1/(Hτ ) for de Sitter space and where the parameter ξ introduced in Eq. (2.35) is now constant (for definiteness, let us assume ξ > 0).In this case, Eq. ( 2.42) has an exact solution in terms of Whittaker functions.The Bunch-Davies boundary condition (2.19) allows to extract a unique solution of the form where W κ,µ (z) is the Whittaker W function.Eq. (2.42) implies that negative-helicity modes A − are not enhanced, because the expression in square brackets is always positive.However, positive-helicity modes A + experience a tachyonic instability when k < 2ξ aH, which leads to their exponential amplification.This can be directly seen from the following approximate expression for A + in terms of elementary functions, which is valid for x ≪ 2ξ [11]: Now, substituting Eq. (2.43) into the expressions for the Chern-Pontryagin density and energy densities, we find where the functions e 0 , g 0 , and b 0 are given in the form of integrals of the Whittaker functions (A.5)-(A.7) in Appendix A. Using Eq. (2.44), we can approximately write (2.46) Therefore, we now define the AS solution as a solution of the system of equations Note that for a generic inflationary model, the solution of Eqs.(2.47)-(2.48) is not a solution of the full system of equations (including Maxwell's equations for the gauge field), since in the full system the assumptions underlying Eqs.(2.47)-(2.48),namely H = const and ξ = const, are not necessarily satisfied.Therefore, we will often refer to this solution as the "enforced" AS (EAS) solution, which is characterized by the fact that we insist on (or "enforce") the specific functional dependence encoded in the functions e 0 , g 0 , and b 0 .We study the stability of the EAS solution in a specific inflationary model in section 4. On the other hand, in the particular case of constant background quantities, considered in detail in Sec. 3, the AS solution is an exact solution of the full system of equations, which can be realized by choosing the right initial conditions. 5Nonetheless, a main result of our analysis will be that, even though the AS solution is an exact solution of Eq. (2.47) for constant H and ξ, it turns out to be unstable against arbitrarily small perturbations and therefore only has a finite survival time, as we will demonstrate in detail.In Table 1, we summarize our notations for the different types of solutions considered in the next sections as well as the methods used to obtain them.

Parameter choice
Let us now specify the model parameter values that we are going to be interested in, i.e., the region in parameter space corresponding to the strong-backreaction regime.To do so, we assume that the system is initialized in phase space either in the AS solution or at least sufficiently close to it, so that we can employ the equations of motion (2.47)-(2.48).
Specifically, we shall impose two conditions.On the one hand, the backreaction in the Klein-Gordon equation must be strong, meaning that the additional friction from the gauge field dominates over the Hubble friction term and hence governs the evolution of the inflaton field.On the other hand, the contribution of the produced gauge fields to the total energy density of the Universe needs to remain small compared to that of the inflaton, so that we can still realize a stage of accelerated (inflationary) expansion with the effective equation of state parameter w = p/ρ < −1/3.In order to give a quantitative meaning to these two conditions, we introduce two parameters, δ KG and δ F , which measure the strength of the backreaction in the Klein-Gordon and Friedmann equations, respectively, ) Then, the conditions determining the desired parameter range can be formulated as These two conditions are independent and satisfied across an extended volume in the threedimensional parameter space spanned by β, ξ, and H.Still, it will be helpful to define a benchmark in the sense of an "optimal parameter choice" determined by the condition This choice gives the central section of the relevant parameter range where the backreaction is strong in the Klein-Gordon equation and, at the same time, small in the Friedmann equation.We expect that even away from the optimal parameter choice, as long as the conditions (2.51) are satisfied, we will obtain qualitatively the same results.The condition (2.52) allows to eliminate one of the parameters, e.g., the Hubble rate, in terms of the two other parameters, (2.53) To obtain an intuition for the analytical dependence on the parameter ξ, we use the approximate expressions in Eq. (2.46), which allow us to write the more explicit expression Moreover, in the range of ξ values 5 ≲ ξ ≲ 10, which is the most interesting for the present study, we find a simple empirical relation, i.e., a fit formula, which reproduces the exact result up to an error of a few percent.This is good enough for us; the condition in Eq. (2.52) is not an exact requirement, anyway.It merely serves the purpose of providing us with guidance as to where in parameter space we can expect the strong-backreaction regime of axion inflation to be realized.Unless specified otherwise, we will therefore use the relation (2.55) in all computations in the remainder of this work.

Constant background quantities
Let us start our discussion of the AS solution considering the simple case of constant Hubble rate and constant inflaton potential gradient, H = const and V ′ (ϕ) = const.In this case, the Universe expands exponentially (de Sitter spacetime), and we effectively disregard the Friedmann equation. 6At the same time, the Klein-Gordon equation admits a solution with constant inflaton velocity.Indeed, setting φ = 0 in Eq. (2.47) and expressing everything in terms of the ξ parameter, which is also constant, we obtain the following equation: For given values of the axion-vector coupling β, Hubble rate H, and potential gradient V ′ , this equation can be solved for the associated constant value of ξ.Thus, in this simple case, the AS solution turns out to be ξ (t) = ξ 0 = const; see Table 1.It is important to note that this is not only the solution of the approximate Eq. (2.47), but also a particular solution of the full system of equations, including the EOM for the gauge-field modes.Since ξ = const, the gauge-field mode functions have the form (2.43) and, therefore, the simple relations in Eq. (2.45) are exact.In particular, this means that a system prepared exactly in this state will remain in it forever.It is interesting, however, to study the stability of this solution and consider the evolution of small perturbations around it.In what follows, we shall denote all quantities in the AS solution by a bar, e.g., ξ.

Linear perturbation theory
In order to study the stability of the AS solution, let us construct a linear perturbation theory for deviations from this solution following the same strategy as in Ref. [1].We write the perturbed quantities as The background quantities evolve according to equations similar to Eqs. (2.18) and (2.10), Subtracting Eqs.(3.3), (3.4) from Eqs. (2.18), (2.10), respectively, and keeping only perturbation terms up to linear order, we obtain the system of EOMs for the perturbations, The solution of Eq. (3.5) can be formally expressed as, where the Green function G λ,k (τ, τ ′ ) is a solution of the differential equation The differential operator acting on the Green function is exactly the same as the differential operator in Eq. (3.3).Ā(τ, k) is hence the solution of the corresponding homogeneous equation.A second linearly independent solution is Ā * (τ, k), since the mode equation has real coefficients.This allows us to construct the retarded Green function as follows: where we used the fact that the mode functions are normalized in such way that their Wronskian equals ∂ Āλ (τ, k) ∂τ Substituting Eq. (3.7) into Eq.(3.6), we get the source term on the right-hand side, Taking into account that Im[. ..] = 0 when τ = τ ′ in this expression, we evaluate the derivative with respect to τ and finally obtain the EOM for the scalar-field perturbation, In Eq. (3.12), we kept the most general form of the mode functions.Next, in order to simplify Eq. (3.12), we assume that the gauge-field mode functions can be represented as where W λ (x) is an arbitrary function of x and λ for the time being.Both the exact solution in Eq. (2.43) and the approximate one in Eq. (2.44) can be represented in this way.Moreover, following Ref.[1], we look for power-law solutions of the EOM (3.12).We therefore choose the following ansatz for the scalar-field perturbation: where C and ζ are constant.Once this ansatz is inserted in Eq. (3.12), and the allowed ζ are found, we exploit the fact that, for any solution ζ, also its complex conjugate ζ * is a solution, to obtain a real inflaton perturbation, where W ′ λ (x) = dW λ (x)/dx.Up to now, we were able to work with a general function W λ (x).However, in order to determine the allowed power-law exponents ζ n , one needs to specify the function W λ (x) and solve Eq. (3.16) numerically.In the case of constant H and ξ, which we consider here, the function W λ (x) can be extracted from Eq. (2.43), (3.17) This equation may be simplified following Ref.[1].Firstly, we take into account only the enhanced λ = + gauge polarization.Secondly, we replace the upper integration limit in the second integral by 2 ξ, which is motivated by the fact that for x ′ > 2 ξ the integrand is no longer enhanced, but it actually becomes a rapidly oscillating function that integrates to a negligible amount.Thirdly, if we use the approximate form of the mode function A + in Eq. (2.44), we obtain an approximate equation for ζ of the form which agrees with Eq. (3.14) in Ref. [1], expressed in the notation of the present paper.

Linearized gradient expansion formalism
The gradient expansion formalism (GEF) introduced in Sec.2.2 allows us to find an exact numerical solution for the system of coupled inflaton and gauge-field EOMs in the strongbackreaction regime.For constant background quantities H and V ′ , the system of equations can be further simplified.First of all, since the inflaton field ϕ itself does not appear in the system and only its derivative φ is involved, it is convenient to use the parameter ξ in Eq. (2.35) as a new field variable.This renders the Klein-Gordon equation a first-order ordinary differential equation, like all the equations of the GEF.Moreover, since H = const, it is more convenient to work with the number of e-folds, N = ln a = Ht as a time variable, instead of physical time t.This leads us to: where v and b are dimensionless parameters accounting for the dragging force caused by the potential gradient and the axion-vector coupling, respectively.This equation needs to be supplemented by the tower of equations that govern the gauge-field evolution: ) ) where E λ , G λ , and B λ are functions of the production parameter ξ given in Eqs.(2.32)-(2.34).In Eqs.(3.20)- (3.22), moreover, we introduced the dimensionless bilinear functions the dimensionless momentum of the horizon-crossing mode and its derivative where θ(x) is the Heaviside unit step function.
Although this system can be directly employed to study the true solution in the strongbackreaction regime, it is also instructive to linearize it for small deviations from the AS solution.Denoting all quantities in the AS solution by bars and small deviations by δ's, we obtain the following system of equations, which define what we shall refer to as the linearized gradient expansion formalism (LGEF): ) ) ) where ) ) This system is also infinite in principle, and, in order to use it in practice, one has to truncate it at some order n cut .The simplest way to do so is to assume that, for all orders larger than n cut , the bilinear functions exactly coincide with the background values in the AS solution, i.e., The advantage of the LGEF compared to other methods is that it leads to a system of linear ODEs with constant coefficients.Its solution can easily be found by methods of linear algebra.In particular, the ansatz δξ ∝ e ζN (the same as in the previous subsection), and similarly for all perturbations of bilinear functions, reformulates the problem from studying the evolution in time to just finding eigenvalues of the matrix of the linear system, i.e., to a purely algebraic task.In practice, this turns out to be the simplest approach.

Results and discussion
We shall now discuss our numerical results obtained for the case of constant Hubble rate H and constant potential gradient V ′ .As noted above, the AS solution is a true solution of the Maxwell and Klein-Gordon equations.This means that, if we perfectly fine-tune the initial conditions for ξ and all gauge-field bilinear quantities to be in the AS solution, the system will remain in this solution indefinitely.Perfectly fine-tuned initial conditions are, however, of little interest.In fact, they are even impossible to achieve in any numerical study with finite numerical precision.In what follows, we will therefore slightly detune the initial conditions and study the dynamical evolution of the system away from the AS solution.As we will find that the AS solution is unstable, we sometimes denote this as the "decay of the AS solution", which also features prominently in the title of this paper.
In this subsection, we apply all three approaches discussed above in order to study the stability of the AS solution with respect to small perturbations.The linear perturbation theory and the LGEF allow us to determine the spectrum of Lyapunov exponents ζ and, thus, to capture all possible scenarios for the evolution of the system at once, however, only in the regime of small perturbations.On the other hand, the full GEF allows us to get an exact numerical solution of the system for some specified initial conditions that is valid also for large deviations from the AS solution.Therefore, these methods are complementary to each other and allow us to study the same system from different angles.

Lyapunov exponents
It is instructive to first work out the spectrum of Lyapunov exponents ζ for our physical system.For definiteness, we perform the numerical analysis in the region of parameter space ξ 0 ∈ [5, 9], β ∈ [10 1.5 , 10 3.5 ], and H determined by Eq. (2.55).For the benchmark point with β = 10 2.5 and ξ 0 = 7 (which is the central point of the specified parameter range), this spectrum is shown in Fig. 1 in the form of a sequence of red dots in the complex plane for ζ.These points have been found by using the LGEF truncated at n cut = 70.We want to compare the results of this numerical integration with the solutions of Eq. (3.16) for ζ.This equation, which follows from the linear perturbation theory, contains integrals of highly oscillatory special functions, and it turns out that finding its solutions is computationally very costly.We have checked that, for the root with the greatest real part, ζ 1 , the numerical solution of Eq. (3.16) is in perfect agreement with the LGEF result presented in the figure.This root is of great importance since, having the greatest real part, it is the one that controls the growth rate of the instability at late times; see Eq. (3.14).A much quicker comparison of the eigenvalues obtained from the LGEF system with the analytical computation can be done in terms of the approximate equation (3.18) (i.e., Eq. (3.14) in Ref. [1]), which is easier to solve.We show the roots of this equation by green empty circles in Fig. 1.The background color of this figure is the density plot of the absolute value of R − 1, the difference between the function on the lefthand side of Eq. (3.18) and unity.As evident from Fig. 1, the approximate equation (3.18) allows us to obtain the spectrum of the Lyapunov exponents ζ with good accuracy.In particular, it also shows that among the ζ values there is at least one with positive real part meaning that the AS solution is unstable; this confirms the findings of Ref. [1].Fig. 1 also shows that there is a one-to-one correspondence between the solutions obtained from the LGEF method and the solutions that we were able to derive in the context of linear perturbation theory and starting from the ansatz in Eq. (3.14).This observation serves as another (a posteriori ) justification for the ansatz in Eq. (3.14) and confirms that we did not overlook any solutions in our linear-perturbation-theory analysis in Sec.3.1.
It is important to note that the Lyapunov exponents are complex numbers.As is clear from Eq. (3.15), a generic complex ζ with positive real and non-vanishing imaginary parts indicates that the deviation from the AS solution shows an oscillatory behavior with a growing amplitude, where Re(ζ) determines the growth rate while Im(ζ) is the angular frequency of the oscillations.Let us now focus on the eigenvalue with the greatest real part, ζ 1 , which corresponds to the fastest-growing mode.This mode also has a nonvanishing imaginary part.In Figs.2(a) and (b), we show the real and imaginary parts of this root as functions of ξ 0 .For each ξ 0 , we actually present a band of values assumed by this root.The different values inside the band are obtained for different values of the axion-vector coupling in the range β ∈ [10 1.5 , 10 3.5 ].The blue bands shown in the two panels are obtained using the LGEF system while the green bands follow from Eq. (3.18).We find that the exact numerical results in blue are in excellent qualitative agreement and good quantitative agreement with the approximate analytical results in green.This observation serves as a validation and refinement of the results presented in Ref. [1] and is one of the main results of the present work.The AS solution is unstable and the fastest-growing perturbation mode is characterized by the growth rate Re(ζ 1 ) and oscillation frequency Im(ζ 1 ) in Fig. 2.
Re(ζ 1 ) is a monotonically increasing function of ξ 0 and a decreasing function of β at fixed ξ 0 .The imaginary part exhibits the opposite behavior.For comparison, in Fig. 2(b), we also show the estimate for the oscillation frequency found in Ref. [41], This expression follows from the fact that the response of the gauge-field Pontryagin density ⟨E • B⟩ to changes in the inflaton velocity is retarded by approximately ∆N ξ ≃ ln(ξ 2 0 /2).The main reason for this delay is that the growth rate of modes that cross the horizon and undergo the tachyonic instability at a certain moment of time is determined by the instantaneous value of ξ; however, these modes are still not dominating the spectrum of ⟨E • B⟩ at this moment of time and will only do so ∆N ξ e-folds later.We point out the good qualitative agreement between the estimate (3.34) and our numerical results.At the quantitative level, the two values agree up to a factor of roughly 1.5.

Decay of the AS solution
Next, let us discuss what we shall refer to as the survival time of the AS solution.This quantity can be naturally defined as the moment of time when the relative deviation of the full numerical solution from the constant AS solution becomes of order unity.To be specific, in this work, we define the survival time as the first moment of time (or the number of e-folds from the beginning) when the production parameter ξ deviates from the initial value in the AS solution, ξ 0 , by more than half an order of magnitude, where the N i denote the moments (in terms of the number of e-folds N ) when the condition on the right-hand side is satisfied, and where we initialize the system at N = 0. We emphasize that, in contrast to the Lyapunov exponents ζ n , which are intrinsic and characteristic properties of the physical system, the survival time N AS depends on the way in which one imposes initial conditions.For the sake of definiteness, we assume that all gauge-field bilinear quantities are in the AS solution corresponding to a certain value of the production parameter ξ 0 , while the initial ξ value is detuned to ξ 0 + δξ 0 by some small amount δξ 0 .Since, for the constant background case considered in this section, the AS solution is an unstable equilibrium solution, the survival number of e-folds depends on the detuning parameter δξ 0 (it is infinite for δξ 0 = 0) as well as on the instability rate controlled by the Lyapunov coefficients studied above.As mentioned above, in the linear regime the quantity ξ (N ) − ξ 0 is the linear superposition of a series of eigenmodes, each characterized by a Lyapunov exponent ζ n .Let us denote this series as where the real coefficients r n and phases φ n depend on how the initial δξ 0 projects on each eigenmode, and where the initial number of e-folds has been set to N = 0. Assuming a non-zero overlapping with the fastest growing mode (namely, the mode whose Lyapunov coefficient, denoted above as ζ 1 , has the greatest real part), and ignoring the initial phase φ 1 , we then have leading to the survival number of e-folds where we ignore the order one term proportional to ln 10 1/2 r 1 in our discussion.We verify the dependence of N AS on the initial conditions numerically in Appendix B and show that it is in a good accordance with the estimate (3.38); see Fig. 9.This provides a nontrivial check of the validity of our numerical scheme.Moreover, we study the dependence of the survival time on the parameters of the axion inflation model, which is also presented in Appendix B. In this appendix, we also provide more details on how we choose the GEF initial conditions.In Fig. 3, we instead present one specific example of the departure from the AS solution.The evolution shown is characterized by β = 10 2.5 , ξ 0 = 7, and an initial fine-tuning of δξ 0 /ξ 0 = 10 −6 .The green solid line in the figure shows the evolution of δξ/ξ 0 as computed using the GEF system.In panel (a) of Fig. 3 The first two eigenvalues, ζ 1 and ζ 2 , have positive real parts, namely, they correspond to unstable departures from AS.The remaining eigenvalues have negative real parts.They correspond to stable departures of decreasing amplitudes.As visible from panel (a), the most unstable mode ζ 1 is by itself able to account for the departure of ξ from its initial value throughout the linearized stage.The inclusion of the other unstable modes, and, particularly, of the stable mode ζ 5 , well reproduces also the initial phase visible in the figure in which δξ decreases. 7We stress that the decreasing stage visible in the figure is due to the fact that our choice of initial conditions is mostly mapped to this stable mode of eigenvalue ζ 5 , and it is not indicating that the AS solution is initially stable.
Finally, let us consider the late-time behavior of the solution, in the regime of large deviations from the AS solution, i.e., for N > N AS .This is the region where neither linear perturbation theory nor the LGEF are applicable.Therefore, we can only use the GEF in order to find the solution of the equations of motion.For the benchmark scenario with ξ 0 = 7 and β = 10 2.5 , the GEF solution is shown by the blue solid lines in Fig. 4. Panel (a) shows the time evolution of the ξ parameter, while panel (b) illustrates the behavior of the produced gauge-field energy densities ρ E = ⟨E 2 ⟩/2, ρ B = ⟨B 2 ⟩/2, and the Chern-Pontryagin density ρ EB = |⟨E • B⟩|/2.The red dashed lines show the corresponding AS solution.As we see from the plots, the time behavior of all quantities becomes almost perfectly periodic showing a sequence of highly oscillatory phases.In the simple case of constant background quantities, these oscillations will last indefinitely.
It is important to note that, for such a complicated and nonmonotonic behavior of the 7 The eigenvalue ζ5 corresponds to the point in the upper half of Fig. 1 that is well separated from the regular sequence of roots in the lower part of the plot.Its imaginary part is significantly greater than that of the first eigenvalues of the lower sequence.Correspondingly, δξ oscillates much faster during the initial decreasing stage than in the following unstable phase.In passing, we also mention that a separated root, such as the ζ5-mode in the present benchmark scenario, is not always present in the spectrum.For instance, fixing β = 10 2.5 as the value considered in the figure, this separated root also exists for ξ0 = 5, but not for ξ0 = 9.The complicated form of the equations that we are solving, even in the simpler approximate form (3.18), does not allow us to determine a priori whether this separate root is present or not.ξ parameter as shown in Fig. 4, the cutoff momentum k h , for which Eq. (2.36) gives is not growing at all times; instead, there is a sequence of plateaus in k h (N ) where it remains constant.During these plateau stages, the truncation condition in Eq. (2.37) may not be accurately satisfied because some of the underlying assumptions [see the paragraph below Eq. (2.37)] are not valid at this time.In particular, the assumption that the spectral densities of E (n) , G (n) , and B (n) are dominated by the mode k h at large n may be violated.These effects can lead to the breakdown of the GEF in a way that we discuss in Appendix C.
In this appendix, we also show one possible solution to this problem, namely, the GEF selfcorrection procedure.

Time-dependent background quantities
In this section we turn to the realistic case in which the Hubble parameter is not a constant, as assumed in the previous section, but it is consistently determined by the Friedmann equation (2.9).For simplicity, we still consider the potential gradient V ′ to be constant, meaning that the inflaton potential is a simple first-order polynomial, Such a potential is still not a realistic choice that could describe the whole stage of inflation over a large range in field space.A linear potential violates, e.g., the constraints on the tensor-to-scalar ratio imposed by CMB observations [7,8], and does not allow for a graceful exit from the inflationary stage.Nevertheless, it serves as a good local approximation for a variety of potentials in restricted regions of field space.In any case, in the following, we shall use the simple potential in Eq. (4.1) primarily for illustrative purposes.
Firstly, let us give a recipe to determine the AS solution in any realistic inflationary model.As discussed in Sec. 2, the AS solution ignores the retardation of the gauge-field response to the changes in the inflaton velocity.This allows to get a closed set of equations (2.47)-(2.48)determining the inflaton evolution.Since in the derivation of those equations we used the expressions for the gauge-field energy densities and the Chern-Pontryagin density for constant ξ and H parameters, the solution of Eqs.(2.47)-(2.48) is not a solution of the full system of equations, which treats the Maxwell equations for the gauge field on the same footing as the Klein-Gordon and Friedmann equations.For this reason, we refer to the solution of Eqs.(2.47)-(2.48)as the "enforced" AS (EAS) solution.
In order to find the numerical EAS solution of Eqs.(2.47)-(2.48), it is more convenient to rewrite the EOMs as a system of two first-order differential equations for the functions ϕ(t) and ξ(t).In order to do so, we first solve Eq. (2.48) with respect to H, Then, the EOM for ϕ follows from the definition of ξ in Eq. (2.35), while the equation for ξ can be derived from Eq. (2.47).Finally, the desired system of equations has the form In order to study the stability of the AS solution in a fully time-dependent background, we have to initialize the system of Eqs.(4.3)-(4.4)and, simultaneously, the full system of Klein-Gordon (2.10) and GEF equations (2.26)-(2.28)using exactly the same initial conditions.We use the following algorithm to impose these initial conditions: (i) Specify the axion-vector coupling β and the desired value of ξ 0 ; propose a Hubble rate H 0 according to Eq. (2.55) that puts the system in the strong-backreaction regime.
(iii) Initialize the EAS system of Eqs.(4.3)-( 4.4) at a slightly smaller value of ξ < ξ 0 and ϕ < 0, to allow the system to cope with any potential inconsistency of the initial conditions.That is, let the system dynamically roll into a smooth and stable solution.
(iv) Determine the moment of time when the ξ parameter in the EAS solution crosses the desired value ξ 0 and use it as the initial moment of time for the GEF.Note that this typically does not happen exactly at the origin of field space.We therefore no longer intend to initialize the system when φ0 = 0 and ϕ 0 = 0.These conditions were used in step (ii) only to obtain some reasonable values for V 0 and V ′ .For given values of V 0 and V ′ , we can now forget about the fact that they were derived assuming φ0 = 0 and ϕ 0 = 0. Instead, we now use the set of consistent values of ϕ, φ, and H that we dynamically reach when the system has rolled into a smooth solution and ξ corresponds to the desired value ξ 0 .8The advantage of this procedure is that it allows us to initialize the GEF with a self-consistent set of input values, including φ and hence the time derivative of ξ.
(v) Use ϕ and φ from the previous step to initialize the Klein-Gordon equation (2.10) for the GEF; compute the initial conditions for the bilinear functions in the GEF by inserting ξ 0 and H into Eqs.(A.5)-(A.7) in Appendix A.
Below, we present the numerical results we obtain using the GEF and compare them to the corresponding results based on the EAS solution.
Figures 5-7 present the results of the evolution for three different realizations of the model.More precisely, in all three figures the same value of the axion-vector coupling β = 10 2.5 and of the Hubble parameter H 0 = 2.7×10 11 GeV are assumed.The three figures differ by a decreasing steepness of the potential V ′ and, consequently, in the correspondingly required value of ξ 0 (the specific values of ξ 0 and of the slow-roll parameter ϵ V,0 assumed in each figure are shown in their titles).Panels (a) of each figure show the evolution of the ξ parameter as a function of number of e-folds N , while panels (b) of each figure show the evolution of the energy densities ρ E , ρ B , and the Chern-Pontryagin density ρ EB of the produced gauge fields.In all figures and panels, the blue solid lines show the exact solution of the full system of equations obtained from the GEF, while the red dashed lines show the corresponding quantities in the EAS solution.Vertical dotted lines of the same color show the moment of time when inflation ends (ä = 0) for the chosen initial conditions, and the gray dashed vertical lines denote as usual the survival time of the AS solution, i.e., the moment when the relative deviation of the blue curve from the red curve for the evolution of ξ reaches the threshold value of 10 −1/2 for the first time.In the following, we discuss a few features that can be read off from these figures.
Firstly, we note that the survival time of the EAS solution cannot be increased to infinity by fine-tuning the initial conditions.This is due to the fact that, once H and φ are allowed to consistently vary, the EAS solution is no longer an exact solution of the full system of equations.This poses an upper limit on the survival time, which we typically found to be around 5 to 7 e-folds.Secondly, the evolution of the exact numerical solution of the system is qualitatively the same as in the simple case of constant H considered in Sec. 3. Indeed, the curve for the ξ parameter initially shows oscillations around the AS solution with a growing amplitude until the deviations becomes eventually of order unity and the growth stops due to nonlinearities in the system.After that, a new phase typically sets in, characterized by quasiperiodic stages of fast oscillations.However, now, due to the continuously changing Hubble parameter, the amplitude and oscillation frequency slowly change in time.Thirdly, the comparison between the different figures shows that the duration of inflation strongly depends on the steepness of the potential.Note that, for the GeV and slow-roll parameter ϵ V,0 = 6.1 (corresponding to ξ 0 = 7.1).The blue solid lines show the true solution of the system found using the GEF, while the red dashed lines correspond to the enforced AS solution.The vertical dashed lines show the moment of time N = N AS when the linear perturbation theory breaks down (relative deviation of ξ from its initial value exceeds 10 −1/2 ).The vertical dotted lines show the end of inflation for the real system (blue) and for the enforced AS solution (red).case of a very steep potential (with initial value of the slow-roll parameter ϵ V,0 ≈ 6.1) shown in Fig. 5, the duration of inflation is just 7 e-folds.The second phase of fast oscillations of ξ does not even start in this case.For a flatter potential with ϵ V,0 ≈ 1.9, shown in Fig. 6, inflation lasts for approximately 13 e-folds, allowing for one stage of fast oscillations.Further flattening the potential with ϵ V,0 ≈ 0.4, as shown in Fig. 7, the duration of inflation is greater than 30 e-folds, and here five periods of fast oscillations appear. 9Note that the fast oscillatory stages in the ξ evolution lead again to the plateaus in the evolution of k h , which complicates the integration of the GEF system.In particular, in order to obtain the numerical results with a controllable accuracy until the end of inflation, we need to apply several self-correction procedures, which we discuss in Appendix C. In Fig. 8, we further elaborate on the comparison of the evolution of the ξ parameter obtained in a realistic inflationary model (the same blue solid line as in Fig. 7) with that obtained under the assumption of constant H (equal to the initial H value in the realistic model; the evolution of ξ in this case is shown by an orange solid line).Both evaluations are shown as a function of the number of e-folds N , with the value N = N AS set when the ξ parameter has increased by half an order of magnitude with respect to the corresponding AS value (which equals ξ 0 = 6.85 in the case of constant background quantities and which is taken from the EAS solution in the case of time-dependent background).The figure clearly demonstrates that the two solutions are very close to each other, not only qualita-tively but even showing a good quantitative agreement.In particular, the amplitude and frequency of the large-amplitude oscillations in the nonlinear-perturbation regime are in a good agreement.We only observe a slow drift of the parameters due to slow-roll corrections for the realistic inflationary model.This result means that the case of constant background quantities considered in great detail in Sec. 3, despite its simplicity, still allows to obtain a good intuition for the time evolution in the realistic case.We attribute this agreement to the fact that, also in the realistic case, H changes only very slowly during inflation.

Conclusions
An axion-like inflaton coupled to a gauge field provides probably one of the simplest and best motivated models where matter production occurs during inflation.Several studies in the last few years have shown that, in the strong-backreaction regime, this model displays a highly nontrivial behavior, significantly different from the steady state AS solution.Such a behavior has been mostly studied via numerical techniques (with the exception of Ref. [1]), and the origin and the fate of this departures are the subject of ongoing research.
In this work, we used the gradient expansion formalism developed in Refs.[2,43,44] to improve our understanding of the system.In the first part of our work, we studied a setup similar to that of Ref. [1], where the velocity of the inflaton is initially close to its AS value, under the assumption of a constant Hubble parameter.The analysis we presented in Sec. 3 leads to results that agree with those of Ref. [1]: among the complex Lyapunov exponents for perturbations around the AS solution, there is always at least one exponent with a positive real part, which results in oscillation of increasing amplitude of the inflaton velocity about the value predicted by the steady state AS solution.While the results of Sec. 3 were obtained for a system identical to that of Ref. [1], the techniques used, and in particular the approximations, are significantly different.For instance, Ref. [1] uses an approximate Green function, which is not required in the solution of Sec. 3, whereas the truncation of the hierarchy of equations in Sec. 3 has no counterpart in the analysis of Ref. [1].For this reason, those results corroborate each other.In the second part of our work, we studied the onset of the instability in the case in which the Hubble parameter is evolving over time.Previous numerical studies in the literature started in the weak-backreaction regime and saw the instability building up as the backreaction becomes strong.One might therefore wonder whether the AS solution might still be stable (even if, possibly, with a very small basin of attraction) if one started directly from the strong-backreaction regime.The analysis of Sec. 3, which assumes Ḣ = 0, does not fully settle this question, since in that case δξ = 0 is a valid (albeit unstable) solution that might in principle be stabilized when Ḣ ̸ = 0.Although it is hard to imagine that a slow-roll variation of H could change this behavior, this remained a logical possibility from the existing literature.Our analysis of Sec. 4 shows that this is not the case.Even if we place our solution on the AS values in the strong-backreaction regime "by hand", the time dependence of the Hubble parameter will destabilize the system, leading in only a few e-folds to the oscillating behavior observed in previous studies of the transition from weak to strong backreaction.
While we believe that our work clarifies, and possibly settles, the questions around the onset of the instability in the axion-vector system, there are still open questions concerning the subsequent evolution and the possible end of the unstable regime.Thus far, these questions have been tackled only with the use of numerical studies.The majority of those studies assumes a spatially uniform inflaton field, which results in a quasi-periodic pattern of oscillations in the inflaton velocity, where deviations from perfectly periodic oscillations are only due to the slow-roll evolution of the inflaton zero mode probing different parts of the potential at different times.So far, only two works have considered the effects of spatial fluctuations of the inflaton field during the inflationary stage (and none has considered the effects of metric perturbations around an FLRW background).These lattice studies, being computationally expensive, covered only a relatively brief time interval.The study of Ref. [45] was able to capture the first oscillation of the inflaton velocity, whose shape agrees with that found in the works with a homogeneous inflaton.More recently, Ref. [46] studied the system for a more extended time range, showing that, around the time of the first oscillation, spatial inhomogeneities in the inflaton field build up very rapidly, and that the subsequent oscillations in the inflaton velocity have a suppressed amplitude.This is a relevant qualitative change with respect to the previous results, which warrants further study.How does this behavior depend on parameters?(Indeed, the pattern of oscillations in the examples shown in [46] changes significantly for very small variations of the axial coupling.)To what extent are the results in Ref. [46] affected by the fact that the simulations are probing only the last ≲ 10 e-folds of inflation?We hope that a (semi-)analytical study extending the formalism presented here might shed more light on these questions; we plan to return to this problem in the future.the survival time N AS ≈ 29 e-folds, according to Fig. 3 in the main text.If we choose other values for the fine-tuning of the initial condition, δξ 0 /ξ 0 , the survival time of the AS solution will change.This dependence is shown in Fig. 9.The blue solid line (green dotted line) corresponds to the initial deviation in the direction of greater (smaller) ξ values.Not surprisingly, the survival time has a clear logarithmic dependence on |δξ 0 | as long as this quantity is in the linear regime, in agreement with Eq. (3.38).This general decreasing trend visible in Fig. 9 is easy to understand: the closer we are to the AS solution initially, the more time it will take the deviation to grow until they become of order unity.This dependence can be simply estimated as in Eq. (3.38), which is shown by the red dashed line in Fig. 9 and nicely reproduces the slope of the exact solution.The constant shift of this line can be explained by the fact that Eq. (3.38) underestimates the survival time since it does not account for the initial decreasing stage. 10Deviations from this dependence (small wiggles on the blue and green lines) occur because of the phase of the cosine function when the solution crosses the threshold value 10 −1/2 .The agreement with the analytical results confirms the robustness of our numerical techniques also for the small departures shown in the figure .In order to see how the survival time of the AS solution depends on the model parameters, the axion-vector coupling constant β and the initial production parameter ξ 0 , we perform a scan over this two-dimensional parameter space and present the results in the form of heatmap plots in Fig. 10.We fix the initial relative deviation in ξ to be (a) δξ 0 /ξ 0 = +10 −6 , (b) −10 −6 , (c) +10 −3 , and (d) −10 −3 , which are shown in the cor- responding panels of Fig. 10.The comparison between the two top and the two bottom panels again confirms the scaling of Eq. (3.38).For Re (ζ 1 ) ≃ 0.55, as indicated by Fig. 2, a variation of 10 3 in δξ 0 /ξ 0 provides a shift N AS ≃ 12.5, in good agreement with the various panels.More interestingly, each panel shows how the survival number of e-folds depends on the model parameters ξ 0 and β.This dependence is characterized by a wavelike pattern, meaning that the survival time changes non-monotonically with the increase of β or ξ 0 .This can be explained by the fact that the time dependence of δξ/ξ 0 is an oscillatory function with increasing amplitude [see Fig. 3(a)].Typically, a slight change in the model parameters leads to a small phase shift in the oscillations, and the curve for |δξ/ξ 0 | crosses the threshold value of 10 −1/2 at a slightly different time.There can, however, be situations where, after a small change in β or ξ 0 , the curve for |δξ/ξ 0 | does not reach the threshold during the same oscillation as before, but it has to evolve approximately half an oscillation period more to do this.This results in jumps in N AS , as can be clearly seen in Fig. 10.
Overall, our results in Figs. 9 and 10 corroborate our understanding of the relation between the growth rate of the fastest-growing mode, Re(ζ 1 ), and the survival time of the AS solution, N AS , and thus serve as another numerical validation of the analysis in Sec. 3.

C Self-correction algorithm for the GEF
At the end of Sec. 3, we discuss a challenge that we encounter when employing the GEF during the period of fast oscillations shown, e.g., in Figs. 4, 6, 7, and 8.For such a fast and non-monotonic change in ξ (or, equivalently, in the inflaton velocity), the cutoff momentum given by Eqs.(2.36) and (3.39) shows a sequence of plateaus in the time evolution.This is clearly seen in Fig. 11(a), plotted for the case β = 10 2.5 and constant background quantities H 0 = 2.7×10 11 GeV and ξ 0 = 6.85,where the red curve shows the expression 2 e N |ξ(N )|, and the green solid line is the upper envelope of this function, which is k h /H.
During these plateau stages, the underlying assumption that the spectral densities of E (n) , G (n) , and B (n) are dominated by the mode k h at large n, which allows us to truncate the GEF system at some finite order n cut , is violated.Indeed, let us consider the mode equation (2.18).In the case of constant H, it implies that, at a given moment of time N , the tachyonic instability occurs for modes with momenta k < 2He N |ξ(N )| . (C.1) During each of the plateau stages, all these momenta are less than k h , which equals the largest value of 2He N |ξ(N )| in all preceding moments of time.This means that the spectral densities are growing for modes with smaller momenta than k h , such that these modes may become of the same importance as k h in the integrals over the spectra.This introduces a numerical error in the last equations of the GEF, where the truncation is performed, which then quickly propagates through the system of equations, finally reaching the zeroth order.Note that increasing n cut does not help to avoid this problem but only postpones it to higher-order bilinear quantities.Therefore, in order to detect the situation where the GEF starts giving inaccurate results, one should always perform a consistency check using the mode-by-mode (MBM) solution.For this, one takes the time dependence of the scale factor a and the inflaton velocity φ from the GEF result [in the case of constant background quantities, one just needs to take the dependence ξ(N )], and solves the mode equation (2.17) or (2.42) for all modes that cross the horizon during the time interval used in the GEF.Then, using Eqs.(A.1)-(A.3),one can compute the bilinear functions that follow from the MBM solution and compare them to those from the GEF result.The relative deviations between them, where in the place of X one may take, e.g., the lowest-order bilinear gauge-field quantities E (0) , G (0) , and B (0) , can be used to estimate the consistency of the GEF solution.Note that this relative deviation is not a true numerical error of the GEF result, because the MBM solution that we use as a reference is not independent from the GEF solution, but it is based on the time dependence of ξ taken from the GEF result.We will, nevertheless, refer to it as an "error" in what follows.
The empty circles in Fig. 11(b) show the typical behavior of the relative error during the plateau in k h : at some point, it starts increasing and exceeds the 1 % threshold (shown by the purple dashed line).If one allows the GEF system to go further in time, the error reaches much greater values.Therefore, in order to control the accuracy of the GEF results, one needs to reinitialize it at the moment when the error exceeds the selected threshold.
To do so, one may use the spectra obtained by the MBM approach in order to compute the bilinear functions according to Eqs. (A.1)-(A.3).This helps to improve the situation and keep the error under control.The vertical gray lines in Fig. 11 mark the times at which the self-correction was performed.The error in the corrected result remains always less than the chosen threshold of 1 % during the whole duration of the simulation.
The algorithm underlying our self-correction procedure can be summarized as follows: (i) Perform numerical runs of the GEF equations for two different values of n cut (sufficiently large so that, for the time intervals without plateaus in k h , the results of both runs coincide)11 up to the time in which the results start to deviate.
(ii) Use the time dependence of ξ and a from the GEF result and solve the mode equation (2.17) in order to obtain the mode spectrum of the produced gauge fields.
(iii) Compute the zeroth-order bilinear quantities from the spectrum using the expressions in Eqs.(A.1)-(A.3)and find the relative error of the GEF result using Eq.(C.2).
(iv) When the error exceeds the set threshold, compute the values of the bilinear functions according to Eqs. (A.1)-(A.3)for n > 0.
(v) Use these new corrected values for the bilinear functions to reinitialize the GEF.
In order to avoid small jumps (by ∼ 1 %) in the zeroth-order quantities, E (0) , G (0) , and B (0) , which may lead to a short spurious stage of relaxation to a smooth solution just after the reinitialization, it is better to reinitialize only the bilinear quantities starting from n = 1, while keeping the old values for the zeroth-order quantities along with the values of ξ and a, which cannot be updated by the MBM approach.Finally, we comment on the choice of the threshold in the error.If one chooses a greater (smaller) threshold, less (more) frequent self-corrections are required.In practice, one therefore needs to find a compromise between the accuracy of the numerical result and the required computation time.
Open Access.This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
.15) where we have introduced the number of e-folds, N ≡ ln a. Writing δϕ in this way, we can identify Re(ζ) with the growth rate and Im(ζ) with the angular oscillation frequency of δϕ as a function of the number of e-folds.In fact, the most general linearized perturbation will be a linear combination of modes of the form(3.15),where each term has one allowed ζ n (with its complex conjugate) and the coefficients C n are obtained from the initial condition for the perturbation.Now, we set V ′′ = 0 (since in our case V ′ = const), substitute Eqs.(3.13) and (3.14) into Eq.(3.12), and perform the change of integration variables τ ′ → x ′ = −kτ ′ and k → x = −kτ .After that, the combination C (−τ ) −ζ−2 appears in all terms on the leftand right-hand sides.Canceling this τ -dependence, we obtain the following equation for ζ:

Figure 1 .
Figure 1.Eigenvalues of the LGEF system with n cut = 70 (red dots) and solutions of Eq. (3.18), which agrees with Eq. (3.14) in Ref. [1] (green circles), in the complex ζ plane for ξ 0 = 7 and β = 10 2.5 .The contour plot in the background shows the absolute deviation from equality in Eq. (3.18).In the white region, |R − 1| > 10 2 , exceeding the scale of the color code.For each solution ζ, also the complex conjugate ζ * is a solution.In this figure, only the solutions with positive imaginary parts are shown.

Figure 2 .
Figure 2. (a) Real and (b) imaginary parts of the eigenvalue ζ 1 , corresponding to the fastestgrowing mode, as functions of ξ 0 .The bands of finite width reflect different values of β ∈ [10 1.5 , 10 3.5].The blue band follows from the LGEF, while the green band shows the solution of Eq.(3.18).The red dashed line in panel (b) is the analytical estimate in Ref.[41]; see Eq.(3.34).
, the evolution is tracked until the deviation becomes of order unity.The survival number of e-folds N AS ≈ 29 is shown by the vertical gray dashed line.In panel (b), only the first 10 e-folds are shown.The blue dashed lines correspond to the LGEF solution (truncated at n cut = 70) taking into account, respectively, only the fastest-growing mode in panel (a), and three eigenmodes ζ 1 , ζ 2 , and ζ 5 in panel (b) (we recall that the eigenmodes are sorted by decreasing real part of ζ).

Figure 3 .
Figure 3.Time evolution of the relative deviation δξ/ξ 0 for an initial ξ value in the AS solution ξ 0 = 7, axion-vector coupling β = 10 2.5 , and initial fine-tuning (initial displacement) δξ 0 /ξ 0 = 10 −6 .(a) Full time interval, (b) zoom-in into the transition region from the initial decaying phase to the growing phase.The green solid lines show the solution obtained by the GEF truncated at order n cut = 151; the blue dashed lines correspond to the result of the LGEF truncated at order n cut = 70.In panel (a), only the fastest-growing mode ζ 1 is included; in panel (b), three modes, ζ 1,2,5 , are included.The vertical gray dashed line in panel (a) indicates the survival time of the AS solution when the relative deviation of ξ from the AS value ξ 0 reaches 10 −1/2 (horizontal dashed line).

Figure 5 .
Figure 5.Time evolution of (a) the parameter ξ and (b) the energy densities for the axion-vector coupling β = 10 2.5 in the realistic inflationary model featuring a steep inflaton potential with initial values of the Hubble parameter H 0 = 2.7 × 10 11 GeV and slow-roll parameter ϵ V,0 = 6.1 (corresponding to ξ 0 = 7.1).The blue solid lines show the true solution of the system found using the GEF, while the red dashed lines correspond to the enforced AS solution.The vertical dashed lines show the moment of time N = N AS when the linear perturbation theory breaks down (relative deviation of ξ from its initial value exceeds 10 −1/2 ).The vertical dotted lines show the end of inflation for the real system (blue) and for the enforced AS solution (red).

Figure 8 .
Figure 8.Time evolution of the parameter ξ in the realistic inflationary model (blue line) compared to the late-time behavior of ξ in the toy model with constant H (orange line) studied in Sec. 3. The curves are shifted in such a way that the moments of time when the condition (ξ − ξ)/ ξ = +10 1/2 is satisfied for the first time coincide (vertical dashed line), where ξ denotes either the value ξ 0 = 6.85 in the case of constant background quantities or the time-dependent ξ value according to the EAS solution in the case of time-dependent background.The green dashed-dotted horizontal line shows the initial value ξ 0 = 6.85 for both curves.

Figure 9 .
Figure 9. Survival time of the AS solution as a function of the initial fine-tuning δξ 0 /ξ 0 for ξ 0 = 7 and axion-vector coupling β = 10 2.5 .The blue solid and green dotted lines correspond to positive and negative initial deviations, respectively.The red dashed line shows the analytical estimate of the survival time in Eq.(3.38), which is based on the growth rate of the fastest-growing mode.

Figure 11 .
Figure 11.Time dependences of (a) the cutoff momentum k h and (b) the relative error of the GEF result for ρ E , ρ B , and ρ EB compared to the mode-by-mode (MBM) solution for β = 10 2.5 , ξ 0 = 6.85, and H = 2.7×1011 GeV.Gray vertical lines show the moments of time when the selfcorrection procedure has been applied.The empty circles show the relative error of the GEF result compared to the MBM solution that one finds when no self-correction procedure is applied.

Table 1 .
.48) Sec.3: Constant gradient V ′ , constant Hubble rate H, Friedmann equation ignored Overview of the solutions studied in Sec. 3 and 4 and the methods used to obtain them.