WKB approach to pair creation in spacetime-dependent fields

Besides tunneling in static potential landscapes, for example, the WKB approach is a powerful non-perturbative approximation tool to study particle creation due to time-dependent background fields, such as cosmological particle production or the Sauter-Schwinger effect, i.e., electron-positron pair creation in a strong electric field. However, our understanding of particle creation processes in background fields depending on both space and time is rather incomplete. Here, we propose a generalization of the WKB method to truly spacetime-dependent fields in order to fill this gap.


I. INTRODUCTION
Particle creation out of the vacuum due to extreme external influences is an intriguing effect and a fundamental prediction of quantum field theory. In the following, we focus on electron-positron pair production in quantum electrodynamics. There are several possibilities for pairproducing external fields. For example, in the Sauter-Schwinger effect [1][2][3], particles are created due to a strong electric field. This is even possible for slowly varying electric fields (as long as they are strong enough). Note that this process is different from pair creation in the (perturbative) multiphoton regime which requires sufficiently fast varying electromagnetic fields; see, e.g., [4]. As another example, cosmological pair production [5,6] occurs in an expanding or contracting universe.
So far, electron-positron pair production has been verified experimentally only in the perturbative (multiphoton) regime [4]. Nonperturbative pair production due to an external field is far more difficult to observe in nature and also not nearly as well understood on the theoretical side. Although these effects were first considered more than half a century ago, our understanding of these effects is still far from complete. This is manifest in the fact that there is still very limited knowledge about the influence of the external field's spacetime dependence. Besides numerical simulations (see, e.g., [7][8][9][10][11][12][13][14][15]), several analytical methods have been used for computing the pair production probability, such as the Wentzel-Kramers-Brillouin (WKB) method [16][17][18][19] or the worldline instanton method [20]. However, most of the studies so far were limited to fields that depend on a single coordinate, e.g., time [16][17][18][19][20][21][22][23][24], a spatial coordinate [20,22,25] or a light-cone coordinate [26][27][28][29][30]; see also [31]. Via the worldline instanton method, there have been a few works on truly spacetime-dependent fields, but these were limited to special cases [32][33][34] or a fully numerical treatment; see, e.g., [35] (see also [36] for a work using the Wigner formalism). Regarding the WKB approach, there have been even less studies for background fields depending on both space and time.
In this article, we present a WKB method based on the eikonal (or Hamilton-Jacobi) equation that promises to overcome this fundamental restriction (see [37] for a previous approach to electron propagation based on the eikonal equation). For the sake of simplicity, we consider the Dirac equation in 1 þ 1 dimensions. However, we believe that the method can be generalized to higher dimensions in a straightforward way as long as the external field only depends on time and a single spatial coordinate. As an important example, we study electron-positron pair creation due to a spacetime-dependent mass mðt; xÞ in the Dirac equation. As one possible motivation, we note that a curved spacetime metric such as in cosmological particle production can be mapped to a spacetime-dependent mass in the Dirac equation [38].
The article is organized as follows: We start by reviewing the conventional WKB method for the time-dependent Dirac equation and use a specific time-dependent mass as an example in Sec. II. In Sec. III, we expand solutions of the Dirac equation using solutions of the eikonal (or Hamilton-Jacobi) equation, giving two linear coupled partial differential equations. In Sec. IV we show that these equations reduce to known results if the electric field (or mass) is either purely time dependent or purely space dependent. Problems that occur while solving the eikonal equation with a truly spacetime-dependent field are discussed in Sec. V. The case of a spacetime-dependent mass is considered in Sec. VI. We calculate approximative solutions to the equations mentioned above for a spacetimedependent mass with a weak space dependence in Sec. VII.

