The trident process in laser pulses

We study the trident process in laser pulses. We provide exact numerical results for all contributions, including the difficult exchange term. We show that all terms are in general important for a short pulse. For a long pulse we identify a term that gives the dominant contribution even if the intensity is only moderately high, $a_0\gtrsim1$, which is an experimentally important regime where the standard locally-constant-field (LCF) approximation cannot be used. We show that the spectrum has a richer structure at $a_0\sim1$, compared to the LCF regime $a_0\gg1$. We study the convergence to LCF as $a_0$ increases and how this convergence depends on the momentum of the initial electron. We also identify the terms that dominate at high energy.


I. INTRODUCTION
The trident process, e − → 2e − + e + , in laser fields has recently attracted renewed interest [1][2][3][4][5][6] following [7][8][9] and the classic experiment at SLAC [10]. While this process was studied already in the early seventies [11,12], it is only with modern high-intensity lasers that it is becoming experimentally important. Although the SLAC experiment is already two decades old and there today exist much stronger lasers, it is still the only experiment directly relevant for this process 1 . This, however, will likely change soon, given the great deal of interest in this process in itself [14][15][16][17] and as a first step towards production of many particles in cascades [18][19][20]. A related second-order process is double nonlinear Compton scattering [21][22][23][24][25][26][27][28], where the electron instead produces two photons.
We use units with c = = m e = 1, where m e is the electron mass, and we absorb a factor the charge e into the background field, so e only appears explicitly via the coupling α to the quantized field. The trident process depends on a 0 = E/ω, where E is the field strength and ω the characteristic frequency of the plane-wave background field, and χ = a 0 b 0 , where b 0 = kp is the product of the wave vector of the plane wave, k µ with k 0 = ω, and the momentum of the initial electron p µ .
The trident process can be separated into one-step and two-step parts, where the latter is obtained from an incoherent product of nonlinear Compton scattering and nonlinear Breit-Wheeler pair production. If the intensity of the background is sufficiently high then one can approximate the trident probability by the two-step where the field is treated as locally constant at the two steps. This locally constant field (LCF) approximation greatly simplifies the calculation and also makes it possible to approximate complicated higher-order processes by a sequence of first-order processes using particle-incell codes [29][30][31][32]. However, the LCF approximation can break down at both low and high energies [33][34][35][36]. Moreover, this also means that one misses much of the additional structure in the trident probability that is not already in the LCF approximation of the two first-order processes. This gives motivation for studying fields with moderately high intensity.
In this paper we consider the trident process at moderately high intensity, i.e. a 0 1. In this regime we can neither make a perturbative expansion in a 0 nor employ the LCF approximation, which is basically an expansion in 1/a 0 1 [1]. As we showed in [1], for a 0 ∼ 1 the onestep terms are in general on the same order of magnitude as the two-step term, and then one in general has to include the exchange term, i.e. the cross-term between the two terms in the amplitude that are related by exchanging the two electrons in the final state. The exchange term is much harder to calculate than the direct terms, i.e. the non-exchange terms 2 . So, it is important to know in which parameter regimes one can neglect the exchange term. We study this here.
While the trident process 3 has been observed in [10], there the laser only had a small a 0 and they found a probability that scaled as a 0 to the power of the number of photons that need to be absorbed to produce a pair. While that observation was already nonlinear in a 0 , for higher intensities the dependence on the field strength becomes more interesting. For a 0 1 and χ 1 one 2 Note "direct" does not mean "instantaneous" or "one-step". We define "direct" as the complement of the exchange term, and it gives one part of the one-step but also the two-step. 3 Note that [10] used different notation for the different contributions to this process. Here trident includes the incoherent product of nonlinear Compton and Breit-Wheeler as well as everything else.

arXiv:1912.11017v1 [hep-ph] 23 Dec 2019
finds a probability that scales as [11,12] P ∼ exp − 16 3χ . (1) This scaling is basically the same as for the nonperturbative Sauter-Schwinger process [37][38][39], cf. [10], where the probability of spontaneous pair production by a constant electric field scales as Since in the trident case we have χ = E in the rest frame of the initial electron, the difference between (1) and (2) is basically just a numerical factor. There is a similar scaling for nonlinear Breit-Wheeler pair production [40][41][42] where χ γ = a 0 kl is proportional to the momentum l µ of the initial photon. However, in this case there is no frame in which χ γ is equal to the field strength, it always contains the frequency of the initial photon. The advantage of Breit-Wheeler, as well as other photon/fieldassisted mechanisms [43,44], is that one has production of matter from an initial state with only photons. The advantage of the trident process is that one finds an exponential scaling with an even closer resemblance to the elusive Sauter-Schwinger process. This provides further motivation for studying the trident process with upcoming high-intensity laser facilities, such as LUXE [15,16] and FACET-II [17]. This can be generalized to inhomogeneous fields. In the trident case we find [1] for a 0 1 and χ 1 P ∼ exp (−f (a 0 )/χ) with some function f which depends on the field shape, and in the Sauter-Schwinger case one finds, for a time-dependent electric field with frequency ω and E 1, P ∼ exp (−F (γ)/E), where the dependence on ω is in this context usually expressed in terms of the Keldysh parameter γ = ω/E, which corresponds to 1/a 0 .
In our previous paper [1] we studied the trident process in plane-wave background fields with the lightfront approach and by integrating over the transverse momentum components as well as summing over all the spins. The advantage of this approach is that the result only depends on relatively few parameters, viz. the field parameters and the longitudinal momentum components. This approach also allowed us to derive relatively compact expressions for the probability for arbitrary field shapes. In this paper we will use these expressions as a starting point for further investigation into the relative importance of the various terms.
This paper is organized as follows. In Sec. II we give some basic definitions. In Sec. III we study numerically the trident spectrum for a pulsed oscillating field. We study the dependence on pulse length, energy parameter b 0 and a 0 . We have studied the contribution from all one-step terms, including the difficult exchange term, for a 0 1 i.e. beyond the LCF regime. In Sec. IV and Appendix A we describe our numerical methods.

