Semiclassical analysis of axion-assisted and axion-driven pair production

We study the pair production of fermions in a time dependent axion background with and without an electric background. We construct the adiabatic mode functions which incorporate the gauge field and the axion velocity dependence of the dispersion relation. The semiclassical approach using this adiabatic basis shows two types of pair production. One is axion-assisted pair production: the presence of the axion velocity gives enhancement and interference effects on the pair production driven by the electric field. The other is axion-driven pair production: the time variation of the axion velocity causes the pair production even though the electric field is absent.


Introduction
The presence of a strong electric background causes the pair production of charged particles [1]. The so-called Schwinger effect is the most well known example of pair production. We here note that some cosmological scenarios treat electromagnetic fields with time evolution of an axion (or an axion-like particle), which is a pseudo-Nambu-Goldstone boson coupled to a photon via an anomaly. As an example, in the scenario of the magnetogenesis by the axion inflation [2][3][4][5][6][7][8], a large electric field as well as a magnetic field is generated by an axion velocity via the Chern-Simons coupling.
Motivated by the application to the above scenario, Domcke et al. studied the pair production of charged fermions in a constant electric field and a constant axion velocity [9]. 1 They found that the presence of the axion velocity enhances the amplitude of the distribution function of produced fermions, and the distribution function shows an oscillatory behavior with the change of the axion velocity. They named this phenomenon axion-assisted Schwinger effect, and empirically proposed an analytical formula which may describe the distribution function of produced fermions.
In this paper, we derive the analytical formula for the axion-assisted Schwinger effect by using the semiclassical approach [10][11][12][13][14]. The fundamental strategy of the semiclassical 1 They also studied the more appropriate situation for the magnetogenesis scenario, where a constant magnetic field is present. The presence of the magnetic field discretizes quantum numbers while does not change the roles of the axion velocity in the axion-assisted Schwinger effect. In this paper, for simplicity, we discuss the case where the magnetic field is absent. approach is as follows. We expand the Dirac field by using the lowest-order WKB solutions (we call them the adiabatic mode functions in this paper) as a basis. The Bogoliubov coefficients are given by the coefficients of the adiabatic basis, and evolve with time as long as the backgrounds included explicitly in the adiabatic basis are time dependent.
The semiclassical approach is used also in [9], while the adiabatic basis adopted in the previous study does not include the axion velocity dependence explicitly. In this paper, we construct the adiabatic mode functions which incorporate the gauge field and the axion velocity dependence of the dispersion relation. Our semiclassical approach, using the axion velocity-modified adiabatic basis, clearly shows the roles of the axion velocity in the axionassisted Schwinger effect, and gives the analytical formula for this effect.
Furthermore, our semiclassical approach shows another type of pair production. If the axion velocity changes with time, the pair production occurs even though the electric field is absent. We here call this phenomenon axion-driven pair production. For the axion-driven pair production, the distribution function of produced fermions is consistent with that found in [15,16]. This mechanism can be particularly important for the kinetic misalignment mechanism proposed in [17,18], where large axion velocity and axion acceleration are induced in the early Universe.
The organization of this paper is as follows. In Sec. 2, we construct the adiabatic mode functions with the gauge field and the axion velocity dependence, and thus derive the equation for the Bogoliubov coefficients based on the axion velocity-modified adiabatic basis. In Sec. 3, we consider the case where the electric field is present and the axion velocity is constant, and study the axion-assisted pair production (the axion-assisted Schwinger effect). In Sec. 4, we consider the case where the electric field is absent and the axion velocity changes with time, and study the axion-driven pair production. Sec. 5 is devoted for conclusion and discussion.

