## 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, Reference Pukhov and Meyer-ter Vehn2002; Malka, Reference Malka2012; Kostyukov and Pukhov, Reference Kostyukov and Pukhov2015). 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* _{0} = *eA* _{0}/(*mc* ^{2}) > 4 and *S* = *n* _{e}/(*a* _{0}*n* _{c}) ≪ 1, where *a* _{0} 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, Reference Gordienko and Pukhov2005).

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, Reference Muggli2016). The bubble is surrounded by an electronic layer from which electrons can be trapped and focused into the bubble center (Kalmykov *et al*., Reference Kalmykov, Yi, Khudik and Shvets2009). 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*., Reference Faure, Rechatin, Norlin, Lifschitz, Glinec and Malka2006; Li *et al*., Reference Li, Hua, Xu, Zhang, Yan, Du, Huang, Chen, Tang, Lu, Joshi, Mori and Gu2013; Tooley *et al*., Reference Tooley, Ersfeld, Yoffe, Noble, Brunetti, Sheng, Islam and Jaroszynski2017). 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., Reference Baxevanis, Hogan, Huang, Litos, O'Shea, Raubenheimer, Frisch, White, Xu and Mori2017; Gonsalves *et al*., Reference Gonsalves, Pollock and Lu2017; Huang *et al*., Reference Huang, Zhou, Li, Wan, Wu, Hua, Pai, Lu, Wang, Deng, Liu, Wang, Zhao, An, Xu, Joshi and Mori2017; Tooley *et al*., Reference Tooley, Ersfeld, Yoffe, Noble, Brunetti, Sheng, Islam and Jaroszynski2017; Wang *et al*., Reference Wang, Feng, Zhu, Li, He, Li, Tan, Ma and Chen2018). The density-down ramp is achieved by longitudinally modulating the plasma density with extremely large gradients (Gonsalves *et al*., Reference Gonsalves, Nakamura, Lin, Panasenko, Shiraishi, Sokollik, Benedett, Schroeder, Geddes, van Tilborg, Osterhoff, Esarey, Toth and Leemans2011; Swanson *et al*., Reference Swanson, Tsai, Barber, Lehe, Mao, Steinke, van Tilborg, Nakamura, Geddes, Schroeder, Esarey and Leemans2017; Xu *et al*., Reference Xu, Li, An, Dalichaouch, Yu, Lu, Joshi and Mori2017*a*, Reference Xu, Buck, Chou, Schmid, Shen, Tajima, Kaluza and Veisz2017*b*). Ionization injection requires a small amount of higher-*Z* gas, added to the gas used for acceleration (Pak *et al*., Reference Pak, Marsh, Martins, Lu, Mori and Joshi2010; Tochitsky *et al*., Reference Tochitsky, Fiuza and Joshi2016). 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*., Reference Hidding, Pretzler, Rosenzweig, Königstein, Schiller and Bruhwiler2012*a*, Reference Hidding, Rosenzweig, Xi, OShea, Andonian, Schiller, Barber, Williams, Pretzler, Königstein, Kleeschulte, Hogan, Litos, Corde, White, Muggli, Bruhwiler and Lotov2012*b*). 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*., Reference Hidding, Pretzler, Rosenzweig, Königstein, Schiller and Bruhwiler2012*a*; Chen *et al*., Reference Chen, Esarey, Geddes, Cormier-Michel, Schroeder, Bulanov, Benedetti, Yu, Rykovanov, Bruhwiler and Leemans2014).

