Transverse instability in the "light sail" ion acceleration
Y. Wan I. A. Andriyash W. Lu W. B. Mori V. Malka


 


[]yang.wan@weizmann.ac.il
Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 7610001, Israel


Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 7610001, Israel



Department of Engineering Physics, Tsinghua University, Beijing 100084, China


University of California Los Angeles, Los Angeles, CA 90095, USA


Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 7610001, Israel
Laboratoire d'Optique Appliquée, Ecole polytechnique - ENSTA - CNRS - Institut Polytechnique de Paris, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France














Acceleration of ultrathin plasma foils by laser radiation pressure promises compact alternatives to the conventional ion accelerators. It was shown, that a major showstopper for such schemes is a strong transverse instability, which develops the surface ripples, and is often attributed to the Rayleigh-Taylor (RT) type. However, simulations indicate, that these perturbations develop the features, that cannot be consistently explained by the RT mechanism. Here we develop a three-dimensional (3D) theory of this instability, which shows that its linear stage is mainly driven by strong electron-ion coupling, while the RT contribution is actually weak. Our model provides the instability spectral structure and its growth rate, that agrees with the large scale 3D particle-in-cell simulations. Numerical modeling shows, that target destruction results from a rapid plasma heating induced by the instability field. Possible paths to instability mitigation are discussed.




Transverse instability in the "light sail" ion acceleration
    V. Malka
    December 30, 2023
===========================================================





§ INTRODUCTION

Interacting with a thin foil, an ultra-intense laser pulse can accelerate ions on a very short distance, which makes such schemes attractive for the applications ranging from radiography of high-energy density targets <cit.> to radiotherapy <cit.>. Among various acceleration mechanisms, one promising option is to push a solid target directly by the radiation pressure force <cit.>. In the case of a very thin nanofoil, it can be pushed as a single sheet, which is known as the “light sail” (LS) acceleration <cit.>. Ideally, in a one-dimensional LS process, efficient and quasi-monoenergetic ion acceleration is predicted by using a normally incident circularly polarized laser, where 𝐉×𝐁 heating is strongly suppressed <cit.>. In the realistic conditions, the surface deformations from a nonuniform laser profile enable more efficient plasma heating <cit.>. Tailoring the laser and target profiles allows to create a uniform radiation pressure distribution, preventing the macroscopic surface bending <cit.>. However, the local surface distortions produced by the plasma instability, turn out to be a greater problem (see Fig. <ref>).

Simulations and experiments indicate a strong transverse instability in LS process, which generates surface and density perturbations of the plasma foil <cit.>. This process is not yet well understood, and the mechanism of Rayleigh-Taylor (RT) instability driven by the light pressure was suggested as a possible  explanation <cit.>. Simulations indicate, that the developing fluctuations have a well pronounced periodic structure, while the property of RT is to grow faster for the shorter wavelengths (in absence of stabilizing mechanisms). In <cit.>, this contradiction was addressed by suggesting, that the scale of the produced surface ripples is imposed by the diffraction patterns of the laser field on the surface, and therefore should be close to the laser wavelength. In its turn, this consideration also contradicts the numerical simulations, where the modulation scales were shown to vary with different laser plasma parameters <cit.>. More recently an alternative explanation suggested, that this instability does not origin from RT process, but is induced by a strong coupling between oscillating electrons and quasi-static ions <cit.>.



In this paper we propose a generalized model that includes all possible effects in the full three-dimensional (3D) geometry. The new theory allows to self-consistently account for the RT development, and compare the contributions of both mechanisms. For the first time we obtain analytically the instability growth rate and its mode structure, which are fully consistent with the 3D particle-in-cell (PIC) simulations. In the second part of the paper we show, that the instability itself does not destroy the target, but saturates when entering nonlinear stage. During this stage, the density fluctuations induce a rapid plasma heating, which triggers the merging of the developed small structures and finally breaks the whole target thus letting the laser penetrate through the plasma. In the final part of the study, we discuss the possibilities to minimize the instability effect following from the model.



§ THREE-DIMENSIONAL INSTABILITY MECHANISM




 §.§ Numerical demonstration


