Ultra-intense laser pulses in near-critical underdense plasmas – radiation reaction and energy partitioning

Abstract

Although, for current laser pulse energies, the weakly nonlinear regime of laser wakefield acceleration is known to be the optimal for reaching the highest possible electron energies, the capabilities of upcoming large laser systems will provide the possibility of running highly nonlinear regimes of laser pulse propagation in underdense or near-critical plasmas. Using an extended particle-in-cell (PIC) model that takes into account all the relevant physics, we show that such regimes can be implemented with external guiding for a relatively long distance of propagation and allow for the stable transformation of laser energy into other types of energy, including the kinetic energy of a large number of high energy electrons and their incoherent emission of photons. This is despite the fact that the high intensity of the laser pulse triggers a number of new mechanisms of energy depletion, which we investigate systematically.

1 Introduction

Laser wakefield acceleration (LWFA) is a means for utilising the extreme fields accessible in plasmas for the purpose of accelerating electrons to very high energies over short distances (Tajima & Dawson 1979). Originally starting out as a basic principle of electron acceleration, this field has matured into one with a manifold of applications (Malka et al. 2008; Corde et al. 2013). As laser technology has developed, so have the attainable laser intensities and repetition rates. A natural question is therefore if the progress in the field of LWFA is likely to continue as we ramp up the power of the laser facilities, or if new physical phenomena may affect the development of laser-driven electron accelerators? Here we try to shed light on this question by invoking large-scale particle-in-cell (PIC) simulations of LWFA. Indeed, we find that there are effects that will alter the acceleration process as we increase the laser intensity of the pulses used to generate the wakefield. Moreover, due to the nature of these processes, we find that there can also be certain benefits, such as an efficient energy transformation from the optical range to the range of extreme ultraviolet (XUV) and gamma rays (Corde et al. 2013). We therefore expect that, guided by our simulations and analytical understanding, new possibilities may open up for future LWFA systems.

Recent, and planned, developments in laser systems (ELI 2017; Vulcan 2017; XCELS 2017) are aimed at reaching new regimes for laser–matter interactions (Mourou, Tajima & Bulanov 2006). The laser radiation of the expected intensities will trigger new semi-classical (or even quantum electrodynamical (QED)) effects (Marklund & Shukla 2006; Heinzl & Ilderton 2009; Di Piazza et al. 2012), that should be taken into account in theoretical analyses and computational models. Furthermore, computational models are inherently limited by their finite resolution, posing limits to their validity. This can in principle lead to faulty estimates when considering certain physical effects. One such effect is that of radiation reaction (RR) (Burton & Noble 2014; Vranic et al. 2016) in PIC simulations, where the problem is due to the limitations in resolving the emitted high frequency radiation (Wallin, Gonoskov & Marklund 2015). The correct inclusion of such an effect can change the output from e.g. a laser wakefield accelerator.

In our simulations, we overcome this problem via the utilization of (i) a radiation reaction module as well as (ii) a synchrotron module. We furthermore include the process of pair production via a Breit–Wheeler process. Certainly, at high enough intensities of the laser radiation, this classical description is not applicable as the energies of electrons and emitted photons become comparable. One then needs a quantum description of the emission process, which implies probabilistic generation of a photon accompanied with a respective recoil. The major differences between these descriptions are the stochastic properties of photon emission and the discrete nature of the photons, which can for example lead to a broadening of the spatial distribution of photon emission. By comparing our numerical results with the ones obtained using the further extended QED–PIC model (Gonoskov et al. 2015) we identify the range of intensities ${<}10^{25}~\text{W}~\text{cm}^{-2}$ , where the outlined effects do not contribute substantially enough to affect our conclusions (see figure 1). In particular, we do not find any significant pair production in these systems for intensities up to ${\sim}10^{26}~\text{W}~\text{cm}^{-2}$ (this is to be expected, as there are very few head-on collisions between the laser field and the electrons when the electron distribution is close to thermal). Throughout the current study we consider this intensity range and therefore use a PIC code extended with the inclusion of models (i) and (ii), leaving the related QED effects out of the scope of the present paper.

