Quantum radiation reaction in head-on laser-electron beam interaction

Near future facilities [1–3] will provide extreme laser intensities ($I\gt {10}^{22}\,{\rm{W}}\,{\mathrm{cm}}^{-2}$), where quantum effects such as electron–positron pair production and discrete photon emission might play a central role in laser-matter interaction [4–12]. Previously, electron–positron pairs have been produced in experiments using a moderately-intense laser of intensity $I\sim {10}^{18}\,{\rm{W}}\,{\mathrm{cm}}^{-2}$ counter-propagating with the ultra-relativistic (46 GeV) SLAC electron beam [13–15]. This setup takes advantage of the ultra-relativistic energy of the particles to observe certain nonlinear quantum effects in electromagnetic fields whose amplitude remains several orders of magnitude below the critical Schwinger field [16], which constitutes usually the threshold to observe pairs spontaneously created in vacuum. As detailed in [17], the field magnitude in the rest frame of the particles will then be of the order of the critical field and the probability of the process becomes optimal. By leveraging on the tremendous progress accomplished in laser technology in the last decades, one can also envisage nowadays to decrease the energy of the relativistic particles and increase proportionally the magnitude of the field of the laser. This explains the recent growing interest [8, 18–23] on configurations where a relativistic electron beam interacts with laser pulses of significantly higher intensity than in the SLAC experiment. The typical electron energy sufficient to diagnose nonlinear quantum effects is around a few GeV and such electron beams can now be generated from an all-optical source in an efficient manner: the current experimental record for self-injected electrons obtained in a laser wakefield accelerator is 4 GeV [24]. This shows that the near future laser facilities can be used to explore this nonlinear Compton scattering configuration without the aid of conventional accelerators.

In our previous work [25] we have studied the radiation reaction for an electron beam in a laser field with the use of the Landau–Lifshitz equation [26] which has been recognised as the best candidate to describe classically the effect of radiation reaction on charged particle orbits [27–32]. In the Landau–Lifshitz equation framework, charged particles emit radiation continuously and the direct effect of this emission can be represented as a continuous drag force in the particle motion equation. As shown in [25], when a GeV electron beam collides head-on with an intense laser ($I\simeq {10}^{21}\,{\rm{W}}\,{\mathrm{cm}}^{-2}$), one of the key effects of classical radiation reaction is to reduce the width of the electron energy distribution function during the interaction [8, 33]. The tendency of classical radiation reaction to shrink the electron energy distribution function has been studied also in fusion plasmas [34], and in ion acceleration with solid targets [35]. If the laser intensity is raised such that the classical description of radiation reaction becomes inapplicable, the quantum effects in radiation reaction induce the opposite behaviour, leading to an increase of the energy spread of the beam spectra [36]. With the advent of quantum electrodynamics (QED) modules incorporated in the traditional particle-in-cell (PIC) algorithm, we are now able to simulate from first principles quantum radiation reaction in laser-plasma interaction and therefore to validate some of the recent theoretical predictions. This paper thus deals with differences between the classical and the QED description in the transition regime where the probabilities for pair creation are still negligible, but quantum effects in the photon emission can already be significant ('moderately quantum regime' defined in [37]). This is of particular relevance since upcoming experiments at several facilities will be able to operate in this regime. We carry out PIC-QED simulations that allow us to evaluate the influence of quantum emission on the electron energy spread and the divergence of the electron beam. Maximum attainable energy spread due to quantum stochasticity as a function of mean electron energy and the laser intensity is computed. This result further allows us to obtain a semi-classical analytical prediction for the electron energy spread after the shutdown of the laser as a function of the initial beam and laser parameters. The prediction is in agreement with fully quantum Monte-Carlo PIC simulations.

This paper is structured as follows. In section 2, we introduce the QED framework to describe photon emission. We then analyse, in section 3, the evolution of the electron energy spectra, predicted analytically, and we compare the analytical results with QED-PIC simulations. In section 4, we study the evolution of the electron beam divergence, another measurable quantity in these scenarios, and in section 5 we state the conclusions of this work.

In QED, radiation is a discrete stochastic process and this impacts the particle trajectory in a distinct manner from the continuous emission associated with the classical radiation processes. The probabilities of the various processes in an electromagnetic plane wave are based on Volkov states [38] where the quantum-transition probability is evaluated taking into account the interaction between the particle and the background wave. In the event of emission, there is a transfer of energy from the electron to the emitted photon; otherwise, the electron momentum and energy remain unaltered. The classical limit corresponds to the case where a large number of photons, whose energy remains small compared to the electron energy, is radiated: the high frequency of the emission events allows the approximation of the trajectory as a classical trajectory with a continuous drag. The main difference between the classical and the QED approach is that QED accounts for the possibility of emitting high-energy photons even in a setup where the cross-section for Compton scattering is small (i.e. the average energy loss of the particle is negligible). In the quantum regime, the stochastic nature of emission becomes noticeable and one may expect a diffusion in energy around the mean value of the energy loss, as it was first reported in [5, 36].

The total probability of radiation emission by a single particle is relativistically invariant and depends on the normalised vector potential ${a}_{0}={eE}/({mc}{\omega }_{0})$ and the quantum invariant parameter χ (${\chi }_{e}$ for electrons and ${\chi }_{\gamma }$ for photons) defined as:

$$\begin{eqnarray}&&{\chi }_{e}=\displaystyle \frac{\sqrt{{({p}_{\mu }{F}^{\mu \nu })}^{2}}}{{E}_{s}\,{mc}},\quad {\chi }_{\gamma }=\displaystyle \frac{\sqrt{{({\hslash }{k}_{\mu }{F}^{\mu \nu })}^{2}}}{{E}_{s}\,{mc}},\end{eqnarray} \tag{ 1 }$$

where p_μ is the particle 4-momentum, k_μ is the photon wave 4-vector, ${F}^{\mu \nu }$ the electromagnetic tensor, ${E}_{s}={m}^{2}{c}^{3}/{\hslash }c$ the Schwinger critical field [16], m is the electron mass, e is the elementary charge, c is the speed of light and ${\omega }_{0}$ is the frequency of the electromagnetic wave. The differential probability rate of photon emission by nonlinear Compton scattering is then given [17, 39–43] by

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}P}{{\rm{d}}t\,{\rm{d}}{\chi }_{\gamma }}=\displaystyle \frac{\alpha {{mc}}^{2}}{\sqrt{3}\pi {\hslash }\gamma {\chi }_{e}}\left[\left(1-\xi +\displaystyle \frac{1}{1-\xi }\right){K}_{2/3}(\widetilde{\chi })-{\int }_{\tilde{\chi }}^{\infty }{\rm{d}}{{xK}}_{1/3}(x)\right],\end{eqnarray} \tag{ 2 }$$

where $\tilde{\chi }=2\xi /(3{\chi }_{e}(1-\xi ))$, $\xi ={\chi }_{\gamma }/{\chi }_{e}$ and $\alpha ={e}^{2}/({\hslash }c)$ is the fine-structure constant.

In order to simulate the emission of photons (and electron–positron pairs), we have added a QED module [44, 45] to OSIRIS [46] which allows real photon emission from an electron or a positron and decay of the photons into pairs (Breit–Wheeler process). The OSIRIS-QED framework accounts for the differential emission probability rate (2) in a similar fashion as other QED-PIC modules [9–11, 47]. The QED algorithm can be summarised as follows: at particle push-time, the probability of radiating a photon is evaluated, and if the event occurs, the radiated photon quantum parameter is selected to obey the distribution given by equation (2); the particle momentum is then updated to account for the momentum of the emitted photon (assumed to be radiated in the direction of the particle motion). For Breit–Wheeler pair production, the procedure is similar but instead of emission, we evaluate the probabilities of photon decay into an electron–positron pair. If the event occurs, we then remove the photon and initialise the new particles. The OSIRIS-QED framework is also equipped with an advanced macroparticle merging algorithm [48].

We first examine the temporal evolution of the energy spectrum of an electron beam as it collides head-on with an intense laser. In particular, we will focus on how the electron beam energy spread is affected by the QED photon emission. To facilitate this analysis, we define the beam width at a time t as the standard deviation in energy over all the particles

$$\begin{eqnarray}&&\sigma (t)=\sqrt{\tfrac{1}{N}\displaystyle \sum _{i=1}^{N}{(\overline{\gamma (t)}-{\gamma }_{i}(t))}^{2}},\end{eqnarray} \tag{ 3 }$$

where N is the total number of particles, $\overline{\gamma (t)}$ is the mean energy of the electron beam at time t, as measured in the laboratory frame, and ${\gamma }_{i}(t)$ is the energy of the particle i at the same time t. The stochastic effects in quantum radiation reaction that are responsible for the spreading of the distribution have been analytically studied in a similar setup in [36] where the Fokker–Planck equation is used to describe the evolution of the electron distribution function in time. This approach is valid for ${\chi }_{\gamma }\ll 1$. The Fokker–Planck equation [49–51] is usually used in laser plasma interaction, for instance, to model kinetically the collisions between species. One can however see the quantum photon emission as a virtual inelastic collision with an electron; as long as the momentum exchange remains small compared with the emitting particle momentum, the Fokker–Planck equation proves to be adequate.

If $w(\vec{p},\vec{q}){{\rm{d}}}^{3}\vec{q}$ denotes the probability per unit time of momentum change $\vec{p}\to \vec{p}-\vec{q}$ of an electron $\vec{p}$, then the transport equation for the electron distribution function $f(t,\vec{p})$ is given by:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f(t,\vec{p})}{\partial t}=\displaystyle \frac{\partial }{\partial {p}_{l}}\left[{A}_{l}f+\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {p}_{k}}({B}_{{lk}}f)\right],\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{A}_{l}=\int {q}_{l}w(\vec{p},\vec{q}){{\rm{d}}}^{3}\vec{q},\quad {B}_{{lk}}=\int {q}_{l}{q}_{k}w(\vec{p},\vec{q}){{\rm{d}}}^{3}\vec{q}\end{eqnarray} \tag{ 5 }$$

represent the drift and diffusion coefficients respectively, and indexes l and k denote different spatial components. Under the assumption that the electron beam is relativistic, and that the photons are radiated in the direction of motion, the problem is reduced to one-dimension.

