2. Theory

Before the theory is introduced, two aspects should be emphasized in advance. The first is that when a charged particle beam is injected into a resistive plasma, we assume the beam is quickly neutralized by plasma electrons. This assumption holds as the beam duration time is much longer than the electron plasma period $2{\rm \pi} /\omega _{\textrm {pe}}$ ($\omega _{\textrm {pe}}$ is the electron plasma frequency) and the beam radius is much larger than the plasma electron skin depth $c/\omega _{\textrm {pe}}$ ($c$ is the light speed) (Kaganovich et al. Reference Kaganovich, Shvets, Startsev and Davidson2001,  Reference Kaganovich, Startsev, Sefkow and Davidson2008; Berdanier, Roy & Kaganovich Reference Berdanier, Roy and Kaganovich2015). The second aspect is that the focusing we are describing here is the transverse compression of the entire beam. As shown in figure 1, after entering the plasma, the beam is in a state of unbalanced force in the transverse direction, which causes its radius to gradually become smaller and at the same time, its density increases, which lead to the focusing. To intuitively explain this phenomenon and determine the focusing condition, we need to figure out the major forces acting on the beam, so we might as well start with a force analysis of the beam. When a charged particle beam propagates in a uniform background plasma, its motion equation can be written as

(2.1)\begin{equation} n_b\frac{\textrm{d}\boldsymbol{p}_b}{\textrm{d}t} ={-}\boldsymbol{\nabla}p_{\textrm{th}}+\boldsymbol{J}_b\times\boldsymbol{B}+n_bq\boldsymbol{E}. \end{equation}

On the left-hand side of the equation, $n_b$ is the density of the beam and $\boldsymbol {p}_b$ is the momentum of the beam particle; on the right-hand side, $p_{\textrm {th}}$ is the beam thermal pressure and can be estimated from the equation of state $p_{\textrm {th}}=n_bkT_b$, where $k$ is Boltzmann's constant and $T_b$ is the beam temperature, $\boldsymbol {J}_b$ is the beam current density, $q$ is the effective charge of a beam particle, $\boldsymbol {B}$ is the magnetic field and $\boldsymbol {E}$ is the electric field. According to the basic Ohm's law, the electric field $\mathrm {E}$ can be expressed as

(2.2)\begin{equation} \boldsymbol{E}=\eta\boldsymbol{J}_e-\boldsymbol{v_e}\times\boldsymbol{B}-\frac{1}{en_e}\boldsymbol{\nabla} p_e, \end{equation}

where $\eta$ is the resistivity of the background plasma, considered to obey Spitzer formulation, $\eta \propto T^{-3/2}$, $\boldsymbol {J}_e$ is the return current density of plasma electrons, $\boldsymbol {v}_e$ is the flow velocity of plasma electrons, $e$ is the elementary charge. Additionally, $n_e$ and $p_e$ are the density and the thermal pressure of plasma electrons, respectively. In Cartesian coordinates $(x, y, z)$, considering the beam moving in the $z$ direction, when the transverse density distribution is assumed to be Gaussian profile, we can express its current density as $\boldsymbol {J}_b=n_bqv_b\boldsymbol {e}_z=n_{b0}qv_b\exp (-(x^{2}+y^{2})/R^{2})\boldsymbol {e}_z$, where $n_{b0}$ is the central density of the beam, $v_b$ is the beam velocity and $R$ is the beam radius. Here, to illustrate our study method, we only take the $y$ direction as an example, and therefore, it is not necessary to consider the density distribution in the $x$ direction. For simplicity, we assume plasma ions to be at rest. Then, with the application of the beam charge and current neutralization, we have $n_e\approx n_p+({q}/{e})n_b$ and $\boldsymbol {J}_e\approx -\boldsymbol {J}_b$, where $n_p$ is the initial density of plasma electrons. We assume the charge density of the beam is much smaller than that of the plasma. Therefore, we can obtain $\boldsymbol {v}_e\ll \boldsymbol {v}_b$. The thermal pressure $p_e$ can also be estimated from the equation of state $p_e=n_ekT_0$, where $T_0$ is the plasma temperature. Then, the electric field can be written as

(2.3)\begin{equation} \boldsymbol{E}\approx{-}\eta\boldsymbol{J}_b-\boldsymbol{v}_e\times\boldsymbol{B}-\frac{qkT_0}{e^{2}n_e}\boldsymbol{\nabla} n_b. \end{equation}

