{'fraction_non_alphanumeric': 0.03317269076305221, 'fraction_numerical': 0.011004016064257029, 'mean_word_length': 4.6502155661447695, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gradient expansion formalism for magnetogenesis in the kinetic coupling model

In order to describe magnetogenesis during inflation in the kinetic coupling model, we utilize a gradient expansion which is based on the fact that only long wavelength (superhorizon) modes undergo amplification. For this purpose, we introduce a set of functions (bilinear combinations of electromagnetic fields with an arbitrary number of spatial curls) satisfying an infinite chain of equations. Apart from the usual mode enhancement due to interaction with the inflaton, these equations also take into account the fact that the number of relevant modes constantly grows during inflation. Truncating this chain we show that even with a relatively small number of equations it is possible to describe the electric and magnetic energy densities with a few percent accuracy during the whole inflation stage. We arrive at this conclusion for different types of coupling functions (increasing, decreasing, and non-monotonic) in the regime with strong backreaction and its absence.


I. INTRODUCTION
Analysis of the gamma-ray observations from distant blazars indicates the presence of magnetic fields in intergalactic voids of our Universe [1][2][3][4][5][6][7]. The very large coherence length of these fields, measured in Megaparsecs, strongly favors their cosmological nature (see Refs. [8][9][10][11][12][13][14] for a review of different magnetogenesis models). As it is well known, one of the most important successes of inflationary models is the explanation of the cosmic microwave background (CMB) power spectra and the observed large-scale structure of the Universe which arose as a result of the evolution of primordial quantum perturbations generated during inflation [15][16][17] (for a review, see Refs. [18,19]). It is natural to assume that primordial magnetic fields were also formed at that epoch. In such a case, primordial magnetic fields are very important observables, because they, in addition to the abundance of chemical elements, the CMB radiation, and the large-scale structure of galaxies and their clusters, can be a unique source of information about physical processes in the early Universe at high energies, which are unattainable in the laboratory.
One of the key requirements necessary for generating magnetic fields during inflation is the violation of conformal invariance in the electromagnetic sector [20]. This can be done through the coupling of the electromagnetic field to the inflaton field or to the gravitational sector [21][22][23][24]. One of the most popular models in the literature is the kinetic coupling model where the inflaton couples to the square of the electromagnetic field tensor. It was first proposed by Ratra [22] and revisited many times in the literature [25][26][27][28][29][30][31][32][33][34][35][36].
In the standard approach to the kinetic coupling model, one works with separate Fourier modes of the electromagnetic field which evolve in the predefined background of the scale factor and the inflaton field (see, e.g., [27,29]). This approximation is valid only in the weak field regime and is based on the two facts: (i) the equation of motion for the mode function is a linear homogeneous equation [see Eq. (24) below] and (ii) the kinetic coupling to the inflaton does not mix different modes. Sufficiently strong generated electromagnetic fields cause the backreaction on the universe expansion and the inflaton evolution [28,35]. This in turn affects the evolution of all electromagnetic modes. Thus, the system enters a strongly nonlinear regime. Moreover, some other strong-field effects may become relevant, e.g. the Schwinger pair production, which also make all electromagnetic modes coupled to each other [37]. In these cases, one has to consider all relevant modes simultaneously. This means that one should solve the system of all mode equations together with the background equations that is numerically very costly [38,39]. Alternatively, one may use a lattice simulation. However, during inflation the spacetime is expanding exponentially and one should take the exponentially small lattice spacing at the beginning in order to get satisfactory accuracy at the end. That is why the lattice simulations are successfully used to describe only the last few e-foldings of inflation and preheating (see, e.g., Ref. [40] for the case of axial coupling). In the present paper we utilize a different approach based on the gradient expansion formalism. We work with vacuum expectation values of an infinite set of bilinear electromagnetic functions with an arbitrary number of spatial curls. Thus, the coordinate dependence of electromagnetic fields is traded for an infinite set of spatially homogeneous quantities. The same method was previously applied in the axial coupling model by three of us in Ref. [41]. Its extension to the case of the kinetic coupling model provides the main motivation for the present paper. This paper is organized as follows. In Sec. II we formulate the problem and apply the gradient expansion in order to derive an infinite set of equations for the bilinear electromagnetic quantities. Further, in Sec. III we consider the evolution of the electromagnetic mode function and separate the physically relevant modes which contribute to the generated electromagnetic field. Then, in Sec. IV we derive the boundary terms which take into account the variable number of physically relevant modes. Finally, in Sec. V we discuss our truncation scheme for the infinite chain of equations and then provide the numerical results obtained from such a truncated system. Section VI is devoted to conclusions. Appendices A and B provide some supplemental mathematical material which is used in the main text for the derivation of the boundary terms. Throughout the work we use the natural units and set = c = 1.

II. GRADIENT EXPANSION FORMALISM
The action describing the electromagnetic field A µ kinetically coupled to the inflaton field φ reads as where g = det g µν , g µν is the metric tensor, V (φ) is the effective potential of the inflaton, f (φ) is the kinetic coupling function, and F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic field tensor. The equation of motion for the inflaton field has the form while the Maxwell equations read as They are supplemented by the Bianchi identities whereF µν = 1 2 η µνλρ F λρ is the dual electromagnetic tensor, η µνλρ = ε µνλρ / √ −g, and ε µνλρ is the totally antisymmetric Levi-Civita symbol with ε 0123 = +1.
Throughout this paper we work in a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe whose metric is expressed in terms of the cosmic time as follows: Here a(t) is the scale factor which describes the universe expansion. It is convenient to use the Coulomb gauge for the electromagnetic field A µ = (0, A) where div A = 0. Then electric and magnetic fields are defined as follows: where overdot denotes the derivative with respect to the cosmic time. The components of the electromagnetic tensors are expressed through the electric and magnetic fields as usual where ε ijk is the three-dimensional Levi-Civita symbol and indices i, j, k on the right-hand side denote the components of 3-vectors. We would like to note that the electric and magnetic fields defined in Eqs. (6) and (7) are physical fields measured by a comoving observer.
Then, the Maxwell equations for the electric and magnetic fields take the forṁ div B = 0, div E = 0.
It is conventional to assume that the inflaton field is spatially homogeneous. In this case we havë The evolution of universe is driven by the total energy density of all fields which can be calculated as the 00component of the stress-energy tensor. It is defined as usual For a spatially homogeneous inflaton field, the energy density reads as Then, the Friedmann equation governing the evolution of the scale factor takes the form Observable classical electromagnetic fields are created during inflation due to enhancement of vacuum fluctuations of the quantum electromagnetic field. As we will see below, when the universe expands quasiexponentially during inflation, the wavelength of electromagnetic modes becomes much larger than the observable region. Such modes look like some long-range fields and can be treated classically. Nevertheless, their quantum nature implies that they are chaotically oriented in different regions of the universe. Therefore, it is not convenient to work with vector quantities such as E and B. Consequently, it is convenient to introduce a set of bilinear functions of the electromagnetic field where . . . denote the vacuum expectation value. Using Eqs. (8)-(9), we derive the system of coupled equations for these quantitiesĖ These equations describe the evolution of quadratic electromagnetic correlators (15)- (17) which contain the contributions from all Fourier modes of the electromagnetic field. As we will see in Sec. IV, a naive computation of the vacuum expectation values gives the UV-divergent results for them. This is a usual situation in quantum field theory, because the quantum fluctuations of all possible modes are always present in physical vacuum and their total energy is infinite. Therefore, in order to define physically meaningful quantities, we should take into account only modes whose evolution differs from that of vacuum fluctuations. The separation of such modes is performed in Sec. III. Their number changes in time during inflation and this leads to an additional time dependence of electromagnetic quantities which can be described by means of boundary terms in the equations of motion. These boundary terms are derived in Sec. IV.

III. EVOLUTION OF MODE FUNCTIONS
In this section, we consider the evolution of mode functions of the electromagnetic field. For this purpose, we decompose the electromagnetic field operator (in the Coulomb gauge) over the set of creation and annihilation operatorŝ where A λ (t, k) is the mode function with the wavevector (momentum) k and polarization λ; ε λ (k) are two circular polarization 3-vectors, which have the following properties: The creation and annihilation operators have the standard commutation relations Our goal now is to describe the evolution of the mode function A λ (t, k). Substituting decomposition (21) into Eq. (8) we find the following equation: or in the conformal time dη = dt/a(t), Here primes denote derivatives with respect to the conformal time. Initially, both polarizations are present equally in vacuum state. Moreover, from Eq. (24) or (25), we conclude that mode functions of different polarizations satisfy the same equations of motion. Consequently, the mode functions of both polarization would be identical at any moment of time. Therefore, in what follows, we will omit the polarization index λ in the mode function. There are two completely different terms in the square brackets in Eq. (25) and depending on their relative magnitude two completely different regimes of evolution can take place. The transition between these two regimes occurs when both terms are equally important. By analogy with the description of primordial perturbations, we will refer to this moment as the horizon crossing.
We assume that the mode k crosses (exits) the horizon (not to be confused with the Hubble horizon) if the equality is achieved for the first time during inflation. For a general coupling function f , the right-hand side of Eq. (26) may be a non-monotonic function of time that is why it is important to take only the very first moment among all the solutions of Eq. (26). The wavenumber of the last mode which has crossed the horizon before a certain moment of time t is determined by the maximal value of the right-hand side of Eq. (26) in the interval [0; t]. This can be represented as follows: where f 1,2 are the first two functions from the set and stands for the upper monotonic envelope of a function F (t). Its definition, properties, and examples are given in Appendix A. For modes deep inside the horizon k 2 ≫ k 2 h , the Bunch-Davies boundary condition [42] is imposed This solution is realized in the absence of the kinetic coupling and describes vacuum fluctuations of the electromagnetic field. The influence of the coupling function becomes important only after the mode exits the horizon. Our aim is to determine the mode function at the moment of horizon crossing. For this purpose, we will use an approach similar to the Wentzel-Kramers-Brillouin method near the classical turning point. Namely, we will approximate the exact time dependence of the second term in brackets in Eq. (25) by some simple ansatz which is accurate for k ∼ k h and admits for an exact analytical solution to Eq. (25). Unknown integration constants will be fixed by matching this solution with the Bunch-Davies mode function (30) deep inside the horizon.
The essence of our approximation is the following. In general, f n are smooth functions and vary much slower than the scale factor during inflation. Indeed, they depend only on the inflaton, its time derivatives and the Hubble parameter which are changing very slowly in the slow-roll regime. In the vicinity of the horizon crossing we exploit this fact and replace f 1,2 (t) by their (constant) values at the moment of horizon crossing f 1,2 (t h ). Moreover, since the conformal time can be approximately expressed as η = −1/(aH) during inflation, Eq. (25) can be written as This equation is accurate in the vicinity of the horizon crossing and admits for an exact analytical solution which is expressed in terms of cylindric functions as follows (see Appendix B for details): Here J α and Y α are the Bessel and Neumann functions, respectively, and the index α equals In order to determine D 1,2 , we have to match the solution (32) with the Bunch-Davies boundary condition. Using the properties of the cylindric functions listed in Appendix B, we finally get the solution valid in the vicinity of the horizon crossing where H (1) α = J α + iY α is the Hankel function of the first kind. For further convenience, let us find also the time derivative of this solution where we used properties of the Hankel function considered in Appendix B. Finally, far outside the horizon k ≪ k h (t), the term k 2 /a 2 in Eq. (24) can be neglected and the latter can be rewritten in the form Its general solution is , where C 1,2 (k) can be determined by matching the solution (37) with Eqs. (34) and (35) valid at the horizon crossing. The resulting expressions are the following: Here t h = t h (k) is the moment of time when the mode with wavenumber k crosses the horizon; the values of f 1,2 and α should be determined also at this moment of time.
In order to illustrate our results, we consider the time evolution of the mode with k = 10 10 M p in the kinetic coupling model with the Ratra coupling function f (φ) = exp(βφ/M p ) for β = 10. In Fig. 1 we plot the square moduli of A(t, k) and D(t, k) multiplied by 2k in panels (a) and (b), respectively. We compare the exact solution of Eq. (24) determined numerically (blue solid lines) with three different approximate solutions derived in this section. Purple dotted lines correspond to the Bunch-Davies solution (30). As expected, they are in accordance with the exact solution only before the horizon crossing (shown by the vertical black dashed line). The approximate solutions (34)- (35) are shown by green dashed lines. They nicely fit the exact solution in the vicinity of the horizon crossing from both sides (inside as well as outside the horizon). Finally, the approximate solution (37) is plotted by red dashed-dotted lines and works well outside the horizon. Thus, in this section we determined the time evolution of the mode function and found relatively simple and general approximate analytic solutions valid at any stage of the mode evolution. We have shown that only the modes which cross the horizon during inflation are influenced by the kinetic coupling and contain the information about generated electromagnetic fields. On the other hand the modes lying inside the horizon correspond to pure vacuum fluctuations and must be excluded from the analysis.
The fact that we take into account only a finite range of modes and that this range changes during inflation leads to an additional time dependence of the quadratic electromagnetic correlators (15)-(17) which is not encountered in system (18)- (20). In order to fix this inconsistency we would have to introduce some additional terms in the evolution equations, known as the boundary terms. They will be derived in the next section.

IV. BOUNDARY TERMS
As we have shown in the previous section, the physical meaning for magnetogenesis have only those electromagnetic modes which cross the horizon during inflation. Moreover, any electromagnetic observable at the moment of time t is determined only by modes which crossed the horizon before this moment. This means that the number of modes contributing to physical observables changes during inflation. In order to take this into account, we need to determine spectral decompositions for all quadratic electromagnetic correlators (15)- (17).
Using the explicit expression for the field operator (21) and substituting it into Eqs. (15)- (17), we obtain the following decompositions of the vacuum expectation values: Formally, the integration is performed over all possible wavenumbers from 0 to infinity. Let us check that taking into account all subhorizon modes leads to divergent results. Deep inside the horizon, mode functions are described by the Bunch-Davies solution (30), which is valid during the whole inflation stage for modes with extremely large wavenumbers k max h ≪ k < +∞. For these modes, the integrands in Eqs.(40)-(42) are given by The UV asymptotical behavior of all three quantities is ∼ k 2n+4 so that all integrals are severely UV divergent. However, as we discussed in the previous section, mode functions with large momenta A(t, k ≫ k h ) are not affected by the presence of the kinetic coupling function f and do not undergo the quantum-to-classical transition during inflation. Therefore, we exclude such modes from the consideration. Thus, we define the physically meaningful quadratic electromagnetic quantities as those containing only the modes crossing the horizon during inflation, i.e., where we expressed the integrands in terms of the mode function A and the auxiliary function D. Despite the fact that the integration is now performed only over a finite range of modes, we will keep the same notations E (2n) , B (2n) , G (2n+1) since in the rest of the paper only the physical quantities will be considered. Another interesting question is whether the lower integration limit in Eqs. (44)-(46) is zero or it should acquire some finite value. The physical meaning of modes which are outside the horizon during the whole inflation stage is still under debate in the literature. Moreover, the duration of the inflation itself is unknown; there is only the minimal number of e-foldings needed to resolve some cosmological problems, typically N min ≃ 60. That is why the lower boundary can be a wavenumber somewhere between 0 and k h (t 60 ) where t 60 is the moment of time at 60 e-foldings before the end of inflation. Fortunately, for the majority of coupling functions, the spectrum of generated electromagnetic fields has a blue tilt (i.e., a positive spectral index) and the exact value of the lower integration boundary is not important for numerical computations. In the rest of the work we will take the zero value of the lower integration limit. Thus, any electromagnetic quantity X can be written in the following form: i.e., as an integral with variable upper boundary k h (t) which was defined by Eq. (27). However, Eqs. (18)-(20) take into account only the evolution of each mode separately and do not take into account the fact that the number of relevant modes changes in time. Therefore, we have to introduce additional terms to describe the modes which cross the horizon and start to contribute to electromagnetic function X: where the lower index "b" means "boundary term". In order to calculate the boundary terms, we need to know the mode function in the vicinity of the horizon crossing. In Sec. III we have already determined the approximate analytical expressions which work well in this region, they are given by Eqs. (34)- (35). Then, using Eq. (48) and spectra (44)-(46), we get the boundary terms for the quadratic functions in the following form: In the derivation of these boundary terms we made two approximations which may introduce numerical errors in the result: (i) we neglected the time dependence of functions f 1,2 near the moment of horizon crossing by an electromagnetic mode; (ii) we approximated the universe expansion close to that moment of time by the de Sitter solution. Both these assumptions work well during inflation for smooth enough coupling functions. However, close to the end of inflation their impact on the result may become significant. We will analyze errors of the numerical results in the next section for different coupling functions and show that typically errors do not exceed a few percent level. Finally, we get the following system of equations: This system of equations is infinite and in the following section we will use some approximations in order to truncate it.

V. NUMERICAL RESULTS
In the previous sections we derived an infinite chain of equations (54)-(56) for the set of bilinear functions of the electromagnetic field (15)- (17). However, their practical usage is impossible unless we find some physically meaningful way to truncate it. In this section, we will show that such a truncation indeed exists during the inflation stage and then we will present some numerical results obtained from the truncated system of equations for different types of coupling functions.

A. Truncation procedure
As we discussed in Sec. IV, at any given moment of time t the shortest mode in the electromagnetic spectra is determined by the horizon size k −1 h (t). This fact will allow us to express the higher-order electromagnetic quantities E , G , B through the lower-order ones. Let us consider, for example, the magnetic function B (2n+2) . From Eq. (46) we conclude that for a generic momentum dependence of the mode function A(t, k) there exists high enough order n for which the contribution to B (2n+2) is dominated by the ultraviolet part of the spectrum near k h . Let us define the spectral index p in this region by the relation |A(t, k k h )| 2 ∝ k p−1 (e.g., p = 0 for the Bunch-Davies solution). Then, we can find whereC(t) does not depend on n. Therefore, we can approximate The last approximation becomes more and more accurate for higher n. However, we will use it even at low orders because the spectral index p is unknown in general. Similar approximations will be used for E and G functions: These relations enable us an opportunity to truncate the infinite chain at certain finite order. Let us denote by N the greatest order we would like to consider. For even N , the system of equations terminates with Eqs. (54)-(55) where the above mentioned approximation for G (N +1) was useḋ For the lowest possible order N = 0, the system can be closed by setting G (1) = 0. Thus, in total, the system consists of (3N + 4)/2 equations for the electromagnetic functions and 2 equations for the scale factor (52) and the inflaton field (53). For odd N , the last equation in the system is Eq. (56) which takes the forṁ In this case the total number of equations is (3N + 3)/2 for the electromagnetic functions plus 2 background equations (52) and (53).
In the following subsections we provide numerical solutions of the truncated system of equations for different values of N from 0 to 15. In order to determine the accuracy of our approximation, we need to perform a direct and independent calculation of the electromagnetic quantities. This indeed can be done in the regime when the backreaction of generated fields can be neglected. In such a situation, background Eqs. (52) and (53) can be solved by setting the electromagnetic field to zero. Then, having determined the time dependences of the inflaton and scale factor, we solve Eq. (24) numerically for all modes crossing the horizon during inflation. Finally, integrating the spectral decompositions (44), (45), and (46) over the relevant range of modes, we get the exact results for electromagnetic quantities which can be used to estimate the accuracy of the approximate solutions of the truncated system of equations.
For definiteness, we consider in the Starobinsky R 2 -model of inflation [43] with the effective potential (in the Einstein frame) Predictions of this model for the spectral index of the primordial scalar perturbations and the tensor-to-scalar power ratio are in good accordance with the recent cosmic microwave background observations [44]. The amplitude of the potential, V 0 , is chosen to provide the correct value of the scalar power spectrum amplitude [44]. For simplicity, we assume that the CMB pivot scale exits the Hubble horizon 60 e-foldings before the end of inflation.

B. The Ratra coupling function
Here we present the numerical solutions of the truncated system of equations for the Ratra coupling function where β is a dimensionless coupling parameter. In order to test our formalism, we consider three different values of β which correspond to different regimes: β = −7 gives an increasing in time coupling function without the backreaction of electromagnetic fields on the background evolution; β = 7 corresponds to the decreasing coupling function again without the backreaction; finally, β = 15 gives the decreasing function with strong backreaction.

Increasing coupling function without backreaction
We start with the case of a monotonically increasing coupling function with β = −7. Since the corresponding coupling constant is always less than unity during inflation, the dynamics of charged fields becomes strongly coupled [29]. Nevertheless, we consider it in order to test the applicability of our approach to growing coupling functions. Figure 2 (a) shows the time dependences of the electric and magnetic energy densities obtained from the truncated system of equations with N = 15. For such a small value of the coupling parameter β, generated fields are weak, their energy density is several orders of magnitude less than that of the inflaton. The most effective generation occurs at the end of inflation when the inflaton gets out of the slow-roll regime. Although the electric and magnetic components are always of the same order of magnitude, the magnetic component is always slightly larger because the coupling function increases with time. The wavenumber of the mode crossing the horizon as a function of time k h (t) is shown by the red solid line in Fig. 2 (b). Since the argument of the envelope in Eq. (27) is a monotonically increasing function (shown by the black dotted line), the envelope is identically equal to the function itself. For comparison, we also plot in the same figure the time dependence of the Hubble horizon scale k H = a(t)H(t) (blue dashed line). For this coupling function, the horizon for electromagnetic modes and the Hubble horizon appear to be very close. However, the biggest deviations occur close to the end of inflation where the generation is the most effective. Therefore, it is important to use the correct definition of the horizon in order to properly extract all modes contributing to the generated electromagnetic field.
In Fig. 3 we compare the solutions of our system of equations truncated at different orders N with the result of exact numerical integration of Eq. (24) for all relevant modes. Since the backreaction is absent in this case, the latter exact solution can be used as a reference point for comparison. Panels (a) and (b) in Fig. 3 show the time dependences of the electric and magnetic energy densities during the final part of inflation. The results show that errors introduced by the truncation of the system and by the approximate form of boundary terms are accumulated and have maximal values at the end of inflation. Nevertheless, with increasing N the deviations from the exact result decrease and starting from N = 7 they become visually indistinguishable.  In order to quantify the error of the approximate result, we calculate the relative deviation from the exact solution at the end of inflation, t = t e , where X denotes any of the electromagnetic quantities, in particular energy densities, X N corresponds to the solution of the truncated system of order N while X exact is the corresponding exact solution. The errors of computation of the electric and magnetic energy densities are shown in Fig. 3(c) by the blue solid line and red dashed line, respectively.
There are a few features we would like to comment on. First of all, the general tendency is that increasing N typically leads to smaller errors. Second, the truncation at some even order gives worse results than at the neighboring odd order. However, for sufficiently large N (in this case, for N ≥ 12), the error stops decreasing and remains almost constant, of order few percent. This residual error cannot be explained by the truncation of the system of equations as it does not depend on the truncation order. Its possible origin is the approximate expressions for the boundary terms [see the discussion after Eq. (51)]. Fortunately, this error is very small and absolutely irrelevant for practical issues.

Decreasing coupling function without backreaction
Now let us consider the case of a decreasing in time Ratra coupling function with β = 7. It does not cause the strong coupling problem during inflation and does not lead to the backreaction. This allows us to use the exact numerical solution of the mode equation (24) to test the accuracy of our approximate solutions. Before doing this, let us discuss the qualitative behavior of generated electromagnetic fields in this model.
The time dependences of the electric and magnetic energy densities are depicted in Fig. 4 (a) by the blue solid and red dashed lines, respectively. The qualitative time dependence of the coupling function is sketched (also in logarithmic scale) by the black dotted line. The energy densities are monotonically increasing functions except a small period of time close to the end of inflation where they exhibit nonmonotonic behavior. The origin of this minimum can be understood from Fig. 4(b) where the time dependence of the horizon scale is shown by the red solid line. For positive values of β, the combination (f 1 + f 2 ), which determines the horizon scale (27), changes its sign during inflation. Indeed, in the slow-roll regime, we get The time derivative of the inflaton is always negative and grows in the absolute value during the Starobinsky inflation. Then, for sufficiently large positive β, the combination |(βφ)/(HM p )| may exceed unity and the combination f 1 + f 2 flips its sign. As a result, the argument of the envelope in Eq. (27)  line in Fig. 4(b). In the vicinity of this moment of time the horizon scale (the envelope) remains constant [see the inset in Fig. 4(b)]. Physically, this means that no new modes cross the horizon in this period of time. As a result, the universe expansion leads to a decay of the energy densities which is clearly seen from Fig. 4 44), (46). The curves corresponding to the lowest N are shown because they can be visually distinguished from the exact result.
Truncating our system at different orders N we get approximate solutions for the electric and magnetic energy densities which are shown in Fig. 5(a) and (b) by the dashed colored curves (see the legend). For comparison, the black solid curves show the exact solutions obtained from the numerical integration of the mode equation (24) for all relevant modes. For a few first values of N the approximate curves significantly differ from the exact result, however, this difference quickly becomes invisible for N ≥ 4. The dependence of the relative error of the approximate solution at the end of inflation (where the error is the largest) is qualitatively the same as in the case β = −7 shown in Fig. 3(c). Starting from N = 6, the relative errors both for the electric and magnetic energy densities become less than 10% and asymptotically tend to residual values of order few percent.

