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

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, a₀, and quantum nonlinearity parameter of the electron, χ_e, exceed unity. Perturbative QED calculations of the interaction fail once a₀ ≪̸ 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, a₀ ≃ 0.4, laser-electron collision experiments have observed Compton scattering (photon emission) and trident electron–positron pair creation [4, 5]. At higher values of a₀, but small χ_e, they have observed photon emission in the classical regime (nonlinear Thomson scattering) [6–9] and at a₀ ≃ 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 ≲ a₀ ≲ 10 at χ_e ≃ 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 a₀ ≪̸ 1. However, neither is an approach based on the locally constant field approximation (LCFA) [1, 14], as this applies only in the opposite limit, a₀ ≫ 1: this approximation underpins the simulation codes [15–17] used to model QED effects in laser–plasma interactions [18–24], which will be explored in the next generation of multi-petawatt laser facilities [25–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 a₀ ≫ 1 and ${a}_{0}^{3}/{\chi }_{\mathrm{e}}\gg 1$; the limitations of doing so have been thoroughly discussed in the literature [29–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 modeling of plane-wave-like fields. While plane-wave rates have already been used in numerical modeling and analysis [13, 34–36], their derivation from strong-field QED has only recently been formalised by Heinzl et al [37], who combine a slowly varying envelope approximation [38–41] with a 'local' expansion in the interference phase [1, 3, 29, 31, 32, 42].

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-plane-wave 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–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 modeling 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: ℏ = c = ₀ = 1. The electron mass is denoted by m. The fine-structure constant α is related to the elementary charge e by α = e²/(4π).

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. [40, 48–52].) 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 A. We define the dimensionless potential a = eA/m,

$$\begin{equation}a(\phi )={a}_{0}f \left(\frac{\phi }{{\Phi}}\right)\left[\varepsilon \enspace \mathrm{cos}\enspace b(\phi )+\beta \enspace \mathrm{sin}\enspace b(\phi )\right],\end{equation} \tag{ 1 }$$

where a₀ is the dimensionless intensity parameter [53] (also called the 'classical nonlinearity', normalised amplitude or the strength parameter) and ɛ, β are orthonormal polarisation vectors obeying ɛ ⋅ɛ = β ⋅ β = −1. Throughout, we use lightfront coordinates ${x}^{\mu }={({x}^{+},{x}^{-},{ \overrightarrow {x}}^{\perp })}^{\mu }$, where x^± = x⁰ ± x³, ${ \overrightarrow {x}}^{\perp }=({x}^{1},{x}^{2})$, x^± = 2x_∓ and ${ \overrightarrow {x}}^{\perp }=-{ \overrightarrow {x}}_{\perp }$. The function f(ϕ/Φ) is the pulse envelope which depends on the lightfront phase ϕ = κ ⋅ x (where ${\kappa }_{\mu }={\delta }_{\mu }^{+}{\kappa }_{+}$ is the background wavevector), and the pulse phase duration, Φ, is related to the number of cycles, N, via Φ = 2N. The function b(ϕ) describes the chirp of the background. For a pulse without chirp, b is linear in ϕ, i.e. b''(ϕ) = 0 for all ϕ. (In the following, we will pick b(ϕ) = ϕ for the unchirped case.)

We use the scattering matrix approach [54] 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:

$$\begin{equation}{\mathsf{S}}_{{r}^{\prime },r;l}=-ie\int \enspace {\mathrm{d}}^{4}x\enspace {\overline{{\Psi}}}_{{p}^{\prime },{r}^{\prime }}(x){\hspace{2pt}{\epsilon}\hspace{-4pt}/\hspace{2pt}}_{k,l}^{{\ast}}{\text{e}}^{\text{i}k\cdot x}{{\Psi}}_{p,r}(x),\end{equation} \tag{ 2 }$$

where ${\overline{){\epsilon}}}_{k,l}^{{\ast}}$ is the polarisation of the emitted photon with four-momentum k and ${{\Psi}}_{p,r}({\overline{{\Psi}}}_{{p}^{\prime },{r}^{\prime }})$ is the Volkov wavefunction [55] of the incoming (outgoing) electron:

$$\begin{equation}{{\Psi}}_{p,r}(x)=\left(1+\frac{m\overline{)\kappa }\overline{)a}}{2\enspace \kappa \cdot p}\right){u}_{p,r}\enspace {\text{e}}^{\mathrm{i}{S}_{p}(x)},\quad {S}_{p}(x)=p\cdot x+{\int }^{\phi } \mathrm{d}y\enspace \frac{2mp\cdot a(y)-{m}^{2}{a}^{2}(y)}{2\kappa \cdot p}.\end{equation} \tag{ 3 }$$

The matrix element can be simplified to:

$$\begin{equation}{\mathsf{S}}_{{r}^{\prime },r;l}=~{C}{\int }_{{\phi }_{i}}^{{\phi }_{f}} \mathrm{d}\phi \enspace {\overline{u}}_{{r}^{\prime }}\left[{\Delta}{\overline{){\epsilon}}}_{k,l}^{{\ast}}+\frac{m}{2\enspace \kappa \cdot p}\left(\frac{\overline{)a}\overline{)\kappa }{\overline{){\epsilon}}}_{k,l}^{{\ast}}}{1-s}+{\overline{){\epsilon}}}_{k,l}^{{\ast}}\overline{)a}\overline{)\kappa }\right)\right]{u}_{r}\enspace \mathrm{exp} \left[\frac{i}{{\eta }_{0}(1-s)}{\int }_{{\phi }_{i}}^{\phi } \mathrm{d}y\enspace \frac{k\cdot \pi (y)}{{m}^{2}}\right]\end{equation} \tag{ 4 }$$

where s = κ ⋅ k/κ ⋅ p is the lightfront momentum fraction of the emitted photon, η₀ = κ ⋅ p/m² is the initial energy parameter of the probe electron, $~{C}$ contains normalisation constants, the instantaneous electron momentum is given by

$$\begin{equation}\pi (y)=p-ma(y)+\kappa \frac{2\enspace mp\cdot a(y)-{m}^{2}{a}^{2}(y)}{2\enspace \kappa \cdot p},\end{equation} \tag{ 5 }$$

and the regularising factor Δ = 1 − k ⋅ π/k ⋅ p incorporates all the contributions from phases outside of the integral. The total probability can be written:

$$\begin{equation}\mathsf{P}=\frac{\alpha }{{\eta }_{0}}\frac{1}{{2}^{4}{\pi }^{2}}\int {\mathrm{d}}^{2}{ \overrightarrow {r}}^{\perp }\mathrm{d}s\enspace \frac{s}{1-s}{\langle \vert {\mathsf{S}}_{{r}^{\prime },r;l}{\vert }^{2}\rangle }_{\mathsf{p}\mathsf{o}\mathsf{l}.},\end{equation} \tag{ 6 }$$

where ${ \overrightarrow {r}}^{\perp }={ \overrightarrow {k}}^{\perp }/(ms)-{ \overrightarrow {p}}^{\perp }/m$ contains the shifted perpendicular momentum. Here '^⊥' indicates directions perpendicular to the background propagation direction and ${\left\langle \cdot \right\rangle }_{\mathsf{p}\mathsf{o}\mathsf{l}.}$ indicates an average over initial and sum over final polarisation states. The numerical results in exact QED are calculated by evaluating equation (6) directly: the matrix element in equation (2) was evaluated using photon polarisation eigenstates of the background [56] and spin states in the Lepage–Brodsky convention [57].

