Wake excited in plasma by an ultra-relativistic point-like bunch

A plasma flow behind a relativistic electron bunch propagating through a cold plasma is found assuming that the transverse and longitudinal dimensions of the bunch are small and the bunch can be treated as a point charge. In addition, the bunch charge is assumed small. A simplified system of equations for the plasma electrons is derived and it is shown that, through a simple rescaling of variables, the bunch charge can be eliminated from the equations. These equations have a unique solution, with an ion cavity formed behind the driver. The equations are solved numerically and the scaling of the cavity dimensions with the driver charge is obtained. A numerical solution for the case of a positively charged driver is also found.


I. INTRODUCTION
Plasma wakefield acceleration (PWFA) in a highly nonlinear "blowout" regime [1] offers a promising path toward high-gradient compact future accelerators for numerous applications. Over the last two decades various aspects of physics involved into the complicated dynamics of plasma excited by a high-charge driver beam moving with a relativistic velocity has been studied. Unfortunately, due to nonlinear nature of the problem, not much can be done analytically beyond the formulation of the governing system of equations for the plasma flow. Usually, these equation are solved numerically, and computer simulations are widely accepted as a main tool to study PWFA [2][3][4].
While computer simulations are indispensable for testing new ideas and making connections to the experiment, there is also a strong need for analytical techniques that can allow for quick estimates and provide scaling relations for the governing parameters of the problem. Another application of the analytical approach would be a calculation of the longitudinal and transverse wakefields for the witness bunch accelerated behind the driver beam and the subsequent analysis of the beam instabilities. The importance of such an analysis for the feasibility of the PWFA concept has been recently emphasized in Ref. [5].
The perceived need to create an approximate analytical description of the PWFA process has lead some of the researchers to developing approximate models that incorporate important features observed in simulations [6,7]. For example, in Ref. [7] an equation for the shape of the plasma bubble is derived based on the assumption that the bubble boundary is formed by a single electron trajectory and the plasma current is localized in a narrow sheath near the boundary. While being a useful complement to numerical simulations, these models usually lack the rigor of the fundamental approach and often have an uncertain range of applicability. In particular, it is not clear whether the model of Ref. [7] can be used for the calculation of wakefields inside the bubble, as assumed in Ref. [5].
An attempt to develop a consistent analytical model for calculations of the plasma flow in the blowout regime is made in this paper. We note that a part of the complexity of the problem is associated with the fact that there are several parameters that control the dynamics of the bubble. Two of these parameters define the size of the driver bunch: σ r k p and σ z k p , where σ r and σ z are the radial and longitudinal dimensions of the driver and k p = ω p /c = √ 4πn 0 r e , with n 0 the plasma density and r e = e 2 /mc 2 the classical electron radius. In our model, we consider the case of a small-size driver beam, σ r k p 1 and σ z k p 1, and formally take the limit σ r k p , σ z k p → 0. This limit eliminates the two parameters from the consideration and simplifies our analysis. While this regime may be away from an optimal setup for achieving the highest accelerating gradient, it can be considered as a first approximation to more realistic situations.
We also assume an ultra-relativistic driver beam moving through the plasma with the speed of light c. This leaves only one external dimensionless parameter-the dimensionless charge of the beam ν-in the problem [8], where N is the number of particles in the bunch. The second equality in this relation shows that within a factor of 1 3 the parameter ν is equal to the ratio of the number of particles in the bunch to the number of plasma particles in the sphere of radius k −1 p . To further simplify the problem we focus this paper on the limit of small charges, ν 1, where, as we will show below, a remarkably simple scaling for the main bubble parameters can be derived.
This paper is organized as follows. In Section II, we recapitulate the main equations that govern the plasma flow and in Section III we formulate initial conditions for these equations. In Section IV, we develop a simple ballistic model that demonstrates how a bubble is developed when the plasma self fields are neglected. In Section V, the general system of equations of Section II is simplified in the limit ν 1, rescaled and then solved numerically. Various properties of the bubble obtained from the numerical solution are discussed and related to the ballistic model of Section IV. In Section VI, we present analytical formulas that remarkably well approximate some of the properties of the numerical solution of Section V. In Section VII, we derive the longitudinal electric field E z on the axis of the bubble. In Section VIII, a numerical solution for the plasma flow behind a positive driver is presented. The main results of the paper are summarized in Section IX.

