Trident pair production in plane waves: Coherence, exchange, and spacetime inhomogeneity

We study the trident process in inhomogeneous plane wave background fields. We obtain compact analytical expressions for all terms in the probability, including the exchange part, for an arbitrarily shaped plane wave. We evaluate the probability numerically using complex deformation of lightfront time integrals and derive various analytical approximations. Our results provide insights into the importance of the one-step and exchange parts of the probability relative to the two-step process, and into the convergence to the locally constant field approximation.


I. INTRODUCTION
The trident process in electromagnetic background fields, e − → 2e − + e + , is a basic process in electron-laser collisions. It was first studied in detail in the 70's for constant crossed fields in [1,2] before the advent of highintensity lasers (see [3] for a review of the development of such lasers). Trident, or at least part of it, was explored in a famous experiment at SLAC [4] in the late 90's, with a laser of still modest intensity. At that time, complete theoretical predictions for trident were still lacking. Since then, the available laser intensities have increased by a couple of orders of magnitude, but on the theory side the situation has been improved thanks to only a few publications [5][6][7]. The importance of trident in high-intensity laser processes therefore motivates further investigations, see also [8].
One part of trident is a two-step process where the initial electron emits a real photon that subsequently decays into an electron-positron pair, i.e. nonlinear Compton scattering followed by nonlinear Breit-Wheeler pair production, and its contribution is given by incoherently gluing together the two corresponding rates, see [5,9]. For a recent comprehensive analysis of the two-step process see [10]. One of the main questions for trident is under what conditions this two-step process is a good approximation for the total trident process. In the experiment at SLAC [4], corrections to the two-step part were estimated using the Weizsäcker-Williams approximation and in that parameter regime were found to be negligible compared to the two-step process 1 . Further, the two-step process is expected to dominate for sufficiently large a 0 = eE/(mω), where E is the field strength and ω the frequency of the field. This is important for particle-in-cell codes, see [11] for a review, where higher-order processes, e.g. cascades, are described as a sequence of first-order processes. We study here in more detail when corrections to the two-step process can be important.
The trident amplitude has two terms due to the exchange between the identical particles in the final state, and the absolute square of the amplitude gives a cross term which is referred to as the exchange part of the probability. In [1, 2] the direct or non-exchange part of the probability was obtained from the imaginary part of a loop diagram, i.e. using the optical theorem. The direct part is expected to dominate, e.g. for χ = a 0 b 0 1, where b 0 = kp/m 2 is the product of a characteristic wave vector of the laser, k µ , and the initial momentum of the incoming electron, p µ . Much less is known about the exchange part. Here we will calculate both the direct and the exchange parts in order to investigate in more detail when the latter can be important.
As in [5,6], we obtain the probability by calculating the amplitude using Volkov solutions for plane-wave backgrounds. In contrast to previous analytical studies [1, 2, 5], we derive compact expressions for nonconstant plane waves. While we recover the results in [1, 2,5] in the locally constant field limit, i.e. for a 0 1, our results also allow us to see how this limit is approached and to study regimes with a 0 ∼ 1. In addition, our analytical results offer a useful alternative starting point for numerical investigations, compared to ones based on e.g. Monte Carlo integration [7,8], as well as insights into the analytical structure of the trident probability for inhomogeneous fields.
Another difference compared to previous studies is that we use the lightfront formalism [12,13], which is particularly convenient when dealing with plane-wave background fields [14][15][16]. The lightfront Hamiltonian has both a "conventional" term, given by jA, and an "instantaneous" term. These two terms suggest another split of the amplitude and the probability based on the lightfront Hamiltonian. As we will show, this lightfront separation is not the same as the standard one-step/two-step separation as in [1, 2,5] (though the total probability and rate are obviously the same for both separations). While we do not propose this lightfront split as something to replace the standard separation, we will show that it is convenient from an analytical and computational point of view.
The rest of this paper is organized as follows. In Sec. II we describe the basic ingredients needed to calculate the trident probability with the lightfront formalism, and compare this approach with the standard one. In Sec. III we present compact analytical expressions for the trident probability in a general inhomogeneous plane wave. In Sec. IV we compare these expressions with the probabilities for nonlinear Compton scattering and Breit-Wheeler pair production, and in Sec. (V) we consider the locally constant field (LCF) approximation and show how the standard one-and two-step terms can be obtained from our lightfront expressions. In particular, we expand in 1/a 0 1 and recover literature results for the direct part. We also calculate the exchange terms in this limit, which are new results. In Sec. VI we go on to consider χ 1 and a 0 1 but for non-constant background fields. Then in Sec. VII and Sec. VIII we consider χ 1 but a 0 ∼ 1 and obtain simple analytical approximations for a pulse and a monochromatic field. In Sec. (X) we explain how to numerically integrate our results from Sec. III. In Sec. IX we consider a 0 1 for up to moderately large χ and study the importance of exchange terms. In Sec. (XI) we consider a pulsed field and study the convergence to the LCF approximation. We conclude in Sec. (XII). In Appendix A we recover known results for the perturbative limit, and in Appendix C we show the agreement between our analytical approximations and numerical results.
We use units with c = 1 and = 1, and measure energies in terms of the electron mass such that m = 1.

II. FORMALISM
We use lightfront coordinates defined for an arbitrary vector v µ by v ± = 2v ∓ = v 0 ± v 3 and v ⊥ = {v 1 , v 2 }, for coordinatesx = {x − , x ⊥ } and for momentap = {p − , p ⊥ }. We consider an arbitrary pulsed plane-wave background field given by f µν = k µ a ν − k ν a µ , where k µ = k + δ + µ is a null wave vector and a ⊥ (kx) is an arbitrary function. We absorb for convenience a factor of the electron charge into the definition of the background field, i.e. ea µ → a µ . Several of our results are conveniently expressed in terms of the Lorentz momentum of an electron in a plane-wave background, which is given by The initial state contains an electron with momentum p µ and spin σ, and the final state contains two electrons with p µ 1,2 and σ 1,2 and a positron with p µ 3 and σ 3 . The trident amplitude, M , is obtained from where U is the evolution operator andδ(. . . ) = (2π) 3 δ −,⊥ (. . . ). The mode operators are normalized according to {b(q, r),b(q , r )} = {d(q, r),d(q , r )} = 2p −δ (q − q )δ rr , and we use dp = θ(p − )dp − d 2 p ⊥ /(2p − (2π) 3 ) to denote the Lorentz-invariant momentum measure. The initial electron is described by a sharply peaked wave packet f (p), where the last equation ensures that in|in = 1. The total probability averaged over the initial spin is given by all spin dp 1 dp 2 dp 3 dp f 1 where one factor of 1/2 comes from averaging over the spin of the initial electron and another 1/2 is due to having identical particles in the final state, andp 3 =p −p 1 −p 2 . The total amplitude can be written M = M 12 − M 21 , where M 21 is obtained from M 12 by swapping the identical particles, i.e. by replacing p 1 ↔ p 2 and σ 1 ↔ σ 2 . This leads to a separation of P into a direct and an exchange part, |M | 2 = |M 12 | 2 + |M 21 | 2 − 2ReM 21 M 12 , where the first two terms give the direct contribution and the third term gives the exchange contribution, i.e.
where the second term in P dir is obtained from the first by replacing p 1 ↔ p 2 and σ 1 ↔ σ 2 .

