# Lattice Simulations of Axion-U(1) Inflation

Angelo Caravano  Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-
Maximilians Universität, Scheinerstr. 1, 81679 München, Germany  Max-Planck-
Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, 80805
Munich, Germany  Eiichiro Komatsu  Max Planck Institute for Astrophysics, Karl
Schwarzschild Str. 1, Garching, 85741, Germany  Kavli Institute for the
Physics and Mathematics of the Universe (Kavli IPMU, WPI), University of
Tokyo, Chiba 277-8582, Japan  Kaloian D. Lozanov  Illinois Center for Advanced
Studies of the Universe & Department of Physics, University of Illinois at
Urbana-Champaign, Urbana, IL 61801, USA.  Jochen Weller  Universitäts-
Sternwarte München, Fakultät für Physik, Ludwig-Maximilians Universität,
Scheinerstr. 1, 81679 München, Germany  Max Planck Institute for
Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany

###### Abstract

We present the first nonlinear lattice simulation of an axion field coupled to
a U(1) gauge field during inflation. We use it to fully characterize the
statistics of the primordial curvature perturbation $\zeta$. We find high-
order statistics to be essential in describing non-Gaussianity of $\zeta$ in
the linear regime of the theory. On the contrary, non-Gaussianity is
suppressed when the dynamics becomes nonlinear. This relaxes bounds from
overproduction of primordial black holes, allowing for an observable
gravitational waves signal at pulsar timing array and interferometers scales.
Our work establishes lattice simulations as a crucial tool to study the
inflationary epoch and its predictions.

##  I Introduction

In the standard paradigm, cosmological inflation, the accelerated expansion of
the very early universe, is driven by a scalar degree of freedom, the so-
called “inflaton” field $\phi$ Guth (1981); Sato (1981); Linde (1987);
Albrecht and Steinhardt (1982); Starobinsky (1980). The quantum vacuum
fluctuations of $\phi$ provide a natural mechanism to generate the observed
anisotropies in the cosmic microwave background (CMB) Starobinsky (1979);
Mukhanov and Chibisov (1981); Hawking (1982); Guth and Pi (1982); Starobinsky
(1982); Abbott and Wise (1984). Although the simplest single-field slow-roll
scenario is compatible with all current observations Komatsu and Bennett
(2014); Akrami  _et al._ (2020), we still lack a complete theoretical
understanding of the inflationary universe. For this reason, non-minimal
models of inflaton have been studied in the literature, involving multiple
scalars or gauge fields. In many interesting cases, the dynamics is affected
by nonlinear physics, invalidating the standard perturbative technquiques.
Nonlinear lattice simulations, developed to study the reheating epoch after
inflation Khlebnikov and Tkachev (1996); Prokopec and Roos (1997); Felder and
Tkachev (2008); Frolov (2008); Huang (2011); Sainio (2012); Child  _et al._
(2013); Easther  _et al._ (2010); Lozanov and Amin (2020); Figueroa  _et al._
(2023), might be an essential tool to compute predictions from non-minimal
inflationary scenarios.

We consider a model where inflation is driven by a pseudoscalar “axionlike”
$\phi$, which is coupled to gauge fields through Chern-Simons interaction
$\phi F\tilde{F}$ Anber and Sorbo (2006, 2010); Barnaby and Peloso (2011);
Barnaby  _et al._ (2011); Anber and Sorbo (2012); Maleknejad and Sheikh-
Jabbari (2011, 2013); Adshead and Wyman (2012); Adshead  _et al._ (2013);
Maleknejad (2021). This system gives rise to unique observational signatures,
like non-Gaussianities and chiral gravitational waves, which might be observed
with next generation experiments Komatsu (2022); Campeti  _et al._ (2021).
However, this system is often characterized by strong backreaction effects,
associated with a breakdown of perturbation theory Ferreira  _et al._ (2016);
Peloso  _et al._ (2016); Papageorgiou  _et al._ (2018); Maleknejad and Komatsu
(2019); Lozanov  _et al._ (2019); Mirzagholi  _et al._ (2020); Papageorgiou
_et al._ (2019); Ishiwata  _et al._ (2022). Although significant effort has
been put into simulating axion-gauge models during the reheating phase of the
universe Figueroa and Shaposhnikov (2018); Cuissa and Figueroa (2019);
Figueroa  _et al._ (2019, 2023); Deskins  _et al._ (2013); Adshead  _et al._
(2015, 2016, 2018a); Braden  _et al._ (2010); Lozanov and Amin (2020); Adshead
_et al._ (2020a, b); Weiner  _et al._ (2021), they have never been simulated
during the inflationary epoch.

In this paper, we use a lattice simulation to study the following axion-gauge
system during inflation:

