On the effect of time-dependent inhomogeneous magnetic fields on the particle momentum spectrum in electron-positron pair production

Electron-positron pair production in spatially and temporally inhomogeneous electric and magnetic fields is studied within the Dirac-Heisenberg-Wigner formalism (quantum kinetic theory) through computing the corresponding Wigner functions. The focus is on discussing the particle momentum spectrum regarding signatures of Schwinger and multiphoton pair production. Special emphasis is put on studying the impact of a strong dynamical magnetic field on the particle distribution functions. As the equal-time Wigner approach is formulated in terms of partial integro-differential equations an entire section of the manuscript is dedicated to present numerical solution techniques applicable to Wigner function approaches in general.


I. INTRODUCTION
The creation of matter out of the vacuum is one of the most exciting concepts of high-energy physics. In particular strong-field quantum electrodynamics (QED) is perfectly suited to study the generation of matter through energy, because it provides a comparatively clean setting [1][2][3].
The main issue is to create the right laboratory conditions as it takes extremely strong field strengths in order to see any form of signal of quantum vacuum non-linearities [4][5][6][7]. However, as laser technology has significantly advanced over the last years probing the quantum vacuum in an earth-based laboratory seems to be in reach [8][9][10]. In fact, the research field has recently gained interest to the extent that upcoming laser facilities already prepare for high-intensity experiments [11][12][13].
In the context of particle creation, it has been shown already that matter can indeed be created in a leptonphoton collider. It was demonstrated that highly energetic electrons, if subjected to huge electric fields, emit photons in the multi-keV regime which, in turn, interact with the light beam itself eventually forming an electron-positron pair [14,15]. On the other hand, the Schwinger effect the depletion of the quantum vacuum in constant background fields, has been predicted 90 years ago [3], but still remains untested. Only recently it was suggested that both mechanism could be utilized together in a multi-beam scenario leading, in theory, to a tremendous increase in the creation rate [16]. For further reviews on pair production in general see Refs. [17][18][19].
Due to the highly complex nature of the subject, only a few analytical results can be found [20][21][22]. Hence, we have to rely on numerical approaches in order to describe pair production accurately. Over the timespan of approximately 20 years we have seen the rise of an abundance of computational methods, all performing very well in certain parameter regions [23][24][25][26]. One common trait, however, was that they * c.kohlfuerst@hzdr.de worked best for homogeneous electric fields as it was argued that the chances to produce particles should be highest in the vicinity of the focus of two crossed beams where the magnetic field vanishes. In particular the so-called quantum kinetic theory (QKT) proved to be very successful under these special circumstances [27][28][29][30]. Only recently, numerical solvers have been improved to the point that it is now possible to take into account the spatial inhomogeneity of laser beams [31][32][33][34][35][36][37][38] beyond any semi-classical approximations [39][40][41]. Performing calculations with these tools revealed additional features of strong field pair production; particle self-bunching [32], ponderomotive effects [42] and spin-field interactions [43]. Furthermore it was shown that disregarding the spatial finiteness of a laser pulse can lead to spurious effects in the particle spectrum [44,45].
In order to fully account for the spatial inhomogeneity of background fields we rely on the so-called Dirac-Heisenberg-Wigner (DHW) formalism [46][47][48]. Its major advantage is its flexibility, because the formalism allows to study the pair production process for any given background field. Once the structure of electric and magnetic fields has been established the formalism acts as a black box automatically incorporating all possible production and interactions mechanisms. Furthermore, as the DHW formalism is deduced a priori from a Lagrangian level no models or assumptions are needed. The only drawback is that a computational solver based on the Wigner formalism is generally very resource-hungry, as one solves for the time evolution of a multi-dimensional partial differential equation. To overcome these issues in order to make simulations feasible novel computational techniques have to be applied.
Although various articles on the Wigner formalism and its applications have been published in recent years [49][50][51][52], a compendium containing methods and techniques to successfully treat the corresponding equations of motion with computational methods is not available. The intention of this article therefore is to introduce various computational solution techniques that are applicable to solving them. In this context, we discuss caveats of the formalism and bring up ways to overcome them. This is done in addition to discussions on the solvers' effectiveness arXiv:1912.09359v1 [hep-ph] 19 Dec 2019 and overall performance including considerations on noise suppression.
We apply all these novel techniques in order to compute pair production in a regime where tunneling as well as absorption effects are important. In particular, we study the impact the magnetic field has on the production yield and thus on the particles' momentum spectrum. Furthermore, as the whole creation process is very sensitive to variations in the background field, e.g., a small change in one parameter can mean a different creation mechanism becomes favored thus heavily altering the spectrum, we rigorously analyze slightly different field configurations eventually improving our understanding of magnetic field effects in particle creation scenarios in general.
In short, this manuscript is organized as follows. We first introduce a quite general model of a background field, which serves as the connection point for the discussion on the characteristics of the various production mechanism, see section II. In section III we introduce the DHW formalism and in section IV we describe how to solve the governing set of equations of motion efficiently. In the main part of this manuscript, section V, the numerical results are discussed. Eventually, a conclusion is given in section VI.
Throughout this paper we use natural units = c = 1 and express all quantities in terms of the electron mass m (= m e + = m e − ).

II. MOTIVATION
A virtual electron-positron pair is formed due to the omnipresence of vacuum fluctuations, where the scale is determined by the Compton time τ C = 1/m e − ≈ 10 −21 s and Compton length of an, in this case, electron λ C = 1/m e − ≈ 10 −12 m.
The main task is to prevent these two particles from recombining. In this way, a virtual pair is turned into two separate particles. In order to do so, a strong electromagnetic field is employed acting on the charge carriers. There are now two different ways to pump energy into the system. Either via providing enough work to allow for tunneling (Schwinger effect) or due to the absorption of energy given by the individual photons forming the background field (multiphoton pair production). These effects are not mutually exclusive, thus depending on the background field any combination of mechanisms is possible.
The probabilities for such an event to happen scale with the critical field strength It is defined as the minimal work that has to be done by the field over one Compton length in order to create one pair 2eE cr λ C = m e − + m e + . The factor of 2 is due to the fact, that electrons and positrons have opposite charge thus they are accelerated in opposite directions, too.

