Chiral anomaly, Schwinger effect, Euler-Heisenberg Lagrangian and application to axion inflation

Particle production in strong electromagnetic fields is a recurring theme in solid state physics, heavy ion collisions, early universe cosmology and formal quantum field theory. In this paper we discuss the Dirac equation in a background of parallel electric and magnetic fields. We review the Schwinger particle production rate, clarify the emergence of the chiral anomaly equation and compute the induced current of charged fermions. We distinguish the contributions from non-perturbative particle production, from the running of the gauge coupling constant and from non-linearities in the effective QED Lagrangian, and clarify how these contributions arise within a single framework. We apply these results to axion inflation. A Chern-Simons coupling between the pseudoscalar particle driving cosmic inflaton and an abelian gauge group induces a dual production of gauge fields and charged fermions. We show that the resulting scalar and gravitational wave power spectra strongly depend on the fermion mass.


Introduction
The production and motion of charged fermions in (strong) electric and magnetic fields is a recurring phenomena in many different areas of physics. Schwinger production describes the pair creation of particles and antiparticles in strong electric fields [1,2], exponentially suppressed by the particle mass. Adding a magnetic field, the particle mass is replaced by the effective magnetic mass, quantized in Landau levels [3,4] (a phenomena well known in solid state physics in the context of the Quantum Hall Effect [5]). In the limit of massless fermions, the lowest (gap-less) energy level induces asymmetric particle production. This in turn is a beautiful realization of the chiral anomaly in quantum field theory (QFT) (Adler Bell Jackiw anomaly [6,7]) in the presence of helical electric and magnetic background JHEP02(2020)055 fields [8][9][10], discussed first in [11] in the context of Weyl fermions in a crystal. These arguments have been extended to non-abelian gauge fields in [12].
Going beyond solid state physics and formal QFT, particle production and transport phenomena associated with the chiral anomaly also play an important role in the study of the quark gluon plasma (QGP) and in cosmology. In heavy ion collision experiments, intense magnetic fields and local P -/CP -violation are expected [13]. In the presence of a magnetic field on top of the chiral imbalance, a separation of electric charge is induced by the chiral magnetic effect [14][15][16][17][18][19], which is a signature of the local P -/CP -violation [20][21][22]. Turning to cosmology, strong helical gauge fields appear in models of axion inflation [23][24][25], are postulated to permeate intergalactic voids [26], may play a crucial role in baryogenesis [27][28][29][30][31][32][33] and also appear in some implementations of the relaxion model [34][35][36][37][38], constructed to address the hierarchy problem in the Standard Model (SM) of particle physics.
All these examples share the same underlying physical principles. The gauge field background leads to a time-dependent dispersion relation for the fermions, inducing particle production. In the presence of a magnetic field, the dispersion relation features discrete Landau levels. The produced fermions are accelerated along the field lines of the background gauge field, leading to an induced current [39,40]. In general, the final expression for the induced current receives contributions from fermionic excitations, from the vacuum (fermion loops) as well as from non-linearities in the gauge field induced by virtual charged particles (described by the Euler-Heisenberg action [1]). If the gauge fields are not external fields but treated dynamically (as in realistic cosmological models), this induced current will inhibit the gauge production.
Various aspects of this picture have been studied in the literature. The production rate of (massive) fermions in parallel, static electric and magnetic fields was derived in [3,4] (see also [41]). The induced current resulting from this particle production was discussed e.g., in [42]. Our work is moreover closely related to ref. [8], which discusses the contribution to the induced current sourced by particle production in the context of heavy ion collisions, using a methodology similar to ours. For a discussion of the Euler-Heisenberg terms see [41]. For a recent analysis of the induced current in scalar quantum electrodynamics (QED) in de Sitter space without magnetic fields, clarifying nicely all the contributions to the induced current, see [43].
In this work, we compute the full expressions for fermion production and the induced current in a background of constant, helical electric and magnetic gauge fields. We keep the fermion mass as a free parameter and work in conformal coordinates, ensuring that our results can be applied both to flat space time and de Sitter space time. We derive and solve the equation of motion for the fermions (and equivalently for the Bogoliubov coefficients), accounting for the non-perturbative result for fermion production. We clarify how the chiral anomaly equation arises [8,10], connecting smoothly to the results obtained in the case of massless fermions [9,11]. Based on these results, we compute all contributions to the induced current: (i) the current from particle production, (ii) the running of the gauge coupling, and (iii) the 1/m suppressed terms of the Euler Heisenberg action. The framework presented here, based on solving the Dirac equation of motion in a background of helical electric and magnetic fields, captures all these different aspects of fermions in abelian gauge theories in a simple and unified description.

JHEP02(2020)055
The resulting expression for the induced current allows us to tackle applications in cosmology. As an example, we study the case of axion inflation. Here, the rolling pseudoscalar inflaton field induces a tachyonic instability in one of the gauge field helicities. The naive exponential gauge field production is significantly inhibited by the induced current for fermion masses smaller than O(100) H, with H denoting the Hubble rate. Consequently, the resulting scalar and gravitational wave power spectrum at small scales become highly sensitive to the fermion mass, providing a unique opportunity to test the microphysical implementation of these inflation models.
The remainder of this paper is organized as follows. In section 2 we discuss the fermion equation of motion, introducing the notion of discretized Landau levels in the dispersion relation and clarifying some subtleties about the definition of particles and antiparticles in a gauge field background. This enables us to study particle production in section 3, including the derivation of the anomaly equation and the derivation of all the contributions to the induced current. In section 4 we turn to axion inflation as an example for the relevance of these computations in early Universe cosmology. We conclude in section 5. Several important but rather technical steps are relegated to the appendices: appendix A clarifies our notation and conventions, appendix B provides some detail on the computations of section 2 whereas appendix C contains some of the necessary computations for section 3. Finally, appendix D contains the weak field expansion necessary to understand the vacuum contributions to the induced current.

Solving the fermion equation of motion
This section is dedicated to the analysis of the Dirac equation in a background of classical (anti-)parallel electric and magnetic fields. In particular, we identify the eigenvalues (energy levels) and eigenstates (particles and anti particles), both in the absence and presence of an electric field. This sets the stage for computing particle production and the different contributions to the induced current in section 3.

Preliminaries
The fermion equation of motion. A Dirac fermion ψ of mass m which is charged under an Abelian gauge group obeys the equation of motion with / D = (∂ µ + igQA µ )γ µ , Q denoting the charge of the fermion ψ under the abelian gauge group with vector potential A µ and gauge coupling g, and a denoting the scale factor of the FRW metric. For more details on our notation and conventions, see appendix A.
It will prove easier to solve the second order differential equation obtained by acting on eq. (2.1) with the differential operator (i / D + ma). 1 After some algebra (see appendix B),

JHEP02(2020)055
we obtain 0 = (i / D + ma)(i / D − ma)ψ (2.2) Here E and B denote co-moving 'electric' and 'magnetic' fields associated with the vector potential A, i.e. E = −∂ 0 A and B = ∇ × A, σ contains the Pauli matrices and we denote the derivative with respect to conformal time with a prime. From this point on we will ignore the expansion of the Universe, a = 0, assuming it to be slow compared to all relevant microphysical processes (see section 4.1 for details). Moreover, motivated by the applications discussed in section 4, we will consider constant (anti-)parallel electric and magnetic fields. This in particular implies a non-vanishing Chern-Simons term, E · B = 0. Without loss of generality, one may take the z-axis to match with the direction of the electric field: Here the sign of the magnetic field, B ≷ 0, encodes the two possible configurations (parallel vs antiparallel), related by a CP transformation. Note that, throughout this paper, we take temporal gauge for convenience. Inserting eq. (2.4) and performing a Fourier transform in the y-and z-direction, we can simplify eq. (2.3) to Here and hereafter, we mostly drop the k y and k z arguments inψ for brevity.
Mode expansion. Let us take a spinor basis which diagonalizes the last term in eq. (2.6): where with σ = ± and H = L, R describes left-and right-handed particles, respectively. Anticipating that the differential equation (2.6) is separable, we can expandψ as follows [9,11,40]:

JHEP02(2020)055
The x-dependent part of eq. (2.6) describes a re-scaled harmonic oscillator with shifted zero-point, where we have defined The solutions of eq. (2.10) can be expressed in terms of Hermite polynomials H n , where have normalized h n (x) such that dx h n (x)h n (x) = δ nn and the introduction of the discrete parameter s ensures that all minus signs are kept track off in the manipulations of eq. (2.10). After inserting eq. (2.9) and (2.12) into (2.6), the t-dependent part of (2.9) is given by the solution of We note the appearance of the discrete energy levels (2n + 1)sgQB labeling the solutions of the harmonic oscillator, see eq. (2.10). The magnetic field confines the trajectories of the charged particles projected onto the x-y plane to circular orbits, resulting in quantized energy levels, commonly referred to as Landau levels.
In the following section 2.2, we first turn off the electric field and study the dispersion relation in the presence of the magnetic field only. For a vanishing electric field, i.e. for constant A z , one may easily solve eq. (2.13) by g H n,σ ∝ e ∓iΩn,σt , where the dispersion relation is given by The lowest energy solution is obtained from n = 0 and σ = s, which is referred to as the lowest Landau level. For the higher energy levels, referred to as higher Landau levels, we note the degeneracy Ω n,s = Ω n−1,−s , for n ≥ 1. Taking into account left-/right-handed particles, H = L, R, and positive/negative energies, one might think that there are four independent solutions for the lowest Landau level and eight for each higher Landau level. However, only a half of the solutions for the second order differential equation (2.3) actually solve the equation of motion (2.1) which is first order, as we will see in section 2.2. Hence, the number of true independent solutions is two for the lowest Landau level and four for each higher Landau level. In the next section 2.2 we first discuss the lowest Landau level and then move on to the higher Landau levels. JHEP02(2020)055

Landau levels
Lowest Landau level. In the absence of an electric field, i.e., for constant A z in eq. (2.4), the operator on the left-hand side of eq. (2.13) becomes independent of the index H. As we have mentioned, the terms proportional to the B-field cancel and the energy is minimized for n = 0 and σ = s: These solutions span an orthonormal basis for the wave functions in the lowest Landau level, i.e., dx u † It is instructive to consider the limit of chiral fermions, i.e. taking the limit m → 0 in eqs. (2.18) and (2.19): where θ denotes the Heaviside function. From this equation, one can see that the dispersion relation for the left-/right-handed fermions (proportional to χ L s , χ R s respectively) is given by ω L/R 0 = ∓sΠ z , reproducing the well known result for the lowest Landau level [9,11]. This is also illustrated in figure 1.
Higher Landau levels. For a constant A z , the operator on the left-hand side of eq. (2.13) becomes independent of the index H, and both g H n,s and g H n−1,−s solve the same equation: The solution can be easily obtained as (2.22) Here we have defined the transverse mass m T for later convenience. We get eight independent solutions for each n, h n g H n,s χ H s and h n−1 g H n−1,−s χ H −s , which solve the second order differential equation (2.6). However, only four particular linear combinations actually solve JHEP02(2020)055  Figure 1. Energy levels for fermions in the presence of background gauge fields. The solid blue lines denote the lowest Landau level n = 0 for a massive fermion, the dotted lines correspond to higher Landau levels n ≥ 1. The gray lines indicate the corresponding energy levels for massless fermions, for which in the lowest Landau level is asymmetric between left-handed (L) and right-handed (R) fermions. Note that all higher Landau levels (dotted curves) have a two-fold degeneracy, i.e., the dispersion relation is identical for left-and right-handed particles. For the purpose of these plots, we use units where g|Q|E = 1. In these units, we set ma = 0.5, 2g|QB| = 0.7 and t = 0 (t = 0.2) in the left (right) panel.
the Dirac equation (2.1): for positive frequency modes; and for negative frequency modes. Note that due to the degeneracy in Ω n , any linear combination of u (1) n and u (2) n (and correspondingly of v (1) n and v (2) n ) solves the Dirac equation.

JHEP02(2020)055
The basis above is chosen such that i) for each higher Landau level n, these solutions span an orthonormal basis for the wave functions, i.e., dx u

Particles and anti-particles
We here clarify the definition of positive and negative energy states without assuming a vanishing electric field. This discussion is helpful, in particular, when we turn on the electric field and discuss particle production in section 3. To define the positive and negative energy states unambiguously, we need to look at the Hamiltonian for the Dirac fermion, which is given by (2.27) Inserting the mode expansion given in eq. (2.9) into this equation, one may express it in the following matrix form: where m B a ≡ √ 2nsgQB. We write the time argument explicitly in Π z to emphasize that this expression is applicable forȦ z = 0. From this equation, it is clear that the n = 0 and n ≥ 1 solutions never mix, and hence we can study them separately.
Let us first consider the first two-by-two matrix. Its eigenvalues are ±Ω 0 with Ω 0 (t) = Π 2 z (t) + m 2 a 2 . (2.30) The positive and negative eigenvalues correspond to particles and anti-particles, respectively. The corresponding eigenvectors are given by for the positive and negative eigenvalues respectively. These eigenvectors form the basis for the wave functions for the positive and negative energy states: (2.32)

JHEP02(2020)055
One can easily see that they coincide with the solutions for the equation of motion for a vanishing electric field as expected, i.e., U 0 = u 0 and V 0 = v 0 forȦ z = 0. This means that the positive and negative energy modes provide distinct solutions for the equation of motion forȦ z = 0, i.e., they do not mix with each other. However, this is no longer true once we turn on the electric field. There the positive and negative frequencies get mixed during the course of evolution due to the time-dependence in Π z (t), resulting in the particle production as we will see in section 3.
Next, we move on to the four-by-four matrix. Its eigenvalues are degenerate, i.e., we find two positive and two negative energy states, ±Ω n with The corresponding eigenvectors for the positive energy state are (2.34) Those for the negative energy state are (2.35) From these eigenvectors, we obtain the following wave functions for the positive and negative energy states:   n forȦ z = 0. Once we turn on the electric field, eq. (2.36) no longer solve the equation of motion, implying the particle production. Now that we know the wave functions for the positive and negative energy states in the presence of the electric field, we can for any time t define creation (and annihilation) operators for particles and anti-particles,B expand the fermion field at any given time t as follows: (2.37) Note that the following rules are implicit to get this formal expression: U Here creation and annihilation operators fulfill the following commutation relations:

Particle production
To study the fermion production induced by the electric field, let us turn on the electric field adiabatically for a finite amount of time. This allows us to define eigenstates with positive/negative frequency at early and at late times, and avoid spurious behavior associated with sudden turn on/off of the electric field. More concretely, we turn on the electric field from t = 0 to t = τ .
Note here that, to avoid spurious singularities at t = 0 and t = τ , we interpolate these discontinuities smoothly as shown in figure 2. JHEP02(2020)055

