# Crystallization dynamics of magnetic skyrmions in a frustrated itinerant
magnet

Kotaro Shimizu  Department of Applied Physics, The University of Tokyo, Tokyo
113-8656, Japan  Department of Physics, University of Virginia,
Charlottesville, Virginia 22904, USA  Gia-Wei Chern  Department of Physics,
University of Virginia, Charlottesville, Virginia 22904, USA

###### Abstract

We investigate the phase ordering kinetics of skyrmion lattice (SkL) in a
metallic magnet. The SkL can be viewed as a superposition of magnetic stripes
whose periods are determined by the quasi-nesting wave vectors of the
underlying Fermi surface. An effective magnetic Hamiltonian that describes the
electron-mediated spin-spin interaction is obtained for a two-dimensional s-d
model with the Rashba spin-orbit coupling. Large-scale Landau-Lifshitz-Gilbert
dynamics simulations based on the effective spin Hamiltonian reveal a two-
stage phase ordering of the SkL phase after a thermal quench. The initial fast
crystallization of skyrmions is followed by a slow relaxation dominated by the
annihilation dynamics of dislocations, which are topological defects of the
constituent magnetic stripe orders. The late-stage phase ordering also
exhibits a dynamical scaling symmetry. We further show that the annihilation
of dislocations follows a power-law time dependence with a logarithmic
correction that depends on magnetic fields. Implications of our results for
SkL phases in magnetic materials are also discussed.

Complex magnetic textures such as vortices and skyrmions are not only of great
fundamental interest in magnetism but also have important implications in the
emerging technology of spintronics [1, 2, 3, 4]. Both vortices and skyrmions
are nano-sized particle-like spin-textures characterized by nontrivial
topological invariants. The presence of such complex patterns in metallic
magnets could give rise to intriguing electronic and transport properties due
to a nontrivial Berry phase acquired by electrons when traversing over closed
loops of noncollinear or noncoplanar spins [5, 6, 7, 1]. The well-studied
topological Hall effects [8, 9, 10, 11, 12, 13] and topological Nernst effects
[14, 15, 16] are some of the representative examples. Also importantly, such
topological electronic responses in metallic magnets can be controlled via the
manipulation of magnetic textures.

In magnetic materials, skyrmions are often stabilized in the form of a
skyrmion lattice (SkL), which is a periodic array of such particle-like
topological spin-textures. Indeed, SkL as a spontaneous ground state was
already predicted in the pioneering work of Bogdanov and Yablonskii that later
triggered the enormous interest in magnetic skyrmions [17, 18]. SkLs have
since been reported in several chiral magnets such as MnSi [19] and other B20
compounds, as well as centrosymmetric materials [20]. While the picture of SkL
as an array of particle-like objects offers an intuitive framework to
understand certain structural and dynamical aspects of skyrmion phases [21],
recent studies have revealed the close connection between SkL and multiple-$Q$
magnetic ordering driven by partial nesting of the electron Fermi surface [22,
23, 24]. Indeed, this mechanism has been conjectured to be ubiquitous for
itinerant spin systems with a wide range of filling fractions.

Despite the huge interest and extensive research on magnetic SkL over the past
decades, the phase-ordering dynamics of SkLs remains an open subject.
Specifically, here one concerns the dynamical evolution and potential
universal behaviors of skyrmion crystallization when a magnet is quenched into
a skyrmion phase. It is known that the kinetics of phase ordering depends
crucially on topological defects of the symmetry-breaking phase. Several
super-universality classes of domain-growth laws have been established over
the years [25, 26]. The fact that both magnetic and translational symmetries
are broken in a skyrmion crystal indicates rich relaxational dynamics of SkL
phases which has yet to be systematically investigated. Understanding the
phase ordering of SkL is also crucial to the engineering and control of
skyrmion phases in real materials.

In this paper, we make an important step toward this goal by investigating the
crystallization dynamics of a SkL in a realistic model of chiral metallic
magnets. A minimum microscopic model for such itinerant spin systems is the
s-d Hamiltonian $\mathcal{H}_{\rm
sd}=\sum_{\mathbf{k},\sigma}\epsilon_{\mathbf{k},\sigma}c^{\dagger}_{\mathbf{k},\sigma}c^{\,}_{\mathbf{k},\sigma}-J_{\rm
sd}\sum_{i}\mathbf{S}_{i}\cdot
c^{\dagger}_{i\sigma}\bm{\sigma}_{\sigma,\sigma^{\prime}}c^{\,}_{i\sigma^{\prime}}$,
where the first term describes electron hopping on a lattice with $t_{ij}$ the
transfer integrals and the second term represents a local coupling between the
itinerant s-electron and the local moments $\mathbf{S}_{i}$ of d-electrons;
$J_{\rm sd}$ represents the coupling strength. Spin-orbit coupling (SOC)
within such single-band models can be described by either Rashba or
Dresselhaus hopping terms. Dynamical simulations based on such electron models
requires solving a disordered electron tight-binding Hamiltonian at every
time-step, which could be prohibitively expensive for large systems. Instead,
here we consider an effective spin Hamiltonian, similar in spirit to the
Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [27, 28, 29], with SOC
properly included [30]. At the leading second-order perturbation, the
effective Hamiltonian has the general form

$\displaystyle\mathcal{H}=-\sum_{ij}\mathbf{S}_{i}\cdot\mathsf{J}(\mathbf{r}_{i}-\mathbf{r}_{j})\cdot\mathbf{S}_{j}-\sum_{i}\mathbf{H}\cdot\mathbf{S}_{i}.$
(1)

Here we have included a Zeeman coupling to an external field ${\bf
H}=(0,0,H)$, and $\mathsf{J}(\mathbf{r})$ represents an effective $3\times 3$
interaction matrix between two spins separated by $\mathbf{r}$. Its Fourier
transform is given by

$\displaystyle\tilde{\mathsf{J}}(\mathbf{q})=\frac{J_{\rm
sd}^{2}}{N}\sum_{\mathbf{k}}\sum_{\sigma,\sigma^{\prime}}\chi^{0}_{\sigma\sigma^{\prime}}(\mathbf{k},\mathbf{q})\left[\mathsf{I}+\mathsf{F}^{\sigma\sigma^{\prime}}(\mathbf{k},\mathbf{q})\right],$
(2)

where
$\chi^{0}_{\sigma\sigma^{\prime}}=[f(\epsilon_{\mathbf{k},\sigma})-f(\epsilon_{\mathbf{k}+\mathbf{q},\sigma^{\prime}})]/(\epsilon_{\mathbf{k}+\mathbf{q},\sigma^{\prime}}-\epsilon_{\mathbf{k},\sigma})$
is the spin-dependent susceptibility, $\epsilon_{\mathbf{k},\sigma}$ is the
electron band energy, $f(\epsilon)$ is the Fermi-Dirac function, $\mathsf{I}$
is the identity matrix, and the dimensionless matrix
$\mathsf{F}^{\sigma\sigma^{\prime}}$ accounts for the anisotropic spin
interaction due to SOC [30, 31].

The interaction matrix $\tilde{\mathsf{J}}(\mathbf{q})$ is often dominated by
a few wave vectors $\mathbf{Q}_{\eta}$ when part of the electron Fermi surface
is connected by them, i.e.
$\epsilon_{\mathbf{k}+\mathbf{Q}_{\eta},\sigma^{\prime}}\approx\epsilon_{\mathbf{k},\sigma}$.
Such partial nesting of the Fermi surface has been shown to be a primary
mechanism for the stabilization of skyrmion or vortex lattices in metallic
magnets [23]. An effective real-space Hamiltonian is then given by the inverse
Fourier transformation of $\tilde{\mathsf{J}}(\mathbf{q})$, which likely can
only be done numerically for general dispersion relation
$\epsilon_{\mathbf{k},\sigma}$. As a first-order approximation, effective spin
Hamiltonians can be obtained by keeping only contribution from the nesting
wave vectors:
$\tilde{\mathsf{J}}(\mathbf{q})\approx\sum_{\eta}\tilde{\mathsf{J}}(\mathbf{Q}_{\eta})\delta(\mathbf{q}-\mathbf{Q}_{\eta})$.
This approach has been employed to investigate complex spin textures in the
ground state of itinerant magnets [23, 24, 32, 33, 34]. However, this approach
gives rise to unrealistic infinite-range spin interactions. A more realistic
analytical approach, while preserving the correct form of the anisotropic
interaction, is to replace the $\delta$-function by a peak-shape function
$h(\mathbf{q})$ of a finite width, giving rise to spin-spin interactions which
decay with distance in real space [35].

For concreteness, here we apply the procedure described above to a square-
lattice s-d model with a Rashba SOC. A square array of skyrmions can be
stabilized by partial nesting with two wave vectors $\mathbf{Q}_{1}=(Q,0)$ and
$\mathbf{Q}_{2}=(0,Q)$. The resultant real-space interaction matrix is given
by

$\displaystyle\mathsf{J}(\mathbf{r})=\sum_{\eta=1,2}g(\mathbf{r})\left[{\rm
Re}\mathsf{J}_{\eta}\cos(\mathbf{Q}_{\eta}\\!\cdot\\!\mathbf{r})-{\rm
Im}\mathsf{J}_{\eta}\sin(\mathbf{Q}_{\eta}\\!\cdot\\!\mathbf{r})\right],\quad$
(3)

where

$\displaystyle\mathsf{J}_{1}=\left(\begin{array}[]{ccc}J^{\perp}&0&0\\\
0&J^{\perp}&-iD\\\
0&iD&J^{zz}\end{array}\right)\\!,\,\,\mathsf{J}_{2}=\left(\begin{array}[]{ccc}J^{\perp}&0&-iD\\\
0&J^{\perp}&0\\\ iD&0&J^{zz}\end{array}\right)\\!.\quad$ (10)

and $J^{\perp}$, $J^{zz}$ describe electron-mediated exchange interactions,
and $D$ represents effective long-ranged Dzyaloshinskii-Moriya interaction
(DMI) induced by SOC; both $J$ and $D$ are of the order of $J_{\rm sd}^{2}/W$,
where $W$ is the electron bandwidth. In the following, we set
$2J^{\perp}+J^{z}=1$ to serve as the unit for energy (and inverse time) and
$D=0.3$. The function $g(\mathbf{r})$ describes the decay of spin-spin
interaction with distance. For simplicity, we assume a Lorentzian function for
$h(\mathbf{q})$, which leads to an exponential decaying
$g(\mathbf{r})=\mathcal{A}\,e^{-\gamma(|x|+|y|)}$, where $\mathcal{A}$ is a
normalization constant [31]. A hard cutoff $r_{c}$ such that
$\mathsf{J}(\mathbf{r})=0$ for $|x|+|y|>r_{c}$ is further introduced for
large-scale simulations of $N=1000^{2}$ spins. For results presented in the
following, parameters $\gamma=0.3$ and $r_{\rm c}=16$ are used.

The dynamical evolution of the magnet is described by the Landau-Lifshitz-
Gilbert (LLG) equation

$\displaystyle\frac{d\mathbf{S}_{i}}{dt}=\frac{1}{1+\alpha^{2}}\left[\frac{\partial\mathcal{H}}{\partial\mathbf{S}_{i}}\times\mathbf{S}_{i}+\alpha\mathbf{S}_{i}\times\left(\mathbf{S}_{i}\times\frac{\partial\mathcal{H}}{\partial\mathbf{S}_{i}}\right)\right],$
(11)

where $\alpha$ is the Gilbert damping coefficient, which is set to 0.1 in the
following simulations. A fourth-order Runge-Kutta method is used to integrate
the LLG equation with a time-step $\Delta t=0.05$. The phase diagram, shown in
Fig. 1(a) in the plane of exchange anisotropy $J^{zz}$ versus field $H$,
includes a 1$Q$-cycloidal order, a double-$Q$ SkL, and a forced ferromagnetic
state at high field. It is worth noting that our results are consistent with
the phase diagram obtained from simulated annealing minimization of the
infinite-range effective spin model [30]. The double-$Q$ magnetic order at
$H=0$, shown in Fig. 1(b), exhibits a non-coplanar Néel-type vortex texture.

Figure 1:  (a) The ground state phase diagram includes a square Skyrmion
lattice (SkL), an 1$Q$ state, and a forced ferromagnetic state (FFM). (b) The
real-space spin texture for $J^{zz}=0$ and $H=0$. The color of the arrows
represents the out-of-plane component of spins.

Importantly, for small $J^{zz}$, the SkL is stable for a wide range of
magnetic field, allowing us to study the crystallization dynamics of skyrmions
and the field effects. To this end, the LLG dynamics was employed to simulate
thermal quenches of a spin system with the effective RKKY interaction in Eq.
(3). An initial state of random spins, corresponding to equilibrium at high
temperatures, is suddenly quenched to zero temperature at time $t=0$. A
typical example of the subsequent evolution of spins is shown in Fig. 2 for
the case of zero-field quench. The color shows the scalar chirality which is
defined as $\chi_{i}^{\rm
sc}=\sum_{\triangle_{i}}\mathbf{S}_{1}\cdot\mathbf{S}_{2}\times\mathbf{S}_{3}/2$,
where the summation is taken over four triplets of nearest spins on a
triangle. As skyrmions, including the Néel vortex texture at zero field, are
characterized by non-coplanar spins that wrap around a sphere, the emergence
of a SkL domain is indicated by the staggered checkerboard-like arrangement of
positive and negative scalar chirality. As shown in Figs. 2(a) and (b), small
patches of skyrmion arrays quickly emerge after the thermal quench. However,
long-range coherence is yet to be established between different patches, and
large areas of vanishing scalar chirality mark the boundaries between skyrmion
crystallites.

Figure 2:  Spatial distribution of the scalar spin chirality $\chi_{i}^{\rm
sc}$ obtained from spin configurations at different times after a thermal
quench. Parameters $J^{zz}=0$, and $H=0$ are used in the LLG simulation of a
$1000^{2}$ system. Shown here are a selected $200\times 200$ region of the
lattice.

As long-range crystallization order further develops, the incoherent regions
quickly contract to particle-like objects of similar size as a skyrmion as
shown in Figs. 2(c) and (d). These “particles”, also characterized by a
vanishing $\chi_{i}^{\rm sc}$, correspond to dislocation defects, which are
topological defects associated with broken translational symmetries. The
dislocations here correspond to a starting point of an extra row or column of
skyrmion lines in the square lattice. These particle-like defects carry a
topological charge corresponding to the so-called Burgers vector. Their
topological nature also manifests itself in the fact that dislocations are
created and annihilated in pairs. The phase-ordering of SkL is now dominated
by the dynamics of dislocation defects.

