Gauge-field production during axion inflation in the gradient expansion formalism

We study the explosive production of gauge fields during axion inflation in a novel gradient expansion formalism that describes the time evolution of a set of bilinear electromagnetic functions in position space. Based on this formalism, we are able to simultaneously account for two important effects that have thus far been mostly treated in isolation: (i) the backreaction of the produced gauge fields on the evolution of the inflaton field and (ii) the Schwinger pair production of charged particles in the strong gauge-field background. This allows us to show that the suppression of the gauge-field production due to the Schwinger effect can prevent the backreaction in scenarios in which it would otherwise be relevant. Moreover, we point out that the induced current, $\boldsymbol{J} = \sigma \boldsymbol{E}$, also dampens the Bunch-Davies vacuum fluctuations deep inside the Hubble horizon. We describe this suppression by a new parameter $\Delta$ that is related to the time integral over the conductivity $\sigma$ which hence renders the description of the entire system inherently nonlocal in time. Finally, we demonstrate how our formalism can be used to construct highly accurate solutions for the mode functions of the gauge field in Fourier space, which serves as a starting point for a wealth of further phenomenological applications, including the phenomenology of primordial perturbations and baryogenesis.


I. INTRODUCTION
Gamma-ray observations from distant blazars provide indirect evidence for the presence of magnetic fields in voids of our Universe [1][2][3][4][5][6][7][8][9][10] (for a recent review, see Ref. [11]). Together with observations of the cosmic microwave background (CMB) [12][13][14][15], this constrains the strength of the fields in the range 10 −17 G B 0 10 −10 G. The very large coherence length of these fields [16], measured in Megaparsecs, implies that they are of cosmological nature, because astrophysical mechanisms are inefficient in voids due to their small matter content. In principle, it is possible to inject magnetic fields into voids through outflows of matter from galaxies and active galactic nuclei [17][18][19][20]; however, it is quite problematic to attain an extremely large correlation length of the observed fields. The possible cosmological origin of magnetic fields in voids has another important advantage because such primordial fields could serve as seed magnetic fields for galaxies and galaxy clusters [21]. These seed magnetic fields could be later enhanced astrophysically through different types of dynamo [22][23][24][25][26] and adiabatic compression [27] to the strength of magnetic fields observed in galaxies and galaxy clusters [27][28][29][30][31][32][33][34][35].
The mechanism of inflationary magnetogenesis proposed in Refs. [36,37] naturally explains the large correlation length of magnetic fields in voids. This means that, in addition to the abundance of chemical elements, the CMB radiation, and the large-scale structure of galaxies and their clusters, the observed magnetic fields in voids could provide a unique source of information about physical processes in the early Universe at very high energies that are unattainable in the laboratory. Since Maxwell's action is conformally invariant, fluctuations of the massless electromagnetic field are not enhanced in a conformally flat inflationary background [38]. Consequently, the conformal invariance should be necessarily broken in order to generate gauge fields during inflation. Such a breaking can be done in may ways, e.g., by introducing an interaction with a scalar or pseudoscalar inflaton field or with spacetime curvature [36,37,39,40]. In our study, we will consider the axial coupling of an Abelian gauge field to a pseudoscalar "axion" inflaton field, which has attracted a lot of attention in the literature , mainly because it produces magnetic fields with nonzero helicity, which facilitates their survival in the primordial plasma after inflation.
In addition to magnetic fields, strong electric fields are always generated in inflationary models. Due to the Schwinger effect [63][64][65], these strong electric fields produce pairs of charged particles and antiparticles and quickly form an ultrarelativistic plasma. Such a plasma efficiently screens electric fields and, therefore, strongly influences the generation and evolution of the gauge fields, especially near the end of inflation and during reheating. The Schwinger pair production during inflation and its impact on magnetogenesis was investigated in Refs. [53,57,59,. For a recent application in the context of cosmological relaxation as a solution to the hierarchy problem, see also Ref. [89].
To handle the complicated dynamics of the inflaton and gauge fields, taking into account the Schwinger effect and the backreaction of the generated gauge fields and primordial plasma on the cosmological evolution, a novel gradient expansion formalism was developed by three of us in Ref. [57]. The formalism was later extended also to the case of the kinetic-coupling model [90]. In the standard approach to magnetogenesis, one works with separate Fourier modes of the gauge field which evolve in a given inflationary background. In contrast, in the gradient expansion formalism, one considers vacuum expectation values of a truncated set of bilinear functions of the electric and magnetic fields in coordinate space that include all physically relevant modes at once. It was shown that, even taking into account a relatively small number of bilinear functions, it is possible to describe the electric and magnetic energy densities with an accuracy of a few percent during the whole stage of inflation. The formalism also takes into account the fact that the number of relevant modes constantly grows during inflation as new modes cross the horizon and undergo the quantum to classical transition. 1 To account for the growth of the number of relevant modes outside the horizon, boundary terms are added to the equations of motion for the electromagnetic bilinear functions.
In the present paper, we will study the problem of generating helical magnetic fields in a model with an axial coupling of the gauge field to the inflaton. Using the gradient expansion formalism of Ref. [57], which takes into account the backreaction of the generated fields and the Schwinger effect, we will develop an approximation scheme that allows us to reach an unprecedented accuracy in the numerical results -the error compared to the corresponding mode by mode solution is always less than one to two percent during inflation -and compare our results to those in the literature. Here, a characteristic feature and advantage of our formalism consists in the fact that we do not rely on an iterative procedure that needs to be repeated over and over again before it converges to a self-consistent result. Instead, our evolution equations only need to be evolved forward in time once in order to generate the desired output.
We will be particularly interested in the case in which the Abelian gauge field coupling to the axion inflaton field is identified with the hypercharge gauge field in the Standard Model. Consequently, we will evaluate the electric conductivity σ of the charged plasma generated during inflation based on the Standard Model particle content, assuming that all fermions remain massless during inflation. Specifically, we will use the expressions for the induced current and electric conductivity in Ref. [59], which account for the presence of both electric and magnetic background fields. On the technical side, we will work out an improved description of the vacuum fluctuations deep inside the horizon that reflects the damping of the gauge-field amplitude by the charged plasma. This will lead us on the one hand to improved boundary terms in our evolution equations and on the other hand to the notion of a damped Bunch-Davies vacuum, which we describe in terms of a new parameter ∆. As we will see, ∆ is related to the time integral over the conductivity σ and hence is nonlocal in time. An important outcome of our analysis therefore is the fact that the state of the system at any given moment in time t does not only depend on the choice of parameter values in the Lagrangian and quantities such as the instantaneous velocity of the inflaton field; it also depends on the entire prehistory leading up to this moment in time, which controls to what extent the vacuum fluctuations of modes inside the horizon have already been damped in the course of inflation at times t ′ ≤ t.
The paper is organized as follows. The axial-coupling model for generating magnetic fields during inflation is described in Sec. II. In Sec. III, the system of equations for the electromagnetic bilinear quantities is derived using the gradient expansion formalism. Numerical solutions of this system of equations are given and discussed in Sec. IV. In this section, we will compare the outcome of our analysis to other methods and results available in the literature whenever possible. This will allow us to validate our formalism; first in the case of no backreaction and no Schwinger effect and then subsequently in the case of relevant backreaction and no Schwinger effect. In the third step, we will then extend our analysis and include the Schwinger effect, which goes beyond existing results in the literature. An important result of this part of our analysis will be that Schwinger pair production can suppress the energy density of the electromagnetic field to such an extent that no backreaction effects occur after all, even though backreaction effects would be absolutely essential in the absence of the Schwinger pair production. In addition, we will present the Fourier spectra of the electric and magnetic fields in this section as well as their relation to the spectra in the (damped) Bunch-Davies vacuum. Section V finally summarizes the findings obtained in this paper and contains an outlook on possible future steps. In this section, we will also comment on possible phenomenological applications of our results, such as the generation of primordial scalar and tensor perturbations and baryogenesis from hypermagnetic fields. Throughout the paper we use natural units and set = c = 1 with the reduced Planck mass equal to M P = (8πG) −1/2 = 2.435 × 10 18 GeV. We assume that the Universe is described on cosmological scales by a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric in terms of the cosmic time, g µν = diag {1, −a 2 (t), −a 2 (t), −a 2 (t)}.

II. AXIAL-COUPLING MODEL
Let us describe the axial-coupling model that we use for the analysis of the generation of primordial magnetic fields. In this model, the conformal invariance of the Maxwell action is broken by means of the axial coupling of the electromagnetic field A µ with the pseudoscalar inflaton field φ. The corresponding action reads as 2 where g = det g µν is the determinant of the spacetime metric, V (φ) is the inflaton potential, I(φ) is the axial-coupling function, F µν = ∂ µ A ν − ∂ ν A µ is the gauge-field strength tensor, and is the corresponding dual tensor; ε µνλρ is the absolutely antisymmetric Levi-Civita symbol with ε 0123 = +1. The last term in Eq. (1) is the gauge-invariant Lagrangian of a generic matter field χ charged under the U (1) gauge group and, therefore, coupled to the electromagnetic four-potential A µ . For the sake of generality, we will not specify this term and assume that it describes all charged fields in the model. Action (1) implies the following Euler-Lagrange equations for the inflaton and electromagnetic fields, where is the electric four-current. Equation (4) should be supplemented by the Bianchi identity for the dual gauge field strength tensor The energy-momentum tensor equals where the last term describes the contribution of the charged matter fields. Now we use the FLRW metric and choose the Coulomb gauge for the vector potential A µ where A µ = (0, −A). Then, the three-vectors of electric E = (E 1 , E 2 , E 3 ) and magnetic B = (B 1 , B 2 , B 3 ) fields can be defined as Note that these are the physically measured fields by a comoving observer; therefore, we included the scale factor in their definition. The gauge field stress tensor and its dual tensor are expressed in terms of the electric and magnetic fields as follows: The cosmic expansion rate (the Hubble parameter H =ȧ/a) is determined by the Friedmann equation 2 In fact, Aµ represents a generic Abelian gauge field axially coupled to the inflaton. However, for simplicity, we will sometimes refer to it as the electromagnetic field. Our numerical analysis will be performed for the specific case of the Standard Model hypercharge field.
where the total energy density ρ is given by the zero-zero component of the energy-momentum tensor (7), Here the two terms in square brackets correspond to the spatially homogeneous inflaton field, the next term gives the gauge-field contribution (angular brackets denote the vacuum expectation value), while the last one is the counterpart for the charged matter fields. The electric current four-vector can be decomposed as We assume that charged particles were absent in the Universe initially and were produced only later in particleantiparticle pairs. Therefore, we set the charge density to zero, ρ c = 0. On the other hand, the current density three-vector J may be nonzero in the presence of the electromagnetic field. Then the equations of motion (3)-(4) and (6) take the following form in the three-vector notation, Finally, in order to close the system of Maxwell's equations, we need to specify the electric current. It was shown in the literature (see, e.g., Refs. [59,66,71]) that the induced current of charged particles produced by the Schwinger effect in a constant electromagnetic field in de Sitter spacetime satisfies Ohm's law, where σ is the generalized conductivity, which depends only on the absolute values of electric and magnetic fields. We will utilize this approximation in our analysis too, even though we will consider time-dependent electric and magnetic fields. The applicability of this approximation will be discussed in detail in Sec. III B.

III. GRADIENT EXPANSION FORMALISM
In this section, we derive the system of equations describing the self-consistent evolution of the inflaton, gauge fields, and charged particles produced by the Schwinger effect during axion inflation.

A. Equations for bilinear electromagnetic functions
Solving the system of Maxwell's equations requires determining the dependence of the gauge field on the spatial coordinates, which makes numerical calculations very demanding. Therefore, as mentioned in the Introduction, it is advantageous to apply the gradient expansion formalism, where one considers vacuum expectation values of bilinear electromagnetic functions in coordinate space that include all physically relevant modes at once. These functions are and contain spatial derivatives only as powers of the curl. This allows us to get a closed system of equations for these quantities, because Maxwell's equations (13) and (14) contain spatial derivatives only in such a form. Using Eqs.(13)-(15), we get the following equations of motion for the bilinear electromagnetic functions, Note that the equation of motion for the nth order function always contains at least one function with the (n + 1)th power of the curl. As a result, all equations are coupled into an infinite chain that needs to be truncated in practice. The right-hand sides of Eqs.(20)-(22) contain the contributions due to boundary terms. Their necessity is dictated by the following consideration. Modes inside the horizon correspond to vacuum fluctuations of the quantum gauge field. They oscillate in time without significant change of their amplitude. Therefore, their contribution to the electromagnetic energy density should be excluded. However, modes outside the horizon do not oscillate and can be treated as the Fourier modes of a classical electromagnetic field [91]. During inflation, modes that were initially inside the horizon cross the horizon and begin to behave classically. Therefore, the number of physically relevant modes (outside the horizon) constantly grows during inflation. To take this fact into account, one needs to introduce the corresponding boundary terms for the bilinear electromagnetic functions. The resulting system of equations forms a closed set of ordinary differential equations that describe the self-consistent evolution of classical observables in the form of quadratic functions of the electric and magnetic fields with an arbitrary power of the curl.
We would like to mention that the boundary terms for the system of Eqs. (20)- (22) were derived for the first time in Ref. [57]. In this study, we calculate the boundary terms more accurately, including the impact of the Schwinger effect. In addition, we propose a new way of truncating the infinite chain of equations of the system that allows us to reach a high accuracy in the numerical results.

B. Schwinger conductivity
Before calculating the boundary terms, let us discuss the Schwinger pair production. In constant and uniform collinear electric and magnetic fields in an expanding Universe, this effect was studied in Ref. [59]. In the case of one Dirac fermion species with mass m and charge Q, the conductivity induced by the Schwinger effect has the form where |B| and |E| are the absolute values of the collinear magnetic and electric fields and e is the U (1) gauge coupling constant. The contribution of a single Weyl fermion is twice smaller. It is important to remember that this expression was derived in the strong-field regime, |eE| ≫ H 2 , which is the most important one for physical applications. This expression is in good agreement with the result for flat Minkowski spacetime (see, e.g., Refs. [92][93][94][95]) if we replace 1/(3H) with the time ∆t during which the external field is switched on. Moreover, for small magnetic field |B| ≪ |E|, it reduces to the expression for the fermionic Schwinger conductivity in a purely electric field in de Sitter spacetime, σ f,0 = (e|Q|) 3 |E|/(6π 3 H) exp(−πm 2 /|eQE|) obtained in Ref. [71]. Using the same logic, it is easy to construct the Schwinger conductivity for scalar charge carriers in the strongfield regime. In flat Minkowski spacetime, the pair production rate for scalars differs from that of fermions with the replacement of coth(π|B|/|E|) by (1/2)cosech(π|B|/|E|), where cosech x = 1/sinh x [93,94]. Then, in de Sitter spacetime, we expect that A consistency of this result can be checked by taking the limit of a vanishingly small magnetic field, |B| ≪ |E|. In this case, we obtain σ s,0 = (e|Q|) 3 |E|/(12π 3 H) exp(−πm 2 /|eQE|), which is the well-known expression for the Schwinger conductivity in a constant electric field in de Sitter spacetime computed in Ref. [66]. We would like to emphasize that expressions (23) and (24) for the Schwinger conductivity have been derived in the approximation of constant collinear electric and magnetic fields in the de Sitter background. Therefore, we assume that the electric and magnetic fields as well as the Hubble parameter evolve sufficiently slowly during inflation so that at any given moment of time the conductivity is determined by the corresponding values of these quantities.
Another our assumption is the collinearity of the electric and magnetic fields. In the axial-coupling model, the generated fields are almost maximally helical (for sufficiently large axial coupling). In such a case, the electric and magnetic fields are indeed nearly collinear. For instance, in the absence of the Schwinger effect and for |ξ| ≡ |(dI/dφ)(dφ/dt)|/(2H) 1, there exist quite simple analytical expressions for the generated fields [42,57], Then the angle between the electric and magnetic fields can be estimated as which is indeed a relatively small angle. This is the minimal possible value of θ which can be reached in the free case without backreaction and the Schwinger effect in the limit |ξ| ≫ 1. For smaller |ξ|, the angle can be greater, while in the presence of the Schwinger effect it can be even smaller than the value in Eq. (26), see Figs. 6(b) and 7(b) in Sec. IV.
In the subsequent analysis, we will assume that the U (1) gauge group is the Standard Model hypercharge group U (1) Y and e = g ′ , Q = Y . We consider the electroweak-symmetric phase where all SM fermions are massless. This requires the SM Higgs field to be stabilized at the origin in field space during inflation, since any nonzero Higgs field value would otherwise break the electroweak symmetry and render the SM fermions massive (except for neutrinos, possibly). One possibility to stabilize the Higgs field consists, e.g., in a nonminimal coupling to the Ricci curvature scalar R, which endows the Higgs field with an effective mass of the order of |R| ≈ 12H 2 times a positive coupling constant during inflation. In the following, we will therefore focus on scenarios in which the Higgs field is very heavy during inflation, m 2 H ≫ g ′ |E|, such that its contribution to the Schwinger current can be neglected. We express the values of the electric and magnetic fields in the conductivity through the bilinear functions of the gauge field, i.e., we use |B| = √ B (0) and |E| = √ E (0) , and finally get the following result (see Appendix A for more details), We take into account the running of the coupling constant g ′ due to all SM particles [96], At the energy scale of the Z-boson mass, m Z ≈ 91.2 GeV, the gauge coupling equals g ′ (m Z ) ≈ 0.35. As a characteristic energy scale µ relevant for the Schwinger pair production, we use The coefficient in front of the logarithmic term in Eq. (28) may be affected by mass threshold effects related to the mass of the Higgs field during inflation. However, the Higgs contribution is 40 times smaller than that of all SM fermions, so these effects can only lead to very small changes, therefore, we will neglect them in the following. Finally, we would like to discuss the evolution of the energy density of the produced charged particles, ρ χ , which affects the cosmic expansion rate. The equation of motion for this quantity can be derived from the energy conservation law. Indeed, multiplying Eq. (12) for the inflaton field byφ, we rewrite it in the forṁ where ρ inf =φ 2 /2 + V (φ) and p inf =φ 2 /2 − V (φ) are the energy density and the pressure of the inflaton field, respectively. Further, the equation describing the evolution of the electromagnetic energy density ρ em = E 2 +B 2 /2 = (E (0) + B (0) )/2 can be obtained from Eqs. (20) and (22) for n = 0, The first term on the right-hand side describes the increase of the electromagnetic energy density due to new modes crossing the horizon during inflation; therefore, it can be thought of as a vacuum source term. Comparing Eqs. (30) and (31), we conclude that the second term on the right-hand side of Eq. (31) describes the energy transfer from the inflaton to the gauge field due to the axial coupling. Finally, the last term in this equation determines the energy loss due to the Schwinger effect. Then we have the following equation for the energy density of produced particles, where we assumed that m 2 ≪ g ′ |E| so that the produced particles are ultrarelativistic; in fact, we assume them to be massless. In the opposite case, the factor 4H should be replaced by 3H as for nonrelativistic particles.

C. Boundary terms
In order to derive the explicit form of the boundary terms, let us consider the quantized gauge field where A λ (t, k) is the mode function, ǫ λ (k) is the polarization three-vector,â k,λ (â † k,λ ) is the annihilation (creation) operator of the electromagnetic mode with momentum k and circular polarization λ = ±, and k = |k|. The polarization vectors satisfy the following properties The creation and annihilation operators have the canonical commutation relations Substituting decomposition (33) into Eqs. (17)- (19) and taking the vacuum expectation values, we obtain where Note that the integration in the above expressions proceeds over wave numbers less than k h (t), which is the wave number of the mode that crosses the horizon at a given moment of time t and whose explicit form will be given below. Clearly, the dependence of k h on time is the underlying reason for the appearance of the boundary terms. Let us forget for a while that the conductivity depends on the electric and magnetic fields and consider σ(t) simply as a function of time. Then, Eqs. (4) and (33) tell us that the equation of motion for A λ (t, k) is given bÿ where σ(t) appears in the frictionlike term σȦ λ (t, k) on the left-hand side in addition to the usual Hubble friction term HȦ λ (t, k). As we will see shortly, it is this additional friction term induced by the Schwinger conductivity that will lead us to the notion of a damped Bunch-Davies vacuum deep inside the horizon. The longer the friction term is active, the stronger will be the damping of the electromagnetic field. In order to find an approximate solution for the mode equation (40), it is convenient to do the following changes of function and variable where η(t) = t dt ′ /a(t ′ ) is the conformal time. If we assume that during inflation the cosmic expansion can be approximated by the de Sitter solution η ≃ −1/(aH), the function f λ satisfies the following equation For further convenience, we introduce the notations Although the Schwinger conductivity depends on electric and magnetic fields, we assume that they vary slowly during inflation. Therefore, the term withσ in square brackets of Eq. (42) can be neglected. Then we get the following simplified equation Deep inside the horizon, when the first term in square brackets dominates, the solution must satisfy the Bunch-Davies boundary condition [97] f We emphasize that Eq. (45) does not fully describe the gauge-field mode function inside the horizon in the presence of finite conductivity. Indeed, the full mode function is given by Eq. (41) and includes an exponential damping factor, where we introduced the new parameter ∆, which is defined in terms of the exponential of the integrated conductivity, This parameter suppresses the gauge-field amplitude on small scales. It is, moreover, nonlocal in time, since it depends on the conductivity at all times t ′ ≤ t. A mode crosses the horizon when the expression in the square brackets in Eq. (44) vanishes for the first time at least for one polarization, which is the case when In the vicinity of the horizon crossing when t ≈ t h , ξ(t) and s(t) in Eq. (44) can be replaced by ξ(t h ) = constant and s(t h ) = constant, respectively. Then the solution of Eq. (44) satisfying the Bunch-Davies boundary condition (45) can be expressed in terms of the Whittaker function (see Appendix B for more details), The derivative of this expression with respect to z can be computed by using Eq. (B4) in Appendix B. Now we are ready to calculate the boundary terms in Eqs. (20)- (22). Let X represent an electromagnetic bilinear function (any of the functions E (n) , G (n) , B (n) ), whose spectral decomposition reads Here, k h (t) is the momentum of the mode that crosses the horizon at time t. It can be found by inverting the function The function X varies in time because (i) the integrand evolves in time and (ii) the upper integration limit is time-dependent. The time evolution of the integrand is described by the left-hand side of Eqs. (20)- (22). The time dependence of the upper integration limit leads to the boundary terms on the right-hand side of these equations. Thus, the boundary term for X equals its time derivative coming from the variable upper integration limit, Using the spectral decompositions (36)- (38) and the approximate solution (49) valid in the vicinity of the horizon crossing, we obtain the following boundary terms where k h (t) is given by Eq. (51) and with r(ξ, s) = |ξ| + ξ 2 + s + s 2 . For large ξ, these expressions can be expanded into series in inverse powers of r(ξ, s), which is convenient for numerical computations. They are listed in Appendix B. We point out again that, in the presence of conductivity (e.g., caused by the Schwinger effect), the evolution of the gauge field becomes nonlocal in time. This nonlocality originates from the exponential factor ∆(t) in the boundary terms, which contains information about the conductivity in all preceding moments of time, see Eq. (47). This implies that the generation of gauge fields at a given moment of time depends not only on the parameters H, ξ, and s at the same moment of time, but also on the entire prehistory of the system. We also emphasize once more that the presence of boundary terms is caused by the cutoff in the spectrum of physically relevant gauge-field modes, which changes in time. Moreover, the explicit form of these terms depends on the choice of the cutoff. However, expression (51) for the cutoff is quite natural, because for any given k > k h (t), the effective frequency of the mode is always positive for all times t ′ ≤ t. Consequently, this mode has not yet experienced any tachyonic instability. This can be clearly seen from the spectra of generated fields. Figure 1 shows the spectral densities of the electric (blue solid lines) and magnetic (red dashed lines) energy densities as well as of the scalar product Finally, we would like to comment on the behavior of the spectral density of E · B , which is shown by the green dashed-dotted lines in Fig. 1. In the case of a free gauge field, whose mode function is given by the Bunch-Davies solution (45), such a spectral density identically vanishes (while the spectral densities of ρ E and ρ B are equal and scale as ∝ k 4 , see the purple dotted lines in Fig. 1). In the axial-coupling model, we obtain a nonvanishing result for the spectrum of E · B for subhorizon modes; however, it is suppressed compared to the spectra of ρ E and ρ B .

D. Truncation of the chain of equations
As we already mentioned above, the system of equations (20)- (22) in general consists of infinitely many equations which form a chain of coupled equations because the equation for the quantity of nth order contains quantities of order (n + 1). However, there is a physical argument which allows us to truncate this chain at some finite order. It can be deduced by considering the spectral decomposition (36)- (38). Any quantity X (n) (X ∈ {E, G, B}) can be represented in the following form where spectral density function X λ (t, k) does not depend on n. For a sufficiently large n, the ultraviolet part of the spectrum in vicinity of the horizon scale k h gives the dominant contribution to the integral. In this region, we can approximate where the factor 1/a was added for convenience and p is a spectral index typically of order unity (we take it to be the same for both polarizations because, for the modes which have just crossed the horizon, mode functions behave similarly to the corresponding vacuum solutions which are the same for both polarizations; however, this is not essential for the following discussion). Then, we can approximately write This expression implies that increasing n by two does not change the second multiplier. Therefore, we have where we assumed n ≫ p = O(1). This gives us an opportunity to express higher-order quantities in terms of the lower-order ones and, thus, to truncate the chain. Indeed, if we truncate the system of equations at some n max , then the equations of motion for E (n) , B (n) , and G (n) with 0 ≤ n ≤ (n max − 1) retain the form of Eqs. (20)- (22) and the last three equations for E (nmax) , B (nmax) , and G (nmax) take the following forṁ IV. NUMERICAL RESULTS The full system of equations describing the generation of gauge fields during axion inflation in the gradient expansion formalism consists of the Friedmann equation (9) for the scale factor, the Klein-Gordon equation (12) for the inflaton, the system of equations (20)- (22) for bilinear electromagnetic quantities, which must be truncated at some order n max , and the equation for the energy density of charged particles produced by the Schwinger effect (32). In this section, we present our numerical solutions to this system of equations, analyze their accuracy, and compare them to existing results in the literature.
We consider the generation of the Standard Model hypercharge gauge field axially coupled to the pseudoscalar inflaton field. The axial-coupling function is taken in the simplest linear form with one dimensionless coupling parameter β. Note that the coefficient in front of φ in the axial coupling function Eq. (66) is often parametrized as α/f , where f is the axion decay constant and α is another dimensionless parameter. However, we consider a generic axionlike inflaton field and do not specify its decay constant f ; therefore, parametrization (66) is more convenient. For numerical analysis, we take the inflaton effective potential in a simple parabolic form with M = 6 × 10 −6 M P . Although such a potential is already discarded by the CMB observations because it predicts very large magnitude of the tensor-to-scalar power ratio [98], it is still worth considering because many other inflaton potentials can be approximated by Eq. (67) close to their minima. This region appears to be the most important for magnetogenesis. Indeed, the most intensive generation of gauge fields occurs during the last few e-foldings of inflation (this is because the generation is determined by the parameter ξ ∝φ and the inflaton velocityφ typically is the largest close to the end of inflation). The initial condition for the inflaton field φ(0) = 15.55 M P was chosen to provide at least 60 e-foldings of inflation; the initial value of the inflaton velocity was determined from the slow-roll approximation,φ(0) = − 2/3M M P . Zero initial values for all electromagnetic quantities E (n) , B (n) , and G (n) as well as for the energy density ρ χ of charged fermions were assumed. Moreover, we suppose that the gauge field was also absent before the initial moment of time; consequently, the initial value of the damping parameter ∆ [see Eq. (47)] equals unity. We evolve our system of equations from the initial time deep in the inflation stage until its end, where preheating starts and the approximation of a homogeneous inflaton background breaks down. For definiteness, we assume that inflation terminates when the accelerated expansion of the Universe ends, i.e., when the conditionä = 0 is satisfied for the first time.
Numerical values of the parameter β typically considered in the literature are of order 10-100. The lower bound is rarely taken to be less than 10 because the gauge field production is very weak in this case. Concerning the upper bound, it is constrained from the requirement that the generated fields do not modify significantly the primordial power spectra at the scales relevant for CMB. The most stringent limit follows from the non-Gaussianity in the scalar power spectrum which results in |ξ| < 2.5 at the time when the modes relevant for CMB cross the horizon during inflation, i.e., δN e = 50 − 60 e-foldings before the end of inflation (see, e.g., Ref. [99]). For the m 2 φ 2 /2 inflationary model in the slow-roll approximation, |ξ| ≃ β/ √ 2 + 4δN e that immediately implies β max = |ξ max | √ 2 + 4δN e = 25−27 for δN e = 50 − 60, respectively. Consequently, in our numerical analysis we will restrict ourselves to the values of β in the range 10 ≤ β ≤ 25.

A. Small coupling and no Schwinger effect
We begin our analysis with the case of small coupling parameter β = 10 in the absence of the Schwinger effect. In such a case, the energy density of the generated gauge field is much less than that of the inflaton; i.e., the gauge field does not backreact on the inflaton evolution and the Universe expansion. This fact allows us to compute the energy densities of the generated gauge fields by an alternative method by using the spectra of generated fields. Indeed, in the absence of the backreaction and Schwinger effect, all Fourier modes of the gauge field evolve independently. Knowing the time dependences of the background inflaton field and the scale factor, we can solve the mode equation (40) for all physically relevant modes and get the power spectra of the generated fields at any moment of time t [see, e.g., panels (a), (d), and (g) in Fig. 1]. Then, integrating these spectra over the range of modes which crossed the horizon before the moment of time t during inflation, we get the time dependences of all electromagnetic quantities. These results are obtained without any approximation and, therefore, we use them as reference solutions in order to estimate the accuracy of the results of our gradient expansion formalism. Figure 2 illustrates the convergence of the solutions of the system of equations (20)-(22) truncated at a certain order n max to the exact mode by mode solution with an increase of n max . Panels (a), (b), and (c) depict the evolution of the electric, magnetic energy densities, and the scalar product | E · B |, respectively, during a few last e-foldings of inflation (its end is marked by the gray vertical lines on each panel). The dashed curves of different colors show the results for six lowest values of n max while the black solid lines correspond to the reference solution. Panel (d) shows maximal relative errors of the approximate results compared to the mode-by-mode solutions as functions of n max . For any electromagnetic quantity X, the error is defined in a usual way: where X ref is the corresponding reference value. The maximal absolute values of these quantities for X = ρ E , ρ B , and E · B achieved during inflation are shown in Fig. 2(d).
First of all, we should mention that there is a general tendency for the error to decrease with increasing n max . However, this decrease is nonmonotonic; for some specific values of n max (e.g., 9,17,19) the error is less than for both neighboring orders; for some other values (e.g., 8,18,20) the situation is opposite. Finally, for n max 25 the error ceases to decrease and stabilizes at the level of 0.6% − 0.8%. This irremovable error is not caused by the truncation procedure as it does not depend on the truncation order. It comes from other approximations used in the derivation of the boundary terms; e.g., from assumptions of constant ξ and exact de Sitter universe expansion in Eq. (44). This can be confirmed by the fact that the residual error increases during the last e-foldings of inflation (see, the bottom plot in Fig. 3) when the above mentioned approximations become less justified. Nevertheless, the achieved accuracy is absolutely satisfactory for all physical applications of the results.
The energy densities of the generated gauge field during the last few e-foldings of inflation compared to the total energy density of the Universe are shown in the top plot of Fig. 3. They were calculated in the gradient expansion formalism with n max = 30, and are in good accordance with the exact results during the whole inflation stage (see the bottom plot of Fig. 3 for the corresponding relative error of the result). This figure confirms that the generated fields are weak (the corresponding energy density is 4 − 5 orders of magnitude less than that of the inflaton) and cannot backreact on the Universe expansion. Close to the end of inflation, the electric and magnetic fields become almost collinear. This follows from the relative position of the curves for ρ E (blue solid line), ρ B (red dashed line), and 1 2 | E · B | (green dashed-dotted line). In the logarithmic scale, the third one is located almost in the middle between the first two; i.e., | E · B | ≈ E 2 B 2 (more precisely, cos θ grows from 0.86 to 0.96 during the last two e-foldings of inflation). This is a consequence of the fact that in the process of enhancement the gauge field becomes more helical and estimate (26) becomes valid.

B. Large coupling and no Schwinger effect
Let us now consider larger values of the coupling parameter, β = 20 and 25, still in the absence of the Schwinger effect. In this case, the generated gauge field backreacts on the background evolution. First of all, the term on the right-hand side of Eq. (12) becomes important and slows down the inflaton rolling. Second, the energy density of the gauge field becomes comparable to that of the inflaton and impacts the expansion rate of the Universe. Consequently, we cannot apply the standard spectral approach and find the solution of the mode equation (40) for each separate mode on a given inflaton background because this background itself depends on the resulting gauge field. Thus, all modes of the gauge field become coupled due to backreaction. One possible solution to this problem is to apply the iterative procedure which converges to a self-consistent configuration of the gauge field and the inflaton [60].
Another possibility is to use the gradient expansion formalism since it takes into account all relevant modes of the gauge field at once. As already emphasized, its main advantage is that it does not require any iterative procedure and allows us to get the results after a single numerical run. First, we find numerically the truncation order n max starting from which the solution of the truncated system of equations (20)- (22) does not change with increasing n max . Thus, we minimize the error caused by the truncation procedure. As a result, we obtain time dependences of all electromagnetic bilinear functions as well as the inflaton and scale factor. Second, we use the last two quantities to solve the mode equation (40) and get the gauge field power spectra [for β = 20 see, e.g., panels (b), (e), and (h) in Fig. 1]. Finally, integrating these spectra over the momenta of all relevant modes, we compute the electromagnetic bilinear functions. These solutions cannot be considered as the exact ones because they are based on the background values of the inflaton and scale factor which have been found by an approximate method. Nevertheless, we can use them as reference solutions in order to check the consistency of the gradient expansion formalism. Quantitatively it can be characterized by the relative deviation between the two solutions given by Eq. (68), where X is the solution found from the gradient-expansion approach and X ref is the corresponding result calculated from the gauge-field spectrum.   For β = 20, this deviation is always less than 0.6%. This is also true for β = 25 until the last e-folding of inflation where the deviation increases to ∼ 2%. For the majority of physical applications such an accuracy of the solution should be absolutely sufficient. However, it can be further increased applying the iterative procedure of Ref. [60].
In order to analyze the evolution of the system in more detail, we show the total electromagnetic energy density (red line), the scalar product | E · B | (blue line) and the parameter ξ, which is responsible for the gauge field generation, in Fig. 5. The dashed lines of the same colors represent the corresponding quantities in the absence of backreaction (when the inflaton background remains unperturbed). At first, the gauge field energy density monotonically increases in time because the parameter ξ constantly grows during inflation in the absence of the backreaction. However, when the impact of generated fields on the inflaton dynamics becomes important, the growth of the energy densities abruptly stops and the gauge field evolution enters the nonlinear stage. In this strong backreaction regime, the energy densities, the scalar product E · B , as well as the parameter ξ exhibit oscillatory features previously reported in the literature [51,60,100] (for easier comparison with the previously reported results, we plotted Fig. 5 for the same values of the parameters and in the same style as Fig. 6 in Ref. [60]). Moreover, since the inflaton rolls much slower, the inflation stage becomes a few e-foldings longer, which is clearly seen from Figs. 4 and 5, where the actual end of inflation [determined from the conditionä(t) = 0] is shown by the pink vertical lines and the end of inflation without the backreaction is marked by the gray vertical lines. Our gradient expansion formalism allows to recover all features of the backreaction regime previously reported in the literature.  (40), presented in Ref. [60], cf. Fig. 6 there.

C. Small and large coupling with the Schwinger effect included
Finally, let us consider the full physical system which includes also charged fermions produced due to the Schwinger effect during axion inflation. In such a case, even in the absence of backreaction the gauge field evolves in the nonlinear regime, because the Schwinger conductivity entering the mode equation (40) depends on the total electromagnetic field; i.e., it couples all relevant modes to each other. This does not allow us to determine the exact solution for the generated gauge field by solving the mode equation (40) separately for each mode as we did in Sec. IV A. This can be done, e.g., by using the iterative approach which still has not been realized in the literature. Instead, we perform the same procedure as discussed in Sec. IV B; i.e., we take the stable solution obtained from the gradient expansion approach and use it to solve the mode equation (40). Electromagnetic bilinear quantities then can be found from the spectra and compared to the results of the gradient expansion formalism. The relative deviation between them, given by Eq. (68), characterizes the consistency of our approach.
We present the results of our gradient expansion approach in Figs. 6 and 7 for β = 10 and 25, respectively. Panels (a) show the dependences of the electric (blue solid lines), magnetic (red dashed lines) energy densities, and the scalar product 1 2 | E · B | (green dashed-dotted line) on the number of e-foldings close to the end of axion inflation. Bottom plots show relative deviations from the reference solutions found from the spectra. As in the previous subsections, these deviations are of order 1% meaning that our gradient expansion formalism allows to get a highly self-consistent result even in the most complicated case when both backreaction and Schwinger effect are taken into account. Moreover, such an accuracy is achieved in a single run of the code without iterative procedure.
The top plots in Figs. 6(b) and 7(b) represent different components of the energy density of the Universe: electric (blue lines), magnetic (red lines), charged particles (green lines), and the total energy density (purple lines). Dashed lines of the respective color show similar dependences in the absence of the Schwinger effect (considered in more detail in previous subsections). Comparing the solid and dashed lines it is easy to conclude that the Schwinger effect significantly suppresses magnetogenesis. In the case of large coupling parameter β = 25, the gauge-field energy density becomes three orders of magnitude less than in the absence of the Schwinger effect and does not cause the backreaction any more. As a result, the inflation stage has the same duration as in the unperturbed case. Therefore, for typical values of the coupling parameter β considered in the literature, the backreaction does not occur in the presence of the Schwinger effect. Unless we take extremely large values of β or larger values of H, we do not expect any backreaction.
Especially strong suppression occurs for the electric component of the energy density, because it is directly affected by the conductivity, see Eq. (20). Although it dominates over the magnetic one during almost the whole inflation stage, at the end of inflation it rapidly decreases and becomes subdominant. This can also be seen from the evolution of spectra in Fig. 1 (c), (f), and (i). The electric energy density is transferred to charged fermions produced by the Schwinger effect. Indeed, the green curve corresponding to charged particles looks like a continuation of the blue curve representing the electric energy density. Already from Figs. 6(a) and 7(a) we can notice that 1 which is a signature of the fact that the electric and magnetic fields are nearly collinear. In order to check this more precisely, we plot the cosine of the angle between the electric and magnetic fields in the bottom plots in Figs. 6(b) and 7(b). They show that during the final part of inflation when the most significant generation of gauge fields occurs the cosine is greater than 0.8 and approaches unity at the end of inflation. This means that the angle between E and B is typically smaller than 35 o and decreases in time. This justifies our use of Eq. (23) for the Schwinger conductivity which was derived in the case of collinear electric and magnetic fields.

V. CONCLUSION
The explosive production of gauge fields during axion inflation is a complicated nonperturbative process. Its theoretical description becomes particularly challenging in the presence of nonlinear effects such as (i) the backreaction of the produced gauge fields on the evolution of the inflaton field and (ii) the Schwinger pair production of charged particles in the strong gauge-field background. In this paper, we have presented a novel gradient expansion formalism that is capable of successfully tackling this theoretical challenge, putting us in the position to qualitatively study and quantitatively describe the dynamics of gauge-field production during axion inflation at an unprecedented accuracy.
The typical approach to particle production during axion inflation, as it can be commonly found in the literature, consists in analyzing the gauge-field mode functions in Fourier space. This approach works well in the absence of backreaction or Schwinger pair production, in which case all Fourier modes evolve independently in time. However, as soon as at least one of these two nonlinear effects becomes relevant, the evolution of the individual Fourier modes can no longer be disentangled. All modes are coupled to each other, which complicates the theoretical and, in particular, numerical description of the system. In this case, backreaction can be accounted for by an iterative procedure that needs to be repeated until the numerical solutions for the gauge-field mode functions converge to a self-consistent result [51,60,100]. The Schwinger effect, possibly in combination with backreaction, can in principle be treated in the same way, although to the best of our knowledge no such self-consistent iterative analysis in Fourier space has thus far appeared in the literature. In contrast to this approach in Fourier space, the main idea behind our gradient expansion formalism is to study the evolution of a set of bilinear electromagnetic functions in position space. These functions are defined in terms of scalar products of two electric or magnetic field vectors with an arbitrary power of the curl operator. Thus by construction, they automatically include all physically relevant gauge-field modes that experience a tachyonic instability and are hence excited above the vacuum level during axion inflation. This characteristic property is extremely helpful in the presence of mode coupling because of nonlinear effects such as backreaction and Schwinger pair production. In this case, our bilinear functions allow us to capture the dynamics of the entire system, without the need for explicitly disentangling the evolution of the individual Fourier modes.
In this paper, we derived and solved a set of evolution equations for the bilinear electromagnetic functions in the gradient expansion formalism. These equations are coupled into an infinite chain, which, however, can be truncated at some point based on well-justified physical arguments. During inflation, the number of modes that leave the horizon, undergo the quantum to classical transition and become tachyonically unstable constantly grows. This effect is accounted for in our evolution equations by nontrivial boundary terms. We derived new explicit expressions for them in terms of the Whittaker function and its derivative as well as accurate approximate expressions that are convenient for numerical applications. A novel feature of the boundary terms, which can be regarded as vacuum source terms in our evolution equations, is the fact that they are proportional to a new parameter ∆. This parameter is related to the time integral over the time-dependent conductivity and describes the exponential damping of gauge fields deep inside the horizon caused by the presence of the conductive plasma. The state of the vacuum fluctuations inside the horizon thus depends on the evolution history at all earlier times, which renders the description of the system nonlocal in time and which led us to the notion of a time-dependent damped Bunch-Davies vacuum. As we were able to demonstrate, our gradient expansion formalism results in highly accurate and self-consistent solutions with a remaining numerical uncertainty of at most one to two percent. Moreover, it does not require an iterative procedure; all results can be generated in a single numerical integration of our differential equations.
While our formalism does not rely on any kind of spectral information in Fourier space, it can certainly be used as a starting point for solving the equations of motion for the gauge-field mode functions in Fourier space and hence derive approximate solutions for the energy spectra of the electric and magnetic fields. To this end, one only needs to use the numerical output of our formalism, specifically, the time evolution of the Hubble rate H, conductivity σ, and inflaton field φ, as input in the gauge-field mode equation (40). As shown in the previous section, this procedure results in highly accurate self-consistent energy spectra (see Fig. 1), which one can use for further phenomenological applications. This includes, e.g., an accurate calculation of the primordial scalar and tensor power spectra generated during axion inflation, notably, in the presence of nonlinear effects such as backreaction and Schwinger pair production. The first analysis along these lines, focusing on the scalar power spectrum and neglecting the Schwinger effect, has been carried out in Ref. [60]. An important outcome of this analysis was that backreaction can cause nontrivial features in the scalar power spectrum, which are related to the nontrivial evolution that we see in Fig. 5. In this paper, on the other hand, we were able to show that, at fixed coupling strength β, backreaction effects can be suppressed because of efficient Schwinger pair production, see Fig. 7. An interesting open question therefore is whether backreaction effects ever have a chance to leave a noticeable imprint in the primordial scalar and tensor power spectra in the presence of efficient Schwinger pair production, if one is willing to consider significantly larger values of β. A numerical study of the strong-coupling regime, β ≫ O(10), is, however, technically challenging, which is why we leave it for future work.
Another possible application of our formalism consists in a refined description of baryogenesis from hypermagnetic fields [52,59,[101][102][103]. This baryogenesis mechanism relies on the observation that helical hypermagnetic fields generated in the early Universe, possibly during axion inflation, can source a primordial baryon asymmetry around the time of the electroweak phase transition by virtue of the chiral anomaly of baryon number in the Standard Model (see also Ref. [104] for an extended scenario, which in addition involves right-handed neutrinos). The most important input quantity for these baryogenesis scenarios is the amount of hypermagnetic helicity generated during primordial magnetogenesis. Based on our formalism, this quantity can now be accurately calculated as a function of the model parameters of axion inflation and in the presence of nonlinear effects. Finally, our formalism is of course useful with regard to the phenomenology of primordial magnetic fields on their own, which we briefly discussed in the Introduction. For a given model of axion inflation, our formalism can be used to determine the electric and magnetic power spectra towards the end of inflation, which can then serve as the starting point for studying the further evolution of these fields during reheating and beyond. In particular, our formalism can generate input spectra for magnetohydrodynamics simulations in the radiation-dominated era, if the impact of reheating on our spectra is assumed to be small.
In the numerical analysis in this paper, we mostly considered a particular toy model in which axion inflation towards the end of inflation can be described by a simple quadratic mass term around the origin in field space. In the next step, it will be interesting to extend this analysis beyond this simple model and study the generation of the electric and magnetic fields during axion inflation in a model-independent way. In the literature, this is typically done in terms of two effective parameters: the inflationary Hubble rate H and the parameter ξ, which quantifies the inflaton velocity during inflation in Hubble units. An important lesson from our analysis, however, is that such a description will no longer be possible in the presence of the Schwinger effect. In this case, one will need to consider the dependence on at least one more parameter, namely, the parameter ∆, which describes the damping of the Bunch-Davies vacuum inside the Hubble horizon. We will return to this important question in future work. Sum over all fermions: 41 6  hypergeometric function U as follows: W κ,µ (y) = e −y/2 y µ+1/2 U (µ − κ + 1/2; 1 + 2µ; y).
Alternatively, the solution to Eq. (44) satisfying the boundary condition (45) can be expressed in terms of the Coulomb wave functions (see Sec. 14 in Ref. [105]) The boundary terms can also be expressed in terms of these functions as follows: where F s = F s (λξ, r(ξ, s)) and all other functions have the same arguments; r(ξ, s) is defined after Eq. (58) and prime denotes the derivative with respect to r, i.e., F ′ s = ∂F s (λξ, r)/∂r. The functions E λ , G λ , and B λ which determine boundary terms are expressed in terms of the Whittaker functions and are given by Eqs. (56)- (58) [or by Eqs. (B6)-(B8) in terms of the Coulomb wave functions]. These expressions, however, are inconvenient for numerical computations for big values of ξ because they require a significant increase of the working precision. Therefore, we found approximate expressions for these functions in terms of elementary functions which provide the sufficient accuracy for large values of |ξ| and all relevant values of parameter s.
We start with the expressions for s = 0. In such a case, there exist power series for the Whittaker functions (or the Coulomb wave functions) in inverse powers of |ξ|. For λ = sign ξ, it is convenient to work with the Coulomb wave functions. The asymptotical expansions for F 0 (|ξ|, 2|ξ|), G 0 (|ξ|, 2|ξ|), and the corresponding derivatives can be found from the integral representations of the Coulomb wave functions by the method discussed in Ref. [106] -for a few first terms in these expansions, see Eqs. (14.5.10) and (14.5.11) in Ref. [105] and also Refs. [107,108]. Applying this method and using Eqs. (B6)-(B8), we finally get the expansions which reproduce the exact result with a relative error less than 10 −4 % for |ξ| ≥ 3, E sign ξ (ξ, 0) = ( 3 2 ) 1/3 Γ 2 ( 2 3 ) π |ξ| 1/3 −  For λ = −sign ξ, asymptotical expansions for the Whittaker functions can be obtained by the method discussed in Sec. 6.13.3 of Ref. [109] and Sec. 10 of Ref. [110]. Approximate expressions for the functions E λ , G λ , and B λ with an error less than 10 −4 % for |ξ| ≥ 3 have the form For s = 0, it is much more difficult to get similar expressions because the second parameter appears in the expansion. We used the results of Ref. [111] and modified the first few terms of our previous expansions for s = 0 in order to obtain the approximate expressions reproducing exact results with an error of less than 0.5% for |ξ| ≥ 4 and all values of the parameter s. For λ = sign ξ, we obtain E sign ξ (ξ, s) ≈ ψ 1/3 r 1/3 where r = r(ξ, s) ≡ |ξ| + ξ 2 + s 2 + s, ψ = ψ(ξ, s) ≡ 2 ξ 2 + s 2 + s |ξ| + ξ 2 + s 2 + s . (B21)