The filamented electron bunch of the bubble regime

Abstract

We present a theory for describing the inner structure of the electron bunch in the bubble regime starting from a random distribution of electrons inside the bubble and subsequently minimizing the system's energy. Consequently, we find a filament-like structure in the direction of propagation that is surrounded by various shells consisting of further electrons. If we specify a two-dimensional (2D) initial structure, we observe a hexagonal structure for a high number of particles, corresponding to the close packing of spheres in two dimensions. The 2D structures are in agreement with the equilibrium slice model.

Introduction

In contrast to conventional accelerators, plasma-based acceleration has the advantage of attaining higher particle energies over shorter acceleration distances due to the higher electric field that can be applied (Pukhov and Meyer-ter Vehn, 2002; Malka, 2012; Kostyukov and Pukhov, 2015). In laser-driven plasma wake-field acceleration (LWFA), a highly non-linear broken wave regime in form of an electronic plasma cavity called the “bubble regime” can occur. The bubble regime arises for a ₀ = eA ₀/(mc ²) > 4 and S = n _e/(a ₀n _c) ≪ 1, where a ₀ is the normalized amplitude of the laser vector potential and S is the similarity parameter. Here, e is the elementary charge, m is the mass of the electron, c is the speed of light, n _e is the electron density, and n _c is the critical density (Gordienko and Pukhov, 2005).

The bubble potential is a nearly harmonic potential with electric fields of more than 100 GV/m. In general, plasma wakefields similar to the bubble can also be excited by dense and high energetic particle beams (Muggli, 2016). The bubble is surrounded by an electronic layer from which electrons can be trapped and focused into the bubble center (Kalmykov et al., 2009). The trapped electrons form the so-called electron bunch (or beam load).

Besides this self-trapping mechanism, there exist several other injection methods, including pre-acceleration, ionization, and density modulation techniques (Faure et al., 2006; Li et al., 2013; Tooley et al., 2017). In all cases, the objective is to create electron bunches with as small beam emittances as possible. The currently most promising methods are the ionization injection and the density down-ramp. Both methods produce electron beams with sub-fs duration, high peak currents in the range of several kA, energy spreads well below 1% and excellent transverse emittances (Baxevanis et al., 2017; Gonsalves et al., 2017; Huang et al., 2017; Tooley et al., 2017; Wang et al., 2018). The density-down ramp is achieved by longitudinally modulating the plasma density with extremely large gradients (Gonsalves et al., 2011; Swanson et al., 2017; Xu et al., 2017a, 2017b). Ionization injection requires a small amount of higher-Z gas, added to the gas used for acceleration (Pak et al., 2010; Tochitsky et al., 2016). In the case of the wakefield being driven by a short electron beam, the Trojan horse regime of underdense photocathode plasma wake-field acceleration is reached (Hidding et al., 2012a, 2012b). It can be used to decouple the electron bunch generation process from the excitation of the accelerating plasma cavity. The combination of the non-relativistic intensities required for tunnel ionization, a localized release volume as small as the laser focus, the greatly minimized transverse momenta, and the rapid acceleration leads to dense phase-space packets. In homogeneous plasma, they can have ultra-low normalized transverse emittance in the bulk of μm mrad and a minimal energy spread in the 0.1% range (Hidding et al., 2012a; Chen et al., 2014).

Some rough descriptions of the bunch's structure already have been made using shadowgraphy or X-ray betatron radiation (Schnell et al., 2012; Saevert et al., 2015). We, however, are interested in the finer sub-structure of the bunch, which is interesting for the field of short wavelength radiation. Here, a counter-propagating laser pulse is scattered back by a relativistic electron bunch, such that spatially incoherent photons of a wide energy spectrum are obtained. If the electron bunch exhibited a regular sub-structure, higher brightness and spatial coherence could be achieved (Apostol and Ganciu, 2011; Petrillo et al., 2012).