Semiclassical approach
We consider a Dirac field ψ which couples to a U (1) gauge field A µ and an axion(-like particle) φ through the following nonlinear terms, where θ 5 = c 5 φ/f a and θ m = c m φ/f a , f a is the axion decay constant, and c 5 , c m are the dimensionless coupling constants. See Appendix A for our notation of the gamma matrices.
We consider an expanding universe with a scale factor a, while the scale factor dependence has been almost eliminated in (2.1) by the conformal transformation a 3/2 ψ → ψ. Unless otherwise stated, the coordinates are given by the conformally flat ones (i.e., x 0 is given by the conformal time τ ), and the Lorentz indices are lowered and raised by η µν = diag(−1, 1, 1, 1) and its inverse η µν , respectively. Although the mass term is time dependent as m = am phys (m phys is the physical scale of the mass), we assume that the pair production occurs at the short time scale where the time evolution of the mass term is negligible. 2 Using the chiral transformation ψ → e −iθmγ 5 ψ, the two axion couplings in (2.1) are unified as follows where θ = θ 5 + θ m . In the subsequent discussion, ψ denotes this conformal and chiral transformed Dirac field. In this paper, the gauge field and the axion are treated at the background level. Specifically, we consider the situation where these background fields are homogeneous, 3) The action (2.2) includes the axion velocity ∂ 0 θ rather than θ itself. We assume that the axion velocity does not change its sign in the time range of our interest. In such a situation, we can set it to be positive, ∂ 0 θ ≥ 0 without loss of generality. After the spatial Fourier transformation, the Dirac equation is given by wherek i = (k 1 , k 2 , k 3 − eA) and k i is the comoving momentum. Inserting the following equation we obtain the Klein-Gordon like equation First, we derive the lowest-order WKB solutions of the Dirac equation (we call them the adiabatic mode functions in this paper). We can derive the adiabatic mode functions by treating A µ and ∂ 0 θ as if they were constant. The Dirac equation (2.4) includes A µ and ∂ 0 θ, but does not include their time derivatives. For the Klein-Gordon like equation (2.6), the last two terms in the left-hand side are neglected in deriving the adiabatic mode functions.
For the positive-frequency adiabatic mode functions, we replace the time derivative with the frequency, (2.7) 2 Specifically, we consider the parameter region −e∂0A or ∂ 2 0 θ m∂0a/a. Here, we compare the time derivatives of the backgrounds which appear linearly in the Dirac equation (2.4).
From (2.7) and (2.8), the Klein-Gordon like equation (2.6) reduces to We thus obtain the explicit forms of the frequencies, where we consider the parameter regionk ≥ ∂ 0 θ. In this parameter region, the positivity of ω 2 k,± − m 2 allows one solution for each sign choice in (2.10). 3 The dispersion relation with the axion velocity dependence can be seen also in [9,15,16]. The presence of the axion velocity resolves the degeneracy, ω k,+ < ω k,− . (2.12) We thus call ω k,+ the lower frequency, and call ω k,− the higher frequency. Recall that we set ∂ 0 θ to be positive. The sign flip of ∂ 0 θ just exchanges the roles of ω k,± . The positive-frequency adiabatic mode functions u k,± are given by (2.14) where k ⊥ = k 2 1 + k 2 2 , e iϕ ⊥ = (k 1 + ik 2 )/k ⊥ , and the time dependent phase factors and the normalization factors are given by For convenience, we introduced the constant phase factor e −iϕ ⊥ /2 in the second lines of (2.13) and (2.14).
In a similar way, where we replace ∂ 0 with iω k,± , the negative-frequency adiabatic mode functions v k,± are constructed as follows whereν k,± are the eigenvectors of the matrix term in (2.6), (2.20) and the normalization factors are given bỹ For convenience, we introduced the constant phase factor e −iϕ ⊥ /2 in the second lines of (2.17) and (2.18). Although we derived (2.11), (2.13), (2.14), (2.17), and (2.18) fork ≥ ∂ 0 θ, the same frequencies and the same adiabatic mode functions hold true regardless of the sign of (k − ∂ 0 θ). See Appendix B for this universality. These adiabatic mode functions satisfy the following equations by definition, and they are orthonormal, The positive and the negative frequency modes are exchanged by the charge conjugate transformation χ C ≡ γ 2 χ * (where χ is an arbitrary four-spinor), where we define that ϕ ⊥ in u k,s shifts by +π while ϕ ⊥ in v k,s shifts by −π under k → −k.
In order to obtain this simple relation, we introduced e −iϕ ⊥ /2 into the adiabatic mode functions. Second, as with [10][11][12][13][14], we expand the Dirac field by using the adiabatic mode functions as a basis. At the initial time τ ini , we introduce the annihilation and the creation operators as follows where we use σ rather than s to emphasize that this label is defined independently of time.
The Dirac field at an arbitrary time can be expressed as follows where the Bogoliubov coefficients α k,sσ , β k,sσ take the following initial vales Unless both A µ and ∂ 0 θ are constant, the adiabatic mode functions do not satisfy the Dirac equation (2.4), and instead satisfy (2.22). Since the left-hand side of (2.28) satisfies the Dirac equation, the Bogoliubov coefficients evolve with time such that the right-hand side satisfies the Dirac equation. In other words, the pair production occurs as long as the electric field or the axion acceleration is nonzero, Let us derive the equation for the Bogoliubov coefficients. From (2.23) and (2.28), the Bogoliubov coefficients are given by Furthermore, using (2.4) and (2.22), their time derivatives are expressed as follows We list the explicit forms of the products of the adiabatic mode functions and their time derivatives in Appendix C. Substituting (C.1) into (2.32), we obtain (2.33) See (C.2) and (C.3) for the explicit forms of C k and D k . Note that the diagonal elements in the right-hand side are zero. The absence properly reflects that the pair production means the mixing of different modes with time, rather than the time evolution of each mode. From (2.33) and (2.29), we can confirm that the following unitary condition holds true, In general, each (|α k,sσ | 2 + |β k,sσ | 2 ) is not conserved, and |β k,sσ | 2 is not equal to the distribution function of produced fermions. This is because the equation (2.33) has the block-off-diagonal elements between the lower and the higher frequency modes. However, as discussed in the next section, we may practically neglect the contribution from the higher-frequency modes including the block-off-diagonal elements.
The equation (2.33) shows two types of pair production. If the axion velocity is constant, the pair production is driven by the electric field. The presence of the axion velocity modifies the electric pair production because the matrix with the coefficient eE depends on ∂ 0 θ (as discussed in the next section, it enhances the amplitude of the distribution function). If the axion velocity changes with time, the pair production occurs even though the electric field is absent; i.e., it is driven by the axion acceleration. We call the former one axion-assisted pair production (we also call it axion-assisted Schwinger effect by following [9]), and call the latter one axion-driven pair production. In the subsequent sections, we study these two types of pair production more specifically.