II. WKB FORMALISM
Let us start by briefly reviewing the standard derivation of the WKB formalism for purely time-dependent fields in 1 þ 1 dimensions (see, e.g., [18,19] for comparison). As we are interested in pair production due to a spacetimedependent mass (or scalar potential) later on, we consider the case of a time-dependent mass and a time-dependent electric field in 1 þ 1 dimensions.
We start with the covariant Dirac equation where A μ are the components of the electromagnetic potential and γ μ are the gamma matrices satisfying the Clifford algebra's anticommutation relation Now consider the Hamiltonian form of the Dirac equation After expanding ψðt; xÞ into Fourier modes ψ p ðtÞ, we get which are orthonormal, i.e., u † AE u AE ¼ 1 and u † AE u ∓ ¼ 0. As usual, this normalization prescription still leaves the phases of the spinors free to choose. Additionally, one can show that We expand ψ p ðtÞ in terms of these eigenvectors, with the time-dependent phase (eikonal) This expansion (7) reflects the main idea of the WKB approach, i.e., the separation of the rapid oscillation of the phases expfAEiφ p ðtÞg from the slow variation of the background in u AE ðp; tÞ as well as αðp; tÞ and βðp; tÞ. Upon inserting this expansion into (4) and projecting onto u AE ðp; tÞ, we get two coupled ordinary differential equations for αðp; tÞ and βðp; tÞ, We define RðtÞ ¼ βðp; tÞ=αðp; tÞ and find a Riccati equation Note that the exact form of the right-hand side depends on the chosen representation and normalization of u þ and u − due to the factor u † þ _ H p u − . Using γ 0 ¼ σ z and γ 1 ¼ iσ y and assuming vanishing phase difference between the spinors u þ and u − , we find On the other hand, for AðtÞ ¼ 0 we find The number of created positrons (or electrons) with momentum p can be calculated using (see the Appendix A) where β out ðpÞ ¼ βðp; t → ∞Þ, R out ¼ Rðt → ∞Þ and we have used the relation jαj 2 þ jβj 2 ¼ 1 in the last equality. Under the assumption that few pairs are created, i.e., R ≪ 1, a linearized form of the Riccati equation (11) is often used: In that case we get N e þ ðpÞ ∝ jR out j 2 . To obtain R out we integrate the linearized Riccati equation (15) over all times, For symmetric electric fields Að−tÞ ¼ −AðtÞ with a constant mass ( _ m ¼ 0), one expects the maximum number of created pairs for p ¼ 0, as the denominator of the integrand is minimal in that case. On the other hand, in the case with only a time-dependent mass [AðtÞ ¼ 0] the right-hand side of (13) immediately reveals that RðtÞ vanishes for p ¼ 0 and so does the number of produced pairs. Furthermore, upon deforming the integration contour for the integral (16) in the complex plane, we see that the integrand is exponentially suppressed in the lower halfplane. Thus, the integral's value is dominated by the value of the exponential at the singularity closest to the real axis. This singularity at t Ã could be a pole of the prefactor Ω p ðt Ã Þ ¼ 0 or a branch point or any other point where the integrand is not analytic anymore, and thus we cannot deform the integration contour further. Then, R out can be approximated as This estimate does not give the correct prefactor but only the exponent. However, due to the linearization of the Riccati equation, one cannot realistically expect to obtain the prefactor from the integral (16) exactly anyway. If there are multiple singularities that are comparably close to the real axis, contributions from all singularities have to be taken into account, which leads to interference effects in the momentum spectrum [11,19,39].
A. Example: Time-dependent mass As an example, we want to calculate the number of produced pairs for a specific time-dependent mass as a toy model. We use a similar functional dependence later in Sec. VII as a spacetime-dependent mass where some of the results derived here will be useful.
We use a mass of the form with fðτÞ ¼ sechτ and a dimensionless parameter γ that controls the amplitude of the pulse. For large γ, the relative change of mðtÞ is small, and we may use perturbation theory to estimate the pair-creation probability (see below). For small γ, however, the change is large, and we need another method, such as the WKB approach.
This parameter γ also controls the adiabaticity, i.e., the applicability of the WKB approximation. A measure for the adiabaticity is the rate of change _ m of the mass compared to the mass itself, i.e., _ m=m 2 , which scales with γω=m 0 . Thus, the interesting region of small γ can be treated via the WKB approach provided that ω ≪ m 0 .
The expression (18) is motivated by the fact that typically the squares of mass (or potential) terms are added. As an example, let us consider the Dirac equation in 2 þ 1 dimensions where the second spatial dimension is compactified, giving rise to a discrete Kaluza-Klein tower of transversal momenta k ⊥ . Then, the effective masses of the 1 þ 1-dimensional Dirac equations would be m 2 1D ¼ k 2 ⊥ þ m 2 2D . As another example, let us consider a scalar field ð□ þ m 2 Þϕ þ V 0 ðϕÞ ¼ 0 with the interaction potential VðϕÞ. Then, linearization ϕ ¼ ϕ 0 þ δϕ around a given background solution ϕ 0 yields the effective mass m 2 eff ¼ V 00 ðϕ 0 Þ þ m 2 for the perturbation δϕ.
The parameter γ plays a role very analogous to the Keldysh parameter γ ¼ mω=ðqEÞ for strong electric fields; see, e.g., [40]. This analogy can be made even more explicit by considering the form of the effective mass of an electron within a laser pulse (see [41][42][43]), even though the resulting pair-creation probability should be derived by using the vector potentials A μ directly. We then can approximate R out using the linearized Riccati equation (15), where For fðτÞ ¼ sechτ the phase integral can be calculated analytically, giving where The integral for R out is dominated by the value of the exponent at the pole where fðτ Ã Þ ¼ AEiγ, For fðτÞ ¼ sechτ we find and thus This result does not depend on γ which at first is a bit surprising. For example, in the limit γ → ∞, mðtÞ ¼ m 0 ¼ const., and thus no pairs should be produced. This apparent inconsistency can be resolved by the observation that our WKB approach breaks down for large enough γ, where we should use perturbation theory instead (see above).
To confirm our result we computed R out numerically from the full Riccati equation (13). Due to the highly oscillatory coefficients in the Riccati equation, we integrated the equation using the TIDES library [44] in conjunction with the arbitrary-precision library MPFR [45]. To parallelize the computation, GNU Parallel [46] has been used. Figure 1 shows the analytical result from (26) and the numerical result for jR out j 2 together for a specific choice of γ and ω. Because the approximation in (26) does not produce the correct prefactor, we assume it to be ap 2 . This is motivated by the form of the integrand's prefactor in (16) which for E ¼ 0 is proportional to the canonical momentum p. The constant a is then chosen to fit the numerical data.
We find very good agreement between the analytical estimate and the numerical calculation. Even without the heuristically determined factor a, the analytic approximation lies within an order of magnitude of the numerical result.
Indeed, if one plots the values of the numerical results' peaks over different values of ω, the points fall nicely on the curve predicted by the maximum of the exponential in (26) (see Fig. 2). On the other hand, if we fix ω and vary γ, the maximum of the numerical data behaves as in Fig. 3. For small γ ≪ 1 the maximum remains constant, while for large γ ≫ 1 the maximum seems to go like γ −4 . This behavior is due to the prefactor in (19) (18) where fðτÞ ¼ sechτ with ω ¼ 0.1m 0 and γ ¼ 0.1. The squares are numerically calculated results, while the solid line represents the analytical estimate from (26). We used ap 2 as the prefactor of the analytical result with a ¼ 5=m 2 0 chosen to fit the height of the peaks in the numerical result; see (16).
for large γ ≫ 1. In between these two regions the value of the maximum fluctuates. This can be attributed to the prefactor as well because the order of magnitude does not change as one would expect if this behavior came from the exponent.

