Axion assisted Schwinger effect

We point out an enhancement of the pair production rate of charged fermions in a strong electric field in the presence of time dependent classical axion-like background field, which we call axion assisted Schwinger effect. While the standard Schwinger production rate is proportional to $\exp(-\pi(m^2+p_T^2 )/E)$, with $m$ and $p_T$ denoting the fermion mass and its momentum transverse to the electric field $E$, the axion assisted Schwinger effect can be enhanced at large momenta to $\exp(-\pi m^2/E)$. The origin of this enhancement is a coupling between the fermion spin and its momentum, induced by the axion velocity. As a non-trivial validation of our result, we show its invariance under field redefinitions associated with a chiral rotation and successfully reproduce the chiral anomaly equation in the presence of helical electric and magnetic fields. We comment on implications of this result for axion cosmology, focussing on axion inflation and axion dark matter detection.


Introduction
Schwinger production [1,2] describes the quantum mechanical, non-perturbative production of pairs of particles and antiparticles in a strong electric field. The production rate is exponentially suppressed by the mass gap in the dispersion relation, particle production ∼ exp − π m 2 + p 2 T g |Q|E , (1.1) for a particle of mass m and charge Q, traveling with a momentum p T transverse to an electric field with magnitude E and gauge coupling g . Not surprisingly, the spontaneous production of particles with large momenta p T is exponentially suppressed.
In this paper, we study the production of charged particles in a strong electric field in the presence of a time-dependent homogeneous background pseudoscalar field φ, with ∇φ = 0 andφ = 0. In analogy to the QCD axion [3][4][5][6], we will refer to φ as axion-like particle, or axion for short. In the presence of such an axion background field, a chiral current coupling (or equivalently a pseudoscalar coupling) induces an additional spin precession operator in the non-relativistic Hamiltonian of a fermion ψ [7]: 1 where the gauge invariant momentum for a fermion is defined as Π = p−gQA with A denoting the vector potential, and σ denoting the Pauli matrices. This spin precession leads to an additional contribution to the fermion energy, with a sign depending on the sign of the axion velocity. For one of the fermion modes, it reduces the mass gap in the dispersion relation, cancelling the contribution from the transverse momentum for suitable values ofφ. Correspondingly, the particle production rate is exponentially enhanced for p T = 0, 2 maximal particle production ∼ exp − πm 2 Somewhat counter-intuitively, for sufficiently large axion velocities, this axion assisted Schwinger effect is thus not exponentially suppressed at large momenta, implying a fundamentally different resulting fermion spectrum. This is the main result of the present paper.
To obtain this result, we solve the Dirac equation in the presence of an electric background field, including the most general shift-symmetric dimension five couplings of an axion-like field to a fermion and gauge field that do not explicitly break C P . The axion velocity and the electric field imply a timedependent background for the fermions, leading to the particle production described by time-dependent Bogoliubov coefficients. We verify the intuition behind key aspects of our computation using the nonrelativistic effective field theory of fermions. Our results are intrinsically invariant under (chiral) fermion rotations, reflecting the basis invariance of all observables.
The non-vanishing axion velocity moreover leads to a magnetic field parallel to the electric field, i.e., a helical gauge field configuration, if the axion couples to the electromagnetic fields. This is, in particular, relevant if we identify the axion φ as the inflaton particle driving cosmic inflation [10][11][12]. Motivated by this, we also study the axion assisted Schwinger effect in a background of parallel electric and magnetic fields. In this case, the fermion dispersion relation is quantized in terms of the Landau levels [13,14], each corresponding to a particular transverse momentum. We confirm that the enhancement mechanism applies also in this case. The helical gauge field background sources a chiral asymmetry as indicated by the chiral anomaly [15,16]. As a non-trivial consistency check, we reproduce the chiral anomaly equation from our result.
Our analysis builds on earlier studies of particle production in (helical) electromagnetic fields, reproducing the chiral anomaly equation and the Schwinger production rate [17][18][19][20][21]. This paper extends these analyses by including a dynamical axion field with general couplings to the fermions and gauge fields, including in particular both operators in Eq. (1.2). The corresponding results of Refs. [17][18][19][20][21] are obtained as a particular limit of the more general results presented here. Moreover our analysis shares some similarities with the so-called 'dynamically assisted Schwinger mechanism' [22,23], which describes an enhancement of the Schwinger pair production rate in electric fields with a non-trivial time dependence.
The remainder of this paper is organized as follows. In Sec. 2, we compute the particle production due to the axion assisted Schwinger effect in the presence of an electric field. We generalize this analysis in Sec. 3 to include an (anti-)parallel magnetic field, in particular verifying the chiral anomaly equation and providing an estimate for the induced fermion current. We comment on implications for axion cosmology, in particular on axion inflation and on axion dark matter detection in Sec. 4 before concluding in Sec. 5.
We have relegated several technical but important details to six appendices. App. A fixes our notation and conventions. Apps. B and C provide the details on deriving the equations of motion for the Bogoliubov coefficients used in Secs. 2 and 3, respectively. App. D is dedicated to our analytical result for particle production due to the axion assisted Schwinger effect. Our results are interpreted in the language of nonrelativistic effective field theory in App. E. Finally, we briefly review the WKB estimate of Schwinger particle production in App. F.