Axion-assisted pair production
We consider the case where the electric field is nonzero while the axion acceleration is zero in the whole time range, The electric field is set to be constant and positive for simplicity. We consider a sizeable axion velocity such that the axion-assisted effect becomes manifest, Note, however, that we impose the additional condition (3.9), which ensures the frequency difference between the adiabatic modes. The additional condition does not allow the trivial limit of (3.2): k ⊥ = m = 0. In this case, the pair production is driven by the first term in the right-hand side of (2.33). It is difficult to solve the differential equation with the 4 × 4 matrix exactly. However, the presence of the axion velocity resolves the degeneracy as seen in (2.12). Since the contribution from the lower-frequency modes may be dominant compared with that from the higher-frequency modes, we focus on the former one: In this approximation, the unitary condition holds true for the lower-frequency modes, and the distribution function for the lower-frequency modes is given by Extracting the upper left block from the 4×4 matrix (i.e., neglecting the block-off-diagonal elements), we obtain the reduced equation for (α k,+ , β k,+ ), In the subsequent paragraphs, we discuss in which parameter region the reduction process is justified, and how to obtain the analytical solution of (3.6).
In this paper, we consider the situation where the semiclassical picture is valid on the real time axis, where X k are the matrix elements without the oscillating parts (we consider the matrix elements which are relevant with the lower-frequency modes). The condition (3.7) can be understood as the adiabaticity condition, which ensures the instantaneous particle picture at a given time. 4 As we see in the subsequent discussion, it leads to the pair production whose distribution function is exponentially small. In such a situation, the turning points of the frequency in the complex time plane determines the Bogoliubov coefficients [19]. More specifically, given the convergence of contour integrals, we need to pick up the turning points in the complex lower-half time plane. In the case (3.1), X k are given by The explicit from of the validity condition (3.7) is given by The first condition comes from the conditions on the block-off-diagonal elements, and the second condition comes from the conditions on the block-diagonal elements. See Appendix D for the detail of the derivation. We already used the first condition in deriving (3.6). If this condition is not satisfied, the block-off-diagonal elements are not negligible. Let us see the turning points in the complex lower-half time plane. In the case (3.1), the lower frequency has two turning points i.e., ω k,+ (τ Near the turning points, the following identity holds true, 4 In short, (3.7) means that the time evolution of the Bogoliubov coefficients is slow compared with that of the adiabatic basis. For scalar fields, the equation for the Bogoliubov coefficients is given by .g., [13]), and thus the validity condition (3.7) reduces to the well known one: ∂ 0 ω k ω 2 k 2 1, where the upper and the lower choices of signs correspond to the behaviors around τ (1) 0 and τ (2) 0 , respectively. The equation (3.6) is thus written as follows We can approximately solve this equation by using the behavior of the frequency near the turning points: The solution is shown below in the form of the distribution function of produced fermions; i.e., n k,+ |β k,+ | 2 . See [12,14] for the detail of the calculation process.
For τ > Re τ (2) 0 , the distribution function is given by which consists of the contributions from the first and the second turning points τ 0 . The minus sign between the first and the second terms comes from the sign difference between the upper and the lower choices in (3.12). For Re τ (2) 0 > τ > Re τ (1) 0 , the distribution function is given by which consists only of the contribution from the first turning point τ 0 , the pair production has not been occurred. We set the initial time as τ ini < Re τ where any real value can be substituted for the lower bound of the integral. The choice of real values does not change the distribution function. The constant electric field was moved outside the integral. Using the imaginary and the real parts of Θ k,+ (τ (i) 0 ), (3.13) is expressed as follows and (3.14) is expressed as follows where the weight and the phase factors are given by As pointed out in [14], if there exist more than one turning points in the complex lower-half time plane, they may give an interference effect on the distribution function. In (3.13), the cross term between the contributions from τ gives the oscillatory function as . This is the second term in the first line of (3.16). More specifically, more than one turning points should have been passed. This is why (3.14), and thus (3.17), have no interference effect.
The same formula (3.16) was proposed as an empirical one in [9]. We gave the derivation of this analytical formula. In contrast to the previous study, the adiabatic mode functions were constructed corresponding to the ∂ 0 θ-modified frequencies (2.11). Using the appropriate adiabatic basis, we derived the equation for the Bogoliubov coefficients (2.33). This equation clearly shows which frequency modes are related to each matrix element. Therefore, we could derive the reduced equation for the lower-frequency modes (3.6). From the reduced equation, we derived (3.16) by picking up the contributions from the turning points.
We here express W k and P k more explicitly. The complex integrals (3.15) are not exactly solvable except in the k ⊥ = 0 limit. Instead, we evaluate them up to O(k 2 ⊥ ), and thus obtain the following expressions of W k and P k , where S ≡ (∂ 0 θ) 2 + m 2 .
As ∂ 0 θ becomes larger than m and k ⊥ , W k and P k approach the following values The large axion velocity eliminates the k ⊥ dependence from the weight factor. In this regard, the presence of the axion velocity enhances the amplitude of the distribution function. We emphasize that it also enhances the amplitudes of integral quantities like the induced current. From (III.20), we can estimate that the momentum integral over k ⊥ induces √ eE · ∂ 0 θ/m (rather than √ eE) per each power of k ⊥ . We also explain the overall factor of the exponential. This factor is given by the square of the number of relevant turning points; i.e., 4 in (3.16) while 1 in (3.17). 5 Furthermore, the distribution function shows the oscillatory behavior with the change of the axion velocity. In particular, at ∂ 0 θ √ nπeE (n ∈ N), the distribution function can take a tiny value beyond the enhancement effect on the amplitude. In evaluating integral quantities, we may move the sine squared outside the k ⊥ integral because P k (∂ 0 θ) 2 /(eE) for k ⊥ √ eE · ∂ 0 θ/m. Although these features were already found in [9], we explicitly evaluated the coefficients of the k 2 ⊥ terms of W k and P k .
We here mention the ∂ 0 θ = 0 limit where the standard Schwinger effect should be reproduced. The distribution function for the standard Schwinger effect cannot be reproduced naively by taking the ∂ 0 θ = 0 limit of (3.16) with (3.20) and (3.21). In deriving this result, we used (3.2) and (3.9), which do not allow the ∂ 0 θ = 0 limit. Let us get back to the equation for the Bogoliubov coefficients (2.33). At first glance, the ∂ 0 θ = 0 limit of (2.33) looks different from the equation for the Bogoliubov coefficients for the standard Schwinger effect. We should recall that the degeneracy of the frequencies recovers in the absence of the axion velocity (see (2.11)). Making use of the rotational transformation which mixes u k,+ and u k,− (and also mixes v k,+ and v k,− ), we can reproduce the equation for the Bogoliubov coefficients for the standard Schwinger effect.
In Fig. 1, we plot the numerical solution of (2.33) under ∂ 2 0 θ = 0 and that of (3.6), fixing ∂ 0 θ, and changing τ . For the full equation, the distribution function for the lowerfrequency modes is given by which reduces to (3.5) in the limit (3.3). The numerical solution of the full equation is so well described by that of the reduced equation that we can barely distinguish between both lines. This figure shows that the pair production occurs on two peaks, and there exist two steps: the first step appears between the two peaks, and the second step appears after the later peak. It is reasonable to expect that (3.17) describes the value of the first step, and (3.16) describes the value of the second step. In Fig. 2, we verify the above expectation: we compare the ∂ 0 θ dependence of the numerical solution of the full equation with that of the analytical solution, at two fixed times. In the upper panel, the time is fixed at the value inside the first step (τ = 0) for the numerical solution, and (3.17) with (3.20) is used as the analytical solution. In the lower panel, the time is fixed at the value inside the second step (τ = 25) for the numerical solution, and (3.16) with (3.20) and (3.21) is used as the analytical solution. We can confirm that the analytical solution reproduces the following features of the numerical solution. In the upper panel, n k,+ is enhanced with the increase of ∂ 0 θ. In the lower panel, the amplitude of n k,+ is enhanced with the increase of ∂ 0 θ, and it shows the oscillatory behavior with the change of ∂ 0 θ.