II. FORMALISM
We use lightfront variables defined by The field is given by a ⊥ (φ). The initial electron has momentum p µ , and the final state electrons and positron have momenta p µ 1,2 and p µ 3 , respectively. The quantities we are interested in here only depend on the longitudinal momentum components of the fermions, which we denote as b 0 = kp and s i = kp i /kp. We also use q i := 1 − s i , i = 1, 2 for the longitudinal momentum of the intermediate photon.
We consider either the integrated probability P or the longitudinal momentum spectrum P(s), which are related via These are obtained by performing the Gaussian integrals over the transverse momenta p 1⊥ and p 2⊥ . We can perform these integrals analytically for any field shape, so it is a natural step in order to reduce the number of independent parameters. From an experimental point of view, it is also worth noting that in an electron-laser collision, the initial electron will in most cases have a very large γ factor. This means that the produced particles will travel almost parallel to the initial momentum. So, although we have integrated over all possible values of the transverse momenta, the dominant contribution comes from a relatively small region around the forward direction, see Appendix B. So, performing these transverse integrals effectively gives what a detector around the forward Gaussian peak would measure. So, for high energies the ratios of longitudinal momenta s i also give us the fraction of the initial electron's momentum/energy given to the produced particles. The amplitude has two terms, M 12 and M 21 , where M 21 is obtained from −M 12 by swapping place of the two identical electrons in the final state. On the probability level |M 12 | 2 and |M 21 | 2 give what we call the direct part P dir , and 2Re M * 21 M 12 gives the exchange part P ex . Note that "direct" does not mean instantaneous. In fact, P dir = P two + P dir one gives the entire two-step P two as well as one part of the one-step P dir one . The exchange term only contributes to the one-step P ex = P ex one . These three terms can be compared with the results from other approaches. In our approach we find it convenient to further split P dir one and P ex one each into three parts, P one dir = P 11 dir + P 12 dir + P 22→1 dir (5) and P ex = P 11 ex + P 12 ex + P 22 ex .
All the contributions to the trident probability are given by the expressions in Appendix B.