Figure 1. Comparing the effects of classical and quantum radiation reaction in a ultra-high intensity regime. Here $I=4\times 10^{24}~\text{W}~\text{cm}^{-2}$ and $N=3\times 10^{20}~\text{cm}^{-3}$ . In (a) we show the case of QED radiation reaction, in (b) classical radiation reaction via the Landau–Lifshitz model and in (c) the case without radiation reaction. All these three cases have self-consistent ion motion included. In (d) we have included classical radiation reaction via the Landau–Lifshitz model but kept the ions (hydrogen) stationary and in case (e) the ions are fixed and the radiation reaction is turned off. As expected, the notion of fixed ions gives results very far from the self-consistent picture. We note a significant effect of radiation reaction on the ion motion by comparison of (a–c), while the difference in (a) and (b) is negligible. Thus, the small differences between (a) and (b) supports the notion that we in the LWFA case can use the classical radiation reaction model with great accuracy even at very high intensities, putting our use of the classical model for radiation reaction for these systems on a firm footing.

1.1 Laser wakefield acceleration

In laser wakefield acceleration (Tajima & Dawson 1979; Sprangle et al. 1988) a short, high intensity pulse propagates through an underdense plasma. Electrons are pushed to the sides of the laser pulse, with some of them ending up as a bunch of electrons in the wake behind the laser. The electron cavity induces a strong electric field (hundreds of $\text{GV}~\text{m}^{-1}$ (Modena et al. 1995; Leemans et al. 2002; Malka et al. 2002)), accelerating the electrons behind the laser. In the bunch, the electrons undergo transverse oscillations (Kiselev, Pukhov & Kostyukov 2004), emitting high-frequency radiation in the process. The output is a highly collimated, focused and quasi-monoenergetic electron beam (Geddes, Toth & Tilborg 2004; Leemans et al. 2006), as well as high-frequency X-ray radiation (Kiselev et al. 2004; Matsuoka et al. 2010; Corde et al. 2013). Apart from maximizing electron energy (Lu et al. 2006a , b , 2007), the studies of LWFA can be important for increasing the total number of accelerated electrons, as well as intensity of the outgoing X-ray radiation. The latter can certainly benefit from highly nonlinear regimes at ultra-high intensities of the laser radiation. Motivated by this, in the present paper we will be analysing the effects of increasing laser intensities on the properties of the electrons and the radiation spectra in LWFA.

In terms of potential applications, one of the main issues is the pulse depletion mechanism, which is of key importance for controlling the rate of energy transformation to other forms, as well as the proportion between them. The process of depletion is known to be difficult in terms of theoretical analysis and mathematical modelling even in the case of low intensities and densities (Bulanov et al. 1992). As a starting point, here we present a phenomenological analysis based on systematic comparison of simulations with and without taking into account the RR effect for different intensities of the driving laser pulse.

1.2 Wakefield acceleration with radiation reaction included

As we increase the intensity of the laser pulse generating the wake field, we expect that the electrons will be further accelerated. However, as the energy loss from synchrotron emission becomes appreciable, it will begin to affect the electron dynamics. This will then alter the characteristics of the wakefield, particle acceleration process, and radiation emission.

When the RR force is properly modelled in a PIC simulation (in the classical regime), the dynamics of the LWFA changes. One of the main principal changes is that when the laser hits the plasma, not all electrons are forced around the laser pulse. Instead, some electrons can enter the high intensity part and start to co-propagate with the laser pulse, being trapped there by the so-called radiation reaction trapping (RRT) effect (Ji et al. 2014). The number of such electrons is dependent on their energy, and thus the laser intensity. For the electrons trapped in the wakefield, the RR force reduces the acceleration of the electrons, and this also changes the spectrum of the emitted radiation. As the RR force acts like friction, we will also see pulse energy depletion at a faster rate than when not including this force. All-in-all, the combination of these, and other, effects makes for rather complex electron behaviour in an already nonlinear regime.

2 Governing equations

2.1 The particle-in-cell method

Particle-in-cell methods (Dawson 1983; Birdsall & Langdon 1985) are now a standard tool used for simulating laser–matter interactions. In this method the plasma is treated as an ensemble of particles moving in an electromagnetic (EM) field defined on a grid. The dynamics of the particles is calculated using the Lorentz force, and the resulting charge and current distributions in turn give rise to EM fields via Maxwell’s equations.