II. FORMULATION OF THE PROBLEM
In this Section, we formulate the equations that describe plasma dynamics behind a point-like driver moving with the velocity of light along the z axis through a cold plasma of constant density n 0 . The charge density of the driver in the cylindrical coordinate system is ρ dr = ∓(Ne/2πr) δ(z−ct) δ(r), where r is the distance from the axis z and the minus sign refers to an electron driver and the plus-to a positively charged one.
Following the standard convention, we normalize time to ω −1 p , length to k −1 p , velocities to the speed of light c, and momenta to mc. We also normalize fields to mcω p /e, the charge density to n 0 e, the plasma density to n 0 , and the current density to en 0 c. Here e is the elementary positive charge. In these dimensionless units, the charge density of the driver beam is given by the following expression: where ξ = ct − z. The uniform plasma at rest occupies the region in front of the beam ξ < 0. Our goal is to find the plasma flow and the electromagnetic field behind the driver, in region ξ > 0. Due to the symmetry of the problem, the plasma flow is axisymmetric. In a steady state, when one can neglect the transients associated with the entering through the plasma boundary, all the fields depend on z and t in combination ct − z, that is through the variable ξ. From Maxwell's equations, we find the following equations for the non-zero components of the electric and magnetic fields E r , E z and B θ in the cylindrical coordinate system: where j is the plasma current density. The equation for ∇ · E reads: These equations are complemented by the continuity equation for the electron flow In this equation we used nu for the flow density of the electrons, which is valid if at every point the electron flow is characterized by a unique velocity u. In general, however, the flow behind the point-like driver may have several streams at a given point, that is there are several values of u at each point in space. In this case, nu in Eq. (5) should be replaced by s n s u s where the summation goes over the different values of u s with the corresponding densities n s . For notational simplicity, we drop the sum sign in what follows, unless it is specifically indicated. The equations of motion for the plasma electrons in the dimensionless variables are where p r and p z are the radial and longitudinal components of the momentum vector. Eqs. (3)-(6) should be supplemented by initial conditions for n and u at ξ = 0. These conditions are derived in the next section.

III. PLASMA CROSSING OF THE DRIVER FIELD
The driver beam propagating through the plasma with the speed of light carries an electromagnetic field that is localized in an infinitesimally thin "pancake" region in the x − y plane ξ = 0. In what follows, we will call the field in this region the "driver field". At a first glance, it might seem that due to the thinness of the driver field region and the fact that it moves with the speed of light, the plasma does not have time to modify it, so the field in this region is the same as in the free space: where This assumption, however, is incorrect. As was shown in Ref. [8], the radial currents induced when electrons cross the driver field modify it in such a way that where K 1 (r) is the modified Bessel function of the second kind. Using the asymptotic expressions for function K 1 (x), for small distances, r 1 (r k −1 p in dimensional units), we recover from (9) the free space expression A(r) ≈ ∓2ν/r; in the opposite limit, r 1, the plasma shields the driver field to exponentially small values, A(r) ∝ ∓e −r / √ r. Plasma electrons from region ξ < 0 encounter the driver field, receive radial and longitudinal kicks and change their momentum from the initial zero at ξ = 0 − to nonzero values p r0 and p z0 at ξ = 0 + . The values p r0 and p z0 can be easily found if one notices that due to an infinitesimally small thickness of the region of the driver field in the ξ direction, at a given radius r, one can locally approximate it by a plane electromagnetic wave. A solution for a point charge motion in a plane electromagnetic field can be found in the literature (see, e.g., [9,10]); from this solution it follows that where γ 0 = 1 + p 2 r0 + p 2 z0 is the electron energy corresponding to the momentum (p r0 , p z0 ). Note the singularities in p r0 and p z0 when r → 0; these singularities are due to our assumption of the point driver (they do not appear in the case of a finite value of the driver radius). As we will see below, they do not cause problems in finding the plasma flow behind the driver.
Integrating the continuity equation (5) through the driver field and noting that both n and v r have finite values at ξ = 0 (although they are discontinuous at this point) we find that the product n(1 − v z ) does not change from ξ = 0 − to ξ = 0 + . This allows us to find the electron density, n * , at ξ = 0 + , immediately after the driver field region, where v z0 = p z0 /γ 0 . In the limit r → 0, we have v z0 → 1 and the density (11) also has a singularity at r = 0. After crossing the driver field at ξ = 0 the plasma flow is governed by a self-consistent field that is determined by Eqs. (3). Eqs. (10) and (11) provide the initial conditions for Eq. (5) and the equations of motion (6). Note that the whole system of equations and the initial conditions have only one dimensionless external parameter, ν. While equations in this and the previous sections are valid for arbitrary ν, in the rest of this paper we will focus on the limiting case of a small driver charge, ν 1. We will see that in this limit a simple rescaling of the variables eliminates ν from the equations and shows a universal pattern of the plasma flow with a bubble whose shape and characteristic parameters are easily found from a simple computational problem.

