Identifying the source of super-high energetic electrons in the presence of pre-plasma in laser–matter interaction at relativistic intensities

The influence of pre-plasma in laser–matter interaction at relativistic intensities has attracted much attention in both experimental and theoretical investigations, because of its significant effects on a number of applications, such as laser driven ion acceleration [1–7], fast ignition [8–11] and bright x/γ ray sources [12, 13], etc. The pre-plasma produced by the intrinsic laser pre-pulse (usually with ns duration) can be as high as $10~\mu \text{m}$ for the energetic main pulses of energies tens of $\text{kJ}$ with a typical contrast ratio 10⁻⁵. In fast ignition related experiments with relatively long pulses (with tens of ps duration), high intensity and high power laser, even though the contrast ratio can be as high as 10⁻⁸, considerable pre-plasma can still build up in front of a solid-density target. This pre-formed plasma always exists in laser–matter interaction at relativistic intensities, thus the laser–pre-plasma interaction is inevitable.

The issue of the hot electron generation due to relativistic intensity laser–solid interaction is typically characterized by Beg's scaling [14], which is based on experiments, or Wilks' numerically modelled ponderomotive scaling [15] and other following works [16–18]. However these scalings do not address the dependence of hot electrons on finite-scale length preformed plasma in front of the target. The first systematic study of the sources of super energetic electrons with energy exceeding ponderomotive energy is the stochastic heating and acceleration model [19, 20]. It is found that when electrons experience the interaction of two counter-propagating laser pulses, efficient acceleration of the electrons in plasma or a vacuum can be realized provided that some threshold amplitudes of the lasers are exceeded. In a real situation, the counter-propagating laser can be produced automatically during the propagation of an intense laser pulse in the plasma through processes such as stimulated Raman backscattering or reflected laser fields from high density regions. However, recently fast electron generation due to relativistic intensity laser–matter interaction influenced by preformed plasma has been addressed in a number of experimental and theoretical studies [21–28], suggesting that the presence of pre-plasma can significantly affect fast electron distributions. The new foundings in experimental and theoretical studies, addressing the influences of pre-plasmas, seem beyond the interpretations of stochastic heating and acceleration models. Both experiments and numerical simulations have reported an increase of fast electron generation efficiency with the increase of pre-plasma scale length. Recent particle-in-cell simulations [25] have observed super-high energetic electrons with cut-off kinetic energy as high as $100~\text{MeV}$ at laser intensity ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ and pre-plasma scale length $10~\mu \text{m}$. Paradkar et al extend the stochastic heating and acceleration model, suggesting that the electron acceleration in under-dense plasma is a stochastic process in the presence of two counter-propagating laser pulses and an electrostatic potential [25]. Kemp and his cooperators also realize that the electrostatic potential caused by low-density preformed plasma plays an important role in electron heating [28]. However the underlying physics, (i) the increase in the generation efficiency of energetic electrons with the increase of pre-plasma scale length and laser intensity, and (ii) the acceleration mechanism of super-high energetic electrons with kinetic energy greatly exceeding the ponderomotive energy, is still unclear. To thoroughly understand the underlying physics, a theoretical scaling law for the cut-off electron energy based on analysing the electrons' motion in the presence of counter-propagating lasers and external electric field in a dense plasma environment should be figured out, which is the aim of our present work.

In order to simulate laser–matter interaction with pico-second duration in the presence of large scale pre-plasma, we choose to use one-dimensional (1D) particle-in-cell (PIC) simulations, because it is computationally cheap. Although multi-dimensional effects, such as laser filamentation and self-focusing [29, 30], might play roles in these processes, they are neglected in the present work. We focus on the role of pre-plasma in energetic electron beam generation, using systematic particle-in-cell simulations. The questions, (i) 'why the generation efficiency of energetic electrons increases with the increase of pre-plasma scale length', and (ii) 'what is the underlying acceleration mechanism of super-high energetic electrons with kinetic energy greatly exceeding the corresponding laser ponderomotive energy', are answered. A two-stage acceleration model is proposed to identify the source of super-high energetic electrons. The first stage is synergetic acceleration by longitudinal electric field and two counter-propagating laser pulses; a scaling law is obtained theoretically, with its efficiency depending on the pre-plasma scale length and laser intensity. The second stage is related to the intense electrostatic potential building in front of the target and the accompanying reflection of electrons by this electrostatic potential barrier; the potential building and electron energy enhancement are figured out analytically and confirmed through electrostatic PIC simulations.

This paper is arranged as follows: the details of numerical modelling and simulation results are demonstrated in section 2. The two-stage acceleration model suggested by analysing the simulation results is proposed in section 3 to explain the impact of pre-plasma and identify the sources of energetic electrons. Conclusions and discussions are given in section 4.

The simulations are performed with 1D PIC code, which is a newly extended version based on LAPINE [31]. In order to simulate laser–matter interactions with large scale pre-plasmas, the weighted particle technique is applied in the numerical simulations, which has proven to be more efficient than uniform weighted particles in large density gradient calculations [32]. In addition, a 4th order particle cloud and 4th order FDTD method are applied in our simulations, because these features make it suitable for simulating laser–solid-density-plasma interactions at relativistic intensities [32]. The laser is of intensity ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ or normalized amplitude a  =  8.54 (with laser wavelength $1~\mu \text{m}$), entering the simulation box from the left boundary, where the laser amplitude rises over 33 fs to a  =  8.54 and then remains constant. To confirm that the rising amplitude is not the dominant factor, we also run simulations with rising duration of 100 fs, keeping other parameters the same; the results show exactly the same physics. The initial plasma density profile is taken as ${{n}_{\text{e}}}={{n}_{\text{solid}}}/\left(1+\exp \left[-2\left(z-{{z}_{0}}\right)/{{L}_{\text{p}}}\right]\right)$, where ${{n}_{\text{solid}}}=50{{n}_{c}}$ is the solid plasma density and L_p is the pre-plasma scale length. The ions are deliberately treated as immobile in the simulation to prevent gradual pre-plasma scale length variation and eliminate effects associated with ion mobility. As the charge separation electric field calculation is of great importance in our following analysis and the non-neutrality of total electric charge in the simulation box would lead to artificial over- estimation of the charge separation electric field, we choose to use a large simulation box instead of the reduced one to guarantee 'all' particles are confined in the simulation box without escaping. Furthermore, choosing a large simulation box can also avoid boundary effects caused by electrons recoiling due to the artificial electrostatic field on the boundaries, which could otherwise interrupt the physics we are studying. In the PIC simulations, the simulation box is of size $400~\mu \text{m}$, which is divided into $40\,000$ cells, with each cell containing 1000 electrons and 1000 ions. In our simulations, the region $0<z<100~\mu \text{m}$ is left as vacuum, L_p varies from ${{L}_{\text{p}}}=1~\mu \text{m}$, $5~\mu \text{m}$, $10~\mu \text{m}$ to $15~\mu \text{m}$, z₀ is fixed as $180~\mu \text{m}$ and the minimum plasma density is set as 0.001n_c for all simulation cases. In order to analyse the electron energy distributions in detail, we have placed two diagnostic planes to temporally record the electrons passing through. As shown in figure 1(a), the first diagnostic plane is located at $z=100~\mu \text{m}$ to record the electron going through in the z-direction, and the other one is located at $z=300~\mu \text{m}$, recording the electron passing through in the z-direction.