A. Lightfront quantization
We have derived our results using two different approaches. In both approaches the plane-wave background field is taken into account exactly using Volkov solutions in the Furry picture. In the first approach we use the combination of the Hamiltonian-based lightfront formalism [12,13] and plane-wave backgrounds [14][15][16]. The evolution in lightfront time, x + , is determined by the lightfront Hamiltonian, and the interaction part of this Hamiltonian has three terms, where the photon field is given by the current is j µ =Ψγ µ Ψ and the spinor field is expressed in terms of Volkov solutions as where ϕ − = ϕ(−p) andK = 1 − / k/ a/(2kp). Note that all the momenta in these mode expansions for A µ and Ψ are on-shell. In particular, in this formalism all photons are on-shell and so one cannot split trident into two parts with one having on-shell and the other off-shell intermediate photons. Only the first two terms in (6), H (1) int and H (2) int , contribute to trident to lowest order. The first term is familiar from ordinary quantization and its contribution to the amplitude is given by a double x + -integral with x + -ordering, The second interaction term, sometimes referred to as instantaneous [12,17], only involves one x + -integral and is given by After some calculation we find where φ = kx,l =p −p 1 and l + = l 2 ⊥ /4l − . To avoid clutter we put the arguments below the function, e.g. uKφ p1 =ū(p 1 )K(p 1 )φ(p 1 ).

B. Relation to standard approach
There are two main differences between the approach just described and the standard approach: The former is a Hamiltonian formalism in the ligthfront gauge. These two approaches should of course give the same results, and we will show this explicitly in this section. In the non-Hamiltonian approach, the amplitude is given by (up to an irrelevant overall phase) where ψ = Kuϕ, ψ − =Kvϕ − , and D νµ is the photon propagator In the Feynman gauge D νµ = g νµ and in lightfront gauge D νµ = L νµ (see [18] for a recent discussion of the photon propagator in the lightfront gauge). Performing the trivial integrals in (12) gives delta functions implyingl =p −p 1 . To see that the two gauges give the same result, consider the contributions from the l µ -terms in L νµ , which involvē ψ(p 1 )/ l ψ(p) andψ(p 2 )/ l ψ − (p 3 ). Consider first the part with p and p 1 . By writing l µ =: π µ (p, φ) − π µ (p 1 , φ) + c(φ)k µ and using K / p = / πK, we findūK(p 1 )/ l Ku(p) =ū(p 1 )/ ku(p) c(φ). Since the φ-dependence of the exponential is given by the integral of c(φ), we hence find a total derivative that vanishes upon integrating over φ. The same holds for the part with p 2 and p 3 . Thus, D νµ = g νµ and D νµ = L νµ give the same result. The next step is to reproduce the amplitude obtained in the previous section. By separating L µν into an on-shell and off-shell part, one finds (c.f. [18]) When performing the l + -integral in the propagator, the first term in (14) gives a lightfront time-ordering step function θ(φ y − φ x ), which one also finds in the Feynman gauge, see Appendix B, while the second term in (14) gives an "instantaneous" δ(φ y − φ x ), which does not appear in the Feynman gauge. (See [19] for a similar separation of the fermion propagator.) We hence find  )) and the two diagrams on the right represent the two lightfront terms in (15). Note that these diagrams represent terms in the amplitude after performing the photon momentum integrals in the photon propagator, so the photon-momentum, lµ, in these diagrams is on-shell. This is why the third diagram appears. The line through the photon line in the third diagram stands for an "instantaneous photon" [12].
where the photon momentum is now on-shell, i.e. l + = l 2 ⊥ /(4l − ). Because of (14), (15) contains two terms that are equivalent to one term in the Feynman gauge (compare (15) with (B1)), as illustrated in Fig. 1. In (15) we have the same M 1 and M 2 as in (11) (up to an irrelevant overall phase), which we obtained using the lightfront Hamiltonian formalism. In particular, the step-function term in (15) int . While (B1) and (15) give the same amplitude, the lightfront separation of M 12 as in (15) can be convenient because of the vector structure of the two individual components, M 12 1 and M 12 2 ; we have for example L µν k ν = L µν l ν = 0 and / kK = / k. Thus, in addition to the direct/exchange separation of the amplitude, M = M 12 − M 21 , we now separate each of these into a "three-point vertex" and a "lightfront instantaneous" part, M 12 = M 12 1 + M 12 2 . This leads to six different contributions to the probability, The sum of these six terms gives the total probability. Given that M 1 comes from a term in the lightfront Hamiltonian that is usually referred to as instantaneous, it might be tempting to associate the corresponding terms in the probability with what one would usually call one-step terms. However, as we will show, all the terms listed in (16) and (17) contribute to the standard one-step term, and in fact "most" of the standard one-step term comes from P 22 (which also contains the entire standard two-step term).

III. EXACT ANALYTICAL RESULTS
Rather than to keep the dependence on all the particle parameters, we content ourselves with the dependence on the longitudinal momenta, which we denote as s i = kp i /kp. The integrals over the transverse momentum components of the final electrons, p 1,⊥ and p 2,⊥ , are Gaussian. We perform these integrals as well as the traces, coming from the spin summations, analytically for arbitrary pulsed plane-wave backgrounds. The resulting expressions contain integrals over the longitudinal momenta plus integrals over lightfront time. We use the following notation for the probability density where s 3 = 1 − s 1 − s 2 . To make the symmetries in the following expressions manifest, we introduce s 0 = 1 in the appropriate places. For the terms coming from the square of the "instantaneous" part of the amplitude, |M 1 | 2 , we find where b 0 = kp, r 1 := 1 s1 − 1 s0 , r 2 := 1 s2 + 1 s3 , q i := 1 − s i , dφ 12 := dφ 1 dφ 2 , θ ij := φ i − φ j , Θ ij := θ ij M 2 ij , and M is an "effective mass" defined via the lightfront time average of the field as [20] where π µ is given by (1). These shorthand notations are very convenient for the more complicated terms below. To make the transverse momentum integrals converge, we have introduced an infinitesimal convergence factor > 0 via φ 2 → φ 2 + i /2 and φ 1 → φ 1 − i /2, which after performing the momentum integrals gives an i -prescription for how to integrate around the singularity at θ 21 = 0.
Next we consider the cross terms between the "lightfront instantaneous" and the "three-point vertex" parts of the amplitude. To make the momentum integrals for P 12 converge, we similarly take φ 2 → φ 2 +i /2 and φ 1,3 → φ 1,3 −i /2. We find where D 12 = ∆ 12 ·∆ 32 and It turns out that the background field only enters the pre-exponential factors and the exponentials via ∆ ij and M 2 ij , respectively; no other field combinations are needed.
For the direct part of the square of the "three-point vertex" amplitude we find where κ ij = (s i /s j ) + (s j /s i ), D 1 = ∆ 12 ·∆ 21 , D 2 = ∆ 34 ·∆ 43 , and where w 1 = ∆ 12 , w 2 = ∆ 21 , w 3 = ∆ 34 and w 4 = ∆ 43 . Note that W ijkl = 0 for linear polarization. These terms have some symmetries that can be understood in terms of the probability diagrams in Fig. 2. Let P 22 dir (s) = P 22 dir (s 1 ; s 2 ) + P 22 dir (s 2 ; s 1 ), where P 22 dir (s 1 ; s 2 ) is everything before the (s 1 ↔ s 2 )-term in (24). P 22 dir (s 1 ; s 2 ) is invariant under {s 0 ↔ −s 1 } as well as {s 2 ↔ s 3 }, which correspond to reflection of one or the other of the two fermion loops in the P 22 dir diagram in Fig. 2. The integrand of P 22 dir (s 1 ; s 2 ), apart from θ(θ 31 )θ(θ 42 ), is also invariant under a total reflection and a 180-degree rotation, i.e. under {φ 1 ↔ φ 3 , φ 2 ↔ φ 4 , s 0 ↔ −s 3 , s 1 ↔ s 2 , q 1 → −q 1 } and , which is everything before the (s 1 ↔ s 2 )-term in (21), we find similar symmetries, except that P 12 dir (s 1 ; s 2 ) changes sign under {s 0 ↔ −s 1 } as well as {s 2 ↔ s 3 }, which implies that after integrating over the longitudinal momenta we find P 12 dir = 0. Under a reflection of the P 12 diagram in Fig. 2 . The corresponding exchange term, P 12 exch (s 1 ; s 2 ), is invariant under the same reflection with the same change of the step function. P 12 exch (s 1 ; s 2 ) is also invariant under {s 0 ↔ −s 1 , s 2 ↔ s 3 }, but not under {s 0 ↔ −s 1 } and {s 2 ↔ s 3 } separately. The direct and exchange terms of P 11 have similar symmetries. Only P 12 dir = 0, so for the total/integrated probability we have five nonzero terms.
These types of symmetry considerations are particularly useful for the exchange part of the square of the "three-point vertex" amplitude, which is the most complicated term. After some lengthy calculation we eventually find Now, the first term in the prefator of (26) is quartic in the field, where the second term in the square brackets vanishes for linear polarization (c.f. (25)) and can also be written as The next two terms in (26) are quadratic in the field, and the last two terms in (26) do not contain the field, The permutation in (28) is a symmetry of the integrand in the following sense: We have P{F 0 , f 0 , f 1 , z 1 , z 2 } = {F 0 , f 0 , f 1 , z 1 , z 2 } and for the exponential in (26) we have P[exp(i...)] = exp(−i...). (Note that the exponent in (26) reduces to that in (21) and (22) for φ 4 = φ 2 .) Let I be the integrand in (26) without the step functions, then P[I] = −I * , but the step functions are different for each of the 4 permutations. We also have a reflection symmetry: . We have R[I] = I * and R leaves the step functions unchanged, so the integral (26) is invariant under R.