IV. BALLISTIC APPROXIMATION FOR TRAJECTORIES
Before we analyze the full picture of the plasma flow in the next section, we consider here a simpler problem of the flow at small distances behind the driver beam, ξ 1. At these distances we can neglect the effect of the plasma self-fields on the trajectories of the plasma electrons. We call this approximation the "ballistic" regime of plasma motion; it assumes that the plasma electrons are moving with constant velocities defined by Eqs. (10). It provides an insight into the formation of the bubble and gives an estimate of its dimensions that will be used in the next section.
In the ballistic approximation, particle trajectories in the r-ξ plane are straight lines starting at ξ = 0 with an offset r 0 and going at an angle arctan[v r0 /(1 − v z0 )] to the horizontal axis, As was mentioned in the previous section, we are interested in the regime ν 1. Assuming, in addition, that the region of interest is limited by ν r 1, we find that one can use Eq. (8) for A, and also A 1. It then follows from Eqs. (10) that the electrons are moving with non-relativistic velocities v r0 1 and v z0 1. Hence and electron trajectories in the ballistic approximation are given by In Eqs. (13) and (14) we have chosen a minus sign specifying an electron driver.
A set ballistic trajectories is plotted in Fig. 1 where the vertical axis r is normalized by this normalization, the family of the trajectories has a universal shape that does not depend on ν. The plot clearly shows the formation of the bubble boundary (shown by the red line) below which all the electrons are evacuated and the electron density is zero. An equation for the boundary, r b (ξ), is found as an envelope for the trajectories with different r 0 , dr/dr 0 = 0: As is seen from Fig. 1, outside of the bubble there are two electron streams at each point r, ξ, corresponding to two different trajectories passing through the point. One of the trajectories touches the bubble boundary after, and the other before, it reaches the point r, ξ.
One can also find the density n(r, ξ) of the plasma in region r > r b (ξ). Derivation of a formula for n(r, ξ) is given in Appendix A. The density is zero below the bubble boundary, r < r b , and has a singularity n ∝ (r − r b ) −1/2 as r → r b (see Eq. (A5)).