For the beam, the forces of the electric field in the transverse direction, i.e. the second and the third terms on the right-hand side of (2.3), can be considered as small quantities compared with the Lorentz force $\boldsymbol {J}_b\times \boldsymbol {B}$ and the beam thermal pressure $\boldsymbol {\nabla } p_{\textrm {th}}$. This is because of $|n_bq\boldsymbol {v}_e\times \boldsymbol {B}|/|\boldsymbol {J}_b\times \boldsymbol {B}|\approx v_e/v_b\ll 1$ and $|n_bq^{2}kT_0\boldsymbol {\nabla } n_b/(e^{2}n_e)|/|\boldsymbol {\nabla } p_{\textrm {th}}|\approx |n_bq^{2}T_0|/|n_ee^{2}T_b|\ll 1$, where we have assumed the beam temperature is not much lower than the plasma temperature, $T_b\lessapprox T_0$. Therefore, we might approximately express the electric field as $\boldsymbol {E}\approx -\eta \boldsymbol {J}_b$. Hence, the motion of the beam in the $y$ direction is constrained by the thermal pressure $p_{\textrm {th}}$ and the Lorentz force $J_bB_x$. The motion equation in the $y$ direction further becomes

(2.4)\begin{equation} n_b\frac{\textrm{d}p_{by}}{\textrm{d}t}={-}\frac{\textrm{d}p_{\textrm{th}}}{\textrm{d}y}+J_bB_x, \end{equation}

where $p_{by}$ is the component of the beam momentum in the $y$ direction. The magnetic field $B_x$ can be derived from Faraday's law:

(2.5)\begin{equation} \frac{\partial \boldsymbol{B}}{\partial t}={-}\boldsymbol{\nabla}\times\boldsymbol{E}=\boldsymbol{\nabla}\times\left(\eta\boldsymbol{J}_b\right). \end{equation}

This model for solving the electric field and the magnetic field has been widely used by many researchers (Davies et al. Reference Davies, Bell, Haines and Guérin1997; Bell, Davies & Guerin Reference Bell, Davies and Guerin1998; Davies Reference Davies2003; Davies, Green & Norreys Reference Davies, Green and Norreys2006; Robinson et al. Reference Robinson, Key and Tabak2012; Norreys et al. Reference Norreys, Batani, Baton, Beg, Kodama, Nilson, Patel, Pérez, Santos and Scott2014; Kim et al. Reference Kim, McGuffey, Qiao, Wei, Grabowski and Beg2016; Curcio & Volpe Reference Curcio and Volpe2019). For example, to simplify calculations and obtain a specific analytical form of the fields, Davies considered a rigid beam model (Davies Reference Davies2003). Using this model, the plasma temperature is increased by Ohmic heating, $\partial T/\partial t=2\eta J_b^{2}/(3kn_e)$ and the resistivity also changes accordingly, $\eta =\eta _0(T/T_0)^{-3/2}$, where we have applied the Spitzer resistivity formulation. Then combining these two relations with (2.5) and doing some integral operations, we can get the magnetic field. Making use of Davies’ conclusion (Davies Reference Davies2003), we can express the magnetic field as

(2.6)\begin{equation} B_x=\frac{\textrm{d}J_b}{\textrm{d}y}\frac{3n_ekT_0}{2J_b^{2}}\left[1-\frac{1}{5}\left(1+\varOmega\right)^{2/5}-\frac{4}{5}\left(1+\varOmega\right)^{{-}3/5}\right], \end{equation}

where

(2.7)\begin{equation} \varOmega=\frac{5}{3}\frac{\eta_0J_b^{2}\tau}{n_ekT_0}, \end{equation}

and $\tau =t-z/v_b$ is the time variable. Here, $\varOmega$ means the ratio of Ohmic heating energy to plasma electrons thermal energy and $\varOmega _{\mathrm {max}}$ is located at $y=0$. For $\varOmega _{\mathrm {max}}\ll 1$, the plasma temperature is slightly increased, which corresponds to weak heating, and for $\varOmega _{\mathrm {max}}\gg 1$, the plasma temperature is significantly increased, which corresponds to strong heating. A schematic diagram of the magnetic field distributions in the cases of weak heating and strong heating is shown in figure 2. For weak heating, the magnetic field has a pure focusing effect; and for strong heating, the case is different: the magnetic field near the centre changes its sign and will push the beam outward, while there exists a magnetic field away from the centre that will push the beam inward. This may lead to a result where the beam becomes hollowed near the centre and focused away from the centre. To determine the position where the magnetic field changes its sign, we can solve the following equation:

(2.8)\begin{equation} 1-\frac{1}{5}\left(1+\varOmega\right)^{2/5}-\frac{4}{5}\left(1+\varOmega\right)^{{-}3/5}=0, \end{equation}

of which the numerical solution is $\varOmega \approx 44$, and for strong heating, the beam focusing can only occur in the range of $0<\varOmega <44$. The beam focusing depends on the relation between the beam thermal pressure $-({\textrm {d}p_{\textrm {th}}}/{\textrm {d}y})$ and the Lorentz force $J_bB_x$. They have the forms as follows:

(2.9)\begin{equation} -\frac{\textrm{d}p_{\textrm{th}}}{\textrm{d}y}={-}kT_b\frac{\textrm{d}n_b}{\textrm{d}y} \end{equation}

and

(2.10)\begin{equation} J_bB_x=\frac{\textrm{d}J_b}{\textrm{d} y}\frac{3n_ekT_0}{2J_b}\left[1-\frac{1}{5}\left(1+\varOmega\right)^{2/5}-\frac{4}{5}\left(1+\varOmega\right)^{{-}3/5}\right]. \end{equation}

The specific focusing condition can be obtained by solving

(2.11)\begin{equation} |J_bB_x|>\left|-\frac{\textrm{d}p_{\textrm{th}}}{\textrm{d}y}\right|. \end{equation}

For weak heating, the Lorentz force can be simplified as

(2.12)\begin{equation} J_bB_x\approx\eta_0J_b\frac{\textrm{d}J_b}{\textrm{d}y}\tau, \end{equation}

and the focusing condition is

(2.13)\begin{equation} \frac{n_bq^{2}v_b^{2}}{kT_b}>\frac{1}{\eta_0\tau}. \end{equation}

For strong heating, in the range of $\varOmega \gg 1$ (at least $\varOmega > 44$), where $|y|\ll \sqrt {\mathrm {ln}\varOmega _{\max }/2}R$, the Lorentz force can be approximately written as

(2.14)\begin{equation} J_bB_x\approx{-}\frac{1}{2}\eta_0J_b\frac{\textrm{d}J_b}{\textrm{d}y}\tau\varOmega^{{-}3/5}, \end{equation}

which, compared with (2.12), has an opposite sign and will cause beam hollowing (Davies Reference Davies2003; Davies et al. Reference Davies, Green and Norreys2006), as discussed above. In the range of $0<\varOmega <44$, there is a strong focusing magnetic field, which potentially can focus the beam. To get an accurate focusing condition in this range, we have to solve inequality (2.11). For the sake of simplicity and ease of use, we can estimate the order of magnitude of the Lorentz force by letting $\varOmega \approx 1$, which corresponds to $|y|\approx \sqrt {\mathrm {ln}\varOmega _{\max }/2}R$. Based on this position and making use of (2.9) and (2.10) and inequality (2.11), we can get an approximate focusing condition for strong heating,

(2.15)\begin{equation} \frac{qv_b}{kT_b}>\sqrt{\frac{1}{\eta_0\tau n_ekT_0}}, \end{equation}

where for simplicity, we have discarded the constant coefficient. The above results in inequalities (2.13) and (2.15) show that in the weak heating regime, increasing either beam density or velocity is good for focusing and in the strong heating regime, only increasing the beam velocity facilitates focusing. It should be mentioned that in the theory, we have neglected the effect of the magnetic diffusion, where the time scale is estimated by

(2.16)\begin{equation} \tau_{md}=\frac{\mu_0R^{2}}{\eta}, \end{equation}

where $\mu _0$ is the vacuum permeability. The resistivity is large for bad conductors, but usually is small for ionized plasmas. For widely used electron and ion beams generated by intense lasers, their duration time is much shorter than the magnetic diffusion time. In these cases, it is reasonable to ignore the magnetic diffusion. Next, we will discuss the effect of the beam radius $R$ on focusing. Adding (2.9) and (2.10), we can express the resultant force $f_y$ as

(2.17)\begin{equation} f_y={-}\frac{\textrm{d}p_{\textrm{th}}}{\textrm{d}y}+J_bB_x=\frac{y}{R^{2}}F_y, \end{equation}

where