III. NUMERICAL STUDY OF THE TRIDENT SPECTRA
In this section we will present a numerical investigation of the trident process for oscillating fields with Gaussian envelope, We may vary the polarization, between linear (ξ = 0) and circular (ξ = π/4), as well as the CEP, φ 0 . We focus on the choices ξ = 0 and φ 0 = 0 for they give more structure, provide a tougher numerical challenge and strengthen the importance of one-step terms. To reveal the importance of pulse length we have considered both an ultra-short pulse (T = π), which maximizes the contribution of the one-step terms, and a longer, realistic pulse (T = 80), motivated by state-of-the-art lasers like the ones to be used in LUXE [15] and FACET-II [17]. It might be useful to consider a somewhat shorter pulse if, as in the SLAC experiment, the electron and laser beams collide at a nonzero angle, which would give a shorter effective pulse length. For a 0 ≥ 1 and χ < 1 we also notice that, unless a pulse with a flat-top envelope is used, the exponential suppression of the rates associated with lower local χ(φ) = a (φ)b 0 values means that only those peaks that are close to the global maximum contribute significantly, hence the field has a shorter effective length. Fig. 1 shows the s 1 = s 2 and s 2 = s 3 cross-sections through the probability density P (s) for the short pulse, detailing all terms that contribute, their LCF approximation and the improved LCF+1 approximation. LCF+1 is obtained by including the next-to-leading order correction in an expansion in 1/a 0 , cf. [1,45]. The values a 0 = (1, 2, 4) and χ = a 0 b 0 = (0.5, 1, 2, 4, 8) are considered. Fig. 2] shows the same for the long pulse, apart from the two-step term which can be close to the total probability density and was removed to reduce clutter. Only the one-step terms are shown explicitely. Fig. 3 shows 3D plots of the spectrum for the long pulse, scaled by dividing by a 2 0 . For Fig. 2 and 3 we go up to χ = 16, as a larger b 0 value is needed to find an important contribution of the one-step terms for a longer pulse. In Fig. 4 we considered subunitary a 0 values for the short pulse, scaled by dividing by a 2 0 . To produce the plots in Fig. 2] we approximated the computationally intensive but small P 22 ex (s) by its LCF approximation. The short pulse results of Fig. 1 and the tests we made for longer pulses suggest that this leads to a small error for the a 0 = 1 plots and a very small one for higher a 0 . Compare this to the larger discrepancy caused by, for instance, neglecting the focusing of the laser or by the imperfect control of the pulse shape allowed by experiments. Also notice that P 22 ex (s) does not posess those high amplitude oscillations which make P 22→1 dir (s) larger in magnitude than its LCF values for a long pulse, even at low b 0 . We feel that the modest loss of precision of the result due to the approximation we made hardly justifies the use of the great computing power needed to precisely plot P 22 ex (s) on the detailed grid, with many s points. Considering a long, linearly polarized, pulse with subunitary a 0 is a computationally intensive challenge for the future, as P 22 ex (s) can be very important in that case and has to be computed accurately for a good precision of the total result. The contribution to the spectrum of all the terms but P 22 ex (s) can be computed precisely at arbitrarily high resolution, through the two/three-tiered approach A 2 , building intermediate tables to be reused again and again. Thus we were able to produce the high resolution 3D plots of Fig. 3 for the long pulse.
We next present the main conclusions we have drawn from our numerical investigations. The "glued-up" term P 22 dir (s) is the most important, clearly dominates for a long pulse unless a 0 is subunitary. Our interest in the gluing techniques comes from the conjecture that this can generalize to other processes, see [27,46], including the ones involving the emission of hard photons. While good in the LCF limit, the division into steps P 22 dir (s) = P 22→1 dir (s)+P 22→2 dir (s) adds spurious oscillations to P 22→1 dir (s) and P 22→2 dir (s) outside LCF. They can be seen in the dependence on both a 0 and s. Integrating over s smooths out the former. The same oscillations were seen for the corresponding "rates" defined in [1]. Hence, as a 0 increases we see much faster convergence for the undivided P 22 dir (s) than for P 22→1 dir (s) and P 22→2 dir (s) or the "rate". Thus, away from the LCF limit it is natural to study P 22 dir (s) as whole. However, there is a way to divide integrals discussed in Appendix A, I = I aa + I a0 + I 0a , which is elegant, efficient and smooth, adding no such oscillations to its terms. It originates in the need to obtain formulas that are regular without the infinitesimal complex shift implemented by i . We prefer to use a more computationally friendly alternative, which only keeps one exponential in the integrals. The more physically meaningful way to do it involves subtracting field-free counterparts from the factors in the integrand coming from each step, allowing us to interpret the three terms as describing an interaction at both steps or just one of them.
The frequency of the oscillations we mentioned increases with a 0 , while their amplitude slowly decreases. Unlike for the short pulse, for the long one they are frequent and very irregular, which can be linked to the many sharp peaks of the derivative introduced by the change of variable θ ij → Θ ij or to the large number of complex saddle points with important contribution. All this structure exists for linear polarization, while for circular polarization (not illustrated), the curves are much smoother and the computation times needed are significantly shorter, as Θ ij is a much smoother function. The rest of the one-step terms converge more quickly to LCF than P 22→1 dir (s), for linear polarization, as they lack ar- tificially introduced oscillations. Looking at Fig. 2 we see how much smoother the oscillation structure is for P (s) than for P 22→1 dir (s). The latter's oscillations are so frequent as to even become hard to render within the resolution of the plot for larger a 0 . We chose linear polarization and a symmetric pulse in order to have a worst case scenario for the consequence of the P 22→1 dir (s) term, as the symmetry makes P 22→2 dir (s) factorize completely into half the product of independent Compton/Breit-Wheeler probability densities, thus conserving intact and combining all their oscillations. This maximizes the importance of P 22→1 dir (s) compared to P 22 dir (s), whose oscillations are partially washed out by the different choice of step functions that sever the Θ ij integrals. Check [27] for the oscillations shown by the Compton spectra for linear polarization.
In the expansion in powers of a −1 0 , whose dominant terms are the LCF results (order a −2 0 for P 2 (s) and a −1 0 for P 1 (s), dashed lines in the plot), there is also an O (1) term coming from P 2 (s). Including the latter we obtained an approximation called LCF+1, cf. [1,45], which we plotted using dash-dotted lines. As seen from Fig. 1 and 2, the LCF+1 approximation is generally larger than the LCF one and, at small b 0 , it provides an obvious improvement. As b 0 increases for fixed a 0 both approximations tend to grow, so they end up greatly overestimating the exact result. Hence, for a brief b 0 range LCF can, in fact, provide a better approximation than LCF+1. This relates to a well-known behavior of asymptotic series, for which adding more terms proves helpful in a closer vicinity of the limit, but counterproductive further away from it. We must keep in mind that, for each of the terms forming P (s), the dependence on b 0 appears through the frequencies found in the exponential, such as r 1 / (2b 0 ). Therefore, the closeness to the LCF limit also depends on s, i.e. on the momenta (cf. [33][34][35][36]). The larger the frequencies, the closer we are to that limit. Looking at them we see that in the vicinity of the point s 3 = 1 we are closer to LCF than near the points s 1 = 1 or s 2 = 1. It is near the latter points that P 22→1 dir (s) starts to first grow as b 0 is increased, eventually even dominating P 22→2 dir (s). As b 0 becomes large, almost all of P (s) comes from the narrow peaks P 22→1 dir (s) has near s 1 = 1 or s 2 = 1, a fact also due to the presence of the factor q −2 1 in the integral. To study the dominant contribution at large b 0 and a 0 ∼ 1 we can neglect both the two-step term P 22→2 dir and the first term in the division I = I aa + I a0 + I 0a mentioned in Appendix A, so we only have to compute two nested integrals, which makes this an easy task. All one-step terms grow relative to the two-step, including the much harder to compute P 22 ex , but P 22→1 dir dominates globally, so at large b 0 we may neglect P 22 ex when computing the probability. The next-to-leading term will be P 11 dir , which is even faster to compute. A negligible twostep term will also occur if b 0 is not large, but a 0 is small. However, in that case we have no reason to neglect P 22 ex , which can be comparable to the other terms, as seen in The P 22→1 dir dominance is also interesting from a conceptual point of view, because, together with the twostep P 22→2 dir it can be obtained from our gluing approach [27,46]. The only difference consists in what step-function combination is put in the lightfront time integrand, see (B6). So the gluing technique provides two things: On the one hand it provides a beyond-LCF generalization of what we call the N-step part (the cascade part) of tree-level diagrams with N final state particles. This is done through a top-down treatment of the spin/polarization of the intermediate particles. On the other hand, we also provide a large part of the rest of the probability. How this compares with the remaining terms depends on the process and on the different parameters involved. Knowing that for the trident process in a long pulse, a 0 ≥ 1 and any b 0 the sum of the two terms, P 22 dir , is dominating the others is also important for future generalizations to higher order processes.
When at least one of the frequencies becomes too large, the result is exponentially suppressed. Due to this the P (s) distribution covers a small region at the centre of the s triangle for low b 0 , which expands more and more to the sides as b 0 grows. For a 0 > 1 the extent of this region is essentially controlled by the value of the χ parameter, as for the LCF approximation. But notice that a small characteristic frequency also leads to a reduction in probability density. This explains the depression found in the middle of the triangular distributions found at the lower left corner of Fig. 3. It may also be a good idea to look slantwise at Fig. 1 and 2, to compare the plots that correspond to the same b 0 values. In this way we do not risk to overestimate the speed of convergence to LCF as the intensity (proportional to a 2 0 ) is increased at fixed collision energy. Fig. 5 illustrates, for the short pulse, the transition form the one-step dominated regime to the two-step dominated one, by varying a 0 at fixed χ for a symmetric distribution of the lightfront momentum between the three particles. Notice, though, that for subunitary a 0 , the comparison with LCF and, thus, the χ parameter lose their relevance. It is best to consider the dependence on a 0 at constant b 0 . If b 0 is large, a small χ no longer implies an exponential suppression of the probability. While at large a 0 and constant χ, , at low a 0 and constant b 0 the opposite is true, as a direct Taylor expansion of the integrand shows that P 2 (s) = O a 4 0 and P 1 (s) = O a 2 0 . Fig. 4 shows us, for the short pulse, how the two-step term loses importance compared to the one-step term, as a 0 decreases. All plots for a given b 0 have the same scale after dividing by a 2 0 , so all the one-step curves and the one for the total result tend to a definite limit, while the two-step curve tends to zero, as O(a 2 0 ). Thus P (s) ≈ a 2 0 P 0 (s) at low laser intensities. The limit P 0 (s) is independent of a 0 . For a monochromatic field there is a threshold at b 0 = 4. For the pulse we consider here P 0 decays exponentially as b 0 < 4 decreases. 4. Section of the spectra for b0 = 4 and T = π. The |P 1 (s) |/P 2 (s) ratio decreases not only with a 0 , but also with pulse length. However, if the polarization is linear, the aforementioned oscillations make this decrease significantly slower for P 22→1 dir (s) than for the other terms making up P 1 (s). Only after one averages/integrates over s, they may become comparable. Pulse length also influences the validity of the aforementioned low a 0 expansions. The bigger it is, the smaller a 0 has to be for such perturbative expansions to work well. In fact, while an experimentalist's biggest challenge may be to compress a pulse as much as possible, to make it both ultrashort and ultra-intense, the computational challenge is made higher by a longer, more diluted pulse of subuni-tary a 0 , where all terms, two-step and one-step alike, are comparable and neither a perturbative nor an LCF approach can help.
The P 22 ex (s) and P 11 dir (s) terms are the closest to LCF, giving a decent order-of-magnitude estimate even at higher b 0 values where P 22 dir is completely off the mark. A heuristic explanation is again related to the characteristic frequencies. For large b 0 , r 1 / (2b 0 ) can be quite small in some relevant region, but is added to the larger r 2 / (2b 0 ) for P 11 dir (s) , while for P 22 ex (s) there are larger frequency oscillations that couple all phase points φ i .
The interference terms P 12 ex (s) and P 12 dir (s) are seen to be relatively small for a 0 ≥ 1, but can be the same order of magnitude with the other terms for smaller a 0 .

