Discrete worldline instantons

The semiclassical approximation of the worldline path integral is a powerful tool to study nonperturbative electron-positron pair creation in spacetime-dependent background fields. Finding solutions of the classical equations of motion, i.e. worldline instantons, is possible analytically only in special cases, and a numerical treatment is nontrivial as well. We introduce a completely general numerical approach based on an approximate evaluation of the discretized path integral that easily and robustly gives the full semiclassical pair production rate in nontrivial multi-dimensional fields, and apply it to some example cases.

The semiclassical approximation of the worldline path integral is a powerful tool to study nonperturbative electron-positron pair creation in spacetime-dependent background fields. Finding solutions of the classical equations of motion, i.e. worldline instantons, is possible analytically only in special cases, and a numerical treatment is nontrivial as well. We introduce a completely general numerical approach based on an approximate evaluation of the discretized path integral that easily and robustly gives the full semiclassical pair production rate in nontrivial multi-dimensional fields, and apply it to some example cases.

I. INTRODUCTION
An as of today still experimentally unconfirmed prediction of quantum electrodynamics is that of nonperturbative electron-positron pair creation in the presence of a strong electric field [1-3]. Schwinger [4] gave the pair production rate per unit volume P e + e − (or more properly the rate of vacuum decay [5]) in a constant, homogeneous electric field E in 3+1 dimensions as ( = c = 1) where q is the elementary charge and m the mass of the electrons and positrons. The generalization to inhomogeneous and time dependent background fields is far from straightforward, since this is a nonperturbative effect (as is visible from the inverse dependence on q and E in the exponent of (1)). Apart from the fundamental interest in this effect as a prototypical example for a nonperturbative phenomenon in quantum field theory, a better understanding is also desirable in view of the various experimental initiatives aimed at reaching ultra-high field strengths [6].
It is in general difficult to obtain the pair production probability for multi-dimensional fields. While there has recently been some progress [7][8][9][10][11][12] in direct numerical computation of the exact probability for multidimensional fields, we will instead focus on an approach using the worldline path integral.
This formulation is an alternative to path integrals over fields to express amplitudes in quantum field theories. The first steps in this direction were pioneered by Fock, who expressed solutions of the Dirac equation via a four dimensional Schrödinger-type equation with space and time parameterized by an additional parameter [13]. After Nambu emphasized how beneficial this representation would be in the path integral approach [14], Feynman derived the Klein-Gordon propagator [15] and Dirac propagator [16] in this worldline formulation. In parallel, Schwinger's famous paper [4] used a similar representation.
It is possible to approximate this worldline path integral for inhomogeneous fields numerically using discretization and Monte Carlo methods [17][18][19][20]. Although our method is based on discretization as well, we use an instanton approach to compute the integrals instead of statistical sampling. Both Feynman and Schwinger mentioned the four dimensional particle's equations of motion in the classical limit, but the first explicit mention of an instanton approximation to the (Euclidean) worldline path integral was given by Affleck, Alvarez and Manton [21]. They derived the pair production rate for a constant homogeneous background field in a way that is very similar to the method used today. The approach was extended to inhomogeneous fields, including the sub-leading fluctuation prefactor [22,23]. An exact analytic treatment is possible in some simple cases [22][23][24] and analytic approximations allows us to study suitable limiting cases [25][26][27], but in general solutions of the instanton equations of motion have to be found numerically. This can be done using, e.g., shooting methods [28], but the highly nonlinear nature of the equations of motion makes this approach very unstable.
After briefly sketching the semiclassical approximation of the worldline path integral in section II, we introduce a different approach to numerically evaluate the path integral by discretization in sections III and IV, and a method to trace families of solutions over a range of field parameters in section V. Finally, we will apply the method to some example cases, both with results known analytically (to assess the accuracy of the approximation) and new examples to demonstrate the scope of the approach in section VI.

