2.1. Historical summary and context

The study of eigenfunctions of collisionless and collisional plasma kinetic equations has a long history. Vlasov was the first to suggest a method to estimate the plasma oscillation frequencies by prescribing that the principal value be taken at the resonant velocity in the dispersion function. Following this, Bohm & Gross (Reference Bohm and Gross1949) side-stepped the resonant term by considering high-enough phase velocities that the distribution function was effectively zero at the pole. Famously, Landau (Reference Landau1946) formally solved the linearized initial-value problem for Vlasov–Poisson dynamics using a Laplace transformation, discovering a discrete set of solutions for the Maxwellian plasma, wherein the electric potential damped in time. This collisionless decay phenomenon is called Landau damping. For the collisionless plasma the decaying Landau modes are not eigenfunctions, meaning that such solutions cannot evolve independently. Yet Landau's analysis also found unstable modes for certain distributions $f_0$. Thus, despite both unstable and dissipative modes sharing a similar phase space structure like $(v-\zeta )^{-1}\partial _vf_{0}{\rm e}^{{\rm i}kx}$, unstable modes evolve as a single analytic function while dissipation spreads across the entire Landau spectrum.

As the eigenvalue problem remained unsolved, Van Kampen (Reference Van Kampen1955) and Case (Reference Case1959) formally solved the problem and determined the spectrum of the linearized kinetic equation to be continuous with a possibly discrete component. Discrete eigenvalues arise only for unstable modes where $\text {Im}(\omega ) > 0$. The continuous part of the spectrum consists of ballistic modes in the form $\varepsilon (k,\zeta )\delta (v - \zeta )$, while the discrete part is in the analytic form $(v-\zeta )^{-1}\partial _v f_0$. Here $\varepsilon (k,\zeta )$ is the dielectric function and $f_0(v)$ the homogeneous equilibrium. The discrete part of the Case–van Kampen spectrum is identical to Landau's unstable modes. Landau's damping modes are represented as an integral expansion over the continuous Case–van Kampen spectrum. As a complete orthogonal system the Case–van Kampen modes are a useful though under-utilized tool. An insightful application was accomplished by P.J. Morrison and colleagues in constructing a linear integral transform, termed a ‘G transform’, to reduce the linearized Vlasov equation to an advection problem by utilizing the Case–van Kampen modes as a basis (Morrison & Pfirsch Reference Morrison and Pfirsch1992; Morrison Reference Morrison2000; Heninger & Morrison Reference Heninger and Morrison2018).

As the Landau damping modes are not collisionless eigenfunctions their status has remained somewhat obscure. Light is shed on this problem by considering weak dissipation in the Vlasov equation, such as the collision operator of Lenard & Bernstein (Reference Lenard and Bernstein1958). The landmark studies of Ng, Bhattacharjee & Skiff (Reference Ng, Bhattacharjee and Skiff1999) and Short & Simon (Reference Short and Simon2002) have shown that as dissipation tends to zero the dissipative eigenfunctions converge to the Landau damping modes, with the conclusion that dissipation is a singular perturbation of the collisionless dynamics. Bratanov et al. (Reference Bratanov, Jenko, Hatch and Brunner2013) numerically confirmed this limit for discrete systems and Ng, Bhattacharjee & Skiff (Reference Ng, Bhattacharjee and Skiff2004) extended and formalized the results for weakly collisional plasma. However, the authors wish to highlight here that the Case–van Kampen modes still play a role in the dissipative picture. Namely, one can show (Bratanov et al. Reference Bratanov, Jenko, Hatch and Brunner2013; Crews Reference Crews2022) that the propagator of the linearized kinetic equation with the Lenard–Bernstein operator limits to the Case–van Kampen modes as the dissipation $\nu \to 0^+$. This fact clarifies the relationship between the continuous Case–van Kampen spectrum and Landau damping modes, as even in the dissipative picture the Landau modes are represented as an integral over the diffusive propagator, and the collisionless Case–van Kampen modes indeed play the role of non-diffusive propagators (Balescu Reference Balescu1997). To understand why Landau modes are easily identified in the initial-value problem, consider that Landau damping originates from phase mixing (Mouhot & Villani Reference Mouhot and Villani2011) and so the modes possess the peculiar property of decaying in both directions of time. Thus, they arise by propagating initial data, or otherwise must be represented as an interference of free streaming modes unmixing from $t\in (-\infty, 0)$ and remixing from $t\in (0, \infty )$.

In summary, in collisionless plasma kinetic theory unstable modes are normal modes and evolve independently, while dissipative modes either occur as a summation of non-orthogonal transient modes or must be represented in the Case–van Kampen spectrum. The Case–van Kampen spectrum has found fruitful application as the basis of constructing an integral transform theory for linearized dynamics, most recently explored in Heninger & Morrison (Reference Heninger and Morrison2018). In weakly collisional dynamics dissipative modes are also eigenfunctions and limit to the collisionless Landau mode spectrum as dissipation tends to zero. As a consequence, it is possible to excite a single-mode instability by initializing its eigenfunction, but one may not excite a lone Landau-damping plane wave under collisionless dynamics. The following discussion illustrates this point.

2.2. Phase space linear response and the dielectric function

The electric susceptibility $\chi$ is the linear response function which relates in the spatiotemporal frequency domain $(\boldsymbol {x}, t) \to (\boldsymbol {k},\omega )$ the polarization $\boldsymbol {\tilde {P}}(\boldsymbol {k},\omega )$ and electric field $\boldsymbol {\tilde {E}}(\boldsymbol {k}, \omega )$ by the constitutive relation $\boldsymbol {\tilde {P}}=\varepsilon _0\chi \boldsymbol {\tilde {E}}$. In the scheme of electrostatic theory one aims to determine the susceptibility and consequently the dielectric permittivity $\varepsilon (\boldsymbol {k},\omega ) = \varepsilon _0(1 + \chi )$ of a plasma with equilibrium distribution $f_0(\boldsymbol {v})$. The susceptibility is determined through the self-consistent particle response $f_1(\boldsymbol {k},\boldsymbol {v},\omega )$ such that $f = f_0 + f_1$, leading to the sought-after modal structures in the charge density. The distribution $f_1$ encodes a linear response of the charges to the electric potential $\varphi$ according to the relation $f_1 = h(\boldsymbol {k}, \boldsymbol {v}, \omega )\varphi$ for some response function $h$. For this reason, we mean by ‘phase space linear response function’ the self-consistent particle response $f_1$ to the potential $\varphi$. The permittivity is often referred to as simply the dielectric function because the permittivity tensor reduces to a scalar for isotropic equilibria.

The following review of the initial-value problem is done in detail for the simplest case in order to build intuition for the electromagnetic and magnetized problems where propagation of the initial data is tedious.Footnote ¹ Using the results, we then consider the evolution of perturbations in a plasma with Maxwellian $f_0$, namely a so-called Maxwellian perturbation (that is, where $f_1(x,v) = A\sin (kx)f_0(v)$) and a special perturbation describing a propagating, damping plane wave. When this plane wave is unstable the special perturbation grows as an eigenmode, i.e. a constant phase space structure with time-dependent amplitude, and when damped its structure evolves in phase space.

Rather than the usual Laplace transform notation, we use a more standard notation for the one-sided temporal Fourier transform pair as

(2.1)\begin{gather} f(\omega) = \int_0^{\infty} f(t){\rm e}^{{\rm i}\omega t}\,{\rm d}t, \end{gather}

(2.2)\begin{gather}f(t) = \frac{1}{2{\rm \pi}}\int_{-\infty + {\rm i}s}^{\infty + {\rm i}s}f(\omega){\rm e}^{-{\rm i}\omega t}\,{\rm d}\omega, \end{gather}

where the factor $s$ keeps the contour above all poles of the integrand in order that (2.1) converges. The contour in (2.2) is closed at infinity around the lower half-plane, and the inverse transform is then given by the residue theorem as

(2.3)\begin{equation} f(t) ={-}{\rm i}\sum\text{Res}(f(\omega)) \end{equation}

with the sum over all poles of the response function $f(\omega )$. The spatial Fourier transform $\boldsymbol {x}\to \boldsymbol {k}$ is defined as usual over all of space by $f(\boldsymbol {k}) = \int \,{\rm d}\boldsymbol {x}{\rm e}^{-{\rm i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}}f(\boldsymbol {x})$. We consider for simplicity the response of a single species of particle charge $q$ and mass $m$ amidst a uniform neutralizing Maxwellian background, so we do not write a species subscript.

The Vlasov–Poisson system linearized by $f(\boldsymbol {x},\boldsymbol {v},t) = f_0 + f_1$ with $f_1\ll f_0$ is given by

(2.4)\begin{gather} \frac{\partial f_{1}}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{1} - \frac{q}{m}(\boldsymbol{\nabla}_x\varphi)\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_{0} = 0, \end{gather}

(2.5)\begin{gather}\nabla^2\varphi ={-}\frac{qn_0}{\varepsilon_0}\int f_{1} \,{\rm d}\boldsymbol{v}. \end{gather}

Fourier transforming for all of space and for time $t\in (0,\infty )$ as described, the phase space linear response is obtained as

(2.6)\begin{equation} f_1(\boldsymbol{k}, \boldsymbol{v}, \omega) = {\rm i}\frac{g(k,v)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} - \frac{q}{m}\varphi\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_{0}}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}, \end{equation}

where $g(k,v)=f_1(t=0,k,v)$ is the initial condition.

The following result is as described in Landau (Reference Landau1946) with some change in notation. Define the Cauchy integral $g^*(k,\zeta ) \equiv \int _{\mathcal {C}}({g(k,v)}/({\zeta - v}))\,{\rm d}v$, with $\mathcal {C}$ Landau's contour (that is, analytically continued into the lower-half $\zeta$-plane). Choose the coordinates such that one axis aligns with the wavevector $\boldsymbol {k}$, so that when computing the zeroth velocity moment of $f_1$, two of the velocity coordinates integrate out. The potential is found to be

(2.7)\begin{equation} \varphi(k, \zeta) = {\rm i}\frac{\sigma}{k^3}\frac{g^*(k,\zeta)}{\varepsilon(k,\zeta)}, \end{equation}

where $\sigma = q/|q|$ and $\varepsilon (k, \zeta ) = 1 + k^{-2}\int _{\mathcal {C}}({1}/({\zeta -v}))({\partial f_0}/{\partial v})\,{\rm d}v$ is the electrostatic dielectric function. Here the wavenumber $k$ is normalized to the Debye length. If each root $\zeta _j$ of $\varepsilon (k,\zeta )=0$ is simple and $g^*(k,\zeta )$ has no poles then inverse transforming in time gives

(2.8)\begin{equation} \varphi(k, t) = \frac{\sigma}{k^3}\sum_j\frac{g^*(k,\zeta_j)}{\varepsilon_\omega(k,\zeta_j)}{\rm e}^{-{\rm i}\omega_jt} \end{equation}

with $\varepsilon _\omega \equiv {\partial \varepsilon }/{\partial \omega }$. Figure 1 shows the typical locations of the infinite set of zeros $\varepsilon (k,\omega )=0$ in the lower-half frequency plane making up the Landau damping mode spectrum. The amplitude of each mode $\zeta _j$ is given by the Cauchy integral of the initial perturbation weighted by that mode's factor $\varepsilon _\omega ^{-1}$. For example, a typical perturbation used in numerical simulations is Maxwellian in velocity space, namely $g(v) = (\sqrt {{\rm \pi} }v_t)^{-1}{\rm e}^{-v^2/v_t^2}$. In this case our Cauchy integral is $g^*(\zeta )= v_t^{-1}Z(\zeta /v_t)$ with $Z(\zeta )$ the plasma dispersion function defined by Fried & Conte (Reference Fried and Conte1961). Normalizing phase velocity to $v_t$, the self-consistent electric potential response for the Maxwellian perturbation is

(2.9)\begin{equation} \varphi(k,t) ={-}\sigma\sum_{j}\frac{Z(\zeta_j)}{Z''(\zeta_j)}{\rm e}^{-{\rm i}\omega_j t} \end{equation}

as $\varepsilon (\zeta )=1-k^{-2}Z'(\zeta )$ and $\varepsilon '(\zeta )= -k^{-2}Z''(\zeta )$.

Figure 1. Locations of solutions to $\varepsilon (k,\omega )=0$ showing contours of $\text {Re}(\varepsilon )=0$ in red and $\text {Im}(\varepsilon )=0$ in green for a Maxwellian $f_0(v)$ at short wavelength, $k\lambda _D=0.5$. Solutions to $\varepsilon (k,\omega )=0$ occur at the intersection of the real and imaginary zero-contours. These damped modes represent the transient Landau spectrum which accompanies any simulation of kinetic instability not initialized with an eigenfunction.

In the theory of continuous dielectrics (Landau & Lifshitz Reference Landau and Lifshitz1946; Nicholson Reference Nicholson1983), wave energy density consists of electric field energy density multiplied by the so-called Brillouin factor $\partial _\omega (\omega \varepsilon _r)$. The denominator of (2.8) evokes the Brillouin factor because for each root $\omega _n=k\zeta _n$, we have $\partial _\omega (\omega \varepsilon ) = \omega _n\varepsilon _\omega (\omega _n)$. This energy factor $\varepsilon _\omega (k, \zeta _j)$ is monotonically increasing towards the higher Landau modes. For this reason we speculate that in Landau damping the lowest-energy state is also the least-damped mode, and that the lower-energy states are also of greater amplitude in general perturbations.

2.3. Electrostatic eigenfunctions and transient responses

Equation (2.8) is the response when the initial data $g^*(k,\zeta )$ is an entire function of $\zeta$. However, the phase space linear response (2.6) itself has a simple pole. Since the initial data supplied to the kinetic equation is supposed to simulate the self-consistent response of the plasma to some perturbation, it follows that appropriate initial data may also have a pole. Consider the initial condition to be the linear response,

(2.10)\begin{equation} g(k, v) ={-}\frac{q}{m}\frac{1}{k^2}\frac{1}{\zeta_n-v}\frac{\partial f_0}{\partial v}, \end{equation}

where $\zeta _n$ is a root of the dielectric function, $\varepsilon (k, \zeta _n) = 0$. The Cauchy transform of (2.10) is the dielectric function with an added residue for $\text {Im}(\zeta )<0$,

