Perturbatively including inhomogeneities in axion inflation

Axion inflation, i.e. an axion-like inflaton coupled to an Abelian gauge field through a Chern-Simons interaction, comes with a rich and testable phenomenology. This is particularly true in the strong backreaction regime, where the gauge field production heavily impacts the axion dynamics. Lattice simulations have recently demonstrated the importance of accounting for inhomogeneities of the axion field in this regime. We propose a perturbative scheme to account for these inhomogeneities while maintaining high computational efficiency. Our goal is to accurately capture deviations from the homogeneous axion field approximation within the perturbative regime as well as self-consistently determine the onset of the non-perturbative regime.


Introduction
Cosmic inflation remains the most attractive theory to explain the precise observations of the cosmic microwave background (CMB) by the Planck satellite [1,2].Among the particle physics models which can lead to a quasi-exponential expansion, considerable attention has been paid to axion-like particles as the driving field of inflation.The axion's angular symmetry is only broken by non-perturbative effects and thus the observationally required flatness of the potential can be ensured naturally [3].Furthermore, the shift symmetry allows for derivative couplings of the axion to the Chern-Simons density F µν F µν of a (dark) gauge field A µ .Such couplings can lead to the exponential production of A µ due to a tachyonic instability of one of the two helicity modes in the equations of motion (EOM), which is solely controlled by the axion velocity and hence generically impacts the later stages of inflation, corresponding to length scales much smaller than those accessible in CMB observations.The consequences of such a large non-thermal A µ population are diverse and include: i) an effective friction term in the axion EOM [4], ii) a strong enhancement of the scalar and tensor perturbations with possible observational consequences, such as the production of primordial black holes [5][6][7][8][9][10][11] and (chiral) stochastic gravitation waves [12][13][14][15][16][17][18] and iii) a mechanism for magnetogenesis [19][20][21][22][23][24] and baryogenesis [23,[25][26][27] if the gauge field is taken to be the Standard Model (SM) hypercharge.
Obtaining accurate predictions for any of these processes requires evolving the highly non-linear system containing the axion and gauge fields, as well as when present, any light fermions.In this paper we shall focus on the case where the gauge field is a dark photon, with no other couplings to the SM (or beyond) other than the Chern-Simons coupling to the axion. 1 In this case, the most important backreaction to consider is the effective friction induced by the gauge fields on the axion [4], and our interest will be in the regime where this backreaction is strong, i.e. typically towards the end of inflation.Changes in the axion velocity impact gauge field modes within the tachyonic instability window, which contribute to the friction force.As a result, the friction term reacts with some time delay to the changes in the axion velocity, leading to a resonantly coupled system with distinct peaks in the axion velocity [31,32].These results have been confirmed in perturbative stability analysis [33,34] as well as using the gradient expansion formalism (GEF) [30,34].The latter provides a remarkably efficient tool for numerical simulations based on expressing the non-linear EOMs in position space as a tower of ordinary differential equations (ODEs) for the 2-point functions of the axion, the gauge fields, and the gradients of the gauge fields, reducing the computation time by orders of magnitude compared to e.g. the iterative procedure used in [32].
All these methods, however, make a crucial assumption: they take the axion field to be homogeneous.Given the rapid growth of the axion perturbations, the significant departure from the standard slow-roll regime and the strong non-linearities involved, it is not surprising that this approximation breaks down in the strong backreaction regime.This has recently been explicitly demonstrated in a lattice simulation [35] using CosmoLattice [36,37], reproducing earlier results when switching off axion gradients but finding significant departures when the axion gradients are taken into account consistently.While accurately dealing with the full non-linear problem including the strong backreaction, the downside of these simulations is that they are extremely costly.Moreover, as can be seen from the results obtained in [35], the highly non-linear dynamics implies that observable quantities (such as the magnitude of the gauge field contribution or the scalar perturbations at a given scale) are not monotonic functions of the axion gauge field coupling.A full exploration of the phenomenology of axion inflation throughout the parameter space thus seems very costly at best based on lattice simulations only.
Our goal in this paper is to leverage the benefits of the highly efficient GEF to perturbatively include axion gradients.To obtain a closed set of ODEs, this requires evolving not only the 2-point functions but also higher p-point functions under the GEF scheme.We are particularly interested in the regime where axion gradients are relevant and impact the evolution of the axion vacuum expectation value and the gauge field distribution, while still allowing for a perturbative treatment.As time evolves, this will typically give way to a regime in which the axion gradients become too large to be treated perturbatively, calling for a lattice simulation.Within our formalism, we self-consistently determine the breakdown of perturbativity, providing a tool to compute initial conditions for lattice simulations, focusing their computational power on the truly non-linear regimes.Our work should be seen as a first step in developing this methodology, and we discuss possible extensions and scalability.In the process, we gain new insights into the application of the GEF to axion inflation, notably regarding the need to go to rather high order in the GEF tower to obtain convergence as well as an improved truncation relation.
The paper is organized as follows.In Sec. 2 we briefly review the GEF and derive an extension including axion gradients.Some technical details, in particular the more lengthy equations for the 3-point functions are given in App.A, while App.B focuses on the limitations of the GEF and in particular the truncation relation.Our results are presented in Sec. 3, and contrasted with results obtained in lattice computations as well as under the assumption of a homogeneous axion field.The final Sec. 4 summarizes and discusses the results.