IV. TWO-STEP AND ONE-STEP
In this section we will show how the above results are related to the two individual processes of nonlinear Compton scattering followed by nonlinear Breit-Wheeler pair production. We first express the corresponding probabilities in a similar way as the expressions for trident above. We assume here a linearly polarized field and an emitted photon with momentum l µ and polarization vector ⊥ = (cos ϑ, sin ϑ), + = ⊥ l ⊥ /(2l − ), − = 0 (so l = k = 0). Averaging and summing over the spins of the incoming and outgoing electrons, respectively, gives us where the i -prescription, θ 21 + i , is left implicit. For pair production by the emitted photon we find after summing over the spins of the electron-positron pair where, again, the singularity is avoided with θ 43 + i . While (34) and (35) are already suitable for comparing with our expressions above for trident, by performing partial integration for the terms proportional to 1/θ 2 we obtain simpler expressions, and By summing (36) over two orthogonal photon polarization vectors, we recover eq. (3) in [21]. The reason it might not be convenient to perform similar partial integration for e.g. (24), is that there are step functions in (24) that would then generate non-vanishing boundary terms. For a constant field, these two expressions reduce to eq. (22) and and (26) in [22] (up to some overall constants due to the difference between probabilities and rates, and ordinary and lightfront volume factors). Without worrying too much about the formation length for the non-constant fields we consider, we follow the procedure for constant fields [5] and glue together the probabilities for nonlinear Compton scattering and Breit-Wheeler pair production and compare the result to the trident probability. By adding the probabilities (34) and (35) for parallel and perpendicular polarization, i.e. P C P BW (ϑ = 0) + P C P BW (ϑ = π/2), and symmetrizing with respect to s 1 ↔ s 2 we find an expression that is identical to our P 22 dir except for the step functions in (24). In the product approach one would also include a step function associated with causality. We are therefore lead to inserting θ(σ 43 − σ 21 ), where σ ij = (φ i + φ j )/2. While this step function only restricts the averages of the φ-variables associated with photon emission and pair production, in the LCF regime we recover the results in [5], which we will show in the next section. So, to separate the one-step and the two-step parts of (24) we rewrite the step functions there as and separate (24) as P 22 and P 22→1 dir are obtained by replacing θ(θ 42 )θ(θ 31 ) in (24) with the first and second term in (38), respectively. Thus, we separate the total probability as P = P 2 + P 1 , where We will refer to R 1,2 as rates. All the six lightfront terms contribute to R 1 (s), though P 12 dir gives zero contribution after integrating over the longitudinal momenta. The two-step rate can be obtained from the probabilities of nonlinear Compton and Breit-Wheeler as The Because of the scaling of these rates, it is natural to consider instead which can be seen as the rates corresponding to the propertime τ , as in a plane wave it is simply related to lightfront time via φ = kx = kpτ = b 0 τ . While we have motivated the definition of P 2 with its relation to the product of P C and P BW , the separation P = P 2 + P 1 might still seem ambiguous (see also the discussion in [23]). However, we will provide further motivation for the separation (39) by making an expansion in 1/a 0 . In particular, the two leading orders agree with literature results for constant fields.

V. a0 1 AND THE LOCALLY CONSTANT FIELD APPROXIMATION
As mentioned, P 22 should not be confused with the standard two-step obtained by gluing together the rates of nonlinear Compton and Breit-Wheeler as described e.g. in [5]. Indeed, P 22 contributes to both P 2 and P 1 in (39) and we will show in this section that it is P 2 and P 1 that to leading order correspond to the standard two-step and one-step parts, respectively. To do so we consider a linearly polarized field a(φ) = a 0 f (φ), where the maximum of f and f are around unity, with a 0 1 and make an expansion in 1/a 0 . As is well known, this limit takes us to the LCF approximation, for which there are analytical results in the literature to compare with. However, since we do not treat the background as constant from the start, our approach also allows us to obtain corrections to the leading order and we can avoid (very large) volume factors. We expand the probability in 1/a 0 as where, as we will demonstrate below, P 2 and P 1 correspond to a two-step and a one-step term, respectively. Each term in this expansion is a nontrivial function of χ, which is treated as independent of a 0 in this expansion. Since a 0 = E/ω the expansion in a 0 is a derivative expansion. The higher orders are usually not considered in the literature, but can be obtained with this approach. In order to provide further motivation for the separation (39) it is important to note that P 2 = a 2 0 P 2 + O(a 0 0 ), so only P 2 contributes to a 2 0 P 2 and only P 1 contributes to a 0 P 1 . In other words, the next-to-leading-order correction to P 2 does not mix with the leading order of P 1 .

A. Two-step
We begin with the two-step part, P 2 = a 2 0 P 2 + O(a 0 0 ). We change variables to σ 43 , θ 43 , σ 21 and θ 21 . In order to expand in 1/a 0 we rescale θ 21 → θ 21 /a 0 and θ 43 → θ 43 /a 0 . To leading order, the entire θ integral can be performed and expressed in terms of Airy functions. We find (note that a 0 /χ = 1/b 0 )   , where χ(σ) = χ|f (σ)| is the local value of χ, and the Airy integral is defined as in [5], (43) is the natural generalization of e.g. eq. (16) and (17) in [5] to non-constant fields. In [5] the trident probability was calculated for a constant field, and the term quadratic in the volume ∆φ was shown to be equal to what one obtains by gluing together the rates for nonlinear Compton and Breit-Wheeler. In [5] it was also explained how this product approach is generalized to the LCF approximation. Our (43) is in fact equal to the LCF expression given in [5]. To see this, we rewrite (43) as where the first and second factors correspond to nonlinear Compton and Breit-Wheeler, respectively, and the sum is over the photon polarization; these two factors are equal to eq. (20) in [5] but evaluated at two different φ, as in eq. (33) in [5]. This agreement is of course not a surprise given that we have already shown that P 2 can be expressed as (40) for non-constant fields and without expanding in 1/a 0 . Note that b 2 0 (43) only depends on the field via χ(σ 21 ) and χ(σ 43 ), which is why it is natural in the LCF regime to consider the proper-time rates as defined in (41).
Here we are mostly interested in regions of parameter space where one has important contributions from one-step terms, to which we now turn.