II. WORLDLINE INSTANTON METHOD
The central object of the method is the effective action Γ M , defined using the vacuum persistence amplitude e iΓM := 0 out |0 in . (2) We take the probability for pair creation to be the complement of the vacuum remaining stable, so the subscript M denoting the physical, Minkowskian quantity. We will work with the Euclidean effective action Γ, related to the Minkowski expression by Γ M = iΓ, so Γ M = Γ [23]. The Euclidean worldline path integral for spinor QED reads (see, e.g., [29,30]) where A µ is the Euclidean potential representing the external electromagnetic field F µν and x µ (τ ) are periodic worldlines in Euclidean space parametrized by the "proper time" τ withẋ µ = dx µ /dτ . There exist a couple of different representations of the spin factor, see [29,31]. We will use with P denoting path ordering, tr the trace over spinor indices and σ µν the commutator of the Dirac matrices For simple fields, the Euclidean potential A µ and field tensor F µν are purely imaginary, so iA µ and iF µν are real. We immediately introduce dimensionless quantities using some reference field strength E, which makes a numerical treatment possible and simplifies the following derivation, and also rescale the integration variablẽ We can now exchange the order of integration, so we can perform theT -integration using Laplace's method. Due to our rescaling, m 2 /qE is singled out as the large parameter of the expansion, while all other quantities are of order unity. We obtain the saddle point , so including the quadratic fluctuation around the saddle we arrive at the approximation with the non-local (due to a[x]) action and the spin factor Note that in (10) we symbolically restored the original parametrization Dx(τ ) in the path integral differential, this will be relevant for the discretization in the next section.
Applying Laplace's method to the path integral, we need to find a pathx µ (u) that satisfies the periodic boundary conditions and extremizes the exponent in (10), so a solution to the Euler-Lagrange equations (the Lorentz force equation in this case) Contracting (13) withẋ µ we see that (due to the antisymmetry ofF µν )ẋ 2 = const. = a 2 , simplifying the instanton equations of motion tö The prefactor of the Laplace approximation is given by the second variation of the action around the classical solution to (14), amounting to an operator determinant. The determinant has to be defined carefully, but we can completely sidestep this complication by instead performing Laplace's method after discretization, when we can calculate the fluctuation prefactor by standard methods of linear algebra.

III. DISCRETIZATION
We approximate (10) by discretizing the trajectories x µ (u) into N d-dimensional points (in general d = 3 + 1, but for simple field configurations it is possible to only consider d = 1 + 1 or d = 2 + 1 dimensions, so we will keep the dimensionality variable): The velocity is then approximated using (forward) finite differencesẋ with the identificationx N µ ≡x 0 µ . Discretizing the path integral requires a normalization factor for eachx i µ -integral. We could find these factors by performing the integration in the free case, however this is not necessary: In the derivation of (4) the path integral arises as an ordinary nonrelativistic transition amplitude, so we can use Feynman's normalization 1/A for each integral (see, e.g., [32]), with Using this normalization and replacing the N × d integrations by the dimensionless versions we arrive at the discretized worldline path integral As we have now expressed everything in terms of the dimensionless variables, we will from now on drop the tilde. We still need discretized expressions for a, A and Φ, the square brackets denoting dependence on all points, instead of a particular choice of indices. The form of discretization of the gauge term is not at all obvious, other choices like having just A ν (x k µ ) or A ν (x k+1 µ ) or evaluating the gauge field between points A ν ((x k µ + x k+1 µ )/2) would yield the same classical continuum limit. That does not mean however that the resulting propagator is the same, see [33][34][35][36]. The midpoint prescription in (20) arises when the path integral representation is derived from the vacuum persistence amplitude using the time slicing procedure, see e.g. [37].
Special care has to be taken to define a discretized expression for the spin factor that obeys path ordering. Instead of approximating the integral by summation and taking the exponential, we employ the product representation of the exponential function which is automatically path ordered (cf. [38]), The finite dimensional integral (18) can now be approximated using Laplace's method as well, by finding To ease notation, we will condense the proper time index l and the spacetime index µ into a single vector so a discrete instantonX has the property Equation (24) describes a system of N ×d nonlinear equations in N ×d unknowns, which can be solved numerically using the Newton-Raphson method or a similar root finding scheme.
In this discretized picture, the fluctuation prefactor is readily computed as well, via the determinant of the Hessian of A giving the full semiclassical approximation of the discretized worldline path integral with a cl . . = a[X]. If the function A[X] were entirely well-behaved we would be done now, we would just need to find solutions of (24) and plug them into (26). The Gaussian integration resulting in the determinant prefactor however is only defined for positive definite matrices in the exponent, which our Hessian H is not.

