Lattice formulation of axion inflation. Application to preheating

TL;DR

This paper develops a lattice formulation for axion-gauge interactions that accurately models early universe dynamics, including inflation and preheating, while preserving key symmetries and enabling detailed numerical simulations.

Contribution

It introduces a gauge-invariant, shift-symmetric lattice formulation with an implicit solver for simulating axion inflation and preheating, accounting for inhomogeneities and backreaction effects.

Findings

Identifies two regimes of axion-gauge coupling affecting energy transfer.

Demonstrates efficient reheating and rapid transition to radiation domination.

Confirms previous results on preheating dynamics.

Abstract

We present a lattice formulation of an interaction ${\phi\over\Lambda} F\tilde F$ between an axion and some $U(1)$ gauge sector with the following properties: it reproduces the continuum theory up to $\mathcal{O}(dx_\mu^2)$ corrections, it preserves exact gauge invariance and shift symmetry on the lattice, and it is suitable for self-consistent expansion of the Universe. The lattice equations of motion can no longer be solved by explicit methods, but we propose an implicit method to overcome this difficulty, which preserves the relevant system constraints down to arbitrary (tunable) precision. As a first application we study, in a comoving grid in $(3+1)$ dimensions, the last efolds of axion-inflation with quadratic potential and the preheating stage following afterwards. We fully account for the inhomogeneity and non-linearity of the system, including the gauge field contribution to the expansion rate of the Universe and its backreaction into the axion dynamics. We characterize in detail, as a function of the coupling, the energy transfer from the axion to the gauge field. Two coupling regimes are identified, sub- and super-critical, depending on whether the final energy fraction stored in the gauge field is below or above $\sim 50\%$ of the total energy. The Universe is very efficiently reheated for super-critical couplings, rapidly entering in a radiation dominated stage. Our results on preheating confirm previously published results.

Figures

Equations

$$S$$

$$F_{\mu \nu }\equiv {\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\nu }\text{, }{\tilde{F}}_{\mu \nu }\equiv \frac{1}{2}{ϵ}_{\mu \nu \rho \sigma }{F}^{\rho \sigma }\text{ , }$$

$${ϵ}^{0123}\equiv \frac{1}{\sqrt{-g}}=\frac{1}{a^3(t)}\text{ .}$$

$${\pi }_{\varphi }\equiv \dot{\varphi },E_i\equiv {\dot{A}}_i-{\partial }_i{A}_0,B_i\equiv {ϵ}_{ijk}{\partial }_j{A}_k,$$

$$S_m$$

$${\partial }_0{\partial }_0\varphi$$

$${\partial }_{\mu }\left(\sqrt{-g}{F}^{\mu \nu }\right)$$

$${\dot{\pi }}_{\varphi }=-3H{\pi }_{\varphi }+\frac{1}{a^2}{\vec{\nabla }}^2\varphi -m^2\varphi +\frac{1}{a^3\Lambda }\vec{E}\cdot \vec{B}\text{,}$$

$$\dot{\vec{E}}=-H\vec{E}-\frac{1}{a^2}\vec{\nabla }\times \vec{B}-\frac{1}{a\Lambda }{\pi }_{\varphi }\vec{B}+\frac{1}{a\Lambda }\vec{\nabla }\varphi \times \vec{E}-\frac{\varphi }{a\Lambda }\underset{\text{=0}}{\underbrace{\left(\dot{\vec{B}}-\vec{\nabla }\times \vec{E}\right)}}\text{, }$$

$$\vec{\nabla }\cdot \vec{E}=-\frac{1}{a\Lambda }\vec{\nabla }\varphi \cdot \vec{B}-\frac{\varphi }{a\Lambda }\underset{\text{=0}}{\underbrace{\vec{\nabla }\cdot \vec{B}}}\text{ , }$$

$${\partial }_{\mu }(\sqrt{-g}{\tilde{F}}^{\mu \nu })=0⟺\{\begin{array}{l}\dot{\vec{B}}-\vec{\nabla }\times \vec{E}=0\text{vspace}\ast 0.25cm\\ \vec{\nabla }\cdot \vec{B}=0\end{array}\text{ . }$$

$$\frac{\ddot{a}}{a}$$

$${\left(\frac{\dot{a}}{a}\right)}^2$$

$$\rho =T_{00},p=\frac{1}{3a^2}\sum _i{T}_{ii}.$$

$$T_{\mu \nu }$$

$$\rho$$

$$p$$

$${ϵ}_H=-\frac{\dot{H}}{H^2},$$

$$\ddot{\varphi }+3H\dot{\varphi }+m^2\varphi =0,\frac{\ddot{a}}{a}=\frac{1}{3m_{pl}^2}\left(m^2{\varphi }^2-{\dot{\varphi }}^2\right).$$

$$N\equiv -\ln \frac{a_{end}}{a}=-\frac{1}{\sqrt{6}{m}_{pl}}{\int }_t^{t_{\mathrm{end}}}d{t}^{\prime }\sqrt{{\dot{\varphi }}^2(t^{\prime })+m^2{\varphi }^2(t^{\prime })}\sim -{\left(\frac{\varphi }{2m_{pl}}\right)}^2,$$

$$\nu =0$$

$$\nu =i$$

$$\vec{A}(\tau ,\mathbf{x})=\int \frac{d^{3}k}{(2\pi {)}^3}\vec{A}(\tau ,\mathbf{k})e^{i\mathbf{k}\cdot \mathbf{x}}=\sum _{\lambda =\pm }\int \frac{d^{3}k}{(2\pi {)}^3}{A}^{\lambda }(\tau ,k){\vec{\epsilon }}^{\lambda }(\hat{\mathbf{k}})e^{i\mathbf{k}\cdot \mathbf{x}}\text{ , }$$

$$k_i{\epsilon }_i^{\pm }(\hat{\mathbf{k}})=0,{ϵ}_{ijk}{k}_j{\epsilon }_k^{\pm }(\hat{\mathbf{k}})=\mp ik{\epsilon }_i^{\pm }(\hat{\mathbf{k}}),{\epsilon }_i^{\pm }(\hat{\mathbf{k}}{)}^{\ast }={\epsilon }_i^{\pm }(-\hat{\mathbf{k}}),{\epsilon }_i^{\lambda }(\hat{\mathbf{k}}){\epsilon }_i^{{\lambda }^{\prime }}(-\hat{\mathbf{k}})={\delta }_{\lambda {\lambda }^{\prime }}.$$

$$\left({\partial }_{\tau }^2+k^2\mp \frac{k\dot{\varphi }}{\Lambda H\tau }\right)A^{\pm }(\tau ,k)=0\text{ , }$$

$$A^{+}(\tau ,k)=\frac{e^{\frac{\pi }{2}\xi }}{\sqrt{2k}}{W}_{-i\xi ,\frac{1}{2}}(2ik\tau )\simeq \frac{1}{\sqrt{2k}}{\left(\frac{k}{2|\xi |aH}\right)}^{1/4}{e}^{\pi |\xi |-2\sqrt{2|\xi |k/aH}}\text{ , }$$

$$\xi \equiv \frac{\dot{\varphi }}{2H\Lambda }\text{ . }$$

$$\frac{1}{2a^4}⟨a^2{E}^2+B^2⟩\simeq 1.4\cdot 10^{-4}\frac{H^4}{|\xi {|}^3}{e}^{2\pi |\xi |},$$

$$\frac{1}{a^3}⟨\vec{E}\cdot \vec{B}⟩\simeq -2.4\cdot 10^{-4}\frac{H^4}{|\xi {|}^4}{e}^{2\pi |\xi |}.$$

$$\ddot{\varphi }$$

$$\frac{\ddot{a}}{a}$$

$${\ddot{A}}_{\pm }(t,k)+\left(\frac{\dot{a}}{a}\right){\dot{A}}_{\pm }(t,k)+\left[\frac{k^2}{a^2}\pm \left(\frac{k\dot{\varphi }}{a\Lambda }\right)\right]A_{\pm }(t,k)=0,$$

$$A_{+}(k_i,t_i)=\frac{{\tilde{A}}_{+}(k_i,t_i)}{\sqrt{a(t_i)}}\simeq \frac{e^{-i{\omega }_i{t}_i}}{\sqrt{2a(t_i){\omega }_i}}\equiv \frac{1}{\sqrt{2k_i}}\left[\cos ({\omega }_i{t}_i)-i\sin ({\omega }_i{t}_i)\right],$$

$$k_i\equiv x_{i}a(t_i)H(t_i),x_i\gg 1,$$

Full text

1 Introduction

Compelling evidence [1] supports the inflationary paradigm as the leading explanation for the initial conditions of the early Universe. The exact scenario responsible for inflation remains however unidentified. Inflation is often assumed to be driven by a singlet scalar field, the $inflaton$, with potential and initial conditions appropriately chosen to sustain a long enough period of accelerated expansion. From a model building point of view, inflationary constructions suffer typically from being very sensitive to unknown ultraviolet (UV) physics, what can be viewed as a challenge to the robustness of the predictions. A simple way to avoid this problem relies on the introduction of some symmetry that naturally forbids certain interactions between the inflationary sector and other physics. In this respect, a promising candidate as an inflaton is a pseudo-scalar field or axion, that enjoys a shift-symmetry. If the inflationary sector is invariant under a constant shift of the field amplitude, the form of possible interactions of the inflaton with other species is severely restricted, “protecting” the inflationary dynamics from unknown UV physics.

Axions-like fields are actually ubiquitous in particle physics, as they represent a possible solution to the strong CP problem [2], potential dark matter candidates [3, 4], or interesting inflaton incarnations [5, 6]. Several implementations of axion-driven inflation scenarios have been proposed, from Natural inflation [5, 7] to $N$-flation [8, 9, 10] and stringy oriented axion monodromy models  [11, 12]. There has been in fact a significant interest in the construction of viable axion inflation models within the framework of string theory [13, 14, 15, 16, 17, 18, 19, 20]. In this paper we will not be concerned, however, with the exact implementation of axion inflation. For simplicity we will consider a quadratic potential, corresponding to the simplest type of chaotic inflation [21]. Even though these scenarios are under pressure by current CMB observations [1], the addition of a small non-minimal gravitational coupling can easily reconcile them with the data [22]. We consider a quadratic potential simply as a working ’arena’ for scenarios with a monomial shaped potential with a minimum after inflation. This choice will actually enable us to compare our results with others in the literature, where the same choice was made. Furthermore, as it will become clear later, a major aim of this paper is to introduce a robust lattice formulation of the interaction of a shift-symmetric field with a gauge field, in the context of an expanding Universe, and for this matter the choice of potential is simply irrelevant.

For any inflationary scenario, a reheating stage must follow after inflation, converting all the energy available into the particle species constituting all the matter and radiation in the Universe. Eventually, when the created particles dominate the total energy budget of the Universe and thermalize, the Universe enters into the thermal era characteristic of the “hot Big Bang” paradigm. If the inflaton has sufficiently large couplings to other species, the initial stages after inflation are typically characterized by non-perturbative particle production phenomena, typically known as preheating. Preheating mechanisms have been investigated in detail in the past: from parametric resonance [23, 24], spinodal instabilities [25, 26] and instant-like effects [27, 28] of scalar field species, to parametric excitation of fermions [29, 30, 31] or non-perturbative excitation of gauge fields [32, 33]. For reviews on preheating see [34, 35]. In this paper we will actually study the post-inflationary preheating dynamics in axion inflation. In general, once we have chosen the inflaton potential, the details of preheating depend directly on the choice of interactions between the inflaton and the fields coupled to the inflaton111Alternatively, the Universe can be reheated even in the absence of direct couplings of the inflaton to other field species, like in the context of gravitational reheating [36] and variants [37, 38] A recent criticism of standard gravitational reheating [39], reinforces however the idea that the Universe can be successfully reheated most naturally when the inflaton is coupled to other species.. In the case of axion inflation, the shift symmetry severely constrains the form of the inflaton coupling to other species.

In particular, the lowest dimensional operator compatible with a shift symmetry is the interaction term between the pseudo-scalar field and the Pontryagin or topological number density $\mathcal{K}\equiv F_{\mu\nu}\tilde{F}^{\mu\nu}$ of some gauge sector222One can also consider a derivative coupling of the axion inflaton to the axial current of the fermions [40, 41, 42, 43, 44], what leads to very interesting phenomenology, including e.g. the successful realization of baryogenesis.. We will focus on an Abelian gauge group, exploring the particle production of the gauge field due to the inflaton oscillations after inflation, when the inflaton and the gauge field interact via $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$. The phenomenology due to the presence of this coupling in axion inflation models is very rich, and it has been extensively studied in the literature. Namely, during slow-roll inflation, an exponential production of helical gauge bosons is expected [45, 46, 47, 48, 49, 50]. The exited gauge fields can then lead to interesting phenomenology such as the production of non-Gaussian density perturbations [51, 52, 20], primordial black holes [53, 20, 54, 50], magnetic fields [45, 55], (chiral) gravitational waves [6, 52, 56, 49], and a very efficient preheating mechanism [57, 58]. These effects are very sensitive to the inflaton velocity, hence becoming larger towards the end of inflation and in particular during the preheating stage after inflation, as the inflaton velocity oscillates, rapidly changing strength and sign. Let us finally note that a later decay of helical gauge fields after inflation could be used to generate the baryon asymmetry of the Universe [59, 60, 61, 62, 63].

Some works have already studied (p)reheating after axion-driven inflation scenarios, considering peturbative decays [17], linearized analyses of an Abelian gauge field coupled to an oscillating scalar field [64, 65, 66], and a full non-linear analysis [57]. Ref. [57] is actually of particular relevance to us, as it has studied in full detail the non-linear dynamics of preheating after axion inflation with a coupling $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$. The aim of our present work is actually twofold: $a)$ to introduce a refined lattice formulation of axion inflation setup’s, particularly suited for describing an interaction of the form $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$, and $b)$ to consider, as a first application of this technique, a review of the preheating stage after axion inflation studied in [57], so that we can compare numerical results.