Even if the evolution of the plasma (and the laser) is solved self-consistently using the Lorentz force and Maxwell’s equations in this model, there is still some important physics not included (Gonoskov et al. 2015; Wallin et al. 2015). The grid and time resolutions place an upper limit on the frequency of the radiation accounted for in the simulations. This becomes a problem for relativistic particles, which can emit synchrotron radiation with a typical frequency of (Jackson 1998)

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{c}={\textstyle \frac{3}{2}}\unicode[STIX]{x1D714}_{H}\unicode[STIX]{x1D6FE}^{3},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{H}$ is the instantaneous cyclotron frequency of the particle and $\unicode[STIX]{x1D6FE}$ the relativistic gamma factor. In a typical PIC simulation the particle motion (thus $\unicode[STIX]{x1D714}_{H}$ ) is resolved, but for a high enough $\unicode[STIX]{x1D6FE}$ the emitted radiation, with typical frequency $\unicode[STIX]{x1D714}_{c}$ , cannot be resolved.

This radiation is typically emitted over a small angle ${\sim}1/\unicode[STIX]{x1D6FE}$ and can be modelled as photons (Gonoskov et al. 2015; Harvey, Ilderton & King 2015; Wallin et al. 2015) instead of fields. As such these will have no direct impact on the particle dynamics of the simulation, but for relativistic particles the energy of this emitted radiation can constitute a large part of the particle energy. To accurately model this we must include the back reaction of the emitted radiation on the emitting particle, which we do via the RR force.

2.2 The classical radiation reaction force

We include the effects of RR by adding a correction term to the equation of motion

(2.2) $$\begin{eqnarray}m\dot{\boldsymbol{v}}=\boldsymbol{F}_{ext}+\boldsymbol{F}_{rad},\end{eqnarray}$$

where the RR term is determined so that the work performed by the force is equal to the emitted energy. This produces the Abraham–Lorentz equation (Jackson 1998) which contains a third-order derivative, enabling unphysical runaway solutions. This can be avoided by approximating the third-order time derivative with the Lorentz force, yielding the Landau–Lifshitz (LL) equation, valid when the RR force is much less than the Lorentz force in the instantaneous rest frame of the particle. The relativistic, covariant expression for the RR force is then given by (Landau & Lifshitz 1975)

(2.3) $$\begin{eqnarray}\displaystyle f^{\unicode[STIX]{x1D707}} & = & \displaystyle \frac{2e}{3}r_{0}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D60D}^{\unicode[STIX]{x1D707}\unicode[STIX]{x1D708}}u_{\unicode[STIX]{x1D708}}u^{\unicode[STIX]{x1D6FE}}+\frac{2}{3}r_{0}^{2}[\unicode[STIX]{x1D60D}^{\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D60D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D708}}u^{\unicode[STIX]{x1D708}}+(\unicode[STIX]{x1D60D}_{\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FC}}u^{\unicode[STIX]{x1D6FC}})(\unicode[STIX]{x1D60D}^{\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FD}}u_{\unicode[STIX]{x1D6FD}})u^{\unicode[STIX]{x1D707}}],\end{eqnarray}$$

where $\unicode[STIX]{x1D60D}^{\unicode[STIX]{x1D707}\unicode[STIX]{x1D708}}$ is the electromagnetic field tensor, $u_{\unicode[STIX]{x1D708}}$ the 4-velocity and $r_{0}=\text{e}^{2}/m_{e}c^{2}$ is the classical electron radius in cgs units. We neglect the first term since it is linear in the field strength and its contribution is found to be negligible. (Indeed, it can be shown that, in cases where classical RR is important, the derivative term is even smaller than the electron spin force and so should be neglected out of consistency (Tamburini et al. 2010).) The three-dimensional form of the equation, without the first term, is (Landau & Lifshitz 1975)

(2.4) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{rad} & = & \displaystyle \frac{2}{3}r_{0}^{2}\left(\boldsymbol{E}\times \boldsymbol{B}+\frac{1}{c}\left[\boldsymbol{B}\times (\boldsymbol{B}\times \boldsymbol{v})+(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{E})\boldsymbol{E}\right]\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{\unicode[STIX]{x1D6FE}^{2}}{c}\left[\left(\boldsymbol{E}+\frac{1}{c}\boldsymbol{v}\times \boldsymbol{B}\right)^{2}-\left(\frac{\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{v}}{c}\right)^{2}\right]\boldsymbol{v}\right),\end{eqnarray}$$