Bogoliubov coefficients
Definition. We consider solutionsũ n for the equation of motion which fulfill the following boundary condition:ũ n | Az=0 at t 0. Note here that we have used the same convention for the lowest Landau level as eq. (2.37), i.e., u Initially, there is no electric field and hence one may expand the field by means of wave functions defined in eqs. (2.18), (2.19) and (2.23) to (2.26) at A z = 0 together with the creation operatorsb where the annihilation operators erase the initial vacuum state, i.e., b n , must be expressed as a linear combination: Here we have defined the Bogoliubov coefficients: Intuitively, the α-coefficient describes the probability of a particle-to-particle transition, whereas the β coefficient encodes a particle-to-antiparticle transition. Complex conjugation exchanges particle and antiparticle states on both sides of these transitions. Note here again that, by convention, for the lowest Landau level, we only have r = 2, i.e., α Note that the Bogoliubov coefficients satisfy |α (r) n | 2 = 1. From the discussion above, it is clear that the annihilation operators at a later time t > 0 do not erase the initial vacuum state. More explicitly, the number densities become non-vanishing for |β (r) n | > 0 implying the particle production: Evolution of Bogoliubov coefficients. Instead of obtainingũ  where the initial condition at t 0 is α (r) n = 1 and β (r) n = 0. One can see that the structure of the equations is the same as those for the ordinary Schwinger effect without B-field by just replacing m B a with k 2 x + k 2 y [44]. For later convenience, we rewrite eq. (3.7) as follows: From a comparison with the ordinary Schwinger effect [44], we obtain the resulting asymptotic, non-perturbative behavior of |β (r) n | 2 at t τ for |k z | ma and |Π z | ma as where we have introduced p = sgn(Q). To avoid this notational complication, we will sometimes take Q > 0 and B > 0 in the following. Be careful that, in this case, all the quantities in the exponent should be regarded as their absolute values if one would like to recover general signs for them. Note that the inherently non-perturbative solution (3.9) can never be obtained by a perturbative expansion in E/m 2 a 2 . Vice versa, the full solutions to eq. (3.8) will contain perturbative corrections of order E/m 2 a 2 to the asymptotic solution (3.9). We will return to these in section 3.3 and appendix D. Figure 3 shows |β (r) n | 2 as a function of k z as a result of a numerical evaluation of eq. (3.8).

Chiral anomaly and chiral asymmetry
Deriving the chiral anomaly. In the presence of a chiral anomaly, a chiral fermion rotation relates the Chern-Simons term and the chiral fermion current [45,46]. For a JHEP02(2020)055 massive fermion, this chiral fermion rotation also transforms the fermion mass term, In the setup at hand, the (anti-)parallel electric and magnetic fields induce a non-vanishing Chern-Simons term, FF = −4E · B. Consequently, the fermion production discussed in the previous section must obey this anomaly equation. In the following we will verify this explicitly. This is not only an important and non-trivial cross-check of the consistency of our treatment, but also clarifies how the anomaly equation for a massive fermion arises (see also ref. [8] as well as refs. [9,11,12] for a derivation for massless fermions). As we shall see in the following, only the lowest Landau level contributes to the anomaly equation. Let us begin by evaluating the left-hand side of eq. (3.10). Taking the spatial average and using the translational invariance of the initial vacuum state, we are left with only the µ = 0 component J 0 5 , which we can identify as the chiral charge measuring the number density of right-handed minus left-handed particles. Here . . . indicates the vacuum expectation value with respect to the initial vacuum at t 0. 2,3 Inserting eqs. (3.2) and (3.4) into (3.11), we obtain (see appendix C for an explicit derivation and the reason why one may drop the regularization): Since contributions from the higher Landau levels cancel out between the two modes, r = 1, 2, only the lowest Landau level contributes to the chiral charge. In the next subsection, we will estimate q 5 for t τ by inserting the solution of the evolution equation for the Bogoliubov coefficients. Here we leave this expression as it is and proceed to the proof of the anomaly equation (3.10).
To evaluate the left-hand-side of the anomaly equation (3.10), let us take a time derivative of q 5 given in eq. (3.12). By using the equation of motion for the Bogoliubov coefficients in eq. (3.8), we obtaiṅ Note that the spatial average is redundant when assuming the initial vacuum state has translational invariance, and hence the chiral charge is essentially q5 = J 0 5 . Nevertheless we mostly keep it to avoid confusions. 3 We have here utilized the antisymmetrized chiral current. This is always possible by adding a total derivative to the action. This expression is useful because it makes a CP -invariant regularization explicit and allows in general a more transparent treatment of the η-invariant [47,48] (see ref. [12] for details). In the case at hand the η-invariant vanishes and hence this is not strictly necessary as we will see, but we will stick with the antisymmetrized version anyway. Recalling the definition of the dispersion relation, Ω 0 = Π 2 z + m 2 a 2 , and changing the integration variable k z → Π z /ma ≡ u, one may integrate the first term in eq. (3.13) analytically, which provides the Chern-Simons term: (3.14) The last step of the proof is to show the second term in Eq (3.13) coincides with the expectation value of the pseudo scalar operator 2ma ψiγ 5 ψ . Again, inserting eqs. (3.2) and (3.4) into 2ma ψiγ 5 ψ , we get which corresponds precisely to the second term in eq. (3.13). Combining eqs. (3.13), (3.14), and (3.15), we finally arrive aṫ which is nothing but the anomaly equation (3.10). Here, for notational brevity, we have dropped the spatial average on 2ma ψiγ 5 ψ . This derivation clarifies that, if one studies the particle production by solving eqs.  pumping up an electric field from t = 0 to t = τ as in eq. (3.1). In particular, we will clarify how each term in the anomaly equation contributes to the final asymmetry in the chiral charge with the help of a numerical evaluation of eq. (3.8). We take Q > 0 for notational brevity of the following discussion. Let us start with the chiral charge, eq. (3.12). To get the final chiral asymmetry, we need to evaluate the chiral charge at a late time, i.e., q 5 (t τ ). For t τ , the electric field is already turned off and hence the Bogoliubov coefficients are no longer evolving. Contrary to the first term in eq. (3.12), one can show that the second term, dk z (ma/Ω 0 ) (α * 0 β 0 e 2i dt Ω 0 ), does not contain a term which grows with τ . We have verified this numerically. Hence using the asymptotic non-perturbative behavior for |β 2 | given in eq. (3.9), we can estimate the resulting chiral charge as which is consistent with the argument given in refs. [8,49]. We also check this estimation by solving eq. (3.8) numerically, as depicted in figure 4. Taking the massless limit m → 0, one recovers the well-known result for chiral fermions in ref. [11]. It is instructive to see how each term in the right-hand-side of eq. (3.16) contributes to this resulting chiral asymmetry. The Chern-Simons term obviously yields To estimate the contribution arising from 2ma ψiγ 5 ψ , we need to evaluate the imaginary part of α * 0 β 0 e 2i dt Ω 0 as eq. (3.15). As can be seen from figure 5, for a long enough τ such that |gQEτ | ma, the time integration of the pseudo-scalar operator, 2ma ψiγ 5 ψ is JHEP02(2020)055 dominated by the plateau from t = 0 to t = τ and hence we get which coincides with the result in refs. [8,10]. In the massless limit, this term vanishes as expected. Now we see that the estimation for the final chiral charge in eq. (3.17) is consistent with the summation of eqs. (3.18) and (3.19). We note in particular the non-trivial "-1" contribution in eq. (3.19), which cancels the FF contribution to the chiral charge, ensuring the overall exponential suppression. Since due to the operator nature of the anomaly equation this cancellation must hold exactly, one can also see how this term emerges explicitly by perturbative computations omitting the non-perturbative exponential behavior in the limit of m 2 a 2 E. As shown in appendix D, the leading order perturbative contribution is It shows that the perturbative contribution of the mass term indeed cancels the FF term at the leading order, which is essential for the anomaly equation as discussed above.