Since the emission probability is given by equation (2) as a function of ${\chi }_{\gamma }$, we proceed to a change of variables using ${\chi }_{\gamma }/{\chi }_{e}\approx {\hslash }k/\gamma {mc}$ which is a consequence of the collinearity of the electrons and the emitted photons. We then get

$$\begin{eqnarray}&&A=\displaystyle \frac{\gamma {mc}}{{\chi }_{e}}{\int }_{0}^{{\chi }_{e}}\displaystyle \frac{{{\rm{d}}}^{2}P}{{\rm{d}}t\,{\rm{d}}{\chi }_{\gamma }}{\chi }_{\gamma }{\rm{d}}{\chi }_{\gamma },\quad B=\displaystyle \frac{{(\gamma {mc})}^{2}}{{\chi }_{e}^{2}}{\int }_{0}^{{\chi }_{e}}\displaystyle \frac{{{\rm{d}}}^{2}P}{{\rm{d}}t\,{\rm{d}}{\chi }_{\gamma }}{\chi }_{\gamma }^{2}{\rm{d}}{\chi }_{\gamma }.\end{eqnarray} \tag{ 6 }$$

After integration, the drift coefficient and the diffusion coefficient become respectively

$$\begin{eqnarray}&&A\approx \displaystyle \frac{2}{3}\displaystyle \frac{\alpha {m}^{2}{c}^{3}}{{\hslash }}{\chi }_{e}^{2},\quad B\approx \displaystyle \frac{55}{24\sqrt{3}}\displaystyle \frac{\alpha {m}^{3}{c}^{4}}{{\hslash }}\,\gamma \,{\chi }_{e}^{3},\end{eqnarray} \tag{ 7 }$$

which were first calculated in [36]. The Fokker–Planck equation (4), which is a special case of the master equation in the continuous limit, is valid for $q\ll p$. In our scenario, this validity condition can be expressed as ${\chi }_{\gamma }\ll {\chi }_{e}$ which is only conceivable for ${\chi }_{e}\ll 1$.

If there is no diffusion (B = 0), the equation of the characteristic in equation (4) is ${\rm{d}}p/{\rm{d}}t\simeq {mcd}\gamma /{\rm{d}}t=-A$. For electrons counter-propagating with a linearly polarised wave ${\chi }_{e}$ is given by ${\chi }_{e}=\sqrt{2}\gamma {a}_{0}{\hslash }{\omega }_{0}/{{mc}}^{2}$, while in a circularly polarised wave ${\chi }_{e}=2\gamma {a}_{0}{\hslash }{\omega }_{0}/{{mc}}^{2}$. This allows us to retrieve the classical result where the photon emission results in the electron relativistic factor γ decrease. The rate of this decrease in a linearly polarised wave is given [20, 25] by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\gamma }{{\rm{d}}t}=-{\alpha }_{\mathrm{rr}}{\gamma }^{2},\quad {\alpha }_{\mathrm{rr}}=\displaystyle \frac{4{e}^{2}{\omega }_{0}^{2}}{3{{mc}}^{3}}{a}_{0}^{2},\end{eqnarray} \tag{ 8 }$$

where ${\alpha }_{\mathrm{rr}}$ is a constant with units of frequency. For a circularly polarised wave ${\alpha }_{\mathrm{rr}}$ needs to be multiplied by a factor of two. By integrating equation (8) with ${\gamma }_{0}$ for initial Lorentz factor, we obtain $\gamma ={\gamma }_{0}/(1+{\alpha }_{\mathrm{rr}}{\gamma }_{0}t)$ in a linearly polarised wave and $\gamma ={\gamma }_{0}/(1+2{\alpha }_{\mathrm{rr}}{\gamma }_{0}t)$ in a circularly polarised wave, in agreement with [30]. By neglecting diffusion and assuming an initial Gaussian distribution for the electrons with initial standard deviation ${\sigma }_{0}$ and initial mean energy $\bar{{\gamma }_{0}}$, the authors in [36] have shown that if ${\sigma }_{0}\ll \bar{{\gamma }_{0}}$, the distribution remains approximatively Gaussian with an effective standard deviation

$$\begin{eqnarray}&&\sigma (t)=\displaystyle \frac{{\sigma }_{0}}{{(1+2{\alpha }_{{rr}}\bar{{\gamma }_{0}}t)}^{2}},\end{eqnarray} \tag{ 9 }$$

which is expressed for a quasi-monoenergetic relativistic electron beam as $\delta {\gamma }_{0}/\delta \gamma ={({\gamma }_{0}/\gamma )}^{2}$ [25]. It is not straightforward to rigorously expand this result to account for the diffusion term contribution. However, if now we assume that the drift is negligible (i.e. the average energy remains constant over a period of time $\bar{\gamma }\simeq \bar{{\gamma }_{0}}$), we obtain the usual diffusion equation, where we have performed the change of variables $p\simeq {mc}\gamma$:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}=\displaystyle \frac{B(t,\bar{{\gamma }_{0}})}{2{m}^{2}{c}^{2}}\displaystyle \frac{{\partial }^{2}f}{\partial {\gamma }^{2}}.\end{eqnarray} \tag{ 10 }$$

