Microstructure formation in radially counterstreaming electron flows

A theoretical model is developed to describe the formation of microstructures due to plasma streaming instabilities in radially convergent geometries. Microstructures in the form of radial spokes are found experimentally in laser wakefield accelerators. The eigenvalues of a set of coupled linear ordinary differential equations are obtained and the complex wavenumbers calculated to give the local growth rates. The predictions are confirmed using particle-in-cell (PIC) simulations carried out for two counter-propagating converging/diverging plasma annuli. The simulations consistently demonstrate unstable growth for interactions between two counterpropagating annuli with different number densities. The growth rates obtained from the PIC simulations agree well with the growth rates from the semi-analytical model. The theory presented in this paper can provide powerful insight into converging plasma beams found in space and laboratory scale plasma.


Introduction
Converging and diverging flows of charged particles are ubiquitous in physics, from extreme astrophysical events [1], such as supernovae [2], gamma ray bursts [3,4] and coronal mass ejections [5][6][7], to emerging technologies, such as inertial confinement fusion [4] and laser wakefield accelerators (LWFAs) [8][9][10][11]. The LWFA accelerating structure comprises an evacuated bubble of ions surrounded by a high-density sheath of electrons that stream backwards around the bubble perimeter before converging and crossing at its rear. While a small fraction of the sheath electrons are injected into the back of the bubble and accelerate, most converge and stream through the confined crossing region at its rear to form the sheath of subsequent bubbles or are lost as so called 'side electrons' [12]. The inevitable crossing of converging/diverging annular particle bunches sets up ideal conditions for strong mutual interactions between counterstreaming electrons. The effect of these interactions can be observed as quasi-regular modulations of the sheath electron momenta and density of the side electrons emitted from the accelerating structure or those forming the sheath of subsequent bubbles. Converging and diverging electrons with opposite transverse momenta will interact with each other as they counterstream, which provides ideal conditions for streaming instabilities to develop (see figures 3 and 4). Recent experiments reveal periodic striations or density 'rays' imprinted on the side-electron beams [12], which we show, is evidence of streaming instabilities in converging-diverging geometries. In this paper, we develop a model that explains the origins of the microstructure evident in the electron streams/jets ejected from the LWFA bubble's sheath after crossing. The background plasma density is estimated to be two orders of magnitude less than that of the sheath [13], which suggests that the influence of the background plasma is negligible.
Instabilities can develop when counterstreaming particle beams interact with each other. Small perturbations lead to out-of-phase oscillations and drive strong mutual interactions and bunching of the particles [4,14,15]. In 1949, Bohm and Gross [16] showed that counterstreaming charged particles beams can lead to longitudinal density modulations due to an instability known as the two-stream instability (TSI). Weibel [17] later showed that modulations can also develop perpendicular to the flow direction because of thermal anisotropies in plasma. In the same year, Fried [18] demonstrated that the physical mechanism behind the Weibel instability could lead to wave growth transverse to the direction of counterstreaming cold beams. Fried's two beam approach is commonly referred to as the current filamentation instability (CFI). Both TSI and CFI require the mutual interaction between the two counterstreaming particle beams, but CFI usually excites a larger range of wavenumbers than TSI [4,[19][20][21]. Previous studies of instabilities have mostly assumed slab geometry and homogeneous particle streams. However, this does not provide an adequate description of instabilities for many geometries, including those usually found in inertial confinement fusion, astrophysical environments and in the LWFA.