Induced current and Euler-Heisenberg action
Once the fermions get generated, they are accelerated in the background gauge field and tend to erase it. On top of this, even if the fermions are so heavy that their production is not efficient, the coupling to the fermions changes the effective action for the gauge field via the running of the gauge coupling and via higher dimensional operators suppressed by the fermion mass. These contributions correspond to the Euler-Heisenberg action [1], which describes the non-linear dynamics of QED induced by one-loop effects.
In this subsection we explicitly show that the expectation value of the current contains all the effects mentioned above. More concretely, we will evaluate Induced current. Let us first discuss the contribution from the particle production. The effect of particle production is imprinted in the non-perturbative part proportional to e −πm 2 a 2 /gQE . As discussed above eq. (3.17), dΠ z (m T a/Ω n ) (α * βe 2i dt Ω ) does not contain any term which grows with τ . As a result, we arrive at the following expression of JHEP02(2020)055 the current induced for a large enough τ by the particle production after pumping up the electric field as eq. (3.1): reproducing [8]. The subscript "ind" means the current induced by the particle production.
We have multiplied a factor gQ in order to make the structure which couples to the gauge field explicit. Note that contrary to the computation of the chiral charge, all Landau levels contribute to the induced current. This result is also consistent with refs. [3,4,41] obtained from evaluating the imaginary part of the effective action after integrating out fermions. See also ref. [44] to see how to relate their results to ours.
Vacuum contribution. As discussed above, the contribution from the particle production is exponentially suppressed for heavy fermions, i.e., e −πm 2 a 2 /g|Q|E . However, this does not mean that those heavy fermions don't contribute to J z . One may anticipate their effect because the low-energy effective action for the gauge field after integrating out heavy fermions should involve the threshold correction to the running of the gauge coupling as well as higher dimensional operators suppressed by the heavy fermion mass. In the following we explicitly compute these effects by means of a perturbative expansion in E/m 2 a 2 and B/m 2 a 2 . Here we ignore terms that contain more than one time derivatives in E. These terms generate higher derivative corrections to the Euler-Heisenberg terms in the low-energy effective action. See also ref. [43] for the case of a vanishing magnetic field.
As discussed in appendix D, in addition to the non-perturbative term from the particle production, we also have perturbative contributions in |β| 2 and α * βe 2i dt Ω . By utilizing integration by parts repeatedly, one may systematically perform perturbative expansions in E/m 2 a 2 and B/m 2 a 2 . After some algebra, we obtain 2 α (r) * n β (r) n e 2i t dt Ωn where the subscript "P" indicates perturbative contributions and the ellipsis represents the higher order terms in E/m 2 a 2 and B/m 2 a 2 . We have dropped terms which do not survive after integration over Π z , i.e., terms odd in Π z for eq. (3.25) and even in Π z for eq. (3.26).

JHEP02(2020)055
See appendix D for the derivation. Inserting eqs. (3.25) and (3.26) into eq. (3.23), we get with u = Π z /ma. The summation over the first term in the square brackets contains a divergence as we will see. Let us make use of the zeta function regularization. We can cope with this summation with the help of the Hurwitz zeta function ζ(p, q): where to evaluate the two terms in eq. (3.27) we will be interested in p = 1 and p = 3.
Here, in the first equation, we have utilized the fact that higher Landau levels (n ≥ 1) have two modes (r = 1, 2) for each n while the lowest Landau level (n = 0) has only one mode. The last "-1" stems from this mismatch. It is known that the Hurwitz zeta function can be extended to a complex value of p via the analytic continuation, which contains a pole at p = 1. By expanding it around p = 1 and also Taylor-expanding in b, we get lim p→1 n,r 1 (bn + 1) p = 2 where the ellipses involve higher order terms in B/m 2 and terms which vanish at p = 1. We also need the behavior at p = 3, which is given by n,r 1 (bn + 1) 3 = 2 where the ellipses imply higher order terms in B/m 2 . Using the relations (3.29) and (3.30), we finally obtain the following form of the current which contains the divergence and finite terms at the leading order in E/m 2 a 2 and B/m 2 a 2 : Recalling that the level, n, represents the transverse momentum squared and the summation is asymptotically N 1/n ∼ ln N for N 1, we replace the divergence at p = 1 with lnΛ 2 /m 2 withΛ being a UV-cutoff. Now we are ready to discuss the physical origin of eq. (3.31). For this purpose, it is convenient to see how this current enters the equation of motion for the gauge field once JHEP02(2020)055 we make the gauge field dynamical. Since our expression for the current eq. (3.31) was based on the vector potential (2.4), the Maxwell equation reduces to 0 =ĖΛ + gΛQ J z vac → 0 Here we have explicitly indicated bare quantities which should be renormalized with the subscriptΛ.
One can easily see that the first parenthesis is nothing but the running gauge coupling evaluated at the fermion mass m and hence it is finite. Moreover, recalling the Ward-Takahashi identity, one finds gΛEΛ = g m E m and gΛBΛ = g m B m where quantities with a subscript m are the renormalized ones evaluated at the fermion mass scale. Noticing that the electromagnetic field appears the combinations gΛEΛ or gΛBΛ in eq. (3.32), we may rewrite the equation in terms of the renormalized quantities 5 where we have suppressed the subscript m for notational brevity. The final task is to identify the origin of the terms suppressed by the power of the fermion mass in eq. (3.34). Such higher dimensional operators arise in the low energy effective action when we integrate out a heavy fermion. Its concrete form is known as the Euler-Heisenberg action. On top of the kinetic term for the gauge field, this yields the following corrections [1,41] where the ellipses imply terms higher powers in E/m 2 a 2 and B/m 2 a 2 . Inserting the specific gauge field configuration (2.4) into eq. (3.35), we find the following terms in the equation of motion for the gauge field after variation with respect to A z : This coincides with the second and third term in eq. (3.34), implying that the origin of (3.34) is nothing but the Euler-Heisenberg action. Before closing this section we would like to emphasize that, even if one takes a large enough fermion mass to exponentially suppress the fermion production, there always exist power-law suppressed terms in the expectation value of the current. However their origin is merely higher dimensional operators which should not to be confused with the particle production. This conclusion holds for any particle production in general.

JHEP02(2020)055 4 Applications
Fermion production in strong helical, abelian gauge fields plays an important role in many different areas of particle physics, solid state physics and cosmology. While the relevance of these questions for the former fields of research has been known for a long time (see e.g., the seminal papers [50] and [1,2,11]), the importance of chiral fermion production for cosmological processes has only fairly recently gained significant attention. In the remainder of this section, we will focus on the dual production of gauge fields and fermions in cosmic inflation, triggered by coupling the inflaton to the Chern-Simons term of an Abelian gauge theory. However, the results presented here also have important implications for a much broader class of processes in the early Universe. For example, this coupling to the Chern-Simons term can play a crucial role in the phase of preheating after cosmic inflation [51][52][53][54] and has been employed in [33,[55][56][57][58][59] to generate large scale magnetic fields which can survive in the intergalactic voids until today. The (spontaneously) CP -violating nature of this coupling moreover makes it a promising candidate to implement baryogenesis, either by a direct coupling of the inflaton (or SM Higgs) to the fermion current [60][61][62][63][64], or by sourcing a baryon asymmetry from the generated CP -violating background of gauge fields [27][28][29][30][31][32][33]. At much lower energy scales, the same coupling to the Chern Simons term can play a crucial role in explaining the smallness of the electroweak scale [34][35][36][37][38].
In this section, we will apply the results derived in the previous sections to analyze particle production in quasi de-Sitter space, i.e., during cosmic inflation. The analysis of gauge field production in this context, sourced by a coupling of the particle driving cosmic inflation (the inflaton) to the Chern-Simons density of an abelian gauge group, dates back to ref. [65]. The coupling to a chiral current of (uncharged) fermion current was first studied in [66] and recently refined in [67]. The simultaneous coupling to both abelian gauge fields and massless charged fermions, as dictated by the chiral anomaly equation, was first studied in [9]. Here we extend these results to include a finite fermion mass. As we will see, the Chern-Simons coupling induces remarkable signatures in the predicted scalar and tensor spectrum at small (sub-CMB) length scales. The presence of fermions, and in particular the presence of a fermion mass term, drastically alters these predictions.

Upper bound on gauge field production during cosmic inflation
We consider a pseudo-scalar singlet φ (referred to as 'inflaton' or 'axion' in the following), coupled to the Chern-Simons term and to the chiral fermion current, 1) see appendix A for details on our notation and conventions. 6 We will leave the scalar potential V (φ) unspecified for now, assuming only that it is flat enough to support slow- 6 Note in particular that here we have taken the φ-dependent? complex phase in the Dirac mass term to be zero. This corresponds to a particular choice of the couplings transporting the spontaneous CP -violation. We leave an investigation of the fully general model parameter space for future work.

