Approximating higher-order nonlinear QED processes with first-order building blocks

Higher-order tree-level processes in strong laser fields, i.e. cascades, are in general extremely difficult to calculate, but in some regimes the dominant contribution comes from a sequence of first-order processes, i.e. nonlinear Compton scattering and nonlinear Breit-Wheeler pair production. At high intensity the field can be treated as locally constant, which is the basis for standard particle-in-cell codes. However, the locally-constant-field (LCF) approximation and these particle-in-cell codes cannot be used when the intensity is only moderately high, which is a regime that is experimentally relevant. We have shown that one can still use a sequence of first-order processes to estimate higher orders at moderate intensities provided the field is sufficiently long. An important aspect of our new ``gluing'' approach is the role of the spin/polarization of intermediate particles, which is more nontrivial compared to the LCF regime.


I. INTRODUCTION
Consider a single high-energy electron that collides with a high-intensity laser field.
Many electrons, positrons and photons can be produced in such collisions [1][2][3]. One can approximate the laser by a pulsed plane wave, and, thanks to the simplicity of the Volkov solution to the Dirac equation in such fields, one can derive compact expressions for e.g. nonlinear Compton scattering e − → e − + γ for arbitrary pulse shapes and parameter values. After the first photon emission the photon can decay via nonlinear Breit-Wheeler pair production [4,5] γ → e − + e + or the electron can emit a second photon. The first gives one part of the nonlinear trident process [6][7][8][9][10][11][12][13][14][15][16][17] e − → e − + e − + e + , and the second gives one part of double nonlinear Compton scattering [18][19][20][21][22][23][24][25] e − → e − + γ + γ. Since all the processes we consider here are in general nonlinear in the interaction with the laser, we will drop "nonlinear" from here on. These two second order processes are significantly harder to calculate than the first-order processes. In fact, even in a plane wave, it is a challenge to calculate all contributions to the total/integrated probability e.g. for a long pulse. So, at third and higher orders one definitely needs some way to approximate the exact result.
One regime that allows for such an approximation is the high-intensity regime. More precisely, the classical nonlinearity parameter a 0 = E/ω should be large 1 . Exactly how large it has to be depends on the values of the other parameters in the system [26][27][28][29], but assuming it is large enough, then one can make an expansion in 1/a 0 . The formation length is small in this regime, so the field can be treated as approximately constant during particle production. So, in for example trident, photon emission would happen at one constant value of the field strength, the photon would propagate a macroscopic distance and then decay into a pair at another constant field strength. We call this the two-step part, and refer to the rest as one-step terms. In the production of N particles we use "N -step" to refer to the corresponding cascade part. This is the basis of particle-in-cell codes [30][31][32][33].
In most of the standard PIC codes so far, whether or not to produce a particle is determined based on probabilities/rates that are summed/averaged over the spin/polarization at each step. However, it is known, see [12,16,17], how the spin/polarization of intermediate particles in trident and double Compton can be included. For constant fields this is achieved with single sums for each intermediate particle. For example, in trident one would have a single sum over two different polarization vectors of the intermediate photon. Recently, spin effects have started to be included in PIC codes [34][35][36][37]. See also [38,39] for further improvements on standard LCF/PIC methods.
However, if a 0 is not large then one cannot use this LCF approach. In this paper we provide a generalization of the LCF approach to a 0 ∼ 1. For a 0 ∼ 1 one can in general expect important corrections to the Nstep. However, if the field has a 0 ∼ 1 and a sufficiently long pulse length then the N -step again gives the dominant contribution, because the intermediate particles can propagate macroscopic distances, while in the corrections to the N -step at least some of the particles are forced to stay close. Note that what we here (with a 0 ∼ 1) mean by the N -step is different from its LCF approximation, and its relation to the product of nonlinear Compton and Breit-Wheeler is more complicated than in the LCF regime and what is put into standard PIC codes. Thus, while things are simpler for a long pulse, we still have new features compared to the LCF regime. For arbitrarily polarized fields the role of the spin/polarization of intermediate particles is significantly different from the LCF case. Our motivation for considering long pulses with moderately high intensity comes from upcoming experiments such as LUXE [40,41] and FACET-II [42]. This paper is organized as follows. In Sec. II we introduce some notation and basic definitions. In Sec. III we first study the relation between the dominant contribution to trident for long pulses and the incoherent product of nonlinear Compton scattering and Breit-Wheeler pair production. Then we study how this generalizes to other higher-order processes. We present two equivalent formulations of the gluing approach. In Sec. IV we show how higher-order processes can be approximated with a saddle-point method, we show in particular how the intricate pattern in the spectrum at low energy can be understood in terms of a large number of complex saddle points.