Figure 3:  Spatial profile of the three components of spin field
$\mathbf{S}_{i}$ at late stage of the phase ordering of SkL. Parameters are
the same as those in Fig. 2.

The magnetic structure of the double-$Q$ SkL can be approximated by an equal
superposition of two spirals:

$\displaystyle\mathbf{S}(\mathbf{r}_{i})\sim\left(\cos\mathcal{Q}_{1},\cos\mathcal{Q}_{2},b(\sin\mathcal{Q}_{1}+\sin\mathcal{Q}_{2})+m_{0}\right),\quad$
(12)

where $\mathcal{Q}_{\eta}=\mathbf{Q}_{\eta}\cdot\mathbf{r}_{i}+$ const, and
the coefficients $b$ and $m_{0}$ depend on model parameters and magnetic
field. Importantly, the $x$\- and $y$-component of the spin field correspond
to simple unidirectional stripes along the $x$ and $y$ directions,
respectively. The fundamental defects of stripe order are also dislocations,
as shown in Figs. 3(a) and (b). We can thus further classify the dislocation
defects of SkL according to whether it is associated with the $S_{x}$ or
$S_{y}$ component. Indeed, our simulations find that pair annihilations are
possible only for dislocations of the same spin component.

To further quantify the phase ordering of SkL, we examine the time-dependent
spin structure factor, defined as
$\mathcal{S}(\mathbf{q},t)=\frac{1}{N^{2}}\braket{\left
missing}{\sum_{i}\mathbf{S}_{i}(t)\exp(i\mathbf{q}\cdot\mathbf{r}_{i})}{{}^{2}}$,
where $\braket{\cdots}$ denotes the average of independent initial conditions.
Examples of the structure factor at the early and late stages of phase
ordering are shown in Figs. 4(a) and (b), respectively. The structure factor
exhibits four broad peaks at $\pm{\bf Q}_{1}$ and $\pm{\bf Q}_{2}$, where
${\bf Q}_{1}=(Q,0)$, ${\bf Q}_{2}=(0,Q)$, and $Q=0.785$, quickly after the
quench; see e.g. Fig. 4(a) for $t=40$. These peaks at the nesting wave vectors
become sharper as the system relaxes toward equilibrium, as shown in Fig.
4(b). Moreover, satellite peaks at $n_{1}{\bf Q}_{1}+n_{2}{\bf Q}_{2}$ with
$n_{1}+n_{2}=2n+1$ ($n$ is an integer) start to emerge at late times,
signaling the onset of higher harmonics of the constituent spiral orders.

An overall order parameter of the SkL phase can be defined as the sum of peak
intensities $\mathcal{M}(t)=\mathcal{S}({\bf Q}_{1},t)+\mathcal{S}({\bf
Q}_{2},t)$. As shown in Fig. 4(c), this SkL order parameter clearly exhibits a
two-stage ordering discussed above: the fast development of quasi-long-ranged
crystalline domains ($t\lesssim 80$), followed by a slow power-law growth
dominated by the annihilation dynamics of dislocations ($t\gtrsim 80$). To
characterize the late-stage phase ordering, Fig. 4(c) also shows the growth of
the correlation length defined as the inverse widths of the Fourier peaks:
$\xi_{{\bf Q}_{\eta}}=1/\Delta Q_{\eta}$, where $\Delta Q_{\eta}=\sum_{{\bf
q}\sim{\bf Q}_{\eta}}|{\bf q}-{\bf
Q}_{\eta}|\mathcal{S}(\mathbf{q},t)/\sum_{{\bf q}\sim{\bf
Q}_{\eta}}\mathcal{S}(\mathbf{q},t)$. Furthermore, the power-law growth of the
correlation lengths is intimately related to that of SkL order parameter
$\mathcal{M}$. Indeed, we find that the late-stage ordering of SkL exhibits a
dynamical symmetry. The structure factor at different times can be described
by a universal function $\mathcal{F}$ with proper rescaling

$\displaystyle\mathcal{S}(\mathbf{q},t)=\xi_{\mathbf{Q}_{\eta}}^{2}\\!(t)\,\mathcal{F}\\!\left(|\mathbf{q}-\mathbf{Q}_{\eta}|\,\xi^{\,}_{\mathbf{Q}_{\eta}}(t)\right).$
(13)

This is illustrated by the excellent data points collapsing shown in Fig.
4(d).

Figure 4:  Spin structure factor $\mathcal{S}(\mathbf{q},t)$ with $J^{\rm
zz}=0$ and $H=0$ at (a) $t=40$ and (b) $t=400$. (c) The SkL order parameter
$\mathcal{M}(t)=\mathcal{S}({\bf Q}_{1},t)+\mathcal{S}({\bf Q}_{2},t)$ and the
correlation length $\xi_{{\bf Q}_{\eta}}$ versus time. Both exhibit a power-
law growth $t^{\alpha}$ with exponent $\alpha_{\mathcal{M}}\approx 1.046$ and
$\alpha_{\xi}\approx 0.476$, respectively. (d) Scaling plot of the structure
factor around the peak at ${\bf Q}_{1}$ along the $q_{x}$ axis.

The above late-time power-law behaviors are related to the dynamics of
dislocations, which are topological defects of emergent square skyrmion arrays
discussed above. As the long-range coherence of a crystalline order is
disrupted by dislocations, the correlation length of SkL can be interpreted as
the average distance $\ell$ between dislocations, which is related to their
number density as $\ell\sim\rho_{\rm d}^{-1/2}$. The power-law growth of
correlation length $\xi\sim t^{\alpha}$ thus implies a power-law decrease of
dislocation density $\rho_{\rm d}\sim t^{-\eta}$ with the exponent
$\eta=2\alpha$. This is indeed confirmed in our large-scale LLG simulations
summarized in Fig. 5, where a naive power-law fitting gives an exponent $\eta$
that depends weakly on magnetic field $H$. For example, the zero-field
exponent $\eta\sim 1.03$ is consistent with the grow exponent of the
correlation length shown in Fig. 4(c).

To understand the power-law annihilation of dislocations, we note that these
topological defects can be mapped to vortices of effective 2D-XY models. As
discussed in Eq. (12), the SkL can be viewed as comprised of two spiral spin
orders. The $x$ and $y$ components of the coarse-grained spin field represent
two orthogonal stripes
$S_{x}(\mathbf{r})\sim\cos[\mathbf{Q}_{1}\\!\cdot\\!\mathbf{r}+\theta_{1}(\mathbf{r})]$
and
$S_{y}(\mathbf{r})\sim\cos[\mathbf{Q}_{2}\\!\cdot\\!\mathbf{r}+\theta_{2}(\mathbf{r})]$.
The dislocations associated with the two stripes, shown in Figs. 3(a) and (b),
correspond to vortex singularities of the phase fields
$\theta_{1,2}(\mathbf{r})$. At the leading-order approximation, our system can
thus be described by two coupled XY models:

$\displaystyle\mathcal{E}=A\int\\{[(\nabla\theta_{1})^{2}+(\nabla\theta_{2})^{2}]+2\kappa(\nabla\theta_{1}\cdot\nabla\theta_{2})\\}d^{2}\mathbf{r},$
(14)

where $A>0$ represents the stiffness of the phase fields, and $\kappa$ denotes
their coupling. A large $\kappa$ could induce a bound state of vortices from
the two XY fields. Nonetheless, it is straightforward to show that the above
action is equivalent to two independent XY fields
$\theta_{\pm}=\theta_{1}\pm\theta_{2}$ with different stiffness:
$A_{\pm}=A(1\pm\kappa)/2$. Importantly, it has been shown that the phase
ordering of the 2D XY model is governed by the annihilation of vortices, whose
number density follows a power-law $\rho_{\rm v}\sim t^{-1}$ time dependence,
up to a logarithmic correction [36, 37]. Our results thus strongly indicate
that the ordering kinetics of the SkL phase belongs to the same dynamical
universality class of 2D XY model. In fact, while the extracted exponent is
weakly dependent on $H$, the fact that all $\eta$ values lie in the vicinity
of unity suggests a universal $\rho_{\rm d}\sim t^{-1}$ behavior, yet with a
field-dependent logarithmic correction. We have checked that the various
curves in Fig. 5 indeed can be well described by the formula $\rho_{\rm
d}/(\log(\rho_{0}/\rho_{\rm d})-1)=a(t-t_{0})^{-1}$, where $\rho_{0}$,
$t_{0}$, and $a$ are fitting parameters [36].

Figure 5:  Real-time dependence of the density of dislocations defects
$\rho_{\rm d}(t)$ with $J^{zz}=0$ for varying magnetic field $H$, obtained by
averaging over 32 independent runs. The positions of defects are detected from
the energy density [38]. The pale lines represent the fitting by the formula
$\rho_{\rm d}(t)/(\log(\rho_{0}/\rho_{\rm d}(t))-1)=a(t-t_{0})^{-1}$.

To conclude, we have presented a comprehensive study on the ordering dynamics
of SkL phases in chiral metallic magnets, which is important to the design and
engineering of skyrmion-based spintronics devices. Fundamentally, magnetic
skyrmion lattices also provide another platform to study crystallization
phenomena in two dimensions, a field which has yet to be systematically
investigated. An intriguing fast crystallization process has recently been
observed in the “square ice” formed by water molecules locked between two
graphene sheets [39]. While our work here sheds new light on this fundamental
subject of 2D crystallization dynamics, several important issues remain to be
addressed. For example, although the long-range nature of electron-mediated
interactions might be crucial to the initial fast crystallization, a detailed
study on the effect of interaction range and its interplay with other factors
is desired. The built-in chirality due to SOC is another crucial component for
the fast establishment of coherent lattices. The recent experimental
observation of square SkL in centrosymmetric magnets such as GdRu2Si2 [40] and
EuAl4 [41] thus calls for theoretical studies of skyrmion crystallization
dynamics in non-chiral itinerant magnets [23, 42]. Yet another interesting
question is the effects of crystal geometry. More complicated annihilation
dynamics of dislocations might arise in the triangular SkL observed in many
skyrmion materials.

###### Acknowledgements.

The authors thank Y. Kato and Y. Motome for fruitful discussions. This work
was supported by JSPS KAKENHI Grant Number No. JP21J20812. K.S. was supported
by the Program for Leading Graduate Schools (MERIT-WINGS). Parts of the
numerical calculations were performed in the supercomputing systems in ISSP,
the University of Tokyo. G.W.C was partially supported by the US Department of
Energy Basic Energy Sciences under Contract No. DE-SC0020330.

