Single, double and higher-order nonlinear Compton scattering

We study single, double and higher-order nonlinear Compton scattering where an electron interacts nonlinearly with a high-intensity laser and emits one, two or more photons. We study, in particular, how double Compton is separated into one-step and two-step parts, where the latter is obtained from an incoherent product of two single-photon emissions. We include all contributions to double Compton and show that the previously neglected exchange term is in general on the same order of magnitude as the other one-step terms. Our approach reveals practically useful similarities between double Compton and the trident process, which allows us to transfer some of our previous results for trident to double Compton scattering. We provide a new gluing approach for obtaining the dominant contribution to higher-order Compton for long laser pulses. Unlike the standard gluing approach, our new approach does not require the intensity parameter $a_0$ to be much larger than one. For `hard' photons we obtain several saddle-point approximations for various field shapes.


I. INTRODUCTION
In [1] we studied the trident process [2][3][4][5][6][7][8], e − → 2e − + e + , in plane-wave background fields, and derived compact expressions for the probability for arbitrary background field shapes. Here we will apply the same methods to another second-order process, namely double nonlinear Compton scattering [9][10][11][12][13][14], where the incoming electron emits two photons, e − → e − + 2γ. This is also a process that one can separate into one-step and two-step parts, where the latter is obtained by incoherently gluing together the probabilities of two single-photon emissions. The two-step is expected to be a good approximation of the total probability for sufficiently high intensities, or more precisely for a 0 = eE/(mω) 1, where E is the field strength and ω a typical/characteristic frequency of the (in general pulsed) background field. This two-step dominance is what makes it possible to use particle-incell (PIC) simulations to study complicated higher-order processes in high-intensity fields [15]. This regime is also associated with the locally-constant-field (LCF) approximation, which entails further simplifications. In this paper we are interested in corrections to this two-step approximation. In particular, the one-step can be separated into (what we call) direct 1 and exchange terms, where the latter comes from the cross-term between the two terms in the amplitude which are related by exchanging the two emitted photons. A similar exchange term appears in the trident case, and in [1] we showed that, while previously neglected, it is in general on the same order of magnitude as the direct part of the one-step. Here we make a * dinu@barutu.fizica.unibuc.ro † greger.torgrimsson@uni-jena.de 1 Note that we do not use "direct" as synonymous to the one-step term. By "direct" we mean instead the non-exchange part. The two-step only has a direct part while the one-step has both direct and exchange parts. similar investigation into the importance of the exchange term in double Compton. That the exchange term can be important e.g. for a 0 ∼ 1 was also found in [12]. For a 0 ∼ 1 the one-step is in general on the same order of magnitude as the two-step. However, if the field is sufficiently long then the probability is again dominated by a term that can be expressed as an incoherent product of two single-photon emissions. If a 0 is not large one should of course not expect this two-step to be the same as the LCF two-step. While spin effects are usually neglected in PIC simulations, to obtain the complete two-step in the LCF regime one has to sum the incoherent product over the spin of the intermediate electron [9,14]. In this paper we identify a term in double Compton that dominates for sufficiently long pulses without assuming a 0 1 or any particular field shape, and then we show that this two-step can be obtained from an appropriate sum of the incoherent product of two single-photon emissions. We do this for an arbitrary background field. For fields with linear polarization one can obtain the twostep by summing over spin in essentially the same way as in the LCF regime. However, for fields that do not have linear polarization things become more nontrivial, because in general one has to take into account the fact there is a spin sum already on the amplitude level, which in general leads to a double spin sum on the probability level. We have found a simple prescription for obtaining the entire two-step from the spin-dependent probability for single Compton. This gluing approach is to the best of our knowledge new and seems promising for studying higher-order processes. We have checked that it gives the correct results for triple and quadruple Compton, where the electron emits three and four photons.
Calculating higher-order processes means performing higher-dimensional integrals. Numerical integration can quickly become challenging. In our approach we integrate analytically over the transverse components of the momenta, and then the longitudinal momentum spectrum is obtained by performing a number of lightfront-time (x + ) integrals. The exponential part of these integrands can in general be expressed in terms of an (x + -dependent) effective mass, and the integrals can be performed with the saddle-point method. In fact, the integrals for double Compton are very similar to the ones in the trident case [1], so we have for example been able to reuse saddle points we found in [1] for double Compton, and the new saddle-point results we provide here can also be translated to the trident case. For certain simple field shapes we can obtain simple analytical approximations, but the saddle-point method can also be useful even if one has to find the saddle points numerically, as it can offer a quick estimate and a check of exact numerical integrations. We show here that the saddle-point method can give a good approximation of even quite small and fast oscillations in the spectrum.
This paper is organized as follows. We focus first on double Compton. In Sec. II we give the necessary definitions. In Sec. III we provide compact expressions for the exact probability for arbitrary field shapes. In Sec. IV we separate the probability into one-step and two-step terms and compare with the incoherent product of two single-photon emissions. This comparison helped us to find a new gluing approach, which we in Sec. V confirm for triple and quadruple Compton scattering. In Sec. VI we derive simple analytical approximations for "hard" photons for various field shapes. In Sec. VII we apply the saddle point method to fields with many oscillations and hence many contributing saddle points, which lead to interference effects in the momentum spectrum. We consider single Compton and compare this saddle point approximation with an exact numerical integration and find very good agreement. In Sec. VIII we consider double Compton in the LCF approximation. We show, in particular, that the exchange term can continue to be on the same order as the direct part of the one-step also for larger χ.

II. DEFINITIONS
We use the same formalism and notation as in [1], which we briefly recall here for convenience. Lightfront coordinates are defined by µ is a light-like wave vector and a ⊥ (kx) a polarization vector with an arbitrary dependence on lightfront time x + . We use units with c = = 1 as well as m e = 1, and absorb the electron charge into the background field, i.e. ea µ → a µ .
We focus first on double Compton, where the final state contains one electron with p µ and σ and two photons with momenta and polarization vectors l µ 1 , l µ 2 and µ 1 , µ 2 . We use lightfront gauge, so in addition to l (l) = 0 we also have k = 0. The amplitude for two-photon emission, M , is defined via the evolution operator U by (3) As in [1], in order to reduce the number of parameters on which the probability depends, we integrate analytically the probability over the Gaussian transverse momentum integrals [16] and sum/average over spins and polarizations, where the factor of 1/4 is due to spin-averaging and the presence of identical particles, andp =p −l 1 −l 2 . We separate the amplitude into M = M 12 + M 21 , where M 21 is obtained from M 12 by replacing l 1 ↔ l 2 and 1 ↔ 2 , which on the probability level gives |M | 2 = |M 12 | 2 + |M 21 | 2 + 2ReM 21 M 12 . We refer to the first two terms as the direct part and the cross-term as the exchange part, i.e.
where (1 ↔ 2) is obtained from the first term by replacing l 1 ↔ l 2 and 1 ↔ 2 , and We have relegated the calculation of the amplitude to the appendix as it only involves standard methods. The important thing to note is that the amplitude contains two terms, M 12 = M 12 1 + M 12 2 , where M 12 1 has one x + integral and M 12 2 has two. These terms are illustrated in Fig. 1. As in [1], this leads to a separation of the direct and the exchange part of the probability into three terms with different number of x + integrals, We integrate over the transverse components of the photon momenta, and define a longitudinal momentum spectrum P(q) as where q i = kl i /kp and s 3 = kp /kp = 1 − q 1 − q 2 . We also define b 0 = kp, s 1 = 1 − q 1 , s 2 = 1 − q 2 and s 0 = 1. We have chosen the definitions of s i in this way to highlight the symmetry with the expressions in [1] for the trident case.

III. EXACT ANALYTICAL RESULTS
The different contributions are illustrated in Fig. 2. For the direct part of the simplest term we find ij , and M is an effective mass given by [17] where the lightfront-time average is and where the Lorentz momentum is given by The exchange part P 11 ex ( 1 , 2 ) depends nontrivially on the polarization vectors, but after summing over polarization vectors we find P 11 ex (q) = 0, in contrast to the trident case [1] where the corresponding term is nonzero. For the terms with three x + integrals we find and where D 12 = ∆ 12 ·∆ 32 and The i factors initially make the transverse momentum integrals converge and at this stage provide a prescription for how to avoid the singularities in the φ integrals. This is equivalent to a shift in the φ-integration contours. From now on we leave these i factors implicit, they can be reinstated by replacing φ 1,3 → φ 1,3 − i /2 and φ 2,4 → φ 2,4 + i /2. For the direct term with four x + integrals we find where κ ij = (s i /s j ) + (s j /s i ), and where w 1 = ∆ 12 , w 2 = ∆ 21 , w 3 = ∆ 34 and w 4 = ∆ 43 . For linear polarization we have W ijkl = 0. In contrast to the trident case, here we have a dot product between the two steps even for linear polarization. Finally, the last term is given by These expressions for P 22 ex look remarkably similar to the corresponding ones in the trident case [1]. In fact, one can show that (20) can be obtained from Eq. . This means that P 22 ex in double Compton scattering has the same symmetries as in the trident case and can be calculated in a similar way.

IV. TWO-STEP AND ONE-STEP
In this section we compare (17) with the product of two single-photon emissions. To treat the electron spin we use the following representation of the Dirac matrices and the following spinor basis (cf. [18]) (29) This spinor basis is particularly convenient for the quantities that we calculate here. An arbitrary spinor can be expressed as a linear combination of these, Instead of ρ and λ we express the spin dependence in terms of the components of the unit vector n that points in the average spin direction for p = 0, i.e.
Now, the probability of single-photon emission, summed over photon polarization and transverse momenta, is given by where n 0 and n 1 are the spin vectors of the initial and final electron, respectively. The first term P gives the probability averaged 3 over initial and final spins, The remaining terms give the spin dependence, and wherek = {0, 0, 1},k X·V =k(X·V) etc., D 1 = ∆ 12 ·∆ 21 and where the Pauli matrix is given as usual by Note that n 1 gives the average spin direction for p 1 = 0 and we have integrated over p 1⊥ with n 1 fixed. Regardless of whether or not this is the most directly relevant quantity for spin-sensitive experiments, we show below that (32) can be very useful for studying multi-photon emission. For a detailed investigation of spin effects in nonlinear Compton scattering see [19].
In evaluating these expressions we can put s 0 = 1. One reason for keeping s 0 explicit is that it helps us to glue together two single-photon emissions, which one might expect to be done according to where one factor of 1/2 comes from averaging over the spin of the initial electron and another factor of 1/2 comes from the symmetrization. We can write this as where the factor of 2 2 is due to the replacement of the sum of two spins with their average for n 1 and n 3 , and we have omitted the arguments of the probability terms (the second factor in each term is obtained by making the appropriate replacements in (33), (34), (35) and (36)). It is easy to show that the P P term gives the QQ-term in (17). The remaining terms are more subtle. We first note that these terms can be expressed as where V 1 and X 1 are given by (37), and V 2 and X 2 are obtained by replacing (37). This should be compared with the corresponding term in P 1 · n 1 n 1 ·P 0 , i.e.
The gluing approach works if (42) gives (41) after summing over n 1 . In (40) we have only used 1 = 1 and n = 0. For linear polarization with a ∝ e 1 we have X·V = 0, and then we can simply sum over n 1 = ±e 2 . For arbitrary polarization we cannot in general obtain (41) from (42) unless we let n 1 depend on both φ 1 and φ 2 (or φ 3 and φ 4 ). However, for arbitrary polarization in the LCF regime we have where σ ij = (φ i + φ j )/2, so then we can neglect the X terms and obtain (41) by choosing the spin direction to be perpendicular to the locally constant field andk, i.e. either n 1 = ±k ×â(σ 21 ) or n 1 = ±k ×â(σ 43 ), wherê a = a/|a|. The reason that the naive gluing approach does not always work is because we actually have a sum over the spin of the intermediate electron already on the amplitude level, so, instead of having on the probability level just one sum over n 1 , one should have one sum for the amplitude and a second sum for its complex conjugate, where the sum is over ±n (or ρ and ρ + π) for some n. While the momentum p 1 is the same in the amplitude and its complex conjugate, the spins n 1 and n 1 need not be the same. Let Compared to (40), one can show that P same = 2( P P + P 1 ·n 1 n 1 ·P 0 ) and These clearly depend on the spin directions ±n 1 one chooses to sum over, but their sum is independent of n 1 , As we saw above, for linear polarization or in the LCF regime we can choose n 1 such that P diff vanishes, but in general we need to include this term. Fortunately, our results suggests a simple cure of the naive gluing approach: Include factors of 2 in the overall prefactor as if we only had one sum over n 1 as above, and then simplify using 1 = 1, n = 0 and importantly nn = 1, where the last ingredient is motivated by the contribution from n 1 = n 1 . We show in the next section that this simple procedure also works for triple and quadruple nonlinear Compton scattering. Note that this improved gluing procedure gives us the dominant term for sufficiently long pulses, for any polarization and field shape, and we can in particular go beyond the usual LCF regime (where gluing first order, albeit spin-averaged, processes is a basic component of PIC codes for a 0 1) and consider a 0 ∼ 1.
In the gluing approach one also has to make sure that the second step happens after the first, which can be done by including a step function θ(σ 43 −σ 21 ). In (17) we have two step functions, which we deal with in the same way as in [1], i.e. we write P 22 where P 22→2 dir and P 22→1 dir are obtained, respectively, from the first and second term in It is P two := P 22→2 dir (rather than P 22 dir ) which we refer to as the two-step. Although it can be obtained from the above gluing approach, we can obtain it without reference to the gluing approach by selecting one part of the exact/total probability. This part scales quadratically in the volume/pulse length and dominates for sufficiently long pulses.

A. Triple Compton
In this section we calculate the three-step part of triple nonlinear Compton scattering, i.e. the part of the probability of three-photon emission that dominates for long pulses, illustrated in Fig. 3. The emission of three photons by an electron colliding with a single photon has been studied in [20], but to the best of our knowledge nonlinear triple Compton has not been studied in the regime we are interested in here. This is in principle a straightforward generalization of our results for the twostep part of double Compton, except that it takes more time to simplify the prefactor. After some simplification we find where (37), and "permutation" is an instruction to sum over all permutations of the emitted photons. Notice that in this section we have chosen the definition of s i such as to have a simple generalization for higher orders, which is different from the previous section. Note also that the exponential part is a simple generalization from single and double Compton. Compare this with the result of the gluing approach described in the previous section, which in this case gives [ 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 ] + permutations = 4 3 P P P + P P 1 · n 2 n 2 ·P 0 + P 1 · n 1 n 1 ·P 0 P + P 1 · n 1 n 1 ·P 01 · n 2 n 2 ·P 0 where the arguments are again suppressed. The factor of 2 3 comes from the (initial) assumption that we are summing over two spin states for n 1 , n 2 and n 3 , and for linear polarization a ∝ e 1 we can obtain (50) by summing over n 1 = ±e 2 and n 2 = ±e 2 . For arbitrary polarization we can obtain (50) from the following procedure: We write an overall factor of 2 N /N ! and replace all sums with ... , and then we simplify with 1 = 1, n = 0 and nn = 1. Note again that it is the replacement nn = 1 that allows us to obtain all terms in the general case. In analogy to (49), we define the three-step P three by replacing the product of step functions in (50) according to

B. Quadruple Compton
We have also checked that the above gluing procedure gives the correct result for quadruple nonlinear Compton scattering, i.e. the emission of four photons, which is illustrated in Fig. 4. Here our gluing approach is not only useful for interpreting the expressions, it is also very useful for simplifying the complicated prefactor. The expressions obtained by a direct calculation without using the gluing procedure are very long and complicated, but we have checked that the results can be expressed neatly  [ 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 ] [ P + n 3 ·P 0 + P 1 ·n 4 + n 3 ·P 01 ·n 4 ] where 1 = 1, n = 0 and nn = 1 for each n i . Even with the help of an advanced symbolic-calculation program such as Mathematica, obtaining or confirming this result by a direct calculation takes a long time. Instead of calculating the prefactor from the trace of a long expression, we replaced all factors / p i + 1 (which would appear in the trace) by sums of uū expressed with a particular spinor representation. This leads to a multi-dimensional sum over all spins and polarizations. We then integrated and simplified each term in this sum before adding them together. The sum was parallelized on several Mathematica Kernels to speed up the calculation. We mention this to stress that our gluing approach is very useful for quickly obtaining these higher-order results. Note again that we only obtain all terms by replacing nn = 1 to account for the terms that would be missing if one replaces the double sums over the spins of the intermediate states, i.e. n 1 , n 2 and n 3 in this case, with single sums as explained above for double Compton.
Although we have not yet proved that this gluing procedure works at arbitrarily high orders, the fact that it does work for double, triple and quadruple Compton suggests that we have a method for obtaining the exact Nstep part for N -Compton for arbitrary N , where the Nstep dominates for sufficiently long pulses. We plan to further study this gluing approach and to generalize it to other higher-order processes involving more than one fermion, like the trident process.

VI. SADDLE POINT APPROXIMATION
In this section we obtain saddle point approximations, which help us to understand the structure and relative importance of the various terms. We can expect these approximations to be good for χ 1 as long as q 1 and q 2 are not too small, so we are in particular outside the infrared region and do not have to worry about IR divergences. We also have to assume that a 0 is not too small. The calculations are very similar to the ones in [1], except that this time, in order to avoid IR divergences, we do not integrate over the longitudinal momenta. We consider linearly polarized fields, a(φ) = a 0 f (φ). In this section we focus on the dominant contribution from a single saddle point located around a single field maximum.

A. Locally constant fields
We consider first the LCF regime where we can expand the probability in 1/a 0 1. For the one-step terms we find (57) and and for the two-step term we find For a constant field for which χ(φ) is zero outside an interval of length ∆φ, we simply have dφ → ∆φ and dσ 1 dσ 2 θ(σ 2 − σ 1 ) → ∆φ 2 /2. For a pulsed field we can also perform the remaining φ integral with the saddle point approximation. Let us for simplicity assume one dominant field maximum with f (0) = 1, f (0) = 0 and f (3) (0) = −ζ, then the results are obtained from the constant field results by replacing for the one-step terms, so for example and ∆φ 2 2 (62) for the two-step term, which simplifies to We see a few things that are similar to the trident case: All terms have the same exponential, and P 11 and P 12 are smaller than P 22→1 by a factor of χ. We also see that the exchange terms are on the same order of magnitude as the direct part of the one-step. In fact, here P 22 ex cancels P 22→1 to leading order, so the χ expansion of the prefactor of P one starts at one order higher than the leading order of the direct part of P one . This also means that P 11 and P 12 contribute to the first nonzero order, in contrast to the trident case. Thus, the exchange term is even more important for double Compton scattering.
In the trident case we could compare the direct terms with previous LCF results. For double Compton, on the other hand, we are not aware of any previous approximations with which we could compare our results. The χ < 1 approximation in e.g. [9] is for the probability integrated over the photon momenta, which has a different form because of the contribution from softer photons.
The exponential part of the above terms can be written where N = 2. Assuming again one dominant field maximum, for triple Compton it follows from (50) and r ij + r jk = r ik that P three ∼ (64) with N = 3. Similarly, for quadruple Compton we find P four ∼ (64) with N = 4. This suggests a simple generalization to the emission of an arbitrary number of photons.

B. Sauter pulse
In the previous section we considered a 0 1 which allows us to consider an arbitrary pulse shape. Here we will consider a particular pulse shape, namely a Sauter pulse a(φ) = a 0 tanh φ, which allows us to obtain explicit analytical expressions also for a 0 1, i.e. to go beyond the LCF approximation. The calculation is very similar to the corresponding one in [1] for the trident case. In particular, we have a saddle point at the same values of the φ i variables as in [1], independently of q i . For the "two-step" we find For a 0 1 we recover (63) to leading order. For the "one-step" terms we find and while P 11 dir and P 12 are again smaller than the above terms by a factor of χ. Notice that these expressions are very similar to the ones in [1] for trident: the dependence on a 0 in the exponent is exactly the same as in [1], and the relation between P 22→1 and P 22→2 is also exactly the same. We also find that the (leading order) exchange term P 22 ex is on the same order of magnitude as the (leading order) direct terms P 22→1 and P 22→2 . Here, though, P 22 ex is not only on the same order of magnitude, but it in fact cancels P 22→1 to leading order in χ; this generalizes the a 0 1 results in the previous section to a 0 1. Note also that the dependence on the momenta remains the same as in the a 0 1 limit.

C. Monochromatic field
For a monochromatic field we can again find saddle point approximations for general a 0 1. For this field there are many saddle points that contribute. We begin in this section with the simplest ones, which are the same as those we studied [1] for the integrated trident probability, (68) These already give a good approximation to the locallyaveraged spectrum. In the next section we include additional saddle points that give oscillations to the spectrum. For the two-step we have saddle points both for n 1 = n 2 and n 1 < n 2 , where the two photons are emitted at the same and different field maxima, respectively. For the contribution from one saddle point with n 1 = n 2 we find We have again the same function of a 0 in the exponent as in the trident case [1]. For a 0 1 we recover the LCF approximation (63) from (69). The contributions from one saddle point (with n 1 = n 2 ) to the dominant one-step terms are given by (70) and The relation (70) is exactly the same as in the trident case [1], and, as for the LCF and Sauter cases, we find that the exchange term cancels the direct part of the onestep to leading order. The other one-step terms, P 11 dir , P 12 dir and P 12 ex , are again smaller by a factor of χ 1, but have to be included if one is interested in the first nontrivial order of the total one-step, since P 22 ex cancels P 22→1 to leading order. These expressions give the contribution from one field maximum with the shape of a sinusoidal field, and for a 0 1 they are one the same order of magnitude. If we have a sinusoidal field with several equivalent field maxima, then the two-step dominates because it also receives contributions from n 1 < n 2 and not only n 1 = n 2 , which means that it scales quadratically in the number of oscillations compared to the linear scaling of the one-step terms. In contrast to the trident case, here the contributions from n 2 = n 1 +2n−1 are different from the ones from n 2 = n 1 + 2n, where and This difference has the same a 0 dependence but a different dependence on the momenta in the prefactor. This difference is due to the (w 2 − w 1 )·(w 4 − w 3 ) term in (17). For a 0 1 we recover the LCF results. For a 0 1 the exponent goes as which is the expected perturbative scaling: Momentum conservation at O(a N 0 ), Thus, the exponent scales as a 2N0 0 , where N 0 is the minimum number of photons from the background field that need to be absorbed in order to emit two photons with longitudinal momenta q 1 and q 2 .
For a Sauter pulse the exponent scales as Since the Sauter pulse has a wide Fourier transform with only exponential decay (which is slow in this context), this scaling agrees with the absorption of a single photon from the background field with (Fourier) frequency N 0 k 0 (cf. [1,21]).

D. General anti-symmetric potential
Both the Sauter pulse and the sinusoidal field considered in the previous two sections fall in the class of fields that have anti-symmetric potentials, a(−φ) = −a(φ). In this section we derive the probability for such fields, assuming for simplicity one dominant field maximum and linear polarization but without choosing a specific field shape. Let a(φ) = a 0 f (φ). We have a saddle point at where f −1 is the inverse of f , and, as before, φ = ϕ = η = 0. We can still perform the integrals with the saddle where and for the one-step terms we find and In deriving these expressions we have assumed that f (iz) > 0 and 0 < 1 − 1 a0zf (iz) < 1, which we will justify below. Note that P 22 ex cancels P 22→1 to leading order independently on the field shape.
To make these expressions more explicit, we consider the class of fields defined implicitly via [22] where each b 4 characterizes a different field shape, see Fig 5. For example, b = 1/2 and b = 1 give us the sinusoidal field (or rather one peak of it) and the Sauter pulse, respectively. For general b the field f (φ) is given implicitly in terms of a hypergeometric function by 4 Not to be confused with b 0 .
For this class of fields we find simple explicit expressions for the probability, using for the exponent and for the prefactor It is easy to check that 0 < 1 − 1 a0zf (iz) < 1 for general a 0 and b. Now everything is explicitly expressed in terms of a 0 and b, which in turn only enter in the arguments of 2 F 1 . For b = 1/2 and b = 1 we recover the results in the previous two sections for a monochromatic field and a Sauter pulse, and for arbitrary b we recover for a 0 1 the LCF results above by expanding in 1/a 0 and using the relation b = ζ/2. The hypergeometric functions also simplify more generally for b = j/2 where j is an integer. For example, for b = 3/2, which corresponds to f (φ) = (1 + φ 2 ) −3/2 , we find a particularly simple prefactor and while for b = 5/2 we find a simple exponent and The prefactors above have been derived under the assumption that a 0 is not too small. The exponents, on the other hand, have the expected perturbative limit for a 0 1: For b > 1/2 the exponent becomes independent of the field strength, In the perturbative regime the minimum energy that needs to be absorbed is N 0 ω, where N 0 is given by (76). For a monochromatic field, N 0 photons have to be absorbed. For b > 1/2, on the other hand, the Fourier transform a(ω f ) has a slow, exponential decay, which (since |a(ω)| 2N0 /|a(N 0 ω)| 2 ∼ a 2(N −1) 0 1) means that the process occurs already at first order, with the absorption of a single photon with ω f = N 0 ω. At ω f ω, the exponential behavior of the Fourier transform is governed by the singularity φ s closest to the real axis, i.e. a(ω f ) ∼ e −|ω f φs/ω| . We find from the |f | → ∞ limit of (84) a singularity at At ω f = N 0 ω this implies |a(ω f )| 2 ∼ (91), so (91) agrees with what one can expect to find in the perturbative limit.

E. Single Compton scattering
While the results in the previous section are for double Compton, it should be clear that the same method can be used to derive similar expressions for other plane-wave processes, like nonlinear Breit-Wheeler or trident pair production. In this subsection we simply give the corresponding result for single Compton. The saddle point approximation is obtained e.g. from (33) in the same way as for the above expressions for double Compton, and we find where s 1 now corresponds to the final electron. For the class of fields defined by (83) we can again obtain explicit expressions using (85) and (86).

VII. SADDLE POINT APPROXIMATION FOR INTERFERENCE EFFECTS
In this section we study fields with many oscillations and with several saddle points that lead to oscillations in the spectrum. We choose the following field Since the exponential part of the integrand for the Nstep part of the N -photon emission probability is a simple generalization of the N = 1 case, we focus here on single Compton. See [23][24][25] for other semi-classical approximations. The saddle points for (33) are determined by where, again, θ ij = φ i − φ j and σ ij = (φ i + φ j )/2. Note that these equations only depend on the field parameters, a 0 and T in our case, but not on the momenta b 0 or q i . To obtain the saddle points for finite T , we first find the saddle points for a monochromatic field (T = ∞) and then use them as starting points for a numerical rootfinding of the corresponding saddle points for finite T . Depending on how large/small T is, one may find it useful to obtain the saddle points by first considering a sequence of T values between T = ∞ and the desired value, and/or by starting with a simple a 0 value and gradually change to a more difficult one, cf. the numerical continuation in [26]. The saddle point equations can be expressed in terms of the "prefactor functions" ∆ (16) as and which imply that all saddle points, for any field shape, are determined by ∆ = ±i. The saddle points are therefore necessarily complex. For the monochromatic field we find saddle points at where n, m = 0, ±1, ±2, . . . . We also have saddle points at where η m can be found numerically by using 2iarcsinh 1 a0 +(2m−1)π as starting points. In Fig. 6 we show saddle points for a pulsed field, which are obtained numerically with the ones in (98) as starting points. For the first set of saddle points (98) we find and for the second set (99) Note that these values of ∆ do not change as we decrease the pulse length from T = ∞ to a finite T . Let now δσ = σ − σ saddle and δθ = θ − θ saddle . The quadratic fluctuation of Θ around any point can be expressed in terms of ∆ and the derivative of the field f , but at the saddle points we can simplify using ∆ = ±i. To leading order we can put δσ, δθ → 0 in the pre-exponential part of the integrand. Having expanded Θ to second order in δσ and δθ, we now have simple Gaussian integrals for each n and m which we perform analytically, i.e. we have for each saddle point where the coefficients c i are in general complex and obtained by finding the saddle points numerically. For a monochromatic field we find with (98) an exponential part given by (cf. (69)) From this we see that the saddle points with m = 0 lead to oscillations in the spectrum around the m = 0 result studied in the previous section. We also see that the frequencies of these oscillations increase with decreasing χ or increasing a 0 . Since this saddle point approximation is good for small χ, these oscillations can be relatively fast and hence contribute less after integrating over the momenta. In Fig. 7, 8 and 9 we compare this approximation with the results obtained by an exact numerical integration. How many saddle points one needs to include depends of course on several parameters. To obtain these results we have summed over the saddle points with |n| ≤ 40 and |m| ≤ 20. These plots show that the saddle point approximation is remarkably good. It can in fact be difficult to see that there are actually two different curves in the large q 1 part. Note that at a 0 = 1 the LCF ap- proximation is not good, not even for an average where the oscillations are neglected. Our non-LCF saddle-point approximation, on the other hand, gives a very good approximation of even the nontrivial oscillations. From these plots we see that the oscillations in the spectrum become smaller and faster as a 0 increases. Fig. 8 shows that already at a 0 = 2 the oscillations are quite small on a log scale. However, by zooming in one can see that our approximation is capable of correctly describing even very fine details in the spectrum. In these figures we also plot the saddle-point approximation obtained by only including the m = 0 saddles from (98). This give a good approximation of a locally averaged spectrum. While the LCF approximation becomes more accurate for increasing a 0 , for a 0 = 2 our approximation, even just the simpler one, is still much better. In Fig. 9 we see that for a 0 = 4 the oscillations are so small that it might be difficult to see them without zooming in, and in this case  the LCF approximation is quite good. Although there are no IR divergences in single Compton for this field shape [27,28], the probability can become larger than one even for some of these nonextreme parameter values. That this can happen is well known [16,29].

VIII. DOUBLE COMPTON LCF
We now return to double Compton in the LCF regime. In the previous sections we showed that for χ 1 the exchange term is on the same order of magnitude and even cancel the direct part of the one-step to leading order. Here we study what happens at larger χ. We need to keep q i sufficiently large as it is known that the LCF approximation is not good for softer photons [30]. In Fig. 10 and 11 we show the one-step as a function of χ for different values of q i . What is actually shown in these figures is the corresponding "rate" R defined by where φ = (φ 1 + φ 2 + φ 3 + φ 4 )/4, in which φ 4 = φ 2 for P 12 and φ 4 = φ 2 , φ 3 = φ 1 for P 11 . Fig. 10 shows that for q 1 = q 2 the direct and exchange terms of the one-step continue to be close to each other also for large χ. The fact that P 22→1 dir and P 22 ex almost cancel each other means that the other one-step contributions, P 11 dir , P 12 dir and P 12 ex , are more important than in the trident case [1]. So, even though P 22→1 dir and P 22 ex are much larger than P 11 dir , P 12 dir and P 12 ex , the size of the total one-step is closer to the latter rather than the former. At least for the parameter values we consider here we find that, in contrast to the trident case, here P 22→1 dir does not change sign, but the total one-step does change from negative to positive as χ increases, similar to the trident case.
We have made a comparison between our numerical results and our saddle-point approximation similar to the one in Appendix C in [1] for the trident case. For sufficiently small χ we again find that each of the first couple of orders give a better agreement. However, here we find that the coefficients in the series in χ increases quite fast. For example, at q 1 = q 2 = 1/3 we find (1 − 11.0χ + 130.5χ 2 − 1847.5χ 3 + . . . ) .
Given that the saddle point approximation can lead to asymptotic series, this growth of the coefficients should not be too surprising, but it does mean that the higher orders are less useful than in the trident case. They only provide an improvement for quite small χ, but there the probability is very small because of the exponential suppression. This is a bit unfortunate if one wants an approximation for the total one-step, because one needs at least the next-to-leading order of P 22→1 dir and P 22 ex since they cancel each other to leading order.
On the other hand, this cancellation also means that neglecting the total one-step compared to the two-step should be a better approximation 5 here than in the trident case. The two-step is shown in Fig. 12 in terms of the following "rate" P two (q) =: dσ 43 dσ 21 b 2 0 R 2 (χ(σ 43 ), χ(σ 21 ), q) . (106)

IX. CONCLUSIONS
In this paper we have studied double nonlinear Compton scattering. By using the same approach as in our previous paper on trident pair production [1], we have showed that many of the results are very similar, which allows us to use the same methods. We have focused on the emission of "hard" photons which makes things more similar to the trident process than if we had included soft photons, we can for example obtain saddle point approximations for χ < 1 that are similar to the ones we obtained in [1]. Focusing on hard photons is also motivated by the fact that they can be more interesting/useful e.g. for subsequent pair production. The saddle point method has not only allowed us to find simple analytical expressions for simple field shapes, we have also considered a more nontrivial, pulsed oscillating field. We then have to obtain the saddle points numerically, but by comparing with the exact numerical result for single Compton scattering we find a very good agreement, even for small and fast oscillations in the spectrum. Since the saddle point approach is much faster it can therefore be a useful method for studying this as well as similar processes. Indeed, since the exponential part of the integrand is very similar for double and higher-order Compton, one can also apply this method to those processes, for which an exact numerical integration would take a long or too long time. We have also made preliminary calculations for trident and found that the same saddle point method can also be used to study oscillations in the momentum spectrum there.
The two-step part of the probability is related to two one-photon emissions. By studying this relation in detail for arbitrary polarization we have discovered a new gluing approach, i.e. a method for obtaining the dominant part of sufficiently long laser pulses. Gluing (spin-averaged) LCF probabilities is an important part of PIC simulations, where using LCF results is motivated by considering a 0 1. Our new gluing approach takes the spin of  the intermediate electron into account and gives the dominant contribution for arbitrary field polarization and for a 0 1. So, this goes beyond the usual gluing approach. We have checked that our approach gives the correct results for triple and quadruple Compton scattering. To the best of our knowledge, these processes have not been studied in this regime before. In this paper we have only presented this gluing approach for intermediate electrons.
Our preliminary results for trident suggest that we will be able to generalize our gluing approach to processes with intermediate photons.

ACKNOWLEDGMENTS
We thank Tom Blackburn, Anthony Hartin and Anton Ilderton for useful discussions about the field shape where P on + = (m 2 + P 2 ⊥ )/(4P − ). The (lightfront) spatial coordinate integrals in (A11) give delta functions imply-ingP =p −l 1 =p +l 2 , which means P − > 0. Upon performing the P + integral, the two terms in (A13) give terms with δ(x + 2 − x + 1 ) and θ(x + 2 − x + 1 ), respectively. We find that the term with δ(x + 2 − x + 1 ) is exactly equal to the term ((A10)) that comes from the instantaneous part of the lightfront Hamiltonian, and the term with θ(x + 2 −x + 1 ) is exactly equal to the term ((A9)) that comes from two vertices with the non-instantaneous part of Hamiltonian.