Laser acceleration of monoenergetic protons via a double layer emerging from an ultra-thin foil

We present theoretical and numerical studies of the acceleration of monoenergetic protons in a double layer formed by the laser irradiation of an ultra-thin film. The ponderomotive force of the laser light pushes the electrons forward, and the induced space charge electric field pulls the ions and makes the thin foil accelerate as a whole. The ions trapped by the combined electric field and inertial force in the accelerated frame, together with the electrons trapped in the well of the ponderomotive and ion electric field, form a stable double layer. The trapped ions are accelerated to monoenergetic energies up to 100 MeV and beyond, making them suitable for cancer treatment. We present an analytic theory for the laser-accelerated ion energy and for the amount of trapped ions as functions of the laser intensity, foil thickness and the plasma number density. We also discuss the underlying physics of the trapped and untrapped ions in a double layer. The analytical results are compared with those obtained from direct Vlasov simulations of the fully nonlinear electron and ion dynamics that is controlled by the laser light.

3 plasma expansion and shock acceleration mechanisms by means of PIC simulations, while Yin et al [39] found a laser 'break-out afterburner' regime, where the laser penetrates the target and accelerates electron beams, which in turn accelerates the ions to hundreds of MeV by a relativistic Buneman instability.
Macchi et al [40] presented 1D and 2D PIC simulation studies where they observed high intensity ion bunches to be generated at the circularly polarized laser-plasma interaction surface and to be moving into the plasma. To obtain nearly monoenergetic ions, Klimo et al [41] Robinson et al [42] and Ryokovanov et al [43] considered the radiation pressure acceleration (RPA) of the thin film as a whole, using circular and linear polarization of the laser light. They demonstrated that circular polarization is crucial for efficient laser acceleration of a thin foil; in the case of linear polarization there is an oscillating component in the ponderomotive force term that leads to strong chaotic heating of the electrons followed by a foil explosion and broadening of the energy spectrum of the ions by the TNSA process. There is an optimal thickness of the target equal to the distance of maximum charge separation, suggested also by Yan et al [44]. Liseykina et al [45] presented a parametric numerical study of optimal target thickness for proton acceleration. The effects of multiple spatial dimensions on the stability of the accelerated foil have been investigated by Klimo et al [41], Robinson et al [42] and Pegoraro and Bulanov [46], showing that the foil can be fractured due to a Rayleigh-Taylor-like instability, leading to the bunching of the foil and the broadening of the ion energy spectrum. Properly tailored laser pulses with a sharp intensity rise may stabilize the foil and make it possible to efficiently accelerate monoenergetic proton beams in multiple dimensions [41,46].
Yet some questions remain: how does the plasma maintain itself in such a seemingly stable state with density higher than a solid even in 1D? What traps the protons, preventing them from overshooting the electron front? Recently, we proposed [47] the self-organized double layer as a stably trapping mechanism in 1D. The electrons are trapped by the potential well formed by the balance of the radiation pressure and the electric force of the protons left behind, while protons are trapped in the accelerating frame by its inertial force while being accelerated by the electric field of the electron layer. This double layer allows the trapped protons to be accelerated as a monoenergetic group as the foil is accelerated by the radiation force as a whole. But a key question remains: what fraction of the protons are trapped and thus accelerated as monoenergetic ones? In this paper, we address these questions as well as the nature of proton trapping and test the role of the target thickness by employing an analytic and numerical study of the 1D Vlasov-Maxwell equations for protons and electrons.