III. EIKONAL FORMALISM
We now want to develop a more general procedure for calculating the pair production probability that, in principle, also works for spacetime-dependent fields. The main idea of the WKB formalism as presented in the last section is to separate fast and slow oscillations in the wave function: The factor of exp½AEiφ p ðtÞ contains the fast oscillations, while the prefactors α and β contain the slow oscillations. We try a similar approach for spacetime-dependent fields.
First, we define two operators with S AE being the two independent solutions of the relativistic eikonal (or Hamilton-Jacobi) equation, The eikonal equation above can be obtained from classical electrodynamics. Thus, it could be derived from the Dirac equation (1) via inserting the WKB ansatz ψ ∼ expðiS AE =ℏÞ and only keeping the lowest-order terms in ℏ. However, here we motivate the WKB expansion by assuming that the mass m is the largest relevant scale in our problem, leading to rapid oscillations of expðiS AE =ℏÞ.
We use the convention that S þ and S − correspond to solutions with positive and negative energy, respectively, Note that this eikonal equation is an immediate generalization of (8). When A μ and m are constant, the solutions S AE correspond to plane wave solutions, that is, S AE ¼∓ p μ x μ . Squaring the operators M AE , we get Thus, the operators M AE both have the two distinct eigenvalues AEm. Let u AE and v AE be their respective eigenvectors defined as follows: Because the operators M AE are self-adjoint in the sense that whereū AE ¼ u † AE γ 0 and analogously forv AE . We normalize the eigenvectors as follows: Although parts of the following derivation can be carried out in a general manner, we want to focus on the case of a 1 þ 1-dimensional spacetime. Then it is sufficient to use 2 × 2 matrices for the gamma matrices, and the M AE will only have one eigenvector each for every eigenvalue. We expand the spinor ψ in terms of these eigenvectors, which is motivated by the expansion (7) of the spinor in the time-dependent case. There, the functions α and β are the Bogoliubov coefficients of the transformation between in-and out-states (see Appendix A), and therefore we sometimes refer to them as Bogoliubov coefficients here as well. Using the expansion (34), the Dirac equation (1) reduces to In terms of the large-m or small-ℏ expansion mentioned after the eikonal equation (28), the leading-order contribution gives Eq. (28) for the exponent S AE , while the subleading order determines the above equation for the Bogoliubov coefficients, compare Eq. (2) in [47]. Multiplying (35) byū þ orv þ from the left, we get two coupled partial differential equations Analogous to the Dirac convention, we choose for the gamma matrices. Thus, and the eigenvectors u AE and v AE can be written as with the normalization constants After calculating all the inner products that appear in (36), we get the following equations for α and β: where Equations (41) are completely equivalent to the Dirac equation (1) but might offer advantages for numerical simulations and for analytical approximations (see below). For the purely time-dependent case, it is known that solving the quantum kinetic equations (see, e.g., [10,36]) or the Riccati equation (see, e.g., [16][17][18][19]) can be more efficient numerically than the original Dirac equation. Thus, we expect that similar advantages could apply here, especially in cases where the functions S AE are available analytically (e.g., within suitable approximations) or can be efficiently implemented numerically.