(2.18)\begin{equation} F_y={-}3n_ekT_0\left[1-\frac{1}{5}\left(1+\varOmega\right)^{2/5}-\frac{4}{5}\left(1+\varOmega\right)^{{-}3/5}\right]+2n_bkT_b. \end{equation}

Once the relative position of the beam is determined, the resultant force $f_y$ is proportional to ${1}/{R}, f_y\propto {1}/{R}$. According to this relation, we are able to draw a conclusion that if a charged particle beam can be focused, a larger radius $R$ will lead to a smaller resultant force $f_y$, a longer focusing time and a longer focusing distance.

Figure 1. A schematic diagram of beam focusing. A beam, whose spot size is 2R, is injected from the left-hand side of the box into a background plasma. During the transport, the beam radius gradually becomes smaller owing to the focusing force and after the focusing distance, the beam finally gets focused.

Figure 2. A schematic diagram of the magnetic field distributions for weak heating (blue solid line) and strong heating (red solid line). The blue and red dashed arrows indicate the force direction of the Lorentz force for weak heating and strong heating, respectively.

Similar focusing conclusions have also been obtained by other researchers (Bennett Reference Bennett1934,  Reference Bennett1955; Kaganovich et al. Reference Kaganovich, Shvets, Startsev and Davidson2001; Dorf et al. Reference Dorf, Kaganovich, Startsev and Davidson2009,  Reference Dorf, Davidson, Kaganovich and Startsev2012). Considering a mixed stream of counterstreaming ions and electrons, Bennett derived the total focusing force by summing over the force of a single particle on another single particle and found that to make the stream focused, the total current should exceed a critical value to overcome transverse spread arising from transverse velocity components (Bennett Reference Bennett1934,  Reference Bennett1955). Based on cold plasma electron fluid equations with Maxwell's equations, assuming an exact charge neutralization of an ion beam by background plasma electrons, Kaganovich et al. (Reference Kaganovich, Shvets, Startsev and Davidson2001) obtained the focusing force which depends on the degree of the current neutralization and it was found that if the radial ion thermal velocity is ignored, the beam can always be focused (Kaganovich et al. Reference Kaganovich, Shvets, Startsev and Davidson2001).

Though to some extent, our conclusions appear to be alike, the physical image we describe and the model we use to calculate the electromagnetic fields in this paper are different. In our model, we treat the plasma as a conductor with a resistivity and assume the beam can be well neutralized by plasma electrons. Because Ohmic heating is considered, the plasma temperature will be increased and plasma resistivity will be decreased. According to the Ohm's law and Faraday's law, the generated magnetic field is not static and will evolve in response to changes in plasma resistivity. Likewise, the focusing force and the focusing condition are also not constant, but vary with the degree of heating.

3. Simulation

To confirm the above theoretical predictions, we have used LAPINS code (Chen et al. Reference Chen, Wu, Ren, Hoffmann and Zhao2020; Ren et al. Reference Ren, Deng, Qi, Chen, Ma, Wang, Yin, Feng, Liu and Xu2020; Wu et al. Reference Wu, He, Yu and Fritzsche2018,  Reference Wu, Yu, Fritzsche and He2020) to make multiple sets of 2D3V simulations of a proton beam propagating in a uniform large-scale hydrogen plasma. As we know for typical particle-in-cell (PIC) simulations, the plasma frequency needs to be resolved, and moreover the grid size must be comparable to the Debye length to minimize artificial grid heating and suppress numerical instabilities. Therefore, high noise levels and high computational requirements arising from the operation on the shortest time and length scales greatly limit the applications of PIC methods to large-scale simulations.