In particular, we propose the use of the lattice formulation of an interaction term $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$ based on the formulation introduced in [67], but adapted to include self-consistent expansion of the Universe, so that both (pseudo-)scalar and gauge fields contribute to the expansion rate. Our lattice technique reproduces the continuum limit of the theory to quadratic order in the lattice spacing, and more importantly it obeys the following properties on the lattice: $i)$ the system is exactly gauge invariant, and $ii)$ shift symmetry of the axion is exact. Property $i)$ implies that physical constraints such as the Gauss law or the Bianchi identities, are exactly verified on the lattice (up to machine precision). Property $ii)$ implies that the lattice formulation has naturally embedded a construction of the topological number density that admits a total (lattice) derivative representation $\mathcal{K}\equiv F_{\mu\nu}\tilde{F}^{\mu\nu}=\Delta_{\mu}^{+}K^{\mu}$, which reproduces the continuum expression $\mathcal{K}=\partial_{\mu}K^{\mu}\propto\vec{E}\cdot\vec{B}$ up to $\mathcal{O}(dx_{\mu}^{2})$ corrections. Without this property the interaction $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$ cannot be interpreted as a derivative coupling, and hence it would not really be shift symmetric. As it is precisely shift symmetry that justifies the functional form of the interaction in first place, it seems therefore relevant to preserve exactly such symmetry at the lattice level. This leads however to an extra technical complication compared to standard formulations of Abelian gauge lattice theories: as we can no longer solve the lattice equations of motion by an explicit method, but rather by means of implicit methods that require numerically iterative procedures at each time step.

Summarizing, the lattice representation we use preserves all the relevant physical properties of the system, i.e. gauge invariance, shift symmetry and the topological nature of $F_{\mu\nu}\tilde{F}^{\mu\nu}$. Hence it represents the best approximation (to quadratic order in the lattice spacing) to the theory in the continuum. Equipped with our lattice formulation, we will study the last efolds of inflation and the preheating stage in axion-inflation, studying in detail the energy transfer from the inflaton to the Abelian gauge field as a function of the coupling strength.

This paper is organized as follows. In Sect. 2 we briefly review the dynamics during inflation, introducing the standard equations of motion and considering both analytical and numerical solutions in the continuum. In Sect. 3 we introduce a non-compact lattice formulation of an Abelian gauge theory with topological term $F_{\mu\nu}\tilde{F}^{\mu\nu}$, suitable for self-consistent expansion of the Universe. We explain in detail the virtues of this formulation and how to iterate the lattice equations of motion. In Sect. 4 we study in the lattice the last efolds of inflation and the preheating stage following afterwards. In Sect. 4.1 we discuss our strategy for initializing the fields on the lattice, determine what is the appropriate initial cutoff for the gauge field fluctuations, and study the consistency of the dynamics during the last efolds of inflation for simulations initiated at different times. In Sect. 4.2 we study preheating after inflation, investigating the transfer of energy from the axion into the Abelian gauge field. We analyse in detail the evolution of the energy densities and the power spectra of the fields. We quantify the efficiency of preheating as a function of the coupling strength. In Sect. 4.3 we demonstrate our numerical ability to preserve the relevant dynamical constraints on the lattice. In Sect. 5 we summarize our results and discuss future prospects and applications of our lattice formalism. In the appendix we present interesting aspects of our lattice implementation, and relevant tests passed by our simulations.

From now on $m_{pl}={1\over\sqrt{8\pi G}}\simeq 2.44\cdot 10^{18}$ GeV is the reduced Planck mass. The expanding Universe is described by a Friedmann-Lemaître-Robertson-Walker (FLRW) background metric $ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}d\vec{x}$, where $a(t)$ is the scale factor, and $t$ represents the cosmic time. Latin letters run over spatial dimensions, and greek letters over space-time dimensions, as usual. We assume summation over repeated indices in the continuum, but not in the lattice.

2 Axion inflation with $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$ coupling

Here we describe the continuum theory formulation of the scalar-gauge system with coupling $\phi F\tilde{F}$. We first introduce the basic field definitions and conventions in Sect. 2.1, and then discuss the field dynamics during axion inflation in Sect. 2.2, reviewing analytical results from the literature, both in the absence and presence of backreaction of the gauge field into the dynamics. In Sect. 2.3 we describe numerical solutions in the absence of backreaction, that will serve as the basis for our approach to initialize the gauge field on the lattice.

2.1 Basic framework

Let us begin by considering the action in a FLRW background for an Abelian gauge field coupled to a pseudo-scalar field $\phi$, which will be referred to as an axion from now on. Given our choice of signature $(-,+,+,+)$, we have

[TABLE]

where $\Lambda$ is a mass scale controlling the scalar-gauge coupling strength. The quadratic potential supports slow-roll inflation at super-Planckian field values, while the mass of the axion is fixed to $m=1.5\cdot 10^{13}$ GeV, in order to explain the amplitude of the CMB temperature anisotropies. The field strength of the gauge field, and its dual counterpart are defined as usual

[TABLE]

however notice that as we are in FLRW, $\epsilon_{\mu\nu\rho\sigma}$ is the completely anti-symmetric tensor in curved space-time

[TABLE]

Action (2.1) is invariant under gauge transformations $A_{\mu}\longrightarrow A_{\mu}+\partial_{\mu}\alpha(x)$, where $\alpha(x)\in\mathbb{R}$ is an arbitrary function, whilst the axion enjoys a shift symmetry [explicitly broken though by the mass term of Eq. (2.1)] as $\phi\rightarrow\phi+c$, with $c\in\mathbb{R}$ an arbitrary real constant. Defining the axion conjugate momenta, and electric and magnetic fields, as usual by

[TABLE]

the action can be re-written in a vectorial form like

[TABLE]

From here it is worth noticing that the interaction term between the axion and the gauge field is not coupled to the space-time background.

Varying the action Eq. (2.1) with respect to the axion and gauge field, leads to the respective equations of motion (EOM)

[TABLE]

where $H=\dot{a}/a$ is the Hubble parameter. The vectorial counterparts are

[TABLE]

where Eq. (2.10) is not dynamical and represents the Gauss constraint. Note that in Eqs. (2.9)-(2.10) we have kept terms that in reality are identically zero due to the Bianchi identities

[TABLE]

We keep these null terms just to remind the reader that, in the lattice, obtaining vanishing Bianchi identities is not guaranteed by default, as it depends on the lattice representation of the terms involving the gauge field in the Lagrangian. The lattice formulation of $\phi F_{\mu\nu}\tilde{F}^{\mu\nu}$ that we will describe in Sect. 3, is based precisely on the ability to obtain vanishing Bianchi identities, and a Pontryagin density of topological nature (i.e. admitting a total derivative representation).

Together with Eqs. (2.8)-(2.10), we also need to take into account Friedmann equations governing the expansion of the Universe,

[TABLE]

where the first equation describes the dynamics of the scale factor, whereas the second represents a constraint to which we will refer to as the Hubble constraint. Here $\rho$ and $p$ are the background energy and pressure densities, defined from the stress-energy tensor of a perfect fluid $T_{\mu\nu}\equiv(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}$, like

[TABLE]

In our problem at hand, the explicit form of these quantities can be obtained from the energy-momentum tensor in curved space-time

[TABLE]

where $\mathcal{L}_{m}$ is the Lagrangian in Eq. (2.1). The energy and pressure densities then read

[TABLE]

2.2 Axion inflation dynamics

The potential $V(\phi)={1\over 2}m^{2}\phi^{2}$ naturally sustains a slow-roll regime of inflation for large field amplitudes. Introducing the slow-roll parameter

[TABLE]

we expect $\epsilon_{H}\ll 1$ for most of the inflationary period, when the energy is dominated by the potential energy of the inflaton. As inflation goes on and the inflaton rolls down its potential, the slow-roll parameter grows (as the Hubble rate $H$ decreases), and eventually it reaches $\epsilon_{H}=1$, what signals the end of inflation.

Deep inside inflation, Eqs. (2.8)-(2.9) can be simplified by neglecting the gradients of the axion, $|\nabla^{2}\phi|\ll m^{2}\phi$ and $|\nabla\phi\times\vec{E}|\ll|\pi_{\phi}\vec{B}|$, as well as the backreaction of the gauge field into the axion dynamics, $|\vec{E}\cdot\vec{B}|\ll a^{3}m^{2}\Lambda\phi$. Moreover, the contribution of the gauge field to the energy and pressure densities appearing in Eqs. (2.12)-(2.13) should also be negligible, $a^{2}{\vec{E}}^{2}+{\vec{B}}^{2}\ll a^{4}m^{2}\phi^{2}$. The expansion of the Universe in such an approximation is dictated by the dynamics of an homogeneous axion,

[TABLE]