## References

  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nat. Nanotechnol. 8, 899 (2013).
  * Everschor-Sitte  _et al._ [2018] K. Everschor-Sitte, J. Masell, R. M. Reeve, and M. Kläui, Perspective: Magnetic skyrmions—Overview of recent progress in an active research field, Journal of Applied Physics 124, 10.1063/1.5048972 (2018), 240901.
  * Back  _et al._ [2020] C. Back, V. Cros, H. Ebert, K. Everschor-Sitte, A. Fert, M. Garst, T. Ma, S. Mankovsky, T. L. Monchesky, M. Mostovoy, N. Nagaosa, S. S. P. Parkin, C. Pfleiderer, N. Reyren, A. Rosch, Y. Taguchi, Y. Tokura, K. von Bergmann, and J. Zang, The 2020 skyrmionics roadmap, J. Phys. D: Appl. Phys. 53, 363001 (2020).
  * Göbel  _et al._ [2021] B. Göbel, I. Mertig, and O. A. Tretiakov, Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles, Physics Reports 895, 1 (2021), beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles.
  * Volovik [1987] G. E. Volovik, Linear momentum in ferromagnets, J. Phys. C Solid State Phys. 20, L83 (1987).
  * Berry [1984] M. V. Berry, Quantal phase factors accompanying adiabatic changes, Proc. R. Soc. Lond. A 392, 45 (1984).
  * Nagaosa and Tokura [2012] N. Nagaosa and Y. Tokura, Emergent electromagnetism in solids, Phys. Scr. T146, 14020 (2012).
  * Loss and Goldbart [1992] D. Loss and P. M. Goldbart, Persistent currents from Berry’s phase in mesoscopic systems, Phys. Rev. B 45, 13544 (1992).
  * Ye  _et al._ [1999] J. Ye, Y. B. Kim, A. J. Millis, B. I. Shraiman, P. Majumdar, and Z. Tešanović, Berry Phase Theory of the Anomalous Hall Effect: Application to Colossal Magnetoresistance Manganites, Phys. Rev. Lett. 83, 3737 (1999).
  * Bruno  _et al._ [2004] P. Bruno, V. K. Dugaev, and M. Taillefumier, Topological Hall Effect and Berry Phase in Magnetic Nanostructures, Phys. Rev. Lett. 93, 096806 (2004).
  * Onoda  _et al._ [2004] M. Onoda, G. Tatara, and N. Nagaosa, Anomalous Hall Effect and Skyrmion Number in Real and Momentum Spaces, J. Phys. Soc. Jpn. 73, 2624 (2004).
  * Binz and Vishwanath [2008] B. Binz and A. Vishwanath, Chirality induced anomalous-Hall effect in helical spin crystals, Physica B: Condensed Matter 403, 1336 (2008).
  * Nakazawa and Kohno [2019] K. Nakazawa and H. Kohno, Weak coupling theory of topological hall effect, Phys. Rev. B 99, 174425 (2019).
  * Shiomi  _et al._ [2013] Y. Shiomi, N. Kanazawa, K. Shibata, Y. Onose, and Y. Tokura, Topological Nernst effect in a three-dimensional skyrmion-lattice phase, Phys. Rev. B 88, 064409 (2013).
  * Mizuta and Ishii [2016] Y. P. Mizuta and F. Ishii, Large anomalous Nernst effect in a skyrmion crystal, Sci. Rep. 6, 28076 (2016).
  * Hirschberger  _et al._ [2020] M. Hirschberger, L. Spitz, T. Nomoto, T. Kurumaji, S. Gao, J. Masell, T. Nakajima, A. Kikkawa, Y. Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, R. Arita, T.-h. Arima, and Y. Tokura, Topological Nernst Effect of the Two-Dimensional Skyrmion Lattice, Phys. Rev. Lett. 125, 076602 (2020).
  * Bogdanov and Yablonskii [1989] A. N. Bogdanov and D. Yablonskii, Thermodynamically stable “vortices” in magnetically ordered crystals. The mixed state of magnets, Zh. Eksp. Teor. Fiz 95, 178 (1989).
  * Roessler  _et al._ [2006] U. K. Roessler, A. Bogdanov, and C. Pfleiderer, Spontaneous skyrmion ground states in magnetic metals, Nature 442, 797 (2006).
  * Mühlbauer  _et al._ [2009] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion Lattice in a Chiral Magnet, Science 323, 915 (2009).
  * Tokura and Kanazawa [2020] Y. Tokura and N. Kanazawa, Magnetic Skyrmion Materials, Chem. Rev. (2020).
  * Reichhardt  _et al._ [2022] C. Reichhardt, C. J. O. Reichhardt, and M. V. Milošević, Statics and dynamics of skyrmions interacting with disorder and nanostructures, Rev. Mod. Phys. 94, 035005 (2022).
  * Ozawa  _et al._ [2016] R. Ozawa, S. Hayami, K. Barros, G. W. Chern, Y. Motome, and C. D. Batista, Vortex crystals with chiral stripes in itinerant magnets, J. Phys. Soc. Jpn. 85, 103703 (2016).
  * Hayami  _et al._ [2017] S. Hayami, R. Ozawa, and Y. Motome, Effective bilinear-biquadratic model for noncoplanar ordering in itinerant magnets, Phys. Rev. B 95, 224424 (2017).
  * Hayami and Motome [2021a] S. Hayami and Y. Motome, Topological spin crystals by itinerant frustration, Journal of Physics: Condensed Matter 33, 443001 (2021a).
  * Bray [1994] A. J. Bray, Theory of phase-ordering kinetics, Advances in Physics 43, 357 (1994).
  * Puri and Wadhawan [2009] S. Puri and V. Wadhawan,  _Kinetics of phase transitions_ (CRC press, 2009).
  * Ruderman and Kittel [1954] M. A. Ruderman and C. Kittel, Indirect Exchange Coupling of Nuclear Magnetic Moments by Conduction Electrons, Phys. Rev. 96, 99 (1954).
  * Kasuya [1956] T. Kasuya, A Theory of Metallic Ferro- and Antiferromagnetism on Zener’s Model, Prog. Theor. Phys. 16, 45 (1956).
  * Yosida [1957] K. Yosida, Magnetic properties of cu-mn alloys, Phys. Rev. 106, 893 (1957).
  * Hayami and Motome [2018] S. Hayami and Y. Motome, Néel- and Bloch-Type Magnetic Vortices in Rashba Metals, Phys. Rev. Lett. 121, 137202 (2018).
  * [31] A detailed discussion of the effective Hamiltonian for the square-lattice s-d model with Rashba SOC is presented in the supplemental material.
  * Okumura  _et al._ [2020] S. Okumura, S. Hayami, Y. Kato, and Y. Motome, Magnetic hedgehog lattices in noncentrosymmetric metals, Phys. Rev. B 101, 144416 (2020).
  * Yambe and Hayami [2021] R. Yambe and S. Hayami, Skyrmion crystals in centrosymmetric itinerant magnets without horizontal mirror plane, Sci. Rep. 11, 11184 (2021).
  * Shimizu  _et al._ [2021] K. Shimizu, S. Okumura, Y. Kato, and Y. Motome, Phase transitions between helices, vortices, and hedgehogs driven by spatial anisotropy in chiral magnets, Phys. Rev. B 103, 054427 (2021).
  * Kato  _et al._ [2021] Y. Kato, S. Hayami, and Y. Motome, Spin excitation spectra in helimagnetic states: Proper-screw, cycloid, vortex-crystal, and hedgehog lattices, Phys. Rev. B 104, 224405 (2021).
  * Yurke  _et al._ [1993] B. Yurke, A. N. Pargellis, T. Kovacs, and D. A. Huse, Coarsening dynamics of the xy model, Phys. Rev. E 47, 1525 (1993).
  * Bray  _et al._ [2000] A. J. Bray, A. J. Briant, and D. K. Jervis, Breakdown of scaling in the nonequilibrium critical dynamics of the two-dimensional $\mathit{XY}$ model, Phys. Rev. Lett. 84, 1503 (2000).
  * [38] A detailed discussion of the way to count the dislocations and the relationship between dislocations in the present system and vortices in the XY model are presented in the supplemental material.
  * Algara-Siller  _et al._ [2015] G. Algara-Siller, O. Lehtinen, F. C. Wang, R. R. Nair, U. Kaiser, H. A. Wu, A. K. Geim, and I. V. Grigorieva, Square ice in graphene nanocapillaries, Nature 519, 443 (2015).
  * Khanh  _et al._ [2020] N. D. Khanh, T. Nakajima, X. Yu, S. Gao, K. Shibata, M. Hirschberger, Y. Yamasaki, H. Sagayama, H. Nakao, L. Peng, K. Nakajima, R. Takagi, T. Arima, Y. Tokura, and S. Seki, Nanometric square skyrmion lattice in a centrosymmetric tetragonal magnet, Nat. Nanotechnol. 15, 444 (2020).
  * Takagi  _et al._ [2022] R. Takagi, N. Matsuyama, V. Ukleev, L. Yu, J. S. White, S. Francoual, J. R. Mardegan, S. Hayami, H. Saito, K. Kaneko, _et al._ , Square and rhombic lattices of magnetic skyrmions in a centrosymmetric binary compound, Nat. Commun. 13, 1472 (2022).
  * Hayami and Motome [2021b] S. Hayami and Y. Motome, Square skyrmion crystal in centrosymmetric itinerant magnets, Phys. Rev. B 103, 024439 (2021b).


# Wakefield-driven filamentation of warm beams in plasma

Erwin Walter<EMAIL_ADDRESS>Max Planck Institute for Plasma Physics,
85748 Garching, Germany  Exzellenzcluster ORIGINS, 85748 Garching, Germany
John P. Farmer<EMAIL_ADDRESS>Max Planck Institute for Physics, 80805
Munich, Germany  Martin S. Weidl  Max Planck Institute for Plasma Physics,
85748 Garching, Germany  Alexander Pukhov  University of Duesseldorf, 40225
Duesseldorf, Germany  Frank Jenko  Max Planck Institute for Plasma Physics,
85748 Garching, Germany

###### Abstract

Charged and quasi-neutral beams propagating through an unmagnetised plasma are
subject to numerous collisionless instabilities on the small scale of the
plasma skin depth. The electrostatic two-stream instability, driven by
longitudinal and transverse wakefields, dominates for dilute beams. This leads
to modulation of the beam along the propagation direction and, for wide beams,
transverse filamentation. A three-dimensional spatiotemporal two-stream theory
for warm beams with a finite extent is developed. Unlike the cold beam limit,
diffusion due to a finite emittance gives rise to a dominant wavenumber, and a
cut-off wavenumber above which filamentation is suppressed. Particle-in-cell
simulations give excellent agreement with the theoretical model. This work
provides deeper insights into the effect of diffusion on filamentation of
finite beams, crucial for comprehending plasma-based accelerators in
laboratory and cosmic settings.

††preprint: APS/123-QED

##  I Introduction

From supernovae in distant galaxies to laboratory-based wakefield
accelerators, the collisionless interaction of relativistic particles with
plasma is relevant to many physical scales. The interactions are often
governed by kinetic micro-instabilities, which result in electrostatic and
electromagnetic fluctuations [1, 2, 3]. This dissipation of a directed
relativistic flow transfers kinetic energy to field energy, which can give
rise to collisionless shocks in the astrophysical regime. In these
collisionless shocks, non-thermal particles accelerated to TeV energies
through Fermi-type processes [4] or Landau resonance [5, 6, 7] emit
synchrotron radiation across a spectrum from radio to gamma-ray frequencies
[8, 9]. Collisionless shocks are observed in active galactic nuclei and
supernovae-remnants [10], or in gamma-ray bursts that occur during merge
events of neutron stars or black holes [11, 12].

Specially designed experimental setups [13, 14, 15] have recently enabled
unprecedented investigations of electromagnetic plasma instabilities relevant
on the astronomical scale. Beam-driven plasma wakefield accelerators (PWFAs)
[16], which can be utilised as $\gamma$-ray sources [17] or to achieve higher
accelerating fields compared to conventional RF accelerators [18, 19, 20], are
also subject to microinstabilities. Furthermore, PWFAs can be adapted to
investigate regimes relevant to astrophysics [21, 22].

The interaction of a relativistic beam with an unmagnetised plasma can be
usually categorised between the electromagnetic Weibel-like current
filamentation instability (CFI) [23, 24], driven by the plasma return current,
or two-stream instabilities [25, 26, 27], driven by the electrostatic plasma
response. In the latter, the beam excites Langmuir plasma waves [28],
conventionally called wakefields in particle accelerators [29], which lead to
the (longitudinal) two-stream instability (TSI) and the transverse two-stream
instability (TTS) [30]. The combination of TSI and TTS is usually referred to
as the oblique instability (OBI) [31, 27] and allows dilute beams to undergo a
similar filamentary behaviour as CFI.

Previous theoretical work on CFI for cold, spatially uniform streams
determined that the temporal growth rate increases with transverse wavenumber
[32]. These studies were extended to warm streams, in which diffusion acts to
suppress small-scale filamentation, and a dominant wavenumber was calculated
[33, 34]. For cold longitudinally bounded streams, CFI was found to exhibit
spatiotemporal growth at the beam head [35].

For two-stream instabilities in cold uniform streams, the growth rate also
increases with transverse wavenumber [31, 27]. It was predicted that diffusion
would suppress the growth of small-scale filaments [25], which was later
studied numerically, and a threshold above which the system is stable was
found analytically [36]. For a localised disturbance in cold bounded systems,
TSI [37, 38] and TTS [39] demonstrate a pulse-shaped spatiotemporal growth.
However, the effect of a finite beam emittance on the spatiotemporal growth of
the filamentation instability has not previously been treated analytically.

This manuscript introduces a fully three-dimensional, spatiotemporal theory
describing filamentation of a warm beam due to wakefield-driven two-stream
instabilities. This allows limits to be set on the beam temperature for
laboratory astrophysics schemes seeking to investigate these instabilities and
PWFA experiments seeking to avoid them. The work is structured as follows:
Wakefield-driven filamentation is introduced in Section II. In Section III, an
analytical expression for the growth is derived for a cold bounded beam with a
transverse profile. The theory is extended to warm beams in Section IV, which
considers the effect of diffusion. This allows the exact value for the
dominant wavenumber to be calculated, as well as the cut-off above which no
filamentation occurs. The analytical predictions are throughout compared to
two and three-dimensional particle-in-cell (PIC) simulations, which show
excellent agreement.

##  II Wakefield-Driven Filamentation

The regimes for the two filamentation instabilities are defined by the current
imbalance in the system. The beam and plasma currents must be comparable for
CFI to dominate. For a relativistic beam propagating in stationary plasma,
relevant to many astrophysical schemes, this requires a dense beam,
$n_{b}\gtrsim n_{p}$, with $n_{b}$ and $n_{p}$ the beam and plasma density
[2]. For a dilute beam, $n_{b}\ll n_{p}$, the plasma current is negligible,
and plasma electrons are mainly deflected by the beam charge. The resulting
wakefield leads to TSI and TTS [2, 40, 41].

Plasma wakefield experiments use a charged beam, which is usually dense and
short, $k_{p}\sigma_{\zeta}<1$, with $\sigma_{\zeta}$ the rms length [42, 43].
Here, $k_{p}=\omega_{p}/c$ is the plasma wavenumber, where $c$ is the speed of
light and $\omega_{p}=[e^{2}n_{p}/(\varepsilon_{0}m_{e})]^{1/2}$ is the plasma
frequency, with $e$ the elementary charge, $m_{e}$ the electron mass, and
$\varepsilon_{0}$ the vacuum permittivity. A dilute and long beam,
$k_{p}\sigma_{\zeta}\gg 1$, is subject to TTS. For narrow beams,
$k_{p}\sigma_{r}\lesssim 1$, with $\sigma_{r}$ the rms width, TTS can take the
form of the axisymmetric self-modulation instability (SMI) modulating the beam
radius [18], or the antisymmetric hosing instability displacing the beam
centroid [44].

Fully modulated, the beam can resonantly drive a quasi-linear wake with an
accelerating field comparable to that driven by a short, dense beam [45].
Wakefield experiments do not utilise wide beams as they may undergo
filamentation due to transverse perturbations [21] and degrade the wakefield.
Experiments that investigate filamentation instabilities may operate with
quasi-neutral beams (equal populations of particles with opposite charge) to
suppress SMI [46].

Figure 1:  Filamentation of a dilute bunch due to wakefields. Transverse and
longitudinal slices of the a) positron and b) electron density of a dilute
quasineutral bunch after propagating $\omega_{\beta}\tau=2.6$ in an initially
uniform plasma. c) The electron density of the corresponding plasma response.
The bunch propagates to the right, with its head at $\zeta=0$. The transverse
and longitudinal slices are taken at $k_{p}\zeta=-18\pi$ and $k_{p}x=0$,
respectively.

This filamentation of a quasi-neutral, dilute bunch and the corresponding
plasma response is shown in Fig. 1 after propagating $2.6/\omega_{\beta}$ in
plasma, where
$\omega_{\beta}=[q^{2}n_{b}/(2\gamma_{b}\varepsilon_{0}m_{b})]^{1/2}$ is the
betatron frequency, with $q$ the charge, $\gamma_{b}$ the Lorentz factor and
$m_{b}$ the mass of the bunch particles. Both the bunch and the plasma
response exhibit roughly equidistant filaments, where positrons and electrons
are oppositely aligned due to the plasma wakefield that drives the
instability. The plasma electrons align with the bunch positrons driven by the
bunch charge. A periodic modulation occurs along the bunch, arising due to the
oscillation of the wakefield.