Rather than direct numerical evaluation, some of the integrals in equation (6) can be evaluated analytically by generalising the LMA [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 equation (1). For the LMA to approximate the emission spectrum well, the envelope function f(ϕ/Φ) should be slowly varying with respect to the carrier frequency, implying that ${{\Phi}}^{-1}\ll \mathrm{min} \left[{b}^{\prime }(\phi )\right]$ (i.e. Φ ≫ 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 [58]). These approximations then allow the squared Kibble mass, μ, which occurs in an exponent, to be integrated over. The Kibble mass takes the form $\mu =1+{\left\langle { \overrightarrow {a}}^{2}\right\rangle }_{\theta }-{\left\langle \overrightarrow {a}\right\rangle }_{\theta }^{2}$, where ${\left\langle f\right\rangle }_{\theta }={\theta }^{-1}{\int }_{\phi -\theta /2}^{\phi +\theta /2}\hspace{2pt}f$ denotes a phase-window average. In the case of a circularly polarised background, the slowly varying (envelope) and rapid (carrier) timescales occur in ${\left\langle \overrightarrow {a}\right\rangle }_{\theta }$. We can demonstrate the approximation by considering a single component of $\overrightarrow {a}$, e.g. $\overrightarrow {\varepsilon }\cdot \overrightarrow {a}$.

$$\begin{equation}{\left\langle \overrightarrow {\varepsilon }\cdot \overrightarrow {a}\right\rangle }_{\theta }=\frac{{a}_{0}}{\theta }{\int }_{\phi -\theta /2}^{\phi +\theta /2} \mathrm{d}\phi \enspace f \left(\frac{\phi }{{\Phi}}\right)\mathrm{cos}\enspace b(\phi ).\end{equation} \tag{ 7 }$$

Now, one can introduce a local frequency scale, ω(φ) = b'(φ) and integrate by parts as in equation (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

$$\begin{equation}\sim \frac{1}{{\Phi}}\frac{{f}^{\prime }(\phi /{\Phi})}{f(\phi /{\Phi})},\qquad \sim \frac{{\omega }^{\prime }(\phi )}{\omega (\phi )},\end{equation} \tag{ 8 }$$

respectively (neglecting a rapidly varying term that appears ∼cot b(ϕ)). 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 $\overrightarrow {\beta }\cdot \overrightarrow {a}$ term, whereas ${\left\langle { \overrightarrow {a}}^{2}\right\rangle }_{\theta }$ 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 A).

Finally, we arrive at ${P}^{\text{LMA}}=\int \mathrm{d}\tau \enspace {W}^{\text{LMA}}$, where:

$$\begin{equation}{W}^{\text{LMA}}=\sum\limits _{n=1}^{\infty }{\int }_{0}^{{s}_{n,{\ast}}(\tau )} \mathrm{d}s\enspace \frac{{\mathrm{d}}^{2}{P}_{n}^{\text{mono}}\left[{a}_{\text{rms}}(\tau ),\eta (\tau )\right]}{\mathrm{d}\tau \enspace \mathrm{d}s},\end{equation} \tag{ 9 }$$

where ${a}_{\text{rms}}^{2}(\tau )={q}^{2}/{m}^{2}-1$ and η(τ) = ω[ϕ(τ)]η₀, with η₀ = κ ⋅ p/m² the unchirped energy parameter. Here $q=\left\langle \pi \right\rangle$ is the quasimomentum, the laser-cycle-average of the instantaneous electron momentum given in equation (5). The appearance of a local wavevector in η(τ) also follows from considering components of the field-strength tensor, F^μν , for the chirped pulse in equation (1), which contain terms ∼κ^μ (ϕ)∂a^ν /∂b, where κ^μ (ϕ) = b'(ϕ)κ^μ . ${\mathsf{P}}_{n}^{\text{mono}}$ is the probability of nonlinear Compton scattering into the nth harmonic in a monochromatic background, τ is the proper time, related to the phase by dτ/dϕ = 1/(mη₀). The approximation is locally monochromatic because the intensity and energy parameter occurring in the monochromatic probability now take the (cycle-averaged) local value at the position of the electron. The integrand is given explicitly by equation (A28) for nonlinear Compton scattering. Unlike the monochromatic case, here the harmonic range is phase-dependent:

$$\begin{equation}{s}_{n,{\ast}}(\tau )=\frac{{s}_{n}(\tau )}{1+{s}_{n}(\tau )},\qquad {s}_{n}(\tau )=\frac{2n\eta (\tau )}{1+{a}_{\text{rms}}^{2}(\tau )},\end{equation} \tag{ 10 }$$

where s_n (τ) 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 γ⁻¹, where γ is the initial relativistic gamma factor of the incident electron [45–47]. Then assuming γ ≫ a₀, for a head-on collision of the electron probe with the focused laser pulse, one can write:

$$\begin{equation}{P}^{2\text{D}}=\int {\mathrm{d}}^{2}{\to {x}}^{\perp }\enspace \rho ({\to {x}}^{\perp })P[{a}_{\text{rms}}({\to {x}}^{\perp }),\eta ({\to {x}}^{\perp })],\end{equation} \tag{ 11 }$$

where ρ is the electron probe areal density and the plane-wave probability, $\mathsf{P}$ from equation (6), now has an intensity parameter which can depend on the perpendicular spatial co-ordinate.

The inclusion of strong-field QED processes in numerical simulations, such as the particle-in-cell [15, 16] or particle-tracking codes [34, 35, 59] 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 [60]. This is motivated by the appearance of the classical kinetic momentum π, equation (5), in the QED scattering probability, via the exponent of the Volkov wavefunction, equation (3). (This occurs because the Volkov solution is identical to the semiclassical solution of the Dirac equation in a plane-wave background.) This permits the probability, equation (6), to be approximated as the integral $\mathsf{P}\simeq \int W\enspace \mathrm{d}\tau$, where W ⩾ 0 is interpreted as a probability rate, which can depend, inter alia, on the instantaneous momentum and field amplitude.

The approximations applied to the probability rate affect what dynamical quantities must be obtained from the classical trajectory. In the LCFA, for example, the rate W = W[χ(τ)], where the quantum nonlinearity parameter $\chi (\tau )=e\left\vert {F}_{\mu \nu }[x(\tau )]{\pi }^{\nu }(\tau )\right\vert /{m}^{3}$ [1]. Furthermore, the conservation of momentum for the scattering may be written such that it constrains the kinetic, rather than asymptotic, momenta. Thus the classical trajectory must be defined in terms of kinetic momentum π, i.e. instantaneously, and obtained from the Lorentz force equation dπ_μ /dτ = −eF_μν π^ν /m and dx^μ /dτ = π^μ /m. This is illustrated on the left-hand side of figure 1: the classical trajectory is well-defined at all timescales, including that of the laser carrier wave. The angular structure of the photon emission arises from two sources: the oscillation of the trajectory (θ ≃ a₀/γ for γ ≫ a₀ ≫ 1) and the intrinsic beaming of the emission around the instantaneous velocity, the latter being of characteristic size θ ∼ 1/γ [61, 62]. The former is the dominant contributor in the regime a₀ ≫ 1, which is consequently where the LCFA is expected to be valid.

Zoom In Zoom Out Reset image size

The rate in the LMA, 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 $q\equiv \left\langle \pi \right\rangle$ [1, 63]. In a plane wave, π = p − ma + κ(2 mp ⋅ a − m² a²)/(2κ ⋅ p) and π² = m², whereas $q=p+\kappa {m}^{2}{a}_{\text{rms}}^{2}/(2\kappa \cdot p)$ and ${q}^{2}={m}^{2}(1+{a}_{\text{rms}}^{2})$, for ${a}_{\text{rms}}^{2}\equiv -\left\langle {a}^{2}\right\rangle$. In contrast to the LCFA case, the rate is a function of two parameters: the normalised amplitude (or intensity parameter), a_rms, and the energy parameter η ≡ κ ⋅ p/m², both locally defined. (The root-mean-square quantum parameter follows as χ_rms = a_rms η.) Both may be obtained from q as follows: ${a}_{\text{rms}}=\sqrt{{(q/m)}^{2}-1}$ and η = κ ⋅ q/m². An equation of motion for the quasimomentum may be obtained by separating the Lorentz force equation (in a focused, pulsed electromagnetic wave) into quickly and slowly varying components and isolating the latter. The result is the relativistic ponderomotive force equation [64]:

$$\begin{equation}\frac{\mathrm{d}\to {q}}{\mathrm{d}t}=-\frac{{m}^{2}}{2{q}^{0}}\frac{\partial {a}_{\text{rms}}^{2}}{\partial \to {r}},\end{equation} \tag{ 12 }$$

where ${q}^{0}={[{m}^{2}(1+{a}_{\text{rms}}^{2})+{\left\vert \overrightarrow {q}\right\vert }^{2}]}^{1/2}$. The slowly varying components of the position are determined by