Decreasing coupling function with backreaction
Finally, let us consider the case with a large coupling parameter β = 15. At first, it has the same qualitative behavior as in the previously discussed case with β = 7. In fact, the coupling function is also decreasing in time, that is why the electric component slightly dominates the magnetic one, see Fig. 6(a). Both energy densities monotonically grow except the short period of time when they exhibit a sharp minimum near t ≈ 0.75 × 10 7 M −1 p . As previously, the decrease happens when the horizon scale becomes frozen for a short time near the point where f 1 + f 2 changes the sign, see Fig. 6(b). However, for β = 15, the production of the electromagnetic field is much more effective. As a result, its energy density becomes comparable to that of the inflaton near t = 10 7 M −1 p and the backreaction becomes important (more precisely, this happens when ρ EM ≃ ǫρ inf , where ǫ = (M 2 p /2)(V ′ /V ) 2 is the slow-roll parameter [35]). This drastically changes the evolution of the system. The electric energy density (the dominant component) becomes almost constant slowly decreasing while the magnetic one (the subdominant component) starts to decay very rapidly. This can be easily understood from the comparison of Eqs. (54) and (55) for n = 0. In fact, the backreaction slows down the inflaton rollingφ so that the Hubble damping term withφ in Eq. (54) is almost exactly compensated by the third term on the left-hand side with the coupling function (the last term 2G (1) as well as the boundary term are much smaller and can be neglected). On the other hand, in Eq. (55), the term with the coupling function enters with the opposite sign and thus helps the Hubble term to damp the magnetic component.  For β = 15, we can compare our approximate results to the exact numerical solution of the mode equation (24) only until the backreaction becomes important. Figure 7 compares the results of the system truncated at N = 0, 3, and 15 (solid lines) with the exact result (dotted lines). All three panels show a very good accordance when the backreaction is not important. As the latter becomes relevant, the exact mode-by-mode solution of Eq. (24) begins to grow rapidly and quickly becomes larger than the inflaton energy density. In this region the comparison is no more possible. Therefore, in Fig. 8, we compare the results of the truncated system for different N with the results for