The simulation in Fig. 1 was carried out using the three-dimensional,
quasistatic PIC code qv3d [47], built on the VLPL platform [48]. The
relativistic, $\gamma_{b}=22.4$ ($u_{b}/c=0.999$), warm electron-positron
bunch has a longitudinally flat-top profile with extent $-20\pi<\zeta<0$, and
along each transverse axis a Gaussian profile with rms width of
$k_{p}\sigma_{r}=3$ and a momentum width of $\sigma_{pr}/(m_{b}c)=0.05$. The
momentum width is related to the normalised beam emittance $\epsilon_{N}$ by
$\sigma_{pr}/(m_{b}c)=\epsilon_{N}/\sigma_{r}$. The peak density of the bunch
positrons and electrons is $n_{b}/2=0.02\,n_{p}$, i.e. $n_{b}$ is the total
peak density of the bunch. The bunch propagation through a uniform plasma is
considered in the co-moving frame $\zeta=z-u_{b}t,\,\tau=z/u_{b}$, with
$\zeta$ the bunch slice, $u_{b}$ the bulk velocity of the bunch and $\tau$ the
propagation time in plasma. The grid size is
$k_{p}\Delta(x,y,\zeta)=(0.01,0.01,0.1)$, the propagation step is $k_{p}\Delta
z=2$, and the bunch and plasma species with stationary modelled plasma ions
are represented by 16 and 4 macroparticles per cell.

From theory, the filamentation growth rate increases with transverse
wavenumber for a cold bunch. In simulations, the finite spatial resolution
limits the maximum wavenumber which can be modeled. This leads to a dominant
wavenumber determined by the cell size. For the finite emittance considered in
Fig. 1, diffusion results in a physical reduction of the growth rate at higher
wavenumbers, yielding a dominant wavenumber well within the resolution limit
of the simulation. In the next section, an analytical model is developed for
wakefield-driven two-stream instabilities.

##  III Filamentation of Cold Beams

The charge density of the bunch drives an electrostatic plasma response,
expressed as the longitudinal $E_{z}$ and transverse $\bm{W}_{\perp}$
wakefield [41]. For a bunch charge density $\rho_{b}=q\delta n_{b}g(x,y)$,
where $\delta n_{b}$ is the amplitude of the density modulation and
$g(x,y)=\tilde{g}(x,y)\cos{(k_{x}x+\varphi_{x})}\cos{(k_{y}y+\varphi_{y})}$ is
the transverse profile, with $\tilde{g}$ a slowly-varying envelope, and
$k_{x,y}$ and $\varphi_{x,y}$ the perturbation wavenumbers and phases along
the transverse axes, the wakefield is given in the linear regime by (Section
A.1)

$\displaystyle E_{z}$ $\displaystyle=\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{k_{e}^{2}g(x,y)}{k_{e}^{2}+k_{r}^{2}}\int_{\zeta}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})\cos{k_{e}(\zeta-\zeta^{\prime})}$
(1) $\displaystyle\bm{W}_{\perp}$ $\displaystyle=\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{k_{e}\nabla_{\perp}g(x,y)}{k_{e}^{2}+k_{r}^{2}}\int_{\zeta}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})\sin{k_{e}(\zeta-\zeta^{\prime})},$

with $k_{e}=k_{p}c/u_{b}$ the electron wavenumber,
$k_{r}=(k_{x}^{2}+k_{y}^{2})^{1/2}$, and $f(\zeta)$ the longitudinal profile.
The linear regime requires $\delta n_{b}\ll n_{b}$ and
$|E_{z}|,|\bm{W}_{\perp}|\ll E_{0}$, with $E_{0}=m_{e}\omega_{p}c/e$ the non-
relativistic wave-breaking field.

The local self-fields of the bunch are neglected and become relevant for low
bunch velocities. The wakefield acts back on the bunch, where particles are
accelerated or decelerated by $E_{z}$ and focussed or defocussed by
$\bm{W}_{\perp}$. The evolution of a cold bunch is described by the fluid
equation [1]

$\partial_{\tau}^{2}\delta
n_{b}=\frac{2\omega_{\beta}^{2}}{q/\varepsilon_{0}}\left(\frac{\partial_{z}E_{z}}{\gamma_{b}^{2}}+\nabla_{\perp}\cdot\bm{W}_{\perp}\right).$
(2)

A Laplace transform solves the Green’s function to Eq. 2 (Section A.2). For a
longitudinal flat-top bunch with the head at $\zeta=0$, the growth of the
modulation amplitude with respect to its initial value $\delta n_{b0}$ is

$\begin{gathered}\Gamma_{\mathrm{TS}}=\frac{\delta n_{b,\mathrm{TS}}}{\delta
n_{b0}}=\left|\sum_{n=0}^{\infty}\frac{\left[i\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2}\right]^{n}}{n!(2n)!}\right|\\\
\eta_{u}=\frac{(c^{2}-u_{b}^{2})k_{p}^{2}+u_{b}^{2}k_{r}^{2}}{c^{2}k_{p}^{2}+u_{b}^{2}k_{r}^{2}}.\end{gathered}$
(3)

The first and second summand in the spectral factor $\eta_{u}$ represent the
respective contribution of the longitudinal and transverse wakefield component
and, therefore, of TSI and TTS. The spectral dependency agrees with the
analytical expression for OBI in [26, 2]. The series can be asymptotically
expressed by $\delta n_{b,\mathrm{TS}}\approx[\delta
n_{b0}/\sqrt{4\pi}]\exp{\Gamma_{\infty}}/\sqrt{\Gamma_{\infty}}$, with
$\Gamma_{\infty}=(3^{3/2}/2^{5/3})[\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2}]^{1/3}$.
In the non-relativistic and ultra-relativistic limit for streams, the
asymptotic form simplifies to previous works [38, 39]. The phase velocity of
the growing wave (Section A.2)

$\displaystyle u_{\psi}$
$\displaystyle=u_{b}\left[1-\frac{2}{3^{3/2}}\frac{\Gamma_{\infty}}{\omega_{p}\tau}\right]$
(4)

reduces relative to the bunch velocity.

The analytical expressions for the growth and phase velocity of TTS are
similar to the axisymmetric SMI [18, 49], with the key distinction being that
the spectral factor $\eta_{u}$ is substituted by Bessel functions. For a
single-species bunch, the growth rate of a transverse modulation within the
bunch exceeds the rate at which the transverse envelope changes for
$\sigma_{r}\gtrsim 3/k_{r}$.

In order to test this analytic description of two-stream instabilities,
comparisons are made to simulations. Two-dimensional simulations are used, in
which relativistic beam particles effectively have one degree of freedom and
$k_{r}=k_{y}$. The two-dimensional simulations are carried out with the
electromagnetic PIC code OSIRIS [50]. The grid size of the simulation is
$k_{p}\Delta(y,z)=(0.02,0.06)$, and the time step is set to $\omega_{p}\Delta
t=0.0172$. The bunch and plasma species are each represented by 384 and 192
macroparticles per cell. The number of particles per cell is significantly
higher than that used in the three-dimensional simulations due to the relative
decrease in the total number of cells in the two-dimensional simulations. The
boundary conditions are open for the macroparticles and electromagnetic
fields. The bunch is initialised in a vacuum and propagates into plasma. The
bunch parameters are equivalent to Fig. 1, but with an initially cold beam.
The bunch charge density is transversely modulated with an amplitude of
$e\delta n_{b}=0.01\sqrt{2}\,en_{p}$ and a wavenumber of $k_{y}/k_{p}=\pi$ to
allow filaments to develop.

Figure 2:  Bunch-driven wakefield immediately after entering the plasma,
$\tau=0$. Bunch charge density and corresponding longitudinal and transverse
wakefield for a a) wide quasi-neutral bunch with transverse modulation and b)
narrow single-species bunch of identical amplitude $\delta n_{b}$ and width
$k_{y}$.

Figure 2a) shows the initial wakefield driven by the transversely modulated
bunch when each bunch slice just entered the plasma, $\tau=0$. The
longitudinal and transverse wakefield exhibit a longitudinal modulation at
$k_{\zeta}=k_{e}$ and a transverse modulation at the seeded wavenumber
$k_{y}=\pi k_{p}$. The transverse wakefield is stronger than the longitudinal
component in agreement with the theoretical ratio
$\tilde{\bm{W}}_{\perp}=\tilde{E}_{z}\bm{k}_{r}u_{b}/(k_{p}c)$ from Eq. 1. For
comparison, the wakefield driven by a narrow single-species bunch is shown in
Fig. 2b). Unlike the wide bunch, the wakefield extends beyond the narrow
bunch. However, in both cases, the transverse wakefield periodically
alternates between focussing and defocussing along the bunch, which gives rise
to TTS for a transversely modulated bunch or SMI for a narrow single-species
bunch.

Figure 3:  The evolution of a modulated bunch propagating in plasma. The a)
bunch charge density, b) electric and c) magnetic field of a filamented bunch
with a transverse modulation, $k_{y}/k_{p}=\pi$, at $\omega_{\beta}\tau=2.6$.
d) The transversely averaged electric field along the plasma length with the
theoretical phase shift superimposed. e) The envelope growth of the electric
field $\tilde{E}_{y}$, along the length of the bunch $\zeta$ at the
propagation times $\omega_{\beta}\tau=\\{0.8,1.6,2.6\\}$, showing the
simulation and theoretical values by solid and dotted lines, respectively. The
black dotted line indicates the seed of the electric field. f) Simulated and
theoretical envelope growth of the electric field along the propagation $\tau$
at equidistant bunch slices $k_{p}\zeta=\\{-6\pi,-12\pi,-18\pi\\}$. g) Bunch
profile from simulation, obtained by fitting the observed growth to Eq. 3,
compared to the initial bunch shape.

The resulting growth of the filamentation instability from the initial plasma
response in Fig. 2a) is illustrated in Fig. 3 at a propagation of
$2.6/\omega_{\beta}$ in plasma. The modulation amplitude of the bunch charge
density in Fig. 3a) increases along the bunch length, and contains a
longitudinal modulation at $k_{\zeta}=k_{e}$ due to the electrostatic plasma
response. The transverse wakefield from Fig. 3b) and c) alternates between
focusing and defocusing, both transversely and along the bunch, resulting in
alternating positron and electron filaments. The magnetic field in Fig. 3c) is
weaker than the electric field by an order of magnitude and is predominantly
due to the local bunch current. For a relativistic bunch, Coulomb repulsion is
compensated by the magnetic field, so the beam evolution is determined
entirely by the plasma wakefield.

The electric field (taken as the average over the range $0<k_{y}y<\pi$) in
Fig. 3d) shows the growth along the bunch length as the bunch propagates in
plasma. The modulation shifts backwards, illustrating that the phase velocity
is lower than the bunch velocity. The superimposed lines represent the
integral of the phase velocity from Eq. 4 over the length of the plasma and
agree well with the phase of the wave.

Figure 3e) and f) show the envelope growth of the electric field (averaged
over the range $-\pi<k_{y}y<\pi$) along the bunch and the plasma length,
respectively. The seed value agrees well with the analytic expression for the
Fourier spectrum $\hat{E}_{y0}=\mathcal{F}_{\perp}\\{E_{y0}\\}=[e\delta
n_{b0}/\varepsilon_{0}]k_{y}/(k_{p}^{2}+k_{y}^{2})$, obtained by solving Eq. 1
for the initial bunch profile. The growth of the electric field is compared
with the semi-analytic solution to Eq. 3, including the first ten terms, and
shows excellent agreement along the bunch up to a propagation time in plasma
of $\sim 2/\omega_{\beta}$.

To demonstrate the effect of a slowly varying transverse envelope on the
growth, Eq. 3 is fitted to the simulation data along the plasma length at
$k_{p}\zeta=-12\pi$ with $\tilde{g}(y)$ as a free parameter. The fit
coefficient agrees well with the Gaussian profile of the bunch in Fig. 3g). In
contrast to a longitudinal extent resulting in an increase of the growth along
the bunch, the growth rate and seed level correlate with the transverse
envelope $\tilde{g}(y)$. The growth rate at a given transverse coordinate can
be treated as a stream with the local bunch density. The curved phase fronts
in the bunch modulation are due to the dependency of the phase velocity on the
transverse envelope, $\tilde{g}(y)^{1/3}$, in Eq. 4.

Simulations show that beyond $\omega_{\beta}\tau=2$, the field growth begins
to decrease relative to the analytical predictions (Fig. 3d,f) while the phase
velocity increases (Fig. 3e). This saturation is due to the beam becoming
fully modulated, with the electron and positron filaments fully separating, as
seen in Fig. 3a) for $k_{p}\zeta\leq-10\pi$.

##  IV Filamentation of Warm Beams

The filamentation of the bunch depicted in Fig. 1 results in a dominant
wavenumber, a behaviour the theory for cold bunches cannot describe. For warm
bunches, diffusion has to be included in the bunch evolution. For non-
relativistic temperatures, $\sigma_{pr}^{2}/(m_{b}c)^{2}\ll 1$, the fluid
equation in Eq. 2 extends to (Section A.3)

$\left(\partial_{\tau}^{2}+\frac{2}{3}\frac{\sigma_{pr}^{2}k_{r}^{2}}{m_{b}^{2}\gamma_{b}^{2}}\right)\delta
n_{b}=\frac{2\omega_{\beta}^{2}}{q/\varepsilon_{0}}\left(\frac{\partial_{z}E_{z}}{\gamma_{b}^{2}}+\nabla_{\perp}\cdot\bm{W}_{\perp}\right),$
(5)

where diffusion acts to damp transverse density modulations [51]. Since all
bunch slices are equally affected by diffusion, damping is purely temporal and
can be treated separately from the spatiotemporal growth of the filamentation
instability. The exponential damping rate $\delta_{D}$ of a transverse
perturbation can be described by

$\displaystyle\delta n_{b,D}$ $\displaystyle=\delta
n_{b}\exp{(-\delta_{D}\tau)}$ (6) $\displaystyle\delta_{D}$
$\displaystyle=\sqrt{\frac{2}{3}}\frac{\sigma_{pr}k_{r}}{\gamma_{b}m_{b}}.$

The total growth rate is, therefore, the sum of the growth rate from two-
stream instabilities with the damping rate from diffusion, expressed by

$\begin{gathered}\Gamma_{\mathrm{tot}}=\delta n_{b}/\delta
n_{b0}=\Gamma_{\mathrm{TS}}\exp{(-\delta_{D}\tau)}\end{gathered}$ (7)

The effect of temperature can only be considered as purely diffusive for
$\sigma_{pr}/(m_{b}c)<[3/2^{10/3}(n_{b}/n_{p})^{1/3}\gamma_{b}^{1/3}(1+\gamma_{b}^{-2})^{2/3}/(1+\gamma_{b}^{-1})^{2}]^{1/2}$
[36]. This corresponds to $\sigma_{pr}/(m_{b}c)<0.2$ for the bunch parameters
in Fig. 1.

