From local to nonlocal: higher fidelity simulations of photon emission in intense laser pulses

State-of-the-art numerical simulations of quantum electrodynamical (QED) processes in strong laser fields rely on a semiclassical combination of classical equations of motion and QED rates, which are calculated in the locally constant field approximation. However, the latter approximation is unreliable if the amplitude of the fields, $a_0$, is comparable to unity. Furthermore, it cannot, by definition, capture interference effects that give rise to harmonic structure. Here we present an alternative numerical approach, which resolves these two issues by combining cycle-averaged equations of motion and QED rates calculated in the locally monochromatic approximation. We demonstrate that it significantly improves the accuracy of simulations of photon emission across the full range of photon energies and laser intensities, in plane-wave, chirped and focused background fields.


I. INTRODUCTION
The collision of multi-GeV electron beams and intense laser pulses is a promising scenario for precision measurements of quantum electrodynamics (QED) in the strong-field regime, where both the normalised amplitude of the laser, 0 , and quantum nonlinearity parameter of the electron, , exceed unity. Perturbative QED calculations of the interaction fail once 0 1 and must be replaced by 'all-order' approaches, which take the interaction with the strong background field into account exactly [1,2]. While the theory for this regime is now several decades old [3], experiments are limited in number. In the weakly multiphoton regime, 0 0.4, laser-electron collision experiments have observed Compton scattering (photon emission) and trident electron-positron pair creation [4,5]. At higher values of 0 , but small , they have observed photon emission in the classical regime (nonlinear Thomson scattering) [6][7][8][9] and at 0 10, radiation reaction (multiple photon emission) in the nonlinear classical [10] and quantum regimes [11]. However, as yet, there are no experimental measurements charting the transition between the perturbative, multiphoton, and nonlinear regimes, 0.1 0 10 at 1. This is likely to change in the near future, as increasing interest in strong-field QED has led to planned experiments that will combine conventional electron accelerators with intense optical lasers [12,13].
The transition regime represents a particular challenge for theory and simulation. A perturbative approach is not sufficient once 0 1. However, neither is an approach based on the locally constant field approximation (LCFA) [1,14], as this applies only in the opposite limit, 0 1: this approximation underpins the simulation codes [15][16][17] used to model QED effects in laserplasma interactions [18][19][20][21][22][23][24], which will be explored in the next generation of multi-petawatt laser facilities [25][26][27][28]. The versatility of the LCFA comes from its local nature and the neglect of interference effects, i.e. the finite size of the spacetime region over which QED processes take place, which requires both 0 1 and 3 0 / 1; the limitations of doing so have been thoroughly discussed in the literature [29][30][31][32][33]. Experiments that aim at precision measurements of strong-field QED demand precision simulations of the interaction. However, in the transition regime, the error made by simulations based on LCFA rates is unacceptably large.
In this paper, we present a simulation framework that overcomes these issues by using the locally monochromatic approximation (LMA) instead. This achieves greater accuracy by taking into account interference effects at the scale of the laser wavelength, which is possible provided that the laser pulse is relatively unchanged by the collision with a probe electron beam. To do this, we combine classical trajectories, defined on a cycle-averaged basis, with probability rates that treat the background 'locally' as a monochromatic plane wave, with an amplitude and frequency that can vary in space and time. As such, we exchange the ability of the LCFA to model an arbitrary electromagnetic field for significantly increased accuracy in the modelling of plane-wave-like fields.
The derivation of the LMA in Heinzl et al. [37] assumes a plane wave, whereas any experimental configuration will employ a focused laser pulse. This makes it essential to consider beyond-planewave field configurations, for which exact theoretical results are limited in number [43,44]. In order to make progress, we consider the case of plane-wave backgrounds that have a nonlinear dependence on phase, or a 'chirp', which results in a localisation of both the wave's amplitude and frequency. By allowing both the amplitude and wavevector to vary in space and time, we gain analytical insight into the case of a focused background, where this would also be the case. We then describe how the LMA may be implemented in numerical simulations of photon emission and benchmark their predictions against strong-field QED for pulsed plane waves (unchirped and chirped) as well as with focusing pulses. For the last of these, we must employ an approximate solution to the Dirac equation [45][46][47], which, to the best of our knowledge, has not previously been compared to a simulation. Our results confirm that simulations based on this framework may be used for precision modelling of experiments, with an accuracy of a few percent in the integrated probability (improving on the accuracy of the LCFA by orders of magnitude in the transition regime), and correct reproduction of harmonic structure in the differential spectrum, which has been identified as an aim of future experiments [13].
In the following, we use a system of units in which the Planck's reduced constant, the speed of light and the vacuum permittivity are all set to unity: ℏ = = 0 = 1. The electron mass is denoted by . The fine-structure constant is related to the elementary charge by = 2 /(4 ).