B. One-step
Next we consider P 1 in (39). All the six lightfront terms in (16) and (17)  ) and all terms in (39) are O(a 0 ). Together these O(a 0 ) terms give what is usually referred to as the one-step part of the probability.
For P 22→1 dir we change variables according to which in terms of the original variables are given by As for the two-step term we rescale θ → θ/a 0 and η → η/a 0 . As the second term in (38) shows, this also forces ϕ to be small, so we also rescale ϕ → ϕ/a 0 . To leading order, i.e. O(a 0 ), the ϕ integrand is constant and we find The leading order of P 11 , P 12 and P 22 ex can be obtained in a similar way. In the previous subsection we showed that P 2 agrees to leading order with literature results for the two-step term for arbitrary values of χ. However, while there are exact analytical expression also for the one-step terms, see [1, 2,9], these are not easy to compare with, so we will instead consider χ 1 and show that our O(a 0 ) terms agree with the results in [1, 2,9] in this regime.

C. Numerical spectrum
Before turning to analytical approximations, though, let us first plot the one-step and two-step rates, R 1 and R 2 , as defined in (39). To leading order in the LCF approximation we have, at constant χ, R LCF 2 ∝ a 2 0 and R LCF 1 ∝ a 0 (recall that R 2 = R LCF 2 + O(a 0 0 ) so the next-to-leading-order correction to R 2 does not mix with the leading order in R 1 ). Even in the LCF limit, R LCF 2 still depends on two different field values, a (σ 21 ) = a (σ 43 ). It is convenient to plot since this only depends on a 0 , b 0 , σ 21 and σ 43 via Fig. 3 and 4 we have plotted R LCF 1 and R LCF 2 as triangular contour plots where each of the three s i values is given by the distance to one of the triangle's sides, for different values of χ (or rather χ(φ) and χ(σ 21 ) = χ(σ 43 )). As expected e.g. from [22], the spectrum is peaked at s 1 = s 2 = s 3 = 1/3 for low χ, and close to the corners, s 1 ≈ 1 or s 2 ≈ 1, for large χ. Note that these plots contain all the one-step terms, both the direct and the exchange ones. The values on top of each plot give the maximum/minimum value, i.e. min < R1 < max, so for χ = 1/2 we have −2.13 * 10 −10 < R1 < 0, and −9.8 * 10 −5 < R1 < 1.0 * 10 −3 for χ = 16. The difference between two neighboring contours is 5% of the larger of |max| and |min|, and purple and red correspond, respectively, to values close to min and max.