Some rough descriptions of the bunch's structure already have been made using shadowgraphy or X-ray betatron radiation (Schnell *et al*., Reference Schnell, Sävert, Landgraf, Reuter, Nicolai, Jäckel, Peth, Thiele, Jansen, Pukhov, Willi, Kaluza and Spielmann2012; Saevert *et al*., Reference Sävert, Mangles, Schnell, Siminos, Cole, Leier, Reuter, Schwab, Möller, Poder, Jäckel, Paulus, Spielmann, Skupin, Najmudin and Kaluza2015). 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, Reference Apostol and Ganciu2011; Petrillo *et al*., Reference Petrillo, Bacci, Zinati, Chaikovska, Curatolo, Ferrario, Maroli, Ronsivalle, Rossi, Serafini, Tomassini, Vaccarezza and Variola2012).

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*. (Reference Thomas, Günther and Pukhov2017) 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, Reference Wigner1934; Crandall and Williams, Reference Crandall and Williams1971; Meissner *et al*., Reference Meissner, Namaizawa and Voss1976; Dubin and O'Neil, Reference Dubin and O'Neil1999; Morfill and Ivlev, Reference Morfill and Ivlev2009; Radzvilavičius and Anisimovas, Reference Radzvilavičius and Anisimovas2011). The advantage of using a Taylor expansion is that the calculation of the implicitly given retarded time

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*., Reference Reichwein, Thomas and Pukhov2018). 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*. (Reference Thomas, Günther and Pukhov2017). In the following section of our paper, we will cover the Lorentz transformed fields which use the approach of Jackson *et al*. (Reference Jackson, Witte and Diestelhorst2013). 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*. (Reference Kostyukov, Pukhov and Kiselev2004). Here, relativistic electrons are accelerated by the normalized force *F* _{z} = −(1 + *V* _{0})ξ/4, where ξ = *z* − *V* _{0}*t* is the longitudinal position inside of the moving bubble with the velocity *V* _{0}. 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

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*. (Reference Thomas, Günther and Pukhov2017), 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

and the magnetic field

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*., Reference Jackson, Witte and Diestelhorst2013). 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 γ^{−2} since the sine term vanishes, whereas for ψ = ±π/2, the field is stronger by a factor of γ.

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

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

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

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*. (Reference Thomas, Günther and Pukhov2017) and Reichwein *et al*. (Reference Reichwein, Thomas and Pukhov2018) 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

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

Therefore, we can rewrite the total force as follows:

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*. (Reference Reichwein, Thomas and Pukhov2018) for the transverse direction, giving us

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

with the normalization of Reichwein *et al*. (Reference Reichwein, Thomas and Pukhov2018). Here, the variable *r* _{e} = 2π*e* ^{2}/(*mc*)^{2} represents the classical electron radius. Therefore, the scaling of the inter-particle distance in the propagation direction is

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:

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

according to Barzilai and Borwein (Reference Barzilai and Borwein1988) 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, Reference Hamacher and Wenzel1999; Metropolis *et al*., Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953), 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).

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.

## 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,

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.

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

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

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

which is in excellent agreement with the ESM.

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, Reference James1998). 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).

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*. (Reference Pyka, Keller, Partner, Nigmatullin, Burgermeister, Meier, Kuhlmann, Retzker, Plenio, del Campo, Zurek and Mehlstaeubler2013). 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, Reference Schiffer1996; Ikegami *et al*., Reference Ikegami, Okamoto and Yuri2006). 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

Thus, for a plasma wavelength of λ_{p} = 10^{5} nm, we get *n* _{ion} ≈ 1.1 × 10^{23} m^{−3}. Looking at the volume of the structure for 20,000 electrons in Figure 5, which roughly amounts to 32 nm^{3}, we would expect a total of approximately 3.52 × 10^{−3} 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*. (Reference Thomas, Günther and Pukhov2017), 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*. (Reference Thomas, Günther and Pukhov2017). 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*., Reference Reichwein, Thomas and Pukhov2018), 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*. (Reference Thomas, Günther and Pukhov2017), 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, Reference Schiffer1996; Ikegami *et al*., Reference Ikegami, Okamoto and Yuri2006). 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.

## Acknowledgments

This work has been supported in parts by DFG (project PU 213/6-1) and by BMBF (project 05K16PFB).