#  Unveiling the nonlinear dynamics of a rolling axion during inflation

Angelo Caravano  Marco Peloso

###### Abstract

A spectator axion-gauge sector, minimally coupled to the inflaton, with the
axion experiencing a momentary stage of fast roll during cosmological
inflation, can generate unique signatures in primordial density fluctuations
and the gravitational wave background. We present the first lattice simulation
of this system using a novel hybrid numerical scheme. This approach solves the
fully nonlinear dynamics of the axion-gauge sector while treating the
gravitational interaction between the axion and inflaton linearly. Initially,
we test the validity of the WKB approximation in the linear regime. We then
investigate strong backreaction dynamics within the axion-gauge sector. Our
findings reveal that backreaction significantly suppresses the growth of the
gauge field and the amplitude of scalar perturbations. The simulation also
allows us to analyze the non-Gaussianity of scalar fluctuations, including
higher-order statistics. We show that, although non-Gaussianity is suppressed
by strong backreaction, it remains higher than in the minimal model where the
axion coincides with the inflaton. Our results highlight the need for
simulations to make robust predictions to test against data from gravitational
wave interferometers and large-scale structure surveys.

##  1 Introduction

Cosmological inflation [1, 2, 3, 4, 5] explains the properties of the universe
at large scales and provides a mechanism for the generation of the primordial
perturbations [6, 7, 8, 9, 10, 11], with excellent agreement with observations
[12, 13, 14]. A large experimental effort is devoted to the search for the
(still elusive) tensor signal from inflation. Under the standard assumption
that this signal is generated by the inflationary expansion, its amplitude
will indicate the energy scale at which inflation took place, which is
probably the most important unknown quantity in this framework. This is no
longer the case if additional sources produce a tensor signal greater than the
one from the expansion. Examples of this from the literature include the
generation of additional gravitational waves (GW) from spectator fields [15,
16, 17], from an effective field theory approach of broken spatial
diffeomorphisms [18], or from particle and string production during inflation
[19, 20, 21, 22]. Such mechanisms are limited by the fact that these sources
also produce scalar perturbations, typically exceeding the generation of
tensor modes. This results in many of these models having a decreased overall
tensor-to-scalar ratio and makes it impossible to observe the produced
gravitational waves (GW) once the limits from scalar production, enforced by
constraints on non-Gaussianity at Cosmic Microwave Background (CMB) scales or
on Primordial Black Hole (PBH) abundance, are taken into account [16, 21, 23].

The situation is better under specific circumstances. For instance, a
decreased ratio of sourced scalar vs. sourced tensor perturbations is produced
if the sourcing fields are coupled only gravitationally to the inflaton and if
they are relativistic (as a nonrelativistic source has a suppressed quadrupole
moment, and thus a reduced GW production [21]). These two features are present
in the model of Refs. [21, 24], in which massless U(1) gauge fields are
sourced by their axial coupling to a spectator axion. This mechanism was
previously studied in the context in which the rolling axion field was the
inflaton, for which, as already mentioned, stronger bounds from the scalar
production are present [25, 26, 27, 28, 29]. 111See also Ref. [30] for a more
recent treatment of this model using the stochastic formalism. In these works,
the axion is assumed to continuously roll all throughout inflation. Ref. [31]
considered instead a spectator axion with a mass of the order of the Hubble
rate, that rolls in a typical cosine potential for only a few $e$-folds
(typically, $\sim 2-5$) of inflation. This results in sourced signals
localized at the scales that left the horizon while the axion was rolling.
This significantly decreases the limits from the CMB non-Gaussianity, so that
this mechanism can produce an observable tensor signal at CMB scales for an
arbitrarily low scale of inflation [31].

From the theory point of view, axions are ubiquitous in string theory, and for
example the recent work [32] studied how this mechanism might be realized in
string compactifications (where in fact a multiplicity of axions might roll at
different stages during inflation, producing a multitude of peaked signals).
From a phenomenological point of view, the mechanism of [31] for the
generation of localized signals has been employed in a number of works that
directly studied it or discussed it in the context of the CMB [33, 34, 35, 13,
36, 37, 38, 39, 40, 41] (including CMB spectral distortions [42]), of PBH dark
matter [43, 44, 45, 46, 47, 48], and of GW generation at smaller scales [49,
50, 51, 52, 53, 54, 55, 56], also in relation to the recent evidence for a
stochastic gravitational wave background (SGWB) from pulsar timing array (PTA)
measurements [57, 58, 59, 60, 61, 62].222These mechanisms can be readily
extended to the SU(2) case. Refs. [63] and [64] studied the non-abelian
counterparts of [21, 24], with gauge field production in a spectator sector
that resembles the inflationary sector of, respectively, chromo-natural
inflation [65] and gauge-flation [66]. The possibility of a localized axion
roll in this context [67] is technically more difficult than in the U(1) case,
and it can be implemented for instance with a varying coupling between the
axion and the gauge field [42]. Notably, in this context, precision astrometry
will allow us to probe the chirality of the SGWB background [68]. This will
allow us to determine whether the PTA signal can be explained by axions during
inflation, as the SGWB generated by axionic interaction is parity-violating.
Moreover, the parity-odd nature of axion-gauge interactions generates a
sizable parity-odd scalar primordial trispectrum [69, 70], which could explain
the recent claim of a parity-odd 4-point correlation function of galaxies from
the BOSS survey [71, 72]. 333A recent study [73] has shown that an alternative
covariance estimation reduces this signal. Consequently, the statistical
significance of this signal is still an open question, as explained for
example in [74].

In light of the phenomenological interest for these models, an extremely
relevant topic of discussion is the backreaction of the produced gauge fields
on the background evolution of the axion, and the validity of the perturbative
computation with which these scalar and tensor signals are computed [75, 76,
77]. The analytical estimates of [76] indicated that backreaction and
perturbativity are under control for $\xi\lesssim 5$, where $\xi$ is the
typical parameter that controls the gauge field amplification, defined in eq.
(2.3) below. This is close to the value required for several phenomenological
applications of this mechanism. Therefore, while the analytic methods of [31]
can provide a reasonable order of magnitude estimate of these signals,
numerical computations appear to be necessary to provide more precise answers.
The present work is the first numerical study of the rolling axion model of
Ref. [31] using nonlinear lattice simulations.

Understanding backreaction is a central focus in the recent literature. In
Ref. [78], Anber and Sorbo (AS) have provided a very simple and minimal
implementation of warm inflation [79]. In their solution, the axion inflaton
has a steady-state speed, for which the friction due to backreacting gauge
fields balances the slope of the axion potential. This friction can be greater
than the usual Hubble friction term, allowing for an inflationary solution in
potentials that would be too steep to sustain inflation in absence of gauge
field production. A number of studies have shown that the AS steady-state
solution is unstable. Most of these works solved, either analytically [80] or
numerically [81, 82, 83, 84, 85, 86, 87, 88] the evolution equations of this
model for a homogeneous axion coupled to gauge fields (namely, under the
working assumption of neglecting the inhmogeneties of the axion). 444See also
Ref. [89] for a recent work where backreaction is computed with a perturbative
scheme. It was understood in [84] (and then confirmed by the analytic kernel
present in the computation of Ref. [80]) that the instability is due the fact
that the amplitude to which each gauge field mode contributes to backreaction
at any time $t$ is sensitive to the values of the axion velocity at earlier
times $t^{\prime}<t$, while that mode was being amplified. This introduces a
memory effect, to which backreaction tries to ‘adapt’ (rather than the
instantaneous balance between the potential slope and the dissipation).
555Instantaneous backreaction can be achieved in the construction of Ref.
[90]. This results in an oscillatory behaviour of the axion velocity about the
AS solution: sudden episodes of faster axion speed result in bursts of gauge
field amplification, and consequent backreaction, that slow down the velocity
of the axion, so that the gauge field production strongly decreases, reducing
backreaction, giving rise to the next episode of fast roll where the gauge
field grows again. This oscillatory behaviour, illustrated in fig. 1, has a
period of $\sim 4-5$ $e$-folds [84, 80], and it could result in correlated
features in the primordial GW spectrum at different scales [91].

Figure 1: Intuitive illustration of the cyclic backreaction process in axion
inflation, explained in section 1.

Recently, lattice simulations have started to improve our knowledge of the
signatures of gauge fields in axion inflation, confirming the validity of
earlier analytic results in the small backreaction regime, and providing new
insight into the complementary regime of strong backreaction. Lattice
simulations are a well-known tool in early-Universe cosmology, and typically
employed to study nonlinear cosmological phenomena such as cosmological phase
transitions and reheating (see, e.g., [92, 93, 94, 95, 96, 97, 98, 99, 100,
101, 102, 103, 104]). More recently, these numerical techniques have been
extended to the inflationary epoch and are emerging as a crucial tool in this
context [105, 106, 107, 108, 109, 110, 111].

Lattice simulations can take backreaction fully into account, relaxing the
assumptions of previous analytic and numerical studies (such as neglecting
fluctuations in the axion field). The first lattice simulation of axion
inflation was presented in Ref. [107]. This work confirmed the instability of
the AS solution and the oscillatory behavior, showing only the presence of the
first oscillatory peak in the axion velocity due to limited UV resolution.
Subsequently, Ref. [110] showed that these oscillations are partially smoothed
out during the final few $e$-folds of inflation. In light of these results,
the validity of a repeated cyclic backreaction, depicted in fig. 1, remains an
open question. Indeed, these studies represent initial steps towards a
comprehensive understanding of backreaction in this system. A central
challenge is the UV sensitivity of the minimal axion inflation setup, where
the axion coincides with the inflaton. In this scenario, the gauge field is
excited across all scales until the end of inflation, necessitating costly
lattice simulations with a large dynamical range. Furthermore, in the regime
near the end of inflation studied in Ref. [110], metric fluctuations, usually
neglected in lattice simulations, might become significant, requiring
numerical relativity techniques to solve the full General Relativity equations
along with nonlinear field equations, as recently explored in [112].666An
alternative approach would be treating gravity linearly, and check afterwards
that metric perturbations remain small. This is precisely the approach that we
use in this work to study the gravitational interaction between the spectator
axion and the inflaton.

Lattice studies have so far focused on the minimal setup where the axion
coincides with the inflaton. As explained above, this scenario requires a
large dynamical range due to the gauge field being enhanced at all scales,
necessitating costly lattice simulations. However, if the axion acts as a
spectator particle, rolling for only a few $e$-folds, the required dynamical
range to study the nonlinear dynamics of this system is significantly reduced.
The fast-rolling model thus presents an ideal environment to investigate
backreaction in axion inflation and its phenomenological implications. In this
work, we present the first lattice study of a spectator axion-gauge system,
coupled to the inflaton through gravity. To achieve this, we introduce a novel
numerical scheme that solves the spectator axion-gauge sector in a fully
nonlinear fashion (i.e. as in Refs. [107, 110]), while simultaneously
addressing the gravitational interaction between the spectator sector and the
axion in a linear way. 777A similar approach was followed in Refs. [113, 114,
115] that computed GW production during gauge preheating by adding tensor
modes linearly to a lattice computation.

This paper is organized as follows. In section 2, we provide a brief review of
the perturbative (analytical) treatment of this model. In section 3, we
introduce our hybrid nonlinear scheme used to solve this system on the
lattice. In section 4 and section 5, we present results from lattice
simulations in regimes of small and strong backreaction, respectively.
Finally, in section 6, we discuss our results and outline future directions.
The work is concluded by appendix A, in which we present the analytic
computation of the axion perturbations (providing results not explictly given
in [31]) and by appendix B, in which we show some tests of our lattice
computations.

##  2 Review of the linear theory

###  2.1 The model