IV. KNOWN LIMITING CASES
We now want to show that Eqs. (41) reproduce the correct results for both a time-dependent electric field with a time-dependent mass and a space-dependent electric field.

A. Time-dependent electric field and mass
We use the temporal gauge where Then the two independent solutions of the eikonal equation (28) are given by with φ p ðtÞ as in Sec. II. We thus find None of the coefficients in (41) depends on x in this case. Thus, if we impose boundary conditions such that α and β are constant initially (i.e., for t → −∞) then ∂ x α ¼ ∂ x β ¼ 0 for all times. Equations (41) for α and β then simplify to We define the ratio RðtÞ ¼ βðtÞ=αðtÞ and, using (47), calculate its time derivative which is a Riccati equation that is up to a factor of i (that can be attributed to a different normalization for the spinors u þ and v þ used here than for the spinors u AE in Sec. II) identical to the one in ordinary time-dependent WKB [compare (11)].

B. Space-dependent electric field
For a purely space-dependent electric field (compare [22,25,48]) we use the gauge In complete analogy to the time-dependent case, we find Thus Again, the coefficients in the equations for α and β in (41) are solely space dependent, and by requiring that α and β are constant left of the barrier, i.e., for x → −∞, we find that ∂ t α ¼ ∂ t β ¼ 0 for all values of x. Then, after introducing the ratio R ¼ β=α we again find a Riccati equation, This case is related to the one-dimensional Schrödinger scattering problem from nonrelativistic quantum mechanics. Again the WKB expansion (34) is motivated by the separation of the rapidly oscillating phase expfiS AE g from the rest, which is slowly varying. This assumes that the local momentum scale P ω ðxÞ is much larger than all other relevant scales, such as P 2 ω ðxÞ ≫ jP 0 ω ðxÞj or equivalently ð∂ x S AE Þ 2 ≫ j∂ 2 x S AE j. Of course, this assumption breaks down at the classical turning points where P ω ðxÞ ¼ 0. Even though the above Riccati equation is, in principle, exact, integrating it becomes problematic at those points. Note that, in contrast to the purely time-dependent case, these classical turning points x can be real for sub-barrier tunneling problems. For quantum reflection above the barrier, they are again complex.

