compat=1.0.0

IFIC/23-45

FTUV-23-1005.0503

UMN-TH-4226/23

FTPI-MINN-23-18

CERN-TH-2023-186

Perturbatively including inhomogeneities
in axion inflation

Valerie Domckea, Yohei Emab,c, and Stefan Sandnerd

a | _Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland_
---|---
b | _William I. Fine Theoretical Physics Institute, School of Physics and Astronomy,_
| _University of Minnesota, Minneapolis, MN 55455, USA_
c | _School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA_
d | _Instituto de Física Corpuscular, Universitat de València and CSIC,_
| _Carrer del Catedrátic José Beltrán Martinez 2, 46980 Paterna, Spain_

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.

###### Contents

  1. 1 Introduction
  2. 2 The gradient expansion formalism including axion inhomogeneities
    1. 2.1 Equations of motion
    2. 2.2 Boundary terms and truncation
    3. 2.3 Validity of the axion gradient expansion
  3. 3 Numerical results
  4. 4 Discussion and outlook
  5. A Equations of motion for $3$-point functions
    1. A.1 Equations of motion up to one spatial derivative
    2. A.2 Equations of motion with two spatial derivatives
    3. A.3 Numerical implementation
  6. B Truncation relation

##  1 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_{\mu\nu}\tilde{F}^{\mu\nu}$ of a (dark) gauge field $A_{\mu}$. Such
couplings can lead to the exponential production of $A_{\mu}$ 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_{\mu}$ 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 [25, 26, 27, 23] 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.111For the
discussion of the SM case including light fermions see [28, 29, 30]. 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 unfeasible 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.

##  2 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,

$\displaystyle S$ $\displaystyle=\int
d^{4}x\sqrt{-g}\left[\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V(\phi)-\frac{1}{4}g^{\mu\rho}g^{\nu\sigma}F_{\mu\nu}F_{\rho\sigma}-\frac{\alpha\phi}{4\pi
f_{a}}\frac{1}{\sqrt{-g}}F_{\mu\nu}\tilde{F}^{\mu\nu}\right]\,,$ (2.1)

