Hypermagnetogenesis from axion inflation: Model-independent estimates

Axion inflation coupled to the Standard Model (SM) hypercharge gauge sector represents an attractive scenario for the generation of primordial hypermagnetic fields. The description of this scenario is, however, complicated by the Schwinger effect, which gives rise to highly nonlinear dynamics. Hypermagnetogenesis during axion inflation in the absence of nonlinear effects is well studied and known to result in a hypermagnetic energy density that scales like $H^4\,e^{2\pi\xi}/\xi^5$, where $\xi$ is proportional to the time derivative of the axion-vector coupling in units of the Hubble rate $H$. In this paper, we generalize this result to the full SM case by consistently taking into account the Schwinger pair production of all SM fermions. To this end, we employ the novel gradient-expansion formalism that we recently developed in [2109.01651], and which is based on a set of vacuum expectation values for bilinear hyperelectromagnetic functions in position space. We parametrize the numerical output of our formalism in terms of three parameters ($\xi$, $H$, and $\Delta$, where the latter accounts for the damping of subhorizon gauge-field modes because of the finite conductivity of the medium) and work out semianalytical fit functions that describe our numerical results with high accuracy. Finally, we validate our results by comparing them to existing estimates in the literature as well as to the explicit numerical results in a specific inflationary model, which leads to good overall agreement. We conclude that the systematic uncertainties in the description of hypermagnetogenesis during axion inflation, which previously spanned up to several orders of magnitude, are now reduced to typically less than 1 order of magnitude, which paves the way for further phenomenological studies.