$\displaystyle\begin{split}S=\int d^{4}x\sqrt{-g}\Biggl{[}&\frac{M^{2}_{\rm
Pl}}{2}R-\frac{(\partial_{\mu}\phi)^{2}}{2}-V(\phi)\\\
&-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{\alpha}{4f}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}\Biggr{]},\end{split}$ (1)

where $\phi$ is the pseudoscalar inflaton field, coupled to a U(1) gauge field
$A_{\mu}$ with strength tensor
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and
$\tilde{F}^{\mu\nu}=(2\sqrt{-g})^{-1}\epsilon^{\mu\nu\mu^{\prime}\nu^{\prime}}F_{\mu^{\prime}\nu^{\prime}}$.
Here, $\epsilon^{\mu\nu\mu^{\prime}\nu^{\prime}}$ is a totally antisymmetric
symbol with $\epsilon^{0123}=1$, $g$ the determinant of the metric tensor, and
$R$ the Ricci scalar, $\alpha$ the dimensionless coupling constant, and $f$
the axion decay constant. $M_{\rm Pl}$ is the reduced Planck mass, that we set
to 1 throughout this paper.

The Chern-Simons interaction $\phi F\tilde{F}$ leads to an abundant production
of gauge field particles Anber and Sorbo (2006, 2010), which act as a source
for inflaton perturbations $\delta\phi$ via the inverse particle decay Anber
and Sorbo (2010); Barnaby and Peloso (2011); Barnaby  _et al._ (2011); Anber
and Sorbo (2012), hence affecting the statistical properties of the comoving
curvature perturbation $\zeta$. The power spectrum of $\zeta$ on super-horizon
scales $k\ll aH$ has been estimated analytically as Anber and Sorbo (2010);
Barnaby and Peloso (2011); Barnaby  _et al._ (2011); Anber and Sorbo (2012):

$\mathcal{P}_{\zeta}(k)\simeq\mathcal{P}_{\rm vac}+\mathcal{P}_{\rm
vac}^{2}f_{2}(\xi)e^{4\pi\xi},\quad\quad\xi=\frac{\alpha\dot{\phi}}{2fH},$ (2)

where $\mathcal{P}_{\rm vac}=H^{4}/(2\pi\dot{{\phi}})^{2}$ is the vacuum
contribution and $f_{2}(\xi)$ is a function that can be found in Ref. Barnaby
_et al._ (2011). This result was derived assuming a constant $\xi$. Moreover,
statistics of $\zeta$ are highly non-Gaussian, as $F\tilde{F}$ is bilinear in
the gauge field. The bispectrum has also been estimated in the $\xi-$constant
approximation Anber and Sorbo (2010); Barnaby and Peloso (2011); Barnaby  _et
al._ (2011); Anber and Sorbo (2012). We avoid reporting the lengthy expression
for the bispectrum, whose value can be found, for example, in Ref. Barnaby
_et al._ (2011).

These analytical estimates are valid when at least two key assumptions are
satisfied. First, the backreaction of $\phi F\tilde{F}$ on the background
inflationary trajectory is small. This translates into the following bound
Anber and Sorbo (2010); Barnaby and Peloso (2011); Barnaby  _et al._ (2011);
Anber and Sorbo (2012):

$\frac{H^{2}}{26\pi|\dot{\phi}|}\xi^{-3/2}e^{\pi\xi}\ll 1.$ (3)

Second, the power spectrum $\mathcal{P}_{\zeta}(k)$ remains small (typically
smaller than $10^{-1}$), to ensure perturbativity. A violation of these
assumptions invalidates the perturbation theory approach and requires
nonlinear tools, such as presented in this paper.

##  II Lattice simulation

Our simulation is based on the methodology developed in Ref. Caravano  _et
al._ (2021a, b), to which we refer for details. We discretize the classical
equations of motion in real space. We choose to work in the Lorenz gauge
$\partial^{\mu}A_{\mu}=0$, in which the equations read Adshead  _et al._
(2015):

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

where $i,j,k\in\\{1,2,3\\}$ and the prime denotes derivatives with respect to
conformal time. The scale factor in eq. 4 is evolved self-consistently with
the second Friedmann equation. As commonly done in the literature Anber and
Sorbo (2010); Barnaby and Peloso (2011); Barnaby  _et al._ (2011, 2012); Anber
and Sorbo (2012); Linde  _et al._ (2013); Domcke  _et al._ (2020); Gorbar  _et
al._ (2021), we neglect the role of metric perturbations because gravitational
interactions are slow-roll suppressed during inflation. Moreover, the field
dynamics of models characterized by large non-Gaussianity is expected to be
decoupled from the gravitational sector Leblond and Pajer (2011).