II. FORMALISM
In order to introduce notation we present the dominant contribution to the trident probability P. The vector potential is given by a ⊥ (φ), where φ = kx = ωx + = ω(x 0 + x 3 ) and a ⊥ = {a 1 , a 2 }. We also use the notation v ± = 2v ∓ = v + ± v 3 . The initial electron has momentum p µ , the two electrons in the final state have momenta p 1µ and p 2µ , and the positron has p 3µ . Because the field only depends on one space-time coordinate, x + , we have a delta function δ 2 , which can be used to perform e.g. the integral over the positron momentum components. In Sec. III we will show how to approximate higher-order processes by incoherently gluing together a sequence of nonlinear Compton and Breit-Wheeler steps. To do that we need the spin, polarization and longitudinal momentum dependence of these first-order building blocks. However, we can integrate the first-order building blocks over the transverse momenta (P ⊥ ) without losing necessary information for building higher orders. This is related to the fact that after integrating over the transverse momenta of the final state particles, the probability no longer depends on the transverse momentum of the initial particle.
So, we perform the integrals over p 1⊥ and p 2⊥ and we are left with the longitudinal momentum spectrum where s i = kp i /kp is the ratio of the longitudinal momentum of particle i and the initial electron. For the initial electron we use b 0 = kp and for the intermediate photon q 1 = 1 − s 1 . In this paper we are interested in the two-step part of P. This comes from a term that on the amplitude level has two lightfront time x + integrals, so on the probability level we have 4 lightfront time integrals. We use φ 1 and φ 2 for the photon emission step and φ 3 and φ 4 for the pair production step. We find κ 01 κ 23 4 where and where s 0 = 1 is inserted for symmetry reasons. The singularities at θ 21 = 0 and θ 43 = 0 are regulated by φ 2,4 → φ 2,4 + i and φ 1,3 → φ 1,3 − i with > 0, or equivalent integration contours. The field enters the exponent in Θ ij := θ ij M 2 ij via the "effective mass" [43] M , which is obtained from the lightfront-time average of the Lorentz momentum as where and The prefactor depends on the field through which enters via D 1 = ∆ 12 ·∆ 21 , D 2 = ∆ 34 ·∆ 43 , and where and where ε is the Levi-Civita tensor. Note that W ij = 0 for linear polarization. The step function combination can be expressed as where σ ij = (φ i + φ j )/2. Replacing θ(θ 42 )θ(θ 31 ) with θ(σ 43 − σ 21 ) gives the same result to leading order in the pulse length or in an 1/a 0 expansion. It is natural to make this replacement because the second term in (9) scales linearly in the volume and is therefore naturally combined with the other one-step terms. Note that we derived (2) in [6] without reference to the first-order processes, nonlinear Compton and Breit-Wheeler. We showed in [6] that for linear polarization (2) can be obtained from the incoherent product of these first-order processes with a single sum over the polarization of the intermediate photon. In Sec. III we will consider arbitrary polarization, where things are more nontrivial.

