Exact WKB analysis of the vacuum pair production by time-dependent electric fields

We study the vacuum pair production by a time-dependent strong electric field based on the exact WKB analysis. We identify the generic structure of a Stokes graph for systems with the vacuum pair production and show that the number of produced pairs is given by a product of connection matrices for Stokes segments connecting pairs of turning points. We derive an explicit formula for the number of produced pairs, assuming the semi-classical limit. The obtained formula can be understood as a generalization of the divergent asymptotic series method by Berry, and is consistent with other semi-classical methods such as the worldline instanton method and the steepest descent evaluation of the Bogoliubov coefficients done by Brezin and Izykson. We also use the formula to discuss effects of time-dependence of the applied strong electric field including the interplay between the perturbative multi-photon pair production and non-peturbative Schwinger mechanism, and the dynamically assisted Schwinger mechanism.

1 Introduction Our vacuum is not a vacant space, and there are full of quantum fluctuations popping in and out of existence. When the vacuum is exposed to a strong field, those quantum fluctuations can interact with the field, leading to non-trivial responses of the vacuum. A prominent example of such a response is pair production of charged particles from the vacuum by a strong electric field (for review, see Refs. [1][2][3]). Since the first proposal by Sauter in 1931 [4], the vacuum pair production has been under intensive investigation not only as a fundamental prediction of quantum-field theory but also to understand actual physics processes under extreme conditions (e.g., the early-stage dynamics of heavy-ion collisions [5][6][7][8][9][10][11], the vacuum decay induced by a superheavy nucleus [12][13][14][15][16], the magnetogenesis in the early Universe [17][18][19][20][21][22][23], and possibility of prohibition of anisotropic inflation [24]) as well as some analogous phenomena with different kinds of fields/forces (e.g., the Hawking radiation by a strong gravitational field [25,26], the dynamical Casimir effect by time-varying boundaries [27][28][29][30][31][32][33][34], and axion production by a time-dependent condensate [35][36][37]). The vacuum pair production has never been verified by experiments yet, as it requires an extremely strong electric field of the order of eE cr = m 2 e ∼ 10 29 W/cm 2 , with m e = 511 keV being electron's mass (cf. the current world record is eE ∼ 10 22 W/cm 2 , achieved by HERCULES laser [38]). Nevertheless, the up-coming intense laser facilities such as Extreme Light Infrastructure (ELI) may reach or go beyond the critical field strength eE cr , providing the very first opportunity to observe the vacuum pair production in laboratory experiments [39].
Quantitative features of the vacuum pair production drastically change depending on the time-dependence of an applied strong electric field [40][41][42][43]. If the electric field is fast (i.e., the typical frequency Ω is large), the field interacts with the vacuum fluctuations incoherently. That is, the electric field behaves like a dynamical photon γ to produce particles via perturbative multi-photon processes nγ → e + e − , with n being the number of photons involved. This is an analog of the photo-absorption effect in materials. The production number N e + e − becomes proportional to powers of the coupling constant e as lim Ω→∞ N e + e − ∝ e 2n . (1.1) Note that we have a factor of 2 in the exponent because the production number N e + e − is the square of the amplitude ∝ e n . On the other hand, for small Ω, the field interacts with the fluctuations coherently, rather than incoherently, and produces particles via quantum tunneling. The production number N e + e − becomes purely non-perturbative with respect to the coupling constant e as [4,44,45] lim Ω→0 N e + e − ∝ e −const.×e −1 . (

1.2)
This production mechanism is called the Schwinger mechanism, named after Schwinger who first derived the exponential formula (1.2) in a fully quantum-field theoretical manner [45], and can be understood as an analog of the electrical breakdown of materials. In the above, we have implicitly assumed that the applied electric field is dominated by a single frequency mode Ω. For a field with several frequency modes Ω 1 , Ω 2 , · · · , the production mechanisms with different frequency modes interfere with each other. The interference leads to, for example, substantial enhancement in the production number (the dynamically assisted Schwinger mechanism [46][47][48][49][50]), characteristic momentum signatures in the spectrum [51][52][53][54][55][56][57][58][59][60]62], and spin-dependences [58,61]. Those time-dependent effects are utilized to "optimize" electromagnetic field profiles, so that the signatures of the vacuum pair production become the most manifest in intense laser experiments [63][64][65][66][67][68]. Thus, getting a deeper understanding of the time-dependent effects is important, so as not only to understand actual physics phenomena correctly but also to design a better experimental setup. The time-dependent effects have been studied with various analytical approaches such as semi-classical methods (e.g., the worldline instanton method [69][70][71], the imaginary-time method [72], and the steepest descent evaluation of the Bogoliubov coefficients by Brezin and Izykson [40]), the standard perturbative calculation [42], and the perturbation theory in the Furry picture [55][56][57][58][59]73]. Those approaches have different regimes of applicability. For example, the semi-classical methods are justified only in the semi-classical limit, in which one takes a formal limit of 1, while the standard perturbative calculation covers the opposite regime 1. The perturbation theory in the Furry picture can be applied to arbitrary values of , but it is valid only if one can clearly separate a given electromagnetic field into a strong field and perturbations on top of it. It is, therefore, desirable to develop a novel method for the vacuum pair production in order to get a deeper understanding of the time-dependent effects and to cover a wider parameter regime for the production.
The purpose of this paper is to use the exact Wentzel-Kramers-Brillouin (WKB) analysis as a first step toward developing a novel method for the vacuum pair production by a time-dependent electric field. The exact WKB analysis, which has been developed mainly in mathematics since the pioneering work by Voros in 1983 [74], is a powerful tool to analyze Stokes phenomena of Schrödinger-type differential equations. We apply the exact WKB analysis to derive a production number formula, using the fact that the vacuum pair production can be regarded as a Stokes phenomenon of a given field equation. By clarifying the generic structure of a Stokes graph for systems with the vacuum pair production, we show that the production number N e + e − is given by a product of connection matrices for Stokes segments connecting pairs of turning points. The connection matrices are then evaluated using the semi-classical approximation. We show that our formula gives a systematic improvement of the divergent asymptotic series method proposed by Berry [75] (see also Dingle's book [76]) and agrees with other semi-classical approaches such as the worldline instanton method [69][70][71] and Brezin-Izykson's steepest descent evaluation [40]. We also use the formula to discuss the time-dependent effects including the interplay between the perturbative multi-photon pair production and non-perturbative Schwinger mechanism and their interference effects. This paper is organized as follows: In Sec. 2, we review the basics of the exact WKB analysis. Section 3 is the main section of the paper, in which we apply the exact WKB analysis to formulate the vacuum pair production by a time-dependent electric field. Section 4 is devoted to summary and discussion. In Appendix A, we discuss the Airy equation by means of the exact WKB analysis. We briefly review other computational methods for the vacuum pair production such as the standard perturbation theory and the perturbation theory in the Furry picture in Appendices B and C, respectively. We also present a detailed analysis of a Sauter-type electric field in Appendix D and of Stokes graphs in the dynamically assisted Schwinger mechanism in Appendix E.

Preliminaries: Exact WKB analysis
To be self-contained, we here review the basics of the exact WKB analysis. The exact WKB analysis, originally proposed by Voros [74] and developed, e.g., by Pham and his collaborators [77][78][79][80] and by Aoki, Koike, and Takei [81][82][83], is an extension of the conventional (J)WKB method, proposed independently by Jeffreys, Wentzel, Kramers, and Brillouin [84][85][86][87]. The central idea of the exact WKB analysis is to apply the Borel resummation technique [88] to the conventional WKB method, which does not only make the conventional WKB method mathematically well-defined but also provides a powerful tool to investigate Stokes phenomena of WKB solutions. The exact WKB analysis can also be regarded as a generalization and a mathematically rigorous justification of the divergent asymptotic series method proposed by Berry [75].
The conventional WKB method is a method to analyze differential equations with a small parameter. To be precise, let us consider a second order differential equation of the form, where t ∈ R is a real variable and > 0 is some small parameter (e.g., the Planck constant ). One may understand that the parameter controls fastness/slowness of the potential Q. Indeed, by introducing τ ≡ t/ , one can rewrite the differential equation (2.1) as 0 = ∂ 2 τ + Q( τ ) φ. Thus, the potential Q effectively becomes slow and fast in the limit of → 0 and ∞, respectively. The starting point of the WKB method is to make the so-called WKB ansatz, where t 0 is an arbitrary point on R and the WKB ansatz is normalized as Note that φ + (φ − ) describes an out-going (in-coming) wave propagating in the positive (negative) time direction, and hence we added the subscript ±. By substituting the WKB ansatz (2.2) into the differential equation (2.1), one obtains an equation for Ω as with which one may iteratively solve Ω n as and Ω n = 0 if n = odd (2.7) because the differential equation (2.4) does not depend on the sign of . Now, it may look that the series (2.5) gives a "solution" of the differential equation (2.1). Unfortunately, it turns out that the series (2.5) is not necessarily convergent and therefore is in general illdefined. For example, see Appendix A, in which we explicitly show that the Airy equation gives a factorially divergent series. This problem can be resolved by the exact WKB analysis, as we describe below. The idea of the exact WKB analysis is to apply the Borel resummation technique to make the conventional WKB method well-defined. To apply the Borel resummation technique, we first reexpress the WKB ansatz (2.2) by expanding the exponential factor as where we have used Ω 0 = √ Q and complexified the variable t as t ∈ R → z ∈ C. One can explicitly compute the series coefficients ψ ±,n from the iterative solution Ω n (2.6) and understands that the series ψ ± (2.8) is ill-defined in general because of the non-convergence of Ω (2.5). Now, we apply the Borel resummation technique to make the series ψ ± , or accordingly φ ± , well-defined. Namely, we introduce a Borel transformationψ ± as which is a well-defined object if the series ψ ± is at most factorially divergent, i.e., for any n there exist some constants A, C such that |ψ ±,n | < AC n n!. Then, the theory of the Borel resummation guarantees that if a Borel sum Ψ ± exists (Borel summable), i.e., one can perform the following Laplace transformation in a well-defined manner, becomes a well-defined solution of the original equation (2.1) having the asymptotic expansion (2.8). The Borel sum (2.11) can be analytically continued to arbitrary values of , while the smallness of 1 is assumed in the conventional WKB method. Note that Berry [75] assumes some truncation of the series (2.8) and only considers the leading order factorial divergence of ψ ± to compute (a correspondence of) the Borel transformation (2.9), and then further employs the saddle point method to evaluate the Borel sum (2.10). Such treatments are unneeded in the exact WKB analysis.
The Borel (non-)summability is determined by the singularity structure of the Borel transform (2.9) in the η-plane. Since the Borel transform (2.9) is dependent not only on η but also on the complexified variable z, the singularity structure in the η-plane changes as z varies in general. To proceed, let us assume for simplicity that (i) the potential Q is an analytic function on the entire complex z-plane and (ii) its Taylor expansion around a zero always starts from a linear term as where z t is a root of the potential such that Q(z t ) = 0 and is called a turning point (or, specifically, a simple turning point because it is of order one) in the language of the exact WKB analysis. Since the potential Q reduces to the Airy potential Q Airy = c(z −z t ) around the turning point z ∼ z t , one can investigate the change of the singularity structure in the η-plane by analyzing the Airy equation (see Appendix A). Mathematically, this procedure is justified by the so-called WKB-theoretic transformation [81,82]. It turns out that the singularity can hit the integration contour of the Laplace transformation (2.10) when z crosses the so-called Stokes line C zt : Therefore, the Borel sum exists unless z is on top of a Stokes line. It should be noted that we implicitly assumed for the moment that (iii) there are no Stokes lines degenerated with other Stokes lines emanating from other turning points (such a degenerated Stokes line is called a Stokes segment), which case shall be discussed later. The Borel summability is important not only to make the conventional WKB method well-defined but also to describe Stokes phenomena of WKB solutions. Consider a Borel sum defined at some point in a Stokes region, which is defined as a region in the z-plane that is separated from other regions by some Stokes lines. The Borel sum can be analytically continued to the entire Stokes region, since it does not hit any Stokes line within the region. However, whenever moving to another Stokes region, z must hit a Stokes line, at which the integration contour of the Laplace transformation hits singularities in the η-plane. Then, the Borel sum experiences a sudden jump due to the integration of the singularity. This is the Stokes phenomenon of WKB solutions. To be precise, let us consider two Borel sums, say Φ ±,I and Φ ±,II , defined on two neighboring Stokes regions, I and II, separated by a Stokes line C zt . For an analytic potential Q having the property (2.12), one can explicitly evaluate the discontinuity by carrying out the integration around the corresponding singularity and finds 1 [74,81] (see also Appendix A) where we have introduced and one chooses + (−) sign if crossing the Stokes line C zt counter-clockwise (clockwise) with respect to the turning point when moving from the region I to II. Note that i R by definition (2.13) and its sign does not change unless the Stokes line C zt emanating from z t hits another turning point or a singularity, which are forbidden by the assumptions. The connection formula (2.14) is the basis to quantify the Stokes phenomenon of WKB solutions within the exact WKB analysis: (1) Draw a Stokes graph (i.e., draw turning points and Stokes lines) for a given potential Q; (2) Draw a path that connects a WKB solution defined at some Stokes region and another one defined at a distinct Stokes region; and (3) The Stokes phenomenon of the two WKB solutions is quantified by successively applying the connection formula (2.14) whenever the path hits a Stokes line.
Remind that we have assumed the following to justify the above procedures: (i) the potential Q has no singularities in the complex z-plane; (ii) all the turning points are simple satisfying the condition (2.12); and (iii) there are no Stokes segments. For (i), if the Stokes regions that the path traverses contain some singularities, the connection formula (2.14) receives some corrections. Such corrections are important in, e.g., analyses of a Fuchsian differential equation [83]. Nevertheless, in most physical applications, the potential Q is obtained by analytically continuing some regular function defined on R and thus Q would not have singularities near the real axis. Therefore, very roughly speaking, as long as one considers a path that goes near the real axis (which is actually a convenient choice for the vacuum pair production; see Sec. 3), singularities of Q sufficiently far from the real axis might not matter. On the other hand, one should be careful about the assumptions (ii) and (iii) if the potential Q has some symmetries. For example, potentials having a symmetry exp ∓ i z z t dz Ω(z ) , one should replace σz t with unity in Eq. (2.14). For a different normalization, one also needs to multiply a proper normalization matrix in addition to the connection matrix, when computing Stokes constants with the procedure (3).
[Q(z)] * = Q(z * ) must have a Stokes segment connecting a turning point z t and its conjugate z * t , which is also a turning point (see Sec. 3.3 for a proof). The existence of a Stokes segment brings some difficulties in the exact WKB analysis. For example, in a case of the Weber potential Q Weber (z) = z 2 − c 2 , there appear the so-called fixed singularities in the η-plane due to a Stokes segment that connects z = +c and z = −c [74, 77-80, 82, 89, 90]. The fixed singularities do not change the structure on the z-plane and may cover the whole right-half of the η-plane. This implies that the integration contour of the Laplace transformation hits the singularities no matter what values of z, i.e., the Borel sum never exists. Therefore, to apply the exact WKB analysis safely, one needs to avoid the existence of a Stokes segment. In Sec. 3.4, we consider an infinitesimally small perturbation on top of a potential Q to break symmetries causing a Stokes segment and derive a connection formula for a Stokes segment by taking the vanishing limit of the perturbation after safely applying the above procedures (1)-(3).

Exact WKB analysis of the vacuum pair production
We discuss the vacuum pair production on the basis of the exact WKB analysis. We begin with clarifying our physics setup as well as working assumptions in Sec. 3.1; in a word, we consider a scalar quantum electrodynamics (QED) in the presence of a timedependent electric field. After explaining in Sec. 3.2 that the vacuum pair production can be understood in terms of a Stokes phenomenon of WKB solutions, we perform the exact WKB analysis to derive the production number formula in Secs. 3.3-3.5, which contain the main results of the present paper. To be specific, we first discuss the generic structure of a Stokes graph for the vacuum pair production and show the existence of Stokes segments in Sec. 3.3. In Sec. 3.4, we derive a connection formula for a Stokes segment by considering an infinitesimally small perturbation and assuming the semi-classical limit, in which we take a formal limit of 1. Using the results obtained in Secs. 3.3 and 3.4, we explicitly derive the production number formula in Sec. 3.5 and discuss the time-dependent effects including the interplay between the perturbative multi-photon pair production and the nonperturbative Schwinger mechanism and their interference effects. In Sec. 3.6, we compare our exact WKB result in the semi-classical limit with the worldline instanton method [69][70][71] and Brezin-Izykson's steepest descent evaluation [40] and show that they are equivalent up to unimportant prefactors.

Setup
We consider a scalar QED 2 in the presence of a time-dependent and spatially homogeneous U(1) electric field. The complex scalar field φ in the momentum space satisfies a Klein- where m > 0, p, e ∈ R, and A are mass, (canonical) momentum, the QED coupling constant, and a U(1) gauge potential for the electric field E ≡ −∂ t A in the temporal gauge A 0 = 0, respectively. Note that the U(1) gauge potential A is real-valued on t ∈ R. We can naturally identify the parameter in the exact WKB analysis (2.1) with the Planck constant in our quantum problem. We assume that (i) the electric field is switched off at the infinite future and past, i.e., The assumption (i) shall be used to define particle states in a well-defined manner (see Sec. 3.2). For simplicity, we also assume that (ii) the gauge potential A, after analytically continued to the complex z-plane t ∈ R → z ∈ C, is an analytic function in the entire complex plane, as we assumed in Sec. 2. To quantitatively compute the production number of the vacuum pair production (see Secs. 3.4 and 3.5), we also consider (iii) the semi-classical limit, in which we neglect higher order terms in appearing in connection matrices by formally taking the limit of 1.

The vacuum pair production as a Stokes phenomenon
We explain how the vacuum pair production can be understood in terms of a Stokes phenomenon of WKB solutions. Namely, based on the canonical quantization formalism of quantum-field theory, we explain that the vacuum pair production is reduced to a scattering problem and that the production number is quantified by a Bogoliubov transformation that connects solutions of a field equation at the infinite past and future (see also Refs. [11,91]). In the standard canonical quantization formalism, one may expand the field operator φ in terms of an annihilation operatorâ for a particle and a creation operatorb † for an anti-particle as where ϕ + and ϕ − are mode functions for positive and negative energy states, respectively, satisfying the mode equation (3.1). One may normalize the mode functions ϕ ± in the same manner as the WKB ansatz (2.3) as Combining with the mode expansion (3.3), one findŝ The non-vanishing canonical commutation relations for the creation/annihilation operators are (3.6) and the other commutations are vanishing.
One needs a special care in identifying "positive" and "negative" energy states in the presence of an external electric field (or a time-dependent potential, in general). Indeed, the time-translational invariance is explicitly broken while the electric field is being turned on (which physically means that the field is supplying energy to the system). Thus, energy is no longer a good quantum number and there is no good quantum number to characterize a "particle." This implies that one can define the notion of a particle in a well-defined manner only after the electric field is turned off, for which the time-translational invariance is restored. Suppose that there are asymptotic regions where Q(t) becomes constant (the electric field is turned off), which we assume to occur at the asymptotic times (3.2) [assumption (i)]. Then, one can naturally identify the positive and negative energy states by plane waves in these asymptotic regions, which are eigenfunctions of the time-translation operator (or the energy operator) +i ∂ t . However, we need to distinguish two different boundary conditions, since there are two separate asymptotic regions at t = −∞ and t = +∞. We, therefore, introduce two separate sets of mode functions ϕ ±,as ("as" denotes "in, out") as full solutions of the field equation (3.1) satisfying two different boundary conditions associated to two different asymptotic regions: One has to distinguish ϕ ±,in and ϕ ±,out because the plane waves cannot be solutions of the mode equation (3.1) for |t| < ∞, during which the electric field is being turned on, and they mix up with each other during the time-evolution so that ϕ ±,in = ϕ ±,out . Since we now have two distinct mode functions, we must distinguish the corresponding annihilation/creation operators as well. Noting Eq. (3.5), it is legitimate to definê Physically,â as andb as correspond to annihilation operators of a particle and an antiparticle at the corresponding asymptotic times, respectively. Since ϕ ±,in = ϕ ±,out ,â in ,b in = a out ,b out follows, and the mismatch betweenâ in ,b in andâ out ,b out is given in terms of that between ϕ ±,in and ϕ ±,out . To quantify the mismatch, we notice that ϕ ±,in and ϕ ±,out satisfy the same second order differential equation (3.1) and that a second order differential equation can have only two independent solutions. Therefore, there exists a 2 × 2 matrix Note that 1 = det U because ϕ ±,in and ϕ ±,out satisfy the same normalization condition (3.4). Plugging this expression into Eq. (3.8), we find which is called a Bogoliubov transformation. The Bogoliubov transformation (3.10) is the essence to describe the vacuum pair production. It is evident from Eq. (3.10) that the vacuum states at t = −∞ such that 0 =â in |vac; in =b in |vac; in (3.11) is no longer annihilated by the annihilation operators at t = +∞ as This means that particles are produced from the vacuum |vac; in , and the production number reads where V is the spatial volume and we have used δ 3 (p = 0) = V /(2π ) 3 . Thus, analyzing the vacuum pair production is equivalent to computing the off-diagonal components of the Bogoliubov transformation (3.10). The Bogoliubov transformation (3.10) is nothing but a Stokes phenomenon of WKB solutions. Indeed, Borel sums defined around t = −∞ and +∞, which we write Φ ±,in and Φ ±,out , respectively, give solutions of the mode equation (3.1) that asymptote the planes waves with t → −∞ and t → +∞. Therefore, the Borel sum Φ ±,as surely satisfies the boundary condition (3.7), and we can identify Φ ±,as = ϕ ±,as . (3.14) In general, the Borel sum Φ ±,in cannot be analytically continued to Φ ±,out because of a Stokes phenomenon of WKB solutions. To be precise, consider a path in the complex z-plane that connects t → z = −∞ to +∞ traversing n Stokes regions. We label the Stokes region containing t → z = −∞ as 1 and successively increase the number 2, 3, · · · as entering the neighboring regions until reaching n, which is the final region containing z = +∞. Whenever analytically continuing Φ ±,in from a Stokes region i to i + 1, Φ ±,in experiences a sudden jump, which is quantified by a 2 × 2 matrix T i ∈ SL(2, C) just like Eq. (2.14). Note that Stokes regions are not necessarily separated by Stokes lines emanating from a simple turning point in general; a Stokes line emanating from a turning point of order more than one or a pole, or a Stokes segment can separate Stokes regions as well. In such cases, the connection matrix T i cannot be identified with Eq. (2.14). As we discuss later, the Stokes regions 1, 2, · · · , n are, indeed, separated by Stokes segments in the case of the vacuum pair production [or potentials Q having a property [Q(z)] * = Q(z * ) in general].
For the moment, we do not have to specify the explicit form of T i , but an important point here is that one may symbolically express the relationship between Φ ±,in and Φ ±,out in terms of T i as Comparing this expression with Eq. (3.9), we find Thus, evaluating the off-diagonal components of the Bogoliubov transformation (3.10) is reduced to evaluating the product of the connection matrices T i 's. Note that Eq. (3.16) is an exact relation for exact T i 's. Approximation enters only when one approximates T i 's. Before closing this subsection, we comment on the relationship between our use of the exact WKB analysis and another common use of it, i.e., application to bound-state problems to get exact quantization conditions (e.g., Ref. [92]). In this work, we apply the exact WKB analysis to a scattering problem of quantum fields, i.e., to investigate how the behavior of a wave function changes between t = −∞ and t = +∞, and the change is described in terms of the Bogoliubov transformation (3.9). In scattering problems, the magnitude of the Borel sums Φ ±,as (or the mode functions ϕ ±,as ) is always bounded |Φ ±,as | < ∞ because the potential is positive definite Q(t ∈ R) > 0, and accordingly the system never gets quantized. On the other hand, the exact WKB analysis has also been applied to bound-state problems in quantum mechanical systems. In bound-state problems, the potential Q is a function of space x, instead of time t. A crucial difference is that the potential Q takes negative values at the asymptotic points x = ±∞ in boundstate problems. Thus, either of the Borel sums Φ + and Φ − diverges exponentially at x = ±∞. So as to have a normalizable wave function, the divergent Borel sum at the asymptotic points x = ±∞ must be vanishing. The Borel sums at the asymptotic points x = ±∞ are related with each other by a Bogoliubov transformation as in the case of our scattering problem, and hence the normalization condition requires some conditions onto the Bogoliubov transformation, which eventually give exact quantization conditions of a given quantum mechanical system. During the above procedures in bound-state problems, the exact WKB analysis is applied to get the Bogoliubov transformation. This is completely the same usage as in our scattering problem. Nonetheless, our analyses presented below cannot be applied directly to bound-state problems. For example, the generic properties of Stokes graphs that we shall discuss in Sec. 3.3 [in particular, the properties (4) and (5), which assume Q > 0] should be modified by the property Q < 0. Accordingly, our path (see Fig. 1) to compute the Bogoliubov transformation is not necessarily suitable in bound-state problems. It is an interesting topic to extend our analyses to bound-state problems, and we leave it as a future work.

Generic properties of Stokes graph
As the first step toward the exact WKB analysis of the vacuum pair production, let us discuss some generic properties of a Stokes graph for the potential Q (3.1): (1) The Stokes graph is symmetric in the upper and lower half complex planes, i.e., symmetric with respect to Im z ↔ −Im z.
Proof. Since Q is a real-valued function on the real axis, the Schwarz reflection (2) Let z t be a turning point of the potential Q. The turning point z t is of order Proof. Around a turning point z ∼ z t , the potential Q behaves as which proves the statement.
(3) A turning point z t and its conjugate z * t are always connected by a Stokes segment. Proof. It is sufficient to show 0 = Im[+i (4) The Stokes segment connecting a pair of z t and z * t crosses the real axis only once. Proof. A Stokes segment connecting z t and z * t should cross the real axis because of the topology. Namely, in order to connect two points having opposite signs of the imaginary parts Im z t = −Im z * t via a continuous line, the line must cross Im z = 0 because of the intermediate value theorem. This proves the first half of the statement. To prove the second half, suppose that there are two (or more) crossings t 1 , t 2 ∈ R between the Stokes segment and the real axis. One can assume t 1 = t 2 because Stokes lines cannot cross each other except at a turning point, which is absent on the real axis since Q(t ∈ R) ≥ m > 0, or at a singularity, which is assumed to be absent by the assumption (ii). Then, 0 = Im[i √ Qdz] under the assumption (ii). However, Q ∈ R on the real axis, and thus 0 = Im[i Qdz], which contradicts with the above. Therefore, the Stokes segment can cross the real axis only once.

Im z
Re z Figure 1. (color online) Generic structure of a Stokes graph for the potential Q (3.1). The red points, the blue lines, and the green wavy lines are representing turning points, Stokes lines, and branch cuts, respectively. Stokes segments connecting a pair of (z t,i , z * t,i ) are represented by the doubled blue lines. The orange arrow represents the path we consider to compute the Bogoliubov transformation between Φ ±,in and Φ ±,out . Note that the Stokes segments are depicted by straight lines. This is just for simplification, and a Stokes segment is not necessarily a straight line in general. Also, there could exist other Stokes segments connecting pairs on either half of the complex z-plane (z t,i , z t,j ) or (z * t,i , z * t,j ). Those Stokes segments do not cross the real axis directly (see also footnote 3 for appearance of multiply degenerated Stokes lines) and thus are omitted here for simplicity.
(5) Stokes lines/segments emanating from a turning point z t , other than the Stokes segment connecting z t and z * t , cannot cross the real axis. Proof. The proof is the same as that for the property (4). Suppose there are two (or more) Stokes lines/segments emanating from the same turning point z t crossing the real axis at, say, t 1 , t 2 ∈ R. Then, one gets 0 = Im[i Qdz], which is a contradiction. Therefore, only one Stokes line, which is nothing but the Stokes segment connecting z t and z * t , can cross the real axis.
Summarizing the above properties (1)-(5), a Stokes graph for the potential Q (3.1) should generically look like Fig. 1. In Fig. 1, we have considered p such that 0 = (p − eA(z t )) · eE(z t ), so that any turning points z t,i 's become of order one [property (2)]. In general, three Stokes lines emanate from a simple turning point. Indeed, the integration of Eq. (3.17) yields and thus Stokes lines around z ∼ z t are emanating in directions where k is an integer. The integral (3.18) is a multi-valued function around a turning point. Accordingly, we inserted cuts for each turning point in such a way that the cuts do not traverse the real axis. Then, one understands that k can take three values when restricted in one Riemann sheet, and only one out of the three values corresponds to the Stokes segment that connects z t and its conjugate z * t . Below, we concentrate on a Riemann sheet such that √ Q > 0 on the real axis. Note that one may insert a cut in other directions, which does not alter our final results if one takes an appropriate path for the new cut. Also, one may play with the exact WKB analysis for p such that 0 = (p − eA(z t )) · eE(z t ) by first considering some perturbation onto p → p + δp, for which 0 = (p − eA(z t )) · eE(z t ), and then taking δp → 0.
In order to compute the Bogoliubov transformation (3.10) based on the exact WKB formula (3.16), we consider a path from z = −∞ to +∞ just along the real axis as shown in Fig. 1, so as to avoid the branch cuts located away from the real axis. The path crosses all the Stokes segments that connect pairs of turning points z t,i and z * t,i , while it does not cross any other Stokes lines because they never traverse the real axis [properties (4) and (5)]. As explained in the previous subsections, the Borel sum Φ ±,in experiences a sudden jump whenever it crosses the Stokes segments. The jump is quantified by connection matrices T i 's, whose explicit form in the semi-classical limit is obtained in the next subsection.

Connection formula at a Stokes segment
We derive the connection matrix T for a Stokes segment connecting a pair of turning points z t and z * t , assuming the semi-classical limit. In the exact WKB analysis, the existence of a Stokes segment can spoil the Borel summability even inside of a Stokes region. To circumvent this difficulty, let us recall that the existence of a Stokes segment is closely related to symmetries of the potential [e.g, [Q(z)] * = Q(z * ); see the proof for the property (3)]. This fact temps us to consider some perturbations Q → Q + δQ to break the symmetries, so as to de-degenerate a Stokes segment into some Stokes lines. Then, one can safely apply the exact WKB analysis and expects that the correct connection matrix T is reproduced by taking a vanishing limit of the perturbations after completing the procedures (1)-(3). We show that this prescription works in the semi-classical limit and write down an explicit form of T .
To be specific, let us consider a Stokes segment C zt;z * t into which two Stokes lines, C zt and C z * t , emanating from a pair of turning points, z t and z * t , are degenerated 3 . Without loss of generality, we can set Im z t > 0. Note that Im z t = 0 is forbidden because Q > 0 on the real axis. Now, consider an infinitesimally small perturbation 4 , and denote the corresponding Borel sum as One may interpret the perturbation (3.21) as the i -prescription in quantum-field theory.
because of the property (1)]. Repeating the same argument, one can generalize the connection matrix for a doubly degenerated Stokes segment (3.29) as where we have dropped e −(Sz t,i +Sz t,i )/ contributions by virtue of the semi-classical approximation. Just for the sake of simplicity, we implicitly assume in the main text that all the Stokes segments are doubly degenerated and simply use the connection matrix (3.29), instead of the generalized one (3.20). One can carry out the same analysis for multiply degenerated Stokes segments in the following and can show that our main result for the production number (3.33) is unchanged. 4 In principle, one can consider any perturbation here, as long as it breaks the symmetry [Q(z)] * = Q(z * ) so that Stokes segments are de-degenerated, and it gives the same result in the semi-classical limit. Note that one may formally deform the Planck constant → e ±i0 + to break the symmetry, though it is equivalent to our perturbation (3.21). This is because, as we mentioned below Eq. (2.1), changing values of is equivalent to that of the argument of the potential Q. Indeed, the Klein-Gordon The perturbation Q → Q (±) shifts (the real part of) the turning points z t , z * t as . Equation (3.23) implies that the turning points z t and z * t are shifted to the right (left) and left (right), respectively, in the complex z-plane by the perturbation Q → Q (+) (Q → Q (−) ); see Fig. 2. The shift of the turning points z t and z * t also shifts the Stokes lines C zt and C z * t , which in turn de-degenerate the Stokes segment C zt;z * t as shown in Fig. 2. The two Stokes lines de-degenerated from the Stokes segment should have the following property: As a proof of Eq. (3.24), we first notice that the de-degenerated Stokes line C zt (C z * t ) crosses the real axis only once, which is a reminiscence of the property (4), and is heading to the lower (upper) half plane at the crossing because Im z t > 0. Thus, by noticing √ Q > 0 on the real axis, one understands that dz at the crossing. On the other hand, the integrals in Eq. (3.24) monotonically increase or decrease as z varies from their initial point to the end. This is because if they suddenly start decreasing/increasing at some pointz on the corresponding Stokes line while they were increasing/decreasing until that point,z must be a singularity or another turning point, which is assumed to be absent. Therefore, the integral along C zt (C z * t ) is increasing +i √ Qdz > 0 (decreasing +i √ Qdz < 0) as z varies. This proves Eq. (3.24), as the integrals are vanishing at the initial points. Now, we consider a path that crosses the de-degenerated Stokes lines from the left to right in the complex z-plane (see Fig. 2). The path first crosses the Stokes line C z * t (C zt ) in the clockwise (counter-clockwise) direction around the turning point, and then crosses C zt (C z * t ) in the counter-clockwise (clockwise) direction for Q (+) (Q (−) ). Therefore, by successively using the connection formula (2.14), one gets where Φ (−) ±,R/L denotes the Borel sum defined at the right/left Stokes region with respect to the de-degenerated Stokes lines (see Fig. 2) and (3.27) where we have used the assumption (ii) that there are no singularities in the complex z-plane to get the second line, and the positivity of S zt follows from Eq. (3.24). Equation (3.26) indicates that because of the diagonal components. It means that one cannot obtain T in the naive δ → 0 limit 5 of T (±) . Nevertheless, if one neglects the exponentially small factors O(e −2Sz t / ) in the diagonal components, which are negligible compared to the off-diagonal ones in the semi-classical limit, the limit δ → 0 becomes well-defined and one may obtain T as Below, we use this approximate connection matrix (3.29) to derive the production number formula for the vacuum pair production.

Production number formula in the semi-classical limit
We derive the production number formula within the exact WKB analysis in the semiclassical limit (which we shall call semi-classical exact WKB analysis for brevity). Let z t,i (i = 1, 2, · · · , n) be the i-th turning point in the upper half plane such that Im z t,i > 0 and Re z t,1 ≤ Re z t,2 ≤ · · · ≤ Re z t,n (see Fig. 1). One may identify T i 's in Eq. (3.16) as the connection matrices for the Stokes segments emanating from the above z t,i 's. To proceed, we decompose Eq. (3.29) as where I 2 is a 2 × 2 unit matrix. Plugging this expression (3.30) into Eq. (3.16), we obtain ± becomes discontinuous at δ = 0, which leads to the discontinuity in the connection matrices T (+) and T (−) . In general, the discontinuity can be quantified by the so-called Voros symbol and is studied well for the Weber potential Q Weber = z 2 −c 2 ; see, for example, Refs. [74, 77-80, 82, 89, 90].
where we have used δT = O(e −Sz t / ) and neglected terms of the order of O(e −(Sz t,i +Sz t,i )/ ), so as to be consistent with the semi-classical approximation that we used in Eq. (3.29). The off-diagonal components read Therefore, the phase-space density of the produced particles (3.13) reads The same amount of particles and anti-particles are produced, indicating that a particle and an anti-particle are always created as a pair because of the gauge invariance. The production number (3.33) is basically controlled by the exponentially small factor e −Sz t,i / , and a pair of turning points having the smallest S z t,i dominates the production. When there are several pairs of turning points equally contributing to the production, those turning points interfere with each other because of the factor e −i Im σz t,i / . This point was missing in the naive imaginary time formalism developed by Popov [72] and was first clarified by Dumlu and Dunne [60,62], who demonstrated that the interference is responsible for the characteristic momentum signatures observed in field profiles with subcycle structures [51][52][53][54][55][56][57][58][59]. Equation (3.33) agrees with other semi-classical methods such as the steepest descent evaluation of the Bogoilubov coefficients by Brezin and Izykson [40] and the worldline instaton method [69][70][71], as we discuss in detail in Sec. 3.6. The same formula was also used in particle production with other types of an external field in, e.g., Refs. [93,94], whose derivation is based on Berry's divergent asymptotic series method [75].
3.5.1 Interplay between the non-perturbative Schwinger mechanism and the perturbative multi-photon pair production process The semi-classical exact WKB formula (3.33) describes the interplay between the nonperturbative Schwinger mechanism for a slow electric field and the perturbative multiphoton pair production process for a fast electric field. From a view point of the exact WKB analysis, the interplay can be understood in terms of the change of the location of dominant turning points. We also confirm that the formula (3.33) is valid in the semiclassical regime and cannot describe processes beyond that regime such as the low-order pair production processes (e.g., one-photon pair production). To see those points in an analytical manner, let us assume that the electric field has a typical frequency Ω > 0 (without assuming some specific field profile). Then, the Fourier transformation of the electric field eE, may be peaked sharply at ω ∼ ±Ω. Using eE, one may approximate the gauge potential eA as characterizes the typical electric field strength [for example, a monochromatic electric field eE = eE 0 cos(Ωt) gives eĒ = eE 0 ]. Then, the potential Q in the complex z-plane reads (3.37) For small Ω, one may expand the exponentials in Eq. (3.37) as is kinetic momentum (at time t = 0). We neglected O(|Ωz| 2 )-terms, which is justified if being the so-called Keldysh parameter [40][41][42][43]. Indeed, Eq. (3.38) has only a single pair of turning points z t and z * t given by where P and P ⊥ denote kinetic momenta parallel and perpendicular to the electric field eĒ, respectively. Thus, |Ωz| 2 1 is guaranteed along the integration contour of S zt (3.27), i.e., z : z t → z * t , as long as 1 |Ωz t | ∼ γ is satisfied. Now, using Eqs. (3.38) and (3.41), one can explicitly evaluate S zt (3.27) as where we have used 2 dz √ 1 + z 2 = z √ 1 + z 2 + ln z + √ 1 + z 2 . Therefore, the semiclassical exact WKB formula (3.33) yields Re eĒ , (3.43) which precisely agrees with Schwinger's result for a constant electric field [45]. On the other hand, for large Ω such that γ 1, one may approximate Q (3.37) by expanding the square as The corresponding turning points z t and z * t read Then, we can explicitly evaluate S zt (3.27) as Therefore, which describes the multi-photon pair production involving n = 2 m 2 + p 2 / Ω 1 photons. Equation (3.47) also indicatess that the low-order perturbative processes such as the one-photon pair production (see Appendix B) cannot be described within the semiclassical approximation, for which higher order corrections O(e −2Sz t / ) must be taken into account. Note that the scalar QED is a derivatively coupled theory, so that N e ± ∼ 0 for p ∼ 0.
We demonstrate the interplay and validity of the semi-classical exact WKB formula (3.33) by considering, as an example, the so-called Sauter electric field [4,42], (3.48) An advantage of this field configuration is that one can derive an exact formula for the production number by analytically solving the mode equation (3.1). The exact result, in comparison with the formula (3.33) and the standard perturbation theory at the lowest order (see Appendix B), is shown in Fig. 3 (see Appendix D for more details such as the derivation of the analytical formulas and the Stokes graph). Qualitative features of Fig. 3 can be understood in terms of two dimensionless parameters [42], i.e., the Keldysh parameter γ  For small Ω, where the semi-classical approximation is justified, the formula (3.33) and the semi-classical exact WKB result coincide. Both results agree with the Schwinger formula for a constant electric field (3.43) at around Ω ∼ 0 (or γ 1, ν 1) and deviate from it with increasing Ω (or γ 1, ν 1), implying that the production mechanism becomes dominated by the perturbative multi-photon pair production process. For larger values of Ω (such that γ 1, ν 1), the semi-classical approximation is invalid. In such a parameter region, the low-order perturbative processes such as the one-photon pair production process dominate the production, and the standard perturbation theory works better than the semi-classical approaches.

The dynamically assisted Schwinger mechanism
The production number can be enhanced significantly if one superimposes a weak fast electric field onto a strong slow electric field (the dynamically assisted Schwinger mechanism [46][47][48][49][50]), which is an analog of the Franz-Keldysh effect in semi-conductor physics [57,[95][96][97][98]. The semi-classical exact WKB formula (3.33) gives a good description of this mechanism in the semi-classical parameter regime.
To get an analytical understanding of the above mechanism within the semi-classical exact WKB formula (3.33), let us assume that the strong slow electric field is sufficiently slow and treat it as a constant field. For simplicity, we also assume that the strong and weak fields are pointing in the same direction ∝ e . The gauge potential may be expressed as where we normalized e as |e | 2 = 1 and set eE s > 0 and A(0) = 0 without losing generality. eE s and eE w (such that |eE w | eE s ) denote the strength of the strong and weak fields, respectively. ε is a book-keeping parameter, inserted so as to make sure the weakness of the weak field. Now that eE w is weak, one speculates that the location of a turning point may not change significantly from that for the strong constant electric field alone eE s : where z 2π e +iωt E w (ω) and we used +1 −1 dz √ 1 − z 2 e +ωz = πI 1 (ω)/ω, with I 1 being the Bessel function of the first kind. Equation (3.53) gives a faithful description as long as the semi-classical approximation is valid and Eq. (3.51) is satisfied, i.e., the weak electric field E w is weak enough that the deviation z t − z (0) t is controlled well only by the bookkeeping parameter ε. To get a better understanding of Eq. (3.53), it may be instructive to expand the Bessel function I 1 in terms of ω as In the static limit E w → const., Eq. (3.54) reproduces the Schwinger formula (3.43) for the total electric field E = E s + εE w up to O(ε 2 ). The particle production by E s occurs at t = − p / eE s , and thus the value of E w at that time becomes relevant. Physically, this value corresponds to the instant of time when the energy threshold ∼ m 2 + (p + eE s t) 2 becomes the minimum. From an exact-WKB point of view, it is the location of the crossing between the Stokes segment by the strong constant field E s and the real axis Re z (0) t = − p / eE s , at which the Stokes phenomenon occurs. Equation (3.54) implies that the production number can be enhanced by the time-dependence of the weak field if eE w (− p / eE s ) < 0, compared to the naive prediction of the Schwinger formula (3.43). Even though eE w is weak, the enhancement can be significant if eE s is sub-critical m 2 / eE s 1, which is the essence of the dynamically assisted Schwinger mechanism. Note that the expanded result (3.54) agrees with the perturbation theory in the Furry picture (see Appendix C) if one neglects high frequency corrections O(γ 4 ). The disagreement originates from the O(e −2Sz t / ) terms neglected by the semi-classical approximation.
As an explicit demonstration, let us consider a situation in which a constant strong electric field is superimposed by a weak monochromatic perturbation, i.e., eE = [eE s + eE w cos Ωt] × e , (3.56) where we have assumed eE s > 0 and eE s |eE w |. We have computed the production number by numerically solving the mode equation (3.1), and compared it with the semi-classical WKB formula (3.33) and the perturbation theory in the Furry picture (see Appendix C). The semi-classical WKB formula (3.33) gives a good description for sufficiently small values of Ω, where the semi-classical approximation is valid. For large Ω, it fails to reproduce the exact result quantitatively. Nevertheless, it qualitatively captures the oscillating behavior of the production number. The oscillating behavior is an analog of the Franz-Keldysh oscillation in semi-conductor physics, whose origin is the modification to the energy spectrum by the strong electric field eE s [57]. Even though the semi-classical approximation cannot fully take into account effects of eE w in the quantum regime where Ω is large, it captures those of eE s correctly and thus it shows the qualitative agreement.

Relation to other semi-classical approaches
We discuss the relationship between our semi-classical exact WKB result (3.33) and other semi-classical methods to compute the vacuum pair production; in particular, Brezin-Izykson's steepest descent evaluation of the Bogoliubov coefficients [40] and the worldline instanton method [69][70][71]. Although they look different at first sight, we show that they are equivalent in the semi-classical regime, if one neglects unimportant prefactors.

Brezin-Izykson's steepest descent evaluation
Brezin and Izykson [40] derived an integral representation of the Bogoliubov coefficient U 12 and evaluated it within the steepest descent method, explicitly for a linearly polarized cosine electric field E ∝ cos ωt×e . Note that the Dykhne-Davis-Pechukas formula [99,100] (see also Ref. [101]), which is widely used in analyses of the Landau-Zener transition in the condensed matter community, is obtained in the same manner. For the sake of simplicity, we here focus on U 12 only, but one can equally perform the same calculation for U 21 . We shall see that the steepest descent method and the exact WKB analysis have a close relationship to each other, as already indicated in mathematics [102], and their leading order results in the semi-classical limit coincide.
First, we derive an integral representation of U 12 . From Eq. (3.9), we have Since U 12 is a constant independent of t, one can evaluate Eq. (3.57) at any convenient time. Taking t → +∞, we obtain Note that u 12 (−∞) = 0. To proceed, let us consider to act ∂ t onto u 12 . We thus get where we have used [ 2 ∂ 2 t + Q]ϕ +,in = 0. Integrating u 12 over t, one finds The production number may be suppressed strongly in the semi-classical limit, so one may simply approximate ϕ +,in by the lowest order WKB solution as Note that this approximation is not so accurate and gives an incorrect prefactor [62]; Nevertheless, it is sufficient to extract the exponentially small main factor in U 12 and is the same approximation level considered by Brezin and Izykson [40]. Substituting the approximation (3.61) into Eq. (3.60), one gets an integral representation of U 12 as 6 One may evaluate the integral representation (3.62) with the steepest descent method, which is a good approximation in the semi-classical regime, in which is formally regarded as a small quantity. After complexifing the integration variable t ∈ [−∞, +∞] → z ∈ C, we define saddle points z s,i as Apparently, z s,i is nothing but the turning point z t,i in the language of the exact WKB analysis. We notice that not all the saddle points but only a part of them are relevant in evaluating the integral (this point shall be clarified in more detail later). By picking up the contributions from the relevant saddle points and noticing Eq. can be evaluated as where Γ is an integration contour running from z = −∞ to +∞ during which it wraps each relevant saddle point z s,i over the angle 4π/3. The integrand in the first equality is singular at the saddle points, and we have picked up the corresponding residues to get the second equality. s i ≡ +1 (−1) if the contour Γ z s,i wraps the saddle point z s,i counterclockwise (clockwise). One may expect that the contour wraps a saddle point clockwise (counter-clockwise) if that point is located in the lower (upper) half plane, i.e., Therefore, Essentially, this is Eq. (42) in Brezin-Izykson's paper [40]. Next, we need to identify which saddle points are relevant. For this, we need to understand the topology concerning how steepest descent/ascent lines are located in the complex z-plane. This can be achieved by making use of the properties (1)-(5) for the Stokes graph in the exact WKB analysis. Indeed, in the steepest descent method (or the Lefschetz thimble method in more general), steepest descent/ascent lines C z s,i emanating from a saddle point z s,i are given by a set of points z ∈ C satisfying 0 = Im −2 i z∈Cz s,i z s,i dz Q(z ) , which are nothing but the Stokes lines (2.13) in the language of the exact WKB analysis. Therefore, we can directly use the properties (1)-(5) to determine the structure of steepest descent/ascent lines. In the language of the steepest descent method, one may rephrase the properties (1)-(5) as:    According to the steepest descent method, only saddles whose steepest ascent line crosses the real axis (or the original integration contour of the integral that one wants to evaluate) contribute to the integral. Therefore, either of z s,i and z * s,i is relevant because of the property (3'). Without loss of generality, we can assume Im z s,i > 0. Then, by repeating the same argument below Eq. (3.24), one can show that the "action," −2 i z∈C * z s,i z * s,i dz Q(z ) ∈ R, increases as it goes from z = z * s,i to z s,i along the steepest line C z * s,i . Therefore, the steepest line emanating from z = z * s,i is the steepest ascent, and thus z = z * s,i is relevant. Thus, we have where n is the number of saddle points in the upper half plane, and we have used . Equation (3.67) reproduces the semi-classical exact WKB result (3.32), putting aside the unimportant prefactor 5π/18, which is inaccurate within the approximation (3.61). This coincidence is not an accident, as the Stokes graph in the exact WKB analysis has precisely the same structure as that of steepest ascent/descent lines in the steepest descent method. In obtaining the exact WKB result (3.32), we neglected higher order terms of the order of O(e −2Sz t,i / ) by virtue of the semi-classical approximation. Similar approximations were implicitly used in the steepest descent method in Eqs. (3.61) and (3.64).

Worldline instanton method
The worldline instation method [69][70][71] computes the one-loop QED effective action Γ 1-loop , whose imaginary part gives the production number, based on the Feynman's worldline path integral representation [105,106]. The steepest descent method is applied to evaluate the worldline path integral, and thus the worldline instation method is valid in the semi-classical regime. We shall see that the semi-classical exact WKB result (3.33) is equivalent to the worldline instation method by showing that worldline instanton actions are nothing but cycle integrals enclosing turning points of the potential Q.
The starting point of the worldline instaton method is Schwinger's proper-time representation of the one-loop effective action Γ 1-loop [45]: where τ is the so-called proper-time parameter andD µ is a covariant derivative operator under an external gauge field A µ such that x|D µ |y = δ 4 (x − y) ∂ µ x + i eA µ (x) . We have also assumed an i -prescription m 2 → m 2 −i0 + and used an identity Ô −1 = i ∞ 0 dτ e − i Ô τ for Im O < 0 to get the second equality. One may interpret the last exponential factor as a "time-translation operator" with respect to the proper-time τ with a "Hamiltonian"Ĥ, (3.69) The time-translation operator can be expressed as a path integral. For a spatially homogeneous gauge potential A µ (t, x) = (0, A(t)) [see Eq. (3.1)] 7 , the HamiltonianĤ does not depend on x and one can explicitly carry out the spatial part of the path integration [107]. The resulting expression is the so-called worldline representation of the one-loop effective action Γ 1-loop : where p is canonical momentum, which is an eigenvector of the spatial derivative −i ∂, and we have introduced an "action" S as with Q[z] = m 2 + (p − eA[z]) 2 as before. Note that all the integration variables z, u, and τ are real at this stage. The path-and proper-time τ -integrations in Eq. (3.70) are not analytically doable, except for a few special cases. In the worldline instanton method, those integrations are evaluated by the steepest descent method, which is justified in the semi-classical regime. To be concrete, we first expand the path-integration variable z and the action S as 8 z ≡ z cl + δz, where z cl (or "instanton" on the worldline) is a solution of the classical equation of motion for S, i.e., (3.73) 7 We here assume a spatially homogeneous gauge potential in order to discuss the momentum distribution and to directly see the relationship with our exact WKB analysis at the semi-classical level. In general, the worldline instanton method can be applied to inhomogeneous gauge potentials depending on several coordinate variables [71]. 8 We here carry out the path integral first and then the τ -integration. Alternatively, one may carry out the τ -integration first and then the path-integral. For this case, the resulting classical equation of motion as well as the classical action S cl might look a bit different, but they are essentially the same and one can discuss the coincidence of our semi-classical exact WKB formula in a similar manner that we discuss below.
Here, we complexified the integration variable z ∈ R → C by virtue of the steepest descent method, and thus the classical solution z cl is in general complex-valued. The one-loop effective action (3.68) now reads where the prefactor C comes from the path-integration of the fluctuation δS. The remaining τ -integration may be evaluated by the steepest descent method, which is again valid in the semi-classical regime. Complexifying the variables τ, u ∈ R → C and assuming that the integral is dominated around τ ∼ τ st at which the classical action S cl becomes stationary, we find where the prefactor C is replaced by C after absorbing contributions from the fluctuations around the stationary point Γ dτ τ e + i (S cl (τ )−S cl (τst)) , with Γ being the steepest descent path, and from the time integral dx 0 = du dz cl (0) du = finite [70]. In the second equality of Eq. (3.76), we have assumed C = O(1) for simplicity. In principle, the prefactor C is calculable numerically/analytically [70,108], but it is quite complicated and not essential in our discussion below. In fact, the dominant factor of the production number comes from the exponential e + i S cl , compared to which the prefactor C is not so important. Note that the stationary condition (3.75) is just a necessary condition for a relevant stationary point, which shall be discussed later. So far, we have implicitly assumed that there is only one classical solution z cl with a single relevant stationary point τ st . Generally, one can have several classical solutions z cl,i and/or stationary points τ st,ij . In such a case, the effective action Γ 1-loop shall be given by a sum of all the relevant classical solutions z cl,i and/or stationary points τ st,ij as (3.77) Having obtained the effective action (3.77), we turn to compute the production number. The effective action (3.77) is related to the vacuum persistence probability P as P ≡ | vac; in|vac; out | 2 = e −2 Im Γ ∼ 1 − 2 Im Γ. Since the vacuum decay occurs due to the particle production, it is natural to assume 9 79) where N e − = N e + is assumed because of the gauge invariance (no spontaneous charge production). Then, using Eq. (3.77), we arrive at (3.80) We show that the worldline instanton method (3.80) agrees with our semi-classical exact WKB formula (3.33) (see also Ref. [107]). To this end, we first integrate the classical equation of motion (3.73) over u to find where s ≡ ±1 is fixed after choosing a Riemann sheet associated with the square root in Eq. (3.81). The integration constant a is independent of u but can be dependent on τ . The arbitrariness of a is essentially related to the number of classical solutions. For our steepest descent evaluation (3.76), it is sufficient to know the value of a at τ = τ st , which is uniquely fixed by the stationary condition (3.75) and thus we just have to consider a single classical solution characterized by the unique value a(τ st ). Noticing Equation (3.84) implies that the classical action S cl can have non-trivial values only when a classical path z cl encloses some singularities of the integrand √ Q as shown in Fig. 5. One can always construct such a classical path by properly choosing a path for the complexified time variable u. Topologically distinct paths u give distinct stationary points τ st , whose contributions should be summed up as in Eq. (3.77). In principle, one may consider multiloops and/or loops containing more than two pairs of turning points, whose contributions are exponentially suppressed in the semi-classical regime by the worldline instanton action O(e + 2i S cl ) [i.e., O(e − 2 Sz t )-contributions in the notation of the exact WKB analysis] and are negligible compared to loops containing a single pair of turning points shown in Fig. 5. For such a classical path enclosing only a single pair of turning points z t,i (such that Im z t,i > 0) and z * t,j , we have in which, without loss of generality, the path is assumed to be on the first (second) Riemann sheet dz cl /du = + Q[z cl ] (− Q[z cl ]) when going from z t,i to z * t,j (coming back from z * t,j to z * t,i ). Before proceeding, we remark that the classical path must be in the counter-clockwise direction, so that Im τ st < 0. To show this, we notice that if τ st is a stationary point, its conjugate τ * st must also be a stationary point because 0 = dS cl /dτ | τ =τst Repeating a similar discussion that we presented in Sec. 3.6, one can show that, in the complex τ -plane, the pair of the stationary points τ st and τ * st is connected by steepest descent and ascent lines and that the lines cross the real axis. One can also identify that the steepest line emanating from the point with negative imaginary part is ascent. Thus, a stationary point which has negative imaginary part is relevant, and its conjugate pair, having positive imaginary part, is always irrelevant in the steepest descent evaluation of the τ -integration. The value of τ st is determined by the cycle integral (3.82). For a classical path in the counter-clockwise direction as shown in Fig. 5, it can be evaluated as where we have deformed the integration contour in the second line and x α ∈ R (α = i, j) is a crossing between the real axis and a Stokes line C connecting the pair of turning points z t,α and z * t,α on which Im +i z∈C zt,α dz/ Q[z] = 0. In a similar manner as in Eq. (3.24), one can show +i z∈C zt,α dz/ Q[z] > 0, and thus Eq. (3.86) surely has negative imaginary part, as we wanted. Note that one can show in a similar manner that Re + i S cl,ij < 0 for the counter-clockwise classical path, and thus the production number (3.80) never becomes exponentially large but is always suppressed exponentially by the worldline instanton action. We also remark that classical paths enclosing (z t,i , z t,j ) or (z * t,i , z * t,j ) do not contribute because they never cross the real axis, and hence the boundary condition z cl (0) = z cl (τ ) = x 0 ∈ R cannot be satisfied. Having explained that only counter-clockwise classical paths enclosing a pair of z t,i and z * t,j are relevant at the leading order in the worldline instatnton action, we substitute Eq. (3.85) into Eq. (3.80) to arrive at with t 0 being an arbitrary point on R. Equation (3.87) agrees with our semi-classical exact WKB formula (3.33).