Let us start the discussion by showing a typical LS scenario in a highly resolved 3D PIC simulation, using a fully relativistic electromagnetic code OSIRIS <cit.>. In this simulation, a circularly polarized (CP) laser pulse with a frequency ω_0 and normalized amplitude a_0=5 is launched at t=0 along z-axis, from the left boundary of the simulation domain with the grid resolutions of 0.04 c/ω_0, 0.04 c/ω_0 and 0.02 c/ω_0 in x, y and z directions, where x, y represent transverse directions and z is the longitudinal axis; c is the speed of light in vacuum. For a more accurate comparison with the following theoretical analysis, we consider a laser with the flat-top temporal profile bounded by 15 c/ω_0 rise and fall ramps, and uniform transverse profile. The foil target is placed at z=5 c/ω_0, and is modeled as a pre-ionized uniform electron-proton plasma with the thickness l_0=c/ω_0, density of 10n_c. Plasma is initialized with equidistantly distributed macro-particles – 64 particles per species per cell, and with electron temperature of 10 eV, which provides initial perturbations to seed the instability.





Fig. <ref>a shows a 3D distribution of proton density within the region of its full width half maximum (FWHM), at t=60 ω_0^-1, when the foil is pushed by the laser pulse. At that time one can see, that both target density and shape are strongly modulated, and the average proton density in the compressed foil reaches, n_p ≃ 23 n_c. The (x-y) density distribution averaged along z is plotted in Fig. <ref>b, demonstrating a quasi-regular structure of the emerging perturbations. This structure can also be analyzed in the spectral space, and in Fig. <ref>c we show the 2D Fourier transformation of Fig. <ref>b. The unstable mode has a ring structure with the central radial wavenumber k_m = 5.3 ω_0/c, and its amplitude can be used to describe instability dynamics, as shown by a black curve in Fig. <ref>a. At early times, t≲ 25 c/ω_0, initial thermal energy of electron plasma provides a quick rise of charge density fluctuations <cit.>. After that, the resonant mode k_m is excited and grows exponentially, δ n∝exp(γ_m t), with the rate γ_m = 0.18 ω_0 (red dashed line in Fig. <ref>a).

The electron energy spread shown with the blue curve in Fig. <ref>a, characterizes plasma heating dynamics. One can see, that during the exponential instability growth, the heating is effectively prevented. Instability gets saturated at t≳ 70 c/ω_0, and from this moment the electron energy spread starts to grow. Rapid heating disperses electrons around the target, and leads to the target breaking at t≃ 100 c/ω_0. This final stage is also observed in proton acceleration dynamics shown in Fig. <ref>b. During the linear and saturation stages, acceleration is uniform, and results in a quasi-monoenergetic proton spectrum. At the final stage, the foil surface is no longer sharp, and does not support the radiation pressure acceleration. At that point, the acceleration is governed by the plasma expansion process, which drastically increases the proton spectral bandwidth, washing off the quasi-monoenergetic feature.



 §.§ Analytic model


Let us now describe the linear stage of instability using the relativistic two-fluid equations of electrons and ions. For the derivation simplicity we consider a reference frame, which follows the accelerated plasma foil. This frame is non-inertial as it experiences acceleration, α_in≡dv_in/dt = P_rad/(m_i n_0 l_0), where P_rad = 2 I_las/c is the radiation pressure produced by a laser with intensity I_las, m_i is the ion rest mass, n_0 and l_0 correspond to the non-perturbed target density and thickness respectively. In this analysis we consider only electrostatic fluctuations along the laser polarizations, and neglect the small contributions of electromagnetic modes <cit.>.

The continuity and transverse motion equations for the cold electron and ion fluids read:


    ∂ n_s∂ t + ∇· (n_s 𝐯_s)=0  , 
       ∂𝐩_⊥ ,s∂ t + 𝐯_s·∇𝐩_⊥ ,s = q_s 𝐄_⊥ ,

where 𝐯_s, 𝐩_s, n_s and q_s denote respectively the velocity, momentum, density and charge of the specie s (either electrons or ions). In the accelerating frame, the longitudinal ion motion is governed by the pressure P_rad and the inertia force m_iα_in,

    ∂ p_z ,i∂ t=P_radn_i l_0-m_iα_in ,

and the charge densities of the ion and electron fluids are coupled via the Poisson equation:

    ∇·𝐄=4π (q_in_i-en_e) .



To simplify further derivations, we make a few assumptions. Firstly, let us consider ions to be non-relativistic, which is well justified at the beginning of the LS process, when the instability is already developing, while particle energies are still relatively low. We also note, that surface distortions produced by the instability have typically much greater size than the thickness of the compressed foil. This allows us to assume plasma to be infinitely thin, and consider all related quantities to only depend on (x, y, t). Following the linear perturbation theory, we divide all field and plasma quantities into the zero- and first-order parts, f=f_0+f̃, and neglect the terms ∝f̃^2 and f̃^3. Electric fields of a CP laser E_x0=E_0sin(ω_0 t) and E_y0=E_0cos(ω_0 t), and the corresponding electron oscillations define the zero-order values of the transverse field and electron velocity v_e0 = v_⊥ e 0 respectively, where ω_0 is the laser frequency in the co-moving frame.