Axion-driven pair production
We consider the case where the electric field is zero in the whole time range while the axion acceleration is nonzero, We eliminate the constant gauge field by shifting the momentum:k i → k i . The axion acceleration is set to be constant and positive for simplicity. Furthermore, we restrict the time range such that the axion velocity is kept to be positive. 6 Physically speaking, we consider the local time range where the axion velocity changes linearly with time.
In this case, the pair production is driven by the second term in the right-hand side of (2.33). Since the second term is block-diagonal, we can separately treat the contribution from the lower-frequency modes and that from the higher-frequency modes: α k,sσ = α k,s δ sσ , β k,sσ = β k,s δ sσ . (4. 2) The unitary condition holds true for each mode, |α k,s | 2 + |β k,s | 2 = 1, (4.3) and the distribution function for each mode is given by In the time range where the axion velocity is positive, the higher-frequency modes do not cause the pair production. We thus consider only the contribution from the lower-frequency modes. The equation for (α k,+ , β k,+ ) is given by  In the upper panel, the numerical solution at τ = 0 is shown, and the analytical solution is given by (3.17) with (3.20). In the lower panel, the numerical solution at τ = 25 is shown, and the analytical solution is given by (3.16) with (3.20) and (3.21). The parameters are expressed in the eE = 1 unit. We here note that the distribution function for the standard Schwinger effect n k = 2 exp[−π(m 2 + k 2 ⊥ )/(eE)] takes the tiny value 9.4 × 10 −15 .
We again consider the situation (3.7) where the semiclassical picture is valid on the real time axis. In the case (4.1), X k is given by (4.6) The explicit form of the validity condition is given by See Appendix D for the derivation. We thus focus on the turning point of the frequency in the complex lower-half time plane.
In the case (4.1), the lower-frequency has one turning point i.e., ω k,+ (τ 0 ) = 0, Im τ 0 < 0. Near the turning point, the following identity holds true, The equation (4.5) is thus written as follows We can approximately solve this equation by using the behavior of the frequency near the turning points: ω k,+ ∝ (τ − τ 0 ) 1 2 . For τ > Re τ 0 (i.e., ∂ 0 θ > k), the distribution function is given by n k,+ (τ ) exp 4Im Θ k,+ (τ 0 ) . (4.11) For τ < Re τ 0 , the pair production has not been occurred. We set the initial time as τ ini < Re τ 0 . We here explain why the higher-frequency modes do not cause the pair production. We can investigate the contribution from the higher-frequency modes just by flipping the sign of the axion velocity: ∂ 0 θ → −∂ 0 θ. The occurrence condition of the pair production is thus given by −∂ 0 θ > k. Obviously, this condition cannot be satisfied as long as ∂ 0 θ > 0. It is convenient to express Θ k,+ (τ 0 ) as the integral over ∂ 0 θ, where any real value can be substituted for the lower bound of the integral. The choice of real values does not change the distribution function. The constant axion acceleration was moved outside the integral. The complex integral is exactly solvable, and thus (4.11) is given by The weight factor has no momentum dependence. We can understand the momentum independence in the following way.
Let us see that the axion-driven pair production is analogous with the Schwinger effect. The dispersion relation for the axion-driven pair production is given by 14) and that for the Schwinger effect is given by Each quantity in (4.14) corresponds to a certain quantity in (4.15): k ↔ k 3 , ∂ 0 θ ↔ eA (i.e., ∂ 2 0 θ ↔ −eE), and m 2 ↔ m 2 + k 2 ⊥ . The first terms in the square roots determine the occurrence condition of the pair production; i.e., the pair production occurs when the first terms become zero, and the remaining terms appear as the numerators of the weight factors. In (4.14), all the momentum dependence is assigned to the first term, and thus it does not appear in the distribution function. 7 Before closing this section, we emphasize that the semiclassical approach does not use the exact solution of the Dirac equation, and thus can be applied to the situation where ∂ 2 0 θ and m are not exactly constant (recall that m is not constant in an expanding universe). We have only to assume that the pair production occurs at the short time scale where their time evolutions are negligible.
Specifically, the turning point is determined as follows where ∂ 0 θ(Re τ 0 ) = k. The general form of the distribution function is thus given by .

