Numerical study of the Schwinger effect in axion inflation

Previous studies demonstrate that the inflaton, when coupled to the hypercharge Chern-Simons density, can source an explosive production of helical hypermagnetic fields. Then, in the absence of fermion production, those fields have the capability of preheating the Universe after inflation and triggering a successful baryogenesis mechanism at the electroweak phase transition. In the presence of fermion production however, we expect a strong damping of the gauge fields production from the fermion backreaction, a phenomenon called Schwinger effect, thus jeopardizing their original capabilities. Using numerical methods we study the backreaction on the generated gauge fields and revisit the processes of gauge preheating and baryogenesis in the presence of the Schwinger effect. We have found that gauge preheating is very unlikely, while still having a sizable window in the parameter space to achieve the baryon asymmetry of the Universe at the electroweak phase transition.


Introduction
Cosmological inflation [1][2][3] is nowadays a well-established paradigm to solve the classical (flatness, horizon, . . . ) problems of the Standard Cosmological Model, and to generate the primordial density perturbations giving rise to the present Universe structure. The achievements of cosmological inflation usually require the presence of one (or several) scalar field -the inflaton-giving rise to physics beyond the Standard Model (SM) of Particle Physics (BSM). In this way, along with the classical problems of the SM (hierarchy problem, baryogenesis, strong CP-problem, dark matter,...), cosmological inflation provides yet another motivation for BSM physics.
Although the existence of a period of cosmological inflation is pretty well stablished by observational cosmological data [4], there is no consensus on a detailed model. An interesting candidate for the inflaton is a pseudoscalar field ϕ, denoted in this paper as axionlike particle 1 , which can then couple to the Chern-Simons density F µνF µν of a U (1) gauge field. In this case, and depending on the size of the coupling of the inflaton to the Chern-Simons term, there can be an explosive production of helical gauge fields at the end of inflation [7][8][9]. This exponential production can dominate the energy density of the Universe during the coherent oscillations of the inflaton around its minimum, a phenomenon dubbed as gauge preheating [10][11][12], and lead to a rapid production of inhomogeneities sourcing a significant gravitational wave background, leading to strong constraints on the inflaton Chern-Simons coupling from the Planck (and future CMB-S4) limits on the net energy density in gravitational waves [5,6].
When the gauge field is identified with the SM hypercharge, Y µ , with strength Y µν , the inflaton coupling to the Chern-Simons density Y µνỸ µν gives rise to the production of helical hypermagnetic fields which can then survive until the electroweak phase transition (crossover), and trigger the baryon asymmetry of the Universe (BAU) [13][14][15][16][17][18]. However, in the presence of strong gauge fields, light fermions charged under the gauge group are produced by the backreaction of gauge fields that source the fermion equations of motion (EoM) [19,20]. The corresponding currents can then, in turn, backreact on the produced gauge fields, a phenomenon called Schwinger effect (see e.g. Ref. [21]). The backreaction of fermion currents on the produced gauge fields acts as a damping force in the explosive production of helical gauge fields, and many of the conclusions from the gauge field production should be revised in the presence of the Schwinger effect 2 , in particular those concerning the preheating capabilities and the baryogenesis mechanism.
In this paper, we will study the effect of the Schwinger particle production on the helical hypermagnetic fields produced at the end of inflation, and in particular its influence on the gauge preheating efficiency and baryogenesis capability. In order to consider the backreaction of the produced gauge fields on the inflationary equations of motion and that of the Schwinger effect on the gauge field production, we will use numerical methods, in particular, the fourth order Runge-Kutta (RK4) algorithm. Our numerical results are validated as they overlap with some recent semianalytical methods and the gradient expansion formalism of Refs. [22][23][24][25]. Our general finding is that the gauge field production is much less explosive than in the absence of the Schwinger effect, which will jeopardize the conclusions concerning the possibility of gauge preheating, although they leave an open window for baryogenesis.
The contents of this paper are as follows. In Sec. 2 we present the general lines of the model and the methods we will consider, including the relevant equations of motion in momentum space and the observable quantities we will compute. In Sec. 3 we will present numerical results for the gauge sector assuming the slow roll conditions in the inflaton equations of motion. In order to check the validity of our approach, numerical results will be compared with some estimates from the literature, namely the backreactionless solution, the Schwinger equilibrium and maximal estimates as well as the gradient expansion formalism where dynamical results are obtained analytically and numerically in configuration space. Some details about the numerical methods will be explained in App. A. In Sec. 4 we will perform the full numerical calculation for two kinds of models that predict cosmological observables in agreement with the observed values: the α-attractor models and the quartic hilltop models. In all cases the gauge preheating efficiency does not seem good enough to ensure complete reheating, which has to be completed by other perturbative or nonperturbative mechanisms. Moreover we have reanalyzed the baryogenesis predictions in the presence of the Schwinger effect and found a sizable window where the BAU is correctly predicted. Again, some details about the numerical methods we used are described in App. B. Finally, our conclusions are presented in Sec. 5.

The model
The model action is given by where ϕ is the pseudoscalar inflaton, V the inflaton potential, and f ϕ provides the inverse coupling of the inflaton to the Chern-Simons term. Y µν is the field strength of the hypercharge gauge field Y µ and Y µν = 1 2 ϵ µνρσ Y ρσ its dual tensor. We also have included the interaction of fermionic currents, corresponding to hypercharge Q Y fermions, with the hypercharge fields (encoded in the covariant derivative D µ ≡ ∂ µ − g ′ Q Y A µ ). All gauge field quantities are U (1) hypercharge fields, i.e. A Y , E Y , B Y , etc. To make the notation lighter, we drop the index Y as there will be no ordinary electromagnetic fields in this work.

