Benchmarking semiclassical approaches to strong-field QED: nonlinear Compton scattering in intense laser pulses

The recoil associated with photon emission is key to the dynamics of ultrarelativistic electrons in strong electromagnetic fields, as are found in high-intensity laser-matter interactions and astrophysical environments such as neutron star magnetospheres. When the energy of the photon becomes comparable to that of the electron, it is necessary to use quantum electrodynamics (QED) to describe the dynamics accurately. However, computing the appropriate scattering matrix element from strong-field QED is not generally possible due to multiparticle effects and the complex structure of the electromagnetic fields. Therefore these interactions are treated semiclassically, coupling probabilistic emission events to classical electrodynamics using rates calculated in the locally constant field approximation. Here we provide comprehensive benchmarking of this approach against the exact QED calculation for nonlinear Compton scattering of electrons in an intense laser pulse. We find agreement at the percentage level between the photon spectra, as well as between the models' predictions of absorption from the background field, for normalized amplitudes $a_0>5$. We discuss possible routes towards improved numerical methods and the implications of our results for the study of QED cascades.


I. INTRODUCTION
Petawatt and multipetawatt laser facilities that reach focussed intensities in excess of $10 22 W cm À2 (Refs. 1 and 2) hold great promise for the study of the interaction of charged particles with electromagnetic fields of unprecedented strength. [3][4][5][6] In these environments, the recoil associated with emission of radiation, known as radiation reaction, can become so large that it dominates the particle dynamics. 7,8 Furthermore, when the energy of individual photons of this radiation becomes comparable to that of the emitting particle, it becomes essential to incorporate quantum effects on this radiation reaction 9 in our modelling of plasmas as sources of high-energy photons, 10 electron-positron pairs, 11,12 or as laboratory analogues of high-field astrophysical environments. 13,14 However, it is not currently possible to use the most general and accurate approach, the theory of strong-field quantum electrodynamics (QED), to model many scenarios of interest, for reasons we will shortly outline. Instead, a semiclassical approach has been widely adopted for use in numerical simulations of laser-plasma and laser-particlebeam interactions. Inherent to this model are a number of assumptions, making it essential that we benchmark its predictions against those from QED. In this work, we do so for photon emission in the collision of ultrarelativistic electrons with intense laser pulses. We focus on the classically nonlinear, moderately quantum regime, motivated not only by progress in the development of the next generation of highintensity laser facilities, 15,16 but also by recent experimental work on radiation reaction in strong fields. [17][18][19] Strong-field QED is not used directly to model these kinds of interactions for a number of reasons. First, a scattering-matrix calculation connects asymptotic free states, thereby requiring ab initio complete knowledge of the spatiotemporal structure of the background electromagnetic field; exact analytical calculations have thus far proven to be possible only for certain field configurations that possess high symmetry, 20 such as plane electromagnetic waves 21 and static magnetic fields. 22 Furthermore, it is generally assumed that back-reaction effects may be neglected. This is especially important when considering QED cascades, 11,23,24 in which the initial state contains a single electron, positron, or photon and the final state many of these, because we expect significant absorption of energy from the background. 25 (See Ref. 26 for analysis beyond this approximation.) Even in the absence of significant depletion, the multiplicity alone makes the calculation of a cascade within strong-field QED extremely challenging. State-of-the-art results are those in which the final state contains two additional particles, such as double Compton scattering [27][28][29] and trident pair creation [30][31][32][33] from single electrons.
These conditions, namely, complex field structure, strong depletion due to back-reaction, and high multiplicity, are ubiquitous in the interaction of high-intensity lasers with particle beams or plasma targets. As such, considerable effort has been devoted to the development of numerical schemes that can model QED cascades as well as selfconsistent plasma dynamics. 34, 35 We characterize these as semiclassical because they factorize the cascade into a product of first-order processes that occur in vanishingly small regions linked by classically determined trajectories; the production rates and spectra are calculated for the equivalent QED process in constant, crossed fields. 36 This is possible because at high intensity (to be defined in Sec. II), the formation length over which a photon is emitted, or an electronpositron pair is created, is much smaller than the length scale over which the field varies. 21 The approximation that emission occurs instantaneously, and therefore that the rate for a constant, crossed field may be employed, is called the locally constant field approximation (LCFA). Monte Carlo implementations of QED processes based on this have found widespread adoption in particle-in-cell (PIC) codes (see Ref. 37 for details). Depletion in these codes is therefore treated classically, as QED processes alter the plasma current density j, which alters the energy density of the self-consistent electric and magnetic fields E and B via the j Á E term in Poynting's theorem; see Refs. [38][39][40] for examples of how this drives laser absorption.
Identifying the parameter regime where this semiclassical picture works, and why, has been the subject of much theoretical work. [41][42][43] However, there has been limited direct benchmarking of numerical and analytical results in regimes of experimental interest. For nonlinear Compton scattering (photon emission by an electron), Harvey et al. 44 compared the frequency and angular spectra predicted by (1) integration of the QED probability rate for a monochromatic plane wave and (2) semiclassical simulation of a 100 fs pulsed plane wave with a super-Gaussian temporal profile, concluding that the neglect of interference effects in the latter caused the harmonic structure to be missed.
In this work, we present systematic comparisons of not only the longitudinal and transverse momentum spectra (Secs. III A and III B) but also the absorption of energy from the background field (Sec. III C). We introduce a normalization framework in Sec. II C that guarantees that we compare precisely the same physical scenario. This permits direct, quantitative benchmarking of semiclassical codes against analytical results from QED in the parameter regime relevant for recent and upcoming experiments.