The model [31] is characterized by an inflaton field $\phi$ and an axion
$\sigma$ that are coupled to each other only gravitationally. The axion has a
typical cosine potential, and a standard pseudo-scalar interaction with a U(1)
field, so that the lagrangian is

${\cal
L}=-\frac{1}{2}\left(\partial\phi\right)^{2}-V\left(\phi\right)-\frac{1}{2}\left(\partial\sigma\right)^{2}-\frac{\Lambda^{4}}{2}\left[\cos\left(\frac{\sigma}{f}\right)+1\right]-\frac{1}{4}F^{2}-\frac{\alpha}{4f}\sigma
F{\tilde{F}},$ (2.1)

where $F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the
field strength of the U(1) field, and
${\tilde{F}}^{\mu\nu}\equiv\frac{\epsilon^{\mu\nu\rho\sigma}}{2\sqrt{-g}}F_{\rho\sigma}$
is its dual. 888The quantity $g$ is the determinant of the metric, while
$\epsilon^{\mu\nu\rho\sigma}$ is totally anti-symmetric, with
$\epsilon^{0123}=1$. We assume a FLRW geometry, with line element
$ds^{2}=a^{2}\left(\tau\right)\left[-d\tau^{2}+d\vec{x}^{2}\right]$, where $a$
is the scale factor and $\tau$ the conformal time. The parameter $f$ is
denoted as the axion decay constant and it has mass dimension one. The
coefficient $\alpha$ controls the coupling to the gauge field and it is
dimensionless. The discussion in this section is independent of the specific
inflaton potential, which is only assumed to be flat enough to lead to slow
roll. The specific potential employed in our simulations is described in
Subsection 3.2.1.

The axion mass is tuned to be of the same order as the Hubble rate $H$ during
inflation,
$\delta\equiv\Lambda^{4}/\left(6H^{2}f^{2}\right)~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle<}}{{\sim}}$}}~{}1$.
Then, for $\delta\ll 3$, and to leading order in slow roll [31]

$\dot{\sigma}=\frac{fH\delta}{\cosh\left[\delta
H\left(t-t_{*}\right)\right]},$ (2.2)

where dot denotes derivative with respect to (wrt) physical time $t$, while
$t_{*}$ is the time at which $\dot{\sigma}$ reaches its maximum value, that,
to leading order in slow-roll, occurs when $\sigma_{*}=\frac{\pi}{2}f$, namely
when the axion is in the steepest part of the potential. Throughout this work,
we use the subscript “$*$” to denote quantities evaluated at the time $t_{*}$.
The expression (2.2) shows that the axion roll is significant only for a
number of e-folds $\Delta N\simeq\frac{1}{\delta}$ around this time. The
motion of the axion can results in a large gauge field amplification, whenever
$\xi>1$, where [31]

$\xi\equiv\frac{\alpha\,\dot{\sigma}}{2Hf}=\frac{\xi_{*}}{\cosh\left[{\delta
H\left(t-t_{*}\right)}\right]}=\frac{2\,\xi_{*}}{\left(\frac{a}{a_{*}}\right)^{\delta}+\left(\frac{a_{*}}{a}\right)^{\delta}},$
(2.3)

where $\xi_{*}=\frac{\alpha\,\delta}{2}$ is the maximum value assumed by $\xi$
at the time $t=t_{*}$. We stress that these analytic results are given to
leading order in slow-roll and neglectinig nackreaction from the gauge fields
amplified by the motion of the axion. In particular, they assume a constant
Hubble rate, and therefore the scale factor evolution $a=-\frac{1}{H\tau}$
during inflation (de Sitter geometry). In Section 4, these analytic results
are compared with the full numerical solutions from the lattice.

###  2.2 Gauge field production

The equation of motion for the gauage field mode function following from eq.
2.1 is [31]

$A^{\prime\prime}_{\pm}+\left(k^{2}\pm\frac{4k\xi_{*}}{\tau\left[\left(\frac{\tau}{\tau_{*}}\right)^{\delta}+\left(\frac{\tau_{*}}{\tau}\right)^{\delta}\right]}\right)=0,$
(2.4)

where the axion is assumed to be homogeneous and eq. 2.3 has been used for the
axion time derivative. Prime denotes differentiation wrt conformal time
$\tau$. Ref. [31] obtained an approximate anlaytic solution for the gauge mode
function using the WKB method:

$A_{+}\simeq\left(\frac{-\tau}{8k\xi}\right)^{1/4}\tilde{A}(\tau,k),\quad
A^{\prime}_{+}\simeq\left(\frac{k\xi}{-2\tau}\right)^{1/4}\tilde{A}(\tau,k),$
(2.5)

where:

$\tilde{A}(\tau,k)\equiv
N[\xi_{*},k,\delta]\exp\left[-\frac{4\sqrt{\xi_{*}}}{1+\delta}\left(\frac{\tau}{\tau_{*}}\right)^{\delta/2}\sqrt{-k\tau}\right],$
(2.6)

and $N[\xi_{*},k,\delta]$ is a time-independent constant obtained, for every
$k$, by matching the WKB solution for $A_{+}$ with the numerical solution of
eq. 2.4 at late times $-k\tau\ll 1$.

The amplified gauge field sources perturbations $\delta\sigma$ of the axion
field by their direct coupling in eq. 2.1. These modes are then
gravitationally coupled to (scalar and tensor) pertrubations of the metric, as
well as to the perturbations of the inflaton field.

###  2.3 Inflationary scalar perturbations

We work in the spatially flat gauge, for which (disregarding vector and tensor
perturbations) the spatial components of the metric are unperturbed,
$g_{ij}=a^{2}\left(\tau\right)\delta_{ij}$. The $00$ and $0i$ components of
the metric are non-dynamical, and they can then be integrated out, leading to
the following equations for the rescaled inflaton and axion modes
$Q_{\phi}\equiv a\,\delta\phi$ and $Q_{\sigma}\equiv a\,\delta\sigma$ [31]

$\displaystyle\left(\frac{\partial^{2}}{\partial\tau^{2}}+k^{2}+\tilde{M}^{2}_{\phi\phi}\right)Q_{\phi}+\tilde{M}^{2}_{\phi\sigma}Q_{\sigma}=0,$
$\displaystyle\left(\frac{\partial^{2}}{\partial\tau^{2}}+k^{2}+\tilde{M}^{2}_{\sigma\sigma}\right)Q_{\sigma}+\tilde{M}^{2}_{\sigma\phi}Q_{\phi}=\frac{\alpha
a^{3}}{f}\int\frac{d^{3}x}{(2\pi)^{3/2}}e^{-i\vec{x}\cdot\vec{k}}\vec{E}\cdot\vec{B},$
(2.7)

where the “mass term” in these relations is given by

$M^{2}_{ij}\equiv-\frac{a^{\prime\prime}}{a}\delta_{ij}+a^{2}V_{,ij}+\left(3-\frac{a^{2}\bar{\phi}^{\prime}_{k}\bar{\phi}^{\prime}_{k}}{2a^{\prime
2}}\right)\bar{\phi}^{\prime}_{i}\bar{\phi}^{\prime}_{j}+\frac{a^{3}}{a^{\prime}}(\bar{\phi}^{\prime}_{i}V_{,j}+\bar{\phi}^{\prime}_{j}V_{,i}),$
(2.8)

where the indices $i$ and $j$ run from $1$ ot $2$, and the value $1$ refers to
$\phi$, while $2$ to $\sigma$ (with $\phi_{1}\equiv\phi$,
$\phi_{2}\equiv\sigma$). We also note that in this relation $V$ refers to the
sum of the infaton and the axion potentials, and a comma on the potential
refers to differentiation wrt the corresponding field.

The last term in section 2.3 represents the direct coupling between gauge and
axion modes. 999As detailed in section 3, only the first equation in (2.3) is
added to the system of equations solved in the lattice simulations, since the
set of equations (3.1) already incorporates the coupling between gauge fields
and the axion. This means that we are neglecting the backreaction of inflaton
and metric perturbations on the axion-gauge field dynamics. Additionally, in
the first equation of (2.3), we disregard the gravitational coupling between
gauge modes and inflaton perturbations, as it is subdominant to the
$\delta\phi-\delta\sigma$ coupling considered, as explained in footnote 10 of
Ref. [31]. To obtain analytical predictions, one must substitute the WKB
ansatz for the gauge field, derived in the previous section, into this set of
equations. In appendix A, we solve the second of eqs. (2.3) and derive
analytical predictions for the power spectrum of $\sigma$. The first of eqs.
(2.3) is analytically solved in Ref. [31], which presents results for the
power spectrum of $\phi$. In sections 4 and 5, we will compare these
analytical calculations with the results from lattice simulations.

###  2.4 Backreaction and perturbativity

The perturbative treatment outlined in this section is valid under certain
assumptions, which we will relax thanks to the lattice simulation. In
particular, if gauge field production is too abundant, it can affect the
background dynamics of the axion field, invalidating the standard perturbative
approach. Backreaction and perturbativity bounds for this model are detailed
in Refs. [39, 76], which provide a comprehensive analysis of these
constraints. A recent motivation for this study is the fact that the parameter
range of interest for this model to explain the recently observed
gravitational wave background using PTA experiments lies in the strong
backreaction regime of the theory [59]. In section 5, we will run simulations
in this regime and quantify the deviations from the predictions of standard
perturbation theory.

##  3 Lattice simulation

To simulate the dynamics of the spectator axion-gauge system and its
interaction with the inflaton, we introduce a novel hybrid numerical approach.
This method addresses the fully nonlinear dynamics of the spectator axion-
gauge system while treating the inflaton sector perturbatively. In this
section, we first present the details of this hybrid scheme and discuss its
validity and applicability. Then, we describe the procedure to generate the
initial conditions for the simulation.

###  3.1 Equations of motion

The lattice approach consists of discretizing continuous space onto a grid of
$N_{\rm pts}^{3}$ points separated by the distance $dx=L/N_{\rm pts}$, where
$L$ is the comoving length of the box. We then numerically solve the equations
of motion in real space for each grid point, using periodic boundary
conditions. The nonlinear equations for the axion-gauge spectator sector are
the following [107]:

$\displaystyle\begin{split}&\sigma^{\prime\prime}+2H{\sigma^{\prime}}-\nabla^{2}\sigma+a^{2}\frac{\partial
V}{\partial\sigma}=-a^{2}\frac{\alpha}{4f}F_{\mu\nu}\tilde{F}^{\mu\nu},\\\
&{A}^{\prime\prime}_{0}-\nabla^{2}A_{0}=\frac{\alpha}{f}\epsilon_{ijk}\partial_{k}\sigma\partial_{i}A_{j},\\\
&{A}^{\prime\prime}_{i}-\nabla^{2}A_{i}=\frac{\alpha}{f}\epsilon_{ijk}\sigma^{\prime}\partial_{j}A_{k}-\frac{\alpha}{f}\epsilon_{ijk}\partial_{j}\sigma(A^{\prime}_{k}-\partial_{k}A_{0}),\end{split}$
(3.1)

where $i,j,k\in\\{1,2,3\\}$ are spatial indices, and
$\nabla^{2}=\partial_{j}\partial_{j}$. We choose to work in the Lorenz gauge
$\partial^{\mu}A_{\mu}=0$. Due to their nonlinear nature, these equations are
able to capture the backreaction of the gauge field on the background dynamics
of the axion. Note that, as already mentioned in section 2.3, we are
neglecting the backreaction of inflaton perturbations on the axion-gauge field
dynamics. This assumption is justified as long as inflaton perturbations
remain small, which is true for all the cases considered in this work.

The gravitational interaction between the inflaton and the axion is determined
by the linear section 2.3, that we rewrite in real space as:

$\left(\frac{\partial^{2}}{\partial\tau^{2}}+\nabla^{2}+\tilde{M}^{2}_{\phi\phi}\right)Q_{\phi}=-\tilde{M}^{2}_{\phi\sigma}a\left(\sigma-\langle\sigma\rangle\right).$
(3.2)

The evolution of the background quantities in eqs. 3.1 and 3.2 is determined
by solving the second Friedmann equation and the Klein-Gordon equation for the
background inflaton, which are, respectively:

$\displaystyle\begin{split}&a^{\prime\prime}=\frac{1}{6}\left(\langle\rho\rangle-3\langle
p\rangle\right)a^{3},\\\
&\bar{\phi}^{\prime\prime}+2\frac{a^{\prime}}{a}\bar{\phi}^{\prime}=-a^{2}\frac{\partial
V(\bar{\phi})}{\partial{{\bar{\phi}}}},\end{split}$ (3.3)

where $\langle\rho\rangle$ and $\langle p\rangle$ are the background energy-
density and pressure, which are obtained as lattice averages of all the energy
and pressure contributions.

The set of eqs. 3.1 and 3.2, evaluated at each lattice point, plus the two
eqs. (3.3), forms a system of $6N_{\rm pts}^{3}+2$ second-order differential
equations. We solve these equations using a 4th-order Runge-Kutta time
integrator and the discretization technique developed in Refs. [105, 106, 107,
108], to which we refer for detailed information on the discretization of eqs.
(3.1). Notably, in this approach, we evolve all four components of the gauge
field $A_{\mu}$ as independent variables. Consequently, the Lorenz constraint
$\partial^{\mu}A_{\mu}=0$ is not automatically preserved on the lattice
(barring infinite accuracy) and must be verified separately. This verification
is addressed in appendix B, where we show that the Lorenz constraint is well
satisfied throughout the evolution in all the cases considered in this work.

The hybrid nonlinear description presented in this section is valid if the
main inflationary sector remains perturbative. Also the right-hand side of eq.
3.2 needs to remain perturbative, requiring axion perturbations to remain
smaller than their background
$\sigma-\langle\sigma\rangle\ll\langle\sigma\rangle$. While this is guaranteed
in the linear regime of the theory, we need to check that it remains true even
in the case of strong backreaction from the gauge field. We will show this in
section 5.

###  3.2 Initial conditions and resolution

The initial conditions on the lattice are set perturbatively. We initialize
all fields as their background values plus some space-dependent fluctuations,
that we discuss separately in the following.

####  3.2.1 Background

In order to set the initial values for the inflaton background, we need to
choose a particular slow-roll potential $V(\phi)$. We assume the following
$\alpha$-attractor potential [116, 117]:

$V(\phi)=\bar{V}\left(1-e^{-\sqrt{\frac{2}{3\alpha}}\phi}\right)^{2},$ (3.4)

where we set $\bar{V}=1.348\times 10^{-11}$ and $\alpha=0.1$. The initial
values of the background inflaton and its velocity are set to
$\bar{\phi}_{0}=2.53$ and $\dot{\bar{\phi}}_{0}=-0.0043\sqrt{\bar{V}}$, so
that there are $N\simeq 52$ $e$-folds left before the end of inflation.
101010We work in reduced Planck mass units $M_{\rm Pl}=1$. This choice
translates to

$n_{s}=0.961,\quad A_{s}=2\times 10^{-9},\quad r=4.5\times 10^{-4}\,,$ (3.5)

for, respectively, the spectral tilt and the amplitude of the scalar
primordial perturbations and for the tensor-to-scalar ratio [12]. We emphasize
that, although we must select a specific potential for the lattice simulation,
the results presented in this work do not crucially depend on the particular
form of the potential, as long as it provides the conditions to sustain slow-
roll inflation.

The background value of the axion field is set to
$(\sigma_{0},\dot{\sigma}_{0})=(0.89\,f,0.20\,f\sqrt{\bar{V}})$. 111111The
initial values for the background inflaton, the background axion, and their
time derivatives are obtained from a numerical integration of about $2$
e-folds of the exact equations for the homogeneous inflaton and axion, and for
the scale factor (without including the gauge fields), that starts from slow-
roll initial conditions. With this choice, there are $N\simeq 1.84$ $e$-folds
left before the axion reaches its maximum velocity. The gauge field background
is set to zero at the beginning of the simulation. The initial value of the
scale factor $a_{0}$ is chosen so that $a=1$ when the axion reaches its
maximum velocity in the simulation.

####  3.2.2 Perturbations

Perturbations of the axion, inflaton and gauge field are initiated in the
Bunch-Davies vacuum at the beginning of the simulation:

$\delta\phi(\vec{k},\tau)=\frac{1}{\sqrt{2k}a}e^{-i{k}\tau},\quad\delta\sigma(\vec{k},\tau)=\frac{1}{\sqrt{2k}a}e^{-i{k}\tau},\quad
A_{\pm}(\vec{k},\tau)=\frac{1}{\sqrt{2k}}e^{-i{k}\tau},$ (3.6)

where $A_{\pm}$ are the gauge field polarizations. The initial fluctuations
are the result of a statistical random process, which sets the power spectrum
equal to the Bunch-Davies one. More details regarding the initialization
procedure, and how it accounts for discretrization and finite-size effects,
can be found in Refs. [108, 106].

The lattice size is set to $N_{\rm pts}^{3}=256^{3}$ number of points and
comoving length of $L=25/(a_{*}H_{*})$. This translates to the following IR
and UV comoving cutoffs for the simulation [105]:

$k_{IR}=\frac{2\pi}{L}\simeq 0.25\,k_{*},\quad\quad k_{UV}\simeq 35\,k_{*},$
(3.7)

where $k_{*}=a_{*}H_{*}$ is the mode exiting the horizon at $t_{*}$, acting as
a reference scale (as just mentioned, $a_{*}=1$). We recall that $t_{*}$
indicates the time at which the axion reaches its maximum velocity in the full
lattice simulation. We run the simulation until $N=\log(a/a_{*})\simeq 4$,
which makes all the modes in our lattice super-horizon at the final time. In
the main text, we show results from simulations with this resolution. In
appendix B we show that our results do not depend on the spatial resolution by
running simulations with different IR and UV cutoffs.

##  4 Weak backreaction

We now present simulation results for a scenario where the backreaction of the
produced gauge quanta on the background dynamics is small, and where linear
theory predictions are expected to be reliable. In this regime, our
simulations can quantitatively test the validity of the WKB approximation
(2.5). Unlike the analytical predictions, our lattice calculation does not
assume a constant $H$ and it does not employ the slow-roll approximation of
eq. 2.2 for the axion dynamics. This allows for a more accurate background
solution, and ultimately for more detailed characterization of the scale-
dependence of the gauge fields amplified by the axion, and for the scalar
fluctuations sourced by the gauge fields. Additionally, we use the simulation
to investigate the role of high-order non-Gaussianity beyond the bispectrum.

For this simulation, we set $\alpha=20.14$, $f=0.1$ and $\delta=0.5$. This
translates into $\Lambda^{4}=6\,\delta H^{2}f^{2}$ and a maxmum value of
$\xi=5$ in the background numerical simulation that neglects the backreaction
from gauge fields.121212To determine the value of $\alpha$ corresponding to
$\xi_{*}=5$, we do not use the slow-roll approximation, which would give
$\alpha=2\,\xi_{*}/\delta=20$. Instead, we set a more precise value by
numerically solving the background dynamics of this model starting from the
slow-roll solution 2 $e$-folds before the simulation starting time, neglecting
backreaction. This is analogous to what we do to set the initial background
values, as explained in section 3.2.1. With this choice of parameters, the
backreaction of the gauge field on the background axion dynamics is expected
to be moderate [76, 39].131313Ref. [39] reports perturbativity bounds for
$\delta=0.6$, that we use as conservative bounds for $\delta=0.5$.

###  4.1 Background dynamics

We begin by presenting the evolution of background quantities, which are
obtained from the simulation as averages over the $N^{3}_{\rm pts}$ lattice
points. In fig. 2 we illustrate the time evolution of the background
quantities $\langle\dot{\sigma}\rangle$ and $\langle\xi\rangle$, and compare
them with the de Sitter approximation given by eq. 2.3, that disregards the
backreaction of the gauge fields. 141414In comparing numerical vs. analytic
results, we set $N=0$ when both the numerical and analytical values of
$\dot{\sigma}$ are at their maximum. Although the discrepancy between the
lattice results and the analytical approximation is small, the axion-gauge
dynamics are exponentially sensitive to slight variations in $\xi$. In the
following section, we will show how this difference impacts the tachyonic
amplification of the gauge field.

$N$ $N$ $\frac{\langle\dot{\sigma}\rangle}{f\sqrt{\bar{V}}}$
$\langle\xi\rangle$ Figure 2: Evolution of the average value of $\dot{\sigma}$
(left) and $\xi$ (right) from the lattice simulation, compared with the
analytical estimate of eq. 2.3.

###  4.2 Gauge field enhancement

In this subsection, we present the results for the growth of the gauge field
on the lattice, focusing on the enhanced polarization $A_{+}$. Specifically,
we compute the following comoving energy density 151515See Refs. [105, 108]
for details on how we compute power spectra on the lattice.

$\rho_{+}(k,\tau)=\langle|A^{\prime}_{+}(k,\tau)|^{2}\rangle+k^{2}\langle|A_{+}(k,\tau)|^{2}\rangle,$
(4.1)

and compare it to the WKB solution given by eq. 2.5, as well as to the linear
solution obtained by numerically solving eq. 2.4. In fig. 3, we illustrate the
time evolution of $\rho_{+}(k,\tau)$ for a characteristic mode
$\tilde{k}\simeq 4k_{*}$. From this plot, we can see that the lattice
correctly reproduces the linear prediction, with small differences due to the
different evolution of $\xi$ (see the previous section).

$\rho_{+}(\tilde{k})$ $N$ Figure 3: Evolution of $\rho_{+}(k,\tau)$, as
defined in eq. 4.1, for a representative mode $\tilde{k}=4k_{*}$ in the case
of weak backreaction. The lattice solution (blue line) is compared with the
WKB solution (green dashed line) and the linear solution (red dashed line).
$\frac{\rho_{+}(k)k^{3}}{2\pi^{2}}$ $k/k_{*}$ $k/k_{*}$
$\frac{\rho_{+}(k)k^{3}}{2\pi^{2}}$ $k/k_{*}$ $k/k_{*}$ Figure 4: Plot of the
excited gauge field polarization $\rho_{+}$, as defined in eq. 4.1, at four
different times during the simulation in the case of weak backreaction. The
lattice results are compared with linear theory and the WKB approximation. The
vertical line indicates $\tilde{k}\simeq 4k_{*}$, the representative mode
shown in fig. 3. $\mathcal{P}_{\sigma/f}(k)$ $k/k_{*}$ Figure 5: Power
spectrum of the axion field $\sigma$ in the case of weak backreaction.
Different colors correspond to different times during the lattice simulation,
as indicated by the legend. The final simulation result is compared with the
analytical prediction from the WKB approximation (red dashed line), that we
computed in appendix A.

In fig. 4, we show the full spectrum of $\rho_{+}(k)$ at four distinct times
during the lattice simulation. We compare these results to the predictions
from linear theory and the WKB solution for all modes. The discrepancy between
the WKB and the numerical solution visible in the last two figures is due to
the fact that the former is valid only after significant amplification of the
gauge field has occurred [31]. In addition, our findings indicate that,
compared to linear theory, the gauge field exhibits less enhancement at large
$k$. This is consistent with the right panel of fig. 2, which shows that $\xi$
is smaller at later times compared to the analytical estimate.

###  4.3 Axion perturbation

#### Power spectrum