$$\begin{equation}\frac{\mathrm{d} \overrightarrow {r}}{\mathrm{d}t}=\frac{ \overrightarrow {q}}{{q}^{0}}.\end{equation} \tag{ 13 }$$

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 figure 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 s, in the limit a₀ ≫ 1, there is a most probable harmonic index $n={a}_{\text{rms}}^{2}s/[\eta (1-s)]$ [65]. Combining this relation with the conservation of quasimomentum, which reads ${k}_{\perp }^{2}/{m}^{2}=2n\eta s(1-s)-{s}^{2}(1+{a}_{\text{rms}}^{2})$ for p_⊥ = 0, one finds that the most probable emission angle is θ ≃ a_rms/γ for γ ≫ a₀ ≫ 1 [65] (see also [63]). 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 1.

Table 1. Overview of the conceptual differences between LCFA- and LMA-based simulations of photon emission in strong laser pulses.

	LCFA	LMA
Rate derived for	Constant, crossed field	Monochromatic plane wave
		(fast quiver motion here)
and controlled by	Instantaneous momentum πμ	Quasimomentum ${q}^{\mu }=\left\langle {\pi }^{\mu }\right\rangle$
	via quantum parameter χe	via arms and η
	Lorentz force:	Ponderomotive force:
equation of motion	$\frac{\mathrm{d}{\pi }_{\mu }}{\mathrm{d}\tau }=-\frac{e{F}_{\mu \nu }{\pi }^{\nu }}{m}$	$\frac{\mathrm{d} \overrightarrow {q}}{\mathrm{d}t}=-\frac{{m}^{2}}{2{q}^{0}}\frac{\partial {a}_{\text{rms}}^{2}}{\partial \overrightarrow {r}}$
	(fast quiver motion here)

The emission of photons, and its effect on this trajectory, is modeled in the following way. At any particular timestep, we have the electron quasimomentum q and position r from the classical equations of motion, as well as the local values of the laser normalised amplitude a_rms(r), wavevector κ(r) and polarisation (taken to be circular throughout). In fact, κ and q are sufficient to determine the properties of the emission, as they define the two invariant parameters, a_rms and η, that control the rate and the conservation of momentum. This is given by

$$\begin{equation}q+n\kappa ={q}^{\prime }+k,\end{equation} \tag{ 14 }$$

