Pair production in temporally and spatially oscillating fields

Electron-positron pair production for inhomogeneous electric and magnetic fields oscillating in space and time is investigated. By employing accurate numerical methods (Furry-picture quantization and quantum kinetic theory), final particle momentum spectra are calculated and analyzed in terms of effective models. Furthermore, criteria for the applicability of approximate methods are derived and discussed. In this context, special focus is placed on the local density approximation, where fields are assumed to be locally homogeneous in space. Eventually, we apply our findings to the multiphoton regime. Special emphasis is on the importance of linear momentum conservation and the effect of its absence in momentum spectra within approximations based on local homogeneity of the fields.


I. INTRODUCTION
Quantum field theory has taught us to see the quantum vacuum as a key factor towards understanding fundamental physical processes. The reason is that the quantum vacuum is far from being empty but can rather be described as a nonlinear medium with fluctuating virtual particles acting as mediators of interactions. In quantum electrodynamics (QED), for example, it is perpetual creation and annihilation of virtual electrons, positrons, and photons that gives rise to new, staggering phenomena such as Sauter-Schwinger effect [1][2][3], vacuum birefringence [4][5][6], and light-by-light scattering [7][8][9][10].
Especially in the context of strong-field QED, understanding these effective interactions and finding ways to probe vacuum nonlinearities have driven the evolution of the research field ever since; for reviews see Refs. [11][12][13][14][15][16][17][18]. However, testing these conjectures experimentally has been proven challenging. Due to the small cross sections of direct light-by-light scattering, it takes extremely high field strengths in order to obtain any signal of underlying nonlinearities. As a consequence, experiments have been limited to atomic fields and highly charged ions (see Refs. [19][20][21][22], where high-energy processes are discussed). New theoretical developments concerning spontaneous pair production in low-energy ion collisions have been recently reported in Ref. [23] (see also references therein).
The famous SLAC experiment on multiphoton pair production marked a turning point in evaluating the state of the research field by demonstrating that state-of-the-art, high-intensity laser systems could perfectly probe the strong-field regime [24,25]. In the years following this seminal observation newly built laser facilities have been recognized as an opportunity to advance the field even further, thus a variety of new theoretical predictions regarding the feasibility of observing similarly astonishing phenomena have been published (see, e.g., Refs. [26][27][28]). Especially creation of matter through high-intensity electric fields, the famous Schwinger effect, and the prospect of finding direct proof for quantum vacuum nonlinearities has gained renewed interest as current laser technology is on the brink of investigating these effects in laboratory.
This development has also huge impact on the theoretical aspect of the research field as experimental conditions are often far away from perfect theoretical settings, cf. Refs. [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48] considering various temporal and spatial pulse shapes. The possibility of working with an imperfect vacuum as well as having imprecise laser and detector equipment plays a huge role in evaluating the chances for detecting, e.g., signal photons [49,50]. Hence, this article is mainly devoted to analyzing currently used numerical methods and approximations in studies on electron-positron pair production. To be more specific, we evaluate the Furry-picture quantization formalism and quantum kinetic approaches in order to determine their strengths and weaknesses. Moreover, we test the applicability of the local dipole approximation, where spatial dependencies are treated as locally homogeneous. This is an important step forward in finding a suitable method to actually perform calculations with respect to pair production for given laboratory conditions.
For this work, we assume to have two counterpropagating laser beams which are focused at the same point creating a highly intense standing wave pattern. Such a setup provides a perfect environment for testing various numerical methods, because the temporal structure allows one to study Schwinger as well as multiphoton effects and the inhomogeneous spatial structure is ideal to see the implications of magnetic fields. Furthermore, the computational techniques applied offer a lot of flexibility allowing us to further investigate the fundamental differences in purely time-dependent electric fields and timedependent, spatially inhomogeneous backgrounds.
The article is structured as follows: In Sec. II we introduce the external laser field and discuss its characteristic traits. Secs. III, IV, and V are the basis of this paper as we individually introduce the different methods step-by-step. To be more specific, Sec. III is about the Furry-picture quantization, in Sec. IV a phase-space formalism is introduced, and in Sec. V we discuss the local dipole approximation (LDA). In the largest section, Sec. VI, we present and discuss our results on the basis of a comparison between the methods for short-pulsed fields, Sec. VI A, as well as for many-cycle pulses, Sec. VI B.
The last section, Sec. VII, contains the conclusions of our study.
Throughout the article, we use = c = 1 and display observables in terms of the electron mass m. The electron charge is e < 0.

II. EXTERNAL FIELD CONFIGURATION
One of our goals is to study the performance of different numerical methods regarding calculating the particle yield for spatially inhomogeneous fields within various setups. Moreover, we are determined to inspect various multiphoton signatures in standing wave patterns which approximate a scenario involving two counterpropagating laser pulses. Hence, we introduce a vector potential of the form This expression represents a standing wave being linearly polarized along the x direction, with a cosine profile along the z axis and a super-Gaussian temporal envelope. Four individual field parameters allow for great flexibility; the peak field strength ε, the pulse length τ and the field frequencies ω and k z . The critical field strength E cr = m 2 /|e| has been factored out for the sake of convenience. The electric and magnetic fields are derived from Eq. (1), In this way a single vector potential is sufficient in order to interpolate from the Schwinger to the multiphoton regime and from few-cycle to many-cycle pulses. In addition, all fields vanish at asymptotic times as A(t → ±∞) → 0. As we will discuss in Sec. VI, neither does ω = k z hold in general nor do the parameters ω and k z always represent physical quantities. Although Maxwell's equations lead to |k z |/ω ! = 1, our toy model allows for a gradual change of this ratio in order to switch on/off the magnetic field component and make the spatial variations more/less pronounced. In what follows, it will provide additional insights into the particle dynamics and its relation to the features of the pair production process. Furthermore, for few-cycle pulses ωτ ≈ O(1), a Fourier analysis of the pulse profiles does not yield one dominant field frequency. In many-cycle pulses ωτ 1, the situation changes drastically in particular close to the critical frequency ω ∼ m. Then, it can be shown that ω indeed corresponds to the photon energy and k z represents the photon momentum.