The influence of diffusion on the filamentation instability is examined for
bunches with different temperatures. Since diffusion has a larger effect at
higher wavenumbers, the parameters are as for the bunch in Fig. 3 but with a
transverse modulation at $k_{y}/k_{p}=2\pi$. The excited electric field is
shown in Fig. 4a) at $2.6/\omega_{\beta}$. The field is lower compared to Fig.
3b) due to the difference in wavenumber, agreeing with $\hat{E}_{y}\sim
k_{y}/(k_{p}^{2}+k_{y}^{2})$ from Eq. 1. For the cold bunch, the seeded
wavenumber continues to dominate along the length of the bunch.

Figure 4:  a) The transverse electric field resulting from a cold and warm
bunches of increasing temperature with a transverse modulation at
$k_{y}/k_{p}=2\pi$. b) The corresponding growth of the spectrum at the seeded
wavenumber $\hat{E}_{y}=|\hat{E}_{y}|(k_{y})$, obtained from simulation and
theory at different slices $\zeta$.

For warm bunches, the phase fronts deviate from the curve given by the bunch
profile. The field reduces with temperature close to the bunch head since the
filamentation instability grows along the bunch while diffusion is spatially
uniform. The transverse modulation shifts from the seeded wavenumber, a change
that becomes evident further away from the bunch head. The growth of the field
spectrum along the plasma length in Fig. 4b) reveals that the seeded
wavenumber is damped proportionally to the bunch temperature. The observation
is in good agreement with the analytical description for the effect of
diffusion on the growth in Eq. 7.

The development of filaments with wavenumbers lower than the seeded wavenumber
indicates a higher growth rate for larger-scale filaments, such that the whole
transverse spectrum of the instability has to be considered. In order to
investigate the variation of the filamentation wavenumber, the electric fields
corresponding to the transverse slice at $k_{p}\zeta=-12\pi$ in Fig. 1 are
shown in Fig. 5a). The transverse component $E_{y}$ is predominantly modulated
along $y$, and $E_{x}$ is predominantly modulated along $x$. However,
transverse modulations occur with a broad range of spatial scales and
orientations in the transverse plane.

Figure 5:  a) The transverse electric fields driven by the filamented bunch
are shown as a transverse slice at $\omega_{\beta}\tau=2.6$,
$k_{p}\zeta=-12\pi$. b) The corresponding 2D power spectrum, and c) the 1D
power spectrum $|\hat{\bm{E}}_{\perp}|(k_{r})$, showing the results from
simulation and theory, as well as the theoretical growth $\Gamma$. The dotted
and dash-dotted lines indicate the theoretical value for the wavenumber of
maximum electric field amplitude $k_{E_{\mathrm{max}}}$, and the cut-off
wavenumber $k_{\mathrm{cut}}$.

For unseeded bunches, the instability grows from fluctuations in the bunch due
to the finite temperature, and the resulting electric field is a superposition
of all growing transverse modulations. The respective contributions of the
wavenumbers can be separated by a Fourier transform. Taking the two-
dimensional Fourier transform of the transverse electric field components and
plotting the absolute amplitude, i.e.
$|\hat{\bm{E}}_{\perp}|=|\hat{E}_{y}|+|\hat{E}_{x}|$, in Fig. 5b) reveals a
wide range of growing transverse wavenumbers. The spectrum is radially
symmetric, showing that growing transverse modulations have no preferred
orientation in the transverse plane. The radial symmetry is in agreement with
the spectral factor in Eq. 3, $\eta_{1}\rightarrow
k_{r}^{2}/(k_{p}^{2}+k_{r}^{2})$, which predicts that the growth rate of the
filamentation instability only depends on the absolute value of the transverse
wavevector. Thus, the filamentation in transverse planes is coupled, and the
transverse modulations in each plane cannot be treated independently.

The spectrum of the electric field grows with transverse wavenumbers up to
$k_{r}/k_{p}\sim 5$ due to the higher growth rate of the filamentation
instability and reduces for higher wavenumbers due to diffusion. This leads to
a transverse wavenumber of maximum growth
$k_{\Gamma_{\mathrm{max}}}(\tau,\zeta)$ and cut-off wavenumber
$k_{\mathrm{cut}}(\tau,\zeta)$ above which the instability is suppressed. In
the relativistic limit, the wavenumbers are numerically obtained by solving
the following expressions for $k_{\Gamma_{\mathrm{max}}}$ or
$k_{\mathrm{cut}}$ (Section A.3)

$\displaystyle 2\Gamma_{\infty}(k_{\Gamma_{\mathrm{max}}})$
$\displaystyle=3(1+k_{\Gamma_{\mathrm{max}}}^{2})\delta_{D}(k_{\Gamma_{\mathrm{max}}})\tau+1$
(8) $\displaystyle\Gamma_{\infty}(k_{\mathrm{cut}})$
$\displaystyle=\delta_{D}(k_{\mathrm{cut}})\tau+\ln{\sqrt{4\pi\Gamma_{\infty}(k_{\mathrm{cut}})}}.$

The wavenumber of maximum growth scales by
$k_{\Gamma_{\mathrm{max}}}\sim\sigma_{pr}^{-1/3}$ and the cut-off wavenumber
scales by $k_{\mathrm{cut}}\sim\sigma_{pr}^{-1}$ with the bunch temperature.
Since the two-stream instability is spatiotemporal, while diffusion is
spatially uniform, the characteristic wavenumbers depend on the propagation
time in plasma and position within the bunch. The spectrum of the electric
field is obtained by multiplying the seed spectrum of the electric field with
the growth spectrum, $\hat{\bm{E}}_{\perp}=\hat{\bm{E}}_{\perp
0}(k_{r})\Gamma_{\mathrm{tot}}(k_{r})$. The seed is assumed to vary inversely
with the square root of the number of particles within a filamentation length,
$\hat{\bm{E}}_{\perp 0}\sim\sqrt{k_{r}}$, for which the wavenumber of maximum
spectral value $k_{E_{\mathrm{max}}}(\tau,\zeta)$ is obtained from

$1+3k_{E_{\mathrm{max}}}^{2}+4\Gamma_{\infty}(k_{E_{\mathrm{max}}})=6(1+k_{E_{\mathrm{max}}}^{2})\delta_{D}(k_{E_{\mathrm{max}}})\tau.$
(9)

Averaging the spectrum of the electric field in Fig. 5b) over all
orientations, $(k_{x},k_{y})\rightarrow k_{r}$, shows the radial spectrum in
Fig. 5c). Excellent agreement is shown between the simulation and the product
of the theoretical growth spectrum with the assumed initial spectrum of the
field, $|\hat{\bm{E}}|_{\perp}\sim\sqrt{k_{r}}\Gamma_{\mathrm{tot}}$. The
absolute value of the spectrum from theory is chosen to align with the
simulation. The predicted wavenumber at which the electric field is maximum,
$k_{E\max}\approx 4.9$, from Eq. 9 aligns well with the simulation data. The
electric field above the calculated cut-off wavenumber, $k_{r}/k_{p}\gtrsim
50$, is attributed to numerical noise.

The whole scope of the introduced theory is compared to two and three-
dimensional simulations of unseeded warm bunches with different temperatures.
Other parameters are as for the bunch in Fig. 1. The growth spectrum from
simulations is obtained by
$\hat{\bm{E}}_{\perp}(k_{r})/(k_{r}\sigma_{pr}^{3})^{1/2}$ and aligned to the
growth spectrum from theory for two-dimensional and three-dimensional
simulations, respectively. Small variations in the field spectrum can occur
when the filamentation instability grows from random fluctuations in the
bunch. Thus, the growth spectrum is averaged over three three-dimensional runs
and five two-dimensional runs for each temperature and compared to the
analytical expression for the total growth in Eq. 7.

Agreement is found for the dependency of the growth spectrum on the
temperature for both two-dimensional and three-dimensional simulations, shown
in Fig. 6. The alignment is better in three dimensions since the total number
of bunch particles is an order of magnitude higher. For cold bunches, theory
predicts that the growth increases with wavenumber due to the filamentation
instability. For warm bunches, the growth increases with wavenumber up to
$k_{\Gamma_{\mathrm{max}}}$ and then decreases as the influence of diffusion
becomes stronger. With higher temperatures, the growth is lower for all
wavenumbers and the wavenumber of maximum growth and cut-off wavenumber shift
to lower values in good agreement with the predicted values from evaluating
Eq. 8. Thus, transverse modulations in the bunch occur at larger scales. The
distance between filaments is inversely related to $k_{E_{\mathrm{max}}}$.
However, this means that the in-plane distance is higher in three-dimensional
simulations with $k_{x}\sim k_{y}\sim k_{E_{\mathrm{max}}}/\sqrt{2}$, compared
to the distance in two-dimensional simulations with $k_{y}\sim
k_{E_{\mathrm{max}}}$.

The analytical expression accurately predicts the dependency of the growth
from the wakefield-driven filamentation instability and the damping from
diffusion. The theory also verifies that the growth of the filamentation
instability can be effectively modelled in two dimensions at a lower in-plane
wavenumber without losing generality.

Figure 6:  The spectral growth dependency on the bunch temperature. The
theoretical growth for different bunch temperatures at
$\omega_{\beta}\tau=2.6$ and $k_{p}\zeta=-12\pi$ (solid lines), compared to 2D
electromagnetic (dashed lines) and 3D quasi-static (dotted lines) PIC
simulations.

The expected distance between filaments,
$\lambda_{f}=2\pi/k_{E_{\mathrm{max}}}$, is shown in Fig. 7 as a function of
the spatiotemporal growth and damping from diffusion. At the back of the
bunch, where filamentation is strongest, the expected distance between
filaments is independent of the bunch length, depending instead on the total
bunch charge, $\omega_{\beta}^{2}\zeta\sim\int
n_{b}\mathop{}\\!\mathrm{d}\zeta$. While the theory developed in this work
considers a longitudinally flat-top beam, this general dependence can readily
be applied to bunches with arbitrary longitudinal profiles.

In experiments carried out with both proton [21] and electron [22] bunches,
the onset of filamentation was studied by varying the plasma density. Taking
the proton bunch parameters 111$400\,\mathrm{GeV}$ proton bunch with a total
charge of $43\,\mathrm{nC}$, an rms width $\sigma_{r}=0.5\,\mathrm{mm}$, a
normalised emittance of $2.5\,\mathrm{mm}\,\mathrm{mrad}$, and a
longitudinally Gaussian profile with $\sigma_{\zeta}/c=163\,\mathrm{ps}$
($\int\omega_{\beta}^{2}\mathop{}\\!\mathrm{d}\zeta=(2\pi)^{1/2}\omega_{\beta}^{2}\sigma_{\zeta}$).
The plasma length was $c\tau=10\,\mathrm{m}$ in [21] and varying the plasma
density gives the dashed line in Fig. 7. Point (a) corresponds to a plasma
density $n_{p}=9.38\times 10^{14}\,\mathrm{cm}^{-3}$, for which filamentation
was observed. The predicted distance between filaments,
$\lambda_{f}=2/k_{p}=0.34\,\mathrm{mm}$, is comparable to the observed
distance of $0.27\,\mathrm{mm}$. Taking the electron bunch parameters
222$0.06\,\mathrm{GeV}$ electron bunch with a total charge of
$1\,\mathrm{nC}$, an rms width $\sigma_{r}=0.065\,\mathrm{mm}$, a normalised
emittance of $6\,\mathrm{mm}\,\mathrm{mrad}$, and an rms length
$\sigma_{\zeta}/c=5\,\mathrm{ps}$. The plasma length was
$c\tau=0.02\,\mathrm{m}$ in [22] and varying the plasma density gives the
dotted line in Fig. 7. Point (b) corresponds to a plasma density
$n_{p}=12\times 10^{16}\,\mathrm{cm}^{-3}$, for which filamentation was
observed. The predicted distance between filaments,
$\lambda_{f}=2.7/k_{p}=0.042\,\mathrm{mm}$, appears to give agreement with the
observed filamentation distance, although the bunch envelope observed in the
experiment was significantly modified through its interaction with the plasma.

The points (c) and (d) correspond to the cases in [21], [22] where a low
plasma density was used, and filamentation was suppressed. For (c),
approximately $50\,$% of shots led to filamentation, suggesting this the
threshold for the instability. For (d), no filamentation was observed. This
threshold for filamentation correlates with the predicted distance between
filaments exceeding the rms width of the bunch. The points ($\alpha$) and
($\beta$) correspond to the cases in [21], [22], where the predicted distance
between filaments is equal to the rms bunch width. Point (c), with a plasma
density of $n_{p}=2.25\times 10^{14}\,\mathrm{cm}^{-3}$, is close to
($\alpha$), with a plasma density of $2.44\times 10^{14}\,\mathrm{cm}^{-3}$.
Point (d), with a plasma density of $n_{p}=1.6\times
10^{16}\,\mathrm{cm}^{-3}$ is well below point ($\beta$), with a plasma
density of $n_{p}=3.4\times 10^{16}\,\mathrm{cm}^{-3}$. The distance between
filaments for the instability cutoff, $2\pi/k_{\mathrm{cut}}$, corresponds to
a plasma density 50–140 times lower than the observed threshold. This
dependence of the instability threshold on $k_{E_{\mathrm{max}}}$ and not
$k_{\mathrm{cut}}$ may be due to competition of the filamentation instability
with SMI of the charged bunches used in these experiments. Further
experimental and numerical studies would allow this prediction for the
instability threshold to be tested across a larger parameter space.

Figure 7: The distance between filaments as a function of bunch and plasma
parameters. The y and x-axis are proportional to the growth of TTS and damping
due to diffusion. The dashed line corresponds to the experimental parameters
in [21] for varying plasma density, while the dotted line corresponds to the
experimental parameters in [22]. Points a, b, c, and d correspond to
individual measurements in [21], [22], with $\alpha$, $\beta$ marking the
point at which the distance between filaments is predicted to reach the rms
bunch width.

##  V Conclusion

A three-dimensional, spatiotemporal theory for the wakefield-driven
filamentation instability is presented for warm bunches of finite size. The
weakly and strongly relativistic regimes, referred to as TSI and TTS, arise
from the longitudinal and transverse wakefield components. In the limit of a
cold stream, the analytical expressions for TSI and TTS simplify to previous
works. The electrostatic plasma response leads to the growth of transverse
filaments with an additional longitudinal modulation. The transverse bunch
profile influences both the growth rate and the seed level, with the growth
rate at a fixed transverse position being equivalent to a stream with the
local bunch density.