(4.18)
In the physical scale and the cosmic time t, it is expressed as follows where t k = Re τ 0 dτ a, H phys =ȧ/a, and˙means the derivative with respect to the cosmic time. This distribution function is consistent with that found in [15,16].

Conclusion and Discussion
We studied the pair production of fermions in a time dependent axion background with and without an electric background. For this study, we developed the semiclassical approach such that the adiabatic mode functions incorporate the gauge field and the axion velocity dependence of the dispersion relation. The equation (2.33) summarizes our semiclassical approach, and shows that the presence of the axion velocity gives enhancement and interference effects on the pair production driven by the electric field (axion-assisted pair production), and the presence of the axion acceleration itself causes the pair production (axion-driven pair production).
For the axion-assisted pair production, we derived the analytical formula (3.16) and (3.17) which can describe the distribution function of produced fermions in the parameter region (3.9) for a sizeable axion velocity (3.2). This result justifies and improves the empirical formula proposed in [9].
Also for the axion-driven pair production, we derived the analytical formula (4.13) which can describe the distribution function of produced fermions in the parameter region (4.7). The consistent results were found in the context of the axion inflation [15], and the cosmological collider [16].
There are some applications of our results to phenomenologically interesting cosmological scenarios. The axion-assisted pair production can be important in the scenario of the magnetogenesis by the axion inflation [2][3][4][5][6][7][8]. In the original paper of the axion-assisted pair production [9], the authors discussed that the axion-assisted pair production results in a further reduction of the gauge field production during the axion inflation. The axiondriven pair production can be important, e.g., in the kinetic misalignment mechanism, which was recently proposed by [17,18]. This mechanism provides large axion velocity and axion acceleration in the early Universe, which may result in significant axion-driven pair production.
Recalling that the Schwinger effect puts an upper bound on the magnitude of the static electric field, we conclude that the axion-driven pair production puts an upper bound on the magnitude of the axion acceleration. If an axion couples to a light fermion, it cannot be accelerated faster than of order the fermion mass squared. This may strongly restrict model constructions in some cosmological scenarios, including the kinetic misalignment mechanism [17,18,[20][21][22][23] and the axion inflation [15,24,25].
For quantitative study of the influence of the pair production, we need to evaluate physical quantities like the induced current and the energy-momentum tensor, which include the Bogoliubov coefficients as their integrands. We expect that the usage of our adiabatic basis enables the simple expression of these quantities; i.e., they are diagonalized with respect to the lower and the higher frequency modes. This expectation has been confirmed at least for the induced current, detailed calculations of which will be published elsewhere.
where I is the identity matrix. In this paper, we use the chiral representation where O is the zero matrix and σ i are the Pauli matrices The fifth gamma matrix is defined as follows and it anti-commutes with all the gamma matrices

B Universality of adiabatic mode functions
Fork ≤ ∂ 0 θ, the frequencies (2.11) are given by solutions of the following equation, Note that unlike fork ≥ ∂ 0 θ, we consider only the case where the sign of the last term in the left-hand side is minus. In this case, both solutions satisfy the positivity of ω 2 k,± − m 2 . If the sign is plus, both solutions do not satisfy the positivity and thus we do not consider that case.

D Validity condition of semiclassical picture
In the case (3.1), (X k /ω k,+ ) 2 and |∂ 0 X k /ω 2 k,+ | can be large whenk and ω k,+ in the denominators take their minimum values. Therefore, the validity condition of the semiclassical picture (3.7) is satisfied as long as the reference quantities are much smaller than unity at k 3 = 0 andk = ∂ 0 θ. Atk 3 = 0, the reference quantities are given by Recall that we consider a sizeable axion velocity ∂ 0 θ k ⊥ , m. Atk = ∂ 0 θ, the reference quantities are given by We thus obtain the explicit form of the validity condition, k ⊥ ∂ 0 θ eE, m 2 eE.
(D. 5) In the case (4.1), (X k /ω k,+ ) 2 and |∂ 0 X k /ω 2 k,+ | can be large when ω k,+ in the denominators takes its minimum value. Therefore, the reference quantities have to be much smaller than unity at k = ∂ 0 θ.
At k = ∂ 0 θ, the reference quantities are given by We thus obtain the explicit form of the validity condition,