C. Nonmonotonic coupling function
In this subsection we consider the nonmonotonic coupling function where β, C 0 , and φ 1 /M p are three dimensionless parameter which determine the steepness, height, and location of the maximum of this function. For C 0 > 0, it is always greater than unity which ensures that the strong coupling problem never occurs during inflation. The values of coupling parameters should be chosen in such a way that (i) the maximum of the coupling function is reached during the simulation time (so that we exploit its nonmonotonic behavior); (ii) the value f ≃ 1 is reached at the end of inflation (electric charges of all fields are not rescaled); (iii) C 0 is large enough to provide a significant generation of electromagnetic fields. We do not claim that such a function has some physical origin, but only use it to test the applicability of our formalism to more complicated (nonmonotonic) functions. We will see that the evolution of the electromagnetic energy density has much more interesting features and our truncated system already can reproduce all of them even for small N .
In the numerical analysis, we use two different sets of parameters which represent the regimes without backreaction and with it.

Nonmonotonic function without backreaction
In order to be able to compare our approximate results with the exact solution of Eq. (24), we first consider the coupling function which does not lead to the backreaction on the inflaton evolution. This is realized, e.g., for the following set of coupling parameters: β = 10, C 0 = 10 19 , φ 1 = 4.5 M p .
The time dependences of the energy densities in this case are shown in Fig. 9(a). As long as the function increases (it is sketched by the black dotted line), the magnetic component is slightly greater than the electric one. As the coupling function starts to decrease, their roles interchange. At first, the curves grow monotonically, but then two sharp minima appear. They happen when the combination f 1 + f 2 changes sign. For the nonmonotonic function (67), this happens twice: once close to its maximum and once during the decreasing stage. At these moments of time the horizon scale k h , which is the envelope of aH|f 1 + f 2 |, is frozen, see Fig. 9 the boundary terms vanish at that time. Therefore, the electromagnetic field undergoes short-term damping due to the universe expansion. The generation becomes more efficient when the inflaton starts rolling faster close to the end of inflation. However, this lasts until the hyperbolic cosine in Eq. (67) becomes greater than √ C 0 . After that the coupling function remains almost constant f ≈ 1, see the final horizontal part of the black dotted line in Fig. 9(a). At this stage, the horizon scale k h remains constant, see Fig. 9(b). As a result both the generation terms ∝ (ḟ /f ) and the boundary terms in Eqs. (54)-(56) vanish. Then, the Hubble expansion leads to a fast damping of generated fields. Let us now compare the results of the approximate system of equations truncated at different orders. The electric and magnetic energy densities are shown in Fig. 10(a) and (b), respectively. The dashed colored curves corresponding to different N are compared to the exact numerical result obtained from the mode-by-mode solution of Eq. (24) which is shown by the black solid line. As in the previously discussed cases with monotonic coupling functions, the dominant (electric) component is satisfactorily described already by the approximation with N = 2. However, the subdominant (magnetic) component starts to converge to the exact solution only after N = 4. The maximal error is, as usually, reached at the end of inflation. We show its dependence on the truncation order N in Fig. 10(c). Starting from N = 8 the errors for both electric and magnetic components is less than 10%. With further increase of N the errors stabilize and remain at the level of order 1%. This residual error is caused not by the truncation but by other approximations we did in the derivation of our system.