III. GLUING APPROACH
In this section we present our new gluing approach, where higher-order processes are approximated by linking together the spin and polarization dependent probabilities of nonlinear Compton scattering and Breit-Wheeler pair production. This is a generalization of the case in [24], where we considered processes with only intermediate electrons. We expect this to give a good approximation for sufficiently long pulses and/or large a 0 . To treat the spin and polarization we use the following basis. For fermions we choose and (cf. [44,45]) (12) These spinors are convenient because of their simple dependence on p ⊥ 2 , so all integrals over p 1,2,⊥ are Gaussian. A general spin state can be expressed as 2 Recall that p− is independent on p⊥ while p+ = (1 + p 2 ⊥ )/(4p−), which also means that factors of p 0 = p− +p+ could lead to more complicated expressions. v = cos We assume that ρ and λ do not depend on the momentum p, and then we integrate over all the transverse momentum components. It turns out that the results can be expressed neatly in terms of the following vector ={cos λ sin ρ, sin λ sin ρ, cos ρ} , where the spin matrix is given by We also have where γ 5 = iγ 0 γ 1 γ 2 γ 3 and the spin 4-vector α µ is a linear combination of three basis vectors α (i) p = 0 with the components of n as coefficients, with (cf. [34]) wherexv = v 1 ,ŷv = v 2 and α (i) α (j) = −δ ij . Thus, n is the Stokes vector with respect to the spin basis given by α (i) . For a photon with momentum l µ we choose a polarization vector with − = 0, + = l ⊥ ⊥ /(2l − ) and ⊥ = cos We again keep ρ and λ constant while integrating over the transverse momenta. Similar to the fermion case, we find that the polarization dependence can be expressed in terms of another three-dimensional unit vector, where the Pauli matrices are as usual given by The meaning of this vector is of course different from the fermion case and it transforms differently under e.g. a rotation of the transverse coordinates. However, the role this vector plays in linking together the first order processes is basically the same as in the fermion case. n is the Stokes vector with respect to the two polarizations given by . Spin and polarization effects have been studied in many papers, see for example [34][35][36][37][46][47][48][49][50][51][52][53][54][55].