where $\boldsymbol{E},\boldsymbol{B},\boldsymbol{v}$ and $c$ are electric field, magnetic field, velocity and speed of light respectively.

We note that although there are alternative formulations in the literature, all RR models result in approximately the same particle dynamics (Kravets, Noble & Jaroszynski 2013; Vranic et al. 2016), and the LL equation has been shown to be consistent with QED (Ilderton & Torgrimsson 2013). For an overview of other models see Burton & Noble (2014).

Figure 2. Comparison of the results for simulations with and without the RR force at $0.53~\text{ps}$ , when the laser has propagated $160~\unicode[STIX]{x03BC}\text{m}$ into the plasma. (a) The electron density (green to black) and the $y$ -component of the electric field (red to blue) for the intensity equal to $8\times 10^{21}~\text{W}~\text{cm}^{-2}$ , for which the two cases are inseparable. (b) The angular plot of the synchrotron emission for the same case. (c–h) The same information (note: with the same scales and axes) for the cases with (a,c,e,g) and without (b,d,f,h) the RR force shown in a more compact form for the intensities $8\times 10^{22}~\text{W}~\text{cm}^{-2}$ , $8\times 10^{23}~\text{W}~\text{cm}^{-2}$ and $4\times 10^{24}~\text{W}~\text{cm}^{-2}$ counting from left to right. (e) The two boxes show the regions in which the number of electrons in the laser pulse and electron bunch regions are measured.

3 Simulations

We use the three-dimensional PIC code ELMIS (Gonoskov 2013) with enabled modules for the classical RR in the form of the LL equation (2.4) and which accounts for the classical synchrotron emission via statistical routines (Wallin et al. 2015). The modules are incorporated via the module development kit interface (Gonoskov et al. 2015). This shared interface (with the PIC code PICADOR (Bastrakov et al. 2012)) enables extensions of the PIC scheme. The simulations are performed with a fixed plasma configuration for a number of different laser intensities, comparing the cases with and without the RR force. The laser pulse is a linearly polarized Gaussian shaped beam with duration $20~\text{fs}$ and diameter $8~\unicode[STIX]{x03BC}\text{m}$ (full width at half maximum (FWHM) in intensity). The wavelength is $1~\unicode[STIX]{x03BC}\text{m}$ and the pulse energy varies from $10^{2}~\text{J}$ to $5\times 10^{4}~\text{J}$ for different simulations, giving an intensity range from $8\times 10^{21}~\text{W}~\text{cm}^{-2}$ to $4\times 10^{24}~\text{W}~\text{cm}^{-2}$ (total energy over FWHM values of duration and diameter squared).

To enable highly nonlinear stable regimes for the whole range of intensity we consider for the target a near-critical density hydrogen plasma with a radially (perpendicular to the pulse propagation direction) parabolic density profile. Electrons are initially inserted into the wakefield via a longitudinal (along the pulse propagation direction) density profile. The density of the plasma along the laser propagation axis is $N_{axis}=3\times 10^{20}~\text{cm}^{-3}$ , or $N_{axis}\approx 0.27N_{cr}$ where $N_{cr}$ is the critical density, with the density increasing radially to $2N_{axis}$ at $40~\unicode[STIX]{x03BC}\text{m}$ (the edge of the simulation box). The longitudinal density profile of the plasma is such that the density increases from $0$ to $3N$ over the first $10~\unicode[STIX]{x03BC}\text{m}$ , then it decreases to $N$ over another $10~\unicode[STIX]{x03BC}\text{m}$ . This density is then kept for the remainder of the plasma. The size of the simulation box is $80~\unicode[STIX]{x03BC}\text{m}\times 80~\unicode[STIX]{x03BC}\text{m}\times 80~\unicode[STIX]{x03BC}\text{m}$ , on $512\times 128\times 128$ cells with $2$ macroparticles per cell. The simulation box is comoving with the laser pulse, though the simulation is still performed in the laboratory frame. In figure 2 pictures of the simulations after $0.53~\text{ps}$ can be seen, with the intensity ranging from $8\times 10^{21}~\text{W}~\text{cm}^{-2}$ to $4\times 10^{24}~\text{W}~\text{cm}^{-2}$ .