D. Constant fields and χ 1
As Fig. 3 and 4 show, for χ 1 the spectrum is peaked at s 1 = s 2 = s 3 = 1/3. We can therefore perform these integrals using the saddle-point method. The lightfront time integrals can also be performed by the saddle-point method as described in the next section for non-constant fields.
For a constant field and χ 1 we find We have already shown that the O(a 2 0 )-term in P 22→2 dir agrees with the literature for arbitrary χ. By comparing P 22→1 dir in (49) with eq. (29) in [2] or eq. (6.57) in [9] (see also eq. (28) in [5] for normalization similar to ours), we see that P 22→1 dir for χ 1 also agrees with previous results. So, P 22 dir already gives the known χ 1 results for both the two-step and the one-step terms, even though P 22 dir is only one out of a total of six terms. The reasons for this are: P 12 dir (s) vanishes upon integrating over s 1 and s 2 , P 11 dir turns out to be smaller than P 22→1 in (49) (see below), and the remaining terms, P 11 ex , P 12 ex and P 22 ex , give the exchange part of the probability, while the previous results just cited only give the direct part. Indeed, to the best of our knowledge, there are no analytical expressions for the exchange terms in this regime to compare with. So, our results for the exchange part below are new.
From (19), (21) and (22) we find where P 22→1 dir is given by the a 0 term in (49). Eq. (50) includes both the direct and the exchange part. In fact, their contribution is on the same order of magnitude: we have P 11 ex = −(1/2)P 11 dir (to leading order in χ), and P 12 dir is identically zero so P 12 = P 12 ex . However, both P 11 and P 12 are smaller than P 22→1 in (49) by a factor of χ 1. In deriving (49) we have already thrown away terms on the order of (50), so the terms in (50) are only part of the higher-order corrections, see Appendix C.
For the last term we find Importantly, P 22 ex is on the same order of magnitude as P 22→1 dir in (49), so to leading order in 1/a 0 1 and χ 1 the one-step part is given by P 1 = a 0 P 1 = P 22→1 dir + P 22 ex , i.e. the exchange and the direct part of P 1 are on the same order of magnitude. Of course, in the regime we consider here, both these contributions to P 1 are small corrections to the two-step P 2 , but at least among the one-step terms, the direct and exchange parts are on the same order of magnitude as long as χ is not too large. We have checked that the analytical approximations in this section agree with our numerical results for sufficiently small χ, see Appendix C. In the next section we will consider a 0 1 and χ 1 for non-constant fields.
VI. PULSED FIELDS WITH a0 1 AND χ 1 In this section we will generalize the results in the previous section and obtain simple analytical approximations for a 0 1 and χ 1 for non-constant fields. We consider a linearly polarized field a(φ) = a 0 f (φ) that has a single maximum at φ = 0, i.e. f (0) = 0, and expand around this point. All the integrals are finite and there are no volume factors. We fix the field strength and the frequency by normalizing f (0) = 1 and f (3) (0) = −ζ < 0. We could set ζ = 1, but this might not always be convenient. For example, for a Sauter pulse it is natural to choose f (φ) = tanh φ and then ζ = 2.
We begin with P 22 dir in (24) and change variables as in (46). We rescale θ → θ/a 0 and η → η/a 0 , for P 22→1 dir we also rescale ϕ → ϕ/a 0 , and then expand the integrand in powers of 1/a 0 . For χ 1 we can perform all the integrals with the saddle-point method. We have a saddle point at The second-order variation in the exponent is for the momentum integrals given by exp − 36 χ (δs 2 1 + δs 1 δs 2 + δs 2 2 ) , where δs i = s i − 1/3, so δs i ∼ √ χ 1. The other variations are also ∼ √ χ 1, so we rescale all the integration variables with √ χ and expand the integrand in a series in χ. The resulting integrals are all elementary, and we thus find where f (n) 0 = ∂ n φ f (0). We have assumed here that 1/(χa 2 0 ) 1 in order to expand a term in the exponent down to the prefactor. The terms proportional to a 2 0 and a 0 correspond, respectively, to what are usually referred to as two-step and one-step terms. The a 2 0 term can hence be obtained by expanding the Airy functions in (43), which in turn can be obtained directly from the literature. The reason that we have included higher orders in χ for the a 2 0 term is that these can be on the same order as the leading-order-in-χ contribution to the a 0 0 term, and similarly for the first term in (54). How many orders in χ one should keep at a given order in a 0 depends of course on the relative size of χ and a 0 ; we have included the higher-order terms in (53) and (54) mainly as examples of this double expansion. In any case, using Mathematica it is straightforward to calculate higher orders in both 1/a 0 and χ. Note that there are no volume factors here (cf. the ∆φ terms in the previous section) because of the damping given by ζ > 0.
For (19) and (22) we find These are both smaller by a factor of χ 1 than the leading-order term in P 22→1 dir , i.e. the first term in (54). This agrees with what we found above for constant fields, see (50).
Next we consider P 22 ex given by (26). In order to expand in 1/a 0 we rescale θ → θ/a 0 , η → η/a 0 and ϕ → ϕ/a 0 . We expand the integrand around the saddle point (52) and perform the resulting integrals. We find As in the constant-field case in the previous section, P 22 ex gives the dominant contribution to the exchange part and it is on the same order of magnitude as the direct one-step term (54). Thus, the one-step terms, a 0 P 1 = a 0 P dir 1 + a 0 P ex 1 , are to leading order in χ 1 given by a 0 P dir 1 = lin a0 P 22 dir and a 0 P ex 1 = lin a0 P 22 ex , where lin a0 refers to the terms that are linear in a 0 . We find, to leading order in χ, the same relation as in (51) between P dir 1 and P ex 1 , i.e. P ex 1 /P dir 1 = 13/18, and both P dir 1 and P ex 1 are negative. The exchange term is also negative for a 0 1, which we can confirm by comparing with the literature, see Appendix A.
VII. SAUTER PULSE WITH a0 ∼ 1 AND χ 1 In this section we will consider a Sauter pulse, a(φ) = a 0 tanh φ. We again assume χ 1, but now we do not assume that a 0 is large. The results in this section for a 0 ∼ 1 therefore go beyond the LCF approximation. Our starting point is the exact expressions in Sec. III. Recall that the field enters these expressions only via M 2 and ∆. For the Sauter pulse the integrals in M 2 and ∆ can be performed analytically and this gives For χ 1 we can perform all the integrals with the saddle-point method. We first change variables as in (46). We have a saddle point at By expanding the integrands around this saddle point we find and As in the previous section, P 11 and P 12 are smaller by a factor of χ 1, so to leading order in χ we have (these expressions can also be obtained without much extra work by starting instead with (B2) and (B3)) and In the limit a 0 1, (60) and (61) reduce to the expansions in (53), (54) and (56) (with ζ = 2 for this field). Eq. (61) shows that the exchange term is on the same order of magnitude as P 22→1 dir , and their ratio, P 22 ex /P 22→1 dir = 13/18, is the same as the O(a 0 ) terms we found in the previous section. In the limit a 0 1 these two terms gives the one-step term, i.e. P 22→1 dir + P 22 ex → a 0 P 1 , while P 22→2 dir → a 2 0 P 2 . However, for a 0 ∼ 1 the "one-step" terms, P 22→1 dir and P 22 ex , are on the same order of magnitude as the "two-step" term, P 22→2 dir . For a 0 = 1 we have P ex /P dir ≈ −0.4, while for a 0 4 the ratio decreases to −0.1 P ex /P dir < 0, see Fig. 5.
The saddle-point calculation of the pre-exponential factors above breaks down when a 0 becomes too small. However, the exponential part of (60) scales for a 0 1 as which is what one would expect for perturbative trident: In terms of the Fourier frequency,ω, of a(φ), the threshold for pair production is given byω th = 4/b 0 . Forω th 1 the Fourier transform of the Sauter pulse is approximately exponential, which gives the exponential in (64) (cf. the treatment of Sauter-like pulses in [24] for Schwinger pair production). Interestingly, the exponent in (60) has the same functional form as the one in photon-stimulated Schwinger pair production in a constant electric field, eq. (5) in [25], but with different parameters.
VIII. MONOCHROMATIC FIELD WITH a0 ∼ 1 AND χ 1 We now consider a monochromatic field, a(φ) = a 0 sin φ. Consider first P 22→2 . We find saddle points at where σ ij = (φ i + φ j )/2 and n 1 and n 2 are integers. The saddle points for σ correspond to maxima and minima of the field. The step function θ(σ 43 − σ 21 ) implies n 2 ≥ n 1 , so either photon emission and pair production happen at the same max/min or the photon is emitted at one max/min and travels to a different max/min where it decays into a pair. All n 1 and n 2 give the same contribution, except the term with n 1 = n 2 which is, due to θ(σ 43 − σ 21 ), a factor of 1/2 smaller than the other terms. With 2N maxima and minima, adding the contribution from all saddle points gives P 22→2 22→2 is the contribution from a single maximum or minimum, which we find to be The argument of the exponential in (67) has the same functional dependence on a 0 as the exponent in [26] (see the equation before eq. 15 in [26]) for the Breit-Wheeler probability. The only difference is an overall factor of 2kl/kp, where l µ is the photon momentum in the Breit-Wheeler case. This factor is expected in the locally constant field limit, i.e. for a 0 1, and now we can see the same relation also for general a 0 . For a large number of periods, N 1, P 22→2 dir gives the dominant contribution, as the other terms only scale as N . For a 0 1 we recover from (67) the a 2 0 and a 0 0 terms in the general expansion (53) (with ζ = 1 for this field). For a 0 1 the exponent in (67) scales as where 4/b 0 can be interpreted as the number of photons that need to be absorbed to produce a pair. Let us compare with the SLAC experiment [4]. The most straightforward result to compare with is their logarithmic plots of the number of detected positrons as a function of 1/χ, which they fit to exp(−c/Υ) 3 , where Υ = χ/ √ 2. According to eq. (8.3) in [27], for an average of η = 0.2, where η = a 0 / √ 2, the SLAC experiment gave c = 2.4 ± 0.1(stat.) +0.2 −0.6 (syst.) For this a 0 our (67) predicts c ≈ 2.46, which is in agreement with (69). However, the errors in (69) are too large for us to really be able to confirm (67). Indeed, in [4,27] c was also shown to roughly agree with an estimate obtained using a result for Schwinger pair production by a purely time-dependent electric field. Also, this value of a 0 is quite small, so the exponent is close to the perturbative one in (68). It would therefore be interesting to compare our (67) with future trident experiments with larger a 0 . If the number of oscillations are not large then we should also consider the other terms. We find P 22→1 These relations are similar to the ones we found in the previous section for a Sauter pulse. For a 0 1 we recover the a 0 and 1/a 0 terms in (54) for P N =1 22→1 , and (56) for P 22,N =1 ex (ζ = 1 for this field). Once again P 11 and P 12 are smaller by a factor of χ and can therefore be neglected to leading order.  We can also obtain these expressions by starting instead with (B2) and (B3). We have again saddle points given by (66): For the exchange term (B3) we have saddle points for n 1 = n 2 , which give P ex = P 22 ex with P 22 ex as in (70). For the direct term (B2) we have saddle points for all n 1 and n 2 , but again because of the step functions only the saddles with n 1 ≤ n 2 contribute. We find where P N =1 22→2 is given by (67). which agrees with what we found with the first approach.
IX. a0 1 AND GENERAL χ In the previous couple of sections we have studied analytically the low-χ regime. In this section we will consider the dependence on χ up to large values of χ. To do so assume a 0 1 and integrate the integrals in the LCF approximation using the numerical methods described in Sec. X. In Fig. 6 we have plotted the five terms that contribute to the total one-step rate R LCF 1 as a function of χ, where the rate is obtained from the probability as in (41). (Recall that there are only five and not six terms because P 12 dir vanishes upon integrating over the longitudinal momenta.) In Fig. 6 we see that the direct part of the one-step, R 1,dir , is negative for small χ (in agreement with the χ 1 approximations), grows in magnitude until it reaches a minimum, then starts to increase and eventually becomes positive. This dependence is nothing new, it can already be found in [5] 4 and the fact that R 1,dir becomes positive at large χ can also be seen from the χ 1 approximations in e.g. [2] 5 . Our results for the exchange part, R 1,ex , on the other hand, are new. We have shown in the previous sections that the exchange terms can be important at low χ. Fig. 6 shows us that the exchange terms continue to be non-negligible even up to quite large χ. In fact, there is a whole range {17 χ 26} of moderately large χ, around the region where R 1,dir changes sign, where |R 1,ex | is even larger than |R 1,dir |. For χ above this interval, R 1,dir eventually becomes much larger than |R 1,ex |, as expected. If we would increase χ even further, then at some point the α-expansion is expected to break down [28,29]. 4 When comparing our results with [5], note that we have plotted R 1 which is related to R 1 /a 0 by a factor of χ. 5 We have checked analytically that the χ 1 limit of P 22→1 dir + P 11 dir agrees with the χ 1 approximation in [2]. Of course, for a 0 1 these one-step terms are corrections to the two-step term R 2 . In Fig. 7 we have plotted R LCF 2 as a function of the two local values of χ, i.e.

