## 1. Introduction

Collimated relativistic electron beams (REBs) with high energy are widely used in astrophysics, inertial confinement fusion (Craxton *et al.* Reference Craxton, Anderson, Boehly and Zuegel2015; Yang *et al.* Reference Yang, Zhuo, Xu and Ma2016), plasma based accelerators (Geddes *et al.* Reference Geddes, Toth, van Tilborg and Leemans2004; Esarey, Schroeder & Leemans Reference Esarey, Schroeder and Leemans2009) and new diagnostic technologies (Barbieri *et al.* Reference Barbieri, Hansen-Ilzhöfer, Ilzhöfer and Holzwarth2000). In these areas, many potential application scenarios such as space active experiments (Marshall *et al.* Reference Marshall, Nicolls, Sanchez, Lehtinen and Neilson2014; Borovsky & Delzanno Reference Borovsky and Delzanno2019; Xue *et al.* Reference Xue, Hao, Zhao and Dong2020) in the middle atmosphere require long-distance transport of electron beams. Moreover, the plasmas produced by the gas ionization and the plasmas of the middle atmosphere have clear effects on REB transport, such as self-modulation (Pukhov *et al.* Reference Pukhov, Kumar, Tückmantel and Shvets2011; Schroeder *et al.* Reference Schroeder, Benedetti, Esarey, Grüner and Leemans2011). Thus, the study of the stability of long-distance transport of REB in plasmas is of great significance to explore the applications of electron beams.

The transport behaviour of charged particle beams in plasmas is mainly affected by the wakefield (Lu *et al.* Reference Lu, Huang, Zhou, Mori and Katsouleas2006; Esarey *et al.* Reference Esarey, Schroeder and Leemans2009; Gonsalves *et al.* Reference Gonsalves, Nakamura, Daniels and Leemans2019; Batsch *et al.* Reference Batsch, Muggli and Agnello2021) driven by these beams. Lotov found that the electron bunch can produce a well-controlled wakefield that is suitable for particle acceleration, while the shaped long bunch is unusable since it is rapidly destroyed by the transverse two-stream instability (Lotov Reference Lotov1998). Pukhov *et al.* found that column-shaped proton beams with TeV energy and steep density gradient can effectively excite a GV m^{−1} wakefield and achieve electron acceleration (Pukhov *et al.* Reference Pukhov, Kumar, Tückmantel and Shvets2011). The results also show that a proton beam with truncated Gaussian profile can excite the wakefield with higher acceleration gradient compared with that of Gaussian profiles (Kumar, Pukhov & Lotov Reference Kumar, Pukhov and Lotov2010; Hui *et al.* Reference Hui, Hu, Zhao, Cheng, Mei and Wang2019). Schroeder *et al.* studied the influence of instability competition (Schroeder *et al.* Reference Schroeder, Benedetti, Esarey, Grüner and Leemans2012) and phase velocity dephasing (Schroeder *et al.* Reference Schroeder, Benedetti, Esarey, Grüner and Leemans2011) during beam transport from the aspect of a self-modulating beam-driven plasmas wave. They show that the development of the hosing instability and phase velocity dephasing in homogeneous plasmas will affect the acceleration of the wakefield and lead to beam defocusing. Fang *et al*. studied the effect of seeding excitation on wakefield acceleration (Fang *et al.* Reference Fang, Yakimenko, Babzien and Muggli2014*b*), and the envelope of the seeds has an obvious effect on the excited wakefield (Fang *et al.* Reference Fang, Vieira, Amorim, Mori and Muggli2014*a*) in electron acceleration experiments. Zhou *et al.* studied the effect of proton beam density distribution on self-modulation instability in the process of a proton beam excited plasma wave. It is revealed that a long pulse proton beam with a smooth density distribution is beneficial for long-distance transport (Zhou *et al.* Reference Zhou, Ma, Tian and Shao2010). Whittum and Vieira *et al.* studied the relationship of the hosing instability and the self-modulation instability (Whittum Reference Whittum1993; Vieira, Mori & Muggli Reference Vieira, Mori and Muggli2014), and found that the hosing instability of long particle beams can be suppressed and stabilized by full self-modulation. Neubert *et al.* studied electron beam transport in the neutral atmosphere by Monte Carlo simulations, and showed that the electron behaviour is affected by electron--neutral collisions and the geomagnetic field (Neubert *et al.* Reference Neubert, Gilchrist, Wilderman, Habash and Wang1996). Generally, most of the related studies mainly discussed the stable accelerating field generated by the REBs in the areas of plasma accelerators and long-distance transport in the neutral atmosphere. However, long-distance REB transport in plasmas is rarely investigated and the influence of the instabilities on the REB transport is still not well understood.