JHEP02(2020)055
roll inflation,φ/H 1,φ/(φH) 1, where H =ȧ/a denotes the Hubble parameter. 7 The resulting equation of motion for the fermion is given in eq. (2.1). The equations of motion for the scalar field and the gauge fields read Here as usual, we have decomposed the inflaton field into a homogeneous background with (small) perturbations, φ(t, x) → φ(t) + δφ(t, x). The backreaction of the gauge fields on the inflaton background is encoded in the spatial average FF in eq. (4.2), whereas the fermion backreaction enters through J ν ψ = ψγ µ ψ. Decomposing the vector potential into longitudinal and transverse modes, A = A L + A T with 0 = ∇ · A T , one can readily see thatφ never affects A L but leads to an instability for A T as discussed below.
Without fermion backreaction. Neglecting for a moment the fermion current J ν ψ , eq. (4.3) can be expressed in Fourier space as where ξ parameterizes the inflaton velocity, and we have introduced the helicity decomposition for the transverse mode with the polarization tensors obeying k· (σ) (k) = 0, (σ) (k) * · (σ ) (k) = δ σσ , and k× (σ) = −σik (σ) (k). Assumingξ/H 1, consistent with the slow-roll approximation, eq. (4.4) is solved by the Whittaker functions W k,m (z). This leads in particular to one exponentially growing mode (in accordance with the negative effective squared mass in eq. (4.4)). After imposing Bunch-Davies initial conditions in the far past, the solution for this mode reads 7 See ref. [68] for a detailed discussions on the role of the shape of the scalar potential.

JHEP02(2020)055
From this, we can compute the expectation values of the electric and magnetic fields as with κ = −kη = k/(aH) and where we have introduced physical electric and magnetic fields defined byÊ = E/a 2 andB = B/a 2 , respectively. See e.g., [32] for a detailed derivation and discussion of these results. Note that in accordance with statistical isotropy, we expect E = 0 = B when averaging over the entire Universe. However, well within any given Hubble patch, the electric and magnetic fields are homogeneous, (anti-)parallel and thus locally select a preferred direction. On average, the magnitude of these fields in given by eqs. (4.8) and (4.9), respectively. Under the assumption that both ξ and H are approximately constant, the magnitudes of these physical fields are constant. Hence, well within a given Hubble patch, the gauge field configuration sourced by eq. (4.1) corresponds precisely to the constant, helical external gauge field configuration studied in sections 2 and 3. We note that the choice of sign ofφ distinguishes between parallel and anti-parallel E and B fields, with λ > 0 corresponding to B > 0 in sections 2 and 3.
The induced current. The analysis above neglected the fermion current. As discussed in section 3, the presence of strong electric and magnetic fields leads to the production and acceleration of charged fermions, resulting in the induced current J ind . In sections 2 and 3 we considered these gauge fields to be classical, external fields. In the context of axion inflation these gauge fields are sourced dynamically as discussed above and hence the generation and acceleration of charged fermions drains energy from these fields. We thus expect the amplitude of the gauge fields to be reduced compared to the expressions (4.8)-(4.10). See ref. [9] for an analogous discussion for the case of massless fermions. In section 3.3 we computed the different contributions to the induced current assuming homogeneous, constant (anti-)parallel electric and magnetic fields and a static Universe. From the discussion above, we note that approximately homogeneous and constant physical (anti-)parallel electric and magnetic fields are sourced during axion inflation. We now turn on the cosmic expansion adiabatically and assume de Sitter spacetime. Note that this adiabatic approximation only holds if all relevant microphysical processes are much faster than the cosmic expansion, see discussion below. By replacing the time with the conformal JHEP02(2020)055 time η, one may rewrite eq. (3.24) as g|Q|Ê , (4.11) Assuming a constant physical electromagnetic field, we can perform the conformal time integral, which gives where we have defined the physical current byĴ = J/a 3 . Taking the mass to zero in this expression, we can recover the known result for the chiral fermions in refs. [9,39]. One may also take B → 0 instead. Then we recover the known result for the Schwinger effect without the magnetic field [69,70]. This indicates that our understanding smoothly connects massive/massless regimes with/without magnetic field. Eq. (4.12) is a key result for the study of cosmological applications in the remainder of this section. While eq. (3.24) is an exact result, the generalization (4.12) to de Sitter space requires some additional assumptions. Let us pause here for a moment to clarify these.
1. As pointed out in ref. [43], the effect of cosmic expansion becomes relevant for weak electric fields, when the acceleration induced by the electric field is comparable or even smaller than the Hubble expansion. The latter case can even lead to an effective motion of the fermions antiparallel to the electric field lines. The validity of eq. (4.12) thus requires gQ|Ê|/(mH) 1. From figure 6, we see immediately that this condition is fulfilled in the entire parameter space where the induced current is relevant.

As we have seen in section 3.3, the induced current contains the Euler-Heisenberg
terms which we ignored in the discussion above. For (mildly) strong field case g|Q|Ê m 2 , the leading Euler-Heisenberg term sets the renormalization scale for the coupling of O(g|Q|Ê), while the real part of the rest is suppressed by m 2 /g|Q|Ê (see [41] as a review). The latter is negligible compared to the particle production as long as g|Q|Ê/(mH) 1. 8 It by chance coincides with the previous condition and is satisfied for the entire parameter space of our interest. Thus the Euler-Heisenberg terms are important only for g|Q|Ê m 2 , which is out of our interest since the particle production and hence the fermion backreaction are anyway negligible. Here we implicitly assumed thatÊ ∼ |B|. It is valid in the case of our interest with of order unity ξ or ξ eff . Note that it is ξ eff that determines the hierarchy betweenÊ and |B| for the equilibrium solution.
3. In the discussion above, we neglected cosmic expansion when computing the fermion production rate. This is valid as long as the relevant microphysical length scale is much smaller than the Hubble horizon, i.e., m T H. From eq. (3.9) we note JHEP02(2020)055 that particle production is most efficient for m 2 T ∼ g|Q|E/π, and hence the flat space expressions for the production rate can be applied if E H 2 , which again holds whenever the induced current is relevant. Note that for massless fermions, this condition becomes equivalent to requiring that the particle production is fast compared to the cosmic expansion rate (see ref. [9]), whereas the additional mass scale of massive fermions relaxes this condition.
The equation of motion (4.3) implies the energy conservation equatioṅ where ρ A = 1 2 (Ê 2 +B 2 ) denotes the physical energy density stored in the gauge fields. To obtain this expression, we have assumed that the average of Ê ·Ĵ is approximated withÊĴ z ind . The second term on the right-hand side describes the increase of the gauge fields sourced by the inflaton motion whereas the third term describes the transfer of gauge field energy into the fermion sector. Assuming dynamical equilibrium between the inflaton, gauge field and fermion sector,ρ A = 0, and inserting the explicit expression for the induced current, eq. (4.12), this leads to an algebraic consistency equation forÊ andB, − 2H(Ê 2 +B 2 ) + 2ξ eff HÊ|B| = 0 , (4.14) with Neglecting the induced current (i.e the second term in eq. (4.15)), eq. (4.14) has two solutions (for any given ξ), which can be depicted as two straight lines in theÊ vsB plane. The analytical solution given by eqs. (4.8), (4.9) corresponds to a particular point on one of these lines. Switching on the induced current, eq. (4.14) describes a closed contour in theÊ vsB plane. Since the equation of motion is highly non-linear, we are no longer able to solve it analytically. However, we can place an upper bound onÊ,B orÊB by extremizing these quantities over the closed contour described by eq. (4.14). 9 This is depicted in the left panel of figure 6 for different values of the fermion mass. For simplicity, we have neglected the running of gauge coupling (see eq. (3.33)), simply setting Q = 1 and g = 1/ ξ breaks down, and the analysis we perform here is no longer valid. When treatingφ and hence ξ as dynamical parameters by solving eq. (4.2), this constraint is automatically fulfilled since the backreaction of the gauge fields on inflaton equation of motion limits the growth ofφ. In practice, this implies that the large values of ξ required for a violation of eq. (4.16) are not reached dynamically, see also section 4.2. Finally, the gray shaded region indicates the regime where thermalization of the produced fermions can no longer be neglected according to the estimate (4.20). Here we only show the condition for the massless case because the bound becomes the most stringent as can be readily seen from eq. (4.20). We also check that the bound relaxes only by an O(1) magnitude for the heaviest mass parameter taken in figure 6. In summary, as expected, increasing the fermion mass reduces the fermion production and hence the induced current, leading to larger values of the electric and magnetic fields.
The upper bounds depicted in the left panel of figure 6 are clearly conservative and here is no compelling reason why they should be saturated. We can proceed by solving eq. (4.14) under an additional assumption: if all the equilibration processes between the inflaton, gauge field and fermion sector are sufficiently fast, we expectÊ andB to be given by eqs. (4.8) to (4.10) with the replacement ξ → ξ eff . This set of equations can be solved selfconsistently, and we depict the results for this 'equilibrium' solution in the right panel of figure 6. Comparing figure 6 with the conditions listed below eq. (4.12), we see that also for the equilibrium solution, these are always fulfilled in the regime of interest.
We emphasize that both panels of figure 6 only provide an estimate for the magnitude of the generated gauge fields. An exact solution of the non-linear equation of motion (4.3) would be highly desirable for a more precise study of the phenomenology of these models, but is unfortunately beyond the scope of this paper. Nevertheless, we very clearly see that the presence of fermions, and in particular their mass, changes the prediction for the gauge field production by several orders of magnitude.