IV. REGULARIZATION OF THE PREFACTOR
We have two problems with the Hessian matrix of the action A. One is that of negative eigenvalues of H. The corresponding direction in the Gaussian integration diverges, and the integral has to be defined by analytic continuation. A single negative mode (which is present for a static electric field) thus turns the determinant negative, and the whole expression (26) imaginary. This could seem troubling at first, as the pair production is given by the real part of the Euclidean effective action. For a field not depending on time we expect a volume factor from the x 4 -integration though, which has to be purely imaginary for a real temporal volume factor V t = −iV x4 .
A more serious technical issue is that of zero modes. One or more zero eigenvalues of H immediately spoil our result, so they have to be removed from the integration in some way. One zero mode that is always present in the worldline path integral is the one corresponding to reparametrization. Due to the periodic boundary conditions we can move every point of the curve along the trajectory without a change in action. We would thus like to separate the integration in this direction (resulting in a "volume factor" of the periodicity, in our rescaled expression just unity) from the other integrations.
We will use the Faddeev-Popov method [39] to perform this separation. While it is commonly used to remove gauge-equivalent configurations from a gauge theory path integral, it can be applied to this simpler scenario as well. We insert a factor of unity into the path integral in terms of the identity where χ(λ) is some function chosen so that χ = 0 fixes the zero mode, λ parametrizes the symmetry and w is the number of times χ(λ) = 0 occurs over the integration interval [40]. The idea is now that the λ-integration can be performed due to the symmetry of the path integral, resulting in the desired volume factor and a Dirac delta that fixes the corresponding mode. This is especially elegant for a discrete numerical evaluation of the semiclassical approximation, as we can use an exponential representation of the delta function where the Gaussian integration over the zero mode produces a factor of √ ε canceling the prefactor, enabling us to simply set ε = 1. We insert the factor of m 2 /qE for convenience, so the action A in (18) just gets an additional term πχ 2 .
To fix the reparametrization mode, we take (cf. [40][41][42]) which is chosen so that at the saddle point. Due to the translation invariance we can set λ u = 0 in the integrand so the λ u integration is equal to one. This means we only need to add the (discretized version of) χ u (0) to the action as in (28), the second derivatives to H and a factor of m 2 /qE from (28) to the prefactor, and just proceed as if no zero mode were present. Other zero eigenvalues appear if the electric background field does not depend on all spacetime coordinates. They are of course easier to deal with, we could just omit the corresponding integrals and add a volume factorL µ (the tilde is to stress that this is in terms of the dimensionless coordinates) per invariant direction x µ . We can, however, treat these just as the reparametrization direction, which simplifies a numerical implementation that supports arbitrary fields. Choosing χ µ to be the average of x µ along the trajectory we obtain the volumẽ L µ , and again a factor of m 2 /qE.
To summarize, our final expression for the semiclassical approximation of the effective action is where the appropriate terms of χ and its derivatives have been added to A and H, N 0 is the number of invariant spacetime directions, and V N0 the corresponding volume factor (with units reinstated). Note that (31) unambiguously contains the full prefactor including spin effects for an arbitrary background field, without having to resort to limiting cases to determine any normalization constants. In addition, the reference field strength E enters only in the combination qE/m 2 in front of the action and in the prefactor, which has two advantages. First, having found an instantonX, we can evaluate (31) for arbitrary values of qE/m 2 without any additional computational effort. Secondly, the accuracy of the discretization does not depend on the field strength, so there are no numerical instabilities for small E. Figure 1 shows how the discretization error scales with the number of points N for a constant, homogeneous electric field. For scalar QED (that is, without the spin factor Φ) the error in the prefactor decreases as N −1 as expected for a first order discretization procedure. As the first variation of the action vanishes for an instanton, the error of the exponent even decreases as N −2 . For spinor QED, on the other hand, the error in the prefactor decreases as N −2 as well. The reason for this is not obvious, as the only difference is an additional, seemingly independent multiplicative spin factor.