V. CAUSTICS
If we consider a truly spacetime-dependent problem, difficulties in solving the eikonal equation (28) may occur. Due to the nonlinear nature of the eikonal equation, it may not be possible to find global solutions in a classical sense; i.e., a solution might not be differentiable everywhere. Note that these singularities of the eikonal equation (28) do not (necessarily) imply that the solutions of the original Dirac equation (1) become singular. They just indicate that the lowest-order WKB approach (34) employed here breaks down. This is very similar to caustics in geometric (ray) optics-e.g., the rainbow effect-where the density of light rays shows a singularity while the full solution of the wave equation remains perfectly regular. Another example is the one-dimensional stationary Schrödinger scattering problem (discussed above) where the WKB approach breaks down at the classical turning points (indicating the onset of tunneling) while the solutions to the original Schrödinger equation remain perfectly regular.
We use the method of characteristics to visualize such situations (see, e.g., [49] or many other standard textbooks on partial differential equations for more details). Using this method any first-order partial differential equation can be cast as a system of ordinary differential equations by finding certain characteristic curves along which the solution of the partial differential equation can be integrated easily. Afterwards, the solutions along many of those curves can be combined into a solution surface. This essentially amounts to going over to another set of coordinates where one coordinate is the parameter to move along the curve and the other coordinates number the curves.
Difficulties appear where two characteristic curves intersect. At such a point the solution is not uniquely defined as we might use the value on either one of the intersecting characteristic curves. Many of these points form a caustic surface.
For example, Fig. 4 shows the spacetime-dependent mass given in (84) together with the (numerically calculated) characteristic curves. We see that such a spacetimedependent mass has a focusing or defocusing effect on the characteristic curves similar to optical lenses on light rays. Indeed we can estimate that the onset of the caustic surface is at time FIG. 4. Projected characteristic curves (dotted) for the mðt; xÞ from (84) together with mðt; xÞ itself (contour) using for p ¼ 0 and m only weakly space dependent; see Appendix B for details. Diagrams like Fig. 4 are well known from geometrical optics. In fact, geometrical optics is just an approximation to wave optics based on the eikonal equation (for optics). That is why it is not too surprising that the above formula (55) for p ¼ 0 is strikingly similar to the formula for the focal length of a thin, biconvex spherical lens [50], where L ∼ 1=ðεωÞ, D ∼ 1=ω and n 2 − n 1 ∼ 1=γ 2 . We see that when the spatial inhomogeneity is weak (i.e., ε is small), the caustics occur far away from the spacetime region in which the mass is nonconstant, i.e., where pairs are produced. Thus, in the case of a purely time-dependent problem no caustics occur and our solution is differentiable everywhere [compare (44)].
In conclusion, assuming that the spacetime region of particle creation is sufficiently localized and that the spatial dependence is weak enough (compared to the temporal variation), the potential problem of caustics (indicating singular solutions of the eikonal equation) occurs far away from the spacetime region where the particles are created and thus does not invalidate our analysis. To cast this statement in a more formal form, there are two options: One option is to choose a finite final time t out which is large enough such that it occurs after all pair-creation processes have taken place, but small enough such that it still occurs before any caustics appear. As another option, one could apply a mild deformation of the mass function mðt; xÞ in this time window which is so slow that the generated pair creation (i.e., mixing of positive and negative frequencies) can be neglected, but it undoes the focusing or defocusing effects and thus avoids caustics.

VI. SPACETIME-DEPENDENT MASS
We now turn to a truly spacetime-dependent problem, namely that of a spacetime-dependent mass mðt; xÞ with no electromagnetic potential, i.e., A μ ¼ 0. This case occurs in a 1 þ 1-dimensional spacetime with curvature: Every 1 þ 1-dimensional spacetime is conformally flat; i.e., its metric can be written as Writing down the Dirac equation in such a spacetime reveals that it is equivalent to the Dirac equation in flat spacetime but with a spacetime-dependent mass mðt; xÞ ¼ ℧ðt; xÞm 0 (see, e.g., [38] for details).
In that case, the eikonal equation (28) is considerably simpler: We may write the two independent solutions S þ and S − using two different functions R and S by splitting S AE into a symmetric and an antisymmetric part, The inverse transformation is given by When the mass is constant, the solutions S AE ¼ ∓ϵ p t þ px correspond to plane-wave solutions. In that case, R ¼ px and S ¼ −ϵ p t. Thus, the case p ¼ 0 is singular in the sense that R vanishes identically. We avoid this case as this leads to problems when using R and S as coordinate transformations (see following subsection).
Using the above definition (60) of R and S in the eikonal equation and computing the sum and difference of the two equations, we find We solve the latter equation for ∂ t R and obtain where we have introduced the abbreviation λ which will be more convenient later on. Inserting this into the first equation in (61) we get Finally, the coefficients in the equations for α and β are