To efficiently investigate the transport of intense charged particles, the significant update to the LAPINS code has recently been made. For our interests, the current density of proton beams is so high that electromagnetic effects need to be considered. However, when compared with the density of the plasma target, the current density of the proton beam is still several orders of magnitude lower. For this particular case, instead of solving the full Maxwell's equations, we couple Ampere's law, Faraday's law and Ohm's law to update the electric and magnetic fields. To be specific, plasma ions and the injected beam particles are treated by the PIC method, while plasma electrons are treated as a fluid, which is solved by Ampere's law, $\boldsymbol {J}_e=(1/2{\rm \pi} )\boldsymbol {\nabla }\times \boldsymbol {B}-(1/2{\rm \pi} )(\partial \boldsymbol {E}/\partial t)-\boldsymbol {J}_b-\boldsymbol {J}_i$, where $\boldsymbol {B}$ is the magnetic field, $\boldsymbol {E}$ is the electric field, $\boldsymbol {J}_b$ is the beam current and $\boldsymbol {J}_i$ is the plasma ion current. The electric field is obtained by Ohm's law, $\boldsymbol {E}=\eta \boldsymbol {J}_e-\boldsymbol {v}_e\times \boldsymbol {B}-\boldsymbol {\nabla } p_e/en_e$, where $\eta$ is the resistivity, $\boldsymbol {v}_e$ is the flow velocity of plasma electrons, $p_e$ is the plasma electron pressure, $n_e$ is the plasma electron density and $e$ is the elementary charge. The magnetic field is obtained by Faraday's law, $\partial \boldsymbol {B}/\partial t=-\boldsymbol {\nabla }\times \boldsymbol {E}$. As only a part of Maxwell's equations needs to be solved, this method is quick, which is useful for large-scale simulations. Furthermore, collisional effects are included based on a Monte Carlo binary collision model, which can deal well with calculating the beam stopping, plasma heating and thermal conduction.

In the simulations, the size of a cell is $0.1\, \mathrm {mm} \times 0.1\, \mathrm {mm}$ and the time step is 0.17 ps, the plasma temperature is 4 eV and density is $10^{5}n_0$, where $n_0=1.1\times 10^{19}\, \mathrm {m}^{-3}$. The beam temperature is 4 eV. The beam moves in the $z$ direction, where the transverse density distribution is a Gaussian profile in the $y$ direction and the radius is 1 mm. The time variable $\tau$ is set to 0.6 ns. In the case with the above parameters kept fixed, we change the beam density $n_p$ and velocity $v_p$, and then investigate their propagations in the plasma. The specific density $n_p$ and velocity $v_p$ of the proton beam vary, as shown in table 1. Under these conditions, the plasma is weakly heated. The length along the $y$ direction is $5\,\mathrm {mm}$. The length along the $z$ direction is several hundred millimetres. It is not fixed because the focusing distance changes as the beam parameters change.

Table 1. In the 30 sets of simulations run by LAPINS code, the values of the proton beam density (in units of $n_0$) and velocity (in units of the light speed $c\approx 3.0\times 10^{8}\, \mathrm {m}\,\mathrm {s}^{-1}$). The first column is the proton beam density parameters, and the second to the sixth columns are the proton beam velocity parameters at different densities.

Figure 3 shows four sets of typical simulation results, from which, we can roughly judge the influences of the beam density, velocity and radius on focusing. In figure 3(a), we see that if the beam density is $1n_0$ and velocity is $0.05c$, it cannot become focused. However, when we increase its density to $10n_0$, it is found in figure 3(b) that the beam can become focused. Similarly, when we increase its velocity to $0.31c$, as shown in figure 3(c), the beam also gets focused. Comparing figures 3(c) and 3(d), we notice that in the case where the other parameters are kept fixed, the focusing distance is different at different radii. In figure 3(c) where the beam radius is 1 mm, the focusing distance is 290 mm, while in figure 3(d) where the beam radius is reduced to 0.5 mm, the focusing distance is shortened to 200 mm.

Figure 3. Beam densities distributions in four sets of simulations under different beam parameters. Initial densities of the beams are $1 n_0$ for (a,c,d) and $10 n_0$ for (b). Initial velocities of the beams are $0.05 c$ for (a,b) and $0.31c$ for (c,d). Initial radii of the beams are $1\,\mathrm {mm}$ for (a–c) and $0.5\,\mathrm {mm}$ for (d). The dashed lines indicate $z=290\, \mathrm {mm}$ for (c) and $z=200\, \mathrm {mm}$ for (d), which are the positions at which the beams eventually become focused.

The simulation results under table 1 parameters and the theoretical results are displayed in figure 4. If $n_p$ and $v_p$ are above the green diamonds, it implies $|J_pB_x |>|-({\textrm {d}p_{\textrm {th}}}/{\textrm {d}y})|$ and the proton beam can be focused. If $n_p$ and $v_p$ are below the green diamonds, it implies $|J_pB_x |<|-({\textrm {d}p_{\textrm {th}}}/{\textrm {d}y})|$ and the proton beam cannot be focused. Qualitatively, as predicted by the theory in § 2, in the weak heating regime, increasing the proton beam density or velocity is beneficial to focusing. The simulation results also conform to this law. For example, we can clearly see in figure 4 that when the proton beam density $n_p$ is set to $0.05n_0$, if the beam velocity $v_p$ is less than or equal to $0.20c$, the beam cannot be focused. With the increase of $v_p$, the proton beam can be focused only when $v_p$ is greater than or equal to $0.31c$. However, when the proton beam velocity $v_p$ is set to $0.10c$, if the beam density $n_p$ is less than or even equal to $0.5n_0$, the beam cannot be focused. Only when $n_p$ is greater than or equal to $1.0n_0$ can the beam be focused.