V. NUMERICAL CONTINUATION
For most fields we are interested in, there is one (or multiple) parameters that we would like to vary, for example the timescale of a pulsed field or the inhomogeneity of a spatially varying field configuration. Let us denote such a parameter γ. In general we are interested in the full family of instantonsX(γ). Methods to numerically map such a solution space are known as numerical continuation algorithms [43,44].
If we know an instanton for a particular value γ i of the parameter (e.g. the limit ω → 0 for a timedependent pulse), we can use it as the starting point for the numerical solution of (24) for a parameter value γ i+1 = γ i + ∆γ, which is the method used in [45]. If we choose a sufficiently small ∆γ, we can expect the root finding procedure to quickly converge. This process is called natural parameter continuation, because we vary a physical parameter of the problem at hand, instead of introducing an artificial variable to blend between an easy and our actual problem (e.g. solving 0 = G(X, γ) . . = γF (X) + (1 − γ)F 0 (X)).
Natural parameter continuation works well if the so-lutionsX(γ) depend on the parameter in a smooth and uniform manner. If, however, the dependence on γ varies strongly, it is difficult to choose appropriate step lengths ∆γ. For some spatially inhomogeneous fields the instantons even grow infinitely large in some limit γ → γ crit , so we need to take ever smaller steps to reach this value. We could, in principle, adaptively adjust the step length when the root finding for the next parameter value converges poorly, but there is an easier method of choosing the increment ∆γ: Natural parameter continuation can be viewed as a predictor-corrector scheme, with the "zeroth-order" predictor step of just taking the last solution as the starting point for the next parameter, and performing the numerical root finding as a corrector step. We can find a better prediction by taking the γ derivative of (24), yielding the Davidenko differential equation [46]: and thus, provided that H is invertible (which it is by our regularization scheme), We can now use (33) in two ways: first, having found an instantonX i for a parameter value γ i , it tells us in which way the instanton for a slightly different value of γ differs from the current one, so we can use it as a much improved predictor in our predictor-corrector scheme, i.e. X i+1 ≈X i + ∆γ dX/dγ. In fact, we could directly integrate (33) to obtain all solutions. Unfortunately, evaluating the Hessian is costly and we can afford a much larger step size by performing the corrector steps. As a compromise it is possible to perform multiple steps according to (33) before starting the root finding routine. Furthermore, we can use the derivative to scale the step ∆γ by instead specifying a maximum (or mean) difference between the points ofX i and the proposed guess forX i+1 , or even a fixed arclength ∆s of the solution curve in R N ×d+1 , ∆s = (∆γ dX/dγ) 2 + (∆γ) 2 ⇔ ∆γ = ∆s (dX/dγ) 2 + 1 .
A situation may be conceivable where it is not possible to parametrize the solution set asX(γ) at all, because such a function would not be single-valued or have infinite slope somewhere. In this case, we can parametrize both the solution and the parameter γ by a new parameter whereH is now an (N d+1)×(N d) matrix, so (35) has to be augmented by an additional condition. This is chosen to be a constraint on the orientation and the "velocity" of the flow 1 = dȲ /du , soȲ (u) is parametrized by arclength, hence the name pseudo-arclength continuation (pseudo because this is only approximately true, as we are taking discrete steps). As long as γ is a suitable parameter, this is equivalent to (34), which is what we will be using in the following.

VI. APPLICATIONS
Let us now apply the method outlined above to some background fields. The strategy in all cases is to start  Imaginary part of the Minkowskian effective action for the spacetime bump profile E = E cosh −2 (ωt) cosh −2 (kz)ez with k = 3ω and E = 0.033m 2 /q. Here (and in the following cases) there are no analytical results to compare with, so we just add the connecting dashed lines as a guide to the eye.
with a limit that is reasonably close to a static, homogeneous field and perform pseudo-arclength continuation to map the solution space for a chosen parameter range. In all figures depicting worldline instantons we color the homogeneous limit (i.e. a circular instanton) purple, and all further instantons proportional to the change in action (blue for a decrease, red for an increase, so blue means more, red less pair production). In all figures that show the full effective action we choose E = 0.033m 2 /q for the reference field strength. This is simply the value we already used in earlier works, and it does not influence the quality of the discretization in any way. We also use N = 500 points in the discretization, which yields good accuracy while it still takes less than thirty seconds to obtain the family of instantons in the cases below, with the exception of the e-dipole pulse.