where q' is the electron quasimomentum after the scattering, k is the momentum of the emitted photon, and n is the harmonic index (the net number of laser photons absorbed). The emission rates themselves control n and subsequently s ≡ κ ⋅ k/κ ⋅ q, the lightfront momentum fraction. Given n, s and q, it is a matter of kinematics to determine k and then q'. Our Monte Carlo algorithm is as follows: (i) advance the electron trajectory by solving equations (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 n with probability W_n /W, where W_n is the partial rate and $W={\sum }_{n=1}^{\infty }{W}_{n}$ is the total rate, (iv) sample s from the partial spectrum (dW_n /ds)/W_n , (v) determine k given n, s and q and (vi) reset the electron quasimomentum from q to q'.

The probability that emission takes place in small interval of lab time Δt is given by $\mathsf{P}=W{\Delta}\tau$ and Δτ = Δt(m/q⁰) is the equivalent interval of proper time. We obtain W by integrating, and then summing, the partial, differential rates of emission W_n , which are given by [37]

$$\begin{equation}\frac{\mathrm{d}{W}_{n}}{\mathrm{d}s}=-\alpha m\left\{{J}_{n}^{2}(z)+\frac{{a}_{\text{rms}}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]\left[2{J}_{n}^{2}(z)-{J}_{n-1}^{2}(z)-{J}_{n+1}^{2}(z)\right]\right\}.\end{equation} \tag{ 15 }$$

The argument z of the Bessel functions J_n (of the first kind [66]) and auxiliary variables are

$$\begin{equation}{z}^{2}=\frac{4{n}^{2}{a}_{\text{rms}}^{2}}{1+{a}_{\text{rms}}^{2}}\frac{s}{{s}_{n}(1-s)}\left[1-\frac{s}{{s}_{n}(1-s)}\right],\qquad {s}_{n}=\frac{2n\eta }{1+{a}_{\text{rms}}^{2}}\end{equation} \tag{ 16 }$$

and the bounds on s are 0 < s < s_n /(1 + s_n ). Note that s_n depends on a_rms and η and is therefore a function of proper time τ, as shown explicitly in equation (10). While the summation should run from n = 1 to infinity, it is sufficient to sum up to a largest value ${n}_{\mathrm{max}}=10(1+{a}_{\text{rms}}^{3})$. In principle, the integration and summation can be done at every timestep, given the particular values of a_rms and η. However, it is significantly faster to obtain W by interpolating from a lookup table, where W(a_rms, η) is precalculated over the domain ${a}_{\text{rms}}^{\mathrm{min}}{< }{a}_{\text{rms}}{< }{a}_{\text{rms}}^{\mathrm{max}}$ and η_min < η < η_max. The upper bounds are fixed by the problem space under consideration; we have taken ${a}_{\text{rms}}^{\mathrm{max}}=10$ and η_max = 2 in our code. The lower bounds are chosen such that alternative sampling strategies may be used.

First, if ${a}_{\text{rms}}{< }{a}_{\text{rms}}^{\mathrm{min}}\enspace \ll 1$, only the first harmonic, n = 1, contributes significantly to the probability. In this limit, the rate may be obtained analytically:

$$\begin{equation}W\simeq {W}_{1}+O({a}_{\text{rms}}^{4}),\qquad {W}_{1}=\frac{\alpha m{a}_{\text{rms}}^{2}}{2\eta }\left[\frac{2+8\eta +9{\eta }^{2}+{\eta }^{3}}{{(1+2\eta )}^{2}}-\frac{2+2\eta -{\eta }^{2}}{2\eta }\enspace \mathrm{ln}(1+2\eta )\right].\end{equation} \tag{ 17 }$$

Second, if η < η_min ≪ 1, we may take the classical limit, whereupon the partial rates become:

$$\begin{equation}\begin{aligned}\frac{\mathrm{d}{W}_{n}}{\mathrm{d}v}& \simeq \frac{\alpha mn\eta }{1+{a}_{\text{rms}}^{2}}[{a}_{\text{rms}}^{2}{J}_{n-1}^{2}(z)+{a}_{\text{rms}}^{2}{J}_{n+1}^{2}(z)-2(1+{a}_{\text{rms}}^{2}){J}_{n}^{2}(z)]+O({\eta }^{2}),\\ {z}^{2}& =\frac{4{a}_{\text{rms}}^{2}{n}^{2}v(1-v)}{1+{a}_{\text{rms}}^{2}},\end{aligned}\end{equation} \tag{ 18 }$$

but where we fix v = s(1 + s_n )/s_n to be 0 < v < 1. Equation (18), integrated over 0 < v < 1 and summed over n = 1 to n_max, is tabulated over the same range ${a}_{\text{rms}}^{\mathrm{min}}{< }{a}_{\text{rms}}{< }{a}_{\text{rms}}^{\mathrm{max}}$. In our implementation, ${a}_{\text{rms}}^{\mathrm{min}}=0.02$ and η_min = 10⁻³. Thus at every timestep, the emission probability $\mathsf{P}=W{\Delta}\tau$ is obtained by interpolating from the appropriate lookup table, or using the limiting analytical expression. Emission is deemed to occur if a pseudorandom number R, drawn from the uniform distribution U(0, 1), satisfies $R{< }\mathsf{P}$.

If emission takes place, the next step is to determine n and s. The former is obtained by solving for n, ${R}^{\prime }={\sum }_{i=1}^{n}{W}_{i}/W$, where R' is another pseudorandom number drawn on the unit interval U(0, 1). In our implementation, the total rate of emission W is already available at this point; however, the sequence of partial rates must be evaluated explicitly, by integrating equation (15) over s. We do this, rather than store a lookup table in n (as well as in a_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 n is fixed, the lightfront momentum fraction transferred, s, is obtained by rejection sampling of equation (15).

The kinematical calculation of k is performed in the zero momentum frame (ZMF), which moves with four-velocity $u=(q+n\kappa )/[m\sqrt{1+{a}_{\text{rms}}^{2}+2n\eta }]$ with respect to the lab frame. In the ZMF, the emitted photon has momentum $\left\vert { \overrightarrow {k}}_{\text{zmf}}^{\prime }\right\vert =mn\eta /\sqrt{1+{a}_{\text{rms}}^{2}+2n\eta }$ and polar scattering angle $\mathrm{cos}\enspace {\theta }_{\text{zmf}}=1-s(1+{a}_{\text{rms}}^{2}+2n\eta )/(n\eta )$. The azimuthal angle φ_zmf, which is arbitrary for circularly polarised backgrounds, is pseudorandomly determined in 0 ⩽ φ_zmf < 2π. Once ${ \overrightarrow {k}}_{\text{zmf}}$ is determined, it may be boosted back to the lab frame, where q' follows from equation (14).

While LMA rates have already been implemented in simulation codes used to study laser–electron interactions [34–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, 67] and proposed improvements in [31–33]. In this section, we compare the results of simulations based on the LMA, as outlined in section 3, with QED theory calculations without the LMA, for photon emission in a pulsed, plane-wave background. We focus on the transition regime a₀ ∼ 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 $\overrightarrow {a}(\phi )={a}_{0}f(\phi )[ \overrightarrow {x}\enspace \mathrm{cos}\enspace b(\phi )+ \overrightarrow {y}\mathrm{sin}\enspace b(\phi )]$, where f(ϕ) =  cos²[ϕ/(2N)] for $\left\vert \phi \right\vert {< }\pi N$. Here N 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 T[fs] ≃ Nλ[μm]/0.8. The function b(ϕ) controls the frequency chirping of the pulse and is initially set to b(ϕ) = ϕ (i.e. unchirped) for the results in section IVA. The electrons counterpropagate head-on to the laser pulse, with initial energy parameter η₀ = 0.1. This is equivalent to an initial Lorentz factor of γ₀ = 1.638 × 10⁴ for a laser wavelength of 0.8 μm.

The theoretical calculations described in section 2 are for single emission only. However, for sufficiently large a₀ or pulse length N, it is possible for the total probability of emission $\mathsf{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 [68–70]. Neglecting other contributions, where the intermediate state does not propagate, is expected to be a good approximation if ${a}_{0}^{2}\enspace {\Delta}\phi \gg 1$ [71], where Δϕ = 2πN is the phase duration of the pulse, which allows simulations to model cascades of photon emission and pair creation [60]. 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 [72]. 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, s, and normalised perpendicular momentum r_⊥ = k_⊥/(ms). These provide proxies for the emitted photon energy ω' and polar scattering angle θ, respectively: s = ω'(1 + cos θ)/p⁻ ≃ ω'/(mγ₀) and r_⊥ = (p⁻/m)tan(θ/2) ≃ γ₀ θ, where p⁻ = m² η₀/ω₀ is the initial lightfront momentum of the electron and γ₀ its Lorentz factor.

4.1. Pulsed plane waves

Figures 2(a) and (e) show photon spectra, double-differential in s and r_⊥, obtained from simulations in the linear and nonlinear regimes (a₀ = 0.5 and 2.5 respectively) for a pulse that is N = 16 cycles in duration. In the former case, radiation emission is dominated by the first harmonic, which displays the expected, characteristic energy–angle correlation. In the latter case, the radiation is composed of a broad range of high harmonics, extending the spectrum to much larger s. The effect of the pulse envelope is evident in the broadening of the first harmonic for small r_⊥: recall that the position of the first Compton edge, ${s}_{1}^{{\ast}}=2\eta /(1+{a}_{\text{rms}}^{2}+2\eta )$, is phase-dependent through a_rms and η. We also see that the higher harmonics are predominantly emitted at r_⊥≃ a₀, as expected in the nonlinear regime, whereas for a₀ = 0.5, the characteristic r_⊥ < a₀.

Zoom In Zoom Out Reset image size

The three plots accompanying each double-differential spectrum compare lineouts at fixed r_⊥ against theoretical results. The simulations capture the position and overall shape of the harmonics well, but miss the subharmonic substructure visible in figures 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 figure 2 integrated over r_⊥, are shown in figure 3. We compare the simulation results with QED for normalised amplitudes a₀ = 0.5 and 2.5 and for pulse durations equivalent to N = 4 and 16 cycles. The agreement is much better for the longer pulse, which we expect because the LMA neglects terms of order 1/N (see equation (8) and [37]). The LMA simulations capture the harmonic structure and correctly reproduce the small-s behavior of the theory, where the spectrum tends to a constant value $\propto {a}_{0}^{2}\int {f}^{2}(\phi )\enspace \mathrm{d}\phi$ [31, 37]. The LCFA simulations are significantly wrong in this region $s{< }{s}_{1}^{{\ast}}$, where we see the characteristic divergence ∝s^−2/3 [1].

Zoom In Zoom Out Reset image size

The intermediate structure, which appears below the first Compton edge for a₀ = 2.5, shown in figure 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 [73]. 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 toward smaller s as N increases and it is eventually lost in the monochromatic limit [37]. Integrating over the s-weighted probability, shown in figures 3(c) and (e), yields the total lightfront momentum transfer from electron to photon. If a₀ > 1, this is dominated by contributions from $s{ >}{s}_{1}^{{\ast}}$, where the LCFA works well [30]. However, it is evident from figure 3(c) that the LCFA fails globally for a₀ < 1.

Finally, we consider the total probability that a photon is emitted, $\mathsf{P}$, and the average lightfront momentum fraction of that photon, $\left\langle s\right\rangle \equiv \int s\frac{\mathrm{d}\mathsf{P}}{\mathrm{d}s}\enspace \mathrm{d}s$, as a function of a₀ for a four-cycle pulse. The values obtained from theory and from LMA and LCFA simulations are shown in figure 4, along with the percentage error made by the simulations. The LMA-based simulations are accurate at the level of a few per cent across the full range of a₀ explored. The improvement over the LCFA is particularly dramatic for the probability, where the error made is larger than 10% even when a₀ = 5. The average lightfront momentum fraction is more sensitive to the contribution of higher harmonics, i.e. large s; as this is where the LCFA works rather well, the accuracy for $\left\langle s\right\rangle$ is better than that for $\mathsf{P}$. However, the LMA simulations are significantly more accurate when a₀ ≲ 1.

Zoom In Zoom Out Reset image size

4.2. 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 a₀ ∼ 1. In this section we benchmark our simulation results against theory for this case.

The first example we consider is that of a linearly chirped laser pulse, which has potential $\overrightarrow {a}(\phi )={a}_{0}f(\phi )[ \overrightarrow {x}\enspace \mathrm{cos}\enspace b(\phi )+ \overrightarrow {y}\enspace \mathrm{sin}\enspace b(\phi )]$, where f(ϕ) =  cos²[ϕ/(2N)] for $\left\vert \phi \right\vert {< }\pi N$ and b(ϕ) = ϕ[1 + cϕ/(2 N)]. The instantaneous frequency, ω(ϕ) = ω₀(1 + cϕ/N) for chirp parameter c, must be positive throughout the pulse, which imposes the restriction c < 1/π. This is consistent with the condition for the chirp to be slowly varying, equation (8), which may be cast as c ≪ N/(1 + πN). Furthermore, for a particular pulse duration, there is an upper bound on the largest chirp that can be obtained [74]. In our notation, this maximum is given by χrp_max ≃ 10/N. 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 a₀ = 0.5, N = 16 and c = 1/(2π) in figure 5. The unchirped results, c = 0, are also shown for reference. The theoretical results are obtained numerically, using equation (6) and the explicit form of the potential $\overrightarrow {a}(\phi )$. For this case, the electron trajectory can be written in a closed form in terms of Fresnel functions. In the simulations, a chirp is included by promoting the frequency of the background κ^μ to be a function of phase κ^μ (ϕ). We find that the simulations capture the softening of the harmonic structure evident in the theory results for the chirped pulse. Lineouts through the theoretical double-differential spectrum at fixed r_⊥ demonstrate that chirping smooths out the subharmonic structure; as a consequence, simulation results appear to be more accurate than in the unchirped case.

Zoom In Zoom Out Reset image size

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 a₀ > 1. In a pulsed plane wave, the position of the first harmonic edge varies from s = 2η₀/(1 + 2η₀) to $s=2{\eta }_{0}/(1+{a}_{0}^{2}+2{\eta }_{0})$ as the cycle-averaged potential a_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 a₀ is not small, it has been proposed to chirp the pulse in a particular way [75–79]. If the instantaneous frequency of the pulse varies as $\omega (\phi )={\omega }_{0}[1+{a}_{\text{rms}}^{2}(\phi )]$, then s = 2η₀/(1 + 2η₀) for all ϕ and the nonlinear redshift is perfectly compensated. Although there are significant obstacles to achieving this in experiment, it is a useful test case for the simulation method we have introduced. We therefore consider a pulse with envelope f(ϕ) =  cos²[ϕ/(2N)] for $\left\vert \phi \right\vert {< }\pi N$ and $b(\phi )=\phi +{a}_{0}^{2}{\int }_{0}^{\phi }{f}^{2}(y)\enspace \mathrm{d}y$. In this case, the chirp may be considered to be slowly varying if $2{a}_{0}^{2}/[N(1+{a}_{0}^{2})]\ll 1$. We show results for a₀ = 1, N = 16 in figure 6. The lightfront momentum spectrum for theory and simulation both show a shift of the edge of the first harmonic from the nonlinear, to the linear position, as expected for this choice of chirp. However, this rather extreme choice of chirp leads to a larger discrepancy in the in the height of the spectra: the simulations underestimate the total yield by a small but not insignificant amount. We have verified that both theory curves tend to the same value in the limit of vanishing s, and that the simulation curves do as well: the limiting value, ${\mathrm{lim}}_{s\to 0}\enspace \frac{\mathrm{d}\mathsf{P}}{\mathrm{d}s}\propto {a}_{0}^{2}\int {f}^{2}(\phi )\enspace \mathrm{d}\phi$, is sensitive only to the pulse envelope (for circular polarization) [31, 37].

Zoom In Zoom Out Reset image size

Theoretical calculations of strong-field QED effects in experimentally relevant scenarios must deal with three-dimensional effects: the nonlinear regime a₀ ≳ 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 γ ≫ a₀ [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 a_rms and wavevector κ to be functions of spatial coordinate as well as phase. For Gaussian focusing, within the paraxial approximation, we have

$$\begin{equation}{a}_{\text{rms}}=\frac{{a}_{0}f(\psi )}{\sqrt{1+{\zeta }^{2}}}\enspace \mathrm{exp} \left(-\frac{{\rho }^{2}}{1+{\zeta }^{2}}\right),\qquad {\rho }^{2}=\frac{{x}^{2}+{y}^{2}}{{w}_{0}^{2}},\qquad \zeta =\frac{z}{{z}_{R}},\end{equation} \tag{ 19 }$$

where w₀ is the beam waist (the radius at which the intensity falls to 1/e² of its central value), ${z}_{\mathrm{R}}=\pi {w}_{0}^{2}/\lambda$ is the Rayleigh range, and the factor f(ψ) is the pulse envelope [80]. The local wavevector κ_μ = ∂_μ ψ, where ψ = ϕ − ρ² ζ/(1 + ζ²) + tan⁻¹  ζ 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 [81]: we include terms up to fourth-order in the diffraction angle ɛ = w₀/z_R in the Gaussian beam, which is then multiplied by a temporal envelope f(ϕ). Electron trajectories are determined by solution of the ponderomotive force equation, equation (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 w₀ = 2 μm, reaching a peak amplitude of a₀ = 10, with Gaussian temporal envelope of (full width at half maximum) duration 30 fs. The electrons have initial energy parameter η₀ = 0.01 (equivalent to γ₀ = 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 r₀ = w₀, such that they encounter a range of peak intensities. We have both a₀ ≫ 1 and ${a}_{0}^{2}/{\eta }_{0}\gg 1$, so the LCFA is expected to be a good approximation. The results presented in figure 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 figure 1.

Zoom In Zoom Out Reset image size

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) b from the laser axis, at ultrarelativistic velocity such that q⁰ ≃ −q³ ≫ q^⊥. We approximate ${a}_{\text{rms}}^{2}\simeq [{a}_{0}\enspace \mathrm{exp}(-{r}^{2}/{w}_{0}^{2})$ ${f(\phi )]}^{2}$ and solve the equation of motion, equation (12), perturbatively in the small parameter ≡ 1/γ₀. The first-order correction to the perpendicular momentum q^⊥ is obtained by substituting into equation (12) q⁰ = mγ₀ and r = b, i.e. assuming the electron is undeflected. The deflection angle follows as θ ≃ q^⊥/q⁰:

$$\begin{equation}{\theta }_{\mathrm{e}}\simeq \frac{{a}_{0}^{2}}{{\gamma }_{0}^{2}}\frac{b{e}^{-2{b}^{2}/{w}_{0}^{2}}}{{z}_{R}}{\int }_{-\infty }^{\infty }{f}^{2}(\phi )\enspace \mathrm{d}\phi .\end{equation} \tag{ 20 }$$

The outer ring in figures 7(a) and (b) corresponds to scattering at b = w₀/2 (shown by the black, dashed line), at which equation (20) is maximised, and the inner ring to scattering at b = w₀ (shown by the black, dotted line), which is the radius of the electron beam.

As discussed in section 3, and shown in figure 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 θ_γ ≃ a₀/γ₀ to the instantaneous quasimomentum. In the latter, it arises from the instantaneous oscillation in the electron kinetic momentum, which has characteristic angle θ_e ≃ a₀/γ₀, 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 figures 7(c) and (d), the characteristic angle is smaller than the expected value θ_γ = a₀/γ₀, 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 a₀: ${a}_{0}^{\text{eff}}(b)\simeq {a}_{0}\enspace \mathrm{exp}(-{b}^{2}/{w}_{0}^{2})$.

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 a₀ ∼ 1, η₀ = 0.1. This also provides an opportunity to crosscheck our LMA simulation results for focused lasers with theory. The latter is obtained using equation (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 ${({a}_{0}/{\gamma }_{0})}^{2}$, which is indeed very small. In this case, the laser amplitude is either a₀ = 0.5 or 2.5, its waist is w₀ = 4 μm, and its temporal envelope (electric-field) is f(ϕ) =  cos²[ϕ/(2N)] with N = 16. The electrons have energy parameter η₀ = 0.1 (equivalent to γ₀ = 1.638 × 10⁴ for a head-on collision with a laser pulse of central wavelength λ = 0.8 μm) and are distributed uniformly over a disk of radius 2w₀.

In figure 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 R_3D which approximately accounts for the effect of focal spot averaging. In the perturbative limit a_rms ≪ 1, the emission rate is proportional to ${a}_{\text{rms}}^{2}$. Thus we expect the overall number of photons, in the 3D case, to be reduced by a factor ${R}_{3\text{D}}\simeq \left(\int {a}_{\text{rms}}^{2}(b)\frac{\mathrm{d}{N}_{\mathrm{e}}}{\mathrm{d}b}\enspace \mathrm{d}b\right)/{a}_{0}^{2}$, where $\frac{\mathrm{d}{N}_{\mathrm{e}}}{\mathrm{d}b}$ is the distribution of electron impact parameters b, and we may take ${a}_{\text{rms}}(b)={a}_{0}\enspace \mathrm{exp}(-{b}^{2}/{w}_{0}^{2})$ for beam waist w₀. For a beam of electrons which are uniformly distributed over a disk of radius 2w₀, we have R_3D = (1 − e⁻⁸)/8 ≃ 0.125. The distribution of photon lightfront momentum fraction s is shown in figures 8(a) and (c) for a₀ = 0.5 and 2.5 respectively. Figures 8(b) and (d) show lineouts through the double-differential spectrum at fixed r_⊥ = a₀/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 figure 8(b) than figure 8(d), i.e. at lower a₀, 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 r_⊥, this width is ${\Delta}s/s={a}_{0}^{2}/(1+2n\eta )$. There is an additional contribution from the non-zero bandwidth of the pulse, which is given approximately by Δs/s ≃ 0.187/(ω₀ τ), where τ is the FWHM duration of the pulse intensity profile: for the cosine-squared envelope under consideration here, ω₀ τ ≃ 0.364 N and Δs/s ≃ 1/(2N). At sufficiently small a₀, it is the latter contribution, from the laser pulse bandwidth, that dominates. Note that in a focused pulse, the effective amplitude at finite impact parameter a_rms(b) < a₀ and so such effects are magnified. Integrating over a finite range of r_⊥ partially mitigates this, which is why the single-differential spectra are in much better agreement.

Zoom In Zoom Out Reset image size

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 a₀ = 2.5. This is because the position of the edge differs for electrons at different impact parameters, as increasing b means reducing the effective a₀. We have repeated the comparison between LMA-based simulations and QED for more tightly focused laser pulses, reducing the waist w₀ to 3λ and 1λ, while holding the peak a₀ and the electron-beam–laser overlap fixed. The detailed results are shown in the supplementary material (https://stacks.iop.org/NJP/23/085008/mmedia): 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 focused laser pulses.

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 modeling 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 ω₀ τ = 2πN, 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 2 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 a₀ and η, the prediction can be wrong by orders of magnitude (as demonstrated by figure 4). The requirement that ω₀ τ ≫ 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 N = 4, which is already much shorter than most typical laser pulses (where the equivalent N ≳ 10): furthermore we expect the accuracy to improve with increasing N. Similarly, the error made by the LCFA in the photon spectrum becomes arbitrarily large as s → 0, even if a₀ ≫ 1, whereas the LMA is guaranteed to obtain the correct limiting value for all a₀. The LMA result is inaccurate only in a small region around $s\simeq {s}_{0}^{{\ast}}=(\eta /N)/(1+{a}_{0}^{2}+\eta /N)$ [73], which is in any case missed by the LCFA: radiation here arises from the slow ponderomotive scattering of the electron.

Table 2. Conditions on the pulse amplitude a₀, phase duration ω₀ τ and electron energy parameter η = κ.p/m² for simulated photon yields and spectra to be accurate under the LMA and LCFA.

Predicted quantity	LCFA is accurate when	LMA is accurate when
Yield, Nγ	${a}_{0}^{2}/\eta \gg 1$, a0 ≫ 1, ω0 τ ≫ 1 a	ω0 τ ≫ 1
Spectrum, $\frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}s}$	$s\hspace{2pt}\gtrsim \hspace{2pt}\frac{2\eta }{1+{a}_{0}^{2}+2\eta }$	$s\hspace{2pt}\gtrsim \hspace{2pt}{s}_{0}^{{\ast}}$, $s\ll {s}_{0}^{{\ast}}$
		where ${s}_{0}^{{\ast}}=\frac{\eta /N}{1+{a}_{0}^{2}+\eta /N}$

^a.The last condition expresses that the photon formation length should be smaller than both the wavelength λ and the pulse length τ. It follows from the standard LCFA applicability conditions ${[({a}_{0}^{2}/\eta )s/(1-s)]}^{1/3}\gg 1$ and a₀ ≫ 1, when one chooses $s={s}_{0}^{\star }$ and $\eta \ll N{a}_{0}^{2}$, corresponding to the mid-IR peak.

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/N ≪ 1. In practice, one is often limited to such chirps by virtue of the available bandwidth [74]. 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 $\overrightarrow {k}$ obeying $\vert { \overrightarrow {k}}^{\prime }\vert /\vert \overrightarrow {k}\vert \ll 1$. This reduces to a condition on the focusing, which may be expressed through the diffraction angle ɛ = λ/(πw₀) as ɛ ≪ 1. (A similar condition applies to the validity of the ponderomotive-force approach to the particle's classical dynamics [64], which confirm for ɛ ≃ 0.12 in figure 8.) Using the high-energy approximation with the focused background in equation (19), we note z/z_R = −ϕ/(2ωz_R), and the relevant change in wave-vector on the focal axis is k'/k ∼ 1/(ωz_R) = ɛ²/2 ≪ 1. We have cross-checked 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 $\vert { \overrightarrow {k}}^{\prime }\vert /\vert \overrightarrow {k}\vert \ll 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 [43, 83, 84].

Motivated by the imminent need for precision simulations of strong-field QED processes in the transition regime a₀ ∼ 1, we have presented here a novel simulation framework which incorporates quantum effects via probability rates calculated within the 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 γ ≫ a₀.

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 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 [85–87].

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 [88–92].

We would like to thank members of the LUXE collaboration for helpful discussions during preparation of this work. We acknowledge funding from the Engineering and Physical Sciences Research Council (Grant EP/S010319/1, BK, AJM). Simulations were performed on resources provided by the Swedish National Infrastructure for Computing at the High Performance Computing Centre North.

The source code for the simulation program described is available at reference [93]. 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 reference [94].

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 a_μ (φ) = eA_μ (φ)/m is

$$\begin{equation}{a}_{\mu }(\phi )={a}_{0}f \left(\frac{\varphi }{{\Phi}}\right)\left[{\varepsilon }_{\mu }\enspace \mathrm{cos}\enspace b(\varphi )+{\beta }_{\mu }\enspace \mathrm{sin}\enspace b(\varphi )\right],\end{equation} \tag{ A1 }$$

and the phase is φ = κ ⋅ x. 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 m and four-momentum p_μ , [55],

$$\begin{equation}{{\Psi}}_{p,r}(x)=\left(1+\frac{m\overline{)\kappa }\overline{)a}(\varphi )}{2\kappa \cdot p}\right){u}_{p,r}{\mathrm{e}}^{-\mathrm{i}{S}_{p}(x)},\end{equation} \tag{ A2 }$$

where u_p are constant spinors. The Volkov phase term is given by,

$$\begin{equation}{S}_{p}(x)=p\cdot x+{\int }_{-\infty }^{\varphi } \mathrm{d}y\enspace \frac{2\enspace mp\cdot a(y)-{m}^{2}{a}^{2}(y)}{2\kappa \cdot p},\end{equation} \tag{ A3 }$$

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 $\mathsf{S}$ in the plane-wave background field only holds for three of the four directions, {−, ⊥}. As such, the scattering amplitude takes the form,

$$\begin{equation}\mathsf{S}={(2\pi )}^{3}{\delta }_{-,\perp }^{3}({p}_{\text{in}}-{p}_{\text{out}})\mathcal{M},\end{equation} \tag{ A4 }$$

where ${\delta }_{-,\perp }^{3}(p)=\delta ({p}_{-})\delta ({p}_{1})\delta ({p}_{2})$, and $\mathcal{M}$ is the invariant amplitude.

Closed form solutions to equation (A3) are not always available. A simple example is the infinite monochromatic plane wave, which is the f(φ/Φ) → 1, b(φ) → φ limit of the background field equation (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–41].

The slowly varying envelope approximation for an arbitrarily chirped plane-wave field was derived in [58], and we follow this approach here. For the circularly polarised background equation (A1), the terms which are quadratic in the field depend only on the slowly varying envelope, ${a}^{2}(\varphi )=-{a}_{0}^{2}{f}^{2}(\varphi /{\Phi})$, while the terms linear in the field contain both slow (through f) and fast (through b) timescales. This gives integrals of the form,

$$\begin{equation}\mathcal{I}={\int }_{-\infty }^{\varphi } \mathrm{d}y\enspace f\left(\frac{y}{{\Phi}}\right)\left[\mathrm{cos}\enspace b(y),\mathrm{sin}\enspace b(y)\right].\end{equation} \tag{ A5 }$$

To deal with these integrals, we first transform the trigonometric functions of b(y) to pull out a factor depending on the inverse of ω(y) = b'(y), where a prime denotes a derivative of the argument:

$$\begin{equation}\mathcal{I}={\int }_{-\infty }^{\varphi } \mathrm{d}y\enspace \frac{f(y/{\Phi})}{\omega (y)}\frac{\mathrm{d}}{\mathrm{d}y}\left[\mathrm{sin}\enspace b(y),-\mathrm{cos}\enspace b(y)\right].\end{equation} \tag{ A6 }$$

The function ω(y) 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

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}y}\frac{f(y/{\Phi})}{\omega (y)}=\frac{1}{{\Phi}}\frac{{f}^{\prime }(y/{\Phi})}{\omega (y)}-\frac{f(y/{\Phi})}{{\omega }^{2}(y)}{\omega }^{\prime }(y).\end{equation} \tag{ A7 }$$

Provided this is a small correction, which is valid for sufficiently long pulses, Φ ≫ 1 and when the derivative of the chirp function satisfies ω'(y) ≪ ω(y), we can neglect these slowly varying terms, and approximate the integrals by,

$$\begin{equation}\mathcal{I}\simeq \frac{f(\varphi /{\Phi})}{\omega (\phi )}\left[\mathrm{sin}\enspace b(\varphi ),-\mathrm{cos}\enspace b(\varphi )\right].\end{equation} \tag{ A8 }$$

Applying these approximations to the classical action S_p in equation (A3) gives,

$$\begin{equation}{S}_{p}(x)=G(\varphi )-z(\varphi )\mathrm{sin} \left[b(\varphi )-\vartheta \right].\end{equation} \tag{ A9 }$$

The function G(φ) contains only slowly varying terms, or terms linear in φ. The function z(φ) depends on the phase only through the slowly varying envelope f(φ/Φ) and local frequency ω(φ), and the angle ϑ is independent of the phase.

The exponential of the trigonometric function in equation (A9) can be expanded into an infinite sum of Bessel functions using the Jacob–Anger expansion,

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}z\enspace \mathrm{sin}(b-\vartheta )}=\sum\limits _{n=-\infty }^{\infty }{\mathrm{e}}^{-\mathrm{i}n(b-\vartheta )}{J}_{n}(z).\end{equation} \tag{ A10 }$$

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 equation (A10), the invariant amplitude, $\mathcal{M}$, in equation (A4), takes on the form,