To solve this system of equations, we associate field values
$\phi_{n_{1},n_{2},n_{3}}$ and ${A}_{\mu,n_{1},n_{2},n_{3}}$ to the $N^{3}$
points of a periodic cubic lattice with comoving volume $L^{3}$. After
defining a discretization scheme for the spatial derivatives, eq. 4 constitute
a set of second order coupled differential equations that we solve numerically
with a Runge-Kutta 4th order integrator. We start the simulation when the
lattice box size satisfies $L\lesssim 1/(aH)$, so that the fields are
approximately in their Bunch-Davies vacuum state at the beginning of the
simulation. The gauge condition $\partial^{\mu}A_{\mu}=0$ is only imposed at
the initial time. Therefore, we need to check by hand that
$\partial^{\mu}A_{\mu}$ vanishes with sufficient precision throughout the
evolution. We find that the dimensionless gauge constraint
${\partial^{\mu}A_{\mu}}/{\sqrt{\sum_{\rho}|\partial^{\rho}A_{\rho}|^{2}}}$ is
always smaller than $3\times 10^{-4}$ for all the simulation runs shown below.

##  III Negligible backreaction

We show the results of the simulation starting from the case when backreaction
is negligible (see eq. 3), and compare them to the known analytical results.
We assume a monodromy potential for the inflaton McAllister  _et al._ (2014)
$V(\phi)=\frac{1}{2}m^{2}\phi^{2}$ with $m=0.51\cdot 10^{-5}$. The system is
initiated far from the end of inflation by setting ${\phi}=-14.5$. We run a
simulation with $N^{3}=256^{3}$ points and comoving size $L=2/m$. We evolve
the system for $N_{e}=6$ e-folds, which makes the simulation box satisfy $L\gg
1/(aH)$ at the end of the simulation. For this run, we set the gauge coupling
$\alpha/f=42$, which is excluded by CMB observations but allows us to better
compare the results of the simulation with the existing analytical estimates.
Below we consider a more realistic value of the coupling. In the left panel of
fig. 1 we show the value of $\xi$ during this simulation, which monotonically
grows following the slow-roll trajectory.

Figure 1: Time evolution of the $\xi$ parameter as a function of the number of
e-folds $N_{e}$, in the case of negligible backreaction (left) and strong
backreaction (right).

In the upper panel of fig. 2 we show the power spectrum of the comoving
curvature perturbation $\zeta\equiv-\delta\phi H/\dot{\phi}$ at different
times during the simulation111This relation for the curvature perturbation is
valid as long as the energy density of the Universe is dominated by the
background $\phi$ field, which remains true for all the cases considered in
this work.. We compare the final power spectrum with eq. 2, which is shown as
a shaded region as $\xi$ varies during the evolution. The black dashed lines
delimiting this region are computed using the initial and final values of
$\xi$. In the bottom panel of fig. 2 we show the bispectrum
$\mathcal{B}_{\zeta}(k)\equiv\langle\zeta({\vec{k}_{1}})\zeta({\vec{k}_{2}})\zeta^{*}({\vec{k}_{1}}+{\vec{k}_{2}})\rangle$
on equilateral configurations
$k\equiv|\vec{k}_{1}|=|\vec{k}_{2}|=|\vec{k}_{1}+\vec{k}_{2}|$ at the final
time, and compare it to the analytical estimate of Ref. Barnaby  _et al._
(2011). We find that both the bispectrum and the power spectrum are in
agreement with the analytical estimates. Note that for the largest modes there
is a drop in the lattice spectra, which is unphysical and it is caused by the
lattice UV cutoff.

Figure 2: (Top) Power spectrum of $\zeta$ in the case of weak backreaction.
The shaded region, delimited by black dashed lines, shows the analytical
prediction of eq. 2. The blue dashed line shows the vacuum contribution
$\mathcal{P}_{\rm vac}$. (Bottom) Equilateral bispectrum of $\zeta$ compared
to the analytical prediction.

Thanks to the lattice approach, we have access to the curvature perturbation
in real space. In the left panel of fig. 3 we show the normalized histograms
of the values of $\zeta$ across the $N^{3}$ points at different times during
the simulation.

Figure 3: (Left) Normalized histograms of $\zeta$ in real space in the case of
weak backreaction. (Right) Time evolution of the correlators defined in eq. 5.
For a more detailed version of this figure, with associated errors, see Ref.
Caravano (2022).

We find that non-Gaussianity manifests as a pronounced exponential tail in the
distribution of $\zeta$. To quantify non-Gaussianity, we compute the cumulants
of the one-point probability density function Bernardeau  _et al._ (2002):

$\kappa_{3}=\frac{\langle\zeta^{3}\rangle}{\sigma^{3}},\quad\kappa_{4}=\frac{\langle\zeta^{4}\rangle-3\sigma^{4}}{\sigma^{4}},\quad\kappa_{5}=\frac{\langle\zeta^{5}\rangle-10\langle\zeta^{3}\rangle\sigma^{2}}{\sigma^{5}},$
(5)