In fig. 5, we show the evolution of the power spectrum of the axion field,
introduced in Fourier space as 161616Details on how we compute power spectra
from the lattice simulation can be found in Refs. [105, 108].

$\left\langle\frac{\sigma(\vec{k})}{f}\,\frac{\sigma(\vec{k}^{\prime})}{f}\right\rangle=\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\sigma/f}(k)\delta_{D}(\vec{k}+\vec{k}^{\prime})\;,$
(4.2)

where $\delta_{D}$ is the Dirac delta and $k=|\vec{k}|$.

The $\propto k^{2}$ profile visible at the earlier times corresponds to the
Bunch-Davies vacuum initial conditions (3.6). The axion field is then enhanced
due to the exponential growth of the gauge field, with the most significant
growth occurring between $N=0$ and $N=2$. The final profile has an enhancement
(a “bump”) corresponding to the scales that left the horizon during the
simulation, namely while the axion had a significant roll. We compare the
final time lattice result with the analytical estimate for the power spectrum,
depicted as a dashed red line. We observe that the analytical estimate
correctly predicts the amplitude of the power spectrum. However, the
analytical solution predicts a different scale dependence. This discrepancy
arises because the analytical calculation neglects the time dependence of $H$,
resulting in a different time dependence for $\xi$. Specifically, scalar modes
enhanced at the beginning of the simulation (small $k$) are expected to have a
larger amplitude due to the higher value of $\xi$ (see fig. 2). Conversely,
modes excited later (large $k$) should exhibit a smaller amplitude,
corresponding to the lower value of $\xi$.

#### Non-Gaussianity

$\delta\sigma/\sqrt{\langle\delta\sigma^{2}\rangle}$ Figure 6: 2D snapshot of
the 3D lattice simulation in the case of negligble backreaction. The snapshots
display the values of axion fluctuations
$\delta\sigma=\sigma-\langle\sigma\rangle$ normalized by their standard
deviation. The four panels correspond to different $e$-folding times $N$. $N$
$N$ $\delta\sigma/\sqrt{\langle\delta\sigma^{2}\rangle}$ Figure 7: Non-
Gaussianity of $\sigma$ in the case of weak backreaction. Left panel:
Ridgeline plot showing the normalized one-point PDF of
$\delta\sigma=\sigma-\langle\sigma\rangle$ at various simulation times,
indicated by the color bar. The vertical axis, representing the probability
(namely, the fraction of grid points having the value indicated in the
horizontal axis), is shown only for the first PDF to enhance clarity, while it
is omitted for the other ones. Right panel: Evolution of the cumulants defined
in eq. 4.3.

The simulation allows us to access the real-space field distribution of the
axion, as shown in fig. 6. These simulation boxes can be used to compute the
non-Gaussianity of scalar statistics. In the left panel of fig. 7 , we present
the 1-point probability density function (PDF) of the axion field at different
times during the simulation. The distributions are calculated as normalized
histograms of the field values across the $N_{\rm pts}^{3}$ lattice points. We
observe that the field distribution deviates significantly from Gaussian
statistics. In particular, we notice an exponential tail in the distribution.
This behavior is similar to what occurs in the most minimal version of this
model, where the inflaton coincides with the axion [107].

To quantify the deviation from Gaussian statistics, we introduce the following
dimensionless cumulants [118]:

$\displaystyle\begin{split}&\kappa^{(\delta\sigma)}_{3}\equiv\frac{\langle\delta\sigma^{3}\rangle}{\langle\delta\sigma^{2}\rangle^{3/2}},\quad\kappa^{(\delta\sigma)}_{4}\equiv\frac{\langle\delta\sigma^{4}\rangle-3\langle\delta\sigma^{2}\rangle^{2}}{\langle\delta\sigma^{2}\rangle^{2}},\quad\kappa^{(\delta\sigma)}_{5}\equiv\frac{\langle\delta\sigma^{5}\rangle-10\langle\delta\sigma^{3}\rangle\langle\delta\sigma^{2}\rangle}{\langle\delta\sigma^{2}\rangle^{5/2}},\\\
&\kappa^{(\delta\sigma)}_{6}\equiv\frac{\langle\delta\sigma^{6}\rangle-10\langle\delta\sigma^{3}\rangle^{2}+30\langle\delta\sigma^{2}\rangle^{3}-15\langle\delta\sigma^{4}\rangle\langle\delta\sigma^{2}\rangle}{\langle\delta\sigma^{2}\rangle^{3}},\end{split}$
(4.3)

that all vanish in the Gaussian case. In the right panel of fig. 7, we show
the evolution of these cumulants during the simulation. We can see that
higher-order cumulants are very large and do not converge. This behavior
reflects the presence of the exponential tail in the PDF, and it is analogous
to what is reported in Ref. [107] for the minimal model in which the axion is
the inflaton. Note that this implies that the bispectrum contains very limited
information about the non-Gaussianity of the axion field.

The exponential tail of the final-time PDF of $\delta\sigma$ can be fitted by
the function $f(x)=e^{-mx}$, where $m\simeq 1.2$. It would be interesting to
investigate how the quantity $m$ depends on theoretical parameters such as
$\xi_{*}$ and $\delta$. We leave this exploration for future work. Note that
the tail of the PDF is different from what we would expect from a $\chi^{2}$
distribution (which is the simplest guess that is reasonable to make for the
PDF that results from the convolution of two gaussian gauge modes [29]), which
also has an exponential behaviour in the tail, but with $f(x)\propto
e^{-x/2}/\sqrt{x}$, for $x\gg 1$ [29].

###  4.4 Curvature perturbation

With the choice of parameters adopted in this section, the curvature
perturbation $\zeta$ sourced by the axion field is of the same order as the
standard vacuum contribution, according to the analytical estimate of Ref.
[31]. More precisely, the sourced contribution estimated analytically is
approximately 40% of the vacuum one at the peak scale. This is consistent with
the simulation results, as shown in fig. 8. The plot shows that the power
spectrum from the simulation remains close to the nearly scale-invariant
prediction with amplitude $\mathcal{P}_{\zeta,\rm vacuum}\simeq 2\times
10^{-9}$. The power spectrum shows only a small peak of about 20% at the peak
scale, roughly consistent with the analytical estimate. The comoving curvature
perturbation is computed from the lattice using the linear relation
$\zeta=-\delta\phi H/\dot{\phi}$.

$\mathcal{P}_{\zeta}(k)$ $k/k_{*}$ Figure 8: Power spectrum of $\zeta$ in the
case of weak backreaction. At the final time, we compare the lattice result
with the analytical prediction computed in [31], depicted as a red dashed
line.

###  4.5 Energy contributions

Before proceeding to the strong backreaction case, we present the evolution of
the different contributions to the mean energy density $\langle\rho\rangle$.
In fig. 9, we plot the evolution of all contributions during the simulation.
We observe that the inflaton potential remains the dominant component.
Additionally, the kinetic energy of the axion, $K_{\sigma}$, remains much
larger than the energy density of the gauge field. This is expected in the
regime with small backreaction. In this figure, the dotted purple and green
lines represent the vacuum contributions to the gauge field energy density and
the gradient energy of the axion field, respectively. These vacuum
contributions diverge in the limit of infinite UV resolution and should be
subtracted. However, they play a negligible role in the evolution of the
system. This is confirmed in appendix B, where we show that changing the UV
cutoff, which corresponds to changing the amount of vacuum energy, does not
affect any of our results.

Finally, we note that once the gauge field is enhanced, it does not
immediately source axion perturbations. Axion particle production occurs
roughly $1.5$ $e$-folds later. This delay is crucial for understanding how
backreaction operates in this system, as discussed in section 5.1.

$N$ Figure 9: Evolution of different components contributing to the total
energy density. For a field $f$, $K_{f}$ and $G_{f}$ represent its kinetic and
gradient energy, respectively, while $V_{f}$ denotes its potential energy. The
purple curve illustrates the total energy density of the gauge field. These
quantities are computed as average energies over the lattice simulation. The
dashed lines represent the vacuum contributions, as explained in section 4.5.

##  5 Strong backreaction

We now turn to the case of strong backreaction. To do this, we maintain the
same parameters as in the previous section and increase the axion-gauge
coupling to $\alpha=24.17$, corresponding to a maximum value $\xi=6$
(according to the background simulation, see the beginning of section 4). This
is the value required to explain the PTA signal, according to the analytical
calculations that neglect backreaction [59]. In this regime, however, we
expect significant deviations from the linear theory and the WKB approximation
[76, 39, 59].

###  5.1 Background evolution

$N$ $N$ $\frac{\langle\dot{\sigma}\rangle}{f\sqrt{\bar{V}}}$
$\langle\xi\rangle$ Figure 10: Evolution of $\dot{\sigma}$ (left panel) and
$\xi$ (right panel) in the case of strong backreaction, compared with the
analytical estimate of eq. 2.3.

In fig. 10 we show the evolution of $\dot{\sigma}$ and $\xi$ during the
lattice simulation and compare it to the linear theory prediction. The
background dynamics is clearly affected by the gauge field growth. The
deviation manifests itself as an oscillatory feature about the background
solution with a half-period of $\Delta N\simeq 1.5$ e-folds. This oscillation
is caused by the fact that, when the gauge field is enhanced, it sources
scalar perturbations with a delay, as we discussed already at the end of
section 4.5. This behavior is similar to what happens in the minimal
axion-U(1) model of inflation, where the axion is the inflaton [81, 82, 83,
119, 107, 110, 80], explained in section 1 (see fig. 1).

###  5.2 Gauge field enhancement

In fig. 11 we show the time evolution of $\rho_{+}$, defined in eq. 4.1, for
one characteristic mode $\tilde{k}\simeq 4k_{*}$. From this plot, we can see
that the lattice solution deviates from the analytical and linear estimates,
giving a much lower value of $\rho_{+}(\tilde{k})$ for $N>2$. In fig. 12, we
plot the power spectrum of $\rho_{+}$ at different times during the
simulation. These plots show that the WKB approximation matches well the
linear computation at late times. However, both approximations overestimate
the amplitude of the power spectrum resulting from the lattice dynamics at
times when the gauge field is excited the most. This is due to the fact that
now the backraction of the gauge fields and the axion dynamics (not included
in the WKB and in the linear computaton) is significant.

$\rho_{+}(\tilde{k})$ $N$ Figure 11: Same as fig. 3, but in the case of strong
backreaction. In this case, the lattice solution (blue) is suppressed compared
to the WKB solution (green dashed) and to the linear solution (red dashed
line). $\frac{\rho_{+}(k)k^{3}}{2\pi^{2}}$ $k/k_{*}$ $k/k_{*}$
$\frac{\rho_{+}(k)k^{3}}{2\pi^{2}}$ $k/k_{*}$ $k/k_{*}$ Figure 12: Same as
fig. 4, but in the case of strong backreaction. The vertical line shows
$\tilde{k}\simeq 4k_{*}$, the representative mode depicted in fig. 11.

These results show that the gauge field growth is significantly suppressed by
the backreaction. This is expected because, as we have seen in the last
section, the growth of the gauge field dampens the velocity of the axion,
resulting in a lower value of $\xi$ during the evolution (see fig. 10).

###  5.3 Axion perturbation

#### Power spectrum

As a result of suppressed gauge field growth, backreaction also suppresses the
power spectrum of the axion field. This effect is illustrated in fig. 13,
which displays the power spectrum of $\sigma$ computed from the lattice
simulation, compared to the analytical prediction at the final time. Note
that, even if we are in the strong backreaction regime, the axion power
spectrum is still perturbative, $\mathcal{P}_{\sigma/f}\ll 1$. Therefore, we
can still trust our hybrid numerical approach in the calculation of
$\delta\phi$, which will be computed in section 5.4.