II. METHODS
The interaction geometry is illustrated in Fig. 1. An electron with initial Lorentz factor c 0 collides head-on with a circularly polarized laser pulse that has dimensionless amplitude a 0 , central frequency x 0 , and invariant duration s. Throughout this work, we set h ¼ c ¼ 1 and denote the elementary charge by e and the electron mass by m. The pulse vector potential eA l ð/Þ ¼ ma 0 gð/Þð0; sin /; cos /; 0Þ, where a 0 ¼ eE 0 /(mx 0 ) for electric field strength E 0 (Ref. 45) and gð/Þ ¼ cos 2 ½/=ð4sÞ for phases j/j < 2ps. In all the results presented here, c 0 ¼ 1000 and x 0 ¼ 1.55 eV (equivalent to a wavelength of k ¼ 0.8 lm). We will consider dimensionless amplitudes in the range of 5 a 0 30, which covers the transition between the weakly and highly nonlinear classical regimes, and restrict the laser duration to be s ¼ 2 or 3 so that the expected number of photons is of order one. This is because our QED calculations are performed for single scattering only, and so that we can gather sufficient statistics in the semiclassical simulations (the fraction of collisions in which only one photon is emitted is exponentially suppressed with increasing a 0 ).
The quantum interaction of charged particles and photons with strong fields is characterized by the invariants v e ¼ e ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ÀðF:pÞ 2 q =m 3 and v c ¼ e ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ÀðF:k 0 Þ 2 q =m 3 , where F is the electromagnetic field tensor and p and k 0 are the fourmomenta of the electron and photon, respectively. 21 v e may be interpreted as a measure of the field strength in the rest frame of the electron relative to that of the critical field of QED E crit ¼ m 2 /e. [46][47][48] It is often referred to as the "quantum nonlinearity parameter" by analogy with a 0 , which is the classical nonlinearity parameter. 21 We have v e $ 0.1 for the interaction parameters under consideration here, and so, quantum effects are non-negligible.

A. QED
The strong-field QED scattering matrix (S-matrix) connects asymptotic free states, evolving the initial state from the distant past to the distant future. The calculation is performed to all orders in the coupling to the background field a 0 , i.e., non-perturbatively, as the number of photons absorbed and reemitted by an electron in an intense laser field is very large. For the lowest order process shown in Fig.  2, the emission of one photon or single nonlinear Compton scattering, 49-52 it reads The delta function ensures the conservation of momentum in the lightfront and transverse directions. By lightfront momentum, we mean p þ k.p/x 0 , which is conserved in a plane wave with wavevector k in the absence of radiation reaction. Other features are the transition operators T j , which are sensitive to the electron spins and the photon polarization, and C j , which are integrals over the laser phase (see the study by Seipt et al. 52 for details) Here, F j are functions of the vector potential A(/) and p(/) in the exponent is the classical kinematic four-momentum of the electron, a solution to the radiation-free Lorentz force equation.
The one-photon emission probability is where we have defined the lightfront momentum transfer fraction f ¼ k 0þ =p þ and normalized transverse photon momentum r ? ¼ k 0 ? =ðfmÞ. The magnitude of the latter r ? ¼ ðp þ =mÞ tanðh=2Þ, where h is the polar angle of the emitted photon, becoming r ? ' ch if c ) 1 and h ( 1. If a 3 0 =v e ) 1 and f is not too small, the phase interval which contributes to the emission of a single photon becomes much smaller than the wavelength of the background field 21 and interference between emission from different formation regions is suppressed. 41 In this case, the field may be treated as constant over the photon formation region. As the photons are emitted into a narrow cone around the direction of the electron's instantaneous momentum, we can integrate over the transverse momenta r ? to obtain the instantaneous probability rate per unit phase and lightfront momentum transfer where z 3=2 ¼ f =½v e ð1 À f Þ and v e v e ð/Þ the local value of the quantum parameter.