II. THEORY BACKGROUND
We begin with an explanation of how the full QED plane-wave results are calculated, as well as a summary of the main details arising from the analytical calculation underpinning the LMA.
(Many papers have investigated the effect of pulse shape on nonlinear Compton scattering, see e.g. [48][49][50][51][52][53].) For concreteness, we specify from the outset that we will be assuming a background that is a circularly polarised, chirped, plane-wave pulse with potential . We define the dimensionless potential = / , where 0 is the dimensionless intensity parameter [54] (also called the "classical nonlinearity", normalised amplitude or the strength parameter) and , are orthonormal polarisation vectors obeying · = · = −1. Throughout, we use lightfront coordinates which depends on the lightfront phase = · (where = + + is the background wavevector), and the pulse phase duration, Φ, is related to the number of cycles, , via Φ = 2 . The function ( ) describes the chirp of the background. For a pulse without chirp, is linear in , i.e.
( ) = 0 for all . (In the following, we will pick ( ) = for the unchirped case.) We use the scattering matrix approach [55] to calculate the probability of single nonlinear Compton scattering from a single incoming electron colliding with a plane-wave background. We can write the scattering matrix element as: where / * , is the polarisation of the emitted photon with 4-momentum and Ψ , (Ψ , ) is the Volkov wavefunction [56] of the incoming (outgoing) electron: The matrix element can be simplified to: where = · / · is the lightfront momentum fraction of the emitted photon, 0 = · / 2 is the initial energy parameter of the probe electron,˜contains normalisation constants, the instantaneous electron momentum is given by and the regularising factor Δ = 1 − · / · incorporates all the contributions from phases outside of the integral. The total probability can be written: where ì ⊥ = ì ⊥ /( ) − ì ⊥ / contains the shifted perpendicular momentum. Here " ⊥ " indicates directions perpendicular to the background propagation direction and · pol. indicates an average over initial and sum over final polarisation states. The numerical results in exact QED are calculated by evaluating eq. (6) directly: the matrix element in eq. (2) was evaluated using photon polarisation eigenstates of the background [57] and spin states in the Lepage-Brodsky convention [58].
Rather than direct numerical evaluation, some of the integrals in eq. (6) can be evaluated analytically by generalising the locally monochromatic approximation [37] to arbitrarily chirped plane-wave pulses. In the following, we present an overview of this approach, and direct the reader to appendix A for details.
The background field is given by eq. (1). For the LMA to approximate the emission spectrum well, the envelope function ( /Φ) should be slowly varying with respect to the carrier frequency, 1 for the unchirped case, which corresponds to a many-cycle pulse). However, in this work, we also include the chirp. Therefore we will also make a "slowly varying chirp" approximation (see e.g. Seipt et al. [59]). These approximations then allow the squared Kibble mass, , which occurs in an exponent, to be integrated over. The Kibble mass takes the form = 1 + ì 2 − ì 2 , where denotes a phase-window average. In the case of a circularly polarised background, the slowly varying (envelope) and rapid (carrier) timescales occur in ì . We can demonstrate the approximation by considering a single component of ì, e.g. ì · ì.
Now, one can introduce a local frequency scale, ( ) = ( ) and integrate by parts as in eq. (A6).
The fast timescale of the cosine term is included exactly. The remaining terms for the envelope and chirp variations have a size, relative to the leading term, of the order of respectively (neglecting a rapidly varying term that appears ∼ cot ( )). As long as the magnitudes of both of these are much less than unity, we should expect the slowly varying approximation to be good. (The same arguments apply to the ì · ì term, whereas ì 2 is not affected by chirp in a circularly polarised background.) Beyond the additional constraints on the chirp, no further modifications to [37] are required in the derivation (more details are given in Appendix appendix A).
Finally, we arrive at P LMA = ∫ LMA , where: where 2 rms ( ) = 2 / 2 − 1 and ( ) = [ ( )] 0 , with 0 = · / 2 the unchirped energy parameter. Here = is the quasimomentum, the laser-cycle-average of the instantaneous electron momentum given in eq. (5). The appearance of a local wavevector in ( ) also follows from considering components of the field-strength tensor, , for the chirped pulse in eq. (1), which contain terms ∼ ( ) / , where ( ) = ( ) . P mono is the probability of nonlinear Compton scattering into the th harmonic in a monochromatic background, is the proper time, related to the phase by / = 1/( 0 ). The approximation is locally monochromatic because the intensity and energy parameter occurring in the monochromatic probability now take the (cycleaveraged) local value at the position of the electron. The integrand is given explicitly by eq. (A28) for nonlinear Compton scattering. Unlike the monochromatic case, here the harmonic range is phase-dependent: where ( ) is the edge of the classical (nonlinear) harmonic range.
To obtain the probability of Compton scattering in a focused laser background, we must use some approximation, as analytical solutions to the Dirac equation in a realistic focused laser background are unavailable (some progress has recently been made in this direction: see e.g. [43,44]). One method is to find an approximate solution to the Dirac equation using a WKB expansion in a small parameter −1 , where is the initial relativistic gamma factor of the incident electron [45][46][47].
Then assuming 0 , for a head-on collision of the electron probe with the focused laser pulse, one can write: where is the electron probe areal density and the plane-wave probability, P from eq. (6), now has an intensity parameter which can depend on the perpendicular spatial co-ordinate.