A. Averaging approach
With these vectors we find that the probability of nonlinear Compton scattering and Breit-Wheeler pair production can be expressed as with two n i 's for the fermions and one for the photon (cf. Sec. 87 in [56] for ordinary Compton scattering). P gives the spin and polarization averaged probability, P i gives the dependence on the spin/polarization of one particle when averaging over the spin/polarization of the other two particles, P ij describes the correlation between the spin/polarization of two particles, and P ijk describes the correlation between the spin/polarization of all three particles.
For Compton scattering we find that the probability that an electron with longitudinal momentum s 0 and spin vector n 0 scatters into a state with s 1 and n 1 by emitting a photon with momentum q 1 and polarization n γ is given by where R C is given, with the same notation as in (23), by and The corresponding expressions for Compton scattering by a positron can be obtained either by 1) replacing a → −a and n → −n for the two spin vectors, or 2) by replacing s 0 ↔ −s 1 (except in the overall factor of 1/s 2 0 in (24)), n 03 ↔ n 13 and n 0⊥ ↔ −n 1⊥ . For Breit-Wheeler we find that the probability that a photon with q 1 and n γ decays into an electron with s 2 and n 2 and a positron with s 3 and n 3 is given by where . These expressions for P BW can be obtained from (24) to (32) by replacing s 0 → −s 3 , s 1 → s 2 , q 1 → −q 1 , n 03 → n 33 , n 0⊥ → −n 3⊥ , n 1 → n 2 , n γ2 → −n γ2 and n γ1,3 → n γ1,3 , and finally multiplying with an overall −1. The sign change for one of the components of n γ can be understood as a consequence of the fact that, when changing an incoming photon to an outgoing one, one takes the complex conjugate of the polarization vector µ , which in (20) corresponds to λ → −λ, which in turn changes the sign of n γ2 in (21). The goal is now to link together these first-order terms to approximate higher-order processes for sufficiently long pulses or large a 0 . It might seem like we have a quite large number of terms compared to the familiar LCF case, but note that the idea is that these are all the terms we need to construct the N -th step for any higherorder process for a 0 1 and arbitrary field polarization.
We start with trident. In this case we only have an intermediate photon, but no intermediate fermions, which means that to obtain the probability summed and averaged over the spins of initial and final state particles we only need to consider four terms, namely (25), (28), (36) and (39). It turns out to be convenient to write all spin and polarization sums as averages, so as our initial ansatz we take where we have a factor of 2 3 because we have replaced sums with averages over the spins for the three particles in the final state, we have similarly included a factor of 2 for the intermediate photon, the factor of 1/2 is because we have two identical particles in the final state and the (1 ↔ 2) term makes the probability symmetric with respect to these two particles. Here and in the following we have, to avoid clutter, omitted the arguments of P C and P BW in . . . , which are of course different. Apart from the different momenta and n we have chosen φ 1 and φ 2 for the first factor/step in . . . , φ 3 and φ 4 for the second factor etc. In this gluing approach we also include step functions, e.g. θ(σ 43 − σ 21 ), so the different steps happen in the right chronological order, but this is done in exactly the same way as in [6] so we leave it implicit in the following. The average/sum over the spins (and polarization in the general case) of initial and final state particles simply corresponds to sums over two antiparallel vectors, e.g. n 1 = ±n r 1 where n r 1 is an arbitrary vector. So, we have 1 = 1 and n = 0. However, things are nontrivial for the spins of intermediate particles. We can still use 1 = 1 and n = 0 which give where n is the 'polarization vector' of the intermediate photon.
The question now is what to do with nn . If one sums over n = ±n (r) , where n (r) is some reference/basis vector, then the matrix nn clearly depends on the choice of n (r) . We want P glue to be equal to P 2 in (2), which we can write as For linear polarization a 2 (φ) = 0 we have P C,BW γ,1 = P C,BW γ,2 = 0 and then we can obtain P 2 by summing over two real polarization vectors with λ = 0 and ρ = 0, π, which correspond to n = ±e 3 . So, for linear polarization there is a choice of n (r) that gives the desired result. However, this does not work for more general field polarization where P C,BW γ,1 , P C,BW γ,2 = 0, because then we would not be able to obtain e.g. the term with κ01 2 W 12 κ23 2 W 43 in P 2 . If one instead tried to sum over circular polarizations, with ρ = π/2 and λ = ±π/2, then one would have n 3 = 0 and would again be missing terms. The fact that we in general do not recover all terms in this way is because we already have a sum over the photon polarization on the amplitude level, so on the probability level we have in general a double sum, over and say, and the naive gluing approach only takes into account the two terms where = .
However, by comparing (45) and (46) we immediately see that by simply replacing nn → 1 we obtain all terms in P 2 for any polarization. We thus propose the following improved gluing approach: Include a factor of 2 for each intermediate particle and then simplify with the rules 1 = 1, n = 0 and nn → 1. We have showed in [24] that this procedure also works for double nonlinear Compton scattering, where we have a vector n for an intermediate electron rather than a photon. We have in fact showed [24] that this procedure also works for triple and quadruple nonlinear Compton scattering, where an electron interacts nonlinearly with the background field and emits three and four photons, respectively. These processes only have intermediate electrons and are hence built from (25), (26), (27) and (29).
As an example of a process that involves both an intermediate photon and an intermediate fermion we consider double Breit-Wheeler pair production, i.e the decay of an initial photon into two electron-positron pairs as illustrated in Fig. 1. There are two different contributions to this process, where the intermediate photon is emitted by either an electron or a positron. We have checked that both can be obtained from where P p C is the probability of positron Compton scattering, we have 2 6 because we have two intermediate particles and we sum over the spins of the final state particles and average over the initial polarization, and we have a factor of 1/2 2 because there are two pairs of identical particles in the final state. The brackets in (47) are calculated using 1 = 1, n = 0 and nn = 1 for all 7 "spin/polarization vectors".
In Fig. 2 we have (one part of) a process with two cascade branches. (The two photons can of course also be emitted by the same fermion, but that is another example of a single-branch cascade.) We have checked that the absolute squared of the diagram in Fig. 2 can be obtained from So, our gluing approach is not restricted to single-branch cascades.