Equations of motion
Variation of the action with respect to ϕ and the hypercharge gauge field A µ = (A 0 , A) leads to the gauge equations of motion in the radiation gauge, A 0 = 0 and ∇ · A = 0, We assume that initially the Universe does not contain any asymmetry of charged particles and that these ones are produced only later in particle-antiparticle pairs. Therefore, we set the charge density to zero, ρ c = 0. Finally, the current J is given by the Ohm's law where σ is the generalized conductivity, which will be defined later.
As it can be seen from the above system, we use cosmic time t for the inflaton dependence and the conformal time τ , defined by g µν = a 2 (τ ) η µν , for the gauge field dependence. We denote the derivative with respect to conformal time τ with a prime and the derivative with respect to the cosmic time t with a dot, e.g. a ′ = da/dτ andȧ = da/dt. The Hubble parameter is defined as H =ȧ(t)/a(t) where a is the scale factor. We assume a homogeneous inflaton with only zero mode, ϕ(t, x) = ϕ(t).
We now quantize the gauge field A in momentum space where λ is the photon polarization, and a λ (k) (a † λ (k)) are annihilation (creation) operators that fulfill the canonical commutation relations The polarization vectors ϵ λ (k) satisfy the conditions 3 where k ≡ |k|. Therefore, the equation of motion for the gauge modes yields In some special cases (σ = 0 and slow-roll inflation), this equation can be solved analytically, and we will do it in Sec. 3.1. In the general case, we will solve it using numerical methods.

Observable quantities
Once we obtain a solution to the modes A λ , we can compute the (hyper)electromagnetic energy densities as The upper integration limit comes because subhorizon modes have an oscillatory behavior and should be regarded as quantum fluctuations. Therefore, such modes do not contribute to the above classical observables and are excluded from the integration. More details and precise value of k c will be given in Sec. 2.3. For the lower integration limit k min , see Eq. (A.5).
In this work, we will also make use of the (hyper)magnetic helicity and its derivative, defined as In the case of one Dirac fermion with mass m and hypercharge Q Y , the conductivity can be written as 4 [19] where g ′ ≃ 0.4 is computed at the characteristic scale µ ≃ (⟨E⟩ 2 + ⟨B⟩ 2 ) 1/4 where the Schwinger effect takes place [24]. This estimation is valid in the case of collinear electric and magnetic fields, an assumption that we have numerically checked by verifying that where θ is the spatial angle between E and B. Moreover, the massless hypercharged fermions that are continuously produced during inflation have an energy density given by Notice that the observable quantities ρ E , ρ B , ρ ψ , H and G are physical 5 , while the fields A, E and B as well as the conductivity σ and current J are comoving.
Concerning the Higgs vacuum expectation value, there are two possibilities during the inflationary period: i) The first possibility, which we will consider throughout this paper, is that ⟨h⟩ = 0, and so the electroweak symmetry is unbroken during the inflationary period. In order to ensure unbroken electroweak symmetry and hence massless SM fermions, which all contribute to the conductivity (2.10), we assume that the SM Higgs field h remains stabilized at the origin in field space by a large mass term throughout the inflationary period. Such a large mass can, e.g., be induced by a nonminimal coupling to the Ricci curvature scalar as L = 1 2 ξh 2 R with ξ > 3/16 (see e.g. Ref. [26]). Hence, we get ii) The second possibility is that the electroweak symmetry is broken during the inflationary period.
In this case after ∆N e-folds of inflation, there is a Gaussian distribution of values of the Higgs field with zero mean and variance ⟨h 2 [26] 6 . In this case, the electroweak symmetry is broken and the hypercharge field strength Y µν in Eq. (2.1) is projected onto the electromagnetic field strength F µν with a coupling to the inflaton given by f ϕ / cos 2 θ W where θ W 4 As the conductivity σ relates J and E in (2.3), it is a comoving quantity, i.e. it scales with the Universe expansion. Our definition differs from the one in [23,24] where the authors used a physical conductivity that we will denoteσ in this paper, their relation being σ = aσ. 5 They relate to the comoving ones ρ c E , ρ c B , H c , and G c by the relations ρ c B,E = a 4 ρ B,E , H c = a 3 H, G c = a 4 G. 6 The SM Higgs potential is still unstable at a value of the Higgs field h = h I ≃ 10 11 GeV and the condition for P (h I , ∆N ) < e −3∆N (so that it is unlikely to find the Higgs away from its EW vacuum in any of the e 3∆N causally disconnected regions formed during inflation) implies H E < 2/3 πh I /∆N , a condition that is not fulfilled by any of the models of inflation we have considered. Therefore this possibility would require stabilization of the Higgs potential by some new physics.
is the electroweak angle. Now the conductivity for the hypermagnetic field in Eq. (2.10) should be replaced by a similar expression for the magnetic field, with the replacement |g ′ Q Y | → |eQ|, where e = gg ′ / g 2 + g ′2 and Q is the fermion electric charge. The condition for a fermion f to contribute to the magnetic conductivity πm 2 f < √ 2ρ E |eQ f | translates into the condition, for the fermion Yukawa coupling, and we have computed all couplings at the characteristic scale µ ≃ (⟨E⟩ 2 + ⟨B⟩ 2 ) 1/4 where the Schwinger effect takes place. If the three generations of fermions satisfy the above condition then the conductivity for the magnetic field is given by Eq. (2.13) with the replacement 41g ′3 72π 2 → e 3 π 2 . We have checked that, in this case, the results for f ϕ ≲ 0.2 M pl are consistent with all three generation fermions contributing to the magnetic conductivity. For f ϕ ≳ 0.2 M pl only the top quark does not contribute. Given that 41g ′3 /72 ≃ 0.37 while e 3 ≃ 0.36, at the scales where the Schwinger effect takes place, we have found that the results in this second case are qualitatively similar to those for the previous case, which will be worked out in detail in this paper.
Considering then the case i) above, the conductivity (2.13) yields a nontrivial integro-differential system as the damping term grows with the magnetic energy and hence backreacts on the amount of produced electric/magnetic fields. We aim to solve this setup of the Schwinger effect numerically. In the next sections we will consider specific cases where this system can be further simplified.