At the onset of inflation the initial value of the inflaton is super-Planckian $\phi_{i}\gg m_{pl}$ and from the slow-roll conditions we can infer that $\dot{\phi}_{i}\simeq-V_{,\phi}/3H_{i}$, $H_{i}\simeq(\phi_{i}/\sqrt{6})({m/m_{pl}})$.

The inflaton slowly rolls its potential reducing its amplitude, until it becomes eventually of order $\phi\sim m_{pl}$, just when $\epsilon_{H}\gtrsim 1$. More precisely, from the condition $\epsilon_{H}(\phi_{\rm end})=1$ we obtain $\phi_{\rm end}=1.039~{}m_{pl}$. The number of efolds before the end of inflation is simply obtained as

[TABLE]

where we use the convention that during inflation $N<0$, and in the last equality we have used properties of slow-roll, so it should only be taken as indicative333We have actually checked that the expression in Eq. (2.20) as a function of the inflaton amplitude, predicts the number efolds before the end of inflaton up the last efold, with only $\sim 3\%$ error.. In order to ensure that the inflationary period lasts $\sim 50-60$ efolds, the value of the inflaton at the moment when CMB fluctuations exited the Hubble radius must have been of the order $\phi_{\rm CMB}\sim 14-15~{}m_{pl}$.

Deep inside inflation, say $\mathcal{O}(10)$ efolds before the end, the homogeneous axion domination is very well justified. However, because of the axionic-coupling $\phi F_{\mu\nu}\tilde{F}^{\mu\nu}$, the gauge field will be excited during inflation. Thus, even if we start with a gauge field amplitude given by quantum fluctuations, as we approach the end of inflation, the gauge field will be largely amplified, and it might – depending on the scalar-gauge coupling – impact on the inflationary dynamics. Analytically, it has been studied that the backreaction of the gauge field into the axions can potentially delay the end of inflation in a significant way, see e.g. [52].

Let us briefly describe the status of the art on the analytical understanding of the gauge field dynamics and its backreaction. We can start by considering the gauge field evolution imposing homogeneity of the axion field. Switching to conformal time $\tau\equiv\int{dt^{\prime}\over a(t^{\prime})}$, Eq. (2.7) becomes

[TABLE]

In order to solve these equations, it is useful to work in Fourier space, where we decompose each mode in a basis of helicity states

[TABLE]

where the helicity vectors $\varepsilon_{i}^{\pm}(\hat{\mathbf{k}})$ are defined in such a way that they satisfy the following properties

[TABLE]

Using this mode decomposition and the properties of the helicity vectors just exposed, the differential equations for the two polarizations of the gauge field read

[TABLE]

where we have used the slow-roll condition $\tau\approx-1/aH$. Because of the different sign in Eq. (2.25), if $\dot{\phi}<0$ ($\dot{\phi}>0$) the polarization $A^{+}$ ($A^{-}$) will be exponentially amplified, while $A^{-}$ ($A^{+}$) will not, and hence can be neglected. A solution for the functional form of the excited polarization (say $A^{+}$ if we choose $\dot{\phi}<0$), normalized to match the form of quantum fluctuations in the ultraviolet limit, is given by [47]

[TABLE]

where the second expression is valid for $k|\tau|\ll 2|\xi|$, with the parameter $\xi$ is defined as

[TABLE]

The expression in Eq. (2.26) is expected to approximate well the real solution to the gauge field mode function in the limit where $\xi=\text{const}$, which is well justified deep inside inflation, but less and less accurate towards the end of it. The second expression in Eq. (2.26) exhibits clearly the fact that the gauge field is mostly amplified around the Hubble scale at each time.

Using Eq. (2.26) one can compute the mean electromagnetic energy density and the topological term, as a function of $\xi$ [47],

[TABLE]

In the Hartree approximation, the dynamics are described by the equations

[TABLE]

so by plugging Eqs. (2.28), (2.29) in Eqs. (2.30), (2.31), one can quantify the backreaction of the gauge field into both the axion and expanding Universe dynamics within such approximation. One can estimate in that way the moment during inflation (for a given coupling $\Lambda^{-1}$) when the terms $\langle\vec{E}\cdot\vec{B}\rangle$ and $\langle a^{2}E^{2}+B^{2}\rangle$ become relevant.

Before we move further into other considerations, let us discuss the allowed range of values for the scale $\Lambda$. In particular, to respect CMB limits on the amount of primordial non-Gaussianity, the function $\xi$ defined by Eq. (2.27) should respect $\xi\lesssim\xi_{\rm CMB}\simeq 2.5$ at CMB scales [68, 69], which translates (given our choice of potential ${1\over 2}m^{2}\phi^{2}$) into $1/\Lambda\lesssim 35\,m_{pl}^{-1}$. However, to ensure that primordial black holes are not over-produced by large density fluctuations, a more stringent constraint actually emerges, requiring $\xi\lesssim\xi_{CMB}\simeq 1.5-1.7$ [53, 54]. This implies that $1/\Lambda\lesssim 22-25~{}m_{pl}^{-1}$. In practice, as the primary focus of the application of our lattice formalism in this paper will be preheating, we will not need to reach so large couplings: as we will show in Sect. 4.2, the preheating efficiency does not improve further as soon as we reach $1/\Lambda\gtrsim 10~{}m_{pl}^{-1}$. Hence, a range like $1/\Lambda\leq 15~{}m_{pl}^{-1}$, which from now on we adopt as a fiducial range, suffices for our purposes.

2.3 Backreaction-less (numerical) solution

The non-linearity of the system makes actually the dynamics very complicated, and probably inaccurate, to be dealt with an analytical approach. Hence the importance of our lattice study (Sect. 4.1, 4.2), which will capture in full detail the non-linearities of the system. However, to be able to address the problem from the lattice point of view, we need to have control over the field configuration that will be introduced in the lattice as an initial condition.

As the solution Eq. (2.26) may not be so precise towards the end of inflation, we will use a different approach in order to introduce an initial condition of the fields in the lattice. In particular, back to cosmic time, but still keeping the axion homogeneous, Eq. (2.25) becomes

[TABLE]

where we have used again the helicity mode decomposition from Eq. (2.23). Once again, given our convention $\dot{\phi}<0$, only $A_{+}$ becomes tachyonic, so we will only worry about this polarization. Note that since we are working in cosmic time, no approximations have been done apart from assuming the homogeneity of the axionic field. Furthermore, since this equation is linear in $k$, we can solve it numerically in a 1-dimensional grid of discretized comoving momenta $\{k_{1},k_{2},...,k_{n}\}$, initializing each mode $k_{i}$ at a time when its wavelength is well inside the Hubble radius. Introducing a conformal transformation of the gauge field, $A_{+}\equiv a^{-1/2}{\tilde{A}}_{+}$, we can get rid off the friction term in Eq. (2.32). We then obtain that modes which are deep inside the Hubble radius $k\gg aH$ during inflation, evolve according to $\ddot{\tilde{A}}_{+}+\omega_{k}^{2}\tilde{A}_{+}\simeq 0$, with $\omega_{k}^{2}\simeq{k^{2}\over a^{2}}\gg H^{2},{\ddot{a}\over a}$. Choosing the Bunch-Davies vacuum, ${\tilde{A}}_{+}(t,k\gg aH)\simeq{1\over\sqrt{2\omega_{k}}}e^{-i\omega_{k}t}$, we can use the following initial condition,

[TABLE]

where $\omega_{i}\equiv k_{i}/a$. Furthermore, the initial time for each mode is obtained from the condition

[TABLE]

with $x_{i}$ some large positive value determining the initial “penetration” of the scale $1/k$ inside the Hubble radius. In practice, we choose $x_{i}=1000$ for all the modes, making sure in this way that we initialize the gauge field sufficiently deep inside the Hubble radius.

In Fig. 1(a) we show the numerical evolution in time of the gauge field power spectrum ${k^{3}\over m^{2}}|A_{+}(k,t)|^{2}$ for two modes, $k/m=2.58\cdot 10^{-5}$ (lower curve) and $k/m=9.50\cdot 10^{-3}$ (upper curve). We decompose each curve in real and imaginary contributions so one can see the mode oscillations deep inside the Hubble scale. As the physical wavelength approaches the Hubble radius, the field mode amplitude becomes excited, growing during few efolds. As it becomes more and more super-Hubble, its amplitude asymptotically freezes to a constant value, as expected from an interaction $\propto k\dot{\phi}$. In Fig. 1(b) we illustrate the evolution in time of gauge field power spectra for different values of the coupling. As expected, the larger the axionic coupling $\Lambda^{-1}$, the more excited the gauge field is. Let us remark that these spectra, built from the mode evolution displayed in the panel (a), are obtained by numerically evolving Eq. (2.32) with initial conditions (2.34), (2.34), considering that the inflationary background dynamics are entirely due to the homogeneous axion, according to Eqs. (2.19). In other words, numerical solutions shown in Fig. 1 ignore the possible backreaction of the gauge field itself in the inflationary dynamics. We will refer from now on to this approximation as the “backreaction-less” solution, which differs more and more from the truly physical solution the closer we get to the end of inflation, or the larger the strength of the coupling $\Lambda^{-1}$. This numerical backreaction-less solution is however very useful for two purposes: $a)$ it will serve as the initial condition on the lattice before backreaction becomes relevant, and $b)$ it shows clearly that the analytical solution fails noticeably towards the end of inflation (see Fig. 1(b)), even if the gauge field backreaction is still negligible (see dicussion in Sect. 4.1.1).

3 Lattice formulation of axion inflation with $\phi\,F_{\mu\nu}\tilde{F}^{\mu\nu}$ coupling

In this section we present a lattice formulation of an axionic-coupling $\phi F\tilde{F}$ for a $U(1)$ gauge group based on [67], but adapted to self-consistent expansion of the Universe. Armed with this formulation, we will derive the equivalent lattice equations to simulate numerically the field dynamics described by Eqs. (2.8)-(2.9), together with Eq. (2.12) describing the expansion of the Universe sourced by (pseudo-)scalar and gauge fields as in Eqs. (2.16)-(2.17), while preserving the constraint equations on the dynamics Eq. (2.10), and on the expansion rate Eq. (2.13). As the resulting lattice equations of motion are more complicated than usual and require implicit solution methods, we also present our methodology to solve them.

3.1 Lattice action in flat space-time

We will start by considering the interaction $\phi F_{\mu\nu}\tilde{F}^{\mu\nu}$ in flat space-time, so we can set up notation and introduce the proper problems arising in the lattice formulation of such term, postponing for later other problems due to the inclusion of the expansion of the Universe in the lattice. We label each lattice point as $n=(t,\vec{n})=(t,n_{1},n_{2},n_{3})$, the lattice spacing and time step as $\Delta x$ and $\Delta t$, and indicate a displacement by a lattice unit in direction $\hat{\mu}$ (spatial or temporal) as $n+\hat{\mu}$, or simply $+\hat{\mu}$. To avoid confusion, we do not consider summation over repeated indices in lattice variables. We use a non-compact formulation of lattice gauge theory, where gauge fields live on the links between two neighbouring sites, $A_{\mu}\equiv A_{\mu}(n+\frac{1}{2}\hat{\mu})$. We define the lattice versions of electric and magnetic fields as

[TABLE]

and, for later convenience, the more elaborated versions

[TABLE]

where in the right hand side we specify the natural site where each field naturally lives.

In flat space-time the action of our problem [c.f. Eq. (2.5) with $a=1$] reads

[TABLE]

Following [67], an appropriate lattice action mimicking Eq. (3.3) is

[TABLE]