III. FURRY-PICTURE QUANTIZATION
The general formalism of quantization within the Furry-picture is described in detail in Ref. [51].
We will first briefly discuss how the rigorous QED approach allows one to extract the necessary pairproduction probabilities from the one-particle solutions of the Dirac equation including the interaction with a classical external field A µ . After that we will show how this approach can be implemented in the case of the standing-wave background (1).
The external field is assumed to vanish outside the time interval t in < t < t out , where in the case of the field configuration (1), t in/out → ∓∞. We introduce two complete sets of the one-particle Hamiltonian eigenfunctions at t = t in and t = t out , respectively: where H e (t) = α[−i∇ − eA(x)] + eA 0 (x) + βm, and the eigenvalues denoted by plus (minus) are positive (negative). These sets are orthonormal and complete with respect to the usual inner product.
In terms of these functions, the field operator can be decomposed as where we have introduced the electron (positron) creation and annihilation operators a † n (b † n ) and a n (b n ), respectively. These operators satisfy the usual anticommutation relations. The Dirac-field Hamiltonian H e (t) = ψ † (x)H e (t)ψ(x)dx is then diagonalized at time instants t in and t out .
We will turn to the Heisenberg representation which is achieved by the unitary evolution operator We perform transformation by means of the operator U e (0, t), so the field operator gains now a temporal dependence: The creation and annihilation operators are transformed according to a n (in) = U e (0, t in )a n (t in )U † e (0, t in ), a n (out) = U e (0, t out )a n (t out )U † e (0, t out ).
The other creation/annihilation operators with indexes (in) and (out) are defined similarly. The anticommutation relations match those taking place in the Schrödinger picture. The in (out) vacuum in the Schrödinger and Heisenberg representation is denoted by |0, t in (|0, t out ) and |0, in (|0, out ), respectively, i.e., |0, in = U e (0, t in )|0, t in and |0, out = U e (0, t out )|0, t out .
One can demonstrate that the evolution of the field operator ψ(x) is governed by the equation it is a solution to the similar Dirac equation as that regarding the timedependent one-particle solutions. From this, it follows that the field operator can be represented as where + ϕ n (x) and − ϕ n (x) are the so-called in solutions of the Dirac equation evolved from the corresponding functions (4). One can also show that where + ϕ n (x) and − ϕ n (x) are the out solutions, which coincide at t = t out with the eigenfunctions (5).
Suppose that the Schrödinger field operator ψ(x) is known, the Heisenberg operator ψ(t out , x) can be constructed in two different ways. First, one can perform the transformation (9) at t = t out .
Second, one can evolve the Heisenberg operator ψ(t in , x) by means of the one-particle propagator. This allows one to express the out operators in terms of the in operators according to where Note that these inner products do not depend on time.
We can now evaluate the number density of electrons produced. The starting point is the following which is obvious in the Schrödinger representation. Using the relation U e (t out , t in ) = U e (t out , 0)U e (0, t in ) and inserting the identity operator 1 = U e (t out , 0)U † e (t out , 0) between a † m (t out ) and a m (t out ), we find With the aid of Eqs. (14) and (15), one obtains The analogous expressions for the number density of positrons have the form It turned out that the full information about the number of particles and their spectrum can be extracted from the (in and out) one-particle solutions of the Dirac equation incorporating the interaction with the external background.
We will now discuss how one can employ this approach in the case of a standing wave (1), which can be represented as A x (t, z) = Q(t) cos k z z. Since the external field depends only on t and z, the in and out solutions of the Dirac equation have the form The time-dependent functions χ satisfy the following conditions: where p 0 = m 2 + p 2 , and the bispinors u p,s and v p,s form an orthonormal and complete set (for given p). Since the external field is monochromatic as a function of z, one can represent the functions ζ χ p,s (t, z) as where ζ = ±. Similar series are introduced for ζ χ p,s (t, z). In terms of the time-dependent Fourier coefficients ζ w j p,s (t), the Dirac equation takes the following form where K = k z e z and j ∈ Z, so one has to solve an infinite system of ordinary differential equations (in practical calculations, it can, of course, be truncated at sufficiently large |j|). The same holds true for the in functions ζ w j p,s (t).
The G matrices can be evaluated according to where The number density of electrons produced reads Propagating the Fourier coefficients + w j p,s (t) backwards in time according to the system (27), one can obtain the electron number density for given quantum numbers p and s. It turns out that in the case of a linearly polarized standing wave (1), the number density does not depend on s, n p,s = n p (see, e.g., Ref. [52]). The individual contributions for given l in Eq. (36) can be interpreted as separate channels of pair production by absorbing external-field quanta of total momentum lK. This aspect will be addressed in Sec. VI B.