Summary and discussion
We have studied the vacuum pair production by a time-dependent strong electric field on the basis of the exact WKB analysis under the semi-classical approximation. First, we have explained that the vacuum pair production can be formulated in terms of a Bogoliubov transformation, which can be regarded as a connection matrix that describes a Stokes phenomenon of WKB solutions at the asymptotic times t = ±∞. To apply the exact WKB analysis, we have identified the generic structure of a Stokes graph for the vacuum pair production (see Fig. 1), assuming that the potential is adiabatic at the infinite times and is an analytic function on the entire complex plane. Then, we have shown that the total connection matrix is given by a product of that for a Stokes segment connecting a pair of turning points, z t and z * t , which we have evaluated in the semi-classical limit [Eq. (3.30)]. From the connection matrix, we have obtained the production number formula [Eq. (3.33)], which is given by a sum of exponential factors e −Sz t / , controlling the magnitude of the production, as well as imaginary factors e − i Im σz t , responsible for interference between different pairs of turning points. The obtained formula is equivalent to other semi-classical approaches such as the steepest descent evaluation by Brezin and Izykson [40] and the worldline instanton method [69][70][71], and generalizes the divergent asymptotic series method by Berry [75]. The time-dependent effects such as the interplay between the perturbative multiphoton pair production and non-peturbative Schwinger mechanism and the interference effects including the dynamically assisted Schwinger mechanism have also been discussed within the obtained formula, and we have found a good agreement with the exact results in the semi-classical regime.
As a future work, it is desirable to include the higher order O(e −2Sz t / ) corrections, which are important for the production in the quantum regime where cannot be regarded as a small quantity and to describe smooth interplay between the multi-photon pair production processes and the low-order ones. For this, it is essential to improve the connection matrix T for a Stokes segment [Eq. (3.30)]. To the best of our knowledge, there does not exist such a mathematical formula applicable in generic situations. For some special cases, one may transform a generic potential Q into a simple potential that can be analyzed exactly. For example, in the case of the so-called merging-turning-points (MTP) equation, it is rigorously proved that the corresponding potential Q can be locally transformed into the Weber potential [82]. The monodromy structure of the Weber potential (or the resulting parabolic cylinder function) is well-known, so that one can compute T without resorting to any approximations. In a physics sense, this mathematical transformation amounts to mapping a generic electric field configuration into a superposition of constant electric fields. It is interesting to study such special cases as a first step toward a complete formula for the vacuum pair production by a time-dependent electric field.
Another interesting direction is to extend our formulation to the vacuum pair production by other kinds of fields/forces other than a strong electric field or to analogous processes such as the Landau-Zener transition in materials. The formulation presented in Sec. 3 is quite general and is not limited to a strong electric field, as we have just assumed that the potential Q is adiabatic at the infinite times and is an analytic function on the entire complex plane. Our formulation gives a powerful and general framework to discuss the time-dependent effects (or similar effects may occur for spatial variation) in the particle production. One of the most interesting examples is the Hawking radiation. One can use the similar Bogoliubov transformation technique to the Hawking radiation. The problem is then reduced to solving a field equation under a potential determined by the background gravitational field (e.g., for a Klein-Gordon equation under the Kerr-Newman background, the potential is given by that for a confluent Heun equation [112]), and thus our exact WKB analysis can be applied straightforwardly. It is also reassuring that the worldline instanton approach has been applied recently to the Hawking radiation [113], which implies the applicability of our semi-classical WKB analysis because of the equivalence between the two approaches that we have shown in this paper. The spacetime dependent effects, which can be conveniently captured with our framework, could play an important role in the Hawking radiation, e.g., by giving rise to non-thermal corrections.
Note added in proof: While preparing the final draft of our paper, a preprint [114] has appeared in arXiv, which applied the exact WKB analysis to cosmological particle production (particularly for Q(t) ∝ t 2 + gt 4 -type potential) and has some overlap with our paper.