JHEP02(2020)055
Thermalization and Pauli blocking. The key quantity in the study of backreaction is the induced current. So far, we have neglected the interaction of fermions in the estimation of induced current. Here we discuss the condition justifying this approximation. As one may see from eqs. (3.24) and (4.12), we have assumed that the acceleration by the electric field continues for 1/H. However, if fermion scatterings are efficient, this approximation is no longer justified.
Let us start our discussion with the "would-be" temperature if the particles were thermalized: T wb 30 π 2 g * n ψ ω 1 4 , ω = g|Q|Ê where ω denotes the energy of the most relevant Landau level after experiencing acceleration in the electric field for a time period min[τ s , H −1 ] with τ s being a typical time scale of scatterings. If we suppose that the plasma were thermalized, the typical interaction rate of particles in the thermal plasma is In most cases of our interest, we find that the typical energy of fermions after production is much larger than the would-be temperature, ω T wb . This leads to an additional suppression of the scattering rate of the generated fermions compared to τ th : For instance, the concrete form including this suppression factor is α 2T wb (T wb /ω) 1/2 [71][72][73][74] if the generated particles are also charged under non-Abelian gauge field as in the SM [9,75,76]. Here, to avoid complications and model-dependent discussions, we just conservatively require 1 < τ th H as a condition to neglect the interactions of the generated fermions. One can readily see that, for 1 < τ th H, even an initially thermal plasma would drop out of thermal equilibrium, indicating how conservative this requirement is. The condition, 1 < τ th H, gives 10 4g|Q|Ê . (4.20) Here we have assumed m 2 g|Q|Ê and g|Q|Ê > H 2 , since otherwise the backreaction from the induced current is anyway negligible. For the parameter choices in figure 6, this condition is always fulfilled. However, we note that for a smaller energy scale of inflation, the upper bound onÊB/H 4 from eq. (4.16) becomes less stringent, potentially opening up some parameter space where thermalization is relevant, in particular for a heavier mass.
It is sometimes argued that the production of fermions is not efficient due to Pauli blocking in the final state. In our computations in section 3 this is already intrinsically JHEP02(2020)055 accounted for, and is reflected e.g., in the condition |α| 2 + |β| 2 = 1. Nevertheless, we find significant particle production. Here we interpret these results in terms of simple quantum mechanical arguments.
Immediately at their generation, the fermions are characterized by a vanishing velocity in z-direction and a finite (for n > 0) velocity in the transverse direction, encoded in the transverse mass m T . From the spacing of the energy levels in eq. (2.22), we see that the quantum mechanical uncertainty in position space is given by dxdy ∼ 2π/(gsQB). With this, we find for the production rate within a box of dimensions dx dy L, where we have employed the asymptotic solution 3.9 to obtainṅ ψ . Consequently the acceleration in z-direction accumulated before the creation of the subsequent fermion is given by Comparing this to the spacing of the p z levels, ∆p z = 2π/L, we conclude that Pauli block is inefficient if which is clearly intrinsically satisfied. In particular, we note that for massless fermions in the lowest Landau level, this bound is saturated (as expected for fermion production from a gap-less dispersion relation), whereas massive fermions and higher Landau levels feature an occupation number significantly below the Pauli limit. This extends and improves the discussion given on this matter in ref. [9].