$\mathcal{P}_{\sigma/f}(k)$ $k/k_{*}$ Figure 13: Power spectrum of the axion
field $\sigma$ in the scenario of strong backreaction, depicted at various
times during the lattice simulation. The final simulation result (blue) is
compared with the prediction from the WKB approximation (red dashed line).

#### Non-Gaussianity

$\delta\sigma/\sqrt{\langle\delta\sigma^{2}\rangle}$ Figure 14: Same as fig.
6, but in the case of strong backreaction, showing values of the axion field
$\sigma$ across a 2D slice of the simulation. $N$ $N$
$\delta\sigma/\sqrt{\langle\delta\sigma^{2}\rangle}$ Figure 15: Same as fig.
7, showing non-Gaussianity of $\sigma$ in the case of strong backreaction. The
right panel includes an inset that provides a detailed view of the region
$N>2.5$.

In fig. 14, we present snapshots of our simulation boxes for the axion field
$\sigma$ in the case of strong backreaction. As we did in section 4.3, we
compute the 1-point PDFs of the axion field across the simulation and analyze
their cumulants. The results are displayed in fig. 15.

We observe that non-Gaussianity is significantly suppressed compared to the
case with weak backreaction. This suppression is consistent with the
interpretation given in Ref. [107], which attributes the reduction in non-
Gaussianity to the large number of excited gauge field modes, leading to
Gaussian statistics due to the central limit theorem. Notably, unlike the
minimal model discussed in Ref. [107], non-Gaussianity remains sizable and
$O(1)$ at the end of the simulation. This is consistent with the central limit
theorem interpretation, as the strong scale dependence in this model (due to
the fact that the axion roll is significant only for a few e-folds) limits the
number of excited gauge field modes compared to the minimal scenario.

###  5.4 Curvature perturbation

We now come to the results for the comoving curvature perturbation, that we
compute from the lattice using the linear relation $\zeta=-\delta\phi
H/\dot{\phi}$.

$\mathcal{P}_{\zeta}(k)$ $k/k_{*}$ Figure 16: Power spectrum of $\zeta$ in the
case of strong backreaction. Different colors correspond to different
simulation times, while the red dashed line shows the analytical calculation.

#### Power spectrum

In fig. 16, we present the power spectrum of $\zeta$ at various times during
the lattice simulation and compare it to the WKB approximation at the final
time. The scalar power spectrum is significantly suppressed relative to the
analytical estimate. This suppression is a consequence of the reduced gauge
field growth due to strong backreaction. The sourced contribution is still
much larger than the standard vacuum contribution $\mathcal{P}_{\zeta,\rm
vacuum}=2\times 10^{-9}$.

#### Non-Gaussianity

$\zeta/\sqrt{\langle\zeta^{2}\rangle}$ Figure 17: Simulation snapshots for
$\zeta$ in the case of strong backreaction. $N$ $N$
$\zeta/\sqrt{\langle\zeta^{2}\rangle}$ Figure 18: Non-Gaussianity of $\zeta$
in the case of strong backreaction. See the caption of fig. 7 for the
explanation of the different panels.

We now discuss the non-Gaussianity of $\zeta$. In fig. 17, we present 2D
snapshots from the simulation. From the full 3D simulation boxes, we extract
the 1-point PDF of $\zeta$ and its cumulants, as shown in fig. 18. Non-
Gaussianity of the observable $\zeta$ on super-horizon scales is large and
non-trivial. In particular, we observe that a large amount of non-Gaussianity
is encoded in cumulants beyond the 3rd order. This implies that the bispectrum
will carry limited information about the statistical properties of the
primordial curvature fluctuation. This has crucial phenomenological
implication, as we discuss in the Conclusion (section 6).

###  5.5 Energy contributions

$N$ Figure 19: Energy contributions in the case of strong backreaction. See
the caption of fig. 9 for the definition of the different components.

Finally, in fig. 19, we show the evolution of the various contributions to the
energy density, which source the evolution of the Universe through the
Friedmann equations. We can see that, even in the strong backreaction regime,
the energy density of the Universe is dominated by the background inflaton
field. However, in this scenario, the mean energy density of the gauge field
becomes comparable to the kinetic energy of the axion field $K_{\sigma}$
around $N=2$. This indicates that the gauge field cannot be neglected in the
background dynamics of the axion field, signaling significant backreaction in
the system. A close inspection to the $K_{\sigma}$ profile shows the
oscillatory behaviour of the axion kinetic energy at the moment in which the
backreaction is more relevant, as we also observed in fig. 10. Moreover, we
observe that the gradient energy of the axion field (which is the cause of the
suppression of the oscillatory behaviour observed in the simulations of [110]
in the model in which the inflaton axion is continuously rolling) remains
small in this simulation.

##  6 Conclusions

Understanding axion-gauge systems during inflation is of utmost importance due
to the theoretical and observational interest in these models. In particular,
it is important to develop a nonlinear and non-perturbative understanding of
these systems, as the observationally relevant regime is often beyond the
validity of standard perturbation theory. Lattice simulations of inflation are
starting to make significant progress in this direction. In this work, we
considered a fast-rolling axion coupled to a gauge field, both minimally
coupled to the main inflaton sector. We performed the first numerical lattice
simulation of this system, employing a novel numerical scheme.

In the regime where standard perturbation theory techniques are valid, we used
the simulation to verify the accuracy of analytical approximations. While we
found that the magnitude of the sourced axion fluctuations is consistent with
analytical predictions, the lattice simulation provides a more precise
characterization of the scale dependence of this signal. This improvement is
due to our ability to relax the assumption of a constant Hubble rate $H$ and
to go beyond the slow-roll approximation for the dynamics of the axion. More
importantly, the simulations allow us to understand non-Gaussianity resulting
from the nonlinear axion-gauge interactions beyond just the three-point
function. We showed that scalar fluctuations exhibit a non-Gaussian
exponential tail, with high-order information being crucial for describing the
statistics. This finding is consistent with results from the minimal version
of this model, where the axion is the inflaton and the sourced contribution
does not have such a strong scale dependence. From this, we learn that the
exponential tail is a general feature of this model and is not linked to the
scale dependence of the signal.

After exploring the linear regime, we investigated the regime strongly
affected by backreaction, where standard perturbation theory is not expected
to provide an accurate picture. This regime is extremely relevant because many
interesting observational consequences, such as a significant gravitational
wave background at PTA scales, lie precisely in this regime. We found that
backreaction strongly suppresses the growth of the gauge field. As a
consequence, the amplitude of sourced scalar perturbations is also suppressed.
Moreover, we showed that the exponential tail in the scalar distribution is
suppressed by backreaction. This is due to the large number of excited gauge
field modes, which tend to ‘Gaussianize’ the system due to the central limit
theorem. However, we found that the suppression is less severe compared to the
minimal model where the axion is the inflaton, due to the strong scale
dependence in this model. Therefore, while the scale dependence does not
affect the presence of an exponential tail, it influences the suppression of
this tail in the strong backreaction regime. As a result, high-order non-
Gaussianity can still play an important role in this system when backreaction
is at play.

These non-perturbative non-Gaussianities have crucial observational
implications. First, high-order information, neglected in previous studies, is
expected to play an important role in shaping the large-scale structure of our
Universe. This will have observational consequences for the statistical
properties of large-scale structures and the CMB (see, for example, [120, 121,
122]). Moreover, the non-trivial non-Gaussian statistics will have important
implications for the formation of primordial black holes on small scales,
which is highly sensitive to the statistical properties of the tail of the
distribution [123, 124, 125, 126, 48].

As already mentioned, this model can contribute to the gravitational wave
signal observed by PTA experiments. However, this requires working in the
strong backreaction regime where analytical calculations cannot be trusted, as
discussed in Ref. [59]. Although we did not calculate gravitational waves
directly, we showed that backreaction significantly suppresses the growth of
the gauge field in this regime, which is the primary source of the
gravitational wave background. Therefore, we conclude that the parameter
choices outlined in Ref. [59] are likely to result in a gravitational wave
signal that is too weak to be detected by PTA experiments, due to the
suppression effect observed in our simulations. Future dedicated work,
employing higher-resolution lattice simulations, will be necessary to identify
the viable parameter space for a detectable PTA signal. In this paper, we
focused on the theoretical understanding of this system. All the
phenomenological and observational prospects outlined in this Conclusion are
left for future work.

## Acknowledgements

We are grateful to Drew Jamieson, Eiichiro Komatsu, Ippei Obata, Sébastien
Renaux-Petel, Caner Ünal, and Denis Werth for insightful discussions and
comments. AC acknowledges funding support from the Initiative Physique des
Infinis (IPI), a research training program of the Idex SUPER at Sorbonne
Université. MP acknowledges support from Istituto Nazionale di Fisica Nucleare
(INFN) through the Theoretical Astroparticle Physics (TAsP) project, and from
the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando
2022 - grant 20228RMX4A.