$$\begin{equation}\mathcal{M}=\int \mathrm{d}\varphi \sum\limits _{n=-\infty }^{\infty } {\mathcal{M}}_{n}^{\prime }(\varphi ).\end{equation} \tag{ A11 }$$

The probability, $\mathsf{P}$, is then found in the usual way by squaring the scattering amplitude equation (A4) and integrating over the Lorentz invariant phase space for the particular process, dΩ_LIPS,

$$\begin{equation}P\propto \iint \mathrm{d}\varphi \enspace \mathrm{d}{\varphi }^{\prime }\sum\limits _{n,{n}^{\text{{'}}}=-\infty }^{\infty }\int \mathrm{d}{{\Omega}}_{\text{LIPS}}\enspace {\mathcal{M}}_{n}^{\prime \enspace {\dagger}}(\varphi ){\mathcal{M}}_{{n}^{\prime }}^{\prime }({\varphi }^{\prime }).\end{equation} \tag{ A12 }$$

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,

$$\begin{equation}\phi =\frac{1}{2}(\varphi +{\varphi }^{\prime }),\qquad \theta =\varphi -{\varphi }^{\prime },\end{equation} \tag{ A13 }$$

and then take the small phase difference approximation θ ≪ 1 to expand the probability in a Taylor series in θ, retaining only the leading-order, O(θ), contributions.

The θ-integral can be performed analytically, leaving the probability in the form,

$$\begin{equation}{\mathsf{P}}^{\text{LMA}}=\int \mathrm{d}\phi \enspace {\mathsf{R}}^{\text{LMA}}(\phi ).\end{equation} \tag{ A14 }$$

The function, ${\mathsf{R}}^{\text{LMA}}(\phi )$, 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 $\mathsf{R}(\phi )$ as a local rate which can be used in simulations. (In the main paper, we instead use a rate W^LMA defined as a probability per unit proper time.) To make this discussion more explicit, we consider the process of nonlinear Compton scattering.

A.1. Nonlinear Compton scattering in a chirped plane-wave pulse

Consider an electron with an initial momentum p_μ interacting with a plane-wave electromagnetic field to produce a photon of momentum k_μ and polarisation ${{\epsilon}}_{k,l}^{{\ast}}$. The scattering amplitude, in terms of the Volkov wave functions equation (A2), is given by,

$$\begin{equation}{\mathsf{S}}_{{r}^{\prime },r;l}=-ie\int {\mathrm{d}}^{4}x\enspace {\overline{{\Psi}}}_{{p}^{\prime },{r}^{\prime }}(x){\overline{){\epsilon}}}_{k,l}^{{\ast}}{\text{e}}^{\text{i}k\cdot x}{{\Psi}}_{p,r}(x).\end{equation} \tag{ A15 }$$

Here we use the Dirac slash notation, ${ \overrightarrow {r}}_{\perp }={ \overrightarrow {k}}_{\perp }/(sm)-{ \overrightarrow {p}}_{\perp }/m$, where γ^μ are the Dirac gamma matrices. The momentum p_μ ' is the momentum of the outgoing electron.

Performing all of the trivial integrations to express the scattering amplitude in the form equation (A4), the invariant amplitude is found to be,

$$\begin{equation}\mathcal{M}=-ie\int \mathrm{d}\varphi \enspace \mathcal{S}(\varphi )\mathrm{exp} \left[i{\int }_{-\infty }^{\varphi } \mathrm{d}y\enspace \frac{k\cdot \pi (y)}{\kappa \cdot (p-k)}\right],\end{equation} \tag{ A16 }$$

where the spin dependent structure is given by,

$$\begin{equation}\mathcal{S}(\varphi )={\overline{u}}_{{p}^{\prime },{r}^{\prime }}\left[1+\frac{m\overline{)a}(\varphi )\overline{)\kappa }}{2\kappa \cdot {p}^{\prime }}\right]{\overline{){\epsilon}}}_{k,l}^{{\ast}}\left[1+\frac{m\overline{)\kappa }\overline{)a}(\varphi )}{2\kappa \cdot p}\right]{u}_{p,r}.\end{equation} \tag{ A17 }$$