Nonmonotonic function with backreaction
Finally, let us discuss the situation when the nonmonotonic coupling function (67) leads to a strong backreaction of generated electromagnetic fields. For example, this happens for the coupling parameters β = 12, C 0 = 10 21 ,  Energy densities of the electric and magnetic fields generated in this case are shown in Fig. 11(a). Their time dependences have several interesting features. At first, when the function is almost constant and close to unity (see the horizontal part of the black dotted curve), generated fields are very weak. Further, the magnetic field slightly dominates over the electric one as the function grows. Then, when the coupling function starts to decrease, their roles change. A rapid decay of the magnetic field is caused by two factors: (i) the horizon scale becomes frozen twice for a short time (again, this happens when f 1 + f 2 changes sign), see Fig. 11(b), and the electromagnetic field is no longer enhanced by new modes crossing the horizon; (ii) the third term on the left-hand side in Eq. (55) leads to a decay of the magnetic field.
A rapid decrease of the coupling function generates very strong electric fields which then induce the growth of the magnetic field. This is due to the Faraday law which is incorporated in Eq. (55) by the last term on the left-hand side. This explains the second increase of the magnetic energy density. Finally, the total electromagnetic energy density reaches the value ρ EM = ǫρ inf and the backreaction occurs. The electric energy density, which is the dominant component, stops growing and remains almost constant (as a result of the dynamical equilibrium between the inflaton and the generated field). The magnetic field at this time decays due to the universe expansion. This analysis explains the "zig-zag" form of the red dashed curve in Fig. 11.
Numerical results for the energy densities obtained from the exact solution of Eq. (24) for all relevant modes naturally does not take into account the backreaction. Therefore, we can use it for comparison with our approximate results only until the backreaction is important. Such a comparison is shown in Fig. 12, where the solid lines correspond to the results for the truncated system of equations while the dotted lines show the mode-by-mode solutions. As it can be seen from this figure, the truncated system reproduces very well the exact solution until the backreaction turns on.
Because of the backreaction we cannot use the exact mode-by-mode solution as a reference for quantifying the error of our approximate results. That is why we use the result of the truncated system with N = 15 for this purpose. The latter is shown by the black solid line in Fig. 13 while the solutions of the system with smaller N are shown by the colored dashed lines (see the legend). Again, the dominant component of the generated energy density (the electric one) is reproduced with high accuracy even by the low-N approximations, while for the magnetic component the convergence starts from N ≥ 4.