IV. NUMERICAL METHODS
To efficiently compute all terms that make up the exact (non-LCF) probability density P(s) we have used two different approaches, benchmarking one against the other when possible. The two classes of numerical methods have already been mentioned in [1], where, however, only class B was used and only to compute rates, not total probabilities. We now add to the B-type methods an ingredient to boost computation speed, which proves vital for the difficult P 22 ex (s) term. We give here a general description of the ideas behind the two classes of methods. In Appendix A we explain in detail the proper implementation of method A for the most important terms. We stress the fact that, in talking about regularization, we do not mean the formulas from [1] are not well defined as such, but only that, for practical computations, an infinitesimal epsilon is not convenient.
A. (A 1 and A 2 ) This approach is well suited for computing all terms apart from P 22 ex (s) and, in particular, it is vital for an efficient calculation of P 22 dir (s) for realistic laser pulses, which have many oscillations. While being a bit involved, its complexity pays off greatly. The main idea is to change variables to the quantities Θ ij , whose linear combination appears in the exponent of the integrands, so as to reduce to a minimum the number of fast oscillating integrals and write the spectrum in the form of a linear combinations of Fourier transforms plus some residual terms involving lower dimension quadratures and special functions.
As already mentioned in [1] the residual terms are computation-friendly analogues of half-free terms that involve interaction with the background field at just one of the two steps. They are due to the Heaviside functions originating in amplitude-level time ordering and do not occur in the two-step term P 22→2 dir (s), where the domain of the Θ ij integrals has been extended to the whole real axis, or in the lightfront instantaneous part P 11 (s), where there is no step function.
For P 22 dir (s), a manipulation of the Heaviside functions allows us to obtain triangular domains for the inner double quadrature in variables (φ 1 ,φ 3 ), which precedes the (Θ 21 , Θ 43 ) Fourier transforms, and have integrands made up of independent factors that can be sampled independently on a grid, thus keeping the numerical complexity of the (φ 1 ,φ 3 ) quadrature close to that required by just one single integral. Then, instead of performing the outer (Θ 21 , Θ 43 ) integrals as they are, the total phase of the complex exponential is used in method A 1 as a new variable for a last, Fourier, integral, such that the quadrature that precedes it does not take oscillations from the complex exponential, so it is faster to perform. A specialized routine is used for the last integral, which only approximates the prefactor and uses exact integration of the resulting polynomial multiplied by trigonometric functions, so very frequent oscillations in the latter can accurately be followed without the need for a huge grid.
An even more efficient variant of this method that we will term A 2 , as opposed to the direct integration A 1 , involves computing intermediate interpolation tables. These are meant to store, for a given pulse shape and intensity, a table of partial derivatives of the inner integrals and then their Fourier transforms, which can be reused to quickly generate results for any combination of initial/final lightfront momenta (b 0 and s). To use an ordinary Fourier transform we take advantage of the fact that for the variable pair (Θ 21 , Θ 43 ), too, the integration domain is triangular and that it can be brought to the upper half-plane by a 45 degree rotation. Moreover, the fact that an equally spaced rectangular grid in the new variables is a subset of a similar grid for the unrotated coordinates means we can produce the interpolation tables for the Fourier integrand efficiently by an independent sampling of the aforementioned factors and their derivatives in variables Θ 21 and Θ 43 , too. All these reasons make it a very efficient method, well suited for computing results for long pulses, where A 1 will be slow even for the computation of just a few points. The region covered by the tables has a limited extent, so an asymptotic expansion of the integrand outside them has to be per-formed and its integral can be computed analytically or by complex deformation for fast convergence.
B. This method uses regularization by deformation of the complex integrals into the complex plane, using the finite epsilon prescription described in [1]. This is essential for P 22 ex (s), for which the complexity of the formula and its singularities make regularization by either subtraction of simplified singular terms or partial integration rather unwieldy for practical purposes. Had the undeformed expression of P 22 ex (s) (B7) been regular or easy to make so, an A-type method involving changing variables to the total phase and thus limiting fast oscillations to only one integral would have been convenient.
At first sight, it seems that to produce even just one P 22 ex (s) point, for a short pulse length, one needs a huge computing power. That is because the formula (B7) involves four nested integrals, with long oscillating tails towards infinity, and the integrand itself has a complex formula. When counting integrals we do not include the averages a and a 2 which we pre-compute and tabulate for all terms and methods.
The factor r d−1 in the Jacobian associated with this change amplifies the asymptotic oscillations and makes integrating to infinity unfeasible, in particular for P 22 ex (s), where d = 4. But, if one implements a cut-off in the radial direction and varies it, one sees that convergence of the total result within a very good precision is quickly reached. The reason is that if one first integrates on a hypersphere, at constant r, due to the many canceling oscillations on its surface, one obtains a result that decreases fast asymptotically, a thing unaffected by the r d−1 -proportional dilation of the hypersphere, as the density of the surface oscillations per unit solid angle also increases. After the integration over the radial direction and θ angular direction and with an appropriate choice of ε, the remaining angular integrals will have a limited amount of oscillations in them and few, if any, sign changes, which makes their computation efficient and accurate.
We want to limit the oscillations in the angular directions to a minimum and avoid small but fast canceling oscillations coming from the boundary of the ball, centered on the pulse, to which we restrict the integral. Therefore, instead of an abrupt cut-off at the boundary of the ball, we use a mollified one, multiplying the integrand by a C ∞ function of r that equals unity in this ball and slowly and smoothly goes to zero between it and a larger, concentric ball. In this way, convergence when varying the cut-off radius is faster and the computation time for a given radius is reduced.
We obtain a fast convergence even for small χ values, such as 1/2, which are in the region of exponential suppression, unless only a 0 is small and not b 0 . In this limit it is a greater challenge to compute the integral of what are initially fast oscillations, whose almost perfect cancelling brings the result down by several orders of magnitude, so the advantages of our method are even more apparent.
Especially for small a 0 , it is very helpful to subtract from the integrand its free-field (a 0 = 0) counterpart, since the latter can be proven to yield a vanishing integral, as expected. We checked that in doing so the change in the result is negligible, providing another good test for our method.
We have also checked that, for P 11 and P 12 , methods A and B give consistent results. We also checked that the exact and LCF trident rates, like the ones shown in [1], can be obtained more efficiently using spherical coordinates.
Outside the LCF approximation the rates are oscillating many times while slowly decaying outside the pulse, so it is preferable to apply method B directly to the probability terms. For LCF calculations the complex deformation introduces rapid decay at infinity, so the radial integrals are fast converging without any cut-off and the use of rates, tabulated in advance as function of χ, is a good idea.