## References

  * [1] A.H. Guth, _Inflationary universe: A possible solution to the horizon and flatness problems_ , _Phys. Rev. D_ 23 (1981) 347.
  * [2] K. Sato, _First Order Phase Transition of a Vacuum and Expansion of the Universe_ , _Mon. Not. Roy. Astron. Soc._ 195 (1981) 467.
  * [3] A.D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Adv. Ser. Astrophys. Cosmol._ 3 (1987) 149.
  * [4] A. Albrecht and P.J. Steinhardt, _Cosmology for grand unified theories with radiatively induced symmetry breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220.
  * [5] A. Starobinsky, _A new type of isotropic cosmological models without singularity_ , _Physics Letters B_ 91 (1980) 99 .
  * [6] A.A. Starobinsky, _Spectrum of relict gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682.
  * [7] V.F. Mukhanov and G.V. Chibisov, _Quantum Fluctuations and a Nonsingular Universe_ , _JETP Lett._ 33 (1981) 532.
  * [8] S. Hawking, _The development of irregularities in a single bubble inflationary universe_ , _Physics Letters B_ 115 (1982) 295 .
  * [9] A.H. Guth and S.-Y. Pi, _Fluctuations in the new inflationary universe_ , _Phys. Rev. Lett._ 49 (1982) 1110.
  * [10] A. Starobinsky, _Dynamics of phase transition in the new inflationary universe scenario and generation of perturbations_ , _Physics Letters B_ 117 (1982) 175 .
  * [11] L. Abbott and M.B. Wise, _Constraints on Generalized Inflationary Cosmologies_ , _Nucl. Phys. B_ 244 (1984) 541.
  * [12] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [13] Planck collaboration, _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9 [1905.05697].
  * [14] BICEP, Keck collaboration, _Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season_ , _Phys. Rev. Lett._ 127 (2021) 151301 [2110.00483].
  * [15] M. Biagetti, M. Fasiello and A. Riotto, _Enhancing Inflationary Tensor Modes through Spectator Fields_ , _Phys. Rev. D_ 88 (2013) 103518 [1305.7241].
  * [16] M. Biagetti, E. Dimastrogiovanni, M. Fasiello and M. Peloso, _Gravitational Waves and Scalar Perturbations from Spectator Fields_ , _JCAP_ 04 (2015) 011 [1411.3029].
  * [17] T. Fujita, J. Yokoyama and S. Yokoyama, _Can a spectator scalar field enhance inflationary tensor mode?_ , _PTEP_ 2015 (2015) 043E01 [1411.3658].
  * [18] D. Cannone, G. Tasinato and D. Wands, _Generalised tensor fluctuations and inflation_ , _JCAP_ 01 (2015) 029 [1409.6568].
  * [19] J.L. Cook and L. Sorbo, _Particle production during inflation and gravitational waves detectable by ground-based interferometers_ , _Phys. Rev. D_ 85 (2012) 023534 [1109.0022].
  * [20] L. Senatore, E. Silverstein and M. Zaldarriaga, _New Sources of Gravitational Waves during Inflation_ , _JCAP_ 08 (2014) 016 [1109.0542].
  * [21] N. Barnaby, J. Moxon, R. Namba, M. Peloso, G. Shiu and P. Zhou, _Gravity waves and non-Gaussian features from particle production in a sector gravitationally coupled to the inflaton_ , _Phys. Rev. D_ 86 (2012) 103508 [1206.6117].
  * [22] D. Carney, W. Fischler, E.D. Kovetz, D. Lorshbough and S. Paban, _Rapid field excursions and the inflationary tensor spectrum_ , _JHEP_ 11 (2012) 042 [1209.3848].
  * [23] A. Papageorgiou, M. Peloso and C. Unal, _Nonlinear perturbations from axion-gauge fields dynamics during inflation_ , _JCAP_ 07 (2019) 004 [1904.01488].
  * [24] S. Mukohyama, R. Namba, M. Peloso and G. Shiu, _Blue Tensor Spectrum from Particle Production during Inflation_ , _JCAP_ 08 (2014) 036 [1405.0346].
  * [25] N. Barnaby and M. Peloso, _Large Nongaussianity in Axion Inflation_ , _Phys. Rev. Lett._ 106 (2011) 181301 [1011.1500].
  * [26] L. Sorbo, _Parity violation in the Cosmic Microwave Background from a pseudoscalar inflaton_ , _JCAP_ 06 (2011) 003 [1101.1525].
  * [27] N. Barnaby, R. Namba and M. Peloso, _Phenomenology of a Pseudo-Scalar Inflaton: Naturally Large Nongaussianity_ , _JCAP_ 04 (2011) 009 [1102.4333].
  * [28] 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. D_ 85 (2012) 023525 [1110.3327].
  * [29] A. Linde, S. Mooij and E. Pajer, _Gauge field production in supergravity inflation: Local non-Gaussianity and primordial black holes_ , _Phys. Rev. D_ 87 (2013) 103506 [1212.1693].
  * [30] T. Fujita, K. Mukaida and Y. Tada, _Stochastic formalism for U(1) gauge fields in axion inflation_ , _JCAP_ 12 (2022) 026 [2206.12218].
  * [31] R. Namba, M. Peloso, M. Shiraishi, L. Sorbo and C. Unal, _Scale-dependent gravitational waves from a rolling axion_ , _JCAP_ 01 (2016) 041 [1509.07521].
  * [32] E. Dimastrogiovanni, M. Fasiello, J.M. Leedom, M. Putti and A. Westphal, _Gravitational Axiverse Spectroscopy: Seeing the Forest for the Axions_ , 2312.13431.
  * [33] O. Özsoy, _On Synthetic Gravitational Waves from Multi-field Inflation_ , _JCAP_ 04 (2018) 062 [1712.01991].
  * [34] NASA PICO collaboration, _PICO: Probe of Inflation and Cosmic Origins_ , 1902.10541.
  * [35] S. Shandera et al., _Probing the origin of our Universe through cosmic microwave background constraints on gravitational waves_ , _Bull. Am. Astron. Soc._ 51 (2019) 338 [1903.04700].
  * [36] P. Campeti, D. Poletti and C. Baccigalupi, _Principal component analysis of the primordial tensor power spectrum_ , _JCAP_ 09 (2019) 055 [1905.08200].
  * [37] K. Abazajian et al., _CMB-S4 Science Case, Reference Design, and Project Plan_ , 1907.04473.
  * [38] E. Komatsu, _New physics from the polarized light of the cosmic microwave background_ , _Nature Rev. Phys._ 4 (2022) 452 [2202.13919].
  * [39] P. Campeti, O. Özsoy, I. Obata and M. Shiraishi, _New constraints on axion-gauge field dynamics during inflation from Planck and BICEP/Keck data sets_ , _JCAP_ 07 (2022) 039 [2203.03401].
  * [40] T. Fujita, Y. Minami, M. Shiraishi and S. Yokoyama, _Can primordial parity violation explain the observed cosmic birefringence?_ , _Phys. Rev. D_ 106 (2022) 103529 [2208.08101].
  * [41] LiteBIRD collaboration, _LiteBIRD science goals and forecasts. A case study of the origin of primordial gravitational waves using large-scale CMB polarization_ , _JCAP_ 06 (2024) 008 [2312.00717].
  * [42] M. Putti, N. Bartolo, S. Bhattacharya and M. Peloso, _CMB spectral distortions from enhanced primordial perturbations: the role of spectator axions_ , 2403.08594.
  * [43] J. Garcia-Bellido, M. Peloso and C. Unal, _Gravitational waves at interferometer scales and primordial black holes in axion inflation_ , _JCAP_ 12 (2016) 031 [1610.03763].
  * [44] J. Garcia-Bellido, M. Peloso and C. Unal, _Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter_ , _JCAP_ 09 (2017) 013 [1707.02441].
  * [45] O. Özsoy and Z. Lalak, _Primordial black holes as dark matter and gravitational waves from bumpy axion inflation_ , _JCAP_ 01 (2021) 040 [2008.07549].
  * [46] A. Talebian, S.A. Hosseini Mansoori and H. Firouzjahi, _Inflation from Multiple Pseudo-scalar Fields: Primordial Black Hole Dark Matter and Gravitational Waves_ , _Astrophys. J._ 948 (2023) 48 [2210.13822].
  * [47] O. Özsoy and G. Tasinato, _Inflation and Primordial Black Holes_ , _Universe_ 9 (2023) 203 [2301.03600].
  * [48] LISA Cosmology Working Group collaboration, _Primordial black holes and their gravitational-wave signatures_ , 2310.19857.
  * [49] C. Caprini, D.G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni et al., _Reconstructing the spectral shape of a stochastic gravitational wave background with LISA_ , _JCAP_ 11 (2019) 017 [1906.09244].
  * [50] W. Giarè and A. Melchiorri, _Probing the inflationary background of gravitational waves from large to small scales_ , _Phys. Lett. B_ 815 (2021) 136137 [2003.04783].
  * [51] P. Campeti, E. Komatsu, D. Poletti and C. Baccigalupi, _Measuring the spectrum of primordial gravitational waves with CMB, PTA and Laser Interferometers_ , _JCAP_ 01 (2021) 012 [2007.04241].
  * [52] D.G. Figueroa, A. Florio, N. Loayza and M. Pieroni, _Spectroscopy of particle couplings with gravitational waves_ , _Phys. Rev. D_ 106 (2022) 063522 [2202.05805].
  * [53] LISA Cosmology Working Group collaboration, _Cosmology with the Laser Interferometer Space Antenna_ , _Living Rev. Rel._ 26 (2023) 5 [2204.05434].
  * [54] W. Giarè, M. Forconi, E. Di Valentino and A. Melchiorri, _Towards a reliable calculation of relic radiation from primordial gravitational waves_ , _Mon. Not. Roy. Astron. Soc._ 520 (2023) 2 [2210.14159].
  * [55] S.P. Corbà and L. Sorbo, _Correlated scalar perturbations and gravitational waves from axion inflation_ , 2403.03338.
  * [56] O. Özsoy, A. Papageorgiou and M. Fasiello, _Scale-dependent chirality as a smoking gun for Abelian gauge fields during inflation_ , 2405.14963.
  * [57] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Search for Signals from New Physics_ , _Astrophys. J. Lett._ 951 (2023) L11 [2306.16219].
  * [58] D.G. Figueroa, M. Pieroni, A. Ricciardone and P. Simakachorn, _Cosmological Background Interpretation of Pulsar Timing Array Data_ , _Phys. Rev. Lett._ 132 (2024) 171002 [2307.02399].
  * [59] C. Unal, A. Papageorgiou and I. Obata, _Axion-Gauge Dynamics During Inflation as the Origin of Pulsar Timing Array Signals and Primordial Black Holes_ , 2307.02322.
  * [60] J.-Q. Jiang, Y. Cai, G. Ye and Y.-S. Piao, _Broken blue-tilted inflationary gravitational waves: a joint analysis of NANOGrav 15-year and BICEP/Keck 2018 data_ , _JCAP_ 05 (2024) 004 [2307.15547].
  * [61] C. Fu, J. Liu, X.-Y. Yang, W.-W. Yu and Y. Zhang, _Explaining pulsar timing array observations with primordial gravitational waves in parity-violating gravity_ , _Phys. Rev. D_ 109 (2024) 063526 [2308.15329].
  * [62] X. Niu and M.H. Rahat, _NANOGrav signal from axion inflation_ , _Phys. Rev. D_ 108 (2023) 115023 [2307.01192].
  * [63] E. Dimastrogiovanni, M. Fasiello and T. Fujita, _Primordial Gravitational Waves from Axion-Gauge Fields Dynamics_ , _JCAP_ 01 (2017) 019 [1608.04216].
  * [64] O. Iarygina and E.I. Sfakianakis, _Gravitational waves from spectator Gauge-flation_ , _JCAP_ 11 (2021) 023 [2105.06972].
  * [65] P. Adshead and M. Wyman, _Chromo-Natural Inflation: Natural inflation on a steep potential with classical non-Abelian gauge fields_ , _Phys. Rev. Lett._ 108 (2012) 261302 [1202.2366].
  * [66] A. Maleknejad and M.M. Sheikh-Jabbari, _Gauge-flation: Inflation From Non-Abelian Gauge Fields_ , _Phys. Lett. B_ 723 (2013) 224 [1102.1513].
  * [67] B. Thorne, T. Fujita, M. Hazumi, N. Katayama, E. Komatsu and M. Shiraishi, _Finding the chiral gravitational wave background of an axion-SU(2) inflationary model using CMB observations and laser interferometers_ , _Phys. Rev. D_ 97 (2018) 043506 [1707.03240].
  * [68] Q. Liang, M.-X. Lin, M. Trodden and S.S.C. Wong, _Probing parity violation in the stochastic gravitational wave background with astrometry_ , _Phys. Rev. D_ 109 (2024) 083028 [2309.16666].
  * [69] T. Fujita, T. Murata, I. Obata and M. Shiraishi, _Parity-violating scalar trispectrum from a rolling axion during inflation_ , _JCAP_ 05 (2024) 127 [2310.03551].
  * [70] X. Niu, M.H. Rahat, K. Srinivasan and W. Xue, _Parity-odd and even trispectrum from axion inflation_ , _JCAP_ 05 (2023) 018 [2211.14324].
  * [71] J. Hou, Z. Slepian and R.N. Cahn, _Measurement of parity-odd modes in the large-scale 4-point correlation function of Sloan Digital Sky Survey Baryon Oscillation Spectroscopic Survey twelfth data release CMASS and LOWZ galaxies_ , _Mon. Not. Roy. Astron. Soc._ 522 (2023) 5701 [2206.03625].
  * [72] O.H.E. Philcox, _Probing parity violation with the four-point correlation function of BOSS galaxies_ , _Phys. Rev. D_ 106 (2022) 063501 [2206.04227].
  * [73] O.H.E. Philcox and J. Ereza, _Could Sample Variance be Responsible for the Parity-Violating Signal Seen in the BOSS Galaxy Survey?_ , 2401.09523.
  * [74] D. Jamieson, A. Caravano, J. Hou, Z. Slepian and E. Komatsu, _Parity-Odd Power Spectra: Concise Statistics for Cosmological Parity Violation_ , 2406.15683.
  * [75] R.Z. Ferreira, J. Ganc, J. Noreña and M.S. Sloth, _On the validity of the perturbative description of axions during inflation_ , _JCAP_ 04 (2016) 039 [1512.06116].
  * [76] M. Peloso, L. Sorbo and C. Unal, _Rolling axions during inflation: perturbativity and signatures_ , _JCAP_ 09 (2016) 001 [1606.00459].
  * [77] E. Dimastrogiovanni, M. Fasiello, M. Michelotti and O. Özsoy, _A universal constraint on axion non-Abelian dynamics during inflation_ , 2405.17411.
  * [78] M.M. Anber and L. Sorbo, _Naturally inflating on steep potentials through electromagnetic dissipation_ , _Phys. Rev. D_ 81 (2010) 043534 [0908.4089].
  * [79] A. Berera, _Warm inflation_ , _Phys. Rev. Lett._ 75 (1995) 3218 [astro-ph/9509049].
  * [80] M. Peloso and L. Sorbo, _Instability in axion inflation with strong backreaction from gauge modes_ , _JCAP_ 01 (2023) 038 [2209.08131].
  * [81] S.-L. Cheng, W. Lee and K.-W. Ng, _Numerical study of pseudoscalar inflation with an axion-gauge field coupling_ , _Phys. Rev. D_ 93 (2016) 063510 [1508.00251].
  * [82] A. Notari and K. Tywoniuk, _Dissipative Axial Inflation_ , _JCAP_ 12 (2016) 038 [1608.06223].
  * [83] G. Dall’Agata, S. González-Martín, A. Papageorgiou and M. Peloso, _Warm dark energy_ , _JCAP_ 08 (2020) 032 [1912.09950].
  * [84] V. Domcke, V. Guidetti, Y. Welling and A. Westphal, _Resonant backreaction in axion inflation_ , _JCAP_ 09 (2020) 009 [2002.02952].
  * [85] 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 (2021) 123504 [2109.01651].
  * [86] R. Durrer, O. Sobol and S. Vilchinskii, _Backreaction from gauge fields produced during inflation_ , _Phys. Rev. D_ 108 (2023) 043540 [2303.04583].
  * [87] 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_ , _JHEP_ 11 (2023) 183 [2309.04254].
  * [88] O. Iarygina, E.I. Sfakianakis, R. Sharma and A. Brandenburg, _Backreaction of axion-SU(2) dynamics during inflation_ , _JCAP_ 04 (2024) 018 [2311.07557].
  * [89] D.C. Galanti, P. Conzinu, G. Marozzi and S. Santos da Costa, _Gauge invariant quantum backreaction in U(1) axion inflation_ , 2406.19960.
  * [90] P. Creminelli, S. Kumar, B. Salehian and L. Santoni, _Dissipative inflation via scalar production_ , _JCAP_ 08 (2023) 076 [2305.07695].
  * [91] 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_ , _JCAP_ 01 (2024) 034 [2303.13425].
  * [92] S.Y. Khlebnikov and I.I. Tkachev, _Classical decay of the inflaton_ , _Physical Review Letters_ 77 (1996) 219–222.
  * [93] T. Prokopec and T.G. Roos, _Lattice study of classical inflaton decay_ , _Physical Review D_ 55 (1997) 3768–3775.
  * [94] G. Felder and I. Tkachev, _Latticeeasy: A program for lattice simulations of scalar fields in an expanding universe_ , _Computer Physics Communications_ 178 (2008) 929?932.
  * [95] J. Garcia-Bellido and D.G. Figueroa, _A stochastic background of gravitational waves from hybrid preheating_ , _Phys. Rev. Lett._ 98 (2007) 061302 [astro-ph/0701014].
  * [96] A.V. Frolov, _Defrost: a new code for simulating preheating after inflation_ , _Journal of Cosmology and Astroparticle Physics_ 2008 (2008) 009.
  * [97] Z. Huang, _Art of lattice and gravity waves from preheating_ , _Physical Review D_ 83 (2011) .
  * [98] J. Sainio, _Pycool — a cosmological object-oriented lattice code written in python_ , _Journal of Cosmology and Astroparticle Physics_ 2012 (2012) 038–038.
  * [99] H.L. Child, J.T. Giblin, R.H. Ribeiro and D. Seery, _Preheating with nonminimal kinetic terms_ , _Physical Review Letters_ 111 (2013) .
  * [100] R. Easther, H. Finkel and N. Roth, _Pspectre: a pseudo-spectral code for (p)reheating_ , _Journal of Cosmology and Astroparticle Physics_ 2010 (2010) 025–025.
  * [101] K.D. Lozanov and M.A. Amin, _Gfire—gauge field integrator for reheating_ , _Journal of Cosmology and Astroparticle Physics_ 2020 (2020) 058–058.
  * [102] D.G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, _CosmoLattice_ , 2102.01031.
  * [103] P. Adshead, J.T. Giblin, T.R. Scully and E.I. Sfakianakis, _Magnetogenesis from axion inflation_ , _JCAP_ 10 (2016) 039 [1606.08474].
  * [104] P. Adshead, J.T. Giblin and Z.J. Weiner, _Non-Abelian gauge preheating_ , _Phys. Rev. D_ 96 (2017) 123512 [1708.02944].
  * [105] A. Caravano, E. Komatsu, K.D. Lozanov and J. Weller, _Lattice simulations of inflation_ , _JCAP_ 12 (2021) 010 [2102.06378].
  * [106] A. Caravano, E. Komatsu, K.D. Lozanov and J. Weller, _Lattice simulations of Abelian gauge fields coupled to axions during inflation_ , _Phys. Rev. D_ 105 (2022) 123530 [2110.10695].
  * [107] A. Caravano, E. Komatsu, K.D. Lozanov and J. Weller, _Lattice Simulations of Axion-U(1) Inflation_ , 2204.12874.
  * [108] A. Caravano, _Simulating the inflationary Universe: from single-field to the axion-U(1) model_ , Ph.D. thesis, Munich U., 7, 2022.  2209.13616.  10.5282/edoc.30905.
  * [109] T. Krajewski and K. Turzyński, _Anatomy of geometrical destabilization of inflation_ , _JCAP_ 10 (2022) 064 [2205.13487].
  * [110] D.G. Figueroa, J. Lizarraga, A. Urio and J. Urrestilla, _Strong Backreaction Regime in Axion Inflation_ , _Phys. Rev. Lett._ 131 (2023) 151003 [2303.17436].
  * [111] A. Caravano, K. Inomata and S. Renaux-Petel, _The Inflationary Butterfly Effect: Non-Perturbative Dynamics From Small-Scale Features_ , 2403.12811.
  * [112] P. Adshead, J.T. Giblin, R. Grutkoski and Z.J. Weiner, _Gauge preheating with full general relativity_ , _JCAP_ 03 (2024) 017 [2311.01504].
  * [113] P. Adshead, J.T. Giblin and Z.J. Weiner, _Gravitational waves from gauge preheating_ , _Phys. Rev. D_ 98 (2018) 043525 [1805.04550].
  * [114] 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 (2020) 171301 [1909.12843].
  * [115] P. Adshead, J.T. Giblin, M. Pieroni and Z.J. Weiner, _Constraining axion inflation with gravitational waves from preheating_ , _Phys. Rev. D_ 101 (2020) 083534 [1909.12842].
  * [116] R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198 [1311.0472].
  * [117] S. Ferrara, R. Kallosh, A. Linde and M. Porrati, _Minimal Supergravity Models of Inflation_ , _Phys. Rev. D_ 88 (2013) 085038 [1307.7696].
  * [118] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, _Large scale structure of the universe and cosmological perturbation theory_ , _Phys. Rept._ 367 (2002) 1 [astro-ph/0112551].
  * [119] V. Domcke, V. Guidetti, Y. Welling and A. Westphal, _Resonant backreaction in axion inflation_ , _JCAP_ 09 (2020) 009 [2002.02952].
  * [120] M. Münchmeyer and K.M. Smith, _Higher N-point function data analysis techniques for heavy particle production and WMAP results_ , _Phys. Rev. D_ 100 (2019) 123511 [1910.00596].
  * [121] O.H.E. Philcox, S. Kumar and J.C. Hill, _Too Hot to Handle: Searching for Inflationary Particle Production in Planck Data_ , 2405.03738.
  * [122] W.R. Coulton, O.H.E. Philcox and F. Villaescusa-Navarro, _The Impact of Non-Gaussian Primordial Tails on Cosmological Observables_ , 2406.15546.
  * [123] J.M. Ezquiaga, J. García-Bellido and V. Vennin, _The exponential tail of inflationary fluctuations: consequences for primordial black holes_ , _JCAP_ 03 (2020) 029 [1912.05399].
  * [124] D.G. Figueroa, S. Raatikainen, S. Rasanen and E. Tomberg, _Non-Gaussian Tail of the Curvature Perturbation in Stochastic Ultraslow-Roll Inflation: Implications for Primordial Black Hole Production_ , _Phys. Rev. Lett._ 127 (2021) 101302 [2012.06551].
  * [125] N. Kitajima, Y. Tada, S. Yokoyama and C.-M. Yoo, _Primordial black holes in peak theory with a non-Gaussian tail_ , _JCAP_ 10 (2021) 053 [2109.00791].
  * [126] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, _The formation probability of primordial black holes_ , _Phys. Lett. B_ 820 (2021) 136602 [2105.07810].
  * [127] J.T. Deskins, J.T. Giblin and R.R. Caldwell, _Gauge field preheating at the end of inflation_ , _Physical Review D_ 88 (2013) .
  * [128] P. Adshead, J.T.G. Jr., T.R. Scully and E.I. Sfakianakis, _Gauge-preheating and the end of axion inflation_ , _Journal of Cosmology and Astroparticle Physics_ 2015 (2015) 034–034.