In this paper, REB transport in plasmas of hundreds of metres is studied by using particle-in-cell simulation. It is shown that the stability of REBs during the transport in plasmas is mainly affected by the self-modulation instability (SMI), hosing instability (HI) and filamentation instability (FI). All these instabilities are influenced by the beam density ratio, the beam energy, the beam initial envelopes and the density profiles. By choosing appropriate parameters, we avoid the development of the instabilities and ensure the beam propagates stably for hundreds of metres, and the beam transport energy loss rate is reduced from 40 % to 10 % while controlling the beam spot radius sufficiently to be small in our simulation conditions. Note that the stable transport of REB here means that the remaining energy is more than 50 % during REB transport over the considered distance.

## 2. Theoretical analysis

When REBs propagate through plasmas, various instabilities can be excited due to the convection and resonance between electron beams and background electrons. These instabilities can be divided into transverse instabilities and longitudinal instabilities. The transverse instabilities mainly include SMI, FI and HI, while the longitudinal instabilities mainly include the two-stream instability. The growth rate of longitudinal instabilities is much smaller than that of the transverse instabilities because of the anisotropy of the relativistic beam response to the longitudinal force and transverse force. Therefore, the transverse instabilities are the dominant factors affecting the beam transport and will be discussed as follows. It should be noted that self-focusing by the beam self-generated fields may also focus the beam, however, by comparing the simulation results, we believe that self-modulation is the main factor due to the fact that the field generated by the self-modulation is a few times higher than that of the self-focusing field and the beam will defocus without the self-modulation field. Generally, the SMI can focus the beam in a plasma wavelength because of its strong transverse wakefield, so utilizing the focusing field of the SMI while suppressing the development of the HI and the FI can lead to a stable transport of REBs.

The main factor of beam energy loss is the defocusing field of the SMI. The instability growth rate can be obtained from the beam-envelope approximation. Beam propagates along the *X* direction. The beam density profile can be described as $\rho (r,\xi ) = {\rho _0}\psi (r)f(\xi )$ where $\xi = {\beta _0}ct - x$ is the commoving ordinate, *c* is the light velocity, *t* is time, $\psi$ is the radial profile function, *f* is the longitudinal profile function, ${\rho _0}$ is the initial beam density and ${\beta _0}$ is normalized longitudinal velocity. Under the assumption of uniform plasmas and motionless ions, the longitudinal wakefield and transverse wakefield formed by the beam can be expressed as (Keinigs & Jones Reference Keinigs and Jones1987; Kumar *et al.* Reference Kumar, Pukhov and Lotov2010)

where ${I_{0(1,2)}},{K_{0(1,2)}}$ are the modified Bessel functions, *B* is the magnetic field, *E* is the longitudinal wakefield, *W* is the transverse wakefield, ${k_p} = {\omega _p}/c$ is the wavenumber of background plasmas and ${\omega _p}$ is the plasma frequency, ${r_ < } = \min (r,r^{\prime})$, ${r_ > } = \max (r,r^{\prime})$. For an electron beam with a steep density profile, the initial radial envelope of the beam can be represented by the Heaviside step function as $\psi (r) = \varTheta ({r_b} - r)$. As the beam radial divergence is mainly affected by the transverse field, the equation for the beam envelope with an arbitrary profile $f(\xi )$ in $\xi$ (after being normalized by ${\tilde{r}_b} = {r_b}/{r_{b0}},\; \; \tilde{\tau } = {\omega _{\beta 0}}\tau ,\; \; \tilde{\xi } = {k_p}\xi$) is written as (Reiser Reference Reiser2008; Pukhov *et al.* Reference Pukhov, Kumar, Tückmantel and Shvets2011)

where $\omega _{\beta 0}^2 = \omega _b^2/2{\gamma _0}$ and $\tilde{\epsilon }_n^2 = {\sigma _\theta }/({k_{\beta 0}}{r_{b0}})$ with ${\sigma _\theta }$, ${r_{b0}}$ being the bunch angular divergence and initial radius. The second term on the left is the effect of divergence on the beam envelope, and the term on the right is the effect of the field on the beam envelope.