III. IMPLEMENTATION IN NUMERICAL SIMULATIONS
The inclusion of strong-field QED processes in numerical simulations, such as the particle-incell [15,16] or particle-tracking codes [34,35,60] used in plasma and beam physics, is based on a semiclassical treatment of particle dynamics, which combines classical trajectories with the use of probability rates [61]. This is motivated by the appearance of the classical kinetic momentum , eq. (5), in the QED scattering probability, via the exponent of the Volkov wavefunction, eq. The rate in the locally monochromatic approximation, by contrast, is derived assuming that the envelope of the potential, rather than the potential itself, is slowly varying. Averaging over the fast timescale, the laser period, means that the quantity that enters the rate, and also the conservation of momentum, is not the kinetic momentum directly, but rather the quasimomentum ≡ [1,64]. In a plane wave, = − + (2 · − 2 2 )/(2 · ) and 2 = 2 , whereas = + 2 2 rms /(2 · ) and 2 = 2 (1 + 2 rms ), for 2 rms ≡ − 2 . In contrast to the LCFA case, the rate is a function of two parameters: the normalised amplitude (or intensity parameter), rms , and the energy parameter ≡ · / 2 , both locally defined.
The trajectory obtained from these two equations does not include the fast oscillation at the timescale of the laser period, as shown on the right-hand side of fig. 1. This does not mean that the physical effect of that oscillation is lost: it is accounted for in the emission rate. To see this more clearly, note that at fixed , in the limit 0 1, there is a most probable harmonic index [66]. Combining this relation with the conservation of quasimomentum, which reads 2 angle is rms / for 0 1 [66] (see also [64]). Thus an equivalent angular structure emerges, provided that the classical trajectory is parametrised in terms of quasimomentum. The conceptual differences between LCFA-and LMA-based simulations are summarized in table I.
The emission of photons, and its effect on this trajectory, is modelled in the following way.
At any particular timestep, we have the electron quasimomentum and position from the classical equations of motion, as well as the local values of the laser normalised amplitude rms ( ), wavevector ( ) and polarisation (taken to be circular throughout). In fact, and are sufficient to determine the properties of the emission, as they define the two invariant parameters, rms and , that control the rate and the conservation of momentum. This is given by where is the electron quasimomentum after the scattering, is the momentum of the emitted photon, and is the harmonic index (the net number of laser photons absorbed). The emission rates themselves control and subsequently ≡ · / · , the lightfront momentum fraction. Given , and , it is a matter of kinematics to determine and then . Our Monte Carlo algorithm is as follows: (i) advance the electron trajectory by solving eqs. (12) and (13), (ii) evaluate, at every timestep, the probability of emission and pseudorandomly decide whether to emit a photon or not, and on those timesteps where emission takes place, (iii) select a harmonic index with probability / , where is the partial rate and = ∞

=1
is the total rate, (iv) sample from the partial spectrum (d /d )/ , (v) determine given , and and (vi) reset the electron quasimomentum from to .
The probability that emission takes place in small interval of lab time Δ is given by P = Δ and Δ = Δ ( / 0 ) is the equivalent interval of proper time. We obtain by integrating, and then summing, the partial, differential rates of emission , which are given by [37] d d The argument of the Bessel functions (of the first kind [67]) and auxiliary variables are (16) and the bounds on are 0 < < /(1 + ). Note that depends on rms and and is therefore a function of proper time , as shown explicitly in eq. (10). While the summation should run from = 1 to infinity, it is sufficient to sum up to a largest value max = 10(1 + 3 rms ). In principle, the integration and summation can be done at every timestep, given the particular values of rms and . However, it is significantly faster to obtain by interpolating from a lookup table, where ( rms , ) is precalculated over the domain min rms < rms < max rms and min < < max . The upper bounds are fixed by the problem space under consideration; we have taken max rms = 10 and max = 2 in our code. The lower bounds are chosen such that alternative sampling strategies may be used.
First, if rms < min rms 1, only the first harmonic, = 1, contributes significantly to the probability. In this limit, the rate may be obtained analytically: Second, if < min 1, we may take the classical limit, whereupon the partial rates become: but where we fix = (1 + )/ to be 0 < < 1. Equation (18), integrated over 0 < < 1 and summed over = 1 to max , is tabulated over the same range min rms < rms < max rms . In our implementation, min rms = 0.02 and min = 10 −3 . Thus at every timestep, the emission probability P = Δ is obtained by interpolating from the appropriate lookup table, or using the limiting analytical expression. Emission is deemed to occur if a pseudorandom number , drawn from the uniform distribution (0, 1), satisfies < P.
If emission takes place, the next step is to determine and . The former is obtained by solving for , = =1 / , where is another pseudorandom number drawn on the unit interval (0, 1). In our implementation, the total rate of emission is already available at this point; however, the sequence of partial rates must be evaluated explicitly, by integrating eq. (15) over .
We do this, rather than store a lookup table in (as well as in rms and ), because unlike the total rate, which is needed at every timestep, the partial rates are only needed on emission, which occurs at infrequent intervals. Once is fixed, the lightfront momentum fraction transferred, , is obtained by rejection sampling of eq. (15).