The gradient expansion formalism including axion inhomogeneities
The gradient expansion formalism developed in Refs.[30,38] for axion inflation (see [39,40] for earlier work in related contexts) provides a computationally efficient way of accounting for the backreaction of the gauge fields on the axion dynamics.The system of interest here is an unbroken, dark U (1) gauge group coupled to the axion via a Chern-Simons interaction, where F µν = ϵ µνρσ F ρσ /2 with ϵ 0123 = 1 and for simplicity we will set V (ϕ) = m 2 ϕ ϕ 2 /2.The GEF re-arranges the resulting coupled partial differential equations governing the dynamics of the electric and magnetic field (i.e.Maxwell's equations in the presence of an axion photon coupling) into a tower of linear ODEs for the 2-point functions with X, Y = {E, B}, the bracket indicating the spatial average, and a indicating the scale factor of the expanding universe.Under the assumptions that the axion field is homogeneous and that its velocity varies only slowly, the EOMs of the 2-point functions form a closed system, and the infinite set of equations can be truncated at finite power n of the curl yielding an efficient and accurate evaluation of the axion and gauge field dynamics [30].However, as we will discuss in more detail below, typically neither of these assumptions are fulfilled once the gauge field backreaction becomes significant.
In the strong backreaction regime, the axion field enters into a oscillatory regime, understood as resonances resulting from the time-delayed friction force exerted by the gauge fields [30][31][32][33].More recently lattice studies have qualitatively confirmed the existence of the oscillatory behavior, finding however quantitatively significant differences (notably a significant damping of these oscillations) when including the inhomogeneities in the axion field [35].
The goal of the present paper is to extend the GEF beyond these two assumptions.Below, we derive an extended version of the GEF including axion inhomogeneities (i.e.axion gradient terms) in a perturbative manner.The resulting non-linearities in the EOMs of the 2-point functions prompt us to include higher p-point functions in order to reduce our equations to ODEs, in a similar spirit as the original GEF.Similarly, we will need to find a suitable procedure to truncate this second expansion series as finite (and low) p.In the process, we will shed light on the role of rapid changes in the axion velocity and the resulting limitations of the GEF.

Equations of motion
We start from the exact EOMs as derived from Eq. (2.1), separating the homogeneous component of the axion ϕ(t) from its inhomogeneous component χ(t, ⃗ x), ) for the matter sector with β = αM P /(πf a ), and for the gravity sector.For χ = 0, multiplying Eqs.(2.5) and (2.6) with {( ⃗ ∇×) n ⃗ E, ( ⃗ ∇×) n ⃗ B} yields a tower of ODEs linear in the (scalar) 2-point functions (2.2).To include the axion gradients, the structure of the last two terms in Eq. (2.5) as well as the last term in the Gauss constraint in Eq. (2.7) suggests the inclusion of 3-point functions including one power of either χ or ⃗ ∇χ.Using the EOMs, the ODEs governing these 3-point functions in turn depend on 4-point functions, etc.The result is a double expansion in gradients (∇×) n and p-point functions.
To first order in this expansion, we will only keep 3-point functions with up to one spatial derivative.In this case, the EOMs of the 2-point electromagnetic functions P with n ≥ 2 are the same as in the usual GEF formalism [30], . (2.12) Here, the boundary terms on the right-hand side account for the change in the number of modes which have been excited from the vacuum, see below for a more detailed discussion.The superscript (n) refers to the number of curls, indexing the GEF tower.The EOMs for 2-point functions for the electromagnetic fields for n = {0, 1} now contain 3-point functions B as anticipated, and EB − 2β φ Ṗ( 1) EB − P B − where we have defined with f = χ, χ and n = {0, 1}.The two superscripts on the 3-point functions refer to the number of curls acting on ⃗ E and ⃗ B, respectively.The EOMs for these 3-point functions can be obtained analogously from Eqs. (2.5) and (2.6) and are given explicitly in App. A. Importantly, we note here the key approximations which enter in their derivation: Firstly, we keep only 3-point functions containing up to one spatial derivative.This is not a fundamental limitation of the proposed method but should rather be seen as the first order in a gradient expansion series.Including higher order terms will lead to more (and more lengthy) equations, with the number of equations scaling at most2 as n3 B with n B the maximal number of derivatives appearing in the 3-point function.However, the overall structure of the equations remains unchanged and we expect this to be numerically tractable.The results derived in this way are reliable as long as the subsequent order in this expansion is sufficiently small, and as discussed below, we will use this as a criterion to self-consistently check the validity of our method.Secondly, we factorize the resulting 4-point functions appearing in the EOMs of the 3-point functions into products of 2-point functions assuming Gaussian distributions for the electromagnetic fields. 3We moreover set 4-point functions that involve one factor of ⃗ ∇χ to vanish since this leaves a spatial index uncontracted within the 2-point function, which would be in violation of statistical isotropy.
Finally the 2-point functions of the axion fluctuations evolve as with