Zoom In Zoom Out Reset image size

To ensure the accuracy of the simulation, as we have done previously [33], we record the energy history of laser flux energy entering the simulation box, electromagnetic field energy in the simulation box, and particle kinetic energy in the simulation box. The total simulation time is set to be 400T₀, to avoid the electron recoil effect. Furthermore, we also loosen the simulation resolution by two, four and eight folds to confirm convergence of the simulation results. We find that the resolution demand is at least 25 cells per wavelength, provided that 4th order particle cloud, 4th order FDTD schemes and 1000 particles per cell are used.

The fast electron energy spectra obtained for different pre-plasma scale lengths (${{L}_{\text{p}}}=1~\mu \text{m}$, $5~\mu \text{m}$, $10~\mu \text{m}$ and $15~\mu \text{m}$) while keeping laser intensity fixed at ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$, are presented in figure 2. Solid lines in figure 2 record the energy spectra of electrons passing through the diagnostic plane located at $z=300~\mu \text{m}$. As we can see, there is a clear relation between cut-off kinetic energy and pre-plasma scale length—the larger the scale length the higher the cut-off kinetic energy—which is in agreement with earlier published works [25]. We have also found that the cut-off electron kinetic energy greatly exceeds the corresponding laser ponderomotive energy, which is about $3.8~~\text{MeV}$ at intensity ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$. For pre-plasma of scale length $1~\mu \text{m}$, $5~\mu \text{m}$, $10~\mu \text{m}$ and $15~\mu \text{m}$, the corresponding cut-off energies are $30~\text{MeV}$, $60~\text{MeV}$, $100~\text{MeV}$ and exceeding $120~\text{MeV}$ respectively. The dashed lines in figure 2 record the energy spectra of electrons passing through the diagnostic plane located at $z=100~\mu \text{m}$. By comparing the two energy spectra recorded by two different diagnostic planes, we can find that the cut-off energy recorded at $z=300~\mu \text{m}$ is several times the value recorded at $z=100~\mu \text{m}$. For pre-plasma of scale length $10~\mu \text{m}$, the cut-off energy recorded at $z=100~\mu \text{m}$ is about $50~\text{MeV}$, while that recorded at $z=300~\mu \text{m}$ is about $100~\text{MeV}$.

Zoom In Zoom Out Reset image size

We have found that the cut-off kinetic energy of electrons increases with the increase of the pre-plasma scale length. Meanwhile, we have also noticed that the cut-off electron kinetic energy recorded by the diagnostic plane located at $z=300~\mu \text{m}$ is several times that recorded at $z=100~\mu \text{m}$. The aim of this work is to uncover the mysteries, (i) the increase in the generation efficiency of energetic electrons with the increase of pre-plasma scale length and (ii) the source of super-high energetic electrons with energy greatly exceeding the corresponding laser ponderomotive energy. In order to understand the underlying physics of the observed phenomena, we now turn to analysing the trace of a particular electron as a function of time and the z-p_z phase-space dynamics. Figure 3 shows the trace of a particular electron initially located at $z=130~\mu \text{m}$, with laser of intensity ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ and pre-plasma of scale length $10~\mu \text{m}$. Figure 4 describes the phase-space patterns of laser–pre-plasma interactions with laser of intensity ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ and pre-plasma of scale length ${{L}_{\text{p}}}=5~\mu \text{m}$ (figure 4(a)) and $10~\mu \text{m}$ (figure 4(b)) respectively. The phase-space density $F\left({{n}_{\text{e}}}\right)$ gives a value proportional to the number of electrons found between z and $z+\text{d}z$ having longitudinal momentum ranged between p_z and ${{p}_{z}}+\text{d}{{p}_{z}}$. The normalized electrostatic potential, $-\text{e}\phi /{{m}_{\text{e}}}{{c}^{2}}$, due to the longitudinal charge separation field E_z, is shown in black curves covered on phase plots. The blue lines are the E_y components of the superposition of incoming and reflected laser pulses. The electron longitudinal momentum p_z is normalized by ${{m}_{\text{e}}}c$ and z is in the units of laser wavelength, which is $1~\mu \text{m}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in figure 3, in the very early stages of laser propagation in under-dense preformed plasma, some electrons are swept away in the forward direction by the laser ponderomotive force, leaving behind immobile ions. The electric field E_z due to charge separation within the under-dense plasma region tries to pull the electrons in the backward direction. When the laser arrives at the critical density surface and is reflected back, the ponderomotive force of the reflected laser pulse can further accelerate the electrons in the backward direction. Actually, the first stage acceleration is due to the synergetic effects by this longitudinal charge separation field E_z and the ponderomotive force of the reflected laser pulse. From the Woodward–Lawson theorem [34], we know that a single electron in vacuum, oscillating coherently with a propagating plane laser pulse would gain zero cycle averaged energy since the electron energy gain in one half cycle is exactly equal to the energy loss in the next half cycle. In the presence of plasma, the phase velocity of the light is decreasing/increasing with the decrease/increase of plasma density. Therefore, electrons can never continuously stay in phase with the laser beam. However, when there exists an external electric field [25, 35–37], even though this field is very week, the Woodward–Lawson theorem can be broken and the electron can obtain non-zero energy from the synergetic effects of the external electric field and the laser pulse. If the extension of the external electric field were infinite, the electron would always stay in phase with the laser and be accelerated to infinite energy.

