Numerical modelling of Compton scattering in ultra-intense laser pulses

We examine current methods of numerically implementing Compton scattering in the context of intense laser-matter interactions. In a recent publication [1] it has been shown that a commonly used approach generates the correct spectra in nearly all cases, except those when the harmonic structure is important. Here we provide an explanation for this using an alternative, classical argument.


Introduction
Recent technological advances have led to the development of lasers with unprecedented powers and intensities. The current record was set by the HERCULES laser in 2008, when it achieved a peak focal intensity of 2 × 10 22 W/cm 2 [2]. This upward trend is expected to continue into the future, with numerous facilities being developed to provide ultra-intense laser fields. These include the Vulcan 10 PW upgrade [3] which is expected to provide peak intensities of 10 23 W/cm 2 . In the longer term we look forward to the European "Extreme Light Infrastructure" (ELI) Facility [4] and the XCELS project [5], both of which may surpass this intensity by a further order of magnitude.
These developments have led to a renewed interest in probing strong-field quantum electrodynamics (QED) using high-intensity laser fields [6,7,8] as well as in a number of applications [9,10]. However, the new generation of facilities will achieve such ultra-intense conditions, in part, by a strong focusing of the beam. So far, the majority of analytical QED calculations are limited to plane-wave models which do not capture the structural features of strongly focussed fields. Additionally, multiple QED processes can take place while the particle is in the field, even leading to runaway cascading. Describing such situations analytically is almost impossible and so one is forced to resort to numerical methods (However, see the recent paper [11]). The standard numerical techniques for modelling QED processes in a laser field involves propagating a particle classically thorough a field and then using statistical event generators to determine various QED transition processes [12,13,14].
If we are to rely on these numerical methods, we must first verify their accuracy and gain a full understanding of their limitations. This is best done by considering problems that can also be solved analytically to which explicit comparisons can be made. Such problems are, however, few in number. In this paper we will study the accuracy of these numerical techniques, restricting ourselves to one such process, the Compton scattering of an electron, for which an analytic solution can be obtained. We outline the differences between the exact analytical calculation of the emission spectra and those obtained via numerics based on discrete time approximations. We also present a classical argument explaining why discrepancies occur.

Outline of numerical techniques
While all of the codes used by the high-intensity laser community have their differences, they are all designed with the same general principles at their heart (see, for example, [15]). For the purposes of this study we consider the single-particle QED program SIMLA [16,17], recently developed by the authors, whose operation we now briefly describe.
Working in natural units where = c = 1, we describe the laser field using the dimensionless intensity parameter a 0 = eE/ωm, where e is the electron charge, m the rest mass, E the peak amplitude of the electric field and ω the laser frequency. The importance of quantum effects can be measured by considering the (invariant) quantum efficiency parameter χ e ≡ (F µ ν p ν ) 2 /m 2 ∼ γE/E cr , where E cr = 1.3 × 10 16 Vcm −1 is the QED 'critical' field ('Sauter-Schwinger' field) [18,19,20]. This parameter is equal to the work done by the laser field on the particle over the distance of a Compton wavelength. Thus, when χ e > ∼ 1 quantum effects dominate and processes such as vacuum pair production occur. For the purposes of this work we will restrict our attention to the regime where a 0 1, χ e ≤ 1, such that quantum effects will play a role, but at the same time we can neglect processes such as pair production. In the high-intensity limit a 0 1 the size of the radiation formation region is of the order λ/a 0 λ, where λ = 2π/ω is the laser wavelength [21]. Thus the laser varies on a scale much larger than the formation region and so can be approximated as locally constant and crossed, allowing us to determine the probability of photon emission using the expression for the constant field rate Γ [21] where α is the fine structure constant, K ν is the modified Bessel function of order ν, η ≡ χ γ /χ e , χ ≡ 2η/ [3χ e (1 − η)], and χ γ ≡ (F µ ν κ ν ) 2 /m 2 is the analogous quantum efficiency parameter for a photon emitted photon with momentum κ ν .
Note that although dΓ diverges at small χ γ , the total differential probability of photon emission (i.e., of any χ γ ), dW = Γdt, where Γ ≡ χe 0 dΓ, is finite (see also [22]). It is common to model the trajectory of the electron through the field classically and then to correct its momentum after every time step using statistical routines to calculate the photon emission rates using the constant field rate. It is important to note that this is not equivalent to what happens in a full, analytical calculation. In normal QED calculations scattering amplitudes are determined using asymptotic 'in' and 'out' states (see e.g., [23]). No assumptions are made about the electron's trajectory in the laser. One simply begins with an in-state momentum at t = −∞ and calculates the out-state momentum at t = ∞ without determining how the electron moves in between. In contrast, in the numerical method we are continuously tracking the electron orbit.
The statistical routines that determine the photon emissions are implemented in the following manner. At each time step a uniform random number r ∈ [0, 1] is generated, and emission deemed to occur if the condition r ≤ Γdt is satisfied, under the requirement Γdt 1. Note that during the simulation dΓ (and thus Γ) is a time-dependent quantity owing to the effect of the temporally varying laser pulse and electron motion. Given an emission event, the photon χ γ is determined as the root of the sampling equation ζ = Γ(t) −1 χγ 0 dΓ(t), where ζ is a uniform random number ζ ∈ [0, 1] (In practice we implement an infra-red cut-off so that the code does not include the emission of large numbers of low energy photons; the integral is performed from a lower limit ε ∼ 10 −5 , rather than zero. The emission of soft photons of energy below this cut off does not appreciably affect the electron dynamics (see e.g., Ref. [24]).). Next, we calculate the photon momentum from χ γ assuming that the emission is in the direction of motion of the electron. This is valid for γ 1, since in reality the emissions will be in a cone of width γ −1 [25,26]. Finally, the photon momentum is subtracted from the electron momentum and the simulation proceeds by propagating the particle via the Lorentz equation to the next time step.