Figure 4. Simulation and theoretical results. The red crosses represent that the beam is not focused and the blue circles represent that the beam is focused. The black solid line represents the numerical results solved by $({n_bq^{2}v_b^{2}})/{kT_b}={1}/{\eta _0\tau }$ (corresponding to focusing condition (2.13)) at $y=R$. The green diamonds are obtained by the original theory, $|J_pB_x |_\mathrm {{theory}}=|-({\textrm {d}p_{\textrm {th}}}/{\textrm {d}y})|$ at $y=R$.

There is a little difference between our theory and simulations. This may be caused by two aspects. First, for simplicity, in the theory, we considered the beam temperature to be a constant, while actually, in the process of the beam being focused, collisions between beam particles will be enhanced, which can increase the beam temperature and as a result, the beam thermal pressure will also be increased. Second, in the theory, the proton beam velocity $v_p$ is considered to be a constant. However in the simulations, owing to collisions and collective electromagnetic effects, the beam is gradually slowed down while propagating (Li & Petrasso Reference Li and Petrasso1993; Brown, Dean & Singleton Reference Brown, Dean and Singleton2005; Zylstra et al. Reference Zylstra, Frenje, Grabowski, Li, Collins, Fitzsimmons, Glenzer, Graziani, Hansen and Hu2015; Clauser & Arista Reference Clauser and Arista2018; Ren et al. Reference Ren, Deng, Qi, Chen, Ma, Wang, Yin, Feng, Liu and Xu2020). Therefore, the real velocity of the proton beam in the simulations will be less than the initial theoretical velocity, $v_\mathrm {sim}< v_\mathrm {theory}$. Then from (2.12), it is easy to prove that the Lorentz force in the simulations will be smaller than the theoretical one, $|J_pB_x |_\mathrm {sim}<|J_pB_x |_\mathrm {theory}$. In this way, the theoretical line $|J_pB_x |_\mathrm {theory}=|-({\textrm {d}p_{\textrm {th}}}/{\textrm {d} y})|$ will be located in the range of $|J_pB_x |_\mathrm {sim}<|-({\textrm {d}p_{\textrm {th}}}/{\textrm {d}y})|$.

From figures 3(c) and 3(d), we can roughly judge the relation between the focusing distance and the beam radius. To further test the influence of the beam radius $R$ on focusing, we have done more simulations using LAPINS code. The parameters of the background plasma remain unchanged, and the parameters of the incident proton beam are a density of $1n_0$, velocity of $0.31c$ and temperature of 4 eV. The value of the beam radius $R$ varies from 0.5 mm to 1.0 mm. The simulation results are shown in figure 5. Obviously, as $R$ increases, the focusing distance of the proton beam gradually becomes longer. This law is also consistent with the previous theoretical predictions.

Figure 5. Variation of the focusing distance as the beam radius $R$ varies.

The double-peak focusing structure predicted by the theory in the strong heating regime has not been observed yet. This may be caused by thermal conduction arising from collisions. Moreover, we found that if we further set the beam radius $R$ to a much larger number, such as 4 mm, as shown in figure 6, the proton beam filamentation phenomenon would occur. There have been many important discoveries and conclusions on electron beam filamentation (Benford Reference Benford1973; Lee & Lampe Reference Lee and Lampe1973; Bret, Firpo & Deutsch Reference Bret, Firpo and Deutsch2005; Allen et al. Reference Allen, Yakimenko, Babzien, Fedurin, Kusche and Muggli2012; Wang, Hu & Wang Reference Wang, Hu and Wang2020), but ion beam filamentation has not been thoroughly studied yet, which should be further investigated in the future.

Figure 6. Densities distributions of the proton beams in the case of $R=4\, \mathrm {mm}$. Initial densities and velocities of the proton beams are (a) $n_p=1n_0,\ v_p=0.31c$ and (b) $n_p=10n_0,\ v_p=0.87c$, respectively.