The gauge vacuum
At very early times, when |aφ| ≪ kf ϕ , the modes are in their Bunch-Davies (BD) vacuum, hence Initially, we can consider all the modes in the BD vacuum (which would be possible by initializing the numerical simulation such that a 0 ≪ k 0 /H 0 ). In that case, since |A + | = |A − |, the fields E and B are plane waves perpendicular to each other, as G = 0 in (2.11) yields cos θ = 0. Therefore, there is no Schwinger effect and σ = 0. It has recently been shown that in the presence of the conductivity σ, the BD vacuum amplitude of the modes that are still in the vacuum get damped by the ones that left it [23]. Indeed, consider we are at a time a * where modes k > k * are still in the BD vacuum, while modes k < k * were amplified by both tachyonic and parametric instabilities from Eq. (2.7). Then, the equation of motion for modes such that |a * φ (τ * )| ≪ kf ϕ does not reduce to a plane wave in the presence of a non-zero σ, but instead to A ′′ λ + σA ′ λ + k 2 A λ = 0, and Eq. (2.15) is not a solution anymore. To derive the generalized BD vacuum, we write the gauge equation of motion (2.7) in cosmic time: where we used the identity a −2 A ′′ λ =Ä λ + HȦ λ , and perform the transformation (2.17) We recall that we have definedσ = σ/a as the physical conductivity in footnote 4. The above equation hence becomes where we used the fact that∆(t) = −σ(t)∆(t). A mode crosses the horizon when the expression in the square brackets vanishes for the first time at least for one polarization, at k = k c . The modes in the vacuum are then characterized by k ≫ k c . This yields the momentum of the mode that crosses the horizon at time t, namely the cutoff of the integrals: Deep inside the horizon, when the first term in square brackets of (2.18) dominates, the solution must satisfy the BD condition (2.15). As we have seen, in the presence of finite conductivity, this equation does not fully describe the gauge-field mode function inside the horizon, as the damped BD condition includes an exponential damping factor The bottom line of this section is that the modes still in their BD vacuum see their amplitudes damped because of the effect of the modes that left their vacuum earlier and participate in the equations of motion (2.2a) and (2.7). The parameter ∆ was first introduced in the context of the gradient expansion formalism in Ref. [23], where it was dynamically solved, while in Ref. [24] it was also considered as a free parameter and validated the corresponding procedure by numerical calculations. In order to compare with results from the gradient expansion formalism in configuration space, we will also both compute ∆ numerically and consider it as a free parameter, although our final results will be based upon the dynamical calculation of ∆.

Slow Roll Analysis
The slow-roll inflation paradigm has been used by many authors to compute the amount of electromagnetic energy density [7,27,28], or baryogenesis through helicity [8,9,16,18] at the end of inflation, with or without taking into account the Schwinger effect. Here we aim to validate our numerical results by comparison with the known analytical results at the end of inflation.
In this section we will take ϕ as a slowly rolling inflaton field such thatφ ≃ 0, 3Hφ ≃ −V ′ (ϕ), and so we can considerφ and H = H E as constant. Doing so, we are neglecting the gauge field backreaction in the right-hand side of Eq. (2.2a), a hypothesis that we have consistently checked a posteriori. The results of this section will be model independent, within the hypothesis of the slow roll approximation.

Absence of Schwinger effect
Here we are assuming there is no Schwinger effect 7 , i.e. σ = 0, hence we can rewrite (2.7) as where, following the slow roll equations, is a constant. Since we are in de Sitter space, we can use the scale factor definition a = −(Hτ ) −1 and solve (3.1) asymptotically. At early time, when |kτ | ≫ 2ξ, the modes are in their BD vacuum given by (2.15), as here ∆ = 1. When |kτ | ∼ 2ξ, one of the modes develops both parametric and tachyonic instabilities leading to exponential growth while the other stays in the vacuum. During the last e-folds of inflation, i.e. |kτ | ≪ 2ξ, the growing mode has the solution [8,27] where a E and H E are, respectively, the scale factor and the Hubble parameter at the end of inflation.
Here, as we assume a slow roll regime, we consider H E constant and we take the convention a E = 1. Using (2.8) and (2.9) we can compute all electromagnetic quantities: These results are only valid when the absence of backreaction on the inflaton equation of motion (2.2a) is guaranteed, hence when |G/V ′ (ϕ)| ≪ f ϕ . This model-dependent condition puts a lower bound on the parameter f ϕ or, equivalently, a higher bound on ξ. Using the slow roll equations and the definition of the slow roll parameters, this parameter can be written as Therefore, at the end of inflation, where by definition ϵ = 1, one has ξ = M pl / √ 2f ϕ , and the no backreaction condition in Eq. (2.2) provides the bound ξ < 5.73 (or equivalently f ϕ > 0.12 M pl ). In Fig. 1 we show, with orange lines, the quantities ρ E , ρ B , H and G evaluated at the end of inflation obtained from the analytical backreactionless solutions from Eqs. (3.4), while the blue dots are the numerical solutions, which correspond to the case σ = 0 (no Schwinger effect) and correspondingly ∆ = 1. We have used a Runge-Kutta method which is explained in App. A.

Presence of Schwinger effect
The Schwinger effect is taken into account by means of the conductivity σ in Eq. (2.7), as given by Eq. (2.10) [19]. The growth of σ with time then backreacts on the gauge field, as the damping term grows in its differential equation. We will compare our numerical calculations with three analytical (or semianalytical) results: the Schwinger maximal and equilibrium estimates [19,24], as well as the gradient expansion formalism [22][23][24]. From the numerical point of view however, we aim to solve Eq. (2.7) with σ computed at each time step using (2.13). The details about the numerics will be displayed in Sec. 3.3.

Schwinger equilibrium estimate
In this case, the backreaction of the chiral fermions on the gauge fields is taken into account by just replacing the parameter ξ with the effective one [19] (3.6) Figure 1: Electric ρE and magnetic ρB energy densities, and the helicity H and its derivative G, at the end of inflation (i.e. for ϵ(aE) = 1), in units of HE, as functions of the coupling f ϕ assuming ∆ constant. We see the plots confirm the result from Fig. 1 of [24]. Here, we also assumed ξ constant.
in the backreactionless solutions (3.4). This amounts to solving which provides the function ξ eq = ξ eq (ξ) that we plug in (3.4) instead of the bare ξ to obtain the quantities ρ E eq , ρ B eq , H eq and G eq . These equilibrium estimates are shown with a purple line in the plots of Fig. 1.

