Longitudinal and transverse cooling of relativistic electron beams in intense laser pulses

With the emergence in the next few years of a new breed of high power laser facilities, it is becoming increasingly important to understand how interacting with intense laser pulses affects the bulk properties of a relativistic electron beam. A detailed analysis of the radiative cooling of electrons indicates that, classically, equal contributions to the phase space contraction occur in the transverse and longitudinal directions. In the weakly quantum regime, in addition to an overall reduction in beam cooling, this symmetry is broken, leading to significantly less cooling in the longitudinal than the transverse directions. By introducing an efficient new technique for studying the evolution of a particle distribution, we demonstrate the quantum reduction in beam cooling, and find that it depends on the distribution of energy in the laser pulse, rather than just the total energy as in the classical case.


Introduction
The emergence over the next few years of a new generation of ultra-high power laser facilities, spearheaded by the Extreme Light Infrastructure (ELI) [1], represents a major advance in the possibilities afforded by laser technology. In addition to important practical applications, these facilities will, for the first time, allow investigation of qualitatively new physical regimes. Among the first effects to be explored will be radiation reaction.
Radiation reaction-the recoil force on an electron due to its emission of radiationremains a contentious area of physics after more than a century of investigation. The standard equation describing radiation reaction (the so-called LAD equation, after its progenitors Lorentz, Abraham, and Dirac [2,3,4]) for a particle of mass m and charge q in an electromagnetic field F reads where f a ext = −qF a bẋ b is the Lorentz force. Here, the constant τ := q 2 /6πm is the 'characteristic time' of the particle ‡ and an overdot denotes differentiation with respect to proper time. Indices are raised and lowered with the metric tensor η = diag(−1, 1, 1, 1), and repeated indices are summed from 0 to 3. Theẋ-orthogonal projection ∆ a b := δ a b +ẋ aẋ b ensures thatẍ is orthogonal toẋ, preserving the normalisation conditionẋ aẋ a = −1 (equivalently the mass shell condition, p a p a = −m 2 , where p a = mẋ a = (γm, p)). We work in Heaviside-Lorentz units with c = 1.
Equation (1) may be unpacked and expressed in terms of the 3-momentum as dp dt = q E + p γm × B + τ γ 1 γ d dt γ dp dt where γ = 1 + p 2 /m 2 and dγ/dt = (p · dp/dt)/γm 2 . Despite numerous independent derivations of equation (1), either on the basis of energy-momentum conservation [4,5] or as the Lorentz force due to the particle's (regularized) self-field [6], it is subject to numerous difficulties; see the recent review [7] for an account of these problems and proposed solutions. The most widely used alternative to LAD is that introduced by Landau and Lifshitz [8], by treating the selfforce as a small perturbation about the applied force and retaining terms to leading order in the small parameter τ : (3) ‡ The characteristic time τ = 2r/3c can be interpreted as the time taken for light to travel across the classical radius of the particle, r = q 2 /4π 0 mc 2 . For an electron, τ e = 6.3 × 10 −24 s, corresponding to r e = 2.8 × 10 −15 m. Since radiation damping is proportional to τ , radiation reaction effects will typically be more prominent for electrons, for which q = −e and m = m e , than for particles with larger mass.
It is often claimed that (3) is valid provided only that quantum effects can be ignored, and though a rigorous demonstration remains elusive there is mounting evidence that this is indeed the case [9,10]. Note that equation (3) can also be presented in terms of the electric and magnetic fields, E and B, as [8,11] where the total time derivative acting on the E, B fields is Under the conditions expected at ELI, the caveat 'provided that quantum effects can be ignored' is pertinent. Quantum effects are typically negligible if the electric field observed by the particle,Ê, is much less than the Sauter-Schwinger field [12,13] typical of QED processes, that is provided where E S = m 2 e c 3 /e = 1.32 × 10 18 V/m is the Sauter-Schwinger critical field. For 1 GeV electrons in a laser pulse of intensity 10 22 W/cm 2 (parameters obtainable at ELI), χ ∼ 0.8 and quantum effects cannot be ignored. A complete QED treatment of radiation reaction is difficult to implement and problematic even to define but, provided χ remains small, a semi-classical modification to (3) should be valid [14].
An important difference between the classical and quantum pictures of radiation emission can be seen in the radiation spectrum. Classically, a charged particle can radiate arbitrarily small amounts of energy at all frequencies. However, in the quantum picture, the particle must radiate entire quanta of energy in the form of photons. Thus, the energy (frequency) of the emitted photons is limited by the energy of the particle. This suppresses emission at high frequencies, and introduces a cut-off in the spectral range of the emitted radiation [14]. As such, it is expected that the effects of radiation reaction are overestimated by classical theories in regimes where quantum effects become important [15], since they consider the particle to be radiating at all frequencies.
In order to account for this reduction in the effects of radiation reaction relative to the Landau-Lifshitz equation of motion, we follow Kirk, Bell and Arka [16] and scale the radiation reaction force by the function g(χ): The full expression for g(χ) involves a non-trivial integral over Bessel functions. To make this tractable, we use an approximation introduced by Thomas et al. [17], g (χ) = 1 + 12χ + 31χ 2 + 3.7χ 3 −4/9 .
It can be clearly seen that, in the classical limit χ → 0, we have g (χ) → 1, recovering the classical equation of motion (3). As we move into a more strongly quantum regime, the quantum nonlinearity parameter χ increases and the scaling function g(χ) decreases, in turn reducing the effects of radiation reaction. The model essentially reduces to a rescaling of the characteristic time of the particle, τ → g(χ)τ , which can also be applied to equation (4). For χ ∼ 1, the stochasticity of quantum emission becomes important, and the semi-classical model is no longer applicable [18]. At this point, g(χ) 0.18, which corresponds to a significant reduction in the effects of radiation reaction.
It is generally accepted that radiation reaction effects will be more readily observed in the behaviour of particles than in the radiation they emit [17,19]. As such, it is important to be able to accurately determine the distribution of a bunch of particles evolving according to (3) or its semi-classical extension (6). Usually this would involve solving a Vlasov equation [20] or following the evolution of very large numbers of particles, either of which is computationally very intensive.
In this paper we investigate beam cooling of a particle bunch due to classical and semi-classical models of radiation reaction. In Section 2 we present a detailed discussion of longitudinal and transverse phase space contraction of the particle distribution, along with an analytical solution of the classical Vlasov equation. The longitudinal particle distribution is introduced. Since the semi-classical Vlasov equation has no analytical solution, in Section 3 we introduce a new method of accurately reconstructing the particle distribution from the trajectories of a relatively small number of particles. Classical predictions using this method are compared to the analytical solution with excellent agreement. The method is then applied in Section 4 in order to compare classical and semi-classical predictions for an electron beam colliding with an intense laser pulse. Finally, we conclude by summarising our findings in Section 5.