The are two approaches to describe the bunch's structure that calculate the prevalent fields in different ways. Originally, the sub-structure was described in Thomas et al. (2017) using a Taylor expansion of the retarded Liénard–Wiechert potentials up to second order in v/c. There, electronic filaments along the propagation direction of the bubble and hexagonal lattices in the transversal were observed. These regular electron structures are similar to Wigner crystals, known from other areas of plasma physics than wake-field acceleration (Wigner, 1934; Crandall and Williams, 1971; Meissner et al., 1976; Dubin and O'Neil, 1999; Morfill and Ivlev, 2009; Radzvilavičius and Anisimovas, 2011). The advantage of using a Taylor expansion is that the calculation of the implicitly given retarded time

(1)$$t_{\rm ret} = t - {1\over c}\vert r_i - r_j\lpar t_{\rm ret}\rpar \vert$$

can be circumvented. Here, index j indicates the radiation of a signal at time t _ret from position r _j, while index i denotes an observer at time t and position r _i receiving the signal. This approach yields incorrect inter-particle distances, as radiative terms are neglected. In the frame of this theory, a phase transition was observed: For emittances below a certain critical value depending on the system parameters, the crystalline structure persists in dynamical simulations. If the threshold is exceeded, the structure becomes a degenerate electron fluid. In the second ansatz, the equilibrium slice model (ESM) uses the full Liénard–Wiechert potentials but only examines two-dimensional (2D) slices transverse to the direction of propagation (Reichwein et al., 2018). This approach leads to hexagonal lattices as well, however with different inter-particle distances since the full Liénard–Wiechert potentials are taken into account. The scaling of these distances regarding particle momentum and plasma wavelength were explained analytically by a heuristic two-particle model. Contrary to the approach via the Taylor expansion, the ESM is restrained to only two-spatial dimensions and a static description of the bunch as calculating the dynamics would require to save the particles' history making it computationally expensive.

In the present paper, we derive a new model for the three-dimensional (3D) structure of the bunch in the static case using a Lorentz transformation of the electromagnetic fields under the assumption that the velocity of the particles is constant. This allows us to avoid the calculation of the retarded times while still describing the structure of the bunch with more precision than in Thomas et al. (2017). In the following section of our paper, we will cover the Lorentz transformed fields which use the approach of Jackson et al. (2013). Using the terms for the focusing force of the bubble potential and the ones for the repulsive Coulomb interaction between the electrons, we can formulate an equilibrium state. This state will represent the structure the system will want to attain. In Section "The 3D equilibrium state", we will cover the numerical algorithm for minimizing the total force of the system and the choice of the step size. Further, we will discuss the dependencies of the mean inter-particle distance since the propagation direction will show different scaling than the transverse direction due to the different strength of the electromagnetic fields in different directions. Finally, we will present the results of our simulations and particularly discuss the scaling regarding the total number of electrons.

Mathematical model and scaling laws

In the following, we derive the 3D inter-particle force in a system of interacting alongside propagating relativistic electrons in an external bubble potential in a moving coordinate system. For the potential, we choose the strongly simplified quasi-static 3D bubble model for electron acceleration in homogeneous plasma from Kostyukov et al. (2004). Here, relativistic electrons are accelerated by the normalized force F _z = −(1 + V ₀)ξ/4, where ξ = z − V ₀t is the longitudinal position inside of the moving bubble with the velocity V ₀. The focusing to the ξ-axis is provided by the force $F_r=-\lpar p_z+{\rm \gamma }\rpar \sqrt {x^2+y^2}/\lpar 4{\rm \gamma} \rpar$. Due to the cylinder symmetric form of the bubble potential, the electrons are also focused in the ξ-direction, namely to the bubble center, where they have maximum energy. Since

(2)$${dv\over dt} \approx {F\over 2m {\rm \gamma}^3}$$

for relativistic particles, we have $\vert \dot {v}\vert \approx 0$ and $\dot {{\rm \gamma }}\approx 0$ in a sufficiently small domain around the bubble origin.