Theory
To model the mutual interaction between converging/diverging beams we simplify the geometry by considering two streams that move radially inwards and outwards, respectively. These represent particles just before and just after crossing at the centre. Gratton and Gnavi have demonstrated TSI for a convergent geometry [22], and Genoni et al [23] have investigated TSI for spherical geometry, where particles converge, cross at a point, and then diverge. In this paper we study CF and oblique instabilities, where the wavevector is not aligned with the flow direction. To obtain a comprehensive understanding of the development of instabilities in this geometry we develop a semi-analytical theory based on a fluid description, and compare its predictions with a numerical modelling using particle-in-cell (PIC) simulations.
The fluid description is based on the continuity and momentum equations applied to converging and diverging electron flows, and Maxwell's equations which self-consistently model the electromagnetic fields: where the ± subscript denotes converging/diverging electrons, n is number density, u is velocity, p is momentum, and e is the electron charge. We proceed by linearising the equations for small perturbations about an equilibrium where two electron populations flow radially inward and outward, respectively. This leads to a set of coupled linear partial differential equations that can be solved as an eigenvalue problem. In homogeneous plasma, if the real and imaginary parts of the eigenvalues exhibit coalescence and bifurcation, then instabilities are present. We would therefore expect the same to occur in inhomogeneous plasma. We assume a cold, relativistic, inhomogeneous plasma, where two electron populations counterstream relative to each other with equal and opposite velocities,ū + = −ū − =ūe r . We denote equilibrium fields by a bar. Particle species are labelled ± corresponding to inward or outward propagation, respectively. Assuming no temporal and azimuthal dependence in the equilibrium densities and velocities, the continuity equation implies that the productnūr is constant, wheren is the equilibrium density, and r is the distance from the centre. For constantū,n is proportional to 1/r. The singularity as r → 0 can be avoided by assuming a small thermal spread in the beams. This central region (r ≈ 0) is ignored in our theory.
The perturbations of the velocities, fields, and densities can be assumed to have a harmonic dependence on the azimuthal angle θ and time t. While in principle we could seek a complex frequency to give the instability growth rate, this would involve a rather intractable eigenvalue problem imposing outward propagating boundary conditions on a system with multiple wave modes. Instead, we consider the more tractable problem of calculating local wavenumbers k r (r), giving a spatial growth rate, with a fixed frequency. Using the WKB approximation, the perturbations (denoted by a tilde) have the form ψ(r, θ, t) = ψ 0 exp(i k r (r )dr + i θ − iωt), where ψ 0 is the amplitude, is the mode number, and ω is the frequency. From equations (1) to (6), one can derive a set of coupled linear ordinary differential equations for the perturbations of the number densities, radial and azimuthal velocities of inward and outward electron flows, and of the electromagnetic fields,Ẽ r ,Ẽ θ , andB z .
The equations are made dimensionless by scaling the position ξ = r/r 0 , with the radial scale length r 0 used as the position, wheren(r 0 ) = n 0 . Similarly, the dimensionless forms of the perturbations are given by: N ± =ñ ± /n 0 ; velocity U ±,i =ũ ±,i /ū; electric fields E i = er 0Ẽi /(mū 2 ); magnetic field B z = er 0Bz /(mū); frequency Ω = r 0 ω/ū; and plasma frequency Ω 2 p = e 2 n 0 r 2 0 /( 0 mū 2 ). Here, u is the velocity; the subscript i = r, θ; m is the electron mass; and e is the electron charge. The dimensionless coupled equations are: d dξ where β =ū/c and γ = (1 − β 2 ) −1/2 . These equations can be written as: where X is the 9D vector comprising the perturbations, X is its derivative with respect to ξ, and Λ is the matrix of coefficients multiplying the corresponding perturbations in equations (7) through (12). Radial wave numbers can be defined as k r r 0 = −iλ, where λ denotes the eigenvalues of matrix Λ, which can be found numerically. As r → 0, the number density diverges, rendering the WKB approximation inadequate at small r. In the spatially homogeneous case, stable regions are characterised by equal and opposite real wavenumber with zero imaginary part. The transition to an unstable region arises when the real values coalesce to zero and equal and opposite imaginary parts appear. In the inhomogeneous case the situation is complicated by the fact that there are spatial variations in wave amplitude due to geometrical effects in addition to the instability. There are regions where there are equal and opposite real parts and identical imaginary parts. These can be interpreted as representing oppositely propagating stable waves with identical amplitude variation resulting from the geometry. Then, there are regions where the real parts coalesce, though generally to a non-zero value, while the imaginary parts bifurcate into unequal values. Here the difference in the spatial variation comes from the existence of an instability, which adds spatial growth to the geometrical variation. The TSI can be described by eliminating azimuthal dependences in equations (7)- (12), where perturbations grow in the velocity vector direction. We find no coalescence or bifurcations of real or imaginary parts of the wave numbers in this case, which indicates no instability. While this appears to disagree with Gratton and Gnavi [22], their calculation includes a singularity in the density at the origin. We believe this singularity may be unrealistic for our case and therefore only consider regions excluding the origin, which may be a reason for this disagreement.
Non-zero azimuthal wave-modes imply oblique instabilities because perturbations vary both along and perpendicular to the direction of propagation, i.e. the corresponding wave vector is oblique to the direction of flow. Oblique instabilities combine aspects of both two-stream and current filamentation instabilities, and maximum linear growth can occur at a particular angle [24], although in the nonlinear evolution the interplay of modes with similar growth rates and propagation angles can be more important than a single 'dominant' mode [25]. Figure 1 shows the set of wave numbers k r , with real and imaginary parts in the left and right plots, respectively. Branches with Re(k r ) > 0 correspond to outgoing waves, branches with Re(k r ) < 0 to incoming waves, Im(k r ) > 0 to amplitude decays (for increasing ξ), and Im(k r ) < 0 to amplitude growth. From figure 1 one can identify a point where, as ξ increases, imaginary parts coalesce, while real parts bifurcate, demonstrating evidence of an instability. The range (shaded blue between the imaginary parts) to the left of this point is identified as a region of growth due to an instability.
Besides the spatial growth (for changing r), the local temporal growth (for increasing time at fixed r) is of interest. The WKB analysis of the spatial dependence for a given time dependence indicates that an instability can be expected, but cannot be related directly to the numerical simulations we describe below, where there is time variation in the amplitudes and the beams interact over a limited spatial region. However, it is possible to relate spatial and temporal growth rates in homogeneous plasma [26] therefore we  look at how these rates are related locally at various points in the inhomogeneous case to obtain an indication of the magnitude of the temporal growth rate. By varying the imaginary part of the frequency and mapping the real against the imaginary parts of the wavenumbers at a fixed position, the maximum temporal growth rate is obtained from a point of coalescence and bifurcation of two different branches [26,27]. Due to the spatial density profile, the temporal growth rate inherits a spatial dependence, as shown in figure 2. For bunches of finite length, this spatial dependence corresponds to the position where they interact.
The growth rate decreases according to a power law as ξ increases. The coalescence/bifurcation of real and imaginary parts of the wavenumbers lie on the imaginary k axis, which implies that the modulations are purely azimuthal and correspond to the CFI.