X. NUMERICAL METHOD
In this section we describe the numerical methods we have used to integrate the different terms in P(s). The terms to be computed involve up to four phase point φ i integrals whose integrand has the general structure of a relatively slowly varying pre-exponential factor multiplying a fast oscillating phase factor. The phase generally grows with φ 2 , φ 4 and decreases with φ 1 , φ 3 . All terms but P 22 ex (s) are built up multiplicatively from independent factors describing each individual process. While lightfront-time ordering prevents the full factorization of integrals, factorization of the integrand can still be exploited to greatly reduce the computational effort by the equivalent of up to two quadrature dimensions. An important issue to tackle is the presence of the i -prescriptions originating from the infinitesimal damping of the transverse momentum integrals required for their asymptotic convergence. These i -prescriptions allow for the avoidance of singularities that are present in the remaining integrals. In all terms except P 22 ex these singularities are found at the origin of the separations between conjugate phase points, i.e. at θ 21 = 0 and θ 43 = 0. Interestingly, in studying P 22 ex one has to deal with a singularity moving with s, corresponding to vanishing d 0 , which appears in the pre-exponential factor and also makes the phase diverge. While taking the 0 limit may be defined in a mathematically sound way, the elimination of the singularity is required for numerical computation. Regularization can be done in several ways. Two classes of methods, both with advantages and disadvantages, that we have successfully tested are: A) Subtraction of vanishing or easy-to-compute quantities sharing the same singular part of the Laurent series (or a partial integration to the same purpose), resulting in a pre-exponential factor without any singularity. This requires more analytical effort in order to produce series expansions near the singularity points, where precision loss obviously prevents a direct computation of the difference. For integrals with only one effective mass in the phase, such as P 11 , P NLC or P BW , regularization may be achieved by a simple subtraction of the vanishing a 0 = 0 integral. For factored terms like P 22 dir this procedure can also be used, but by now expressing the product of two factors forming the integrand where the subscript 0 indicates a field-free (a 0 = 0) counterpart. Apart from the last, vanishing, fully free term, and the first term, describing interaction with the field at both steps, there are two additional terms which we can think of as "semi-virtual", in the sense that one of the processes takes place outside the field with no absorption of photons from the field. Note that they only survive theta integration due to the cut-off in the (θ 21 , θ 43 ) integrals arising from lightfront time ordering, and hence they only contribute to P 1 .
A similar procedure, which is numerically preferable as it does not introduce different phase factors (corresponding to a 0 = 0 and a 0 = 0), is achieved by replacing θ 21 , θ 43 with Θ 21 , Θ 43 as variables, apart from φ 1 , φ 3 , and noticing that we can fortunately rewrite the intial θ ij -dependent step functions as a combination of Θ ij -dependent ones, such that the integration domain remains triangular in shape in the new variables. Then the same division of AB is performed where now A 0 and B 0 have the same phase factor as A and B, and since the φ 1 , φ 3 integrals in the first, largest dimension integral do not involve the phase factor, they cost us the resources of only one integral with about as many oscillations as the pulse has. Once these integrals are computed, in a second step one can generate the whole s spectrum (using adequate integration routines that interpolate only the pre-exponential factor in Fourier integrals) and (breaking the integrand into a linear combination of terms factored into products of contributions depending on just one pair of φ values, with s-dependent coefficients) even compute s-integrated quantities such as the total probability or energy averages in a very fast way, through mere two-dimensional quadratures (the s integrals are fast to compute compute by part analytical, part numerical integration). The additional analytical effort and the numerical effort to change variable to Θ ij are thus greatly rewarded, by a huge speed-up.
Unfortunately the change of variables works for all terms but P 22 ex (s), where the more complex structure of the phase prevents this important additional simplification. Also regularisation by subtraction is very laborious to apply in this case. However, the fact that the dependence of the phase of P 22 ex (s) on the field only arises through effective masses and is monotonic can be quite useful, as it may allow here, too, for the fast-oscillating exponential to be separated into only one Fourier integral of a function that involves slowly varying integrals. B) Complex contour deformation generated by φ 2,4 → φ 2,4 + i /2 and φ 1,3 → φ 1,3 − i /2 with a finite . This is essentially the i -prescription we used to regularize the transverse momentum integrals, except that now we let be finite instead of infinitesimal, which means potentially appears in all functions of φ i and not just in the denominators of pre-exponential factors (like the 1/(θ 21 + i ) 2 factor in (19)). Note, though, that this complex deformation does not affect the arguments of the step functions, which hence do not spoil the analyticity. If the analytic properties of the field allow for this contour deformation, then this method provides an elegant alternative to A).
As our formulas are valid not only for pulses and periodic (e.g. monochromatic) fields but also for constant fields, method B) can also be useful in the LCF regime. A suitable choice of for the complex deformation makes the convergence of LCF terms fast, as wildly oscillating integrands are converted into fast-converging ones with only a few oscillations, making the drawing of high-resolution plots like Figs. 3 and 4 a short matter, even on the older-generation personal computer that we have used. When considering pulses, this deformation may not seem like a general enough method as it assumes a certain analyticity of the field, but remember that one can always approximate a given field by a suitable analytic one. For instance we could truncate the pulse's expansion in a basis of Hermite polynomials.
For a pulse, one finds an important advantage brought by contour deformation, apart from the great help it brings in offering a comfortable means of regularization for all terms, P 22 ex (s) included. We refer to the fact that it amplifies the integrand in the interaction region and greatly diminishes the "tails" seen in the dependence on the separations θ ij from the asymptotics of the effective mass that follow us even outside the pulse. Still, the low amplitude oscillating tails stay there and their decay is not as fast as for the LCF approximation, so for a quick and precise integration they are best separated and integrated as such, by a mix of analytic integration and a 90 degrees complex contour deformation, turning oscillating exponentials into decaying ones. The intermingling of all four phase points makes applying this optimization procedure less straightforward for P 22 ex . If one is interested in total quantities, such as the total probability, average energy of the pair and so on, method A) works efficiently for all terms, but P 22 ex . For P 22 ex one may find it simpler to use a complex contour deformation of the s variables. We found such a deformation, defined piecewise with several cases depending on the relative values of φ i , that makes s integration fast, with few oscillations. One still has then to regularize the remaining singularities appearing when some of the phase points φ i merge, through either subtraction or complex deformation of φ.