A. Temporal Sauter pulse
First let us consider simple, one-dimensional inhomogeneities where we can compare to analytic results. As an example, we choose the Euclidean four-potential iA 3 = tan(γ ω x 4 )/γ ω describing the (physical) field E = E cosh −2 (ωt)e z with the Keldysh parameter γ ω = mω/qE [47]. Since the field does not depend on any spatial coordinates, we have N 0 = 3 translational zero modes that need to be held fixed.
The analytical worldline instanton result for this field is [23] Γ Sauter ω In Figure 2 the first plot shows a family of instantons in the range of 0 < γ ω < 3.5, and the top panel in Figure 3 compares the numerical result (31) in this parameter range to the analytical expression (36), showing near-perfect agreement.

B. Spatial Sauter pulse
We can also consider the spatially inhomogeneous profile iA 4 = tanh(γ k x 3 )/γ k describing the (physical) field E = E cosh −2 (kz)e z with (the spatial analog of) the Keldysh parameter γ k = mk/qE. The analytical result is related to (36) where the instanton is now confined in x 3 -direction and we obtain a "temporal volume factor" V t instead. The worldline instantons in this field for the range 0 < γ k < 1 are depicted in the middle of Figure 2, and the comparison of the numeric result and the analytic expression (37) in the bottom panel of Figure 3.

C. Space-time Sauter pulse
As a simple example of a both space-and timedependent background we choose the product of the preceding profiles with γ := γ ω = γ k /3, i.e. iA 3 = cosh −2 (3γx 3 ) tan(γx 4 )/γ. The resulting worldline instantons in the range 0 < γ < 2.5 are shown in the bottom plot of Figure 2 and the resulting pair production rate in Figure 5. With the chosen ratio γ k /γ ω = 3 the spatial inhomogeneity dominates for small γ, giving larger instantons, increased action and lower pair production, while above γ ≈ 1 the time dependence takes over and produces smaller instantons, reduced action and increased pair production.

D. Multidimensional instantons
In [28] multidimensional instantons were found for background fields that depend on multiple spatial coordinates using the shooting method. We can obtain instantons for these fields using discretization as well. Consider the potential tanh(kx 1 + kx 2 ) 1 + (kx 1 ) 2 + 10(kx 2 ) 2 from Figure 1 in [28] (with the factor of 1/ √ 2 added so the peak strength is 1). This yields three dimensional FIG. 6. Worldline instantons for the four-potential (38) from [28]. As before, stronger inhomogeneity stretches the instantons (in a more complicated way than for the one dimensional fields) and increases the action. instantons (in x 1 -, x 2 -and x 4 -direction). Figure 6 depicts a family of instantons in a three dimensional plot, while Figure 8 shows all two dimensional projections of the same trajectories. The resulting pair production rate is given in Figure 7.

E. Plane wave plus electric field
In [48] we already applied the discrete worldline instanton method to calculate the pair creation rate for the superposition of a weak propagating plane wave and a constant field, a variant of dynamically assisted pair production [25]. Different pulse shapes have been considered for the weak field before [49], however a plane wave is special in that it cannot produce pairs on its own, so the process is fully nonperturbative for all frequencies.
In the case of parallel polarization (the plane wave and the constant field point in the same direction, but perpendicular to the propagation direction) this combination can be represented by the four-potential The method can handle the perpendicularly polarized case just as well, however that leads to four-dimensional Pair production rate for a constant field (EStrong = 0.033m 2 /q) with superimposed plane wave (E Weak = 10 −2 EStrong). The temporal volume factor Vt arises from the number of instantons, one per oscillation of the wave at a fixed spatial point. The dashed line is added as a guide to the eye.
instantons that are cumbersome to visualize.
In contrast to the examples considered before, the field (39) leads to complex instantons, in particular purely real x 3 (u), x 4 (u) and purely imaginary x 1 (u). A family of instantons is shown in Figures 9 and 11, while the full pair production rate is given in Figure 10.  Figure 9, projected onto the coordinate planes.