IV. DHW FORMALISM
In heavy contrast to the Furry-picture formalism, Sec. III, where one solves for the wavefunctions of one-particle solutions to the Dirac equation, kinetic theories are formulated on the basis of transport properties of a quantum plasma. The advantage of such a description is mainly given by the fact that the governing set of equations of motions have only to be calculated once in order to determine the full particle spectrum. The downside of this approach is the presence of non-local operators, required to properly account for dynamical pair production. Furthermore, the DHW formalism accesses the full particle phase-space. A speciality of kinetic approaches, that can become numerically challenging.
Similarly to Sec. III for the Furry-picture formalism, we will only state the key points in the derivation of the DHW (Dirac-Heisenberg-Wigner) approach, cf. Ref. [56] for a complete derivation of the governing equations of motion or Refs. [57,58] for a more in-depth look at the features of the formalism.
We begin the derivation of the governing quantum kinetic equations of motion stating the QED and, consequently, the Dirac equation where D µ = (∂ µ + ieA µ ) and D † µ = ∂ µ − ieA µ are the covariant derivatives with a vector potential A µ that vanishes at asymptotic times, and γ µ are the gamma matrices. The fundamental quantity within the DHW approach is given by the gauge-invariant density operator with the center of mass coordinate r, the relative coordinate s and the Wilson line factor U (A, r, s) = exp ie In order to obtain a proper phase-space formalism, we perform a Fourier transform in s leading to the covariant Wigner operatorŴ which is properly defined in terms of four-position r and four-momentum p coordinates. Combining Dirac and adjoint Dirac equation (38)- (39) with derivatives of this operator yields two coupled operator where the derivatives ∂ µ are replaced by nonlocal, pseudo-differential operatorŝ In a crucial step towards obtaining a numerically feasible formalism based on distribution functions, we then take the vacuum expectation value of Eqs. (43)- (44). At this point we introduce a Hartree-type approximation of the form This basically transforms the field strength tensorF µν to a fixed c-number valued function F µν .
The key implication of this approximation becomes apparent when considering terms of the form of F µν (r)Ĉ (r, s) , which simplify tremendously The significance of Eq. (48) is given by the fact that we end up with a closed expression of only oneparticle correlations instead of an infinite BBGKY hierarchy of n-particle correlation functions. This, in turn, allows us to obtain the equations of motion for the covariant Wigner function A decomposition into Dirac bilinears turns the occurring matrix-valued equations into a set of equations for the C-number valued Wigner coefficients S, P, V µ , A µ , and T µν . Up to this point everything is formulated covariantly, which, however, poses the problem of having to know the solution to the governing equations of motion for all times. To avoid this, we project on equal times dp 0 2π thus essentially reformulating everything in terms of an initial-value problem. As all covariant Wigner coefficients are transformed analogously w (t, x, p) = dp 0 2π W (r, p) we eventually end up with a set of coupled integro-differential equations for the equal-time Wigner coefficients The vectors t 1 and t 2 are defined as t 1 = 2t i0 e i and t 2 = ijk t jk e i . The pseudo-differential operators are given by To mimic the vacuum-to-matter formalism introduced in Sec. III, we further employ vacuum initial At this point we can already interpret the functions of the different Wigner coefficients [58]. Most notably, s gives the mass density, v 0 yields the charge density and v detail the current density. Moreover, we can construct observables, valid when evaluated at asymptotic times, e.g., the particle distribution function as well as the particle momentum spectrum As we are employing periodic boundary conditions in x and are mainly interested in particle densities, the spatial integral in (64) requires proper normalization, As exact computations taking into account spatiotemporal inhomogeneities of complex external backgrounds may be extremely time-consuming, it is important to employ approximate techniques or design effective models to advance computations. If, for example, the external field only slowly varies in space and time, one can assume it to be locally constant, which provides a huge advantage since the Heisenberg-Euler effective Lagrangian [2] has a well-known nonperturbative closed-form expression in the case of a constant background. A general idea here is to calculate the necessary physical quantities in the presence of a constant external field and then perform integration over space and time by summing the individual (local) contributions. In the context of particle production, however, applicability of the local constant field approximation is still not perfectly well understood especially when evaluating the momentum spectra of particles produced in the presence of general space-time-dependent fields.
The main obstacle here is the fact that the particle created is being accelerated in the field, so that its dynamics can become overwhelmingly complicated in a general background. Furthermore, the local constant field approximation disregards photon absorption effects entirely thus restricting its validity severely. To overcome these issues, we abandon the strict concept of locally homogeneous fields and discuss a local density approximation instead. Our suggestion is to treat the external field as locally constant in space incorporating the temporal dependence exactly. Since approximating the field by a spatially uniform time-dependent background is often referred to as the dipole approximation [59], we will refer to this technique as the local dipole approximation (LDA). In what follows we will discuss how the LDA is implemented within the methods described in Secs. III and IV.
A. LDA in the Furry-picture approach The external field is approximated by a spatially homogeneous background whose temporal dependence coincides with that of the original field considered at a given position x. In this case, we perform exact computations of the momentum spectra by solving the Dirac equation and employing the formalism described in Sec. III. Varying then x and summing the results over the spatial coordinates, one obtains approximate momentum distributions of particles created. This approach is discussed in detail in Ref. [54]. Although it was already benchmarked against the exact predictions considering several different space-time-dependent field configurations [54], it was not assessed in the case of spatially oscillating fields. In the present study, we carry out the analysis of the LDA applying it to a standingwave background. In this case, instead of integrating over the whole z axis, one should average the results over the spatial period, i.e., z is being varied within the region [0, 2π/k z ].
Within the LDA computations, the external field does not depend on the spatial coordinates, which means that the coordinate part of the one-particle solutions of the Dirac equation is always given by exp(±ipx), so one has to construct only the time-dependent bispinor part being a solution to Eq. (27) for j = 0 in the absence of the coupling terms j ± 1, which do not appear within the LDA as the momentum transfer does not take place here. Instead of an infinite ODE system, one deals now with only one 4-component equation.
Let us point out several important features concerning the validity of the LDA. First, we note that such approximate calculations completely neglect the magnetic field component, so one may expect this approach to perform well only for sufficiently small k z and large ω since the characteristic ratio between the electric field strength (2) and the magnetic field strength (3) is k z /ω. Moreover, the LDA approximates the electric field component as well, so one can expect that averaging the results over z is reasonable only when the external field slowly varies in space, i.e., k z is again sufficiently small. Note that the LDA predictions, i.e., the average pair-production probabilities, do not depend on k z at all.
Thus, there should be a characteristic value k (LDA) z allowing one to employ the LDA for k z k In Sec. VI, we will investigate how k (LDA) z depends on the external field parameters and reveal several important features of this dependence.

