Higher fidelity simulations of nonlinear Breit–Wheeler pair creation in intense laser pulses

When a photon collides with a laser pulse, an electron-positron pair can be produced via the nonlinear Breit–Wheeler process. A simulation framework has been developed to calculate this process, which is based on a ponderomotive approach that includes strong-field quantum electrodynamical effects via the locally monochromatic approximation (LMA). Here we compare simulation predictions for a variety of observables, in different physical regimes, with numerical evaluation of exact analytical results from theory. For the case of a focussed laser background, we also compare simulation with a high-energy theory approximation. These comparisons are used to quantify the accuracy of the simulation approach in calculating harmonic structure, which appears in the lightfront momentum and angular spectra of outgoing particles, and the transition from multi-photon to all-order pair creation. Calculation of the total yield of pairs over a range of intensity parameters is also used to assess the accuracy of the locally constant field approximation (LCFA).


Introduction
A photon propagating through a laser pulse can decay to an electron-positron pair. The phenomenology of the interaction between the photon and the laser pulse is dependent on the centre of mass energy of the collision, the strength of the laser field, and the bandwidth of the pulse. At weak field strengths, the leading-order contribution is determined by the number of photons that must be taken from the laser pulse to overcome the mass threshold. If only one is necessary, the process is referred to as the linear Breit-Wheeler process [1]. At lower centre of mass energies, the leading-order kinematicallyallowed process may involve several laser photons, which is the 'multi-photon' nonlinear Breit-Wheeler process. For a e-mail: tom.blackburn@physics.gu.se b e-mail: b.king@plymouth.ac.uk (corresponding author) higher field strengths, the 'threshold' number of laser photons required for the photon to decay may not be the most probable channel, and all orders of interaction between the laser pulse and the photon must be taken into account.
To measure nonlinear Breit-Wheeler pair creation  in the 'all-order' regime, a sufficiently powerful laser is required: peak intensities > 10 22 Wcm −2 are now accessible with current laser technology [27][28][29][30]. Typical proposals involve providing a particle beam via laser-wakefield acceleration of electrons, which can produce energies of the order of several GeV [31][32][33][34]. This means that, in order for pair creation to be kinematically allowed, many photons are required. If the laser pulse is sufficiently intense, nonlinear Breit-Wheeler becomes sufficiently probable as to be measurable. In this regime of low energy and high field strength, the 'locally constant field approximation' (LCFA) [21,[35][36][37][38][39][40] of approximating the probability for the process as a sum of interactions with 'instantaneously constant' phase slices of the laser pulse, is expected to be sufficiently accurate. The LCFA is the standard method for including strong-field QED effects in numerical simulation of laser-matter interactions [41][42][43][44][45][46][47][48].
Recently, complementary 'high-energy' experiments such as E320 at SLAC and LUXE [49] at DESY have been suggested to measure all-order strong-field QED effects such as nonlinear Breit-Wheeler, Compton and the nonlinear trident process [50][51][52][53][54][55][56][57]. These experiments will collide particle beams accelerated using conventional radiofrequency cavities, with strong laser pulses. The particle beams can reach higher energies (13 GeV for E320 and 11.5-16.5 GeV for LUXE), lower emittances, higher repetition rate and energy stability than current laser-wakefield accelerated beams have been achieved. (LUXE will also collide the driving electron beam with a solid target to generate a source of bremsstrahlung radiation that will collide with the intense optical laser and form a photon-photon experiment [20,49,[58][59][60].) These types of experiments therefore allow for: (i) higher precision measurements of total and differential yields; (ii) access to strong-field QED processes, such as nonlinear Breit-Wheeler, at lower intensities where the transition from perturbative to non-perturbative physics occurs. The field strengths where this transition takes place, are outside the region of applicability of the LCFA. Therefore, a new simulation framework must be developed to correctly model the physics at these experiments.
So far, the only measurement of the nonlinear Breit-Wheeler process was in the landmark E144 experiment, where, in the multiphoton regime [50], both nonlinear Compton scattering [61] and Breit-Wheeler [62] were measured. The experiment was modelled using an 'instantaneously monochromatic' approximation [63], which has since been used in the simulation codes CAIN [64] and IPStrong [65]. This approximation has recently been formalised, by being derived directly from strong-field QED in a plane-wave background [66]. The resulting form, the 'locally monochromatic approximation' (LMA), assumes a pulse envelope that varies much slower than the carrier wavelength, and employs a local phase expansion, which includes interference effects between processes taking place within the same laser wavelength. (The slowly varying envelope approximation in the LMA has been used by other authors in first [67][68][69][70], and second-order processes [71].) Unlike the LCFA, the LMA is not restricted by its use in a particular intensity or energy regime, but it does assume that the interacting electromagnetic field is well-approximated by a plane electromagnetic wave. Using a ponderomotive scattering approach, the LMA has been realised in the simulation code, Ptarmigan [72], which is being used to model the interaction point physics of the LUXE experiment [49]. A comprehensive benchmarking of the LMA with exact expressions from QED, was recently performed for the process of nonlinear Compton scattering [73], and an analysis of pulse shape effects beyond the LMA performed in [74]. The LMA has also recently been used to approximate the two-vertex process of trident electronpositron pair creation [49,75].
Further motivation for improving the accuracy of modelling strong-field QED processes at the LUXE experiment, is provided by the search for new physics. The 'LUXE New Physics search with Optical Dump' LUXE-NPOD [76] is based on a secondary production mechanism that utilises the high-energy inverse Compton-scattered photons [77] produced in the electron beam-laser collision, to produce axionlike particles (ALPs) in a beam dump further downstream of the experiment. Such a setup will be sensitive to ALPs with masses O(MeV)-O(GeV), which is a range that has attracted much attention in recent years [78][79][80][81][82][83][84][85][86]. It is important, therefore, that the modelling of the photon source is as accurate as possible. The strong-field QED source of photons in LUXE-NPOD is modelled using Ptarmigan, which we further benchmark in the current paper. The approach we use to model strong-field QED in particle-laser collisions can be adapted to model the generation of ALPs at the interaction point itself [87][88][89][90][91][92][93].
In the current paper, we calculate totally inclusive and differential probabilities for nonlinear Breit-Wheeler pair creation, using the LMA in a numerical simulation framework, and compare them to the prediction from strong-field QED employing plane-wave Volkov states [94]. Through a series of benchmarks for a typical pulse shape, for a range of parameters and observables, we acquire a measure of the accuracy of the LMA. Our comparison covers a regime soon to be explored in upcoming experiments, but also extends to much higher energies, where harmonic structure appears in particle spectra. We highlight pulse envelope effects that are beyond the LMA, and establish parameters where the LCFA can also be employed.
The paper is organised as follows. In Sect. 2, the phenomenology of nonlinear Breit-Wheeler is recapped, and the expressions used in the calculation directly from theory are stated; in Sect. 3 the implementation of the LMA in the simulation framework is explained; in Sect. 4 the results of benchmarking between the direct theory calculation and numerical simulation are presented. In Sect. 5 the results are discussed and the paper is concluded. Unless otherwise stated, c and have been set to unity.

Theory background
To aid the understanding of the rest of the paper, some of the main phenomenology associated with pair creation in a quasimonochromatic field will be recapped. (Reviews of strongfield QED in a laser background can be found in [35,95].) We define the rescaled vector potential a = e A, as: and a = 0 otherwise, where N is the number of cycles, ϕ = κ · x is the phase, κ is the laser wavevector (satisfying κ · a = 0), ξ is the classical nonlinearity parameter (also referred to as the 'intensity parameter') and m and e > 0 are the positron mass and charge, respectively. We note that, using a vector potential that is non-zero only on a finite lightfront interval will mean having a background with a wide bandwidth. This will lead to an enhancement of pulse-envelope interference effects, which we will comment on when they arise. The momentum contributed by the background field can be written in terms of harmonics n, of the central laser wavevector. In the LMA, conservation of momentum in non-linear Breit-Wheeler can be written: where q and q are the electron and positron quasimomenta, q = p − (a 2 /2κ · p)κ and p is the free electron momentum (and analogously p for the positron), k is the momentum of the photon and we allown to be a real number. Solving Eq.
(2) forn, we find: where s = κ · p /κ · k is the lightfront momentum fraction, η = κ · k/m 2 is the (photon) energy parameter and p ⊥ is the electron momentum transverse to the laser propagation direction (and analogously for the positron, p ⊥ ). We have also defined ξ 2 (ϕ) = −a 2 , which, for the here-considered vector potential, Eq. (1), becomes and ξ(ϕ) = 0 otherwise. From Eq. (2), it can be seen s ∈ [0, 1] and since p ⊥ , and p ⊥ are integrated over, there is no upper bound ton. But there is a lower bound, at p ⊥ = p ⊥ = 0 and s = 1/2, giving: The 'threshold' harmonic for pair creation to occur is then n = n . When ξ 1, the probability of this harmonic contributing to pair creation scales as ∼ ξ 2n , which is perturbative, because the next highest harmonic is a factor ξ 2 1 compared to the leading harmonic, but n can be large and so the leading order generally depends nonlinearly on the background intensity, ξ 2 . This is the 'multi-photon regime'. As ξ increases to a value where it no longer satisfies ξ 1, although one can still define a threshold harmonic, it is no longer the case that the standard harmonic hierarchy holds and that each successive higher harmonic is suppressed. Furthermore, the quasimomentum of the produced pair increases, and so it becomes more difficult to create a pair, and the curve of the probability arches down and away from the multiphoton scaling. If ξ is increased sufficiently, the LMA tends to the LCFA (locally constant field approximation), which depends only on the combination of parameters, χ = ξη. At small χ (but large ξ ) this allows for a tunneling dependency on χ , whereas for large χ , the probability for the tree-level process of nonlinear Breit-Wheeler scales as P ∼ χ 2/3 . (This non-perturbative dependency on the field strength has recently been discussed in the context of the Ritus-Narozhny conjecture [96][97][98][99][100][101][102][103].) In a plane wave background, the total probability for pair creation, P, (or equivalently: yield of pairs per photon), will depend on the intensity parameter ξ , the probe particle energy , where at fixed η, P ∝ ξ 2n , where n is the threshold harmonic (for η = 0.2, n = 11). Inset right: the tunneling regime (ξ 1, χ 1), in which the probability is exponentially suppressed as P ∼ χ exp(−8/3χ) parameter η = κ · k/m 2 and the number of laser cycles N . We illustrate the various phenomenological regimes of P(ξ, η, N ) in Fig. 1.
We calculate the probability for pair creation directly from QED in a plane-wave background. This involves using exact solutions to the Dirac equation in a plane wave potential (Volkov states [94]) to include the charge-laser coupling to all orders, and then perturbatively expanding in vertices of radiated particles. We evaluate the amplitude for each polarisation channel, mod-square and then integrate over outgoing momenta. We outline the formulas used here for completeness and to introduce notation, but an analogous approach can be found in the literature.
The total unpolarised probability P can be written as: where r ⊥ = p ⊥ − sk ⊥ and σ k, p, p are the polarisation-state indices of the photon, electron and positron respectively. The (scaled) amplitude is given by: whereū σ p (v σ p ) are the outgoing electron (positron) bispinors, = 1 − k.π p /k. p is the regularisation factor, π p is the classical positron kinetic momentum: We will refer to the numerical evaluation of the above as the 'QED' result. Although the formula is exact for a plane wave, the numerical evaluation will have a finite error. It will be useful at times to compare our results with the probability of the linear Breit Wheeler process, P BW , which is equal to the leading-order expansion in ξ 2 of Eq. (6): where: is the Fourier transform of the vector potential, a(ψ) = dϕ a(ϕ) exp(iϕψ).

Particle dynamics
Numerical simulations of strong-field QED interactions are fundamentally semiclassical, in that probability rates for the processes under consideration are evaluated along the particles' classical trajectories [104,105]. The derivation of a 'rate' from QED generally involves some level of approximation, as it must be locally well-defined. Two such frameworks are available in our Monte Carlo simulation code, ptarmigan [72]: the locally monochromatic [66] and locally constant field approximations [35]. In this work we focus on the LMA-based approach: for reviews of LCFA-based simulations, see [45,106]. We have already described the simulation concept in [73], so we present only a summary here.
In the LMA approach, the classical trajectory is defined by the particle quasimomentum, which is the cycle-averaged value of the kinetic momentum. For a photon, which propagates ballistically through the strong-field region, the quasiand kinetic momenta coincide. For an electron (or positron), the quasimomentum q (q ), is the cycle average of the kinetic momentum (for the positron given in Eq. (11)), and evolves according to the following equation of motion [107]: where τ is the proper time and a rms is the component of the potential that varies slowly with respect to the wavelength. For a pulsed wave, we have a rms ( i.e. the square root of Eq. (4). The quasimomentum predicted by Eq. (13) for an electron in a pulsed plane wave, where p is the asymptotic momentum, is a slowly varying function of phase. The same result can be obtained by calculating a phase-window average of the kinetic momentum π , which is obtained from Eq. (14). The cycle-averaged position X μ , i.e. the component of the worldline that is slowly varying with respect to the laser wavelength, follows from dX μ /dτ = q μ /(m 1 + a 2 rms ). The local value of q controls the probability rates of photon emission and electron-positron pair creation through the parameters a rms and η: for an electron or positron, these are a rms = q 2 /m 2 − 1 and η = κ.q/m 2 .
Using the LMA means that the fast oscillating component of the trajectory is included at the level of the probability rates, which incorporate the conservation of quasimomentum [73]. In the LCFA framework, it is instead encoded in the trajectory itself, via the kinetic momentum, which evolves according to the Lorentz force equation: The rates of QED processes are then controlled by a single quantum parameter χ = e −(F μν π ν ) 2 /m 3 , which is defined instantaneously along the trajectory x μ (τ ), and dx μ /dτ = π μ /m. As is done in the particle-in-cell approach to kinetic plasma simulations, the code uses ensembles of 'macroparticles' to model real particle beams. The number of real particles each macroparticle represents is called the 'weight'.

Event generation
Under the LMA. The pair creation rate per unit lab time, W , of a photon with momentum k, embedded in a circularly polarised (CP) plane wave of normalized, root-mean-square amplitude a rms and wavevector κ, is given by [66] dW n ds = αm 2 where the argument of the Bessel functions, z > 0 fulfils: and the auxiliary variables are The lightfront momentum fraction s = κ.q/κ.k is restricted to and therefore the rate is non-zero only for s n > 4, or equivalently, harmonic orders n > n = 2(1 + a 2 rms )/η . A pair creation event occurs if, in a small interval of time t along the trajectory (where t is the lab time), the probability W t satisfies U < W t, where U is a pseudorandom number uniformly distributed in (0, 1). In principle, the total rate W can be obtained directly by integrating Eq. (15) over all s and then summing all relevant n: however, it is much faster to precalculate W as a function of a rms and η and implement the evaluation as a table lookup. We restrict the sum to n ≤ n ≤ n max , where n is defined by Eq. (5). The cutoff n max is defined to be the lowest harmonic order that satisfies W n < 10 −4 n i=n W i . In order to obtain a simple, analytical approximation for n max as a function of a rms and η, we obtain it numerically for a range of (a rms , η) points in the region a rms < 10 and η < 1 and then fit a trial function to the data obtained. Our result is We use the following procedure to obtain the quasimomentum of the created electron (equivalently, positron). When pair creation occurs, the harmonic index n and lightfront momentum transfer s are pseudorandomly sampled from the emission rates, Eq. (15). The former is obtained by solving U = n i=n W i / n max i=n W i , where U is a pseudorandom number drawn on the unit interval and W n is the nth partial rate, i.e. Eq. (15) integrated over all s. The cutoff, n max , is the same that used when precalculating the total rate. The lightfront momentum transfer s is obtained by rejection sampling of Eq. (15).
The momentum q is fixed by k, n and s. In the zero momentum frame (ZMF), the electron (positron) is created with energy and momentum (22) and scattering angle The azimuthal angle ϕ zmf is pseudorandomly determined in 0 ≤ ϕ zmf < 2π . The momentum so obtained is then transformed back to the lab frame, using the fact that the fourvelocity of the ZMF is u zmf = (k + nκ)/ (k + nκ) 2 . The quasimomentum of the other particle, q , follows from quasimomentum conservation, q = k + nκ − q. Under the LCFA. We also compare the theoretical results to LCFA-based simulations. The relevant rate, resolved in both positron energy ε = f k 0 and scattering angle ϑ, is [108] (see also [39]) where the photon quantum parameter If an electron-positron pair is created, the (kinetic) momenta are determined by conserving threemomentum k = π + π [109], where the components of π follow from the energy and angle pseudorandomly sampled from Eq. (24). Previous implementations have generally used the angularly integrated form of Eq. (24) and assumed that the electron, positron and photon momenta are all collinear, as ϑ is of order m/k 0 .

Biasing
The extreme rarity of pair creation events in certain regions of parameter space, which may be seen in Fig. 1, is a challenge for Monte Carlo simulations. As events are generated pseudorandomly, the number of macrophotons must be significantly larger than 1/P 1 to resolve the positron yield (and larger still, to resolve differential quantities such as the spectrum). In order to overcome this, our simulations implement a simple form of "event biasing", wherein the macrophoton partially decays. The pair-creation rate, W , is artificially increased by a large factor R ↑ 1, while the weights of any macroelectrons and macropositrons that are created are reduced by the same factor. Thus, if the weight of the photon before pair creation is w γ , the electron and positron are created with identical weights w ± = w γ /R ↑ and the photon weight is changed to w γ = w γ (1 − 1/R ↑ ). As the code checks for pair creation in a single timestep t by comparing In the absence of event biasing, we would need > 10 8 macrophotons in order to resolve the yield with satisfactory accuracy. Instead we use 10 4 macrophotons, each with weight w γ = 10 −4 , and increase R ↑ from 10 to 10 10 . The probability, which is equivalent to the sum of the macropositron weights, is given in Fig. 2a and the raw number of macropositrons in Fig. 2b. As may be expected, if R ↑ 10 4 P 10 3 , no macropositrons are generated. When R ↑ increases above this value, the probability abruptly increases, fluctuates, and then converges to the expected value. Convergence is driven by the increasing number of macropositrons, which grows linearly with R ↑ : the size of the fluctuations therefore scales as 1/ √ R ↑ . For the largest values of the rate increase, the number of macropositrons deviates from a linear scaling because of the upper bound that is automatically placed.
The simulation results we present in this manuscript are all generated using R ↑ > 1. Choosing the value of the rate increase is a matter of balancing two competing concerns: on the one hand, if it is too small, the number of macropositrons is very small and the probability is poorly estimated; on the other hand, if it is too large, many low-weight particles must be advanced through the simulation domain and checked for, e.g., secondary photon emission. Thus the code can spend an inordinate amount of time dealing with particles that have little overall importance to downstream uses such as detector simulation. We have found that a reasonable balance is achieved by setting R ↑ to the ratio between the photon-emission and pair-creation rates, evaluated at the same values of ξ and η.
We conclude by noting that, in situations where the incident photon beam has a broad energy spectrum, increasing R ↑ cannot compensate for poor statistics in the high-energy tail, which generally provides the dominant contribution to the positron yield. In this case, the number of macrophotons must be increased as well. In the present work, we deal exclusively with monoenergetic photon beams, so this is not a concern.

Benchmarking
In this section, we compare the predictions of numerical simulation, with direct evaluation of the QED expressions. Whilst the numerical simulations employ the LMA, the QED expressions do not. The aim is to understand how accurate a simulation employing the LMA can be, and in what parameter region it is accurate when calculating total yields and differential spectra. We expect agreement between the direct QED and LMA results everywhere that the LMA holds, with some disagreement for very short pulses (i.e. N 1). In order to achieve this aim, we will first benchmark in a plane wave background, where we know the direct QED expression to be accurate.
To benchmark in a plane wave background, we will compare the total yield, in each of the three variables: intensity, ξ ; lightfront momentum η and number of cycles, N . We will then benchmark some of the differential spectra: in positron lightfront momentum fraction, s and in positron transverse momentum, r ⊥ (recalling that r ⊥ = p ⊥ −sk ⊥ ) and k ⊥ is set to zero. As an example demonstrating a difference between QED and the LMA, we will consider a chirped laser pulse with more extreme parameters. Then to end this section, we will consider a focused background, which involves comparing two approximations with each other.
To quantify the accuracy of the simulation, we will define a measure of the error E, given by: where P lma is the probability as calculated from simulation using the LMA. (This form of the error was chosen so that it becomes large in the linear Breit-Wheeler region, where P lma underestimates P.) We will also calculate E lcfa (ξ, η, N ), which is the equivalent accuracy measure, but where the locally constant field approximation is used in numerical simulation instead.
To make the comparison as faithful as possible, the simulation results have been generated with the recoil of electron/positron momenta due to photon emission, disabled. Thus electrons and positrons that are created do emit secondary radiation, but their momenta are unchanged when they do so. While recoil would not affect the yield, it would shift the momentum spectrum to smaller values of s. Including these higher order effects, i.e. radiation reaction [110], theoretically would involve a calculation of at least the 'phototrident' process [111,112], which is beyond the scope of the present work.
The intensity parameter range has been chosen to include where the LMA predicts something different to the standard approximation, the LCFA. Because it has been shown that in the region ξ 1, n 1, the LMA tends to the LCFA in a plane wave background [66], and the LCFA has already been benchmarked in several works in the literature [36][37][38][39][40]113], we will concentrate on intensities that span the perturbative ξ 1 to the 'intermediate' intensity regime ξ ∼ O(1). It is known from these literature works that the LCFA should be accurate when the condition ξ 1 and ξ 2 (ϕ)/η 1 are fulfilled. It is therefore of particular interest, to ascertain what parameter values obey these conditions.
The energy parameter range has been chosen to span the range from those that are just beyond experimentally accessible with state-of-the-art strong-field QED laser wakefield acceleration experiments [114,115], η = 0.05 (approximately 4 GeV for a head-on collision with optical laser photons), up to values where the harmonic structure becomes evident in spectra, η ∼ O(1).
The pulse duration range has been chosen to reflect typically available laser pulse durations, with N = 16 cycles (corresponding to a full-width-at-half-maximum duration of 20 fs for a wavelength of 1 µm [fs]) being a standard choice.
The pulse envelope shape is the one given in Eq. (1): half a period of a squared cosine. Being a pulse with finite support has the advantage that it is more straightforward to benchmark with QED. However, such a pulse has a much wider Fourier spectrum than a pulse used in experiment (high and low energy tails of the spectrum would not be transmitted through all the optical elements). Therefore, we will see in what follows, that pulse envelope effects are significantly enhanced, but this allows us to better understand where the LMA may potentially fail.

Yield
We refer to the totally inclusive probability as the yield. First, we begin with benchmarking the yield as a function of the intensity parameter, P = P(ξ ; η, N ) where η and N are constant values. In Fig. 3 we present results for N = 16 which for a 1 µm wavelength (photon energy 1.24 eV) laser background corresponds to a full-width-athalf-maximum duration of 20 fs, and choose energy parameters η ∈ {0.05, 0.1, 0.2, 1}, which correspond, for a head-on collision of a high-energy probe photon with the propagat- ing wave, to photon energies of the order of 4, 8, 16, 80 GeV respectively. We make the following observations.
For most of the intensity range, the LMA agrees very well with the direct QED result, with an approximate error of the order of E ≈ 2-4%. We contrast this with the expansion parameter in the slowly-varying-envelopeapproximation used in the LMA, which is of the order of 1/ = 3.125% for pulse duration = 2N .
The LCFA is inaccurate for ξ < 1, which is expected, but we point out that due to the strong scaling with intensity in the multi-photon regime, the LCFA is incorrect by orders of magnitude. As is already known, the LCFA is more accurate for pair creation at higher photon energies, and this is reflected in the benchmarking. We also note the shallower gradient in how the accuracy of the LCFA increases as ξ is increased: for example, an accuracy of 50% is achieved at η = 1 already when ξ ≈ 0.8, but to achieve an accuracy of 5% requires ξ ≈ 2.5.
Below some value of the intensity parameter, ξ , the QED signal is dominated by the linear Breit-Wheeler process. A comparison is made in Fig. 3 between the QED result, and the linear Breit-Wheeler process from Eq. (12). The enhancement of linear Breit-Wheeler is emphasised in the pulse form we have chosen. This is because the pulse is only non-zero for |ϕ| < N π , and hence is effectively multiplied by a flattop envelope of width 2N π . Since the linear Breit-Wheeler signal is proportional to the square of the Fourier transform of the potential, and since the flat-top envelope has a large bandwidth, this allows the linear Breit-Wheeler process to occur (a) (b) Fig. 4 a Pair creation probability P for a photon with energy parameter η = 0.2 colliding with a plane-wave, CP laser pulse of peak amplitude ξ and duration equivalent to a number of wavelengths N : results from QED for ξ = 0.35 (blue), ξ = 0.5 (green) and ξ = 1.0 (orange) and LMA-based simulations (black, dashed). b Percentage error made by the simulations. Vertical, dashed lines give the N at which the linear Breit-Wheeler contribution from the pulse envelope surpasses the multi-photon scaling from the LMA with a higher probability than e.g. in a Gaussian envelope of a comparable width. proceed, compared to, e.g. a Gaussian envelope. Since this bandwidth effect, which is central to the linear Breit-Wheeler process, is due to pulse envelope interference effects, and since the LMA only includes carrier frequency interference, the QED results deviate from the LMA at the linear Breit-Wheeler ξ 2 'floor' shown in Fig. 3.
To illustrate the dependency of the yield on the number of laser cycles, N , P = P(N ; ξ, η), we fix η = 0.2 and choose ξ from {0.35, 0.5, 1.0}. As we saw in Fig. 3, for N = 16, these ξ values span a range including the linear Breit-Wheeler, multiphoton nonlinear Breit-Wheeler and nonperturbative nonlinear Breit-Wheeler physics. Since the bandwidth is proportional to 1/N , we would expect that, as N is reduced, the contribution from linear Breit Wheeler should increase. In Fig. 4, this is indeed what we find. This is contrasted with the prediction from the LMA. The LMA probability depends on the pulse envelope, which is a function of ϕ/2N and so a simple change of integration variables shows that the LMA simply scales linearly with N , which is also demonstrated in Fig. 4.
One can estimate the pulse parameters where the linear contribution dominates by equating the perturbative result Eq. (12), with the prediction from the LMA. The predicted value of N where linear Breit Wheeler is of equal probability to the LMA prediction, is indicated by the vertical, dashed lines in Fig. 4, which shows good agreement with the QED results.
To complete this yield section, we analyse the dependency on lightfront momentum, η. In Fig. 5, we plot the function P = P(η; ξ, N ), fixing the field parameter at ξ 0 = 0.2, and comparing three different fixed number of laser cycles, N 0 ∈ {4, 16, 32}. We recall that the 'threshold harmonic', n = n (withn given by Eq. (5)), is a function of phase, and therefore in a given laser pulse, different parts of the pulse can access different threshold harmonics. However, we also recall that the higher the field parameter, ξ , the more probable the Breit-Wheeler process is. When considering the yield of pairs the relevant threshold harmonic, is the one given by the maximum value of ξ(ϕ) in the pulse, namely ξ . In Fig. 5, we note that the η values of the first four harmonics, correspond to 'steps' in the 'staircase'-like dependency of the yield on η. For the short pulse (N = 4), the full QED result gives a smooth increase at each harmonic, and as the pulse increases in duration, the QED results tends towards the LMA prediction from simulation. Since the LMA involves application of the 'slowly-varying-envelope-approximation' and neglects terms in the Kibble mass of order ∼ 1/N , we expect better agreement with QED as N is increased, which is indeed what we find.
In Fig. 5b, we plot the error function, E(η; ξ, N ). Whilst the general trend is that a longer pulse (larger N ), leads to a smaller value of E, if the parameters are such that P is close to a channel opening, the error can increase. This is particularly the case for shorter pulses, where the slowlyvarying-envelope-approximation is already predicted to lead to a lower accuracy. However, we see that if the pulse is sufficiently long, E becomes relatively insensitive to the channelopening effect and the error saturates.
Finally, we note how Fig. 5 shows that if η is reduced enough, the probability due to the carrier frequency drops sufficiently, that the finite bandwidth linear Breit-Wheeler probability, contained in the QED but not the LMA, dominates again.

Lightfront momentum spectra
Here we plot the single differential probability, dP(s; ξ, η, N )/ds, which, when integrated over s ∈ [0, 1], gives the total yield, P, as in the previous section. We pick two cases to investigate: (i) ξ = 0.2, N = 16, and various energy parameters η ∈ [1.9, 2.4], which correspond to the opening of the n = 1 harmonic channel; and (ii) η = 0.2, N = 16, and various intensity parameters around the intermediate range, ξ ∈ [0.5, 2.5]. The parameters in (i) demonstrate how the harmonic channel-opening phenomenon is approximated by the LMA; those in (ii) correspond to those that will be used in upcoming experiments [49], but for which there is less structure in the spectra as the threshold harmonic is much larger than one. Due to the symmetry of the spectrum around s = 1/2, we will only plot the 1/2 ≤ s < 1 part of the spectra.
The plot of the yield in Fig. 5 revealed a staircase-like structure due to the opening of harmonic channels as η is varied. This phenomenon can also be revealed in the lightfront momentum spectrum. In Fig. 6, we demonstrate this, by picking parameters that best display the effect, but are beyond what has been considered for near-future experiments. The parameters ξ = 0.2 and N = 16 are set as constant, and various constant values of η are plotted. For these parameters, the first harmonic channel 'opens' when η = 2.08. Therefore in Fig. 6, we vary η parameters across this value, to show how channel-opening manifests in the lightfront momentum spectrum. We find that, as η is increased, the centre of the distribution rises as a peak, demonstrating channel opening first, and as η is further increased, the width of the peak broadens. Recalling the introductory discussion about momentum conservation and the results Eq. (3) and Eq. (5), we note that the 'threshold' harmonic corresponds to a condition fulfilled by the most probable lightfront momentum, which is at s = 1/2. The 'threshold' then increases as s is increased/decreased from the central value. The result in Fig.  6 reflects this behaviour.
Comparing the QED with the LMA result, we see that whereas QED describes a smooth transition as the harmonic range is crossed over, the LMA result displays a much clearer 'jump'. For η < 2.08, there is no discernible change in the LMA spectrum, which reflects the fact that the LMA is only including interference effects over a wavelength of the pulse, and therefore the transition is not 'softened' by the finite spectrum of the pulse envelope. This 'jump' of the LMA, is demonstrated more clearly in Fig. 6, where the simulation and QED results are compared with η as a parameter as well as s. The sudden jump of the simulation spectrum can clearly be seen, although at s = 1/2, some smoothness can be discerned. The error in the LMA is largest at η = 2 because the full QED result can already access the first harmonic (due to the finite bandwidth of the envelope), whereas for η < 2.08, the LMA is essentially 'blind' to the first harmonic. (See also [116], which analyses pulse envelope effects on harmonic transitions in more detail.) Also visible in Fig. 6 is the second harmonic, which is sharper in the simulation results than in QED. Also evident, is the Monte-Carlo noise on the simulation results.
We illustrate the spectra for the soon experimentallyaccessible case, in Fig. 7. To aid comparison, the spectra have been normalised by their maximum values. We note  Fig. 7 a Lightfront momentum spectra dP/ds for a photon with energy parameter η = 0.2 colliding with a plane-wave, CP laser pulse of peak amplitude ξ and duration equivalent to a number of wavelengths N = 16: results from QED for ξ = 0.5 (blue), ξ = 1.0 (orange), ξ = 1.5 (green), ξ = 2.0 (red) and ξ = 2.5 (grey); and LMA-based simulations (black, dashed). b Percentage error made by the simulations the behaviour, that as the intensity increases, the lightfrontmomentum spectrum becomes wider. Calculation of the first differential of the error function, dE/ds reveals an accuracy of around the 2 − 5% level for all the intensity values, with the exception of ξ = 0.5. We recall from the study of the yield in Fig. 3, that (ξ, η, N ) ≈ (0.5, 0.2, 16) is the point at which the LMA started to diverge from QED due to the linear Breit-Wheeler effect for the parameters we have chosen. Since linear Breit-Wheeler has a wider momentum spectrum [116], and since at ξ = 0.5, the spectrum is heavily suppressed at values of momentum in the 'tails', if s is far enough from the centre at s = 1/2, the linear Breit-Wheeler contribution will dominate. This explains the large increase in the error in Fig. 6 for s > 0.7 and ξ = 0.5.

Angular spectra
In this section, the double differential spectrum in the positron's perpendicular momentum r (recalling that r = |r ⊥ |/m) and the the lightfront momentum, s, is presented. In calculating d 2 P(r, s; ξ, η, N )/drds, η = 0.2 and N = 16 are chosen constant, and three cases of constant ξ are plotted. ξ = 0.2 has been chosen as an example of the linear Breit-Wheeler signal, ξ = 0.5 is in the multiphoton regime, and ξ = 2.5 is in the all-order regime.
In the ξ = 0.2 plot in Fig. 8 (top row), the double-differential spectrum has a fine, highly-oscillating structure. The period of oscillation becomes shorter as r is increased from r = 0. In the figure, above r ≈ 0.3, the oscillation is no longer resolved by the density of data points that was calculated, and any structure along the r -axis that can be discerned above this value, is just an aliasing effect. A lineout at s = 0.5 is presented for comparison. When the corresponding spectrum is calculated using the linear Breit-Wheeler formula Eq. (12), we see a good agreement when r < 0.3, but where the full QED spectrum is slightly shifted with respect to the perturbative result. Although r 0.3 is difficult to resolve, the overall shape and magnitude of the rest of the spectra show excellent agreement. For comparison in the middle plot, the prediction of the LMA is presented, but as this only includes carrier-frequency and not pulse-envelope interference, misses the linear Breit-Wheeler contribution entirely, as expected.
The ξ = 0.5 plots in Fig. 8 (middle row), illustrate the double-differential spectrum in the multi-photon case. The threshold harmonic is n = 13, but already at ξ = 0.5, the threshold harmonic does not give the leading contribution, rather this occurs at some orders above the threshold. In the figure it is the harmonics from n = 15 to n = 23 harmonics that are visible in the spectrum. In the lineout at s = 0.5 (right-hand plot), the main contribution from the QED result agrees very well with the LMA from numerical simulation. Sub-harmonics between the main peaks, which are due to interference on longer length scales than a wavelength (and hence are beyond the LMA), can also be observed in the QED result.
The ξ = 2.5 plots in Fig. 8 (final row) demonstrate the angular spectra in the 'all-order' regime, where many harmonic orders contribute to pair creation. Although the threshold harmonic here is n = 73, by the number of peaks in the plots, it can be seen that the harmonic order that contributes most, is much larger than n = 73. In the parameter region n 1, ξ 1, it is known that predictions using the LMA tend to those using the LCFA [66]. Confirmation of this is seen in the right-hand panel lineout at s = 0.5, there the LCFA shows good agreement with both the QED and the LMA curves. (This agrees with what we expect from the yield plot in Fig. 3.)

Chirped pulses
The LMA can be used with a spacetime-dependent wavevector. An example of this is a chirped plane wave. In this section, we consider pair creation in a symmetrically chirped pulse with scaled vector potential defined over the domain |ϕ| < N π as: with the chirp function: and a = 0 otherwise. This would give a 'local' energy parameter of the formη(ϕ) = η(1 − bϕ/N ), where η is the 'unchirped' value. At the centre of the pulse,η = η and the chirp is effectively zero, but towards the leading edge of the pulseη takes its maximum value. In Fig. 9, a comparison is made between QED and simulation, for the calculation of the transverse momentum spectrum, dP(r ; ξ, η, N , b)/dr, with ξ = 0.2, η = 1.5, N = 16 and the chirp parameter b varied up to its highest, physical value b = 1/π (at which point, at the trailing of the pulse,η(ϕ) = 0). This set of parameters is a rather extreme case: in the middle of the pulse,η(0) = 1.5, which is almost midway between the energy values corresponding to the n = 2 and n = 1 harmonics (as can be seen from Fig. 5). The leading edge of the pulse is at a frequency which is below the n = 1 harmonic threshold Fig. 9 for b = 1/4π and b = 1/3π (first row), and above the harmonic threshold for b = 1/2π and b = 1/π (second row). Therefore the parameters describe a channel opening due to chirp effects, where, in different parts of the pulse, the probe photon can access different threshold harmonics. From studying channel opening by varying the energy of the photon in Fig. 6, we expect the errors of the LMA to be larger in this case because it predicts a more sudden channel opening than in QED, as η is varied. In Fig. 9 we indeed see that, at b = 1/4π , the QED result, which demonstrates sub-threshold pair creation, is noticeably different to the LMA prediction. However, as the chirp is increased, and the LMA describes pair creation above the n * = 1 threshold, the agreement between QED and simulation improves.

Focused lasers
In any real experiment, finite-size effects such as the laser pulse waist and focusing will be important in determining the yield of pair-creation events and in influencing particle spectra. Although no analytical form of an electron wavefunction in a typical focused laser background is yet known (but see e.g. [117][118][119]), we can employ a high-energy approximation [120][121][122] that averages over the probability for the process in a plane-wave background: where x ⊥ is the spatial co-ordinate transverse to the laser propagation direction, ρ(x ⊥ ) is the photon probe areal density and the plane-wave probability P is calculated from Eq. (6). The expression in Eq. (27) is then taken as the 'QED result' (although it also involves an approximation) and compared with simulation employing the LMA. The paraxial Gaussian beam [123] is taken to model the focused laser pulse, which has a potential of the form: where g(ϕ) is a pulse envelope, ς = z/z r , the Rayleigh range z r = ωw 2 0 /2, the carrier frequency of the laser pulse is ω, the beam waist is w 0 and the phase dependency is: . In order to use the high-energy approximation Eq. (27) for the QED calculation with the paraxial Gaussian beam, we employ the infinite Rayleigh-length approximation [124,125] (no such approximation is required in the simulation), which sends ς → 0, meaning that a ⊥ can be written as a ⊥ = exp[−(x ⊥ ) 2 /w 2 0 ]a ⊥ pw where a ⊥ pw is the plane-wave potential. We pick the same pulse envelope as in previous sections, i.e. g(ϕ) = cos 2 ϕ 2N for |ϕ| < N π and g(ϕ) = 0 otherwise.
As an example of the agreement between simulation and QED for the case of a focused background, a scenario is calculated where photons with energy parameter η = 0.2, arranged in a flat disc of radius 2 λ, collide head-on with a paraxial Gaussian beam with N = 16 and a waist of w 0 = 5λ. In Fig. 10, agreement in the energy spectrum, dP/ds is illustrated for two values of intensity, ξ = 0.5 and ξ = 2.5.
In the multi-photon regime of ξ = 0.5 (first column of Fig. 10), the centre of the spectrum (s = 1/2) agrees with the QED result to give an error dP/ds ≈ 3%(≈ 1/2N ). However, moving away from the centre of the spectrum, the disagreement between theory and simulation increases dramatically. One reason for this, which can be seen in the logarithmic plot (middle row), is that the QED result includes the linear Breit-Wheeler signal, which, being due to interference on the length scale of the pulse envelope, is missed by the LMA. The linear Breit-Wheeler signal has the effect of suppressing the rapid decay of the spectrum at larger values of s [116]. This disagreement between simulation and theory, whilst large, is occurring in a region which does not contribute significantly to the total probability of the process.
In the all-order regime of ξ = 2.5 (second column of Fig.  10), there is virtually no contribution from pulse-envelope interference, and the agreement between simulation and QED is good. For the 1D (plane wave) case, the agreement is dP/ds ≈ 3%(≈ 1/2N ), and for the focused case, the error is slightly larger, of the order of 5%. (The comparison is made up to s = 0.8 due to the finite number of harmonics included in the simulation rate, which leads to the rapid decay of the LMA away from the QED result, which can be seen in the logarithmic plot of the middle panel of the figure.) For comparison, the LCFA is also plotted, which performs worse in the 1D (plane-wave) comparison, but has a comparable error to the LMA for the 3D focused case at ξ = 2.5.

Conclusion
The process of nonlinear Breit-Wheeler pair creation exhibits a perturbative, but generally nonlinear (multiphoton), dependency on the field strength when the intensity parameter ξ 1, but an all-order (non-perturbative) dependency on the field strength when ξ 1. To model this process for experiments using a field with intermediate intensity ξ ∼ O(1), a simulation framework is required that is accurate across these regimes. The standard method of including strongfield QED effects in numerical simulations, by employing the locally-constant field approximation (LCFA), is insufficient for this purpose, as the errors of the LCFA grow rapidly as ξ is reduced from large values to values where ξ ∼ O(1) and smaller.
By directly comparing with numerical evaluation of exact QED expressions, the accuracy of simulations employing the locally monochromatic approximation (LMA) in calculating a range of observables and in different field configurations has been assessed. The accuracy can be quantified using an error function, E, which is defined to be the relative difference between the QED and simulation results, and has the same functional dependence on parameters as the full or differential probability. The error function is useful for defining a 'theory error' on simulation results.
Typically, in prediction of the total yield of pairs, the lightfront (energy) spectrum and the angular spectrum, E took values of the order of 1/ (where is the pulse duration), which is the expansion parameter of the LMA. Errors increase around harmonic channel openings, particularly around the n = 1 and n = 2 harmonics (which, for optical laser frequencies, require a probe photon energy of at least around 80 GeV and 40 GeV respectively, to access). Errors also increase when pulse-length interference effects begin to dominate, which can happen when it becomes favourable to create pairs via the linear Breit-Wheeler process using photons from a wide pulse envelope bandwidth. This effect is increased by: i) having a pulse shape with a wide bandwidth (in this paper the envelope was a squared cosine multiplied by a flattop); ii) having lower probe photon energies (this means the yield of pairs is suppressed more strongly when ξ is reduced in the multi-photon regime). In experiment, pulse envelope interference effects are likely to be minimal, because higher energy photons from the laser pulse will not be transmitted by various optical elements. (The pulse envelope chosen in the current paper has a particularly wide bandwidth.) A further source of error was found in the wings of pair energy spectra. When ξ is increased above ξ = 1, an increasing number of harmonics must be calculated in order that the decay of the spectrum in the wings is correctly predicted. Since only a finite number of harmonics can be tabulated in the numerical calculation, inevitably there are parts of the wings where the error can be large. However, this error occurs in regions of the spectrum that are orders of magnitude smaller than in the centre. Furthermore, as ξ is increased, the LCFA becomes more accurate, and could potentially be used in a hybrid approach (perhaps along with an improved LCFA approximation [38][39][40]), if resolution of these high-energy tails were crucial to an experiment. A comparison was also made for pair creation in a focused laser pulse, where the baseline error was found to increase to around 5% in regions where pulse-envelope and harmonic sampling effects were negligible.
Therefore, although parameters can be found that test the accuracy of simulations calculating strong-field QED effects using a locally monochromatic approach, these parameters are for laser pulses or probe particle energies that are currently beyond any experimental realisation.
One can compare the current study of calculating nonlinear Breit-Wheeler with a similar study [73] of nonlinear Compton. Whereas for Compton, harmonic structure in particle spectra can be accessed at low intensities and easily achievable electron energies, for Breit-Wheeler, clearly discernible harmonic structure requires very high probe photon energies, beyond anything currently planned in experiment. At intensity parameters ξ < 1, the LCFA is often cited as becoming inaccurate for Compton, because of the importance of the lowest harmonics, which are missed by the LCFA. However, for the Breit-Wheeler process at ξ < 1, the LCFA is much less accurate, and for experimentally relevant parameters, the threshold harmonic order is rather large. Thus for Breit-Wheeler, the inaccuracy of the LCFA is not due to missing harmonic structure, but rather the influence on the total probability due to interference effects on the length scale of the laser wavelength.
To conclude, in calculating how the yield of pairs depends on the intensity parameter ξ being varied from the multiphoton (ξ 1) to the all-order (ξ > 1) regime, simulations based on the LMA sustain an accuracy level of around 5%. This demonstrates the suitability of the approach for mod-elling high-energy experiments such as LUXE [49] at DESY and E320 at SLAC.