I. INTRODUCTION
Baryonic matter in the Universe mostly exists in the form of plasma. Being composed of free-streaming charged particles, plasma very efficiently screens electric fields. On the other hand, its large electric conductivity keeps magnetic fields frozen for a long time. Therefore, it is not surprising that magnetic fields are observed everywhere in the Universe, namely, in stars, galaxies, and clusters of galaxies [1][2][3][4][5][6][7][8][9]. There exist several astrophysical mechanisms that could be responsible for the generation of magnetic fields on these length scales. In contrast, the evidence for magnetic fields in voids with coherent lengths of the order of megaparsecs based on the observation of blazars [10][11][12][13][14][15][16][17][18][19][20] is quite unexpected and fascinating. Indeed, the small matter content in voids makes the direct generation of magnetic fields in voids impossible. Although these fields could be induced by outflows of magnetized matter from galaxies [21][22][23][24], such outflows would need to be strong and coherent over tens of millions of years, which appears implausible. As a consequence, a cosmological origin of these magnetic fields emerges as an interesting and realistic possibility.
Inflationary magnetogenesis [25,26] naturally addresses the large coherence length of magnetic fields observed in voids. In addition, cosmological magnetic fields generated in the early Universe provide the necessary seeds for magnetic fields in protogalaxies, whose amplification through adiabatic compression [1] and different types of dynamo processes [27][28][29][30][31] could easily explain the magnetic fields observed in galaxies and clusters of galaxies today [32]. One of the most attractive features of inflation for cosmology is that it results in an isotropic and homogeneous Universe consistent with the smallness of the temperature fluctuations in the cosmic microwave background. Since magnetic fields are not enhanced in a conformally flat inflationary background [33], this means that the conformal symmetry of Maxwell's equations should be broken to ensure the possibility of inflationary magnetogenesis. Although this breaking can be done in many ways (see, e.g., Refs. [25,26,34,35] or Ref. [36] for a recent effective field-theoretical analysis), we consider in this work the axial coupling of the Standard Model (SM) hyperelectromagnetic field to a pseudoscalar "axion" inflaton field . This scenario results in helical hypermagnetic fields, which enhances the chance of their survival in the primordial plasma.
The production of hypercharge gauge fields during axion inflation is subject to several nonlinear effects, which highly complicates the theoretical analysis. These effects include (i) the backreaction of the produced gauge fields on the evolution of the inflation field [44,48,56,65] as well as (ii) the Schwinger pair production of hypercharged matter degrees of freedom [66][67][68]. The produced pairs of particles and antiparticles quickly form an ultrarelativistic plasma, which efficiently screens the electric field. This strongly affects the generation and evolution of electromagnetic fields, especially near the end of inflation and during reheating [49,53,55,. A quantitative description of hypermagnetogenesis during axion inflation is therefore theoretically challenging, which is why up to now only some order-of-magnitude estimates of its efficiency have been worked out in the literature [49,54].
The goal of this paper is to improve on this situation, leveraging the quantitative accuracy of the gradient-expansion formalism [53] that we recently successfully applied to axion inflation coupled to the SM hypercharge gauge field in Ref. [92]. As demonstrated in Ref. [92], this novel gradient-expansion formalism allows us to consistently account for the above-mentioned nonlinear effects, which enables us to evaluate the efficiency of gauge-field and fermion production during axion inflation at unprecedented accuracy. The basic idea behind the formalism is to consider vacuum expectation values of a truncated set of bilinear electromagnetic functions in coordinate space rather than momentum space. Solving the equations underlying our formalism in a single numerical run, we are able to describe the evolution of the electric and magnetic energy densities at percent-level accuracy during the whole inflation stage without the need for an iterative procedure. The formalism also takes into account the fact that the number of relevant gauge-field modes constantly grows during inflation as new modes become tachyonically unstable by adding appropriate boundary terms in the equations of motion for the bilinear electromagnetic functions.
In Ref. [92], we provided a detailed description of the gradient-expansion formalism and confirmed its validity by comparing its numerical output to existing results in the literature for specific model and parameter benchmark scenarios. In this paper, we shall continue our investigation of hypermagnetogenesis during axion inflation based on the gradient-expansion formalism, now turning to a model-independent analysis. In the following, we will study the efficiency of gauge-field production during axion inflation in terms of a minimal number of parameters (the gaugefield production parameter ξ, Hubble rate H, and damping factor ∆; see Sec. II for the precise definition of these quantities), which will provide us with numerical results that are applicable across a large range of models based on different types of scalar potentials. In fact, as we will show, our numerical results will always provide a good estimate of the efficiency of gauge-field production whenever the three parameters ξ, H, and ∆ vary only very slowly during axion inflation, such that their time dependence can be approximately neglected. To facilitate the application of our numerical results in future studies, we will also present semianalytical fit functions that reproduce our full numerical results to very good accuracy. These fit functions are compact and ready to use, which means that, in future studies, it will not be necessary to implement our full gradient-expansion formalism and redo the entire numerical analysis.
Finally, we will also compare the model-independent estimates in Refs. [49,54] to our new model-independent estimates and present fit functions for the estimates in Refs. [49,54]. An important outcome of this exercise will be that, while the estimates in Refs. [49,54] span several orders of magnitude, our new estimates are capable of reducing the uncertainty in the description of hypermagnetogenesis (without specifying a concrete model and solving the equations of the gradient-expansion formalism explicitly) down to roughly less than 1 order of magnitude. This becomes particularly apparent when comparing the explicit outcome of a specific inflationary model to three available model-independent estimates (i.e., the two estimates in Refs. [49,54] as well as our new estimate).
The paper is organized as follows. In the next section, we will review the gradient-expansion formalism that was first developed in Ref. [53] and then further refined in Ref. [92]. In particular, we will slightly adapt our notation compared to our earlier work so as to account for the fact that we are now dealing with constant values of the three parameters ξ, H, and ∆. In Sec. III, we will then represent the numerical output of our formalism after scanning over the three-dimensional parameter space of our model. We will specifically construct fit functions for our results as well as for the estimates in Refs. [49,54] and validate our approach by comparing it to the explicit results in a specific model: a simple m 2 φ 2 /2 model for three different values of the axion-vector coupling constant. The good agreement between our model-independent estimates and the explicit numerical results in this model implies that the results presented in this paper provide a good description of all scenarios of axion-driven hypermagnetogenesis that close to the origin in field space are characterized by a simple quadratic mass term. Finally, we will summarize our findings and conclude in Sec. IV. In the Appendix, we collect a number of numerical fit coefficients that enter the constructions of our fit functions. 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 by a spatially flat Friedmann-Lemaître-Robertson-Walker metric in terms of cosmic time, g µν = diag {1, −a 2 (t), −a 2 (t), −a 2 (t)}.

II. GRADIENT EXPANSION FORMALISM
Let us consider the Abelian gauge field A µ (which we will identify with the SM hypercharge gauge field shortly) axially coupled to the pseudoscalar axion inflaton field φ. The corresponding action has the form where g = det g µν is the determinant of the spacetime metric, I(φ) is the axial-coupling function, is the corresponding dual tensor; ε µνλρ is the absolutely antisymmetric Levi-Civita symbol with ε 0123 = +1. The last term in Eq. (1) corresponds to all matter fields χ a charged under the U (1) gauge group and, therefore, coupled to A µ . For the sake of generality, we do not specify the inflationary model and the axial-coupling function I(φ); we assume that the inflaton dynamics is known and consider the generation of gauge fields on this background. (This can be done only in the absence of the backreaction of generated fields which will be discussed below.) The equation of motion for the gauge field following from action (1) reads as where is the electric 4-current induced by the Schwinger effect. In addition, the Bianchi identity for the dual gauge-field strength tensor must be satisfied: To switch to a 3-vector notation, we introduce the electric E = (E 1 , E 2 , E 3 ) and magnetic B = (B 1 , B 2 , B 3 ) fields as follows: Moreover, the electric current 4-vector can be represented as where we assumed a vanishing charge density because of the quasineutrality of the plasma produced due to the Schwinger effect. Note that all 3-vectors represent physical quantities measured by a comoving observer. Then, Maxwell's equations readĖ where the dot over a symbol denotes its derivative with respect to time t, H ≡ H(t) =ȧ (t) a(t) is the Hubble parameter, and ∇ ph = ∂/∂x ph = (1/a)∂/∂x is the spatial gradient operator in physical coordinates x ph = ax. We also introduced the dimensionless parameter which controls the efficiency of gauge-field production due to the axion-vector coupling. To close the system of Maxwell's equations, we will assume that the induced current of charged particles produced by the Schwinger effect satisfies Ohm's law, where σ is the generalized conductivity, which depends only on the absolute values of electric and magnetic fields. In the case of one Dirac fermion of mass m and hypercharge Q, the Schwinger conductivity reads where g ′ is the SM U (1) Y gauge coupling constant, |E| = E 2 , |B| = B 2 , and . . . denotes the vacuum expectation value.
This expression was derived in the case of constant and collinear electric and magnetic fields in de Sitter spacetime (see, e.g., Ref. [55]). We will utilize this approximation in our analysis, too, assuming that the electric and magnetic fields change adiabatically slowly. To be more precise, we employ Eq. (13) in order to estimate the hyperelectric conductivity of the SM plasma in the presence of a hyperelectromagnetic background field, where the factor of 2H is moved to the left-hand side in order to obtain a dimensionless quantity, which will become convenient later on (see below). The expression in Eq. (14) only accounts for the production of massless SM fermions during axion inflation. In principle, the SM Higgs boson, which also interacts with the hypercharge gauge field, could be produced during axion inflation as well. However, to ensure unbroken electroweak symmetry and hence massless SM fermions, we assume that the SM Higgs field remains stabilized at the origin in field space by a large mass term throughout inflation. Such a large mass can, e.g., be induced by a nonminimal coupling to the Ricci curvature scalar R. In the numerical evaluation of Eq. (14), specifically in the evaluation of the numerical coefficient a SM , we also take into account the energy dependence of the hypercharge gauge coupling constant [93], Here, we use the full SM beta function of g ′ ; threshold effects because of the large Higgs mass during inflation are model dependent and numerically negligible. At the energy scale of the Z-boson mass, m Z ≈ 91.2 GeV, the gauge coupling equals g ′ (m Z ) ≈ 0.35. For a characteristic energy scale µ relevant for Schwinger pair production, we use In any specific model of inflation, the Hubble rate H and ξ are functions of time. However, if they change adiabatically slowly (which is consistent with the slow-roll regime during inflation), some order-of-magnitude estimates for the generated gauge fields can be obtained by considering the case of H = const and ξ = const. These estimates can be then used in any other model where the same values of H and ξ are realized. The main purpose of the present work is to derive such model-independent estimates for a wide range of constant H and ξ.
Handling vector quantities in position space makes the numerical analysis very demanding. That is why we will utilize the gradient-expansion formalism developed in Ref. [92] for the description of hypermagnetogenesis during axion inflation. It employs the vacuum expectation values of scalar products of the electric and/or magnetic field vectors with an arbitrary number of spatial curls acting on them. In this work, it will be more convenient for us to slightly adapt our notation and introduce the following set of dimensionless bilinear electromagnetic functions: Now, using Maxwell's equations (8) and (9), we obtain the system of equations for these functions, where the prime denotes the derivative with respect to dimensionless time τ = Ht and s ≡ s(τ ) = σ(τ ) 2H is the dimensionless conductivity. Terms on the right-hand sides of these equations are the boundary terms, which take into account that the number of physically relevant gauge-field modes outside the horizon continuously grows in time during inflation. Indeed, in Ref. [92], it is shown that the momentum k h of the mode crossing the horizon [defined in such a way that all modes with k < k h (t) have already experienced the tachyonic instability at time t] changes in time as Since a(t) exponentially grows during inflation, new modes cross the horizon, undergo the quantum-to-classical transition and start contributing to the classical gauge field. The explicit expressions for the boundary terms were derived in Ref. [92]. For the dimensionless quantities introduced in Eqs. (17)- (19), they take the form where G λ (ξ, s) = e πλξ r(ξ, s) ℜe W iλξ, 1 2 +s (2ir(ξ, s))W 1−iλξ, 1 2 +s (−2ir(ξ, s)) − s W −iλξ, 1 2 +s (−2ir(ξ, s)) Here, W −iλξ, 1 2 +s is the Whittaker function, and r(ξ, s) = |ξ| + ξ 2 + s + s 2 is the dimensionless physical momentum of the horizon-crossing mode.
We would also like to highlight the parameter ∆ in Eqs. (24)- (26), which was recently discussed for the first time in Ref. [92] and which modulates the magnitude of the boundary terms. It originates from the fact that gauge-field vacuum fluctuations inside the horizon are damped in the conducting medium. Indeed, as is shown in Ref. [92], the mode function deep inside the horizon (i.e., for k ≫ k h ) is represented by the damped Bunch-Davies vacuum where η = t dt ′ /a(t ′ ) is the conformal time and is the damping factor which depends on the conductivity at times t ′ ≤ t and thus makes the gauge-field evolution inherently nonlocal in time. Note that, in Eq. (31), we integrate over t ′ up to the infinite past t ′ → −∞. This implies that nonzero conductivity at some t ′ results in a suppression of all subhorizon gauge-field modes up to arbitrarily large k values, even modes that are deep inside the Hubble horizon at time t ′ and which only become tachyonically unstable at times much later than t ′ . To some extent, this is an approximation and technical simplification, as we expect the range of k values affected by the electric conductivity on subhorizon scales to be finite. Gauge-field modes with momenta much larger than the momenta of the charged fermions in the plasma should, e.g., not suffer from the damping induced by the nonvanishing conductivity. In principle, the lower integration boundary in Eq. (31) should therefore be replaced by a finite k-dependent cutoff t UV (k) ensuring that ∆ only receives contributions from times t ′ ≥ t UV (k) when the gauge-field mode with momentum k has a spatial extent larger than some UV length scale. However, at present, no exact expression for t UV (k), is known. In the following, we will therefore stick to the standard approach in the literature and simply treat all k modes on an equal footing. In particular, we do not attempt to determine the UV cutoff scale that is expected to enter Eqs. (30) and (31) at some point. This task requires further investigation and is left for future work. For the purposes of this work, it suffices to note that the conductivity σ is often a monotonically (even exponentially) increasing function of time that reaches its largest value toward the end of inflation. In realistic scenarios, the integral in Eq. (31) is therefore typically dominated by the upper integration boundary, which drastically reduces the sensitivity to the lower integration boundary. Even without a more sophisticated treatment of the lower integration boundary in Eq. (31), the dependence of the boundary terms on ∆ complicates our model-independent analysis for several reasons. For instance, if two models exhibit the same values of H and ξ, they could still have different values of ∆ because of a different prehistory. Moreover, even though the parameters H and ξ can be constant (or changing adiabatically slowly) in a real situation, the parameter ∆ always decreases in time unless the conductivity vanishes. Finally, in contrast to H and ξ, a priori it is difficult to estimate the value of the parameter ∆ without solving the full self-consistent system of equations for the inflaton, gauge fields, and charged particles. Nevertheless, we will show in the subsequent section that the magnitude of the generated field, in the case when the Schwinger effect is important, depends only weakly on ∆. Therefore, one does not necessarily need to know the exact value of ∆; a rough estimate suffices. This serves as an other reason why we defer a more detailed investigation of the lower integration boundary in Eq. (31) to future work.
The system of equations (20)- (22) for the gradient-expansion formalism is infinite by construction, since the equation for the quantity of order n contains quantities of order (n + 1). However, there exists a simple approximation, allowing us to truncate this chain at some finite order. Indeed, for large enough n, the dominant contributions to the quantities E (n) , B (n) , and G (n) correspond to the shortest modes in their spectra, i.e., modes in the vicinity of the horizon crossing mode k h . This allows us to express the higher-order quantities in terms of the lower-order ones [92], e.g., and so on. Applying these relations for some n = n max , one can truncate the infinite system of equations (20)- (22).

III. NUMERICAL ANALYSIS
The gradient-expansion formalism outlined in the previous section now allows us to determine model-independent estimates for the efficiency of gauge-field production during axion inflation. To this end, we need to fix the values of the parameters H, ξ, and ∆ and look for the stationary solution of the system of equations (20)- (22). We repeat this analysis for a large set of parameter points in the three-dimensional parameter space spanned by H, ξ, and ∆ and present our numerical results for the electric and magnetic energy densities ρ E and ρ B as well as for | E · B | in Fig. 1. In the following, we will compare our numerical results to existing estimates of ρ E , ρ B , and | E · B | in the literature [49,54] (see Secs. III A and III B) and construct semianalytical fit functions for these old estimates as well as for our own new results, i.e., for all functions shown in the first row of Fig. 1 (see Sec. III C). Rows 2 to 7 in Fig. 1 show the differences between these fit functions and the corresponding exact numerical results in row 1, which clearly demonstrates the high accuracy of our fit functions. Based on these fit functions, it is therefore now possible to reconstruct and utilize all existing estimates of gauge-field production during axion inflation, including the estimates in Refs. [49,54] as well as our new estimates based on the gradient-expansion formalism, without any further numerical analysis. Our fit functions provide all the necessary information in a compact and ready-to-use form for future applications. Finally, to validate our results, we will compare all estimates of ρ E , ρ B , and | E · B | to the exact numerical results in a specific inflationary model, namely, m 2 φ 2 /2 inflation, in Sec. III D. This will lead us to the conclusion that our new model-independent results can reproduce the exact results in a given model within an order of magnitude or so as well as that our new estimates typically improve over the existing estimates. The Hubble rate is varied between 10 8 and 10 14 GeV; the dependence on H is shown in the form of bands whose upper (lower) edges correspond to H = 10 8 GeV (H = 10 14 GeV). The gray and purple bands in the first row, respectively, correspond to the maximal and equilibrium estimates derived in Refs. [49,54] (see Secs. III A and III B). All other colorful bands in the first row correspond to our new results based on the gradient-expansion formalism. For all three types of estimates, we derive fit functions; in rows 2 to 7, we compare these fit functions to the exact numerical results in row 1.

A. Maximal estimate
The authors of Refs. [49,54] account for the effect of fermion production during axion inflation in terms of an effective gauge-field production parameter ξ eff . This is possible in the limit of perfectly parallel or antiparallel electric and magnetic fields, in which the induced current J can also be expressed in terms of the magnetic field, where the relation between sgn (E · B) and sgn (ξ) follows from our sign convention in the axion-vector coupling in Eq. (1). Making use of this relation, Ampère's law in the axion-inflationary background [see Eq. (8)] readṡ with the effective gauge-field production parameter ξ eff being defined as The same effective parameter also appears in the equation that describes the time evolution of the energy density stored in the hyperelectromagnetic field [92], which, in the limit of parallel or antiparallel electric and magnetic fields, turns intȯ Based on Eqs. (34) and (37), the authors of Refs. [49,54] now construct two estimates for the electric and magnetic fields generated during axion inflation. We will first discuss the estimate based on Eq. (37), which we will refer to as the "maximal estimate" in the following, for reasons that will become clear shortly.
Let us consider the idealized situation of a stationary de Sitter background with constant values of ξ and H. In such a background, we expect the electric and magnetic field strengths to reach a stationary attractor solution that is solely described by the values of ξ and H. 1 This solution will therefore be independent of time, which allows us to setρ em = 0 in Eq. (37). Moreover, if we momentarily interpret ρ em as the energy density stored in the entire hyperelectromagnetic field on super-and subhorizon scales and not only stored in the gauge-field modes that have already become tachyonically unstable, we can also drop the boundary term in Eq. (37), This relation represents a consistency condition, based on the requirement of energy conservation, that applies to any stationary solution for the electric and magnetic fields. In the context of our gradient-expansion formalism and introducing the shorthand notation which can be solved for |G| ≡ | E · B | /H 4 as a function of the ratio of the electric and magnetic field strengths, R, where R takes values in the range Outside this interval, no solution of the consistency condition in Eq. (39) exists. Any attractor solution withρ em = 0 satisfies the relation between |G| and R in Eq. (40), which means in particular that any attractor solution for |G| is always bounded from above by the maximal value that |G| R can obtain as a function of R, Because of the nontrivial functional dependence of |G| R on R, it is unfortunately not possible to write down a closed analytical expression for |G| max . Instead, we need to maximize |G| R , for fixed |ξ|, over the admissible range of R values numerically. Moreover, once we know the value R max that maximizes |G| R for a given |ξ|, we can use it to determine the corresponding values of E and B, We caution that E max and B max do not correspond to the maximal values of E and B that are consistent with the condition in Eq. (39); they rather correspond to the pair of |G| and R values that maximize |G| R . Together, E max , B max , and |G| max represent what we will refer to as the maximal estimate in the following. Our numerical results for the maximal estimate, based on the maximization described in Eq. (42), are shown by the gray bands in the first row in Fig. 1. The dependence on the Hubble rate in these results enters through the running of the hypercharge gauge coupling constant inside the factor a SM , which we determine in a self-consistent manner so as to satisfy the relations in Eqs. (15) and (16). As can be seen in Fig. 1, the gray bands extend down to |ξ| ∼ 4 − 5. At lower values of |ξ|, the numerical results for E max , B max , and |G| max begin to exceed the corresponding quantities in the free case without any backreaction or fermion production, which is physically not motivated. The effect of fermion production should always suppress the efficiency of gauge-field production and not enhance it. At low |ξ| values, we therefore show the free solutions for E, B, and |G| rather than the numerical output of Eqs. (42) and (43).
In the following, we will now discuss simple and novel fit functions that manage to describe our exact numerical results for E max , B max , and |G| max to high accuracy. The starting point of our construction is going to be the value of |G| R evaluated at R = 1, which corresponds to equal amounts of energy in the electric and magnetic fields, It then turns out that our numerical data are well described by the fit functions In addition to |ξ|, these three functions also depend on H through the energy dependence of the hypercharge gauge coupling constant. To the first approximation, this effect may be neglected, and g ′ may be simply set to a characteristic value around g ′ ∼ 0.4. However, to obtain even more accurate estimates, we are actually able to determine the selfconsistent value of g ′ analytically. To this end, we need to express the renormalization scale µ in Eq. (16) in terms of the semianalytical expressions for E max and B max and factor out the dependence on the gauge coupling constant g ′ , where the rescaled quantityμ max ,μ is defined in terms of rescaled versions of E max and B max that no longer depend on g ′ , because they are no longer proportional to a SM but instead proportional tō Making use of these definitions, the hypercharge gauge coupling constant in Eq. (15) can be written as The self-consistent solution for g ′ at the one-loop level thus follows from solving this relation for g ′ , where the subscript "max" again indicates that this expression for g ′ corresponds to the self-consistent solution in the case of the maximal estimate. W 0 denotes the principal branch of the Lambert W function or product logarithm, and the argument of the Lambert W function, z max , and the numerical coefficient b are given as follows: which makes the dependence on the value of the Hubble rate H explicit. If we now use the result in Eq. (54) in order to evaluate a SM in Eqs. (44), (45), (46), and (47), we obtain semianalytical fit functions for E max , B max , and |G| max that self-consistently account for the running of the gauge coupling constant. In row 2 of Fig. 1, we compare these fit functions to the exact numerical result indicated by the three gray bands in row 1 and find excellent agreement. Finally, before we move on to the next estimate, which we will refer to as the "equilibrium estimate", we mention that g ′ max in Eq. (54) is also well approximated by the following, much simpler expression: g ′ max ≃ 0.4162 + 0.0068 log 10 H 12 + 0.0030 ln |ξ| 10 , Using this expression in Eqs. (44), (45), (46), and (47) results in nearly equally accurate fit functions.

B. Equilibrium estimate
Next, we turn to the second estimate proposed by the authors of Refs. [49,54], the equilibrium estimate, which is based on the modified version of Ampère's law in Eq. (34). The maximal estimate constructed in the previous section represents an upper bound on all attractor solutions in a stationary de Sitter background with constant ξ and H. The equilibrium estimate, by contrast, aims at actually constructing an explicit attractor solution. The basic idea is that, once the electric and magnetic fields have reached the attractor solution, Ampère's law in Eq. (34) will have the same solution as in the free case without any backreaction or fermion production, the only difference being that the parameter ξ in this solution needs to be replaced by ξ eff . If this is the case, the equations of motion for the gauge-field modes in Fourier space will also be solved by the usual Whittaker functions, with ξ → ξ eff , such that Thus, an underlying assumption of this approach is that, after a sufficiently long time, the system reaches an attractor solution that resembles the free solution (modulo the replacement ξ → ξ eff ) not only at the level of the integrated quantities E, B, and |G| but also at the level of the Fourier mode spectrum. For large values of the effective gauge-field production parameter, |ξ eff | ≫ 3, the quantities in Eqs. (58), (60), and (59), for n = 0, can in particular be written as where the numerical coefficients are roughly given by C E ≃ 2.6×10 −4 , C B ≃ 3.0×10 −4 , and C G ≃ √ C E C B ≃ 2.8×10 −4 . In the following, we will explicitly work with the relation C G ≃ √ C E C B , which is valid in the limit of parallel or antiparallel electric and magnetic fields, i.e., the limit that has been used in the derivation of the induced current.
Equation (61) implicitly defines the equilibrium estimate for E, B, and |G|. The evaluation of this estimate is, however, complicated by the fact that the effective parameter ξ eff in Eq. (61) also depends on the electric and magnetic field strengths. Again, we are thus not able to write down a closed analytical solution but have to resort to a numerical approach. To this end, we first note that Eq. (61) results in a simple relation between R and |ξ eff |, Using the definition of the effective gauge-field production parameter |ξ eff | in Eq. (35), we are therefore able to write which can be solved for E eq as a function of |ξ eff | = C B /C E R eq . Together with the relations [see Eq. (61)] we thus find By comparing these expressions with the expressions in Eq. (61), we obtain a single consistency condition for |ξ eff |, Therefore, to evaluate the equilibrium estimate, we need to numerically solve this condition for |ξ eff |, for fixed values of ξ and H, while making sure that the relations in Eqs. (15) and (16) are self-consistently satisfied. The numerical result for |ξ eff | that we obtain in this way can then be used in Eq. (61). This procedure defines our numerical results for E eq , B eq , and |G| eq , which are shown in the form of purple bands in the first row of Fig. 1.
Similarly as in the case of the maximal estimate, we are again able to describe our numerical results in terms of fit functions. This time, the entire relevant information can be encoded in the fit function for |ξ eff |, |ξ eff | eq ≃ a eq ln (|ξ| + b eq ) + c eq (69) with coefficients a eq ≃ 0.3679 − 0.0004 log 10 H 12 , (70) b eq ≃ −3.3668 + 0.0099 log 10 H 12 , (71) c eq ≃ 3.7012 − 0.0152 log 10 H 12 , where H 12 denotes again the Hubble rate in units of 10 12 GeV [see Eq. (57)]. In the third row of Fig. 1, we compare the approximate results for E eq , B eq , and |G| eq based on this fit function to the exact numerical results in the first row. Again, we find excellent agreement, in the regime where the expressions in Eq. (61) are valid, i.e., for |ξ eff | ≫ 3. By construction, the fit function in Eq. (69) already takes into account the running of the gauge coupling constant. It is therefore not necessary to work out an independent fit function for g ′ . This differs from the case of the maximal estimate, where we were able to solve the maximization condition in Eq. (42) without specifying the coefficient a SM . For completeness, we, however, note that the same strategy that eventually led to g ′ max in Eq. (54) can be applied in order to determine the self-consistent solution for the gauge coupling constant in the case of the equilibrium estimate, where the rescaled renormalization scaleμ eq is now given bȳ Equation (73), together with |ξ eff | eq in Eq. (69), results in excellent agreement (at the level of 10 −3 %) with our numerical results for g ′ eq . Alternatively, we can simply solve the consistency condition in Eq. (68) for g ′ , This is an exact expression for g ′ eq , which, however, is more sensitive to deviations of the fit function in Eq. (69) from the exact numerical result. The combination of Eqs. (69) and (77) still results in a good approximation of the exact numerical result for g ′ eq , with the numerical deviations mostly remaining below the percent level. Finally, we are also able to describe the exact numerical result for g ′ eq with the simple fit function [see also Eq. (57)] g ′ eq ≃ 0.4131 + 0.0067 log 10 H 12 + 0.0025 ln |ξ| 10 , which reproduces the exact result up to deviations at the level of around 0.1%.

C. Semianalytical fit functions
In the previous two sections, we constructed novel fit functions for the two estimates of the efficiency of gauge-field production that had originally been proposed in Refs. [49,54]. Now, we turn to our own numerical results based on the gradient-expansion formalism, i.e., the colorful bands for ∆ = 10 −6 , 10 −4 , 10 −2 , and 1 in the first row of Fig. 1. To fit our numerical results, we make an ansatz for X = E, B, and |G| of the form where X ∆=1 is supposed to describe our data for ∆ = 1 and the function S X accounts for the suppression of the quantity X if the parameter ∆ is smaller than unity. In fact, it turns out convenient to write S X as a power of ∆, We furthermore approximate X ∆=1 by two different expressions at small and large values of |ξ|, where the threshold value |ξ| X at which we switch from one expression to the other is chosen as for E, B, and |G|, respectively. At small |ξ|, we relate our results to the free solution without any backreaction or fermion production, while at large |ξ|, we express our results in relation to the maximal estimate defined in Sec. III A, Here, X free (X = E, B, |G|) is given by the three integral expressions in Eqs. (58), (60), and (59) for n = 0 and after undoing the replacement |ξ| → |ξ eff |; X max (X = E, B, |G|) corresponds to our three fit functions for the maximal estimate in Eqs. (45), (46), and (47) in combination with our result for g ′ max in Eq. (54).
The nontrivial information contained in our numerical results is thus captured by the three functions P X , T X , and U X , for each of the three quantities E, B, and |G|, in Eqs. (80) and (83). For each of these functions, we make a particular (purely phenomenological) ansatz that turns out to describe our numerical data to very good accuracy, For each X, we hence need to determine 11 fit coefficients: Our results for these, in total, 33 coefficients, which depend on the logarithm of H as well as partially on the logarithm of ∆, are listed in the Appendix. In rows 4 to 7 of Fig. 1, we compare the resulting fit functions for E, B, and |G| with the exact numerical results shown in row 1. As before, we find excellent agreement. On a logarithmic scale, our fit functions begin to deviate from the exact numerical results typically only in the third significant digit. In summary, we therefore conclude that the semianalytical fit functions constructed in Secs. III A and III B as well as in the present section are capable of reproducing all relevant numerical results with very good accuracy. In future work, it will no longer be necessary to repeat the numerical analysis that originally led to these fit functions.

D. Validation in a specific model
In this section, we test the accuracy of our model-independent approach by comparing it to the exact results for a specific inflationary model. For this purpose, we consider the simple model with a quadratic inflaton potential, This quadratic dependence is universal for a wide class of inflaton potentials close to their minima; therefore, validating our formalism in this model will allow us to make more general conclusions. Indeed, since the most efficient gauge-field production occurs close to the end of inflation, the behavior of the inflaton potential far from its minimum is not of great importance for hypermagnetogenesis. In our numerical analysis, we will set m = 6 × 10 −6 M P , for concreteness, which is the same value that we used in Ref. [92]. However, because of the one-to-one relation between m and the Hubble rate [see Eq. (91)], we stress that other values of the inflaton mass will only lead to logarithmic corrections to our results. We therefore expect that the following analysis applies, in fact, to a large range of m values. We take the axial-coupling function in a linear form, with a dimensionless coupling parameter β. To obtain the exact numerical results for the generated gauge fields in this model for a given value of β, we apply the gradient-expansion formalism developed by us in Ref. [92]. It is worth noting that the gradient-expansion formalism is itself an approximate method; however, it was shown in Ref. [92] that its error compared to the exact mode-by-mode solution can be made as small as 1% − 2% during the whole stage of inflation. As we will see, such an accuracy is much better than that of the model-independent approach; therefore, we can use the gradient-expansion result as a reference solution.
In practice, we proceed as follows. For a given value of β, we first apply the gradient-expansion formalism of Ref. [92] and obtain the hyperelectric and hypermagnetic energy densities, ρ E and ρ B , as well as the Chern-Pontryagin density | E · B | as functions of the number of e-folds N e until the end of inflation. In addition, we also compute the corresponding values of the gauge-field production parameter ξ, the Hubble rate H, and the damping factor ∆. Knowing the latter three parameters at the same moment of time then allows us to apply our model-independent approach and find the predictions for the generated gauge fields at this moment. Comparing these predictions with the reference solutions, we are able to draw conclusions concerning the accuracy of the model-independent approach.
One may argue that such a usage of the model-independent approach has no advantage compared to the full gradient-expansion formalism because we have to launch the latter method in any case in order to obtain the set of parameters (ξ, H, ∆) for the former method. However, this is done only for the purpose of comparing the two methods. Normally, to arrive at model-independent predictions, it suffices to determine the parameters (ξ, H, ∆) from some other, much simpler consideration. For instance, the values of ξ and H can be estimated by considering the inflaton dynamics neglecting the backreaction of the generated gauge fields (it is shown in Ref. [92] that this is a reasonable approximation for a wide range of parameters in the presence of the Schwinger effect). However, the parameter ∆ cannot be determined by a simple method. Therefore, it is interesting to check whether one can simply use the fixed value ∆ = 1 in the model-independent approach. This value is well motivated by the following arguments. For small ξ, when the dependence of the generated gauge fields on ∆ is strong (see Fig. 1), the gauge fields are rather weak; therefore, the Schwinger conductivity is small compared to the Hubble parameter, and ∆ is indeed close to unity (unless the system had a nontrivial prehistory including a period with very high conductivity). In the opposite case of large ξ, the parameter ∆ can be much less than unity; however, the generated fields exhibit only a very weak dependence on ∆, and their values do not differ much from those with ∆ = 1; see Fig. 1. Therefore, to check the validity of this approximation, we will in addition also apply the model-independent approach to the same values of (ξ, H) as before in combination with a fixed value of ∆ = 1.
Finally, for given ξ and H, we also compute the maximal and equilibrium estimates discussed in the previous sections. Thus, for a fixed value of β, we obtain exact numerical results for the generated gauge fields and four different approximate results. We compare them in Figs Let us now analyze and comment on the obtained numerical results. First of all, we mention that the modelindependent approach presented in this paper allows us to estimate the magnitude of the generated gauge field within at most 1 order of magnitude in comparison to the exact numerical results. The maximal deviation occurs very close to the end of inflation (during the last one to two e-folds), where the change in the parameters ξ, H, and ∆ cannot be considered to be adiabatically slow (see Fig. 5 and its discussion below). This deviation can be easily understood since the model-independent approach assumes that the above-mentioned parameters are constant. Second, the model-independent approach with a fixed value of ∆ = 1 typically results in a slightly worse agreement (although it is still comparable with the model-independent approach with the exact value of ∆). However, for larger values of β, at the very end of inflation, its predictions are accidentally even in better accordance with the exact result (compare the red and blue curves in Figs. 3 and 4). Third, the model-independent approach typically gives more accurate results than the equilibrium or maximal estimates. Only for the hyperelectric energy density, it sometimes happens that the latter estimates have a comparable accuracy with the model-independent approach.
Finally, let us discuss the reasons for the deviation of the approximate methods considered above from the exact numerical solution. As we already pointed out, all approximate methods rely on the fact that their input parameters ξ, H, and ∆ are constant. Such an approximation would be reasonable if these parameters were changing adiabatically slowly, i.e., their change during the Hubble time was much smaller than the absolute value of the parameter. This condition can be characterized by the adiabaticity parameter ǫ P defined for any P = {ξ, H, ∆} in the following way: In particular, the parameter ǫ H is the well-known slow-roll parameter that controls when inflation terminates. Some preliminary estimates for the parameters ǫ H and ǫ ξ can be obtained from the slow-roll analysis. For the inflationary model with the scalar potential (87) in the absence of any backreaction from the generated gauge fields, it is possible to analytically find the dependence of the inflaton field on the number of e-folds before the end of inflation, from which one immediately obtains the slow-roll expressions for ǫ H , ǫ ξ , given the coupling function I(φ) = βφ/M P , |ξ| ≃ βM P 2 Concerning the parameter ǫ ∆ , it follows from Eq. (31) that Although we cannot derive any slow-roll estimates for this parameter, it is clear that it increases when the generated gauge field becomes stronger. Figure 5 shows the parameters ξ, H, and ∆ during the last 15 e-folds of inflation (row 1) and the corresponding parameters ǫ ξ , ǫ H , and ǫ ∆ (row 2). In panel a1 of Fig. 5, we show the parameter ξ for β = 10 (red line), β = 15 (green line), and β = 20 (blue line). The solid lines correspond to the exact numerical result, while the dashed lines of the same colors give the slow-roll estimates according to Eq. (92). Since ξ ∝ β, the parameter ǫ ξ does not depend on β. The same also holds for H and ǫ H , if there is no backreaction from the generated gauge fields on the inflaton dynamics (this is indeed the case in our model). Therefore, in panels a2, b1, and b2, we show exact results for ǫ ξ , H, and ǫ H , respectively, by the blue solid lines and the corresponding slow-roll estimates by the red dashed lines. It is worth noting that in the slow-roll approximation ǫ ξ = ǫ H . As expected, ǫ H and ǫ ξ tend to unity when inflation ends. Panel c1 of Fig. 5 shows the dependence of the parameter ∆ on the number of e-folds for β = 10 (red dotted line), β = 15 (green dashed line), and β = 20 (blue solid line). Far from the end of inflation, this parameter is very close to unity because the gauge fields are weak and, consequently, the Schwinger conductivity is small (see panel c2, where the corresponding values of ǫ ∆ = σ/H are shown). However, when inflation ends, it becomes exponentially small. The corresponding adiabaticity parameter exceeds unity (this happens earlier for larger β), and the adiabatic approximation fails. The behavior illustrated in Fig. 5 thus explains why the approximate model-independent results begin to differ from the exact solution close to the end of inflation.

IV. CONCLUSION
The description of gauge-field production during axion inflation is relevant for a variety of phenomena, ranging from primordial magnetic fields over baryogenesis to primordial gravitational waves and black holes. In the presence of the Schwinger effect, this process becomes highly nonlinear and typical approaches dealing with separate gauge-field modes in Fourier space become inapplicable. To overcome this difficulty, we proposed in Ref. [92] a novel gradient-expansion formalism that operates with a set of bilinear gauge-field quantities in coordinate space and that allows us to build a complete and self-consistent system of equations for studying hypermagnetogenesis during axion inflation in the presence of nonlinear effects, such as the backreaction of the generated fields on the inflaton field and the Schwinger effect. Although this system of equations can be explicitly solved in a given inflationary model, such an analysis is quite complicated and requires some computational efforts. However, for many practical purposes, it would be desirable and convenient to have some model-independent predictions allowing one to estimate the magnitude of the generated gauge fields without carrying out a complicated numerical analysis. In the free case involving no nonlinear effects at all, such model-independent estimates are well known in the literature: the electric energy density, magnetic energy density, and Chern-Pontryagin density, respectively, scale like H 4 e 2πξ /ξ 3 , H 4 e 2πξ /ξ 5 , and H 4 e 2πξ /ξ 4 for large ξ, if all nonlinear effects can be neglected. The aim of the present work was to generalize these model-independent estimates to the full SM case in the presence of Schwinger pair production of all SM fermions.
The fact that gauge-field production during axion inflation is controlled by ξ and H has been known for a long time. We recall once more that the parameter ξ characterizes the velocity of the inflaton field and H denotes the Hubble expansion rate. Both parameters can be easily determined from the standard slow-roll analysis in a concrete inflationary model. In the presence of Schwinger pair production, there appears one additional parameter ∆, which describes the damping of the vacuum gauge-field fluctuations due to the finite conductivity of the Universe on subhorizon scales [92]. This parameter depends on the prehistory of the system [see Eq. (31)] and thus makes its evolution nonlocal in time, which complicates the model-independent analysis of gauge-field production to some degree.
The main assumption behind our analysis was that gauge-field production at some moment of time during axion inflation is determined by the momentary values of the parameters ξ, H, and ∆. This assumption is based on the fact that the exponentially fast cosmic expansion quickly dilutes the gauge fields generated at earlier times and thus the dominating contribution to the magnitude of the gauge fields at some moment of time originates from the modes that crossed the horizon just shortly before. A convenient criterion for the validity of this assumption is the adiabaticity of evolution of all three parameters: their change during the Hubble time must be much less than the absolute value of the parameter itself. For the parameters ξ and H, this adiabatically slow evolution indeed holds during slow-roll inflation. As for the parameter ∆, it changes slowly only in the weak-field regime when the corresponding Schwinger conductivity is much less than the Hubble parameter, σ ≪ H; see Fig. 5 for an explicit example.
To derive our new model-independent results, we employed the gradient-expansion formalism presented in Ref. [92], for which we used constant values of ξ, H, and ∆ as input parameters. The main difference compared to our earlier model-dependent analysis in Ref. [92] is that, in the full system of equations, these parameters are computed selfconsistently by considering the evolution of the inflaton field and scale factor together with the gauge fields. However, as we have shown in this paper, the approximation of constant ξ, H, and ∆ indeed works well when their respective time variation in a concrete model remains adiabatically slow. In our numerical analysis, we fixed the parameters ξ, H, and ∆ and looked for the stationary solution of our system of equations, which provided us with the prediction for the generated gauge field. Then, we scanned over wide ranges of parameter values that are sufficient for most physical applications, namely, 1 < ξ < 15, 10 8 GeV < H < 10 14 GeV, and 10 −6 < ∆ < 1. Our numerical results for the hyperelectric and hypermagnetic energy densities, ρ E and ρ B , as well as the Chern-Pontryagin density | E · B | are summarized in Fig. 1, which is the main result of our study. There, we also compare our predictions to estimates that had previously been obtained in the literature, in particular, the maximal and equilibrium estimates derived in Refs. [49,54]. They were derived as upper and lower constraints on | E · B | without taking into account the damping of vacuum fluctuations by the parameter ∆, i.e., for ∆ = 1. Our predictions for | E · B | in the case ∆ = 1 lie between the two above-mentioned estimates, thus being in good accordance with them.
For small values of ∆, the results change dramatically only for small values of ξ, when the gauge fields are weak. Indeed, in such a case, each new mode crossing the horizon makes an important contribution to the total energy density. Damping of these new modes thus significantly changes the resulting gauge field. On the contrary, for the case of strong gauge fields, the contributions of new modes crossing the horizon are small compared to those that are already outside the horizon and are enhanced due to the axion coupling. Therefore, damping of these new modes by the ∆ parameter makes a small impact on the generated field. These features are clearly seen from Fig. 1.
By construction, our new model-independent estimates do not account for possible backreaction effects, as they are based on the assumption of adiabatically slowly varying values of ξ and H. This, however, does not limit the range of applicability of our results as severely as one may naively think. In Ref. [92], we showed that Schwinger pair production often suppresses backreaction effects in scenarios in which it would otherwise be relevant. In the presence of the Schwinger effect, backreaction therefore only occurs in extreme regions of parameter space. As long as it is negligible, our results are applicable and can be considered as the straightforward generalization of the corresponding expressions in the free case (i.e., the expressions proportional to H 4 e 2πξ /ξ 3 , H 4 e 2πξ /ξ 5 , and H 4 e 2πξ /ξ 4 ). To make our results more accessible and easier to use, we provide semianalytical fit functions that describe our entire numerical data with very high accuracy across the full range of parameter values that we considered in this paper.
To validate our model-independent results, we considered a concrete inflationary model with potential V (φ) = m 2 φ 2 /2. Although such a potential is already discarded by cosmic microwave background observations, it is still worth considering because many other inflaton potentials can be approximated by m 2 φ 2 /2 close to their minima, and this region appears to be the most important for the generation of gauge fields, which occurs during the last few e-folds 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). We implemented the inflationary model with potential V (φ) = m 2 φ 2 /2 in the full gradient-expansion formalism. Then, we used the exact values of the parameters ξ, H, and ∆ at a sequence of moments of time close to the end of inflation and launched the model-independent approach for these values of parameters (again treating them as constants). By comparing these approximate results with the results of the self-consistent gradient-expansion formalism, we conclude that the model-independent results indeed can be used to estimate the magnitude of the produced gauge fields with an error typically less than 1 order of magnitude. This is a significant improvement over previous estimates, which spanned several orders of magnitude. The main accomplishment of the present paper is therefore a significant reduction in the theoretical error in the description of hypermagnetogenesis during axion inflation. The largest error is reached close to the end of inflation, where the adiabaticity conditions for ξ, H, and ∆ break down, while far from the end of inflation, the accuracy of the modelindependent result is much better. Moreover, we show that one can even use the fixed value ∆ = 1 during the whole stage of inflation and the accuracy of our model-independent predictions remains of the same order, which facilitates the usage of our results.
We stress again that our model-independent results are particularly well suited to estimate the efficiency of gaugefield production during inflation. Toward the end of inflation, where the slow-roll approximation breaks down, only estimates within 1 order of magnitude are possible (which, however, still improves on earlier estimates). In future work, we will therefore turn to a dedicated and more precise description of the initial conditions for reheating after the end of inflation. This analysis will then allow one to connect the dynamics of hypermagnetogenesis during axion inflation to the subsequent evolution during reheating and the radiation-dominated stage after axion inflation.