2 EV formation and the proton acceleration

Figure 1 shows the simulation set-up. A linearly (p-)polarised laser beam with normalised amplitude $a_{0}=30$ (intensity $I=1.2\times 10^{21}~\text{W}~\text{cm}^{-2}$ ), is incident at a grazing angle of $80^{\circ }$ on a plasma foil. Here $a_{0}=eE_{0}/m_{\text{e}}c\unicode[STIX]{x1D714}_{0}$ , $E_{0}$ , $\unicode[STIX]{x1D714}_{0}$ and $\unicode[STIX]{x1D706}_{0}=1~\unicode[STIX]{x03BC}\text{m}$ are the peak electric field, frequency and wavelength of the laser pulse, while $c$ , $e$ and $m_{\text{e}}$ denote the speed of light in vacuum, the elementary charge and electron mass, respectively. The focal spot and the full width at half maximum (FWHM) duration of the main laser pulse are $w_{0}=3\unicode[STIX]{x1D706}_{0}$ and $5.4T_{0}=18~\text{fs}$ , respectively, where $T_{0}$ is the laser period. The main laser pulse propagates to the right with a grazing angle with respect to the target surface. It first penetrates the underdense plasma region until it reaches the simulation box (marked by the black-dashed line in figure 1, there is a 2 $\unicode[STIX]{x03BC}$ m vacuum gap between the plasma and the left boundary where laser is injected), the density at the incident point is $n_{0}\approx 0.4n_{c}$ (in reality this is the position where the density becomes strongly inhomogeneous). The density profile within the simulation box is initialised with a density decreasing exponentially in the $y$ direction as $n_{\text{e}}(y)=100n_{\text{c}}\exp [-(y-y_{\text{t}})/\unicode[STIX]{x1D70E}]$ , where $n_{\text{c}}=m_{\text{e}}\unicode[STIX]{x1D714}_{0}^{2}/4\unicode[STIX]{x03C0}e^{2}$ is the critical density, $y_{\text{t}}=-16.9\unicode[STIX]{x1D706}_{0}$ is the target position and $\unicode[STIX]{x1D70E}=3\unicode[STIX]{x1D706}_{0}$ is the scale length. The density profiles in the tangential directions are assumed to be uniform since the scale lengths in those directions are much larger than $\unicode[STIX]{x1D70E}$ . Note that it is possible to achieve such density profile experimentally, as reported by Gauthier et al. (Reference Gauthier, Lévy, d’Humières, Glesser, Albertazzi, Beaucourt, Breil, Chen, Dervieux and Feugeas2014). The proton density is $n_{\text{p}}(y)=n_{\text{e}}(y)/55$ . The proton species is assumed to consist only a small fraction of the total positive charge, the majority of ions are considered to be immobile in this simulation. The effect of heavy ion motion will be addressed in the discussion section. The PIC simulations are performed with EPOCH (Arber et al. Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz and Bell2015). The simulation window $x\times y=200\unicode[STIX]{x1D706}_{0}\times 20\unicode[STIX]{x1D706}_{0}$ is sampled by $4000\times 400$ cells with 10 macro electrons and 20 macro protons in each cell.

Figure 1(b,c) shows the evolution of the current density ( $j_{x}$ ) and the out-of-plane magnetic field ( $B_{z}$ ) during the interaction. The first column ( $t=50T_{0}$ ) shows the laser channelling in the NCD plasma, and the associated strong currents inside the channel due to the laser–plasma interaction. In the NCD plasma, the laser effectively transfers all its energy to the electrons within a distance of a few tens of microns (Sylla et al. Reference Sylla, Flacco, Kahaly, Veltcheva, Lifschitz, Malka, d’Humières, Andriyash and Tikhonchuk2013). In the second column ( $t=100T_{0}$ ), the laser energy is completely depleted and the electron dynamics is dominated by self-generated fields. The electrons gyrate in the self-generated magnetic fields, forming a stable vortex structure. The two columns on the right ( $t=150T_{0}$ and $200T_{0}$ ) show that the EV continues to propagate along $x$ with a drift velocity ${\sim}0.28c$ . In the meantime, the background protons are accelerated in the $x$ -direction by the EV as shown in figure 1(d).