B. LDA in the DHW formalism
Due to the fact that evaluating the pseudo-differential operators in Eqs. (59)-(61) is numerically challenging, the LDA provides a good opportunity to significantly reduce numerical costs. The phasespace formalism is inherently formulated as an initial value problem for transport properties, e.g., mass or charge density. Hence, it is impossible to perform independent calculations for separate mode functions.
Within the LDA, however, the governing equations of motion simplify drastically, as a Taylor expansion in x of the pseudo-differential operators terminates at first order and, as the spatial coordinates x i are fixed, the spatial derivative operators ∇ x , and with it the magnetic fields, vanish entirely leading to a much simpler system of equations, cf. Ref. [60]. To be more specific, within the LDA the set of partial differential equations decouples into a set of ordinary differential equations. The particle momentum spectrum is then obtained by simply summing the individual local density functions n (x i , p) over all instances of x i . See Refs. [61][62][63][64] for more details on a phase-space formalism for purely time-dependent electric fields.
In all our simulations, the LDA results obtained within the Furry-picture approach and DHW formalism were found to be identical (see Appendix A for more details).

VI. NUMERICAL RESULTS
We vary the external field parameters within the following domain: In what follows, we will separately consider the cases of short and long laser pulses, i.e., we will employ small and large values of the pulse duration τ .

A. Short pulses
In this subsection we choose τ = 5m −1 and vary ε, k z , and ω. For such a small value of the pulse duration, the number of carrier cycles N ∼ ωτ /(2π) is always less than 1 in our computations, so the carrier frequency is not well-localized in the Fourier spectrum of the pulse.
Validity of the LDA.-In order to determine k (LDA) z , we compared the exact momentum distributions obtained by means of the techniques discussed in Secs. III and IV to those predicted by the LDA varying first k z for given ε and ω and turning then to other values of these two parameters. In Fig. 1 we present the momentum spectra along the magnetic field direction y which were obtained exactly for various k z and with the aid of the LDA. We observe that the LDA provides indeed very accurate predictions in the case of small k z which become worse with increasing k z . Note that the discrepancy between the exact and approximate results appears to be independent of ω, i.e., it is fully determined by k z . This fact is illustrated in Fig. 1 and was confirmed in our computations employing other different field parameters. It means that k  results and indicate that the ratio k z /ω does not play any important role for given k z . The spectra computed within our study proved this ratio to be an irrelevant parameter in the case of such short laser pulses.
The fact that the validity of the LDA strongly depends on k z and does not depend on ω can also be accounted for by the requirement that "formation length" = 2m/|eεE cr | must be much smaller than the spatial wavelength λ = 2π/k z . In other words, the e + e − pair must be formed within a sufficiently narrow spatial region. This condition is equivalent to k z /(πmε) 1, which is indeed independent of ω and also explains why smaller values of k z are preferable. However, as will be seen below, this requirement is not the only criterion regarding the LDA justification.
Next we will examine how k (LDA) z depends on ε. In Fig. 3 we present the number density of electrons produced with zero momentum as a function of ε for various values of k z (ω = 0.2m). We observe that the discrepancy between the exact curves and the LDA predictions is not reduced for larger values of the field amplitude, so the results suggest that the concept of the formation length does not always provide a relevant parameter for justifying such kind of an approximate technique.
As was demonstrated in Ref. [54], in order to properly justify the LDA, one should also make sure that the classical particle trajectories are well localized within the spatial region where the external field does not change much. This condition arises due to the fact that the LDA treats different positions in space independently, and thus does not appropriately take into account the particle post-creation dynamics leading to inaccurate momentum spectra. This drawback becomes crucial when the particle trajectory stretches over a large part of the spatial period of the external standing-wave background.
To carefully examine this issue, we solve the equations of motion in the case of a relativistic classical particle interacting with the standing electromagnetic wave of the form of Eq. (1). We assume that the particle is initially at rest [p(t in ) = 0] at the position z(t in ) = z 0 . Since the external field does not depend on x and y, two momentum components have obvious temporal dependences: p y (t) = 0.
In order to evaluate the function z(t), we solve the following ODE system: As the particle is likely to be produced in the vicinity of the electric field maxima z = πn/k z (n ∈ Z), we choose z 0 to be close to the origin z = 0. To take the magnetic field component into account properly, we use a small nonzero value z 0 = λ/20 = π/(10k z ). We then propagate the solution z(t) from t = t in to t = 1.5τ , where the external field is about to vanish, and evaluate the ratio In this context, η is a measure for the particles' maximal displacement within the strong-field region of the background field. The criterion for applicability of the LDA is given by Trajectory analysis.-For the field configuration under investigation in Fig. 4 we find that only when k z is sufficiently small, particle positions only fluctuate locally., i.e., η 1. Starting with k z ≈ 0.4m, the value of η exceeds 0.1, so we expect a noticeable discrepancy between the LDA and exact simulations. A result nicely illustrated in Fig. 1, where calculations for k z = 0.8m deviate considerably from the curves obtained through the LDA. More important, calculating η as a function of ε, the discrepancies displayed in Fig. 3 can be well understood considering that for strong laser pulses (ε ≥ 0.4) particle trajectories become poorly localized. On the other hand, for weak fields (ε < 0.4) the first criterion based on the particle formation length is not fulfilled as l/λ > 0.24. Accordingly, the LDA is inaccurate for fields where k z ≥ 0.3m independent of ε. For slowly varying fields with k z ∼ 0.1m, on the other hand, η remains small for all ε ≤ 1. Hence, applying the LDA and averaging over z is completely justified according to this criterion.  As a similar argument can actually be made when η is considered on its own, it shows that it is the combination of both criteria, /λ 1 and η 1, that determines the validity of the LDA.
Magnetic fields.-Finally, let us note that the LDA leads to identical spectra in the y and z directions, whereas the exact results are anisotropic due to the presence of the magnetic field. In Fig. 5 we display this anisotropy for two different values of k z . We observe again that the accuracy of the LDA results strongly depends on k z .