We also note, that the transverse force in Eq. (<ref>), corresponds to a case of the ideally flat foil. Adding the surface distortions z̃, produces the transverse force, -P_rad∇_⊥z̃, in the first order component of Eq. (<ref>) for electrons:

    ∂𝐩̃_⊥ ,e∂ t + e 𝐄̃_⊥  = - 𝐯_e0·∇𝐩̃_⊥ ,e  - P_rad∇_⊥z̃ .

The two terms on the right hand side of Eq. (<ref>) represent the contributions of the electron-ion coupled parametric instability <cit.>, and the Rayleigh-Taylor (RT) instability driven by the the radiation pressure coupled to the foil distortions <cit.> respectively.

After all mentioned considerations, a linear set of equations can be obtained for the first-order terms for both fluids. Eliminating the terms ñ_i, 𝐯̃_i, 𝐩̃_i, 𝐄̃_⊥ and z̃, the equations for electron density and momenta fluctuations can be obtained. Assuming further, that all quantities are ∝exp[i(k_x x + k_y y - ω t)], we can linearize the derivatives, and finally obtain the equations for the Fourier components of ñ_e, p̃_x ,e and p̃_y ,e (in following we suppress the subscript e):

    Ω_g ñ_ω+ω_0 + Ω_g^* ñ_ω-ω_0 +  n_0 κ( k_x p̃_x, ω+k_y p̃_y ,ω) = ωñ_ω ,
       Ω_gp̃_x, ω+ω_0 + Ω_g^*p̃_x, ω-ω_0 + k_x ñ_ωχ(ω) =ωp̃_x, ω ,
       Ω_g p̃_y, ω+ω_0  + Ω_g^* p̃_y, ω-ω_0 + k_y ñ_ωχ(ω)= ωp̃_y, ω ,

where we denote Ω_g=v_os(k_x-ik_y)/2, κ=(2-v_os^2)/2γ_0, with v_os=P_os/γ_0 being the electron quiver velocity amplitude and γ_0 is the electron's zero-order Lorentz factor. We have also defined here a dispersion function:

    χ(ω) = α_in^2 m_i ω_pi^2n_e0ω^2 (ω_pi^2-ω^2) - 4πω^2(k_x^2+k_y^2)(ω_pi^2-ω^2) ,

where ω_pi is the ion plasma frequency. The first term in Eq. (<ref>) is responsible for the RT process, and the second one defines the electron-ion coupled instability. In the expressions above, mass is in units of electron rest mass m_e, velocity in units of c and momentum in units of m_e c.

Equations (<ref>) are written for a mode ω, and they couple it with the up- and downshifted modes ω±ω_0. Following the standard parametric analysis <cit.>, we write these equations for the modes ω±ω_0, and, neglecting the off-resonant terms ω± 2ω_0, obtain a closed set of equations for ñ, p̃_x and p̃_y for the modes ω, ω+ω_0 and ω-ω_0. Since these equations are homogeneous (have no constant terms), the non-trivial solution exists only for the coefficient matrix with a zero determinant, which gives a dispersion relation between ω and k. More details on the derivation can be found in the Supplementary Material <cit.>.

The resulting dispersion equation can be solved numerically, and in Fig. <ref>d we show distribution of the imaginary part of ω (i.e. instability growth rate) in (k_x, k_y) space. For the coefficients in Eqs. (<ref>) we use the values of the electron and ion densities n̅_e=55 n_c and n̅_i=23 n_c extracted from the simulation and averaged in the FWHM region, and γ̅_0=2.2 is calculated from the simulated electron spectrum. Figures  <ref>c and <ref>d clearly demonstrate the similar ring structures, and the simulated radius value 5.3 ω_0/c has a good agreement with the theoretical value 5.8 ω_0/c. The growth rate given by the maximal imaginary part of this numerical solution is γ_m = 0.21 ω_0, also in good agreement with the simulated value (see Fig. <ref>a). Considering separately two terms in Eq. (<ref>), one can estimate the contributions of each effect in the instability processes: γ_ei= 0.2 ω_0 is a growth rate of electron-ion coupling instability and γ_RT=0.06 ω_0 corresponds to the RT process.

Making a few more assumptions for the numerical solution of the dispersion equation (similar to ones in section III of <cit.>), we obtain simpler expressions for the maximal growth rate:


    γ_m≈√(γ_ei^2+γ_RT^2) ,