which we normalized by powers of $\sigma^{2}=\langle\zeta^{2}\rangle$ to make
them dimensionless. In the right panel of fig. 3 we show the evolution of the
cumulants during the simulation. We find that
$\kappa_{6}>\kappa_{5}>\kappa_{4}>\kappa_{3}$ at late times. For illustrative
purpose, we avoid showing the evolution of $\kappa_{6}$, whose final value is
$\kappa_{6}\simeq 10.3$. This result shows that higher-order statistics are
essential to characterize non-Gaussianity of $\zeta$. This has important
observational consequences, as discussed below.

##  IV Strong backreaction

We now turn to the case of strong backreaction. We set the gauge coupling to
$\alpha/f=25$, so that the imprints of the Chern-Simons coupling on $\zeta$
are unobservable at CMB scales Anber and Sorbo (2010); Barnaby and Peloso
(2011); Barnaby  _et al._ (2011); Anber and Sorbo (2012). Later during
inflation, however, $\xi$ increases and the universe eventually enters a
nonlinear phase.

We start the simulation when $\phi=-5.5$. With this choice, the universe is
still in the weak backreaction phase at the beginning of the simulation. Then,
after roughly 2 e-folds, the system enters a strong backreaction phase where
the bound of eq. 3 is violated and eq. 2 gives $\mathcal{P}_{\zeta}\sim 0.1$,
which indicates a breakdown of perturbativity. We show results from a run with
$(N,L)=(256,1.5/m)$, but we tested our simulation also with other values of
$(N,L)$ to ensure that our results are physical and do not depend on the
spatial resolution. Moreover, we ensured the stability of the time integration
by checking energy conservation and time-step convergence.

In the right panel of fig. 1 we show the evolution of $\xi$ during the
simulation. We find the departure from the slow-roll trajectory as an
oscillatory behavior in $\xi$. This is intuitive, as one can see from eq. 4
that a strong $F\tilde{F}$ leads to a depletion of the inflaton velocity; this
lowers the value of $\xi$ and reduces the backreaction, bringing the system
momentarily back to the slow-roll trajectory. Oscillations of similar period
and size were already predicted by previous studies Cheng  _et al._ (2016);
Notari and Tywoniuk (2016); Dall’Agata  _et al._ (2020); Domcke  _et al._
(2020), which explored backreaction effects using semi-analytical tools.
Another consequence of the backreaction is that, after $6.5$ e-folds of
evolution, the background inflaton value is $\phi=-3.02$. This value would be
reached after $5.4$ e-folds of evolution if the backreaction were negligible,
which means that the backreaction significantly delays the background
dynamics.

In fig. 4 we show the histograms of $\zeta$ and the evolution of the cumulants
$\kappa_{i}$. These plots show that the non-Gaussianity of $\zeta$
substantially decreases during the strong backreaction phase. At late times,
it is mainly described by a (small) negative $\kappa_{4}$, while the other
cumulants are negligible. Moreover, $\kappa_{5}$ shows oscillations. The
suppression of non-Gaussianity in this regime is a consequence of the central
limit theorem, and it is caused by the fact that the number of excited gauge
field modes grows with $\xi$. To understand this, we expand the source term
$F\tilde{F}$ in Fourier space as follows:

$\left(F_{\mu\nu}\tilde{F}^{\mu\nu}\right)(k)=\sum_{k^{\prime}}F_{\mu\nu}(k^{\prime})\,\,\tilde{F}^{\mu\nu}(k-k^{\prime}).$
(6)

This shows that each Fourier mode of $F\tilde{F}$ is the sum of several non-
Gaussian quantities. For $\xi\sim 1$, there are few elements contributing to
this sum due to the small number of excited gauge field modes. For $\xi\gg 1$,
the number of statistically independent elements in this sum is large, and
$F\tilde{F}$ converges to a Gaussian distribution due to the central limit
theorem. Therefore, it sources a Gaussian $\zeta$. This is analogous to what
happens when $\phi$ is coupled to fermionic fields Adshead  _et al._ (2018b).

Parametrizing non-Gaussianity in the local type form Komatsu and Spergel
(2001):

$\zeta(\vec{x})=\zeta_{g}(\vec{x})+f_{\rm
NL}[\zeta_{g}^{2}(\vec{x})-\langle\zeta^{2}_{g}(\vec{x})\rangle],$ (7)

where $\zeta_{g}$ is a Gaussian field and $f_{\rm NL}$ a real number, one can
use the linear results of Ref. Barnaby and Peloso (2011) to show that $f_{\rm
NL}\propto e^{-2\pi\xi}$. This shows that the suppression of non-Gaussianity
for large $\xi$ can be guessed from the linear regime of the theory in the
$\xi$-constant approximation. Although the linear theory is not reliable in
this regime, and we find that non-Gaussianity is not of the local type, the
simulation confirms this intuition.