From the previous work in Thomas et al. (2017), we can expect that the equilibrium configuration will be located in or at least near the bubble center. If we want to circumvent computing a Taylor series of the retarded electromagnetic potentials, we have to find a different way evaluating retarded times. In principle, this is impossible for accelerated particles. Thus, we use the sensible approximation $\vert \dot {v}_i\vert =0$ and v _i = v for all particles i, allowing for a Lorentz transformation of the electromagnetic fields from the rest frame to the laboratory frame. In this case, the Coulomb interaction between all electrons is determined by the electric field

(3)$$E = {q r\over r^2 {\rm \gamma} \lpar 1- {\rm \beta}^2 \sin\lpar {\rm \psi}\rpar \rpar ^{3/2}}$$

and the magnetic field

(4)$$B = {\rm \beta} \times E,$$

where β = |β| is the velocity v normalized with the speed of light c and n is the unit vector pointing from the charge moving, respectively, to the resting observer (Jackson et al., 2013). Further, we have ${\rm \psi }=\arccos \lpar {\bf n} \cdot \hat {v}\rpar$ with $\hat {v}=v/\vert v\vert$ (Fig. 1). For all cases but β = 0, we see that the electric field is anisotropic. For angles ψ = 0, π, we see a weaker field by a factor of γ⁻² since the sine term vanishes, whereas for ψ = ±π/2, the field is stronger by a factor of γ.

Fig. 1. Depiction of the setting for the Lorentz transformation: the charged particle moves along the x′₁ axis and is seen by an observer at position P (Jackson et al., 2013).

Having calculated the electromagnetic fields, we can calculate the forces affecting a given particle. The total force onto the ith electron is

(5)$$F_i = F_{{\rm ext}\comma i} + F_{{\rm C}\comma i}\comma $$

where F _ext,i is the force exerted by the bubble potential and F _C,i is the sum over all Coulomb forces between electrons i and j, such that

(6)$$F_{{\rm C}\comma i} = \sum_{i=1}^{N} F_{{\rm C}\comma ij}.$$

The forces can be calculated by using the electromagnetic fields and the equation for the Lorentz force

(7)$$F_{\rm L} = q\lpar E + {\rm \beta} \times B\rpar .$$

Minimizing the total forces F _i for all particles i in the system, we find its energetic minimum and thereby the corresponding structure. The equivalence of a Hamiltonian approach as in Thomas et al. (2017) and Reichwein et al. (2018) to our force balance here can be seen by writing the forces as the gradient of the Lagrangian L from the references above such that

(8)$${d\over dt}p_i = {d\over dt} \nabla_{{v}_i} {\rm L} = \nabla_{{r}_i} {\rm L}.$$

If we split up the momentum and the Lagrangian into external and Coulombic parts in a similar fashion to the forces, we have

(9)$${d\over dt}{\,p}_{{\rm ext}\comma i} = \nabla_{{r}_i} {\rm L}_{{\rm ext}\comma i} + \nabla_{{r}_i} {\rm L}_{{\rm C}\comma i} - {d\over dt}{\,p}_{{\rm C}\comma i}\comma $$

(10)$${d\over dt}{\,p}_{{\rm C}\comma i} = \nabla_{{v}_i} {\rm L}_{{\rm ext}}.$$

Therefore, we can rewrite the total force as follows:

(11)$$F_i = {d\over dt}{\,p}_{{\rm ext}\comma i} - {q_i\over c} {d\over dt}{A}_{{r}_i}+\left(\nabla_{{r}_i} - {d\over dt} \nabla_{{v}_i}\right){\rm L}_{\rm C}\comma $$

finally leading to our equation for the force (5), showing that the force balance is an equivalent way of describing the energy minimization.

Before conducting the simulations, we already can estimate the scaling of the mean inter-particle distance regarding the parameters momentum p and plasma wavelength λ_pe. For this, we are able to use the heuristic approach of Reichwein et al. (2018) for the transverse direction, giving us

(12)$$\Delta r = \sqrt [3] {{r_e\over 2{\rm \pi}^3}}\left({{\rm \lambda}_{{\rm pe}}\over \sqrt{{\rm \gamma}}}\right)^{2/3}\propto p^{-1/3}{\rm \lambda}_{{\rm pe}}^{2/3}$$