where $\Delta_{\mu}^{\pm}\varphi\equiv\frac{\pm 1}{dx}\left(\varphi_{\pm\mu}-\varphi\right)$ are standard forward and backward lattice derivatives. We refer the interested reader to [67] for details on the derivation of this action. Here we simply repeat its virtues. In particular, Eq. (3.4) ensures that:

Gauge invariance is manifest with the $U(1)$ lattice transformation $A_{\mu}\rightarrow A_{\mu}+\Delta_{\mu}^{+}\alpha$, where $\alpha(n)$ an arbitrary real function. 2. 2.

All terms reproduce the continuum limit to order $\mathcal{O}(dx^{2})$. Variation of the action will therefore lead to lattice equations of motion that reproduce their continuum counterparts also to order $\mathcal{O}(dx^{2})$, see below Eqs. (3.5)-(3.7). 3. 3.

Variation with respect to the gauge field leads to $\sum_{j,k}\epsilon_{ijk}(\Delta^{+}_{j}+\Delta^{-}_{j})(E_{k}^{(4)}+E_{k,+k}^{(4)})=(\Delta_{o}^{+}+\Delta_{o}^{-})(B_{i}^{(4)}+B_{i,+i}^{(4)})$ and $\sum_{i}\Delta_{i}^{-}(B_{i}^{(4)}+B_{i,+i}^{(4)})=0$, i.e. vanishing Bianchi identities on the lattice. 4. 4.

The topological term $\mathcal{K_{L}}\equiv(F_{\mu\nu}\tilde{F}^{\mu\nu})_{L}\equiv{4}\sum_{i}{1\over 2}E_{i}^{(2)}(B_{i}^{(4)}+B_{i,+0}^{(4)})$ admits a discrete total derivative representation $\mathcal{K}_{L}=\sum_{\mu}\Delta_{\mu}^{+}K_{L}^{\mu}$, with $K_{L}^{0}=-K_{0}^{L}\equiv{2}\sum_{i}A_{i}^{(2)}B_{i}^{(4)}$, $K^{i}_{L}=K_{i}^{L}\equiv-\sum_{j,k}\epsilon_{ijk}\left(E_{j}^{(2)}A_{k,-i}^{(2)}+E_{j,-i}^{(2)}A_{k}^{(2)}\right)$.

The equations of motion describing the real time dynamics on a lattice follow from variation of action Eq. (3.4). Fixing the gauge to $A_{0}=0$, we obtain [67]

[TABLE]

where the electric field $E_{i,+\frac{\hat{0}}{2}}\equiv\Delta_{0}^{+}A_{i}$ lives at semi-integer time steps, whereas the axion conjugate field $\pi_{\phi}\equiv\Delta_{0}^{-}\phi$ lives at integer times, as the axion field itself $\phi_{+\frac{\hat{0}}{2}}$ lives at semi-integer times, due to its pseudo-scalar nature. Eq. (3.7) represents the Gauss constraint of the system. It can be shown that the set of Eqs. (3.5)-(3.7) reproduces correctly the continuum equations of motion to order $\mathcal{O}(dx^{2})$.

3.1.1 Implicit (iterative) solution for the electric fields

Standard lattice equations of a gauge system, e.g. an Abelian-Higgs model, typically admit an explicit scheme for their solution: field amplitudes at a given time can be obtained from the corresponding values at the previous time step, i.e

[TABLE]

and similarly the evolution of conjugate momenta can be obtained from previous ones as

[TABLE]

with $\mathcal{Q}_{+{\hat{0}}/{2}}$ and $\mathcal{P}_{i}$ characteristic kernels of the system of equations. This simple leap-frog scheme assumes that $\mathcal{Q}_{+{\hat{0}}/{2}}$ and $\mathcal{P}_{i}$ are functions only of the field amplitudes. The set of Eqs. (3.5)-(3.7) does not admit however such explicit scheme of evolution (unless the axion is forced to remain homogeneous). Inspecting Eq. (3.6), we observe that on the right hand side there are terms containing $E_{k,\pm j+{\hat{0}}/{2}}$. This means that to solve for the electric field $E_{i,+{\hat{0}}/{2}}$ on the left hand side, one needs to know the neighbouring electric fields $E_{k,\pm j+{\hat{0}}/{2}}$ at the same evolved time. This makes it impossible to solve explicitly for electric fields.

A potential solution to this problem is to obtain an explicit scheme by using an approximated version of the lattice equations of motion. Making the following approximation on the right hand side of Eq. (3.6),

[TABLE]

leads to an explicit, but approximated, equation of motion for the electric fields

[TABLE]

The error caused in Eq. (3.11) is of order $\mathcal{O}(\text{dt})$. Therefore, the set of coupled equations Eqs. (3.5), (3.11) will only reproduce the continuum limit of the dynamics up to order $\mathcal{O}(dx_{\mu})$, instead of $\mathcal{O}(dx_{\mu}^{2})$. Furthermore, as the left and right hand sides of Eq. (3.7) are not well balanced anymore, Eq. (3.11) leads into a violation of the Gauss law (which accumulates in time), and hence the whole scheme may exhibit inaccuracies in the numerical solution, which are out of control. Solving the dynamics with Eq. (3.11) instead of Eq. (3.6), does not represent therefore an appealing solution.

In order to preserve an order $\mathcal{O}(dx^{2})$ in the accuracy of the equations of motion, and to verify the Gauss law to machine precision, [67] proposed an implicit method for solving the electric fields iteratively. In practice the solution E_{i,+\frac{\hat{0}}{2}}{\big{|}}_{1} of Eq. (3.11) is simply viewed as an intermediate value of the updated electric fields at each lattice site. Then, by inserting this solution into the original Eq. (3.6), we obtain a new and more precise value of it, which will be called E_{i,+\frac{\hat{0}}{2}}{\big{|}}_{2}. Then, in a following iteration E_{i,+\frac{\hat{0}}{2}}{\big{|}}_{2} acts as an intermediate solution to compute E_{i,+\frac{\hat{0}}{2}}{\big{|}}_{3}, and so on. The iteration can be done $n$ times, obtaining a solution

[TABLE]

where E_{i,+\frac{\hat{0}}{2}}{\big{|}}_{0}\equiv E_{i,-\frac{\hat{0}}{2}}. In the limit $n\rightarrow\infty$ this solution approaches the exact solution $E_{i,+\frac{\hat{0}}{2}}$ to Eq. (3.6). To reproduce the continuum theory up to order $\mathcal{O}(dx^{2})$, it should be enough to iterate two times. In practice, to verify the Gauss constraint Eq. (3.7) to machine precision, we need a solution as close as possible to the exact solution to Eq. (3.6), and hence it is not clear a priori how many iterations are required for such purpose. Approximating the solution of the electric fields by E_{i,+\frac{\hat{0}}{2}}\approx E_{i,+\frac{\hat{0}}{2}}{\big{|}}_{n} will lead to a degree of violation of the Gauss constraint, the smaller the larger the number of iterations $n$. We can therefore exploit this fact and determine (by trial an error in a given simulation) the number of iterations needed to verify the Gauss law to certain degree of violation. In practice, we observe that $n=10$ iterations is typically enough to verify the Gauss law with machine precision, as presented in Sect. 4.3.

3.2 Adding expansion of the Universe in the lattice

Having laid the foundations of an interaction $\phi F_{\mu\nu}\tilde{F}^{\mu\nu}$ on the lattice, let us implement now the expansion of the Universe. We build a lattice action that reproduces Eq. (2.1) in the continuum using the same axionic coupling as in Eq. (3.3),

[TABLE]

where we choose that the scale factor $a_{+\frac{\hat{0}}{2}}$ lives at semi-integer times, so that at integeter times we write

[TABLE]

From here, we can derive the discrete equations of motion in an expanding background. Fixing the gauge to $A_{0}=0$, we obtain

[TABLE]

where Eq. (3.18) represents the Gauss law, and the conjugate momenta are

[TABLE]

By expanding each term in Eqs. (3.16)-(3.18) around their natural site, it can be shown that this set of equations reproduce to order $\mathcal{O}(dx_{\mu}^{2})$ the equations describing the continuum dynamics, Eqs. (2.8)-(2.10).

Let us consider now the expansion of the Universe. In light of the continuum equations Eqs. (2.12), (2.13), and given our choice of a scalar factor living at semi-integer times, lattice counterpart of Friedmann equations can then be written as

[TABLE]

where Eq. (3.20) is the lattice version of the Hubble constraint. From Eqs. (2.16), (2.17) we can identify the mean value (volume average) of the energy and pressure densities in the lattice as

[TABLE]

with

[TABLE]

so that $\bar{H}^{kin},\bar{H}^{grad},\bar{H}^{pot}$ are the terms associated with the axion’s kinetic, gradient and potential energy densities, and $\bar{H}^{E},\bar{H}^{B}$ are those associated with the gauge field electric and magnetic energy densities. Note that $\bar{H}^{grad}$, $\bar{H}^{pot}$ and $\bar{H}^{E}$ naturally live at semi-integer times, whereas $\bar{H}^{kin}$ and $\bar{H}^{B}$ live at integer times. Thus, we choose the semi-sums above so that $lhs$ and $rhs$ in each of the Eqs. (3.20), (3.21) are consistently living at the same time (given our choice of scale factor living at semi-integer times), and reproduce the continuum Eqs. (2.12), (2.13) to order $\mathcal{O}(dx_{\mu}^{2})$.

Let us emphasize that it is only Eq. (3.21) that controls the evolution of the scalar factor; Eq. (3.20) instead, is not a dynamical equation, but rather a constraint that should be conserved by the real time evolution of the field. Its violation indicates the degree of violation of “energy conservation” in an expanding Universe.

3.2.1 Implicit (iterative) solution for the scale factor

Eqs. (3.21), (3.22) and (3.23) show that there is no explicit scheme for evolution of the scale factor. To see this, let us label $b$ the lattice time derivative of the scale factor $a_{+0/2}$, i.e.

[TABLE]

By inspecting Eq. (3.16) we observe that to evolve the axion’s conjugate momentum

[TABLE]

we need $a_{+\hat{0}}=\frac{1}{2}(a_{+\frac{\hat{0}}{2}}+a_{+\frac{3}{2}\hat{0}})$, which contains $a_{+\frac{3}{2}\hat{0}}$, i.e. an unknown quantity which has not been computed yet at that moment. We could find $a_{+\frac{3}{2}\hat{0}}$ by using Eq. (3.24) once $b_{+\hat{0}}$ is known, but to calculate $b_{+\hat{0}}$ we need $\pi_{\phi,+\hat{0}}$. Therefore, there is no consistent scheme of evolution.

It is worth stressing that if we had chosen a scale factor living instead at integer times, an analogous problem would still appear but related with the evolution of $E_{i,+\frac{3}{2}\hat{0}}$, instead of $\pi_{\phi,+\hat{0}}$. This can be traced back to the fact that in our system, the pseudo-scalar field is forced to live at semi-integer times due to the nature of the axionic-coupling interaction, see [67] for further discussion on this.

In order to solve for the scale factor, we will then use an implicit method similar to the one used in Sect. 3.1.1 to solve for the electric field. In particular, we first make the following approximation,

[TABLE]

from where we compute a first approximated value of the evolved conjugate momenta as

[TABLE]

An approximate value of the parameter b_{+\hat{0}}{\big{|}}_{1} then follows as

[TABLE]

from where we can then also compute a first approximation of the evolved scale factor

[TABLE]