The acceleration of protons is due to the charge separation caused by the magnetic pressure in the EV. As shown by figure 2(a,b), the electron density inside the EV is lower than the ambient value. A high-density return current layer forms on the lower boundary of the EV. Components of the electrostatic field induced by the charge separation are presented in figure 2(c,d).

Figure 2. Features of the relativistic EV. (a) The electron density distribution, where (b) presents the one-dimensional cut of electron density along the horizontal (red) and vertical (blue) axes (marked as dashed in a). Electric field components in the directions (c) parallel and (d) normal to the target surface. (e) Proton population in momentum space $P_{x}$ – $P_{y}$ , and (f) energy spectrum of protons within opening angle $10^{\circ }$ (between blue dashed lines in e). The inset in (f) shows the phase space map $x$ – $P_{x}$ of these protons. All quantities are shown at simulation time $t=200T_{0}$ .

As the EV drifts, the protons are accelerated in the leading half of the structure. Those protons that gain enough energy can keep up with the EV (they are captured), otherwise they lose energy due to the decelerating field in the second half of the EV. The captured protons can be continuously accelerated until they overrun the vortex with approximately twice its drift velocity.

The electric field in the direction normal to the target surface ( $E_{y}$ ) is defocusing, therefore only a fraction of the captured protons can achieve the maximum energy gain, others are scattered out of the acceleration phase before they can do so. This makes the energy distribution depend strongly on the pitch angle as indicated by figure 2(e), where the slope of the dots shows the pitch angle $\unicode[STIX]{x1D703}_{\text{p}}=P_{y}/P_{x}$ of the corresponding protons. The energy spread of the accelerated protons can then be adjusted by selecting different pitch angles (e.g. with a collimator): at large $\unicode[STIX]{x1D703}_{\text{p}}$ a wide spectrum is obtained with a sudden cutoff, while protons at low $\unicode[STIX]{x1D703}_{\text{p}}$ exhibit a quasi-monoenergetic feature. External focusing with various well-established approaches (Toncian et al. Reference Toncian, Borghesi, Fuchs, d’Humières, Antici, Audebert, Brambrink, Cecchetti, Pipahl and Romagnani2006; Albertazzi et al. Reference Albertazzi, d’Humières, Lancia, Dervieux, Antici, Böcker, Bonlie, Breil, Cauble and Chen2015; Zhai et al. Reference Zhai, Shen, Borghesi, Wang, Zhang, Kar, Ahmed, Li, Li and Zhang2019) can be reliably employed in order to obtain a well-collimated beam for applications such as cancer therapy.

Figure 2(f) presents the energy spectrum of the protons within a $10^{\circ }$ opening angle in the forward direction. The proton beam has a mean energy ${\sim}140$  MeV and an energy spread of approximately 10 %. In this quasi-monoenergetic beam the proton number per 1 % energy range and solid angle increment $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}$ is estimated to be  ${\sim}10^{6}/\text{msr}/1\,\%\unicode[STIX]{x1D716}_{\text{p}}$ (assuming the proton beam to be axisymmetric with respect to the $x$ axis), which is comparable to contemporary laser-driven ion acceleration schemes at their cutoff energy (Schreiber et al. Reference Schreiber, Bolton and Parodi2016), and this can be further improved by increasing the proton-to-heavy-ion ratio in the plasma.

In order to obtain a deeper insight into the electron dynamics, we track 2000 electrons throughout the simulation. These electrons are chosen to be located inside the EV at $t=150T_{0}$ , with half of them in the forward-going current distributed throughout the structure (hereafter laser-induced current) and the other half in the return current. Their positions at $t=150T_{0}$ and $t=400T_{0}$ are represented by blue and red circles in figure 3(a), respectively. Most of the electrons in the laser-induced current stay within the EV, but the electrons in the return current are lost when they reach the end of the EV. This can also be seen from the representative orbits shown by the black curves.