Assuming that the small perturbation $\delta r_b$ can be written as $\delta {r_b} = \delta {\hat{r}_b}exp ( - \textrm{i}\delta \omega \tilde{\tau } + \textrm{i}k\tilde{\xi })$, by introducing $\delta \omega ^{\prime} = \delta \omega - \nu k,\nu = \tilde{\xi }/\tilde{\tau }$, the dispersion equation $\partial D(\delta \omega ^{\prime},k)/\partial k = 0$, $D(\delta \omega ^{\prime},k) = 0$, can be written as

To simplify the solution, we assume that $\tilde{\epsilon }_n^2$ is small enough to be omitted and $k = 1 + \delta k$, $\delta k\ll 1$, thus, in dimensional variables, the instability growth rate can be write as

where $n_b$ is the beam density and $n_e$ is the background plasma density.

The growth rate ratio of the HI to the SMI is ${\varGamma _h}/{\varGamma _{\textrm{SM}}} = {(\mu /2\nu )^{1/3}}\sim 1$, where $\mu = 2{I_1}({r_0}){K_1}({r_0}),\;\nu = 4{I_2}({r_0}){K_2}({r_0})$ (Schroeder *et al.* Reference Schroeder, Benedetti, Esarey, Grüner and Leemans2012). It should be noted that the ${\varGamma _{\textrm{SM}}}$ from Kumar is obtained from the beam-envelope function and the ${\varGamma _{\textrm{SM}}}$ from Schroeder is obtained from the beam centroid evolution. The two growth rates have a similar form, but the influence of ${r_0}$ was omitted during our above solution, so we use the results from Schroeder as a supplement. The HI growth rate is comparable to the SMI growth rate and the dominant relationship of these instabilities may transform with different instability seeds, but generally the SMI is more likely to dominant as the radius ${r_0}$ increases according to the growth rate ratio of the HI to the SMI. The FI is also related to the radius ${r_0}$, and decreasing the radius can reduce beam loss from filamentation (Bret, Firpo & Deutsch Reference Bret, Firpo and Deutsch2005). Therefore, without additional instability seeds, the SMI and HI can be generated by beam misalignments, beam transport errors or any asymmetries in the beam creation and delivery, and by decreasing the electron density ratio ${n_b}/{n_e}$, increasing the relativistic Lorentz factor ${\gamma _b}$ and adjusting the initial beam radius ${r_0}$, the SMI will dominate the main instability and inhibit the HI and the FI, thus realizing stable transport of the REBs.

## 3. Simulation results and discussions

### 3.1. Transport of REB in plasmas

The simulations are performed by using the full three-dimensional particle-in-cell simulation code LAPINE (Xu *et al.* Reference Xu, Yu, Yu, Wong and Sheng2012; Yang *et al.* Reference Yang, Dieckmann, Sarri and Borghesi2012). According to the relevant data in the atmosphere at an altitude of 60 km in the IRI (International Reference Ionosphere) ionospheric model, the background electron density is set to ${n_e} = {10^{19}}\,{\textrm{m}^{ - 3}}$ and the initial temperature is 0.1 eV in this work. We select a cylindrical beam as the standard set-up to discuss. In the simulation, the initial beam is a cylindrical beam with ${n_b} = {10^{15}}\,{\textrm{m}^{ - 3}}$, ${\gamma _b} = 40$, $L = 2\,\textrm{cm}$ and ${r_0} = 0.5\,\textrm{cm}$ which has a steep density profile, and the beam propagates along the *X* direction. The density profile of the cylindrical beam is $\rho (r,\xi ) = {\rho _0}\varTheta ({r_b} - r)\varTheta (\xi )$. The simulation box is $60 \times 80 \times 80\,\textrm{mm}$ with the cell size being $0.2 \times 0.2 \times 0.2\,\textrm{mm}$. The current injected is ${\sim}10$ A at the initial time. In order to decrease numerical heating, a third-order interpolation scheme is used to evaluate the particles and a fourth-order finite difference scheme is used to evaluate the field. The beam electron macroparticles per cell is set to 8 and the background plasma (including N^{+} and electrons) macroparticles per cell is set to 5. For both the transverse and longitudinal boundaries, simple outflow boundary conditions are used for the fields and particles. A moving window is used in the simulation and the velocity of the moving window is ${v_x} = 2.997 \times {10^8}\,\textrm{m}\,{\textrm{s}^{ - 1}}$. Since the impact of the collision and the initial divergence of the REBs on the beam radius during the transport are less than 4 %, we do not consider these factors in this work for simplicity.