An iterative scheme can then be put forward as

[TABLE]

For a sufficiently large number of iterations $n$, we expect to reproduce the continuum evolution up to order $\mathcal{O}(dx^{2})$, ensuring a sufficiently good accuracy in the conservation of the Hubble constraint Eq. (3.20).

3.3 Lattice power spectrum of fields

Let us discuss here briefly the lattice formulation of the power spectrum of a generic field $f$, with ensemble average $\langle f^{2}\rangle$ in the continuum defined as

[TABLE]

In the lattice the ensemble average is substituted with a volume average,

[TABLE]

Using the following convention for the discrete Fourier transform,

[TABLE]

one obtains

[TABLE]

where $\langle(...)\rangle\equiv\frac{1}{4\pi|\tilde{n}|^{2}}\sum_{\tilde{n}^{\prime}\in R(\tilde{n})}(...)$ is an angular average over the spherical shell of radius $\tilde{n}^{\prime}\in[|\tilde{n}|,|\tilde{n}+\Delta\tilde{n}|]$, with $\Delta\tilde{n}$ a given radial binning, and we used the following lattice definitions: $\Delta\log k(\tilde{n})\equiv\frac{k_{IR}}{k(\tilde{n})}\,,~{}\mathbf{k}(\tilde{n})\equiv k_{IR}\tilde{n}$ and $k_{IR}\equiv\frac{2\pi}{L}$. By identification with Eq. (3.32), we conclude that

[TABLE]

The lattice spectra of the axion, gauge field amplitude and its associated electric and magnetic fields, are given by

[TABLE]

Given our choice of quadratic potential $V(\phi)={1\over 2}m^{2}\phi^{2}$, we note that we measure all observables in units of $m$, as this is a natural choice in our problem, given that $a)$ the mass of the axion is the only dimensionful scale in the problem while the gauge field backreaction remains negligible, and $b)$ it will determine the frequency of oscillations of the axion field around the minimum of its potential during the preheating stage. Hence we measure e.g. the gauge field amplitude or its energy in the corresponding powers of $m$, $\Large[\Large\langle{\vec{A}}^{\,2}\Large\rangle\Large]=m^{2}$, $\Large[\Large\langle{\vec{E}}^{2}\Large\rangle\Large]=\Large[\Large\langle{\vec{B}}^{2}\Large\rangle\Large]=m^{4}$. Correspondingly their logarithmic spectra have identical dimensions, $\Large[\mathcal{P}_{A}\Large]=m^{2}$, $\Large[\mathcal{P}_{E}\Large]=\Large[\mathcal{P}_{B}\Large]=m^{4}$, etc.

4 Dynamics on the lattice

In this section we study the evolution of the fields during the last efolds of inflation (Sect. 4.1), and during the preheating stage that follows afterwards (Sect. 4.2).

4.1 Part I. Inflation (last efolds)

Here we concentrate first in the study of the last efolds of inflation. As long as excitation of the gauge field remains “small” – we will quantify what this means in a precise manner –, we can ignore the backreaction of the gauge field in the inflationary dynamics and safely use the initial conditions described in Sect. 2.3. Whereas this should represent a safe procedure deep inside inflation, it is expected to become less and less accurate the closer we are to the end of inflation, as the excitation of the gauge field is larger and the axion may also develop some degree of inhomogeneity. However, due to the natural limitations of using a finite volume, we cannot start our simulation too many efolds before the end of inflation, say $-N\gg 1$. In fact, the initial modes captured in the lattice would be red-shifted too fast out of the box and we would miss the successive relevant ones within our limited dynamical range. It is clear that we can only aim to start our simulations some $-N\sim O(1)$ number of efolds before the end of inflation.

4.1.1 Initial conditions

To determine, for a given coupling, when the backreaction-less solution described in Sect. 2.3 becomes a bad approximation to the real dynamics, we need to determine the moment when the contribution of the gauge field to the expansion of the Universe, as well as its backreaction on the axion dynamics, become non-negligible. In Fig. 2(a) we show the evolution of the electromagnetic energy density built from the gauge field backreaction-less numerical solution, together with the energy density components of the axion field, during the last efolds of inflation. In Fig. 2(b), we plot the different terms involving the gauge field and the homogeneous axion in Eq. (2.8). For the range of couplings exhibited, $6\leq m_{pl}/\Lambda\leq 15$, we see that the axionic energy components always dominate until the end of inflation. The ratio of the electromagnetic to the axion energy density is always very suppressed, for instance, at $N=-1$ the ratio ranges between $\sim\mathcal{O}(10^{-8})-\mathcal{O}(10^{-4})$. Even at the end of inflation the suppression still persists as $\sim\mathcal{O}(10^{-4})$ for the lowest coupling and $\sim\mathcal{O}(10^{-1})$ for $1/\Lambda=15~{}m_{pl}^{-1}$. At the same time, we see that the terms $m^{2}\phi$ and $3H\dot{\phi}$ in the axion field equation of motion dominate over $\ddot{\phi}$ (slow-roll condition) all the way until $\epsilon_{H}=1$, whereas the electromagnetic backreaction term $\frac{1}{a^{3}\Lambda}\vec{E}\cdot\vec{B}$ starts to be significant for the largest couplings only at $-N\lesssim 0.5$ efolds. We see therefore that choosing as initial time for our lattice simulations any moment between $N\sim-2$ and $N\sim-1$ efolds before the end of inflation guarantees that both the energy density and backreaction of the gauge field have a negligible contribution (orders or magnitude suppressed) into the inflationary dynamics for all the relevant coupling range we considered.

Gauge field backreaction effects will be automatically accounted for, as small as they might be, once we plug in the backreaction-less solution as an initial condition in our lattice simulations. In fact, our code will solve the full fields dynamics in a $3+1$ space-time according to the discretized version of the equations of motion, naturally accounting for the gauge field contribution to the inflationary dynamics. Furthermore, let us remind that the above discussion about the weight of gauge field effects on the inflationary dynamics, based on the backreaction-less solution, also assumes that the gradient contributions from the axion field are negligible. That is a very reasonable assumption as the axion is only expected to develop inhomogeneities whenever the backreaction of the gauge field becomes noticeable. Our code will account for possible effects due to the development of axion gradients, allowing us to quantify in a self-consistent manner the validity of neglecting axion inhomogeneities in first place.

Finally, we note that we have chosen to initialize the electric field to zero, $E_{i}(t_{i})=0$. This, together with the fact that we introduce the axion initially as homogeneous, $\phi(t_{i},\vec{x})=\phi_{0}(N_{i})$, guarantees that the lattice Gauss constraint given by Eq. (3.18) is automatically satisfied by the initial conditions. Once the Gauss constraint is verified initially, it will remain conserved through all the running time of the simulations. Depending on the number of iterations used in the implicit method to solve for the electric field (recall Sect. 3.1.1), we will preserve the Gauss law with worse or better accuracy, even down to computer precision if desired. In Sect. 4.1.2 we will elaborate further on the consistency of taking initial vanishing electric field. In Sect. 4.3 we will elaborate on our numerical ability to preserve the Gauss law Eq. (3.18) and the Friedmann constraint Eq. (3.20) [c.f. Eqs. (2.10), (2.13) for the continuum counterparts.]

4.1.2 Overlapping trajectories and fluctuations’ cutoff

In the previous section we have concluded that introducing the backreaction-less solution as an initial condition $-N\sim 1-2$ efolds before the end of inflation seems consistent with the idea that all gauge field contributions to $i)$ the axion dynamics, $ii)$ to the expansion of the Universe, and $iii)$ to the induction of axion gradients, are negligible. Of course, the deeper inside inflation, the less excited the gauge field is, and hence the better such consideration becomes. From that point of view, the larger the number of efolds before the end of inflation to introduce the initial condition in the code, the better we should capture the real dynamics. However, in practice, starting at an arbitrary number of efolds before the end of inflation, say at $N<-2$, does not allow to capture well the whole dynamical range of mode excitation. Modes are red-shifted away out of our simulation box very fast during inflation, so to capture the whole dynamical range all the way from the beginning of the simulation until the end of preheating, is only possible if we start sufficiently close to the end of inflation. Actually, given that we are limited to lattice sizes of $\leq 256^{3}$ sites, starting our simulation at $-N=2$ is already barely acceptable, as we will show in a moment.

To quantify the appropriated choice of the initial number of efolds $-N$ before the end of inflation to introduce the initial condition in our code, we can do the following self-consistency exercise: for each given coupling we introduce initial conditions based on the backreaction-less solution at two different moments. We shall identify these moments $t_{2}$ and $t_{1}$, which are naturally associated to $-N_{2}$ and $-N_{1}$ efolds before the end of inflation. Let us suppose that $t_{2}<t_{1}$ (equivalently $-N_{2}>-N_{1}$). Then we evolve the lattice equations of motion for each initial condition and we label each dynamical solution for the fields as $\{\phi^{(2)}(t,\bar{x}),A_{\mu}^{(2)}(t,\bar{x})\}$ and $\{\phi^{(1)}(t,\bar{x}),A_{\mu}^{(1)}(t,\bar{x})\}$. We shall obtain each solution until the condition $\epsilon_{H}(t_{\rm end})=1$ is reached numerically in each case. If the backreaction is truly negligible at both initial times, essentially we should obtain numerically that $\{\phi^{(2)}(t,\bar{x}),A_{\mu}^{(2)}(t,\bar{x})\}=\{\phi^{(1)}(t,\bar{x}),A_{\mu}^{(1)}(t,\bar{x})\}$ for all the overlapping time range $t_{1}\leq t\leq t_{end}$, where $t_{end}$ signals the end of inflation. We will refer to this situation as correct “overlapping dynamical histories” between two initial conditions started at $t_{2}$ and $t_{1}$. If on the contrary, between $t_{2}$ and $t_{1}$ the backreaction of the gauge field becomes noticeable, the solution initialized at $t_{2}$ will not coincide later at $t_{1}$, with the one just initialized in that moment, i.e. $\{\phi^{(2)}(t_{1},\bar{x}),A_{\mu}^{(2)}(t_{1},\bar{x})\}\neq\{\phi^{(1)}(t_{1},\bar{x}),A_{\mu}^{(1)}(t_{1},\bar{x})\}$. In this case, we then have to initialize the fields at yet another earlier time $t_{3}<t_{2}$, and check again whether the dynamical field trajectories overlap or not for $t_{2}\leq t\leq t_{end}$. Nevertheless, let us remember that we cannot push $t_{3}$ too deep inside inflation because of the finite size of our lattice box. Hence, we expect the trajectories to fail overlapping also when the initialization of fields happens too many efolds before the end of inflation.

In order to perform the “overlapping trajectories” test, we need first to determine what is the typical range of excitation of the gauge field fluctuations during inflation. Based on the backreaction-less solution, we can already anticipate that within the allowed range of coupling values $m_{pl}/\Lambda\lesssim 22$, we only expect to excite the gauge field up to modes $k/m\lesssim 20$ up to the end of inflation, see Fig. 1. However, for our application to preheating, see Sect. 4.1, we will only be interested in the coupling range $6\leq m_{pl}^{-1}/\Lambda\leq 15$, and hence only modes up to $k/m\lesssim 10$ will be excited towards the end of inflation. During preheating, however, higher modes will be excited, and thus we need to guarantee that the natural UV cutoff $k_{UV}=\sqrt{3}{\pi\over dx}$ of the lattice, is large enough to encompass the physically excited modes.