for the distance between two nearest neighbors in the 2D lattice. Due to the structure of the fields, the particles sense a 1/γ times weaker interaction force in the propagation direction, such that for a balance of the bubble force with the interaction term, we have

(13)$${\Delta {\rm \xi}\over 4}= {r_e\over {\rm \lambda}_{{\rm pe}}}{1\over {\rm \gamma}^2}{1\over \lpar \Delta {\rm \xi}\rpar ^2}$$

with the normalization of Reichwein et al. (2018). Here, the variable r _e = 2πe ²/(mc)² represents the classical electron radius. Therefore, the scaling of the inter-particle distance in the propagation direction is

(14)$$\Delta {\rm \xi} \propto p^{-2/3} {\rm \lambda}_{{\rm pe}}^{2/3}.$$

The heuristic analytical model cannot explain the scaling regarding the total number of particles N at the moment. We will, however, look at this dependency numerically and give some ideas to what influences this behavior in the Discussion section.

The 3D equilibrium state

In order to find the equilibrium structure of the electron bunch, we use the so-called steepest descent method. At a given position ${\bf X}^k=\lpar {r}_1^k\comma\; \ldots \comma\; {r}_N^k\rpar$, we calculate the gradient $\nabla f\lpar {\bf X}^k\rpar$ of the function f(X) that is to be minimized (in our case, the magnitude of the total forces F _i). The gradient always points in the direction of steepest ascent, so going in the opposite direction brings us closer to the structure with minimal energy where $\lpar \nabla_{\bf X} f\rpar \lsqb {\bf X}_0\rsqb = 0$. Therefore, we can write our iterative algorithm as follows:

(15)$${\bf X}^{k+1} = {\bf X}^k - {\rm \alpha}_k \nabla_k f\lpar {\bf X}^k\rpar \comma $$

where α_k is a parameter for the step size at iteration k. The choice of this step size is crucial for our algorithm to converge since large steps lead to jumping over the position of the minimum, while too small step sizes lead to slow convergence. In our case, the choice of

(16)$${\rm \alpha}_k = {\Delta x \cdot \Delta g\over \Delta g \cdot \Delta g}$$

according to Barzilai and Borwein (1988) is sufficient. Here, Δx and Δg are the differences in position and gradient between the iterations k and k − 1, respectively. Using the steepest descent method, we find a local minimum. Generally, this does not need to be a global minimum as well. However, using the technique of stochastic tunneling (Hamacher and Wenzel, 1999; Metropolis et al., 1953), we are able to show that the structures we obtain actually are the global minima of the system. In these methods, a random vector is added onto the position X^k, such that valleys in the potential landscape, that normally would have been hidden from the gradient descent method due to surrounding hills, can be reached.

We distribute a fixed number of N particles randomly inside a spherical volume with a given momentum p and plasma wavelength λ_pe. The main structure we observe is a central filament on the ξ-axis (Fig. 2). For a sufficiently high number of particles surrounding elliptic shells can form (Fig. 3b).

Fig. 2. Formation of the central filament for increasing gradient descent iterationsit with decreasing error err from left to right. Notice the particles twirling onto the ξ-axis. A similar behavior can be seen for too many particles that are being forced into one filament; they try to escape into the x–y plane.

Fig. 3. Cross section of the 3D equilibrium structure for N = 20,000 electrons. (a) Transverse cross section in the plane where ξ = 0. (b) Longitudinal cross section for x = 0.

In order to generate a strongly simplified 3D depiction of the electron distribution at hand, we need to classify the various shells and filaments. To do so, we look at the transverse cross section in Figure 3a and plot the corresponding radial density profile (Fig. 4). Then, after having fitted a multi-Gaussian (red curve) to the distribution, we define: A shell is the set of all electrons inside the $\sqrt {2}{\rm \sigma}$environment of one Gaussian distribution. A final visualization of the structure in Figure 3 is shown in Figure 5.

Fig. 4. Histogram of the number of particles in the final distribution depending on the radius R of the total distribution. The different peaks represent the occurring shells with a certain thickness that are fitted using a multi-Gaussian.