XI. CONVERGENCE TO THE LCF LIMIT
In this paper we do not present numerical results for total quantities for a pulsed plane wave, for which the methods of type A) and the one just mentioned have great potential; we reserve such investigations for a future, more numerically oriented paper. Here we will instead apply these numerical methods to local quantities, namely the rates, and study how the LCF approximation is reached for pulsed plane waves.
In particular, we consider a short pulse and compare the rates with their LCF approximations, for a set of increasing a 0 /decreasing b 0 values at constant a 0 b 0 . Both variations contribute to the reduction of what is commonly known as formation length and hence are expected to ensure convergence towards the a 0 = ∞ limit provided by the LCF approximation. We choose as pulse model a linearly polarized plane wave with Gaussian envelope:  We choose an ultra-short pulse length, given by T = π or 2π. The reason for this is not only that an extremely short pulse maximizes the ratio between the one-step and two-step contributions, but also that it allows us to better see tiny features in the plots detailing the convergence to the LCF limit. In Fig. 8 and Fig. 9 we plot two contributions to the one-step rate R 1 , namely R 22→1 dir and R 22 ex . In Figs. 10, 11, 12 and 13 we illustrate the effect on the rates R 2 and R 2 (s) produced by raising a 0 at constant χ = 1, for short Gaussian pulses. We find the remarkable result that coherence may work constructively on this term, increasing the rate on average and therefore the probabilities, unlike what can be seen in [21]. Apart from this increase, at spectrum level, we notice superimposed patterns of oscillations with s and a 0 , getting more frequent but smaller in amplitude as a 0 increases. For T = 2π, the oscillations show beats, unlike for T = π, where they are regular. In the s-integrated rate (shown at the peak of the pulse in Fig. 10, as function of a 0 , and all over the pulse in Fig. 13 for a 0 = 1) we notice how oscillations have been smoothed out and how convergence to LCF is, therefore, much faster. Coherence is not constructive everywhere for the rate. In the ripples seen between the peaks of Fig. 13 (left) there are points where the non-local rate is even slightly negative, unlike the LCF one. However, the global effect we see is a stark increase of the probability relative to the LCF approximation as a 0 decreases towards unity and coherence increases. Another thing we notice is the asymmetry between the two processes, already noted in the LCF plot of Fig. 7. The rate decreases faster with χ BW than with χ C . Conventional wisdom tells us that we should expect the LCF rates to offer good approximations for the exact ones when the scale at which the field varies significantly is considerably larger than a so-called coherence length. This is expected to appear in our integrals as the size of a finite range of values of the phase differences θ ij that give a significant contribution. We mean, for a start, those phase differences that appear in the phase factor of the integrand and can be seen united by a fermion line in Fig. 2. Of the four diagrams in Fig. 2, P 22 dir stands out as the one in which not all points are interconnected by fermion lines, so the contribution of the two pairs φ 1 , φ 2 and φ 4 , φ 3 is not suppressed as they move away from each other, except for the suppression due to the pulse length. However, in breaking the P 22→1 dir term away from P 22 dir we have added a step function that imposes an upper limit on the range of separations |σ 43 − σ 21 |, adding P 22→1 dir to the list of one-step terms, for which all the phase points φ i must be close to one another within the limits of the coherence length. This length depends on a 0 and b 0 but is unrelated to the pulse length, which constitutes a third scale, larger than the period that gives the scale for the pulse's variation. As this third scale becomes larger, P 2 will increasingly dominate P 1 , which is why we find short pulses most interesting.
Before we explain how this suppression outside the coherence length comes into play, let us mention that previous results tell us not to hold the blind belief that convergence to the LCF limit must be uniform or even true for all quantities, see e.g. [30,31]. For instance, [21] shows that for nonlinear Compton scattering this logic is confirmed for the average energy-momentum (as well as its higher moments) but not for the total probability. The more important contribution of high-wavelength photons to the latter justifies our approximating the integrand to the leading-order term to its asymptotic, high θ 21 , expansion (see equation (23) in [21]), rather than the low θ 21 one that emphasizes the dominant role of small wavelengths. For trident, however, the existence of a low energy threshold makes soft photons irrelevant, allowing for the LCF approximation to apply also at probability level.
Another issue is that, in general, when looking at totally integrated quantities we turn the oscillating phase factor present in e.g. (34) and (35) into a decaying function of θ ij (not oscillating for Compton/oscillating for BW). This is likely to make convergence to LCF much faster than before s integration.
As explained before, after proper regularizations of prefactors and manipulations of step functions, all terms but R 22 ex (s) can be written as Fourier transforms of some function of one or two Θ ij values. This is achieved by a change of variable, θ ij → Θ ij , which for large a 0 is highly non-uniform, with sudden leaps at particular points. At these points the function to be Fourier transformed has sharp variation with Θ ij . When the frequencies r k 2b0 are high, which happens at constant χ as we increase a 0 , only these points give a significant contribution. Of the aforementioned sharp variations the most significant is the peak one at the origin of Θ ij , but there are other such points for linear polarization, see the appendix of [16]. In Fig. 14 we see the derivative associated with this change of variable as function of the new variables for two very short Gaussian pulses, with T = π, 2π and a 0 = 4. We notice a central ridge, located at at the origin of Θ ij and some sharp knife-blade-like peaks, located around the points where a i = a j = a ij . Had we studied a monochromatic field, these peaks would have formed an infinite set, extending indefinitely in both directions [16]. For a pulse only a finite number of peaks are visible, lining up into a finite number of ridges, which increases with the pulse length. The height of the ridges between the peaks increases and eventually converges to unity as σ ij goes outside the pulse; hence for the short pulses in Fig. 14, only a few peaks are distinguishable from the corresponding ridges. As the pulse length increases, more and more peaks appear in both Θ ij and σ ij directions. The non-central ridges line up into two diverging bundles, oblique to the σ ij axis, forming the X-shaped image seen in Fig. 14. They correspond to the case where one of the pair φ i , φ j stays outside the pulse and the other, say φ j , is inside the pulse but the condition a j = a ij = 0 is approximately fulfilled. Thus, the number of ridges equals the number of zeros of a j in the pulse and, therefore, increases linearly with pulse length.
When both phase points fall outside the pulse, ∂Θ ij /∂θ ij forms a plateau at unity, into which the central ridge merges as σ ij leaves the pulse. The other ridges extend to infinity (Θ ij → ±∞ when |σ ij | → ∞). Of course, we must not forget that we still have to multiply with the initial integrand, which will bring asymptotic decay in both directions, through its pre-exponential factor. If we increase a 0 we notice that, for a given σ ij , the comb-like dependence on Θ ij spreads away (it has to be scaled by 1/a 2 0 to keep the plot about the same size) and becomes extremely sharp at the tip in the center of the bundle, even for moderately large a 0 . Taking into account just the central ridge gives us the semi-classical LCF approximation. Including the effect of the other features will add oscillatory terms, as seen in e.g. Fig. 8 and Fig. 9. One notices two types of oscillations in the rates: a) those at fixed φ or (σ 21 , σ 43 ) when s or a 0 is varied, b) those at fixed s and a 0 as φ or (σ 21 , σ 43 ) is varied.
Say the pulse has finite length, lasting for φ ∈ (−L/2, L/2) and let c ∞ 2 = ∞ −∞ a 2 . The bundles are located between the intervals ±Θ ij ∈ (2 |σ ij | − L, 2 |σ ij | + L + c ∞ 2 ), as σ ij moves away from the pulse. They are at the origin of the regularly oscillating "tail" seen in the rates as their arguments, φ or (σ 21 , σ 43 ), distance themselves from the pulse center. This tail's amplitude exhibits polynomial decay, introduced by the pre-exponential factor combined with the ridges' obliqueness. In addition, we notice localized type a) oscillations, coming from the peaks, such as the one at the origin of φ, seen in Fig. 8 or 10.
All oscillations described decrease in amplitude as their frequency increases with a 0 , since in unscaled coordinates their Θ ij position moves away from the central ridge. This is mathematically interesting, as the LCF limit is reached in a non-uniform way, and thus convergence may not translate to the rate's derivatives.
For a pulse, superimposed (same σ ij ) peaks and ridges will contribute to the rate with terms of different frequencies, adding up to complicated oscillations in the spectrum after σ ij integration, cf. [30]. Compare the regular oscillations in Fig. 10, (T = π, one peak at σ ij = 0) with the beats exhibited in Fig. 11 (T = 2π, two peaks). When considering a totally integrated quantity or less sharply defined momenta, the non-local oscillations of the rates wash out. In Fig. 13 we have after s-integration only smooth ripples left to remind us of the non-local, not always positive, character of the rate R 2 .
In conclusion, as soon as we increase a 0 above unity to even a moderately high value, we expect LCF to work better on average, but not at spectrum level. The transition from the fully coherent regime of a 0 = b 0 = 1 to LCF passes by an intermediate regime, where quantum correlations generate important oscillatory terms in the rates and spectra, stemming from resonant structures originating in the effective mass. This is true for all terms but P 22 ex (s), where the oscillations with s or a 0 are much reduced even at spectrum level, due to the way the exponent mixes up s-dependent quantities and effective masses, so the integrand cannot be just written as a linear combination of Fourier transforms of s-independent terms with s-dependent frequencies.
If instead of keeping χ constant, we increase the frequencies r k 2b0 at constant a 0 , the result decays exponentially as the corresponding period decreases below the minimum "coarseness" scale at which the prefactor expressed in the variables Θ ij varies significantly. This makes the distribution concentrate near the central point s 1 = s 2 = s 3 and decay exponentially as a whole with the decrease of b 0 , as seen in the LCF plots and the low χ approximations.
A useful step for the σ ij /φ integration of terms in the coherent regime would be to separate oscillatory parts that extend outside a finite pulse. This can be done noticing that they arise from the integration regions where at least one of the points is outside the pulse. As the simplest example, for a 2D term like P C (s), it is natural to perform a piecewise change of variable from θ ij to the one of φ i or φ j that is inside the pulse. Say φ i is inside the pulse. Then the phase factor looks like A similar procedure can be applied to the regularized pre-exponential factor and, at large σ distances, the condition |σ| T > |φ i | allows us to write the integral as an expansion in powers of 1/σ ij times oscillations of frequency r k b0 .