Let us emphasize here that whenever we say we start a simulation at a given number $-N_{i}$ of efolds before the end of inflation, this refers to the would be efolds in a Universe dominated only by the axion. In reality, for the largest a couplings, the contribution of the gauge field to the expansion rate towards the end of inflation may become noticeable, and hence the number of efolds till we reach the condition $\epsilon_{H}=1$, will actually be larger than $-N_{i}$. Having this in mind, in Fig. 3(a) we show the evolution of the magnetic energy density power spectra for $m_{pl}/\Lambda=15$ from a lattice simulation started at $N=-2$ efolds before the end of inflation (in purple) up to the end of inflation (in red). It can clearly be seen that the spectrum shifts towards the ultraviolet and grows in amplitude, as expected. In Fig. 3(b) we illustrate the overlap of magnetic power spectra at $N=-1,0$, superimposing the evolution of simulations initialized at $N_{i}=-2.0,-1.5,-1.0,-0.5,~{}0.0$, and plotting also also spectra given by initial conditions in dashed lines. From this figure we can notice how the dynamical range varies according to the initialization moment $N_{i}$. In fact, the deeper inside inflation we initialize the fields, the more red-shifted the lattice cutoffs $k_{IR}$ and $k_{UV}$ are. This is done to ensure that we are able to capture whole range of relevant modes at the initial condition and during the simulation. It is clear that, because of the ultraviolet shift of excitations, if we were to use a fixed dynamical range for all initializations, we would not capture either the initial spectra deep in inflation, or the excited modes towards the end of it and afterwards. In light of these considerations, we decided to choose as an IR lattice cutoff adapted to the initialization moment $N_{i}$ by fixing $k_{IR}=0.5~{}a(N_{i})H(N_{i})$. This automatically leads to a lattice cutoff $k_{UV}=\sqrt{3}Nk_{IR}/2$ which, for the coupling range considered, suffices to later encompass all the physically excited UV modes during preheating.

Another important aspect concerns the modes which are in vacuum at the moment of initialization and how to deal with them. Since we are interested in the classical simulation of the system, these modes should be eliminated (set to zero) above an initial cutoff $k\geq k_{*}$. Whereas at $N=-2$ the gauge field spectrum shows fluctuations excited up to $k\lesssim 1\,m$, at the end of inflation the cutoff has grown to $k\lesssim 7-10\,m$. However, we cannot introduce as initial cutoff $k_{*}\Large|_{{}_{N_{i}}}$ the one inferred from the backreaction-less initial condition at the initial moment $N_{i}$, because there will be always harder modes $k\gtrsim k_{*}\Large|_{{}_{N_{i}}}$ in vacuum at that initial time, that will be successively excited later on as we approach the end of inflation. This is precisely the meaning of the ultraviolet displacement of the spectra observed in Fig. 3. Therefore we need to introduce a sufficiently large cutoff, otherwise modes in vacuum that are meant to be excited during the last efolds before the end of inflation will not be excited. On the other hand, the initial cutoff cannot be arbitrary large, since the energy density obtained from integrating over the logarithmic spectrum until such cutoff should be dominated by the excited part of the spectrum. The contribution to the energy density of vacuum modes should be sub-dominant because they are not physical in our classical simulation. In practice, we will choose to introduce as initial cutoff the inferred one from the spectra at the end of inflation based on the excitation range of the backreaction-less solution (for each coupling $\Lambda^{-1}$). This way we guarantee to cover the dynamical range of gauge field excitation, which is the most important criteria for the correct evolution of the system. However this choice implies that deeper inside inflation there will be a larger and larger contribution to the energy densities coming from the ultraviolet vacuum tail. This effect becomes more clear when looking at the overlap of dynamical trajectories started at successive earlier times. In Fig. 4 we plot the evolution in time of the magnetic (a) and electric (b) energy density components of the gauge field, for $1/\Lambda=9.5~{}m_{pl}^{-1}$. There, it is actually manifest the failure of the overlapping of trajectories during inflation started at different initial times if the natural UV lattice cutoff is the one considered for the initial fluctuation spectrum (i.e. if we do not introduce any initial cutoff). When instead a shorter UV cutoff is introduced according to the prescription explained before, we see that the trajectories overlap well, as they should. However, it is noticeable an unnatural decay in time from the onset of each simulation. This is a reflection of the fact that some vacuum modes are still contributing to the total energy, in particular those between the physical cutoff and our initial cutoff. This effect could be removed by means of an adaptative cutoff, that dynamically shifts towards the ultraviolet at the same speed as the real cutoff of the physical spectrum, hence capturing only the range of physically excited modes, leaving out the undesired vacuum ones. This approach is shown in dashed lines in both panels of Fig. 4 and there we see how the trajectories correctly grow during inflation and successively overlap. Nonetheless, we notice that after inflation ends, all trajectories overlap relatively well independently of whether we chose the prescribed initial cutoff or an adaptative one, as the UV vacuum tail contribution becomes irrelevant once the spectrum is dominated by the physically excited modes.

Note that in Fig. 4(b), the evolution in time of the electric energy density goes up abruptly initially, and then follows a smooth overlapping trajectory (when the right cutoff is chosen). This is because, as mentioned in Sect. 4.1.1, in order to have the Gauss law initially verified, we set the electric field to zero as initial condition, even though the amplitude of the gauge field is introduced according to the backreaction-less solution. The electric field is however rapidly generated from the equation of motion of the gauge field, as it is not possible to sustain the magnetic field given by the initial condition of the gauge field, without restoring the electric counterparts. Therefore for few time steps after the onset of the initial condition, the evolution of both the gauge and electric field amplitudes do not follow the real physical dynamics, since the electric field need first to catch up with its physical value. However, this occurs very fast and, since we expect the gauge field dynamics to be initially linear, very soon the spectrum of gauge field fluctuations starts shifting towards the UV. Thus the system will completely lose memory of this artefact after a very short period of adjustment. However, thanks to this artefact, we are able to preserve the Gauss law from the onset to the end of our simulation. Once the restoration is completed by the dynamics, the electric field energy density enters into a correct regime of overlapping dynamical trajectories initiated at different times, as clearly exhibited in Fig. 4(b). In conclusion, the initially vanishing electric field trick has no impact in the dynamics, however allowing for a simple initialization of the fields preserving the Gauss law from the beginning.

4.1.3 Choice of initialization moment

In Sect.4.1.2 we have shown that during inflation the trajectories of electromagnetic energy densities overlap consistently for initialization moments going from $N_{i}=-2.0$ to $N_{i}=0.0$ efolds. However, from Fig. 3(b) we notice a fall in power for the simulation that started at $N_{i}=-2.0$ efolds, particularly in the UV tail of the spectrum. This is a symptom of lack of dynamical range. Now we will show how well the trajectories overlap during preheating, and from there decide how much we can initialize the simulations inside inflation without losing the dynamical range of the system.

In Fig. 5(a) we can see evolution and overlap of electromagnetic energy density during the onset of preheating for the couplings $1/\Lambda=9.5~{}m_{pl}^{-1}$ and $1/\Lambda=15~{}m_{pl}^{-1}$. We notice that during the last stages of inflation the trajectories overlap very well, but eventually some discrepancies emerge, in particular for the simulations starting at $N_{i}=-2.0$ and $N_{i}=-1.5$ efolds. In Fig. 5(b) we have zoomed in the overlapping trajectories for the coupling $1/\Lambda=9.5~{}m_{pl}^{-1}$, from where we can see more in detail the discrepancy between the $N_{i}=-1.0,-0.5,~{}0.0$ efolds and the $N_{i}=-2.0,-1.5$ efolds initializations. We notice that the field dynamics starting at $N_{i}=-2.0$ exhibit a clear discrepancy after $mt\approx 1$, while the $N_{i}=-1.5$ case presents differences between $mt\approx 2$ and $\approx 5$, but eventually it catches up the other trajectories. From Figs. 5(c) and 5(d) we can inspect the magnetic power spectra at $mt=20$ for $1/\Lambda=9.5~{}m_{pl}^{-1}$ and $1/\Lambda=15~{}m_{pl}^{-1}$. As notice before, the $N_{i}=-2.0$ case manifests a decrease of power in the UV modes for both couplings. Moreover we can see the presence of extra wiggles which are not present or not very noticeable in the other spectra. From these two plots we can also convince ourselves that the $N_{i}=-1.5$ case is borderline since for IR modes we have a good overlap but the UV modes start to suffer from the same fall in power encountered for the $N_{i}=-2.0$ case.

The reason for this is very simple: when we initialize the fields as early as $N_{i}=-2.0$ efolds before the end of inflation, we need to make larger the volume of the lattice. As we actually take as our lattice IR cutoff as $k_{\rm IR}=0.5aH$, the earlier we start a simulation the more we reduce our lattice IR cutoff to capture the initial spectrum of gauge modes at that time. This reduces automatically the lattice UV cutoff, which is always given by $k_{\rm UV}=\frac{\sqrt{3}N}{2}k_{\rm IR}$. Therefore, the problem with starting our simulation at $N_{i}=-2.0$ efolds is that we barely capture well the UV range that is really needed later on for describing the gauge field excitation during preheating. As a result, the simulations started at $N_{i}=-2.0$ efolds, may eventually differ from others started later, specially in the case of relatively large couplings $m_{pl}/\Lambda\gtrsim 9$, where the enhancement of gauge field is not negligible. Whereas $N_{i}=-2.0$ as a starting time is barely acceptable, we have seen that the dynamically evolved solutions starting at earlier times, say $N_{i}<-2.0$, exhibit clear differences which are only due to the lack of a good dynamical coverage. Therefore, cannot consider them.

Based on the discussion presented so far, we could take as a typical initializing time of our simulations $N_{i}=-1.5$ or $N_{i}=-1.0$ efolds before the end of inflation. For the largest couplings simulations starting at $N_{i}=-1.5$, the results are however borderline acceptable, as they still exhibit some lack of power in UV. Because of this, we have finally decided to take as our canonical starting time as $N_{i}=-1.0$. This choice guarantees for all the couplings that the initial backreaction is negligible, based on Fig. 2(b), and at the same time allows to use a sufficiently large dynamical range in our simulations, so that all the relevant excited modes from the start to the end of the simulation, are well captured. Once again, let us notice that $N_{i}=-1.0$ refers to a Universe dominated by an axion sector only; for the largest couplings however the physical number of efolds till the end of inflation can be larger, as the gauge field backreaction in the inflationary dynamics is not negligible.

4.2 Part II. Preheating

In this section we investigate the preheating stage after inflation. In particular, we analyse the time evolution of the different energy density components and of the gauge field power spectra, and we quantify the efficiency of energy transfer from the axion to the gauge field, as a function of the coupling $1/\Lambda$. All simulations presented here start according to the prescription discussed in Sect. 4.1.1, using a grid of $N^{3}=256^{3}$ points, and initialized at our canonical choice of $N_{i}=-1.0$ efolds before end of inflation (as viewed from the axion-only dominated dynamics).

4.2.1 Energy densities and power spectra