In the case of a Gaussian distribution, the standard deviation evolves as

$$\begin{eqnarray}&&\sigma (t)={\sigma }_{0}{\left(1+\displaystyle \frac{1}{{\sigma }_{0}^{2}{m}^{2}{c}^{2}}{\int }_{0}^{t}B(t^{\prime} ){\rm{d}}t^{\prime} \right)}^{1/2}.\end{eqnarray} \tag{ 11 }$$

It is therefore clear from equations (9) and (11) that there is a competition between the drift-like term that tends to compress the distribution width whereas the diffusion-like term tends to increase it. For an infinitesimally short period of time dt, the change of the distribution width at a time t due to the drift is given by differentiating equation (9) yielding

$$\begin{eqnarray}&&{\rm{d}}{\sigma }_{1}=-\sigma (t)\,4{\alpha }_{\mathrm{rr}}\overline{\gamma (t)}{\rm{d}}t\end{eqnarray} \tag{ 12 }$$

and the change due to the diffusion is obtained in a similar manner by differentiating equation (11)

$$\begin{eqnarray}&&{\rm{d}}{\sigma }_{2}=\sigma (t)\displaystyle \frac{B(t)}{2\sigma {(t)}^{2}{m}^{2}{c}^{2}}{\rm{d}}t.\end{eqnarray} \tag{ 13 }$$

We can then compute the total change of the electron distribution width within an interval dt:

$$\begin{eqnarray}&&{\rm{d}}\sigma =\sigma (t)\left[\displaystyle \frac{B(t)}{2\sigma {(t)}^{2}{m}^{2}{c}^{2}}-4{\alpha }_{\mathrm{rr}}\overline{\gamma (t)}\right]{\rm{d}}t.\end{eqnarray} \tag{ 14 }$$

A direct integration of the equation (14) is not possible because the variables cannot be separated. The expression from [36] can be retrieved by approximating $\sigma (t)={\sigma }_{0}$ in the term within the squared brackets in equation (14) and then integrating in time. The authors in [36] have shown that their expression is valid for relatively short laser interaction times τ such that $\alpha {a}_{0}{(\bar{{\gamma }_{0}}/{\sigma }_{0})}^{2}{\chi }_{e}^{2}\tau {\omega }_{0}\ll 1$.

Considering that the width of the distribution can change significantly throughout the interaction, it turns out impossible to simplify equation (14) and obtain an explicit form for the energy spread evolution. Nevertheless, equation (14) can still provide an insight into the changing features of the electron distribution function. Depending on the specific values of the initial parameters of the beam ${\rm{d}}\sigma$ can be positive (the diffusion wins over the drift) or negative (the drift wins over the diffusion). If we start from a narrow momentum spread like in our simulations, the distribution width first tends to increase, and later shrink. The 'turning point', when classical-like drift starts winning over the QED-induced diffusion can be defined by solving ${\rm{d}}\sigma /{\rm{d}}t=0$. For a circularly polarised laser, we obtain

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\sigma }_{T}}{{\gamma }_{T}}\right)}^{2}\approx \displaystyle \frac{55}{32\sqrt{3}}\displaystyle \frac{{\hslash }{\omega }_{0}}{{{mc}}^{2}}{\gamma }_{T}{a}_{0}=\displaystyle \frac{2.4}{{\lambda }_{0}[\mu {\rm{m}}]}\times {10}^{-6}{\gamma }_{T}{a}_{0},\end{eqnarray} \tag{ 15 }$$

and for a linearly polarised laser

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\sigma }_{T}}{{\gamma }_{T}}\right)}^{2}\approx \displaystyle \frac{55}{32\sqrt{6}}\displaystyle \frac{{\hslash }{\omega }_{0}}{{{mc}}^{2}}{\gamma }_{T}{a}_{0}=\displaystyle \frac{1.7}{{\lambda }_{0}[\mu {\rm{m}}]}\times {10}^{-6}{\gamma }_{T}{a}_{0},\end{eqnarray} \tag{ 16 }$$

where ${\gamma }_{T}$ and ${\sigma }_{T}$ represent the electron average energy and energy spread at the 'turning point', ${\sigma }_{T}/{\gamma }_{T}$ is the relative energy spread and ${\lambda }_{0}$ is the laser wavelength. For a given ${\gamma }_{T}$ and a₀, if $\sigma \lt {\sigma }_{T}$ the energy spread increases, but if $\sigma \gt {\sigma }_{T}$, the width of the energy distribution function decreases. In other words, equations (15) and (16), allow us to estimate what is the maximum attainable energy spread through diffusion depending on the laser vector potential and the average energy of the electron beam. The evolution of the electron energy distribution function can be also interpreted in terms of entropy. Neitz and Di Piazza have shown that the entropy of the electron beam increases when the quantum stochasticity dominates, and decreases when the energy spread is decreasing [36].

We will now compare the above findings with simulation results obtained using the OSIRIS-QED framework. The simulation setup is depicted in figure 1 where an electron bunch is colliding head-on with a circularly polarised laser and emitting photons. The electrons are presented in their transverse momentum space as a function of the longitudinal spatial coordinate. The main characteristic of the quantum radiation emission is already visible in this figure: even though there is an average trend to emit and to lose more energy further into the laser pulse, the energy of an individual electron is subject to fluctuations due to stochastic nature of the quantum photon emission.