Figure 3. The dynamics of the trapped electrons in the EV. (a) The upper and lower plots show the positions of 1000 tracked electrons (each) in the laser-driven current and return current, selected at $t=150T_{0}$ . The red and blue circles show their position at $t=150T_{0}$ and $t=400T_{0}$ , respectively. A representative electron trajectory in each case is plotted with the black curve, i and f mark the initial and final points of the trajectory. (b) A comparison of the average drift velocity of electrons (blue line) and the prediction from (2.1) (red line), the black dashed line shows the $\boldsymbol{E}\times \boldsymbol{B}$ drift velocity (first term in (2.1)). All quantities are taken from the simulation at $t=200T_{0}$ , along $x=83.5\unicode[STIX]{x1D706}_{0}$ .

An individual electron drifts in the electric and magnetic fields self-generated by the vortex, where the drift velocity can be obtained by the sum of the $\boldsymbol{E}\times \boldsymbol{B}$ drift and grad- $B$ drift. (Note that only the electrons are magnetised, so the protons and ions are not subject to $\boldsymbol{E}\times \boldsymbol{B}$ drift in the time scale of interest)

(2.1) $$\begin{eqnarray}\displaystyle \frac{\boldsymbol{v}_{\text{d}}}{c}=\frac{\boldsymbol{E}\times \boldsymbol{B}}{B^{2}}-\left(\unicode[STIX]{x1D6FE}_{\text{e}}-\frac{1}{\unicode[STIX]{x1D6FE}_{\text{e}}}\right)\frac{m_{\text{e}}c^{2}\boldsymbol{B}\times \unicode[STIX]{x1D735}B}{2eB^{3}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{\text{e}}$ is the relativistic factor of electrons. By assuming $\unicode[STIX]{x1D6FE}_{\text{e}}=\unicode[STIX]{x1D6FE}_{a}\equiv \sqrt{1+a_{0}^{2}/2}$ in the laser-induced current and $\unicode[STIX]{x1D6FE}_{\text{e}}\sim 1.4$ in the relatively cold return current, the prediction of (2.1) matches the PIC simulation data as shown by figure 3(b). Moreover, it can be seen that the $\boldsymbol{E}\times \boldsymbol{B}$ drift is the dominating term in (2.1); the grad- $B$ drift only becomes important in the upper boundary of the EV, where the gyro-radii of electrons are not negligible compared to the variation scale of magnetic fields.

3 Model of EV velocity in strongly non-uniform plasmas

Since the laser-induced electrons propagate with the EV (see figure 3(a) and supplementary movie, available at https://doi.org/10.1017/S0022377819000485), their average drift velocity must equal the EV propagation speed. Due to the approximate anti-symmetry of $E_{x}$ , shown in figure 2, the drifts in $y$ direction must average to zero and the EV moves in a direction consistent with $\boldsymbol{B}\times \unicode[STIX]{x1D735}n$ , as found in weakly inhomogeneous plasmas (Richardson et al. Reference Richardson, Angus, Swanekamp, Ottinger and Schumer2013; Angus et al. Reference Angus, Richardson, Ottinger, Swanekamp and Schumer2014).

We model the EV as two oppositely travelling electron streams: the electrons going backwards are in a thin layer with high density, and the ones going forward are in a larger region but with significantly lower density, as shown in figure 4(a). For simplicity, the electron and ion densities are taken to be uniform in the return current layer: $n_{\text{e}}=n_{\text{e}1}$ and $n_{\text{i}}=n_{\text{i}1}$ , and the $y$ coordinate is set to be zero at the bottom of the EV. The interface between the two oppositely travelling currents is at $y=y_{1}$ and the upper boundary of the EV is at $y=y_{2}$ (where the gyro-radii of the electrons become comparable to the size of EV). As the proton density is negligible, the positive charge density is $Zen_{\text{i1}}\exp [-(y-y_{1})/\unicode[STIX]{x1D70E}]$ , where $Z$ is the charge number and $n_{\text{i}1}$ is the immobile ion density at $y=y_{1}$ . Furthermore, we assume the shape of EV to be elongated in the $x$ direction (which is increasingly well satisfied for high $a_{0}$ ), effectively reducing the problem to one dimension ( $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x\ll \unicode[STIX]{x2202}/\unicode[STIX]{x2202}y$ ). The relevant electric and magnetic fields at $y_{1}$ are $E_{y}=E_{1}$ , $B_{z}=B_{1}$ , and $E_{2}$ , $B_{2}$ at $y_{2}$ , respectively, where $|B_{2}|\ll |B_{1}|$ and $|E_{2}|\ll |E_{1}|$ .

Figure 4. Analytical model to calculate the EV drift velocity. (a) The one-dimensional EV model with two oppositely travelling electron streams. The drift velocity of EV and the maximum velocity of the accelerated protons (black squares) are plotted as functions of (b) $a_{0}$ , (c) $\unicode[STIX]{x1D703}_{0}$ and (d) $\unicode[STIX]{x1D70E}$ , respectively. These parameters are set to their default values $a_{0}=30$ , $\unicode[STIX]{x1D703}_{0}=80^{\circ }$ and $\unicode[STIX]{x1D70E}=3\unicode[STIX]{x1D706}_{0}$ , when not being scanned. The red diamonds and dashed lines in (b–d) represent the drift velocities in PIC simulations and the fit suggested by (3.2).

The electric and magnetic fields satisfy Maxwell’s equations $\unicode[STIX]{x2202}E_{y}/\unicode[STIX]{x2202}y=4\unicode[STIX]{x03C0}e(Zn_{\text{i}}-n_{\text{e}})$ and $\unicode[STIX]{x2202}B_{z}/\unicode[STIX]{x2202}y=-4\unicode[STIX]{x03C0}j_{x}/c\approx -4\unicode[STIX]{x03C0}en_{\text{e}}v_{\text{d}x}/c$ . Integrating these two equations in the region of laser-induced current gives an estimate for the drift velocity of the EV

(3.1) $$\begin{eqnarray}\displaystyle U_{\text{d}}=\frac{\displaystyle \int _{y_{1}}^{y_{2}}v_{\text{d}}n_{\text{e}}\text{d}y}{\displaystyle \int _{y_{1}}^{y_{2}}n_{\text{e}}\text{d}y}\approx \left(\frac{E_{1}}{B_{1}}+\frac{4\unicode[STIX]{x03C0}eN_{\text{i}}}{B_{1}}\right)^{-1}c, & & \displaystyle\end{eqnarray}$$

where $N_{\text{i}}=Zn_{\text{i}1}\unicode[STIX]{x1D70E}(1-\exp [(y_{1}-y_{2})/\unicode[STIX]{x1D70E}])$ is the ion areal density inside the EV. Thus, the drift velocity of an EV depends on $N_{i}$ as well as the electric and magnetic fields at $y_{1}$ .

The first term in the bracket on the right-hand side of (3.1) can be obtained by integrating Maxwell’s equations in the return current layer (note that the $\boldsymbol{E}\times \boldsymbol{B}$ drift is the dominating term in this region, as illustrated in figure 3(b), so that $v_{\text{d}}$ can be estimated by $cE_{y}/B_{z}$ ), which gives $E_{1}/B_{1}=-\sqrt{1-Zn_{\text{i}1}/n_{\text{e}1}}$ . From this it follows that $|E_{1}/B_{1}|$ varies over a range of ${\sim}$ 0.1, that is typically an order of magnitude smaller than the variation of the second term in (3.1). As a representative value, we assume $Zn_{\text{i}1}\approx n_{\text{e}1}/2$ , yielding $E_{1}/B_{1}\approx -0.7$ .

An estimate of the second term in (3.1) relies on the fact that the amplitude of $B_{1}$ in the EV approximately equals the magnetic field in the laser created channel (see figure 1). Consider the case where all the electrons in the channel move forward with velocity ${\sim}c$ , and the bottom of the channel is at the position $y=y_{1}$ . Then $B_{1}\approx 4\unicode[STIX]{x03C0}eZn_{\text{i}1}\unicode[STIX]{x1D70E}[1-\exp (-R_{\text{ch}}/\unicode[STIX]{x1D70E})]$ , where $R_{\text{ch}}=\unicode[STIX]{x1D706}_{0}\sqrt{\unicode[STIX]{x1D6FE}_{a}n_{\text{c}}/4\unicode[STIX]{x03C0}^{2}n_{\text{m}}}$ is the characteristic radius of the laser-driven channel (Borisov et al. Reference Borisov, Shiryaev, McPherson, Boyer and Rhodes1995; Pukhov, Sheng & Meyer-ter Vehn Reference Pukhov, Sheng and Meyer-ter Vehn1999), and $n_{\text{m}}$ is the local plasma density corresponding to the maximum laser penetration depth in the NCD plasma.

From Snell’s law, we obtain $\unicode[STIX]{x1D702}\sin \unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D702}_{0}\sin \unicode[STIX]{x1D703}_{0}$ , where $\unicode[STIX]{x1D702}\approx 1-n_{\text{e}}/(2\unicode[STIX]{x1D6FE}_{a})$ is the refractive index of the plasma [the initial value is $\unicode[STIX]{x1D702}_{0}\approx 1-n_{0}/(2\unicode[STIX]{x1D6FE}_{a})$ , and $\unicode[STIX]{x1D6FC}$ is the angle between the laser propagation direction and $y$ axis. Thus, by setting $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}/2$ , the density corresponding to the maximum penetration of the laser is obtained as $n_{\text{m}}/n_{\text{c}}=2\unicode[STIX]{x1D6FE}_{a}-(2\unicode[STIX]{x1D6FE}_{a}-n_{0}/n_{\text{c}})\sin \unicode[STIX]{x1D703}_{0}$ .

For large EVs ( $N_{\text{i}}\approx Zn_{\text{i}1}\unicode[STIX]{x1D70E}$ ), the drift velocity in terms of laser–plasma parameters can then be approximated as

(3.2) $$\begin{eqnarray}\displaystyle U_{\text{d}}\approx c\left[\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D706}_{0}}\sqrt{2-\left(2-\frac{n_{0}/n_{\text{c}}}{\sqrt{1+a_{0}^{2}/2}}\right)\sin \unicode[STIX]{x1D703}_{0}}-0.7\right]^{-1}. & & \displaystyle\end{eqnarray}$$