Fig. 5. Simplified schematic depiction of the resulting structure for N = 20,000 electrons. Notice the central filament surrounded by several ellipsoid shells. Depending on the number of particles, the main filament (here shown as a continuous red line) is broken up into little pieces, and some of its electrons are assigned to the surrounding shells.

Numerical scaling

For our first series of simulations, we vary the momentum between 50 and 500 MeV/c for N = 1000 electrons and a plasma with λ_pe = 100 μm. The resulting structure is a single filament in the propagation direction (Fig. 2). We obtain a dependence according to our previously calculated scaling laws, that is,

(17)$$\Delta {\rm \xi} \propto p^{-2/3}.$$

The inter-particle distances are in the region of some picometers in the longitudinal direction (see Fig. 6). These distances are well under the diameter of an atom. Thus, we need to consider if quantum effects could play any role in this regime, even though electrons can be considered as point-like particles. The ratio between inter-particle distance and de Broglie wavelength λ_dB gives some indication regarding that. The wavelength is given by ${\rm \lambda }_{{\rm dB}} = 2{\rm \pi} \hbar /\lpar pc\rpar$, where $\hbar$ is the reduced Planck constant and p = γmc is the relativistic momentum of the electron. For a Lorentz factor γ = 100, we obtain λ_dB ≈ 24 fm, which is orders of magnitude below our simulation results for the inter-particle distance. Therefore, we can neglect quantum effects here.

Fig. 6. Dependence of the mean inter-particle distance for a constant number of N = 1000 electrons and λ_pe = 10⁵ nm in the propagation direction (Δξ) and the transverse direction (Δr). The circles represent the simulation data, while the lines show the power fit.

Scaling regarding λ_pe (Fig. 7) is in agreement with our analytic results as well and therefore yields

(18)$$\Delta {\rm \xi} \propto {\rm \lambda}_{{\rm pe}}^{2/3}.$$

Fig. 7. Scaling of the mean inter-particle distance with the plasma wavelength λ_pe for N = 1000 electrons and p = 100 MeV/c. The circles represent the simulation data, while the lines show the power fit.

Lastly, we keep p and λ_pe fixed but vary the number of electrons for each simulation. For a sufficiently high number of particles, we observe additional shells surrounding the main filament we have seen before (Fig. 3). The central filament now is discontinuous as some of the electrons go to the various shells. The resulting dependence of the inter-particle distance on N is

(19)$$\Delta {\rm \xi} \propto N^{-0.75}\comma\; $$

which can be seen in Figure 8. We further specify an initial distribution on a 2D slice (such that ξ = constant for all electrons) and embedding it in the 3D model. As a result, the 2D structure persists and hexagonal lattices are observed. The scaling of the mean inter-particle distance Δr in the slice regarding the different parameters is given by

(20)$$\Delta r \propto p^{-1/3}{\rm \lambda}_{{\rm pe}}^{2/3}N^{-0.14}\comma $$

which is in excellent agreement with the ESM.

Fig. 8. Simulation results for λ_pe = 10⁵ nm and p = 100 MeV/c for a varying number of particles. The circles represent the simulation data, while the lines show the power fit.

The exponents for the dependency on N cannot be explained by our two-particle model, in fact scaling laws regarding the number of particles are a problem in further fields of physics (James, 1998). We can, however, explain the behavior phenomenologically to some extent.

Discussion

Considering the case of the 2D structure, for N = 2 we see a straight line for the equilibrium structure, and for N = 3, an equilateral triangle since the repelling Coulomb force causes the electrons to be apart as far as possible from each other (see Fig. 9). Opposing to that the parabolic bubble potential confines the electrons and focuses them to its center. The interplay of these to effects leads to the close packing of spheres in two dimensions with an hcp lattice. For N = 7, we have one particle in the middle surrounded by one full shell of six further electrons. Adding another particle to the densest packing, we break the symmetry, meaning that we now see different distances for the electrons while before, every electron had the same distance to one another. Only if a sufficient amount of further electrons are supplied, we can fill up the next shell, such that the maximum symmetry is restored. For higher shells, a lot more particles are needed than the six that make up the first shell (see the case for N = 20 in Fig. 9).