The high density (near-critical) plasma results in a high level of interaction between the particles and the pulse, depleting the pulse after a moderate amount of time. For the laser intensity of $4\times 10^{23}~\text{W}~\text{cm}^{-2}$ the pulse loses half of its energy in $1.7~\text{ps}$ , having travelled $L\approx 500\unicode[STIX]{x1D706}$ . This enables us to study the effect of the pulse depletion, for not too long simulation times.

Figure 3. Ratio of peak electron energies (left axis) and pulse depletion times (right axis) for simulations without to those with RR included. The pulse depletion time is calculated as the time taken for the laser to lose half its initial energy. The ratio of the peak electron energies in the two cases is calculated throughout the simulations and is generally found to increase with time. Here we plot the maximum value which typically occurs when the pulse is depleted. Note the different scales on the vertical axes.

4 Results

We compare, for a number of different intensities, how the inclusion of the RR force affects the time it takes for the laser pulse to be depleted. The result can be seen in figure 3, where we plot the ratio of the pulse depletion times (i.e. the time for the pulse to lose half its energy) without and with the RR force included ( $T_{NoRR}/T_{RR}$ ). It can be seen that when the laser intensity surpasses $10^{23}~\text{W}~\text{cm}^{-2}$ this ratio begins to differ from $1$ , indicating that RR effects are starting to reduce the propagation distance of the laser. If the pulse loses too much of its energy this effect could prevent the application of wakefield acceleration schemes where long plasmas are used. The inclusion of the RR force is equivalent to restoring the possibility of the electrons emitting high energy photons. Thus the faster pulse depletion can be understood as being a result of the extra work the field does in reaccelerating these electrons. This decrease in energy of the electrons can also be seen in figure 3, where the ratio of the maximum electron energies without and with RR included is shown. For intensities above $10^{22}~\text{W}~\text{cm}^{-2}$ this starts to be noticeable and for even higher intensities there is a great difference in the maximum electron energy due to the inclusion of RR. An example of the electron energy distribution, after 1.06 ps can be seen in figure 4. In figure 5 an example of the transformation of the initial pulse energy to different types of particle energy can be seen, for the two cases with and without RR for the intensity $4\times 10^{24}~\text{W}~\text{cm}^{-2}$ . It is shown that the pulse is drained faster in the case where RR is included, and one can also see how the total electron energy is overestimated when it is not. Also, one can see that a great amount of energy actually goes to the ions. This is regardless of whether RR is included or not.

Figure 4. Electron energy distribution at 1.06 ps. The inclusion of RR reduces the peak energy of the electrons.

Figure 5. Redistribution of laser pulse energy to electrons, ions (hydrogen) and photons over time for the intensity $4\times 10^{24}~\text{W}~\text{cm}^{-2}$ . The energy of the electrons is clearly overestimated without RR. Also, a substantial amount of energy goes to the ions in both cases.

Our regime of simulation is similar to LWFA, but with some important differences. The field intensity is sufficiently strong to perturb the hydrogen ions, and their motion changes the dynamics of the laser–plasma interaction (Capdessus, d’Humières & Tikhonchuk 2013). One effect of this is that the particle densities do not close behind the bubble as in conventional LWFA; instead there is a column after the bubble region containing substantially fewer particles, both electrons and ions. There seem to be two effects on the ions as they are hit by the laser, (i) some ions are forced around the pulse by the radiation pressure and (ii) some ions pass through the laser pulse region and interact strongly with the electron bunch. The result is a column behind the bubble region void of ions, except a filament in the centre. The void region does not attract either particle species and the particle densities do not close behind the bubble. In spite of this, the bubble is still formed and electrons stream in and become trapped behind the laser just as in conventional LWFA, as seen in figure 2. The ion density in the laser pulse region is relatively unchanged and there is thus still a strong electric field, induced by the charge separation, which can accelerate the electrons.

Figure 6. Electron density in the laser pulse region, measured in a box as shown in figure 2(e). The range shown is until the laser pulse starts to break down and the wakefield collapses allowing particles flow freely through the box. For the lowest intensity, the electrons do not flow as much around the pulse, but pass through it and there is no difference between the RR and no RR cases. For higher intensities the electrons are forced around the pulse, but the effect of RR counteracts this. More electrons flow through the pulse due to RR, and there is also a trapping effect, where for the higher intensity cases the electrons spend more time in the high intensity region. The thickness of the box is $2~\unicode[STIX]{x03BC}\text{m}$ , the same as the electron density slice in figure 2