Figure 4: Histograms of $\zeta$ (left) and the cumulants $\kappa_{i}$ (right)
in the case of strong backreaction.

In fig. 5 we show the power spectrum from the simulation. Although the
analytical estimates are not reliable in this regime, we still compare it with
eq. 2 using the initial and final values of $\xi$ from the simulation. As
$\kappa_{3}\ll 1$, we are not able to compute the bispectrum because it is
below the noise of our bispectrum estimator.

Figure 5: Power spectrum from the simulation in the case of strong
backreaction.

##  V Discussion

We presented the first lattice simulation of nonlinear perturbations generated
during inflation. For the first time, we obtain the full primordial density
fluctuation predicted by inflation. We used the simulation to fully
characterize the statistics of $\zeta$ from the axion-U(1) model.

In the weak backreaction regime, power spectrum and bispectrum agree with
analytical estimates. Nevertheless, high-order information is crucial in
describing the statistics, showing that $n$-point correlators are not
efficient to fully characterize this signal. Current large-scale bounds,
constraining $\alpha/f\lesssim 32$ for the quadratic potential, are derived
using power spectrum and bispectrum of the CMB, ignoring higher-order
statistics Anber and Sorbo (2010); Barnaby and Peloso (2011); Barnaby  _et
al._ (2011); Anber and Sorbo (2012).

In light of our results, we expect the information beyond power spectrum and
bispectrum to play a key role. The output of the simulation can be used as
truly ab initio initial condition for the cosmological simulations of
structure formation, allowing to test inflation using the full information
contained in the density field. This opens a new possibility for the field of
cosmological simulations, that we plan to explore in future work.

In the case of strong backreaction, typically relevant for small scales, the
system enters an instability phase, characterized by oscillations in the
inflation background velocity. In this regime, non-Gaussianity of $\zeta$ is
strongly suppressed. This is a consequence of the central limit theorem, and
it is caused by the large number of excited gauge field modes contributing to
the source term $F\tilde{F}$. We expect this suppression to be a general
feature of models where matter fields are coupled linearly to the inflaton
$\mathcal{L}\supset\phi f(X)$, with $f(X)$ being a quadratic function of a
generic matter field $X$, that could be for example a scalar $X=\psi$ or a
gauge field $X=A^{a}_{\mu}$. If $X$ is copiously produced during inflation via
some mechanism, its contribution to the statistics of $\zeta$ is expected to
be Gaussian for the same reason. Due to the simplicity of this argument, the
same conclusion could hold for more generic nonlinear functions $f(X)$,
although this needs to be investigated in future studies.

This result relaxes the bounds from the overclosure of the universe due to
overproduction of primordial black holes (PBH), constraining $\alpha/f\lesssim
23$ for the quadratic potential Linde  _et al._ (2013); Garcia-Bellido  _et
al._ (2016). These bounds strongly rely on the assumption that $\zeta$ can be
approximated by a (non-Gaussian) $\chi^{2}$ distribution during the strong
backreaction phase Linde  _et al._ (2013); Garcia-Bellido  _et al._ (2016),
which corresponds to the $f_{\rm NL}\gg 1$ limit of eq. 7. We find a nearly
Gaussian $\zeta$, requiring a much larger power spectrum at small scales to
efficiently produce PBH Linde  _et al._ (2013). Estimating PBH production
requires a more extensive and detailed study of the final e-folds of
inflation, for two main reasons. First, as we show in this paper, backreaction
significantly delays the end of inflation, making it problematic to identify
the range of modes relevant for PBH production. Second, as we show in fig. 4,
there is still a small remnant non-Gaussianity at the end of the simulation,
to which the production of PBH is extremely sensitive.

We conclude that the most stringent bounds on $\alpha/f$ are the ones from the
statistics of $\zeta$ at large scales, discussed above. This allows for an
inflationary gravitational waves (GW) signal within reach of LISA Amaro-Seoane
_et al._ (2017); Bartolo  _et al._ (2016), advanced LIGO Abbott  _et al._
(2016) and PTA-SKA Foster and Backer (1990); Dewdney  _et al._ (2009); Kramer
_et al._ (2004) experiments. Indeed, the gauge field acts as a source for GW
Anber and Sorbo (2010); Barnaby and Peloso (2011); Barnaby  _et al._ (2011);
Anber and Sorbo (2012), and the signal can be above the projected sensitivity
of all these experiments in the parameter range compatible with current CMB
constraints Garcia-Bellido  _et al._ (2016).