B. Matrix approach
Consider higher-order nonlinear Compton scattering (emission of more than one photon with nonlinear interaction with the background field). For this single-branch cascade we can replace the ... "operator" with a matrix formulation. From the 3D spin unit vector n we define a 4D vector Summing over the photon polarization we find that the probability of single Compton scattering can be expressed as with N 0 and N 1 for the initial and final state electron, respectively, and where P is defined in [24]. The last step defines a 4 × 4 matrix M in terms of P 0 , P 1 and P 01 .
As an illustration of this matrix formulation, consider triple nonlinear Compton scattering. The ... approach described in [24] can be expressed in terms of the 4D approach according to [ P + n 0 ·P 0 + P 1 ·n 1 + n 0 ·P 01 ·n 1 ] [ P + n 1 ·P 0 + P 1 ·n 2 + n 1 ·P 01 ·n 2 ] [ P + n 2 ·P 0 + P 1 ·n 3 + n 2 ·P 01 ·n 3 ] ={1, 0}·M 10 (3)·M 10 (2)·M 10 (1)·{1, 0} , where N 0 = N 1 = {1, 0} means that we are averaging and summing over the spin of the initial and final state electron, and M 10 (i) depends on the lightfront time integration variables and longitudinal momenta associated with the emission of photon i. This seems like a simple formulation for a single-branch cascade, and it is similar to the Müller-matrix formulation of the evolution of Stokes vectors in optics, see also [49] for perturbative QED with no background field. However, in cascades with more than one branch, a 4 × 4 matrix would not describe the most general spin and polarization dependence, instead the dimensionality increases with the number of branches.
We will now demonstrate that these two formulations are equivalent. The spin and polarization dependent parts factorize because for a fermion propagator we can express / p + 1 = ρ=ρ0,ρ0+π uū (52) and where the spinors u and v are given by (13) and (14). And for the photon propagator we have where the polarization vector µ is given by (20). Consider a Compton scattering step, where a photon is emitted by an electron. From the amplitude we have, regardless of whether or not the particles come from the initial, some intermediate or the final state, and From the complex conjugate of the amplitude we have where ρ c i and λ c i can be different from ρ i and λ i for intermediate particles. We will integrate X Y over the transverse momenta. Assume first that we are dealing with a final-state step, where no more particles are produced by these particles (but more particles can be produced at a later lightfront time by a different cascade branch). We have three cases: 1) We use the overall momentum conservation delta function to perform the p ⊥ n integral, which means p ⊥ n = p ⊥ m − l ⊥ , where p ⊥ m depends on the other, independent momenta. Only this step depends on l ⊥ , so we have where N is a normalization factor, which we will come back to. 2) If we have already used the overall delta function for a different step, then we have integrals over both l ⊥ and p ⊥ n . These momenta also appear in some other step(s) because of the momentum conservation, but only via their sum p ⊥ m = p ⊥ n + l ⊥ . So, we change variable from p ⊥ n to p ⊥ m , and then l ⊥ only appears in this step and we again have (60). 3) If we initially have an integral over p ⊥ n and where the photon momentum is fixed by momentum conservation l ⊥ = p ⊥ m − p ⊥ n , then we just change variable from p ⊥ n to l ⊥ . So, in all cases we have (60). After performing the integral over l ⊥ we find that the result is independent of all the other transverse momenta. This means that this step factors out and we can treat the step that produced the p m electron as if it were a final step. After performing all the transverse momentum integrals we find that all steps have factorized. The different steps are linked via matrix multiplication. For an initial electron we have u α (ρ, λ)ū β (ρ, λ), with u from the amplitude andū from its complex conjugate, which we can write in terms of the basis spinors as (omitting the spinor indices) where (62) For an outgoing electron we haveū α (ρ, λ)u β (ρ, λ), with u from the amplitude and u from its complex conjugate, soū Similarly, for an outgoing photon we have¯ µ ν , with the first term from the amplitude and the second term from its complex conjugate, where µ ↑ is given by (ρ = 0, λ = 0) and µ ↓ by (ρ = π, λ = 0). For an intermediate particle we a double spin sum, one for the amplitude and a second for its complex conjugate. We can without loss of generality sum over the spin basis with λ = 0, i.e. r,s=↑,↓ u rūs . We can make this a single sum by summing over (u ↑ū↑ , u ↑ū↓ , u ↓ū↑ , u ↓ū↓ ). We order the sums for the other two particles also as (↑↑, ↑↓, ↓↑, ↓↓). The most general case is thus given by a 4 × 4 × 4 matrix Z i l imin , where i j = 1, .., 4 and the {first, second, third} index for the {photon, incoming electron, outgoing electron}. We want to express the spin and polarization of the particles in terms of their Stokes vectors (as above) rather than N (0) , and for this it is natural to transform the indices using which obeys 2T T † = 1 (T † is the conjugate transpose). Using this matrix we define If one of these particles is in the initial or final state, then we have 2T † N (0) = (1, cos λ sin ρ, sin λ sin ρ, cos ρ) = N (67) andN (0) 2T = N , which is the Stokes vector as defined above. For single Compton scattering we have where N l , N m and N n are the Stokes vectors for the photon and the initial and final state electron, respectively. So, the building block M for constructing the gluing estimate of higher orders is naturally interpreted in terms of the Stokes-vector dependent first-order process. For double Compton scattering, for example, we would have We can find similar expressions for photon emission by a positron and pair production. Because an incoming positron hasv instead of v in the amplitude, for the positron spin sums we order the terms according to (↑↑, ↓↑, ↑↓, ↓↓) instead of (↑↑, ↑↓, ↓↑, ↓↓), so that we have the same N (0) and T for all particles.
To show that this matrix formulation is equivalent to the average ... approach, start with the gluing ansatz where N is the Stokes vector of some intermediate particle. Using the prescription 1 = 1, n = 0 and nn = 1 we also have N α N β = δ αβ , so the ... approach reduces to the matrix formulation. We now return to the normalization factor N in (60). We use the same normalization as in [6]. So, in particular, the amplitude is given by where b, d and a are the mode operators (with the momentum and spin arguments suppressed) for electrons, positrons and photons, respectively, and U is the evolution operator. The delta function is given byδ(P ) = (2π) 3 δ(P − )δ 2 (P ⊥ ). P in is the momentum of the initial particle and P out is the sum of the (− and ⊥ components of the) momenta of all the final state particles. The initial state is given by where B is the mode operator for an electron, a positron or a photon, and dP in = θ(P in The mode operators are normalized such that in | in = 1 implies dP in |f | 2 = 1 .
We assume for simplicity that the wave packet is sharply peaked, which means dp where p is the momentum of one of the final-state particles. For each of the rest of the outgoing particles we have a momentum integral dP . For each intermediate particle we have i d 4 xd 4 P = i dφ dP+ 2πk+ d 3 x −,⊥ d 3 P −,⊥ . The integral over x −,⊥ gives a delta function which we use to perform the integral over P −,⊥ . The P + integral is elementary and independent of the field, where ∆x + is the difference between two x + variables, and (1, P + ) means that the numerator is either independent of or linear in P + . Using P 2 = 4P − P + −P 2 ⊥ , the integral either gives a term with an "instantaneous" δ(∆x + ) or one with a time-ordering θ(∆x + ), cf. [21]. Only the step-function term contributes to the cascade/gluing estimate. Apart from the terms that we have already included in X and Y, each propagator gives a factor of 1/(2kP ). There is one x µ integral more than there are propagators and, since d 4 x = dx + dx − d 2 x ⊥ /2, it gives an overall factor of 1/4, which we combine with 1/(kP in kp ) in (74). Thus, each external and intermediate particle gives a factor of 1/(2kP ). So, in order for P = dkP i θ(kp ) dφ i "time-ordering" all steps Z i (76) to give the probability with the correct normalization, where kP i = b 0 s i are all the independent longitudinal momentum variables of the outgoing particles and "timeordering" denote the product of step functions that give lightfront-time ordering, we need .
The fundamental matrix for a Compton scattering step can be expressed compactly as with X and V obtained from (34) by replacing w 1 → w m and w 2 → w n . For single Compton scattering, N l i l N n in M C i l imin N m im is equivalent to (23), (24) to (32). The corresponding matrices for positron Compton and Breit-Wheeler, M pC and M BW , can be obtained from M C using the same replacement rules as explained for (24) to (32).

C. Spin dependence
The spin/polarization treatment above is needed in order to sum over the intermediate particles. The same treatment can also be used in order to study the dependence on the spin/polarization of initial and final state particles. In this section we consider the dependence of the two-step part of trident on the spin of the initial electron, described by the vector n as in Sec. III. Spin effects in single and double nonlinear Compton scattering have been studied in e.g. [23,34,35]. By either a direct calculation or by simply omitting the average over n in P C P BW we find where the first term gives the spin average, which we calculated in [6], and the second term gives the spin dependence, where P = P ⊥ + P , and where X = 1 2 (w 2 + w 1 ), Y = 1 2 (w 2 − w 1 ), and W 34 = w 3 · iσ 2 · w 4 . Note that only the photon emission part of the integrand depends on the spin. Consider a constant field a µ = δ 1 µ a 0 φ. To obtain the dominant contribution we replace θ(θ 42 )θ(θ 31 ) → θ(σ 43 − σ 21 ). We find where . So the maximum and minimum probability is obtained with spin orthogonal to the field and the propagation direction. For χ 1 we find so the spin dependence is smaller than the average by a factor of χ/27 1 in this regime. We recognize this factor from Eq. (24) in [17].

IV. SADDLE-POINT APPROXIMATION
In this section we will present a method that can be used for a first approximation of the higher-order processes considered in the previous section. In [24] we used a saddle-point approach to obtain an approximation for χ < 1 of the spectrum for Compton scattering with a remarkably good agreement with the exact numerical result, including small and fast oscillations. In this section we use the same method for trident. Saddle-point/semiclassical methods are of course often used for various strong-field processes, see [57,58] for two recent studies of Breit-Wheeler pair production. For the field we consider here we can perform the integrals in Θ ij and ∆ analytically in terms of error functions. We consider linear polarization which is simpler in this approach and also leads to a richer spectrum. The saddle points are determined by these two equations and similar for σ 43 and θ 43 . Note that these are exactly the same saddle-point equations as in [24] for Compton scattering. So, we can reuse the saddle points we already have, and we refer to [24] for the details on how to obtain them. For a field with many oscillations there is a large number of saddle points to include, which all lie in the complex plane. For one-dimensional integrals one can deform the original integration contour to a sum over steepest descent contours that go through some saddle points, but not necessarily all saddle points. For multidimensional integrals like the ones we have here, it is much more nontrivial to construct the higher-dimensional version of steepest descent contours and it is in general a nontrivial question which saddle points one should actually include. (These questions are also considered in Monte-Carlo/Lefschetz thimbles approaches to e.g. the sign problem in QCD [59].) We should of course not include the saddle points that give exponentially large contributions, but by including the other saddle points we have found a good agreement with the exact numerical result. For a pulsed oscillating field we in general have to obtain these saddle points numerically. To do so we need starting points. As explained in [24] we obtain the saddle points by using the corresponding ones for a monochromatic field (which are easier to find) as starting points. The saddle points move continuously through the complex plane as we decrease the pulse length T (or change a 0 ), and in some cases it can be useful to consider a couple of intermediate values of T between the monochromatic case T = ∞ and the actual value of T . As explained in [24], each pair of σ and θ have a set of saddle points which are characterized by two integers, n and m, where increasing n and m correspond, respectively, to increasing Re σ and Re θ (see [24] for the exact definition of n and m). There are two sets of saddle points. Here the dominant contribution comes from the ones that are continuously connected to in the monochromatic limit T → ∞.
A new aspect compared to the first order processes is that now we have step functions for the φ integrals, which lead to restrictions on which saddle points to include. For P 22→2 dir we have θ(σ 43 − σ 21 ) and therefore we should only include saddle points with n 43 ≥ n 21 . For n 43 = n 21 the step function removes one half of a Gaussian integral, which gives an overall factor of 1/2 compared to the cases with n 43 > n 21 . For P 22 dir we have a more complicated step function, θ σ 43 − σ 21 − |θ43−θ21| where the factor in front of the integral in the second line gives the relative factor compared to the case without a step function.
Note that to find the saddle points numerically, we only have to specify a 0 and T , so we have P(s) =α 2 1 + [1 − κ 01 ][1 + κ 23 ] q 2 1 r 1 r 2 n21,m21,n43,m43 "prefactor" where "prefactor" and Θ ij depend on a 0 and T via the numerically obtained saddle points (the Gaussian integrals, obtained from expanding to second order around the saddle points, are performed analytically and then evaluated by inserting the numerical saddle points), but not on the momenta s i or b 0 .
The saddle-point approximation can be expected to be good when χ is sufficiently small. As can be seen by comparing Fig. 3 with the corresponding exact results in [60], this saddle-point approximation captures many of the features of the exact result even for χ = 0.5, which is not particularly small. The most easily seen difference is the slightly higher peak in the saddle-point approximation. For χ = 1 we start to see a bit larger differences also in the shape of the spectrum, but the saddle-point result still gives a good first approximation. Producing the data for the saddle-point approximation in Fig. 3 is of course much faster than obtaining the exact results. So, the saddle-point method can be used to quickly test and find interesting parameter values, for which one can then produce more accurate results with an exact integration, e.g. with the methods described in [60].

V. CONCLUSIONS
For sufficiently large a 0 or sufficiently long pulse, one can expect the dominant contribution to trident to come from the incoherent product of nonlinear Compton scattering and Breit-Wheeler pair production. The nontrivial problem is how to treat the polarization of the intermediate photon. For constant-crossed fields it was already known [12,16,17] that the two-step can be obtained by summing the incoherent product over two suitable polarizations vectors. In [6] we showed that this simple sum can be generalized to inhomogeneous fields for a 0 ∼ 1 (the two-step in LCF/constant-crossed fields corresponds to the leading order term in an expansion in 1/a 0 1). In this paper we have studied how to generalize this to general (e.g. circular) field polarization. This turns out to be nontrivial. We have managed to find such a generalization, and it is not simply a (single) sum over the polarization of the intermediate photon.
This gluing generalization involves a 3D unit (Stokes) vector describing the photon polarization. In [24] we provided a similar gluing generalization to double nonlinear Compton scattering, where the intermediate particle is an electron instead of a photon. In that case the spin of the intermediate electron enters via another 3D unit vector (and the two-step is again not simply a sum over two independent spins). Interestingly, although the 3D unit vectors for the photon polarization and electron spin are different objects and transforms differently, they enter the construction of gluing estimate in basically the same way. Now, if this gluing approach only worked for the second-order processes, trident and double Compton, it would perhaps not have been so useful, because we have calculated all contributions to the probability independently of this gluing approximation and so we anyway know what the two-step is and what the gluing approach has to give. However, in this paper we have generalized to higher orders and showed that the N -step part of N'th-order processes with N > 2 (which can be obtained in a spin/polarization basis independent way with Dirac traces of / p+1 etc.) can be obtained with the same gluing approach. ACKNOWLEDGMENTS 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.