A. Coordinate transformation
Somewhat similar to the method of characteristics mentioned in the previous section, we want to introduce new coordinates which simplify the evolution equations (41) for the Bogoliubov coefficients. The rapidly oscillating exponential contains the difference of the phases S ¼ ðS þ − S − Þ=2, and hence we choose one coordinate (the new time coordinate) in this direction. In order to have the same dimension as time, we define the new time coordinate s via s ¼ Sðt; xÞ=m 0 where m 0 ¼ lim t→−∞ m is the asymptotic value of the mass. To simplify scalar products, the new spatial coordinate r should be locally orthogonal to s. Inspecting the equations above, we find that this is automatically satisfied if we define r ¼ Rðt; xÞ=m 0 .
Then we have and thus The components of the inverse metric tensor in r − s coordinates are then given by where we see explicitly that the coordinates r and s are indeed locally orthogonal. Finally, the components of the Levi-Civita tensor are Additionally, we need to introduce the covariant derivative Furthermore, we rescale α and β according to Again, we assume nonvanishing p ≠ 0 as this would be singular for p ¼ 0 because λ ∝ p. Finally, after several manipulations and simplifications, we get, as equations for α and β, with the abbreviations Equations (71) are still exact, but they have several advantages in comparison to the original Dirac equation (1). First, as in the purely time-dependent case, the rapidly oscillating phase e AE2iS is a function of the new time coordinate s only. Thus, they might also be advantageous for numerical simulations, especially when the transformation from ðx; tÞ to ðr; sÞ coordinates can be implemented efficiently. Second, if λ is small enough (see below) such that we may neglect terms of order λ 2 , Eqs. (71) can be approximated by Third, in the relevant case of α ≫ β, we see thatα does not evolve with time s but stays nearly constant,α ¼αðrÞ, which fits the picture of the characteristics. This suggests a wave packet αðrÞe im 0 rþim 0 s moving along curves of constant r (i.e., in the s direction) whose shape is given by αðrÞ.
Going back to Cartesian coordinates t and x, this corresponds to a wave packet traveling at varying speeds, with the form of the wave packet changing over time (i.e., becoming wider or narrower). Then, we may solve the evolution equation for β by integrating over s for fixed values of r. For each value of r, we then have the same situation as in the purely time-dependent case; i.e., the paircreation exponent will be determined by the complex value of S at the first relevant singularity in the complex s plane. Note that this requires rewriting all functions of t and x as functions of s and r. Then, for all fixed (real) values of r, one should analytically continue in s and find the singularities in the complex s plane. Since this procedure can only be applied fully analytically to special cases, we develop a suitable approximation scheme based on weak spatial dependencies in the following.
Another approach that could be considered is an inverse one (see also [51]): If solutions R and S are given, one can calculate the associated mass m from Eqs. (61). These solutions could be obtained by choosing S such that R can be calculated easily from (61).

VII. WEAKLY SPACE-DEPENDENT MASS
Consider a spacetime-dependent mass where the space dependence is much weaker than the time dependence, i.e., m ¼ mðt; εxÞ ¼ mðt; ξÞ with ε ≪ 1. As before, we use the initial condition S AE ðt in → −∞; xÞ ¼ px ¼ pξ=ε. We can then expand the solutions of the eikonal equation (61) in a power series for small ε, where R n and S n , n ¼ 0; 1; 2; …, are functions of t and ξ.
Because only squares of the derivatives of R and S appear in (61), every second term in the expansions of R and S vanishes, i.e., R 2nþ1 ¼ S 2nþ1 ¼ 0, n ¼ 0; 1; 2; …. To lowest order, we find These are exactly the same expressions as in the purely time-dependent case, with the only change being that the mass m now also depends on x (or ξ). The next nonvanishing terms are given by To simplify this further we assume that p ¼ Oðε 2 Þ, i.e., p ¼ ε 2p wherep ¼ Oð1Þ. Note that our WKB approximation is based on the assumption that the temporal oscillations of expfiS AE g are fast (of order m), and the spatial variation (and thus the momentum p) can be small. In fact, pair creation is expected to be suppressed for large momenta p. Inserting p ¼ Oðε 2 Þ, we obtain Using this approximation in (74) we get It should be noted that in the strict sense the square root ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m 2 þ p 2 p should be expanded in a power series in ε as well. However, we assume that keeping this expression as it is will only enhance the accuracy of our approximation. Inserting these expansions into the definition of λ, we find Hence, if we only keep terms up to order ε 2 in (71), the equations for α and β are where λ 0 ¼ −p= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m 2 þ p 2 p is the leading-order term of λ. Again assuming the dominance of the positive frequency part α ≫ β (i.e., that only a few pairs are created), we find α ≈αðrÞ. Thenβ out can be obtained from the second equation in (80) by integrating over all s. While performing that integral the other coordinate r ¼ px=m 0 þ Oðε 3 Þ has to be held constant. Fortunately, if we only keep terms up to order ε 2 , holding r constant is the same as holding x constant.
The integral is dominated by the singularity closest to the imaginary axis at s Ã ¼ Sðt Ã ; xÞ. Typically, this will occur where λ 0 diverges, i.e., where Here we assume that the function m 2 ðt Ã ; xÞ itself does not possess singularities which are even closer to the real axis. (This could be the case for dynamically assisted pair creation; see, e.g., [18,39,52].) In this case, they would determine t Ã . Thus, we expect β out to behave like β out ðxÞ ∝ e 2iSðt Ã ;xÞ : The density of produced pairs will then be [see (A18)] jβ out ðxÞj 2 ∝ e −2ℑSðt Ã ;xÞ : ð83Þ In complete analogy to the purely time-dependent case, we do not expect this method to yield the correct prefactor due to the approximations made.