Simulations
To verify the analytic studies described above we have undertaken PIC simulations, using the EPOCH [28] code, replicating the convergent geometry at the back of an LWFA bubble, with a planar domain representing a section of the sheath crossing region (transverse to the plasma bubble propagation). The cells have a width and length of 30 nm and each contain 40 particles; and simple outflow boundaries have been used. To avoid the singularity at r = 0, ingoing and outgoing electron bunches are shaped into annuli, representing a segment of the converging/diverging electron streams. The annuli counter-propagate radially and interact where they pass through each other. To replicate the conditions at the back of the LWFA bubble, the bunches are given initial radial velocities corresponding to a Lorentz factor γ = 2 and initial densities ofn + (r 0 ) = 10 19 cm −3 (inner) andn − = A +n+ /A − (outer), where r 0 = 1.7 μm and A + , A − are the respective cross sections of the bunches, so their electron number densities are equal when they fully  overlap [13]. They have equal initial widths of 2.6 μm, and are separated by 5.2 μm. A sinusoidal azimuthal density modulation of 20% is added to the outer annulus to seed the instability. Mode number = 10 was found to yield the highest growth rate for the densities and resolution used. To avoid expansion of the bunches due to space charge effects, a co-propagating ion bunch is added to each electron bunch. The initial electric and magnetic fields of the modulated bunch are explicitly included in the simulations to mitigate artificial growth.
The electron number density is radially integrated to obtain the angular density, which is then normalized to the total number of electrons in each bunch. The spectrum (spectral amplitude) is then calculated using a Fourier transform of the angular density. Initially, while the bunches have no spatial overlap, interaction is negligible. The spectrum shows a peak at = 10 for bunch 1 and noise for bunch 2, as shown in figure 3. The series of small peaks present in the angular number density of bunch 2 are artefacts arising from the Cartesian grid. The increase and decrease of the average number density of both bunches, approximately proportional to 1/r, is a geometry effect. The overlap and interact of the bunches results in a change in the spectrum. The effects of interaction are evident when the bunches fully overlap, and corresponding modulation, with the same mode number develops in the initially unmodulated annulus, but π out of phase, as shown in figure 4. Figure 5 shows the evolution of the signal of both bunches with and without interaction. The bunch 1 signal does not vary until later in the simulation, where space charge effects begin to dilute the modulation depth. Soon after the beginning of the interaction, bunch 1 experiences a slowing of the rate of decrease in signal until the bunches begin to fully overlap when compared to the non-interacting case. At this time, bunch 2 experiences signal growth and a modulation due to the instability that is one fifth of the initial modulation of bunch 1, but at a rate that decreases as the bunches move away from each other and the overlap region diminishes. The temperatures of the annular bunches are calculated using expressions from Schroeder and Esarey [29]; as their averages do not exceed ∼30 meV the bunches can still be considered cold, which implies that a coupling between 'pure' Weibel and current filamentation instabilities [30] is unlikely.