Boundary terms and truncation
Two subtleties in the derivation above deserve a more detailed discussion: the boundary terms in Eqs.(2.13) to (2.18) and the truncation of the gradient expansion series for the 2-point functions at finite n.Both of these are a priori not related to the inclusion of the axion inhomogeneities as they appear only in the EOMs of the 2-point functions.In fact, for a homogeneous axion field, these are the only two approximations by which the GEF prescription deviates from an exact solution and this is where the approximation of a slowly varying axion velocity enters.However, while for a homogeneous axion field the gauge field power spectra P (0) E,B,EB have been shown to be robustly reproducing the results found solving the mode equations of the gauge field [30] and are moreover in very good agreement with lattice results after setting the axion gradient terms to zero [35], we find the impact of these approximations on the higher orders in the GEF tower P (n) E,B,EB to be significant (see also [34]).We will show that these difficulties can be mitigated by an improved truncation relation.Including axion gradients, we observe that this ensures sufficient stability in the algorithm within the perturbative regime of axion inhomogeneities.
Boundary terms.In terms of the mode functions of the vector potential A, the gauge field 2-point functions are given as (2.25) where σ encodes the two gauge field polarizations and the Heaviside function ensures the vacuum subtraction, i.e. that only modes with k < k h which have been exited out of their vacuum state contribute to the regularized integral.To determine k h , we refer to the EOM for A k (τ ), which encounters a tachyonic instability for the polarization σ = λ for Here we assume that the axion inhomogeneity does not affect the gauge boson dynamics at the subhorizon scales.Since the gauge bosons are excited only after exiting the horizon, the axion inhomogeneity produced from the gauge bosons is also expected to be outside the horizon, justifying our assumption.Taking this into account leads to boundary terms in the EOMs for the 2-point functions, which are explicitly given as To evaluate these, we follow the prescription given in [30] which is based on the solutions for the mode functions for constant ξ.The boundary terms thus arise from the need to regularize the vacuum contribution, and in this implementation, rely on an at most slowly varying axion velocity.The final results are not particularly sensitive to the value of the cut-off k h since the change of the integration range in Eqs.(2.24)-(2.26)by changing k h is compensated by the change of the boundary term.For example, reducing k h by 25% induces changes which are smaller than 10% (and largely smaller than 1%) in ξ.
Truncation relation.The infinite tower of equations (2.10) to (2.12) needs to be truncated at some finite n to implement the GEF numerically.The truncation relation employed in [30] is where X = E, B or EB.This can be derived from Eqs. (2.29) to (2.31) assuming a power law spectrum for |A σ (k)| 2 around k h , as proposed in Ref. [30].However, both results from lattice simulations [35] as well as the results obtained using solving the gauge field mode equations show that there can be a lot of structure in the spectrum at this scale, leading to significant deviations in the ratio P (n+1) X /P (n−1) X , as shown in App.B. However, as we show in that appendix, Eq. (2.32) can also be obtained (without invoking the power law approximation of the spectrum) as the asymptotic large-n limit under the assumption that ξ is constant.As we show there, a value of n ∼ 55 is necessary for the relation (2.32) to hold up to about 5%.To our understanding, this explains why values of n ∼ 100 are necessary to achieve convergence in the GEF, whereas n ≫ 1 should have been sufficient based on the assumption of power law spectrum around k h .
This in turn indicates that rapid changes in ξ will induce errors in this truncation relation, which propagate through the coupled system of equations down to low n.This is supported by our observations in the numerical studies shown in App.B and was also recently observed in Ref. [34].When the axion velocity drops rapidly, the P (n) X take (unphysically) large values at large n, which over time propagate down to lower n modes.If the phase of rapidly changing (in particular dropping) axion velocity is sufficiently long, this can in principle impact the observables, i.e. the P (0) X power spectra.This calls for an improvement of the truncation relation (2.32) in the further development of the GEF formalism.We show in App.B results obtained using not only where P(n) X /H 4 0 (k h /a) n .For L = 1 we trivially recover Eq. (2.32).As we show in App.B for L = 4 and L = 10 this improves the stability of the system to a point which is sufficient for the study of the regime of perturbative axion inhomogeneities.