For beams with finite emittance, diffusion acts to damp small-scale
filamentation. The dependency of the growth spectrum on the temperature is
identified for dilute bunches. Theory and simulations show that the
filamentation growth rate depends on $k_{r}=(k_{x}^{2}+k_{y}^{2})^{1/2}$. Two-
dimensional simulations reproduce the behaviour of three-dimensional
simulations in the linear regime, with the caveat that $k_{r}=k_{y}$ in this
reduced geometry, resulting in filaments that are more tightly clustered.
Explicit expressions for the dominant and cut-off wavenumber are calculated
and depend on the propagation time in plasma and position within the bunch.
This arises as diffusion is spatially uniform while the filamentation
instability grows along the bunch length. Remarkable agreement is found
between theory and PIC simulations.

Although the analytical treatment developed here considers a longitudinally
flat-top beam, a general dependence on the expected distance between filaments
is found for bunches with arbitrary profile. The predicted distance between
filaments gives good agreement with previously published experimental results.
For single-species beams, filamentation appears to be suppressed when the
predicted distance between filaments is larger than the rms beam width. These
findings provide a crucial basis for designing laboratory astrophysics
experiments investigating filamentation instabilities and for PWFA experiments
seeking to avoid them.

###### Acknowledgements.

We thank Patric Muggli for helpful discussions related to this work. Funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany´s Excellence Strategy – EXC 2094 – 390783311.

##  Appendix A Derivation of Wakefield-Driven Two-Stream Growth for Warm
Bunches

###  A.1 Wakefield Induced by a Modulated Bunch

A dilute bunch propagating in the $+z$ direction through an unmagnetised
plasma leads to an electrostatic plasma response [41, 40, 30, 27]. The
associated fields are $E_{z}$,
$\bm{W}_{\perp}=\bm{E}_{\perp}+u_{b}\bm{\hat{z}}\times\bm{B}_{\perp}$, with
$\bm{\hat{z}}$ the unit vector along z and $u_{b}$ the bunch velocity. Only
the oscillatory plasma current $\bm{j}_{p}$ is considered, and the small bulk
return current for underdense bunches is neglected. Therefore, Ohm’s law
reduces to $\mu_{0}\partial_{t}\bm{j}_{p}=-k_{p}^{2}\bm{E}$, with $\mu_{0}$
the vacuum permeability and $k_{p}$ the plasma wavenumber. The fields can then
be described by [40]

$\displaystyle\left(\nabla^{2}-\partial_{t}^{2}/c^{2}-k_{p}^{2}\right)\bm{E}$
$\displaystyle=\mu_{0}\partial_{t}\bm{j}_{b}+\nabla(\rho_{b}+\delta\rho_{p})/\varepsilon_{0}$
(10)
$\displaystyle\left(\nabla^{2}-\partial_{t}^{2}/c^{2}-k_{p}^{2}\right)\bm{B}$
$\displaystyle=-\mu_{0}\nabla\times\bm{j}_{b},$

with $\rho_{b}$ the bunch charge density, $\delta\rho_{p}$ the charge density
of the plasma perturbation, $\bm{j}_{b}=\rho_{b}u_{b}\bm{\hat{z}}$ the bunch
current density, $c$ the velocity of light and $\varepsilon_{0}$ the vacuum
permittivity. The plasma perturbation connects to the bunch charge density via
its fluid equation,
$m_{e}\partial_{t}^{2}\delta\rho_{p}=-e^{2}\nabla\bm{E}=-e^{2}(\rho_{b}+\delta\rho_{p})/\varepsilon_{0}$,
where $m_{e}$ is the electron mass.

With bunch slice $\zeta$ and propagation time in plasma $\tau$, the Lagrangian
frame of the bunch is defined by $\zeta=z-u_{b}t,\,\tau=z/u_{b}$. The partial
derivatives transform in the bunch frame to
$\partial_{t}=-u_{b}\partial_{\zeta}$ and
$\partial_{z}=\partial_{\zeta}+\partial_{\tau}/u_{b}$. Assuming that the bunch
evolution along its propagation $z$ is significantly slower than the response
of the plasma electrons along the bunch $\zeta$, the quasi-static
approximation for the plasma quantities and wakefield can be assumed.
Therefore,
$|\partial_{\zeta}\delta\rho_{p}|\gg|\partial_{\tau}\delta\rho_{p}/u_{b}|$ and
$(\partial_{z}^{2}-\partial_{t}^{2}/c^{2})\delta\rho_{p}\to(1-u_{b}^{2}/c^{2})\partial_{\zeta}\delta\rho_{p}=\partial_{\zeta}^{2}\delta\rho_{p}/\gamma_{b}^{2}$,
with $\gamma_{b}=(1-u_{b}^{2})^{-1/2}$ the Lorentz factor. The same applies
for $\bm{E}$ and $\bm{B}$.

Defining the 3D Fourier transform

$\displaystyle\hat{\rho}_{b}$ $\displaystyle=\mathcal{F}_{\zeta
xy}\\{\rho_{b}\\}(k_{\zeta},k_{x},k_{y})$ (11)
$\displaystyle=\iiint_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}\zeta\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y\rho_{b}\exp\left(-{i\mkern
1.0mu}k_{\zeta}\zeta-{i\mkern 1.0mu}k_{x}x-{i\mkern 1.0mu}k_{y}y\right),$

the spectral form of the plasma fluid in the bunch frame is given by
$\delta\hat{\rho}_{p}=-k_{e}^{2}\hat{\rho}_{b}/(k_{\zeta}^{2}+k_{e}^{2})$,
with $k_{e}=ck_{p}/u_{b}$. The field components transform to

$\displaystyle\hat{E}_{z}$ $\displaystyle=-\frac{{i\mkern
1.0mu}}{\varepsilon_{0}\sqrt{2\pi}}\frac{k_{\zeta}(k_{\zeta}^{2}/\gamma_{b}^{2}+k_{p}^{2})\mathcal{F}_{\zeta
xy}\\{\rho_{b}\\}}{(k_{\zeta}^{2}-k_{e}^{2})(k_{\zeta}^{2}/\gamma_{b}^{2}+k_{p}^{2}+k_{r}^{2})}$
(12) $\displaystyle\bm{\hat{E}}_{\perp}$
$\displaystyle=\frac{1}{\varepsilon_{0}\sqrt{2\pi}}\frac{k_{\zeta}^{2}\mathcal{F}_{\zeta
xy}\\{\nabla_{\perp}\rho_{b}\\}}{(k_{\zeta}^{2}-k_{e}^{2})(k_{\zeta}^{2}/\gamma_{b}^{2}+k_{p}^{2}+k_{r}^{2})}$
$\displaystyle\bm{\hat{B}}_{\perp}$
$\displaystyle=\frac{u_{b}/c^{2}}{\varepsilon_{0}\sqrt{2\pi}}\frac{\mathcal{F}_{\zeta
xy}\\{\nabla^{\perp}\rho_{b}\\}}{k_{\zeta}^{2}/\gamma_{b}^{2}+k_{p}^{2}+k_{r}^{2}}.$

with $k_{r}=(k_{x}^{2}+k_{y}^{2})^{1/2}$,
$\nabla_{\perp}=(\partial_{x},\partial_{y})$ and
$\nabla^{\perp}=(-\partial_{y},\partial_{x})$.

To obtain the fields for a small-scale perturbation in 3D configuration space,
a quasi-neutral bunch (equal populations of particles with opposite charge) is
superimposed by a non-neutral transverse modulation
$g(x,y)=\tilde{g}(x,y)\cos{(k_{x}x+\varphi_{x})}\cos{(k_{y}y+\varphi_{y})}$,
with $\tilde{g}$ the slowly varying transverse envelope, i.e.
$|\partial_{y}\tilde{g}|\ll k_{y}|\tilde{g}|$, and $k_{x,y}$ and
$\varphi_{x,y}$ the respective modulation wavenumbers and phases. The positron
density may be given by
$n_{bp}=[n_{b}/2]\tilde{f}(\zeta)\tilde{g}(x,y)+[\delta
n_{b}/\sqrt{2}]f(\zeta)\tilde{g}(x,y)\cos{(k_{x}x)}\cos{(k_{y}y-\pi/4)}$ and
respectively for the electron density
$n_{be}=[n_{b}/2]\tilde{f}(\zeta)\tilde{g}(x,y)+[\delta
n_{b}/\sqrt{2}]f(\zeta)\tilde{g}(x,y)\cos{(k_{x}x)}\sin{(k_{y}y-\pi/4)}$, with
$n_{b}$ the total density amplitude of the bunch and $\delta n_{b}$ the
amplitude of the density perturbation. The longitudinal bunch shape and its
slowly varying envelope are given by $f(\zeta)$ and $\tilde{f}(\zeta)$. The
net charge density of the bunch, $\rho_{b}=q\delta n_{b}f(\zeta)g(x,y)$, with
$q$ the charge, serves as the source for the fields. The inverse Fourier
transforms for $\zeta<0$ are

$\displaystyle\hat{E}_{z}$ $\displaystyle=\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{\mathcal{F}_{xy}\\{g(x,y)\\}}{k_{e}^{2}+k_{r}^{2}}\int_{-\infty}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})\left[k_{e}^{2}\cos{k_{e}(\zeta-\zeta^{\prime})}+k_{r}^{2}\exp\left(-\gamma_{b}\sqrt{k_{p}^{2}+k_{r}^{2}}|\zeta-\zeta^{\prime}|\right)\right]$
(13) $\displaystyle\bm{\hat{E}}_{\perp}$ $\displaystyle=\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{\mathcal{F}_{xy}\\{\nabla_{\perp}g(x,y)\\}}{k_{e}^{2}+k_{r}^{2}}\int_{-\infty}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})\left[k_{e}\sin{k_{e}(\zeta-\zeta^{\prime})}-\gamma_{b}\sqrt{k_{p}^{2}+k_{r}^{2}}\exp\left(-\gamma_{b}\sqrt{k_{p}^{2}+k_{r}^{2}}|\zeta-\zeta^{\prime}|\right)\right]$
$\displaystyle\bm{\hat{B}}_{\perp}$
$\displaystyle=\frac{\gamma_{b}u_{b}}{c^{2}}\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{\mathcal{F}_{xy}\\{\nabla^{\perp}g(x,y)\\}}{\sqrt{k_{p}^{2}+k_{r}^{2}}}\int_{-\infty}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})\exp\left(-\gamma_{b}\sqrt{k_{p}^{2}+k_{r}^{2}}|\zeta-\zeta^{\prime}|\right).$

Neglecting the small spectral broadening due to $\tilde{g}(x,y)$, the
transverse inverse Fourier transform for the transverse component of the
electric field gives

$\bm{E}_{\perp}\sim\mathcal{F}_{xy}^{-1}\left\\{\frac{\mathcal{F}_{xy}\left\\{\nabla_{\perp}g(x,y)\right\\}}{k_{e}^{2}+k_{r}^{2}}\right\\}\approx\frac{\nabla_{\perp}g(x,y)}{k_{e}^{2}+k_{r}^{2}}$
(14)

and for the $z$ component gives $E_{z}\sim-k_{e}/(k_{e}^{2}+k_{r}^{2})g(x,y)$.

The second electromagnetic summand in the integral can be split into the
contribution of the local bunch slice and the inductive, purely decaying
fields due to a change in bunch shape. The latter can be safely ignored if the
plasma is non-diffusive [40]. The electromagnetic terms simplify to

$\int_{\zeta}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})\exp\left(-\gamma_{b}\sqrt{k_{p}^{2}+k_{r}^{2}}|\zeta-\zeta^{\prime}|\right)\approx\frac{f(\zeta)}{\gamma_{b}\sqrt{k_{p}^{2}+k_{r}^{2}}}.$
(15)

Without any limitation on the longitudinal shape, the fields can be expressed
by

$\displaystyle E_{z}$ $\displaystyle=\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{k_{e}g(x,y)}{k_{e}^{2}+k_{r}^{2}}\int_{\zeta}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})k_{e}\cos{k_{e}(\zeta-\zeta^{\prime})}$
(16) $\displaystyle\bm{E}_{\perp}$ $\displaystyle=\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{\nabla_{\perp}g(x,y)}{k_{e}^{2}+k_{r}^{2}}$
$\displaystyle\times\left[\int_{\zeta}^{0}\mathop{}\\!\mathrm{d}\zeta^{\prime}f(\zeta^{\prime})k_{e}\sin{k_{e}(\zeta-\zeta^{\prime})}-f(\zeta)\right]$
$\displaystyle\bm{B}_{\perp}$ $\displaystyle=-\frac{u_{b}}{c^{2}}\frac{q\delta
n_{b}}{\varepsilon_{0}}\frac{\nabla^{\perp}g(x,y)}{k_{p}^{2}+k_{r}^{2}}f(\zeta).$

For relativistic bunches $k_{e}\approx k_{p}$, the latter charge-repulsion
term in $\bm{E}_{\perp}$ and the magnetic field $\bm{B}_{\perp}$ approximate
to the local bunch contribution $W_{f}\sim
f(\zeta)(1-u_{b}^{2})=f(\zeta)/\gamma_{b}^{2}$ and are usually neglected for
$\gamma_{b}\gg 1$. Each bunch slice $f(\zeta)$ drives a wakefield with an
amplitude proportional to the transverse bunch shape. These contributions sum
up along the bunch.

###  A.2 Growth of Two-Stream Filamentation

The excited fields act on the bunch. Assuming a cold bunch with a longitudinal
momentum much larger than its transverse momentum, the linearised fluid
equation gives

$\partial_{\tau}^{2}\delta
n_{b}=\left(\partial_{t}+u_{b}\partial_{z}\right)^{2}\delta
n_{b}=\frac{2\omega_{\beta}^{2}}{q/\varepsilon_{0}}\left(\frac{\partial_{z}E_{z}}{\gamma_{b}^{2}}+\nabla_{\perp}\cdot\bm{W}_{\perp}\right),$
(17)

with $\omega_{\beta}=[q^{2}n_{b}/(2\gamma_{b}\varepsilon_{0}m_{b})]^{1/2}$ the
betatron frequency and $m_{b}$ the mass of bunch particles.