##  Appendix A Analytic computations of the axion perturbations

This appendix presents the analytic computation of the power spectrum of the
axion field perturbations sourced by the gauge field. The method follows that
of Ref. [31], which studied the power spectrum of the inflaton field (sourced
by the axion) but not that of the axion.

Decomposing the axion field into a classical homogeneous value plus a quantum
operator

$\sigma\left(\tau,\,\vec{x}\right)=\sigma\left(t\right)+\delta\sigma\left(\tau,\,\vec{x}\right)=\sigma\left(t\right)+\int\frac{d^{3}k}{\left(2\pi\right)^{3/2}}\,{\rm
e}^{i\vec{k}\cdot\vec{x}}\,\frac{Q_{\sigma}\left(\tau,\,\vec{k}\right)}{a\left(\tau\right)},$
(A.1)

and analogously for the inflaton, the operator
$Q_{\sigma}\left(\tau,\,\vec{k}\right)$ satisfies the second of eqs. (2.3).
Disreegarding the backreaction from the inflaton perturbations (namely,
neglecing the term proportional to $Q_{\phi}$ in that equation) and to leading
order in slow roll, this equation becomes (cf. eq. (3.10) of [31] for details)

$\left(\frac{\partial^{2}}{\partial\tau^{2}}+k^{2}-\frac{2}{\tau^{2}}\right)Q_{\sigma}\left(\tau,\,\vec{k}\right)\simeq\frac{\alpha
a^{3}}{f}\int\frac{d^{3}x}{\left(2\pi\right)^{3/2}}{\rm
e}^{-i\vec{k}\cdot\vec{x}}\vec{E}\cdot\vec{B}\equiv
S_{\sigma}\left(\tau,\,\vec{k}\right).$ (A.2)

where the explicit form of the source can be found in eq. (C.2) of [31].

The axion field in turn sources inflaton perturbations through the first of
eqs. (2.3). Ref. [31] computed the power spectrum and the bispectrum of the
sourced inflaton perturbatons, and fitted their frequency-dependence amplitude
as a function of the model parameter $\xi_{*}$ (for the two choices
$\delta=0.2,\,0.5$). Here we present the results for the computation of the
axion power spectrum, which was not given in [31]. We are interested in the
spectrum at super-horizon scales, which we parametrize as