B. Semiclassical
In the semiclassical interpretation of the collision process, the electron follows a (radiation-free) classical trajectory between point-like, probabilistically determined, QED events. These events are implemented using the standard Monte Carlo algorithm, 36,37 with rates calculated in the LCFA, i.e., Eq. (4) (see also Refs. 21, 22, and 53). We use circe, a particle-tracking code that simulates photon and positron production by high-energy electrons (and photons) that collide with laser pulses that have a 0 ) 1. In one spatial dimension, the external field is assumed to be a plane electromagnetic wave, i.e., the fields are determined by a single parameter / ¼ x 0 n.x, where x 0 is the wave frequency, n ¼ ð1;kÞ for the direction of propagationk, and x is the four-position of the electron. Between emissions, the electron dynamics are determined by the Lorentz force alone. The spatial components of the four-momentum p ¼ (cm, p) that are perpendicular to the wavevector are obtained by integrating dp where E ? is the electric field at phase /. The remaining components of p are determined by the two conditions p þ ¼ const and p 2 ¼ m 2 . The four-position is determined by Photon emission is implemented as follows: Each electron is assigned an optical depth against emission s ¼ Àlogð1 À RÞ for pseudorandom 0 R < 1, which evolves as where W is the instantaneous probability rate of emission given by Eq. (4), until the point where it falls below zero. Then, the lightfront momentum transfer f ¼ v c /v e is pseudorandomly sampled from the differential rate and s is reset.
Assuming that emission occurs in the direction parallel to the initial momentum, as the electron emits into a narrow cone of opening angle 1/c, the momenta of the electron and the photon after the scattering are As discussed in Ref. 34, this leads to an error in energy conservation of which is small for ultrarelativistic particles.

C. Comparison basis
In this work, we present quantitative, as well as qualitative, comparisons of electron and photon spectra predicted by the exact QED and semiclassical methods. We discuss here how the normalization may be set consistently, but independently by the two methods.
The final result of the QED calculation is the probability P 1 that a single photon is emitted in collision of an electron with a pulsed electromagnetic plane wave. However, even for the short pulses under consideration here, the fact that a 0 > 1 makes it possible for P 1 > 1. Where this occurs, the probability is generally interpreted as the mean number of emitted photons, i.e., P 1 ! N c;QED , as this quantity can certainly exceed unity. 9,44,54 We emphasise that this interpretation is exact only in the classical limit, where recoil can be neglected. The true probability for single scattering is given by the renormalized quantity P 1 =ð1 þ P 1 n¼1 P n Þ. To determine this would require the calculation of the scattering probability to a state containing an arbitrary number of photons n. Efforts to characterize such multiphoton interactions analytically have been limited due to the complexity of the calculations; at present, all the results in the literature are restricted to n 2. For these reasons, we will define the QED "number of photons" as In the semiclassical calculation, multiphoton emission is accounted for by factorisation of the multiphoton emission into a product of first-order processes. Localizing emission in this way allows us to determine the branching fraction to a final state containing an arbitrary number of particles, thereby guaranteeing that P 1 < 1. In the classical limit (i.e., negligible recoil per photon), one emission event is independent of any other and the probability that n photons are emitted in a given collision P n ¼ k n e Àk =n!, where k N c;sc , the mean number of photons in the semiclassical case. However, we consider here collisions where v $ 0.1 and recoil is not negligible. As the emission rate (at fixed field strength) decreases with increasing particle energy, emitting a photon and so losing energy make it more probable that further photons are emitted. As such, the numbers of photons emitted in two non-overlapping intervals are not independent and the probability P n ceases to be Poisson-distributed. In summary, radiation reaction, the recoil due to photon emission, affects the average number of photons because "the emission of each photon modifies the electron state and, consequently, the next emissions." 9 Since it is not possible, as yet, to determine the renormalization factor by which the QED results should be scaled, we propose this alternative. The QED results from Eq. (9) are not scaled. Equivalent semiclassical spectra are obtained statistically, by generating a large set of collision data, accepting only those collisions in which exactly one photon is emitted, and rescaling such that the spectra have integral N c,sc . The mean number of photons N c,sc is determined by considering the entire set of collision data, i.e., where N i is the number of simulated collisions in which exactly i photons are emitted. This definition ensures that only collisions with a single emission contribute to the shape of the spectrum and that its integral may be interpreted in a manner consistent with the QED result. From now on, all differential spectra will be given in terms of the "number of photons" defined by Eqs. (9) and (10). We note that while this post-facto selection criterion lets us compare the same physical scenario as the QED approach, multiphoton and recoil corrections are still present because N c,sc is affected by radiation reaction. We do not necessarily expect the QED "probability" P 1 to satisfy P 1 ! N c;QED ¼ N c;sc for this reason.