Particle distribution and phase space contraction
The evolution of a particle beam can be described by the Vlasov equation for the particle distribution F (x, u), where u a = (γ, u) is the 4-velocity. Position and velocity are considered as independent phase space variables. The Vlasov equation for F can be expressed as d ds where V is the phase-space volume element and β s (x, u) describes the rate of change (i.e. expansion or contraction) of V with proper time s. (Technically, β s is the phase-space divergence of the vector field X = u a ∂/∂x a + A I ∂/∂u I (where A is the acceleration) associated with the flow d/ds, given by the Lie derivative L X V = β s V, see [20].) Capital Latin indices take three values. Unlike the Liouville equation (or the case with no radiation reaction) the phase-space volume element is not preserved by the flow, β s = 0.
To facilitate investigation of the interaction of a particle bunch with a laser pulse, we introduce the (null) wavevector k such that the phase of the pulse is The orthogonal (transverse) vectors , λ satisfying together with k and the null vector (defined to satisfy · = · λ = 0 and k · = −1) form a basis. In addition, the coordinates are also defined, along with the corresponding velocities u φ , u ξ , u σ and u ψ . However, u ψ is not independent and may be found from the normalisation condition u a u a = u 2 ξ + u 2 σ − 2u φ u ψ = −1. We note that Greek subscripts are used only as labels and are not free indices.
For a plane wave with arbitrary polarisation, the electromagnetic field tensor F depends on spacetime only through the phase φ, and takes the form where the functions a ,λ (φ) are dimensionless measures of the electric field strength in the ,λ direction. The corresponding electric and magnetic fields are E = (mω/q)[a (φ)ˆ + a λ (φ)λ] and B = k×E/ω, where the orthogonal unit 3-vectorsˆ ,λ satisfy k·ˆ = k·λ = 0.
In a similar manner, we assume that the particle distribution also depends on spacetime only through the phase φ, such that F (x, u) = F (φ, u φ , u ξ , u σ ). The Vlasov equation is then written where u I ∈ {u φ , u ξ , u σ }, and the accelerations A I ∈ {A φ , A ξ , A σ } follow from the singleparticle equations of motion. Dividing through by u φ , we have The quantity β is responsible for any phase space contraction (β < 0) or expansion (β > 0) of the particle distribution, and the associated change in electron entropy [21]. For a highly relativistic particle beam colliding with a laser pulse (the scenario in which radiation reaction effects are most prominent), we are mainly interested in the dependence of F on φ and u φ . An advantage of the coordinate system (9)-(11) is that it decouples the longitudinal from the transverse velocity in the Lorentz invariant measure, d 3ẋ /γ = du ξ du σ du φ /u φ . Hence we can define the longitudinal distribution which satisfies the reduced Vlasov equation Here, β describes the longitudinal phase space contraction. The transverse contribution is then Note that this reduction to the longitudinal distribution is purely a consequence of the coordinate system, and does not rely on the plane wave assumption. It is at this point that a decision must be made as to the appropriate single-particle equations of motion. While there are many classical models for radiation reaction [7], we start by considering the Landau-Lifshitz equation given by (3), before moving on to the semi-classical extension (6). This is in part motivated by the existence of an analytical solution to the single-particle Landau-Lifshitz equation [22]. In our coordinates, the Landau-Lifshitz equations in the plane wave (12) arê where a 2 (φ) = a 2 (φ) + a 2 λ (φ) and prime denotes differentiation with respect to φ. Inserting these equations into (16) and (17), we find for the classical casê It is immediately apparent that half the contraction of the distribution occurs in the longitudinal and half in the transverse directions. The semi-classical equations of motion are just (18) with the replacement τ → g(χ)τ . However, since χ(φ, u φ ) = 3τ a(φ)u φ /2α (where α is the fine structure constant) depends on u φ (but not on the transverse velocities) we pick up an additional contribution to the longitudinal phase space contraction: whereas, the transverse contraction is simply scaled by g(χ): Thus, as quantum effects become more important and g(χ) decreases, the semiclassical model predicts a reduction in both the longitudinal and transverse phase space contraction (reduced beam cooling). As well as this scaling of the classical contraction by g(χ), there is an additional longitudinal heating given by where measures the strength of the longitudinal compared to the transverse phase space contraction. This is shown in Fig. 1(a) for the interval χ ∈ [0, 1]. Even for the weakly quantum regime in which the semi-classical model remains valid, we observe a significant reduction in longitudinal beam cooling. This is especially clear when comparison is made with the classical resultβ as shown in Fig. 1(b). We see that where χ = 0.2 there is nearly a 60% reduction in the longitudinal contraction experienced compared to the Landau-Lifshitz model. For the case of the classical Landau-Lifshitz theory in a plane wave, the Vlasov equation (14) may be solved analytically for the particle distribution: where φ 0 , u 0 φ , u 0 ξ , u 0 σ are the initial phase and velocities of a particle with {u φ , u ξ , u σ } at phase φ. In a similar manner, the longitudinal distribution is found to be The contraction/expansion of phase space is contained in the function Solutions to equation (18) [22] can then be used to rewrite φ 0 , u 0 φ , u 0 ξ , u 0 σ in terms of the independent variables {φ, u φ , u ξ , u σ }: where the functions with i ∈ { , λ} depend only on the properties of the laser pulse.
Using equations (27) we can express u φ (ϑ) in equation (26) in terms of the independent variable u φ , such that the distribution becomes or for the longitudinal distribution This latter result is in agreement with observations made by Neitz and Di Piazza [23], and we see that the longitudinal distribution is only sensitive to the properties of the laser pulse through the function G(φ). After the pulse has passed, G becomes constant and is proportional to the fluence of the pulse. Final-state properties of the longitudinal distribution therefore depend only on the total energy contained in the laser pulse, and are insensitive to how that energy is distributed within the pulse. The full distribution F , on the other hand, depends additionally on the integrals W i given in equation (28). Although the reduced Vlasov solution (31) is somewhat simpler than (30), and captures the key features of the electron beam itself, the solution is not sufficient to calculate the transverse current density, and hence cannot be coupled to Maxwell's equations to determine the radiation produced by the electron beam. However, if the transverse momentum spread is sufficiently small, we can approximate the full distribution by where the δ-functions restrict the transverse velocities to the submanifold X i . Then, in addition to (16), equation (14) yields Note that (33) indicate that the distribution is concentrated on a surface in phase space that itself satisfies the Landau-Lifshitz equation. Given solutions to the reduced Vlasov equation (31) and the transverse Landau-Lifshitz equation (33), the current can be written with j a ⊥ and evaluated as We could also calculate j a directly, but it follows more straightforwardly from charge conservation, ∂ a j a = 0. In the following, we restrict our attention to the longitudinal distribution f (φ, u φ ), and longitudinal beam cooling, as this is more readily measurable in experiments than the transverse cooling. However, the transverse cooling, which can be considerably greater, can be determined from equation (23).