Schwinger maximal estimate
In this case, we assume the exponential behaviors of the backreactionless solutions to be valid until they saturate the maximal value that we will display hereafter. We numerically determine the value of crossing, which happens for ξ ≃ 4.4-4.7 depending on each quantity.
The maximum helicity density can be estimated as the solution of [19] |E| 2 + |B| 2 = ξ eff |E| |B|. This replacement yields an equation relating the |E| and |B| fields that can be solved analytically. We then choose, as definition of our maximal estimate, the solution (|E|, |B|) of (3.8) that maximizes the product |E| · |B| 8 . This yields for ξ ≫ 1 The maximal estimates for the quantities ρ E max , ρ B max , H max and G max are shown with a pink line in the plots of Fig. 1.

Gradient expansion formalism
This method was introduced in Refs. [22][23][24] and transforms the EoM for the vector field A into EoM for observable quantities, in particular the electric E and magnetic B fields. As the spatial gradients in the EoM do always appear as rot E and rot B, the EoM can be written as an infinite series in terms of the bilinears E (n) = ⟨E · rot n E⟩/a n , G (n) = ⟨E · rot n B⟩/a n and B (n) = ⟨B · rot n B⟩/a n , with n = 0, 1, . . . . In this way the coupled system of EoM for the fields E and B transforms into a system of coupled differential equations for the quantities E (n) , B (n) and G (n) . This system is not block diagonal in the space of the n index so that the system has to be truncated to find solutions.
Moreover, the parameter ∆(t) in Eq. (2.17), which suppresses the gauge-field amplitude on small scales depends on the conductivity at all times t ′ < t. So, a precise determination of ∆(t) would require a complete analytical solution of the infinite-dimensional system of equations. While ∆ was dynamically computed in Ref. [23], for the sake of simplicity and generality, it was considered as a free parameter in Ref. [24] and fixed to the values ∆ = 1, 10 −2 , 10 −4 , 10 −6 . In our numerical approach we will consider ∆ as a function of the conductivity σ, as the initial condition for E and B are plane waves, such that E · B = 0 and therefore initially σ = 0 and so ∆ = 1. However, as time is evolving E and B will become collinear, and a nonvanishing conductivity will appear, as well as the function ∆(t) < 1. In order to compare our numerical results with those from Ref. [24], we also will eventually enforce ∆ to be a constant in our code. Upon considering a constant value of ∆, our results will agree pretty well with those obtained in the gradient expansion formalism, see Fig. 1. In the more realistic cases where we just compute the value of ∆(t), we will see that at the beginning t = t 0 , just very deep inside the inflationary period, ∆(t 0 ) = 1, while the value of ∆ will decrease very fast and at the end of inflation t = t E , ∆(t E ) ≪ 1.

Numerical results at the end of inflation
We will find it more convenient to change the variable from the time t to the scale factor a. The gauge field equation of motion (2.7) then becomes We recall that, as we are considering the slow roll regime in this section, we do not need to solve the equation of motion for ϕ.
The Bunch-Davis solutions can now be written as The technical details of the numerical simulations for solving Eq. (3.10), subject to the boundary conditions (3.11), can be found in App. A. We display in Fig. 2 the spectra of all the observable quantities in order to see how the BD vacuum is dominating the spectra for large k and how the cutoff k c (a), given by (2.19), efficiently removes that part of the integration. The difference between the BD vacuum and the damped BD vacuum is also clear, as the first goes like k 3 whereas the second goes like ∆(a)k 3 with ∆ decreasing with time. Hence the asymptotic behaviors are not superimposed since ∆ changes. Finally, we also see explicitly how the growth of ρ E and ρ B with the scale factor a is due to the increase in amplitude of the spectrum hump and its shift to larger values of k. For this illustrative purpose we used a constant value of ξ.
Here we have fixed f ϕ = 0.1 M pl , while for other values of this parameter the plots are similar. Before moving to the full numerical results, we will compare our slow roll based inflaton numerical results with the recent literature on the subject.

Constant ∆ and ξ approximation
We will first assume that the parameters ∆ and ξ are constants. As we already mentioned, the parameter ∆ was fixed to constant values in Ref. [24] while ξ, as defined in Eq. (3.2), is often considered as a constant in the slow roll approximation. In Fig. 1, we displayed several results already present in the literature that we successfully reproduced with our numerical method. First the backreactionless case, where there is no conductivity, by simply enforcing σ = 0 (therefore ∆ = 1) in the code. The data set are displayed in blue and match the corresponding analytical value given by Eq. (3.4). Then, in order to reproduce results from [24], we considered a non-zero conductivity given by (2.13) while assuming ∆ constant during inflation, thus making it a free parameter. In Fig. 1, we plot the quantities ρ B , ρ E , H and G at the end of inflation for chosen values of ∆. We can see that the results agree well with those using the gradient expansion formalism in Ref. [24].

Variable ∆ and ξ
The benefit of the slow roll approximation is that the results look "model independent". However, the tradeoff comes with the need of having a constant parameter ξ as the slow roll regime implies an approximately constantφ. Besides, we know that this parameter can also be expressed in terms of the slow roll parameter ϵ (see Eq. (3.5)), which is indeed small and constant during inflation but then quickly becomes unity during the last e-folds. We also know that the modes produced during the last e-folds are the ones that contribute the most to the integrals (2.8) and (2.9), as all the modes previously generated get washed out by the Universe expansion.
All these observations lead us to conclude that the most important contribution to the quantities ρ E , ρ B , H and G is taking place during an epoch when the constant ξ approximation loses its relevance. Hence, in this section, we will instead specify an inflation model, namely the Starobinsky potential, and make its study in the slow roll regime with a function ξ(a) that can be obtained from the model. We have chosen in this section the Starobinsky potential as it provides a realistic model of inflation, and will be a particular case of a more general class of models we will consider to make predictions using the full solution of the system. The purpose of this section will thus be to assess the goodness of the slow roll approximation when computing the full solution to the system (2.2).
The Starobinsky potential is given by (3.13) Using the slow roll regime, the inflaton field ϕ is given by where W n is the nth branch of the Lambert function. The value of the function ξ is then given by In Fig. 3, we display in blue results for the Starobinsky model, for various values of ϵ, when σ and ∆ vary dynamically. Although the slow roll approximation loses its relevance for values of ϵ closer to 1 (an issue we address in the next section), we already see a difference with the plots in Fig. 1. This is because, no matter the value of the initial time, the function ∆(a) rapidly goes to extremely small values, thus killing the BD modes that would have been amplified at the very end of inflation and that would have contributed the most to the integrals (2.8) and (2.9). With a constant ∆, this suppression is less effective and the tachyonic amplification yields higher energy densities and helicity. Figure 3: Comparison between the slow roll approximation and the full solution for the Starobinsky model. The analytical estimates are given for ϵ = 1. As expected, the slow roll computation diverges from the full solution as inflation is nearing the end, since the slow roll approximation is only valid in the regime ϵ ≪ 1. Hence the slow roll computation overshoots the value of all quantities, closer to the value given by the Schwinger equilibrium estimate for f ϕ ≲ 0.05 M pl . As expected, we also have compared both analysis, slow roll and full solution, for values of a such that ϵ(a) ≪ 1 (in particular ϵ = 10 −1 , 10 −2 , 10 −3 ) and found good agreement.