and the classical action in the exponent is expressed in terms of the kinetic, or local, momentum of the incoming electron,

$$\begin{equation}{\pi }^{\mu }(\varphi )={p}^{\mu }-m{a}^{\mu }(\varphi )+\frac{2\enspace mp\cdot a(\varphi )-{m}^{2}{a}^{2}(\varphi )}{2\kappa \cdot p}{\kappa }^{\mu }.\end{equation} \tag{ A18 }$$

After applying the slowly varying approximation, as detailed above, to the classical action in the exponent, the invariant amplitude equation (A16) can be expressed as

$$\begin{equation}\mathcal{M}=-ie\int \mathrm{d}\varphi \enspace \mathcal{S}(\varphi ){\mathrm{e}}^{\mathrm{i}G(\varphi )-\mathrm{i}z(\varphi )\mathrm{sin}[b(\varphi )-\vartheta ]}.\end{equation} \tag{ A19 }$$

The function G(φ) has the explicit form,

$$\begin{equation}G(\varphi )=\frac{1}{2s\kappa \cdot p(1-s)}{\int }_{-\infty }^{\varphi } \mathrm{d}y\left\{\vert { \overrightarrow {k}}_{\perp }-s{ \overrightarrow {p}}_{\perp }{\vert }^{2}+{s}^{2}{m}^{2}\left[1+{a}_{0}^{2}{f}^{2} \left(\frac{y}{{\Phi}}\right)\right]\right\},\end{equation} \tag{ A20 }$$