Discussion
To compare with the growth rate predicted by the semi-analytical theory, we estimate the growth rate from the simulations using the simplifying assumption that the modulations in both bunches mutually drive each other. Adding the respective differential equations yields where Γ is the growth rate and S is the signal of the bunch at mode = 10. Figure 6 shows the growth inferred from the simulations compared with the growth predicted by the semi-analytical approach for n(1.7 μm) = 10 19 cm −3 , and the local growth rate , with A = c 2 k 2 + 2ω 2 p /γ 3 , derived for a cold relativistic homogeneous plasma [31], but where the plasma frequency and wavenumber are replaced by their local values, ω p = e 2n (r)/(m 0 ) and k = 2π /r, respectively. The growth rates at different locations are determined by performing further simulations at different complete overlapping positions, while scaling the initial number density according to the 1/r dependence. They fit the expression Γ sim ∝ ξ −0.88 , which should be compared with the local growth rate Γ local ∝ ξ −0.75 and the semi-analytical expression using r 0 = 5.6 μm Γ s−a ∝ ξ −0.70 , as shown in figure 6. This estimate of temporal growth rates is obtained from a mapping of spatial to temporal growth, which is valid in a homogeneous system. Since spatial inhomogeneity is expected to reduce the growth rate, this should be regarded as an upper limit. Furthermore, the theoretical maximum growth rate corresponds to a continuous particle stream, while the simulations consider two finite propagating particle bunches, which contributes to discrepancies between the two growth rates. Also, the model underlying equation (14) does not take inertial effects in the interaction into account, whereas figure 5 shows a delay between overlap of the bunches and the response. Nonetheless, the growth rates determined from simulations and analytic theory agree reasonably, to within an order of magnitude, and exhibit similar scaling dependencies. Figure 6. The predicted local maximum temporal growth rate (dashed) compared with the growth rates of a simulation of full overlap at r 0 = 5.6 μm (solid), the maximum growth rates of simulations with overlap at different positions (triangles) and the growth rate from the homogeneous approach for varying number density (dots). The dot dashed line represents a fit to the growth rates from simulations at different complete overlap positions.

Conclusion
We have developed a semi-analytical theory to explain microstructure formation by instabilities arising from the interaction of two counterstreaming plasmas in a convergent, disk-like geometry. We confirm that the CFI occurs for a radial convergent geometry, and have obtained its temporal growth rate. This has been verified by PIC simulations, which demonstrate microstructure formation due to the CFI in the same radially convergent geometry. The temporal growth rates determined from the simulations compare well with our analytical model. The inevitable crossing at the back of the wakefield bubble gives rise to a fine structure, which is observed as azimuthal modulations. The modulations arising from these instabilities are observable in experiments investigating electrons ejected from the bubble's sheath after crossing [12]. In addition, counterstreaming and convergent beams of charged particles are also found in astrophysical scenarios [1] such as in gamma ray bursts [3,4], supernovae [2] and coronal mass ejections [5][6][7], as well as in laboratory plasmas such as inertial confinement fusion [32,33].