Zoom In Zoom Out Reset image size

To illustrate these features, we first present a set of simulations using two different electron bunches with mean energies of 0.5 and 0.85 GeV. The bunches are initialised with a very small thermal momentum spread, equal in all transverse directions (the initial beam divergence is ${p}_{\perp }/{p}_{\parallel }\,\sim \,0.2$ mrad). The laser is modelled as a transverse plane wave with a temporal envelope function f(t). The laser rise and fall sections have the same shape and duration ${\tau }_{\mathrm{rise}}={\tau }_{\mathrm{fall}}=50.0\ {\omega }_{0}^{-1}$, while the duration of the flat part ${\tau }_{\mathrm{flat}}$ is varied between 0.0 and 300.0 ${\omega }_{0}^{-1}$ with a step of 50.0 ${\omega }_{0}^{-1}$ (seven different total pulse durations $\tau ={\tau }_{\mathrm{flat}}+({\tau }_{\mathrm{rise}}+{\tau }_{\mathrm{fall}})/2$). The slope of the envelope function for $t\lt {\tau }_{\mathrm{rise}}$ is defined as $f(t)=10{(t/{\tau }_{\mathrm{rise}})}^{3}-15{(t/{\tau }_{\mathrm{rise}})}^{4}+6{(t/{\tau }_{\mathrm{rise}})}^{5}$, where ${\tau }_{\mathrm{rise}}=50.0\ {\omega }_{0}^{-1}=26.6\ \mathrm{fs}$ and ${\omega }_{0}=1.88\times {10}^{15}$. The variable length middle section of the pulse had the laser vector potential always at the maximum value (${a}_{0}=27$). We shall stress that the interaction between the particles of the beam is negligible and since the laser field amplitude is uniform in the transverse directions, all the particles are subject to the same conditions. This corresponds to the approximation of beam transverse size much smaller than the laser spot size. Therefore, the large amount of PIC particles (about 1 million) provides us with a good statistical sample to study the evolution of the energy spectra of the electrons. The simulations are performed in 2D, where the box size was $500\times 20\ {c}^{2}/{\omega }_{0}^{2}$, resolved with 5000 × 200 cells and the timestep ${\rm{d}}t=0.04\ {\omega }_{0}^{-1}$ using 16 particles per cell. The simulation timestep is chosen such that ${\rm{d}}t\ll {W}_{\mathrm{rad}}^{-1}$, where W_rad is the photon emission rate.

The results of the first set of simulations are summarised in figure 2 which shows the electron energy spectra after the interaction with the laser. All the spectra, as expected, are wider than the quasi-monoenergetic initial distribution. The striking fact is that after interaction with longer lasers, the final energy spread of the electron beam is narrower than after interacting with shorter lasers. This hints that in case of longer interaction, a 'turning point' predicted by the theory with properties given by equation (15) must exist. After the initial increase of the width of the spectrum due to the quantum nature of the radiation process, the 'turning point' indicates the time at which the width reduces anew as predicted by classical radiation reaction [25].

Zoom In Zoom Out Reset image size

To investigate this further, the evolution of the beam energy with the spread $\overline{\gamma (t)}\pm \sigma (t)/2$ as a function of time for several examples is shown in figure 3. Without quantum stochasticity or radiation reaction, both the electron energy and $\sigma (t)$ would remain the same throughout the whole interaction. With only classical radiation reaction, both the average energy and $\sigma (t)$ would reduce with time. The analysis of the beam spectral width evolution during the time of the interaction confirms that the first quantum effect is to broaden the spectrum due to quantum stochasticity. If the laser is short enough, the spectrum stays broad (this is in agreement with [36]). However, if the laser is longer, then there is a specific point in time where the spread starts decreasing due to the classical drift of the electron energy distribution function.

Zoom In Zoom Out Reset image size

A second set of simulations is performed by varying the electron beam initial energy, using a laser pulse similar to the ones described previously (${a}_{0}=27$, ${\tau }_{\mathrm{rise}}={\tau }_{\mathrm{fall}}=50\,{\omega }_{0}^{-1}$, ${\tau }_{\mathrm{flat}}=600\,{\omega }_{0}^{-1}$), with the same simulation box and resolution. We would like to compare the predictions of equation (14) with the simulation results in a regime with ${\chi }_{e}\ll 1$, in a wave with a constant amplitude that allows for direct integration of equation (14). The black line in figure 4(a) shows the evolution of the energy spread from a simulation with $\gamma =100$, where ${\chi }_{e}\simeq 0.01$. We are interested only in the constant amplitude section, so we select $t=27\,{\omega }_{0}^{-1}$ as the new 'initial time'(${\tau }_{\mathrm{rise}}/2+{\tau }_{\mathrm{beam}}$, where ${\tau }_{\mathrm{beam}}=2\ {\omega }_{0}^{-1}$ is the electron beam 'duration'). Therefore, ${\sigma }_{0}=\sigma (t=27\,{\omega }_{0}^{-1})$ and ${\gamma }_{0}=\gamma (t=27\,{\omega }_{0}^{-1})$. The numerical integration of equation (14) is computed from this new 'initial time' in order to compare with the simulation results. The panel (a) of figure 4 shows that equation (14) gives a result in good agreement with the simulation. Even in the case of a higher ${\gamma }_{0}=1700$, which corresponds to ${\chi }_{e}\simeq 0.2$, the integration of equation (14) provides a reasonable agreement (figure 4(b)).