and the wavenumber of the fastest growing mode:

    k_xm^2+k_ym^2≃2κω_pe^2/v_os^2 ,


where the contributions of two instabilities to the growth rate are:

    γ_ei≃ 2 (ω_pi^2ω_pe)^1/3(κ m_e / m_i )^1/6 ,
          γ_RT≃( √(κ/2) α_inω_pe / v_os)^1/2 ,

Equation (<ref>) defines an isotropic mode structure in the (k_x, k_y) space in the form of a ring with the radius of k_r=√(2κ)ω_pe/v_os. In the 1D case, k_y=0, Eq. (<ref>) is reduced to the expression for k_m obtained in Ref <cit.>. In a case γ_0≫1, Eq. (<ref>) can be also simplified as k_r≃ω_pe/c√(γ_0), corresponding to the relativistic plasma skin depth. In most cases, γ_ei is much larger than γ_RT, indicating the former dominates the whole instability development. For a low laser intensity (a_0<1), n_e is close to the value of n_i, and the expression of γ_ei is simplified as 2ω_pi.

Note, that for non-relativistic interaction, the expressions Eqs. (<ref>) and (<ref>) can be directly applied in the laboratory reference frame, while for the higher ion energies these estimates should be Lorentz-transformed.



§ NONLINEAR INSTABILITY EVOLUTION


In Fig. <ref>, we see that after the linear instability phase its exponential growth saturates. Meanwhile, electron heating starts, and until certain level it does not affect acceleration of ions and their spectral bandwidth. As electron temperature grows, particles spread further from the target, and density of electron plasma decreases. At that point laser penetrates through the target, and the LS process can no longer continue. This moment we call the target breaking, and for the optimal experiment design it is important to understand this process better.



To explore the induced electron heating, we perform a 2D PIC simulation with the same physical parameters as before, but finer resolutions Δ y=Δ z=0.02 c/ω_0 and 80 particles per species in each cell. In Fig. <ref>a we look at the electron transverse phase space (y-p_y) at the time when instability saturates. Here we see the periodic transverse modulations produced by the instability, and at the peaks of this mode one can see the particles escaping the modulated structure, which corresponds to the wave-breaking process. These higher energy electrons travel between the modulation nodes, and their dynamics can be observed by following a sample electron trajectory, as shown in Fig. <ref>b. During the linear instability phase (blue part), electron oscillates around its initial position in the laser field with a transverse momentum p_y≃ 2 m_e c, defined by the partially screened laser field. At the saturation stage (green part), the instability field becomes strong enough to force electron to drift away from its initial position along the surface. Finally, near the target breaking (red part) laser field fully penetrates through the plasma, and electron can reach the momentum up to 6 m_e c.

As we see, the instability saturation time is related to the modulation wave-breaking. At the linear stage, the  first-order perturbations of electron motion are confined in the local oscillating electric field with a typical wave number k_m. These perturbations can be expressed as p̃≃Δ psinω_0tsink_m y, where Δ p is the amplitude of the first-order momentum. As instability develops, Δ p grows exponentially, and so does the coordinate oscillation amplitude Δ y∼Δ p/m_eκω_0, until it reaches the plasma potential width, k_mΔ y∼1. At that point, electrons start to escape the unstable plasma wave, which corresponds to the wave-breaking process. Combining this relation with Eqs. (<ref>), one can estimate the saturated density perturbations Δ n_sat, as

    Δ n_sat/n_0≃Δ p_sat/p_0≃2 ω_0 γ_0^7/2/ω_pe (γ_0^2+1)^3/2 .

Eq. (<ref>) predicts that plasma heating onsets, when Δ n_e/n_e≃ 0.3, which compares reasonably well with the simulation value of 0.2 at saturation. We note, that at this time, density fluctuations are still small and the target is opaque to the laser, such that the LS process is maintained well.

Let us now consider more carefully the development of the saturation stage. When an electron obtains enough energy to escape the unstable plasma wave, it travels along the target surface. The fields of the nodes may slightly differ from each other (see Figs. <ref>b and <ref>a), hence, escaping one node electron may still be trapped by another one with a higher field. This stochastic process produces electron flow from smaller to bigger wave buckets, and lead to the nodes merging. In the spectral space, this process corresponds to a shift of the developed unstable plasma mode to the longer wavelengths. This spectral down-shift driven by thermalization is known in physics of turbulent plasma instabilities <cit.>, and it can also be observed in our numerical simulations. In Fig. <ref>, we show the proton density maps and their Fourier images at two different times of 75 ω_0^-1 and 100 ω_0^-1. One may clearly see, that during the saturation phase, the sizes and the spacings between the density perturbations increase. In the spectral space the long-wave fluctuations develop within the initial ring structure, with the unstable k value falling from 5 ω_0/c to about 1.5 ω_0/c. In this simulation, the maximum acceleration time t_acc is around 100 ω_0^-1, corresponding to about 40 fs for a 800nm laser.