VI. CONCLUSION
In this work, we extended the gradient expansion formalism (proposed in Ref. [41] for the axion coupling model) to the kinetic coupling model of the generation of electromagnetic fields during inflation. In contrast to the standard approach operating with separate Fourier modes of the electromagnetic field, we work with physical quantities in the coordinate space which include all physically relevant modes at once. Although, in principle, there are infinitely many of such quantities, at practical level it appears to be sufficient to use the finite number of them to satisfactorily describe the magnetogenesis process during inflation. We found that such a formalism works very efficiently and has an excellent convergence.
An important feature of inflationary magnetogenesis is the fact that the number of relevant modes grows in time as more and more new modes cross the horizon and undergo the quantum-to-classical transition. We take this into account by introducing the boundary terms in the equations of motion for the electromagnetic quantities. Their origin can be described as follows. Any of the electromagnetic quantities (15)-(17) is a result of integration of the corresponding spectral densities over the range of physically relevant modes. Since new modes cross the horizon, the higher integration boundary continuously grows in time during inflation and the electromagnetic quantities acquire additional time dependence. The latter is described by an additional term in the corresponding equation of motion. Thus, we realize the idea of the adaptive UV cutoff of the electromagnetic power spectrum.
The central result of our paper is the system of coupled equations (54)-(56) for the set of bilinear electromagnetic quantities E (n) , G (n) , B (n) which contain n spatial curls acting on the electric or magnetic field. The equation of motion for the n-th order quantity always contains at least one quantity with n + 1 curl. As a result, all equations are coupled into an infinite chain. There is, however, a physical argument allowing us to truncate this chain at some finite order. Indeed, for the quantities with a large number of curls the dominating contribution is made by the UV edge of the spectrum. The shortest relevant mode is determined by the horizon scale (27) and can be easily determined. This provides an opportunity to express the higher-order quantity in terms of the lower-order one and, thus, to truncate the chain.
As a testbed for our formalism we used the Starobinsky model of inflation which is one of the most favored by the cosmic microwave background observations [44]. We considered two types of the coupling functions: the Ratra function (64) and the nonmonotonic inverse cosine function (67). Choosing different sets of the coupling parameters we studied qualitatively different cases of the increasing, decreasing, and nonmonotonic coupling functions in the absence of the backreaction and taking it into account. We compared the numerical results for the electric and magnetic energy densities obtained from the approximate system of equations truncated at different orders N = 0 ÷ 15. For all types of the coupling function we found a convergence of the numerical results to some limiting solution which is independent on N . The residual error of this solution [compared to the exact result obtained from the numerical integration of the mode equation (24) for all physically relevant modes] has the maximal value of order few percent at the end of inflation. This error is caused not by the truncation of the system but by other approximations which were adopted in its derivation (an approximate solution of the mode equation near the moment of horizon crossing, the slow-roll approximation etc.). We observed that the approximate solution has a less than 10% error already for the truncation order N ≃ 10. At this order one has to solve ∼ 20 coupled ordinary differential equations which is a simple task even for an ordinary personal computer.
The nonmonotonic coupling function (67) leads to much more complicated time behavior of the energy densities than in the case of a monotonic Ratra function. At first, when the function grows in time, the magnetic component is slightly greater than the electric one. Further, when the coupling function starts to decrease, the roles of these components interchange. However, this does not have any important effect on the final value of the magnetic field. The reason for this is pretty clear: the most intensive generation occurs close to the end of inflation when the inflaton rolls faster. At that time the coupling function is already a decreasing function and, therefore, the generation occurs qualitatively in the same way as for a function monotonically decreasing from the very beginning. Nonmonotonic coupling functions were studied previously in Refs. [30,31]. The authors considered the piecewise smooth functions composed of 2 or 3 pieces with a power-like dependence on the conformal time f ∝ η α with different α. For a nonmonotonic function with one maximum, their conclusions are similar to ours: they do not observe a significant enhancement compared to the case of a monotonic function.
Our method fails to work after the end of inflation. First of all, this happens because the horizon scale stops increasing and no new modes exit the horizon. As a result, the enhancement occurs only for modes which crossed the horizon during inflation and the spectrum of such modes strongly depends on time. In such a case there is no simple relation between the quantities of different order n any more and we cannot truncate the chain of equations. There are, however, even more fundamental problems in the description of preheating. Firstly, at this stage, the inflaton cannot be treated as spatially homogeneous [40]. Secondly, one has to take into account the inflaton damping due to the perturbative and nonperturbative (e.g., parametric resonance) particle production during preheating [45,46]. The last but not the least, the Schwinger effect may become extremely important after the end of inflation. All these phenomena can be taken into account only in a complicated lattice simulation. Fortunately, during preheating the scale factor grows in time slowly (powerlike time dependence in contrast to the exponential one during inflation) and the lattice spacing remains of the same order of magnitude during the simulation. Indeed, this method has been already applied to the description of preheating, see, e.g., Ref. [40]. The lattice simulations require the initial conditions for all physical quantities imposed at the end of inflation (or a few e-foldings before its end). Our approach can be used to generate such initial conditions.
Another interesting question is about the role of the Schwinger effect during inflation. In the present work we completely neglected this phenomenon. However, it is an important problem to take it into consideration. In particular, it may be particularly important in the backreaction regime when the Schwinger process decreases the electric energy density and helps to end the inflation stage [35]. This issue deserves a separate investigation and will be addressed elsewhere. Their asymptotical behavior for large values of the argument is the following: From this expression we conclude that H (1) (x) matches the correct asymptotics of the Bunch-Davies solution (30) for −kη → +∞. Comparing the prefactors in Eqs. (30), (32), and (B3), we finally get solution (34). Below we list two properties of cylindric functions which were used in algebraic calculations in the main text. They are valid for any kind of cylindric functions Z α = {J α , Y α , H