Zoom In Zoom Out Reset image size

In figure 4(c) we show ${\sigma }_{T}^{2}$ as a function of ${\gamma }_{T}^{3}$ at the 'turning point' for several simulations starting at different average electron energies. All the simulations are performed with ${a}_{0}=27$ and the 'turning point' is located within the constant amplitude section of the pulse. For particles with lower energies (and therefore with lower ${\chi }_{e}$), the 'turning point' is well identified by equation (15). For higher energies, the ${\chi }_{e}$ parameter is close to one and the electron energy spread is high, which makes the simulation results depart from the prediction of equation (15). However, the value obtained in the simulations is always lower than the predicted value. This confirms that the upper limit on the electron energy spread increase through diffusion as a function of a₀ and γ can be estimated using equations (15) and (16). The predictions of equation (15) are shown in figure 4(d) for a range of different values of the laser intensities and electron energies.

Let us comment on the underlying physics involved in this behaviour. The differential probability rate of photon emission given by equation (2) depends on the ${\chi }_{e}$ parameter in such manner that electrons with higher ${\chi }_{e}$ emit on average a larger fraction of their energy than the electrons with low ${\chi }_{e}$. This is what leads to the classical-like shrinking of the electron beam energy distribution in addition to the average energy drift towards a lower value. Moreover, the photons in the nonlinear Compton regime are emitted according to a distribution, such that electrons in identical conditions can radiate photons of different energy; this leads to a diffusion in the electron distribution function. These two tendencies compete, and the drift effect becomes dominant if the energy spread is wide enough ($\sigma \gtrsim {\sigma }_{T}$). On the contrary, if the initial electron energy spread is very low, the diffusion process dominates. This is illustrated in figure 5 that shows the temporal evolution of electron energy spectra in a Gaussian laser pulse with duration $\tau =150\,\mathrm{fs}$ and peak vector potential ${a}_{0}=27$ (a pulse like this will be available, for instance, in ELI Beamlines [1]). An electron beam is initialised with ${\gamma }_{0}=1700$, and we varied the value of ${\sigma }_{0}$. The first value ${\sigma }_{0}=2$ (figure 5(a)) corresponds to an electron beam energy spread of 0.1% that is still a challenge for experiments. The second ${\sigma }_{0}=100$ in figure 5(b) corresponds to a 6% energy spread, that has already been achieved with GeV-class electron beams in state-of-the-art laser facilities [24]. Third value ${\sigma }_{0}=200$ in figure 5(c) corresponds to 12% energy spread, and this is routinely achievable with laser wakefield acceleration in the present day laboratory conditions.

Zoom In Zoom Out Reset image size

It is clear from figure 5 that starting from a very narrow distribution, the diffusion appears to be faster than for an initial wider one. All the three examples finally converge to the same electron energy distribution function. We can take advantage of this to calculate the expected electron energy spread at the end of the interaction. We perform the numerical integration of equation (14) assuming that the initial energy spread ${\sigma }_{0}$ is equal to the maximum attainable energy spread through diffusion ${\sigma }_{T}$ determined by equation (15) for a given a₀ and $\bar{{\gamma }_{0}}$. The result obtained in this way is valid for all ${\sigma }_{0}\lt {\sigma }_{T}$, as long as the total time of interaction is much longer than the typical emission time ${T}_{\mathrm{total}}\gg {W}_{\mathrm{rad}}^{-1}$. The expected final energy spread obtained through numerical integration of equation (14) is compared with the simulation results in figures 6(a) and (b) for all the different laser pulse durations considered in the first set of simulations. As expected, the agreement is better for longer laser pulses and lower average electron energies, but the order of the expected energy spread is well predicted in all cases (maximum error is about 30%).

Zoom In Zoom Out Reset image size

For guiding future experiments, it would be beneficial to have an explicit expression where one could insert the initial electron and laser parameters, and estimate the final energy spread of the electron beam that could be measured directly on a spectrometer. It is possible to attain such an approximate expression for laser pulses with a symmetrical, Gaussian-like temporal profile. We assume that the balance point ${\rm{d}}\sigma =0$ has been reached before the centre of the pulse. This assumption is reasonable for beams with ${\sigma }_{0}\sim {\sigma }_{T}$, and for beams with ${\sigma }_{0}\ll {\sigma }_{T}$ provided that there is enough time for the diffusion to act and increase ${\sigma }_{0}$ to the same order as ${\sigma }_{T}$. This is verified when the pulse duration satisfies ${\tau }_{\mathrm{fwhm}}\gt 1/(8{\gamma }_{0}{\alpha }_{\mathrm{rr}})$. If the pulse is shorter, the stochastic effects dominate over the drift, and the final energy spread ${\sigma }_{{\rm{F}}}$ can be directly estimated using equation (11). Similarly, for ${\sigma }_{0}\gg {\sigma }_{T}$, the drift dominates and we can approximate ${\sigma }_{{\rm{F}}}$ through equation (9). The most difficult case is, therefore, when ${\sigma }_{0}\sim {\sigma }_{T}$, as all the terms in equation (14) are of the same order, which renders the equation unintegrable. However, we can estimate an upper boundary for ${\sigma }_{{\rm{F}}}$ by assuming that at the central point of the laser (at the point of peak intensity) the electron beam is close to the balance between the drift and diffusion, i.e. ${\sigma }_{M}^{2}\approx (2.4/\lambda [\mu {\rm{m}}])\times {10}^{-6}{\gamma }_{M}^{3}{a}_{0}$, where ${\gamma }_{M}$ is the average Lorentz factor of the electron beam in the central laser point. As ${\gamma }_{M}$ is easy to calculate (see [20, 25, 52–54]), we can retrieve an explicit expression for ${\sigma }_{M}$ as a function of laser intensity and duration, and initial electron energy. Beyond this point, the energy spread slowly decreases, and the final electron energy spread ${\sigma }_{{\rm{F}}}$ is smaller than ${\sigma }_{M}$. This yields