The most important prediction of (3.2) is that the drift velocity, and so the maximum proton energy, does not depend crucially on the laser intensity, since the plasma density at the incident point is typically $n_{0}/n_{\text{c}}\ll a_{0}$ . Therefore, controlling the laser incident angle and the plasma density gradient can be effectively used to adjust the proton energy. This agrees well with the PIC simulations as shown by figure 4(b–d), $U_{\text{d}}$ depends strongly on the laser incident angle and plasma scale length, and weakly on the laser intensity (the scaling flattens further as $n_{0}$ decreases).

The maximum proton speed obtained in the PIC simulations are overlaid on the plot in figure 4(b–d) with black squares, and they are approximately twice the corresponding drift velocity. There are cases where no accelerated proton energy is plotted above $U_{\text{d}}$ . This indicates that protons cannot be captured by the EV for those parameters, either because the acceleration field $E_{x}$ is too weak (in figure 4(b) when $a_{0}=10$ ) or because the EV is drifting too fast (for the cases when $U_{\text{d}}>0.4c$ ). In this case one could increase the values of $a_{0}$ , $\unicode[STIX]{x1D70E}$ , or reduce the incidence angle $\unicode[STIX]{x1D703}_{0}$ as suggested in figure 4(b–d).

Unlike most other proton acceleration mechanisms, which require high laser intensity to reach high ion energies, the scheme proposed here provides accessible degrees of freedom to control the acceleration process, allowing accelerating protons to hundreds of MeV with a narrow energy spread. According to figure 4 the laser parameters required to produce a 100–200 MeV monoenergetic, self-injected proton beam are within reach of existing laser facilities.