Validity of the axion gradient expansion
In the scheme described above we included only terms with up to one power of the axion gradient.To ensure the consistency of this approach, we monitor the second order axion derivatives sourced in this manner (without including their backreaction on the truncated system described above).More precisely, we monitor the ratio of the gradient to the kinetic energy, As long as this quantity is small, it is justified to drop higher powers of the axion gradients, whereas R χ above O(0.1 -1) indicates that the axion gradients can no longer be treated perturbatively, calling for a full lattice simulation.To compute (∇χ) 2 = − χ∇ 2 χ , the relevant EOMs for the axion 2-point functions to second order in the gradient expansion are χ + 2m 2 ϕ P (2)  Evolution of ξ assuming a homogeneous axion field (dashed black) and perturbatively including axion gradients (red) for different values of β.The light (dark) gray region indicates that the gradient energy of the axion exceeds 1 % (50 %) of the kinetic energy, while the gray vertical line corresponds to 5 %.Wherever possible we compare to the result of the lattice simulation [35].
which in turn require the evaluation of the 3-point function B (2;0,0) f ;EB , whose definitions and equations are again given in the App. A. The breakdown of this perturbative expansion scheme R χ > 0.5 is indicated by the dark gray region in the figures below.

Numerical results
Figs. 1 and 2 show the results of the methodology described in Sec. 2 for values of the coupling β ranging from 15 to 25.The six panels of Fig. 1 focus on the evolution of the parameter ξ = β| φ|/(2HM P ).In all panels, the dotted black lines show the result for the GEF assuming a homogeneous axion and the solid red lines are new results including perturbatively the axion gradients to first order.Where available, we show for comparison the lattice results obtained in Ref. [35] in dashed blue.The gray regions and horizontal lines indicate the quality of the perturbative expansion.In the dark gray region the axion gradient energy exceeds 50% of the axion kinetic energy, R χ > 0.5, indicating the non-perturbative regime.For all panels, we show only 0.25 e-folds of this non-perturbative regime, though we stress that lattice results have shown that inflation can last several e-folds longer.Defining ∆N = N end − N (R χ = 0.5) leads to ∆N ≃ 1.5 for β = 15, ∆N ≃ 3.5 for β = 18, and ∆N ≃ 4.5 for β = 20 [35].The light gray region indicates the regime in which the gradient energy is sizeable, i.e. 0.01 < R χ < 0.5, but our perturbative treatment is still valid.The vertical gray line indicates R χ = 0.05.We observe that below this value, we recover the full lattice results to good agreement (while the deviation from the homogeneous approximation is already significant).Above this value the leading order correction implemented here is insufficient to fully reproduce the lattice results, but since R χ ≪ 1 a systematic expansion to higher orders in the axion gradients might conceivably achieve this (within the light gray region).
Fig. 2 shows the evolution of the kinetic (ρ K ), potential (ρ V ), electromagnetic (ρ EM ) and axion  gradient energy (ρ G ) defined as for the same values of β.Here we include ρ G in the total energy density, ρ tot = ρ K + ρ V + ρ EM + ρ G , although we do not include ρ G in the Friedmann equation as we explained above. 4This figure gives a more detailed view of the importance of the axion gradient energy, showing that in many cases (e.g.β = 20, 22, 25) there is a prolonged regime in which the axion gradient energy remains at the percent level or below.The opposite is observed for β = 15, where once it becomes relevant, the axion gradient energy rapidly becomes comparable to the other components.The efficiency of the method discussed in Sec. 2 allows us to study a wide range of couplings β, confirming that diverse and complex dynamics occurring in the strong backreaction regime.In particular, possible observable signatures related to the different energy components depend on the coupling β in a non-monotonic way.
As one example of possible observable consequences we show the scalar power spectrum sourced by the axion fluctuations, ∆ 2 ζ = (H/ φ) 2 P (0) χ in Fig. 3. 5 Note that this contains only the contribution from the axion fluctuations, and not the contribution from the gauge field energy density fluctuations or metric fluctuations.We note that even within the perturbative regime of axion gradients we significantly exceed the threshold for primordial black hole (PBH) formation which is estimated to lie around ∆ 2 ζ ∼ 10 −2 for a Gaussian distribution and around ∆ 2 ζ ∼ 10 −4 for the non-Gaussian spectrum assumed in Ref. [5] for axion inflation.PBHs generated in the last few e-folds of inflation are very light and decay rapidly through Hawking radiation leaving no traces in our cosmological history.However, as we 4 Including this term, despite spoiling the consistency of the Friedmann equation, does not affect our result in the regime where we can trust our perturbative treatment of the axion gradient. 5Strictly speaking P (0) χ is an integrated quantity and we do not have access to the power spectrum for each different momentum k in our formalism (unless we compute the axion 2-point functions with sufficiently high powers of the spatial derivatives).However, we expect that at each given time the modes with k ∼ k h dominate the integral and the time evolution of P (0) X roughly corresponds to the spectrum.see from Fig. 3, even for a moderate coupling of β = 15 the non-perturbative regime of axion gradients is reached before the end of inflation, making a precise prediction of the PBH spectrum impossible with perturbative methods.It is nevertheless interesting to note that even for such small values of β (which, in particular, do not produce excessive gravitational waves during preheating [41,42]) our results indicate a phase of PBH formation.Increasing the coupling β will typically extend this phase, leading to more massive PBHs with longer lifetimes.Interestingly, for some parameter choices such as β = 20, displayed in the right panel, we find indication for scalar perturbations above the PBH formation threshold for a significant range of e-folds within the perturbative regime.