$$\begin{eqnarray}&&{\sigma }_{{\rm{F}}}^{2}\lesssim 1.455\times {10}^{-4}\sqrt{{I}_{22}}\displaystyle \frac{{\gamma }_{0}^{3}}{{(1+6.12\times {10}^{-5}{\gamma }_{0}{I}_{22}{\tau }_{0}[\mathrm{fs}])}^{3}},\end{eqnarray} \tag{ 17 }$$

where ${I}_{22}=I[{10}^{22}\,{\rm{W}}\,{\mathrm{cm}}^{-2}]$ and ${a}_{0}=0.855\sqrt{I[{10}^{18}\,{\rm{W}}\,{\mathrm{cm}}^{-2}]}\lambda [\mu {\rm{m}}]$ for linear polarisation and ${a}_{0}=0.855\sqrt{I[{10}^{18}\,{\rm{W}}\,{\mathrm{cm}}^{-2}]}\lambda [\mu {\rm{m}}]/\sqrt{2}$ for circular polarisation. It is worth noting that the result presented in equation (17) does not depend on the laser polarisation, but solely on intensity and duration.

Figures 6(a) and (b) shows the estimate given by equation (17) compared with the simulation results. Even though the lasers in our simulations are not Gaussian, we obtain a satisfactory agreement for the same ${\tau }_{\mathrm{fwhm}}$. Panels (c) and (d) show the predictions for the final energy spread according to equation (17) for electron beams starting at different initial energies after interacting with a 30 and a 100 fs laser of $2\times {10}^{21}\,{\rm{W}}\,{\mathrm{cm}}^{-2}$ intensity. These laser durations are to be available in the near-future laser facilities such as ELI [1], so there is a possibility to verify this model in the next few years.

In addition to the electron energy spread, we can also evaluate the impact of the laser interaction on the electron beam divergence. We define the weighted average of the deflection angle from the main propagation direction as

$$\begin{eqnarray}&&\overline{\tan \theta }=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{q}_{i}{\left(\displaystyle \frac{{p}_{\perp }}{{p}_{\parallel }}\right)}_{i}}{{\displaystyle \sum }_{i=1}^{N}{q}_{i}},\end{eqnarray} \tag{ 18 }$$

where N is the total number of simulation particles, q_i is the charge weight of the ith particle, and ${({p}_{\perp }/{p}_{\parallel })}_{i}$ is the ratio of the transverse to the longitudinal momentum with respect to the direction of laser propagation. For small angles, $\tan \theta \simeq \theta$, and the average divergence shown in figure 7 is determined with this approximation (the error is less than 1 mrad).

Zoom In Zoom Out Reset image size

Figure 7 shows the evolution of the electron bunch divergence as time progresses. Classically, the radiation reaction leads to momentum phasespace contraction proportionally in transverse and longitudinal direction. According to the analytical solution for trajectory of a relativistic electron in an intense plane wave [30], on average, all momentum components and electron energy are reduced by a same factor due to radiation reaction. The angle between the particle momentum and the laser propagation direction is therefore approximately the same before and after the interaction, provided that the laser has a slowly varying temporal envelope compared to the laser period. However, during the interaction with the laser, the electron has an additional oscillatory component of the transverse momentum, whose amplitude depends only on the laser intensity (${p}_{\perp }^{* }\,\simeq \,{a}_{0}{mc}$). It is worth noting that the oscillatory component of the ${p}_{\perp }$ given by the laser is the same with and without radiation reaction, as it depends only on the normalised vector potential of the wave a₀ (see supplementary material [55]).

In our case, the particles are counter-propagating with the laser, so the initial ${p}_{\perp 0}=0$. The only transverse momentum of such a particle is therefore ${p}_{\perp }^{* }$. In the field of a circularly polarised wave, the transverse momentum vector rotates in the plane perpendicular to the laser propagation direction, while its magnitude remains constant ${p}_{\perp }\,=\,{p}_{\perp }^{* }\,\simeq \,{a}_{0}{mc}$.