Note that, both for simplicity and to better compare with previous studies, we
considered a quadratic potential for the inflaton
$V(\phi)=\frac{1}{2}m^{2}\phi^{2}$, which is disfavored by the latest Planck-
BICEP/Keck results Akrami  _et al._ (2020); Ade  _et al._ (2021). The
particular choice of slow-roll potential, however, only affects the
quantitative bounds on $\alpha/f$ given above. Our findings about the
statistics of $\zeta$, that are the main original result of this work, do not
depend on this choice.

The lattice simulation presented in this paper allowed to reveal unknown
aspects of the axion-U(1) system that are beyond the regime of validity of
perturbation theory. At the same time, it allowed to compute inflationary
observables within this model with a precision that exceeds state-of-the-art
analytical and semi-analytical computations, both in the weak and strong
backreaction regimes. There several other cases where a lattice simulation
could be crucial, like models involving non-abelian SU(2) gauge fields
Maleknejad and Sheikh-Jabbari (2011, 2013); Adshead and Wyman (2012); Adshead
_et al._ (2013); Maleknejad (2021), or scalar fields models with a strong turn
in field space Fumagalli  _et al._ (2021). Our work shows that lattice
simulations could be an essential tool to understand the predictions of these
models, and more generically inflationary scenarios characterized by nonlinear
physics.

###### Acknowledgements.

This work is supported in part by the Excellence Cluster ORIGINS which is
funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311 (AC,
EK, JW), and JSPS KAKENHI Grant Number JP20H05859 (EK). The Kavli IPMU is
supported by World Premier International Research Center Initiative (WPI),
MEXT, Japan. The work of KL was supported in part by the US Department of
Energy through grant DE-SC0015655.