(2.11)\begin{equation} g^*(k,\zeta) = \frac{\sigma}{\zeta-\zeta_n}\left(\varepsilon(k,\zeta) + \begin{cases} 0, & \text{Im}(\zeta_n) > 0,\\ \left.\dfrac{2{\rm \pi} {\rm i}}{k^2}\dfrac{\partial f_0}{\partial v}\right|_{v=\zeta_n}, & \text{Im}(\zeta_n) < 0. \end{cases}\right) \end{equation}

Consider first $\text {Im}(\zeta )> 0$. Combining (2.7) and (2.11) and performing an inverse Fourier transform with (2.2) gives the potential $\varphi =k^{-2}{\rm e}^{{\rm i}(kx-\omega _nt)}.$ The phase space structure is given by (2.6), and again observing the form of (2.10) and combining with the expression for the potential gives

(2.12)\begin{equation} f_1(x,v,t) = g(k,v){\rm e}^{{\rm i}(kx - \omega_n t)}. \end{equation}

Therefore, (2.12) is a linear eigenfunction of the Vlasov equation. The phase space fluctuation grows in time and there is no phase mixing. On the other hand, in the case of $\text {Im}(\zeta )<0$ the spectral potential contains a residue,

(2.13)\begin{equation} \varphi(k,\zeta) = \frac{{\rm i}}{k^3}\frac{1}{\zeta-\zeta_n} - \frac{2{\rm \pi}}{k^3}\frac{1}{(\zeta-\zeta_n)\varepsilon(k, \zeta)}\left.\frac{1}{k^2}\frac{\partial f_0}{\partial v}\right|_{v=\zeta_n}. \end{equation}

Equation (2.13) has a double pole at $\zeta =\zeta _n$ in the second term, and a simple pole at all other roots $\varepsilon (k, \zeta _j)=0$ with $j\neq n$. Inverting the solution obtains the expression

(2.14)\begin{align} \varphi(k,t) & = \left\{1 + \left.\frac{{\rm i}{\rm \pi}}{k^2}\frac{\varepsilon_{\omega\omega}(k,\zeta_n)}{\varepsilon_\omega^2(k,\zeta_n)} \frac{\partial f_0}{\partial v}\right|_{v=\zeta_n}\right\}\frac{{\rm e}^{-{\rm i}\omega_n t}}{k^2}\nonumber\\ & \quad + \sum_{j\neq n}\left(\left.\frac{1}{(\omega_n-\omega_j)\varepsilon_\omega(k,\zeta_j)} \frac{2{\rm \pi}}{k^2}\frac{\partial f_0}{\partial v}\right|_{v=\zeta_n}\right)\frac{{\rm e}^{-{\rm i}\omega_jt}}{k^2}. \end{align}

Although the kinetic mode with frequency $\omega _n$ is preferentially excited by this perturbation, it is evident that all Landau modes are also necessarily involved. The form of (2.14) demonstrates phase mixing and decay at the Landau damping frequencies. However, at long wavelength the mode propagates decoupled from the others to $O(k^{-2})$.

So far the results are independent of the specific form of $f_0(v)$ other than its spatial uniformity. Now, to illustrate the partition of energy into the damped modes, the relative response amplitudes to (i) a Maxwellian perturbation and (ii) one with a single pole, namely (2.10), are computed for a Maxwellian equilibrium distribution $f_0(v)$. In case (i) the Maxwellian perturbation evolves with the potential

(2.15)\begin{equation} \varphi(k, t) ={-}2\sum_j\frac{Z(\tilde{\zeta}_j)}{Z''(\tilde{\zeta}_j)}{\rm e}^{-{\rm i}\omega_j t}, \end{equation}

where $\tilde {\zeta }_j = \zeta _j/v_t$, while case (ii) evolves with the potential