Comparisons
We now present a comparison of the numerical photon spectra produced by the SIMLA program with that obtained via the full analytic calculation. The emission spectra can be calculated exactly for an electron in an infinite plane-wave field, details of the calculation can be found in Refs. [27,28,26]. Although we can't run a numerical simulation for an infinite plane-wave field, we can closely approximate such a field using a plane wave with a super-Gaussian temporal profile. The emission spectra in such a field closely resembles that of an infinite plane wave, retaining most of the structural features [29]. In figure 1 we show the emission spectrum calculated analytically (for an infinite plane wave) together with that calculated numerically for an electron in a degree-8 super-Gaussian field. The laser intensity a 0 = 20 with wavelength λ = 0.8 µm and the initial electron γ 0 = 9000. The parameters have been deliberately chosen to emphasis certain structural features in the spectrum. It can be seen that while in general the numerical spectrum matches the analytical spectrum very well, it fails to reproduce certain structural features (caused here by the first and second harmonics). It turns out to be true in general that the numerical method reproduces the correct spectra except that it averages over the structural features. A detailed analysis for a wide range of parameter regimes is presented in Ref [1].

Classical analysis
In these proceedings we explain the discrepancy between the full analytically calculated spectra and that obtained using the numerical method by employing a classical argument.
The expression describing the classical radiation emission of a particle can be found in most textbooks on electrodynamics (see, e.g., [25]). Here we adopt the less common covariant formulation which can be found, for example, in [30]. The four-momentum of the emitted radiation can be written where sgn is the sign function, k µ = ω (1, sin θ cos φ, sin θ sin φ, cos θ) is the vector of the emitted radiation, where θ and φ are the polar angles that determine the direction of the photon, and is the Fourier transform of the electron current, where u and x are the four-velocity and fourposition of the electron, respectively. However, the boundary conditions can result in an infra-red divergence which means it is advisable to regulate the integral: see [31] for a fuller discussion. The resulting expression is Performing an integration by parts we find which is (4) with the addition of some boundary terms. Now, returning to our calculation, the radiated energy is given by the zero component of (2) where is the classical spectral density, which gives the 'number of photons' N γ radiated per unit frequency per unit solid angle. It is possible to formulate a classical equivalent of the QED numerical scheme by dividing the particle orbit into discrete segments, evaluating (6) for each segment, and then adding the resulting spectra together to obtain the total spectrum. Proceeding this way, the electron four-current (4) becomes Clearly, dividing the integral up into segments and summing the result from each segment as such will not alter the final result of the current. However, when it comes to calculating the spectral density (7) this is not the case, as, of course, since, by calculating the spectral density in stages using each discrete segment of the four-current at a time, we are in effect throwing away the cross terms. The evaluation of (6) is a very numerically expensive procedure, especially for large a 0 . We are aided in this respect by the fact that the program SIMLA uses a nonlinear time grid. (Other techniques to improve efficiency can be found in refs. [32,33].) However, we do not have the same constraints in the classical modelling as we do in the QED case. This is because we have the full orbit (and, by implication, the full properties of the field) over each segment, so it is not necessary for us to make any constant field assumptions. We just require that the classical 'formation length' (see, e.g., [34,32]) is shorter than the size of the discretisation. Thus we can perform our analysis with a much lower a 0 than is valid for the QED method.
In figure 2 we show the classical emission spectra, calculated by dividing the orbit into various numbers of segments. It can be seen very clearly that, as we divide the orbit, the neglect of the cross terms in (10) causes a smoothing of the structural features of the spectrum, just as we saw in the QED case. It appears that the discretisation causes an averaging of the emitted energy. To prove this we consider the total radiated energy as given by Larmor's formula [25] P = − 2 3 e 2 m 2 dτ dp µ dτ dp µ dτ .
Since (11) is a function of the acceleration only, it is possible to split the integral up, evaluate each segment separately, and add the parts together without changing the result. This means that the total radiated power is conserved by the discretisation, even though the spectral density is changed. Hence the discretisation averages over the structural features of the spectrum.

Summary
We have analysed the commonly used numerical technique based upon the idea of splitting a classical electron orbit into discrete segments, and then correcting the momentum over each of these segments using a QED-derived photon emission rate. In a previous paper [1] it has been shown that this method reproduces the correct spectra (calculated analytically) reasonably well, except that it fails to reproduce certain structural features of the spectra. By adopting an analogous classical approach we have shown that, in dividing the orbit up into discrete segments, we neglect cross terms which provide information about the field. The shorter these segments become the more information is thrown away and the flatter the spectrum becomes. However, the total radiated energy remains constant, meaning that the result of discretisation is an averaging over the structural peaks and troughs of the spectrum.