Adiabatic sheath model for beam-driven blowout plasma channels

Beam-driven plasma wakefield accelerators (PWFAs) have the capacity to accelerate particles to tens of GeV energies in meter-scale structures, making them competitive candidates for future high-energy accelerators that are used for next-generation particle colliders [1, 2, 3, 4, 5]. In a PWFA, when a high current, relativistic electron beam, called the drive beam, interacts with an underdense plasma, plasma electrons are completely radially expelled while the ions remain immobile due to their much larger mass, forming an electron-free ion channel along the propagation axis. The physics of plasma response in this regime is strongly nonlinear and the regime has been referred to as the blowout or bubble regime. In this regime, the acceleration (or deceleration) field is uniform and the focusing field is linear in regard to the radial offset from the axis [6, 7, 8, 9, 10].

On the boundary of the bubble, the expelled electrons form a narrow sheath, where the electron density and current profiles steepen with a large value and decay in a small thickness [11, 12, 13]. The sheath structure depends on the properties of the drive beam and is strongly correlated with the creation of the plasma blowout, the structure of wakefield inside the blowout, and the processes of particle injection into the wakefield [14, 15, 16].

Existing theories have used rectangular or exponential distributions to simplify the sheath shape [17, 18, 19, 20, 21, 22]. Dalichaouch et al. [23] proposed a multisheath model, which employs two sheath layers to describe the wakefield more accurately. However, these models require the prior knowledge of the sheath thickness, which cannot be obtained self-consistently. Recently, Golovanov et al. [24] have developed a blowout theory from the energy conservation point of view. They have assumed a $\delta$-function distribution for the blowout sheath, and found the evolution function for the blowout channel radius $r_{\delta}$ as

where $A=1+r_{\delta}^{2}/4$, $B=1/2+1/r_{\delta}^{2}$ and $C=1/4$. $\xi=ct-z$ is the longitudinal co-moving coordinate, $c$ is the speed of light in vacuum, $t$ is time, and $z$ is the longitudinal coordinate. The drive term of the equation can be written as $\Lambda=2I/I_{A}$, where $I$ is the instant current of the beam, $I_{A}=4\pi\epsilon_{0}m_{e}c^{3}/e\approx 17\ \rm kA$ is the Alfvén current, $\epsilon_{0}$ is the vacuum permittivity, $m_{e}$ is the electron rest mass, and $e$ is the elementary charge [25, 26, 27]. If derivatives of $r_{\delta}$ are assumed to be zero, one can obtain the balancing radius of the channel based on the $\delta$-sheath theory as

By contrast, the charge neutralization radius,

has been widely accepted as the balancing radius from the electrostatic point of view [28, 10, 26].

In this work, we develop a self-consistent model for the channel sheath using the adiabatic assumption. The sheath equations are examined in two cases with the help of particle-in-cell (PIC) simulations. We demonstrate that our model provides a more accurate description of blowout channel balancing, as its predicted channel radius falls between the electrostatic neutralization radius and that predicted by the $\delta$-sheath theory.

Throughout the paper, we adopt the plasma normalization units, in which charge is normalized to $e$, mass to $m_{e}$, velocity to $c$, density to plasma density $n_{p}$, time to $\omega_{p}^{-1}$, length to $c/\omega_{p}$, electric field to $m_{e}c\omega_{p}/e$, magnetic field to $m_{e}\omega_{p}/e$, electrostatic potential to $m_{e}c^{2}/e$ and vector potential to $m_{e}c/e$, where $\omega_{p}=\sqrt{e^{2}n_{p}/\epsilon_{0}m_{e}}$ is the plasma frequency.

We consider a highly relativistic drive electron beam with density $n_{b}\left(\xi,r\right)$ that expels the plasma electrons to form an ion channel with a radius $r_{c}\left(\xi\right)$, where $r$ is the radial coordinate, and cylindrical symmetry is assumed. We assume that the drive beam is narrower than the channel, so that there is no drive charge outside the channel, or $n_{b}=0$ for $r>r_{c}$. Under these assumptions, the drive term can also be written as $\Lambda\left(\xi\right)=\int_{0}^{\infty}n_{b}rdr$, which has the physical meaning of the (normalized) linear density of the drive beam.

Furthermore, we assume the radius of the ion channel changes adiabatically, which means $\left(dr_{c}/r_{c}\right)/d\xi\ll 1$, $\left(d\Lambda/\Lambda\right)/d\xi\ll 1$, and $\partial_{\xi}\ll\partial_{r}$ for all field variables, as illustrated in Fig. 1. Under the cylindrical symmetry and the quasi-static approximation [29, 30], the pseudo-potential of the wakefield obeys the following Poisson-like equation [31, 10]

where $\psi=\varphi-A_{z}$ is the pseudo-potential, $\varphi$ is the electrostatic potential, $A_{z}$ is the longitudinal component of the vector potential,