Figure 1 shows the evolution of the electron beam transport behaviour in a plasma background with time. It can be seen that, after the REBs enter the plasma, the plasma wave is excited first, and the plasma wavelength is given by ${\lambda _p} = 2\mathrm{\pi }c/{\omega _p},\;\omega _p^2 = 4\mathrm{\pi }{n_e}{e^2}/{m_e}$ (i.e. ${\lambda _p} \approx 1\,\textrm{cm}$ when ${n_e} = {10^{19}}\,{\textrm{m}^{ - 3}}$). The plasma wave and beam interaction excite the SMI, causing changes in the beam envelope. Under the influence of the fields, the beam head is not significantly affected. The beam neck is compressed while the waist is split into filaments (see figure 4 for details) and defocuses gradually, therefore, the main loss of the beam happens in the waist part because of the defocusing field. The head of the next cycle appears in the tail of beam. Significant beam centroid displacement is observed during the transport, leading to the HI during the transport.

To analyse the physical process of transport, a two-dimensional illustration of the wakefield at the time $t = 90\,\textrm{ns}$ is presented in figure 2. The red solid line in figure 2(*a*) represents the change of ${n_e}$ at the position of the dashed line in figure 2(*d*). The other lines are the same. It can be seen that the fields only take large values in the central axis and decay significantly along the radial direction. The fields in this case are ${W_{ \bot \max }} \approx {E_{x\max }} \approx 2 \times {10^5}\,\textrm{V}\,{\textrm{m}^{ - 1}}$. Affected by the SMI, the beam-driven transverse wakefield ${W_ \bot }$ and the beam-driven longitudinal electric field ${E_x}$ will cause the beam to focus (defocus) and accelerate (decelerate), respectively. From the density distribution, the beam is periodically focused and defocused by the influence of the ${W_ \bot }$ field. The change of the beam envelope is highly consistent with the field change.

Since there are no distinct acceleration and deceleration peaks affected by the ${E_x}$ field, for the energy spectrum in figure 3(*a*), the acceleration of the beam is not obvious. There is significant velocity broadening. With the increase of the transport distance, the energy of the electron beam is transferred to the fields, which leads to the deceleration of the whole beam. Figure 3(*b*) shows the evolution of the beam shape at different times, the beam shapes are made from the two-dimensional density contour of ${n_e} = 3 \times {10^{13}}\,{\textrm{m}^{ - 3}}$ at $Z = 0$. It is shown that, with the increase of transport distance, the beam segments into two pieces and the beam centroid moves back gradually.

### 3.2. The effects of beam parameters on the REB transport

According to the instability growth rate (2.6) and the growth rate ratio of the HI to the SMI, one can see that the beam transport is obviously affected by the beam-plasma electron density ratio ${n_b}/{n_e}$, the relativistic Lorentz factor ${\gamma _b}$ and the initial beam radius ${r_0}$. In addition, the beam loss rate is closely related to the electron beam length *L* according to previous simulations. Thus, varying these parameters may improve the beam stability and reduce the beam loss during transport.

In order to discuss the influence of the beam envelope, the beam length *L* and beam radius ${r_0}$ are analysed firstly. Figure 4 shows the three-dimensional density distribution of the beam with $L = 2\,\textrm{cm},\;{r_0} = 1\,\textrm{cm}$ (*a*) and $L = 2\,\textrm{cm},\;{r_0} = 0.5\,\textrm{cm}$ (*b*) at *t* = 50 ns, respectively. It can be seen that, at *t* = 50 ns, the beam current excites three obvious instabilities during transport, namely, the transverse SMI, the HI and the FI in figure 4(*a*). The SMI leads to the beam segmenting into pieces; the HI leads to the destruction of beam symmetry of the beam density distribution along the *X* axis and the FI leads to the beam splitting into filaments at the head and waist of the beam (Vieira *et al.* Reference Vieira, Fang, Mori, Silva and Muggli2012). Compared with the beam with smaller radius, the beam head splits into filaments and the middle part of the beam becomes slightly twisted. As the beam transports further, the SMI increases, while the HI and FI gradually disappear.