When the incident laser arrives at the critical density surface and is reflected back, due to the formation of the steep interface of electron density [38], a strong delta-like charge separation field or step-like electrostatic potential, as shown in figure 4, is built up therein, and is strong enough to drive electrons to very high velocity within very short time and short distance. Imagine we are standing on the frame of a backward propagating electron, we will find that the incident laser pulse is oscillating very fast, and its only contribution to the motion of the electron is to increase its mass by a factor $\gamma ={{\left(1+{{a}^{2}}/2\right)}^{1/2}}$ in an average way (appendix A); however the reflected laser pulse is oscillating so slowly that this electron can be captured and continually accelerated backward by its ponderomotive force. Actually the first stage acceleration strongly depends on the pre-plasma scale length. As clearly demonstrated in figure 4(b) ($10~\mu \text{m}$), the first stage acceleration is stronger than that in figure 4(a) ($5~\mu \text{m}$). According to the Woodward–Lawson theorem [34], a single electron cannot gain non-zero cycle-averaged energy from one plane wave. However, in our case, there exists an external electric field E_z due to the charge separation in the under-dense pre-plasma region. Actually, as we shall analyse in the next section, the pre-plasma scale length determines the space extension of E_z, which eventually determines the maximum possible electron energy gain in this synergetic acceleration process.

The energetic electrons pre-accelerated in the first stage continuously propagate backward and expand freely, building up an intense electrostatic potential barrier therein, as shown by the black curves in figure 4. Actually the peak value of the electrostatic potential barrier is several times as large as the kinetic energy of these electrons pre-accelerated. As we know, for an electron with kinetic energy ${{E}_{k\text{in}}}$ initially located at position with zero electrostatic potential energy, it is impossible for it to arrive at the position with potential energy ${{U}_{p}}>{{E}_{k\text{in}}}$ without any external forces. However, for continuously emitting electron beams or separated multi electron bunches, we find that some electrons can arrive at positions where the potential energies are several times as large as their initial kinetic energies. When these electrons are reflected back to their original positions, the kinetic energies obtained will increase to ${{E}_{k\text{f}\max}}=N\times {{E}_{k\text{in}}}$, with N  >  1. Although it seems impossible, this process conserves the total energy of the system, and ${\sum}^{}{{n}_{\text{in}}}{{E}_{k\text{in}}}={\sum}^{}{{n}_{\text{f}}}{{E}_{k\text{f}}}$ is always satisfied, with ${{E}_{k\text{f}}}$ obeying ${{E}_{k\text{f}\min}}<{{E}_{k\text{in}}}<{{E}_{k\text{f}\max}}$. In the next section, solid interpretations are presented, including mathematical analysis and electrostatic numerical simulations, for the building process of electrostatic potential and the accompanying electron kinetic enhancement by reflection off this potential barrier.

The synergetic acceleration by longitudinal electric field and laser ponderomotive force—We consider the relativistic electron dynamics in the presence of two counter-propagating plane laser waves with vector potential a₊ and a₋ and longitudinal field E_z, as shown in figure 5. a₊ means the laser pulse is propagating in the same direction as an electron in the presence of electric field E_z. Note that in this section—analysing synergetic accelerations—for simplicity the z-direction is taken along the backward direction, i.e. the opposite direction to the PIC simulation configurations. Considering the electron propagates with high velocity along the z-direction, the only contribution of the incident wave a₋ is to increase the electron mass in an averaged way (appendix A). The z-momentum and energy equation, in normalized units, can be written as

$\begin{eqnarray}&&\frac{\text{d}\left(\gamma {{v}_{z}}\right)}{\text{d}t}=\frac{-1}{2\gamma}\frac{\partial a_{+}^{2}}{\partial z}+{{E}_{z}},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\frac{\text{d}\gamma}{\text{d}t}=\frac{1}{2\gamma}\frac{\partial a_{+}^{2}}{\partial t}+{{E}_{z}}{{v}_{z}},\end{eqnarray} \tag{ 2 }$

where v_z is the electron velocity component along the z-direction and the relativistic factor γ defined as $\gamma ={{\gamma}_{A}}{{\gamma}_{z}}$ with ${{\gamma}_{A}}={{\left(1+{{a}^{2}}/2+a_{+}^{2}\right)}^{1/2}}$, a²/2 is the average mass increase due to the incident laser wave of the form ${{a}_{-}}=a\sin (t+x)$, and ${{\gamma}_{z}}=1/{{\left(1-v_{z}^{2}\right)}^{1/2}}$.

Zoom In Zoom Out Reset image size

For a reflecting plane wave of the form ${{a}_{+}}={{a}_{+}}\sin (t-z)$, from equations (1) and (2), we find

$\begin{eqnarray}&&\frac{\text{d}}{\text{d}t}{{\gamma}_{A}}{{\gamma}_{z}}\left(1-{{v}_{z}}\right)=-{{E}_{z}}\left(1-{{v}_{z}}\right).\end{eqnarray} \tag{ 3 }$

Assuming electric field E_z to be constant, equation (3) can be integrated and we have

$\begin{eqnarray}&&{{\gamma}_{A}}{{\gamma}_{z}}\left(1-{{v}_{z}}\right)={{\sigma}_{{{\tau}_{0}}}}-{{E}_{z}}\left(t-{{t}_{0}}-z-{{z}_{0}}\right),\end{eqnarray} \tag{ 4 }$