V. SMALL-CHARGE REGIME OF THE BUBBLE
There are several important simplifications that can be carried out in the limit ν 1 as indicated by the analysis of the ballistic approximation in the previous section. First, the region of the nonlinear plasma flow (the bubble and the adjacent region of a large density perturbation), as it turns out, is localized at r 1. This allows us to use Eqs. (10) with A given by (8); in other words, due to a small transverse size of the flow region, we can neglect the shielding effect in the driver field. Second, a typical value of A in this regime is also small, A 1, and plasma velocities in the dominant part of the flow are non-relativistic. Using these observations, we replace the initial conditions (10) by the following ones: We can obtain a crude estimate of the maximal bubble radius, r bm , and the position ξ m in the r − ξ plane where it is attained from the following consideration. The radius r bm is where the ballistically streaming trajectories start to bend toward the axis and finally collapse under the influence of the ion focusing field. A plasma electron traveling from ξ = 0 at r ∼ r bm will receive a radial kick from the ion electric field E r comparable to its initial radial momentum p r0 . In dimensional units, E r ∼ en 0 r bm , and an electron traveling distance ξ m with the speed c − v z ∼ c changes its momentum by e 2 n 0 r bm ξ m /c. Equating this to p r0 we obtain where on the right-hand side we used the initial radial momentum from (10). On the other hand, for an estimate, we can use the relation (15) from the ballistic approximation in which r b is replaced by r bm and ξ by ξ m . From these two equations and Eq. (8) for A we find (where we now return to the dimensionless variables) and ξ m ∼ 1. Because we assume ν 1, we see that the transverse size of the bubble is much smaller than its longitudinal extension, r m ξ m ,-the bubble has a spindle-like shape. Note also that A from (8) at r ∼ r m is on the order of |A| ∼ √ ν 1, and hence the plasma electron motion is indeed nonrelativistic, as we have assumed above.
Due to the smallness of the electron velocity we can neglect the magnetic force in Eq. (6) and replace it by a simpler one where we used d/dt = (1−v z )d/dξ ≈ d/dξ. Particle orbits r(ξ, r 0 ) are determined by dr/dξ = v r . The electric field in this approximation can be found from Eq. (4) in which the term with E z is neglected as of higher order in parameter √ ν, Eqs. (19) and (20) supplemented by the continuity equation (5) (in which v z in the first term can now be neglected) and the initial conditions (16) constitute a full set of equations. An important feature of these equations is that the only external parameter of the problem, ν, can be eliminated by a rescaling of the variables. Indeed, it is easy to see that changing the variables: does not change the equations but eliminates ν (replaces it by 1) in Eq. (16).
We solved numerically the rescaled equations by launching a large number N tr of trajectories (up to N tr = 5000) at ξ = 0 with a small radial step ∆r 0 uniformly distributed fromr 0 = ∆r 0 tõ r 0 = N tr ∆r 0 (typical ∆r 0 ≈ 0.003). The equations of motion for macroparticles were approximated as d dξṽ with the initial conditionsr where the macroparticle number i = 1, 2 . . . N tr , and η is the Heaviside function, η(x) = 0 at x < 0 and η(x) = 1 at x > 0. The radial electric field was calculated using Eq. (20), by integrating the electron density over r, with the electron density n found from the conservation of the number of macro particles. The equations of motion were solved numerically using an adaptive Runge-Kutta integrator. bubble radius (the maximal value or r on the boundary) is reached at ξ m = 1.57 and is r bm = 2.82 √ ν.
The total bubble length (from ξ = 0 to the point where the bubble boundary collapses on the axis) is ξ b = 3.8 and does not depend on the parameter ν. Note that there are exactly two electron trajectories passing though each point (r, ξ) outside of the bubble which means a two-stream flow at each point. In this regard, the plasma flow is similar to the ballistic model. In the region of not very large values of ξ, ξ 0.5, the trajectories resemble the ballistic ones from Fig. 1. From the analysis of Fig. 2, it follows that the bubble boundary before its maximum r bm (that is ξ < ξ m ) is comprised of many trajectories for which the boundary is an envelope-the same mechanism of the boundary formation as in the case of the ballistic model from Section IV. After the maximum, the boundary is represented by a single trajectory of an electron launched initially at r 0 = 1.44 √ ν. Fig. 3 shows a small number of trajectories that clearly demonstrate a transition from an envelope structure to a single trajectory boundary (it also shows the boundary calculated in the ballistic approximation). After the maximum, it is represented by a single trajectory bent toward the axis by the focusing electric field of the ions.
It is interesting to note that the envelope boundary is different from the assumption made in Ref. [7] where a differential equation for the bubble boundary was derived assuming that the boundary coincides with an electron trajectory. Our solution in Fig. 2 shows that the derivation of [7] is not applicable to our case of a small driver charge and dimensions and hence has a limited range of validity. Fig. 4 shows plots of the radial electron density distribution at several values of the longitudinal coordinate ξ. One can see a zero density inside the bubble with a density jump at the bubble boundary and a gradual fall off to unity as r → ∞. Before the maximum of r b (curves 1 and 2 in Fig. 4) the electron density at the boundary has an infinite value, with the singularity of the same type, n ∝ (r − r b ) −1/2 , as in the ballistic model. After the maximum (curves 3, 4 and 5), the density value on the boundary is finite. Note that the density at the boundary increases from curve 4 to 5 due to the decreasing radius r b ; it becomes infinite at ξ = ξ b = 3.8 where r b = 0.
The fact that the density perturbation in the limit r → ∞ is relatively small, n − 1 1, indicates that the electron flow at large radial coordinates can be described in linear approximation. We will use this observation in Section VII in calculation of the longitudinal electric field inside the bubble.