Discussion and outlook
The key result of this work is the extension of the gradient expansion formalism (GEF) to perturbatively include axion gradients when evaluating the dynamics of axion inflation.This involves applying the GEF not only to the 2-point functions but also to higher p-point functions, resulting in an expansion in (i) the number of derivatives n acting on the gauge fields in the 2-point function as in the original GEF, (ii) the order of p of the highest correlation function which is not factorized into lower p-point functions and (iii) the number n p of derivatives in those p-point functions.Here we consider the leading order correction due to axion gradients, which requires p = 3 and This captures several important aspects of the dynamics of the system.For very small values of the expansion parameter R χ we recover the results of the homogeneous GEF approximation (i.e.assuming the axion field to be perfectly homogeneous but allowing for gradients in the gauge fields) which in this regime agree with the lattice results capturing the full non-linear dynamics.For R χ ≳ 1%, the homogeneous GEF results and the lattice results start to diverge, with our perturbative expansion recovering the lattice results.For R χ ≳ 5% the leading order corrections included here no longer suffice to accurately track the system, though given the smallness of R χ a systematic inclusion of higher order terms might conceivably achieve this.We see no fundamental obstruction to this in this regime.In particular, the number of the 3-point functions scales as n 3 B (or milder by using e.g. the integration by parts) and this eventually limits the inclusion of higher order terms.However, this happens only when n B is sufficiently large and the computational time does not change drastically by a mild increase of n B (note that the number of the 2-point functions we numerically follow is ∼ 750 for our choice of n max = 250).Finally, for R χ ≳ 0.5 perturbativity is violated and a full lattice simulation seems unavoidable.Our method allows to very efficiently identify the regions and provide accurate initial conditions for them, thereby enabling future lattice simulations to focus their computational power on these non-perturbative regimes.
An accurate calculation of the strong backreaction regime in axion inflation is crucial to exploit its rich phenomenology and conclusively test this inflation model.This includes the remaining duration of inflation once the system enters the strong backreaction regime (see Ref. [35]), the magnitude of peaks in the scalar power spectrum as well their position with regards to the end of inflation which are crucial to determine if a sizable amount of primordial black hole are formed and their mass distribution (see Ref. [5]) as well as the anisotropic component of the gauge field energy momentum tensor which will determine the gravitational wave spectrum (see Ref. [18]).Such results will need to be contrasted with model constraints coming from the gravitational wave production in the preheating era [41,42].Assuming a simple shape of the scalar potential throughout the inflation and preheating era, these impose stringent bounds on the axion gauge field coupling, which however still allow for an inflationary phase within the strong backreaction regime.Currently, only costly lattice simulations can give quantitative answers to all these questions.
The method proposed here provides a first step towards an efficient way of studying these questions across the parameter space of axion inflation.Several extensions of the scheme proposed here deserve further study, e.g. the inclusion of higher order p-point functions and/or derivatives therein as well as including a correction algorithm as proposed in [34].We hope that this work will trigger some of these developments.The left (right) panel depicts the exemplary case of ξ = 4 (ξ = 6), the horizontal line indicates the truncation relation assumed in the GEF in [30].
following Ref.[35].The initial Hubble parameter H 0 is computed as Throughout this paper, we follow the electromagnetic 2-point function up to n = n max = 250, and express those with n = n max + 1 by the truncation relation.