A. Example: Hyperbolic secant pulse
As an example, for an only weakly space-dependent mass, we consider which is similar to (18) but with an additional spacedependent function gðχÞ. In complete analogy to Eq. (18), we assume ω ≪ m 0 and sufficiently small γ in order to be in the WKB regime; the limit of large γ corresponds to the perturbative regime.
We again use fðτÞ ¼ sechτ. Solutions to (81) here are the same as in the time-dependent case (25), with the only difference being that nowγ depends on ξ (or, equivalently, x). Comparing S 0 from (75) with φ p from the time-dependent case (20), we see that they are equal up to an overall sign, i.e., S 0 ¼ −φ p , and thus the lowest-order contribution to the exponent of the number of produced pairs is exactly the same as in the time-dependent case. For the next-order contributions we have to calculate with the dimensionless function h depending onγ only, Note that becauseγ ¼γðξ; pÞ this still depends on the momentum p and the spatial coordinate ξ. This integral cannot be solved exactly in terms of elementary functions, but we may obtain the asymptotics. If we expand hðγÞ in a series for smallγ, we find For largeγ ≫ 1, the integrand (89) decays with 4=γ −3 . Note, however, that this limit corresponds to the perturbative regime, where the WKB eventually breaks down.
To test this behavior, we calculate the function hðγÞ numerically (see Fig. 5).
Because hðγÞ > 0 for all values ofγ, the next-order contribution always decreases the pair-creation exponent; i.e., its absolute value increases, thus reducing the number of produced pairs. This is qualitatively consistent with the numerical results from [7] using the worldline formalism. There it was found that the locally constant field approximation overestimates the true pair production probability, at least in the case of a Sauter potential.
Consequently, we see that, to this order of approximation, the density of produced pairs will be at its maximum where g 0 ðωξÞ vanishes. Thus, both minima and maxima of the pulse may give significant contributions to the number FIG. 5. Log-log plot of the numerically calculated function hðγÞ as given in (89). For small γ ≪ 1, hðγÞ approaches π 2 =γ, whereas for large γ ≫ 1 it behaves like 4=γ 3 . of produced pairs (compare [7]) as both are saddle points of the spatial integral in (A18). However, the exact contribution depends on the prefactor in β out ðxÞ which we have not calculated here.
Qualitatively, the momentum dependence of the total number NðpÞ of produced pairs will be the same as in the purely time-dependent case, i.e., quadratically NðpÞ ∼ p 2 for small p and exponentially suppressed for large p. The main effect of the spatial dependence of mðt; xÞ will be an overall reduction of the total amount of NðpÞ, due to the reduced pair-creation volume or length and the correction (88) to the exponent.

B. Higher momenta
After Eq. (76), we used the low-momentum approximation p ¼ Oðε 2 Þ in order to simplify the subsequent expressions. This was sufficient for calculating the lowest-order correction (88) to the pair-creation exponent which shows that the spatial dependence tends to decrease the pair-creation probability. However, as the mass varies on length scales on the order Oð1=εÞ, one might expect that further interesting effects occur on momentum scales of the order p ¼ OðεÞ. Thus, let us briefly discuss this case. According to Eq. (76), R 2 can no longer be neglected, which implies that the coordinates R and x are no longer equivalent. This complicates the analytical continuation because fixed and real values of R do not correspond to fixed and real values of x (for complex t). Furthermore, λ is now less suppressed, λ ¼ OðεÞ, which implies that reaching the desired accuracy of Oðε 2 Þ, one should keep the quadratic terms Oðλ 2 Þ in the evolution equations (71), which adds further complications. Of course, these more complicated equations can also be solved within a consistent expansion in ε, but the resulting expressions will be much more involved than those presented here.
For very large momenta p, on the other hand, one would expect that the results simplify again because the locally homogeneous field approximation along the particle's worldline should become a good approximation.