XII. CONCLUSIONS
We have studied the trident process in pulsed, constant and oscillating plane wave backgrounds. We have derived compact expressions for the exact probability for general pulse shapes. We have used these expressions to obtain various analytical approximations that go beyond the locally constant field (LCF) approximation. The formulas presented in this paper also offer a great numerical advantage due to the reduction of the number of successive quadratures needed, achieved through partial analytical integration. Their simple analytic structure thus not only allows for a more insightful view of the process and for an easy comparison between its components and with their various approximations, but also provides a good starting point for applying efficient methods to reduce the numerical complexity, as we have explained.
The trident probability in a plane wave is separated into direct and exchange terms as well as into terms characterized by the number of lightfront time (x + ) integrals. For a constant field, such x + -integrals give volume factors, ∆x + , and then the terms proportional to (∆x + ) 2 are referred to as two-step terms, while the ones proportional to ∆x + are referred to as one-step terms [1, 2,5]. For general field shapes, our lightfront separation of the probability leads to three direct and three exchange terms, which have either two, three or four x + -integrals. These terms come from squaring an amplitude with two x + -integrals and with a photon propagator in lightfront gauge that leads to one term with θ(x + 2 − x + 1 ) and another one with δ(x + 2 − x + 1 ). In the lightfront quantization formalism, the term with δ(x + 2 − x + 1 ) comes from an "instantaneous" term in the lightfront Hamiltonian. The word "instantaneous" might suggest that the corresponding terms in the probability should be related to the standard one-step terms, at least for constant fields. However, this is not the case; we found instead that all six lightfront terms contribute to the standard one-step term. We therefore group together the six lightfront terms into P 2 and P 1 , where P 2 is simply related to the product of the probabilities of nonlinear Compton scattering and Breit-Wheeler pair production. In the limit of a constant field, the direct parts of P 2 and P 1 agree with the literature results for the two-step and one-step terms, respectively, while our results for the exchange part are new as it has previously been neglected (see [7] though).
In addition to checking our results in various limits by comparing with the literature, we have derived our results using both the Hamiltonian-based lightfront formalism as well as the standard covariant formalism. While these two formulations are expected to give the same results, this equivalence might not always be trivial or obvious, see e.g. [18]. So, our results provide one more explicit example of this equivalence.
In addition to recovering previous analytical results for the direct terms for constant fields, our approach has also allowed us to go beyond these known results: We have obtained various simple analytical approximations for both the direct and the exchange terms for non-constant fields. By considering non-constant fields all the integrals are finite and well behaved, and so we avoid large volume factors ∆x + . The terms that for constant fields would be proportional to (∆x + ) 2 and ∆x + , are instead distinguished as the terms that scale as a 2 0 and a 0 , respectively, for a 0 1 and constant χ. These can be seen as the first two terms in a derivative expansion, and our approach allows us to calculate higher-order corrections.
For large a 0 , the dominant contribution comes from the two-step term, which is simply obtained by gluing together the individual probabilities of nonlinear Compton scattering and Breit-Wheeler pair production, and then the exchange and one-step parts of the probability only give small corrections. This is of course what makes PIC codes based on three-level processes useful in describing e.g. cascades in high-intensity lasers. However, the trident process itself might be most interesting in regimes where one has to take into account corrections to the two-step part, e.g. from the exchange terms. For sufficiently large χ one might expect the exchange terms to be small. However, for small χ and large a 0 we have shown analytically that the exchange part is in general on the same order as the direct part of the one-step P 1 . By considering some simple field shapes we have also shown that for χ 1 and a 0 ∼ 1 the exchange part can even be on the same order of magnitude as the total probability. Further, by numerical integration using complex deformation of x + -integrals we have also shown that the exchange terms continue to be important for P 1 even for quite large χ.
We have also studied how the exact probability converges to the LCF approximation in the limit of large a 0 . In the rate we found oscillations around the LCF approximation with decreasing amplitude but increasing frequency.
In a follow-up paper we plan on exploring trident for more general fields and parameter regimes, in particular by applying the numerical methods described here to the total probability. We believe that the methods that we have used here could also be useful for other higher-order processes in strong laser fields, such as nonlinear double Compton scattering [19,32].

XIII. ACKNOWLEDGEMENT
We thank Anton Ilderton for many valuable and extensive discussions. We also thank Ben King for useful discussions. For this work V. Dinu has been supported by the 111 project 29 ELI RO financed by the Institute of Atomic Physics A. G. Torgrimsson is supported by the Alexander von Humboldt foundation.
agrees exactly with what is stated in [33] for the cross section for single-photon trident. To recover also the overall coefficient for the cross section, σ, we replace the field according to (recall that a factor of e has been absorbed in our definition of the field) a(ω) → e µ 2πδ(ω − ω )/ √ 2ωV 3 and divide the probability by the temporal volume and the initial flux density, which is given by 1/V 3 in this case. We obtain which is exactly the result in [34][35][36].
In (C4) we have included two orders more than for the other terms; these extra terms have not been used when comparing with the numerics, but have been included in order to illustrate the growth of the series coefficients for this term. In Fig. 15 we have compared these approximations with our numerical results. These plots show that adding higher orders improves the approximation, and also demonstrate the accuracy of our numerical method. For the (total) direct term, the approximations are good all the way up to χ ∼ 1 and even by only including the first two orders: At χ = 1 the relative error is still only |R 1,dir,small χ /R 1,dir − 1| = 0.1, 0.03 for the leading order and the leading order plus the next-to-leading-order correction, respectively. The corresponding values for the exchange term are |R 1,ex,small χ /R 1,ex − 1| = 0.3, 0.06. The higher-order corrections for the exchange terms intersect at χ ∼ 0.5 where the relative error is |R 1,ex,small χ /R 1,ex − 1| ≈ 0.01. We conclude that even the leading-order approximations give good order-of-magnitude estimates even for χ ∼ 1.