## References

  * Guth (1981) A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * Sato (1981) K. Sato, Mon. Not. Roy. Astron. Soc. 195, 467 (1981).
  * Linde (1987) A. D. Linde, Adv. Ser. Astrophys. Cosmol. 3, 149 (1987).
  * Albrecht and Steinhardt (1982) A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982).
  * Starobinsky (1980) A. Starobinsky, Physics Letters B 91, 99 (1980).
  * Starobinsky (1979) A. A. Starobinsky, JETP Lett. 30, 682 (1979).
  * Mukhanov and Chibisov (1981) V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33, 532 (1981).
  * Hawking (1982) S. Hawking, Physics Letters B 115, 295 (1982).
  * Guth and Pi (1982) A. H. Guth and S.-Y. Pi, Phys. Rev. Lett. 49, 1110 (1982).
  * Starobinsky (1982) A. Starobinsky, Physics Letters B 117, 175 (1982).
  * Abbott and Wise (1984) L. Abbott and M. B. Wise, Nucl. Phys. B 244, 541 (1984).
  * Komatsu and Bennett (2014) E. Komatsu and C. L. Bennett (WMAP Science Team), PTEP 2014, 06B102 (2014), arXiv:1404.5415 [astro-ph.CO] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Khlebnikov and Tkachev (1996) S. Y. Khlebnikov and I. I. Tkachev, Physical Review Letters 77, 219–222 (1996).
  * Prokopec and Roos (1997) T. Prokopec and T. G. Roos, Physical Review D 55, 3768–3775 (1997).
  * Felder and Tkachev (2008) G. Felder and I. Tkachev, Computer Physics Communications 178, 929?932 (2008).
  * Frolov (2008) A. V. Frolov, Journal of Cosmology and Astroparticle Physics 2008 (11), 009.
  * Huang (2011) Z. Huang, Physical Review D 83, 10.1103/physrevd.83.123509 (2011).
  * Sainio (2012) J. Sainio, Journal of Cosmology and Astroparticle Physics 2012 (04), 038–038.
  * Child  _et al._ (2013) H. L. Child, J. T. Giblin, R. H. Ribeiro, and D. Seery, Physical Review Letters 111, 10.1103/physrevlett.111.051301 (2013).
  * Easther  _et al._ (2010) R. Easther, H. Finkel, and N. Roth, Journal of Cosmology and Astroparticle Physics 2010 (10), 025–025.
  * Lozanov and Amin (2020) K. D. Lozanov and M. A. Amin, Journal of Cosmology and Astroparticle Physics 2020 (04), 058–058.
  * Figueroa  _et al._ (2023) D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, Comput. Phys. Commun. 283, 108586 (2023), arXiv:2102.01031 [astro-ph.CO] .
  * Anber and Sorbo (2006) M. M. Anber and L. Sorbo, Journal of Cosmology and Astroparticle Physics 2006 (10), 018–018.
  * Anber and Sorbo (2010) M. M. Anber and L. Sorbo, Physical Review D 81, 10.1103/physrevd.81.043534 (2010).
  * Barnaby and Peloso (2011) N. Barnaby and M. Peloso, Physical Review Letters 106, 10.1103/physrevlett.106.181301 (2011).
  * Barnaby  _et al._ (2011) N. Barnaby, R. Namba, and M. Peloso, Journal of Cosmology and Astroparticle Physics 2011 (04), 009–009.
  * Anber and Sorbo (2012) M. M. Anber and L. Sorbo, Physical Review D 85, 10.1103/physrevd.85.123537 (2012).
  * Maleknejad and Sheikh-Jabbari (2011) A. Maleknejad and M. M. Sheikh-Jabbari, Physical Review D 84, 10.1103/physrevd.84.043515 (2011).
  * Maleknejad and Sheikh-Jabbari (2013) A. Maleknejad and M. M. Sheikh-Jabbari, Phys. Lett. B 723, 224 (2013), arXiv:1102.1513 [hep-ph] .
  * Adshead and Wyman (2012) P. Adshead and M. Wyman, Physical Review Letters 108, 10.1103/physrevlett.108.261302 (2012).
  * Adshead  _et al._ (2013) P. Adshead, E. Martinec, and M. Wyman, Journal of High Energy Physics 2013, 10.1007/jhep09(2013)087 (2013).
  * Maleknejad (2021) A. Maleknejad, Phys. Rev. D 104, 083518 (2021), arXiv:2012.11516 [hep-ph] .
  * Komatsu (2022) E. Komatsu, Nature Rev. Phys. 4, 452 (2022), arXiv:2202.13919 [astro-ph.CO] .
  * Campeti  _et al._ (2021) P. Campeti, E. Komatsu, D. Poletti, and C. Baccigalupi, JCAP 01, 012, arXiv:2007.04241 [astro-ph.CO] .
  * Ferreira  _et al._ (2016) R. Z. Ferreira, J. Ganc, J. Noreña, and M. S. Sloth, JCAP 04, 039, [Erratum: JCAP 10, E01 (2016)], arXiv:1512.06116 [astro-ph.CO] .
  * Peloso  _et al._ (2016) M. Peloso, L. Sorbo, and C. Unal, Journal of Cosmology and Astroparticle Physics 2016 (09), 001–001.
  * Papageorgiou  _et al._ (2018) A. Papageorgiou, M. Peloso, and C. Unal, JCAP 09, 030, arXiv:1806.08313 [astro-ph.CO] .
  * Maleknejad and Komatsu (2019) A. Maleknejad and E. Komatsu, Journal of High Energy Physics 2019, 10.1007/jhep05(2019)174 (2019).
  * Lozanov  _et al._ (2019) K. D. Lozanov, A. Maleknejad, and E. Komatsu, Journal of High Energy Physics 2019, 10.1007/jhep02(2019)041 (2019).
  * Mirzagholi  _et al._ (2020) L. Mirzagholi, A. Maleknejad, and K. D. Lozanov, Phys. Rev. D 101, 083528 (2020), arXiv:1905.09258 [hep-th] .
  * Papageorgiou  _et al._ (2019) A. Papageorgiou, M. Peloso, and C. Unal, JCAP 07, 004, arXiv:1904.01488 [astro-ph.CO] .
  * Ishiwata  _et al._ (2022) K. Ishiwata, E. Komatsu, and I. Obata, JCAP 03 (03), 010, arXiv:2111.14429 [hep-ph] .
  * Figueroa and Shaposhnikov (2018) D. G. Figueroa and M. Shaposhnikov, Nuclear Physics B 926, 544–569 (2018).
  * Cuissa and Figueroa (2019) J. R. C. Cuissa and D. G. Figueroa, Journal of Cosmology and Astroparticle Physics 2019 (06), 002–002.
  * Figueroa  _et al._ (2019) D. G. Figueroa, A. Florio, and M. Shaposhnikov, Journal of High Energy Physics 2019, 10.1007/jhep10(2019)142 (2019).
  * Deskins  _et al._ (2013) J. T. Deskins, J. T. Giblin, and R. R. Caldwell, Physical Review D 88, 10.1103/physrevd.88.063530 (2013).
  * Adshead  _et al._ (2015) P. Adshead, J. T. G. Jr., T. R. Scully, and E. I. Sfakianakis, Journal of Cosmology and Astroparticle Physics 2015 (12), 034–034.
  * Adshead  _et al._ (2016) P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis, Journal of Cosmology and Astroparticle Physics 2016 (10), 039–039.
  * Adshead  _et al._ (2018a) P. Adshead, J. T. Giblin, and Z. J. Weiner, Physical Review D 98, 10.1103/physrevd.98.043525 (2018a).
  * Braden  _et al._ (2010) J. Braden, L. Kofman, and N. Barnaby, Journal of Cosmology and Astroparticle Physics 2010 (07), 016–016.
  * Adshead  _et al._ (2020a) P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, Physical Review D 101, 10.1103/physrevd.101.083534 (2020a).
  * Adshead  _et al._ (2020b) P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, Physical Review Letters 124, 10.1103/physrevlett.124.171301 (2020b).
  * Weiner  _et al._ (2021) Z. J. Weiner, P. Adshead, and J. T. Giblin, Physical Review D 103, 10.1103/physrevd.103.l021301 (2021).
  * Caravano  _et al._ (2021a) A. Caravano, E. Komatsu, K. D. Lozanov, and J. Weller, Journal of Cosmology and Astroparticle Physics 2021 (12), 010.
  * Caravano  _et al._ (2021b) A. Caravano, E. Komatsu, K. D. Lozanov, and J. Weller, (2021b), arXiv:2110.10695 [astro-ph.CO] .
  * Barnaby  _et al._ (2012) N. Barnaby, E. Pajer, and M. Peloso, Phys. Rev. D 85, 023525 (2012).
  * Linde  _et al._ (2013) A. Linde, S. Mooij, and E. Pajer, Phys. Rev. D 87, 103506 (2013).
  * Domcke  _et al._ (2020) V. Domcke, V. Guidetti, Y. Welling, and A. Westphal, JCAP 09, 009, arXiv:2002.02952 [astro-ph.CO] .
  * Gorbar  _et al._ (2021) E. V. Gorbar, K. Schmitz, O. O. Sobol, and S. I. Vilchinskii, Phys. Rev. D 104, 123504 (2021), arXiv:2109.01651 [hep-ph] .
  * Leblond and Pajer (2011) L. Leblond and E. Pajer, Journal of Cosmology and Astroparticle Physics 2011 (01), 035.
  * McAllister  _et al._ (2014) L. McAllister, E. Silverstein, A. Westphal, and T. Wrase, JHEP 09, 123, arXiv:1405.3652 [hep-th] .
  * Caravano (2022) A. Caravano,  _Simulating the inflationary Universe: from single-field to the axion-U(1) model_ , Ph.D. thesis, LMU Munich (2022), arXiv:2209.13616 [astro-ph.CO] .
  * Bernardeau  _et al._ (2002) F. Bernardeau, S. Colombi, E. Gaztanaga, and R. Scoccimarro, Phys. Rept. 367, 1 (2002), arXiv:astro-ph/0112551 .
  * Cheng  _et al._ (2016) S.-L. Cheng, W. Lee, and K.-W. Ng, Phys. Rev. D 93, 063510 (2016), arXiv:1508.00251 [astro-ph.CO] .
  * Notari and Tywoniuk (2016) A. Notari and K. Tywoniuk, JCAP 12, 038, arXiv:1608.06223 [hep-th] .
  * Dall’Agata  _et al._ (2020) G. Dall’Agata, S. González-Martín, A. Papageorgiou, and M. Peloso, JCAP 08, 032, arXiv:1912.09950 [hep-th] .
  * Adshead  _et al._ (2018b) P. Adshead, L. Pearce, M. Peloso, M. A. Roberts, and L. Sorbo, JCAP 06, 020, arXiv:1803.04501 [astro-ph.CO] .
  * Komatsu and Spergel (2001) E. Komatsu and D. N. Spergel, Phys. Rev. D 63, 063002 (2001), arXiv:astro-ph/0005036 .
  * Garcia-Bellido  _et al._ (2016) J. Garcia-Bellido, M. Peloso, and C. Unal, JCAP 12, 031, arXiv:1610.03763 [astro-ph.CO] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ 10.48550/arXiv.1702.00786 (2017).
  * Bartolo  _et al._ (2016) N. Bartolo _et al._ , JCAP 12, 026, arXiv:1610.06481 [astro-ph.CO] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 131102 (2016), arXiv:1602.03847 [gr-qc] .
  * Foster and Backer (1990) R. S. Foster and D. C. Backer, Astrophys. J. 361, 300 (1990).
  * Dewdney  _et al._ (2009) P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. W. Lazio, Proceedings of the IEEE 97, 1482 (2009).
  * Kramer  _et al._ (2004) M. Kramer, D. C. Backer, J. M. Cordes, T. J. W. Lazio, B. W. Stappers, and S. Johnston, New Astron. Rev. 48, 993 (2004), arXiv:astro-ph/0409379 .
  * Ade  _et al._ (2021) P. A. R. Ade _et al._ (BICEP, Keck), Phys. Rev. Lett. 127, 151301 (2021), arXiv:2110.00483 [astro-ph.CO] .
  * Fumagalli  _et al._ (2021) J. Fumagalli, S. Renaux-Petel, and L. T. Witkowski, JCAP 08, 030, arXiv:2012.02761 [astro-ph.CO] .