Besides the typical simulation shown above, we have also systematically scanned over a wide range of laser plasma parameters to verify the validity of the theoretical model. Moreover, since the derivation makes assumption of a homogeneous radiation pressure, it is also important to see how instability development is affected by a more realistic laser profile. For this, we have also performed a 3D PIC simulation with a transverse Gaussian profile, that exhibits similar ring structures in the central region within the laser focal spot. For the case of a tightly focused laser, we have observed additional electron heating arousing at the target periphery, which may spread and destroy the target even before the instability fully develops. More details of these additional simulations are presented in the Supplementary Material <cit.>.



§ DISCUSSION


As we have shown, the plasma light-sail instability triggers the target heating and indeed presents a major obstacle for efficient ion acceleration. It also follows, that this instability always develop in the laser-driven thin foils, and cannot be avoided except by preventing the coupling between the oscillating electrons and ions. One such possibility is to physically separate electrons and ions using, for example, a tweezer effect of two colliding lasers with the different frequencies <cit.>. In this case, however, acceleration is no longer provided by the radiation pressure in a normal manner. Considering only LS regime, we will look to minimize the instability growth rate, and for this it is crucial to identify the proper interaction parameters.

Let us consider the acceleration efficiency as a function of the target design, i.e. the ion compound, initial density or temperature. As mentioned before, in the optimal LS process the foil thickness, density and laser amplitude are coupled by the balance condition 2a_0≃ l_0 ω_0 n_0/c n_c <cit.>. Therefore, for the given laser characteristics an optimum can be found by varying a single parameter – either foil thickness or density. From Eq. (<ref>) follows, that instability growth correlates with the plasma density, which in its turn correlates with the initial target density. This means, that for the lower target density the acceleration time can be extended resulting in the higher proton energies. This statement can be verified numerically, and here we perform a series of 2D PIC simulations for a_0=5 and varying foil density, while its thickness is adjusted with respect to the balance condition. In Fig. <ref> we show the times before foil breaking and ion energy corresponding to the peak in their quasi-monoenergetic spectra. This result clearly confirms our hypothesis, as when the target density decreases from 100 n_c to 3 n_c, maximum acceleration time is extended from 70 ω_0^-1 to around 250 ω_0^-1, and the final proton energy increases from 12 MeV approaching 90 MeV. We should note here, that in 2D case, instability has less degrees of freedom and develops slightly slower, while the behavior with respect to the parameters remains equivalent. 



One more option to inhibit the instability, is to consider an initial electron temperature T_e high enough to provide a thermal smoothing of the perturbations <cit.>. The smoothing scale is defined by the electron Debye length λ_D=√(T_e/4π n_0e^2), and to affect the instability it needs to be of the order of the instability scale length, i.e. k_mλ_D∼ 1.  In the case of an intense laser, this condition defines temperature T_e∼γ_0, which can be as high as several MeV and leads to a rapid plasma expansion. An additional numerical study of the electron temperature effect have demonstrated, that in all cases it leads to an earlier target breaking than in a case of the initially cold plasma.

Among other possibilities to improve the LS performance, it is worth to mention the multi-species targets, where heavy ions provide a co-moving electrostatic field for protons even after the plasma is disrupted by the instability <cit.>. Another proposed alternative is to use an extremely intense single-cycle laser to provide acceleration on a time before the instability develops <cit.>.

In conclusion, we have presented and explored a three-dimensional theory of the transverse instability in the laser-driven plasma light sail. It is shown, that both strong electron-ion coupling effect and nonuniform laser pressure contribute to the instability at early linear stage, while the former dominates. We have identified the unstable mode structure and estimated its growth rate, which was found in good agreement with 3D PIC simulations in a wide range of parameters. The nonlinear instability stage was studied numerically, and it was shown that it saturates with the wave-breaking effect, and leads to a strong plasma heating. The following target breaking is accompanied by merging of the smaller density structures into the larger ones. Possible solutions to suppress the instability are discussed and summarized.

This work was supported by NSFC Grant No. 11425521, No. 11535006, No. 11475101 and No. 11775125, the Thousand Young Talents Program, Gerry Schwartz and Heather Reisman, Israel Science Foundation, and VATAT supports. Simulations were performed on Edison and Cori clusters at NERSC.