where $\tilde{F}^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}/2$ with
$\epsilon^{0123}=1$ and for simplicity we will set
$V(\phi)=m_{\phi}^{2}\phi^{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

$\displaystyle\mathcal{P}_{X}^{(n)}$
$\displaystyle=\frac{1}{a^{n}}\left\langle\vec{X}\cdot(\vec{\nabla}\times)^{n}\vec{X}\right\rangle\,,\quad\mathcal{P}_{XY}^{(n)}=-\frac{1}{a^{n}}\left\langle\vec{X}\cdot(\vec{\nabla}\times)^{n}\vec{Y}\right\rangle\,.$
(2.2)

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 [31, 32, 33, 30]. More recently
lattice studies have qualitatively confirmed this result, 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.

###  2.1 Equations of motion

We start from the exact EOMs as derived from Eq. (2.1), separating the
homogeneous component of the axion $\phi(t)$ from its inhomogeneous component
$\chi(t,\vec{x})$,

$\displaystyle 0$
$\displaystyle=\ddot{\phi}+3H\dot{\phi}+m_{\phi}^{2}\phi-\frac{\beta}{M_{P}}\left\langle{\vec{E}\cdot\vec{B}}\right\rangle\,,$
(2.3) $\displaystyle 0$
$\displaystyle=\ddot{\chi}+3H\dot{\chi}-\frac{\nabla^{2}\chi}{a^{2}}+m_{\phi}^{2}\chi-\frac{\beta}{M_{P}}\left(\vec{E}\cdot\vec{B}-\left\langle{\vec{E}\cdot\vec{B}}\right\rangle\right)\,,$
(2.4) $\displaystyle 0$
$\displaystyle=\dot{\vec{E}}+2H\vec{E}-\frac{1}{a}\vec{\nabla}\times\vec{B}+\frac{\beta}{M_{P}}\left(\dot{\phi}+\dot{\chi}\right)\vec{B}+\frac{\beta}{M_{P}}\frac{1}{a}\vec{\nabla}\chi\times\vec{E}\,,$
(2.5) $\displaystyle 0$
$\displaystyle=\dot{\vec{B}}+2H\vec{B}+\frac{1}{a}\vec{\nabla}\times\vec{{E}}\,,$
(2.6) $\displaystyle 0$
$\displaystyle=\vec{\nabla}\cdot\vec{E}+\frac{\beta}{M_{P}}\vec{\nabla}\chi\cdot\vec{B}\,,\quad
0=\vec{\nabla}\cdot\vec{B}\,,$ (2.7)

for the matter sector with $\beta=\alpha M_{P}/(\pi f_{a})$, and

$\displaystyle H^{2}$
$\displaystyle=\frac{1}{3M_{P}^{2}}\left\langle\frac{1}{2}\left(\dot{\phi}^{2}+\dot{\chi}^{2}\right)+\frac{(\partial_{i}\chi)^{2}}{2a^{2}}+\frac{m_{\phi}^{2}}{2}\left(\phi^{2}+\chi^{2}\right)+\frac{1}{2}\left(|{\vec{E}}|^{2}+|{\vec{B}}|^{2}\right)\right\rangle\,,$
(2.8) $\displaystyle\dot{H}$ $\displaystyle=-\frac{1}{6M_{P}^{2}}\left\langle
3\left(\dot{\phi}^{2}+\dot{\chi}^{2}\right)+\frac{(\partial_{i}\chi)^{2}}{a^{2}}+2\left(|{\vec{E}}|^{2}+|{\vec{B}}|^{2}\right)\right\rangle\,,$
(2.9)

for the gravity sector. For $\chi=0$, multiplying Eqs. (2.5) and (2.6) with
$\\{(\vec{\nabla}\times)^{n}\vec{E},(\vec{\nabla}\times)^{n}\vec{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 $\dot{\chi}$ or
$\vec{\nabla}\chi$. 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 $(\nabla\times)^{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 ${\cal P}$ with $n\geq 2$ are the same as in the
usual GEF formalism [30],

$\displaystyle\dot{\mathcal{P}}_{E}^{(n)}+(n+4)H\mathcal{P}_{E}^{(n)}-\frac{2\beta\dot{\phi}}{M_{P}}\mathcal{P}_{EB}^{(n)}+2\mathcal{P}_{EB}^{(n+1)}=\left[\dot{\mathcal{P}}_{E}^{(n)}\right]_{b}\,,$
(2.10)
$\displaystyle\dot{\mathcal{P}}_{B}^{(n)}+(n+4)H\mathcal{P}_{B}^{(n)}-2\mathcal{P}_{EB}^{(n+1)}=\left[\dot{\mathcal{P}}_{B}^{(n)}\right]_{b}\,,$
(2.11)
$\displaystyle\dot{\mathcal{P}}_{EB}^{(n)}+(n+4)H\mathcal{P}_{EB}^{(n)}-\mathcal{P}_{E}^{(n+1)}+\mathcal{P}_{B}^{(n+1)}-\frac{\beta\dot{\phi}}{M_{P}}\mathcal{P}_{B}^{(n)}=\left[\dot{\mathcal{P}}_{EB}^{(n)}\right]_{b}\,.$
(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 ${\cal B}$ as anticipated,

$\displaystyle\dot{\mathcal{P}}_{E}^{(0)}+4H\mathcal{P}_{E}^{(0)}+2\mathcal{P}_{EB}^{(1)}-\frac{2\beta\dot{\phi}}{M_{P}}\mathcal{P}_{EB}^{(0)}-\frac{2\beta}{M_{P}}\mathcal{B}_{\dot{\chi};EB}^{(0)}=\left[\dot{\mathcal{P}}_{E}^{(0)}\right]_{b}\,,$
(2.13)
$\displaystyle\dot{\mathcal{P}}_{B}^{(0)}+4H\mathcal{P}_{B}^{(0)}-2\mathcal{P}_{EB}^{(1)}=\left[\dot{\mathcal{P}}_{B}^{(0)}\right]_{b}\,,$
(2.14)
$\displaystyle\dot{\mathcal{P}}_{EB}^{(0)}+4H\mathcal{P}_{EB}^{(0)}-\mathcal{P}_{E}^{(1)}+\mathcal{P}_{B}^{(1)}-\frac{\beta\dot{\phi}}{M_{P}}\mathcal{P}_{B}^{(0)}-\frac{\beta}{M_{P}}\mathcal{B}_{\dot{\chi};B}^{(0)}-\frac{\beta}{M_{P}}\left(\mathcal{B}_{\chi;EB}^{(1,0)}-\mathcal{B}_{\chi;EB}^{(0,1)}\right)=\left[\dot{\mathcal{P}}_{EB}^{(0)}\right]_{b}\,,$
(2.15)

and

$\displaystyle\dot{\mathcal{P}}_{E}^{(1)}+5H\mathcal{P}_{E}^{(1)}+2\mathcal{P}_{EB}^{(2)}-\frac{2\beta\dot{\phi}}{M_{P}}\mathcal{P}_{EB}^{(1)}-\frac{2\beta}{M_{P}}\mathcal{B}_{\dot{\chi};EB}^{(1,0)}=\left[\dot{\mathcal{P}}_{E}^{(1)}\right]_{b}\,,$
(2.16)
$\displaystyle\dot{\mathcal{P}}_{B}^{(1)}+5H\mathcal{P}_{B}^{(1)}-2\mathcal{P}_{EB}^{(2)}=\left[\dot{\mathcal{P}}_{B}^{(1)}\right]_{b}\,,$
(2.17)
$\displaystyle\dot{\mathcal{P}}_{EB}^{(1)}+5H\mathcal{P}_{EB}^{(1)}-\mathcal{P}_{E}^{(2)}+\mathcal{P}_{B}^{(2)}-\frac{\beta\dot{\phi}}{M_{P}}\mathcal{P}_{B}^{(1)}-\frac{\beta}{M_{P}}\mathcal{B}_{\dot{\chi};B}^{(1)}=\left[\dot{\mathcal{P}}_{EB}^{(1)}\right]_{b}\,,$
(2.18)

where we have defined

$\displaystyle\mathcal{B}_{f;E}^{(n)}$
$\displaystyle=\frac{1}{a^{n}}\left\langle{f\left(\left(\vec{\nabla}\times\right)^{n}\vec{E}\right)\cdot\vec{E}}\right\rangle\,,\quad\mathcal{B}_{f;B}^{(n)}=\frac{1}{a^{n}}\left\langle{f\left(\left(\vec{\nabla}\times\right)^{n}\vec{B}\right)\cdot\vec{B}}\right\rangle\,,\quad\mathcal{B}_{f;EB}^{(0)}=-\left\langle{f\vec{E}\cdot\vec{B}}\right\rangle\,,$
$\displaystyle\mathcal{B}_{f;EB}^{(1,0)}$
$\displaystyle=-\frac{1}{a}\left\langle{f\left(\vec{\nabla}\times\vec{E}\right)\cdot\vec{B}}\right\rangle\,,\quad\mathcal{B}_{f;EB}^{(0,1)}=-\frac{1}{a}\left\langle{f\vec{E}\cdot\left(\vec{\nabla}\times\vec{B}\right)}\right\rangle\,,$
(2.19)

with $f=\chi,\dot{\chi}$ and $n=\\{0,1\\}$. The two superscripts on the
3-point functions refer to the number of curls acting on $\vec{E}$ and
$\vec{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 most222The scaling will be
milder with suitable optimizing of the equations, using in particular
integration by parts.  as $n_{\cal B}^{3}$ with $n_{\cal 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.333
Note that in the presence of axion gradients, the electric field is no longer
divergence free. Defining $\vec{D}=\vec{E}+(\beta/M_{P})\chi\vec{B}$, we can
rewrite our expressions in terms of a divergence free quantity,
$\vec{\nabla}\cdot\vec{D}=0$. The price to pay is the appearance of $5$-point
functions in the EOMs of the $3$-point functions, whose factorization requires
assumptions not only on the gaussianity of the gauge fields but also of the
axion perturbation $\chi$. While for the former, this is in good agreement
with results found in lattice simulations (we thank Dani Figueroa and Ander
Urio Garmendia for providing valuable input and cross-checks on this point),
the axion fluctuations are expected to be highly non-gaussian due to the non-
linear source terms. Therefore, under the approximations we employ we find it
more convenient to work with original electric field, taking into account that
it is not divergence free.  We moreover set $4$-point functions that involve
one factor of $\vec{\nabla}\chi$ 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

$\displaystyle\dot{\mathcal{P}}^{(0)}_{\chi}-2\mathcal{P}^{(0)}_{\chi\dot{\chi}}=0\,,$
(2.20)
$\displaystyle\dot{\mathcal{P}}^{(0)}_{\chi\dot{\chi}}+3H\mathcal{P}_{\chi\dot{\chi}}+m_{\phi}^{2}\mathcal{P}_{\chi}^{(0)}+\frac{\beta}{M_{P}}\mathcal{B}_{\chi;EB}^{(0)}-\mathcal{P}_{\dot{\chi}}^{(0)}=0\,,$
(2.21)
$\displaystyle\dot{\mathcal{P}}^{(0)}_{\dot{\chi}}+6H\mathcal{P}^{(0)}_{\dot{\chi}}+2m_{\phi}^{2}\mathcal{P}^{(0)}_{\chi\dot{\chi}}+\frac{2\beta}{M_{P}}\mathcal{B}_{\dot{\chi};EB}^{(0)}=0\,,$
(2.22)

with

$\displaystyle\mathcal{P}_{\chi}^{(2n)}=\frac{1}{a^{2n}}\left\langle{\chi\nabla^{2n}\chi}\right\rangle\,,\quad\mathcal{P}_{\dot{\chi}}^{(2n)}=\frac{1}{a^{2n}}\left\langle{\dot{\chi}\nabla^{2n}\dot{\chi}}\right\rangle\,,\quad\mathcal{P}_{\chi\dot{\chi}}^{(2n)}=\frac{1}{a^{2n}}\left\langle{\chi\nabla^{2n}\dot{\chi}}\right\rangle\,.$
(2.23)

###  2.2 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 ${\cal
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 ${\cal 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

$\displaystyle\mathcal{P}_{E}^{(n)}$
$\displaystyle=\frac{1}{a^{n+4}}\int\frac{d^{3}k}{(2\pi)^{3}}\theta(k_{h}(t)-k)\sum_{\sigma}(\sigma
k)^{n}\left|{\frac{dA_{\sigma}}{d\tau}}\right|^{2}\,,$ (2.24)
$\displaystyle\mathcal{P}_{B}^{(n)}$
$\displaystyle=\frac{1}{a^{n+4}}\int\frac{d^{3}k}{(2\pi)^{3}}\theta(k_{h}(t)-k)\sum_{\sigma}(\sigma
k)^{n+2}\left|{A_{\sigma}}\right|^{2}\,,$ (2.25)
$\displaystyle\mathcal{P}_{EB}^{(n)}$
$\displaystyle=\frac{1}{2a^{n+4}}\int\frac{d^{3}k}{(2\pi)^{3}}\theta(k_{h}(t)-k)\sum_{\sigma}(\sigma
k)^{n+1}\frac{d}{d\tau}\left|{A_{\sigma}}\right|^{2}\,,$ (2.26)

where $\sigma$ 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}(\tau)$,

$\displaystyle\frac{d^{2}A_{\sigma}}{d\tau^{2}}+k\left(k-\lambda\sigma 2\xi
aH\right)A_{\sigma}=0\,,\quad\lambda=\text{sign}(\dot{\phi})\,,$ (2.27)

which encounters a tachyonic instability for the polarization $\sigma=\lambda$
for

$\displaystyle k_{h}(t)=\underset{t^{\prime}\leq
t}{\mathrm{max}}\left[2\xi(t^{\prime})a(t^{\prime})H(t^{\prime})\right]\quad\text{with
}\xi=\beta|\dot{\phi}|/(2HM_{P})\,.$ (2.28)

Taking this into account leads to boundary terms in the EOMs for the $2$-point
functions, which are explicitly given as

$\displaystyle\left[\dot{\mathcal{P}}_{E}^{(n)}\right]_{b}$
$\displaystyle=\frac{\dot{k}_{h}}{a^{n+4}}\frac{k_{h}^{2}}{2\pi^{2}}\sum_{\sigma}(\sigma
k_{h})^{n}\left|{\frac{dA_{\sigma}}{d\tau}}\right|_{k=k_{h}}^{2}\,,$ (2.29)
$\displaystyle\left[\dot{\mathcal{P}}_{B}^{(n)}\right]_{b}$
$\displaystyle=\frac{\dot{k}_{h}}{a^{n+4}}\frac{k_{h}^{2}}{2\pi^{2}}\sum_{\sigma}(\sigma
k_{h})^{n+2}\left|{A_{\sigma}}\right|_{k=k_{h}}^{2}\,,$ (2.30)
$\displaystyle\left[\dot{\mathcal{P}}_{EB}^{(n)}\right]_{b}$
$\displaystyle=\frac{\dot{k}_{h}}{a^{n+4}}\frac{k_{h}^{2}}{2\pi^{2}}\sum_{\sigma}(\sigma
k_{h})^{n+1}\mathrm{Re}\left[A_{\sigma}^{*}\frac{dA_{\sigma}}{d\tau}\right]_{k=k_{h}}\,.$
(2.31)

To evaluate these, we follow the prescription given in [30] which is based on
the solutions for the mode functions for constant $\xi$. 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 $\xi$.

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

$\displaystyle\mathcal{P}_{X}^{(n_{\mathrm{max}}+1)}$
$\displaystyle\simeq\left(\frac{k_{h}}{a}\right)^{2}\mathcal{P}_{X}^{(n_{\mathrm{max}}-1)}\,,$
(2.32)

where $X=E,B$ or $EB$. This can be derived from Eqs (2.29) to (2.31) assuming
a power law spectrum for $|A_{\sigma}(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
$\mathcal{P}_{X}^{(n+1)}/\mathcal{P}_{X}^{(n-1)}$, 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 $\xi$ is constant. As we show there,
a value of $n\sim 55$ is necessary for the relation (2.32) to hold up to about
$5\%$. To our understanding, this explains why values of $n\sim 100$ are
necessary to achieve convergence in the GEF, whereas $n\gg 1$ should have been
sufficient based on the assumption of power law spectrum around $k_{h}$.

This in turn indicates that rapid changes in $\xi$ 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 $\mathcal{P}_{X}^{(n)}$ 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
$\mathcal{P}_{X}^{(0)}$ 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
$\mathcal{P}_{X}^{(n_{\text{max}}-1)}$ but instead a series of
$\mathcal{P}_{X}^{(n_{\text{max}}+1-2l)}$ to determine
$\mathcal{P}_{X}^{(n_{\text{max}}+1)}$,

$\displaystyle\bar{\mathcal{P}}_{X}^{(n_{\mathrm{max}}+1)}=\sum_{l=1}^{L}(-1)^{l-1}\begin{pmatrix}L\\\
l\end{pmatrix}\bar{\mathcal{P}}_{X}^{(n_{\mathrm{max}}+1-2l)}\,,$ (2.33)

where
$\bar{\mathcal{P}}_{X}^{(n)}=\mathcal{P}_{X}^{(n)}/H_{0}^{4}(k_{h}/a)^{n}$.
For $L=1$ we trivially recover Eq. (2.32). As we show in 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.

###  2.3 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,

$\displaystyle
R_{\chi}\equiv\left|\frac{\left\langle{(\nabla\chi)^{2}}\right\rangle}{\dot{\phi}^{2}+\left\langle{\dot{\chi}^{2}}\right\rangle}\right|\,.$
(2.34)

As long as this quantity is small, it is justified to drop higher powers of
the axion gradients, whereas $R_{\chi}$ above ${\cal O}(0.1\,\mathchar
45\relax\mathchar 45\relax\,1)$ indicates that the axion gradients can no
longer be treated perturbatively, calling for a full lattice simulation. To
compute
$\left\langle{(\nabla\chi)^{2}}\right\rangle=-\left\langle{\chi\nabla^{2}\chi}\right\rangle$,
the relevant EOMs for the axion $2$-point functions to second order in the
gradient expansion are

$\displaystyle\dot{\mathcal{P}}_{\chi}^{(2)}+2H\mathcal{P}_{\chi}^{(2)}-2\mathcal{P}_{\chi\dot{\chi}}^{(2)}=0\,,$
(2.35)
$\displaystyle\dot{\mathcal{P}}_{\chi\dot{\chi}}^{(2)}+5H\mathcal{P}_{\chi\dot{\chi}}^{(2)}+m_{\phi}^{2}\mathcal{P}_{\chi}^{(2)}+\frac{\beta}{M_{P}}\mathcal{B}_{\chi;EB}^{(2;0,0)}-\mathcal{P}_{\dot{\chi}}^{(2)}=0\,,$
(2.36)
$\displaystyle\dot{\mathcal{P}}_{\dot{\chi}}^{(2)}+8H\mathcal{P}_{\dot{\chi}}^{(2)}+2m_{\phi}^{2}\mathcal{P}_{\chi\dot{\chi}}^{(2)}+\frac{2\beta}{M_{P}}\mathcal{B}_{\dot{\chi};EB}^{(2;0,0)}=0\,,$
(2.37)

which in turn require the evaluation of the 3-point function
$\mathcal{B}_{f;EB}^{(2;0,0)}$, whose definitions and equations are again
given in the App. A. The breakdown of this perturbative expansion scheme
$R_{\chi}>0.5$ is indicated by the dark gray region in the figures below.

##  3 Numerical results

Figure 1: Evolution of $\xi$ assuming a homogeneous axion field (dashed black)
and perturbatively including axion gradients (red) for different values of
$\beta$. 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].

Figs. 1 and 2 show the results of the methodology described in Sec. 2 for
values of the coupling $\beta$ ranging from 15 to 25. The six panels of Fig. 1
focus on the evolution of the parameter $\xi=\beta|\dot{\phi}|/(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_{\chi}>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 [35]. The light gray region indicates the regime
in which the gradient energy is sizeable, i.e. $0.01<R_{\chi}<0.5$, but our
perturbative treatment is still valid. The vertical gray line indicates
$R_{\chi}=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_{\chi}\ll 1$ a systematic expansion to higher orders in
the axion gradients might conceivably achieve this (within the light gray
region).

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

Fig. 2 shows the evolution of the kinetic ($\rho_{\mathrm{K}}$), potential
($\rho_{\mathrm{V}}$), electromagnetic ($\rho_{\mathrm{EM}}$) and axion
gradient energy ($\rho_{\mathrm{G}}$) defined as

$\displaystyle\rho_{\mathrm{K}}=\frac{1}{2}\left(\dot{\phi}^{2}+\mathcal{P}_{\dot{\chi}}^{(0)}\right),\quad\rho_{\mathrm{V}}=\frac{m_{\phi}^{2}}{2}\left(\phi^{2}+\mathcal{P}_{\chi}^{(0)}\right),\quad\rho_{\mathrm{EM}}=\frac{1}{2}\left(\mathcal{P}_{E}^{(0)}+\mathcal{P}_{B}^{(0)}\right),\quad\rho_{\mathrm{G}}=-\frac{1}{2}\mathcal{P}_{\chi}^{(2)},$
(3.1)

for the same values of $\beta$. Here we include $\rho_{\mathrm{G}}$ in the
total energy density,
$\rho_{\mathrm{tot}}=\rho_{\mathrm{K}}+\rho_{\mathrm{V}}+\rho_{\mathrm{EM}}+\rho_{\mathrm{G}}$,
although we do not include $\rho_{\mathrm{G}}$ in the Friedmann equation as we
explained above. This figure gives a more detailed view of the importance of
the axion gradient energy, showing that in many cases (e.g. $\beta=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 $\beta=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 $\beta$, 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 $\beta$ in a non-monotonic way.

|
---|---
Figure 3: Power spectrum of axion perturbations,
$\Delta_{\zeta}^{2}=(H/\dot{\phi})^{2}\mathcal{P}_{\chi}^{(0)}$, for
$\beta=15$ (left) and $\beta=20$ (right). Large values imply PBH formation,
see text for details. Gray bands as in Fig. 1.

As one example of possible observable consequences we show the scalar power
spectrum sourced by the axion fluctuations,
$\Delta_{\zeta}^{2}=(H/\dot{\phi})^{2}\mathcal{P}_{\chi}^{(0)}$ in Fig. 3.444
Strictly speaking $\mathcal{P}_{\chi}^{(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\sim k_{h}$ dominate the integral and the time evolution
of $\mathcal{P}_{X}^{(0)}$ roughly corresponds to the spectrum.  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 $\Delta_{\zeta}^{2}\sim
10^{-2}$ for a Gaussian distribution and around $\Delta_{\zeta}^{2}\sim
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 see from Fig. 3, even for a moderate
coupling of $\beta=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 $\beta$ (which, in particular, do
not produce excessive gravitational waves during preheating [41, 42]) our
results indicate a phase of PBH formation. Increasing the coupling $\beta$
will typically extend this phase, leading to more massive PBHs with longer
lifetimes. Interestingly, for some parameter choices such as $\beta=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.

##  4 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 $n_{3}=n_{\cal B}=1$.

This captures several important aspects of the dynamics of the system. For
very small values of the expansion parameter $R_{\chi}$ 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_{\chi}\gtrsim 1\%$, the homogeneous GEF results and the
lattice results start to diverge, with our perturbative expansion recovering
the lattice results. For $R_{\chi}\gtrsim 5\%$ the leading order corrections
included here no longer suffice to accurately track the system, though given
the smallness of $R_{\chi}$ a systematic inclusion of higher order terms might
conceivably achieve this. We see no fundamental obstruction to this in this
regime. Finally, for $R_{\chi}\gtrsim 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.

#### Acknowledgements

We thank Kai Schmitz, Oleksandr Sobol and Richard Voneckardstein for fruitful
discussions and the very helpful comparisons of our numerical codes, as well
as Dani Figueroa and Ander Urio Garmendia for insightful cross-checks from
lattice simulations. We moreover thank Ryo Namba for helpful discussions, and
Kyohei Mukaida for collaboration at the initial stages of this project. Y.E.
and S.S. would like to thank the CERN Theory Department for their hospitality
during crucial parts of this project. The work of Y.E. is supported in part by
DOE grant DE-SC0011842. The work of S.S. received the support of a fellowship
from “la Caixa” Foundation (ID 100010434) with fellowship code
LCF/BQ/DI19/11730034 and by the Generalitat Valenciana grant
PROMETEO/2021/083.

##  Appendix A Equations of motion for $3$-point functions

The main goal of this paper is to extend the GEF and include the axion
inhomogeneity $\chi(t,\vec{x})$. Once we include $\chi$, the EOMs become non-
linear in terms of the inhomogeneous quantities, $\chi$, $\vec{E}$ and
$\vec{B}$. As a result, the EOMs of the electromagnetic $2$-point functions no
longer form a closed system and we need to include the evolution of the
$3$-point functions. The EOM of the $3$-point functions then depends on
$4$-point functions and a similar structure persists for higher-point
functions. To truncate this tower of $p$-point functions, we factorize the
$4$-point functions as products of the $2$-point functions. Moreover, for
simplicity, we include the $3$-point functions with only up to one spatial
derivative in the full system. The latter approximation is expected to be
valid when the axion gradient energy is suppressed, and to check the quality
of our approximation, we compute the axion gradient energy by treating it as a
perturbation and monitor its size.

In the following, we list the time evolution equations of the $2$-point and
$3$-point functions that we need in our numerical computation. These are
derived by the repeated use of Eqs. (2.4)–(2.6). In this appendix we use the
following notation for 3-point functions in the GEF,

$\displaystyle\mathcal{B}_{f;X}^{(2l;n,m)}$
$\displaystyle=\frac{1}{a^{n+m+2l}}\left\langle{\nabla^{2l}f(\vec{\nabla}\times)^{n}\vec{X}\cdot(\vec{\nabla}\times)^{m}\vec{X}}\right\rangle\,,\quad\mathcal{B}_{f;XY}^{(2l;n,m)}=-\frac{1}{a^{n+m+2l}}\left\langle{\nabla^{2l}f(\vec{\nabla}\times)^{n}\vec{X}\cdot(\vec{\nabla}\times)^{m}\vec{Y}}\right\rangle\,,$
(A.1)

with $X,Y=\\{E,B\\}$ and $f=\\{\chi,\dot{\chi}\\}$. These are related to the
simplified notation used in the main text (which does not include axion
gradients) as

$\displaystyle\mathcal{B}_{f;E}^{(n)}=\mathcal{B}_{f;E}^{(0;n,0)}\,,\quad\mathcal{B}_{f;B}^{(n)}=\mathcal{B}_{f;B}^{(0;n,0)}\,,\quad\mathcal{B}_{f;EB}^{(0)}=\mathcal{B}_{f;EB}^{(0;0,0)}\,,\quad\mathcal{B}_{f;EB}^{(1,0)}=\mathcal{B}_{f;EB}^{(0;1,0)}\,,\quad\mathcal{B}_{f;EB}^{(0,1)}=\mathcal{B}_{f;EB}^{(0;0,1)}\,.$
(A.2)

Notice that

$\displaystyle\mathcal{B}_{f;E}^{(2l;n,m)}=\mathcal{B}_{f;E}^{(2l;m,n)}\,,\quad\mathcal{B}_{f;B}^{(2l;n,m)}=\mathcal{B}_{f;B}^{(2l;m,n)}\,,$
(A.3)

by definition, but in general

$\displaystyle\mathcal{B}_{f;EB}^{(2l;n,m)}\neq\mathcal{B}_{f;EB}^{(2l;m,n)}\,,$
(A.4)

for $n\neq m$.

###  A.1 Equations of motion up to one spatial derivative

As we stated above, we include the $3$-point functions with only up to one
spatial derivative. With this in mind, the electromagnetic $2$-point functions
evolve as

$\displaystyle\dot{\mathcal{P}}_{E}^{(n)}+(n+4)H\mathcal{P}_{E}^{(n)}+2\mathcal{P}_{EB}^{(n+1)}-\frac{2\beta}{M_{P}}\left(\dot{\phi}\mathcal{P}_{EB}^{(n)}+\mathcal{B}_{\dot{\chi};EB}^{(0;n,0)}\right)+\frac{2\beta}{M_{P}}\left(\mathcal{B}_{\chi;E}^{(0;n+1,0)}-\mathcal{B}_{\chi;E}^{(0;n,1)}\right)=\left[\dot{\mathcal{P}}_{E}^{(n)}\right]_{b}\,,$
$\displaystyle\dot{\mathcal{P}}_{B}^{(n)}+(n+4)H\mathcal{P}_{B}^{(n)}-2\mathcal{P}_{EB}^{(n+1)}=\left[\dot{\mathcal{P}}_{B}^{(n)}\right]_{b}\,,$
(A.5)
$\displaystyle\dot{\mathcal{P}}_{EB}^{(n)}+(n+4)H\mathcal{P}_{EB}^{(n)}-\mathcal{P}_{E}^{(n+1)}+\mathcal{P}_{B}^{(n+1)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{P}_{B}^{(n)}+\mathcal{B}_{\dot{\chi};B}^{(0;n,0)}\right)-\frac{\beta}{M_{P}}\left(\mathcal{B}_{\chi;EB}^{(0;1,n)}-\mathcal{B}_{\chi;EB}^{(0;0,n+1)}\right)=\left[\dot{\mathcal{P}}_{EB}^{(n)}\right]_{b}\,,$

where it is understood that

$\displaystyle\mathcal{B}_{f;XY}^{(2l;n,m)}=0,~{}~{}\mathrm{for}~{}~{}n+m+2l>1\,,$
(A.6)

in these equations. The axion $2$-point functions with no spatial derivatives
are given by

$\displaystyle\dot{\mathcal{P}}_{\chi}^{(0)}-2\mathcal{P}_{\chi\dot{\chi}}^{(0)}=0\,,$
$\displaystyle\dot{\mathcal{P}}_{\chi\dot{\chi}}^{(0)}+3H\mathcal{P}_{\chi\dot{\chi}}^{(0)}+m_{\phi}^{2}\mathcal{P}_{\chi}^{(0)}+\frac{\beta}{M_{P}}\mathcal{B}_{\chi;EB}^{(0;0,0)}-\mathcal{P}_{\dot{\chi}}^{(0)}=0\,,$
$\displaystyle\dot{\mathcal{P}}_{\dot{\chi}}^{(0)}+6H\mathcal{P}_{\dot{\chi}}^{(0)}+2m_{\phi}^{2}\mathcal{P}_{\chi\dot{\chi}}^{(0)}+\frac{2\beta}{M_{P}}\mathcal{B}_{\dot{\chi};EB}^{(0;0,0)}=0\,.$
(A.7)

The time evolution equations of the $3$-point functions are given as follows.
For $(2l;n,m)=(0;0,0)$, we obtain

$\displaystyle\dot{\mathcal{B}}_{\chi;E}^{(0;0,0)}+4H\mathcal{B}_{\chi;E}^{(0;0,0)}+2\mathcal{B}_{\chi;EB}^{(0;0,1)}-\frac{2\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\chi;EB}^{(0;0,0)}+\mathcal{P}_{\chi\dot{\chi}}^{(0)}\mathcal{P}_{EB}^{(0)}\right)-\mathcal{B}_{\dot{\chi};E}^{(0;0,0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\chi;B}^{(0;0,0)}+4H\mathcal{B}_{\chi;B}^{(0;0,0)}-2\mathcal{B}_{\chi;EB}^{(0;1,0)}-\mathcal{B}_{\dot{\chi};B}^{(0;0,0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\chi;EB}^{(0;0,0)}+4H\mathcal{B}_{\chi;EB}^{(0;0,0)}-\mathcal{B}_{\chi;E}^{(0;1,0)}+\mathcal{B}_{\chi;B}^{(0;1,0)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\chi;B}^{(0;0,0)}+\mathcal{P}_{\chi\dot{\chi}}^{(0)}\mathcal{P}_{B}^{(0)}\right)-\mathcal{B}_{\dot{\chi};EB}^{(0;0,0)}=0\,,$
(A.8)

and

$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};E}^{(0;0,0)}+7H\mathcal{B}_{\dot{\chi};E}^{(0;0,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;E}^{(0;0,0)}+2\mathcal{B}_{\dot{\chi};EB}^{(0;0,1)}-\frac{2\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\dot{\chi};EB}^{(0;0,0)}+\mathcal{P}_{\dot{\chi}}^{(0)}\mathcal{P}_{EB}^{(0)}\right)+\frac{2\beta}{3M_{P}}\mathcal{P}_{EB}^{(0)}\mathcal{P}_{E}^{(0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};B}^{(0;0,0)}+7H\mathcal{B}_{\dot{\chi};B}^{(0;0,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;B}^{(0;0,0)}-2\mathcal{B}_{\dot{\chi};EB}^{(0;1,0)}+\frac{2\beta}{3M_{P}}\mathcal{P}_{EB}^{(0)}\mathcal{P}_{B}^{(0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};EB}^{(0;0,0)}+7H\mathcal{B}_{\dot{\chi};EB}^{(0;0,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;EB}^{(0;0,0)}-\mathcal{B}_{\dot{\chi};E}^{(0;1,0)}+\mathcal{B}_{\dot{\chi};B}^{(0;1,0)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\dot{\chi};B}^{(0;0,0)}+\mathcal{P}_{\dot{\chi}}^{(0)}\mathcal{P}_{B}^{(0)}\right)+\frac{\beta}{3M_{P}}\left[\left(\mathcal{P}_{EB}^{(0)}\right)^{2}+\mathcal{P}_{E}^{(0)}\mathcal{P}_{B}^{(0)}\right]=0\,.$
(A.9)

Notice, in particular, that the equations contain the products of the
electromagnetic $2$-point functions. This is a result of the factorization of
the $4$-point functions. Assuming the isotropy and Gaussianity of the
electromagnetic fields, we can factorize e.g.

$\displaystyle\left\langle{(\vec{E}\cdot\vec{B})^{2}}\right\rangle\simeq\frac{4}{3}\left(\mathcal{P}_{EB}^{(0)}\right)^{2}+\frac{1}{3}\mathcal{P}_{E}^{(0)}\mathcal{P}_{B}^{(0)}\,,$
(A.10)

and so on. For $(2l;n,m)=(0;1,0),(0;0,1)$, we obtain

$\displaystyle\dot{\mathcal{B}}_{\chi;E}^{(0;1,0)}+5H\mathcal{B}_{\chi;E}^{(0;1,0)}-\frac{\beta\dot{\phi}}{M_{P}}\left(\mathcal{B}_{\chi;EB}^{(0;1,0)}+\mathcal{B}_{\chi;EB}^{(0;0,1)}\right)-\frac{2\beta}{M_{P}}\mathcal{P}_{\chi\dot{\chi}}^{(0)}\mathcal{P}_{EB}^{(1)}-\mathcal{B}_{\dot{\chi};E}^{(0;1,0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\chi;B}^{(0;1,0)}+5H\mathcal{B}_{\chi;B}^{(0;1,0)}-\mathcal{B}_{\dot{\chi};B}^{(0;1,0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\chi;EB}^{(0;1,0)}+5H\mathcal{B}_{\chi;EB}^{(0;1,0)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\chi;B}^{(0;1,0)}+\mathcal{P}_{\chi\dot{\chi}}^{(0)}\mathcal{P}_{B}^{(1)}\right)-\mathcal{B}_{\dot{\chi};EB}^{(0;1,0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\chi;EB}^{(0;0,1)}+5H\mathcal{B}_{\chi;EB}^{(0;0,1)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\chi;B}^{(0;1,0)}+\mathcal{P}_{\chi\dot{\chi}}^{(0)}\mathcal{P}_{B}^{(1)}\right)-\mathcal{B}_{\dot{\chi};EB}^{(0;0,1)}=0\,,$
(A.11)

and

$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};E}^{(0;1,0)}+8H\mathcal{B}_{\dot{\chi};E}^{(0;1,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;E}^{(0;1,0)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{\beta}{M_{P}}\left(\dot{\phi}\left(\mathcal{B}_{\dot{\chi};EB}^{(0;1,0)}+\mathcal{B}_{\dot{\chi};EB}^{(0;0,1)}\right)+2\mathcal{P}_{\dot{\chi}}^{(0)}\mathcal{P}_{EB}^{(1)}\right)+\frac{\beta}{3M_{P}}\left(\mathcal{P}_{EB}^{(0)}\mathcal{P}_{E}^{(1)}+\mathcal{P}_{E}^{(0)}\mathcal{P}_{EB}^{(1)}\right)=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};B}^{(0;1,0)}+8H\mathcal{B}_{\dot{\chi};B}^{(0;1,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;B}^{(0;1,0)}+\frac{\beta}{3M_{P}}\left(\mathcal{P}_{B}^{(0)}\mathcal{P}_{EB}^{(1)}+\mathcal{P}_{EB}^{(0)}\mathcal{P}_{B}^{(1)}\right)=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};EB}^{(0;1,0)}+8H\mathcal{B}_{\dot{\chi};EB}^{(0;1,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;EB}^{(0;1,0)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\dot{\chi};B}^{(0;1,0)}+\mathcal{P}_{\dot{\chi}}^{(0)}\mathcal{P}_{B}^{(1)}\right)+\frac{\beta}{3M_{P}}\left(\mathcal{P}_{B}^{(0)}\mathcal{P}_{E}^{(1)}+\mathcal{P}_{EB}^{(0)}\mathcal{P}_{EB}^{(1)}\right)=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};EB}^{(0;0,1)}+8H\mathcal{B}_{\dot{\chi};EB}^{(0;0,1)}+m_{\phi}^{2}\mathcal{B}_{\chi;EB}^{(0;0,1)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\dot{\chi};B}^{(0;1,0)}+\mathcal{P}_{\dot{\chi}}^{(0)}\mathcal{P}_{B}^{(1)}\right)+\frac{\beta}{3M_{P}}\left(\mathcal{P}_{E}^{(0)}\mathcal{P}_{B}^{(1)}+\mathcal{P}_{EB}^{(0)}\mathcal{P}_{EB}^{(1)}\right)=0\,.$
(A.12)

Finally, the background EOM is given by

$\displaystyle
0=\ddot{\phi}+3H\dot{\phi}+m_{\phi}^{2}\phi+\frac{\beta}{M_{P}}\mathcal{P}_{EB}^{(0)}\,,$
$\displaystyle
H^{2}=\frac{1}{6M_{P}^{2}}\left[\dot{\phi}^{2}+\mathcal{P}_{\dot{\chi}}^{(0)}+m_{\phi}^{2}\left(\phi^{2}+\mathcal{P}_{\chi}^{(0)}\right)+\mathcal{P}_{E}^{(0)}+\mathcal{P}_{B}^{(0)}\right]\,,$
$\displaystyle\dot{H}=-\frac{1}{6M_{P}^{2}}\left[3\dot{\phi}^{2}+3\mathcal{P}_{\dot{\chi}}^{(0)}+2\left(\mathcal{P}_{E}^{(0)}+\mathcal{P}_{B}^{(0)}\right)\right]\,.$
(A.13)

These equations are consistent, even after our approximations, in the sense
that we can derive the last equation from the others.

###  A.2 Equations of motion with two spatial derivatives

To monitor the size of the axion inhomogeneity, we follow the time evolution
of the axion gradient energy, i.e., $\mathcal{P}_{\chi}^{(2)}$. For this
purpose, we need to follow some of the $3$-point functions with two spatial
derivatives. We treat these quantities as perturbations, by treating the
quantities with only up to one spatial derivative as the source term and
ignoring the backreaction of the axion gradient energy. The relevant equations
are given by

$\displaystyle\dot{\mathcal{P}}_{\chi}^{(2)}+2H\mathcal{P}_{\chi}^{(2)}-2\mathcal{P}_{\chi\dot{\chi}}^{(2)}=0\,,$
$\displaystyle\dot{\mathcal{P}}_{\chi\dot{\chi}}^{(2)}+5H\mathcal{P}_{\chi\dot{\chi}}^{(2)}+m_{\phi}^{2}\mathcal{P}_{\chi}^{(2)}+\frac{\beta}{M_{P}}\mathcal{B}_{\chi;EB}^{(2;0,0)}-\mathcal{P}_{\dot{\chi}}^{(2)}=0\,,$
$\displaystyle\dot{\mathcal{P}}_{\dot{\chi}}^{(2)}+8H\mathcal{P}_{\dot{\chi}}^{(2)}+2m_{\phi}^{2}\mathcal{P}_{\chi\dot{\chi}}^{(2)}+\frac{2\beta}{M_{P}}\mathcal{B}_{\dot{\chi};EB}^{(2;0,0)}=0\,.$
(A.14)

Obviously we need to compute $\mathcal{B}_{f;EB}^{(2;0,0)}$. For this, it is
enough to consider only the following equations:

$\displaystyle\dot{\mathcal{B}}_{\chi;B}^{(2;0,0)}+6H\mathcal{B}_{\chi;B}^{(2;0,0)}-\mathcal{B}_{\dot{\chi};B}^{(2;0,0)}=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\chi;EB}^{(2;0,0)}+6H\mathcal{B}_{\chi;EB}^{(2;0,0)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\chi;B}^{(2;0,0)}+\mathcal{P}_{\chi\dot{\chi}}^{(2)}\mathcal{P}_{B}^{(0)}\right)-\mathcal{B}_{\dot{\chi};EB}^{(2;0,0)}=0\,,$
(A.15)

and

$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};B}^{(2;0,0)}+9H\mathcal{B}_{\dot{\chi};B}^{(2;0,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;B}^{(2;0,0)}-\frac{2\beta}{3M_{P}}\left(\mathcal{P}_{EB}^{(0)}\mathcal{P}_{B}^{(2)}+\mathcal{P}_{B}^{(0)}\mathcal{P}_{EB}^{(2)}-\mathcal{P}_{EB}^{(1)}\mathcal{P}_{B}^{(1)}\right)=0\,,$
$\displaystyle\dot{\mathcal{B}}_{\dot{\chi};EB}^{(2;0,0)}+9H\mathcal{B}_{\dot{\chi};EB}^{(2;0,0)}+m_{\phi}^{2}\mathcal{B}_{\chi;EB}^{(2;0,0)}-\frac{\beta}{M_{P}}\left(\dot{\phi}\mathcal{B}_{\dot{\chi};B}^{(2;0,0)}+\mathcal{P}_{\dot{\chi}}^{(2)}\mathcal{P}_{B}^{(0)}\right)$
(A.16)
$\displaystyle~{}~{}~{}~{}~{}-\frac{\beta}{3M_{P}}\left[\mathcal{P}_{B}^{(0)}\mathcal{P}_{E}^{(2)}+\mathcal{P}_{E}^{(0)}\mathcal{P}_{B}^{(2)}+2\mathcal{P}_{EB}^{(0)}\mathcal{P}_{EB}^{(2)}-\mathcal{P}_{B}^{(1)}\mathcal{P}_{E}^{(1)}-\left(\mathcal{P}_{EB}^{(1)}\right)^{2}-\frac{\beta^{2}}{3M_{P}^{2}}\mathcal{P}_{\chi}^{(2)}\left(\mathcal{P}_{B}^{(0)}\right)^{2}\right]=0\,,$

where we used the factorization

$\displaystyle\left\langle{\partial_{i}(\vec{E}\cdot\vec{B})\partial_{i}(\vec{X}\cdot\vec{Y})}\right\rangle$
$\displaystyle\simeq\frac{1}{3}\left(\left\langle{\vec{E}^{(1)}\cdot\vec{X}^{(1)}}\right\rangle+\left\langle{(\vec{\nabla}\cdot\vec{E})(\vec{\nabla}\cdot\vec{X})}\right\rangle\right)\left\langle{\vec{B}\cdot\vec{Y}}\right\rangle-\frac{1}{6}\left\langle{\vec{E}^{(1)}\cdot\vec{Y}}\right\rangle\left\langle{\vec{B}\cdot\vec{X}^{(1)}}\right\rangle$
$\displaystyle+(X\leftrightarrow Y)+(E\leftrightarrow B)+(X\leftrightarrow
Y,E\leftrightarrow B)\,.$ (A.17)

Here we use the short-hand notation
$\vec{X}^{(n)}=(\vec{\nabla}\times)^{n}\vec{X}$. These are closed by
themselves, and hence we do not need to compute e.g.
$\mathcal{B}_{f;E}^{(2;0,0)}$ nor $\mathcal{B}_{f;EB}^{(0;2,0)}$.

###  A.3 Numerical implementation

The system of the equations described above has to be solved numerically. It
is then convenient to apply the following variable re-definition

$\displaystyle\begin{split}&\bar{\phi}=\frac{\phi}{M_{p}}\,,\,\,\bar{\mathcal{P}}_{X}^{(n)}=\frac{\mathcal{P}_{X}^{(n)}}{H_{0}^{4}\left(k_{h}/a\right)^{n}}\,,\,\,\bar{t}=H_{0}t\,,\,\,\bar{H}=\frac{H}{H_{0}}\,,\,\,\bar{m}=\frac{m}{H_{0}}\,,\,\,\frac{\bar{k}_{h}}{\bar{a}}=\frac{k_{h}}{aH_{0}}\,,\,\,\\\
&\bar{\mathcal{P}}_{\chi}^{(n)}=\frac{\mathcal{P}_{\chi}^{(n)}}{M_{p}^{2}\left(k_{h}/a\right)^{n}}\,,\,\,\bar{\mathcal{P}}_{\dot{\chi}}^{(n)}=\frac{\mathcal{P}_{\dot{\chi}}^{(n)}}{M_{p}^{2}H_{0}^{2}\left(k_{h}/a\right)^{n}}\,,\,\,\bar{\mathcal{P}}_{\chi\dot{\chi}}^{(n)}=\frac{\mathcal{P}_{\chi\dot{\chi}}^{(n)}}{M_{p}^{2}H_{0}\left(k_{h}/a\right)^{n}}\,,\,\,\\\
&\bar{\mathcal{B}}_{\chi;X}^{(2l;n,m)}=\frac{\mathcal{B}_{\chi;X}^{(2l;n,m)}}{M_{p}H_{0}^{4}\left(k_{h}/a\right)^{2l+n+m}}\,,\,\,\bar{\mathcal{B}}_{\dot{\chi};X}^{(2l;n,m)}=\frac{\mathcal{B}_{\dot{\chi};X}^{(2l;n,m)}}{M_{p}H_{0}^{5}\left(k_{h}/a\right)^{2l+n+m}}\,,\end{split}$
(A.18)

where $H_{0}$ is the Hubble parameter at the beginning of the simulation. The
normalization of the power spectra and bispectra by respective powers of
$(k_{h}/a)$ is found to be crucial to numerically evolve the tower of the
$2$-point functions of Eqs. (2.29)–(2.31) to high $n$. Indeed, we empirically
know that $(k_{h}/aH_{0})\sim\mathcal{O}(0.01)$ at the end of inflation, and
hence a different normalization (e.g. by $H_{0}$ instead of $k_{h}/a$) can
cause an extremely tiny numerical value for high enough $n$. Our normalization
ensures that each term of the gradient expansion of a given $p$-point function
is of the same order, resulting in a more numerically stable system.

We use the multistep Adams method of the GNU scientific library implemented in
C++ which integrates the full system in $\mathcal{O}(1\,\mathrm{s})$. For all
plots in this paper we take the initial conditions and parameters following
Ref. [30], i.e.,

$\displaystyle\phi=-15.55\,M_{P}\,,\quad\dot{\phi}=\sqrt{\frac{2}{3}}m_{\phi}M_{P}\,,\quad
a=1\,,\quad\mathcal{P}_{X}^{(n)}=0\,,\quad\mathcal{B}_{f;X}^{(2l;n,m)}=0\,,$
(A.19)

at the beginning of our simulation, with the reduced Planck mass
$M_{p}=2.435\times 10^{18}\,\mathrm{GeV}$. We take the axion mass as

$\displaystyle\frac{m_{\phi}}{M_{P}}=6.16\times 10^{-6}\,,$ (A.20)

following Ref. [35]. The initial Hubble parameter $H_{0}$ is computed as

$\displaystyle
H_{0}=\sqrt{\frac{\dot{\phi}^{2}+m_{\phi}^{2}\phi^{2}}{6M_{P}^{2}}}\,.$ (A.21)

Throughout this paper, we follow the electromagnetic 2-point function up to
$n=n_{\mathrm{max}}=250$, and express those with $n=n_{\mathrm{max}}+1$ by the
truncation relation.

##  Appendix B Truncation relation

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

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 $\xi$. 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={\cal O}(100)$ are required to obtain accurate results in the GEF. It also
illustrates how when $\xi$ 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 $\xi$. In
the slow-roll approximation with constant $\xi$ the gauge field $A_{k}(\tau)$
mode function can be determined analytically by solving

$\displaystyle\left(\partial_{x}^{2}+1\pm\frac{2\xi}{x}\right)\omega_{\mp}=0\,,$
(B.1)

where we defined $\omega_{\pm}=\sqrt{2k}A_{k}^{\pm}(\tau)$ and $x=-k\tau$, see
Eq. (2.27). Matching to the Bunch-Davies vacuum fixes

$\displaystyle
A_{k}^{+}(\tau)=\frac{1}{\sqrt{2k}}e^{\pi\xi/2}W_{-i\xi,1/2}(2ik\tau)\,.$
(B.2)

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\to 0$), which dominates the contribution to
the $n$-point functions, is given by

$\displaystyle
A_{k}^{+}(\tau)\sim\frac{1}{\sqrt{2k}}\frac{e^{\pi\xi/2}}{\Gamma(1+i\xi)}e^{-\xi\sqrt{-k\tau}}\,.$
(B.3)

For the $n$-point function of the $E$-field, see Eq. (2.24), the relevant
quantity is

$\displaystyle\left|\frac{\mathrm{d}A_{k}^{+}(\tau)}{\mathrm{d}\tau}\right|=\left|\frac{1}{\Gamma(1+i\xi)}\right|\frac{1}{2\sqrt{2}}\frac{\xi
e^{\pi\xi/2}}{\sqrt{-\tau}}e^{-\xi\sqrt{-k\tau}}\,.$ (B.4)

The integral of Eq. (2.24) can then be done analytically and evaluates to

$\displaystyle\mathcal{P}_{E}^{(n)}=\frac{4^{-4-n}H^{4+n}}{2\pi^{2}|\Gamma(1+i\xi)|^{2}}e^{\pi\xi}\xi^{-2(2+n)}\left(\Gamma(6+2n)-\Gamma(6+2n,2\xi\sqrt{-k\tau})\right)\,.$
(B.5)

The ratio of $\mathcal{P}_{E}^{(n+2)}/\mathcal{P}_{E}^{(n)}$ is thus

$\displaystyle\frac{\mathcal{P}_{E}^{(n+2)}}{\mathcal{P}_{E}^{(n)}}\frac{a^{2}}{k_{h}^{2}}=\frac{\bar{\mathcal{P}}_{E}^{(n+2)}}{\bar{\mathcal{P}}_{E}^{(n)}}=\frac{1}{(2\xi)^{6}}\frac{\Gamma(10+2n)-\Gamma(10+2n,2\sqrt{2}\xi^{3/2})}{\Gamma(6+2n)-\Gamma(6+2n,2\sqrt{2}\xi^{3/2})}\,,$
(B.6)

where as in the main text, we have introduced the notation
$\bar{\mathcal{P}}_{X}^{(n)}\equiv\mathcal{P}_{X}^{(n)}/H_{0}^{4}(k_{h}/a)^{n}$.
For large $n$, and any $\xi$, the ratio of Eq. (B.6) asymptotically reaches
unity,

$\displaystyle\mathrm{lim}_{n\to\infty}\frac{\mathcal{\bar{P}}_{E}^{(n+2)}}{\mathcal{\bar{P}}_{E}^{(n)}}=1\,,$
(B.7)

with a deviation of $5\,\%$ for $n=45$ ($n=55$) for $\xi=4$ ($\xi=6$), see
Fig. 4. This, on the one hand confirms the truncation relation chosen in [30]
in the constant $\xi$ 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
$\mathcal{\bar{P}}_{B}^{(n)}$ and $\mathcal{\bar{P}}_{EB}^{(n)}$.

The truncation relation remains a good approximation for slowly varying $\xi$,
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
$\mathcal{\bar{P}}_{E}^{(n+1)}/\mathcal{\bar{P}}_{E}^{(n-1)}$ for fixed e-fold
$\mathcal{N}$ and $\beta=15$. The left panel shows an example from the weak
backreaction regime, where $\xi$ 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\gtrsim 50$, although at
$n\sim n_{\mathrm{max}}=250$ the onset of a violation of this relation can be
observed.

|
---|---
Figure 5: Truncation relation (2.32) evaluated numerically for fixed $\beta=15$, showing good agreement for $n\gtrsim 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: GEF tower for rapidly changing $\xi$ 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.

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 which shows ${\cal\bar{P}}_{E}^{(2n)}$ at a slightly later point in time,
${\cal N}=62$ for different values of $n$. For $50\lesssim 2n\lesssim 150$ we
observe a plateau in ${\cal\bar{P}}_{E}^{(2n)}$, 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$. 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 $\xi$ 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 $\mathcal{\bar{P}}_{E}^{(n+1-2l)}$ to determine the
truncation relation for $\mathcal{\bar{P}}_{E}^{(n+1)}$,

$\displaystyle\bar{\mathcal{P}}_{X}^{(n_{\mathrm{max}}+1)}=\sum_{l=1}^{L}(-1)^{l-1}\begin{pmatrix}L\\\
l\end{pmatrix}\bar{\mathcal{P}}_{X}^{(n_{\mathrm{max}}+1-2l)}$

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,
$\mathcal{P}_{E}^{(2n)}$ should always be positive, while the ‘wiggles’
contain negative valued $\mathcal{P}_{E}^{(2n)}$. Similar observations hold
for the other $2$-point functions $\mathcal{\bar{P}}^{(n)}_{B}$ and
$\mathcal{\bar{P}}^{(n)}_{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 Fig. 7 indicate the regions in which (even with the use of Eq. (2.33)
with $L=10$) we obtain
$|\bar{\mathcal{P}}_{X}^{(2n+2)}/\bar{\mathcal{P}}_{X}^{(2n)}-1|>0.1$, i.e. a
significant violation of the truncation relation, for some values of $50\leq
n\leq 75$ with $X=\\{E,B\\}$. As can be seen, this typically happens after a
sharp drop in $\xi$. For a milder evolution of $\xi$, e.g. for $\beta=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 $\beta$
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_{\chi}>0.5$, for all values of $\beta$
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 $\xi$. This sheds some light on the limitations of the
GEF formalism for rapidly varying $\xi$, 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 7:  Evolution of $\xi$ for $beta=15$ (left) and $\beta=22$ (right),
assuming a homogeneous axion field. The gray regions indicate a violation of
the truncation relation above $10\%$, see text for details.

## References

  * [1] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO].
  * [2] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [3] K. Freese, J. A. Frieman, and A. V. Olinto, “Natural inflation with pseudo - Nambu-Goldstone bosons,” Phys. Rev. Lett. 65 (1990) 3233–3236.
  * [4] M. M. Anber and L. Sorbo, “Naturally inflating on steep potentials through electromagnetic dissipation,” Phys. Rev. D 81 (2010) 043534, arXiv:0908.4089 [hep-th].
  * [5] A. Linde, S. Mooij, and E. Pajer, “Gauge field production in supergravity inflation: Local non-Gaussianity and primordial black holes,” Phys. Rev. D 87 no. 10, (2013) 103506, arXiv:1212.1693 [hep-th].
  * [6] E. Bugaev and P. Klimai, “Axion inflation with gauge field production and primordial black holes,” Phys. Rev. D 90 no. 10, (2014) 103501, arXiv:1312.7435 [astro-ph.CO].
  * [7] S.-L. Cheng, W. Lee, and K.-W. Ng, “Numerical study of pseudoscalar inflation with an axion-gauge field coupling,” Phys. Rev. D 93 no. 6, (2016) 063510, arXiv:1508.00251 [astro-ph.CO].
  * [8] J. Garcia-Bellido, M. Peloso, and C. Unal, “Gravitational waves at interferometer scales and primordial black holes in axion inflation,” JCAP 12 (2016) 031, arXiv:1610.03763 [astro-ph.CO].
  * [9] V. Domcke, F. Muia, M. Pieroni, and L. T. Witkowski, “PBH dark matter from axion inflation,” JCAP 07 (2017) 048, arXiv:1704.03464 [astro-ph.CO].
  * [10] J. Garcia-Bellido, M. Peloso, and C. Unal, “Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter,” JCAP 09 (2017) 013, arXiv:1707.02441 [astro-ph.CO].
  * [11] S.-L. Cheng, W. Lee, and K.-W. Ng, “Primordial black holes and associated gravitational waves in axion monodromy inflation,” JCAP 07 (2018) 001, arXiv:1801.09050 [astro-ph.CO].
  * [12] L. Sorbo, “Parity violation in the Cosmic Microwave Background from a pseudoscalar inflaton,” JCAP 06 (2011) 003, arXiv:1101.1525 [astro-ph.CO].
  * [13] J. L. Cook and L. Sorbo, “Particle production during inflation and gravitational waves detectable by ground-based interferometers,” Phys. Rev. D85 (2012) 023534, arXiv:1109.0022 [astro-ph.CO]. [Erratum: Phys. Rev.D86,069901(2012)].
  * [14] N. Barnaby, E. Pajer, and M. Peloso, “Gauge Field Production in Axion Inflation: Consequences for Monodromy, non-Gaussianity in the CMB, and Gravitational Waves at Interferometers,” Phys. Rev. D85 (2012) 023525, arXiv:1110.3327 [astro-ph.CO].
  * [15] N. Barnaby, R. Namba, and M. Peloso, “Phenomenology of a Pseudo-Scalar Inflaton: Naturally Large Nongaussianity,” JCAP 1104 (2011) 009, arXiv:1102.4333 [astro-ph.CO].
  * [16] M. M. Anber and L. Sorbo, “Non-Gaussianities and chiral gravitational waves in natural steep inflation,” Phys. Rev. D85 (2012) 123537, arXiv:1203.5849 [astro-ph.CO].
  * [17] V. Domcke, M. Pieroni, and P. Binétruy, “Primordial gravitational waves for universality classes of pseudoscalar inflation,” JCAP 1606 (2016) 031, arXiv:1603.01287 [astro-ph.CO].
  * [18] J. Garcia-Bellido, A. Papageorgiou, M. Peloso, and L. Sorbo, “A flashing beacon in axion inflation: recurring bursts of gravitational waves in the strong backreaction regime,” arXiv:2303.13425 [astro-ph.CO].
  * [19] W. D. Garretson, G. B. Field, and S. M. Carroll, “Primordial magnetic fields from pseudoGoldstone bosons,” Phys. Rev. D 46 (1992) 5346–5351, arXiv:hep-ph/9209238.
  * [20] M. M. Anber and L. Sorbo, “N-flationary magnetic fields,” JCAP 10 (2006) 018, arXiv:astro-ph/0606534.
  * [21] C. Caprini and L. Sorbo, “Adding helicity to inflationary magnetogenesis,” JCAP 10 (2014) 056, arXiv:1407.2809 [astro-ph.CO].
  * [22] P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis, “Magnetogenesis from axion inflation,” JCAP 10 (2016) 039, arXiv:1606.08474 [astro-ph.CO].
  * [23] D. Jiménez, K. Kamada, K. Schmitz, and X.-J. Xu, “Baryon asymmetry and gravitational waves from pseudoscalar inflation,” JCAP 12 (2017) 011, arXiv:1707.07943 [hep-ph].
  * [24] R. Durrer, O. Sobol, and S. Vilchinskii, “Backreaction from gauge fields produced during inflation,” Phys. Rev. D 108 no. 4, (2023) 043540, arXiv:2303.04583 [gr-qc].
  * [25] M. M. Anber and E. Sabancilar, “Hypermagnetic Fields and Baryon Asymmetry from Pseudoscalar Inflation,” Phys. Rev. D 92 no. 10, (2015) 101501, arXiv:1507.00744 [hep-th].
  * [26] V. Domcke, B. von Harling, E. Morgante, and K. Mukaida, “Baryogenesis from axion inflation,” JCAP 10 (2019) 032, arXiv:1905.13318 [hep-ph].
  * [27] V. Domcke, K. Kamada, K. Mukaida, K. Schmitz, and M. Yamada, “Wash-in leptogenesis after axion inflation,” JHEP 01 (2023) 053, arXiv:2210.06412 [hep-ph].
  * [28] V. Domcke and K. Mukaida, “Gauge Field and Fermion Production during Axion Inflation,” JCAP 11 (2018) 020, arXiv:1806.08769 [hep-ph].
  * [29] V. Domcke, Y. Ema, and K. Mukaida, “Chiral Anomaly, Schwinger Effect, Euler-Heisenberg Lagrangian, and application to axion inflation,” JHEP 02 (2020) 055, arXiv:1910.01205 [hep-ph].
  * [30] E. V. Gorbar, K. Schmitz, O. O. Sobol, and S. I. Vilchinskii, “Gauge-field production during axion inflation in the gradient expansion formalism,” Phys. Rev. D 104 no. 12, (2021) 123504, arXiv:2109.01651 [hep-ph].
  * [31] G. Dall’Agata, S. González-Martín, A. Papageorgiou, and M. Peloso, “Warm dark energy,” JCAP 08 (2020) 032, arXiv:1912.09950 [hep-th].
  * [32] V. Domcke, V. Guidetti, Y. Welling, and A. Westphal, “Resonant backreaction in axion inflation,” JCAP 09 (2020) 009, arXiv:2002.02952 [astro-ph.CO].
  * [33] M. Peloso and L. Sorbo, “Instability in axion inflation with strong backreaction from gauge modes,” JCAP 01 (2023) 038, arXiv:2209.08131 [astro-ph.CO].
  * [34] R. von Eckardstein, M. Peloso, K. Schmitz, O. Sobol, and L. Sorbo, “Axion inflation in the strong-backreaction regime: decay of the Anber-Sorbo solution,” arXiv:2309.04254 [hep-ph].
  * [35] D. G. Figueroa, J. Lizarraga, A. Urio, and J. Urrestilla, “The strong backreaction regime in axion inflation,” arXiv:2303.17436 [astro-ph.CO].
  * [36] D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, “The art of simulating the early Universe – Part I,” JCAP 04 (2021) 035, arXiv:2006.15122 [astro-ph.CO].
  * [37] D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, “CosmoLattice: A modern code for lattice simulations of scalar and gauge field dynamics in an expanding universe,” Comput. Phys. Commun. 283 (2023) 108586, arXiv:2102.01031 [astro-ph.CO].
  * [38] E. V. Gorbar, K. Schmitz, O. O. Sobol, and S. I. Vilchinskii, “Hypermagnetogenesis from axion inflation: Model-independent estimates,” Phys. Rev. D 105 no. 4, (2022) 043530, arXiv:2111.04712 [hep-ph].
  * [39] O. O. Sobol, E. V. Gorbar, and S. I. Vilchinskii, “Backreaction of electromagnetic fields and the Schwinger effect in pseudoscalar inflation magnetogenesis,” Phys. Rev. D 100 no. 6, (2019) 063523, arXiv:1907.10443 [astro-ph.CO].
  * [40] O. O. Sobol, A. V. Lysenko, E. V. Gorbar, and S. I. Vilchinskii, “Gradient expansion formalism for magnetogenesis in the kinetic coupling model,” Phys. Rev. D 102 no. 12, (2020) 123512, arXiv:2010.13587 [astro-ph.CO].
  * [41] P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, “Constraining axion inflation with gravitational waves from preheating,” Phys. Rev. D 101 no. 8, (2020) 083534, arXiv:1909.12842 [astro-ph.CO].
  * [42] P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, “Constraining Axion Inflation with Gravitational Waves across 29 Decades in Frequency,” Phys. Rev. Lett. 124 no. 17, (2020) 171301, arXiv:1909.12843 [astro-ph.CO].