The fields result in positive feedback, which gives rise to spatiotemporal
growth. The growth due to the wakefield, given by the integral terms in
$E_{z}$ and $\bm{E}_{\perp}$, can be tracked by applying the spatial
derivative along $\zeta$ to Eq. 17. In the strongly coupled regime,
$\omega_{\beta}\tau\ll k_{e}\zeta$, the bunch perturbation can be described by
$\delta n_{b}g(x,y)\tilde{f}(\zeta)[\exp{({i\mkern
1.0mu}k_{e}\zeta)}/2+\mathrm{c.c.}]$, considering the longitudinal wavenumber
of the wakefields at $k_{\zeta}=k_{e}$. The integral along $\zeta$ from Eq. 16
reduces to

$\partial_{\zeta}\int_{\zeta}^{0}\tilde{f}(\zeta^{\prime})\frac{\exp{({i\mkern
1.0mu}k_{e}\zeta^{\prime})}}{2}\sin{k_{e}(\zeta-\zeta^{\prime})}\approx\frac{{i\mkern
1.0mu}}{2}\tilde{f}(\zeta)\exp{({i\mkern 1.0mu}k_{e}\zeta)}.$ (18)

For a flat-top bunch, $\tilde{f}(\zeta)=\Theta(-\zeta)$, the initial
perturbation is given by $\delta n_{b}(\tau=0,\zeta)=\delta
n_{b0}\Theta(-\zeta)$. The local terms in Eq. 16, which only act within a
bunch slice $E_{x,y}\sim f(\zeta)$ and $B_{x,y}\sim u_{b}f(\zeta)$, are
negligible compared to the growing wakefield term. For a slowly varying
transverse envelope, the transverse gradient simplifies to
$\nabla_{\perp}^{2}g(x,y)\approx-k_{r}^{2}g(x,y)$ and the perturbation
amplitude follows [18]

$\begin{gathered}\left[\partial_{\zeta}\partial_{\tau}^{2}+{i\mkern
1.0mu}\eta_{u}k_{e}\omega_{\beta}^{2}\tilde{g}(x,y)\right]\delta
n_{b}(\tau,\zeta)=0\\\
\eta_{u}=\frac{(c^{2}-u_{b}^{2})k_{p}^{2}+u_{b}^{2}k_{r}^{2}}{c^{2}k_{p}^{2}+u_{b}^{2}k_{r}^{2}}.\end{gathered}$
(19)

The spectral parameter $\eta_{u}$ includes the dependency on the bunch
velocity, representing the relative contribution of the (longitudinal) two-
stream and transverse two-stream instability. The Green’s function can be
solved by a double Laplace transform [18]

$\displaystyle\mathcal{L_{\zeta\tau}}\\{\delta
n_{b}\\}(k_{\zeta},k_{\tau})=\iint_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}\tau\mathop{}\\!\mathrm{d}\zeta\delta
n_{b}\exp\left(-k_{\zeta}\zeta-k_{\tau}\tau\right)$ (20)
$\displaystyle=\frac{k_{\zeta}\mathcal{L}_{\zeta}\left\\{(k_{\tau}+\partial_{\tau})\delta
n_{b}(\tau=0,\zeta)\right\\}+\mathcal{L}_{\tau}\left\\{\partial_{\tau}^{2}\delta
n_{b}(\tau,\zeta=0)\right\\}}{k_{\zeta}k_{\tau}^{2}+{i\mkern
1.0mu}\eta_{u}k_{e}\omega_{\beta}^{2}\tilde{g}(x,y)}.$

Assuming a sharp plasma boundary at $\tau=0$ sets the initial condition
$\delta n_{b}(\tau,\zeta=0)=\delta n_{b0}$, which results in
$\partial_{\tau}\delta n_{b}(\tau=0,\zeta)=\partial_{\tau}^{2}\delta
n_{b}(\tau,\zeta=0)=0$ and $\mathcal{L}_{\zeta}\\{\delta
n_{b}(\tau=0,\zeta)\\}=\delta n_{b0}/k_{\zeta}$. Using the Residue theorem for
the inverse Laplace transform in $\zeta$ and the relation
$\mathcal{L}_{\tau}^{-1}\\{\tau^{-2n-1}\\}=t^{2n}/(2n)!$ [54] for the inverse
transform in $\tau$ gives the solution to LABEL:eq:app_lapl_dn_obi as a
complex power series for $\tau\geq 0,\zeta\leq 0$

$\delta n_{b,\mathrm{TS}}=\delta
n_{b0}\sum_{n=0}^{\infty}\frac{\left[i\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2}\right]^{n}}{n!(2n)!}.$
(21)

The solution contains a growing imaginary and oscillatory real term, which can
be obtained by the absolute value $\Gamma_{\mathrm{TS}}=|\delta n_{b}/\delta
n_{b0}|$ and the phase $\psi(\delta n_{b})$. The asymptotic expansion,
$\tau\rightarrow\infty$, to Eq. 21 gives

$\displaystyle\delta n_{b,\mathrm{TS}}$ $\displaystyle\approx\frac{\delta
n_{b0}}{\sqrt{4\pi}}\frac{\exp\\{(3/2^{2/3})[{i\mkern
1.0mu}\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2}]^{1/3}\\}}{\sqrt{(3/2^{2/3})[i\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2}]^{1/3}}}.$
(22)

The growth of the bunch perturbation due to the combined two-stream
instabilities is

$\displaystyle\delta n_{b,\mathrm{TS}}$ $\displaystyle\approx\frac{\delta
n_{b0}}{\sqrt{4\pi}}\frac{\exp{\Gamma_{\infty}}}{\sqrt{\Gamma_{\infty}}}$ (23)
$\displaystyle\Gamma_{\infty}$
$\displaystyle=\frac{3^{3/2}}{2^{5/3}}\left[\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2}\right]^{1/3}.$

The oscillatory term yields a phase $\psi$ and corresponding phase velocity
$u_{\psi}=-\partial_{t}\psi/\partial_{z}\psi$ of the growing wave

$\displaystyle\psi$
$\displaystyle=\frac{\pi}{4}-k_{e}|\zeta|-\frac{3}{2^{5/3}}(\eta_{u}\tilde{g}(x,y)k_{e}|\zeta|\omega_{\beta}^{2}\tau^{2})^{1/3}$
(24) $\displaystyle u_{\psi}$
$\displaystyle=u_{b}\left[1-\frac{1}{2^{2/3}}\left(\eta_{u}\tilde{g}(x,y)\frac{\omega_{\beta}^{2}}{\omega_{p}^{2}}\frac{|\zeta|}{c\tau}\right)^{1/3}\right],$

comparable to the phase velocity from SMI [18, 49].

At early propagation times, the initial bunch perturbation and, consequently,
the plasma perturbation is much larger than the exponential growth of the two-
stream instability. For short times, the growth evolves as

$\delta n_{b,S}=\delta
n_{b0}\left[\eta_{u}\tilde{g}(x,y)\omega_{\beta}^{2}\tau^{2}+1\right].$ (25)

This can be observed in Fig. 8, where the bunch and plasma perturbation are
not purely exponential. The same initial field dominates their initial growth,
and the exponentially growing term only dominates after the bunch has
propagated for some time. However, the transverse electric field, being the
difference between plasma and bunch charge density

$\frac{\bm{E}_{\perp}}{E_{0}}=\frac{\bm{k}_{r}}{k_{e}^{2}+k_{r}^{2}}\frac{e\delta
n_{p}-q\delta n_{b}}{n_{p}}$ (26)

exhibits exponential growth even at early times.

Figure 8:  The growth of the spectral a) plasma perturbation, b) bunch
perturbation and c) electric field at $k_{p}\zeta=-12\pi$ from simulation and
theory for the cold bunch in Fig. 4. The yellow and green dashed lines show
the asymptotic (23) and semi-analytic (21) solution, and the orange line
includes the growth from the initial fields. The dotted lines indicate the
seed for the two-stream instability.

###  A.3 Influence of Diffusion

Extending Eq. 17 for warm bunches with a thermal spread $\sigma_{pr}$ requires
the pressure term $\mathcal{P}$ to be included in Eq. 17

$\displaystyle\partial_{\tau}^{2}\delta
n_{b}=\frac{2\omega_{\beta}}{q/\varepsilon_{0}}\left(\frac{\partial_{z}E_{z}}{\gamma_{b}^{2}}+\nabla_{\perp}\cdot\bm{W}_{\perp}\right)+\frac{\nabla^{2}\mathcal{P}}{\gamma_{b}m_{b}},$
(27)

where the pressure can be described by $\mathcal{P}=(2/3)\sigma_{pr}^{2}\delta
n_{b}/(\gamma_{b}m_{b})$ for non-relativistic temperatures,
$\sigma_{pr}^{2}/(m_{b}c)^{2}\ll 1$ [51]. The thermal spread can be related to
the normalised emittance $\epsilon_{N}$ by
$\sigma_{pr}/(m_{b}c)=\epsilon_{N}/\sigma_{r}$, where $\sigma_{r}$ is the rms
width for a Gaussian bunch. The effect of emittance-related diffusion is
purely temporal. It can be considered separately from the wakefield-driven
two-stream instability since all bunch slices are equally affected by the
bunch divergence, and the fluid equation reduces to

$\left[\partial_{\tau}^{2}-\frac{2}{3}\frac{\sigma_{pr}^{2}}{m_{b}^{2}\gamma_{b}^{2}}\nabla_{\perp}^{2}\right]\delta
n_{b}g(x,y)=0.$ (28)

Considering the slowly varying envelope in $g(x,y)$, the damping of the
perturbation amplitude is described by

$\left[\partial_{\tau}^{2}+\frac{2}{3}\frac{\sigma_{pr}^{2}k_{r}^{2}}{m_{b}^{2}\gamma_{b}^{2}}\right]\delta
n_{b}=0.$ (29)

The Green’s function is readily obtained by a Fourier transform for $\tau\geq
0$ to

$\begin{gathered}\delta n_{b,D}=\delta n_{b}\exp{(-\delta_{D}\tau)}\\\
\delta_{D}=\sqrt{\frac{2}{3}}\frac{\sigma_{pr}k_{r}}{\gamma_{b}m_{b}}=\sqrt{\frac{2}{3}}\frac{\sigma_{pr}/(m_{b}c)}{\gamma_{b}}\frac{k_{r}}{k_{p}}\omega_{p}.\end{gathered}$
(30)

Consequently, the total growth rate of the bunch perturbation is a sum of the
growth rate from the two-stream instability with the damping rate from
diffusion,
$\Gamma_{\mathrm{tot}}=\Gamma_{\mathrm{TS}}\exp{(-\delta_{D}\tau)}$.

The growth of the two-stream instability is larger for higher wavenumbers, as
seen by the spectral parameter $\eta_{u}$ in Eq. 19. However, these
wavenumbers are more strongly damped by diffusion. This gives rise to a finite
wavenumber $k_{\Gamma_{\mathrm{max}}}(\tau,\zeta)$ for which the growth is
largest, which can be derived in the asymptotic limit by
$\partial_{k_{r}}\exp{(\Gamma_{\infty}-\delta_{D}\tau)}/(4\pi\Gamma_{\infty})^{1/2}=0$.
Further, a cut-off wavenumber $k_{\mathrm{cut}}(\tau,\zeta)$ exists at which
the growth and damping rates are equal,
$\exp{(\Gamma_{\infty}-\delta_{D}\tau)}/(4\pi\Gamma_{\infty})^{1/2}=0$. For
higher wavenumbers, an initial bunch perturbation will be damped. Their
respective values can be numerically obtained by

$\displaystyle 2\Gamma_{\infty}(k_{\Gamma_{\mathrm{max}}})$
$\displaystyle=3(1+k_{\Gamma_{\mathrm{max}}}^{2})\delta_{D}(k_{\Gamma_{\mathrm{max}}})\tau+1$
(31) $\displaystyle\Gamma_{\infty}(k_{\mathrm{cut}})$
$\displaystyle=\delta_{D}(k_{\mathrm{cut}})\tau+\ln{\sqrt{4\pi\Gamma_{\infty}(k_{\mathrm{cut}})}}$

###  A.4 Transition from TTS to TSI

To qualitatively compare the dominant regime of the (longitudinal) two-stream
and transverse two-stream instability, respectively referred to as TSI and
TTS, the spectral parameter from Eq. 19 can be rewritten to
$\eta_{u}=\eta_{\mathrm{TSI}}+\eta_{\mathrm{TTS}}$. The longitudinal and
transverse contributions are provided by

$\eta_{\mathrm{TSI}}=\frac{(c^{2}-u_{b}^{2})k_{p}^{2}}{c^{2}k_{p}^{2}+u_{b}^{2}k_{r}^{2}},\quad\eta_{\mathrm{TTS}}=\frac{u_{b}^{2}k_{r}^{2}}{c^{2}k_{p}^{2}+u_{b}^{2}k_{r}^{2}},$
(32)

and shown in Fig. 9a) and b).

Figure 9:  The spectral dependency $\eta_{u}$ for a) TSI and b) TTS on bunch
velocity $u_{b}$ and wavenumber $k_{r}$. The purple line indicates where TSI
and TTS equally contribute.

As expected, TSI is dominant for non-relativistic bunches, and the
longitudinal wakefield component predominantly modulates the bunch. However,
for transverse perturbations with a long scale, $k_{r}/k_{p}<1$, TSI remains
dominant even in mildly relativistic regimes. This is a consequence of the
transverse electric field scaling
$\tilde{\bm{E}}_{\perp}=\tilde{E_{z}}\bm{k}_{r}u_{b}/(k_{p}c)$ to the
longitudinal field from Eq. 16. TTS is dominant for highly relativistic
bunches or high transverse wavenumbers in mildly relativistic bunches, such
that the transverse wakefield predominantly modulates the bunch. Given a
negligible energy spread of the bunch, the longitudinal wavenumber of the two-
stream instability uniformly equals $k_{\zeta}=k_{e}=ck_{p}/u_{b}$. The
combined influence of TSI and TTS is generally referred to as oblique
instability (OBI) [25, 26, 55, 56, 2, 57, 46, 39]. However, the current
filamentation instability (CFI), which becomes dominant for overdense beams,
represents a different longitudinal wavenumber ($k_{\zeta}=0$) and growth
scaling as discussed in [58, 35].

Figure 10:  The longitudinal $E_{z}$ (green) and transverse $E_{y}$ (red)
electric field envelope for $u_{b}=0.25\,c$ and $0.999\,c$ as a function of a)
propagation time and b) position along the bunch. The solid vertical lines
indicate the theoretical ratio between the field components. The dashed lines
indicate the respective slice taken.