is the source term, $\rho$ is the charge density, $J_{z}$ is the longitudinal component of the current density, $n_{e}$ is the plasma electron density, and $v_{z}$ is the longitudinal velocity of the plasma electrons. Because $n_{e}=0$ for $r<r_{c}$, we can write the form of the pseudo-potential inside the blowout

where

is to be determined. It can be derived that $\psi_{c}=0$ with the $\delta$-function distribution, i.e. $S=1-H\left(r-r_{c}\right)-r_{c}\delta\left(r-r_{c}\right)/2$, where $H$ is the Heaviside step function, and $\delta$ is the Dirac delta function. But in this work, $S$ is finite everywhere, so that $\psi_{c}$ is non-zero, and $\partial\psi/\partial r$ is continuous everywhere because of Eq. (4), which means

For a plasma electron which is at rest before the drive beam arrives, there is a constant of motion [30]

where $\gamma=1/\sqrt{1-v_{z}^{2}-v_{r}^{2}}$ is the Lorentz factor, $v_{z}$ and $v_{r}$ are the longitudinal and radial components of the velocity of the electron, respectively. This equation can be solved under the adiabatic assumption, which means $v_{r}\approx 0$, as

Based on the adiabatic assumption, the transverse electromagnetic field can be written using the Gauss’ Law and Stokes’ Theorem

The transverse force $F_{r}=-E_{r}+v_{z}B_{\theta}$ should be balanced for $r>r_{c}$, which means

The channel radius is determined by letting $r=r_{c}$ in Eq. (13), as

where

We see that $r_{c}$ reduces to $r_{n}$ if $\psi_{c}=0$, and to $r_{\delta 0}$ if $\psi_{c}\rightarrow\infty$. Then, we know $r_{n}<r_{c}<r_{\delta 0}$ if $\psi_{c}$ is non-zero and finite, because $r_{c}$ monotonically increases with $\psi_{c}$.

The sheath equations defined by Eqs. (4), (10) and (13) are solvable using the left boundary condition Eq. (8) and the right boundary condition

which is a natural requirement that the plasma disturbance is local.

We use the shooting method to numerically solve the equations. In other words, the numerical integral is performed from $r=r_{c}$ to a sufficient large value of $r$, while $\psi_{c}$ is predicted and corrected in iterations, so that Eq. (16) is satisfied to a certain accuracy. Practically, there is a critical value of $\psi_{c}$ for each $r_{c}$, which prevents divergence at large $r$. The Python script for the numerical solution is attached as a Supplementary Material, and is maintained on GitHub [32]. The numerical solution of $\psi_{c}$ with different $r_{c}$ is shown in Fig. 2. The polynomial fit,

is a good estimate for $r_{c}\lesssim 8$ as shown in Fig. 2(a).

After $\psi_{c}$ is obtained, the electron density at the sheath boundary can be predicted. By taking derivative of Eq. (13), we obtain

At $r=r_{c}$, the integral is zero, the derivative of $v_{z}$ can be determined by Eqs. (8) and (10), and $\Lambda$ can be written as the function of $r_{c}$ and $\psi_{c}$ according to Eq. (14). As a result,

where $\psi_{*}=1+\psi_{c}$.

Next, we study the behavior of $\psi$ and $\psi_{c}$ in the limit $\psi_{c}\ll 1$. This also means $r_{c}\ll 1$, $\Lambda\ll 1$, and $\psi\ll 1$. By expanding Eqs. (10) and (18) and keeping only the 1st order terms of $\psi$, we have $v_{z}\approx-\psi$ and $n_{e}\approx 1$. Thus, the linear form of Eq. (4) is obtained

The solution satisfying the right boundary condition Eq. (16) is $\psi=R\left(r_{c}\right)K_{0}\left(r\right)$, where $R\left(r_{c}\right)$ is a constant of $r$ to be determined, and $K_{0}\left(r\right)$ is the modified Bessel function of the second kind of order zero [33]. By also considering the left boundary condition Eq. (8), we may get the expression of $R\left(r_{c}\right)$ to the lowest order of $r_{c}$, and find

Thus,

or $\psi_{c}\approx-\frac{r_{c}^{2}}{2}\left(\ln\frac{r_{c}}{2}+0.577...\right)$, where the second summand is the Euler-Mascheroni constant. We can see in Fig. 2(b) that Eq. (22) is satisfactory with the numerical results for $r_{c}\lesssim 0.28$, while Eq. (17) is better for $r_{c}\gtrsim 0.28$.

Once $\psi_{c}$ is determined, the one-to-one mapping between $\Lambda$ and $r_{c}$ can be obtained according to Eq. (14). For comparison, PIC simulations have been performed using the quasi-static code QuickPIC [34, 35]. In the simulations, we consider a 5 GeV drive electron beam, transversely positioned at $x=y=0$ (i.e. the center of the simulation domain). It has a sufficiently small transverse size, while its longitudinal size is the same as the entire simulation domain. Two sets of simulations with different drive beam current profiles have been conducted. One is called the near-adiabatic profile, where $\Lambda$ increases linearly with a sufficiently small slope. The other is called the near-stationary profile, where $\Lambda$ increases from 0 to a plateau, and the increasing profile and length are designed so that the oscillation of the channel radius in the plateau is not severe.