Physics of thin foil acceleration
The acceleration of protons in a thin foil by intense laser light is based on the following physics. The foil is irradiated by an intense laser light that completely ionizes the foil to form a slab of plasma, and pushes the electrons forward towards the front end of the plasma slab by the radiation pressure. Einstein [48] first derived the radiation pressure for arbitrary incidence on a moving mirror using Lorentz transformation and the invariance of Maxwell's equations. This work has been extended to arbitrary moving objects by Censor [49], and the subject also has been treated in textbooks [50,51]. Debye [52] considered radiation pressure in his work on the physics of comets, while Marx [53] proposed that interstellar vehicles could be propelled by terrestrial laser beams. For our purposes, we will consider the radiation pressure of normal 4 incidence on a perfectly reflecting moving mirror [47] which takes into account the decrease of the radiation pressure due to the Doppler shift of the reflected light. Here I 0 is the intensity of the laser beam, v is the speed of the foil and c is the speed of light in vacuum. The ions are accelerated by the space-charge electric field induced by the separation of the electrons and ions in plasmas. This leads to an acceleration of the plasma slab as a whole by the radiation pressure. The RPA has been discussed earlier in [37,42,53]. We assume a thin plasma slab of an equal number N 0 of the ions and electrons, so that its total rest mass density is N 0 m i , where m i is the ion rest mass, in a 1D geometry along the z-axis. The equation of motion for the plasma slab is then where T L is the temporal laser period, γ = (1 − v 2 /c 2 ) −1/2 is the relativistic gamma factor, and we have denoted the dimensionless parameter Equation (2) can be integrated to give the speed [47] } and h(t) = 6Pt/T L + 4. The position z(t) of the plasma slab (center of mass position) is found from dz/dt = v, which combined with (2) yields Integration of equation (5) gives where v is given by (4). The kinetic energy of a typical ion in the slab is These are the ions trapped in a self-organized double layer that is formed by the balance of the electrostatic attractive potential of the electrons and the inertial force of the accelerating frame. These trapped ions form a blob of monoenergetic ions in phase space. Now we wish to examine the conditions for the validity of the assumption of plasma slab accelerated as a whole, pushed by the radiation pressure, since the slab is in fact a plasma. The rigidity is due to the formation of a self-arranged accelerated double layer, which traps the ions and electrons. The formation of the electron layer is due to the balance between the laser ponderomotive force pushing the electrons forward and the electrostatic force by the ions behind, pulling the electrons back. For the ion layer, it is the balance between the electrostatic 5 force of the electrons accelerating them and the inertial force of the acceleration holding them back. The effective potential well of the double layer traps the ions and moves them together as a whole system, which is accelerated by the radiation pressure.
Due to the acceleration, there will also be a group of ions that are untrapped and left behind the accelerated layer, making the plasma slab negatively charged. An ion will be untrapped if the acceleration by the space charge electric field is weaker than the mean acceleration of the whole plasma slab. It is interesting to consider the forces on an ion in the moving frame of the foil accelerated by the laser radiation. In this frame, there will be an accelerating force by the electrostatic field and a decelerating inertial force. The balance between these two forces determines whether the ions will be trapped or untrapped. Hence, we can define the effective potential as a sum of the electrostatic and inertial potentials, φ (z, t) = φ(z, t) + (F rad /eN 0 )z. In a frame moving with the speed v(t) of the foil, the equation of motion for the ion is then given by where E = −∂φ/∂z is the electric field. For applications of the ion acceleration, it is interesting to estimate how large a portion of the ions is trapped in the plasma slab and are accelerated to monoenergetic energies and how many ions are untrapped and spread out in energy. The loss of the ions will stop when the electric field behind the plasma slab is large enough, say due to untrapped ions, so that the acceleration of an ion by this electric field is as large as the mean acceleration of the plasma slab. We assume a thin plasma slab of initially an equal number of ions and electrons, in a 1D geometry along the z-axis. Then the electric field behind the accelerated layer is obtained from the Poisson equation ∂ E/∂z = en i / 0 as where N un = z 0 0 n i dz is the integrated (along the z-axis) number of untrapped ions and z 0 is the boundary between trapped and untrapped ions. Using E = E un in (8), we have If the right-hand side of (10) is larger than zero, then the ion is falling into the potential well and can be considered as trapped, while if it is less than zero, it can be considered as untrapped.
The number of trapped ions can be found when m i d( The fraction of the number of trapped ions to the total number of ions is In particular, we note that all ions are untrapped (N un = N 0 ) when 6 and for F rad > e 2 N 2 0 / 0 , one would have a complete separation between the electrons and ions. In this case, the electrons will be accelerated in front of the ions by the radiation pressure, while the left-behind ions will experience a Coulomb explosion due to the electrostatic repulsion. This will result in a large spread in ion energies and will decrease the usefulness of these ions in applications where monoenergetic ions are needed.

Comparison between theory and numerical modeling of the acceleration of a thin plasma slab
To assess the validity of the present theory, we have performed a series of Vlasov simulations and have compared the numerical results with those deduced from analytical formulae presented in section 2. While many simulation studies of proton acceleration have been carried out using PIC simulations, we believe that this is the first time a grid-based Vlasov solver has been used for this purpose. Therefore, it is worthwhile to discuss the difference between the two methods. In the PIC method [54], the Vlasov equation is solved by following the trajectories of a set of statistically distributed super-particles, which resolves the particle distribution functions in phase space. Each super-particle represents a large number of real particles. PIC simulations have proven to be extremely successful due to their relative simplicity and adaptivity, especially in problems involving large-amplitude waves and beams, and in multiple spatial dimensions. However, the statistical noise of PIC simulations sometimes overshadows the physical results, and for some problems, the low-density velocity tail of the particle distribution cannot be resolved with high enough accuracy by the super-particles. As a contrast, the grid-based Vlasov solver treats the particle distribution function as a continuous phase fluid that is represented on a grid in both space and velocity (or momentum) space. The advantage with grid-based Vlasov solvers is that there is no statistical noise in the simulations, and that the dynamical range is determined by the number system of the computer rather than super particles. Hence, the low-density velocity tail of the particle distribution can be resolved much more accurately by grid-based Vlasov solvers, which makes them suitable for certain types of problems. Simulations in higher dimensions using grid-based Vlasov solvers are very challenging since the full phase-space has to be represented on a grid, which makes both the storage of the data in the computer's memory, and the numerical calculations, extremely demanding. One also has to be aware of the tendency of the distribution function to become oscillatory in velocity space, which can lead to unphysical noise and recurrence effects in the numerical solution unless proper smoothing or viscosity is used in velocity space [55]. In our simulations, we used proper numerical viscosity in configuration space and momentum space to minimize these effects.
Let us present the relevant nonlinear kinetic model describing the proton acceleration resulting from the nonlinear interaction between intense laser light and a collisionless plasma. A more detailed derivation of the governing equations is given in the appendix. The electromagnetic wave gives rise to the relativistic electron mass increase m e γ e and a relativistic ponderomotive force [1], F e = −m e c 2 ∂γ e /∂z, where m e is the electron rest mass, and γ e = (1 + p 2 z /m 2 e c 2 + e 2 |A| 2 /m 2 e c 2 ) 1/2 is the relativistic electron gamma factor and p z is the z component of the electron momentum. The relativistic ion mass is m i γ i , and the ponderomotive force acting on the ions reads F i = −m i c 2 ∂γ i /∂z, where m i is the ion rest mass and γ i = (1 + p 2 z /m 2 i c 2 + e 2 |A| 2 /m 2 i c 2 ) 1/2 is the relativistic ion gamma factor. For the ions, the relativistic effects due to the quivering of the ions in the radiation field (proportional to |A| 2 inside the ion gamma factor) can usually be neglected due to the large mass of the ions compared to that of the electrons. The complex envelope of the perpendicular (to the z-axis) vector potential of the circularly or linearly polarized electromagnetic wave is given by the wave equation is the local plasma frequency, accounting for the relativistic mass increase and plasma density variations, and ω 0 is the laser angular frequency. For circularly polarized light, we will have |A| 2 = |A| 2 , where we note that the oscillatory part of the light cancel exactly. On the other hand, for linear polarization, we would instead have |A| 2 = |A| 2 [1 + cos(2ω 0 t − 2φ)]/2, where φ is the complex phase of A. We note that in this case, |A| 2 has an oscillatory part with a frequency twice that of the laser light, which will enter into the ponderomotive force for the electrons [42,43]. We will consider circularly polarized light below unless stated otherwise.
The parallel electric field is obtained from the Maxwell equation with initial conditions obtained from Gauss' law where we have introduced the electrostatic potential φ and is the electric charge density. The electromagnetic field is coupled nonlinearly with the ions and electrons, whose dynamics is governed by the ion and electron Vlasov equations in the lab frame and respectively. Equations (14), (16), (19) and (20) form a closed set for our purposes. Similar models have been used to study electron phase-space structures [56]- [59]. In our simulations, we used a simulation box in z space of size 200 c/ω 0 , which was resolved with 2000 grid points. For the electrons, we used a momentum space spanning ±10 m e c, resolved with 60 grid points, while the ion momentum space was taken from −30 m e c to +1470 m e c, resolved with 6000 grid points. We used absorbing boundary conditions for the Vlasov equations such that the electrons or ions hitting the right or left boundary were absorbed and removed from the simulation. For the electromagnetic wave, we used an outflow boundary condition on the right boundary, the electromagnetic wave of amplitude A 0 was 8 injected at the left boundary while absorbing outgoing waves. The initial conditions for the ions and electrons were taken to be a Maxwellian distribution of the form f i,e (z, p z ) t=0 = {n i,e (z)/[(2π ) 1/2 α i,e ]}exp[ − p 2 z /2α 2 i,e ] with α i = 0.8 m e c and α e = 0.4 m e c, and with the initial ion and electron number density n i,e = n 0 in the interval −L/2 < z < L/2 and n i,e = 0 elsewhere. In our computer simulations, we used a compact Padé scheme [60] to approximate the p z and x derivatives, while we used the standard fourth-order Runge-Kutta scheme for the time-stepping.
The laser intensity is related to the laser amplitude through where a 0 = e|A 0 |/m e c is the normalized amplitude of the injected laser light. We assume a thin plasma slab with the electron number density n 0 and the width L, so that the total number of ions is N 0 = n 0 L. In terms of these parameters, the dimensionless parameter P in (3) is where ω pe = (n 0 e 2 / 0 m e ) 1/2 is the plasma frequency of the overdense plasma slab, λ L = 2π c/ω 0 is the laser wavelength in vacuum, ω 0 is the laser angular frequency and m i /m e = 1836 is the proton to electron mass ratio. In terms of our parameters, the expression for the fraction of the trapped ions (12) can be written as We note that the fraction of the trapped ions is smallest initially when the speed of the foil v is small and increases for larger values of v. Hence, a fraction of the ions initially untrapped may later catch up and again become trapped. In our simulations, we used a plasma slab whose density n 0 was 10 times larger than the critical density n c = 0 m 2 e ω 2 0 /e 2 of the laser light, so that ω 2 pe /ω 2 0 = 10, and a laser amplitude of a = 5. We made three runs with different initial thicknesses of the plasma slab, with thickness L = 0.1λ L , L = 0.2λ L and L = 0.4λ L . We note from (23) that if no ions will be initially trapped and the theory predicts a separation between the electron and ion layers. For ω 2 pe /ω 2 0 = 10, the electron and ion layers will be separated if a 0 > 4.5 for L = 0.1 λ L , a 0 > 9 for L = 0.2 λ L , and a 0 > 18 for L = 0.4 λ L . Hence, for the chosen value a 0 = 5, we are above the threshold for complete separation of the electron and ion layers for the thinnest slab L = 0.1 λ L , while we are below the threshold for L = 0.2 λ L and L = 0.4 λ L .
We first discuss the case L = 0.2λ L . This is close to the optimal thickness L = (ω 2 0 /ω 2 pe )a 0 λ L /π ≈ 0.16λ L for the monoenergetic ion acceleration [47]. In figures 1-4, we have presented details of the simulations at different times. At time t = 3.4 T L , we see in figure 1(a) that the laser amplitude has set up a standing wave pattern behind the plasma layer, with an amplitude of a ≈ 10 of its first maximum, i.e. 2 times that of the injected wave. A small portion of the laser light, with an amplitude of a ≈ 1, has also tunneled through to the right side of the plasma layer. The electrons, seen in figure 1(b), are pushed forward by the laser light but are kept back by the space charge electric field. The charge density in figure 1(e) shows a bipolar structure where the left-hand side is positively charged and the right-hand side negatively charged, giving rise to the localized positive electric field in figure 1(e) and the associate doublelayer structure of the potential in figure 1(f). The effective trapping potential for the ions given in equation (8), which is derived from the sum of the electrostatic and inertial forces, is depicted in figure 1(c). The latter shows a local minimum at z ≈ 0.1λ L and a maximum at z ≈ 0. The ions that are to the left of the potential maximum are starting to lag behind the accelerated plasma layer, and can be considered as untrapped. In figure 1(c), we see a bunch of monoenergetic ions at z = 0, which coincide with the left edge of the electron layer in figure 1(c). This ion bunch is accelerated uniformly by the space charge electric field. At t = 16 T L , displayed in figure 2, the effective potential has deepened and broadened further with a minimum at z = 1.5 λ L and the local maximum at z = 1 λ L , making the trapped ions more deeply trapped. The untrapped population of the ions, clearly seen in figure 2(d) between z = 0 and z = 1.5 λ L , form an almost homogeneous positive charge density seen in figure 2(e) and associated linearly increasing positive electric field and concave potential profile, seen in figures 2(f) and (g). Note that most of the ions are trapped in the bottom of the effective potential well at z = 1.5 λ L . They have reached a momentum of ∼ 0.2 m i c, corresponding to an energy of ∼ 20 MeV.
Turning to the results at t = 35 T L , shown in figure 3, we observe that the effective potential has deepened and broadened further. Here, the trapped protons have been accelerated to momenta of 0.35 m i c corresponding to a proton energy of ∼ 55 MeV. We see in figure 3 that a portion of the previously untrapped ions begin to fall into the potential and to overtake the trapped ion population at z ≈ 6.5 λ L . This is due to the decrease of the radiation pressure at higher speeds of the plasma layer by the Doppler shift of the reflected light, giving rise to the factor (1 − v/c)/(1 + v/c) in the expression for the radiation pressure (1) and in the equation of motion (2). Also here the untrapped population of the ions, seen in figure 3(d) gives rise to   For the thinnest plasma slab, L = 0.1λ L , illustrated in figure 5, we observe that a few laser periods after the laser light has reached the plasma slab, the electrons have not been separated from the ions, but instead the laser light is tunneling through the plasma barrier. The normalized amplitude a = e|A|/m e c of the laser light is of the order a ≈ 5 inside the electron slab, so that the relativistic electron mass increase is a factor ≈ 5, making the plasma less overdense and the thin layer almost completely optically transparent. At this point, the proton acceleration stops and the fastest protons have gained energy of only ∼ 0.2 MeV.
Finally, we display the results with the thicker slab L = 0.4 λ L at t = 73 T L in figure 6. We see that the accelerated ions are less coherent with a larger spread in momentum (and ion energy) than for the thinner foil in figure 4. Initially, we, similarly to Macchi et al [40], observed 'pulsed' acceleration and the production of ion bunches directed towards the plasma slab. The acceleration is slower due to the larger inertial of the thicker slab in figure 6. Here, we could measure that about 95% of the ions are trapped in the effective potential well.
In the RPA discussed here, the laser beam accelerates the whole plasma slab (ions + electrons) as a whole. For this to work, we clearly need the electron pressure to be sufficiently low so that the plasma slab does not expand by the electron thermal pressure to make the foil transparent, and thus circularly polarized light is preferred since this minimizes the electron heating. Hence, the radiation pressure should be larger than the electron thermal pressure, 2I 0 /c > n 0 T e . From the electron distribution function (panel b in figures 1-4 and 6), we can estimate the electron thermal energy to be of the order T e ∼ 1m e c 2 , which should be less than (using the expression (21) for I 0 ), 2I 0 /(n 0 c) = 2 0 ω 2 0 m 2 e c 3 a 2 0 /(e 2 n 0 c) = 2(ω 2 0 /ω 2 pe )a 2 0 m e c 2 = 5m e c 2 , for our parameters ω 2 0 /ω 2 pe = 1/10 and a 0 = 5. Hence, the radiation pressure is here ∼ 5 times larger than the electron thermal pressure and the foil is prevented from undergoing thermal expansion. For the electrons, we instead have a balance between the radiation pressure and the electrostatic force as discussed above. This is approximately what is predicted by equation (23) at small speeds (v c). As the accelerated layer increases in speed, the efficiency of the radiation pressure decreases due to the Doppler shift effect. This means that a portion of the untrapped ions starts to overtake the accelerated layer and to be re-trapped as discussed in figures 4 and 6. However, the re-trapping of the untrapped ions is a relatively slow process and we see only a slight increase of the trapped ions in the simulations at later times, where the formula (23) tends to overestimate the number of the trapped ions at later times.

Summary and discussion
In summary, we have presented theoretical and numerical studies of ion acceleration in an ultra-thin plasma foil by intense laser light. In the RPA, the radiation pressure pushes the electrons forward into the foil until they are stopped by the induced space-charge electric field. On the other hand, the electric field accelerates the ions, which are held back by the inertial force in the accelerated frame. As a result, both the electrons and ions are trapped, and the thin foil is accelerated as a whole. We have derived analytic formulae for the ion energy of the trapped ions, which depend on the laser intensity and its duration, foil thickness and the plasma number density. Also discussed is the physics of trapped and untrapped ions in the double layer, where, due to the acceleration, some ions are untrapped and spread out in energy. The obtained theoretical results show excellent agreement when compared against direct Vlasov-Maxwell simulations of the ion and electron dynamics in the laser field. Our numerical Vlasov simulation results are generally in agreement with previous works using PIC simulations. The present investigation indicates that monoenergetic ions can be accelerated to energies of the order of 100 MeV in 100 laser periods, making them suitable for medical treatment. Even though investigations in multiple dimensions show that the foil is unstable to a Rayleigh-Taylor-like instability [41,42,46], properly tailored laser pulses may stabilize the foil long enough to accelerate it efficiently [46]. Klimo's work [41] also show that the acceleration of ∼ 100 MeV quasi-monoenergetic proton beams remains a possibility in multiple dimensions. Further investigations are underway and will be published elsewhere.
BE acknowledges the support and hospitality of the University of Maryland where this work was carried out. We will consider a 1D geometry along the z-axis. In this geometry, it is convenient to introduce the scalar and vector potentials, φ and A, respectively, through their relations to the electric and magnetic fields where A is the vector potential of the electromagnetic wave, which in the Coulomb gauge (∇ · A = 0) is given by A =xA x +ŷA y . Herex,ŷ andẑ are the unit vectors along the x, y and z directions in a Cartesian coordinate system. We then have from the Gauss law respectively.
In the limit of linearly polarized light, e.g. A = (1/2)A(z, t)x exp(−iω 0 t) + complex conjugate, equation (A.17) would still be valid, but we would have |A| 2 = |A| 2 [1 + cos(2ω 0 t − 2φ)]/2, where φ is the complex phase of A. Hence, the ponderomotive force term will in this case have an oscillatory part with a frequency twice that of the laser light, which may lead to chaotic heating of the electrons [42,43].