A Exact WKB analysis of the Airy equation
We discuss the exact WKB analysis of the Airy equation. The Airy potential is defined as a potential having one simple turning point at some point z = z t , where ξ, z, z t ∈ C. For the Airy potential (A.1), one can analytically compute the coefficients ψ ±,n in Eq. (2.8). The solution reads Notice that d n is factorially divergent, so is ψ ±,n . In turn, one can explicitly compute the Borel transformationψ ± (2.9) as Therefore, the Borel sum (2.10) reads Since the hypergeometric function 2 F 1 (a, b; c; z) has a branch cut on z ∈ [1, ∞], the integrand becomes singular on Therefore, the integration contour of the Laplace transformation η : 0 → ∞ can hit the singularity if 0 = Im iξ 1/2 (z − z t ) 3/2 = Im i z zt dz Q(z ) , which is nothing but the Stokes line (2.13).
For the Airy potential (A.1), there exist three Stokes lines emanating from the turning point z t . The Stokes lines are straight lines emanating with angle arg(z −z t )+(1/3) arg ξ = −π/3, +π/3, π (mod 2π) and eventually flow into the infinity. For convenience, let us insert a cut on arg(z −z t )+ 1 3 arg ξ = π (mod 2π) and restrict ourselves on the first Riemann sheet such that −π < arg(z − z t ) + 1 3 arg ξ < π. Then, we have three Stokes regions A = I, II, III in the z-plane that are separated by the Stokes lines as . (A.8) In each Stokes region A = I, II, III the Borel sum Φ ±,A is well-defined, but the Borel sums in the different Stokes regions are not necessarily identical Φ ±,A = Φ ±,B if A = B because of the Stokes phenomenon of WKB solutions. Namely, suppose z is initially located in a Stokes region A and smoothly moved to another region B by crossing the Stoke line separating the regions A and B. In the η-plane, the singularity can hit the integration contour of the Laplace transformation as changing z. Then, the Borel sum in the region B gets an additional contribution from the singularity compared to that in the region A as where the second term describes the contribution for the singularity and Γ is a path that wraps the singularity of the Borel transformationψ ± . The second term can be evaluated with the help of 10) and the result reads This is nothing but Eq. (2.14), as Re i z zt dz Q(z ) < 0 on the Stokes line between I and II and > 0 between III and I.