Axion assisted Schwinger effect in an electric field 2.1 Dirac equation with classical background fields
Let us consider an axion-like particle φ, coupled through shift-symmetric dimension five operators to an Abelian gauge boson A µ and a Dirac fermion ψ. Our goal is to study non-perturbative fermion production by explicitly solving the Dirac equation, see App. A for more details on our notation and conventions.
Without loss of generality, we choose the direction of the electric field along the z-axis, i.e., we work in a classical background field configuration described by where the electric field is given by we can write the equation of motion (2.1) in Fourier space as We now decompose the fermion modes as with u λ (v λ ) denoting the (anti-)particle eigenfunctions of the fermionic part of the Hamiltonian for constant A z and φ, respectively, see App. B for details. Inserting this into Eq. (2.4) yields the dispersion relation with the effective mass m T = p 2 x + p 2 y + m 2 now including the transverse momentum. Note that, although the dispersion relation (2.6) is independent of θ 5 and θ m , the eigenvectors u λ , v λ depend on θ 5+m ≡ θ 5 +θ m . The appearance of only this linear combination can be traced back to the symmetry of our setup under chiral fermion rotations. This in particular implies that in the limit m → 0 any dependence on θ 5 and θ m must drop out, since it can be removed from Eq. (2.1) by a suitable fermion rotation.
Bogoliubov coefficients. We now turn on the time-dependence of A z and φ, thus promoting the coefficients α λ and β λ to time dependent Bogoliubov coefficients. Then Eq. (2.4) requires that (2.7) 3 In an FRW background with scale factor a(t ), the conformal mass m is related to the physical mass m ph as m = m ph a.
By exploiting the orthogonality among the eigenvectors, we obtaiṅ Using the relations for the scalar products of the eigenvectors u λ and v λ and their time derivatives given in App. B, we obtain the equations of motion for the Bogoliubov coefficients, where we denote Θ = dt Ω. Here and henceforth we ignore the time dependence of the scale factor, which is justified as long as particle production rate, driven by the time dependence of A z and φ, is much larger than the Hubble expansion rate.
Quantization and initial conditions. In terms of these Bogoliubov coefficients, the quantized fermionic mode function is given by with b λ , b † λ , d λ and d † λ denoting the usual fermionic annihilation and creation operators. Here the minus signs arise from the observation that ± −β * 1 , β * 2 , α * 1 , −α * 2 T satisfies the same equation of motion as α 1 , α 2 , β 1 , β 2 T , and the superscript indicates the initial conditions for α and β, that is, at the initial time t = −∞. In other words, there are two types of particles that are the remnants of the two helicities, and we take both particles absent for the initial state. In our final results we need to sum over both particles since they both contribute to physical quantities. Note that, due to the degeneracy in the eigenvalue Ω, the labelling of the states λ = {1, 2} is arbitrary, and any two orthogonal linear combinations would lead to the same final results. In particular, (non-perturbative) particle production corresponds to a non-zero value of λ,λ |β (λ ) λ | 2 arising from the time-dependence in the background A z and φ. The quantized fermionic mode function (2.11) can be equally expressed in terms of the fermionic annihilation and creation operators B λ , B † λ , D λ and D † λ at any given time t , defined as (2.14) They satisfy the standard anti-commutation relations, where here we have made the dependence on the momentum explicit. These relations indicate that the Bogoliubov coefficients satisfy independently of the time t . We show that these relations follow from Eqs. (2.10) and (2.12) in App. D.