Fig. 9. Formation of the different shells for higher number of particles N. Notice the symmetry breaking for the transition N = 7 → N = 8.

If we move on to the 3D case, we now have two competing effects. At first, only the central filament is being filled for a low number of electrons due to the different strength of the electric field in the different spatial directions. Increasing the number of particles is accompanied by reducing the inter-particle distance. If this distance cannot be further reduced without sacrificing minimal energy, the filament starts to curl into a helix-like structure (similar to the modus operandi of the algorithm seen in Fig. 2). This is comparable to the 2D zigzag structure observed in Pyka et al. (2013). Even higher number of particles break up this structure; the main filament is broken up into various pieces that cannot be clearly assigned to one single shell. This is the transition to additional shells surrounding the one at the center: the structure still does not have enough particles to completely fill those new shells but starts to build up the hcp structure in the transverse direction. These two competing minima, one being the central filament structure, the other being the surrounding shells, leads to the scaling of Δξ ∝ N ^−0.75. The same structural behavior but at different length scales has also been observed in the field of circular accelerators with storage rings (Schiffer, 1996; Ikegami et al., 2006). The main differences here are the vastly different length scales of milimeters instead of our findings of some nanometers or even picometers. Furthermore, these two structures occur at different time scales due to their methods of acceleration.

In the considered density regime, we can neglect the effect of resting single ions onto the lattice structure since the ion density scales as

(21)$$n_{{\rm ion}} = {1\over {\rm \lambda}_{\,p}^2} {m\over {\rm \pi} e^2}.$$

Thus, for a plasma wavelength of λ_p = 10⁵ nm, we get n _ion ≈ 1.1 × 10²³ m⁻³. Looking at the volume of the structure for 20,000 electrons in Figure 5, which roughly amounts to 32 nm³, we would expect a total of approximately 3.52 × 10⁻³ ions to be inside this volume. For larger bunches, the higher number of ions inside the volume would be compensated by an even larger amount of electrons since the density ratio remains constant. Therefore, we are able to further neglect any larger effects by ions onto the lattice structure.

In comparison to the model of Thomas et al. (2017), where the Taylor expansion in v/c was used, we now see one intact filaments throughout the whole length of the structure, whereas then many short filaments could be seen. This is due to the incorporation of more relativistic effects in our model, also leading to smaller distances in the range of picometers rather than nanometers like in Thomas et al. (2017). We do however still see those hexagonal lattices, albeit with smaller inter-particle distances. As we have already seen in the ESM (Reichwein et al., 2018), this again is due to retardation effects. Instead of those additional filaments, we now observe the formation of shells that exhibit some patterning on their surface for a high precision of the steepest descent algorithm.

Conclusion

We have presented a theory for describing the 3D structure of the electron bunch in the bubble regime. The basis of our model is the Lorentz transformation of the electromagnetic fields, allowing us to avoid the calculation of implicit retarded times. Our model uses a quasi-static picture and considers the electron bunch in equilibrium being around the center of the bubble. The electrons used are perfectly mono-energetic. This approach leads to the observation of electronic filaments in the propagation direction of the bubble. For a low number of particles, instead of many fragmentary filaments as in Thomas et al. (2017), one main filament containing all electrons of the bunch can be seen. A higher number of particles leads to the breaking of this filament and finally various surrounding shells, similar to structures previously found in simulations for circular accelerators (Schiffer, 1996; Ikegami et al., 2006). The formation of these additional shells corresponds to the genesis of the outer shells of the ESM . Scaling laws regarding the dependence of the mean inter-particle distance on momentum and plasma wavelength are derived by a heuristic two-particle model. The distances in the sub-nanometer regime in the transverse direction or even tens of picometers in the propagation direction are smaller than previously observed due to the higher incorporation of relativistic effects. Since, however, the de Broglie wavelength is in the range of some femtometers, we can neglect quantum effects at this point in time.