Full Analysis
In this section, we are not using the slow roll hypothesis for the inflaton equation of motion and consider the full solution to the system (2.2) in specific models of inflation. We will choose a set of inflationary models that are well known to be in agreement with all cosmological constraints. Also, we do not assume any peculiar geometry of the Universe.
The equations to be solved during inflation are the system (2.2) written in terms of the variable a. Unlike in the previous section, the current change of variables must take into account that the Hubble parameter is not constant, but moreover we have da dt =ȧ = aH, and we will define the auxiliary quantity F as We will relate it to the Friedmann equations where the total energy density and pressure are Hence we have (4.5) and the system (2.2) becomes The Hubble parameter can be computed from the Friedmann equation (4.2a), where ρ is given by (4.4a). This way, we can compute the value of H and F at each time step recursively to feed the equations of motion, like we already do for σ and G. The BD vacuum modes are identical to the previous case, see Eqs. (3.11). Finally, for comparison purposes, we can define a generalized time dependent instability parameter ξ(a) as such that it corresponds to the definition (3.2). The simulations show that this parameter, obtained from full solution computation, significantly differs from the slow roll one at the very end of inflation.

Full numerical results at the end of inflation
In this subsection we will compare our results at the end of inflation, where we are making a full numerical analysis of the EoM, with those obtained using the slow roll approximation for the inflationary potential. For the sake of comparison we will concentrate on the Starobinsky model given by (3.13).
In this current framework, we see in Fig. 3 that the four studied quantities, namely ρ B , ρ E , H and G, are much closer to the Schwinger equilibrium estimate at the end of inflation. We present in Fig. 3 the values of the physical observables evaluated at various stages of inflation, i.e. various values of the scale factor a, from ϵ(a) = 10 −3 to ϵ(a) = 1, as a function of the coupling f ϕ for the Starobinsky model. We superimpose the analytical results from Secs. 3.1 and 3.2, and hence the backreactionless solution as well as the Schwinger maximal and equilibrium estimates. From the plots we see that for f ϕ ≲ 0.05 M pl the equilibrium estimate is a good approximation, especially for ρ E where the predictions of maximal and equilibrium estimates merge. We also verify that cos θ ≃ 1 hence satisfying the assumption on parallel electric and magnetic fields leading to the conductivity definition (2.10).
In this setup, our numerical code is computing a value of the conductivity σ and ∆ for each time step, hence we got the functions σ(a) and ∆(a). The variation and presence of ∆(a) is not without effect on the final results. Indeed, the smallest (k ≫ H E ) modes are the ones that most contribute to the integrals (2.8) and (2.9). Without the Schwinger effect, these modes are produced last, just at the end of inflation, and only briefly leave the horizon. They therefore should have a significant impact on preheating. When the Schwinger effect prevents their generation, by reducing them by a ≪ 1 factor, while they are still in the BD vacuum, we can ask ourselves about the effectiveness of gauge preheating. It was shown in previous studies of gauge preheating [12] that its efficiency mainly depends on the electromagnetic energy fraction available at the end of inflation ρ EM /ρ tot . To shed light on the last point, we will extend, in the next section, our numerical results beyond the end of inflation when the inflaton is coherently oscillationg around its potential minimum. We will do that in a set of particularly interesting phenomenological models that we describe in the next section.

Inflationary models
We will here introduce two classes of models that all satisfy the cosmological constraints. They should be considered as a sample of possible models, and they are just chosen for illustrative purposes, as they do not exhaust by any means the allowed inflationary models.

α-attractor models
The α-attractor potential is given by [29] V α (ϕ) = Λ 4 (4.8) Setting α = 1 yields the R 2 model or Starobinsky potential (3.13). To make the comparison interesting, we choose to have 1 ≤ α ≤ 100, where cosmological observables are correctly reproduced. In the slow roll approximation, the field value at the end of inflation is We can readily compute ϕ * , and evaluate the slow roll parameters N * = 60 e-folds before the end of inflation. The slow roll parameters and the cosmic observables are in agreement with the cosmological contraints for the range 1 ⩽ α ≲ 100. Using the observed value of A s from Ref. [4], A obs s = 2.2 · 10 −9 , we fix the vacuum energy. The result depends on α and is approximately given by Λ α ≃ 3.4 · 10 −3 α 1/5 M pl . We then obtain the values Λ 1 = 3.152 · 10 −3 M pl and Λ 100 = 8.313 · 10 −3 M pl . Figure 4: Inflaton kinetic and potential energy density, as well as electric, magnetic and fermion energy density ratios to the initial total energy density of the Universe for the α-attractor models with α = 1 (upper panels) and α = 100 (lower panels). The vertical gray lines display the value a for which ϵ(a) = 1 and the dashed line shows the expected scaling of the dominant sector.