B Standard perturbation theory (at the lowest order)
The standard perturbative treatment with respect to the applied electric field eA gives a faithful description when the field varies very fast compared to the mass scale; or, physically, when the particle production mechanism is dominated by the one-photon (or a fewphoton) pair production process γ → e + e − [42]. In such a regime, the electric field interacts with particles incoherently rather than coherently, for which the semi-classical treatment is not justifiable. Thus, the standard perturbation theory is applicable to different parameter regimes that the semi-classical exact WKB analysis as well as other semi-classical approaches covers, and vice versa.
For later convenience, we here use the Green function technique to derive the production number formula within the standard perturbation theory. We first rewrite the mode equation [i.e., the Klein-Gordon equation (3.1) in terms of the mode function ϕ ±,as ] as Below, we assume smallness of eA, or V , which is justified when the applied electric field is very fast 11 . Indeed, for an electric field with a monochromatic frequency ω, i.e., E = E(ωt), the corresponding gauge potential decreases linearly with ω as |A| ∝ ω −1 . Then, one may solve Eq. (B.1) perturbatively using the Green function technique to obtain where G (0) is a retarded Green function satisfying ϕ ±,(0) is a plane wave solution of 0 = 2 ∂ 2 t + Q 0 ϕ ±,(0) , i.e., ϕ ±,(0) = 1 and we imposed a boundary condition onto ϕ ±,as as [the same as Eq. Using one can express Eq. (B.3) as To be precise, the lowest order perturbation works for ν ≡ ( ω) 2 /| eE| 1 and γ ≡ m ω/| eE| 1 [42]. Notice that the largeness of the Keldysh parameter alone γ 1 is not enough to justify the lowest order perturbation theory.
Taking t → ∞ and comparing with Eq. (3.9), one understands Using the analytical expression for the plane wave ϕ ±,(0) (B.5), one can explicitly evaluate the off-diagonal component of U (B.9b) and arrives at (B.10) where E(ω) ≡ +∞ −∞ dt e −iωt E(t) is the Fourier transformation of E(t), as in the main text. The Fourier transformation E(ω) may be peaked at some characteristic frequency ω ∼ ω c of the applied field E. Then, Eq. (B.10) implies that the production occurs only when ω c ∼ 2Q 1/2 0 = 2 m 2 + p 2 , reflecting the threshold nature of the one-photon pair production.

C Perturbation theory in the Furry picture
The perturbation theory in the Furry picture is an improved version of the standard perturbation theory (see Appendix B) by using a dressed propagator and wavefunction instead of the bare ones G (0) (B.4) and ϕ ±,(0) (B.5). This theory works quite well when there is a clear scale separation in the applied electric field as where E w is some perturbation on top of a strong field E s . In particular, when the strong field E s is sufficiently slow such that it is well approximated by a constant electric field (such a situation may be realized in the dynamically assisted Schwinger mechanism), one can carry out all the calculations exactly, finding an analytical production number formula for E w with arbitrary time-dependence. The derivation of the formula can be done in a parallel manner as the standard perturbation theory (we resort details of the derivation to Ref. [57]). The only difference is that the bare propagator and wavefunction, G (0) and ϕ ±,(0) , respectively, are replaced by the dressed ones G (0) → G (d) and ϕ ±,(0) → ϕ ±,as,(d) under the strong field E s such that 0 = G (d) (t − t < 0), δ(t − t ) = 2 ∂ 2 t + m 2 + (p − eA s (t)) 2 G (d) (t, t ), 0 = 2 ∂ 2 t + m 2 + (p − eA s (t)) 2 ϕ ±,as,(d) (t), where A s is the gauge potential for the strong field E s and we required that ϕ ±,as,(d) asymptotes a plane wave at t → ±∞ as in Eq. (3.7). Treating the remaining field E w as a perturbation and repeating a similar calculation that we explained in Appendix B, one can obtain d 6 N e ± dx 3 dp 3 = 1 (2π ) 3  with A w being the gauge potential for the perturbation E w . One can explicitly evaluate Eq. (C.3) for a constant strong electric field [34,57], which is quite powerful to discuss the dynamically assisted Schwinger mechanism. Namely, we consider E s = E s × e with eE s > 0. (C.5) For this case, one can derive an analytical expression for the dressed wavefunction ϕ ±,as,(d) :

D Analysis of the Sauter electric field
As a supplement to Fig. 3, we here describe the details of the analysis of the Sauter electric field (3.48).

D.3 Standard perturbation theory at the lowest order
The Fourier transformation of the Sauter electric field (3.48) is given by (D.13) Thus, the formula (B.10) is expressed as (2π ) 3 d 6 N e ± dp 3 dx 3 = = ν −4 absent in the above. This implies that the standard perturbation theory is valid outside of the semi-classical regime, or when the applied field is fast enough compared to the electric field strength, such that ν 1.

E Stokes graphs in the dynamically assisted Schwinger mechanism
As a supplement to Sec. 3.5.2, i.e., the analysis of the dynamically assisted Schwinger mechanism with a constant strong electric field superimposed by a weak cosine perturbation (3.56), we here discuss the structure of Stokes graphs. We numerically solved Q(z t ) = 0 and the condition (2.13) to get turning points and Stokes lines, respectively. Figure 8 displays the obtained Stokes graphs for various frequencies of the perturbation Ω.
For small Ω (the leftmost panel in Fig. 8), the particle production is dominated by the pair of turning points closest to the real axis around the origin of the plot. The location of the dominant pair (z t,dom , z * t,dom ) is approximated well by that for a constant electric field (3.52) as (z t,dom , z * t,dom ) ∼ (z of the Stokes graphs changes at intermediate values of Ω (the middle panel of Fig. 8). Namely, we had some Stokes segments passing the infinity before crossing the real axis for small Ω as in the leftmost panel, but those Stokes segments change their topology as increasing Ω and eventually cease to passing the infinity. For example, the four turning points Im z t ∼ 15 in the leftmost panel go down to the real axis as increasing Ω. At some point they cross a Stokes line from the dominant turning point z t,dom , after which the Stokes segments emanating from those turning points cease to passing the infinity. This is intuitively because any Stokes lines cannot cross each other except at a turning point or a pole for finite |z| < ∞, and thus if a Stokes segment eventually crossing the real axis is (not) separated from the real axis by some Stokes lines, it must (needs not) pass the infinity.
For sufficiently large values of Ω (the rightmost panel of Fig. 8), the above-mentioned change in the topology of the Stokes graphs finishes, and all the Stokes segments cross the real axis without passing the infinity. Contributions from each Stokes segment are roughly equal, as the distances between the turning points and the real axis are almost the same and so are the actions S zt 's. The structure of the Stokes graph is essentially the same as that for the cosine perturbation (3.56) alone and is almost unaffected by the presence of the strong constant electric field. This is reasonable since the particle production in this parameter regime is dominated by the perturbative multi-photon pair production processes by the perturbation.