Figure 5 shows the evolution of the electron beam transport behaviour in the plasma background with time. Here, the initial beam is a cylindrical beam with ${n_b} = {10^{15}}\,{\textrm{m}^{ - 3}}$, ${\gamma _b} = 40$, $L = 1\,\textrm{cm}$ and ${r_0} = 0.5\,\textrm{cm}$ which has a steep density profile. It can be found that the beam with $L\sim {\lambda _p}$ can also excite the SMI, and the beam loss of the short beam is much smaller than the beam with $L \ge {\lambda _p}$ because less of the beam is located in the defocusing field. Since the first period of the plasma wave is always larger than the rest, most of the beam is located in the focusing field, which occupies approximately half the plasma wave. At the same time, because of the small beam initial radius ${r_0}$, the FI of the beam is effectively controlled, and the beam transport is relatively stable.

The effects of the beam density and relativistic Lorentz factor on the beam transport are presented in figure 6. It can be found that, as the beam density ${n_b}$ increases, the amplitude of the beam in the defocusing field expands fast and the instability develops rapidly, which is attributed to the fact that the transverse wakefield ${W_ \bot }$ is enhanced. The highest beam density of ${n_{b\max }} = {10^{19}}\,{\textrm{m}^{ - 3}}$ is considered here. For such densities, the transverse field ${W_ \bot }$ will strongly disturb the background electrons and the HI will cause the beam symmetry to be destroyed (Krall & Joyce Reference Krall and Joyce1995). Due to the focusing effect of the ${W_ \bot }$ field, the beam will also collapse rapidly as part of the beam density in the focusing field is close to that of the background electrons. The beam collimation gets better and the instability develops more slowly with the increase of the electron Lorentz factor. One can see that the simulations share the same trend with the theoretical results, and are in good agreement in the low beam density and low Lorentz factor cases. The influence of the longitudinal deceleration that causes the ${\mathrm{\gamma }_b}$ drop on the beam transport gets more obvious as ${n_b}/{n_e}$ and ${\mathrm{\gamma }_b}$ become larger, leading to the deviation between the simulation and theoretical results gradually growing.

### 3.3. Improvement of the beam transport by control the initial beam profile

Since the beam envelope is mainly affected by the wakefields, a smaller wakefield amplitude would be beneficial for the propagation of the beam, which can be achieved with a smooth density profile. Here, we will control the beam transport by modulating the initial beam profile to a profile with a smooth density gradient such as a Gaussian profile and other smooth density profile.

Figure 7 shows the evolution of the electron beam transport behaviour in a plasma background with time. The initial beam is a Gaussian beam with ${n_{b\max }} = {10^{15}}\,{\textrm{m}^{ - 3}}$, ${\gamma _b} = 40$, $L = 2\,\textrm{cm}$ and ${r_0} = 0.5\,\textrm{cm}$. The density profile is $\rho (r,\xi ) = {\rho _0}{\rm exp} ( - {\xi ^2}/\sigma _\xi ^2){\rm exp} ( - {r^2}/\sigma _r^2)$. It can be seen that, unlike the cylindrical beam, the Gaussian beam gradually becomes a conical structure during the transport process, and can maintain its shape for a long distance without obvious HI, FI and even SMI. The segmentation is not obvious under the Gaussian profile and the beam defocusing only happens in the tail. The beam centroid displacement is not obvious in this case, which means the HI may be suppressed.

To analyse the physical process compared with the cylindrical beam, a two-dimensional illustration of the wakefield at the time $t = 90\,\textrm{ns}$ is presented in figure 8. The cylindrical beam and the Gaussian beam have a similar field wavelength under the same plasma conditions. The Gaussian profile obviously has two different features compared with those shown in figure 2. Firstly, the field strength of Gaussian beams is only 20 % of the cylindrical beam, which is only ${W_{ \bot \max }} \approx {E_{x\max }} \approx 4 \times {10^4}\,\textrm{V}\,{\textrm{m}^{ - 1}}$. Secondly, the field strength of the beam head is similar to the other parts and has no significant phase lag, which may cause the beam change to a cone. Furthermore, transport of a density-uniform ellipsoidal beam with a steep density profile is also studied (not shown for brevity). The ellipsoidal beam with ${r_{\textrm{long}}} = 1\,\textrm{cm},\;{r_{\textrm{short}}} = 0.5\,\textrm{cm}$ also has an obvious segment like the cylindrical beam during its transport. That is, the segmentation can be inhibited by the smooth density profile. It is worth mentioning that REB transport with initial seeds of the HI is also performed, where a beam density profile given by $\rho (r,\xi ) = {\rho _0}{\rm exp} ( - {\xi ^2}/\sigma _\xi ^2){\rm exp} ( - ({x^2} + {(y - {\delta _{\textrm{HI}}}\xi )^2})/\sigma _r^2)$ (Vieira *et al.* Reference Vieira, Mori and Muggli2014) is applied. By comparing the results of ${\delta _{\textrm{HI}}} = 0$, ${\delta _{\textrm{HI}}} = 0.001$ and ${\delta _{\textrm{HI}}} = 0.01$, we find no obvious difference between these cases, which is different to that of Vieira's results (not shown for brevity). Here, the HI is suppressed by the saturated SMI.