Hilltop quartic models
The hilltop model potential is given by [30] (4.13) The case p = 4 can be compatible with the Planck measurements. There are two ways for the field to relax to the minimum at ϕ = µ, with different initial conditions: 1. ϕ * > ϕ E : In this case the field ϕ > µ is relaxing in a potential region that can be approximated by V h ∼ ϕ 8 , and thus, the slow roll conditions are not met, as chaotic inflation is ruled out.
2. ϕ * < ϕ E : In this case the field ϕ < µ is relaxing in a flat potential region and the model predicts correct inflationary observables for a large range of the parameter. In this work, we will study this option.
The slow roll parameters and the cosmic observables are in agreement with the contraints for the range 10M pl ≲ µ ≲ 50M pl . (4.14) We fix the vacuum energy from the constraint on the amplitude of scalar fluctuations. The result depends on µ and is approximately Λ h ≃ 6 · 10 −4 µ 2/3 M

Numerical results beyond the end of inflation
Now that we have established a method to numerically compute the quantities ρ E , ρ B , ρ ψ , H and G, we aim to study the system evolution past ϵ = 1, and the onset of reheating. Indeed, the system (4.6) describes the most general interaction of the zero mode of both hypercharge gauge and inflaton fields. In particular, no assumption was made on the Universe geometry, hence there is no specific reason to stop its numerical computation at the end of inflation. We will also find it convenient to present some numerical results using as the variable the number of e-folds before the end of inflation N , instead of the scale factor a, and related to it by N = − log a E a (4. 16) such that N = 0 corresponds to the time a E when ϵ(a E ) = 1. We show the postinflationary energy breakdown for selected value of f ϕ , for the α-attractor models in Fig. 4, α = 1 (upper panels) and α = 100 (lower panels), and the hilltop models of inflation in Fig. 5, with µ = 10 M pl (upper panels) and µ = 50 M pl (lower panels). From the inflaton behavior, we see that the Universe enters a matter domination era as ρ ϕ ∼ a −3 . For high enough values of f ϕ , i.e. f ϕ ≳ 0.1 M pl , we reproduce the results shown in Ref. [12], whereas for f ϕ ≲ 0.1 M pl the electric and magnetic fields exhibit a different behavior: the former decays faster than the latter while oscillating. This is due to the fact, already mentioned in Ref. [23], that the energy density for the electric component E = −A ′ is much more sensitive to the Schwinger effect than the magnetic component B, because it directly couples to the conductivity in the gauge field equation of motion (2.7). On the other hand, the magnetic component reflects spatial effects, as it is defined by B = ∇∧A. In this work, we do not consider the inflaton spatial effects, ∇ϕ, because this would require one to implement real fermion interactions in a lattice simulation. Hence, for low values of f ϕ , when the Schwinger effect is strongly affecting the system, the behavior of ρ B is expected to be subject to changes when the spatial effects are enabled; namely we expect to see a faster decay, like that of ρ E . As also observed in Ref. [23], the electric field, which is dominant during inflation, becomes subdominant afterwards. Finally, we can see that for low values of f ϕ the fermion energy density dominates the radiation energy density at the end of inflation as already highlighted in Ref. [23].
The authors of Ref. [12] quote a sufficient criterion for gauge preheating to happen, namely that at least an 80% fraction of the total energy density of the Universe is electromagnetic energy. In the absence of the Schwinger effect, they found that this criterion is satisfied for values f ϕ ≲ 0.1 M pl . However, as expected, the Schwinger effect significantly reduces the share of electromagnetic energy, as shown on Fig. 6 for the considered models, which displays the ratio ρ EM /ρ total for the four previous considered cases. We can see that the maximum is attained with a value ∼ 10 −3 , which precludes any gauge preheating, at least for f ϕ ≳ 0.01 M pl . Another conclusion from Ref. [12] is that the spatial effects of the inflaton become relevant for sufficiently low values of f ϕ and contribute to preheating. Since we are neglecting them in our simplified calculation, any negative statement concerning the possibility of gauge preheating due to the lack of enough electromagnetic energy should be a conservative one. Figure 7: Maximum value of the electromagnetic to total energy fraction as a function of f ϕ for the four considered models: α-attractor models, with α = 1, 100, and hilltop models, µ = 10, 50 M pl . Preheating seems unlikely to occur.
The final results from our analysis can be summarized in Fig. 7, where we plot the maximum value of the electromagnetic to total energy fraction as a function of f ϕ (preheating efficiency) for the Starobinsky model, the α-attractor model with α = 100 and the hilltop models with µ/M pl = 10, 50. For f ϕ ≳ 0.01 M pl , we obtain ρ EM ρ tot ≲ 0.01, (4.17) which seems to prevent gauge preheating as its efficiency is far from the value of ∼ 0.8 established in the numerical analysis of Ref. [12].

End of reheating
If gauge preheating does not occur, the inflaton will eventually decay by perturbative processes which depend on the inflaton total decay width Γ ϕ . Therefore at the time t rh ∼ 1/Γ ϕ , the inflaton has completely decayed and the radiation domination era starts. Results from last sections have shown that shortly after inflation ends, the Universe is dominated by matter, hence we can approximate the Hubble parameter by where is the end value after reheating by inflaton perturbative decays. Of course a rh is a model-dependent quantity, which depends on the value of Γ ϕ , which in turn, depends on the couplings of the inflaton to the matter. In particular, the coupling 1/f ϕ of the inflaton to the hypercharge Chern-Simons density provides a channel for the perturbative decay of the inflaton into a pair of hyperphotons A, as ϕ → AA. This decay has a width given by [10] where m ϕ is the inflaton mass given by For the α-attractor (hilltop quartic) model, we have ϕ min,α = 0 (ϕ min, h = µ) and In the presence of extra couplings of the inflaton to matter, the predictions for the inflaton decay width, Eqs. (4.23) and (4.24), and the reheating temperature, Eqs. (4.26) and (4.27), will change in a model-dependent way, as well as the model predictions concerning the generation of the baryon asymmetry.
Of course in the hypothetical case where the explosive production of gauge fields should have prevailed over the perturbative inflaton decays, gauge preheating would have taken place over a few e-folds after the end of inflation. As we see from the previous results, this is never the case and gauge preheating is never strong enough to reheat the Universe after the period of cosmological inflation. This result does not preclude that, in the presence of a strong coupling λ of the inflaton with some other field, e.g. a scalar (or a fermion), there could exist an explosive production of that scalar (or fermion), triggering preheating of the Universe after inflation [31].