A. Lightfront momentum
The symmetries of a plane electromagnetic wave make lightfront momentum u þ n.u a natural choice of parametrization for the differential scattering probability because the conservation of momentum for (single) nonlinear Compton scattering reads =m are the normalized lightfront momenta of the scattered electron and photon, respectively. We will plot differential spectra in terms of the transfer fraction In the back-scattering limit, f ' x 0 =ðc 0 mÞ, where x 0 is the energy of the scattered photon.
We compare the analytical and simulation predictions for the number of photons emitted in the head-on collision of an electron with initial energy c 0 ¼ 1000 and a two-or threecycle laser pulse in Fig. 3(a), with examples of the differential spectra in Fig. 3(c). The percentage difference between the results of the two methods is given for the total number of photons in Fig. 3(b). We find that the semiclassical method systematically overestimates the number of photons but that the fractional discrepancy diminishes with increasing a 0 , falling below 10% when a 0 ! 20. N c scales approximately linearly with a 0 as expected in the strong-field regime; the growing discrepancy towards the lower end of the plotted range is an indication of the transition to the perturbative regime where N c / a 2 0 instead. The origin of the discrepancy is shown in Fig. 3(c). While there is very good agreement for large f, i.e., high energy, the semiclassical method strongly overestimates the number of photons with vanishing f. This is because the underlying LCFA rate contains an integrable singularity / f À2=3 absent in the exact QED calculation. 55 In the latter case, the probability tends to a finite value 43 lim f !0 where g(/) is the envelope function given in Sec. II. It is not surprising that the semiclassical spectra do not reproduce this limit because the LCFA arises from an expansion in the parameter v e =a 3 0 ( 1. 41 The physical meaning of this parameter is that it compares the formation length of the photon to the length scale over which the field varies. If this is sufficiently small, we can assume that emission occurs instantaneously and thereby neglect interference effects. As discussed in the study by Harvey et al., 44 this means that Monte Carlo implementations of localized rates cannot reproduce the harmonic structure that becomes visible in the emission spectrum at small f. A simple way to estimate the smallest f for which the LCFA should be valid is to recall that in a monochromatic plane wave, emission over a complete phase oscillation contributes to photons at the first nonlinear Compton edge, for which the transfer fraction f C ' 2v e =a 3 0 . The requirement that the formation length be smaller than the laser wavelength is then equivalent to having f տ f C . This limit is consistent with the results shown in Fig. 3(c) and with a more detailed calculation performed by Di Piazza et al. 43 In fact, if we cut off the QED and semiclassical spectra below f ¼ f C , the percentage discrepancy in the number of photons falls from 17% to 5% at a 0 ¼ 10 and from 5% to 2% at a 0 ¼ 30.
While it is important to capture the number spectrum accurately, the dynamically significant quantity is the spectrum weighted by u þ c , as this gives the momentum change of the electron or radiation reaction. 56,57 When we compare the total radiated lightfront momentum in Fig. 4(a), we find much better agreement between the two methods, with a relative discrepancy below 4% even for a 0 ¼ 5. This is because the contribution to the electron recoil from the low-energy tail, i.e., the part of the spectrum where the LCFA fails, is very small and diminishes with increasing a 0 . This can also be seen in Figs. 4(c) and 4(d) where we show the intensity spectra without log scaling on the vertical axis.

B. Perpendicular momentum
We parametrize the perpendicular momentum spectrum using the scaled quantity r ? u ? c =f . For c 0 ) a 0 and c 0 ) 1, as we have here, r ? ' c 0 h, where h is the photon scattering angle. The comparison between the analytical and simulation results, shown in Fig. 5, is for r ? scaled by the laser strength parameter a 0 for the following reason.
Analysis of nonlinear Compton scattering in a monochromatic, circularly polarized plane wave only in terms of the number of laser photons absorbed has been shown to reproduce the classical result that the photons are typically emitted along the direction of the instantaneous momentum of the electron in the electromagnetic field. 25 Assuming that the electron and the laser were initially counterpropagating and that both the electron Lorentz factor c 0 and the laser amplitude a 0 ) 1, we have that the most probable angle of emission is tan h ¼ 4a 0 c 0 = ð4c 2 0 À a 2 0 Þ ' a 0 =c 0 . In the semiclassical calculation, we capture the electron's transverse oscillation directly by solving the classical equations of motion and rely on relativistic beaming to justify setting the photon's emission direction to be parallel to the electron's instantaneous momentum. We may derive a scaling relation for the average r ? for the pulsed plane waves under consideration here within the framework of the LCFA. The instantaneous angle between the electron momentum and the laser axis is hð/Þ ' a 0 gð/Þ=c 0 for c 0 ) a 0 ) 1, where g(/) is the pulse envelope described in Sec. II. Assuming that the photon is emitted parallel to the electron momentum, we have that the mean value of r ?
where W(/) is the emission rate [Eq. (4)] integrated over all f, and we obtain an analytical result by working in the classical limit v e ( 1. In this expression, the mean r ? is normalized to the number of photons. Consequently, if v is not too large, we expect hr ? i=a 0 predicted by semiclassical simulation to be independent of a 0 and the pulse duration s. This is indeed what is shown in Fig. 5(a). The exact QED results are generally larger but tend towards the semiclassical results as a 0 is increased. This is because photons are emitted into a broader range of angles in the QED calculation, which can be seen by the fact that the standard deviations of the spectra (normalized to the mean) shown in Fig. 5(c), i.e., the widths of the distribution, are larger in the analytical case. These integral comparisons lead us to expect important qualitative differences between the r ? spectra predicted by QED and by semiclassical simulation. Indeed, Fig. 5(e) shows that, unlike the former, the latter exhibit a universal shape when plotted as a function of r ? /a 0 . Furthermore, they diverge as r ? ! a 0 . No photons are emitted with r ? > a 0 . The range of accessible r ? is much larger for the QED spectra although we note that as a 0 increases, the peak of the spectrum (i.e., the most probable r ? ) tends towards a 0 and the spectra generally become narrower. As we have scaled r ? down by a 0 , this indicates that the characteristic width of the spectrum is approximately constant for all a 0 .
The cause of these differences is the assumption of collinear emission in the semiclassical simulations, i.e., assigning the final momenta according to Eq. (7). Relativistic beaming means that the photon is emitted forward into a cone of half-angle $1/c 0 , corresponding to a width in r ? of r r $ 1. We have neglected this extra angular divergence, which is why the semiclassical results have a sharp edge at r ? ¼ a 0 . The range of angles at which a photon can be emitted is bound by the angle between the electron's instantaneous momentum and the laser wavevector, which is at most a 0 /c 0 for the circularly polarized pulses under consideration here. A more subtle discrepancy may be seen for small r ? in Fig. 5(e). Whereas the analytical results tend to zero as r ? does, exhibiting a pronounced shoulder as they do so, all the semiclassical spectra tend to a finite value of approximately 0.022. (This is consistent with a classical calculation of the angular spectrum, which gives dN c =dr ? ' 5a= ffiffi ffi 3 p ' 0:021 for r ? ¼ 0.) The difference is a consequence of the LCFA: recall that in the semiclassical approach, photons are emitted parallel to the electron's instantaneous momentum. Therefore, photons with small r ? originate from the leading and trailing edges of the pulse, where the local field strength is small and so too is the angle between the electron trajectory and the laser wavevector. In these regions, the effective a 0 is small enough that interference effects become important, suppressing photon emission.

C. Absorption
Energy-momentum conservation demands that the emission of a photon by an electron in a strong background field be accompanied by the absorption of energy from that background field. As the background under consideration here is an electromagnetic wave, this can be interpreted as the absorption of a certain number ' of photons. Seipt et al. 25 have shown that in a circularly polarized, monochromatic plane wave with strength parameter a 0 , the emission of a photon with quantum parameter f ¼ v c =v e is associated with the absorption of ' ¼ s laser photons, where For the short pulses in this work, the probability that ' photons are absorbed is determined numerically.
In the semiclassical method, the background field is treated entirely classically. Nevertheless, we may define an equivalent number of absorbed photons by dividing the classical work done on the moving charge by the laser frequency Note that for a plane wave, the above integral would be identically zero in the absence of radiation. (The same result holds in the QED calculation: if no photon is emitted, ' ¼ 0.) circe computes Eq. (15) for each test electron, integrating the work done across the entire trajectory of the electron in the pulse. A comparison between the total number of absorbed photons N L Ð ' dN c , as computed by the two methods, is shown in Fig. 6(a). We find that the semiclassical method systematically underestimates the absorption but that this difference occurs at the level of a few percent and decreases with increasing a 0 . This is in contrast to what we found for the number of photons and the total radiated lightfront momentum, where the semiclassical result was generally larger than the QED result. In all three cases, we expect errors to arise due to the finite size of the formation length and the associated interference; however, here, our results imply that there is some "missing" absorption.
In their analysis of electron-positron pair creation by a photon in a strong laser field, Meuren et al. 42 divided the absorbed energy into "classical" and "quantum" parts, with the former being the work done accelerating the daughter particles out of the laser pulse and the latter being the absorption of photons over the formation length. They showed that the classical part scales approximately like a 3 0 =v e and the quantum part like a 0 /v e , concluding that the classical absorption should be dominant at high intensity. This is consistent with the results presented here, in that we capture the acceleration of the electron post-emission but not any absorption over the formation length, which is assumed to be vanishingly small. Recall that in the semiclassical simulations, there is an error in energy conservation due to the assumption of collinear emission (see Sec. II B). Equation (8) predicts that the magnitude of this error is 2.8% at a 0 ¼ 10 and 0.3% at a 0 ¼ 30. This is comparable to the discrepancy shown in Fig. 6(b) but for the largest a 0 , where recoil corrections to N c,sc take effect.
Equation (14) indicates that the larger the lightfront momentum transfer f, the more the photons are absorbed from the external field. Both the lightfront momentum transferred to an individual photon and the number of emissions increase with a 0 , and so, we expect the absorption to increase as well. Integrating Eq. (14) weighted by the emission rate Eq. (4) over all f, we find that N L $ a 2 0 I þ . Here, I þ is the total radiated lightfront momentum given by Eq. (12), which scales like a 2 0 in the classical limit. Then, we expect N L $ a 4 0 , which agrees reasonably well with a power-law fit to the data shown in Fig. 6(a); we find N L / a 3:7 0 for both the QED and semiclassical results. The true scaling is weaker than a 4 0 because of quantum corrections that reduce the radiated power. 22 We showed in Sec. III A that the semiclassical method predicts the large-f part of the emission spectrum accurately even for a 0 ¼ 5. As this part of the spectrum is associated with the largest ', the agreement between the QED and semiclassical results should be best for ' ) 1. Four examples of the spectrum of probable ' are shown in Fig. 6(c). The agreement between the two is excellent for ' > 10, but the semiclassical method fails to capture the small-' part of the spectrum accurately. This is because it localizes emission, thereby neglecting interference effects; these suppress the emission probability for small ' and give rise to the harmonic structure that can be seen in the QED spectra.

D. Exemplary case
Finally, we present a comparison for a specific set of collision parameters, drawing on the systematic results we have so far, to discuss the role of multiple emissions. Figure  7 shows the full set of double-and single-differential photon spectra for lightfront momentum, perpendicular momentum, and absorption for a collision between an electron with c 0 ¼ 1000 and a laser pulse with a 0 ¼ 20 and s ¼ 2. The average number of photons is N c ¼ 2.36 for these parameters, and therefore, multiphoton effects should be taken into account. However, as the QED calculation is performed only for single scattering, we filter the semiclassical collision data to ensure that the same physical scenario is being compared. Now, we can show the effect of this filtering on the semiclassical results.
The double differential spectrum obtained when all emissions are taken into account is shown in Fig. 7(b); when only single scattering events are binned, we obtain the spectrum shown in Fig. 7(c). As discussed in Sec. III B, in the latter case, we find a sharp cutoff at r ? ¼ a 0 as this is the largest angle between the electron momentum and the laser wavevector, and we assume that photons are emitted in the collinear direction. Then, photons with r ? > a 0 can only come from secondary scattering. Notice that while the QED result [ Fig. 7(a)] is generally broader in the vertical direction, the probability that r ? > a 0 diminishes with increasing f, which apparently justifies the assumption of collinear emission for f $ 1. The QED result is also smoother as it is free from the numerical noise inevitable in Monte Carlo simulations. The single scattering spectrum appears noisier as it represents only 20% of the collision data. Figure 7(d) shows that the effect of the filtering on the lightfront intensity spectrum is rather small. Nevertheless, the agreement is better when only single scattering is included. The spectrum in this case is slightly harder, as secondary photons tend to be emitted with smaller f. Taken in isolation, Fig. 7(e) appears to suggest that the agreement between the QED and semiclassical spectra is better if we do not filter the collision data. The shape of the peak, if not its maximum value, is actually captured better when all emissions are included. This is coincidental. Recalling that we have assumed collinear emission in the semiclassical approach, we would expect the QED result for double Compton scattering to be even broader in r ? . In fact, the most significant difference is found when we compare the probability that ' photons are absorbed from the laser pulse in Fig. 7(f). The two semiclassical spectra have the same integral (by construction, see Sec. II C), and therefore, small values of ' must be suppressed when multiple emissions are included. The probability that ' ¼ 2 Â 10 4 , for example, is 6Â larger when all emissions are accounted for. Without the selection procedure we have introduced, it would not be possible to compare these against the QED result in a consistent way.

IV. DISCUSSION
It is generally assumed that the semiclassical approach used in modelling high-intensity laser interactions is valid when a 0 ) 1 and a 3 0 =v e ) 1. The precise value of a 0 for which these conditions are satisfied depends, however, on the particular quantity that is being calculated. We have shown that for a 0 as low as 5, semiclassical codes accurately capture the part of the emission spectrum for which f $ 1. On the other hand, there is still a significant discrepancy in the total number of photons even when a 0 ¼ 30. The error is concentrated in the low-f part of the spectrum, i.e., photons with low energy and large angle; clearly, a semiclassical code should not be used to predict the result of an experimental measurement in the spectral region f Շ 2v e =a 3 0 . Improving these codes could be accomplished by calculating this threshold and replacing the LCFA rate for photons with smaller f (Ref. 43) although this does require the external field to be treated as a slowly varying plane electromagnetic wave.
Alternatively, it might be possible to bypass this problem by using a PIC code, in which a hybrid description of the electromagnetic field is employed. The Nyquist frequency associated with the finite spacing of the grid naturally separates radiation into two components: lower frequencies are resolved on the grid, i.e., classically, and higher frequencies are treated as "photons," just as we outlined in Sec. II B. This is why many implementations of QED processes in PIC codes include the option of a low-frequency cutoff below which macrophotons are not emitted. 58 It would be interesting to compare the QED results in this work with the predictions of a semiclassical code in which the radiation spectrum below a certain cutoff is obtained by Fourier analysis of the Li enard-Wiechart potentials in the far field. This would ensure that the formation length is resolved at low f, thereby capturing interference effects. It is reasonable to expect a classical description to be appropriate because both the quantum corrections and the electron recoil should be negligible for photons with f ( 1. While it is important to make these improvements at low f, this part of the spectrum contributes negligibly to the momentum change of the electron, which is dominated by photons with large f. As the agreement between the QED and semiclassical spectra is much better here, it is not surprising that we find the average lightfront momentum loss predicted semiclassically to be within a few percent of the QED value even at a 0 ¼ 5. At a 0 ¼ 10, for example, the error in the total number of photons is 16% for both s ¼ 2 and 3, whereas in the total radiated momentum I þ , it is 1.5%, an order of magnitude better. For the experimental parameters in the study by Cole et al. 17 (a 0 ' 10; v e ' 0:1, and c 0 ' 1000), the lower limit on f is equivalent to x Շ 100 keV; this part of spectrum represents approximately 16% of the total number of photons but only 0.04% of the total radiated energy, using their parametrization of the spectrum and the measured critical energy of 30 MeV. This is encouraging for semiclassical or PIC-based modelling of radiation reaction of an electron population. In laser-beam interactions, the particle number density is generally low enough that the radiation spectrum can be obtained by incoherent summation over all photons emitted by the individual particles; this allows the treatment of a beam with a spectrum of energies and non-zero divergence. Although these two effects will wash out, for example, detailed harmonic structure in the momentum spectrum, 59 the overestimate at low f will survive. Nevertheless, as radiation reaction is an intrinsically multiphoton process, 9 proper benchmarking requires the calculation of the higher-order diagrams shown in Fig. 2. The selection and scaling scheme we have presented here could easily be extended to comparisons with QED calculations of double, triple, etc., nonlinear Compton scattering.
Perhaps more important for the case of multiple emissions are the comparisons we present for the perpendicular momentum spectra. These test the assumption of collinear emission, which is distinct from the LCFA. While the peak at r ? ' a 0 is common to both QED and semiclassical results, the width of the distribution around this point is not captured semiclassically. The angle at which the photon and the electron travel after the scattering affects their quantum parameter, and so, the rates at which secondary processes occur. Even if the change to the rates is small, the cumulative error could become large if the multiplicity is high. Accurate modelling of the angular spectrum is important because, for example, the transverse broadening of an electron beam has been proposed as a signature of quantum effects on radiation reaction. 60 This broadening would be in addition to that from the finite beaming of the radiation and any initial divergence of the beam (a few milliradians for the laser-wakefield-accelerated electron beams reported in Refs. 17 and 19). It would also be important for the study of QED avalanches, in which even a single electron accelerated by counterpropagating lasers can seed the creation of a critically dense electronpositron pair plasma. One possible approach would be to implement an angularly resolved LCFA rate that includes the finite 1/c beaming of the scattered photon (see Ref. 53 for example).
We have also found that the absorption of energy from the background is dominated by the "classical" component, i.e., j Á E work done by the external field in accelerating the scattered electron. This is assumed to be the case in PIC modelling of QED avalanches and suggests that a classical treatment of backreaction is reasonably accurate at high intensity. (Recall that in these codes, the fields and currents are evolved self-consistently but classically.) There remains the question of the "missing" absorption we discussed in Sec. III C. On the one hand, the fraction of the total depletion diminishes with increasing a 0 ; however, if this error does arise on a "per-emission" basis, the increased multiplicity at high a 0 could mean that it becomes significant. It is reasonable to expect a causal relationship between the absorption discrepancy and the assumption of collinear emission. This is something we will consider in future work.

V. CONCLUSIONS
In this work, we have presented benchmarking of semiclassical simulations against exact QED results for nonlinear Compton scattering in an intense laser pulse, using a method that guarantees that we compare precisely the same physical scenarios.
The differential spectra agree both qualitatively and quantitatively in the dynamically important region f տ 2v e =a 3 0 that dominates the electron recoil and absorption from the laser fields. We find that the lightfront momentum loss and the number of absorbed photons from semiclassical simulations are within a few percent of the exact QED results for a 0 > 5. However, improvements are clearly called for at low f, where the LCFA breaks down, and in the angular distribution, where the agreement is only qualitative due to the assumption of collinear emission in the simulations.
It remains to be seen whether improving these will lead to significant differences in the results of simulations of laser-plasma interactions. In deciding what is most important, we should be guided by further comparison with QED calculations that include multiple emissions. These will place more stringent limits on the validity of the approximations that underpin the semiclassical approach. The fact that experimental exploration of the strong-field, multiphoton regime (a 0 $ 10, v e $ 0.1, N c $ 10) is already underway makes this an urgent question.