Figure 10 shows the growth of an initial perturbation $k_{r}/k_{p}=\pi$ for
two different bunch velocities. As can be seen, the growth scales with
$\omega_{\beta}\tau$ and $k_{e}\tau$, in agreement with Eq. 23, as the
spectral parameter $\eta_{u}$ remains roughly constant between non-
relativistic and relativistic bunches for $k_{r}/k_{p}\gtrsim 3$. For a
constant wavenumber, $E_{y}$ is weaker in the non-relativistic limit, given by
the theoretical ratio.

Bunches with reduced mass are often used to lower the computational overhead
of simulations. It should be noted that the two-stream instability growth
scales with $\omega_{\beta}\sim m_{b}^{-1/2}$ along the propagation time while
damping scales with $\sigma_{pr}/m_{b}$. Therefore, when scaling the bunch
mass, the bunch thermal spread should be scaled by a factor of
$(m_{b}/m_{\mathrm{reduced}})^{1/2}$ to maintain the ratio of the growth and
diffusion rate.

## References

  * Chen [2016] F. F. Chen, in  _Introduction to Plasma Physics and Controlled Fusion_ (Springer Interational Publishing Switzerland 2016, 2016) Chap. 7.3, pp. 222–224, 3rd ed.
  * Bret  _et al._ [2010a] A. Bret, L. Gremillet, and M. E. Dieckmann, Physics of Plasmas 17, 120501 (2010a), https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.3514586/16019035/120501_1_online.pdf .
  * Michno and Schlickeiser [2010] M. J. Michno and R. Schlickeiser, The Astrophysical Journal 714, 868 (2010).
  * Spitkovsky [2008] A. Spitkovsky, The Astrophysical Journal 682, L5 (2008).
  * Landau [1946] L. Landau, J. Phys. 10, 25 (1946).
  * Iwamoto  _et al._ [2019] M. Iwamoto, T. Amano, M. Hoshino, Y. Matsumoto, J. Niemiec, A. Ligorini, O. Kobzar, and M. Pohl, The Astrophysical Journal Letters 883, L35 (2019).
  * Tajima  _et al._ [2020] T. Tajima, X. Q. Yan, and T. Ebisuzaki, Reviews of Modern Plasma Physics 4, 10.1007/s41614-020-0043-z (2020).
  * Hillas [1984] A. M. Hillas, Annual Review of Astronomy and Astrophysics 22, 425 (1984), https://doi.org/10.1146/annurev.aa.22.090184.002233 .
  * Piron [2016] F. Piron, Comptes Rendus Physique 17, 617 (2016), gamma-ray astronomy / Astronomie des rayons gamma - Volume 2\.
  * Bohdan  _et al._ [2021] A. Bohdan, M. Pohl, J. Niemiec, P. J. Morris, Y. Matsumoto, T. Amano, M. Hoshino, and A. Sulaiman, Phys. Rev. Lett. 126, 095101 (2021).
  * Medvedev and Loeb [1999] M. V. Medvedev and A. Loeb, The Astrophysical Journal 526, 697 (1999).
  * Perna  _et al._ [2016] R. Perna, D. Lazzati, and B. Giacomazzo, The Astrophysical Journal Letters 821, L18 (2016).
  * Fiuza  _et al._ [2020] F. Fiuza, G. Swadling, A. Grassi, and et al., Nat. Phys. 16, 916 (2020).
  * Arrowsmith  _et al._ [2021] C. D. Arrowsmith, N. Shukla, N. Charitonidis, R. Boni, H. Chen, T. Davenne, A. Dyson, D. H. Froula, J. T. Gudmundsson, B. T. Huffman, Y. Kadi, B. Reville, S. Richardson, S. Sarkar, J. L. Shaw, L. O. Silva, P. Simon, R. M. G. M. Trines, R. Bingham, and G. Gregori, Phys. Rev. Res. 3, 023103 (2021).
  * Zhang  _et al._ [2022] C. Zhang, Y. Wu, M. Sinclair, A. Farrell, K. A. Marsh, I. Petrushina, N. Vafaei-Najafabadi, A. Gaikwad, R. Kupfer, K. Kusche, M. Fedurin, I. Pogorelsky, M. Polyanskiy, C.-K. Huang, J. Hua, W. Lu, W. B. Mori, and C. Joshi, Proceedings of the National Academy of Sciences 119, e2211713119 (2022), https://www.pnas.org/doi/pdf/10.1073/pnas.2211713119 .
  * Chen  _et al._ [1985] P. Chen, J. M. Dawson, R. W. Huff, and T. Katsouleas, Phys. Rev. Lett. 55, 1537 (1985).
  * Macchi and Pegoraro [2018] A. Macchi and F. Pegoraro, Nature Photonics 12, 314 (2018).
  * Schroeder  _et al._ [2011] C. B. Schroeder, C. Benedetti, E. Esarey, F. J. Grüner, and W. P. Leemans, Phys. Rev. Lett. 107, 145002 (2011).
  * Caldwell and Lotov [2011] A. Caldwell and K. V. Lotov, Physics of Plasmas 18, 103101 (2011), https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.3641973/13613329/103101_1_online.pdf .
  * AWAKE Collaboration [2018] AWAKE Collaboration, Nature 561, 363 (2018).
  * Verra  _et al._ [2024] L. Verra, P. Muggli, _et al._ (AWAKE Collaboration), Phys. Rev. E 109, 055203 (2024).
  * Allen  _et al._ [2012] B. Allen, V. Yakimenko, M. Babzien, M. Fedurin, K. Kusche, and P. Muggli, Phys. Rev. Lett. 109, 185007 (2012).
  * Weibel [1959] E. S. Weibel, Phys. Rev. Lett. 2, 83 (1959).
  * Fried [1959] B. D. Fried, The Physics of Fluids 2, 337 (1959), https://pubs.aip.org/aip/pfl/article-pdf/2/3/337/12401439/337_1_online.pdf .
  * Bludman  _et al._ [1960] S. A. Bludman, K. M. Watson, and M. N. Rosenbluth, The Physics of Fluids 3, 747 (1960), https://pubs.aip.org/aip/pfl/article-pdf/3/5/747/12499127/747_1_online.pdf .
  * Faǐnberg  _et al._ [1969] Y. B. Faǐnberg, V. D. Shapiro, and V. I. Shevchenko, Soviet Journal of Experimental and Theoretical Physics 30, 528 (1969).
  * Bret  _et al._ [2004] A. Bret, M.-C. Firpo, and C. Deutsch, Phys. Rev. E 70, 046401 (2004).
  * Tonks and Langmuir [1929] L. Tonks and I. Langmuir, Phys. Rev. 33, 195 (1929).
  * Dawson [1959] J. M. Dawson, Phys. Rev. 113, 383 (1959).
  * Lawson and Lawson [1977] J. Lawson and J. Lawson,  _The Physics of Charged-particle Beams_, International series of monographs on physics (Clarendon Press, 1977).
  * Watson  _et al._ [1960] K. M. Watson, S. A. Bludman, and M. N. Rosenbluth, The Physics of Fluids 3, 741 (1960), https://pubs.aip.org/aip/pfl/article-pdf/3/5/741/12499098/741_1_online.pdf .
  * Davidson  _et al._ [1972] R. C. Davidson, D. A. Hammer, I. Haber, and C. E. Wagner, The Physics of Fluids 15, 317 (1972), https://pubs.aip.org/aip/pfl/article-pdf/15/2/317/12743024/317_1_online.pdf .
  * Silva  _et al._ [2002] L. O. Silva, R. A. Fonseca, J. W. Tonge, W. B. Mori, and J. M. Dawson, Physics of Plasmas 9, 2458 (2002), https://pubs.aip.org/aip/pop/article-pdf/9/6/2458/19097722/2458_1_online.pdf .
  * Jia  _et al._ [2013] Q. Jia, H.-b. Cai, W.-w. Wang, S.-p. Zhu, Z. M. Sheng, and X. T. He, Physics of Plasmas 20, 032113 (2013), https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4796052/14795548/032113_1_online.pdf .
  * Pathak  _et al._ [2015] V. B. Pathak, T. Grismayer, A. Stockem, R. A. Fonseca, and L. O. Silva, New Journal of Physics 17, 043049 (2015).
  * Bret  _et al._ [2010b] A. Bret, L. Gremillet, and D. Bénisti, Phys. Rev. E 81, 036402 (2010b).
  * Bers [1983] A. Bers, Handbook of plasma physics. vol. 1 (North-Holland Publ., Amsterdam, 1983) Chap. 3.2 Space-time evolution of plasma instabilities - Absolute and convective, pp. 451–517.
  * Jones  _et al._ [1983] M. E. Jones, D. S. Lemons, and M. A. Mostrom, The Physics of Fluids 26, 2784 (1983), https://pubs.aip.org/aip/pfl/article-pdf/26/10/2784/12495716/2784_1_online.pdf .
  * San Miguel Claveria  _et al._ [2022] P. San Miguel Claveria, X. Davoine, J. R. Peterson, M. Gilljohann, I. Andriyash, R. Ariniello, C. Clarke, H. Ekerfelt, C. Emma, J. Faure, S. Gessner, M. J. Hogan, C. Joshi, C. H. Keitel, A. Knetsch, O. Kononenko, M. Litos, Y. Mankovska, K. Marsh, A. Matheron, Z. Nie, B. O’Shea, D. Storey, N. Vafaei-Najafabadi, Y. Wu, X. Xu, J. Yan, C. Zhang, M. Tamburini, F. Fiuza, L. Gremillet, and S. Corde, Phys. Rev. Res. 4, 023085 (2022).
  * Keinigs and Jones [1987] R. Keinigs and M. E. Jones, The Physics of Fluids 30, 252 (1987), https://pubs.aip.org/aip/pfl/article-pdf/30/1/252/12365569/252_1_online.pdf .
  * Katsouleas  _et al._ [1987] T. C. Katsouleas, S. Wilks, P. Chen, J. M. Dawson, and J. J. Su, Part. Accel. 22, 81 (1987).
  * Clayton  _et al._ [2016] C. E. Clayton, E. Adli, J. Allen, W. An, C. I. Clarke, S. Corde, J. Frederico, S. Gessner, S. Z. Green, M. J. Hogan, C. Joshi, M. Litos, W. Lu, K. A. Marsh, W. B. Mori, N. Vafaei-Najafabadi, X. Xu, and V. Yakimenko, Nature Communications 7, 12483 (2016).
  * Albert  _et al._ [2021] F. Albert, M. E. Couprie, A. Debus, M. C. Downer, J. Faure, A. Flacco, L. A. Gizzi, T. Grismayer, A. Huebl, C. Joshi, M. Labat, W. P. Leemans, A. R. Maier, S. P. D. Mangles, P. Mason, F. Mathieu, P. Muggli, M. Nishiuchi, J. Osterhoff, P. P. Rajeev, U. Schramm, J. Schreiber, A. G. R. Thomas, J.-L. Vay, M. Vranic, and K. Zeil, New Journal of Physics 23, 031101 (2021).
  * Moreira  _et al._ [2023] M. Moreira, P. Muggli, and J. Vieira, Phys. Rev. Lett. 130, 115001 (2023).
  * Gschwendtner  _et al._ [2022] E. Gschwendtner, K. Lotov, P. Muggli, M. Wing, _et al._ (AWAKE Collaboration), Symmetry 14, 10.3390/sym14081680 (2022).
  * Shukla  _et al._ [2018] N. Shukla, J. Vieira, P. Muggli, G. Sarri, R. Fonseca, and L. O. Silva, Journal of Plasma Physics 84, 905840302 (2018).
  * Pukhov [2016] A. Pukhov, CERN Yellow Reports 10.5170/CERN-2016-001.181 (2016).
  * Pukhov [1999] A. Pukhov, Journal of Plasma Physics 61, 425–433 (1999).
  * Pukhov  _et al._ [2011] A. Pukhov, N. Kumar, T. Tückmantel, A. Upadhyay, K. Lotov, P. Muggli, V. Khudik, C. Siemon, and G. Shvets, Phys. Rev. Lett. 107, 145003 (2011).
  * Fonseca  _et al._ [2002] R. Fonseca, L. Silva, F. Tsung, V. Decyk, W. Lu, C. Ren, W. Mori, S. Deng, S. Lee, T. Katsouleas, and J. Adam, eds.,  _OSIRIS: A Three-Dimensional, Fully Relativistic Particle in Cell Code for Modeling Plasma Based Accelerators_, Computational Science - ICCS 2002 No. 2331 (Springer Berlin Heidelberg, 2002).
  * Bret and Deutsch [2006] A. Bret and C. Deutsch, Physics of Plasmas 13, 042106 (2006), https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.2196876/16011061/042106_1_online.pdf .
  * Note [1] $400\,\mathrm{GeV}$ proton bunch with a total charge of $43\,\mathrm{nC}$, an rms width $\sigma_{r}=0.5\,\mathrm{mm}$, a normalised emittance of $2.5\,\mathrm{mm}\,\mathrm{mrad}$, and a longitudinally Gaussian profile with $\sigma_{\zeta}/c=163\,\mathrm{ps}$ ($\intop\ilimits@\omega_{\beta}^{2}\mathop{}\\!\mathrm{d}\zeta=(2\pi)^{1/2}\omega_{\beta}^{2}\sigma_{\zeta}$). The plasma length was $c\tau=10\,\mathrm{m}$.
  * Note [2] $0.06\,\mathrm{GeV}$ electron bunch with a total charge of $1\,\mathrm{nC}$, an rms width $\sigma_{r}=0.065\,\mathrm{mm}$, a normalised emittance of $6\,\mathrm{mm}\,\mathrm{mrad}$, and an rms length $\sigma_{\zeta}/c=5\,\mathrm{ps}$. The plasma length was $c\tau=0.02\,\mathrm{m}$.
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York, 1964).
  * Thode [1976] L. E. Thode, The Physics of Fluids 19, 305 (1976), https://pubs.aip.org/aip/pfl/article-pdf/19/2/305/12638730/305_1_online.pdf .
  * Califano  _et al._ [1998] F. Califano, R. Prandi, F. Pegoraro, and S. V. Bulanov, Phys. Rev. E 58, 7837 (1998).
  * Chang  _et al._ [2016] P. Chang, A. E. Broderick, C. Pfrommer, E. Puchwein, A. Lamberts, M. Shalaby, and G. Vasil, The Astrophysical Journal 833, 118 (2016).
  * Bret  _et al._ [2008] A. Bret, L. Gremillet, D. Bénisti, and E. Lefebvre, Phys. Rev. Lett. 100, 205008 (2008).