B Truncation relation
One of the key approximations in the gradient expansion formalism is the truncation of the tower of equations (2.10)-(2.12)through the relation (2.32).In this appendix we demonstrate that this relation holds asymptotically for large n and approximately constant ξ.In contrary to the derivation given in Ref. [30], this does not rely on approximating the gauge field spectra by a power law around k h , and explains why values of n = O(100) are required to obtain accurate results in the GEF.It also illustrates how when ξ changes rapidly, imposing this truncation relation introduces sizable errors, as observed also in Ref. [34].Finally, we will introduce an improved truncation relation which mitigates some of these effects.For simplicity, the discussion in this appendix is focused on the application of the GEF assuming a homogeneous axion field, as in Ref. [30].The arguments can immediately be extended to include axion gradients as described in Sec. 2. Let us begin by studying the truncation relation assuming constant ξ.In the slow-roll approximation with constant ξ the gauge field A k (τ ) mode function can be determined analytically by solving where we defined ω ± = √ 2kA ± k (τ ) and x = −kτ , see Eq. (2.27).Matching to the Bunch-Davies vacuum fixes The mode A − k does not experience a tachyonic enhancement and can thus be safely neglected.A reasonable approximation of the mode function A + k in the IR regime (x → 0), which dominates the contribution to the n-point functions, is given by For the n-point function of the E-field, see Eq. (2.24), the relevant quantity is The integral of Eq. (2.24) can then be done analytically and evaluates to The ratio of E is thus where as in the main text, we have introduced the notation P(n) X ≡ P X /H 4 0 (k h /a) n .For large n, and any ξ, the ratio of Eq. (B.6) asymptotically reaches unity, with a deviation of 5 % for n = 45 (n = 55) for ξ = 4 (ξ = 6), see Fig. 4. This, on the one hand confirms the truncation relation chosen in [30] in the constant ξ limit, while on the other hand explains why large values of n are required in order for this truncation relation to become relatively accurate.Similar conclusion can also be drawn for P(n) B and P(n) EB .The truncation relation remains a good approximation for slowly varying ξ, as can be seen in Fig. 5 which displays the result of a fully numerical solution of the GEF.In particular, we show the ratio of P(n+1) E / P(n−1)