IV. BENCHMARKING
While LMA rates have already been implemented in simulation codes used to study laserelectron interactions [34][35][36], the accuracy of these simulations has not been thoroughly benchmarked against the underlying theory. Placing quantitative bounds on the error made, is essential for experiments that aim for precision characterisation of strong-field QED processes [13]. These analyses have been performed for LCFA-based simulations, however: see [29,30,68] and proposed improvements in [31][32][33]. In this section, we compare the results of simulations based on the LMA, as outlined in section III, with QED theory calculations without the LMA, for photon emission in a pulsed, plane-wave background. We focus on the transition regime 0 ∼ 1, where currently existing approaches based on the LCFA are likely to fail. The laser pulses we consider are circularly polarised with a cosine-squared temporal envelope: the potential ì( Here is the number of cycles corresponding to the total duration of the pulse. One may estimate the (intensity) full-width-at-half-maximum duration of this pulse as [fs] [μm]/0.8. The function ( ) controls the frequency chirping of the pulse and is initially set to ( ) = (i.e., unchirped) for the results in section IV A. The electrons counterpropagate head-on to the laser pulse, with initial energy parameter 0 = 0.1. This is equivalent to an initial Lorentz factor of 0 = 1.638 × 10 4 for a laser wavelength of 0.8 μm.
The theoretical calculations described in section II are for single emission only. However, for sufficiently large 0 or pulse length , it is possible for the total probability of emission P to exceed unity. This indicates that higher order processes, including the emission of multiple photons by a single electron, become important. Simulations model multiple emissions as the incoherent combination of single-vertex processes, transporting the electron classically between emission events. This is motivated by theoretical calculations of higher order processes which show that part of the probability can be factorised into a product over polarised, first-order processes [69][70][71].
Neglecting other contributions, where the intermediate state does not propagate, is expected to be a good approximation if 2 0 Δ 1 [72], where Δ = 2 is the phase duration of the pulse, which allows simulations to model cascades of photon emission and pair creation [61]. In the present case, we consider only the comparison for single photon emission results. Therefore, the probability obtained theoretically is interpreted as the average number of emitted photons [73]. As our simulations allow for an arbitrary number of emission events per electron, we obtain equivalent results by artificially disabling recoil, i.e. the electron momentum is not changed self-consistently when a photon is emitted. The number of emitted photons therefore scales exactly linearly with pulse duration. This does not apply to the theoretical results.
The symmetries of a plane wave suggest that the photon spectrum is best characterised in terms of the lightfront momentum fraction, , and normalised perpendicular momentum ⊥ = ⊥ /( ).

A. Pulsed plane waves
The three plots accompanying each double-differential spectrum compare lineouts at fixed ⊥ against theoretical results. The simulations capture the position and overall shape of the harmonics well, but miss the subharmonic substructure visible in fig. 2(f) and (g) in particular. This structure arises from interference effects at the scale of the pulse envelope, whereas the LMA accounts only for interference effects at the scale of the wavelength. The LCFA, by contrast, captures neither, which causes the spectra to be smeared between the clear peaks seen in both the theory and LMA simulation results [29].
Single-differential spectra, i.e. the results from fig. 2 integrated over ⊥ , are shown in fig. 3.
We compare the simulation results with QED for normalised amplitudes 0 = 0.5 and 2.5 and for pulse durations equivalent to = 4 and 16 cycles. The agreement is much better for the longer pulse, which we expect because the LMA neglects terms of order 1/ (see eq. (8) and [37]). The  fig. 3(e), is ponderomotive in origin: it is radiation from the slow decrease and increase of the electron momentum caused by gradients in the intensity profile [74]. While this is accounted for at the level of the classical trajectory in the simulations, its contribution to the emission spectrum is neglected. The peak moves towards smaller as increases and it is eventually lost in the monochromatic limit [37]. Integrating over the -weighted probability, shown in fig. 3(c) and (e), yields the total lightfront momentum transfer from electron to photon. If 0 > 1, this is dominated by contributions from > * 1 , where the LCFA works well [30]. However, it is evident from fig. 3(c) that the LCFA fails globally for 0 < 1. is particularly dramatic for the probability, where the error made is larger than 10% even when 0 = 5. The average lightfront momentum fraction is more sensitive to the contribution of higher harmonics, i.e. large ; as this is where the LCFA works rather well, the accuracy for is better than that for P. However, the LMA simulations are significantly more accurate when 0 1.

B. Chirped pulses
In Heinzl et al. [37], the LMA is derived for a pulse in which the amplitude is slowly varying.
However, a monochromatic plane wave is defined by both an amplitude and a frequency. By extending the LMA to the situation where both may vary with phase, it becomes possible to simulate radiation generation in chirped laser pulses in the transition regime 0 ∼ 1. In this section we benchmark our simulation results against theory for this case.
Furthermore, for a particular pulse duration, there is an upper bound on the largest chirp that can be obtained [75]. In our notation, this maximum is given by max 10/ . We note that chirping a pulse, which is accomplished by introducing a frequency-dependent phase shift, also changes its duration and peak amplitude; we neglect these such that the only difference between the chirped and unchirped case is the variation of the instantaneous frequency.
We compare the photon spectra obtained from theory and LMA-based simulations for 0 = 0.5, double-differential spectrum at fixed ⊥ demonstrate that chirping smooths out the subharmonic structure; as a consequence, simulation results appear to be more accurate than in the unchirped case.
The second example we present is that of a highly nonlinear chirp, where the instantaneous frequency varies in such a way as to compensate for the classical broadening of the photon spectrum at 0 > 1. In a pulsed plane wave, the position of the first harmonic edge varies from = 2 0 /(1 + 2 0 ) to = 2 0 /(1 + 2 0 + 2 0 ) as the cycle-averaged potential rms ( ) sweeps up and down. As such, the on-axis emission is broadband unless the intensity is rather low. In order to overcome this, and obtain a narrowband source of Compton rays even when 0 is not small, it has been proposed to chirp the pulse in a particular way [76][77][78][79][80] sensitive only to the pulse envelope (for circular polarization) [31,37].

V. FOCUSED LASERS
Theoretical calculations of strong-field QED effects in experimentally relevant scenarios must deal with three-dimensional effects: the nonlinear regime 0 1 is reached by focusing laser light to a spot of small, even diffraction-limited, size, so the laser pulse will differ significantly from a plane wave; the electron beam that probes the laser will also have finite size and temporal duration.
Theoretical results build upon analytical solutions of the Dirac equation in a background field and are therefore only available for plane waves, focusing models of very high symmetry [43,44], or under a high-energy approximation 0 [45,47]. In this section, we discuss the application of simulations, based on LMA emission rates, to model the interaction of electron beams with focused laser pulses.
Within the LMA, the field is treated locally as a monochromatic plane wave. In order to model a focused laser pulse, we therefore promote the cycle-averaged amplitude rms and wavevector to be functions of spatial coordinate as well as phase. For Gaussian focusing, within the paraxial approximation, we have where 0 is the beam waist (the radius at which the intensity falls to 1/ 2 of its central value), = 2 0 / is the Rayleigh range, and the factor ( ) is the pulse envelope [81]. The local wavevector = , where = − 2 /(1 + 2 ) + tan −1 is the total phase. However, in what follows we neglect the wavefront curvature and Gouy phase so that = and takes its usual, plane-wave value. We compare the results so obtained with simulations based on the LCFA, which is a more standard approach [15,16]. In the LCFA simulations, the laser pulse is defined using the paraxial solution for the fields given in [82]: we include terms up to fourth-order in the diffraction angle = 0 / in the Gaussian beam, which is then multiplied by a temporal envelope ( ).
Electron trajectories are determined by solution of the ponderomotive force equation, eq. (12), for the quasimomentum, or the Lorentz force for the kinetic momentum, as appropriate.
First, we verify that LMA and LCFA simulations yield consistent results in a regime where they are expected to do so. We consider a laser pulse that is focused to a spot size 0 = 2 μm, reaching a peak amplitude of 0 = 10, with Gaussian temporal envelope of (full width at half maximum) duration 30 fs. The electrons have initial energy parameter 0 = 0.01 (equivalent to 0 = 1638, given a laser wavelength of 0.8 μm) and are initially counterpropagating, with zero initial divergence. Their initial positions are distributed over a disk of radius 0 = 0 , such that they encounter a range of peak intensities. We have both 0 1 and 2 0 / 0 1, so the LCFA is expected to be a good approximation. The results presented in fig. 7 are obtained from simulations of this scenario using the LMA and LCFA, with recoil on photon emission artificially disabled.
This means that the electron trajectory is determined solely by the action of the laser fields, allowing us to confirm the equivalence between the LMA and LCFA at the level of the electron dynamics, illustrated in fig. 1. Figure 7 shows the angular distributions of the electrons and emitted photons, after the collision has taken place. We see immediately that the LMA and LCFA simulations yield almost identical results. In order to explain the double ring structure evident in the electron distributions, we derive an approximate, analytical prediction for the expected ponderomotive scattering angle .
Consider an electron that is initially counterpropagating, with no initial transverse momentum, at radial distance (impact parameter) from the laser axis, at ultrarelativistic velocity such that 0 − 3 ⊥ . We approximate 2 rms [ 0 exp(− 2 / 2 0 ) ( )] 2 and solve the equation of motion, eq. (12), perturbatively in the small parameter ≡ 1/ 0 . The first-order correction to the perpendicular momentum ⊥ is obtained by substituting into eq. (12) 0 = 0 and = , i.e. assuming the electron is undeflected. The deflection angle follows as ⊥ / 0 : The outer ring in fig. 7(a) and (b) corresponds to scattering at = 0 /2 (shown by the black, dashed line), at which eq. (20) is maximised, and the inner ring to scattering at = 0 (shown by the black, dotted line), which is the radius of the electron beam.
As discussed in section III, and shown in fig. 1, angular structure in the photons emerges differently in the LMA and LCFA simulations. In the former, it is the emission rate and the conservation of quasimomentum that ensures that photons are most probably emitted at angles 0 / 0 to the instantaneous quasimomentum. In the latter, it arises from the instantaneous Analytical predictions for the scattering angle are also given in [83], but these are derived under the assumptions that the laser transverse intensity profile is flat up to a radius equal to the waist, and that the pulse duration is infinitely long. Neither condition applies here. oscillation in the electron kinetic momentum, which has characteristic angle 0 / 0 , and the fact that the radiation is beamed parallel to this. The azimuthal symmetry of a circularly polarised laser means that the radiation angular profile is annular in shape: while this is evident in fig. 7(c) and (d), the characteristic angle is smaller than the expected value = 0 / 0 , which is shown by the black, dashed line. This is caused by the fact that the electrons are distributed over a range of impact parameters and therefore encounter lower effective values of 0 : eff 0 ( ) 0 exp(− 2 / 2 0 ).
Focal spot averaging not only lowers the yield of photons, as compared to a plane wave with the same peak amplitude, it also reduces the clarity of signatures of strong-field QED effects. We demonstrate this in particular for the position of the first nonlinear Compton edge, at 0 ∼ 1, This also provides an opportunity to crosscheck our LMA simulation results for focused lasers with theory. The latter is obtained using eq. (11), i.e. under the high-energy approximation that the electron is undeflected during its passage through the laser pulse. We have already shown that the total deflection angle scales as ( 0 / 0 ) 2 , which is indeed very small. In this case, the laser amplitude is either 0 = 0.5 or 2.5, its waist is 0 = 4 μm, and its temporal envelope (electric-field) and are distributed uniformly over a disk of radius 2 0 .
In fig. 8, we compare the theory and simulation results with those obtained for a plane wave with the same peak amplitude. As the total yield is reduced in the former case, we scale the plane-wave results by a factor 3D which approximately accounts for the effect of focal spot averaging. In the perturbative limit rms 1, the emission rate is proportional to 2 rms . Thus we expect the overall number of photons, in the 3D case, to be reduced by a factor 3D  fig. 8(a) and (c) for 0 = 0.5 and 2.5 respectively. Figure 8(b) and (d) show lineouts through the double-differential spectrum at fixed ⊥ = 0 /2. The agreement between theory and simulation is reasonably good. The detailed structure in the lineouts is not resolved, because the LMA misses interference effects at the scale of the pulse envelope. This is more evident in fig. 8(b) than (d), i.e. at lower 0 , for the following reason. In the LMA, the only contribution to the bandwidth of an individual harmonic is the variation in the mass shift over the pulse duration: at fixed ⊥ , this width is Δ / = 2 0 /(1 + 2 ). There is an additional contribution from the non-zero bandwidth of the pulse, which is given approximately by where is the FWHM duration of the pulse intensity profile: for the cosine-squared envelope under consideration here, 0 0.364 and Δ / 1/(2 ). At sufficiently small 0 , it is the latter contribution, from the laser pulse bandwidth, that dominates. Note that in a focussed pulse, the effective amplitude at finite impact parameter rms ( ) < 0 and so such effects are magnified.
Integrating over a finite range of ⊥ partially mitigates this, which is why the single-differential spectra are in much better agreement.
The difference between the 1D and 3D cases, evident in the theory, is captured very well by the simulations. We see that the first nonlinear edge is smeared out by focal spot averaging, particularly for 0 = 2.5. This is because the position of the edge differs for electrons at different impact parameters, as increasing means reducing the effective 0 . We have repeated the comparison between LMA-based simulations and QED for more tightly focussed laser pulses, reducing the waist 0 to 3 and 1 , while holding the peak 0 and the electron-beam-laser overlap fixed. The detailed results are shown in the Supplementary Material: we find that the spectra are barely affected by the reduction and the agreement between simulations and theory is consistently good.
This supports our expectations that LMA-based simulations are accurate even for focussed laser pulses.

VI. DISCUSSION
The focus of this paper has been incorporating the LMA into numerical simulation and providing the first benchmarks for a range of parameters and observables with direct calculation from QED.
As part of this work, we compare the LMA to the LCFA, which is the standard scheme for including QED effects in the modelling of intense laser interactions. The power of the LCFA is, in part, due to its versatility. It can be used when the strong electromagnetic field is not known a priori, which is a particular advantage when dealing with a laser-plasma interaction. However, in situations where the shape of the intense laser pulse is well-known, and unchanged in the interaction, the LMA can be used to attain a higher precision than the LCFA. The demand for the precision is acute if the field strength and particle energies in question are outside the region of validity of the LCFA, as is the case in some upcoming high-energy experiments [12,13].
Using a plane-wave pulse with phase duration 0 = 2 , we can give some indication for parameter regimes where the approximations can be used. We designate the region of accuracy as being when particle and field parameters are far away from values where it is known that the approximation is in doubt. We summarize our findings in table II for the photon yield and spectrum separately, as the conditions depend on the quantity that is to be measured. Note that, for these two quantities, the 'penalty' from violating the validity conditions is not equal in each case: if the LCFA is used to calculate the yield outside of the range given by 0 and , the prediction can be wrong by orders of magnitude (as demonstrated by fig. 4). The requirement that 0 1 for the LMA to be accurate is softer. Figure 4 shows that the simulations are accurate to within a few per cent even if = 4, which is already much shorter than most typical laser pulses (where the equivalent 10): furthermore we expect the accuracy to improve with increasing . Similarly, the error made by the LCFA in the photon spectrum becomes arbitrarily large as → 0, even if 0 1, whereas the LMA is guaranteed to obtain the correct limiting value for all 0 . The LMA result is inaccurate only in a small region around * 0 = / 1+ 2 0 + / [74], which is in any case missed by the LCFA: radiation here arises from the slow ponderomotive scattering of the electron.
If the plane wave contains a chirp, then the condition that the LMA is still a good approximation is found to be | / | 1, which is a 'slowly varying' condition of the same nature as 1/ 1. In practice, one is often limited to such chirps by virtue of the available bandwidth [75]. For focused laser pulses, we may expect that the applicability regime apparent when a chirp was introduced, should be adaptable for the change in wave vector ì obeying | ì |/| ì | 1. This reduces to a condition on the focusing, which may be expressed through the diffraction angle = /( 0 ) as 1. (A similar condition applies to the validity of the ponderomotive-force approach to the particle's classical dynamics [65], which confirm for 0.12 in fig. 8.) Using the high-energy approximation with the focused background in eq. (19), we note / = − / (2 ), and the relevant change in wave-vector on the focal axis is / ∼ 1/( ) = 2 /2 1. We have crosschecked our simulation results against a theoretical calculation which uses a WKB approximation, suitable for high-energy electrons. This integrates the plane-wave probability over a transverse distribution of intensity and particle flux and therefore works in a similar way to our simulations.
We find good agreement even for diffraction-limited focal spots (see Supplementary Material), when | ì |/| ì | 1 still holds. To benchmark the LMA fully would require solutions to the Dirac equation in backgrounds beyond a plane wave, which is still an active area of research [84][85][86].

VII. SUMMARY
Motivated by the imminent need for precision simulations of strong-field QED processes in the transition regime 0 ∼ 1, we have presented here a novel simulation framework which incorporates quantum effects via probability rates calculated within the locally monochromatic approximation (LMA) [37]. From the theory perspective, the formalisation of the LMA from the plane-wave model has been extended to include chirped pulses, under a "slowly varying chirp" approximation.
We have also adapted the LMA to model focused laser backgrounds, under the approximation that the incident electron has a relativistic factor satisfying 0 .
The emission rates so derived are embedded within a classical simulation framework that assumes a definite particle trajectory. In contrast to simulations based on the locally constant field approximation (LCFA), the electron quasimomentum (the cycle-averaged kinetic momentum) plays the essential role here, appearing in the classical equations of motion and the conservation of momentum. The fast oscillation of the particle momentum, at the timescale of the laser frequency, is nevertheless included, but at the level of the emission rates. This simulation framework therefore has conceptual similarities to the "envelope solvers" used to model laser-wakefield acceleration [87][88][89].
In benchmarking the simulations against QED results, we have found excellent agreement for a variety of background field configurations. Furthermore, we obtain significant reductions in the relative error when compared to the use of the LCFA in the transition regime. While we have focused, in this work, on the specific example of nonlinear Compton scattering in a circularly polarised background, our results can be extended to other processes, such as electron-positron pair creation [1,37], and to include spin-and polarisation-dependence [90][91][92][93][94].

ACKNOWLEDGMENTS
We would like to thank members of the LUXE collaboration for helpful discussions during Centre North.

DATA AVAILABILITY
The source code for the simulation program described is available at Ref. [95]. Version 0.6.0, which is used in this work, the input configurations necessary to reproduce the simulation results, and the analytical results, are archived at Ref. [96].
In [37], the LMA was derived from plane-wave QED for a simple plane-wave pulse. A plane wave is a highly idealised model of a laser field, which does not take into account some of the important characteristics of pulses in a real experiment. Here we extend the LMA to the case of a plane-wave pulse which includes an arbitrary chirp. We begin with a general overview of the LMA for a plane-wave field with a general chirp term.
For concreteness, we use a circularly polarised pulse with an arbitrary chirp, where the dimensionless gauge potential ( ) = ( )/ is and the phase is = · . In the derivation of the LMA, it is more natural to work with functions of the phase variable , than the proper time , which is used in the main text, and so in what follows we work with . The discussion here can be generalised to linearly or elliptically polarised backgrounds (see [37] for more details on the subtleties involved in the LMA for a linear, unchirped, plane-wave pulse).
We follow the standard approach of defining the scattering amplitude for our process in terms of the Volkov wavefunctions for the background dressed fermions of mass and 4-momentum , [56], where are constant spinors. The Volkov phase term is given by, which is just the classical action for an electron in a plane-wave background field. The nontrivial dependence of the Volkov wavefunctions on the phase means that overall momentum conservation for an arbitrary scattering amplitude S in the plane-wave background field only holds for three of the four directions, {−, ⊥}. As such, the scattering amplitude takes the form, where 3 −,⊥ ( ) = ( − ) ( 1 ) ( 2 ), and M is the invariant amplitude.
Closed form solutions to eq. (A3) are not always available. A simple example is the infinite monochromatic plane wave, which is the ( /Φ) → 1, ( ) → limit of the background field eq. (A1). However, one can separate the fast and slow dynamics of the background field in such a way that the field dependent terms in the exponent can by integrated by parts, and simplified by neglecting derivative corrections. This technique is known as the slowly varying envelope approximation [37][38][39][40][41].
The slowly varying envelope approximation for an arbitrarily chirped plane-wave field was derived in [59], and we follow this approach here. For the circularly polarised background eq. (A1), the terms which are quadratic in the field depend only on the slowly varying envelope, 2 ( ) = − 2 0 2 ( /Φ), while the terms linear in the field contain both slow (through ) and fast (through ) timescales. This gives integrals of the form, To deal with these integrals, we first transform the trigonometric functions of ( ) to pull out a factor depending on the inverse of ( ) = ( ), where a prime denotes a derivative of the argument: The function ( ) is taken to define a local frequency scale. Each term can then be readily integrated by parts, giving two contributions: a boundary term and a term proportional to d d Provided this is a small correction, which is valid for sufficiently long pulses, Φ 1 and when the derivative of the chirp function satisfies ( ) ( ), we can neglect these slowly varying terms, and approximate the integrals by, Applying these approximations to the classical action in eq. (A3) gives, The function ( ) contains only slowly varying terms, or terms linear in . The function ( ) depends on the phase only through the slowly varying envelope ( /Φ) and local frequency ( ), and the angle is independent of the phase.
The exponential of the trigonometric function in eq. (A9) can be expanded into an infinite sum of Bessel functions using the Jacob-Anger expansion, For the case of a one vertex process, such as nonlinear Compton scattering or Breit-Wheeler pair production, once the oscillating phase term has been expanded by eq. (A10), the invariant amplitude, M, in eq. (A4), takes on the form, The probability, P, is then found in the usual way by squaring the scattering amplitude eq. (A4) and integrating over the Lorentz invariant phase space for the particular process, dΩ LIPS , There are now two phase integrals, and what distinguishes the LMA from the slowly varying approximation (which is all we have applied so far) is performing a local expansion in the phase variables. To achieve this we introduce the sum and difference variables, and then take the small phase difference approximation 1 to expand the probability in a Taylor series in , retaining only the leading-order, ( ), contributions.
The -integral can be performed analytically, leaving the probability in the form, The function, R LMA ( ), contains summations over the Bessel harmonics and integrations over the final states, but crucially only depends on one phase variable. This allows us to interpret R( ) as a local rate which can be used in simulations. (In the main paper, we instead use a rate LMA defined as a probability per unit proper time.) To make this discussion more explicit, we consider the process of nonlinear Compton scattering.

Nonlinear Compton scattering in a chirped plane-wave pulse
Consider an electron with an initial momentum interacting with a plane-wave electromagnetic field to produce a photon of momentum and polarisation * , . The scattering amplitude, in terms of the Volkov wave functions eq. (A2), is given by, Here we use the Dirac slash notation, / = , where are the Dirac gamma matrices. The momentum is the momentum of the outgoing electron.
Performing all of the trivial integrations to express the scattering amplitude in the form eq. (A4), the invariant amplitude is found to be, where the spin dependent structure is given by, The function ( ) has the explicit form, where we have defined the lightfront momentum fraction = · / · . As stated above, this only has dependence on the phase through either linear or slowly varying terms.
The term ( ) is and so the only dependence on the phase comes through the ratio of the slowly varying pulse envelope and the local frequency. The angle is defined through the relationship, and so can be interpreted as the angle between the components of the 4-vector − projected onto the directions of background field polarisation.
We skip now to the explicit form of the probability. Expanding into Bessel harmonics according to eq. (A10), the probability eq. (A12) becomes (A23) The probability in this form contains two infinite sums over the Bessel harmonics and integrals over the outgoing photon momentum. Note the exponential dependence on the chirp function, ( ), and the angle . If we consider the definitions eq. (A20)-eq. (A22), we notice that the only dependence on the transverse photon momentum is through the combination ì ⊥ = ì ⊥ /( )− ì ⊥ / .
We can then shift the integration variables in eq. (A23), and using eq. (A22) express the integration measure in polar coordinates, The only dependence of the probability on the angle is then through the exponential factor exp(+ ( − ) ). The integration over the angle sets = . This allows the probability to be well approximated by, Following through with the local expansion, using eq. (A13) and 1, the integral over d can be performed, which gives a -function: where we have defined 0 = · / 2 . The probability only has support when the argument of the -function satisfies: which (upon adapting the notation) is found to be exactly the stationary phase condition which is evaluated in [59] (see eq. (25) of [59]). In that work, the stationary phase approximation is carried out at the level of the amplitude for nonlinear Compton scattering in the slowly varying envelope approximation. Here we have shown that the exact same kinematic relationship reappears at the probability level after the explicit application of a local expansion.
The integral over the remaining perpendicular momentum dependence can be trivially carried out using the -function in eq. (A26), which gives the relatively concise expression (suppressing explicit dependence on ) where the argument of the Bessel functions is now and we have defined the cycle-averaged potential rms = 0 ( /Φ) and the upper bound on the integration over is (A30) Thus, when compared with the expressions found for the LMA in a non-chirped pulse [37], the chirp function, ( ), contributes an effective rescaling of the lightfront energy parameter, 0 → 0 ( ), inside the argument of the Bessel functions. In eq. (10) we have redefined and , * by absorbing