In the near-adiabatic simulations, we have chosen the maximum $\Lambda_{\max}=4$ for the large $r_{c}$-scale shown in Fig. 2(a) and in Fig. 3(a), and $\Lambda_{\max}=0.4$ for the small $r_{c}$-scale shown in Fig. 2(b). The simulation domain has a size of 200 in $\xi$-direction, which is the propagation direction, and 40 (10) in $x$- and $y$-direction for the $\Lambda_{\max}=4$ (0.4) case. The numbers of cells is chosen as 2048, 1024 and 1024, respectively. The near-adiabatic simulation results for $\psi_{c}$ vs. $r_{c}$ are shown in Fig. 2, which are in good agreements with our model.

In the near-stationary simulations, $\Lambda$ increases from 0 to the plateau value $\Lambda_{0}$ within $0\leq\xi\leq 25$, and keeps the plateau value in the remaining of the simulation domain, as shown in Fig. 3(c). Cubic splines have been designed to smooth the transitions near $\xi=0$ and 25. Here, the simulation domain has a longitudinal size of 48, while the transverse size is adjusted in accordance with $\Lambda_{0}$, such that $\psi\rightarrow 0$ at the simulation boundary is guaranteed ($\partial\psi/\partial r\approx 0$ at the boundary), and the transverse resolution is adequate to resolve the sheath. Again, the number of cells is chosen to be $2048\times 1024\times 1024$. $\Lambda_{0}$ has been scanned from $2.5\times 10^{-5}$ to 6.25 in the simulations. The results of $\psi_{c}$ vs. $r_{c}$ are shown in Fig. 2, which are also in good agreements with our model.

Plasma density slices at $y=0$ for two of the above mentioned simulations are shown in Fig. 3(b) and 3(d). We compare the simulations with $r_{\delta 0}$ obtained by Eq. (2), $r_{n}$ obtained by Eq. (3) and $r_{c}$ obtained by Eq. (14). We see that the channel radius in simulations best matches with $r_{c}$, while $r_{n}$ is an underestimate and $r_{\delta 0}$ is an overestimate.

Furthermore, the channel radius slightly oscillates around $r_{c}$ in the near-stationary case because of the imperfect balancing, as shown in Fig. 3(d). This oscillation behavior can be numerically solved using Lu’s model [10, 17], which is similar to Eq. (1), but the parameters are $A=1+\left[1/4+\beta/2+\left(r_{c}/8\right)\left(d\beta/dr_{c}\right)\right]r_{c}^{2}$, $B=1/2+\left(3/4\right)\beta+\left(3r_{c}/4\right)\left(d\beta/dr_{c}\right)+\left(r_{c}^{2}/8\right)\left(d^{2}\beta/dr_{c}^{2}\right)$ and $C=\left[1+1/\left(1+\beta r_{c}^{2}/4\right)^{2}\right]/4$. The variable $\beta$ needs prior knowledge to the sheath thickness, and was recommended to be $\beta=1/r_{c}+0.1$. Although Lu’s model is in good agreement with the simulation in the short drive cases (drive beam is shorter than or comparable to the plasma wavelength), it does not predict the oscillation of $r_{c}$ very well in the long drive cases as shown in Fig. 3(d).

To improve Lu’s model, we replace $\beta$ by

where $\psi_{c}$ is expressed in a piecewise manner by Eq. (22) for $r_{c}<0.28$ and by Eq. (17) for $r_{c}>0.28$, and call it the combined model. The combined model has been applied to the near-stationary cases, and its numerical solutions have the best agreement with the simulations compared with other models, with one example shown as the violet dotted line in Fig. 3(d).

In conclusion, we have introduced a self-consistent model for the sheath structure of the blowout plasma channel driven by an electron beam based on the adiabatic assumption. The model consists of the equation for the pseudo-potential Eq. (4), the longitudinal velocity of the sheath electrons Eq. (10), and the transverse force balancing condition Eq. (13), which is solvable with the boundary conditions Eqs. (8) and (16). An analytical solution is obtained in the small-blowout limit as Eq. (21), and a general solution is retrieved numerically by the shooting method. The PIC simulations show that the radius obtained by Eq. (14) is a better estimate of the channel balancing radius than that from previous neutralization or $\delta$-sheath models.

It should be noted that the accuracy of our model is limited in short drive beam cases. Nevertheless, our model is very suitable for long drive beams. Especially, by combining our model with Lu’s model, we have obtained the most correct description of the oscillation of the channel radius. The better understanding of the pseudo-potential at the sheath boundary and inside the blowout for long drive cases is crucial for high transformer ratio and future high-energy PWFA studies [5, 36, 37].