Numerical (re)construction of the particle distribution
The motion of a single charged particle colliding with a laser pulse, including radiation reaction, has been extensively studied [9,22,24,25]. As shown in Section 2, the Vlasov equation with classical readiation reaction can be solved analytically. However, this is not the case for the semi-classical model (6) or for stochastic models of radiation reaction in the quantum regime. Instead of attempting to solve a Vlasov-type equation on the phase space numerically, which would require significant computing resources, we propose an innovative method which allows for the dynamics of a particle distribution to be explored using single-particle equations of motion such that the distribution can be efficiently reconstructed. While this approach is quite general and could be used for a variety of systems, here we consider a distribution of particles subject to equation (3) and its semi-classical extension (6), without particle-particle interactions §.
Assuming that the laser pulse can be approximated by a plane wave with compact longitudinal support , any spatial spread in the initial particle distribution would only determine the moment when each particular particle enters the pulse. For simplicity, we therefore take all particles to originate from the same point. This is reasonable as we are primarily interested in the longitudinal momentum distribution.
Since our pulse is modelled by a plane wave and we focus on the longitudinal properties of the distribution, we consider the initial momenta to be strongly peaked about zero in the transverse directions. As such, the initial distribution can be taken to be a Maxwellian distribution for the (longitudinal) momentum p (in units of mc) with φ = ωt − k · x the phase, θ the variance of the distribution, and N P the number of particles. The momentum p is related to our velocities of Section 2 by We stress that this initial distribution is chosen for its simplicity; alternative distributions could be used where appropriate (such as Maxwell-Jüttner).
Typically, one would sample the distribution at random, which would require a large number of particles to accurately represent the distribution. Instead, since the particle number is simply we determine the momentum spacing δp between the particles from the initial distribution by truncating the integral in (38) so that the particle number increases by unity in the given momentum interval: (39) § For a highly relativistic particle bunch, these interactions can be neglected on the time scale of the laser interaction.
A function has compact support if it is zero outside a finite interval.
This leads to a set of N P = 2N c + 1 initial momenta V (0) = p i (0) for i ∈ [−N c , N c ], with the p i generated iteratively from p 0 =p and p ±1 =p ± 1/f (0,p) using The momentum space is not sampled uniformly, instead more particles are located in regions where the distribution function is large.
As the evolution proceeds, this procedure is applied in reverse to reconstruct the distribution. The set of momenta V (φ) is ordered such that p i+1 ≥ p i and used to find Reconstruction of a distribution from a particle sample can be problematic, but in our formalism it becomes quite natural. This is achieved by using the momentum spacing between particles to determine the value of the distribution such that equation (39) is satisfied for all φ (i.e. integration over each of the measured momentum spacings always contributes a single particle to the total particle number). The closer the measured momenta are together, the 'more likely' one is to have a particle in that momentum range, resulting in a larger value for the distribution. The definition (41) allows properties of the distribution to be calculated directly from the momenta of the individual particles. The mean is simply evaluated as In a similar manner, higher-order moments of the distribution X n may be calculated straightforwardly as: For example, the variance θ(φ) = X 2 (φ) and skewness S(φ) = X 3 (φ)/X 3/2 2 (φ). We note that angular brackets · · · will be used to denote an average over the bunch, not time.
A major advantage of this new methodology is that the distribution can be accurately represented and reconstructed using far fewer particles than would be required using random sampling. Figure 2(a) shows a simple Gaussian distribution with zero mean and unit variance, A(z) = (1/ √ 2π) exp(−z 2 /2), reconstructed using two methods. First, the distribution A(z) was sampled according to the iterative relation  (40) and (41) with N P = 401 particles (red, solid) and compared to the standard approach of random sampling, using both N z = 401 (blue, dot-dash) and N z = 100, 000 (black, dashed) particles. Part (b): Variation of the measurement of the initial relative momentum spreadσ i as the number of particles N P is increased, compared to the desired value,σ d .
(40) with N P = 401 and then reconstructed using (41) (-). Next, N z = 100, 000 random numbers were generated from A(z) using a Mersenne-Twister pseudo-random number generator, from which the distribution was found using 100 fixed-width bins (---). Even with 100,000 random samples, the standard approach does not describe the Gaussian perfectly, while the new method presented in this paper does an excellent job using only 401 points. For comparison, a sample of N z = 401 was also reconstructed using 50 fixed-width bins (-· -). The advantage of the new method with such a small sample size is clear. Figure 2(a) also shows how sampling the distribution according to equation (40) cuts off the tails of the distribution, which will have an effect on the measured properties of the reconstructed distribution. The finite number of particles used to represent the distribution causes the measured relative momentum spread (defined by equation (44) below) to be less than that specified when defining the initial distribution. Essentially, it comes down to the ' ' in equation (39), compared to the definition given by equation (41). As N P is increased, the approximation in equation (39) improves, and the measured value approaches the desired value. In addition, with more particles the distribution is sampled further into the tails. Figure 2(b) shows the measured initial spread,σ i , as a fraction of the desired spread,σ d , when the particle number is varied from as few as 11 up to N P = 2001. We see that the approximation quickly improves as N P is increased up to about 500. The value N P = 401 is chosen to give less than 0.5% error in the initial measured momentum spread. In practice, good agreement can be found with lower N P , with the caveat that properties sensitive to the tails of the distribution (such as the skewness) may be strongly affected. (This is confirmed in Fig.  3, where the distribution and its statistics calculated from the analytical solution to the classical Vlasov equation are compared to those obtained with this new method.)