Axion inflation
In the previous section we discussed the gauge field production induced by the nonvanishing velocity of a pseudo-scalar field. This is a key ingredient in axion inflation (see e.g., [77] for a review). In this context, the gauge field production has shown to lead to striking and observable signatures both in the scalar and tensor power spectrum [68,77,[77][78][79][79][80][81][82][83][84][85][86]. In this section, we revisit these predictions accounting for the presence of massive fermions. For the purpose of this discussion, we will choose a linear scalar potential in eq. (4.2), motivated by models of axion monodromy [87] with the parameter µ determined by the observed amplitude of the scalar power spectrum at CMB scales. For the resulting phenomenology it will be mainly relevant that this leads to an inflaton velocity which grows monotonously during inflation, for a discussion of the impact of the choice of scalar potential see [68]. In the figures below, we will further choose α/(πf a ) = (0.02 M P ) −1 , which is the maximal value consistent with current constraints on scalar non-gaussianities in the CMB. With these choices, we can numerically solve the inflaton equation of motion (4.2) inserting the results forÊB(ξ) depicted in figure 6. Figure 7 illustrates how the fermion masses impact the evolution of inflaton velocity. For massless fermions (red curve) we reproduce the result of ref. [9]: due to the efficient generation of fermions, the gauge field abundance and hence the backreaction on the inflaton field are strongly suppressed, leading to ξ ∝φ/H ∝ N −1/2 . This is the usual result for inflation in monomial potentials without a gauge field background. As the fermion mass increases, the energy conversion from gauge fields into fermions becomes less efficient. This results in a stronger gauge field background and hence in a stronger backreaction on the inflaton field, reducing the inflaton velocity towards the end of inflation. Defining the end of inflation as (φ) 2 /(2H 2 M 2 P ) = 1, this additional friction prolongs the duration of inflation, explaining why the x-axis in figure 7 extends to negative values (here N = 0 is defined as the end of inflation in the absence of gauge fields). This is particular visible for m/H = 100 (purple curve) in the left panel of figure 7, which essentially reproduces the result of axion inflation in the absence of fermions, see e.g., ref. [77]. In the right panel of figure 7 we show the corresponding evolution for ξ eff , defined in eq. (4.15). In this case, the induced current is included in the effective parameter ξ eff . Lighter fermions correspond to a smaller value of ξ eff and hence a less efficient gauge field production. In summary, light fermions lead to a strong induced current which makes gauge field production more difficult, whereas sufficiently heavy fermions effectively decouple. Scalar power spectrum. In slow-roll approximation, the equation of motion for the inhomogeneous perturbations δφ of the inflaton field is given by [81] δφ where the first term is the standard vacuum contribution and the second term in sourced by the Chern-Simons charge. We note that both the homogeneous equation of motion for the inflaton as well as the scalar power spectrum for its fluctuations are sensitive to the gauge field and fermion production only through the magnitude ofÊB, which in turn is a function of ξ ∝φ. Consequently we can convert the estimates depicted in figure 6 into estimates of the scalar power spectrum, taking into account both the backreaction of the gauge fields in eq. (4.2) as well as the additional source term in eq. (4.26). The resulting scalar power spectrum for different values of the fermion mass is shown in figure 8. At CMB scales (N 55), the parameter ξ is small and both gauge field and fermion production are inefficient. At these scales, the model thus closely resembles a standard single-field slow roll inflation model (indicated by the dashed black line). At smaller scales, towards the end of inflation, ξ increases and the scalar power spectrum is strongly enhanced. Lighter fermion masses lead to a smaller scaler power spectrum, since the induced current inhibits the gauge field production. The slightly earlier rise in the scalar power spectrum for heavy fermions is due to the definition of N = 0 asφ 2 /(2H 2 M 2 P ) = 1, which occurs a few e-folds later when efficient gauge field production induced a strong friction in the inflaton equation of motion which delays the end of inflation.
For strong gauge fields, β 1, the scalar power spectrum can be approximated as ∆ 2 s 1/(2πξ) 2 . An upper bound on the gauge fields implies an upper bound on the amount of gauge friction and thus translates to a lower bound on ξ. In this sense, the 11 Here we are assuming that the estimate of the variance ofÊB is given by the expression found in [81] in the absence of fermions. In other words, we are assuming that the fermions only change the magnitude of the generated gauge fields but not their spectral shape. It will be crucial to scrutinize this assumption using dedicated numerical simulations in the future.

JHEP02(2020)055
curves in the left panel of figure 8 can be viewed as upper bounds for the scalar power spectrum for different fermion masses. The equilibrium solution predicts overall smaller gauge fields compared to the upper bound, resulting in a smaller scalar power spectrum, in particular for light fermions.
The horizontal lines in figure 8 indicate the CMB normalization and the primordial black hole bound [81], respectively. Above the latter, the large scalar perturbations lead to a too large probability of producing primordial black holes, in contradiction with observations. This can be remedied by considering couplings to multiple abelian gauge groups [68]. However, given the significant uncertainties in our estimates ofÊ andB as well as the limits of perturbation theory for large values of ξ [88,89], we shall for the moment ignore this problem. The main purpose of this section is not to derive a precise upper bound on the fermion mass which would circumvent the primordial black hole bound, but rather to highlight the significant impact of (massive) fermions on the scalar power spectrum of axion inflation.
Gravitational wave spectrum. The computation of the tensor power spectrum requires the knowledge of the anisotropic contribution to the energy momentum tensor, sourced by vacuum fluctuations as well as by the gauge field and fermion sector. The contribution from the vacuum perturbation is given by the usual expression, Ω vac GW = Ω r /12(H/(πM P )) 2 with Ω r = 9.1 · 10 −5 denoting the fraction of radiation energy density today. In the absence of fermions, the gauge field contribution can be estimated as (see e.g. [77,79]) The contribution from uncharged chiral fermions, in the absence of gauge fields, was computed in [90] and found small compared to the vacuum contribution. A full computation of the gravitational wave spectrum resulting from the interplay of the fermion and gauge field sector requires the solution of the inhomogeneous non-linear gauge field equations taking into account the fermion backreaction, which is beyond the scope of this paper. However, for the equilibrium solution (requiring a self-consistent solution of eq. (4.8) and (4.9) with ξ → ξ eff ), we can perform an estimate for the gravitational wave spectrum by inserting the resulting ξ eff into eq. (4.27), adding the vacuum contribution but neglecting the fermion contribution.
The resulting spectra are shown in figure 9. As expected, the presence of light fermions reduces the gauge field abundance and hence the expected gravitational wave signal. This is a striking example of the ability of gravitational wave detectors to probe the properties of particle physics models.
Finally, we note that the equilibrium solution assumes a fast backreaction of the generated fermions on the gauge field abundance. If this backreaction is slower than, e.g., the change in ξ, we would expect the gauge fields to 'overshoot' their equilibrium value. Moreover, the gravitational wave spectrum depicted in figure 9 neglects any contribution from the fermion sector. Given that we found the fermion production to be fast compared to the Hubble expansion rate, we expect a larger fermion abundance than found from the gravitational production in ref. [90], and hence the resulting contribution to the gravita- tional wave spectrum may be relevant. In view of this, the results depicted in figure 9 should be viewed as lower bounds on the gravitational wave production in axion inflation. We leave a more detailed investigation to future work.

Baryogenesis from axion inflation
As pointed out in refs. [31][32][33], axion inflation with the abelian gauge group identified as the SM hypercharge U(1) Y can provide the initial conditions required for successful baryogenesis from decaying helical hypermagnetic fields [27][28][29][30]. If the SM Higgs expectation value is stabilized at zero during inflation, axion inflation leads to a dual production of helical hyper gauge fields and massless chiral fermions. The latter thermalize once Yukawa interactions and sphaleron processes become relevant in the thermal plasma after inflation. However, under certain conditions, depending also on the chiral charge in the fermion sector, the helical hypermagnetic fields can survive until the electroweak phase transition. During the electroweak phase transition the hypermagnetic fields are converted into electromagnetic fields, thereby generating a baryon asymmetry which can explain the matter-antimatter asymmetry observed in our Universe [33].
If the SM Higgs expectation value is stabilized at non-zero vev during inflation, axion inflation will lead to a dual production of electromagnetic gauge fields and massive chiral fermions. Depending on the details of the Higgs potential and the properties of thermal plasma, the electroweak symmetry may be restored after inflation, implying a conversion of the electromagnetic fields to hyper gauge fields. The magnitude of these hyper gauge fields, together with chiral charge stored in the fermion sector at this point in time, determines the efficiency of the subsequent baryogenesis process. We leave a detailed study of this (model-dependent) symmetry restoration process and the subsequently generated baryon asymmetry to future work.

Conclusions
The production and transport of charged particles in the presence of electromagnetic fields are ubiquitous in nature. In this paper, we study the production of charged massive Dirac fermions in parallel electric and magnetic fields. By solving the equation of motion for fermions in a background gauge field, we compute various quantities semi-analytically in terms of Bogoliubov coefficients. We analytically show how the chiral anomaly arises from the equation of motion [eqs. (3.13), (3.14), and (3.15)], and also numerically evaluate the chiral charge and the pseudo-scalar operator yielding the correct chiral anomaly equation [eqs. (3.17) and (3.19)], which guarantees a smooth connection between the case of massless and massive fermions. In addition, we evaluate the induced current and explicitly show how three contributions appear: (i) the current from fermion production [eq. (3.24)], (ii) the running of the gauge coupling constant [eq. (3.33)], and (iii) the Euler-Heisenberg terms in a 1/m-expansion [eq. (3.31)]. This result clarifies that, while particle production is exponentially suppressed in the limit of heavy fermions, the expectation value of any operator contains terms suppressed in powers of 1/m, which is nothing but higher dimensional operators from integrating out heavy particles and should not be confused with the particle production. We moreover demonstrate explicitly how the Pauli exclusion principle manifests itself, showing that the occupation number per phase space volume is saturated only for massless fermions in the lowest Landau level.
Equipped with these results, we discuss a pseudo-scalar inflaton which couples to a U(1) gauge theory solely via the Chern-Simons coupling. The velocity of the inflaton increases over the course of inflation, which triggers the efficient production of helical gauge fields. We demonstrate how the production of charged fermions backreacts on the gauge field equation of motion and suppresses the production of helical gauge fields. As expected, this backreaction becomes more significant for lighter fermions, resulting in a suppression of the sourced scalar/tensor perturbations. We find that depending on the fermion mass, the stringent constraints from primordial black hole overproduction may be avoided, while still implying a gravitational wave spectrum which is well above the standard vacuum contribution. Our discussion neglects the interactions among fermions after production. We clarify the condition for this approximation to be valid (see also [9]), i.e. the condition that ensures that thermalization of the fermion population during inflation is inefficient. To go beyond this limitation, one would need to treat the evolution of the phase space density of the fermions dynamically, as recently discussed in ref. [91] in the context of an f 2 (φ)F F inflation model. Extending their analysis to our case is surely worth pursuing in future.
Although we restrict ourselves to the Chern-Simons coupling as a starting point for cosmological applications, one may also consider more general shift-symmetric couplings such as (∂φ) · J 5 and ψe iγ 5 φ/fa ψ. Up to a field redefinition associated with a chiral rotation, one of these three couplings is redundant and we are left with two independent couplings, i.e. the Chern-Simons coupling and φ-dependent mass coupling. This field redefinition never changes the physical result as it should be, which is however has been explicitly shown only in two limited cases: the case without the gauge coupling (uncharged fermions) in refs. [66,67] and the case without mass in ref. [9]. Including this basis independence, a smooth JHEP02(2020)055 connection between two regimes is therefore desirable for more complete understanding of the system at hand.
For simplicity we assume a constant mass for the Dirac fermions throughout this paper. However, if we would like to identify this gauge group as part of the SM gauge group, we need to recast the mass in terms of the Higgs field. As recently pointed out in ref. [92], once we treat the Higgs field dynamically, the asymmetric production of fermions backreacts onto the Higgs, destabilizing the Higgs effective potential in contrast to the usual thermal mass, which can lead to an interesting signature in the bispectrum of the scalar perturbations generated during inflation. Although ref. [92] focuses on the production via the φ-dependent mass coupling, it would be interesting to extend their analysis to include the Chern-Simons coupling which then involves the production of helical gauge fields, the production of fermions obeying the anomaly equation, and also the Higgs production via the Schwinger effect. We hope that our work will stimulate further studies on this direction, including the impact of these effects on the dynamics of the Higgs field and their observational signatures.

Acknowledgments
It is a pleasure to thank Tomohiro Fujita for helpful discussions related to this project. This work was funded by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306.

A Notation and conventions
We consider a pseudo-scalar singlet φ together with a Dirac fermion ψ of mass m which is charged under an Abelian gauge group. Taking into account all interaction terms (up to dimension 5) which respect the shift-symmetry of φ, the most general action (which does not explicitly break parity) reads with / D = (∂ µ + igQA µ )γ µ , α = g 2 /4π and Q denoting the charge of the fermion ψ under the abelian gauge group with vector potential A µ . Here c A , c 5 and c m are real coupling constants, f a indicates the cut-off scale of the effective theory and g µν denotes the FRW metric. Performing a chiral fermion rotation, we can eliminate (for example) the term proportional to (∂ µ φ)J µ 5 where J µ 5 =ψγ µ γ 5 ψ. In this paper, we will moreover for simplicity drop the φ-dependent phase of the mass term. This leads to the fermion equation of motion (2.1).
In eq. (A.1) we have introduced the comoving quantities ψ, A µ and g µν , related to the corresponding physical quantities (indicated by a hat) as

B Solving the fermion equation of motion
Second order equation of motion. This appendix provides some additional details on the results derived in section 2. We start be deriving the second order equation of motion (2.3) starting from the first order equation (2.1): Here the relative minus sign in the spatial part of the gauge co-variant derivative arises from our convention ∂ µ = (∂ 0 , ∇) and A µ = (A 0 , −A). In the third term we have introduced the anti-symmetric operator In temporal gauge, A 0 = 0, we can immediately evaluate explicitly the relevant components of σ µν and F µν , With this, the third term of eq. (2.1) can be expressed as
For the higher Landau levels, an analogous procedure leads to eqs. (2.23) to (2.26).

JHEP02(2020)055
cancel out between r = 1 and r = 2. Inserting eq. (2.5) and more specifically expanding in terms of the creation and annihilation operators defined at some time t > 0 yields Using the explicit expression for the eigenvectors U 0 and V 0 in eq. (2.32) we find ExpandingB 0 andD 0 in terms of the creation and annihilation operators of the original vacuum at t < 0 yields The factor δ (2) (0) arises since the creation and annihilation operators contain the same wave vector k y and k z , as enforced by the delta function in the second line of eq. (C.7). Plugging these results back into eq. (C.7), we obtain To treat the vacuum contribution unambiguously, we have explicitly inserted the regulator RΛ withΛ denoting a UV-cutoff. By changing the integration variable k z → Π z , one can see that the vacuum contribution proportional to "−1" in the first bracket vanishes because the integrand including the regulator is an odd function of Π z . This confirms that, in the case at hand, the vacuum contribution originating from the η-invariant vanishes. In contrast, the η-invariant plays a crucial role in the case discussed in ref. [12]. In summary, we obtain the following expression for the chiral charge: q 5 (t) = gQB 2π dk z 2π Π z Ω 0 2 |β 0 | 2 + s ma Ω 0 e 2i t dt Ω 0 (t ) α * 0 β 0 + H.c. . Since the equations of motion are the same for r = 1, 2, we omit the index r here and hereafter for notational simplicity. We solve these equations in the weak field expansion in the following.
In the weak field limit, one can clearly separate terms that contribute to the Bogoliubov coefficients into two classes. The first are those that accompany the exponential suppression factor e − πm 2 a 2 g|Q|E , which corresponds to the particle production due to the electric field. We call it a "non-perturbative" contribution since it is not obtained from the weak field expansion. The second are those that do not accompany such a suppression factor, which we call a "perturbative" contribution. It is not understood as particle production, but it nevertheless has effects on, e.g., the induced current of the fermion and hence the equation of motion of the gauge field. Indeed, the gauge coupling renormalization and the Euler-Heisenberg Lagrangian originate from the latter class of terms, as shown in section 3.3. For this reason, we sometimes refer to the latter as a vacuum contribution.
In the following, we focus on the perturbative contribution. We first describe how the weak field expansion proceeds with the help of integration by parts. We then compute the perturbative contribution up to third order in the weak field expansion. The results are used in sections 3.2 and 3.3.
Weak field expansion as integration by parts. The weak field expansion can be systematically performed by integration by parts. Eq. (D.1) is formally solved as where we have performed the integration by parts to get the second line. We ignore the time dependence of the scale factor here and henceforth. Note that E and its derivatives have to vanish at the initial time for the initial state to be defined unambiguously. Since the time derivatives of Ω n and |β n | 2 contain additional E orĖ, the second term in eq. (D.3) is sub-leading compared to the first term for the perturbative contribution. In this sense, the above integration by parts is equivalent to the weak field expansion. By repeating the JHEP02(2020)055 same procedure, we obtain the following formula: α * n β n (t)e 2i t dt Ωn where the subscript "P" stands for "perturbative". The corresponding |β n | 2 P is obtained by integrating eq. (D.2). In the following, we will recursively solve these equations to derive the explicit forms of (α * n β n ) P and |β n | 2 P in the weak field expansion.
Leading order contribution. It is easy to compute the leading order contribution to α * n β n . From eq. (D.4), it is given by α * n β n e 2i t dt Ωn P ism T a gQE Higher order perturbative contribution. Now we compute higher order perturbative terms in the weak field expansion. In particular, we focus on [α * n β n e 2i t dt Ωn ] P and |β n | 2 P that are necessary to derive the Euler-Heisenberg term. We compute terms up to third order in the weak field expansion that correspond to the lowest order terms of the Euler-Heisenberg Lagrangian. If one goes beyond the third order terms, one obtains higher order Euler-Heisenberg terms. Both [α * n β n e 2i t dt Ωn ] P and |β n | 2 P vanishes at the zero-th order since α * n β n e 2i t dt Ωn where we have ignored |β n | 2 P in the second line since it vanishes at the zero-th order. Plugging it into the equation for |β n | 2 , we obtain at the first order gQE Ω 3 n + · · · = m 2 T a 2 (gQ) 2 16Ω 6 n E 2 + · · · . (D.7) Now we move to the higher order terms. The second order terms vanish since the second order term of α * n β n e 2i t dt Ωn P is again purely imaginary. Up to the third order, we n E 2Ė 3m 2 T a 2 − 38Ω 2 n + 248Π 2 z + (Odd in Π z ) + · · · , (D.8) where we have inserted eq. (D.7) in the second line, and ignored terms higher than the third order. We have also ignored terms of O ( ... E) that correspond to higher derivative terms in the effective action for the gauge field. Correspondingly |β n | 2 P is given by gQE Ω 2 n α * n β n e 2i t dt Ωn P = 7 32 m 2 T a 2 (gQ) 3 Ω 10 n E 2Ė Π z + (Even in Π z ) + · · · . (D.9) These expressions are used in section 3.3 to derive the Euler-Heisenberg term. Note that only lower order terms are necessary to obtain terms at a specific order (the third order in the above case), and hence this procedure can be systematically extended to higher order terms.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.