A. Model for the fields
The ultimate goal of this article is to study the consequences of applying an additional magnetic field to the particle creation process and investigate how the different creation mechanism are affected by it. As we want to smoothly interpolate between Schwinger and multiphoton pair production we employ a toy model of the form of (2) and derive electric and magnetic fields accordingly Here, ε gives the peak field strength in terms of the critical field strength m 2 /e, ω states the field frequency and τ describes the pulse length. Moreover, we have decided to use a variable envelope function, exp − t j τ j , in order to have better control over the number of cycles in the field. This choice makes it indeed very easy to switch from few-cycle to prolonged many-cycle pulses, c.f. Fig. 1. The parameter λ, on the other hand, controls the extent of the peak in spatial direction and indirectly also defines the peak strength of the magnetic field, where a smaller λ accounts for higher magnetic field strengths. Technically, Eqs. (2)-(4) describe a head-on collision of two identical, linearly polarized laser beams within a quasi-dipole approximation. To be more specific, we have two non-propagating pulses e ±ikzz ≈ 1 + O (k z ) each consisting of left-and right-handed waves that form a temporally and spatially localized interference pattern mimicking the region of highest intensity in a realistic setup. In this context, we can therefore interpret the background field as incoming spherical waves with angular momentum j = 1 (electric dipole). Moreover, the electric field (3) still accounts for the individual photon energies in the form of the field frequency ω.
Naturally, one were also to assume that photons would transfer linear momentum to the created particles. Our model (2), however, does not take into account the photon's linear momentum as we are only interested in basic particle behavior in presence of a strong magnetic field. In order to allow for a non-zero photon momentum one would have to go beyond the dipole approximation and examine pair production in, e.g., a standing wave One of the reasons we have decided to not use such a model is that a standing wave pattern creates the opportunity for simultaneous pair production at a multitude of locations. This would render any analysis on specific details of the production process unclear, because in a phase-space approach every bunch of particles would be stacked on top of each other leading to a confusing particle spectrum. Given the definition of the background field (3)-(4) we can already discriminate between the Schwinger effect and multiphoton pair production using the Keldysh parameter [53] γ = mω eε 1 Schwinger 1 Multiphoton (6) As γ was first introduced for infinitely long fields it does not take the temporal envelope into account. In our case, this means that in addition to the choice of ε and ω the shape of the envelope function is a major contributing factor.

B. Regimes of Pair Production
In the following, we will briefly introduce the different regimes of pair production (Schwinger, multiphoton and intermediate) and discuss their characteristic signatures in the momentum spectrum.

The Schwinger effect
The Schwinger regime applies when a strong background field is employed over a sufficiently long period of time in order to allow the particle to tunnel through the Coulomb barrier. As we work exclusively with linearly polarized fields we can adopt the notion in Refs. [54,55] (describing the conceptually similar qq-string breaking) yielding where a = 1 The key aspects of the probability rate P (t, z, p z ) are, that (i) particles are produced at times when the electric field is strong and (ii) that these particles have zero initial parallel momentum p x . A non-vanishing initial transversal momentum p z is possible. However, it comes with a penalty in terms of likeliness, see the exponential term m 2 + p 2 z . After creation, these particles are exposed to strong forces due to the extreme strength of the background field. As we are primarily interested in the momentum spectrum at final times (when the electromagnetic field is switched off) we have to account for these fieldparticle interactions. The simplest way to achieve such a description is by applying a single-trajectory formalism, where electrons/positrons are considered to be semiclassical point-like particles, which follow classical paths [56]. Excluding any kind of interaction between generated particles, thus naturally ignoring quantum interferences and particle collisions, these trajectories can then be tracked and evaluated.
We can already deduce that in case of a quasihomogeneous configuration and therefore a vanishing magnetic field the particle trajectory can be easily (11) p z (t) = p z,0 = const, (12) where t 0 gives the particle's time of creation.
If the temporal structure of the electric field exhibits multiple peaks it essentially acts as a double-slit in time reflecting the particle's quantum statistics, see Refs. [60,61] for a more profound discussion. Alternatively, one can interpret the final distribution as the sum over trajectories with the particles initial time of existence as starting points [43]. Quantum interferences are then obtained through summing up the quantum phases these particles have picked up when exposed to the background field [62].
If a strong magnetic field is present, however, these trajectories are altered significantly. Speaking in terms of Eqs. (8)-(10), as soon as a particle has acquired a parallel momentum p x , the magnetic force automatically converts it to transversal momentum p z . As such the particle is pushed into direction ±z. This happens the quicker the stronger the magnetic field is. Interestingly, if the particle leaves the strong-field region at times of high acceleration, it cannot be pushed back once it reaches regions of low field strengths. As a result this particle is basically unaffected by variations of the background field at later times.
The spin force term F S = s ∂ z eB (t, z(t)) in Eq. (10) is an interesting addition to the classical Lorentz force equations, because it allows for symmetry-breaking phenomena. In the vein of Ref. [43], we have decided to work in 2+1 dimensions in order to keep the volume of the phase-space on a numerically manageable level. As a side effect, there are multiple ways to define QED (reducible or irreducible representation). To keep the calculations as simple as possible we have opted for one of the 2−spinor representation. As a result, the spin in our calculations is fixed, thus a spatially rapidly varying magnetic field can indeed break the symmetry in transversal direction.

Multiphoton pair production
Contrary to the Schwinger effect, multiphoton pair production relies on a sufficient energy transfer from photons to virtual pairs. Moreover, electrons and positrons quiver within an oscillating background field resulting in an increased potential energy. Only if the amount of transfered energy exceeds the pairs rest energy, a particle and an antiparticle is formed. As energy has to be conserved, the surplus of photon energy is directly converted to the particle's initial kinetic momentum. All in all, we can define the effective energy of a particle as [63] E (p x , p z ) = ω 2π Correspondingly, energy conservation reads 2E (p x , p z ) = nω with n being the number of absorbed photons.