Particle production
For concreteness, let us consider thatφ andȦ z are turned on and off adiabatically at t min and t max respectively, i.e.,θ Hereθ and E denote the constant amplitude of these functions for t min t t max and T denotes the characteristic time-scale for switching theφ andȦ z on and off. We can now proceed to solve Eq. (2.10) numerically. The result, depicted in Fig. 1, displays a remarkable property: The exponential suppression of the non-perturbative Schwinger particle production is not given by exp(−πm 2 T /(g |Q|E )), as one would expect from the gap in the dispersion relation (2.6), but instead the suppression is governed by the bare mass m, This corresponds to an exponential enhancement of the Schwinger production rate for p x,y = 0, possible only in the presence of a suitableθ = 0, i.e., in the presence of a moving axion background field. We call this exponential enhancement of the particle production the axion assisted Schwinger effect. The remainder of this section is dedicated to explaining this result.
Numerical results. We first explain our numerical results displayed in Fig. 1 in more detail. In the upper panels, we plot the spectrum of the produced particles for several different values ofθ with two sets of m, p x and p y . We take t min = 0, t max = 100/ g |Q| E and T = 50/(g |Q| E ), and evaluated the spectrum at t = 150/ g |Q| E . The resulting spectrum features an approximate plateau for −τ < p z /gQE < 0 with τ = t max −t min . The width of the plateau thus depends linearly on the duration of the non-zero electric field τ which implies that for large values of τ, any transient effects at t min, max become irrelevant. The height and the width of the plateau in the caseθ = 0 can be intuitively understood as follows. When the electric field is imposed, the particle is accelerated and Π z increases with time. The gap between the particle and antiparticle energy levels is minimized when Π z = 0, and hence particle production is most efficient at this point. Therefore, only the modes which cross Π z = 0, i.e. −τ < p z /gQE < 0, are efficiently enhanced, corresponding to the plateau of the spectrum. The gap at Π z = 0 is given by m T = p 2 T + m 2 , and hence the height of the plateau is suppressed by exp(−πm 2 T /g |Q| E ). 4 As one can see in the upper panels of Fig. 1, once we turn onθ 5+m , the height of the plateau depends strongly on this parameter. This dependence is shown explicitly in the respective lower panels, where we plot the occupation number for p z / g |Q| E = −50, corresponding to the center of the plateau, as a func-tion ofθ. Two important features stand out. First, the particle production is drastically enhanced asθ increases, and the envelope of the suppression factor asymptotically approaches to exp(−πm 2 /g |Q| E ). In other words, the part of the gap in the dispersion relation due to p T is overcome and the spectrum is correspondingly enhanced by exp(πp 2 T /g |Q| E ). Second, on top of the exponential enhancement, the height of the plateau oscillates withθ. As we explain below, we may interpret this oscillation as an interference effect of two saddle points. We note in passing that we have checked that the particle production is independent ofθ 5+m if the fermion is massless. This should be the case since θ 5+m can be rotated away if the fermion is massless and hence θ 5+m is physical only when m = 0.
Non-relativistic limit. We now proceed to interpret our numerical results. For this purpose, we study the dispersion relation of the Bogoliubov coefficients. As we saw in Sec. 2.1, the eigenvalues Ω of the Hamiltonian are independent ofθ 5+m . However, since the equation of motion of the Bogoliubov coefficients depends onθ 5+m , their time evolution is not simply governed by Ω. Indeed, by taking the non-relativistic limit, we show in App. E that the non-zero axion velocity induces the following operators (in the particle sector): where η is a two-component spinor corresponding to the positive frequency part, σ is the Pauli matrix and Π = p x , p y , Π z . Thus, the axion induces a coupling between the spin and momentum. Accordingly we obtain for the time-dependent eigenfrequency, by diagonalizing the equation of motion of η. 5 These eigenvalues have an interesting property. For definiteness, let us takeθ 5+m > 0 and consider the minimum value ofΩ − NR with respect to Π z . The minimum value isΩ − NR = p 2 T /(2m) at Π z = 0 forθ 5+m = 0. Onceθ 5+m is turned on, however, the minimum value of Ω − NR is smaller than p 2 T /(2m), and eventually becomes zero whenθ 5+m > p T . In other words, the spinmomentum interaction induced by the axion velocity compensates the gap from the transverse momentum. We thus naturally expect that the exponential suppression from p T is compensated by the axion velocity, which is exactly what we find in Fig. 1.
Semi-analytical results. The above observation relies on the non-relativistic limit, but this limit is not essential. Indeed, as we see in App. D, the full relativistic equation of motion of the Bogoliubov coefficients has the following eigenvalues (among others), which reduce to Eq. (2.22) in the non-relativistic limit. As an empirical proof that these modes play the essential role, we find that the following integral approximates the numerical results well for m, p T g |Q| E : 6 whereθ 5+m is replaced byθ inΩ − in this expression, and In short, we interpret our numerical results as follows. The axion motion induces the spin-momentum interaction. It modifies the dispersion relation and in particular can compensate the energy gap from the transverse momentum p T (see Eq. (2.23)). The remaining minimal gap is given solely by the intrinsic mass m in the largeθ 5+m limit, and hence the particle production is suppressed only by exp(−πm 2 /g |Q| E ). In other words, the particle production is exponentially enhanced by exp(πp 2 T /g |Q| E ). The UV cut-off of this process is set byθ 5+m , as will be discussed in more detail in Sec. 3.3.

Dirac equation with classical background fields
This section generalizes the analyses performed in Sec. 2 by including a constant magnetic field, aligned (anti-)parallel to the electric field. Such a configuration is generated dynamically through a tachyonic instability in the axion gauge field coupling, leading to an exponential enhancement of one of the two gauge field helicities [10][11][12]. We will thus consider the background field configuration where the magnetic field B is constant and as before the electric field is given by −Ȧ z . Here, our starting point will be the most general action coupling axions, fermions and gauge bosons through dimension five operators, while preserving the shift-symmetry of the axion and C P -invariance, This in particular extends the discussion of Ref. [21], which assumed for simplicity θ 5 + θ m = 0.
The introduction of the magnetic field induces a spiralling motion of the charged fermions. This spiralling motion is quantum-mechanically quantized and the dispersion relations are consequently described by discrete Landau levels (see Ref. [21] for an explicit construction). The equation of motion for the lowest Landau level reads where s = sign(QB ). We note that in this lowest Landau level, the effect of θ 5 is degenerate with Π z . In particular, ifφ is constant, we can absorb it by shifting the momentum p z . Hence, since the computation of particle production involves an integral over p z , we can anticipate that the particle production in the lowest Landau level will not be affected by the axion motion.
Similarly, introducing the magnetic mass m 2 B = 2ng |QB |, the equation of motion for the higher Landau levels reads where n = 1, 2, .. labels the Landau levels and ψ n is a vector containing the four fermionic modes, i.e., (anti-)particles of both helicities, which mix for a given Landau level, see Ref. [21] for details.
As in Sec. 2.1 we may now expand the fermionic mode functions as for the lowest Landau level, and with the effective mass labelling the Landau level, m T = m 2 B + m 2 .
Bogoliubov coefficients. We again turn on the time-dependence of A z and φ, thus promoting the coefficients α λ and β λ to time dependent Bogoliubov coefficients. Proceeding as in Sec. 2.1, the equations of motion are given byα for the lowest Landau level, and for the higher Landau levels, where Θ 0 = dt Ω 0 and Θ n = dt Ω n . In particular, we see that Eq. (3.10) is equivalent to its counterpart without the magnetic field (2.10) after replacing p T with m B . The only difference is that the transverse momentum is quantized due to the magnetic field in Eq. (3.10), while it is continuous in Eq. (2.10).
Quantization and initial conditions. In terms of these Bogoliubov coefficients, the quantized mode function is given by for the lowest Landau level, and for the higher Landau levels, where b 0 , d 0 , b n,λ , d n,λ are the fermion creation and annihilation operators.
The initial conditions are given by for the lowest Landau level and for the higher Landau level, respectively, at the initial time t = −∞.
As in the case without a magnetic field, one can define the fermionic annihilation and creation operators at any given time t as for the lowest Landau level, and for the higher Landau levels. They satisfy the standard anti-commutation relations, indicating that for the lowest Landau level, and for the highest Landau level, for all times t . The former trivially follows from Eqs. (3.8) and (3.9), while the latter is shown to be satisfied in App. D.