VI. ANALYTICAL CONSIDERATION
It seems unlikely that the numerical solution of the previous section can be also obtained analytically, however, we found approximate formulas that remarkably well agree with the numerical results. These formulas are obtained if we replace the ballistic trajectories (14) of Section IV with electron orbits derived in linear theory. In linear approximation, electrons oscillate with the plasma frequency, and their orbits, in our dimensionless variables, are given by the following equations: where for v r0 we use Eq. (16), v r0 = 2ν/r 0 . Strictly speaking, we expect these equations to be valid at a large distance from the axis, r 1, but we will now assume that we can use them within a quarter of the plasma period, 0 < ξ < π/2, for all values of r 0 . In this way, replacing the straight orbits with the oscillating ones, we overcome the deficiency of the ballistic approximation that completely neglects the plasma self field. Note that in the limit ξ 1 Eqs. (24) reduce to Eqs. (14). The envelope of trajectories (24) is found from the equation dr/dr 0 = 0 which gives Using this equation in the interval 0 < ξ < π/2 we find for the maximal radius of the bubble and the location of the maximum: in a remarkable agreement with the numerical values for ξ m = 1.57 and r bm = 2.82 √ ν of the previous section.
As was mentioned in the previous section, in region ξ > ξ m = π/2, the bubble boundary is a single electron trajectory. This trajectory is easily found from Eqs. (19), (20) in which we now set n = 0 (because the electron density is zero inside the bubble). The electric field on this trajectory is E r = 1 2 r which gives for the orbit the following equation Its solution in the region ξ > ξ m that matches the initial value r b (ξ m ) = r bm with an addition condition r b (ξ m ) = 0 is where we have used Eq. (26). As follows from this formula, the bubble boundary collapses on the axis at ξ = π 2 (1 + √ 2) = 3.79, again in a remarkable agreement with the numerical value for ξ b from the previous section.
Comparison of the analytical expressions for the bubble boundary (25) and (28) with the numerical solution is shown in Fig. 5-they agree very well in the whole interval 0 < ξ < ξ b .

VII. LONGITUDINAL ELECTRIC FIELD IN THE BUBBLE
We will now discuss how to calculate the longitudinal electric field E z (ξ) inside the bubble, on the axis r = 0. We first integrate Eq. (3c), where we chose the integration constant from the requirement that E z → 0 when r → ∞. The sum over index s in this equation goes over the two electron streams at a given point r, ξ, as discussed in Section V. Unfortunately, an attempt to apply this equation to the numerical solution of the previous section shows that the integral logarithmically diverges at the upper limit, because v r at large distances falls off as 1/r. This happens because we neglected the screening effect of the plasma at r 1. To include the screening effect into calculations, we split the integration interval in (29) into two parts choosing some value r * such that r bm r * 1: the integral from r to r * is computed using the numerical solution, and the contribution from the region r > r * is calculated analytically. At large distances, r > r * , we can use a linear theory of plasma oscillations because the initial velocity perturbation v r0 is small. In the linear theory, v r oscillates with the plasma frequency (cf. Eq. (24)), v r = v r0 cos ξ ≈ 2νK 1 (r) cos ξ, where we have used in the expression for the velocity the value of A given by Eq. (9), thus including the effect of the screening. Since the velocity v r is small, we can neglect the density perturbation and replace n by unity which gives for the contribution to E z where K 0 is the modified Bessel function of the second kind, γ E ≈ 0.577 is Euler's constant and we have used an asymptotic expression for K 0 in the limit r * 1. Hence we obtain We now introduce R = r * / √ ν; since r bm = 2.82 √ ν, choosing R 1 makes r * r bm . We then write where and The contribution E z1 is computed by solving the plasma flow in the region 0 < r < R √ ν and numerically calculating the integral in (34) using this solution. As follows from the derivation, E z1 actually does no depend on R: the change in the integral due to the variation of R is compensated by the logarithm in the second term in (34). Moreover, looking at the scalings (21) we conclude that E z1 is linearly proportional to ν. Hence, it can only be computed once, and then used for different values of ν by a simple rescaling. The second term, E z2 , is nonlinear in ν, but has an explicit simple analytical form.
To illustrate the calculation of E z we show in Fig. 6 the plots of E z1 /ν and E z /ν calculated with the help of Eqs. (33)-(35). Note that E z1 /ν does not depend on ν while E z /ν has a part (35) that logarithmically depends on ν; the plot of E z /ν in Fig. 6 corresponds to ν = 0.1. The black line shows E z calculated with R = 40, and the red dots indicate the result of the calculation of E z with R = 20. The agreement between these two calculations corroborate the statement that Eq. (34) uniquely defines E z1 independent of the value of R.
One can see from Fig. 6 that the field E z changes sign at ξ = 1.69, close to, but somewhat behind, the location of maximum value of the bubble radius. The sign of E z is such that electrons inside the bubble would be decelerated at ξ < 1.69 and accelerated in the region 1.69 < ξ < ξ b = 3.8.