As the initial electron energy is on the order of a 0.5–1.0 GeV, we take ${\gamma }_{0}\gg {a}_{0}$ and the parallel momentum becomes ${p}_{\parallel }\,=\,{mc}\sqrt{{\gamma }^{2}-{a}_{0}^{2}-1}\simeq \gamma {mc}$. As a result, the average angle that a single electron makes with the direction of laser propagation is approximately $\theta \simeq {a}_{0}/\gamma$. Without radiation reaction, this angle would stay the same during the constant amplitude section of the laser, as there would be no change in γ. However, equation (8) shows that with radiation reaction the electron Lorentz factor decreases through $\gamma ={\gamma }_{0}/(1+2{\alpha }_{\mathrm{rr}}{\gamma }_{0}t)$ and the angle is then expected to increase linearly as a function of time:

$$\begin{eqnarray}&&\theta \simeq \displaystyle \frac{{a}_{0}}{{\gamma }_{0}}(1+2{\alpha }_{\mathrm{rr}}{\gamma }_{0}t).\end{eqnarray} \tag{ 19 }$$

The dashed lines in figure 7 show the average expected angle ${a}_{0}/\gamma$ during the constant amplitude part of the pulse, where γ is taken as the average relativistic factor of the electron bunch and ${a}_{0}=27$. We observe a similar trend with the simulation data, which indicates that the average divergence increase due to radiation emission in the constant amplitude region of the laser envelope is well explained by the semi-classical approach. However, there is a slight difference between the simulation data and the expected ${a}_{0}/\gamma$ which increases over time.

After the interaction has finished, the electron beam has a residual divergence on the order of ${\theta }_{{\rm{F}}}\sim 10\,\mathrm{mrad}$ which is larger than the initial divergence on the 0.2 mrad level. Semi-classical radiation reaction predicts the final divergence to be approximately equal to its initial value. In QED, we expect the final $\bar{{p}_{\perp }}=0$, but the divergence defined by equation (18) can have a larger value than initial if there is a wider particle distribution function in transverse momentum space.

To ascertain the origin of this effect, we have examined the transverse momentum space at different times (see figure 8). Initially the electron beam has a narrow momentum spread. During the plane wave stage the average transverse momentum is indeed around the predicted value ${p}_{\perp }\,\simeq \,{a}_{0}{mc}$. Howbeit, the QED simulations show the existence of a momentum spread around the average value that increases with time. This is consistent with the transverse momentum spreading previously reported in [56]. The variation around the average angle as defined in equation (18) during the interaction with the constant-amplitude section of the laser can be approximately related to the variation in energy:

$$\begin{eqnarray}&&{\rm{\Delta }}\theta \simeq \displaystyle \frac{{a}_{0}}{{\gamma }^{2}}{\rm{\Delta }}\gamma .\end{eqnarray} \tag{ 20 }$$

This variation persists and finally becomes the net beam divergence when the laser shuts down: ${\theta }_{{\rm{F}}}\simeq \sqrt{2/\pi }({a}_{0}/{\gamma }_{{\rm{F}}}^{2}){\sigma }_{{\rm{F}}}$. As ${\gamma }_{{\rm{F}}}$ converges to a lower value for a longer interaction time, and σ is from equation (15) approximately proportional to ${\gamma }^{3/2}$, we can then conclude that the width of the final angular spread increases slowly with the length of the interaction (as seen in figure 7).

Zoom In Zoom Out Reset image size

In classical radiation reaction, the energy loss of a single electron depends on its initial energy. For an electron beam, the main effects are the decrease in its mean energy and reduction of the energy distribution width. When QED effects are taken into account, the intrinsic stochastic nature of photon emission leads to diffusion in the energy distribution around the mean value which would translate into the increase of the energy spread of the beam. Therefore, in the general scenario, there is a competition between these two tendencies. If we allow a long enough interaction time, there is a point when the diffusion in momentum, intrinsically quantum, is balanced by the energy width reduction, associated with the classical regime. Beyond this point, the energy spread only decreases. This allows for estimating the maximal attainable energy spread through diffusion for a set of initial parameters ${\gamma }_{0}$ and ${\sigma }_{0}$. We have estimated this limit (and confirmed it with numerical simulations), which has further allowed us to predict analytically the final electron energy spread, which is the relevant quantity to be measured in experiments.

The average divergence of the electron beam during the laser interaction is well-described by the classical radiation reaction. However, we have observed that the electron distribution function in momentum space has a certain spread around the average value that increases with the interaction time. This spread persists after the interaction is shut down and leads to a residual divergence of the electron beam that can be estimated analytically through its connection with the electron energy distribution function.

The control of beam properties is of relevance for all near future laser facilities that will operate at high intensities, regardless if they are aimed at optimising particle acceleration, radiation sources or fundamental research. As the quantum spreading might discriminate between the measurable effects and those whose signatures are too small to be observed due to the width of the final distribution function, our findings are vital for the design of upcoming experiments. They are also valuable for numerous applications with specific beam quality requirements.

This work is supported by the European Research Council (ERC-2010-AdG Grant 267841) and FCT (Portugal) Grants SFRH/BD/62137/2009 and SFRH/IF/01780/2013. Simulations were performed at Supermuc (Germany) under a PRACE Grant and at the Accelerates cluster (Lisbon, Portugal).