Figure 7. Electron density in the electron bunch region, measured in a box around the bunch as shown in figure 2(e). At higher intensities the inclusion of RR increases the amount of charge in the bunch up to the point were the laser is too depleted to drive the bunch. For the case without RR the bunch electrons, as well as the laser pulse, have lost less energy and the bunch can be sustained for a longer time. For the lower intensities the pulse is quickly depleted, after which only the background plasma is measured.

The inclusion of the RR force also changes the interaction between the laser pulse and the plasma. One such effect is that electrons can become trapped and pass through the high intensity pulse region of the laser as it enters the plasma. This can be seen in figure 6 where the electron density in the high intensity part of the pulse is plotted as a function of time. For the lower intensity there is no visible effect of the RR. Some electrons pass through the pulse region as this is not sufficiently strong to force them around. As the laser gets more intense the electrons are forced around the pulse, and the electrons in the high intensity region decrease. The effect of including RR is to change this, allowing electrons to go through the high intensity pulse region, as well as to spending a longer time there. In figure 2 the electron density and the laser pulse is shown for the different intensity cases. The trapped particles passing through the intense pulse region can be seen for the intense cases when RR is included (figure 2 g). One can also note the difference in the size of the electron bunch, with a larger bunch for the cases when RR is included. This is further shown in figure 7 where electron density in the bubble region is plotted as a function of time.

The photon emission calculated from the electron motion is vastly affected by the inclusion of the RR force for these intensities, as expected. In figure 8 the frequency spectra of the emitted photon radiation is shown. The inclusion of RR prevents the very highest electron energies, and their high-frequency contribution is removed from the spectra. This is also visible in figures 2 and 9 where the angular spectra of the emitted radiation is shown. The high-frequency radiation removed by the inclusion of RR was emitted very close to the propagation direction of the pulse, which can be seen in the broader spectra for the case with RR included.

Figure 8. Radiation spectra for the four different intensity cases. For the lowest intensity there is little difference between the RR and no RR cases, but this grows very notable for the higher intensities where the inclusion of RR removes the high-frequency contribution, as the peak energy of the electrons is lowered.

Figure 9. Angular distribution of radiated energy as a function of the angle from the propagation direction of the laser pulse. Note that the energy scale is normalized.

5 Conclusions

Including a radiation reaction force in the equation of motion for electrons in a PIC simulation of a laser induced plasma wakefield can alter the dynamics of the electrons. Our simulations show that it can affect the pulse depletion distance for laser intensities from $10^{23}~\text{W}~\text{cm}^{-2}$ and the energy of the electrons in the electron bunch from intensities of $10^{22}~\text{W}~\text{cm}^{-2}$ , as well as affecting how the electrons interact with the laser.

The RR force models the radiation from the electrons which is of too high frequency/short wavelength (thus of large energy) to be resolved by the time step and the grid. The inclusion of this restores an important channel of energy loss for the electrons, with the result that the electrons in the bunch become less energetic than they would otherwise be. There is then also an effect where the electrons take more energy from the laser pulse, shortening the pulse depletion distance. The energy lost due to the RR force is emitted as photons. In the considered cases these have little impact on the particle dynamics, but allow for a wakefield set-up to be used as a gamma radiation source.

Our simulations reveal a highly nonlinear regime of laser pulse propagation in underdense or near-critical plasmas, similar to LWFA but with sufficiently high intensity for ion motion to be of importance. This affect the dynamics of the laser–plasma interaction, but still allows for electron acceleration.

For intensities approaching $10^{24}~\text{W}~\text{cm}^{-2}$ and beyond, it is necessary to approach the photon emission and RR using QED methods (Di Piazza et al. 2012). We have implemented such a stochastic QED-based routine (Gonoskov et al. 2015) for emission of photons and also pair production, and used for comparison for intensities up to $4\times 10^{24}~\text{W}~\text{cm}^{-2}$ . The results are similar to the case of the classical RR force as presented here and the main results are unchanged, even if the QED runs indicate that the energy loss of the electrons is slightly overestimated by the classical RR force.

For even higher intensities than those presented here, the stochastic nature of quantum photon emission, and the resulting electron recoil, is likely to change the dynamics of the wakefield significantly. In this case, due to the large number of hard photons emitted, Breit–Wheeler pair production will also have to be included (Gonoskov et al. 2015). However, we leave this issue for future studies.