VIII. POSITIVE DRIVER BEAM
Equations of Section V can be also used to simulate a positive driver beam of small dimensions with ν 1 by changing the sign of v r0 in Eq. (16). All the subsequent equations in that section remain valid, including the scaling of the variables (21) that eliminates ν from the equations.
We solved numerically the rescaled equations for the positive driver by the same method as for the electron driver-launching a large number of trajectories uniformly distributed over the initial radial position r 0 at ξ = 0 and then tracing them with account of the electric field of the plasma. An additional complication in numerical algorithm in comparison with the electron driver was caused by the fact that many trajectories cross the axis of the system generating a density singularity on the axis. In the tracking algorithm these trajectories were mirror reflected from the horizontal axis in the r, ξ plane. Calculated electron orbits in plasma are shown in Fig. 7. Again, as in the case of the electron driver, there are two electron streams at each point r, ξ with different values of density n and velocity v r . While there is no bubble in this case where electron density is equal to zero, one can see that there is a rarefaction in plasma trajectories in the region around ξ ≈ 5, r 1. This is more clearly seen in the radial density plots shown in Fig. 8 for several values of ξ. The plasma density in the region ξ ≈ 3 − 5 drops considerably below the equilibrium value n 0 = 1 before it reaches a spike on the axis.

IX. SUMMARY
In this paper we studied the plasma flow behind a relativistic electron bunch propagating through a cold plasma with the velocity of light. We assumed that the dimensions of the bunch are small so that the bunch can be treated as a point charge. In this model, the governing system of equations for the plasma contains only one parameter, the dimensionless charge of the driver ν, defined by Eq. (1). Assuming additionally that ν 1, we showed that this parameter can be eliminated from the equations, which then have a unique solution. In this solution, for an electron driver, a bubble is formed in the plasma from which electrons are evacuated. We showed that the first part of the boundary of the bubble, before it reaches its maximal radius, is an envelope of many electron trajectories each or which touches the boundary at one point. The electron density has a singularity at this part of the boundary. The second part of the boundary consists of a single electron trajectory, and the electron density is finite at the boundary. Simple analytical expressions were obtained in Section VI that describe the shape of the boundary with high precision. We found how the shape of the boundary scales with the driver charge Q: its maximal radius r bm scales as r bm ∝ √ Q, and its length does not depend on Q.
We also found the longitudinal electric field on the axis of the bubble and showed that it consists of two parts which have different scalings with Q: the first one is proportional to Q, and the second one scales as Q ln Q.
Finally, a numerical solution of the plasma equations was obtained for the case of a positively charged point driver. This solution does not have a bubble, but shows a region of rarefied electron density at the distance of approximately 5k −1 p behind the driver.
The two solutions correspond to two trajectories that arrive from different initial radii r 0 to a given point r, ξ. One of these trajectories arrives before and the other one after it touches the envelope. Correspondingly, at a given ξ, r, we need to sum the two densities for both trajectories. Substituting (A2) into (A1) we obtain n ± (r, ξ) = 1 2 For the total density, after some simplifications, we find n(r, ξ) = n + (r, ξ) + n − (r, ξ) = 1 2 This density has a square root singularity at the boundary of the bubble.