where t₀ is the time at which the electron crosses z  =  z₀ and ${{\sigma}_{{{\tau}_{0}}}}={{\gamma}_{A}}{{\gamma}_{z}}\left(1-{{v}_{z}}\right){{|}_{t={{t}_{0}},z={{z}_{0}}}}$. Note for the highly relativistic case, we have ${{\sigma}_{{{\tau}_{0}}}}\sim (1/2)\left({{\gamma}_{A}}/{{\gamma}_{z}}\right)\ll 1$.

The trajectory of the electron z can be found by introducing a local time $\tau =t-z$, in which $\text{d}\tau /\text{d}\tau =\text{d}t/\text{d}\tau -\text{d}z/\text{d}\tau$ and $\text{d}t/\text{d}\tau =\left(\text{d}z/\text{d}t\right)\left(\text{d}t/\text{d}\tau \right)/{{v}_{z}}$, as

$\begin{eqnarray}&&\frac{\text{d}z}{\text{d}\tau}=\frac{{{v}_{z}}}{1-{{v}_{z}}}.\end{eqnarray} \tag{ 5 }$

Using v_z from equation (4), $\text{d}z/\text{d}\tau$ can be found to be

$\begin{eqnarray}&&\frac{\text{d}z}{\text{d}\tau}=\frac{1}{2}\left[\,{{f}^{2}}(\tau )-1\right],\end{eqnarray} \tag{ 6 }$

where $f(\tau )={{\gamma}_{A}}\left(\tau +{{\tau}_{0}}\right)/\left({{\sigma}_{{{\tau}_{0}}}}-{{E}_{z}}\tau \right)$.

The change in the electron energy only due to the contribution of laser waves is given by $\Delta \varepsilon (\tau )={{\gamma}_{A}}\left(\tau +{{\tau}_{0}}\right){{\gamma}_{z}}\left(\tau +{{\tau}_{0}}\right)-$ ${{\gamma}_{A}}\left({{\tau}_{0}}\right){{\gamma}_{z}}\left({{\tau}_{0}}\right)-{{E}_{z}}\left[z\left(\tau +{{\tau}_{0}}\right)-z\left({{\tau}_{0}}\right)\right]$.

Following the equation (4) and making use of the inequality (${{\sigma}_{{{\tau}_{0}}}}\ll 1$, ${{\sigma}_{\tau +{{\tau}_{0}}}}\ll 1$ and ${{E}_{z}}\tau \ll 1$), $\Delta \varepsilon (\tau )$ can then be rewritten as

$\begin{eqnarray} &&\Delta \varepsilon (\tau )=\frac{1}{2}{\int}_{0}^{\tau}\frac{\text{d}\gamma _{A}^{2}\left(\tau +{{\tau}_{0}}\right)/\text{d}\tau}{{{\sigma}_{{{\tau}_{0}}}}-{{E}_{z}}\tau}\text{d}\tau.\end{eqnarray} \tag{ 7 }$

Through equations (6) and (7), we can find the maximal-possible energy gain within the limited longitudinal scale length L and the maximal in-phase time $\tau =\pi /2$,

$\begin{eqnarray}&&L=\frac{1}{2E_{z}^{2}}\left[\frac{\gamma _{A}^{2}\left(\pi /2+{{\tau}_{0}}\right)}{{{\sigma}_{E}}-\pi /2}-\frac{\gamma _{A}^{2}\left({{\tau}_{0}}\right)}{{{\sigma}_{E}}}-a_{+}^{2}f\left({{\sigma}_{E}}\right)\right]-\frac{\pi}{4},\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray} &&\Delta \varepsilon (\pi /2)=\frac{a_{+}^{2}}{2{{E}_{z}}}f\left({{\sigma}_{E}}\right),\end{eqnarray} \tag{ 9 }$

where we define ${{\sigma}_{E}}={{\sigma}_{{{\tau}_{0}}}}/{{E}_{z}}\geqslant \pi /2$, and

$\begin{eqnarray}&&f\left({{\sigma}_{E}}\right)={\int}_{0}^{\pi /2}\frac{\sin (2x)}{{{\sigma}_{E}}-x}\text{d}x.\end{eqnarray} \tag{ 10 }$

As ${{\tau}_{0}}$ is just an arbitrary initial local time, for simplicity we set ${{\tau}_{0}}=0$ in the following expressions. Assuming $a\gg 1$, $L\gg 1$ and $a_{+}^{2}=R{{a}^{2}}$, where R is the reflection rate, based on equation (8) we can obtain,

$\begin{eqnarray}&&{{E}_{z}}=\frac{a}{{{L}^{1/2}}}{{\left[\frac{R{{\sigma}_{E}}+\pi /4}{2{{\sigma}_{E}}\left({{\sigma}_{E}}-\pi /2\right)}-\frac{R}{2}f\left({{\sigma}_{E}}\right)\right]}^{1/2}}.\end{eqnarray} \tag{ 11 }$

Combining equation (9) and equation (11), the maximal-possible electron kinetic energy gain within the limited longitudinal length L from the laser of incident amplitude a and reflection rate R can be expressed as

$\begin{eqnarray} &&\Delta \varepsilon =\eta a{{L}^{1/2}}=\frac{R}{2}\frac{f\left({{\sigma}_{E}}\right)}{g\left({{\sigma}_{E}}\right)}a{{L}^{1/2}},\end{eqnarray} \tag{ 12 }$

with ${{g}^{2}}\left({{\sigma}_{E}}\right)=\left(R{{\sigma}_{E}}+\pi /4\right)/\left[2{{\sigma}_{E}}\left({{\sigma}_{E}}-\pi /2\right)\right]-Rf\left({{\sigma}_{E}}\right)/2$.