Interaction of a particle bunch with high-fluence laser pulses
The analytical solution to the Vlasov equation including radiation reaction according to the Landau-Lifshitz theory derived in Section 2 predicts that the collision of an energetic electron beam with a high-intensity laser pulse leads to a significant contraction of the particle phase space, resulting in a reduction in the relative momentum spread. Agreement between this analytical solution and numerical results obtained using the approach discussed above is shown below to be excellent. As previously observed, classical beam cooling depends only on the total fluence of the pulse, rather than its duration or peak intensity independently [23,26]. However, for the semi-classical extension, the Vlasov equation is no longer tractable. This highlights the value of our approach and, to demonstrate the use of our proposed method in such a case, we consider the importance of quantum effects in the interaction of an electron bunch with a high-intensity laser pulse.
To establish the impact of quantum effects on the evolution of the particle distribution subject to radiation reaction, we introduce the relative momentum spread and the momentum skewness (calculated from the meanp and variance θ): The former gives a measure of the beam quality, while the latter indicates how symmetric the distribution is about its mean. We restrict our attention to a linearly polarised N -cycle plane wave pulse (12), modulated by a sin 2 -envelope [9], with a λ = 0 and a = a(φ), where where a 0 is the dimensionless (peak) intensity parameter (the so-called "normalised vector potential") and L = 2πN is the pulse length ¶. This pulse shape offers compact support, allowing the particles to begin and end in vacuum. The total fluence (energy per unit area) of the pulse is proportional to In this work, E is kept constant, which fixes a 0 for each N . It has been shown [23,26] (see also Section 2) that the classical Landau-Lifshitz prediction for the final state of a particle distribution emerging from the pulse is completely determined by the fluence, whereas quantum effects are expected to depend on the value of a 0 itself. We are then able to explore the impact of the reduced emission in the quantum model with varying a 0 while maintaining the same classical prediction. This allows us to explore the relative importance of quantum effects.
To motivate this study, parameters have been chosen to be relevant at the forthcoming ELI facility. We have chosen to consider N a 2 0 = 9.248 × 10 3 which, for N = 20 with a wavelength of λ = 800 nm, represents a full-width half-maximum pulse duration of 27 fs with peak + intensity 2 × 10 21 W/cm 2 . We have investigated pulses of length N ∈ [5,200] cycles (together with their corresponding a 0 ) counter-propagating relative to a bunch of N P = 401 particles, with an initial momentum spread of 20% around 1 +p 2 = 2 × 10 3 . This corresponds to an average particle energy of just over 1 GeV, which should be well within the capabilities of the laser-plasma wakefield accelerator at ELI.
Before comparing predictions of the classical and semi-classical models, we briefly confirm the validity of our method by comparing numerical results with the analytical solution (31) obtained in Section 2. Figure 3 shows the interaction of a 1 GeV electron beam with a plane-wave laser for two pulse lengths, N = 20 (a 0 22) in part (a) and + Peak intensity is obtained from I peak = (4π 2 m 2 e c 5 0 /e 2 λ 2 )a 2 0 2.74 × 10 10 (a 0 /λ) 2 W/m 2 . N = 5 (a 0 43) in part (b). The left-hand panels show the numerical results obtained using the new method described in Section 3, while the right-hand panels show the solution (31) using the initial distribution f (0, u φ ) corresponding to f (0, p) given by (36) with p = 1 2 (u φ /ω − ω/u φ ). The momentum p is evaluated during the evolution using (37), along with the solutions (27) for u ξ , u σ when u 0 ξ = u 0 σ = 0. The agreement is excellent. Note that the measured values for the initial and final momentum spread also agree, while the skewness is underestimated by the numerical method (as discussed in Section 3). Figure 4 shows the variation of the particle distribution on the (φ, p) phase space. As can clearly be seen in moving from the classical Landau-Lifshitz theory (left) to the quantum model (right), there are noticeable differences in the meanp, spreadσ, and skewness S of the distribution. We first note that the deficit in measuring the initialσ i = 19.9% < 20% is due to the finite number of particles used to represent the distribution (as discussed at the end of Section 3 and illustrated in Fig. 2(b)).
For the classical theory, the final distribution only depends on the fluence of the pulse, though this does not prevent the system from taking different routes along the way. As the number of cycles is decreased, very different intermediate behaviour is observed in Fig. 4, yet the measured properties of the final distribution support this prediction: in each case, we measure the mean momentump f = 1197.7 with a relative spreadσ f = 12.5%. This represents a significant contraction of the phase space, where the average energy of the particle bunch decreases significantly, as does its thermal spread (beam cooling), and the distribution becomes more sharply peaked. In addition, we find the development of a negatively-skewed distribution with S f = −0.46. In the classical model this is readily understood, since the higher a particle's momentum the more it radiates. This causes particles in the positive tail of the distribution to be slowed down more than those in the negative tail.
The introduction of a semi-classical model in which the effect of radiation reaction is reduced by the function g(χ) given by equation (7) results in a reduction in the amount of phase space contraction. Figure 4(a) for N = 200 clearly demonstrates this, with the final average momentump f = 1301.1 only slightly higher than the classical case. The final relative momentum spread is now 14.4%, showing that the final distribution is less sharply peaked. While remaining negative, the skewness reduces in magnitude to −0.280, because it is precisely the higher-energy particles (which were classically most affected by radiation reaction) that now have this damping suppressed due to larger χ (smaller g(χ)).
These changes become more pronounced as we move to higher intensities (by reducing the number of cycles). For N = 20, as shown in Fig. 4(b), we find that p f = 1451.8 andσ f = 16.6% have both increased, with the skewness also increasing to S f = −0.129. This trend continues to N = 5 as shown in Fig. 4(c). In this case, very little beam cooling occurs for the quantum model, with the final relative momentum spread taking the valueσ f = 18.0% aroundp f = 1581.3. The profile also remains much more Gaussian, with S f = −0.0597. The reduction in phase space contraction (beam cooling) observed here is in agreement with previous predictions [27], which provides further validation of our method for reconstructing the particle distribution. Using this new method, it has been possible to investigate the effects of the semi-classical model on a distribution of particles. The distributions are nicely reconstructed and do not feature any artefacts, in contrast to other approaches [17], which emphasises the power of our method. Figure 4 also shows how the difference between the classical and quantum results increases as the intensity is increased (or as N is decreased). It is therefore interesting to consider the difference δσ f =σ qm f −σ cl f as a fraction of the total (constant) classical change in momentum spread, ∆σ cl =σ i −σ cl f . This can be found in Fig. 5(a), where we see that for N = 5 the two predictions differ by about 75%. As N is increased, this ratio is reduced because the average quantum parameter χ becomes smaller and radiation reaction is not so heavily suppressed. It would be expected that the two models converge as N → ∞.
In cases where there is a large discrepancy between the predictions of the two theories, it is especially important to be confident in the validity of the model. As a semiclassical model, we expect it to remain valid into the weakly quantum regime, such that particles experience instantaneous values χ 2 1. Figure 5(b) shows the evolution of the bunch-average χ as the bunch moves through a laser pulse with N = 20 according to the classical and semi-classical models. Initially, there is good agreement between the two models, until the bunch approaches the centre of the pulse, where the intensity becomes higher and the models are significantly different. For completeness, we note that our highest intensity case with N = 5 satisfied χ 2 < 0.22.