E
for fixed e-fold N and β = 15.The left panel shows an example from the weak backreaction regime, where ξ changes only slowly and hence the truncation relation is relatively accurate at large n.On the other hand, the right panel shows the situation when the inflaton field velocity begins to undergo a more rapid change.In this case, we observe that the truncation relation remains a good approximation for a large range of n ≳ 50, although at n ∼ n max = 250 the onset of a violation of this relation can be observed.
This violation becomes more dramatic when inflaton velocity changes, and in particular drops, rapidly.This is displayed by the dotted black curve in Fig. 6  If left unchecked, this will finally affect the n = 0 components of the GEF tower.Similar effects were recently observed in Ref. [34], and similarly to those results we find that these difficulties occur when ξ drops rapidly and k h features a plateau, see Eq. (2.28), and the boundary terms are zero.To counter this issue, we propose the improved truncation relation (2.33).Averaging over several P(n+1−2l) E to determine the truncation relation for P(n+1) gives a more robust procedure, as shown by the colored curves in Fig. 6.Compared to the original truncation relation (dotted black), we see that the striking and unphysical feature ('bump') at large n has largely disappeared, at no additional computational cost.Nevertheless, smaller unphysical features ('wiggles') remain, which over time can numerically destabilize the system (see right panel).In particular, P (2n) E should always be positive, while the 'wiggles' contain negative valued P EB as well as for odd powers of gradients (although in the latter two cases positivity is not necessarily required).We were not able to conclusively determine the origin of these 'wiggles' nor find a strategy to remove them.To illustrate their relevance, the grey shaded regions  As can be seen, this typically happens after a sharp drop in ξ.For a milder evolution of ξ, e.g. for β = 18 or 20, this issue does not arise.In this context, it may be interesting to further study the proposal given in Ref. [34] based on re-initializing the GEF through the mode-by-mode method (which in turn only requires the input of the n = 0 mode).Since it takes some time for these unphysical effects to propagate to the n = 0 mode, this procedure might further improve the situation.To our understanding, this method was initially proposed to remove the 'bump' feature, which Eq. (2.33) achieves more efficiently.However, a similar method may prove useful to address the remaining issue of the 'wiggles'.
So far, the discussion in this appendix has focused on the homogeneous axion case.Despite the difficulties mentioned above, in practice, for the quadratic scalar potential studied here, the numerical issues in the higher orders of the GEF tower do not propagate to the lowest orders for any coupling β considered here before the end of inflation (see however [34] for different background dynamics).Moreover, including the axion gradients, the remaining difficulties with maintaining a stable truncation relation only occur in the non-perturbative regime, i.e. for R χ > 0.5, for all values of β considered, and do hence not seem to pose a problem within the region of validity of the perturbative method proposed here.In this sense, the methods employed here ensure sufficient stability of the GEF scheme to study the perturbative regime of axion gradients which is the main goal of this work.
In summary, this appendix clarifies the origin of the truncation relation of the GEF as the asymptotic limit of the GEF tower for large n and approximately constant ξ.This sheds some light on the limitations of the GEF formalism for rapidly varying ξ, and prompts the introduction of an improved truncation relation.This increases the stability and range of validity of the GEF approach.We point out a remaining instability which however in practice does not impact the results of this work as it does not occur within the perturbative regime of axion inhomogeneities.

Figure 1 :
Figure1: Evolution of ξ assuming a homogeneous axion field (dashed black) and perturbatively including axion gradients (red) for different values of β.The light (dark) gray region indicates that the gradient energy of the axion exceeds 1 % (50 %) of the kinetic energy, while the gray vertical line corresponds to 5 %.Wherever possible we compare to the result of the lattice simulation[35].

Figure 2 :
Figure 2: Evolution of the energy densities for different values of β.Gray bands as in Fig. 1. Results from lattice simulations for β = 15, 18 and 20 from Ref. [35] are shown in dashed.

Figure 4 :
Figure 4: For constant values of ξ the truncation relation (2.32) holds asymptotically for large n, see Eq. (B.6).The left (right) panel depicts the exemplary case of ξ = 4 (ξ = 6), the horizontal line indicates the truncation relation assumed in the GEF in[30].

Figure 5 :
Figure 5: Truncation relation (2.32) evaluated numerically for fixed β = 15, showing good agreement for n ≳ 100 when the backreaction is weak (left panel) and mild (right panel).The horizontal line indicates the truncation relation assumed in the GEF.

Figure 6 :
Figure 6: GEF tower for rapidly changing ξ using the truncation relation (2.32) (dotted black) and the improved truncation relation (2.33) with involving the highest four (solid blue) and ten (dashed red) even powers of the GEF tower.

Fig. 7
indicate the regions in which

Figure 7 :−
Figure 7: Evolution of ξ for β = 15 (left) and β = 22 (right), assuming a homogeneous axion field.The gray regions indicate a violation of the truncation relation above 10%, see text for details.
, indicating that the truncation relation (2.32) is well satisfied.However, at larger values of n we see very large deviations from the truncation relation, propagating over time down to lower values of n. E