In equation (12), the coefficient η is a function of R and ${{\sigma}_{E}}$. From figure 6, for the typical reflection rate R  =  0.9, α almost saturates at 0.5 for a large range of ${{\sigma}_{E}}$. Finally, we give a scaling law which describes the maximal-possible electron energy gain for the synergetic acceleration process, where the laser intensity I is normalized by $1.37\times {{10}^{18}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ and the longitudinal length $L\sim \beta {{L}_{\text{p}}}$ is normalized by $\mu \text{m}$,

$\begin{eqnarray}&&\varepsilon ~\left[\text{MeV}\right]=0.64\times {{\beta}^{1/2}}\times {{I}^{1/2}}\times L_{\text{p}}^{1/2}.\end{eqnarray} \tag{ 13 }$

Zoom In Zoom Out Reset image size

In equation (13), we assume that the longitudinal length is determined by the pre-plasma scale length with $L\sim \beta {{L}_{\text{p}}}$. Coefficient β can be obtained by the direct comparisons of equation (13) and values from PIC simulations, where the pre-plasma is assumed to be of the form $1/\left[1+\exp \left(-2z/{{L}_{\text{p}}}\right)\right]$ in the PIC simulations. According to the scaling law of equation (13), we can see that the first stage acceleration, or the synergetic acceleration by longitudinal electric field E_z and the ponderomotive force of the reflected laser, depends on both the incident laser intensity and the pre-plasma scale length.

In order to validate the PIC results with theory, a direct comparison of maximum electron energy observed in PIC simulation and theoretical scaling law is present in figure 7. The electron energy spectra are collected by the diagnostic planes located at $100~\mu \text{m}$, which could record the maximum electron energy accelerated in the first stage. Figure 7(a) corresponds to the cases with fixed pre-plasma scale length ${{L}_{\text{p}}}=10~\mu \text{m}$ and varying incident laser intensity, where the black line is for incident laser of intensity ${{10}^{19}}~\text{W}~\text{c}{{\text{m}}^{-2}}$, blue line is for ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$, and red line for ${{10}^{21}}~\text{W}~\text{c}{{\text{m}}^{-2}}$. Figure 7(c) corresponds to the cases with fixed incident laser intensity ${{10}^{20}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ and varying pre-plasma scale length, where black line is for pre-plasma scale length of ${{L}_{\text{p}}}=1~\mu \text{m}$, blue line is for ${{L}_{\text{p}}}=5~\mu \text{m}$, and red line for ${{L}_{\text{p}}}=10~\mu \text{m}$. Figures 7(b) and (d) are the comparisons between the theoretical scaling law and PIC simulation results, where black lines are the plot of equation (13) and red squares are the simulation results. The comparisons are of good fitness, provided the coefficient is set to be ${{\beta}^{1/2}}=2.5$.

Zoom In Zoom Out Reset image size

Electrostatic potential building and the accompanying electron reflection—To get insights on both (i) the possibility of the formation of the electrostatic potential barrier with the maximal value significantly larger than electron kinetic energy, and (ii) the role of the potential barrier in electron acceleration, let us consider a 1D model problem. Assume that at t  =  0 we have a bunch of electrons with density n_b and momentum p₀  >  0 occupying region 0  <  z  <  z_b (${{z}_{b}}\ll {{\lambda}_{\text{De}}}$) and a bunch of immobile ions, located at z  <  0, such that total electron and ion charges compensate each other. We consider dynamics of electron bunch expansion assuming that the electrons, which come back to their original positions, do not move any more. Since we are considering the 1D geometry, the electric field acting on an electron solely depends on its original position at t  =  0 and does not vary in time. Therefore, for an electron having $z(t=0)={{z}_{0}}<{{z}_{b}}$ we have the following equation of motion,

$\begin{eqnarray}&&\frac{\text{d}}{\text{d}t}\frac{p}{\sqrt{1-{{p}^{2}}}}=-{{E}_{z}}\left({{z}_{0}}\right),\end{eqnarray} \tag{ 14 }$

where ${{E}_{z}}\left({{z}_{0}}\right)$ is the original electric field, normalized by $e/{{m}_{\text{e}}}c$. From equation (14) we find the time dependence of the position z(t,z₀) of the electron initially located at z₀ as

$\begin{eqnarray}&&z\left(t,{{z}_{0}}\right)={{z}_{0}}+{\int}_{0}^{t}\frac{{{p}_{0}}-{{E}_{z}}\left({{z}_{0}}\right){{t}^{\prime}}}{\sqrt{1+{{\left[\,{{p}_{0}}-{{E}_{z}}\left({{z}_{0}}\right){{t}^{\prime}}\right]}^{2}}}}\text{d}{{t}^{\prime}}\\ &&={{z}_{0}}-\frac{1}{{{E}_{z}}\left({{z}_{0}}\right)}\left\{\sqrt{1+{{\left[\,{{p}_{0}}-{{E}_{z}}\left({{z}_{0}}\right)t\right]}^{2}}}-\sqrt{1+p_{0}^{2}}\right\},\end{eqnarray} \tag{ 15 }$

where p₀  =  p(t  =  0). From equations (14) and (15) one can easily see that within the setting of the problem the electron coming back to its original position has p  =  −p₀ and, therefore, returns to the original energy.

The original increase of the normalized electrostatic potential within the electron bunch, $\delta {{\phi}_{0}}$, can be easily found from Poisson equation,

$\begin{eqnarray}&&\delta {{\phi}_{0}}=\frac{1}{2}{{\left(\frac{{{\omega}_{\text{pe}}}{{z}_{b}}}{c}\right)}^{2}},\end{eqnarray} \tag{ 16 }$

where $\omega _{\text{pe}}^{2}=4\pi {{\text{e}}^{2}}{{n}_{b}}/{{m}_{\text{e}}}$. Now we will analyse time variation of the electrostatic potential at relatively large time $t>{{p}_{0}}/{{E}_{z}}\left({{z}_{0}}\right)$, when the majority of electrons have already returned to their original positions. Estimating the magnitude of ${{E}_{z}}\left({{z}_{0}}\right)$ from the Poisson equation, we can re-write this inequality as

$\begin{eqnarray}&&t>{{\tau}_{b}}=\frac{{{p}_{0}}c}{\omega _{\text{pe}}^{2}{{z}_{b}}}.\end{eqnarray} \tag{ 17 }$

Then the difference of the normalized electrostatic potential, $\Delta \phi (t)$, between the head of the expanding electron bunch, ${{z}_{h}}(t)=z\left(t,{{z}_{b}}\right)$, and the coordinate z_r(t) with ${{z}_{r}}=z\left(t,{{z}_{r}}\right)$ of electrons returning to their original position at time t, can be written as follows,

$\begin{eqnarray} &&\Delta \phi (t)={\int}_{{{z}_{r}}(t)}^{{{z}_{h}}(t)}{{E}_{z}}(z)\text{d}z,\end{eqnarray} \tag{ 18 }$

or

$\begin{eqnarray} &&\Delta \phi (t)=-{\int}_{0}^{{{E}_{z}}\left[{{z}_{r}}(t)\right]}{{E}_{z}}\left({{z}_{0}}\right)\frac{\text{d}z\left(t,{{z}_{0}}\right)}{\text{d}{{E}_{z}}\left({{z}_{0}}\right)}\text{d}{{E}_{z}}\left({{z}_{0}}\right)\\ &&=-\frac{1}{2}{{\left(\frac{{{\omega}_{\text{pe}}}}{c}\right)}^{2}}\left[z_{b}^{2}-z_{r}^{2}(t)\right]+{\int}_{0}^{2{{p}_{0}}}\frac{\sqrt{1+p_{0}^{2}}-\sqrt{1+{{\left(\,{{p}_{0}}-\xi \right)}^{2}}}}{\xi}\text{d}\xi.\end{eqnarray} \tag{ 19 }$

Since we are considering the time $t\gg {{\tau}_{b}}$ where ${{z}_{r}}(t)\to {{z}_{b}}$, we find the following asymptotic expression, $\Delta {{\phi}_{\infty}}= \Delta \phi (t\to \infty )$,

$\begin{eqnarray} &&\Delta {{\phi}_{\infty}}={\int}_{0}^{2{{p}_{0}}}\left[\sqrt{1+p_{0}^{2}}-\sqrt{1+{{\left(\,{{p}_{0}}-\xi \right)}^{2}}}\right]\frac{d\xi}{\xi}.\end{eqnarray} \tag{ 20 }$

From equation (20) we derive $\Delta {{\phi}_{\infty}}\sim p_{0}^{2}$ for ${{p}_{0}}\ll 1$ and $\Delta {{\phi}_{\infty}}\sim 2\ln (2)p_{0}^{2}$ for ${{p}_{0}}\gg 1$. In other words, for the non- relativistic case $\Delta {{\phi}_{\infty}}$ is twice the initial electron kinetic energy ${{E}_{k\text{in}}}$, while for a super-relativistic case $\Delta {{\phi}_{\infty}}\sim 2\ln (2){{E}_{k\text{in}}}\sim 1.4{{E}_{k\text{in}}}$.

As we mentioned before, electrons, being finally reflected back by potential, will come to their original positions and obtain their original kinetic energy. Thus in the process of launching just one electron bunch, there is no possibility of increasing electron energy. However, the situation changes drastically when we launch a few electron bunches separated by a dwell time ${{\tau}_{\text{dw}}}$. To get an insight into the electron acceleration mechanism, consider the case of two bunches. The first bunch, launched at t  =  0 will expand as was discussed before. At time $t={{\tau}_{\text{dw}}}>{{\tau}_{b}}$, the second bunch starts launching. At that time, the first bunch has formed a 'potential barrier' between the head of the first bunch and the launch point, with ${{\phi}_{\text{bar}}}= \Delta {{\phi}_{\infty}}$. However, almost all electrons of the first bunch have already come back to their original position and the electric field within the 'potential barrier' becomes very small, with $E\sim \Delta {{\phi}_{\infty}}/{{p}_{0}}t\ll {{E}_{z}}\left({{z}_{0}}\right)$. As a result, the second bunch also expands virtually into vacuum and at time $t=2{{\tau}_{\text{dw}}}$, the cumulative contribution of the first and second bunches will create a 'potential barrier' with ${{\phi}_{\text{bar}}}=2\times \Delta {{\phi}_{\infty}}$. In addition, a relatively small number of electrons at the head of the first bunch turning back after the expansion of the second bunch will finally acquire not only their initial kinetic energy but also potential energy created by the second bunch. As a result, their total kinetic energy as they reach the launching location will double. Their additional energy comes at the expense of electron energy from the second bunch, which are decelerated somewhat as the bunches pass through each other.

We can consider the injection of many identical electron bunches with the dwell time between them such that the previous bunches do not impact the dynamics of later ones. One can easily find that the number of such bunches is limited by ${{N}_{b}}\sim \ln \left(t/{{\tau}_{b}}\right)$. Therefore, maximum kinetic energy, acquired by the returned electrons of the very first bunch, after being accelerated by the electric field of all bunches can be estimated as ${{E}_{k\max}}\sim {{N}_{b}}\times {{E}_{k\text{in}}}$, which, nonetheless, can be significantly larger than ${{E}_{k\text{in}}}$.

We can consider also continuous injection of electrons into a half-space taking the time-dependent distribution function of launching electrons f(t, v). Considering the non-relativistic case we take

$\begin{eqnarray}&&f(t,v)=\frac{{{n}_{0}}\delta \left(v-{{v}_{0}}\right)}{1-\alpha t{{\omega}_{\text{pe}}}\left({{n}_{0}}\right)},\end{eqnarray} \tag{ 21 }$

where $\alpha \ll 1$. This temporal evolution of electron launch, limited by $\alpha t{{\omega}_{\text{pe}}}\left({{n}_{0}}\right)<1$, resembles the rate of bunch launches. The final energy by electric field of all bunches can be estimated as ${{N}_{b}}\times {{E}_{k\text{in}}}$, which, nonetheless, can also be significantly larger than ${{E}_{k\text{in}}}$.

In order to confirm the above theoretical analysis, we also run a serious of 1D electrostatic PIC simulations, which are solved by an energy conserving method (appendix B). The electrostatic PIC simulations solve the following equations,

$\begin{eqnarray}&&\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial z}+\frac{\text{e}{{E}_{z}}}{{{m}_{\text{e}}}}\frac{\partial f}{\partial v}=0,\end{eqnarray} \tag{ 22 }$

$\begin{eqnarray}&&\frac{\partial {{E}_{z}}}{\partial z}=4\pi \text{e}{\int}^{}f\text{d}v,\end{eqnarray} \tag{ 23 }$

$\begin{eqnarray}&&f(t=0)={{n}_{\text{e}}}\delta \left(v-{{v}_{0}}\right),\end{eqnarray} \tag{ 24 }$

with ${{\omega}_{pe0}}=4\pi {{n}_{0}}{{e}^{2}}/{{m}_{\text{e}}}$, $v=\bar{v}[c]$, $t=\bar{t}\left[1/{{\omega}_{pe0}}\right]$, $z=\bar{z}\left[c/{{\omega}_{pe0}}\right]$, ${{E}_{z}}=\overline{E}_{z}\left[c{{\omega}_{pe0}}{{m}_{\text{e}}}/e\right]$, $\phi =\bar{\phi}\left[{{m}_{\text{e}}}{{c}^{2}}/e\right]$, ${{\omega}_{\text{pe}}}=\overline{\omega}_{\text{pe}}\left[{{\omega}_{pe0}}\right]$, ${{n}_{\text{e}}}= {\overline{\omega}}_{{\text{pe}}^{2}}\left[{{n}_{0}}\right]$ and $f=\bar{f}\left[{{n}_{0}}/c\right]$. We define a reference density n₀, corresponding to a reference plasma frequency ${{\omega}_{pe0}}$. $1/{{\omega}_{pe0}}$ define the time scale in simulation, $c/{{\omega}_{pe0}}$ define the length scale and c is the speed of light. We can change the plasma density in simulation by adjusting $\overline{\omega}_{\text{pe}}$. If $\overline{\omega}_{\text{pe}}=1$, the plasma density used in simulation is exactly n₀, if $\overline{\omega}_{\text{pe}}=0.5$, the corresponding plasma density in simulation is $0.5\times 0.5\times {{n}_{0}}$.

Figure 8 shows the simulation results, in which an electron bunch of velocity v₀  =  0.7, thickness L₀  =  0.2 and plasma frequency ${{\omega}_{\text{pe}}}=0.5$ is emitted from the surface z  =  0. Figures 8(a)–(c) show the time-snap of z-v_z phase-space, electric field and potential profile at t  =  0.5, t  =  28 and t  =  80, which clearly demonstrates that at t  =  28, the electrons in the rear start returning to the emitting point at z  =  0, well consistent with the theoretical analysis, ${{\tau}_{b}}=\left(2/\omega _{\text{pe}}^{2}\right)\left({{v}_{0}}/{{L}_{0}}\right)=28$. In our simulation, we include a numerical friction mechanism to stopping electrons when re-entering into the emitting point. Figure 8(d) shows the maximal electric field and potential evolution with time, and we find that the maximal potential almost keeps constant even when the back edge of the bunch returns to the emitting point, which is also consistent with theoretical prediction. As expected by theoretical analysis, the maximal electric field decrease with time as ${{\tau}_{b}}/t$ when $t>{{\tau}_{b}}$. The kinetic energy of the returning electron is exactly equal to the initial value, having v  =  −v₀, which is, nonetheless, consistent with the theoretical prediction.

Zoom In Zoom Out Reset image size

Let us consider the situation of emitting multiple bunches. Figure 9(a) and (b) show the two bunch cases with the dwell time ${{\tau}_{\text{dw}}}=280$ greatly larger than ${{\tau}_{b}}=28$. We noticed that the maximal potential energy can be further increased by the emission of the second bunch, finally reaching four times the original kinetic energy. The velocity of the returned electron can be as high as v  =  −0.99 compared with the initial value v₀  =  0.7, confirming the theoretical prediction that the kinetic energy of the returning electron is doubled. Figure 9(c)–(f) are cases of three (${{\tau}_{\text{dw}1}}=200$ and ${{\tau}_{\text{dw}2}}=100$) and four (${{\tau}_{\text{dw}1}}=200$, ${{\tau}_{\text{dw}2}}=100$ and ${{\tau}_{\text{dw}3}}=50$) bunches, the maximal potential and the returning electron kinetic energy can be further increased as expected. Limited to the computational ability of our simulation, if the dwell time is long enough, the finial maximal potential energy will be close to the theoretically predicted value ${{E}_{k\max}}\sim \ln \left(t/{{\tau}_{b}}\right){{E}_{k\text{in}}}$.

Zoom In Zoom Out Reset image size

As shown in figure 4, the emission of electrons is a continuous process. Here in figure 10, we show the simulation results for a continuous electron beam with constant velocity v₀  =  0.4 and density profile $\omega _{\text{pe}}^{2}\exp \left(t{{\omega}_{\text{pe}}}\right)$, where ${{\omega}_{\text{pe}}}=0.125$. The simulation results, as shown in figures 10(a) and (b), indicate that the maximal potential energy is more than three times as large as the initial kinetic energy at t  =  40 and is still increasing gradually with time. Note that the oscillation of maximal potential energy, with the oscillation frequency increasing with time, comes from the plasma intrinsic oscillations, with its frequency determined by the density of the emitted electron beam. With the increase of electron beam density, the maximal potential energy and oscillation frequency are also increasing with time. Figure 10(c) records the velocity spectra of the returning electrons collected at the emission point, indicating that the returning electrons actually span a large velocity range, from  −0.3 to  −0.7. This spanning of the velocity range is also consistent with the theoretical prediction, with some of the electrons having velocity higher than the initial value 0.4, and some having velocity smaller than 0.4. The cut-off kinetic energy of the returning electrons can be about three times as large as their initial value.

Zoom In Zoom Out Reset image size

The generation of super-high energetic electrons influenced by pre-plasma in relativistic intensity laser–matter interaction is studied in a one-dimensional slab approximation with particle-in-cell simulations. Different pre-plasma scale lengths and incident laser intensity are considered, showing an increase in both particle number and cut-off energy of energetic electrons with the increase of the pre-plasma scale length and laser intensity, and the obtained cut-off energy of electrons greatly exceeding the corresponding laser ponderomotive energy. The two questions, (i) 'why is the generation efficiency of energetic electrons increasing with the increase of pre-plasma scale length and laser intensity', and (ii) 'what is the underlying acceleration mechanism of super-high energetic electrons with kinetic energy greatly exceeding the ponderomotive energy', are answered in this work.

Furthermore, a two-stage electron acceleration model is proposed to explain the underlying physics in detail. The first stage is attributed to the synergetic acceleration by the longitudinal charge separation electric field E_z and the ponderomotive force of the laser beams. The efficiency of the first stage acceleration depends on both the pre-plasma scale length and the laser intensity. The maximal possible energy gain during the first stage acceleration is analysed, and a scaling law is obtained by solving the relativistic electron motions in the presence of two counter-propagating plane laser waves and the external electric field due to the charge separation within limited space extension, which is on the order of the pre-plasma scale length. The maximal-possible energy gain in the first stage is estimated to be $\varepsilon ~\left(\text{MeV}\right)=1.6\times {{I}^{1/2}}\times L_{\text{p}}^{1/2}$, where I is laser intensity normalized by $1.37\times {{10}^{18}}~\text{W}~\text{c}{{\text{m}}^{-2}}$, and L_p is pre-plasma scale length normalized by $\mu \text{m}$. The scaling law indicates that with the increase of pre-plasma scale length and incident laser intensity, the maximal-possible electron energy also increases, which agrees well with the simulation results.

The energetic electrons pre-accelerated in the first stage could expand freely and build up an intense electrostatic potential barrier in front of the target, with potential energy several times as large as the electron kinetic energy. Some of the energetic electrons could be reflected by this potential and return back to the target, obtaining finial kinetic energies several times as large as the initial values. The processes of potential building and the accompanying electron kinetic enhancement by this potential barrier are analysed theoretically and confirmed by electrostatic PIC simulations, in which the theoretical predictions and electrostatic PIC simulations are also in good agreement.

Note that the modulational (M) and Raman (RS) instabilities may also play roles in this process. It was shown that [39–41], for a relativistic laser with laser amplitude a  >  1 and homogeneous density close to relativistic critical, the non-linear stage of the instabilities will result in a strong heating of the electron distribution function. The existing analyses are restricted to homogeneous plasmas. While in our case, plasma is rather strongly inhomogeneous, their inhomogeneous thresholds will differ a lot [42], so that if the existing conclusions match our cases need further investigation.

In experiment, reflected light does not always retrace the path of the incident light. However, the backward going electrons do not necessarily have an additional acceleration by reflected light. For those electrons which fail to get accelerated by the first stage, their initial energy acquired at the critical surface is enough to build up large potential, provided that we have a large number of electron bunches or sustained emission of electron beam, as shown in figures 9 and 10.

Furthermore, other multidimensional effects of laser propagation through under-dense plasmas, such as beam self-focusing [43] and self-generated electromagnetic fields [44, 45], could also enhance the electron energy by breaking the Woodward–Lawson theorem. These multidimensional effects are ignored in the present work. The extension of our work to multi-dimensional configurations and addressing the effects of self-focusing and self-generated electromagnetic fields shall be presented in a following separate paper.

This work was supported by the National Natural Science Foundation of China (11304331, 11174303, 61221064), the National Basic Research Program of China (2013CBA01504, 2011CB808104) and USDOE Grant DENA0001858 at UCSD.

We have studied the motion of a single electron in the field of a₊ , a₋ and E_z by numerically solving the 1D-3V electron equation of motion with the standard Boris algorithm. Figure A1(a) shows the motion of a single electron in the fields of only a₊ and E_z. It indicates that when the Woodward–Lawson theorem is broken, electrons will be continuously accelerated forward and the final kinetic energy increases with increasing acceleration length. Figure A1(b) shows when there exists two counter-propagating laser pulses, i.e. a₊ and a₋, the dynamics of the electron initially at rest is quite complicated, resulting in stochastic-like motions.

Zoom In Zoom Out Reset image size

However, when an electron with an initial large momentum p_z enters the fields of two counter-propagating laser waves and longitudinal electric field, the influence of the incident laser a₋ can be simplified. The only contribution of the incident laser wave a₋ is to increase the electron mass in an averaged way. In figure A1(c), black line shows the full dynamics of the electron under a₊ , a₋ and E_z, and the red line shows the dynamics of electron under only a₊ and E_z but replacing ${{\gamma}_{a}}={{\left(1+a_{+}^{2}+a_{-}^{2}\right)}^{1/2}}$ by ${{\gamma}_{a}}={{\left(1+a_{+}^{2}+{{a}^{2}}/2\right)}^{1/2}}$. The results of full dynamics and reduced model are well fitted, which confirms our assumption in the article well.

The majority of PIC schemes—including LAPINE, as we are using, have the property of exactly conserving the system total momentum, while not conserving the system total energy exactly. In fact, the typical PIC methods, which use explicit differentiation in time, i.e. explicit PIC schemes, tend to increase the total energy of the system by numerical heating. Conversely, the other categories of PIC methods, which use implicit differentiation in time, i.e. the so-called implicit PIC schemes, tend to decrease the total energy of the system by numerical cooling.

For large scale simulations, with $L\gg {{\lambda}_{d}}$, both explicit PIC and implicit PIC can be applied; this does not affect the overall simulation picture, provided that the energy conservation is satisfied within an acceptable error, like $\delta E/E\sim 1 \%$. However in our case, where we need to resolve significantly fine structures, with ${{L}_{0}}\ll {{\lambda}_{d}}$, the demand for energy conservation is of extremely high level.

To overcome this issue, we refer to the method introduced by Markidis and Lapenta [46], where a new PIC method, which conserves energy exactly, is used for the electrostatic PIC simulations, shown in figures 8–10. The equations of motion of particles and the Maxwell's equations are differenced implicitly in time by the mid-point rule and solved concurrently by a Jacobian-free Newton Krylov (JFNK) solver. The particle average velocities and the electrostatic field are calculated self-consistently by the JFNK solver to preserve the total energy of the system. The shortcomings of the method introduced are their difficulties in parallel implementation and extension to relativistic dynamics.