B. Long-pulsed, high-frequency fields
External fields that exhibit a high frequency-duration product, ωτ 1, open up an opportunity for multiphoton pair production. Especially for rapidly oscillating fields, photon absorption easily becomes the main source for particle creation. If we further assume that high-intensity light sources in the hard x-ray regime operate at a specific carrier frequency, setups involving collisions of multiple pulses give rise to distinctive and well pronounced multiphoton signatures in the momentum spectrum.
As the focus in this section of the manuscript is primarily on evaluating the importance of the incident photons' momenta, we rely again on a comparison between full-scale simulations and the LDA. In this way, we can perfectly isolate the impact that momentum transfer has on the particle spectra, as in the LDA the "energy packets" do not carry any linear momentum.
To be more specific, we take advantage of the fact that we have access to various performant codes that excel in different areas. To gain an overview over the full particle spectrum, we solve the equations of motion within the DHW formalism and display the results in terms of 2D density plots.
High-precision computations, which are required to perform an in-depth analysis, are done within the Furry picture approach. In order to understand the outcome of these simulations in simple terms, we introduce a model based on energy conservation laws. This absorption model is capable of predicting the appearance of resonances in the particle distribution and thus perfectly suited to examine spectral signatures in the multiphoton regime.
To illustrate our goal, we display in Fig. 6 an example of the momentum distribution in direction parallel to the employed electric field (pulse duration τ = 40m −1 , peak field strength ε = 0.2, and temporal frequency ω = 0.8m). We have chosen this direction specifically because it serves as a good starting point for the following discussion on more complex, fully developed 2D spectra. In particular, because although the background fields' spatial frequency still plays a vital role in creating distinctive patterns in the particle distribution, the net momentum transfer in direction of p x is always zero leading to an easier to interpret distribution function.
In the left panel of Fig. 6, we present the results obtained within the LDA. For comparison, Absorption model.-The key to understanding this apparent mismatch and multiphoton pair production in general lies in proper application of conservation laws; here for energy and linear momentum (it was already successfully employed, e.g., in Refs. [43,65,66]). To derive a conclusive picture, we start stating full 4-momentum conservation for any single multiphoton creation process where p µ e + and p µ e − give the 4-momentum of electron and positron, respectively. The photons are characterized by the 4-vectors k µ + and k µ − , where we used a subscript to distinguish photons with positive and negative momentum with respect to the z axis, c.f. the definition of the vector potential (1). In particular, we have As one can see, the only relevant parameters in our analysis are the fields' temporal and spatial frequency (which coincide with the photon energy ω and photon momentum k z ). 1 We see immediately, that linear momentum conservation for two components are trivial; p x,e + = −p x,e − and p y,e + = −p y,e − . The momentum p z , however, yields a fascinating connection between electrons, positrons, and photons. Formally, we have In a fully symmetric process, n + = n − , the particle-antiparticle-pair is created with zero net momentum as the linear momentum of electron and positron cancel out each other p z,e + = −p z,e − . If, however, n + = n − , particle and antiparticle are accelerated in opposite directions minus an offset q 0 due to a surplus of momentum stemming from the absorbed photons, The interesting aspect of this relation is that although for the latter p z,e + + p z,e − = 0 holds, the whole systems' total linear momentum P tot = particles p is still conserved. This is due to the fact that, because in order to obtain the full picture we always have to take into account all multiphoton channels, in particular the mirrored case n + ↔ n − . As the sum over all particle momenta in the combined situation vanishes the total momentum is preserved (due to vacuum conditions, we have initially P tot = 0). Nevertheless, to keep the discussion concise, we derive all expressions for a single scattering channel ignoring the mirrored cases.
Particle and antiparticle momenta in the z direction are therefor given by and related via q z = q z,e + = −q z,e − . Furthermore, in the weak field limit, ε 1, we can assume that the particle energies are given by The final particle momenta for a specific setup can thus be fully calculated by finding a solution of the equation If, on the other hand, a channel is open, particles are created with a specific kinetic energy. As energy is conserved throughout the creation process any excessive photon energy is converted into higher particle momenta by forming above-threshold peaks. In order to analyze these peaks in the particle spectrum we derive simple expressions within the LDA (k z → 0) and the on-shell limit (k z = ω), respectively.
Particle spectra.-In order to allow for a gentle start, we begin the discussion on the basis of one-dimensional spectra. For the special case where p y = p z = 0 we obtain in case of a vanishing spatial frequency k z recovering the simple energy-momentum relation. Consequently, within the LDA the positions of the above-threshold peaks are determined solely by the total number of photons absorbed n = n + + n − as there is no distinction between n + and n − due to the fact that one cannot distinguish between photons propagating in direction +z and photons propagating in direction −z.
As Fig. 6 has been the showcase for this section, we apply our analysis to these plots first. From the data we retrieve peak positions at p x = 0.65m, 1.24m and 1.73m for the LDA. Employing our absorption model (78) we find that these peaks correspond to our model's prediction regarding 3-, 4-, and 5-photon pair production, respectively. Hence, we label the peak positions p x,3 = 0.66m, p x,4 = 1.25m, and p x,5 = 1.73m.
In the same vein we can perform an analysis for on-shell photons. To be more specific, for k z = ω we obtain finally allowing us to analyze also the right panel in Fig. 6. From Eq. (79) we immediately see, that the total number of absorbed photons n = n + + n − is not sufficient anymore to interpret the particle spectrum. In turn, this automatically explains the rather involved spectrum in Fig. 6 as the various channels become fully distinct. One consequence is that only signatures at high momenta become fully distinguishable. This includes the symmetric channel n + = n − = 2, which is still located at the same position p x,2−2 = p x,4 = 1.25m due to a vanishing linear momentum surplus. Additionally, we have the channel n + = 3, n − = 2 at p x,3−2 = 1.64m and the channel 4 − 2 at p x,4−2 = 1.88m. The broad particle distribution at low momenta p x , however, is the result of summing over various overlapping contributions. In Fig. 8 we have illustrated the particle spectrum with a focus on the net photon number at p x < 0.5m although 2-photon pair production is well below the threshold 2 × 0.8m < m + m.
Furthermore, the right hand side in Eq. (79) yields only imaginary, and thus unphysical values, for one-sided photon absorption. In other words, multiphoton pair production is generally forbidden if the photon is on-shell ω = k z and either n + or n − vanishes. This is in accordance with the fact that an individual plane-wave pulse cannot produce pairs. In particular, one-photon pair production is therefore only possible if the photon is off-shell [44].
On a similar note we can derive analytical expressions for the two different limits in the z component of the final particle momentum, too. For on-shell photons, k z = ω, our model yields which was also found in Ref. [43]. Similarly to Eq. (79), the right hand side diverges for one-sided photon absorption.
On the contrary, in the homogeneous limit, k z = 0, photons only carry energy which, in turn, enables, e.g., one-photon pair production. Analyzing the expression we can further deduce that this includes not only one-photon pair production but also channels where only photons from one beam are absorbed. Furthermore, as (n + + n − ) is equal to the total number of absorbed photons n we have recovered the following simple energy conservation relation given in, e.g., Ref. [67]: The ring pattern revealed in Fig. 7 can be described if we keep both p x and p z in Eq. (77). In particular, we obtain the following relation: This equation indicates that each of the "resonance rings" is, in fact, an ellipse with semi-minor axis b = n + n − ω 2 − m 2 corresponding to the p x direction and semi-major axis a = b(n + + n − )/(2 √ n + n − ) ≥ b regarding the p z direction. The center of the ellipse is located at p x = 0, p z = q 0 . These predictions of our model are in perfect agreement with the full-simulation results.
Resonances.-With the aid of the model introduced above, we can also efficiently search for resonance effects in the spectrum. A similar study has already been performed on the basis of twolevel systems exploiting Rabi oscillations in order to determine resonance frequencies with maximal transition probability between the two states [40,68,69] (in the case of more complex space-timedependent setups, it was done in Refs. [43,52,65]). In terms of our model, such a study can be equivalently formulated as a search for optima in the particle distribution for particles at rest p x = p y = p z = 0. The corresponding optimal frequencies are found to be We point out that such resonances precisely at p = 0 appear only for odd values of n = n + + n − due to charge-conjugation symmetry [40,43,68,69].
To demonstrate conclusively how big the impact of a nonzero spatial frequency is, we illustrate Eq. (84) in Fig. 9. By varying the spatial frequency k z , we can easily calculate peak positions for any Figure 9. Map of the resonance peaks as a function of temporal ω and spatial frequency k z for odd values of the total number n = n + + n − of photons absorbed. The peaks are color-coded according to n (n = 3, 5, and 7).
The thick, solid light-blue line is an indicator for k z = ω.
given absorption channel. The resulting 2D map eventually reveals a highly complex structure of peak locations in the multiphoton regime. Furthermore, it makes it very clear that taking into account the spatial dependence of the background field completely changes the landscape of few-photon pair production. Most important, however, it shows spectacularly that within the LDA we obtain nonphysical contributions. Simply because none of the solid, colored lines in Fig. 9 intersect with the on-shell limiter (light-blue thick line), pair creation, for which either n + or n − is zero, is impossible for on-shell photons. On a similar note, for k z = 0 above-threshold peaks obtained though multiple, different channels can easily be located at the same position in the spectrum leading to nonphysical enhancements.
Such an effect can already be observed in Fig. 6, where the 4-photon channel is misleadingly large in the LDA.
For the sake of completeness, the on-shell limit for this type of resonance is given by The LDA yields local pair production rate. They do not yield any information regarding the total particle number N .
More specifically, only a tiny fraction of electrons/positrons are created at rest at these frequencies as most of the particles are going to be produced with nonvanishing momentum. Moreover, the frequencies ω 0 are not directly linked to the resonance frequencies ω * and thus to the sudden increases in the production rate as the latter are expected to happen when a new channel opens p x = p y = q z = 0.
In order to find these resonance frequencies ω * , we again rely on our model. In this case, energy conservation gives A resonance in the particle yield is then obtained for frequencies and, for the sake of completeness, the resonance condition for on-shell photons is given by Additionally, for light that does not carry momentum, Eq. (88) takes on the more familiar form Hence, only within the LDA the appearance of a resonance in the particle spectrum for particles at rest coincides with the resonance condition for the particle yield, cf. Eq. (86) and Eq. (90).
In Fig. 10 we display this connection for a series of setups. Most notably, in Fig. 10d three different channels are observable (not counting the two mirrored cases). The channels 2 − 3 and 3 − 2 qualify for resonance in the zero-frequency ω 0 . However, the majority of particles created within this channel is ejected in a ∼ 30 • angle from the center. Additionally, although ω 0 exhibits a local optimum due to the 3 − 2 channel, the actual resonance in the particle yield N (ω * ) is expected to come from the 2 − 1 channel. In Fig. 10d at p z ≈ ±0.335m one can already see a new peak to form.
Lastly, one fascinating feature of our setup is that it allows one to study channel closing directly.
Solving Eq. (88) for k z we obtain the resonance condition for the spatial frequency The channels 3 − 2 and 2 − 3 contribute heavily to the particle yield at ω = k z = 0.5m ellipses with center at p = (0, 0.25m) . Their impact diminishes for higher frequencies, because the higher the photon energy the more the produced particles are accelerated and, in turn, the smaller the local rate becomes. At ω = k z = 0.67m the channels 2 − 2 (ring with center at the origin) as well as the channels 1 − 3 and 3 − 1 are clearly visible.
Additionally, around p z = ±0.335m, signatures of 3-photon pair production are already observable.
In Fig. 11 we have displayed the full particle momentum spectrum for a configuration with ε = 0.2, τ = 40m −1 , ω = 0.8m, and various k z . This series of plots strikingly demonstrates the transition from the LDA results to the on-shell limit up to the value of k z where the resonance condition is met.
Technically, for vanishing k z we cannot distinguish between the 2 − 1 and the 1 − 2 channels.
Consequently, both form the same structure in momentum space centered at p = 0 even constructively interfering with each other. We also see, that an increase in the spatial frequency automatically leads to an increase in the net momentum transfer as n + = n − . Closing of the channels 2 − 1 and 1 − 2 is predicted at k * ≈ 1.33m explaining the tremendous decrease in the distribution function in Fig. 11f.
The drop-off is apparently so significant that it leads to the 2 − 2 channel being dominant.
Beyond the weak-field approximation.-In this whole section, we have never mentioned shifts in the final particle spectrum due to the oscillatory energy of the created particles. The reason is, that for the exemplary setups presented here such a field-dependent effect amounts to a ∼ 1% correction only. Nevertheless, in order to build a bridge to the established results given in the literature we have included a discussion on the effective mass effects in Appendix B.  Figure 11. Density plots of the particle momentum spectrum n (p x , p z ) (p y = 0) for various spatial frequencies k z . For otherwise identical field configurations (ε = 0.2, τ = 40m −1 , ω = 0.8m), a higher frequency in spatial direction k z results in a shift towards higher particle momenta. Going from low to high frequencies the signatures also become smaller in diameter. At a photon momentum of k z = 1.0m, a great amount of particles is created at rest, p = 0. In the last plot (k z = 1.33m), the channels 2 − 1 and 1 − 2 are close to the resonance condition.
Additionally, the interference pattern per shell structure strongly varies with k z .
Electron-positron pair production within high-intensity electromagnetic background fields is a highly complex process, thus adequately precise theoretical predictions for realistic scenarios are difficult to obtain. Specifically, the fact that phenomena concerning the particle distributions occur within laser fields requires resolving effects on two vastly different scales (1/m ∼ 10 −21 s and 1 fs ∼ 10 −15 s), so a complete simulation is out of reach for the foreseeable future. In this regard and especially with a view to upcoming strong-field laser facilities, it is thus of utter importance to develop a formalism that can accurately describe such strong-field experiments.
In this work we have therefor inspected the applicability of the local dipole approximation (LDA) in the context of short-pulsed fields and in the regime of few-photon pair production. Within our study, we have found a second criterion for the applicability of the LDA based on the expected particle movement within the background field; the first involving the formation length of the particle pair.
To be more specific, in order to have the LDA to hold, the particle trajectories after creation must be confined to an area smaller than the fields' wavelength. Note, however, that neither the formation length nor the particle movement make hard, exclusive statements as we have still found setups where only one of these two criteria was fulfilled.
In the second part of the manuscript, we have made a comparison between the LDA and full simulation in the regime of multiphoton pair production on the basis of the particle momentum spectrum.
We have reproduced the spectra given in the literature for the LDA confirming that these momentum distribution functions indeed carry signatures of particle absorption. However, upon further inspection, these distributions fail to give the right above-threshold peak positions (except for the symmetric case, where an equal amount of photons was absorbed from each direction) and even show unphysical contributions. Moreover, quantum resonances show up at wrong frequencies and momenta and, worst, the amplitudes in the spectrum are wrong by at least one order of magnitude. The latter being true even for symmetric, few-photon signals, where enhancements stemming from quantum interferences of different scattering channels are rare.
All limitations of the LDA in this regime can be traced back to its inability to take into account the photon's linear momentum. As a matter of fact we have developed an effective model based upon the conservation of energy and linear momentum in order to predict peak locations and to obtain the right resonance conditions. While this model does not yield any information on the signal strength, it can give an idea about the landscape in the multiphoton regime. This includes estimates for quantum resonances, channel-resolved peak locations as well as zero-momentum resonance frequencies. Although the model presented here in this manuscript is based on monochromatic, linearly polarized fields, the basic concept is universal and can be readily applied to more involved field configurations.
be different. Within this manuscript, we have employed one solution technique based on Furry-picture quantization uniquely suited to calculate particle distributions to an unprecedented level of accuracy.
The alternative way is given by the DHW formalism, a relativistic quantum kinetic approach, that automatically yields, by its design, a complete overview of the particle phase-space and therefor the full momentum spectrum.
As was demonstrated in the manuscript multiple times, the strengths of these two significantly different formalism lie in totally different areas; one is about high-precision while the other aims at completeness. It is therefor challenging to show that both approaches indeed produce the same result. In this section of the article, we want to provide some insight into the technical details of our computations.
To do so, we employ a test scenario that allows us to examine in particular the accuracy of the results displayed in the main text.
The test configuration is given by a long-pulsed field with the pulse length τ = 40 −1 m, peak field strength ε = 0.2, and temporal frequency ω = 0.8m. For the spatial frequency we opted for two different scenarios. In the first setup, we studied on-shell photons, k z = ω = 0.8m, ultimately proving trustworthiness of both approaches, see Fig. 12. In the second scenario, c.f. Fig. 13, we have averaged over z allowing for inspection within the LDA.
The Furry-picture formalism is formulated in terms of solving many first-order, ordinary differential equations in time. As memory is not an issue, the only concern is about resolving any given oscillation in the mode functions. As the approach is used sparsely in the manuscript mainly in order to determine distribution functions locally, there is no need to use coarse step functions in order to save computation time. Consequently, the usage of a small step size in time t leads to tremendously accurate solutions in the particle momentum spectra. In Figs. 12 and 13 the actual limiter was given by the precision of long double variables ∼ O (10 −15 ).
The DHW approach, on the other hand, operates intrinsically with partial differential equations.
Hence, accuracy and precision of such an approach cannot compete with more streamlined formalism using, e.g., mode decomposition. However, given sufficiently advanced computational techniques and sophisticated parallelization schemes, phase-space approaches can make up for their slowness in calculating distributions at a specific point in the domain by the huge number of output they can produce. On top of that, although kinetic approaches are not on the same level as other formalism in terms of accuracy, they still provide reasonably precise results, in Fig. 12 the orange curve evens out at ∼ O (10 −9 ).
Within the LDA the situation changes significantly. Due to the fact that within this approximation, Figure 12. Comparing the particle distribution function n p as a function of the momentum p x (p y = p z = 0) obtained through the Furry-picture quantization approach (blue line) and DHW formalism (orange markers). The signal-to-noise ratio for the DHW formalism is ∼ 10 −5 , while the Furry-picture technique is accurate to more than 8 orders of magnitude. At peak value around p x = 0.32m, both approaches deviate by roughly 3%. The configuration tested here is: τ = 40 −1 m, ε = 0.2, ω = k z = 0.8m.
spatial inhomogeneities are singled out altogether, the coupled system of partial differential equations governing the particle dynamics in the DHW formalism is converted into a set of ordinary differential equations. Consequently, the time integration of this noncoupled system behaves similarly to the mode equations within the Furry-picture formalism leading to highly-accurate, easily calculable distribution functions. In Fig. 13 a comparison between the Furry-picture formalism and DHW approach in the LDA is displayed.