Baryon asymmetry
Before concluding this paper we wish to make a small comment on the baryogenesis issue at the electroweak phase transition. In Ref. [18], we presented a model of inflation that leads to a successful BAU. The effective potential for the inflaton, labeled therein as χ, was the Starobinsky potential 9 , and we did consider the Schwinger equilibrium and maximal estimates. Hence it is straightforward, using our numerical analysis in this paper, to make an update of the final results for the BAU for inflation driven by the α-attractor models with α = 1.
As all details are explained in Secs. 6 and 7 of Ref. [18], we skip them here and go straight to the final result. First of all we show in Fig. 8 the analogous plot to Fig. 9 of Ref. [18], namely the parameter space that provides a successful BAU. In particular, we display in blue the region where  9 In fact, we used in Ref. [18] an scalar field ϕ non-minimally coupled with gravity as L = − 1 2 gϕ 2 R + . . . , which yields for the canonically normalized field χ in the Einstein frame an α-attractor potential with α = 1 + 1 6g ∈ [4.3, 17.6], where the lower bound was coming from imposing the naive unitarity bound gϕ 2 < M 2 pl . As the dependence in α (hence in g) is tiny, we choose to show in the present paper the result for α = 1, hence for the Starobinsky potential (which would correspond in Ref. [18] to the limit g ≫ 1). the asymmetry parameter meets its observational value given by where we have imposed the observed value [32] in the right-hand side. Following Refs. [16,17] we define the parameter f θ W , which encodes all the details of the EW phase transition and its uncertainties, as In addition to their dependence on the gauge sector observables, the quantities used in this section vary according to the ratio of the reheating over the instant reheating temperature. This parameter hence adds to f ϕ in the parameter space. The reheating temperature is computed as where g * = 106.75 is the SM number of relativistic degrees of freedom, and we define T ins rh as a reference temperature given by the above equation with Γ ϕ ≃ H E , which is obtained from the simulation. It would correspond to the reheating temperature for instant reheating, and takes the value T ins rh ≃ 2.87 · 10 15 GeV. Using Eq. (4.26) it is possible to link the reheat temperature to the parameter f ϕ . The corresponding plot is shown in Fig. 8 which shows that it provides a wide window for baryogenesis.
Second, we display in orange the region where the magnetic Reynold's number at reheating R rh m is bigger than one, hence ensuring that the required magnetohydrodynamical conditions are fulfilled for the (hyper)magnetic fields to survive until the electroweak crossover. As we are in the viscous regime, it can be computed as [18] R rh m ≈ 5.9 · 10 −6 ρ B ℓ 2 where ℓ B is the physical correlation length of the magnetic field given by which can be numerically computed during the simulation in the same way as the other observables. Third, and last 10 , we show in green the condition on the chiral plasma instability (CPI) temperature, ensuring that the CPI time scale is long enough to allow all right-handed fermionic states to come into chemical equilibrium with the left-handed ones via Yukawa coupling interactions (so that sphalerons can erase their corresponding asymmetries in particle number densities) before CPI can happen. The estimated temperature at which CPI takes place is (4.33) The constraint T CPI ≲ 10 5 GeV (the temperature at which e R comes into chemical equilibrium) guarantees that the CPI cannot occur before the smallest Yukawa coupling reaches equilibrium and 10 Besides, we checked that the generation of baryon isocurvature perturbation provides no constraint.
all particle number density asymmetries are erased, preventing thus the cancellation of the helicity generated at the reheating temperature. Therefore, as we can see from Fig. 8, the resulting baryogenesis window for the Starobinsky potential is close to the Schwinger equilibrium estimate for f ϕ ≲ 0.06 M pl , just as the corresponding results on helicity and magnetic energy density suggest (see the green dots of Fig. 3). However, for f ϕ ≳ 0.06 M pl there is no space for the BAU, as the production of gauge fields is too weak, unlike in the previous results from Ref. [18]. In addition to this we have seen that the reheating temperature is constrained by the model, see Eq. (4.26), as we can see from Fig. 8 and compatibility of the model reheating temperature with the baryogenesis results translates into the baryogenesis region on the parameter f ϕ f ϕ ≲ 0.03 .M pl (4.34) Finally, one of the results of this paper is then that baryogenesis at the electroweak phase transition is favored by low reheating temperatures, in the range 10 −6 T ins rh ≲ T rh ≲ 10 −3 T ins rh .