V. CONCLUSIONS
We have studied the trident process in laser pulses with a 0 1, which mean that we have gone beyond the LCF approximation. We have found that for a 0 1 the (longitudinal) momentum spectrum has a richer structure compared to LCF. The corrections to the two-step, which we call one-step terms, are difficult to compute for long pulses, but they are anyway expected to be small for sufficiently long pulses. For shorter pulses they can be important though. We have studied all contributions for a short pulse, and found that the one-step terms can even become larger than the two-step for large electron momentum b 0 . That the one-step can become dominant is expected from the perturbative limit, because only the one-step contributes to leading order O(a 2 0 ). This is also consistent with the fact that the convergence to the LCF approximation needs larger a 0 for larger b 0 .
For a field with longer pulse length we have explicitly demonstrated that, unless b 0 is large or a 0 small, the one-step terms become negligible and so the two-step can indeed be used to approximate the probability. This is evidence in support of using our new gluing approximation [27,46] for studying other higher-order processes in long laser pulses with moderately high intensity a 0 1.

ACKNOWLEDGMENTS
We thank Christian Kohlfürst and Andre Sternbeck for advice on the use of computer clusters. We thank Sebastian Meuren for discussions about the planned experiments at FACET-II. We thank Tom Blackburn, Anton Ilderton and Anthony Hartin for useful discussions about the field shape. V. Dinu is supported by the 111 project 29 ELI RO financed by the Institute of Atomic Physics A, and G. Torgrimsson was supported by the Alexander von Humboldt foundation during most of the work and writing of this paper.