Intermediate regime
For few-cycle pulses with rapidly varying timedependency or for background fields exhibiting a Keldysh parameter of γ ≈ 1 neither Schwinger nor multiphoton effects dominate. In such a case either no distinct signatures of Schwinger nor multiphoton effects are visible or the particle spectrum shows characteristics of both. Conceptually, it is, for example, easily possible, that a virtual particle pair absorbs n photons and is simultaneously pulled apart by the strong background field resulting in interferences showing up in the particle spectrum. Interactions of this kind are difficult to calculate, thus they make the perfect case study to demonstrate the power of kinetic approaches.
On a side note, the so-called dynamically assisted mechanism, where photon absorption is used to effectively lower the threshold for pair production, is quite similar. The main difference is that the enhancement effect is much more pronounced, because multiple fields are fine-tuned in order to give a strong signal in the spectrum [16,[64][65][66].

III. DHW FORMALISM
The DHW formalism is a very general method that allows to study particle creation within arbitrary electromagnetic fields. In the context of this article, this means we only have to develop one single numerical algorithm to study the spectra produced through fewcycle as well as many-cycle pulses.
As providing a detailed introduction into phase-space methods in general is not the goal of this manuscript, we only state the most important steps in deriving the DHW formalism here. We recommend Ref. [46] for a more profound introduction into quantum transport theories. Additional information on quantum kinetic theories in the context of particle creation can be found in Refs. [27,29]. A complete derivation of the equations of motion has been performed in Ref. [67].
To save computational costs later, but without disregarding magnetic fields, we introduce the QED Lagrangian in 2 + 1 dimensions where D µ = (∂ µ + ieA µ ) and Additionally, we have the vector potential A µ , the electromagnetic field strength tensor F µν = ∂ µ A ν − ∂ ν A µ and the spinor fields Ψ andΨ. In an attempt to further reduce the computational overhead we use only one 2spinor basis 2 effectively cutting numerical costs in half. Varying with respect to ψ (ψ) we obtain the (adjoint) Dirac equation The backbone of the formalism is the density operator where we introduced the center-of-mass coordinate r and the relative coordinate s. The Wilson line factor is essential to guarantee gauge invariance. The covariant Wigner operator is then obtained via a Fourier transform with respect to the relative coordinate ŝ In the following we will drop the indices inŴ to improve readability.
In a rather lengthy calculation, see Ref. [46,67] for details, we can combine the Wigner operator (20) with Eqs. (16) and (17) to generate the operator equations 2 We perform all calculations in the xz-plane, thus we label the third matrix γ 3 .
with the pseudo-differential operators Due to the operator nature of the equations above, it is difficult to use them in a numerical simulation.
In order to obtain computational feasible equations of motion we therefore introduce a mean-field (Hartree) approximation Consequently, the operator-valued electromagnetic field strength tensorF µν is treated as a C-number field F µν . In fact, a direct evaluation of Φ|F µν (r)Ĉ (r, s) |Φ would introduce couplings to an arbitrary number of n-body terms (BBGKY hierarchy). The Hartree approximation is in principle a truncation at one-body level effectively turning an infinite hierarchy of coupled differential equations into a closed system, c.f. Ref. [47] for an in-depth analysis.
We proceed by taking the vacuum expectation value of Eqs. (21) and (22). As a result of the Hartree approximation Φ|F µν (r)Ĉ (r, s) |Φ ≈ F µν (r) Φ|Ĉ (r, s) |Φ , (26) Eqs. (21) and (22) are transformed into equations of motion for the covariant Wigner function Working directly with the covariant Wigner function is problematic, because it requires solving the equations of motion at all points in time at once. This is certainly not feasible in a real-time simulation, thus we choose a definite frame to project on equal-time effectively reformulating everything to an initial-value problem. In this way, we obtain the equal-time Wigner function which we decompose into Dirac bilinears This makes it easy to interpret the quantities s as mass density, v 0 as charge density and v 1 , v 3 as the current density.
Finally, we obtain a coupled set of differential equations for the equal-time Wigner coefficients, cf. Refs. [67], with the pseudo-differential operators The initial conditions are given by the vacuum solution

IV. NUMERICAL SOLUTION TECHNIQUES
One of the aims of this manuscript is to demonstrate how to solve the equations of motion (30)-(33) effectively. Although a brute-force calculation might be successful, running such simulations is computationally expensive and numerically quite unstable. To combat these issues we have developed various solution techniques that certainly help to speed-up the calculations while reducing the computational load and still maintaining a high accuracy. It is worth mentioning that all methods and numerical recipes introduced in this manuscript can be easily applied to other system (QED 1+1 , . . .), too.

Preconditioning
In order to solve the differential equations (30)-(33) for a configuration featuring fields of the form of Eqs. (3) and (4) multiple preconditioning steps are in order. At first we introduce modified Wigner components [32] w v = w − w vac (40) such that we obtain a set of inhomogeneous partial differential equations with vanishing initial conditions. In this way, we pro-actively ensure that all components fall off at the boundary already minimizing truncation errors.
Numerical feasibility can then be significantly improved by establishing the relation where A(t, z) is the vector potential defined previously (2). In this way, we take into account so-called "minimal coupling". More specifically, we switch from kinetic momenta p to canonical momenta q. Naturally, the derivative terms have to be transformed, too The most prominent advantage of this transformation is the change in the time-development operator As is immediately obvious from a Taylor expansion of the derivative operator where E (2n) denotes the 2n-th derivative with respect to z, the leading term in Eqs. (44) is of the order of O E (t, z) . Hence, quasi-homogeneous quantum kinetics is already incorporated in the transformed equations. In particular, in the limit of a locally constant vector potential, E (t, z + iξ∂ qz ) ≈ E(t, z) and B (t, z + iξ∂ qz ) ≈ 0, the momentum derivatives vanish, thus one is left with a partial differential equation in t and z. If one employs a local constant field approximation (LCFA) the derivatives with respect to z vanish, too [32]. Eventually, we change the geometry of the domain we work in. The equations (30)- (33) are defined on an open domain z ∈ (−∞, +∞). Due to the change of variables (40) and the fact, that we use localized fields (3) the region of interest is clearly finite. Hence, we truncate the length of the domain z ∈ [−L z , +L z ] to minimize memory requirements. The same holds for the momentum variables q x and q z .

Background fields
In the code, we only specify the time-dependence of the electric field analytically. Time dependence of the magnetic field and thus of the vector potential are calculated numerically by adding one term of the form of dA dt = −E(t) to the system of differential equations.
In previous studies [32,43] an analytical form of the differential operators has been used. Such an approach is highly inflexible and one has to deal with cancellation effects in the anti-derivative. A better way is to perform the integration with respect to ξ numerically and store the result in a lookup table.

Spectral solver
In order to solve the equations of motion a combination of a spectral solver (phase-space) and a finite-difference solver (time dependency) is used. To be more specific, we employ pseudo-spectral methods at every time step to achieve optimal accuracy at minimal grid size in the phasespace domain. The general idea of using a pseudo-spectral solver is to solve the system of equations in phase-space, but evaluate derivatives via multiplications in Fourier space, c.f. Fig. 2. In this way, exponential accuracy can be achieved [68]. Time integration is done using a Dormand-Prince Runge-Kutta integrator of order 8(5,3) [69].
A generic derivative is calculated in the following way wheref (k) denotes a Fourier transformed quantity. Derivatives with respect to z and q x are evaluated in this way. Pseudo-differential operators (34)- (38) and therefore derivatives with respect to q z are evaluated differently. For a generic function G (z + iξ∂ qz , t) = G t (t) G x (z + iξ∂ qz ) we write In contrast to previous studies [42,43] this technique has been applied to evaluate derivatives of the inhomogeneities, too. The use of fast Fourier transforms (FFTs) accelerates the simulations tremendously [68,70]. Only downside is that we have to have Fourier basis functions namely sin (kx) and cos (kx). This choice of basis functions requires (i) an equidistant discretization of the domain and (ii) periodic boundary conditions. Alternatively, Chebychef polynomials can be used, which have the advantage that they can be used together with a nonperiodic grid. However, they demand a discretization where the point density is highest at the boundaries. As we wanted to have a high resolution at the center of the domain to capture as many physical effects as possible we chose equidistant point sampling.
A further improvement on lattice point distribution is given by the fact, that we can in principle introduce any bijective transformation of the Eqs. (34)- (38) and  then solve a similar set of equations of motion. The solution to the original problem is then simply given by re-transforming the results once the simulation has finished. In our case, we introduce the transformations where L q (L z ) gives the length of the domain in q x (z)-direction and α q (α z ) controls the strength of the deformation. In this way, we can (i) achieve a much better convergence rate and (ii) we are still allowed to use a customized rasterization of the domain. The derivative operators transform accordingly Keep in mind, however, that the grid inz andq x is regular, the grid in the physical domain (z, q x ) is not.
Regarding periodic boundary conditions, we observe that due to the transformation to modified Wigner components (40) all quantities fall off for larger values ofz,q x or q z . We can thus truncate the domain at points where the coefficients supposedly vanish yielding a finite domain of size Lz × Lq x × L qz . Then we setz 1 =z Nz ,q x,1 =q x,Nq x and q z,1 = q z,Nq z , where N z , N qx and N qz give the total number of grid points inz,q x and q z , respectively. While momentum and spatial derivatives have to be evaluated at every single time step (often multiple times), the overall time evolution can be calculated in a straightforward way, see Fig. 2. First, a starting point t vac of the simulation has to be determined. Ideally, the vector potential as well as the fields are still zero at this point such that one is left with the pure vacuum state. However, background fields may have infinite support making a truncation of the time domain necessary. To minimize the error due to such a truncation, while not needlessly increasing computing time, we have set the initial times to t * = −6τ (Gaussian envelope) and t * = −2τ (super-Gaussian).
The actual time stepping t i+1 = t i + ∆t might be done through explicit or implicit methods. We have chosen a higher-order explicit stepper, because it seems to be the ideal compromise between accuracy and run time. Specifically, we neither have to allocate additional memory due to the overhead in terms of inverting the derivative matrix nor do we have to solve the corresponding algebraic equations. On the other hand, we still obtain high accuracy using an 8th-order method with error-correction.
As we are only interested in the particle's momentum spectrum, simulations terminate once the background fields vanish, here t * = −t * .

Filters
Due to the fact, that we want to apply the procedure introduced in Eq. (47) only trivial transformations in q z are feasible. Consequently, we have to use an equidistant grid in q z , thus we are automatically limited to a small domain L qz . In turn, we have to account for the fact that the Wigner coefficients may not fall off sufficiently fast, thus the requirement for periodic boundary conditions could be in conflict with domain truncation. More specific, there could be a significant discrepancy between coefficients defined at the now neighboring points q z,Nq z −1 , q z,1 and q z,2 resulting in Gibbs phenomenon induced effects, e.g., an artificial growth of the coefficients w.
To avoid the creation of spurious patterns or a wrong sampling of high-frequency modes so-called anti-aliasing methods are employed, see Refs. [68,71] for a more detailed analysis. First, we introduce an artificial dampening factor at every time step, which is fairly easy to implement as we only have to modify the prefactor of the ∂ qz -derivative to with N w ∈ N.
Secondly, to resolve any additional problems regarding high-order modes we cut the upper part of the Fourier spectrum at the end of a simulation (at the last time step t f , but before evaluating any observables). The combination of applying these filters proved to work excellently, as it was possible to avoid the appearance of any non-physical terms in the final distribution functions. On a side note, in previous calculations [43] the highest order modes have been terminated at each time step. However, that technique seemed to fail when applied to spatially strongly inhomogeneous, long pulsed fields as spurious artifacts started to show up at long run times.
Unfortunately, for configurations that require a long runtime, none of these methods work sufficiently well. In such a case it seems as if the only even remotely successful way is to sacrifice some resolution in direction of q z and expand the total domain size. As the source of these artificial oscillations is the jump at the boundary, any artifacts associated with it will start at ±L qz and then propagate towards the center. If the domain, however, is large enough, the time it takes for these oscillations to reach the real particle distributions is so long that we can safely stop the simulation at times t f , where we still have an intact particle signal. Of course, before evaluating any observables the Wigner components at all higher transversal momenta have to be artificially set to zero.

Observables
The DHW formalism only allows to solve for the Wigner coefficients w. As we are, however, mainly interested in observables we have to rearrange these coefficients to obtain particle distribution functions. To keep the discussion simple, we only state the particle number density at asymptotic times t f (A(t → ∞, z) → 0) Accordingly, the particles momentum spectrum is given by as well as The total particle number is then given by Note, that we have already transformed the data obtained from the simulation to ordinary spatial coordinates and canonical momenta {z, q x , q z }, because it is easier to discuss the results in terms of a physical basis. The crucial point in evaluating Eqs. (53)-(57) is given by the fact, that the particle signal is generally only acquired after summing up all relevant Wigner components. The reason is, that each of s v , v v,1 and v v,3 not only holds the information of the particle spectrum but also the information for additional observables, c.f. Ref. [46]. As a consequence, the fraction of the particle spectrum on a Wigner component might be lower than 1%. Hence, when calculating particle distributions one technically extracts sub-leading contributions, which turns this final step into a delicate process.
As the particle yield is naturally higher the further the electric field is extended we have to normalize the rates to allow for a fair comparison between results with different λ. The problem is that in principal Schwinger effect, multiphoton pair production and assistance mechanism all demand different normalization. The compromise was to normalize parallel f (p x ) as well as transversal f (p z ) distribution functions in terms of semi-classical expectation values with a = 1 2 |B 2 − E 2 | − (B 2 − E 2 ) due to its simplicity and its capability of taking into account the electric field's spatial finiteness as well as the magnetic field's attribute to suppress pairs to form. In the limit of pure Schwinger pair production Eq. (58) is even exact. Moreover, signatures of multiphoton particle creation generally extent non-linearly in transversal direction thus we display such results in 2d contour plots, where no particular normalization is needed.

V. RESULTS
One big goal of this manuscript is to provide some insight into particle creation and, consequently, particle dynamics in high-intensity electromagnetic background fields. We discuss the final particle momentum spectra with respect to the different creation mechanisms and show how to interpret certain features by simple means. On a side note, the full set of data in terms of contour plots of the particle distribution function f (p x , p z ) is attached at the end of the manuscript, c.f. section VII.

A. Tunneling-dominated pair production
To obtain a clear qualitative picture of particle creation via tunneling we employ slowly varying, few-cycle fields with a peak field strength of eε = 0.5m 2 as well as eε = 0.2m 2 , respectively.
Both field configurations have in common that particles are predominantly produced via tunneling. Reason is that, although the Keldysh parameter (γ = 1, γ = 0.4) might indicate that both particle creation mechanism are important, the total pulse length is quite short τ = 25m −1 and the envelope function does not allow for many subcycles. As a result multiphoton effects are suppressed, thus the spectra can be nicely described through a semi-classical single-trajectory picture. Nevertheless, the temporal variation in the field is sufficiently strong such that one cannot observe the "pure" Schwinger effect, see the deviations in the normalized total yield in Fig. 3.
The simple man's model introduced in Ref. [43] states that particles are produced in regions where E(t, z) 2 − B(t, z) 2 > 0 holds (if E·B = 0) and the higher the effective field strength E(t, z) 2 − B(t, z) 2 the higher the chance for pair production. For a field with a Gaussian envelope, pulse length τ = 25m −1 and frequency ω = 0.2m only three peaks are really capable of producing a sizable amount of particles. More specifically, given that z = 0 the field is strongest at t = 0, while at t = ±13.7m −1 it reaches ∼ 75% of its maximum value. Due to the Schwinger effect's exponential suppression by the effective field strength every other peak can be considered as minor, thus they can be safely neglected in the further discussion.
In a quasi-homogeneous calculation (λ = 1000m −1 ), see the solid blue curves in Fig. 3, we obtain a broad peak superposed by oscillations in f (p x ) and an exponentially declining distribution in p z . Assuming, that for this configuration we can neglect the magnetic field and consider the electric field as homogeneous in z, we can solve the set of equations (30)-(33) analytically. As particles in the simple man's model (7)-(10) are created with zero initial parallel momentum p x,i the final particle momentum p x,f solely depends on the strength of the vector potential at the time of creation t 0 . Hence, we can easily calculate the reference points for the final distribution. Unsurprisingly, pair creation at t 0 = 0 translates into a peak at p x,f = 0. At t = ±3.25m −1 the electric field still has 75% of its maximal strength translating into a 90 % smaller chance for pair production (7). Consequently, if a particle was created at t 0 = 3.25m −1 it would obtain a final momentum of p x,f = 2.98 ε/E cr m, thus roughly determining the point where the particle density has fallen off by 90 %. Particles created at the two side maxima of the electric field (t = ±13.7m −1 ) acquire a final momentum of p x,f = ±1.44 ε/E cr m. Hence, particle bunches stemming from main and side peaks in E(t, z) are clearly overlapping. Moreover, one would expect these particles to carry different phase information due to their different times of creation, c.f. Refs. [39] for a quantum field theoretical explanation and Refs. [62,72] for discussions in the context of atomic ionization. Adding up the individual contributions then automatically results in quantum interferences.

FIG. 3.
Particle distribution functions f (px) (top) and f (pz) (bottom) normalized by the expected production rate in local constant field approximation. Due to the fact, that the field is polarized along ex quantum interferences show up in parallel direction px only. The smaller λ the higher the magnetic peak field strength and, in turn, the stronger particles are accelerated along pz-direction. As for λ = 3m −1 particle bunches do not occupy the same region in phase-space any more, the interference pattern vanishes. Both simulations feature a temporal pulse length of τ = 25m −1 and a frequency of ω = 0.2m.
As λ decreases the magnetic field strength rises and thus the overall impact of the magnetic field increases. Pictorially speaking, the magnetic field between the three main peaks in E(t, z) acts as an accelerator in direction of p z . However, as B(t, z) is oscillating in space and time, particles created at subsequent peaks are accelerated in opposite directions. As a result, those particles do not share the same phase space at final times any more. Averaging over particle paths with similar phase information, however, does not result in quantum interference ultimately leading to a smooth distribution function, c.f. the grey dotted curve in Fig. 3.
The distribution function f (p z ) can be very well understood in the same way assuming that (i) particle creation is exponentially suppressed with exp − m 2 + p 2 z eε , see Eq. (7), and (ii) it takes a strong magnetic field to deflect the particles significantly. In case of a vanishing magnetic field the transversal particle spectrum f (p z ) shows a Gaussian distribution, see Fig. 3. The absence of quantum interferences is given due to the fact that although all particles have picked up a phase, the distribution of the phase information only varies in direction of p x . Hence, when summing over the phases in direction of p z no interference pattern appears.
For strong magnetic fields the peak at p z = 0 splits into two weakly pronounced peaks at p z ≈ ±0.6m (eε = 0.2m 2 ) and p z ≈ ±1.25m (eε = 0.5m 2 ). These two peaks are not equal in height clearly favoring the ones at p z > 0, respectively. The reason for this asymmetry lies in our choice of representation as this has intrinsically fixed the particle spin, see Sec. III. Hence, despite the fact that A (t, z) is symmetric in z performing all calculations for only one 2-spinor basis and, therefore, neglecting half of all electrons and positrons automatically results in an uneven particle distribution in phase-space. As a consequence the spatially varying magnetic field introduces a net force in one transversal momentum direction (8)-(10). In summary, smooth particle distributions superposed by quantum interferences are clear signatures of the Schwinger regime. Furthermore, the momentum spectrum can be understood under the assumptions, that the effective field strength E(t, z) 2 − B(t, z) 2 determines the chances for particle creation and particles, once created, follow semi-classical trajectories.

B. Absorption-dominated pair production
In a multi-cycle field the virtual pair can obtain energy from the background field via absorbing photons. If the total energy gain is higher than the production threshold an electron-positron pair is created. The decisive quantities are the number and energy of the absorbed photons nω as well as the particle's rest energy plus a modification factor due to the oscillatory motion of the particle, see Ref. [73]. Furthermore, if an l-photon process could occur there is also the chance for l + sphoton absorption. In such a case, particles are created with higher kinetic momenta and a different angular momentum profile [74]. As a result, the final particle distribution is given by multiple ellipses of fixed energy E n , where n gives the number of photons initially involved in the process.
Employing a background field with a field strength of eε = 0.5m 2 and a field frequency of ω = 0.5m we find many characteristic traits of photon absorption in the particle distribution function, see Fig. 4. These structures are thus best discussed in terms of energy conservation. The only energy source is given by the absorption of photons nω. The total absorption energy must be equal the total particle energy E n , which is given by the particles' rest energy plus their ponderomotive energy. The former can be stated easily as E = m 2 + p 2 x + p 2 z , while the latter depends on the fields' strength, frequency and polarization direction, c.f. Ref. [63] for a complete derivation. Here, the oscillations in the electric field only lead to a modification of the particles' parallel momentum component As a result the various absorption channels form ellipses in the particle phase-space instead of circles as one would naively expect.
For the sake of a better understanding of the interference pattern we display the particle distribution along the ellipses in Fig. 5. Note, that we have used linear interpolation techniques to illustrate the distribution function f (p x , p z ) as a function of the polar angle ϑ, where ϑ is defined as the ejection angle with respect to the fields' polarization direction. We immediately see, that lines with an odd number of photons are strongest at ϑ = 0 and ϑ = π (p z = 0) and show side maxima at ϑ = ±π/2 (p x = 0). In turn, an even photon number corresponds to minima at ϑ = −π/2 and ϑ = +π/2 (dotted orange line).
This behavior can be very well understood considering that the toy model given in Eq. (2) still describes fields within the dipole approximation. As parity as well as charge parity has to be conserved, the dipole approximation allows us to discuss the particles' angular distribution by simply counting the number of possible FIG. 4. Non-linearly scaled particle momentum spectrum f (px, pz) for large (top) and small (bottom) spatial extent λ. Electric field strength eε = 0.5m 2 , pulse length τ = 25m −1 and field frequency ω = 0.5m are fixed. Due to the presence of a strong field and the high photon energies particles are predominantly created along ellipses. For a vanishing magnetic field (top) multiphoton peaks as well as a pronounced interference pattern are clearly visible. A strong magnetic field disturbs these patterns and additionally applies a strong force in perpendicular direction ±pz.
final quantum states. To be more specific, the intrinsic parity of an electron-positron pair is (−1). The particles' orbital momentum contributes by an additional factor (−1) L . Charge conjugation symmetry gives (−1) L+S depending on the particles' spin orientation (S = 0 or S = 1, cf. Ref. [56]). For the incoming photons we have C-parity (−1) n and parity (−1) due to the fact that only electric dipole transitions are possible. Nevertheless, upon absorption every photon transports a unit of angular momentum to the pair ∆L = ±1. The change in the magnetic quantum number is zero though, because Eq.
An n-photon process therefor requires a parity of (+1) for an even number of photons and (−1) for an odd number. Consequently, if n is even L has to be odd. Performing a partial wave analysis [56] we find that the FIG. 5. Normalized radial distribution function f (ϑ) for various n-photon channels for a configuration with field frequency ω = 0.5m, field strength eε = 0.5m 2 and pulse length τ = 25m −1 . The higher the photon number n the more minima/maxima in the spectrum are visible. In particular at ϑ = ±π/2 we have alternating extrema. The variable ϑ rotates clockwise with starting point (px,n, 0). final particle state can be conveniently written in terms of Legendre polynomials In this case, L is odd thus the sum in Eq. (60) is over the Polynomials with L odd only. At vanishing parallel momentum we have ϑ = ±π/2 for which all these remaining Legendre polynomials vanish. As a result, in a pure multiphoton absorption process with an even number of photons, neither electrons nor positrons can be emitted in a 90 degree angle. Besides, these structures in the spectrum are very sensitive to an external magnetic field. With increasing magnetic field strength the particles are accelerated in transversal direction similarly to the results displayed in section V A. The most notable difference is the intact symmetry in p z even for extremely strong magnetic fields. The reason is that in this case the particle formation time is much longer, thus it is impossible to attribute a peak in E(t, z) with a peak in f (p x , p z ). Consequently, particle trajectories are not unevenly separated and as a result the symmetry in p z is not broken.

C. Multi-mechanism pair production
In the following, we discuss field configurations that combine multi-cycle pulses with high field strengths to enhance Schwinger pair production via absorption effects and vice versa. In order to do so we employ a super-Gaussian envelope function to ensure that the background operates close to its maximum value for multiple field cycles. The consequences are twofold. At the one hand, chances for Schwinger pair production increase, because also the "side peaks" can produce a sizable amount of particles. On the other hand, a higher amount of significant field oscillations also increases the likelihood of n-photon absorption processes.

Quasi-homogeneous fields
As one can see in Fig. 6 the calculated particle spectrum indeed displays a mixture of exponentially decaying as well as elliptical structures. At vanishing particle momentum we obtain an interference pattern typical for tunnelingenhanced pair production. To be more specific, the field configuration under consideration (super-Gaussian envelope with eε = 0.2m 2 and ω = 0.2m) is at the edge of seeing the 12-particle channel directly. Nevertheless, the absorption of 12 photons creates a highly excited state. As this is the equivalent of lowering the threshold by the same amount subsequent tunneling is barely suppressed. In this way the onset of a new channel can be seen very well even below the threshold. On a side note, the particle distribution falls off exponentially in transversal direction, which is a clear sign for tunneling.
We analyze this area using the simple man's model. As the model's output is a smooth particle density peaked at vanishing momenta we can easily obtain the equipotential lines, where the distribution function holds 10% of its maximal value. In case of eε = 0.2m 2 we obtain the shape parameters a = 0.375 and b = 0.3. This is extremely close to the values obtained from the DHW calculation: p * x = 0.34 and p * z = 0.27. Similar holds for the widespread particle distribution in p x for strong fields, see Fig. 7. In terms of the simple man's model particles are produced around the main peaks of the electric field t k = kπ ω with k = −2, −1, . . . 2. In a quasi-homogeneous setup, a particle created exactly at one of these peaks acquires nearly no net momentum due to the almost flat envelope. However, the high field strength of eε = 0.5m 2 allows for easy particle production at the slopes, too. For example, att = 2.25m −1 the electric field shows a local field strength of E(t, 0) = 0.4E cr . As the vector potential is highly nonzero att, particles created at this instant in time are strongly accelerated up to a final momentum of p x,f ≈ 1.09m. As a result, a sizable amount of particles can be found even at large parallel momenta, e.g. f (p x ≈ 1.1m) ≈ 1 4 f (0). In addition, elliptical structures appear at higher particle momenta in Figs. 6-7, which can be interpreted as 13+ and 20+ photon absorption processes, respectively. The corresponding analysis using fits to determine the corresponding effective energy of the elliptic particle distributions in comparison with the predictions from the effective energy model is given in Tabs. I-II. The most surprising result is that simulation and model deviate slightly in assessing the particles final transversal momentum (∆p z ∼ 5%). Predictions for the parallel FIG. 6. Particle distribution f (px, pz) for a background field featuring a field strength of eε = 0.2m 2 , a super-Gaussian envelope with pulse length τ = 75m −1 , a field frequency of ω = 0.2m and a spatial extent of λ = 1000m −1 (top) as well as λ = 20m −1 (bottom). The spectrum shows a Schwinger-like distribution (within white dashed area) as well as multiphotonlike patterns (ellipses). momenta perfectly agree with the simulation for all abovethreshold peaks, while for each of the lowest ellipses the model underestimates the final value by up to ∼ 9%. It is very likely that for such low energies remnants of tunneling influence the final particle momentum.
Similarly to the previously presented case for highly oscillating fields, we display the multiphoton structures obtained for the strong field configuration in terms of radial distribution functions, see Fig. 8. Due to the large photon count n the maximal orbital angular momentum particles can acquire is much higher than in the previous case. As a result, these functions show a large number of side peaks. Moreover, multiphoton and tunneling distributions interfere making an evaluation based on conservation laws difficult. The wild pattern around ϑ = π are a remnant of this superposition. Nevertheless, the considerations on parity and C-parity still hold FIG. 7. Particle distribution f (px, pz) for a background field featuring a field strength of eε = 0.5m 2 , a super-Gaussian envelope with pulse length τ = 75m −1 , a field frequency of ω = 0.2m and a spatial extent of λ = 1000m −1 (top) as well as λ = 20m −1 (bottom). The main structure (cigar shaped area) is superposed by multiple ellipses stemming from multiphoton pair production giving rise to additional quantum interferences.

Strong magnetic fields
In multi-cycle fields the particles' phase-space occupancy is much more involved, thus quantum interferences form easily. Nevertheless, when exposed to strong forces due to the magnetic field, particle bunches are accelerated in direction of p z . However, they are boosted in such a way that their relative quantum phases hardly change. In fact, all the individual peaks in the particle distribution function for quasi-homogeneous fields can still be linked to the peaks observable in the spectrum for strong magnetic fields, c.f. Fig. 7. What changes are the positions of these peaks in the spectrum as well as their relative size.
FIG. 8. Normalized radial distribution function f (ϑ) for various above-threshold signals in momentum space for a background field with peak field strength eε = 0.5m 2 , pulse length τ = 75m −1 (super-Gaussian envelope) and frequency ω = 0.2m. Due to the high intensities the threshold for pair production is increased thus particles acquire less kinetic energy. Here, the photon numbers n correspond to effective energies of E20 = 0.65m, E21 = 1.06m, E22 = 1.36m and E23 = 1.60m.
The distribution function for parallel momenta f (p x ), see the blue solid line in Fig. 9, indicates that the Schwinger effect is the main source of particle production. In section V A we have already established, that a smooth distribution function superposed by quantum interferences are typical signs of tunneling pair production. The difference here is that we observe an irregular pattern on top of a broad spectrum. We interpret the data such that particles mostly tunnel through the Coulomb barrier, thus also the smooth decrease in f (p z ) for small momenta p z , c.f. Fig. 10. The irregular oscillations that superpose the smooth tunneling spectrum, c.f. Fig. 9, are caused by multiphoton processes which includes assisted tunneling. These patterns are still visible if the spatial extent is chosen to be small λ = 20m −1 . Only at extreme values, λ = 5m −1 , where pair production in general starts to break down, see Tab. III, the interferences fade away.  Table showing the effective energy of the particles created through the absorption mechanism. In order to compare the outcome of the simulation with the effective energy model we fit ellipses to the particle spectrum. This gives us the shape parameters a and b. The parameters p * x and p * z determine the ellipses of equal effective energy E (p * x , p * z ). Background field: eε = 0.2m 2 , τ = 75m −1 (super-Gaussian envelope), ω = 0.2m. We complete this section by comparing the total particle yield N for the configuration eε = 0.5m 2 , τ = 75m −1 (super-Gaussian envelope), ω = 0.2m for different values of λ. Similar investigations have already been performed for short, high-intensity fields [75]. However, as the particle spectrum in this setup shows clear signatures of Schwinger as well as multiphoton pair production, c.f. Fig. 7, it is an ideal candidate for a study of the impact of magnetic fields on the general creation rate. Naively, one would expect a linear dependence on the spatial extent λ, given the particle yield scales linearly with the volume. This assumption is indeed true for wide fields λ ≥ 10m −1 .
For strongly focused fields, one might expect a fasterthan-linear decrease due to the fact that the regions of significantly high effective field strength E(t, z) 2 −B(t, z) 2 shrink non-linearly. However, studying the normalized yield N/N cl in Tab. III we find that starting with λ ≈ 5m −1 the particle yield decreases much faster than TABLE II. Comparison of the effective energy model with the outcome of the simulation for particles created through the absorption mechanism. We fit ellipses to the particle spectrum in order to obtain the shape parameters a and b. The parameters p * x and p * z determine the ellipses of equal effective energy E (p * x , p * z ). Background field: eε = 0.5m 2 , τ = 75m −1 (super-Gaussian envelope), ω = 0.2m. FIG. 10. Normalized particle distribution function f (pz) for multiple values of λ (smaller λ corresponds to higher magnetic field strength). The change in the distributions with increased magnetic field strength can be understood as a two-step process. For weak magnetic fields particles that have already been accelerated by the electric field are deflected by the magnetic field in direction of pz. The stronger the magnetic field becomes the stronger the effect. At a certain limit (λ ∼ 10m −1 ) the particle rate drops considerably and a strong magnetic field prevents a clear signal to form. Additionally, for high magnetic field strengths the symmetry in pz is broken. Field parameters: eε = 0.5m 2 , τ = 75m −1 (super-Gaussian envelope), ω = 0.2m.
TABLE III. Table of total particle number N as well as the normalized particle number N = N/N cl as a function of the spatial extent λ. All other field parameters have been kept fixed at eε = 0.5m 2 , τ = 75m −1 (super-Gaussian envelope) and ω = 0.2m. For comparison, for every value of λ the deviation to an assumed linear dependency is shown ∆N lin . expected. This might be a hint towards a critical point in a time-dependent, spatially inhomogeneous, high-intensity field [76]. Such an investigation, however, is beyond the scope of this article, and will be addressed elsewhere.

VI. CONCLUSION & OUTLOOK
This study on pair production in inhomogeneous electromagnetic fields mark a significant step forward in demonstrating the capabilities of the Wigner formalism in general. Especially the possibility to perform calculations for almost arbitrarily focused backgrounds is a clear sign that quantum kinetic approaches are a valuable asset towards understanding quantum field theories in general.
To be more specific, we have adopted the Dirac-Heisenberg-Wigner formalism for large-scale computations by taking into account novel numerical methods. In this way, it was possible to calculate the particle creation rates as well as their momentum spectra in inhomogeneous electromagnetic fields. In this regard, we have significantly expanded the capabilities of quantum kinetic approaches such that we could finally study the impact of a strongly inhomogeneous, time-dependent magnetic field on pair production processes even in long-pulsed fields.
By thoroughly analyzing the so acquired data we were able to identify signatures of Schwinger as well as multiphoton effects in an intermediate regime, where no creation mechanism is dominant.
We further demonstrated how much impact the temporal envelope has on the final particle distributions as we could easily enhance and suppress certain signatures by simply switching from a Gaussian to a flat-top envelope.
In the course of this study we also observed symmetrybreaking due to spin-field interactions in tunnelingdominated regions as well as vanishing above-threshold peaks in absorption-dominated areas. For strong, multicycle field we were able to show that quantum interference patterns are preserved even in the presence of strong magnetic fields. Only if the spatial extent of the electric field is of the order of the Compton wavelength of the pair, these interferences vanish as the pair production process breaks down independent of the regime.

VII. APPENDIX
Collection of figures not suited for publication in the main body of the manuscript. Nevertheless, they carry interesting information, thus we decided to present all data available.