Conclusions
The next few years will see the emergence of a number of new high-power laser facilities operating at unprecedented field strengths, providing access to fundamentally new physical regimes. This will allow us to experimentally probe previously untested areas of physics, such as the long-standing question of radiation reaction.
In this paper, we have analysed the transverse and longitudinal cooling of a relativistic electron beam as it interacts with an intense laser pulse, according to classical and semi-classical theories of radiation reaction. In the classical theory, we have found these two contributions to be equal, but quantum effects break this symmetry, leading to significantly less cooling in the longitudinal than the transverse directions.
To facilitate evaluation of the longitudinal beam cooling effects, we have introduced an innovative method to efficiently and accurately calculate the distribution function for an electron beam interacting with an intense laser pulse. This has been validated by comparison with an analytical solution to the Vlasov equation in the classical case, and used to compare classical and quantum predictions of radiative cooling. We have found that quantum effects can significantly alter the beam properties and, unlike the classical case, can be influenced by the shape of the laser pulse, not just its energy.
As we move into the quantum regime where final-state electron beam properties become sensitive to pulse shape, it is becoming increasingly important to have an efficient method in order to investigate the full parameter space. The approach developed here to facilitate this study of beam dynamics provides a powerful tool with wide-ranging application within the discipline.
The results presented in this paper are limited to the semi-classical case χ 2 1. However, it should be noted that, for the longitudinal beam cooling, this restriction is due to the use of a deterministic equation of motion, and not the method of sampling and reconstructing the distribution. There should be no obstruction to exploring more strongly quantum regimes (such as higher initial beam energies ∼ 5 GeV available at ELI) using this approach with a stochastic equation where photon emission probabilities are determined by strong field QED, as in [28,29]. This will be addressed in future work, along with an investigation of stochastic transverse beam cooling.