Appendix B: Effective mass effects
To put it simple, the sole idea of having an effective mass is based on the concept of summarizing the electrons' and positrons' collective interactions with the environment and absorbing it in a single variable. If successful, this approach leads to a massive simplification of the process one wants to describe. It is therefor not surprising that the concept of an effective mass is very often used in the Figure 13. Comparison between Furry-picture quantization and DHW formalism in the LDA on the basis of the particle distribution function n px,py=0,pz=0 . Both approaches produce results that are accurate to at least 8 orders.
literature when discussing physics in the multiphoton regime. In its most common form, it appears as a field-dependent modification m → m * (E) with primary focus on the particles' oscillatory energy within a high-frequency field.
The huge advantage of this build is its simplicity. In its most uncomplicated form, the effective mass is easily derived through a single averaging over the temporal shape of the vector potential. Reference [67] is an excellent example of such an approach, where the field-dependent mass term is found to be For the vector potential in this manuscript (1), this would amount to m * ≈ m 1 + ε 2 m 2 2ω 2 (B2) completely neglecting any spatial dependency.
A slightly different approach is given in Refs. [43,65,75], where the effective mass is not calcu- (ring with center "x" at origin), 1 − 3 and 3 − 1 (ellipses with center "+" at p z = ±0.67m) as well as 1 − 2 and 2 − 1 (ellipses with center "o" at p z = ±0.335m). Field parameters: ε = 0.2, τ = 60m −1 , and ω = k z = 0.8m. lated directly but through an effective particle energy defined via This technique is more accurate than the previous one and applicable even for very strong fields ε ∼ O (1). The downside, however, is given by the fact that in order to calculate peak locations (78)-(81) or resonance frequencies, Eqs. (84) and (88), one has to solve an implicit equation.

Appendix C: Supplements
We could not put all data we have generated in the main body of the manuscript without seeming to be repetitive. The utility of the following attachments is too narrow to warrant putting it into the main body of the manuscript. Nevertheless, they might still turn out to be of good use. Table I. Supplemental material for the figures displayed in Sec. VI B; List of peak locations in the particle spectrum according to the exact numerical simulation p x and our model p x,n + −n − , where n + (n − ) gives the number of photons absorbed with positive (negative) momentum. Effective mass effects were not taken into account in Sec. VI B, thus the peak location for the channel 2 − 1 is slightly off. a The background field configuration is determined by a temporal length of τ = 40m −1 , a peak field strength of ε = 0.2, and a frequency of a If effective mass effects were taken into account, we would get p x = 0.325m.
b This configuration exhibits a plateau rather than a peak.