(2.16)\begin{equation} \varphi(k, t) = \frac{1}{k^2}\left\{\left(1 - 2{\rm \pi} {\rm i} \frac{Z'''(\tilde{\zeta}_n)}{(Z''(\tilde{\zeta}_n))^2} \frac{\partial f_0}{\partial v}\right){\rm e}^{-{\rm i}\omega_n t} - 4\sqrt{2}{\rm \pi}\frac{\partial f_0}{\partial v}\sum_{j\neq n}\frac{{\rm e}^{-{\rm i}\omega_jt}}{(\zeta_n-\zeta_j)Z''(\tilde{\zeta}_j)}\right\}. \end{equation}

Figure 2 compares these relative potential amplitudes for the two cases (i) and (ii), where figure 2(b) represents a rightward-propagating Langmuir wave. While mostly the primary plasma oscillation mode is excited, the higher modes make a substantial contribution.

Figure 2. Relative potential amplitudes (normalized to the highest magnitude) of the Landau damping spectrum for the first ${\pm }10$ frequencies $\omega _j$ given a Maxwellian equilibrium distribution $f_0(v)$ at $k\lambda _D=0.5$ for (a) the Maxwellian perturbation with $g^*(\zeta )=Z(\zeta )$, and (b) initial data as (2.10) with a single complex pole. The modes are ordered by the real part of their frequency. Higher modes are rapidly damped, as seen in figure 1 which is computed for the same equilibrium Maxwellian distribution. Panel (b) shows that while (2.10) can excite a propagating, damping plane wave it is necessarily accompanied by transient oscillations in its phase space structure with initial amplitudes up to 20 % of the primary oscillation.

In summary, Landau damping modes are not eigenfunctions of the Vlasov equation. If they are initialized and time is run either forward or backward they damp through phase mixing in either direction of time. However, their phase space structure is essentially the same as that of the unstable eigenfunctions, namely the plasma part of the plasma-field configuration occurring in a plasma wave. On the other hand, unstable modes are true eigenfunctions whose phase space structures do not change in time.

2.4. Visualizing the electrostatic phase space eigenfunctions

Let us visually explore the phase space structure just discussed mathematically. Given a solution $\zeta _n$ to $\varepsilon (\zeta _n, k_n)=0$ for particular $k_n$ (for instability, $\text {Im}(k_n\zeta _n) > 0$) the perturbed distribution is given by (2.12) in complex-conjugate pair. Examining the real part gives

(2.17)\begin{gather} f_1(x, v, t_0) = \alpha\text{Re}(\psi)\frac{\partial f_0}{\partial v}, \end{gather}

(2.18)\begin{gather}\psi(x, v) \equiv \frac{{\rm e}^{{\rm i}k_n (x - \zeta_n t_0)}}{\zeta_n - v} \end{gather}

where $\alpha$ is the perturbation amplitude. Provided that no solution to $\varepsilon (\zeta, k)=0$ has $\text {Im}(\zeta ) = 0$, the denominator of $\psi$ does not vanish and the function is well-defined. The complex function $\psi$ is defined for convenience to account for all phase information in the perturbation. The real part is chosen arbitrarily as the linear modes come in conjugate pairs. Note that one can think of the mode as an instantaneous Cauchy transform of the distribution gradient.

For example, consider the unstable two-stream mode of two drifting Maxwellians of drift velocities $v_d = \pm 4 v_t$ with wavenumber $k_0\lambda _D = 0.1$, giving a growth rate $\omega _i \approx 0.28\omega _p$. Figure 3 visualizes the phase space of the corresponding mode given by (2.17)–(2.18). The perturbation is a two-dimensional oscillatory structure in phase space, yet the zeroth moment is a pure sine wave resulting in an initial electric potential $\varphi (x) = \varphi _0\sin (x)$. The non-separability of the perturbation is evident in figure 3.

Figure 3. Depicted is the phase space eigenfunction of an unstable two-stream mode on two drifting Maxwellians of drift velocities $v_d = \pm 4v_t$, which appears as coupled plasma waves on each drifting distribution in the specific sense that the phase space structure on each beam is visually and mathematically similar to a Landau-damped mode of a thermal Maxwellian (i.e. the phase space structure of (2.10)). However, these plasma wave structures only occur as a normal mode, or eigenfunction, when unstable. The mode has structure in both $x$ and $v$, but its zeroth velocity moment is a pure sine wave and its first moment is zero.

2.5. Kinetic eigenfunctions applied to nonlinear initial-value problems

The physically correct initial data is usually discussed in the context of the thermal fluctuation spectrum (Ichimaru Reference Ichimaru1992). Linear eigenfunctions grow from spontaneous thermal fluctuations until nonlinear saturation at some significant fraction of the thermal energy (Yoon Reference Yoon2007; Crews & Shumlak Reference Crews and Shumlak2022). The Vlasov model does not resolve thermal fluctuations, but this is acceptable for typical plasmas as the magnitude of such fluctuations is much less than the thermal energy.

The eigenfunction part of a general perturbation amplifies its energy while the non-eigenfunction part decays with the same time scale $\omega _p^{-1}$. Clearly, with sufficiently small initial amplitude the non-eigenfunction part of a general perturbation does not participate in nonlinear saturation, so that sufficiently low-amplitude general perturbations are physically correct. Yet in the same way, an eigenfunction perturbation of initially large amplitude compared with the thermal fluctuation level is also physically correct. For this reason, eigenfunction perturbations yield a physically meaningful computational cost-savings when initialized with amplitudes just below nonlinear levels, while high-amplitude general perturbations introduce nonlinear Landau damping. Initialization at high amplitude translates to considerable computational savings for high-dimensional, computationally intensive continuum kinetic problems.

We mention an application of eigenfunction perturbations to small-amplitude perturbation problems. Sometimes linear instability growth rates are measured for verification of model implementation (Ho, Datta & Shumlak Reference Ho, Datta and Shumlak2018; Einkemmer Reference Einkemmer2019). Small-amplitude eigenfunction perturbations allow linear instability growth rates to be deduced from data with basically arbitrary precision because there is a complete absence of Landau damping.

2.5.1. Phase space eigenfunctions applied to the two-stream instability problem

Here we refer to perturbations as ‘separable’ when they factor as $g(x,v) = h(x)f_0(v)$ with $h(x)$ representing the desired density perturbation. Given the preceding discussion, it is illustrative to compare energy traces of fully nonlinear Vlasov–Poisson simulations initialized both with general separable perturbations and eigenfunction perturbations. Consider, for example, the two-stream unstable distribution

(2.19)\begin{equation} f_0(v) = \frac{1}{2\sqrt{2{\rm \pi}}v_t}\left(\exp\left\{-\frac{(v-v_d)^2}{2v_t^2}\right\} + \exp\left\{-\frac{(v+v_d)^2}{2v_t^2}\right\}\right). \end{equation}

Initialization with a separable Maxwellian perturbation, namely $g(x,v) = \alpha f_0(v){\rm e}^{{\rm i}kx}$ with $\alpha$ a scalar amplitude, leads to the linear solution

(2.20)\begin{equation} \varphi(k, t) ={-}2\alpha \sum_j \frac{Z(\zeta_{+,j}) + Z(\zeta_{-,j})}{Z''(\zeta_{+,j}) + Z''(\zeta_{-,j})}{\rm e}^{-{\rm i}\omega_j t}, \end{equation}

where $\zeta _{\pm,j} \equiv (\zeta _j \pm v_d)/\sqrt {2}v_t$ are the beam-shifted phase velocities. Table 1 lists the greatest amplitudes of (2.20) for drifts $v_d=5v_t$ at wavenumber $k\lambda _D=0.126$, and shows that the unstable mode is only the fifth-largest amplitude. Langmuir waves influence dynamics by masking the growing mode or through nonlinear Landau damping.

Table 1. Relative amplitudes of electric potential for the first six linear mode pairs in a separable perturbation of two counterstreaming thermal plasmas of drift velocities $v_d= 5v_t$ and at wavenumber $k\lambda _D=0.126$. In this case, for each non-zero frequency $\omega$ there is a pair of solutions. Due to these mode pairs the unstable mode (here of growth rate $\gamma /\omega _p=0.335$) has only the fifth-largest relative amplitude of electric potential in the initial condition (counting the modes twice, as they come in pairs). The significance is that the unstable mode contains only a small fraction of the initialized energy, as all modes occur at the same wavelength. The amplitude and growth rate matching is here only a coincidence.

Figure 4 compares the energy traces of nonlinear simulations of the two-stream instability initialized by both separable and eigenfunction perturbations. When initialized at small amplitude the type of perturbation does not make a difference to saturation dynamics. On the other hand, a perturbation of large amplitude reaches saturation much faster. The large amplitude Maxwellian perturbation in figure 4(a) introduces nonlinearities and changes the energy trace from the desired evolution. Observe that the simulations initialized with the eigenfunction perturbations undergo a pure growth.

Figure 4. Comparison of energy traces for nonlinear simulations of a single wavelength two-stream instability of drift velocities $v_d = 5v_t$ at wavenumber $k\lambda _D=0.126$ showing (a) separable perturbation and (b) eigenfunction perturbation. Two subcases are shown. The black trace represents a small-amplitude perturbation and the blue trace a large-amplitude one. The large-amplitude Maxwellian perturbation introduces nonlinearity, polluting the solution, while the large-amplitude eigenmode saturates equivalently as if seeded from small amplitude. In any case solutions initialized by eigenfunction perturbation undergo pure growth in the linear phase.

2.6. Multidimensional dispersion function for the two-stream instability

Electrostatic turbulence at the Debye length scale generated by streaming instability of electron beams is a ubiquitous plasma phenomenon (Rudakov & Tsytovich Reference Rudakov and Tsytovich1978; Che Reference Che2016), and is inherently three-dimensional. This section considers the two-stream instability in the computationally tractable two-dimensional configuration space as an example of the methodology used to compute electrostatic phase space eigenfunctions in multiple dimensions. Recall that the electrostatic dielectric function is determined by

(2.21)\begin{equation} \varepsilon(\omega, \boldsymbol{k}) = 1 + \frac{\omega_p^2}{k^2}\int \frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_0}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\,{\rm d}\boldsymbol{v} = 0. \end{equation}

Now consider a thermal two-stream distribution with equal temperatures on each beam,

(2.22)\begin{equation} f_0(u,v,w) = \frac{{\rm e}^{-(v^2+w^2)/(2u_t^2)}}{2(2{\rm \pi})^{3/2}u_t^3} \left(\exp\left(-\frac{(u-u_d)^2}{2u_t^2}\right) + \exp\left(-\frac{(u+u_d)^2}{2u_t^2}\right)\right), \end{equation}

which differs from (2.19) in retaining three components of velocity. Having assumed an isotropic thermal velocity $u_t$ greatly simplifies analysis; otherwise complications arise due to the elliptical level-sets of $f_0$. Let the wavevector lie in the $(x,y)$-plane and consider (2.21). The $w$-component integrates out immediately, while the $(x,y)$-directed velocities must be rotated into the frame of the wavevector. Rotating through an angle $\varphi$ to coordinates $(v_\parallel, v_\perp )$, the distribution function is $f_0(v_\parallel, v_\perp ) = (f_+ + f_-)/2$ with

(2.23)\begin{equation} f_{{\pm}} \equiv \frac{1}{2{\rm \pi} u_t^2}\exp\left(-\frac{(v_\parallel{\mp} \cos(\varphi)u_d)^2}{2u_t^2}\right) \exp\left(-\frac{(v_\perp{\pm} \sin(\varphi)u_d)^2}{2u_t^2}\right). \end{equation}

Evaluating the integral for each drifting component $f_\pm$ as $-Z'(\zeta _\pm )/2u_t^2$ as in Skoutnev et al. (Reference Skoutnev, Hakim, Juno and TenBarge2019) with drift velocity-shifted phase velocity $\zeta _\pm \equiv (\zeta \mp \cos \varphi u_d)/\sqrt {2}u_t$ gives

(2.24)\begin{equation} \varepsilon(k, \zeta, \varphi) = 1 - \frac{1}{(k\lambda_D)^2}\frac{Z'(\zeta_+) + Z'(\zeta_-)}{2} = 0. \end{equation}

The wavevector having transverse components to the drift axis decreases the effective drift by the cosine of $\varphi$, leading to maximum growth rate of longitudinal waves parallel to the streaming velocity. However, the growth of these transverse-axis components seeds a multidimensional turbulence, depending on the configuration space dimensionality.

2.7. Nonlinear simulation of two-stream instability in two spatial dimensions

Although the fastest growing mode of the two-stream problem has a beam-axis-aligned wavevector, eigenmodes with a long-wavelength transverse part $k_\perp \ll k_\parallel$ grow comparably to the fastest mode, as illustrated in figure 5 by the dielectric function $\varepsilon (k_\parallel, k_\perp )$ computed using (2.24). In this way, streaming in unmagnetized plasma generically produces multidimensional Langmuir turbulence. In practice this is a three-dimensional phenomenon, but for computational tractability a streaming simulation is presented here with two space and two velocity dimensions (2D2V). We speculate that 2D2V nonlinear phase mixing is quite similar to three space and three velocity dimensions (3D3V) because the ‘total dimensionality’ of the phase space turbulence is greater than two. On the other hand, nonlinear spatial dynamics of the saturated state are likely quite different in 2D2V versus full 3D3V. A dedicated study of this issue with sufficient computer capability would be welcome.

Figure 5. Growth rates of the electron two-stream instability as $\text {Im}(\omega )/\omega _{pe}$ for two counter-streaming Maxwellians with drift speeds relative to the thermal speed of $u_d/u_t = \pm 3$, in terms of parallel and perpendicular wavenumbers relative to the beam axis vector $\hat {u}_d$. The red contour indicates the line of marginal stability where $\text {Im}(\omega ) = 0$. While the fastest growing mode is aligned with the beam axis (occurring at $(k\lambda _D)\boldsymbol {\cdot }\hat {u}_d\approx 0.2$), modes of comparable growth rates occur with transverse wavelengths roughly five-to-ten times the unstable wavelength on-axis.

2.7.1. Initialization of the two-dimensional two-stream simulations

Our numerical method is summarized in Appendix A. The domain used is periodic and set by fundamental wavenumbers $k_x\lambda _D=0.05$ and $k_y\lambda _D=0.02$. The $x$-axis is divided into 40 evenly spaced collocation nodes and the $y$-axis into 50 nodes. Velocity space is truncated at $v_{\max } = \pm 11.5 v_t$, and each axis divided into 14 finite elements each of a seventh-order Legendre–Gauss–Lobatto polynomial basis. A non-uniform velocity grid is used; 10 elements are linearly clustered between $v\in (-7.5, 7.5)v_t$ and two elements into $v\in \pm (7.5, 11.5)v_t$. The drift velocity in (2.22) is set to $u_d=3v_t$. Finally, a spatial hyperviscosity $\nu \nabla _x^4 f$ with $\nu =10$ is added to the kinetic equation to mitigate spectral blocking with this low spatial resolution, as in Crews & Shumlak (Reference Crews and Shumlak2022). Many modes of comparable growth rates are initialized using the eigenfunction perturbations

(2.25)\begin{equation} f_1(x, y, u, v) = \frac{\alpha}{k}\text{Re}\left( \frac{k_x \partial_u f + k_y\partial_v f}{\omega(k_x, k_y) - k_x u - k_y v} \exp({\rm i}(k_x x + k_y y + \theta))\right) \end{equation}

with $\alpha$ an amplitude scalar, $\omega (k_x, k_y)$ solution to the dispersion relation for the mode $(k_x, k_y)$, and of random phases $\theta$. A total of 33 modes are excited, each with the amplitude $\alpha =0.01$; for each harmonic of the $x$-fundamental $k_{n,x} = n k_{1,x}$ with $n=2$, $3$ and $4$, 11 harmonics of the $y$-fundamental are excited with $k_{m,y} = m k_{1,y}$ with $m = 0$, ${\pm }1$, ${\pm }2$, ${\pm }3$, ${\pm }4$ and ${\pm }5$. There is no symmetry in the $y$-direction as different phases $\theta$ are used for the modes ${\pm }m$. Mode $n=1$ has small growth rate and is not initialized.

2.7.2. Two-dimensional nonlinear simulations of the streaming instability

The simulation is run to a stop time of $t\omega _p=30$. Figure 6 shows the simulation's electric field energy trace. Hyperviscosity with this spatial resolution leads to a domain energy loss of $O(10^{-3})$ while electric energy saturates at $O(10^{-1})$. Due to the use of eigenfunction perturbations there are no oscillations in the electric field energy trace. Therefore, the simulation was initialized with a perturbation energy just two orders below saturation. In this case, this saves approximately $10\omega _p^{-1}$ of simulation time compared with, for example, a perturbation of initial energy $10^{-6}$.

Figure 6. Electric energy trace of the 2D2V two-stream instability with thirty-three excited domain modes. Initialization with eigenmodes enables a perturbation just two decades below saturation energy. The aggressive spatial hyperviscosity of $\nu =10$ leads to an energy loss of $O(10^{-3})$ by the stop time of this simulation. This artificial dissipation would not be necessary with greater spatial resolution, but the loss is well beneath the electric energy at the stop time.

Figure 7 plots evolution of electric potential. Beam-axis wave energy initially predominates as electron holes form with wavenumber $k_x\gg k_y$, yet the transverse perturbations with comparable growth rates to the maximal beam-axis part lead to two-dimensional structures. Following nonlinear saturation, wave energy significantly increases in the transverse direction as the holes tilt, consolidate and isotropize. Figure 8 visualizes the domain-averaged (coarse-grained) distribution and demonstrates that isotropization associated with beam-driven electrostatic turbulence distributes coarse-grained energy into the beam-transverse directions. That is, beam-transverse temperature increases significantly in the resulting marginally stable double-humped distribution (Penrose Reference Penrose1960), as heat transfers from coarse-grained wave energy to the coarse-grained distribution (Nicholson Reference Nicholson1983). These results support the notion that the continuum of electron hole solutions is key to strongly driven plasma transport (Schamel Reference Schamel2023).

Figure 7. Electric potential $\varphi (x,y)$ of the 2D2V electron two-stream instability for three times: (a) $t\omega _p=8$ (linear phase); (b) $t\omega _p=18$ (nonlinear phase); (c) $t\omega _p=27$ (isotropizing). The evolution begins with the formation of one-dimensional electron holes, or phase space vortex lines, transverse to the streaming axis. These vortex lines then break up after saturation into two-dimensional hole structures with more complex orbits which maintain connection to one another by the Vlasov-dynamical conservation of phase space circulation. The lines of potential in this simulation may be understood as a projection of these phase space vortex tubes.

Figure 8. Domain-averaged (coarse-grained) distribution $\langle \,f\rangle _{(x,y)}(v_x,v_y)$ showing relaxation to a Penrose-stable distribution on averaged scales due to multidimensional Langmuir turbulence, for two times: (a) $t\omega _p=0$, the unstable initial condition; (b) $t\omega _p=30$, the stop time of the simulation. At the end of the simulation the average distribution is double-humped, having heated significantly in the direction transverse to the beam axis compared with the initial state.

2.8. Quasilinear kinetic simulation with phase space eigenfunctions

Quasilinear theory is the name given to the simplest closure in the hierarchy of equations resulting from separating the variables of a turbulent system into fluctuating and mean components (Vedenov Reference Vedenov1963). The scheme of the theory is as follows: a suitable method of averaging is defined, typically temporal, spatial or ensemble averaging; the dynamical equation is averaged and the mean subtracted from the original equation to obtain the mean and fluctuating components of the system; lastly a closure hypothesis is made by neglecting the ‘second fluctuation’. Under this procedure, known as ‘classic QLT’, the equation for the fluctuation becomes quasilinear and can be solved by spectral methods (Crews & Shumlak Reference Crews and Shumlak2022). Substitution into the equation for the mean gives a diffusion equation in velocity space. Beyond classic QLT, a formulation more applicable to inhomogeneous plasmas is presented in Dodin (Reference Dodin2022).

Diffusion equations are numerically stiff when diffusivity is large. A drawback of QLT posed as a diffusion problem is that the diffusivity is asymptotically singular in the relaxed state of an unstable system. The singularity arises as $\text {Im}(\zeta )\to 0$ around the purely real frequencies of the relaxed state (Crews & Shumlak Reference Crews and Shumlak2022). This singularity is side-stepped by solving the equations of QLT as an initial-value problem for a system of first-order equations, resolving the linear kinetic eigenfunctions, linear Landau damping and the asymptotic ($t\to \infty$) saturation of the distribution function. There is no dimensionality reduction like in the diffusion theory, but significant advantage remains over the fully nonlinear theory because the turbulent nonlinear cascade does not form and only the unstable scales need be resolved. Further, as a first-order system there is no need to solve the dielectric function in a quasilinear simulation. To prevent spurious Landau damping it is wise to utilize kinetic eigenfunction perturbations as in § 2.5.

To demonstrate we consider the kinetic equation for electrons in a neutralizing background. Splitting the distribution function $f = \langle \,f\rangle _L + \delta f$ where $\langle \boldsymbol {\cdot }\rangle _L$ is a spatial average, the quasilinear system in normalized units is (Crews & Shumlak Reference Crews and Shumlak2022)

(2.26)\begin{gather} \frac{\partial\langle \,f\rangle_L}{\partial t} = \frac{\partial}{\partial v}\langle E\delta f\rangle_L, \end{gather}

(2.27)\begin{gather}\frac{{\rm d}(\delta f)}{{\rm d}t} = \frac{\partial}{\partial v}E\langle \,f\rangle_L, \end{gather}

with ${{\rm d}}/{{\rm d}t}=\partial _t + v\partial _x$ the change along a zero-order trajectory. Consider a finite periodic domain $x\in (0,L)$ and the expansion of the distribution in Fourier series,

(2.28)\begin{equation} f(x,v,t) = f_0 + \sum_{n=1}^{\infty}\left(f_n {\rm e}^{{\rm i}k_n x} + f_n^*{\rm e}^{-{\rm i}k_{n}x}\right). \end{equation}

The $n=0$ component of (2.28) is the average distribution, $f_0=\langle f\rangle _L$, and the remaining Fourier coefficients make up the Fourier spectrum of the fluctuation with $f_n^* = -f_n$. Thus, dropping the symbols $\delta ({\cdot })$ and $\langle \cdot \rangle _L$, (2.26) and (2.27) are

(2.29)\begin{gather} \frac{\partial f_0}{\partial t} = \frac{\partial}{\partial v}\langle Ef\rangle_L, \end{gather}

(2.30)\begin{gather}\frac{\partial f_n}{\partial t} ={-}{\rm i}k_{n}vf_n + \frac{\partial}{\partial v}\left(E_n f_0\right),\quad n\geq 1, \end{gather}

(2.31)\begin{gather}E_n = {\rm i}k_n^{{-}1}\int_{-\infty}^{\infty} f_{n}\,{\rm d}v,\quad n\geq 1, \end{gather}

where (2.31) is obtained from Gauss's law in the Fourier basis.

2.8.1. Numerical method for the initial-value problem for the quasilinear equations

Equations (2.29) and (2.30) are discretized in velocity space and (2.31) is applied as a constraint. First, the Fourier series is truncated at a chosen mode number (Galerkin projection) to resolve the range of instability, with the corresponding spatial grid identified as the evenly spaced collocation nodes of the frequency range in a standard manner through the fast Fourier transform. The velocity axis is to be discretized by the discontinuous Galerkin (DG) method similarly to the method in Appendix A. We consider the two velocity fluxes in (2.29) and (2.30),

(2.32)\begin{gather} \mathcal{T}(v) \equiv \langle Ef\rangle_L = L^{{-}1}\int_0^L E\delta f(x,v) \,{{\rm d}\kern0.7pt x}, \end{gather}

(2.33)\begin{gather}\mathcal{M}_n(v) \equiv E_n f_0 = \left({\rm i}k_n^{{-}1}\int_{-\infty}^\infty f_n(v')\,{\rm d}v'\right)f_0(v). \end{gather}

We evaluate the flux $\mathcal {T}(v)$ by the trapezoidal rule because of its ideal trigonometric quadrature properties (Boyd Reference Boyd2001) using the inverse fast Fourier transform of the spectra $E_n$ and $f_n(v)$. On the other hand, the fluxes $\mathcal {M}_n(v)$ depend only on the local field mode $E_n$ and the mean distribution $f_0(v)$, so this quantity is simply computed by quadrature in $v$. Both fluxes $\mathcal {T}$ and $\mathcal {M}_n$ are then utilized in the DG method as a nonlinear flux. The linear translation operator ${\rm i}k_n vf_n$ is discretized by quadrature and the system is integrated in time semi-implicitly by Strang splitting, as outlined in Appendix A.

2.8.2. Simulation of the bump-on-tail instability using phase space eigenfunctions

We repeat the calculation of Crews & Shumlak (Reference Crews and Shumlak2022) for the nonlinear and quasilinear evolutions of the bump-on-tail instability, with the difference that here QLT is solved as an initial-value problem in phase space instead of as a diffusion problem for the reduced distribution, and using the numerical method described in § 2.8.1. See Crews & Shumlak (Reference Crews and Shumlak2022) for the details of initialization including the domain and the zero-order distribution. In Crews & Shumlak (Reference Crews and Shumlak2022) the perturbation is constructed using the phase space eigenfunctions given by (2.10). Figure 9 compares the energy traces of the nonlinear and quasilinear simulations where QLT is solved as an initial-value problem, while figure 10 compares the spectral phase space at saturation and demonstrates the absence of phase space cascade in QLT. In addition, there is an absence of Landau damping as the perturbation is constructed from eigenfunctions.

Figure 9. Electric energy traces are compared for (a) the quasilinear evolution, and (b) the fully nonlinear Vlasov–Poisson simulation (from the data of Crews & Shumlak (Reference Crews and Shumlak2022)) from identical initial conditions and perturbations. The quasilinear system approaches the marginally stable state asymptotically, while the fully nonlinear system saturates in a finite time due to nonlinear particle trapping. Both simulations are initialized at high though subnonlinear amplitude without subsequent Landau damping oscillations because they utilize kinetic eigenfunction perturbations, while non-eigenfunction perturbations would oscillate significantly.

Figure 10. Spectral $(k\lambda _D, v/v_t)$ phase space is compared at saturation $t=150\omega _p^{-1}$ between (a) the quasilinear system, and (b) the corresponding fully nonlinear Vlasov–Poisson simulation. The spectral phase of the quasilinear system is grossly over-resolved in $k\lambda _D$ space only for the purpose of comparison; energy remains spatially localized in QLT because the system is unable to cascade through turbulent mixing of phase space eddies, a fully nonlinear phenomenon. The phase space structures in a simulation of QLT always remain linear eigenfunctions of the form of (2.10). In this way sub-Debye length scales need not be resolved in simulation of QLT.

3.1. Review of the Harris dispersion relation

Our treatment here follows Gurnett & Bhattacharjee (Reference Gurnett and Bhattacharjee2017, § 10.2.1) with a particular focus on the phase space perturbation as a kinetic eigenfunction. Coordinates are chosen such that $\boldsymbol {B} = B\hat {z}$, $\boldsymbol {k} = k_\perp \hat {x} + k_\parallel \hat {z}$ and the velocity space is expressed in cylindrical coordinates as $\boldsymbol {v} = v_\perp \cos (\phi )\hat {v}_x + v_\perp \sin (\phi )\hat {v}_y + v_\parallel \hat {v}_z$. In these coordinates the Vlasov equation linearizes around $f_0(v_\perp, v_\parallel )$ as, with $\omega _c=qB/m$,

(3.1)\begin{equation} \frac{\partial f_1}{\partial t} + v_{{\perp}}\cos\phi\frac{\partial f_1}{\partial x} + v_{{\parallel}}\frac{\partial f_1}{\partial z} - \omega_c\frac{\partial f_1}{\partial\phi} = \frac{q}{m}\left(E_x\cos\phi\frac{\partial f_0}{\partial v_\perp} + E_z\frac{\partial f_0}{\partial v_\parallel}\right). \end{equation}

The linearized equation is then Fourier transformed $(x,z,t)\to (k_\perp,k_\parallel,\omega )$ to yield

(3.2)\begin{equation} \frac{{\rm d}f_1}{{\rm d}\phi} + {\rm i}\frac{\omega - k_\perp v_\perp \cos\phi - k_\parallel v_\parallel}{\omega_c}f_1 ={-}{\rm i}\frac{q\varphi}{m}\left(k_\perp\cos\phi \frac{\partial f_0}{\partial v_\perp} + k_\parallel\frac{\partial f_0}{\partial v_\parallel}\right). \end{equation}

Equation (3.2) is a first-order inhomogeneous equation in the cylindrical velocity-space angle $\phi$ and can be solved by the usual methods. Integrating the inhomogeneous term along the solution of the homogeneous equation yields

(3.3)\begin{equation} f_1(v_\perp, \phi, v_\parallel) ={-}\frac{q\varphi(k)}{m}\exp({\rm i}k_\perp v_y/\omega_c)\left[ \frac{1}{v_\perp}\frac{\partial f_0}{\partial v_\perp}\varUpsilon_1 + k_\parallel\frac{\partial f_0}{\partial v_\parallel}\varLambda_1\right], \end{equation}

where the terms $\varUpsilon _1$ and $\varLambda _1$ are auxiliary wavefunctions defined as the polar Fourier series

(3.4)\begin{gather} \varUpsilon_1 = \sum_{n={-}\infty}^{\infty}\frac{n\omega_c}{\omega-k_\parallel v_\parallel{-} n\omega_c} {\rm J}_n(k_\perp v_\perp{/}\omega_c){\rm e}^{{\rm i}n\phi}, \end{gather}

(3.5)\begin{gather}\varLambda_1 = \sum_{n={-}\infty}^{\infty}\frac{\omega_c}{\omega-k_\parallel v_\parallel{-} n\omega_c} {\rm J}_n(k_\perp v_\perp{/}\omega_c){\rm e}^{{\rm i}n\phi}, \end{gather}

with ${\rm J}_n(z)$ a Bessel function of the first kind. The terms of $\varUpsilon _1$ associated with perpendicular propagation decay one order slower in $n$ than $\varLambda _1$, meaning that high-order resonances are more important for perpendicular propagation. Equation (3.3) is the phase space linear response for electrostatic fluctuations in a strongly magnetized plasma.

The self-consistent spectrum consists of all pairs $(\omega, k_\perp, k_\parallel )$ such that the zeroth moment of $f_1(\omega, k_\perp, k_\parallel, v_\perp, \phi, v_\parallel )$ results in an electric potential mode of wavenumber $(k_\perp, k_\parallel )$. Integration of the phase space fluctuation gives the density fluctuation as

(3.6)\begin{equation} n_1(\omega, k_\perp, k_\parallel) ={-}{\rm i}\frac{q\varphi}{m}\int_{-\infty}^{\infty}\int_0^\infty \left[ \frac{1}{v_\perp}\frac{\partial f_0}{\partial v_\perp}\varUpsilon_2 + k_\parallel\frac{\partial f_0}{\partial v_\parallel}\varLambda_2\right]2{\rm \pi} v_\perp \,{\rm d}v_{{\perp}}\,{\rm d}v_{{\parallel}} \end{equation}

where a set of additional series, analogues of (3.4) and (3.5), are defined as

(3.7)\begin{gather} \varUpsilon_2 = \sum_{n={-}\infty}^{\infty}\frac{n\omega_c}{\omega - k_\parallel v_\parallel{-} n\omega_c}{\rm J}_n^2(k_\perp v_\perp{/}\omega_c), \end{gather}

(3.8)\begin{gather}\varLambda_2 = \sum_{n={-}\infty}^{\infty}\frac{\omega_c}{\omega - k_\parallel v_\parallel{-} n\omega_c} {\rm J}_n^2(k_\perp v_\perp{/}\omega_c). \end{gather}

Substitution of $n_1$ into Poisson's equation gives Harris's dispersion relation

(3.9)\begin{equation} \varepsilon(\omega, k_\perp, k_\parallel) \equiv 1 - \left(\frac{\omega_p}{\omega_c}\right)^2\frac{1}{(k\lambda_D)^2} \int_{-\infty}^{\infty}\left(\mathbb{V}_\perp f_\parallel{+} k_\parallel \mathbb{V}_\parallel \frac{\partial f_\parallel}{\partial v_\parallel}\right) \,{\rm d}v_{{\parallel}} = 0, \end{equation}

where the integration over perpendicular velocities is broken out into the two quantities

(3.10)\begin{gather} \mathbb{V}_{{\perp}} = \int_0^\infty\frac{1}{v_\perp}\frac{\partial f_\perp}{\partial v_\perp}\varUpsilon_2(v_\perp, v_\parallel)2{\rm \pi} v_{{\perp}}\,{\rm d}v_{{\perp}}, \end{gather}

(3.11)\begin{gather}\mathbb{V}_{{\parallel}} = \int_0^\infty f_\perp \varLambda_2(v_\perp, v_\parallel)2{\rm \pi} v_{{\perp}}\,{\rm d}v_{{\perp}}, \end{gather}

and separability of the background $f_0(v_\perp, v_\parallel ) = f_\perp (v_\perp )f_\parallel (v_\parallel )$ has been assumed.

3.2. Amplitude limitation of linearization around zero-order cyclotron orbits

Given a zero-order spatially uniform magnetic field $\boldsymbol {B} = \boldsymbol {B}_0$ and first-order electric field perturbation $\boldsymbol {E} = \boldsymbol {E}_1$, the zero- and first-order kinetic equations are

(3.12)\begin{gather} (\boldsymbol{v}\times\boldsymbol{B}_0)\boldsymbol{\cdot}\boldsymbol{\nabla}_{v}f_0 = 0, \end{gather}

(3.13)\begin{gather}\partial_{t}f_1 + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}_{x}f_1 + \frac{q}{m}(\boldsymbol{v}\times\boldsymbol{B}_0)\boldsymbol{\cdot}\boldsymbol{\nabla}_{v}f_1 + \frac{q}{m}\boldsymbol{E}_1\boldsymbol{\cdot}\boldsymbol{\nabla}_{v}f_0 = 0 \end{gather}

assuming a homogeneous zero-order distribution $f_0 = f_0(\boldsymbol {v})$. Equation (3.12) indicates gyrotropy of $f_0$. This ordering is valid when the zero-order cyclotron acceleration is much greater than the electrostatic acceleration of a typical particle. Validity translates to an amplitude restriction on electric potential and the density fluctuation. Assuming $k_\parallel =0$ and comparing terms proportional to $\boldsymbol {\nabla }_{v}f_0$ in (3.12) and (3.13) for a thermal particle gives $E_1 \ll v_{th}B_0$. Estimating the field $E_1$ of wavenumber $k$ by Gauss's law gives $E_1 = e\delta n / k\varepsilon _0$ for density fluctuation $\delta n$. Combining these estimates results in equivalent conditions on amplitude as measured by $\delta n$ or $\varphi$,

(3.14)\begin{gather} \frac{\delta n}{n_0} \ll \left(\frac{\omega_c}{\omega_p}\right)^2(kr_L), \end{gather}

(3.15)\begin{gather}\frac{e\varphi}{kT} \ll \frac{1}{kr_L}, \end{gather}

for Larmor radius $r_L = v_{th}/\omega _c$. Amplitudes which exceed these inequalities are subject to electrostatic Landau damping as the dielectric function of (3.9) is not valid. Typical cyclotron instabilities have $k_\perp r_L \approx 1$ which limits the amplitudes of the linear modes considered in this section to amplitudes $\delta n/n_0 \ll (\omega _c/\omega _p)^2$ and $e\varphi \ll kT$.

3.3. The dielectric function for ring ($\chi$-) distributions

The linear mode spectrum depends on the background distribution (the zero-order equilibrium). Plasma theory textbooks consider Maxwellian plasmas by expansion in the cyclotron harmonics (Gurnett & Bhattacharjee Reference Gurnett and Bhattacharjee2017, § 10.2.3). Of course, in ideal collisionless plasmas with plasma parameter $\varLambda \to \infty$ distributions are expected to be observed only close enough to Maxwellian such that the Penrose criterion is satisfied.

Recent analytical work on non-Maxwellian distributions focuses on the kappa distributions (Mace & Hellberg Reference Mace and Hellberg2009) to model observations in space plasma (Pierrard & Lazar Reference Pierrard and Lazar2010). Kappa ($\kappa$-) distributions, also called $q$-Gaussians, are motivated by recent advances in entropy methods (Livadiotis & McComas Reference Livadiotis and McComas2023; Zhdankin Reference Zhdankin2023). It is thought that such entropy methods may facilitate the extension of maximum entropy principles to the prediction of metastable equilibria such as non-Maxwellian velocity distributions or self-organized equilibria in magnetic confinement including tokamaks (Dyabilin & Razumova Reference Dyabilin and Razumova2015) and Z pinches (Crews et al. Reference Crews, Datta, Meier and Shumlak2024). Ewart et al. (Reference Ewart, Brown, Adkins and Schekochihin2022) is a significant recent advance with a lucid description of collisionless relaxation.

Spatially uniform strongly magnetized plasmas must have zero-order gyrotropy so that non-Maxwellian features in perpendicular velocity space are typically ring-shaped. Ring distributions commonly arise from the loss-cone mechanism of magnetic traps or planetary magnetospheres. Early identifications of velocity-space instability in ring-distributed plasmas were made by Dory et al. (Reference Dory, Guest and Harris1965). Dory's ring distribution, known in the mathematics literature as a $\chi$-distribution, is a type of maximum entropy distribution subject to two constraints on variance. The studies of Tataronis & Crawford (Reference Tataronis and Crawford1970a) and Tataronis & Crawford (Reference Tataronis and Crawford1970b) extended the theory to oblique propagation, showing maximal growth rates for near-perpendicular propagation, though analytical work was performed only with singular ring distributions. Around the turn of the millennium $q$-analogues of Dory's analytic ring distributions were introduced by Leubner & Schupfer (Reference Leubner and Schupfer2001) and extended in Pokhotelov et al. (Reference Pokhotelov, Treumann, Sagdeev, Balikhin, Onishchenko, Pavlenko and Sandberg2002), motivated by the successful use of $\kappa$-distributions as $q$-deformations of Maxwell–Boltzmann statistics. Dory's $\chi$-distribution is the $q\to 1$ limit of Leubner's $\kappa$-like ring distributions in the same way that the Maxwell–Boltzmann distribution is the $q\to 1$ limit of the $q$-Gaussian (or $\kappa$-) distributions.

For this reason, in this section we focus on Dory's ring distribution and analyse the dielectric function for such rings assuming separability of the zero-order distribution as

(3.16)\begin{equation} f_0(v_\parallel, v_\perp) = f_{{\parallel}}(v_\parallel)f_{\gamma}(v_\perp), \end{equation}

where $f_\parallel (v_\parallel )$ is a Maxwellian of thermal velocity $v_t$ and $f_\gamma (v_\perp )$ is Dory's ring function

(3.17)\begin{equation} f_{\gamma}(v) = \frac{1}{{\rm \pi}\alpha^2}\frac{1}{\varGamma(\gamma+1)}\left(\frac{v^2}{\alpha^2}\right)^{\gamma}\exp\left(-\frac{v^2}{\alpha^2}\right) \end{equation}

of thermal velocity $\alpha$ and parameter $\gamma$, where $\varGamma (z)$ is the Gamma function. The function $f_\gamma$ is a polar distribution,

(3.18)\begin{equation} \int_0^\infty f_\gamma(v)2{\rm \pi} v\,{\rm d}v = 1,\quad\text{for } \text{Re}(\gamma) >{-}1 \end{equation}

yet is bounded at zero only for $\gamma \geq 0$. Equation (3.17) is also known as a $\chi$-distribution and is ‘two-temperature’. That is, the ring distribution $f_\gamma (v_\perp )$ is the maximum entropy distribution for $v_\perp \in [0,\infty )$ subject to the two constraints $\langle v_\perp ^2\rangle = (\gamma +1)\alpha ^2$ and $\langle (v_\perp - \langle v_\perp \rangle )^2\rangle = \alpha ^2(1 + \gamma - ({\varGamma (\gamma +3/2)}/{\varGamma (\gamma +1)})^2)$ or $\langle (v_\perp - \langle v_\perp \rangle )^2\rangle = \tfrac{1}{4}\alpha ^2 + O(\gamma^{-1})$ such that the temperature in the gyrating frame is independent of $\gamma$ as $\gamma \to \infty$. Thus, the physical meaning of $\alpha$ is the thermal velocity in the gyrating frame, while $\gamma$ is the boost to thermal energy in the laboratory frame (and need not be an integer). In this sense the distribution has two temperatures.

The ring distribution satisfies the recurrence (and $f_{-1}=0$),

(3.19)\begin{equation} \frac{1}{v}\frac{{\rm d}f_\gamma}{{\rm d}v} = \frac{1}{\alpha^2}\frac{f_{\gamma-1} - f_\gamma}{2}. \end{equation}

In the case of ring distributions of the form of (3.17), the integrals in both quantities $\mathbb {V}_\perp$ and $\mathbb {V}_\parallel$ involve only $f_\gamma$ due to the recurrence Eq. (3.19), and suggests defining

(3.20)\begin{equation} \mathbb{F}_{n,\gamma}(k) \equiv \int_0^{\infty}f_{\gamma}(v){\rm J}_n^2(kv)2{\rm \pi} v\,{\rm d}v \end{equation}

which, as shown in Appendix B, is a type-$_{2}F_2$ hypergeometric function with series representation (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2015)

(3.21)\begin{equation} \mathbb{F}_{n,\gamma}(x) = \frac{1}{\varGamma(\gamma+1)} \sum_{\ell=0}^{\infty}\frac{\varGamma(2n+2\gamma+2\ell+1)} {\varGamma^2(n+\ell+1)\varGamma(2n+\ell+1)}\frac{({-}1)^\ell}{\ell!}x^{2(\ell+n)}. \end{equation}

We now proceed with integrating this power series over the parallel velocities. First, we can make a note on an alternative possibility. Rather than integrating the auxiliaries (3.7) and (3.8) in their summation form, it is possible to first close the summation with the Lerche–Newberger summation theorem (Newberger Reference Newberger1982), and to determine the perpendicular velocity integrals in (3.10) and (3.11) in closed form as hypergeometric functions, as in Appendix C, for arbitrary $k_\parallel$. However, the integration over parallel velocities must then proceed by a series expansion around the poles of these hypergeometric functions, making a power series approach inevitable. On the other hand, the power series developed in this section maintains separability of terms containing the parallel velocity.

Proceeding to the integration over parallel velocities, with the parallel distribution $f_{\parallel }(v_\parallel )$ taken as a Maxwellian distribution, the two integrals are

(3.22)\begin{gather} \int_{-\infty}^\infty \frac{n\omega_c}{\omega - k_\parallel v_\parallel{-} n\omega_c}f_\parallel \,{\rm d}v_\parallel{=}{-}\frac{n}{k_\parallel}Z(\zeta_n), \end{gather}

(3.23)\begin{gather}\int_{-\infty}^{\infty} \frac{\omega_c}{\omega - k_\parallel v_\parallel{-} n\omega_c} \frac{\partial f_\parallel}{\partial v_\parallel}\,{\rm d}v_{{\parallel}} = \frac{1}{2k_\parallel}Z'(\zeta_n), \end{gather}

where $Z(\zeta )$ is the plasma dispersion function and the cyclotron harmonic-shifted phase velocity is defined as $\zeta _n \equiv ({\omega -n\omega _c})/{\sqrt {2}k_\parallel v_t}$. Thus, the dielectric function for loss-cones is

(3.24)\begin{align} \varepsilon(\omega, k) & = 1 - \left(\frac{\omega_p}{\omega_c}\right)^2\frac{1}{(kr_L)^2} \sum_{n={-}\infty}^{\infty}\left[\vphantom{\frac{n}{k_\parallel r_L}}Z'(\zeta_n)\mathbb{F}_{n,\gamma}(k_\perp r_L) \right.\nonumber\\ & \left.\quad +\,\frac{n}{k_\parallel r_L}Z(\zeta_n) \left\{\mathbb{F}_{n,\gamma-1}(k_\perp r_L) - \mathbb{F}_{n,\gamma}(k_\perp r_L)\right\}\right]. \end{align}

Equation (3.24) reduces to the standard series for a Maxwellian ($\gamma =0$) by the identity $\mathbb {F}_{n,0}(k) = {\rm e}^{-k^2}{\rm I}_n(k^2)$ where ${\rm I}_n(k)$ is the modified Bessel function of the first kind.

3.4. Perpendicular propagation in ring distributions

In the limit of perpendicular propagation, $k_\parallel \to 0$, the series (3.24) in the cyclotron harmonics simplifies as the terms proportional to $Z(\zeta _n)$ vanish. Further, by use of the Lerche–Newberger summation theorem two alternatives to the power series for the perpendicular cyclotron wave dielectric function may be developed which incorporate the contributions from the cyclotron harmonics to all orders. The derivation of these expressions may be found in Appendix C. The first is a closed form, a hypergeometric function with complex poles at the cyclotron resonances,

(3.25)\begin{align} \varepsilon(\omega, k_\perp) & = 1 + \left(\frac{\omega_p}{\omega_c}\right)^2\frac{1}{(k_\perp r_L)^2} \left\{ {}_{2}F_{2}\genfrac[]{0pt}{}{\frac{1}{2}, \qquad\gamma+1}{1+\omega/\omega_c, \ 1-\omega/\omega_c}({-}2(k_\perp r_L)^2)\right.\nonumber\\ & \quad\left.-\,{}_{2}F_{2}\genfrac[]{0pt}{}{\frac{1}{2},\qquad \gamma}{1+\omega/\omega_c,\ 1-\omega/\omega_c}({-}2(k_\perp r_L)^2) \right\} = 0 \end{align}

and the second form is a representation of (3.25) as a trigonometric integral generalizing that used in Tataronis & Crawford (Reference Tataronis and Crawford1970a), Vogman, Colella & Shumlak (Reference Vogman, Colella and Shumlak2014) and Datta, Crews & Shumlak (Reference Datta, Crews and Shumlak2021),

(3.26)\begin{equation} \varepsilon(\omega, k) = 1 + \left(\frac{\omega_p}{\omega_c}\right)^2\csc({\rm \pi}\omega/\omega_c) \int_0^{\rm \pi}\sin(\theta)\sin(\theta\omega/\omega_c)L_\gamma(\beta)\exp(-\beta)\,{\rm d}\theta \end{equation}

with $\beta \equiv 2\cos ^2(\theta /2)(k_\perp r_L)^2$ and $L_\gamma (\beta )$ the Laguerre polynomial of order $\gamma$. Equation (3.26) is particularly suited to numerical calculation by quadrature. Hypergeometrics similar to (3.25) have been obtained for the Maxwellian plasma and for $\kappa$-distributions (Mace Reference Mace2004; Mace & Hellberg Reference Mace and Hellberg2009), and reduce to the Maxwellian result for the parameter $\gamma = 0$. Observe that (3.25) and (3.26) are functions of frequency and not of phase velocity as there is no ballistic contribution to the zero-order motion. It is expected that (3.24)–(3.26) will be the basic dispersion functions from which to build the $q$-analogues for $\chi$-distributions proposed by Leubner & Schupfer (Reference Leubner and Schupfer2001).

3.5. Visualization of the dispersion function and phase space eigenfunctions

Figure 11 plots the electrostatic dispersion function in the complex frequency plane for a ring distribution of $\gamma =6$ for magnetization $\omega _c/\omega _p=0.1$ at wavenumber $kr_L=0.9$. A plethora of solutions to the complex dispersion function $\varepsilon =0$ is illustrated by the intersection of the zeros of the real and imaginary parts, which take place at both solutions and poles. Complex poles at the cyclotron harmonics occur only for the case $k_\parallel \to 0$, as evident from the cosecant function $\csc ({\rm \pi} \omega /\omega _c)$ in (3.26). In this way solutions and simple poles can be clearly distinguished. It is clear that, given zero-order cyclotron orbits, oblique modes are Landau damped but perpendicular modes are not. When wave amplitudes violate the inequality of (3.14) nonlinear phenomena occur, and perpendicular waves are also Landau damped. In this situation one may see streaming instabilities in numerical experiments in the magnetization transition regime $\omega _c \approx \omega _p$.

Figure 11. Contours of $\text {Re}(\varepsilon )=0$ (green) and $\text {Im}(\varepsilon )=0$ (red) for the electrostatic dispersion function of a loss-cone-distributed plasma with parameters $\omega _p = 10\omega _c$, $\gamma =6$, and $k_\perp r_L=0.9$ for case (a) $90^\circ$ propagation and case (b) $85^\circ$ propagation. Either a solution or a pole occurs at an intersection of these two curve families. Solutions to the dispersion relation $\varepsilon (\omega, k)=0$ are labelled $\mathbb {S}$ while simple poles are labelled $\mathbb {P}$. There are no poles for oblique $(< 90^\circ )$ propagation. Solutions correspond to the first few electron Bernstein modes. The dispersion function shows Gaussian-like responses in the lower-half plane for the oblique wave, an indicator of dissipative wave–particle resonance near the cyclotron harmonics.

Therefore, perpendicularly propagating linear modes do not experience Landau resonance so that all such modes are true eigenfunctions. The phase space structure associated with these cyclotron waves (that is, (3.3)) consists of helical modes in the perpendicular velocity space, since the primary phase component is $\exp ({\rm i}(kx + n\phi - \omega t))$ for a mode with frequency $\omega \approx n\omega _c$, with $\phi$ the cylindrical velocity space angle. The simplest example of such helical phase space modes are the electron Bernstein modes. Figure 12 visualizes the eigenfunctions of the first and second cyclotron harmonics for a Maxwellian background distribution $f_0(v_\perp )$ in the phase space $(x, u, v)$ with $(u, v)$ the perpendicular velocity space and $x$ the propagation coordinate perpendicular to the background magnetic field. With non-zero real frequency these helical modes propagate through phase space.

Figure 12. Phase space eigenfunctions of (a) the first and (b) the second cyclotron-harmonic electron Bernstein waves propagating orthogonal to $\boldsymbol {B}_0$, visualized in the phase space $(x, u, v)$. The order of the cyclotron harmonic determines the number of helical strands of a given sign. The phase space eigenfunction propagates in the laboratory frame, but when the velocity space is observed at a fixed spatial coordinate the eigenfunction rigidly rotates in the perpendicular velocity space. Balanced counterpropagating modes produce a stationary wave.

3.6. Simulation of perpendicular electron cyclotron loss-cone instability

Here kinetic eigenfunction initialization is illustrated for the instability of an electron loss cone to perpendicular cyclotron waves in a neutralizing background, as studied in Vogman et al. (Reference Vogman, Colella and Shumlak2014). Normalizing to the Debye length, plasma frequency and thermal velocity, and the fields by $E_0 = v_{th}B_0$, the Vlasov–Poisson equations are

(3.27)\begin{gather} \partial_{t}f + F^j\partial_j f = 0, \end{gather}

(3.28)\begin{gather} \frac{{\rm d}^2\varphi}{{{\rm d}\kern0.7pt x}^2} = \int_{-\infty}^\infty\int_{-\infty}^\infty f(x,u,v) \,{\rm d}u\,{\rm d}v - 1,\end{gather}

(3.29)\begin{gather} F = \begin{bmatrix} u, & \dfrac{{\rm d}\varphi}{{\rm d}\kern0.7pt x} - vB_{\text{ext}}, & uB_{\text{ext}} \end{bmatrix}^{\rm T} \end{gather}

with coordinates $(x,u,v)$. The external magnetic field $B_{\text {ext}}$ is set such that the magnetization parameter $\omega _c / \omega _p = B_\text {ext}$ is in normalized units. Two single-mode simulations termed $A$ and $B$ are performed for (3.27)–(3.29) in the highly unstable over-dense parameter regime $\omega _p = 10\omega _c$ using as eigenvalues two solutions to $\varepsilon (\omega,k)=0$, namely $(k_A\lambda _D,\omega _A /\omega _c)\approx (0.886, 0.349{\rm i})$ and $(k_B\lambda _D, \omega _B/\omega _c)\approx (1.4, 1.182 + 0.131{\rm i})$ found using the integral form of (3.26) with 50-point Gauss–Legendre quadrature.

In this situation, case $A$ corresponds to a stationary mode with wavelength long compared with the thermal Larmor radius, similar to the two-stream instability studied in the unmagnetized case, while case $B$ corresponds to a destabilized propagating Bernstein-like mode at the first cyclotron harmonic with more significant finite Larmor radius effect. The spatial domain is set to $L = 2{\rm \pi} / k$ in each case and the velocity boundaries to $u_{\max }, v_{\max }= \pm 8.5$. We perturb these nonlinear simulations using the phase space eigenfunctions corresponding to the eigenvalue pairs $(k_A, \omega _A)$, $(k_B, \omega _B)$.

Figure 13 shows isocontours of the phase space eigenfunctions used as the initial perturbations. Case $A$ has a phase space structure similar to the basic plasma wave seen in the two-stream instability, while case $B$ has a helical structure as a cyclotron mode with $\omega _r\approx \omega _c$. We reiterate here that the eigenfunction perturbation allows arbitrary perturbation amplitude and still produces the same nonlinear phenomena, namely the saturated state or mode coupling/conversion. However, with zero-order $\boldsymbol {B}_0$ the initial amplitude must not exceed the limit of (3.14) or nonlinear phenomena will develop as the perturbation electric force is not first order compared with the thermal magnetic force.

Figure 13. Shown here are the perturbations $f_1 \equiv f - f_0$ at time $t=0$ for (a) simulation $A$ with $\omega _r = 0$ and (b) simulation $B$ with $\omega _r\neq 0$, each with isosurfaces at $30\,\%$ of the minimum (green) and maximum (yellow). The considered eigenfunctions $f_1(x,u,v)$ consist of twisting islands in the phase space, capturing the combined physics of translation, electric acceleration and magnetic gyration. Translating modes ($\omega _r\neq 0$) are associated with a helical structure.

3.6.1. Numerical method for loss-cone simulations

The problem is evolved numerically using the DG method described in Appendix A and in Crews & Shumlak (Reference Crews and Shumlak2022), with the difference that the spatial coordinate is not Fourier transformed but also discretized by DG method. We use an element resolution $(N_x, N_u, N_v) = (25, 50, 50)$ and nodal basis of $n=8$ Legendre–Gauss–Lobatto nodes per dimension, while the Shu–Osher SSPRK3 method is used to integrate the semidiscrete equation in time. These instabilities grow on a slow time scale relative to the plasma frequency; that is, they grow at a fraction of the cyclotron time scale $\omega _c^{-1}$, while time $t$ is normalized to the plasma frequency $\omega _p^{-1}$. Thus, these instabilities take many plasma periods to reach nonlinear saturation beginning from amplitudes below the limit of (3.14). Simulation $A$ reaches saturation around $t=100$ and runs to $t=175$ while simulation $B$ saturates at around $t=175$ and stops at $t=200$.

Three-dimensional isosurface plots were produced using PyVista, a Python package for VTK (visualization toolkit). To prepare the data, an average is taken of nodal values lying on element boundaries for smoothness, and the eight-nodes per element are resampled to $25$ linearly spaced points per axis and per element onto the basis functions of the DG method. These isocontours are shown for simulations $A$ and $B$ in figures 14 and 15, respectively.

Figure 14. Phase space view $(x,u,v)$ of simulation $A$ focused on $(x,u)$-plane as isocontours at $15\,\%$ of $\max (f)$ shown in yellow, at times (a) $t=0$, (b) $t=80$ and (c) $t=120$. The domain within the original ring is shown with $(u,v)\in (-3.5, 3.5)$ to focus on the trapping dynamics. The trapping structure consists of a ribbon winding around a separatrix, while the outer bulk ring distribution maintains passing trajectories. This saturation geometry is typical for electrostatic potentials in a magnetic field as the magnetic force depends on the sign of the transverse velocity.

Figure 15. Phase space view looking on $(-x,v)$-plane of simulation $B$ at $15\,\%$ isocontours of $\max (f)$ (yellow), with (a) the nonlinear mode developing at $t=100$, (b) the developed vortex at $t=160$ and (c) the saturated vortex translating at $t=180$. The mode is seen to be a growing, translating electrostatic potential $\varphi (x)$ of positive phase velocity $v_{\varphi } = \omega _r / k$ with an underlying phase space vortex structure centred at $(u,v)=0$. The vortex shape is explained by considering the trajectory of a test particle in the wave. That is, particles with a velocity close to that of the wave see a stationary potential and are accelerated to a high $u$-velocity. They then translate towards positive $x$ while their velocity vector is rotated by the Lorentz force to $-u$ at a rate close to the wave frequency (as $\omega _r \approx 1.2\omega _c$) and repeats the cycle.

3.6.2. Fully nonlinear single-mode loss-cone simulation results

Figure 16 shows the electric potentials $\varphi (x)$ in simulations $A$ and $B$, computed using the numerical method described in § 3.6.1. In simulation $A$ the wave potential $\varphi (x)$ is stationary with a weakly fluctuating boundary, so that part of the density $f(x,v)$ within the potential well executes trapped orbits. This results in a trapping structure with orbits tracing a nonlinear potential similar to the characteristic pendulum-like cats-eye separatrix of the single-mode electrostatic two-stream instability. In this case the electrons are magnetized and execute zero-order cyclotron motion so that the trapping separatrix in simulation $A$ is instead in the form shown by the isosurfaces. The saturated state of simulation $A$ is a lattice of one-dimensional strongly magnetized electron holes. The continuum of such hole solutions is key to strongly driven plasma transport physics (Schamel Reference Schamel2023). Two-dimensional axisymmetric magnetized holes are the focus of analytical work in Hutchinson (Reference Hutchinson2020) and Hutchinson (Reference Hutchinson2021).

Figure 16. Electric potentials $\varphi (x)$ at saturation of the two studied cases, for (a) simulation $A$ at $t=120$ and (b) simulation $B$ at $t=180$. The potential of $A$ is stationary while that of $B$ is translating to the right. The negative of potential $-\varphi (x)$ is shown in order to account for the electron's negative charge. In both cases electron holes develop in the potential wells of $-\varphi (x)$.

The saturated wave potential of simulation $B$, on the other hand, translates with positive phase velocity $\zeta \approx \omega _r / k$. The region of particle interaction translates along with the wave potential and forms a vortex structure in the phase space density $f(x, u, v)$. The centre of this kink continues to tighten as the simulation progresses, leading to progressively finer structures just as in simulation $A$. This effect is in agreement with the filamentation phenomenon and introduces simulation error as the structures lead to large gradients on the grid scale where discreteness produces dispersion error. For this reason the simulation is stopped at $t = 200$. This solution is perhaps somewhat artificial as it is obtained by single-mode initialization via its phase space eigenfunction, thereby not disturbing modes of greater growth rate. This is demonstrated through the energy traces in figure 17, showing that simulation $A$ saturates with a greater proportion of the plasma thermal energy than simulation $B$. The solution of simulation $B$ has the significance of a propagating train of nonlinear electrostatic cyclotron waves with associated electron holes, driven unstable by the free energy of a ring distribution.

Figure 17. Domain-integrated electric field energy traces for simulations (a) $A$ and (b) $B$. The thermal energy of the zero-order distribution is $0.25$ per unit length with $\gamma =6$ and $\alpha = 1$. In simulation $A$ this corresponds to a domain-integrated thermal energy of $E_A\approx 17.8$ and in simulation $B$ to $E_B\approx 11.2$. Therefore, in both cases the instability saturates with an electric energy a few per cent of the thermal energy, approximately $6\,\%$ in $A$ and $2.5\,\%$ in $B$.

3.6.3. Experimental consequences of cyclotron loss-cone instabilities

Magnetic mirror trapping requires the maintenance of a loss-cone distribution in the confined plasma. Simulations such as these, and QLTs, maintain that kinetic instabilities lead to a relaxation of the distribution function on macroscopic scales towards near-Maxwellian distributions, along with long-lived vortical structures in the phase space. By examination of the dispersion functions one finds that phase space instability may be suppressed when $\omega _p \ll \omega _c$. Assuming equal electron and ion temperatures and densities, we may write for the plasma beta

(3.30)\begin{equation} \beta = \left(\frac{v_{ti}}{c}\right)^2\left(\frac{\omega_{pi}}{\omega_{ci}}\right)^2 = \left(\frac{r_{Li}}{\lambda_i}\right)^2 \end{equation}

with $\lambda _i = c/\omega _{pi}$ the ion skin depth. Thus, non-Maxwellian features such as ring distributions, as $\chi$-distributions or their $\kappa$-analogues (Pokhotelov et al. Reference Pokhotelov, Treumann, Sagdeev, Balikhin, Onishchenko, Pavlenko and Sandberg2002; Leubner Reference Leubner2004), are expected to persist in very low-$\beta$ plasma, in much the same way in which for weakly magnetized plasma the $\kappa$-distributions persist in the collisionless regime when the plasma parameter $\varLambda \gg 1$. Further discussion on the consequences of mirror instability in high-$\beta$ space plasma can be found in Pokhotelov et al. (Reference Pokhotelov, Sagdeev, Balikhin and Treumann2004).

4.1. Tensor components for arbitrary Cartesian coordinates

In Cartesian coordinates $(x,y,z)$ with wavevector $\boldsymbol {k} = k_x\hat {x} + k_y\hat {y} + k_z\hat {z}$, the dielectric tensor (Skoutnev et al. Reference Skoutnev, Hakim, Juno and TenBarge2019) is obtained from combination of (4.6) and (4.3) as

(4.7)\begin{equation} \left. \begin{gathered} \varepsilon_{ij} E_j = 0,\\ \varepsilon_{ij} = \delta_{ij} - \frac{c^2k^2}{\omega^2}\left(\delta_{ij} - \frac{k_ik_j}{k^2}\right) - \frac{\omega_p^2}{\omega^2}\left(\delta_{ij} - \int v_i v_j \frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_0}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\,{\rm d}\boldsymbol{v}\right). \end{gathered} \right\} \end{equation}

In the formal initial-value problem an initial-value vector results in the system $\varepsilon _{ij}E_j=g_i$. Typically in Cartesian form the moment integrals are inseparable because of the resonant denominator $\omega - \boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}$, yet are separable when the frame is chosen with one coordinate aligned with the wavevector such that the resonant denominator appears as $\omega -kv_\parallel$.

4.2. Electromagnetic susceptibility and the eigenvalue problem

Reformulation of (4.7) as an eigenvalue problem for the phase velocity $\zeta$ allows calculation of electric field eigenfunctions $\boldsymbol {E}$ naturally consistent with the corresponding phase space eigenfunction $f_{1}$ given by (4.5). This facilitates construction of initial conditions for simulation of Vlasov–Maxwell instabilities. Multiply by $\omega ^2$ and pull out the diagonal tensor $\omega ^2 \delta _{ij}$. Define the integrals encoding resonant wave–particle interaction into the self-consistent perturbation current as

(4.8)\begin{equation} {\rm I}_{ij} \equiv \int_{\mathcal{C}} v_i v_j\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_0}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\,{\rm d}\boldsymbol{v}, \end{equation}

where the integral is evaluated on the Landau contour $\mathcal {C}$, i.e. analytically continued to the lower-half complex $\omega$-plane. The dielectric tensor system may then be rewritten as

(4.9)\begin{equation} \left\{\delta_{ij} - {\rm I}_{ij} + (k\lambda_c)^2\left(\delta_{ij} - \frac{k_ik_j}{k^2}\right)\right\}E_j = \left(\frac{\omega}{\omega_p}\right)^2 E_i, \end{equation}

where $\lambda _c = c/\omega _p$ is the inertial length. The resonant integrals are naturally functions of the phase velocity $\zeta =\omega /k$, so one can also express the system as

(4.10)\begin{equation} \frac{1}{(k\lambda_D)^2}\left\{\delta_{ij} - {\rm I}_{ij} + (k\lambda_c)^2\left(\delta_{ij} - \frac{k_ik_j}{k^2}\right)\right\}E_j = \left(\frac{\zeta}{v_t}\right)^2 E_i. \end{equation}

Equation (4.10) casts the problem in eigenvalue form as the system of integral equations

(4.11)\begin{equation} \chi(k, \zeta)\boldsymbol{E} = \zeta^2\boldsymbol{E}. \end{equation}

As in the scalar Poisson problem there is a spectrum of solutions $\zeta$ for a given $k$, obtained by determining a root of the characteristic function $D(k, \zeta ) \equiv \textrm {det}(\chi - \zeta ^2 I) = 0$. As in the scalar problem only unstable solutions satisfying $\textrm {Im}(\zeta ) > 0$ constitute normal modes of oscillation. In unmagnetized spatially uniform plasma these modes are either streaming instabilities (two-stream, Buneman, ion-acoustic) or generalized Weibel instabilities. In either case their effect is thermalizing on coarse-grained scales by reducing relative velocities far from equilibrium. Casting the dielectric tensor for zero-order cyclotron motion (Gurnett & Bhattacharjee Reference Gurnett and Bhattacharjee2017, p. 408) into an eigenvalue problem proceeds in the same manner; the main difference is that the integrals $\textrm {I}_{ij}$ are sums over Doppler-shifted cyclotron resonances and their calculation, though methodical, is lengthy.

Having determined a particular eigenvalue $\zeta _n^2$ such that $D(k, \zeta _n)=0$, the corresponding electric field eigenfunction $\boldsymbol {E}_n$ is found by solving for the eigenvector of the matrix $\chi _n \equiv \chi (\zeta _n, k)$ with eigenvalue $\zeta _n^2$. The other two eigenvalues of the matrix $\chi _n$ are spurious as they do not correspond to solutions of (4.11). The magnetic field eigenfunction is obtained through $\boldsymbol {B}_n = \zeta _n^{-1}\hat {\boldsymbol {k}}\times \boldsymbol {E}_n$, and the phase space eigenfunction $f_{1,n}$ from (4.5).

4.3. Dielectric tensor components for the anisotropic Maxwellian

In order to illustrate the Weibel instability due to anisotropy in a plasma with zero-order ballistic trajectories it is useful to consider the anisotropic Maxwellian, or multiple-temperature, zero-order distribution (Davidson et al. Reference Davidson, Hammer, Haber and Wagner1972). Anisotropy is a key element of kinetic equilibrium in flowing magnetized plasmas (Mahajan & Hazeltine Reference Mahajan and Hazeltine2000), and the anisotropic Maxwellian is the simplest anisotropic model distribution to investigate pure Weibel instability. Let $(x,y,z)$ and $(u,v,w)$ be Cartesian coordinates in configuration space and velocity space, respectively. Consider a three-temperature Maxwellian distribution given by

(4.12)\begin{equation} f_0 = \frac{1}{(2{\rm \pi})^{3/2}\theta_u\theta_v\theta_w}\exp\left\{-\frac{1}{2}\left(\frac{u^2}{\theta_u^2} + \frac{v^2}{\theta_v^2} + \frac{w^2}{\theta_w^2}\right)\right\}. \end{equation}

For purpose of illustration we obtain from (4.12) one- and two-dimensional model problems of Weibel instability by letting the wavevector lie in the $(x,y)$-plane and the $\hat {z}$-direction be the wave binormal. These model problems have the advantage of zero-order electrostatic stability such that all eigenfunctions are fully electromagnetic.

With the wavevector in the $(x,y)$-plane the off-diagonal integrals $\textrm {I}_{13}$, $\textrm {I}_{23}$ vanish as $\langle w\rangle =0$, so that the $\hat {z}$-directed perturbation does not contribute to the $(x,y)$-plane's perturbation current, decoupling the binormal from the longitudinal and transverse components (Sharma & Bhatnagar Reference Sharma and Bhatnagar1976; Datta et al. Reference Datta, Crews and Shumlak2021). A fully general wavevector would couple all three components of the perturbation. Equation (4.7) simplifies to

(4.13)\begin{align} \begin{Bmatrix} \omega^2 - \omega_p^2(1 - {\rm I}_{11}) & \omega_p^2 {\rm I}_{12} & 0\\ \omega_p^2 {\rm I}_{12} & \omega^2 - c^{2}k^2 - \omega_p^2(1 - {\rm I}_{22}) & 0\\ 0 & 0 & \omega^2 - c^2 k ^2 - \omega_p^2 (1 - {\rm I}_{33}) \end{Bmatrix} \begin{Bmatrix} E_1\\ E_2\\ E_3 \end{Bmatrix} = 0. \end{align}

Focusing on the $(x,y)$-plane we consider a reduced distribution, the bi-Maxwellian $f_0(u,v)=\int f_0(u,v,w)\,\textrm {d}w$ whose level sets form ellipses in the $(u,v)$ velocity plane. Take $\theta _v > \theta _u$ so that the semimajor axis of each ellipse is $v$-directed and the characteristic eccentricity is $e^2=1 - \theta _u^2/\theta _v^2$. To ensure separability of the resonant integrals $\textrm {I}_{ij}$, the $(x,y)$-plane is rotated through an angle $\varphi$, transforming velocities $(u,v)\to (v_\parallel,v_\perp )$ such that the resonant denominator is $\omega - k v_\parallel$. By completing the square on $v_\perp$, the anisotropic Maxwellian of (4.12) in the reduced coordinates $(u,v)\to (v_\parallel,v_\perp )$ is

(4.14)\begin{equation} f_0(v_\parallel, v_\perp) = \frac{1}{2{\rm \pi}\theta_\parallel\theta_\perp} \exp\left(-\frac{v_\parallel^2}{2\theta_\parallel^2}\right)\exp\left(-\frac{(v_\perp{-} \alpha v_\parallel)^2}{2\theta_\perp^2}\right), \end{equation}

where the rotated thermal and mean velocities are defined as

(4.15a-c)\begin{equation} \theta_{{\parallel}}^2 \equiv \frac{(1 - e^2)\theta_u^2}{1-e^2\cos^2(\varphi)},\quad \theta_{{\perp}}^2 \equiv \frac{\theta_u^2}{1 - e^2\cos^2(\varphi)},\quad \alpha \equiv \frac{e^2\sin(\varphi)\cos(\varphi)}{1-e^2\cos^2(\varphi)}. \end{equation}

Equation (4.14) shows that in the rotated frame the distribution is a shifted Maxwellian. That is, in the frame of a resonant particle the distribution has a mean velocity $\alpha v_\parallel$ transverse to the wavevector, as illustrated in figure 18. The off-diagonal integrals are non-zero, $\textrm {I}_{ij}\neq 0$, for anisotropic distributions as the apparent transverse current in the resonant particle frame couples the longitudinal and transverse wave components for propagation not aligned with the principal axes, even though there is no net current in the laboratory frame. In general all three components are coupled for non-principal propagation. This coupling is characteristic of anisotropy and contrasts with isotropic distributions for which the longitudinal plasma wave and the two transverse electromagnetic wave components have fully independent dispersion relations.

Figure 18. In the frame of reference of a particle propagating in a direction not aligned with the principal axes of an anisotropic distribution there is a mean drift velocity in the direction transverse to the particle's motion which couples together the electromagnetic wave's longitudinal and transverse components. In other words, in a frame where the velocity coordinates $(u,v)$ are rotated through an angle $\varphi$ into coordinates $(v_\parallel, v_\perp )$, for a particular value of $v_\parallel$ there is a non-zero mean value of $v_\perp$ such that the off-diagonal integrals $\textrm {I}_{ij}\neq 0$ of (4.8).

To evaluate the integrals $\textrm {I}_{11}$, $\textrm {I}_{12}$ and $\textrm {I}_{22}$ for the distribution in (4.14) in the rotated coordinates, note the following integrals related to the plasma dispersion function $Z(\zeta )$:

(4.16)\begin{gather} \frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^{\infty}\frac{x}{\zeta-x}{\rm e}^{{-}x^2/2a^2}\,{{\rm d}\kern0.7pt x} = \frac{a}{2}Z'(\tilde{\zeta}), \end{gather}

(4.17)\begin{gather}\frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^{\infty}\frac{x^2}{\zeta-x}{\rm e}^{{-}x^2/2a^2}\,{{\rm d}\kern0.7pt x} ={-}\frac{a^2}{2}(Z''(\tilde{\zeta}) + 2Z(\tilde{\zeta})), \end{gather}

(4.18)\begin{gather}\frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^{\infty}\frac{x^3}{\zeta-x}{\rm e}^{{-}x^2/2a^2}\,{{\rm d}\kern0.7pt x} = \frac{a^3}{2^{3/2}}(Z'''(\tilde{\zeta}) + 6Z'(\tilde{\zeta})), \end{gather}

where $\tilde {\zeta }=\zeta / \sqrt {2} a$. These identities are found through integration by parts using the Hermite relation $\psi _n(x)=({\textrm {d}^n}/{{\textrm {d} x}^n})\textrm {e}^{-x^2/2}$ and the identity $Z'(\zeta )=-2(1+\zeta Z(\zeta ))$. With the gradient ${\partial f}/{\partial v_\parallel } = (-v_\parallel /\theta _\parallel ^2 + \alpha (v_\perp -\alpha v_\parallel )/\theta _\perp ^2)f$ the integrals $\textrm {I}_{ij}$ in (4.13) work out to

(4.19)\begin{gather} {\rm I}_{11} ={-}\frac{1}{4}(Z'''(\tilde{\zeta}) + 6Z'(\tilde{\zeta})), \end{gather}

(4.20)\begin{gather}{\rm I}_{12} ={-}\frac{\alpha}{4}(Z'''(\tilde{\zeta}) + 6Z'(\tilde{\zeta})), \end{gather}

(4.21)\begin{gather}{\rm I}_{22} ={-}\frac{\alpha^2}{4}(Z'''(\tilde{\zeta}) + 2Z'(\tilde{\zeta})) - \frac{\theta_\perp^2}{\theta_\parallel^2}\frac{Z'(\tilde{\zeta})}{2}, \end{gather}

(4.22)\begin{gather}{\rm I}_{33} ={-}\frac{\theta_w^2}{\theta_\parallel^2}\frac{Z'(\tilde{\zeta})}{2}, \end{gather}

where $\tilde {\zeta }\equiv \zeta /\sqrt {2}\theta _\parallel$. Let the angular eccentricity be $\gamma = \arcsin (e)$. The characteristic anisotropies are then $A_1\equiv \theta _\perp ^2/\theta _\parallel ^2-1={e^2}/({1-e^2})=\tan ^2(\gamma )$ for the in-plane anisotropy and $A_2\equiv \theta _w^2/\theta _\parallel ^2-1$ for the binormal anisotropy. These two anisotropy parameters induce the electromagnetic Weibel instability to relax the anisotropy of their respective dimensions. For propagation along the principal axes, $\varphi \to 0, {\rm \pi}/2$ the parameter $\alpha \to 0$, decoupling the $\hat {x}$- and $\hat {y}$-components as $\textrm {I}_{12}\to 0$ and $\textrm {I}_{22} \to -({\theta _v^2}/{\theta _u^2})({Z'(\zeta /\sqrt {2}\theta _u)}/{2})$ as in (4.22).

4.3.1. Dispersion function for the classic Weibel instability

Focusing on the decoupled component in (4.13), namely $\varepsilon _{33}E_3=0$, leads to

(4.23)\begin{equation} 1 - \left(\frac{\zeta}{c}\right)^2 + \frac{1}{(k\lambda_c)^2}\left(1 + \frac{Z'(\tilde{\zeta})}{2} + A_2\frac{Z'(\tilde{\zeta})}{2}\right) = 0 \end{equation}

with $A_2$ the binormal anisotropy. Equation (4.23) is valid for all propagation angles. For principal propagation at $\varphi =0$ (the $\hat {x}$-direction) the principal anisotropies $A_1=A_2$ such that (4.23) describes both the transverse and binormal components. It is instructive to observe that isotropy reduces (4.23) to the ordinary wave kinetic dispersion relation. The coupled branch of solutions to (4.13) is described by

(4.24)\begin{equation} \left(-\left(\frac{\zeta}{c}\right)^2 + \frac{1}{(k\lambda_c)^2}(1 - {\rm I}_{11})\right) \left(1 - \left(\frac{\zeta}{c}\right)^2 + \frac{1}{(k\lambda_c)^2}(1 - {\rm I}_{22})\right) - \frac{1}{(k\lambda_c)^4}{\rm I}_{12}^2 = 0 \end{equation}

with transverse magnetic field out-of-plane and two components of electric field in-plane.

4.4. Single-mode saturation of Weibel instability in one spatial dimension

Just as electrostatic instability saturates by nonlinear trapping of near-resonant particles in an electric potential energy well, electromagnetic instability saturates by magnetic trapping. As a supplement to this section, Appendix D builds an analytic model of the phase portraits associated with ideal magnetic trapping and a conceptual model of magnetic trapping as a magnetic potential momentum well. Simulation of a single-mode Weibel instability by initialization with kinetic eigenfunctions allows one to self-consistently visualize the saturated phase space structures. The foundational simulations of an electromagnetic instability due to anisotropy considered Weibel instability evolving from the anisotropic Maxwellian using the particle-in-cell method (Davidson et al. Reference Davidson, Hammer, Haber and Wagner1972). Later, Califano, Pegoraro & Bulanov (Reference Califano, Pegoraro and Bulanov1997) conducted one- and two-dimensional simulations in an inhomogeneous plasma using anisotropy induced by streaming beams. Building on this, Cagas et al. (Reference Cagas, Hakim, Scales and Srinivasan2017) simulated single-mode Weibel saturation with a continuum-kinetic method using a zero-order counter-streaming electron beam distribution. In Cagas et al. (Reference Cagas, Hakim, Scales and Srinivasan2017) an electrostatic streaming instability was proposed to explain the growth of a beam-axis directed electric field close to nonlinear saturation.

Zero-order beam distributions are often used because Weibel instability is induced in the laboratory by colliding high velocity plasmas (Hill et al. Reference Hill, Key, Hatchett and Freeman2005; Fox et al. Reference Fox, Fiksel, Bhattacharjee, Chang, Germaschewski, Hu and Nilson2013; Huntington et al. Reference Huntington2015; Shukla et al. Reference Shukla, Vieira, Muggli, Sarri, Fonseca and Silva2018). While zero-order beam distributions are inherently anisotropic they are also possibly unstable to electrostatic streaming instability. The possible introduction of electrostatic streaming instability can confuse and complicate an attempt to isolate Weibel instability. On the other hand, there is no possibility of streaming instability when the Weibel instability is induced by a zero-order anisotropic Maxwellian distribution. For this reason the anisotropic Maxwellian is considered here using the continuum kinetic method, as originally considered with particle-in-cell method by Davidson et al. (Reference Davidson, Hammer, Haber and Wagner1972).

4.4.2. Initialization with field and phase space eigenfunctions

The characteristic parameters are chosen such that $v_t/c=0.3$ and anisotropy $A = 3$ (or ratio $\theta _{v}/\theta _{u}=2$) with the zero-order distribution given by (4.12) and the direction of propagation set to $\varphi =0$. This is equivalent to using (4.23) for the dispersion function. Velocities are normalized to $\theta _{u}$, time to $\omega _p^{-1}$ and lengths to $\lambda _D=\theta _{u}/\omega _p$. The domain length is then specified by the chosen wavenumber $k_x\lambda _D = 0.1$. Solution of (4.23) gives the eigenvalue of the problem as the phase velocity $\zeta /v_t = 1.23\textrm {i}$. The phase space perturbation is constructed using (4.5) for the phase space eigenfunction in the form

(4.25)\begin{equation} f_1(x, v_x, v_y) = \text{Re}\left\{ {\rm i}\frac{\mathcal{A}}{k}\left(\frac{v_y}{\zeta - v_x}\frac{\partial f_0}{\partial v_x} + \frac{\partial f_0}{\partial v_y}\right) \exp({\rm i}kx)\right\}, \end{equation}

longitudinal field $E_{x1}=0$ and the initial transverse electrodynamic fields by

(4.26)\begin{gather} B_{z1} = \text{Re}\left(\mathcal{A} \exp({\rm i}kx)\right), \end{gather}

(4.27)\begin{gather}E_{y1} = \text{Re}\left(\zeta \mathcal{A} \exp({\rm i}kx)\right), \end{gather}

where the amplitude is set to $\mathcal {A} = 10^{-3}$. This initial condition is consistent in the sense that the charge density is zero and the current density satisfies Ampere's law.

4.4.3. Fully nonlinear simulation results for single-mode Weibel saturation

The simulation is run until $t\omega _p=100$ using the numerical method of § 4.4.1 with time step $\Delta t = 2.0\times 10^{-3}$. The change in domain-integrated wave energies is shown in figure 19. Of note is the oscillating transverse electric field at saturation, and how longitudinal electric energy corrects from zero to trend with magnetic energy.

Figure 19. Growth of domain-integrated wave energies towards saturation of Weibel instability in one spatial dimension, namely (a) the magnetic energy, (b) the transverse electric energy of $E_y$ and (c) the longitudinal electric energy of $E_x$, or energy along the axis of the wavevector. A nonlinear phase is reached at $t\omega _p=40$, with peak magnetic energy around $t\omega _p=50$. Longitudinal electric energy is observed to grow with a similar trend to the magnetic energy.

Figure 20 shows the time evolution of magnetic field and electron density in increments of $t\omega _p=20$. As magnetic energy grows, electrons are progressively magnetically trapped as described analytically in Appendix D. With the initial free energy released, the distribution function attempts to evolve towards a function of the constants of motion, namely the energy $H = \tfrac {1}{2}m(v_x^2 + v_y^2) - e\varPhi (x)$ and canonical momenta $P_y = mv_y - e A_y(x)$, $P_x = mv_x$ where $\varPhi (x)$ is the electric potential and $\boldsymbol {A}=A_y(x)\hat {y}$ is the magnetic vector potential. Specifically, trapped and passing phase space trajectories are determined by the equation (Morse & Nielson Reference Morse and Nielson1971)

(4.28)\begin{equation} mv_x = \sqrt{2m(H + e\varphi(x)) - (P_y + eA_y(x))^2}, \end{equation}

for a particle's energy $H=\frac {1}{2}m(v_{x0}^2 + v_{y0}^2) - e\varphi (x_0)$ and momentum $P_y = mv_{y0} - eA_{y0}(x_0)$ constants. Figure 21 visualizes the phase space structure at nonlinear saturation by an isocontour of the distribution at $10\,\%$ of its maximum. Trapped and passing trajectories are seen at the right- and left-hand of the figure, respectively, for $v_y>0$. Trapped trajectories circulate within a phase space vortex while passing trajectories execute motion in the vicinity of the separatrix. Upon inversion of the transverse velocity, $v_y\to -v_y$, the positions of the trapped trajectories and passing trajectories are inverted, $x\to -x$.

Figure 20. Evolution of (a) magnetic field $B_z(x)$ and (b) electron density $n_e(x)$, to nonlinear saturation of a single unstable Weibel mode. The mode saturates around $t\omega _p\approx 45$. Prior to saturation the magnetic field has a spectrum consisting of only even mode numbers with an apparent power law in logarithmic amplitudes. Since the function is clearly analytic this spectrum is consistent with an elliptic cosine function until nonlinear saturation. It is interesting that this higher-order phenomenon arises even from a single-mode eigenfunction perturbation. The electron holes at saturation are bounded by the maxima of magnetic energy $B^2/2\mu _0$, or equivalently bounded by the potential wells $A_y(x)$.

Figure 21. Phase space vortex shown by contour at $0.1\max (f)$ in the coordinates $(x, v_x, v_y)$ at saturation of a single Weibel mode. The mid-domain $x$-coordinate corresponds to $x=0$. The magnetic trapping phase space vortex is distinguished from the electrostatic vortex as velocity-dependent since the momentum $p_y$ in (4.28) is linear in $v_y$. In this way the single-mode phase space vortex is antisymmetric, with trapped and passing orbits swapping under $v_y\to -v_y$. For this reason there is no magnetic trapping along the $v_y=0$ plane.

Linear electrostatic instability is not possible due to the zero-order anisotropic Maxwellian, so here we explain the growth of electrostatic energy observed both here and in Cagas et al. (Reference Cagas, Hakim, Scales and Srinivasan2017) as a second-order phenomenon arising from space-charge-generating filamentation (that is, the magnetic-trapping electron holes of figure 20). While saturated filaments can be understood intuitively as electron holes, the progressive development of longitudinal electric energy in the linear phase can be understood as mode coupling of the longitudinal field to the transverse field at second order in the transverse dynamic field (Taggart et al. Reference Taggart, Godrey, Rhoades and Ives1972). Appendix D presents a visual description of the density filamentation phenomenon as resulting from a bifurcation in phase space topology as the thermal trapping parameter $mv_{th}/eA_{\max }$ passes through unity. When $eA_{\max } \ll mv_{th}$ the primary phase space vortex structures are antisymmetric across the $v_y$-plane, producing a perturbation in current density without a perturbation in charge density. Once $eA_{\max } \approx mv_{th}$ particle orbits with low $P_y$ become more symmetric across the $v_y$-plane and produce a coherent perturbation in the density.

Thus, both electric and magnetic field trapping are associated with local variations in the electron density which manifests as space charge. In the case of magnetic trapping the charge density is a higher-order effect, and the first-order effect is production of electric currents to sustain the magnetic mode. The development of space charge from magnetic pressure is anticipated in Morse & Nielson (Reference Morse and Nielson1971), and the longitudinal field is explained in Taggart et al. (Reference Taggart, Godrey, Rhoades and Ives1972) to arise at second order from the coupling of two magnetic modes. Dynamic space charge, or filamentation, has been observed in modelling to modify growth rates, in both early and more recent studies (Taggart et al. Reference Taggart, Godrey, Rhoades and Ives1972; Tzoufras et al. Reference Tzoufras, Ren, Tsung, Tonge, Mori, Fiore, Fonseca and Silva2006), pointing to the importance of higher-order effects prior to saturation. Since the fraction of trapped electrons is proportional to the magnetic energy, with saturation when the characteristic magnetic bounce frequency reaches the growth rate (Davidson et al. Reference Davidson, Hammer, Haber and Wagner1972), it follows that longitudinal electric field should trend nonlinearly with the magnetic field. Finally, the theory of Appendix D interprets the density perturbation as arising from particles with $v_y\approx 0$ and energies $H<{(eA_{\max })^2}/{2m}$.

4.5. Saturation of many unstable Weibel modes in two spatial dimensions

Weibel instability of a homogeneous unmagnetized plasma is inherently multidimensional as wavevectors oblique to the principal anisotropy axes have comparable growth rates to the principal axes (similar to the multidimensional streaming instability studied in § 2.7). We consider for example the branch of the dispersion function for the anisotropic Maxwellian given by (4.24) with transverse magnetic field out-of-plane and the longitudinal and transverse electric fields in-plane. The out-of-plane magnetic field allows a reduced 2D2V phase space geometry for tractable continuum-kinetic simulation (Skoutnev et al. Reference Skoutnev, Hakim, Juno and TenBarge2019). It should be kept in mind that the instability dynamics are truly three-dimensional just as in the multidimensional Langmuir turbulence simulation of § 2.7. Figure 22 plots the growth rate as a function of the wavevector $\boldsymbol {k}=k_x\hat {x} + k_y\hat {y}$ of the unstable branch of (4.24) for a zero-order bi-Maxwellian of anisotropy $A=3$ and with $v_t/c=0.3$. If the destabilized ordinary wave studied in § 4.4 were included in the simulation, the magnetic field would also take in-plane components, necessitating a third velocity dimension. We point out the same caveats for this simulation as for the two-dimensional two-stream simulation; the nonlinear phase space dynamics are likely similar to the three-dimensional problem, but the nonlinear dynamics of the saturated state (that is, the turbulence statistics) will likely be different. The essence of the problem lies in restricting the cylindrical symmetry of the near-maximal-growth-rate modes about the principal anisotropy axis to a plane. Fully three-dimensional simulations are necessary to clarify these issues.

Figure 22. Growth rates of multidimensional Weibel instability as $\textrm {Im}(\omega )/\omega _{pe}$ for an anisotropic bi-Maxwellian of anisotropy $A=3$ where the electric field is assumed to lie in the $(x,y)$-plane and the first-order transverse magnetic field to be out-of-plane. The higher-temperature direction is assumed to lie in the $\hat {y}$-direction assuming $v_t/c=0.3$. The red line identifies marginal stability.

4.5.1. Numerical methods for the two-dimensional multimode Weibel instability

The simulations are conducted with the methods of Appendix A and Crews & Shumlak (Reference Crews and Shumlak2022) with a few key differences. Specifically, Ampere's and Faraday's laws are used for the field equations which are time-integrated in Fourier spectral space with the same third-order Adams–Bashforth method as the kinetic equation. There is also an important difference in initialization of the unstable modes. In the Vlasov–Poisson system the field part of the kinetic eigenfunctions is described by the scalar potential, yet in the multidimensional Vlasov–Maxwell system the kinetic eigenfunctions consist of a vector self-consistent field-plasma configuration. Initialization with kinetic eigenfunctions in an electrodynamic problem necessitates solving the eigenvalue problem $\chi \boldsymbol {E} = \zeta ^2\boldsymbol {E}$ in the wavevector frame, as discussed in § 4.2. In our implementation, the eigenfunctions are computed in the rotated wavevector frame and then antirotated back into the $(x,y)$-plane with components $\boldsymbol {E} = E_x\hat {x} + E_y\hat {y}$. Thus, for each desired pair of wavenumbers $(k_x, k_y)$ a perturbation is applied as

(4.29)\begin{equation} f_1 = \mathcal{A}\text{Re}\left[\frac{q}{{\rm i}\omega}\left\{E_x\left(\frac{\partial f_0}{\partial u} + \frac{v_x(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_0)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\right) + E_y\left(\frac{\partial f_0}{\partial v} + \frac{v_y (\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}_v f_0)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\right) \right\} {\rm e}^{{\rm i}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x} + \tilde{\varphi})}\right] \end{equation}

with $\tilde {\varphi }$ a randomly chosen phase shift per wavevector and the amplitude $\mathcal {A}=2\times 10^{-3}$. The magnetic field is then initialized as $B_z = \zeta ^{-1}\hat {k}\times \boldsymbol {E}$ where $\zeta$ is the eigenvalue. A higher amplitude than the one-dimensional problem is chosen to reduce time to saturation.

4.5.3. Nonlinear simulations of the two-dimensional Weibel instability

The problem described in § 4.5.2 was simulated using the methods of § 4.5.1. Figure 23 shows the energy traces of the electrodynamic field during the linear phase and beyond instability saturation. Nonlinear saturation occurs around $t\omega _p=17$, as gauged by the $\hat {y}$-directed transverse electric field. The $\hat {x}$-directed electric field energy follows a similar trend as the one-dimensional problem, its growth composed of two effects: to first-order from the initialized oblique modes, and to second-order in the magnetic field as discussed in § 4.4. At nonlinear saturation, magnetic-trapping electron holes form between the counter-streaming mean flows.

Figure 23. Evolution of domain-integrated energy for (a) the out-of-plane magnetic field $B_z$, (b) the electric field transverse to the maximum growth-rate axis and (c) the electric field parallel to that axis (formerly the longitudinal field of the one-dimensional simulation). A difference from the one-dimensional simulation is a steadily decreasing transverse electric energy rather than a coherent oscillation after saturation time. The $\hat {x}$-directed electric energy increases at a rate faster than the growth of the linear modes due to the same space charge effects as in the single-mode simulation discussed in § 4.4.

Figure 24 illustrates the dynamics of the multidimensional instability through streamlines of the current density $\boldsymbol {j}(x,y)$ and filled contours of its magnitude for three times in the simulation output. Trigonometric interpolation is used to visualize the current density by zero-padding the spectrum and inverse Fourier transformation, as the spectrum is properly resolved up to the chosen spectral cutoff. In figure 24(b) one can observe spiral current streamlines indicating the local production of space charge as the electron density filaments into two-dimensional analogues of the magnetic-trapping electron holes studied in the single-mode problem, also observed in Califano et al. (Reference Califano, Pegoraro and Bulanov1997). By the simulation's end, nonlinear mixing has caused some of the filaments to rotate, as in the multidimensional electrostatic problem considered in § 2.7, indicating a similar isotropization on averaged scales. Inspection of figure 25, which plots the spatially averaged distribution function $\langle f\rangle _{(x,y)}(v_x, v_y)$ at $t\omega _p=0$ and $t\omega _p=40$ in figure 25(a,b) and in figure 25(c) the trace of the anisotropy parameter $A$, shows the anisotropic Maxwellian to relax towards isotropy through Weibel-induced turbulence. The initial anisotropy $A=3$ is observed to decrease monotonically until at saturation $A\approx 1$, meaning that the saturated state is weakly anisotropic. This persistent anisotropy coexists with the fluctuating filamentation currents, consistent with the observations of Davidson et al. (Reference Davidson, Hammer, Haber and Wagner1972). Indeed, persistent anisotropy accompanying sheared flows in the saturated state is expected for collisionless dynamics (Del Sarto & Pegoraro Reference Del Sarto and Pegoraro2017).

Figure 24. Evolution of current at times (a) $t\omega _p=0$, (b) $t\omega _p = 25$ and (c) $t\omega _p = 39.5$. Plotted are streamlines of $\boldsymbol {j}$ and its magnitude $|\boldsymbol {j}|$ as filled contours. Current tends to form closed paths at long wavelength in $y$. The indicated spiral vortices in the streamlines demonstrate rapid local space charge production ($\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {j}\neq 0$) as expected from the filamentation and trapping dynamics of the one-dimensional problem. The large circulating configuration space electron current vortices manifest on small scales as counter-propagating electron streams which sustain a train of magnetic-trapping electron phase space vortices. In this unmagnetized model problem, these structures isotropize from four-dimensional phase space dynamics.

Figure 25. The relaxing spatially averaged distribution function $\langle f\rangle _{(x,y)}(v_x, v_y)$ shown at two times: (a) the initial condition $t\omega _p=0$, (b) the stop-time $t\omega _p=40$ and in addition (c) the time evolution of the anisotropy parameter $A = \langle v_y^2\rangle / \langle v_x^2\rangle - 1$. The turbulence of magnetic trapping phase space vortices in the turbulent currents of the saturated instability isotropizes the distribution function. Note that $A=0$ would indicate isotropy, so that the saturated state consists of persistent anisotropy coexisting with the wave field (Davidson et al. Reference Davidson, Hammer, Haber and Wagner1972).