We start by studying the evolution in time of the volume averaged energy density components. We will consider three different coupling values, representative of three distinct regimes, characterized by the fraction of energy transferred from the axion into the gauge field. In Fig. 6(a) we see the evolution of the axion and gauge field energy density components for $1/\Lambda=6~{}m_{pl}^{-1}$. This value is representative of the “sub-critical” coupling regime, for which the energy density of the gauge field remains always sub-dominant with respect to the axionic one through the whole evolution (during and after inflation). For sub-critical couplings the gradient energy density is tiny compared to the other components, so we do not display it in Fig. 6(a). We observe that the kinetic and potential energy densities of the axion behave as those of an harmonic oscillator with frequency $m$ in an expanding Universe. This is expected as in this weak coupling regime the axionic sector dominates the energy budget of the Universe (and hence the expansion rate). In other words, the backreaction of the gauge field remains mostly negligible. As expected, the oscillation-averaged kinetic and potential energy densities correctly scale as $\rho~{}\propto~{}a^{3}$ (say the envelope of the oscillation maxima), as if the Universe was matter dominated with $a\propto~{}t^{2/3}$. Concerning the gauge field, we see that the evolution of electric and magnetic energy densities are characterized by wiggles around a common trajectory. These wiggles appear after the end of inflation during the axion oscillation regime, and decrease in amplitude towards late times, as the velocity of the axion also becomes gradually smaller due to the expansion of the Universe. After few axion oscillations (approximately $\sim$ 1-2 efolds after the end of inflation) the dilution of the gauge field due to the expansion of the Universe becomes stronger than the tachyonic instability from the axionic coupling, and hence on the overall the electromagnetic energy density starts decaying in time.

The energy densities for $1/\Lambda=9.5~{}m_{pl}^{-1}$ are shown in Fig. 6(b). This case corresponds to the critical coupling value, as the gauge field is well excited and eventually, from $mt\sim 20$ ($N\sim 2$ efolds after inflation), its mean energy density becomes comparable to the axion energy density. Furthermore, we notice an enhancement of the gradient mean energy density, which reflects a larger degree of inhomogeneity of the axion field. This is due to the backreaction of the gauge field into the scalar sector, which is not completely negligible anymore. Once the total energy density becomes almost equivalently split between the axion and the gauge sector, we notice a change in the scaling of the kinetic and potential mean energy densities, as expected. In particular, the observed scaling is neither characteristic of matter- or radiation-domination, as the Universe is neither axion or gauge field dominated, but rather dominated by an even mixture of the two. We note nonetheless that the axion field mean value still oscillates around the minimum of the potential, as the axion inhomogeneities still represent only a small perturbation over the homogeneous background configuration.

Finally in Fig. 6(c), we present the energy densities for a super-critical coupling, $1/\Lambda=15~{}m_{pl}^{-1}$. In this case the gauge field mean energy densities rapidly become dominant. We see that the oscillations of the mean kinetic and potential energy densities have disappeared and that the gradient energy also becomes rapidly comparable to the others. This shows how for a super-critical coupling the gauge field backreaction, and correspondingly the inhomogeneity of the axionic field, are very noticeable, clearly impacting on the evolution of the system dynamics. We see that for this large coupling, the wiggles that appeared during the evolution of electric and magnetic energy densities for sub-critical and critical couplings, have now disappeared. This can be explained by the fact that the gauge field are not being excited periodically by the axion since the latter cannot be considered anymore an oscillating homogeneous condensate. We also note that at late times the gauge sector reaches equipartition, with both the electric and magnetic energy densities correctly scaling as $\rho~{}\propto~{}a^{-4}$, while the scale factor scales as $a~{}\propto~{}t^{1/2}$. This is expected, as it corresponds to a radiation dominated Universe with expansion rate dictated by a dominating massless gauge field.

For further clarity, we plot in Fig. 7 the different energy components of the system separately, comparing in each panel the evolution of a given component for the same three couplings $1/\Lambda$. The amount of gauge field excitation and corresponding impact in the system dynamics, comparing sub-critical, critical and super-critical couplings, is clearly appreciated in the different panels, in an component by component basis.

For a better understanding of the dynamics of preheating, we analyse now the time evolution of the power spectra of the gauge field. We will actually consider only the power spectra of the “comoving” electric and magnetic energy densities, removing the appropriate scale factor powers from the physical expressions, so that the dilution due to the expansion of the Universe is not visible, and we can focus only in the gauge field excitation due to its coupling to the axion. In Fig. 8 the evolution of the magnetic energy power spectra for the previous three representative couplings, while the same analysis for the electric power spectrum is shown in Fig. 9. We will focus our discussion about power spectra only based on the magnetic energy spectrum, as the details on the electric counterpart are rather similar.

In Figs. 8(a)-8(c) we see that the growth of the three coupling cases is, as expected, visibly distinct. This fact is highlighted by Fig. 8(d), where the line-styles represent different stages of the lattice evolution. We see that for the sub-critical case $1/\Lambda=6~{}m_{pl}^{-1}$ plot in Fig. 8(a), only infrared modes $k<k_{*}=5~{}m$ are excited, whereas higher modes are barely excited and don not play any role in the dynamics444Ultraviolet modes $k>k_{*}=5~{}m$ in the sub-critical coupling example are interpreted as remaining in vacuum and hence should not be considered in a classical simulation. Their energy density however is totally sub-dominant compared to the infrared modes, and their presence in our simulation is harmless..

The picture changes for the critical couplings $1/\Lambda=9.5m_{pl}^{-1}$ shown in Fig. 8(b), where modes above the initial cutoff are excited through the simulation, due to the non-linearity of the system. Let us recall that during inflation higher modes were successively excited as they approached the Hubble scale from initially sub-horizon scales, what corresponds to a linear regime of the system, where modes of different wavelength are not influenced by each other. However, after inflation, when the axion starts oscillating, the backreaction between the gauge field and the axion becomes more noticeable, so the dynamics become non-linear as both fields influence each other, leading to a spread of power into higher momenta. We also note that for both sub-critical and critical couplings, the power spectra grow towards a stationary configuration, where the modes start to oscillate until they relax. These oscillations are due to the change of sign and amplitude of the axion momenta $\pi_{\phi}$, which is still oscillating around the bottom of its potential.

Finally we can take a look at the super-critical case $1/\Lambda=15~{}m_{pl}^{-1}$ shown in Fig. 8(c), where the excitation of ultraviolet-modes above the initial cutoff takes place very rapidly, immediately after the end of inflation. The spread of power to higher modes occurs from the beginning of preheating (before the axion can even reach its minimum), as the growth of the mode amplitudes is actually significantly larger compared to the critical and sub-critical cases. This is precisely the spectral counterpart of what we observed in terms of the energy densities [c.f. Figs. 6(c), 7(f)], as in supercritical cases the energy is transferred so efficiently to the gauge field that from the very end of inflation the systems becomes non-linear, resulting in an immediate excitation of ultraviolet modes well above the initial cutoff. In summary, the gauge field excitation encompasses now a broader range of momenta, and it is actually significantly larger in amplitude, as it can be immediately appreciated by comparing the range of scales in both $x$ and $y$ axis of Figs. 8(a), 8(b) and 8(c). Most of the excited modes in the dynamical range of supercritical coupling, have actually amplitudes within a range of $3$ orders of magnitude. An interesting aspect appreciated in Fig. 8(c) is that the oscillations of the modes around the final stable configuration are not present, except for a small range around $k/m\approx 1$.

In Fig. 8(d) we present finally a direct comparison of the power spectra for the three chosen couplings, showing them at the beginning of the simulation $\sim 1$ efold before the end of inflation (continuous lines), at the end of inflation (dashed lines), and at the end of the simulation when the system has reached a stationary configuration (dotted lines). The difference in the amplitudes and in the range of excited modes are quite evident. However, we would like also to point out two aspects that the three cases have in common. The first one is that the peak of the power spectra is around $k/m\approx 3$, and the second is that the highest rate of excitation happens during the very first stages after the end of inflation. The transfer of energy with an axionic coupling is therefore extremely efficient, and will lead to very efficient preheating for the highest couplings.

4.2.2 Preheating efficiency

We analyze now the efficiency of preheating as a function of the coupling $1/\Lambda$. In particular we monitor the evolution in time of the ratio between the electromagnetic and total energy densities $\rho_{EM}/\rho_{tot}$. This is what is shown in Fig. 10 for $1/\Lambda$ ranging from $6~{}m_{pl}^{-1}$ to $15~{}m_{pl}^{-1}$, sampled in $\Delta\Lambda^{-1}=0.25~{}m_{pl}^{-1}$ intervals. We notice that, for small couplings $1/\Lambda<9~{}m_{pl}^{-1}$, the preheating efficiency is typically smaller than $\sim 10~{}\%$, while for large couplings couplings $1/\Lambda>10~{}m_{pl}^{-1}$, the energy stored in the electromagnetic field is above $\sim 80~{}\%$, reaching an almost perfectly reheated Universe as soon as the couplings reach $1/\Lambda\gtrsim 14~{}m_{pl}^{-1}$. The transition between inefficient ($\lesssim 50\%$) and efficient ($\gtrsim 50\%$) preheating occurs precisely around the critical value $1/\Lambda_{c}\approx 9.5~{}m_{pl}^{-1}$.

Looking in detail at Fig. 10(a) we can appreciate different patterns of the preheating efficiency. A first pattern concerns the sub-critical coupling regime $1/\Lambda\lesssim 9.5~{}m_{pl}^{-1}$, for which the gauge energy density growth is dictated in first place by the tachyonic instability, but once the scalar sector starts oscillating, the growth suddenly decreases drastically and wiggles appear. We can see that for such small couplings the growth continues at a very slow rate until $mt\approx 20-30$. From that moment onward the electromagnetic energy density starts decreasing in time, as the instability due to the axionic coupling cannot compensate the dilution due to the expansion of the Universe (which effectively behaves as matter domination after $mt\approx 20-30$).

For supercritical couplings $1/\Lambda\gtrsim 9.5~{}m_{pl}^{-1}$, the transfer of energy into the gauge field is very efficient immediately after the end of inflation, the more efficient the larger the coupling $\Lambda^{-1}$. However, for the largest couplings $1/\Lambda\gtrsim 12~{}m_{pl}^{-1}$, the transfer is so efficient that the growth of the gauge energy reaches rapidly a local maximum around $mt\sim 1-2$, and the growth then ceases. This pattern can be clearly appreciated in Fig. 10(b). The larger the coupling $\Lambda^{-1}$ the sooner the maximum in $\rho_{EM}/\rho_{\rm tot}$ is reached, and the smaller its amplitude is. The reason for this is that for larger couplings in the supercritical regime, the earlier the kinetic energy density of the axion drops below the electromagnetic energy density. When the kinetic axion energy becomes very sub-dominant (typically $\sim 10\%$ of the electromagnetic energy), the tachyonic instability looses efficiency and the excitation of gauge field stops, until the kinetic energy is eventually restored, inducing again the tachyonic instability. In order words, after reaching the local maximum, the ratio $\rho_{EM}/\rho_{tot}$ falls for short while, but begins to grow again soon afterwards, as the tachyonic instability becomes again active. This pattern can be observed for instance in Fig. 6(c), and we certify that we have always observed it for sufficiently large couplings. It is not clear to us, however, whether the gradients of the axion field are also playing any role in this pathological behaviour, as the term $\propto\vec{\nabla}\phi$ is present in the equation of motion of the gauge field Eq. (2.9).