Anomaly equation
The anomaly equation relates the chiral fermion current with the Chern-Simons term as [15,16] An important property of the system in this section is that the Chern-Simons density is non-vanishing, FF = −4E B = 0. Therefore a non-trivial consistency check of our computation is to correctly reproduce the anomaly equation. Ref. [21] shows that the anomaly equation holds in the case θ 5 + θ m = 0, and we now generalize this result including a non-vanishing θ 5 + θ m .
In the following we focus on the spatially averaged version of the anomaly equation, which readṡ with the expectation value 〈· · · 〉 taken with respect to the initial vacuum state. We consider contributions from the lowest and higher Landau levels separately. In particular, we will see that only the lowest Landau level contributes to the Chern-Simons term. See App. C for some useful relations that are used in the computation below.
Lowest Landau level. The chiral charge from the lowest Landau level is given by 7 where the subscript "0" indicates the contribution of the lowest Landau level. With Eqs. (3.8) and (3.9), its time derivative readsq The first term is easily integrated, and the second term corresponds to the mass term, where "LLL" stands for the lowest Landau level. As a result, we obtaiṅ Thus, the Chern-Simons term is supplied by the lowest Landau level. Below we confirm that the higher Landau levels do not induce additional contributions to the Chern-Simons term.
Higher Landau level. The chiral charge from the higher Landau levels is given by where the subscript "n" indicates the contribution of the nth Landau level for n ≥ 1. By taking the time derivative and using Eq. (3.10), we obtaiṅ It is straightforward to show that the right-hand-side corresponds to the mass term, where "HLL" stands for the higher Landau levels. Therefore the higher Landau levels do not contribute to the Chern-Simons term. This completes the proof of the anomaly equation.

Particle production
We now study the particle production for the lowest and higher Landau levels separately, and estimate the induced current. We impose the electric field and the axion dynamics as in Eqs. (2.18) and (2.19). 7 We dropped a regulator in this expression. One can show, as in Ref. [21], that the results do not depend on the choice of a regulator function as long as it depends only on Ω 0 .  Lowest Landau level. We have checked numerically that the spectrum of the lowest Landau level does not depend onθ. This result is easily understood based on our discussion in Sec. 2.2 as follows. The lowest Landau level corresponds to the mode that moves parallel to the magnetic field, or equivalently has a vanishing transverse momentum. In this case, the minimum of Ω − is not affected by the presence oḟ θ 5+m , and hence no exponential enhancement of the particle production is expected to occur. The nonvanishing axion velocity simply leads to the replacement p z by p z −θ 5+m in the eigenvalue, and this is absorbed by a constant shift of p z for a constantθ. As a result there is no effect on the particle production of the lowest Landau level.
Higher Landau levels. In Fig. 2, we plot our numerical results for the particle production for the higher Landau levels. It shows that the spectrum is again exponentially enhanced whenθ is large. As we noted above, the equation of motion of the Bogoliubov coefficients for the higher Landau levels are equivalent to Eq. (2.10) after replacing p T with m B . Therefore, we can interpret this result in exactly the same way as we did in Sec. 2.2. In particular, the enhancement of the particle production is well approximated by Eq. (2.24) after replacing p T with m B as we show in the right panel. The axion-induced spin-momentum interaction does not care whether the transverse momentum is continuous or discretized, and the axion assisted Schwinger effect is at work for both cases.