where we have defined the lightfront momentum fraction s = κ ⋅ k/κ ⋅ p. As stated above, this only has dependence on the phase through either linear or slowly varying terms.

The term z(φ) is

$$\begin{equation}z(\varphi )=\frac{m{a}_{0}}{\kappa \cdot p(1-s)}\frac{\vert f(\varphi /{\Phi})\vert }{\vert \omega (\varphi )\vert }\sqrt{{\left\vert { \overrightarrow {k}}_{\perp }-s{ \overrightarrow {p}}_{\perp }\right\vert }^{2}},\end{equation} \tag{ A21 }$$

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,

$$\begin{equation}\vartheta =\mathrm{arctan} \left[\frac{(k-sp)\cdot \varepsilon }{(k-sp)\cdot \beta }\right],\end{equation} \tag{ A22 }$$

and so can be interpreted as the angle between the components of the four-vector k_μ − sp_μ projected onto the directions of background field polarisation.

We skip now to the explicit form of the probability. Expanding into Bessel harmonics according to equation (A10), the probability equation (A12) becomes

$$\begin{align}{\mathsf{P}}^{\text{LMA}}=& -\frac{\alpha {m}^{2}}{4{\pi }^{2}{(\kappa \cdot p)}^{2}}\iint \mathrm{d}\varphi \enspace \mathrm{d}{\varphi }^{\prime } \sum\limits _{n,{n}^{\prime }=-\infty }^{\infty }\int \frac{\mathrm{d}s}{s(1-s)}\int {\mathrm{d}}^{2}{ \overrightarrow {k}}_{\perp }{\mathrm{e}}^{\mathrm{i}G(\varphi )-\mathrm{i}G({\varphi }^{\prime })-\mathrm{i}nb(\varphi )+\mathrm{i}{n}^{\prime }b({\varphi }^{\prime })+\mathrm{i}(n-{n}^{\prime })\vartheta }\\ & {\times}\left(\left\{1+\frac{{a}_{0}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]\left[{f}^{2}\left(\frac{\varphi }{{\Phi}}\right)+{f}^{2}\left(\frac{{\varphi }^{\prime }}{{\Phi}}\right)\right]\right\}{J}_{n}(z(\varphi )){J}_{{n}^{\prime }}(z({\varphi }^{\prime }))\right.\\ & \left.-\frac{{a}_{0}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]f\left(\frac{\varphi }{{\Phi}}\right)f\left(\frac{{\varphi }^{\prime }}{{\Phi}}\right)\left[{J}_{n+1}(z(\varphi )){J}_{{n}^{\prime }+1}(z({\varphi }^{\prime }))+{J}_{n-1}(z(\varphi )){J}_{{n}^{\prime }-1}(z({\varphi }^{\prime }))\right]\right).\end{align} \tag{ 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, b(φ), and the angle ϑ. If we consider the definitions equations (A20)–(A22), we notice that the only dependence on the transverse photon momentum is through the combination ${ \overrightarrow {r}}_{\perp }={ \overrightarrow {k}}_{\perp }/(sm)-{ \overrightarrow {p}}_{\perp }/m$. We can then shift the integration variables in equation (A23), and using equation (A22) express the integration measure in polar coordinates,

$$\begin{equation}\int {\mathrm{d}}^{2}{ \overrightarrow {k}}_{\perp }\to {s}^{2}{m}^{2}\int {\mathrm{d}}^{2}{ \overrightarrow {r}}_{\perp }=\frac{{s}^{2}{m}^{2}}{2}{\int }_{0}^{2\pi }\mathrm{d}\vartheta \int \mathrm{d}\vert { \overrightarrow {r}}_{\perp }{\vert }^{2}.\end{equation} \tag{ A24 }$$

The only dependence of the probability on the angle ϑ is then through the exponential factor exp(+i(n − n')ϑ). The integration over the angle ϑ sets n = n'. This allows the probability to be well approximated by,

$$\begin{align}\mathsf{P}\simeq & -\frac{\alpha {m}^{4}}{4\pi {(\kappa \cdot p)}^{2}}\int \int \mathrm{d}\varphi \enspace \mathrm{d}{\varphi }^{\prime } \sum\limits _{n=-\infty }^{\infty }\int \mathrm{d}s\frac{s}{(1-s)}\int \mathrm{d}\vert { \overrightarrow {r}}_{\perp }{\vert }^{2}{\mathrm{e}}^{+\mathrm{i}G(\varphi )-\mathrm{i}G({\varphi }^{\prime })}{\mathrm{e}}^{-\mathrm{i}n(b(\varphi )-b({\varphi }^{\prime }))}\\ & {\times}\left(\left\{1+\frac{{a}_{0}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]\left[{f}^{2}\left(\frac{\varphi }{{\Phi}}\right)+{f}^{2}\left(\frac{{\varphi }^{\prime }}{{\Phi}}\right)\right]\right\}{J}_{n}(z(\varphi )){J}_{n}(z({\varphi }^{\prime }))\right.\\ & \left.-\frac{{a}_{0}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]f\left(\frac{\varphi }{{\Phi}}\right)f\left(\frac{{\varphi }^{\prime }}{{\Phi}}\right)\left[{J}_{n+1}(z(\varphi )){J}_{n+1}(z({\varphi }^{\prime }))+{J}_{n-1}(z(\varphi )){J}_{n-1}(z({\varphi }^{\prime }))\right]\right).\end{align} \tag{ A25 }$$

Following through with the local expansion, using equation (A13) and θ ≪ 1, the integral over dθ can be performed, which gives a δ-function:

$$\begin{align}\mathsf{P}\simeq & -\frac{\alpha }{{\eta }_{0}}\int \mathrm{d}\phi \hspace{2pt}\sum\limits _{n=1}^{\infty }\int \mathrm{d}s\int \mathrm{d}\vert { \overrightarrow {r}}_{\perp }{\vert }^{2}\delta \left[\vert { \overrightarrow {r}}_{\perp }{\vert }^{2}+1+{a}_{0}^{2}{f}^{2}\left(\frac{\phi }{{\Phi}}\right)-\frac{2{\eta }_{0}n\omega (\phi )(1-s)}{s}\right]\\ & {\times}\left\{{J}_{n}^{2}(z(\phi ))+\frac{{a}_{0}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]{f}^{2}\left(\frac{\phi }{{\Phi}}\right)\left[2{J}_{n}^{2}(z(\phi ))-{J}_{n+1}^{2}(z(\phi ))-{J}_{n-1}^{2}(z(\phi ))\right]\right\},\end{align} \tag{ A26 }$$

where we have defined η₀ = κ ⋅ p/m². The probability only has support when the argument of the δ-function satisfies:

$$\begin{equation}\vert { \overrightarrow {r}}_{\perp }{\vert }^{2}+1+{a}_{0}^{2}{f}^{2}\left(\frac{\phi }{{\Phi}}\right)-\frac{2{\eta }_{0}n\omega (\phi )(1-s)}{s}=0,\end{equation} \tag{ A27 }$$

which (upon adapting the notation) is found to be exactly the stationary phase condition which is evaluated in [58] (see equation (25) of [58]). 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 equation (A26), which gives the relatively concise expression (suppressing explicit dependence on ϕ)

$$\begin{equation}\mathsf{P}\simeq -\frac{\alpha }{{\eta }_{0}}\int \mathrm{d}\phi \hspace{2pt}\sum\limits _{n=1}^{\infty }{\int }_{0}^{{s}_{n,{\ast}}(\phi )} \mathrm{d}s\left\{{J}_{n}^{2}({z}_{n})+\frac{{a}_{\text{rms}}^{2}}{2}\left[1+\frac{{s}^{2}}{2(1-s)}\right]\left[2{J}_{n}^{2}({z}_{n})-{J}_{n+1}^{2}({z}_{n})-{J}_{n-1}^{2}({z}_{n})\right]\right\},\end{equation} \tag{ A28 }$$

where the argument of the Bessel functions is now

$$\begin{equation}{z}_{n}(\phi )=\frac{2n{a}_{\text{rms}}}{\sqrt{1+{a}_{\text{rms}}^{2}}}\sqrt{\frac{1}{\omega (\phi )}\frac{1}{{s}_{n}(\phi )}\frac{s}{1-s}\left[1-\frac{1}{\omega (\phi )}\frac{1}{{s}_{n}(\phi )}\frac{s}{1-s}\right]},\end{equation} \tag{ A29 }$$

and we have defined the cycle-averaged potential a_rms = a₀ f(ϕ/Φ) and the upper bound on the integration over s is

$$\begin{equation}{s}_{n,{\ast}}(\phi )=\frac{{s}_{n}(\phi )\omega (\phi )}{1+{s}_{n}(\phi )\omega (\phi )},\qquad {s}_{n}(\phi )=\frac{2n{\eta }_{0}}{1+{a}_{\text{rms}}^{2}}.\end{equation} \tag{ A30 }$$

Thus, when compared with the expressions found for the LMA in a non-chirped pulse [37], the chirp function, b(ϕ), contributes an effective rescaling of the lightfront energy parameter, η₀ → η₀ ω(ϕ), inside the argument of the Bessel functions. In equation (10) we have redefined s_n and s_n,* by absorbing the local frequency, ω (where ω = ω(ϕ)), into the definition of the local energy parameter, η = η₀ ω (where η = η(ϕ)).