VIII. CONCLUSIONS & OUTLOOK
Calculating the creation of particle pairs by truly spacetime-dependent external fields (such as gravitational or electromagnetic fields) in the nonperturbative regime is a very challenging task. For purely time-dependent fields, a very powerful method to estimate the pair-creation exponent is the WKB approximation. In this work, we propose a generalization of this approach to truly spacetimedependent background fields, which is based on solutions of the relativistic eikonal equation (28). For fields that only depend on either time or a spatial coordinate, our method reproduces the known results (see Sec. IV).
One of the first obstacles we encounter is the problem of caustics. They indicate that the eikonal equation (28) in truly spacetime-dependent background fields does not have globally differentiable solutions in general, in contrast to the purely time-dependent case. However, if the spatial dependence is sufficiently weak compared to the temporal variation of the background, these caustics are well separated from the spacetime region of particle creation and thus do not spoil our approach (see Sec. V).
Then, via a transformation to adapted coordinates r and s, the Dirac equation in the presence of a spacetimedependent mass mðt; xÞ can be mapped exactly to Eqs. (71) for the Bogoliubov coefficients. These equations have several advantages and could also be suitable for improved numerical simulation schemes. In the low-momentum approximation λ ≪ 1, they simplify to (73). Then, via the usual assumption that the positive frequency part dominates, α ≫ β, we may estimate the Bogoliubov coefficient β associated with pair creation via a simple integral over the new time coordinate s in complete analogy to the purely time-dependent case. Thus, as in the purely time-dependent case, the pair-creation exponent is determined by the first singularity in the complex s plane.
Finally, consistent with our assumption to avoid caustics, we consider the case in which the spatial dependence is much weaker than the temporal variation and employ an expansion in terms of the relative strength ε of the spatial dependence in Sec. VII. To leading order, we obtain a result which is analogous to the locally constant field approximation: At each point x in space, we simply have to integrate the evolution equation for βðt; xÞ over time-in complete analogy to the purely time-dependent case (as if we had a spatially homogeneous background). In analogy to the locally constant field approximation, this leading order could be referred to as a locally homogeneous field approximation.
Calculating the next-to-leading order correction (88) to the pair-creation exponent (for our example), we find that the spatial dependence tends to decrease the pair-creation probability-which is qualitatively consistent with the behavior for the Sauter-Schwinger effect in an inhomogeneous electric field (see, e.g., [20]). Note that this nextto-leading order correction vanishes at maxima and minima of the pulse, where g 0 ðΩξÞ is zero.
We expect that other field configurations, where the dependence on one spacetime coordinate is weak, can be treated similarly (e.g., tunneling through a weakly timedependent barrier or a light-front field pulse depending on x þ plus a pulse only weakly dependent on x − ). In the presence of an electromagnetic field A μ , one formally obtains the same equations (41) for α and β, but the subsequent steps such as the transformation to new coordinates r and s are more involved. It is still possible to use S=m 0 and R=m 0 as coordinates, but they are not locally orthogonal anymore. Alternatively, one can obtain the coordinate s ¼ S=m 0 in a similar way as before and then construct another locally orthogonal coordinate, but the equations for the Bogoliubov coefficients α and β become more sophisticated nevertheless [53]. However, the main strategy should also be applicable in this case.
The curves are deflected when they reach the region of the pulse. This deflection is manifest in a change of a curve's slope dx=dt after passing the region of nonconstant mass. The above equations imply that the change of the slope with the parameter τ is d dτ To approximate the change in the slope, we use the initial form of the characteristic curves (B3) in (B4). We expect this to be a good approximation if m is only weakly space dependent; i.e., the timescale on which the value of the mass changes is much smaller than its length scale. For p ¼ 0 this approximation yields d dτ Thus, the slope after passing the region of nonconstant mass is approximately Hence, a characteristic starting at x ¼ x 0 will have the form after passing the pulse. The intersection of this characteristic and the one starting at x 0 þ δ is at For δ → 0 this gives the intersection of two neighboring curves Consequently, the focal point or onset of the caustic surface is where this is at its minimum with respect to x 0 . For a weakly space-dependent mass of the form (84), we get which immediately leads to the proportionality given in (55). For fðτÞ ¼ sechτ and gðχÞ ¼ sechχ we find the minimum to be