Conclusions
In this paper, we have studied by means of numerical computations the effect of the Schwinger particle production on the helical hypermagnetic fields produced at the end of inflation. The inflaton field ϕ can decay, through its coupling to the Chern-Simons density ϕ 4f ϕ Y µνỸ µν , into helical hypermagnetic fields in a nonperturbative process. When exiting the vacuum, the gauge modes are strong enough to create particle/antiparticle pairs of light fermions, which contribute to the electrical conductivity of the plasma. The backreaction of fermion currents on the produced gauge fields acts as a damping force in the explosive production of helical gauge fields. This effect, called Schwinger effect, was already considered in numerous studies of inflation and/or baryogenesis, where some analytical and numerical estimates were computed, mainly in configuration space while our calculation is done in momentum space.
The equations of motion are in fact a nontrivial integro-differential system. It was solved numerically by using a fourth order Runge-Kutta method, with details being displayed in the Appendices. The computed observables of interest are the electric and magnetic energy density, the helicity as well as the helicity time derivative. We assumed a homogeneous inflaton with only zero mode, hence we did not treat any spatial effects. Besides, we also ensured the convergence of the algorithm and its invariance to the initial conditions.
First of all we have checked that we recover previous results in the slow roll inflation regime by making the same approximations required by an analytical resolution. In this way, we validate our code, i.e. we verify that our code produces the right results in known cases such as the backreactionless case, where the Schwinger effect is turned off, and the gradient expansion formalism, where the Bunch-Davies parameter ∆ was first introduced.
In a second step, still in the slow roll regime, we considered a specific model of inflation, namely the Starobinsky potential, in order to account for the instability parameter as a function, ξ(a), instead of the constant imposed by the analytical approximations. That way, we could also implement the effects of a function ∆(a) obtained from the plasma evolution on the gauge production itself.
We then simulated, in a third step, the full system, where neither the slow roll conditions nor the Universe geometry (e.g. de Sitter) are imposed. In order words, the inflaton equation of motion was computed alongside with the gauge one, taking the backreaction of the latter to the former into account along with the Schwinger effect. We compare our result to the previous setup and found perfect agreement as long as the slow roll conditions are met. When inflation is near its end, the full solution diverges from the slow roll results and produces, as expected, less energy density and helicity.
Finally we will comment on the implication about two related topics: gauge preheating and baryogenesis. As our code is free from any geometrical issues, and only requires a model of inflation, we let the simulations run until the onset of reheating to compute the electromagnetic to total energy density ratio. We choose two well-known classes of models that satisfy the cosmological constraints as illustrative examples. Previous studies have quoted a sufficient criterion for gauge preheating to happen, namely that this fraction should be at least ≳ 80% [12]. However, our numerical estimates suggest that the Schwinger effect significantly reduces the share of electromagnetic energy for the considered models and preheating is unlikely to occur. Moreover, since we are neglecting all spatial effects, any negative statement concerning the possibility of gauge preheating due to the lack of electromagnetic energy should be a conservative one. On the other hand our results do apply to the considered class of inflationary models. They show a certain degree of model dependence, so we cannot exclude a qualitatively different result for models of inflation other than the considered ones.
On the other hand, as a successful baryogenesis does depend on a delicate equilibrium between the amount of helicity, magnetic energy density, and magnetic correlation length, damped fields do not necessarily mean no baryon asymmetry in the late Universe. Actually, as a result of our numerical calculation, we have found there is still a window in the parameter space for baryogenesis to happen as long as f ϕ ≲ 0.05 M pl , while consistency from the perturbative decay channel of the inflaton into hypergauge bosons implies the bound f ϕ ≲ 0.03 M pl . Moreover, baryogenesis is favored for low enough values of the reheating temperature T rh ≲ 10 −3 T ins rh . Of course, the baryogenesis predictions should, to some extent, depend on the model of inflation. In this way, our result here is restricted to the Starobinsky model and should be considered just as a "proof of existence" for baryogenesis in the presence of the Schwinger effect.
These two comments should be viewed as hints for future studies that address the production of gauge fields at the end of inflation. Of course, a full lattice simulation of the Schwinger effect involving fermions remains to be done.
Eq. (3.10) becomes the following system: To perform each time step ∆a, we use the fourth order Runge-Kutta (RK4) algorithm: Note that x is complex, hence we solve the above system for both real and imaginary parts but with their specific initial conditions. These are mode-dependent as it takes longer for modes with bigger wave number to leave the BD vacuum. Therefore, we choose as initial condition for each mode where we choose the factor x BD in order to make sure that we initialize the gauge field sufficiently deep inside the Hubble radius. Its exact value is subject to analysis and is discussed later. As we can see from (2.8) and (2.9), high values of k are dominating the integral hence large modes are negligible compared to small ones. This makes us to choose a lower bound on the k range such that the initial time of the simulation is a 0 = k min x BD . (A.5) In that way, at a 0 we make sure that all the modes are in their respective vacua, which implies σ = 0 as explained above.
In practice, this means that the modes with k > x BD a are given by the following relations so that the discretization is the same for each order of magnitude. This means ∆a grows exponentially with a. The advantage of this method is that there is a refinement of the grid for small values of a, at the beginning of inflation. The same is done for the discretization in k.
We explored the numerical convergence of the solution, both in the number of a i 's, labeled as N a , and in the number of k j 's, labeled as N k . Provided that N a > 2000 and N k > 200, the simulations are very stable and the output does not depend on the discretization. For big values of f ϕ , f ϕ ≳ 0.1, we can even lower the number of time steps needed.
Besides, we must choose the BD penetration factor x BD such that it produces trustable results. We have done a numerical analysis and conclude that depending on the value of N a , a range 20 < x BD < 50 yields trustable results. We hence choose throughout this work the following values x BD = 20, N a = 500, 1000, 2000, N k = 300. (A.8) At each time step, we compute the electric and magnetic energy density as where we choose such that we cut off the spectra to retain only modes outside the horizon. The helicity (2.9a) and its derivative (2.9b) become In the numerics, these integrals are performed numerically over the range of k that takes N k discrete values. If the Schwinger effect is taken into account, we turn on the possibility of having σ computed at each time step a i of the numerical computation with σ i+1 = 41 g ′3 72 π 2 a i 2ρ i B coth π ρ i B ρ i E (A. 12) and injected into the calculation of the next step. Otherwise, we keep it zero. Last, the fermion energy density is computed as Re(x λ i )Re(y λ i ) + Im(x λ i )Im(y λ i ) .

B Numerical method: full analysis
The numerical implementation follows from the previous case. Defining the variables w = ϕ, x = dϕ da , y λ = A λ , z λ = dA λ da (B.1) we transform the above coupled system of differential equations (4.6) into the system which is equivalent to writing dx da = f (a, x).

(B.3)
We recall that w, x ∈ R and y λ , z λ ∈ C. Similarly to the previous calculation with the slow roll approximation, we use the RK4 algorithm (A.3) with the values of H, σ, F and G computed at each time step. Inflaton initial condition could be set to However, the number of e-folds sets the initial time as a 0 = e −|N * | ∼ 10 −26 , which is too small a number for the numerical implementation. We then proceed as follows. For a ≲ k min /x BD , and sufficiently low k min , the gauge field modes stay in their vacuum and the total contribution to ∆(a) is negligible. Hence we do not need to perform the numerical simulation before that time, as the inflaton is the main player, so we can solve its equation of motion analytically. Instead, we fix the start of the simulation like before, at a 0 = k min /x BD and we compute the corresponding number of e-folds N which leads us to the corresponding value of ϕ(N ). Therefore, the initial condition must be set to w 0 such that and, usingφ ≃ − V ′ (ϕ) 3H which is valid at the early stages of inflation, As for the gauge field, initial conditions are set in the same way as in the slow roll approximation, see App. A.