$\lim_{k\tau\to
0}\left\langle\delta\sigma\left(\tau,\,\vec{k}\right)\delta\sigma\left(\tau,\,\vec{k}^{\prime}\right)\right\rangle\equiv
f^{2}\frac{2\pi^{2}}{k^{3}}{\cal
P}_{\delta\sigma/f}\left(k\right)\delta^{(3)}\left(\vec{k}+\vec{k}^{\prime}\right).$
(A.3)

Eq. (A.2) is formally solved by

$Q_{\sigma,{\rm
sources}}\left(\tau,\,\vec{k}\right)=\int^{\tau}G_{k}\left(\tau,\,\tau^{\prime}\right)S_{\sigma}\left(\tau,\,\vec{k}\right),$
(A.4)

where the (de Sitter) Green function in the super-horizon limit is

$\lim_{k\tau\to
0}G_{k}\left(\tau,\,\tau^{\prime}\right)=\frac{k\tau^{\prime}\,\cos\left(k\tau^{\prime}\right)-\sin\left(k\tau^{\prime}\right)}{k^{3}\tau\tau^{\prime}}\theta\left(\tau-\tau^{\prime}\right).$
(A.5)

Inserting the solution (A.4) into eq. A.3, and evaluating the expectation
value as done in appendix C of [31] for the power spectrum of the inflaton, we
arrive to

$\displaystyle{\cal P}_{\delta\sigma/f}\left(k\right)=$
$\displaystyle\frac{\alpha^{2}H^{4}}{16\pi^{2}f^{4}}\int\frac{d^{3}{\tilde{p}}}{\left(2\pi\right)^{3}}\left|\epsilon_{i}^{(+)}\left(\hat{\tilde{p}}\right)\epsilon_{i}^{(+)}\left(\frac{{\hat{k}}-\vec{\tilde{p}}}{\left|{\hat{k}}-\vec{\tilde{p}}\right|}\right)\right|^{2}{\tilde{p}}^{1/2}\left|{\hat{k}}-\vec{\tilde{p}}\right|^{1/2}\left({\tilde{p}}^{1/2}+\left|{\hat{k}}-\vec{\tilde{p}}\right|^{1/2}\right)^{2}$
$\displaystyle
N^{2}\left[\xi_{*},\,{\tilde{p}}x_{*},\,\delta\right]N^{2}\left[\xi_{*},\,\left|{\hat{k}}-\vec{\tilde{p}}\right|x_{*},\,\delta\right]{\cal
T}_{\sigma}^{2}\left[\xi_{*},,x_{*},\,\delta,\,\sqrt{\tilde{p}}+\sqrt{\left|{\hat{k}}-\vec{\tilde{p}}\right|}\right]\;,$
(A.6)

where we have introduced the rescaled internal momentum
$\vec{\tilde{p}}\equiv\frac{\vec{p}}{k}$ and the rescaled time
$x_{*}\equiv-k\tau_{*}$ at which the axion roll is fastest. We recall that
$N[\xi_{*},k,\delta]$ is the normalizaton of the gauge mode functions defined
in eq. (2.6), while

${\cal
T}_{\sigma}\left[\xi_{*},\,x_{*},\,\delta,\,Q\right]\equiv\int_{0}^{\infty}dx\left[\sin
x-x\,\cos x\right]{\rm
exp}\left[-\frac{4\xi_{*}^{1/2}}{1+\delta}\,\frac{x^{\left(1+\delta\right)/2}}{x_{*}^{\delta/2}}\,Q\right]\;,$
(A.7)

arises from the time integration of the Green function times the gauge mode
function. Finally, $\vec{\epsilon}^{(+)}$ is the polarization operator of the
enhanced gauge mode, and the product entering in this expression is given
after eq. (C.17) of [31].

One angular integration (namely, the rotation of $\vec{p}$ around $\vec{k}$)
is trivial; for the remaining two integrations, we introduce the variables
$x\equiv{\tilde{p}}+\left|{\hat{k}}-\vec{\tilde{p}}\right|$ and
$y\equiv{\tilde{p}}-\left|{\hat{k}}-\vec{\tilde{p}}\right|$, and note that the
integrand is even in $y$, so that

$\displaystyle{\cal P}_{\delta\sigma/f}\left(k\right)=$
$\displaystyle\frac{\alpha^{2}H^{4}}{256\pi^{4}f^{4}}\int_{1}^{\infty}dx\int_{-1}^{1}dy\,\frac{\left(x^{2}-1\right)^{2}\left(\sqrt{x+y}+\sqrt{x-y}\right)^{2}}{\sqrt{x+y}\,\sqrt{x-y}}$
$\displaystyle
N^{2}\left[\xi_{*},\,\frac{x+y}{2}\,x_{*},\,\delta\right]N^{2}\left[\xi_{*},\,\frac{x-y}{2}\,x_{*},\,\delta\right]{\cal
T}_{\sigma}^{2}\left[\xi_{*},,x_{*},\,\delta,\,\frac{\sqrt{x+y}+\sqrt{x-y}}{\sqrt{2}}\right]\;,$
(A.8)

As in [31], we consider the two values $\delta=0.2,\,0.5$. We then consider
the set of values $\left\\{3,\,3.5,\,4,\,\dots,\,7\right\\}$ for $\xi_{*}$.
For each choice of these two parameters we evaluated the rescaled power
spectrum

$\tilde{\cal
P}_{\delta\sigma/f}\left(\frac{k}{k_{*}},\,\xi_{*},\,\delta\right)\equiv{\cal
P}_{\delta\sigma/f}\left(k\right)\Big{/}\frac{\alpha^{2}H^{4}}{f^{4}}\,.$
(A.9)

where $k_{*}\equiv a\left(\tau_{*}\right)H\left(\tau_{*}\right)$ is the
wavenumber of the mode that exited the horizon at the time $\tau_{*}$, when
the roll of the axion is fastest.

Figure 20: Bump of the (rescaled) power spectrum of the source perturbations
of the axion field, for $\xi_{*}=5$ and for $\delta=0.2$ (left panel) and
$\delta=0.5$. The exact integration of (A.8) is compared with the fitting
function (A.10).

In Figure 20 we show the rescaled power for $\xi_{*}=5$ and for the two
choices $\delta=0.2$ (left panel) and $\delta=0.5$ (right panel). The exact
numerical integration of eq. (A.8), shown with green dots, is compared with
the fitting function

$\tilde{\cal
P}_{\delta\sigma/f}\left(\frac{k}{k_{*}},\,\xi_{*},\,\delta\right)\simeq
f_{2,\sigma}^{c}\left[\xi_{*},\,\delta\right]\,{\rm
exp}\left[-\frac{1}{2\sigma_{2,\sigma}^{2}\left[\xi_{*},\,\delta\right]}\,\ln^{2}\left(\frac{k}{k_{*}\,x_{2,\sigma}^{c}\left[\xi_{*},\,\delta\right]}\right)\right]\,,$
(A.10)

shown with a red dashed line in the figure. We verifued that the fitting
function provides and equally good fit of the exct result for all the other
values of $\xi_{*}$ that we have considered. The function (A.10) is a log-
normal bump centered at wavenumbers close to $k_{*}$, and characterized by
three parameters $f_{2,\sigma}^{c},\,x_{2,\sigma}^{c}$, and
$\sigma_{2,\sigma}^{2}$, which control, respectively, the height, the central
position, and the width of the bump.

Figure 21: Dependence on $\xi_{*}$ of the fitting parameters $\ln
f_{2,\sigma}^{c},\,x_{2,\sigma}^{c}$, and $\sigma_{2,\sigma}^{2}$ entering in
eq. (A.10) for the choice $\delta=0.2$. The red dots denotes the values
obtained by evaluating and fitting the power spectrum at fixed $\xi_{*}$. The
solid lines are the polynomial fits reported in Table 1.

The three fitting parameters are in turn functions of the model parameters
$\xi$ and $\delta$. In Figure 21 we show the value of the parameters for
$\delta=0.2$ and for the varoius $\xi_{*}$ that we have considered. The
behaviour of these values with $\xi_{*}$ is compared with a quadratic
polynomial fit. The fit works equally well also for the other choice
$\delta=0.5$ that we have studied. We provide the fitting functions in Table
1.

| $\ln f_{2,\sigma}^{c}$ | $x_{2,\sigma}^{c}$ | $\sigma_{2,\sigma}$
---|---|---|---
$\delta=0.2$ | $-19.42+9.313\,\xi_{*}+0.1349\,\xi_{*}^{2}$ | $5.097+0.047\,\xi_{*}+0.04193\,\xi_{*}^{2}$ | $1.428-0.1785\,\xi_{*}+0.01096\,\xi_{*}^{2}$
$\delta=0.5$ | $-18.66+8.6457\,\xi_{*}+0.08085\,\xi_{*}^{2}$ | $1.666+0.67\,\xi_{*}+0.00332\,\xi_{*}^{2}$ | $0.732-0.0672\,\xi_{*}+0.00428\,\xi_{*}^{2}$
Table 1: Dependence on $\xi_{*}$ of the parameters in (A.10), for the two
values of $\delta$ that we have studied.

##  Appendix B Numerical convergence and lattice tests

In this appendix, we demonstrate the numerical stability of our results.

###  B.1 Resolution dependence tests

In this section, we show that our results do not depend on the resolution by
running simulations with different numbers of points $N_{\rm pts}$ and
comoving length $L$, corresponding to different IR and UV resolutions.

In fig. 22 we show the final power spectrum of $\sigma$ in the case of weak
backreaction for different lattice resolutions. Similarly, fig. 23 displays
the power spectrum of $\zeta$ in the case of strong backreaction. In both
cases, we see that the size of the power spectrum does not change with the
resolution, indicating that removing IR and UV modes does not affect our
results.

In fig. 24, we also plot the evolution of $\xi$ in the case of strong
backreaction, showing that the amount of backreaction does not depend on the
resolution.

$\mathcal{P}_{\sigma/f}$ $k/k_{*}$ Figure 22: Final time power spectrum of
$\sigma/f$ in the case of small backreaction ($\xi=5$). Different colors
correspond to different resolutions, as given by the legend.
$\mathcal{P}_{\zeta}$ $k/k_{*}$ Figure 23: Final time power spectrum of
$\zeta$ in the case of strong backreaction ($\xi=6$). Different colors
correspond to different resolutions, as given by the legend. $N$
$\langle\xi\rangle$ Figure 24: Time evolution of $\langle\xi\rangle$ in the
case of strong backreaction from simulations with different resolutions, as
indicated by the legend.

###  B.2 Gauge constraint

As already mentioned in section 3.1, the Lorenz constraint
$\partial^{\mu}A_{\mu}=0$ is not automatically enforced, as we evolve all four
components of the gauge field independently. While an exact solution of the
system would respect this constraint, our numerical integration unavoidably
introduces small numerical errors that break it. Therefore, we need to check
that this breaking of the Lorenz gauge is sufficiently small all throughout
the evolution, to ensure that we do not evolve any appreciable spurious
effect.

To this aim, we introduce the following dimensionless quantity [127, 128,
108]:

$\mathcal{G}\equiv\frac{\partial^{\mu}A_{\mu}}{\sqrt{\sum_{\rho}|\partial^{\rho}A_{\rho}|^{2}}},$
(B.1)

which quantifies the violation of the Lorenz constraint in a dimensionless
way. In fig. 25, we show the evolution of this quantity in the cases of strong
and weak backreaction. From this plot, we see that the Lorenz constraint is
preserved in both cases, and the size of the violation is consistent with what
is observed in the minimal model [107, 108].

$\mathcal{G}$ $N$ Figure 25: Time evolution of $\mathcal{G}$, defined in eq.
B.1, quantifying violation of the Lorenz constraint.