The beam spectrum is relatively narrow and the beam energy is more concentrated in figure 9(*a*) compared with those in figure 3, which may be due to the small wakefield ${E_x}$ and ${W_ \bot }$ excited by the Gaussian beam. Figure 9(*b*) shows the evolution of the beam shapes at different times, the beam shapes are extracted from the two-dimensional density contour of ${n_e} = 3 \times {10^{13}}\,{\textrm{m}^{ - 3}}$ at $Z = 0$. It is shown that, with the increase of the transport distance, the beam head becomes to a cone-like distribution and the beam tail defocuses, while the beam centroid has no significant change. Note that these contours are relatively larger than the contours of the cylindrical beam due to the Gaussian distribution of REB.

The beam energy loss of different initial parameters are displayed in figure 10(*a*), it is shown that the energy loss turning point of a Gaussian beam is behind that of a cylindrical beam under the same conditions $(L = 2\,\textrm{cm,}\;{r_0} = 0.5\,\textrm{cm})$. The beam energy loss can also decrease significantly as the beam length decreases $(L = 1\,\textrm{cm,}\;{r_0} = 0.5\,\textrm{cm})$, indicating that the short Gaussian beam has more advantages in beam stable propagation. The energy loss caused by the increase of beam radius is relatively smooth, and the increase of initial energy can further reduce the energy loss rate. The remaining beam energy and charge at $t = 200\,\textrm{ns}$ are shown in figure 10(*b*). The energy and charge are calculated within areas of different radii around the beam propagation axis. It can be seen that more than 70 % of the energy and the charge are concentrated within the circle of $R = 1\,\textrm{cm}$ around the centre. The remaining energy and charge have a similar trend. With these improvements, we eventually decrease the beam transport energy loss rate from approximately 40 % to 10 % with a beam spot radius that is sufficiently small. The results meet our expected requirements for stable transport.

## 4. Conclusion

The long-distance stable transport of a REB in plasma is studied by theoretical analysis and numerical simulations. To ensure stable transport, generating a stable SMI while suppressing the HI and the FI are of great importance. There are five main factors that can affect the beam stable transport, i.e. the electron density ratio ${n_b}/{n_e}$, the relativistic Lorentz factor ${\gamma _b}$, the beam envelope $(L,{r_0})$ and the beam profile. The decrease of the electron density ratio ${n_b}/{n_e}$ and increase of the relativistic Lorentz factor can decrease the growth rate of SMI; the adjustment of the initial radius ${r_0}$ can inhibit the development of FI and HI; the decrease of beam length *L* can reduce the loss of current and a beam with a Gaussian profile may reduce the strength of the defocusing field and further ensure beam transport. Therefore, relativistic bunches having a smooth density profile and a length of several plasma wave periods can suppress the beam-plasma instabilities and propagate in plasmas over long distances with small energy losses. The results provide a good reference for REB applications that require long-distance and stable transport.

## Acknowledgements

*Editor Luís O. Silva thanks the referees for their advice in evaluating this article*.

## Declaration of interests

The authors report no conflict of interest.

## Funding

This work was supported by the National Natural Science Foundation of China (Grant Nos. 11775305, 11975308, 12005297 and 12075317), the Strategic Priority Research Program of Chinese Academy of Science (Grant No. XDA25050200), the fund of Science Challenge Project (No. TZ2018001) and the State Key Laboratory of Laser Interaction with Matter (No. SKLLIM1908). X.H.Y. also acknowledges the financial support from Fund for NUDT Young Innovator Awards (No. 20180104).

## Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.