Finally, let us note that the precise value of the critical coupling $\Lambda_{c}^{-1}\simeq 9.5m_{pl}^{-1}$ can be actually inferred from Fig. 11, where we show the maximum values that the ratio $\rho_{EM}/\rho_{tot}$ reaches for each simulation. Looking at Fig. 11 it becomes very clear why we speak in first place about sub-critical $\Lambda^{-1}\lesssim\Lambda_{c}^{-1}$ and supercritical $\Lambda^{-1}\gtrsim\Lambda_{c}^{-1}$ coupling values. For couplings below the critical value the preheating efficiency decreases actually exponentially fast. We also notice that only when the system transfers at least $\sim 80~{}\%$ of its energy into the gauge field, then it reaches eventually a stationary configuration where the energy fractions remain constant. Since in the supercritical case the Universe energy budget is dominated by the massless gauge field, the energy density of the gauge fields scales as radiation. The Universe is therefore fully reheated, entering into a radiation dominated epoch very soon after the end of inflation, for couplings $\Lambda^{-1}\gtrsim 14~{}m_{pl}^{-1}$ (see zoomed plot in Fig. 11).

4.3 Verification of dynamical constraints

In order to test the correctness of our lattice approach, we are going to study the evolution of the Gauss and Hubble constraints of the system, which are given by Eqs. (3.18), (3.20). To be more precise, we will obtain the violation of these two constraints at each time step by computing the following quantity,

[TABLE]

where $lhs$ and $rhs$ refer to the left and right hand sides of Eqs. (3.18), (3.20) respectively. Then a total violation of the constraint corresponds to $\Delta_{G,H}=1$, while the exact conservation would imply values of $\Delta_{G,H}$ close to machine precision.

We start by studying the dependence of the Hubble constraint on the number of scale factor iterations, which will be denoted by $I_{S}$, while keeping fixed the number of electric iterations, $I_{E}$, at $2$. In Fig. 12(a) we can inspect its violation for a simulation with coupling $1/\Lambda=15~{}m_{pl}^{-1}$. We notice that for $I_{S}=0$ the system reaches the complete violation of the constraint, which corresponds to the development of an instability. Indeed in this case the dynamics of the fields are completely spoiled as the values quickly diverge towards infinities. The end of the orange line, $I_{S}=0$, corresponds to the appearance of a non-numerical value in the outputs of the simulation. This implies that the explicit approximation of the scale factor, Eq. (3.26), makes the system of equations unstable. Thus the results will be totally spoiled if the iterative method presented in Sect. 3.2.1 is not used. Concerning the other three cases we notice that the conservation of the Hubble constraint has almost the same evolution. More precisely, the case $I_{S}=1$ presents a small deviation compared to $I_{S}=2,3$, which instead are indistinguishable. This behaviour repeats itself also for $I_{E}\neq 2$. From these considerations we conclude that $I_{S}=2$ is the minimal number of iterations of the scale factor in order to have the best possible conservation of the Hubble constraint.

Having fixed the number of scale factor iterations, we can study the evolution of the Gauss constraint as a function of the number of electric iterations. This dependence is shown in Fig. 12(b), for $I_{E}$ going one by one from [math] (purple) to $11$ (red) and the coupling $1/\Lambda=15~{}m_{pl}^{-1}$. We see that each extra iteration reduces the violation of the Gauss constraint by a factor $\sim 10$, until a lower limit is reached. In fact the lines given by $I_{E}=10$ and $I_{E}=11$ are superimposed and higher iterations give the same result. This demonstrates that, for a sufficiently high number of electric iterations, our lattice approach conserves the Gauss constraint up to machine precision. The question now is to know if all these iterations and the consequently increased precision in the Gauss constraint are really necessary to capture the dynamics of the problem. This is important because, for each iteration, we are adding loops of computations over the whole lattice grid, hence increasing the computational cost. In fact we have found out that each iteration increases of $\sim 7~{}\%$ the total time of simulation, as shown in Appendix A.2. The conclusion is that the dynamics are almost insensitive to the number of electric iterations, as opposed to what happens with $I_{S}$. In fact we only noticed a small departure from the physical trajectory for $I_{E}=0$, while for $I_{E}\geq 1$ substantially we have the same evolution through the whole simulation. We then decided to set both electric and scale iterations to $2$ during our simulations.

5 Discussion

In this paper we have introduced a lattice formulation of an axionic coupling ${\phi\over\Lambda}F\tilde{F}$, generalizing the original flat-space formulation from Ref. [67] to curved space with self-consistent expansion of the Universe, so that both (pseudo-)scalar and gauge field degrees of freedom contribute to the expansion rate. This formulation preserves (at the lattice level) all the properties of the continuum theory: gauge invariance, shift symmetry and the topological nature of $F_{\mu\nu}\tilde{F}^{\mu\nu}=\partial_{\mu}K^{\mu}$. We have proposed iterative methods to deal with the fact that we can no longer solve the lattice equations of motion by an explicit method. We have demonstrated that despite this technical complication, our numerical evolution is capable of preserving the Gauss constraint down to machine precision, and the Hubble expansion constraint down to the usual $\sim 10^{-4}$ accuracy in standard leap-frog evolution algorithms. Our lattice formulation represents therefore an optimal tool for studying the dynamics of axion-driven scenarios with a derivative coupling to an Abelian gauge field, as we approximate correctly the continuum theory up to quadratic order corrections in the lattice spacing, while preserving all of its relevant properties on the lattice.

Using the proposed lattice formulation we have studied the last efolds of inflation and the onset of preheating in an axion-inflation model with quadratic potential. We have characterized in detail the energy transfer from the inflaton to the Abelian gauge field as a function of the strength of the coupling $1/\Lambda$. We have first observed that the analytical solution provided in the literature (for $\xi$ constant), fails significantly towards the end of inflation. Hence, in order to set the initial condition on the lattice, it is better to solve numerically in Fourier space the evolution of the gauge field deep inside inflation, when the backreaction is still negligible. Once the lattice code is initialized, full account of the backreaction of the gauge field is automatically considered in our simulations, including the influence of the gauge field on the axion dynamics, on the expansion rate, and even on the possible development of axion inhomogeneities.

We have quantified in detail for quadratic axion-inflation the efficiency of preheating, i.e. the ability of the system to transfer a given percentage of energy from the axion into the gauge field, as a function of the strength of the coupling $1/\Lambda$. Two coupling regimes are clearly identified, sub-critical $\Lambda^{-1}\lesssim\Lambda_{c}^{-1}$ and super-critical $\Lambda^{-1}\gtrsim\Lambda_{c}^{-1}$, depending on whether the final energy fraction stored in the gauge field is below or above $\sim 50\%$. The approximated value of the critical coupling $\Lambda_{c}^{-1}\simeq 9.5~{}m_{pl}^{-1}$ can be extracted from Fig. 11, which summarizes well the preheating efficiency in this scenario. For sub-critical couplings the efficiency of preheating drops exponentially when reducing the coupling, as the tachyonic growth due to the axionic coupling can never compete with the dilution due to the expansion of the Universe. For super-critical couplings the efficiency grows very fast when increasing the value of the coupling, reheating the Universe very efficiently (with more than $\sim 95\%$ of the total energy in the gauge field) for couplings $\Lambda^{-1}\gtrsim 14\,m_{pl}^{-1}$.

Despite of all the care put in the construction of a lattice formulation that respects all the necessary symmetries of the continuum theory, our numerical results on preheating are however very similar to those presented already in Ref. [57], where lattice simulations of the same system were also explored with a less refined lattice formulation. Actually, a flat-space version of our lattice formulation has been used before in the context of the anomalous non-conservation of fermion number in Abelian gauge theories at finite temperature [70], where the statistical study of the Chern-Simons number $Q(t)\propto\int_{0}^{t}dt^{\prime}\int d^{3}xF\tilde{F}=N_{\rm CS}(t)-N_{\rm CS}(0)$, with $N_{\rm CS}(t)\propto\int d^{3}x\vec{A}\vec{B}$, proved to be very sensitivity to the lattice formulation used. It seems therefore that depending on the observable, and perhaps on the context, the introduction of a lattice representation of $F\tilde{F}$ that can be expressed as a total derivative (hence preserving exactly the shift symmetry on the lattice), matters or not. Clearly, in the context where we have applied our formalism here – preheating after axion inflation –, we are forced to conclude that the precise conservation on the lattice of gauge invariance and shift symmetry, seems not to be particularly relevant. We consider our analysis however, the demonstration of this fact; something that, admittedly, we were not expecting.

Our lattice formulation represents in any case a powerful tool for studying the non-linear dynamics of any system composed by an axion-like field (with arbitrary potential) and an Abelian gauge field coupled to the latter via a Chern-Simons interaction, considering self-consistent expansion of the Universe, and preserving exactly gauge and shift symmetries on the lattice. We expect therefore that our present paper will be the first one of a series, where we plan to apply our lattice formulation to review and investigate further phenomenological aspects of axion-driven inflation scenarios with a $\phi F\tilde{F}$ coupling, particularly in the non-linear regime where analytical studies cannot be precise.

Acknowledgments

We are very grateful to M. Shaposhnikov and V. Domcke for discussions on the project. The work of DGF was supported partially by the ERC-AdG-2015 grant 694896 and partially by the Swiss National Science Foundation (SNSF).

Appendix A Appendix

A.1 Dynamical range study

An important aspect which can be analysed from the power spectra of the fields is how well we are covering the dynamical range of momentum of our system. In order to do this we run simulations with slightly different infrared lattice cutoff which we define by $k_{IR}=x_{k}aH$, where $x_{k}<1$ is the arbitrary parameter and $aH=a(t_{i})H(t_{i})$ is defined by the initial conditions. We know that once $k_{IR}$ is defined, also the UV cutoff $k_{UV}$ is fixed by the relation $k_{UV}=\frac{\sqrt{3}N}{2}k_{IR}$. Then, by comparing the resulting power spectra we will be able to know how different values of $x_{k}$ affect the dynamics. If no evident distortion of the spectra appear, it means that our momentum coverage is sufficient to capture the dynamics of the system. The results of this study are presented in Fig. 13, where we show the magnetic energy density power spectra for three representative couplings at $mt=50$. We see that the range of momentum is shifted for each different $x_{k}$, accordingly to our definition of $k_{IR}$. Nevertheless, the picture is broadly the same and no particular difference can be spotted, even though $x_{k}=0.3$ is borderline in the UV. We then conclude that with our choice of $x_{k}=0.5$ we are well capturing the dynamical range of the system.

A.2 Time-demand of iterative methods

In this appendix we present a study on the time-demand of iterative methods. More precisely we would like to quantify how much costs, from the computational point of view, each extra loop imposed by these methods. In order to do this we will compute, for each combination of iterations, the ratio $\frac{t_{i,j}}{t_{0,0}}$, where $t_{i,j}$ is the time needed for a code implementing $I_{E}=i$, $I_{S}=j$ iterations to run the simulation.

In figure 14(a) we show the results for $I_{E}=0,1,\dots,15$ and $I_{S}=0,1,2,3,4$. From it we see that the increase on time is proportional to both number of electric and scale factor iterations. This is not surprising since each iteration adds always the same amount of computations to the simulation. From fitting the data, we find that the mean incremental percentage for each extra electric iteration is

[TABLE]

Moreover, we figured out that for each extra scale factor iteration the computational time increases by $\sim 8\%$.

A.3 Statistical fluctuations

Here we present a statistical study to understand how much the fluctuations due to randomness in the initial conditions affect the dynamics of the system. In fact, each initialization of gauge field follows a random realization of a Rayleigh distribution, which rms value coincides with the input value of the initial spectra computed following the prescription of Sect. 2.3. In Fig. 14(b) we see that basically the dynamical evolution in insensitive to the statistical realization, as the standard deviation at each time step is essentially negligible.