F. E-dipole pulse
An especially interesting, highly non-trivial example is that of an e-dipole pulse. It is a solution to Maxwell's equations in vacuum that represents a localized pulse of finite energy [50]. It saturates the theoretical upper bound of peak field strength for given laser power [51] and is thus in a sense the optimal (and at the same time physically viable) configuration to study pair creation [52]. Its name stems from the structural similarity to dipole radiation, however it does not suffer from the strong singularities at the origin for a simple non-stationary dipole.
The electromagnetic field of the e-dipole pulse can be given in terms of a driving function g using the vector Z [52] Z = e z d |r| g(t + |r|) − g(t − |r|) , We choose the function g(t) = t 4ω 2 e −ω 2 t 2 + √ π 8ω 3 (1 + 2ω 2 t 2 ) erf(ωt) (41) and the virtual dipole moment d = 3E/4, so that at the origin We cannot immediately apply the instanton approach to this field since it is not given in terms of a fourpotential. It is however possible to obtain an expression for the potential in coordinate gauge A(x) · x = 0 from the field tensor [53], For the field (40) this gives a lengthy expression, which can now be used to obtain worldline instantons.  Figure 12 shows the result. The top plot compares the instanton action for the e-dipole pulse to the action for a field with Gaussian time depencence only. Due to the additional spatial inhomogeneity in the e-dipole field the action is slightly larger (and thus pair production slightly lower) than for the purely time dependent pulse. We can also compare the full imaginary part of the effective action with the locally constant field approximation (in the bottom plot of Figure 12), which can be calculated using the saddle point method for E below the critical field strength, giving As expected, the worldline instanton result tends to the locally constant field approximation for small values of γ, while it is exponentially larger for higher γ. For the parameters considered in [52] the adiabaticity is very small, with γ < 10 −3 , so the locally constant field approximation is accurate. For high frequency pulses however, the pair production rate is higher than the constant field estimate.

G. Transversal standing wave
Let us now briefly consider a purely transversal inhomogeneity. Two counterpropagating laser beams create a standing wave pattern, i.e. E = cos(ωt) cos(kx)e z with k = ω. In [10] the authors find that omitting the spatial inhomogeneity leads to qualitatively incorrect results in the high frequency regime. In [8] strong deviations in the momentum spectrum have been found in the homogeneous approximation as well.
In the semiclassical regime and for the total integrated rate however we can now check that approximating the standing wave by an oscillating homogeneous field works well. It is easy to see that the transversal inhomogeneity does not change the instanton and thus the action [49], but the effect on the prefactor is not as obvious. Calculating the full effective action using the discrete instantons shows that while the prefactor does change, the difference from the homogeneous result is small and barely visible, see Figure 13. Note, however, that the momentum spectrum could still display noticeable differences between the standing wave and the purely time dependent field.

H. Constant electric and magnetic fields
In all examples up to now, the spin factor had only a small impact, apart from the trivial factor of two in the pair production probability. Let us finally consider a simple example where there is a large, qualitative difference between scalar and spinor QED, a parallel superposition of constant electric and magnetic fields of strength E and B respectively.
The (first term of) the effective action for this combination is given by (see e.g. [54] and references therein) Figure 14 depicts the prefactors of these expressions, so the B/E dependence, together with the discrete instanton result, showing perfect agreement.

VII. SUMMARY AND CONCLUSION
We have introduced a new approach to numerically implement the worldline instanton method for electronpositron pair creation. We use a discretization scheme that turns the infinite-dimensional path integral into a finite dimensional integration that we can then perform using Laplace's method. Crucially, this also means that the fluctuation prefactor is simply given by a finite dimensional determinant that can be computed without the great care that is needed for a properly normalized treatment of the functional determinant.
After having implemented the necessary root finding and continuation steps outlined in sections III, IV and V, full pair production results for arbitrary background fields can be obtained in minutes. Section VI gives a (by no means exhaustive) sample of such applications.
Although we used a frequency or inhomogeneity scale as the continuation parameter in all examples, we could have also chosen a different field parameter like the polarization direction or the ellipticity of the field, or even an entirely synthetic parameter to slowly transition to an especially complicated field configuration.
In this paper we have only considered cases for which there is one dominant instanton, which is continuously connected to a circular one in the constant field limit. It would be interesting for future studies to consider cases where there are more than one instanton, and where some of them might have a nontrivial topology.