Appendix A: Numerical methods
Let us study the application of methods A 1 and A 2 to the integral I appearing in the expression: In order to make the changes of variables θ 21 → Θ 21 , θ 43 → Θ 43 we notice the very useful fact that, since Θ ij always grows with φ i and decreases with φ j , we have the implications: so we can express the time-ordering step functions as: Therefore, we break the integral into two pieces, for which we choose the new variables to be (φ 1 , Θ 21 , φ 3 , Θ 43 ) and (φ 2 , Θ 21 , φ 4 , Θ 43 ), respectively. At first sight, the pieces are complex conjugates, so it seems we only have to compute the first one. The integration domain determined by the step functions is now triangular for both the pairs (φ 1 , φ 3 ) and (Θ 21 , Θ 43 ). Hence, we would like to be able to put the integral into the form: We will however see later on that for some terms the two apparently complex conjugate integrals need to be kept separate and the way their boundary values extend to infinity must be correlated in a particular way, corresponding to a symmetric integration in the initial variables. To regularize the expression, i.e. be able to obtain a finite integrand J without further need of the epsilon prescription, we keep in mind that J is a linear combination of factorized terms of the form and Only F 0 and F 1 are singular functions at Θ = 0, from whom we plan to subtract simpler, analytically integrable functions F 0 i possessing the same singular part, using the decomposition and making sure that the last term in A8 vanishes when integrated. As explained in [1], doing this before the change of variables to Θ ij and using for F 0 i the free field (a 0 = 0) version of F i , provides an interpretation of the second and third terms in A8 as half-virtual, that is involving an interaction with the field at only one of the two steps. From a computational point of view, though, it is preferable to have only one Fourier exponential as a common factor, so as to regularize directly the pre-exponential factor. Hence we choose: and write all formulas in terms of the regularized factors Successively plugging into the linear combination J just one of the first three terms in A8 produces the decomposition: We detail them below, starting with the first, four dimensional quadrature term: where and T ij are obtained from double integrations which can be done efficiently, as they are defined on a triangular domain and the oscillations of the integrand only come from the pulse shape: where and g 6 = g 1 + g 2 . A special routine was written for computing S ij . For method A 1 it first uses a preliminary adaptive Simpson method to sample the two functions g i and g j on suitable nonuniform grids. Then follows a φ 1 −integration of the product of the piecewise polynomial approximations generated for g i (φ 1 , Θ 21 ) and φ1 dφ 3 g j (φ 3 , Θ 43 ). In this way we only sample the two functions separately on a line, as for a one-dimensional integral. For method A 2 each of the functions g i is accompanied in the above procedure by a few of its Θ ij derivatives and, thus, a whole matrix of partial derivatives of S ij with respect to (Θ 21 , Θ 43 ) is produced. This is because method A 1 is aimed at producing just one point value, while method A 2 is meant to produce a whole set of spectra, for any b 0 and (s 1 , s 2 ) values.
For method A 1 it is best to change variable to the exponent Θ = ω 1 Θ 21 + ω 2 Θ 43 and leave it for the last Fourier integral.
The procedure for A 2 is more complex, as we want to create tables of values for A13 and we start by tabulating T ij together with a matrix of partial derivatives, obtained as described above. To turn the integration domain into the upper half-plane, we make a 45 degree rotation by changing the variables to Z 1 = Θ43+Θ21 2 and Z 2 = Θ43−Θ21 2 . We then compute Fourier transforms with new frequencies o 1/2 = ω 1 ± ω 2 for each element of the matrix of partial derivatives of each term T ij and a dense set of o 1 and o 2 values so as to create up to 9 interpolation tables, for general polarization.
The factor 1/q 2 1 in A1 can ruin precision in the vicinity of s 1 = 1, but this can be avoided by a double partial integration, writing: are linear combinations of . The replacement of T ij byT ij means that, for a given maximum total order of the Z partial derivatives, we now need to add two more orders in the tables that preceed the 45 degree rotations, whose elements are linear combinations of: i (φ, Θ) = ∂ p g i (φ, Θ) /∂Θ p . We tabulate the set of partial derivatives on a grid with step h, so that This can be done highly efficiently in the following manner: In order to reduce memory storage, we concentrate at one time on just one small piece of the φ 1 integrals, covering an interval of length l, adding to the contribution of the previous intervals. Say we use derivatives up to order 4. For a grid of i 1 and i 2 values we determine, using an adaptive Simpson method, piecewise polynomial approximations for g on two different grids of φ 1 values, for p, q = 0, 1, ..., 6 and for a square n × n matrix of (Θ 21 , Θ 43 ) values, which cover only a small square part of the domain we want covered.
If n is large enough, the computation of the functions and derivatives, the adaptive Simpson procedure and the subsequent sorting of the grid points cost little compared to the n 2 φ 1 integrals that are performed, using a custom integration tool for the product of piecewise polynomial interpolations. Then the procedure is repeated for all relevant φ 1 intervals and, after implementing the rotation, we save to an external memory the values ofT ij and its derivatives on a square matrix of equally spaced (Z 1 , Z 2 ) values, one out of many that are to be computed and stored for the next step. If the grid spacing h is small enough, by a Taylor expansion, we can find with good precisionT ij at any point inside the square block where the (Z 1 , Z 2 ) values lie.
The second step involves a custom code for implementing the double Fourier transform, using, for each h interval centered at a grid point, analytical formulas for the integral of an oscillating exponential multiplied by the aforementioned Taylor series. Since we cannot compute an infinite number of blocks, we cover a finite rectangle with these blocks and use forT ij an asymptotic extrapolation as a polynomial decay times oscillation, which can be used to extend the integration to infinity, with good approximation. In this asymptotic region the integral can be expressed in terms of trigonometric integrals, or a complex rotation of the integration path may be used to turn the oscillating exponential into a decaying one.
The double Fourier transform of eachT ij is computed and stored as a table for a dense set of o 1 , o 2 values, with |o 2 | < o 1 < o max , where the o max is chosen large enough to provide all the nonnegligible part of the spectrum for the minimum χ value we have in mind. Then, for any b 0 and s values, I aa is produced from interpolations of the data in the tables. Next we describe the residual terms, choosing I a0 for more detailed explanations.
Letg i be the function obtained from g i by replacing φ 1 by φ 2 as variable independent of Θ 21 . where r2 d 2 , and d 1 = d 2 + 1. We have to write the expression like this, before going to just one integral like in A4, for an essential ingredient is finding out how L 1 and L 2 are to be related for a symmetric integration around the average phase σ ij , instead of a discontinuous choice of the variables that are independent from Θ ij . By breaking the integration domain into two parts, each with a different definition of the variables independent of (Θ 21 , Θ 43 ), we effectively severed the Θ 21 and Θ 43 integrals and reattached them with a displacement, which has to be taken into account. Hence we cannot just keep L 1 and L 2 equal as they tend to infinity, for this would give a value for the probability that can even be negative at low a 0 . Instead, since for a fixed Θ 43 value, at large φ, Θ 43 = θ 43 , to keep the integration symmetric, we should choose L 1 = L − Θ 43 /2, L 2 = L + Θ 43 /2, with L → ∞. We have checked that this A-type method gives the same results as what we obtain by integrating the rates in [1], which were obtained with a B-type method.
We introduce the following special functions: We use the above integrals, then in the second term we swap the two indices of Θ 21 , after which we change the sign of Θ 21 to keep the same exponent.
We now use the symmetric role of the two φ values, which implies thatg 0 (φ 1 , Using the notation W 0 (x) = W (x) e −ix for what is a non-oscillating function, easy to be produced from an interpolation table and series expansions near the origin and at large values, we obtain the final formulas for the residual terms: where G i = φḡ i (φ, Θ) dφ,Ǧ i = ǧ i (φ, Θ) dφ and for i = 0. We have checked analytically that for a 0 1 we find agreement between I 0a and the high-energy limit of perturbative O(a 2 0 ) trident. Since both the Re... and Im... terms in I 0a contribute to leading order, this is a nontrivial check of the choice of the boundaries L 1 and L 2 above.
Similarly to what has already been described before, method A 1 uses direct integration, while method A 2 involves the tabulation of the functions G i ,Ǧ i and a few of their derivatives, for a large region of Θ values, so as to be used repeatedly to compute the integrals for any b 0 and s i values. Here, too, a double integration by parts is very useful to prevent precision loss. We can extend the tables enough as to not need asymptotic expansions, as convergence is fast enough.
Next we want to calculate (B2) and (B3): To make the formulas very similar to the ones before, we interchange even and odd indices, so φ 1 ↔ φ 2 , φ 3 → φ 4 and the three phase points are φ 1 , φ 2 , φ 4 . The direct term needs no regularization and, therefore, it has no residual terms, I 12 dir = I 12aa dir , while the exchange one has the decomposition I ex = I 12aa where g a (φ, Θ) = w2/θ (∂Θ/∂θ) φ , and Notice that all residual terms vanish in the LCF approximation, but they can be essential at low a 0 /large b 0 , where they dominate.
The procedure to compute just the two-step term, P 22→2 dir (s), is similar to the one for P 22 dir (s), except that we use the average coordinates σ 21 = φ1+φ2 2 and σ 43 = φ3+φ4 2 as independent variables and there are no residual terms. We will not present the details, as they are easy to infer, but mention that we need only up to five tables here, for general polarization, if we make the partial integration used in Eqs. (34)-(37) in [1]. Again, a double partial integration can bring out a factor of ω −2 2 , to mitigate the singularity q −2 1 present in the formula when s 1 is close to unity. Asymptotic expansions outside the region covered by the tables are essential for the two-step terms, which have slower asymptotic decay than for the integrals encountered before and even more important for the LCF approximation to the two-step probability that can also be computed by the same method. To compute P 11 (s) we need similar interpolation tables for the Fourier integrand and its transform, but these are one-dimensional, so require fewer resources. Here, too, the integrand is approximated by an expansion outside the tables.

Appendix B: Formulas
For the readers' convenience we collect the formulas from [1] in this appendix. The total probability is given by the sum of the following 6 terms: The prefactor of (B7) can be expressed in terms of a permutation defined by We have obtained these results by performing the Gaussian integrals over p 1⊥ and p 2⊥ . For the direct terms we find that the Gaussian peak is located at and which means It is experimentally relevant to assume that γ = p 0 is much larger than all the other parameters, and then we have to leading order The peak for the exchange terms is in general more complicated, but reduces (B17) to leading order. Thus, in this regime all particles in the final state have momenta approximately parallel to the initial momentum. s i is by definition the fraction of the initial longitudinal momentum given to particle i. From (B17) we see that it also gives the fraction of the total momentum to leading order. s i are therefore natural variables to focus on. After shifting and diagonalizing the transverse momentum variables, the exponential becomes e −c1p 2 1⊥ −c2p 2 2⊥ . The coefficients c 1 and c 2 determine the width of the Gaussian peak. c i can be expressed in terms of s i and φ i and scale with a 0 and b 0 , but do depend on γ separately (b 0 is much smaller than γ). The width of the Gaussian peak is therefore much smaller than its position.