Induced current.
We finally estimate the current of the produced fermions. As discussed in Sec. 4, this is a key quantity to determine the backreaction of the fermion production on the background gauge fields in a dynamical system. In our setup, only the z-component of the induced current is non-vanishing, and is given as We focus on the first term since it is proportional to the duration τ of the non-zero electric field, and hence dominates for large τ. As we saw before, the spectrum develops a plateau approximated as where we denote the height of the plateau as β n 2 , Θ is the Heaviside theta function and sgn(Q) = Q/ |Q|.
Therefore the induced current is estimated as where we assumed τ m T /(g |Q| E ).
If there is no coupling to the axion,θ 5+m = 0, the height is given by and hence the induced current is estimated as reproducing our previous result [21].
Once we turn on the axion coupling, the spectrum for the higher Landau level is exponentially enhanced and β n 2 is estimated by Eq. (2.24) after replacing p T with m B . Unfortunately, Eq. (2.24) is still complicated enough so that we could not obtain the induced current analytically. Therefore, we just make a crude estimation of the induced current by further simplifying Eq. (2.24). First, we approximate the oscillatory behavior of β n 2 with respect toθ by simply inserting half the envelope of the oscillation. We also focus on the modes that satisfy m B g |Q| E since otherwise the axion assisted Schwinger effect is not so drastic (see App. D for a discussion of the m g |Q|E limit). With these simplifications, we find As shown in Fig. 2, this formula describes the envelope well in the asymptotic regime. We see that the suppression factor approaches exp(−πm 2 /(g |Q| E )), i.e., exp(−πm 2 m 2 B /(2g |Q|Eθ 2 ) 1, when where we simply count the number of the modes that satisfy Eq. (3.38) and introduce the "max" function so that it reduces to the previous result whenθ 5+m is small. One can see that the induced current is enhanced asθ 5+m increases. In particular, for E ∼ |B |, the axion assisted Schwinger effect leads to an enhancement of the induced currently by a factor of roughlyθ 2 /m 2 compared the standard result in the absence of the axion field. Here we emphasize again that our estimation above is quite rough, and we leave a more precise estimation to future work.

Implications for axion cosmology
Axion inflation. Identifying the inflaton, i.e., the particle responsible for driving cosmic inflation, with an axion-like particle with shift-symmetric dimension five couplings to gauge fields and fermions, naturally ensures a sufficiently flat scalar potential as required for slow-roll inflation [24,25]. The coupling to gauge fields induces a tachyonic instability in one of the helicities of the vector potential, leading to the production of a strong, large-scale helical gauge field configuration during inflation, driven by the kinetic energy of the axion field [10][11][12]. 8 The production of fermions in this helical gauge field background is well described by the analysis of Sec. 3 as long as the axion velocity varies only slowly, implying approximately constant physical electric and magnetic field strengths. 9 In the absence of charged fermions, the exponential gauge field production in axion inflation leads to striking signatures [36], including the generation of gravitational waves in the range of ground-and spacebased interferometers [37,38] and of primordial black holes [39][40][41]. The dual production of helical gauge fields and charged fermions was first studied in Ref. [19] for massless fermions and extended in Ref. [21] to massive fermions, for the particular parameter choice of c m + c 5 = 0, corresponding to the absence of the last term in Eq. (3.2). The fermion production and the resulting induced current lead to the formation of electric and magnetic fields anti-aligned to the background fields, and thus to a reduction of the net gauge field background generated in axion inflation by several orders of magnitude. This dramatically changes the predictions of axion inflation.
Our new estimate for the induced current, Eq. (3.39), indicates that for suitable values ofθ, the induced 8 As in the rest of this paper, we focus on Abelian gauge fields here. Non-Abelian gauge field configurations can also be sourced by a non-vanishing axion velocity [26][27][28][29], however due to the self-interactions of the non-Abelian gauge fields, the fermion backreaction originating from the induced fermion current is less relevant in this case [30]. Moreover, for a discussion of the gravitational production of neutral fermions in axion inflation, see Ref. [31]. 9 For very strong gauge field backgrounds with correspondingly strong backreaction effects on the axion equation of motion, this approximation becomes invalid due to resonance effects in the coupled axion gauge field system [32][33][34][35]. However, in the presence of light fermions, the gauge field production is inhibited, and hence the resonance effects are expected to be less rele- In axion inflation, we typically expect H 2 /E ∼ 10 −4..−2 [21] and (M P / f ) 2 10 3 [36], indicating a potentially sizable enhancement of the induced current. This is in particular true towards the end of inflation, when ∼ 1. Consequently, we expect a further reduction of the gauge fields production compared to the analysis of Ref. [21].
Moreover, due to the slow variation ofφ over the course of slow-roll inflation, we expect to scan the oscillatory features of Fig. 2. For moderate values ofθ when only a small number of Landau levels contribute significantly to the induced current, this may lead to an oscillation in the fermion backreaction which could induce oscillations in gauge field background, and thus, through the gauge friction effect, in the inflaton velocity. Since the approximation of an adiabatically varying inflaton velocity and gauge field strength may fail in this case, a detailed analysis of this system requires non-linear methods such as lattice simulations. Qualitatively, resonance effects similar to Refs. [32][33][34][35] may occur, leading to characteristic 'spikes' in the spectra of scalar and tensor perturbations. We leave a thorough investigation to future work.
Axion dark matter. Axion-like particles are intriguing dark matter candidates. Here we are particularly interested in ultralight (sub-eV) axions, which are stable on cosmological time scales and can be described by a coherently oscillating axion background field. 10 In the rest-frame of an earth-based laboratory, this leads to an 'axion wind', with characteristic frequencies associated with the motion with respect to DM halo of our galaxy as well as with the axion mass, m φ . For an overview of axion searches exploiting this effect, see Ref. [42].
In principle, earth-based experiments aimed at detecting Schwinger production in strong electric fields may thus be sensitive to the axion-electron coupling through a distortion of the high-momentum tail of the electron spectrum, resulting from the enhanced electron production for m T > m e depicted in Fig. 1, with m e denoting the electron mass. However, sinceθ < m φ m e ∼ g E , the enhancement is only mild (see Fig. 1). Since moreover Schwinger production has not yet been successfully observed in the laboratory due to the experimental challenges involved [43], this conceptually interesting observation seems of little practical use in the immediate future.

Discussion and conclusions
In this paper, we study Schwinger production of charged fermions in the background of a homogeneous axion field with a non-vanishing velocity. By numerically solving the Dirac equation in background gauge and axion fields, we obtain time-dependent Bogoliubov coefficients describing the non-perturbative particle production. We find that the Schwinger production rate is exponentially enhanced when the axion velocity is sufficiently large (1.3) and the transverse momentum of the produced particle is non-zero, which we dub axion assisted Schwinger effect. We also provide a semi-analytic expression for the number densities of produced particles, which well explains our numerical results [see Eq. (2.24) and below].
Throughout this paper, we allow for general dimension-five couplings which preserve the axion shift symmetry. By means of a chiral rotation, some couplings can be expressed by the others, implying that unphysical. For this reason, the axion assisted Schwinger effect is more pronounced for a larger mass, while we should bear in mind that the overall production rate gets more suppressed as Eq. (1.3). As a consistency check, we also confirm that our result reproduces the earlier studies [17][18][19][20][21] in the limit of a vanishing θ 5+m or massless fermions (see also Figs. 1, 2, and 5).
Based on these results, we discuss phenomenological implications of the axion assisted Schwinger effect. The enhanced production rate results in an enhancement of the induced current, which backreacts on the gauge field equation of motion. This is, in particular, relevant for axion inflation which predicts the production of helical gauge fields. Since the backreaction is enhanced by the axion assisted Schwinger effect, we expect a reduced helical gauge field production in parameter regimes where the axion assisted Schwinger production is relevant. Moreover, the enhanced production rate predicts a distortion of the high-momentum tail of the electron spectrum in laboratory experiments aiming at measuring the Schwinger mechanism, such as the X-ray laser XFEL and the extreme-light infrastructure ELI [43][44][45]. We briefly discuss its implication in the case of axion dark matter, but the enhancement is small for the parameters of our interest.
Acknowledgements It is a pleasure to thank Ben Mares for helpful discussions related to this project.
This work was partly funded by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306.

A Notation and conventions
In Eq. (3.2) we have introduced the comoving quantities ψ, A µ and g µν , related to the corresponding phys- Using the chiral representation of the γ matrices, (γ µ ) = (γ 0 , γ) with the left(right-)handed component of the four-spinor ψ = (ψ L , ψ R ) is projected out by the projection operator P L/R = (1 ∓ γ 5 )/2. The (dual) field strength tensor of the gauge field is given by with 0123 = +1.

B Particles and antiparticles
Hamiltonian. To introduce the notion of particles and anti-particles, we derive the Hamiltonian density of our system. The conjugate momentum are given from Eq. (3.2) with a = 1 as The Hamiltonian density is then given by For our purpose the fermion part is important. Imposing the equation of motion for ψ, it is given by Wave function. In the absence of any magnetic field, we decompose the fermion mode function in with Here the superscript (σ) labels left-and righthanded particles in the m → 0 limit. The meaning of the subscript s will be more transparent once we include magnetic fields, for the moment this just labels to two linearly independent modes associated with each value of λ.
Positive and negative frequency modes. Inserting Eq. (B.6) into Eq. (B.5) gives where we have used the notation ψ s ≡ σ ψ (σ) s χ (σ) s . We note that the time derivative is shifted by iθ 5 γ 5 , and consequently the eigenstatesψ of the Hamiltonian are given by extracting this phase factor,ψ = exp(−i γ 5 θ 5 )ψ. In this basis, the equation of motion (2.4) simplifies to Inserting the ansatzψ ∝ exp(−i Ωt ), we find which in particular does not depend onθ 5 orθ m .
We now proceed to express our wave function as a decomposition of positive and negative frequency states, where u λ , v λ are eigenvectors of the Hamiltonian (B.5) for constant A z : with the normalization factor Here we define the transverse momentum and choose the phases of u 2 and v 2 such that the analogy to the case with the magnetic field becomes transparent. Note that the eigenvalues Ω associated with u 1 and u 2 (and correspondingly v 1 and v 2 ) are degenerate. Here we have chosen u n,2 and v n,2 so as to match u 0 , v 0 for p T → 0.

Some useful relations.
Computing the equations of motion for the coefficients α and β requires the evaluation of inner products among the vectors u λ , v λ and their time derivatives in the presence of timedependent A z and θ 5+m . We give some useful relations in the following:

C Particles and antiparticles -with magnetic field
Wave function. In the presence of (anti-)parallel electric and magnetic fields pointing in the z-direction, we can expand the fermion wave function as where we have performed a Fourier-transform with respect to the y and z direction,x s = g |QB | x − s p y g |QB | and h n is related to the Hermite function that satisfieŝ with the ladder operatorsâ Positive and negative frequency modes. As before, we obtain the eigenstatesψ of the Hamiltonian by extracting the θ 5 phase factor,ψ = exp(−i γ 5 θ 5 )ψ. In this basis, the equations of motion for the lowest and higher Landau levels, Eqs. (3.3) and (3.4), simplify to where we used the short-hand notatioñ Inserting the ansatzψ ∝ exp(−i Ωt ), we find which is in particular does not depend onθ 5 orθ m .
We now proceed to express our wave function as a decomposition of positive and negative frequency states,ψ α n,λ u n,λ exp (−i Ω n t ) + β n,λ v n,λ exp (+i Ω n t ) , n ≥ 1 , (C.9) where u 0 , v 0 and u n,λ , v n,λ are eigenvectors of the Hamiltonian (B.5) for constant A z . For the lowest Landau level, this yields 2Ω 0 (Ω 0 + sΠ z ) whereas for the higher Landau levels we find where the normalization factor is Note that the eigenvalues Ω n associated with u n,1 and u n,2 (and correspondingly v n,1 and v n,2 ) are degenerate. Here we have chosen u n,2 and v n,2 so as to match u 0 , v 0 for m B → 0 for s > 0.

Some useful relations.
Computing the equations of motion for the Bogoliubov coefficients requires the evaluation of inner products among the vectors u 0 , v 0 and u n,λ , v n,λ . We give some useful relations in the following (dropping for notational simplicity the index n for the higher Landau levels): for the lowest Landau level, and correspondingly for the higher Landau levels.
In addition, the following relations are useful to show the anomaly equation. For the lowest Landau level, and for the higher Landau levels, are useful relations for the evaluation of the chiral charge. Moreover, for the lowest Landau level, and for the higher Landau levels, with i , j ∈ {1, 2}, are useful for evaluating the mass term.

D Bilinear forms of Bogoliubov coefficients and particle production
In this appendix, we consider products of the Bogoliubov coefficients. The purpose of this appendix is twofold. First, we show the existence of the conserved quantities (Eqs. (2.16) and (2.17) in the case with only an electric field). Second, we motivate the analytical formula that we use to estimate the spectrum (2.24).
Although we could not derive this formula rigorously, we outline our computation that leads us to this formula, with the hope that one may find the argument there useful to derive a more complete analytical formula in the future. We consider the case with only an electric field in the following, but the results equally apply to the case with both an electric and magnetic field after replacing p T by m B .
For our purpose, it is convenient to treat the Bogoliubov coefficients associated with the two sets of initial conditions (see Eq. (2.12)) in a unified way. We thus define matrices as where the quantities in the left-hand-side are now understood as 2 × 2 matrices. These matrices satisfẏ where σ i is the standard Pauli matrix, and their initial conditions are given by We consider products of these matrices. They are again 2×2 matrices and hence can be expanded in terms of the Pauli matrix as Note that R µ α and R µ β are real while C µ are complex. The initial condition now reads R 0 α = 1 and R i α = R µ β = C µ = 0 at the initial time.
Conserved quantities. The conserved quantities (2.16) and (2.17) are expressed as These relations can be shown as follows. The equations of motion of R One can see thatṘ 0 α +Ṙ 0 β = 0, (D. 20) and it follows from the initial condition that R 0 α + R 0 β = 1. One can also see that the equations of motion of where we omit the subscript 5 + m here and here only for notational simplicity. We can formally solve this equation as follows. We define the transfer function as The parameters are shown in the unit g |Q| E = 1.
Empirical analytical formula. Unfortunately, we could not evaluate the expression (D.37) analytically in the case of our interest. Nevertheless we can learn properties of our system from this expression. In particular, we see that R 0 β is obtained after convoluting the source term with the phase factor e 2i dt λ n (t ) . If |λ n | is large, the integral is suppressed since the integrand oscillates rapidly. Thus, we expect that the smallest eigenvalue is the most important (see also App. F). One can see that the smallest eigenvalue for a sizableθ 5+m is given by if we ignore the terms proportional toΠ z in M , where we take Ω − forθ 5+m > 0 and Ω + forθ 5+m < 0. If this mode is the most important, we may expect that R 0 β is related to the following integral: with some function f (t ), where we assumeθ 5+m > 0 to be specific. This integral can be evaluated by the saddle point approximation. The saddle points of the phase factor are given by where σ, σ = ±. We take the saddle points in the upper half plane, and then the integral is given by where we changed the integral variable from t to Π z , and we ignored all the prefactors. Since R 0 β is positive definite, we may expect that it is given by the integral squared as where we chose the minus sign in the interference term simply because it describes the numerical results well. 12 We show in Figs. 1 and 2 that this formula reproduces the full numerical results extraordinary well for several model parameters. Since our setup only contains fairly few parameters, we can empirically gain confidence in the expression (D.42) by systematically varying all these parameters and comparing with the full numerical result. For large m and p T , we find excellent agreement for the height of the plateau of the spectrum, despite the rather heuristic approach, as depicted in Fig. 3 for a variety of model parameters.
Eq. (D.42) deviates from the full results when m and/or p T is small as shown in Fig. 4. We note that the case of small p T is however not particularly relevant for our discussion since the axion assisted Schwinger effect is most prominent when p T is large. It would be certainly interesting if one could derive an analytical formula that works for both small and large values of m and p T , which we leave as a future work. Small mass/transverse momentum limit. We finally comment on the limit m → 0. Our numerical result shows that the enhancement becomes less significant and eventually dies out as m gets smaller, as shown in the left panel of Fig. 5. In particular, we observe a scaling well approximated as λλ |β (λ ) λ | 2 ∝ m 2 /(g |Q|E ) for m 2 g |Q|E . Since θ 5+m is unphysical when m = 0, this result is consistent with a continuous m → limit.
We also comment on the limit p T → 0. As depicted in the right panel of Fig. 5, the axion assisted Schwinger effect again becomes less significant and eventually dies out in this limit. Although our formula (D.42) does not reproduce the numerical results accurately in this limit as we mentioned above, the deviation is at most factor two simply because the suppression due to p T becomes irrelevant in this limit. 12 See also Ref. [23], which traces the minus sign in the interference term of back to the fermionic nature of the produced particles in the context of the dynamically assisted Schwinger mechanism.
Effective action. As demonstrated in the main text, Schwinger pair-production sourced byφ is enhanced when the mass term cannot be neglected. This implies that it should be also understood in non-relativistic effective field theory of fermions. Let us take the Dirac basis for the γ matrices in this section (and this section only) as they are convenient in the non-relativistic limit: In this basis, the upper and lower components of fermion fields stand for particle and anti-particle: We have factored out the fast oscillation coming from the mass term and η (ξ) represents a slowly varying field for the (anti-)particle. One may derive the non-relativistic effective action of Eq. (3.2) in terms of η and ξ by performing the expansion of 1/m systematically, which reads Here we have assumed that the axion field is homogeneous and slowly varying compared to 1/m. Notice that the effective action solely depends on the particular combination θ 5+m signaling the basis independence under rotations of the fermion fields.
The second (and correspondingly the fourth) term can be rewritten as follows 1 2m where the gauge invariant momentum is defined by Π = −i D = −i ∇ − gQA. The first two terms are the ordinary kinetic term and the magnetic dipole interaction, respectively. The other terms depend onθ 5+m .
In particular, the third term provides the spin-momentum interaction induced byθ 5+m . This operator reduces (increases) the energy of a particle for a given Π = |Π| if its spin is (anti-)parallel to Π. The same argument holds for the anti-particle ξ. As a result, we have two modes in total that are produced more efficiently than in the case withoutθ 5+m .
Dispersion relation. Let us first derive the dispersion relation explicitly without the magnetic field, i.e., A µ = (0, 0, 0, A z ). The case with the magnetic field parallel to the electric field will be mentioned shortly after. The equation of motion in this case reads with ϕ p = arctan(p y /p x ). Here we expand the mode as whereη represents a field in the new basis. For a constant Π z , one may solve this equation with two frequencies:Ω Now it is clear that the production of the mode withΩ ± is enhanced forθ 5+m ≶ 0, since the gap in the corresponding dispersion relation is reduced. It is straightforward to obtain the same dispersion relation for the anti-particle.
Including the magnetic field. Once we turn on a magnetic field parallel to electric field, the momentum transverse to the electromagnetic field becomes discretized as the Landau levels. The equation of motion is obtained by just replacing m T with m B = 2ng |QB | as Note that the mode expansion is also modified as η p y ,p z (x) = n,σ η σ n h n S z σ . (E.10) The dispersion relation is now given bỹ (E.11) In the same way as the case without the magnetic field, the production of the mode with Ω ± n is enhanced forθ 5+m ≶ 0 respectively.

F Phase integral method
The exponential suppression of Schwinger particle production can be elegantly understood using the phase integral method for quantum mechanical scattering problems (see e.g. [47]), based on the WKB approximation. Here we briefly review the essence of this method, as it provides useful intuition for understanding the role of the gap in the dispersion relations between particles and antiparticles, as well as for understanding interference effects between several (complex) zero-points of the dispersion relation. See Refs. [23,48,49] and references therein for derivations and applications to particle production in scalar and spinor QED.
The computational problem at the heart of particle production in electromagnetic fields can be reduced to solving a harmonic oscillator equation with a time-dependent frequency, ψ(t ) + Ω 2 (t ) ψ(t ) = 0 . (F.1) Starting from the Dirac equation, a second order differential equation of the type of Eq. (F.1) is obtained by acting with the conjugate differential operator. In the WKB approximation, the solution is given by with the Bogoliubov coefficients α and β. Determining the asymptotic value β(t → ∞) for the initial condition β(t → −∞) = 0 thus requires solving an integral of the type with the coefficient C (t ) determined by the coupled first-order ODEs governing the evolution of the Bogoliubov coefficients. This integral can be evaluated by choosing a convenient contour in the complex plane, noting that for Ω(t ) real (which is e.g. typically the case along the real axis), the integrand of the outer integral becomes a highly oscillatory function and hence the contribution to I ± is negligible. Moreover, for i Ω is real and positive, the integrand of I + (I − ) is exponentially suppressed. This qualitatively explains why to good approximation, Eq. (F.3) can be evaluated by integrating along the Stokes line 13 connecting the pair of complex zeroes of Ω(t ), where for a constant electric fieldΠ z = gQE . For gapped dispersion relation, i.e. Ω(p) > 0 for all p ∈ R, the distance 2Im(Π z,0 ) between the pair of complex zeroes corresponds to the minimal energy gap which needs to be overcome for the production of a particle from the Dirac sea.
For the case c 5 + c m = 0 considered in the main part of this paper, an additional subtlety is that the eigenvalues of the Dirac equation (in the basis that diagonalizes the Hamiltonian (B.5)) do not coincide with the eigenvalues of the equations of motion for the Bogoliubov coefficients derived in App. D. However, once we have arrived at the set of coupled equations of motion (D.22), we can again interpret them as coupled oscillators with time-dependent frequencies, similar to Eq. (F.1). Consequently, the particle production is governed by an expression similar to Eq. (F.4), with Ω replaced by the smallest eigenfrequency of the system.

Interference effects.
In general (and in particular in the situation in the main part of this paper), the function Ω(Π z ) can have multiple complex zero points. Intuitively, this can be understood by recasting 13 The Stokes line is a path in the complex t plane along which Ωdt is imaginary. In practice, as long as on pays careful attention to potential branch cuts, on may equally well integrate on a contour parallel to the imaginary axis.
Eq. (F.1) as a quantum mechanical scattering problem over a time-dependent potential barrier with a potentially complicated shape [23]. Partial reflections of different parts of this barrier can lead to interference effects. Eq. (F.4) is then modified to [48], with the index i labelling the pairs of complex zeroes and θ i =Π −1 z Re(Π i z,0 ) Re(Π 1 z,0 ) Ω(Π z )d Π z accounting for the phase accumulated by the integration along the real axis. This leads to interference effects in the final result for |β| 2 .