# On the validity of quasilinear theory applied to the electron bump-on-tail
instability

D.W. Crews  U. Shumlak  Computational Plasma Dynamics Lab, Aerospace and
Energetics Research Program,
University of Washington, Seattle, Washington 98195-2400 USA.

###### Abstract

The accuracy of quasilinear theory applied to the electron bump-on-tail
instability, a classic model problem, is explored with conservative high-order
discontinuous Galerkin methods applied to both the quasilinear equations and
to a direct simulation of the Vlasov-Poisson equations. The initial condition
is chosen in the regime of beam parameters for which quasilinear theory should
be applicable. Quasilinear diffusion is initially in good agreement with the
direct simulation but later underestimates the turbulent momentum flux. The
direct simulation corrects from quasilinear evolution by quenching the
instability in a finite time and producing a robust state of oscillation. Flux
enhancement above quasilinear levels occurs as the phase space eddy turnover
time in the largest amplitude wavepackets becomes comparable to the transit
time of resonant phase fluid through wavepacket potentials. In this regime
eddies effectively turn over during wavepacket transit so that phase fluid
predominantly disperses by eddy phase mixing rather than by randomly phased
waves. The enhanced turbulent flux of resonant phase fluid leads in turn,
through energy conservation, to an increase in non-resonant turbulent flux and
thus to an enhanced heating of the main thermal body above quasilinear
predictions. These findings shed light on the kinetic turbulence fluctuation
spectrum and support the theory that collisionless momentum diffusion beyond
the quasilinear approximation can be understood through the dynamics of phase
space eddies (or clumps and granulations).

##  I Introduction

One-dimensional electrostatic turbulence is among the oldest model problems in
plasma physics and stimulated the development of the quasilinear theory of
plasma diffusion [1, 2, 3]. Quasilinear diffusion theories are used to predict
the statistical properties of turbulent energy and the plasma distribution
function in the saturated state of high-dimensional dynamic processes such as
velocity space instability with a reduced low-dimensional model. Quasilinear
methods have been widely applied to model turbulent diffusion in both physical
and velocity space (e.g. gyrokinetic modeling [4], rotating Couette flow [5],
lower-hybrid drift instability [6], firehose instability [7], electron
interaction with whistlers [8], etc.), often reproducing experimental
phenomena beyond the theory’s formal region of validity [9].

Quasilinear theory is concerned with the dynamics of macroscopic observables
such as spatially-averaged quantities and the fluctuation power spectrum. It
is applicable for weak turbulence, i.e. a broad spectrum of small-amplitude
non-self-correlated linear waves in a homogeneous medium [10, 11]. In order
that phase interference, or wave beating, not lead to nonlinearity it is also
necessary to suppose random phasing of the waves [12]. However the random
phase approximation is somewhat fictional because the spectral width of
unstable wavenumbers is finite. Interference of finite spectral width produces
a spatial distribution of wave energy (wavepackets) even for random phases.
For example the two-species numerical study in [13] saw the late-time
development of Langmuir cavitons from an initial ensemble of randomly-phased
linear waves. The high amplitudes of wavepacket potentials lead to deviations
from quasilinear predictions as linear waves grow from instability [14].

Nevertheless quasilinear theory makes accurate predictions of the saturated
state of instability despite the breakdown of the random phase approximation.
This article demonstrates by a direct simulation of an electron bump-on-tail
instability for which quasilinear theory should be applicable that weak
nonlinearity due to spatial structure in the distribution of wave energy
enhances the rate of turbulent flux above quasilinear levels. However the
effect of this discrepancy is to reach a similar saturated state as the
quasilinear prediction, but in a finite time rather than asymptotically with
corresponding corrections to the transient linear growth rates and with
significant fluctuations in the probability distribution function atop the
quasilinear equilibrium. The success of quasilinear theory is in prediction of
the greatest amplitude part of the field spectrum as this portion is due to
linear growth. On the other hand, the quasilinear prediction does not capture
some aspects of the weak kinetic turbulence spectrum such as plasma wave
harmonics and the eddy turbulence of resonant particles, i.e. phase space
granulations.

Many numerical studies have been done comparing quasilinear analysis to
solutions of the Vlasov equation using particle-in-cell (PIC) method [15, 16,
17]. Yet it is well-known that PIC methods are prone to errors due to
statistically sampling the sensitive trajectories of the continuous
distribution [18, 19, 20], so it is worthwhile to explore alternative kinetic
modeling methods. In one such alternative the equation for the continuous
distribution is solved by a demonstrably convergent and conservative finite
element discretization of phase space, using for example discontinuous
Galerkin methods [21, 22, 23]. This work uses a high-order discontinuous
Galerkin method to sufficiently resolve the detailed phase flow up to and past
saturation of the kinetic instability. From this point of view the turbulent
nature of resonant electrons, and the distinct behavior of resonant and non-
resonant electrons, can be clearly seen.

The structure of the article is a review of the theory of the bump-on-tail
instability with an emphasis on the aspects necessary for a numerical study,
followed by novel methods of discretization, and finally an analysis of the
simulation results. Linear electrostatic theory is reviewed in Section II in
order to identify the kinetic eigenmodes used for the perturbation of the
distribution function, followed by a review of the quasilinear theory of the
bump-on-tail in Section III with comments on the role of resonant and non-
resonant parts of the distribution function and on the hierarchy of timescales
required for validity of the approximation. The details of the discretization
methods are discussed in Section IV. The simulation results are then explored
in Section V by an analysis of phase space structures in the Vlasov-Poisson
simulation and a comparison of quasilinear and Vlasov-Poisson levels of
turbulent flux.

##  II Linear theory of bump-on-tail instability

This section reviews electrostatic plasma instability with an emphasis on the
linear response associated with unstable modes. One-dimensional collisionless
electron dynamics in a uniform neutralizing background are governed by the
Vlasov-Poisson equations [24],

$\displaystyle\frac{\partial f}{\partial t}$ $\displaystyle+v\frac{\partial
f}{\partial x}-E(x)\frac{\partial f}{\partial v}=0,$ (1)
$\displaystyle\frac{dE}{dx}$
$\displaystyle=1-\int_{-\infty}^{\infty}f(x,v,t)dv,$ (2)

which describes the evolution of the probability distribution function
$f(x,v,t)$ through phase space as influenced by the self-consistent electric
field $E(x)$. Velocities are normalized to the thermal velocity $v_{t}$,
lengths normalized to the Debye screening length $\lambda_{D}$, and electric
field normalized to the thermal energy per Debye length. Static ions and
collisionless trajectories are valid for timescales $\tau$ less than the ion
plasma period and the electron-electron collision time, i.e.
$\tau\ll\omega_{pi}^{-1}$ and
$\tau\ll\frac{\Lambda}{\log\Lambda}\omega_{pe}^{-1}$ where $\Lambda\gg 1$ is
the inverse plasma parameter. This article simulates the Vlasov-Poisson
system, and the result will be referred to as either the Vlasov-Poisson
simulation or just the Vlasov simulation.

###  II.1 General linear theory of electrostatic instability

The linearized Vlasov-Poisson equations about a homogeneous equilibrium
$f_{0}(v)$ are

$\displaystyle\frac{\partial f_{1}}{\partial t}$
$\displaystyle+v\frac{\partial f_{1}}{\partial x}+\nabla\varphi\frac{\partial
f_{0}}{\partial v}=0,$ (3) $\displaystyle\nabla^{2}\varphi$
$\displaystyle=\int_{-\infty}^{\infty}f_{1}(x,v,t)dv$ (4)

with $\varphi$ the potential. Following Landau [25], spatial Fourier transform
and one-sided temporal Fourier transform (or Laplace transform with
$s\to-i\omega$) leads, if the initial condition is denoted $g(k,v)\equiv
f_{1}(k,v,t=0)$, to a solution for electric potential,

$\displaystyle g^{*}(k,\zeta)$
$\displaystyle\equiv\int_{\mathcal{C}}\frac{g(k,v)}{\zeta-v}dv,$ (5)
$\displaystyle\varphi(k,\zeta)$
$\displaystyle=-\frac{i}{k^{3}}\frac{g^{*}(k,\zeta)}{\varepsilon(k,\zeta)}$
(6)

with $\zeta=\frac{\omega}{k}$ the complex phase velocity where the complex
dielectric function is

$\varepsilon(k,\zeta)=1+\frac{1}{k^{2}}\int_{\mathcal{C}}\frac{1}{\zeta-v}\frac{\partial
f_{0}}{\partial v}dv$ (7)

and $\mathcal{C}$ the Landau contour. Inversion of Eq. 6 by the residue
theorem suggests that the response to a general perturbation $g(k,\zeta)$ will
include all solutions of $\varepsilon(k,\zeta)=0$. However initial conditions
of the linear response form, with a single complex pole at an unstable
frequency, excite only that mode. That is, suppose the initial condition is
chosen as

$g(k,v)=\frac{1}{k^{2}}\frac{1}{\zeta_{n}-v}\frac{\partial f_{0}}{\partial v}$
(8)

where $\text{Im}(\zeta_{n})>0$ solves the dispersion relation in the upper
half-plane, i.e. $\varepsilon(k,\zeta_{n})=0$. Then its transform along the
Landau contour is the dielectric function with a pole at $\zeta=\zeta_{n}$,

$g^{*}(k,\zeta)=-\frac{\varepsilon(k,\zeta)}{\zeta-\zeta_{n}}.$ (9)

The potential response is
$\varphi(k,\zeta)=\frac{i}{k^{3}}\frac{1}{\zeta-\zeta_{n}}$, inverting to
$\varphi(k,t)=\frac{1}{k^{3}}e^{-i\omega_{n}t}$. An initial condition in the
form of Eq. 8 represents the discrete part of the spectrum of the linearized
Vlasov-Poisson integral operator, as found by Case and van Kampen [26, 27].
Equation 8 is used in this work as the perturbation to excite the bump-on-tail
instability.

###  II.2 Initialization of many discrete linear modes

An unstable distribution which should meet the applicability conditions of
quasilinear theory consists of a hot drifting Maxwellian through a main
thermal population

$\begin{split}f_{0}(v)=\frac{1}{(1+\chi)\sqrt{2\pi}}\Big{\\{}\frac{1}{v_{t0}}&\exp\Big{(}-\frac{v^{2}}{2v_{t0}^{2}}\Big{)}+\\\
&\frac{\chi}{v_{tb}}\exp\Big{(}-\frac{(v-v_{b})^{2}}{2v_{tb}^{2}}\Big{)}\Big{\\}}.\end{split}$
(10)

The bump-on-tail instability is the textbook example of a distribution with
quasilinear evolution. According to Ref. [28], reasoning by the Bohm-Gross
dispersion relation, instability development should be in the weak turbulence
regime given a bump fraction $\chi=0.05$, bump spread
$v_{tb}=\chi^{1/3}v_{b}$, and beam velocity $v_{b}=5v_{t0}$. Velocities are
then normalized to $v_{t0}=1$. Figure 1 shows the unstable branch of
$\varepsilon(k,\zeta)=0$ corresponding to the coupling of right-going
oscillations of the main body and left-going oscillations on the drifting
beam.

Figure 1: Bump-on-tail dispersion relation showing: a) wave velocities; and b)
frequencies. The bump-on-tail distribution has an unstable branch of solutions
with phase velocities bounded by the regions of $\frac{\partial
f_{0}}{\partial v}>0$ in the tail, in this case for $v\in(3,5)v_{t}$. The
unstable solution is an electron-acoustic wave at the beam velocity for long
wavelengths and an acoustic wave propagating at the bulk thermal velocity for
short wavelengths. This unstable coupled acoustic wave is damped at both
extremes and has maximum instability growth where the frequency
$\omega\approx\omega_{p}$.

The Vlasov-Poisson simulation is initialized as a sum over the eigenmodes,

$g(x,v)=f_{0}(v)+\sum_{n}\frac{\alpha_{n}}{k_{n}}\text{Re}\Big{(}\frac{1}{\zeta_{n}-v}e^{i(k_{n}x+\theta_{n})}\Big{)}\frac{\partial
f_{0}}{\partial v}$ (11)

where each $\zeta_{n}$ with $\text{Im}(\zeta_{n})>0$ is the solution to
$\varepsilon(k_{n},\zeta_{n})=0$, the factor $\theta_{n}$ is a random phase
shift, and $\alpha_{n}$ is the field amplitude. The factor $k_{n}^{-1}$ sets
each initialized mode’s field energy to $|E_{n}|^{2}=\alpha_{n}^{2}$. Figure 2
shows the perturbation on a domain of $L=5000\lambda_{D}$ initialized in this
manner. The initial condition is chosen in this way to avoid exciting the
other damped branches of $\varepsilon(k,\zeta)=0$, as the quasilinear theory
studies the unstable wave branch.

Figure 2: a) The perturbation is shown in phase space representation as a
subset of the full domain, which has $5000\lambda_{D}$ length and spans
velocities $v\in(-25,25)$. The main thermal body ($v\in(-3,3)$) supports
coherent plasma oscillations while the modes in the resonant band
($v\in(3,5)$) appear randomly phased. The unstable modes are rightward-
propagating so the non-resonant distribution of positive velocity
$(v\in(0,3))$ oscillates with greater energy than the negative velocity part.
Wavepackets in the initial condition are seen as regions of greater amplitude.
b) The Fourier amplitudes as $\log(1+|f|)$ of the perturbation. Spectral
energy is peaked along a coherent crescent shape bridging the band of unstable
velocities.

##  III Quasilinear theory of the bump-on-tail

The quasilinear approximation is the first approximation to velocity-space
diffusion in electrostatic turbulence and is a reduction from the Vlasov-
Poisson theory [1, 2, 11, 24, 28]. It is mathematically similar to theories of
Reynolds-averaging and eddy-viscosity in neutral fluid turbulence [29].
Consider a spatially-periodic domain of length $L$. Let $\langle
f\rangle_{L}\equiv L^{-1}\int_{0}^{L}fdx$ be the spatially-averaged
distribution and $\delta f\equiv f-\langle f\rangle_{L}$ the fluctuation.
Equation 1 is then spatially averaged and the resulting mean equation is
subtracted out to obtain an equation for the fluctuation. The two equations
for the mean and the fluctuation are then

$\displaystyle\frac{\partial\langle f\rangle_{L}}{\partial t}$
$\displaystyle=\frac{\partial}{\partial v}\langle E\delta f\rangle_{L}$ (12)
$\displaystyle\frac{\partial(\delta f)}{\partial t}$
$\displaystyle+v\frac{\partial(\delta f)}{\partial x}=\frac{\partial}{\partial
v}\Big{(}E\langle f\rangle_{L}+E\delta f-\langle E\delta f\rangle_{L}\Big{)}.$
(13)

Note that $E=\delta E$ for internal fields. The momentum flux of the
fluctuation $\delta f$ is composed of two terms: the flux from the mean
distribution $E\langle f\rangle_{L}$ and the fluctuation in the fluctuating
flux $E\delta f-\langle E\delta f\rangle_{L}$. The quasilinear closure assumes

$E\delta f-\langle E\delta f\rangle_{L}\ll E\langle f\rangle_{L},$ (14)

transforming Eq. 13 into a quasilinear PDE (one linear in its highest
derivatives). The quasilinear equations are

$\displaystyle\frac{\partial\langle f\rangle_{L}}{\partial t}$
$\displaystyle=\frac{\partial}{\partial v}\Big{\langle}E\delta
f\Big{\rangle}_{L}$ (15) $\displaystyle\frac{d(\delta f)}{dt}$
$\displaystyle=\frac{\partial}{\partial v}\Big{(}E\langle f\rangle_{L}\Big{)}$
(16)

with $\frac{d}{dt}=\partial_{t}+v\partial_{x}$ the change along a zero-order
trajectory. The momentum flux of the background distribution $\langle E\delta
f\rangle_{L}$ is statistically the un-normalized field-particle correlation
coefficient [30] and is a tool increasingly used independent of quasilinear
theories in the analysis of energy transfer in collisionless plasmas [31, 32,
33, 34]. The quantity $\langle E\delta f\rangle_{L}$ also has a basic
interpretation as the turbulent momentum flux, in analogy to a similar
quantity arising in the averaged equation for advection of a passive scalar in
a fluctuating velocity field [29]. These interpretations of $\langle E\delta
f\rangle_{L}$ as a correlation coefficient and as a mean turbulent flux are
independent of any closure assumptions.

The quasilinear closure amounts to assuming that the field-particle
correlation has no variance, or in other words negligible structure inside the
averaging box. This assumption is often stated as the random phase
approximation as there are no wavepackets, or structure, in white noise. The
fluctuating distribution and field are expanded in Fourier series to solve the
linear equation for $\delta f$ as

$\delta f_{n}=\frac{iE_{n}}{\omega(k_{n})-k_{n}v}\frac{\partial\langle
f\rangle_{L}}{\partial v}.$ (17)

The covariance $\langle E\delta f\rangle_{L}$ is converted to a spectral sum
by Plancheral’s theorem,

$\langle E\delta f\rangle_{L}=\frac{1}{L}\int_{0}^{L}E\delta
fdx=\sum_{n=-\infty}^{\infty}\delta f_{n}E_{n}^{*}.$ (18)

Substituting Eq. 17 into Eq. 18 gives the turbulent momentum flux in the
quasilinear closure,

$\langle E\delta
f\rangle_{L}=\Big{(}\sum_{n=-\infty}^{\infty}\frac{i}{\omega(k_{n})-k_{n}v}|E_{n}|^{2}\Big{)}\frac{\partial\langle
f\rangle_{L}}{\partial v}.$ (19)

As the left-hand side is the mean momentum flux, and the right-hand side
proportional to the mean gradient, Eq. 19 defines a diffusivity as $\langle
E\delta f\rangle_{L}=D(v)\partial_{v}\langle f\rangle_{L}$. In the sense that
the quasilinear closure sets the turbulent flux as proportional to the mean
gradient this specification of diffusivity is closely related to the eddy
viscosity theory in physical space dating to Boussinesq [35]. Now, the
simplest equation for the evolution of the spectrum $\mathcal{E}(k)$ is
exponential growth (assuming no modal interactions and slow modulations),
where the dispersion relation $\omega(k)=\omega_{r}(k)+i\omega_{i}(k)$ is
taken as the solution to $\varepsilon(\omega,k)=0$. From all this one arrives
at the electrostatic quasilinear diffusion equations [1, 2, 11, 28, 36],

$\displaystyle\frac{\partial\langle f\rangle_{L}}{\partial t}$
$\displaystyle=\frac{\partial}{\partial v}\Big{(}D(v)\frac{\partial\langle
f\rangle_{L}}{\partial v}\Big{)},$ (20) $\displaystyle D(v)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{2|\omega_{i}|}{(\omega_{r}-k_{n}v)^{2}+\omega_{i}^{2}}|E_{n}|^{2},$
(21) $\displaystyle\frac{d|E_{n}|^{2}}{dt}$
$\displaystyle=2\omega_{i}|E_{n}|^{2},\quad\varepsilon(\omega,k_{n})=0,$ (22)

with the absolute value $|\omega_{i}|$ following a one-sided time analysis,
i.e. causality [36]. The equations are intrinsically nonlinear with a self-
consistent determination of the dispersion relation key to the system’s energy
and momentum conservation properties [24]. The analytic solution of
$\varepsilon(\omega,k)=0$ as a complex root can be numerically expensive as
$\langle f\rangle_{L}$ evolves from a sum of Maxwellians to a general form.
The dispersion relation can be solved approximately on the real line using the
well-known small growth-rate approximation, valid for
$\omega_{i}\ll\omega_{r}$ [28],

$\displaystyle\varepsilon_{r}(\omega_{r},k)$ $\displaystyle\equiv
1+\frac{1}{k^{2}}\mathcal{P}\int_{-\infty}^{\infty}\frac{1}{\zeta_{r}-v}\frac{\partial\langle
f\rangle_{L}}{\partial v}dv=0,$ (23) $\displaystyle\omega_{i}$
$\displaystyle=\pi\frac{\partial\langle f\rangle_{L}}{\partial
v}\Big{(}\frac{\partial\varepsilon_{r}}{\partial\omega_{r}}\Big{|}_{\omega=\omega_{r}}\Big{)}^{-1}$
(24)

where $\zeta_{r}=\omega_{r}/k$ and $\mathcal{P}$ is the Cauchy principal value
(P.V.), calculable by Hilbert transform.

###  III.1 The continuous spectrum limit

By defining a spectral density $\mathcal{E}(k)=k_{0}^{-1}|E(k_{n})|^{2}$ with
$k_{0}=\frac{2\pi}{L}$ and considering the limit as $L\to\infty$, the
quasilinear equations can be modeled with a continuous energy spectrum [28]

$\displaystyle\frac{\partial f}{\partial t}$
$\displaystyle=\frac{\partial}{\partial v}\Big{(}D(v)\frac{\partial
f}{\partial v}\Big{)},$ (25) $\displaystyle D(v)$
$\displaystyle=\int_{0}^{\infty}\frac{2|\omega_{i}|}{(\omega_{r}-kv)^{2}+\omega_{i}^{2}}\mathcal{E}(k)dk,$
(26) $\displaystyle\frac{d\mathcal{E}(k)}{dt}$
$\displaystyle=2\omega_{i}\mathcal{E}(k),\quad\varepsilon(\omega,k)=0.$ (27)

The continuous-spectrum approximation is appropriate for large domains and for
analysis of the resonant and non-resonant contributions to the diffusivity. A
numerical advantage of using a continuous spectrum is that a high-order
polynomial representation may be used for the spectral density
$\mathcal{E}(k)$. The solution of the dielectric function on the spectral
points $k_{n}$ is typically the most expensive part of a numerical solution.
With a high-order representation fewer interpolation nodes are required than
the evenly-spaced lattice frequencies of the finite-interval problem. Fewer
evaluations of $\varepsilon(\omega,k)$ speed the simulations considerably.

####  III.1.1 Resonant and non-resonant parts of the diffusivity

It is common to split the diffusion coefficient into resonant and non-resonant
parts in the limit of slow modulation $\omega_{i}\to 0$ [11, 28, 12]. Applying
the Plemelj relation to the linear response

$f_{1,n}\to\mathcal{P}\frac{i}{\omega_{r}-k_{n}v}E_{n}\frac{\partial
f_{0}}{\partial v}+\pi\delta(\omega-k_{n}v)E_{n}\frac{\partial f_{0}}{\partial
v}$ (28)

suggests it has two principal terms. The first term represents the non-
resonant contribution, the response of the main body thermal plasma particles
in sustaining the wave motion (seen in Fig. 2(a) as the checkerboard pattern
for $v\in(-3,3)$), while the second term accounts for resonant energy transfer
and is the primary contribution to the integration [36, 37]. Having
consistently taken the small growth-rate limits, the quasilinear system is
usually written as

$\displaystyle\frac{\partial f}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
v}\Big{(}(D^{r}(v)+D^{nr}(v))\frac{\partial f}{\partial v}\Big{)},$ (29)
$\displaystyle\frac{d\mathcal{E}}{dt}$ $\displaystyle=2\omega_{i}\mathcal{E},$
(30) $\displaystyle D^{r}(v)$
$\displaystyle=\pi\int_{-\infty}^{\infty}\delta(\omega_{r}-kv)\mathcal{E}(k)dk,$
(31) $\displaystyle D^{nr}(v)$
$\displaystyle=``\mathcal{P}"\int_{-\infty}^{\infty}\frac{\omega_{i}}{(\omega_{r}-kv)^{2}}\mathcal{E}(k)dk$
(32)

with $``\mathcal{P}"$ the principal value operator. Quotes are used because
the principal value operation on the denominator, a term quadratic in a
nonlinear function $\omega_{r}(k)$, is not necessarily well-defined. While it
may be tempting to discard the non-resonant diffusivity and evolve only the
resonant distribution as it is smaller than the resonant diffusivity by
approximately $\omega_{i}/\omega_{r}$, it was well-phrased by Kadomtsev that
the population of nonresonant particles is greater by this same factor and
$D^{nr}$ cannot be neglected [36]. The difficulty in a numerical
implementation of the resonantly-split quasilinear equations lies in
approximating the non-resonant diffusivity. For example, a typical
approximation for $D^{nr}$ is constant in velocity [11, 28]. For the problem
under consideration there is no need to make this approximation. The
quasilinear problem solved here numerically is the unreduced system of Eqs. 20
and 21.

####  III.1.2 On analytic prediction of the saturated spectrum

Note that as the growth-rate $\omega_{i}(k_{n})\to 0$ the diffusivity
corresponding to mode $k_{n}$ becomes singular, i.e.
$D(v=\zeta_{r}(k_{n}))=\frac{2}{|\omega_{i}|}|E_{n}|^{2}\to\infty$. As
$\omega_{i}\propto\frac{\partial f_{0}}{\partial v}$ the dynamic equation
approaches a singular diffusion equation in the region of resonant velocities,
similar in form to

$\frac{\partial f}{\partial
t}\to\pi^{-1}\sum_{k_{n}}\partial_{v}\Big{(}\varepsilon_{r}^{\prime}|E_{k_{n}}|^{2}\frac{\partial_{v}f}{|\partial_{v}f|_{v=v_{\phi,{k_{n}}}}}\Big{)}.$
(33)

Such singular diffusion equations instantaneously flatten the diffusing
variable [38], but the singular problem is approached asymptotically in QL
theory. The singular behavior of the turbulent diffusivity does not require
having taken the small-growth rate limit $\omega_{i}\to 0$ or having
considered a resonant/non-resonant splitting of the diffusivity. From a
numerical standpoint the transition from a regular diffusion equation with
smooth diffusivity to a near-singular equation with delta-valued diffusivity
peaks poses a serious challenge for a numerical time-dependent solution of the
QL equations. Yet this is not necessarily a problem for analysis.

A theory is sometimes discussed which analytically predicts the asymptotic
state assuming dispersionless waves $\varepsilon=1-\omega_{p}^{2}/\omega^{2}$
and only resonant interaction, i.e. $D^{nr}=0$ [36]. From this one can show
that, with $v_{1}$ the lower limit of the asymptotically flat region,

$|E_{k}|^{2}(t_{\infty})=\zeta^{3}\int_{v_{1}}^{\zeta}f(t_{\infty},v^{\prime})-f(t_{0},v^{\prime}))dv^{\prime}.$
(34)

The analysis leading to Eq. 34 must be cautioned as approximate because the
neglect of non-resonant diffusivity means that the growth rate at the edges of
the diffusing region is significantly overestimated, as illustrated in Fig. 3
for the considered bump-on-tail problem. Sharp gradients develop in the
distribution near the “notch” and “hill” of the bump as seen in Fig. 3(a)
which would be otherwise filled in by non-resonant diffusion. These sharp
gradients correspond to near-singular instability growth rates, as shown in
Fig. 3(b).

Figure 3: Solution of the quasilinear system using only the resonant diffusion
coefficient and assuming $\omega=\omega_{p}$, showing: a) flattening of the
distribution function in the region of $\partial_{v}f>0$ and b) modal growth
rates. Use of only resonant diffusion results in near-singular growth rates at
the edges of the diffusing region as sharp gradients develop in the
distribution. This means the growth rate is greatly overestimated, and the
saturated spectrum is not correctly predicted by Eq. 34 across the full
breadth of wavenumbers, specifically for non-initially-resonant phase
velocities.

###  III.2 Conditions for validity of the QL approximation

A considerable amount of research has been conducted on the conditions under
which QL theory should be valid [39, 40, 41], aptly summarized in Ref. [11]
and briefly repeated here. The validity conditions center around whether the
condition given by Eq. 14 is justified, so that the fluctuation $\delta f$
should satisfy a quasilinear equation. As noted, the condition for a small
divergence of the fluctuating flux is unimportance of the wavefield’s spatial
structure within the averaging box. Since some distribution of wavepackets
from spectral superposition is inevitable, the condition is generally
rephrased as to whether particles undergo random walks from the wavepacket
potentials [11]. The quantities of importance for this are the electron bounce
frequency $\omega_{b}\approx\sqrt{\rho_{c}}$ with $\rho_{c}$ the normalized
charge density, the wavepacket autocorrelation frequency
$\tau_{ac}^{-1}\equiv|(\zeta-v_{g})\Delta k|$ describing the timescale on
which phase waves remain within a wavepacket, and the distribution relaxation
frequency $\tau_{r}^{-1}\equiv f_{0}^{-1}\partial_{t}f_{0}$.

With these quantities, the conditions for validity are that i) the bounce
frequency be less than the autocorrelation frequency
$\omega_{b}\ll\tau_{ac}^{-1}$, and ii) the growth rates lie between the
relaxation and autocorrelation frequencies
$\tau_{r}^{-1}<\omega_{i}<\tau_{ac}^{-1}$. Since the growth rates
$\omega_{i}\to 0$, this second condition suggests that the QL theory should
saturate asymptotically with $\tau_{r}\to\infty$. Instability saturation in a
finite time violates this condition.

##  IV Numerical simulation methods

The 1D1V Vlasov-Poisson system is a two-dimensional hyperbolic-elliptic
system, while the QL theory is a 0D1V diffusion equation with nonlinear time-
dependent diffusivity. Thus each model requires its own method of numerical
solution. This section describes the high-order discontinuous Galerkin methods
[42, 43, 44, 45] developed and applied to each of the two systems.

###  IV.1 Discretization of the quasilinear equations

The QL diffusion equation is discretized in velocity using the local
discontinuous Galerkin (LDG) method as high resolution and accuracy is needed
in the region of resonant velocities when the diffusivity becomes highly
peaked. A small modification in the LDG method is made to account for the
nonlinearity of the system appearing as, defining $q\equiv\partial_{v}f$, the
flux $F=D(v)q(v)$ being a product of two interpolated functions. The
modification from the usual linear LDG method is that the Galerkin matrix
elements on the right-hand side of the semi-discrete equation are computed by
an exact quadrature of the products of the interpolation polynomial basis for
both $D(v)$ and $q(v)$.

The grid nodes of high-order DG methods are not evenly spaced, so the Hilbert
transform used to compute the principal-value integral in the small-growth
rate approximation for $\varepsilon(k,\zeta)$ cannot be done using FFTs or
most other methods discussed in the literature. Instead the principal value
integral is determined by a quadrature approximation of the velocity Fourier
coefficients to form the analytic signal, and the P.V. integral extracted as
its imaginary part. This is not an optimal approximation but it works
sufficiently well on the existing grid provided that the distribution is
represented to high-enough velocity in the tail regions to be effectively
zero. For this reason velocity-space is resolved from $v=(-20,30)$ with one
hundred finite elements. Each element contains a ninth-order polynomial basis
on Lobatto quadrature node interpolants. Only twenty elements are used in the
tails, with eighty elements spanning $v\in(2,10)$. A Fourier-by-quadrature
evaluation of the Hilbert transform is not computationally time-limiting as
the P.V. integral needs to be evaluated only once per timestep. Rather, the
nonlinear solution to $\varepsilon(k,\zeta)=0$ for each lattice frequency $k$
is typically the computational time-limiting factor.

With the P.V. integral determined, the dielectric function
$\varepsilon(k,\zeta)$ is evaluated as a function of arbitrary phase velocity
$\zeta$ via the Vandermonde matrix of the quadrature grid’s interpolation
polynomials. After solution with a nonlinear solver to obtain the phase
velocity root to $\varepsilon(k,\zeta)=0$ corresponding to the unstable wave
branch, the gradient $\varepsilon^{\prime}(k,\zeta)$ on the interpolated point
is found using the local derivative matrix. Further computing the interpolated
gradient of the distribution function, this determines the approximate linear
growth rates $\omega_{i}(k)$ by Eq. 24. Since the velocity space is highly
resolved the discretization is limited in accuracy by the
$\mathcal{O}\big{(}\frac{\omega_{i}}{\omega_{r}}\big{)}$ expansion of
$\varepsilon$.

As the distribution flattens and $\omega_{i}\to 0$ the diffusivity becomes
highly peaked. Use of an implicit time integration scheme is necessary to
avoid restrictive CFL-limited timestepping and avoid dispersive errors, as
such oscillations feedback into the growth rate calculation. A second-order
implicit midpoint method is used to evolve the diffusion equation. As the
diffusivity nonlinearly depends on the distribution, a nonlinear global solve
is avoided by taking an approximate explicit half-step in wave energy to the
midpoint with a sufficiently small time-step. The diffusion equation is then
solved by the implicit midpoint method. That is, with the semi-discrete QL
system as

$\displaystyle\frac{df}{dt}$ $\displaystyle=A(E)f$ (35)
$\displaystyle\frac{dE}{dt}$ $\displaystyle=2\omega_{i}(\partial_{v}f)E$ (36)

where $A$ represents the discretized diffusion operator and $E=|E_{k}|^{2}$
the spectral energy,

$\displaystyle E_{n+1/2}$
$\displaystyle=E_{n}\exp(h\omega_{i}(\partial_{v}f_{n})),$ (37) $\displaystyle
A_{n+1/2}$ $\displaystyle\equiv A(E_{n+1/2}),$ (38) $\displaystyle f_{n+1}$
$\displaystyle=(I-\frac{h}{2}A_{n+1/2})^{-1}(I+\frac{h}{2}A_{n+1/2})f_{n}$
(39)

where $h$ is the time-step. The calculations use $h=0.1\omega_{pe}^{-1}$ for
the time-advance.

###  IV.2 Discretization of the Vlasov-Poisson system

The Vlasov-Poisson system is solved using a mixed Fourier/DG spectral method.
The spatial coordinate is represented on a grid of evenly spaced nodes and the
variables $f(x,v,t)$, $E(x)$ are projected onto a truncated Fourier series
[46]. The Fourier coefficients each satisfy

$\displaystyle\frac{\partial f_{n}}{\partial t}$
$\displaystyle+ik_{n}vf_{n}+\partial_{v}F_{n}=0,$ (40) $\displaystyle
F_{n}(v,t)$ $\displaystyle=-(E*f)_{n},$ (41) $\displaystyle E_{n}$
$\displaystyle=ik_{n}^{-1}\int_{-\infty}^{\infty}f_{n}dv,\quad E_{0}=0.$ (42)

An advantage of this formulation is that Gauss’s law is solved exactly and
computationally is only a constraint associated with the system, a familiar
approach in incompressible flows. The velocity flux $F_{n}(v,t)$ is computed
pseudospectrally using zero-padded FFTs according to Orszag’s two-thirds rule
and is thus de-aliased [46]. Velocity space is divided into finite elements
and the coupled system of first-order PDEs for $f_{n}(v,t)$ is discretized
with the discontinuous Galerkin method using interpolants on Legendre-Gauss-
Lobatto (LGL) quadrature nodes [42]. As the flux is delocalized in spectral
space by the convolution product the Lax-Friedrichs flux is used rather than
upwinding [47]. The variable-coefficient term $ik_{n}vf_{n}$ is integrated
exactly up to the order of the basis polynomials and is thus also alias-free.

As both spatial and velocity fluxes are alias-free the method has good
conservation properties. For example, domain-integrated particle number is
conserved to machine precision. Energy, with no artificial hyperviscosity, is
conserved to $\mathcal{O}(10^{-10})$ for the resolution used. Spectral
blocking, when the turbulent cascade hits the Nyquist frequency of the spatial
grid, breaks this energy conservation. For this reason a hyperviscosity
$\nu\nabla_{x}^{4}f$ with $\nu=1.0\times 10^{-2}$ for $k\lambda_{D}>1$ is used
to introduce an artificial cutoff scale at the Debye length. This artificial
dissipation results in an energy loss $\mathcal{O}(10^{-7})$ by saturation
while electric energy saturates at $\mathcal{O}(10^{-3})$, meaning the
important dynamics are not adversely affected.

A non-uniform discretization is used, as in the bump-on-tail instability most
of the action occurs in the resonant band of velocities with $v\in(3,6)$ but
velocity space must be truncated at large $v$ for the above-discussed
particle, momentum, and energy conservation. Thus only a few elements are used
in the tails while most elements are clustered around the main bodies of the
two Maxwellian distributions. In this way sufficient resolution is obtained
for fine-scale velocity features without wasted resolution in the tails. A
semi-implicit time-integrator is used called the AB3CN method [46]. In this
use of the method, the linear advection term $ik_{n}vf_{n}$ is treated
implicitly with the second-order Crank-Nicholson method while the nonlinear
velocity flux is evolved explicitly by the third-order Adams-Bashforth
multistep method. This semi-implicit method is CFL-limited by the spectral
velocity flux rather than particle advection, enabling arbitrarily large
velocity domains.

In this work $N_{v}=50$ elements are used to discretize velocity space, each
with a ninth-order interpolant polynomial basis. Velocity space is truncated
at $v_{\text{max}}=\pm 25$ thermal velocities in order that total density be
conserved to machine precision. Thirty elements are clustered into the range
of approximately resonant velocities from $v\in(3,8)$ where high accuracy is
needed to resolve the fine trapping trajectories, with the remainder of
elements distributed in the main thermal body and tails. The physical domain
is taken as a periodic box of length of $L=5000\lambda_{D}$ and is divided
into $N_{x}=5000$ evenly-spaced nodes.

When the perturbation spectrum is initially peaked at
$\max(|E_{k}|^{2})=10^{-6}$ the instability saturates at $t\omega_{p}\approx
130$. The simulation is run to $t_{\text{max}}=170$ with a $\Delta t=2.5\times
10^{-3}$, requiring about fifteen minutes on a single RTX 3090 GPU. For a
given time-advance the numerical results are well converged, such that
discretization errors are below machine precision as verified by varying the
resolution. However some phase error is introduced by the implicit time
integration of spectral velocity advection as the very smallest scales are
slowed to the chosen time-advance [46]. As these smallest modes are damped by
the hyperviscosity no significant difference is seen in domain averages such
as $\langle E\delta f\rangle_{L}$ when compared to explicitly time integrated
simulations. This error does mean that care should be taken with the choice of
time-step as small scales are slowed in addition to large velocities.

##  V Analysis of simulation results

Here the Vlasov-Poisson and quasilinear simulations are analyzed, and the
discrepancy between the two simulations explored. Both results use the initial
background distribution specified in Eq. 10 and suppose a domain length of
$L=5000\lambda_{D}$. Section V.1 makes some comments on the choice of initial
condition. A comparison of the turbulent flux is made in Section V.2, along
with a discussion on the fluctuations in the field-particle correlation
$\langle\delta fE\rangle_{L}$ responsible for deviation from quasilinear
diffusion. Section V.3 analyzes the role of phase space structures in the
solution of the Vlasov equation followed by an analysis of the validity
conditions as discussed in Section III.2.

Figure 4: Domain-integrated energy traces for the Vlasov-Poisson simulation
of: a) electric field energy and b) change in kinetic energy. The bump relaxes
as electrons release free kinetic energy. A typical feature of nonlinear
saturation is an overshoot of equilibrium and slow nonlinear evolution, seen
here as an increase in kinetic energy past saturation at $t\omega_{p}\approx
125$. An artificial hyperviscosity is applied for wavenumbers
$k\lambda_{D}>1$, but total energy per unit length is conserved to
$\mathcal{O}(10^{-7})$.

###  V.1 On the choice of initial energy spectrum

The dynamic behavior of the quasilinear system as an initial value problem
depends on the choice of initial condition, because the spatially-averaged
distribution $\langle f\rangle_{L}$ and also the initial energy spectrum
$|E_{k}|^{2}$ are free to be specified. All modes must be energized as change
takes place only through linear growth or damping. A variety of choices exist
for the initial distribution of spectral field energy $|E_{k}|^{2}$. A
possible choice is a constant energy in all modes $E_{k}$ with random phases,
e.g. a white noise. Yet one cannot put energy into damped “eigenmodes” of the
Vlasov-Poisson system without exciting the other branches of the kinetic
dispersion relation, as Landau-damped modes are not true eigenfunctions [26].
If initialized, damped kinetic modes can couple into the dynamics of growing
modes. For example, Fig. 4 shows the temporal change of domain-integrated
electric and kinetic energies in the development of the instability. If
perturbations other than eigenmodes are initialized then the domain-integrated
energy will oscillate as it grows in time. The Vlasov-Poisson simulation is
perturbed only by unstable modes for better comparison to quasilinear theory.

Initialization of constant energy across all modes results in a cusped
diffusivity profile, nearly discontinuous at the edges of the resonant
velocities, because some wavenumbers have near-zero growth rates. Because of
this cusp effect the energy spectrum is initialized in both simulations to
match a scaled version of the positive growth rates with the maximum
$\alpha_{n}^{2}=10^{-6}$ in Eq. 11. In the quasilinear system the damped part
of the spectrum is initialized to the minimum positive growth rate in order
that the initial spectrum is continuous.

Figure 5: Early-time evolution of diffusivity $D(v)$ in: (a-b) the region of
resonant velocities; and (c-d) for non-resonant velocities in the main thermal
body. The Vlasov-Poisson diffusivity in (a) estimated with Eq. 43. The
diffusivities agree well until around $t\omega_{p}=40$, where they diverge
significantly. A key difference is the width of the diffusing region, as the
Vlasov simulation diffusivity at the resonant edges ($v=3.5$ and $v=5$) is
twice as large as the quasilinear prediction. Note that wave energies, plotted
in Fig. 4, are not yet at nonlinear levels. By $t\omega_{p}=50$ non-resonant
diffusivity of the Vlasov-Poisson simulation in (c) is roughly $30\%$ greater
than the diffusivity in (d) predicted by quasilinear theory for the
$L=5000\lambda_{D}$ interval.

###  V.2 Comparisons of turbulent diffusivity and momentum flux

Primary quantities of interest in phase space turbulence are the velocity-
space diffusivity and the turbulent momentum flux $\langle E\delta
f\rangle_{L}$, in particular its profile and magnitude, in the regions of
resonant and non-resonant velocities. Based on Eq. 19, a turbulent diffusivity
can be estimated from the Vlasov-Poisson simulation by decomposing the
distribution into a component that is averaged over physical space $\langle
f\rangle_{L}$ and a component fluctuating about this average, $\delta
f=f-\langle f\rangle_{L}$. Then diffusivity is estimated as

$D_{\text{estimated}}(v)\equiv\frac{\langle E\delta
f\rangle_{L}}{\partial_{v}\langle f\rangle_{L}}$ (43)

as the turbulent momentum flux normalized to the gradient of the background
distribution.

Figure 5 shows the Vlasov-Poisson and quasilinear velocity-space diffusivities
in the resonant (a-b) and non-resonant regions (c-d), demonstrating that
quasilinear theory correctly predicts the shape and early development of the
profile of turbulent diffusivity, but later presents an underestimate. By
$t\omega_{p}=50$ the velocity-space diffusivity is underestimated by between
30–100%. Regarding the non-resonant diffusivity, Fig. 5(c-d) shows that this
value varies by an order of magnitude across the positive velocities of the
main thermal body, becoming smaller for negative velocities. The variation of
diffusivity with velocity is important for the change in $\langle
f\rangle_{L}$ in the tail region, with the two simulations compared in Fig. 6.
The main thermal body is non-resonantly heated primarily in positive
velocities. The “notch” or depression between the main thermal body and the
bump in $\langle f\rangle_{L}$ is in this way filled by diffusion of the
distribution from near-resonant, but not fully resonant, velocities.

Figure 6: Diffusion of the distribution function in the tail region. Up to
this time the quasilinear solution agrees closely with the Vlasov solution,
though by $t\omega_{p}=50$ the Vlasov-Poisson simulation has a noticeably
smoother profile as the diffusivity has grown significantly greater than the
quasilinear value. By $t\omega_{p}=90$ the Vlasov simulation has completely
filled in the low-velocity notch, while the quasilinear simulation maintains a
steep notch area.

For later times in the Vlasov simulation fluctuations of $\langle
f\rangle_{L}$ make the diffusivity estimate unreliable as its gradient
oscillates. This happens to occur as the two models begin to diverge and the
diffusive closure breaks down. A smoother quantity to compare is the
covariance, or turbulent flux, $\langle E\delta f\rangle_{L}$. Figure 7
compares the quasilinear turbulent flux to that of the Vlasov-Poisson
simulation. This figure shows that the flux in the region of non-resonant
velocities close to the resonant region around $v\approx[3,6]$ is
underestimated by quasilinear theory. This enhanced turbulent flux flattens
the extreme regions of the unstable distribution, meaning the “notch” and
“bump” of the evolving tail. These two regions are flattened asymptotically in
the quasilinear model. The effect of this enhanced diffusion is to quench the
instability in a finite time, in contrast to quasilinear theory’s asymptotic
approach to quenching. This is seen as the saturation of electric energy in
Fig. 4(a).

Figure 7: Time-evolution of the field-particle correlation, or equivalently
turbulent momentum flux of the averaged distribution function, for (a, c) the
Vlasov-Poisson solution and (b, d) the quasilinear solution. The turbulent
flux of both systems peaks around $t\omega_{p}=50$ for this initial condition,
yet the observed flux of the Vlasov simulation is significantly greater than
the quasilinear value for all times beyond $t\omega_{p}=30$.

Another tool to quantify the quasilinear closure is the mean-square of the
fluctuation in the turbulent flux, i.e. the variance $\langle\langle E\delta
f\rangle\rangle_{L}$ plotted in Fig. 8. Recall the closure assumption Eq. 14
comparing a 0D1V-dimensional quantity to a 1D1V-dimensional one. This
comparison may be valid in some areas and invalid in others. The variance
quantifies the fluctuations in the turbulent flux and points out which
velocities contribute most to breaking the assumed closure. Most variance
occurs in the energy transfer corresponding to the non-resonant distribution.
This supports the evidence from Fig. 7 that a modification of the turbulent
diffusivity due to spatial structure, meaning wavepacket potentials, is
primarily responsible for the deviation of the turbulence evolution from the
quasilinear prediction.

Figure 8: Evolution of the variance of the turbulent flux, computed as the
mean-square in the fluctuating flux $E\delta f-\langle E\delta f\rangle_{L}$,
showing significant values in non-resonant velocities only. The variance of
the turbulent momentum flux for non-resonant particles by two orders of
magnitude over resonant ones supports the mechanism for the disagreement with
the quasilinear theory as due to spatial structure in the electric field, i.e.
beat waves.

###  V.3 Analysis of phase space structures in the Vlasov-Poisson simulation

Section III reviewed quasilinear theory and highlighted parallels to models of
eddy viscosity in the turbulence of an incompressible fluid. In such eddy
viscosity models a passive scalar is mixed by its partial participation in an
ensemble of fluid eddies, resulting in a diffusive character to the scalar’s
averaged profile. By observing the phase flow within the averaging box a
similar structure of eddies is observed due to phase fluid resonantly
interacting with phase waves.

Section V.2 shows that the Vlasov simulation evolves in two phases, at first
in agreement with quasilinear diffusion and later diverging as wave energy
increases. Figure 9 shows the development of the fluctuation $\delta f$ at an
early time corresponding to quasilinear evolution for a domain of
$L=1000\lambda_{D}$. Phase fluid mixes without the formation of closed eddies.
A shorter domain was chosen here to present a complete picture which can be
reasonably viewed as a whole. The simulation in Fig. 9 uses $5000$ spatial
nodes and does not use spatial hyperviscosity. Note that the checkerboard
pattern of the plasma oscillations is similar to the eigenmodes discussed in
Section II.1 and visualized in Fig. 2.

The Vlasov simulation begins to disagree with quasilinear theory as phase
space eddies form. These eddies are also referred to as clumps or granulations
and are a central ingredient of the theory of phase space structures [48, 49,
50]. These structures can be roughly divided into coherent (eigenmode-like)
and incoherent (eddy-like) parts. This section discusses the dynamics of the
eddy-like part of the fluctuation. An extensive textbook discussion of
granulations and their role in phase space turbulence can be found in Chapter
8 of Ref. [11]. Of course the analogy to eddy-viscosity models is only partial
since phase space eddies account only for resonant diffusion; it is well-known
from theory that coherent oscillations of the main thermal body are
responsible for energy transfer at non-resonant velocities [37, 51].

Figure 9: The fluctuation in the distribution function, $\delta f\equiv
f-\langle f\rangle_{L}$, of a Vlasov simulation is plotted during the
quasilinear phase of evolution at $t\omega_{p}=10$. A short domain
($L=1000\lambda_{D}$) is shown in order that the whole flow can be observable
at once. A clear distinction can be made between resonant and non-resonant
parts of the fluctuation. Non-resonant velocities undergo coherent
oscillations (seen as the striped pattern) while phase fluid at resonant
velocities ($v\in(3,5)$) is accelerated by randomly-phased waves. These
structures are consistent with quasilinear evolution.

####  V.3.1 Eddy analysis of the direct Vlasov simulation

One might expect each linear mode to evolve independently until terms
nonlinear in the wave amplitude dominate the dynamics. Yet as seen in the
Fourier spectrum of the fluctuation in Fig. 10(b), instability development
leads to a broadening of the initialized coherent crescent shape in the
distribution spectrum $f(k,v)$ and the gradual development of a power law in
these Fourier coefficients around the resonant band of velocities. Initial
broadening of the Fourier spectrum in the linear phase of the instability,
representing evolution of fine structures in the band of resonant velocities
from reaction on the mean momentum flux, is the essence of the quasilinear
effect. However the power law seen in Fig. 10(a) represents turbulent features
not captured by the quasilinear approximation.

Figure 10: a) The field spectrum $|E_{k}|^{2}$ of the Vlasov simulation
averaged from $t\omega_{p}=140$–$170$, consisting of two peaks at
$k\lambda_{D}\approx 0.25$ and $0.5$, as well an inertial range-like power law
for $k\lambda_{D}\gtrsim 0.6$. The primary peak is due to the saturated
instability and the secondary peak to its first harmonic. The spectral knee
and dissipation around $k\lambda_{D}=2$ is due to the artificial
hyperviscosity.
b) Fourier-decomposition of the distribution function, shown in logarithmic
view as $\log(1+|f|)$, demonstrating that the high-$k$ power law is due to
eddy turbulence of resonant electrons, while the energy-containing
oscillations are present in the population of non-resonant electrons.

An analysis of phase space eddies reveals how quasilinear closure breaks down
as this resonant-velocity power law from turbulent eddy mixing develops.
Eddies consist of electron phase fluid oscillating in a local potential well.
The typical electron bounce frequency at saturation, measured as the rms
charge density fluctuation, is $\omega_{b}\approx 0.1\omega_{p}$, extending to
bounce frequencies of $\omega_{b}\approx 0.2$–$0.4\omega_{p}$ in the largest
amplitude wavepackets. This gives a range of eddy turn-over times from
$\tau_{b}\omega_{p}\approx 15$–$60$. Resonant particles do not experience a
full period of bounce motion, as phase fluid transits each potential well in a
time

$\tau_{\text{transit}}\equiv\frac{\lambda}{v-v_{g}}$ (44)

with $v$ the local phase fluid velocity, $v_{g}$ the wavepacket group
velocity, and $\lambda$ the width of a wavepacket potential trough. In the
simulation considered the width of a wavepacket potential trough is roughly
$\lambda\approx 10\lambda_{D}$, the group velocity of wavepackets roughly one
thermal velocity, and the velocities of resonant particles
$v\approx(3,6)v_{t}$ (here estimated as greater than the original
$v\in(3,5)v_{t}$ as the resonant region expands with time). This gives transit
times of $\tau_{\text{transit}}\omega_{p}=2$–$5$. Thus groups of resonant
particles in the mixing layer participate in between $3$–$30\%$ of a typical
phase space eddy motion per trough. However, since the potential profile is
similar in each wavepacket trough, and the clump enters the trough
approximately as it left, a rotating clump can turn over in the process of
passing through multiple potential troughs. In this way a clump of phase fluid
is progessively distorted by its passage through the potential profile of a
wavepacket and trajectories randomized by the process of phase mixing in a
nonlinear potential.

A recurrence effect allows this progressive eddy turnover to be visualized in
a single-time snapshot of the phase flow. To understand the recurrence effect,
consider the fluctuation’s autocorrelation function $\langle\delta
f(x,v,t)\delta f(x-v_{g}\tau,v,t-\tau)\rangle_{L}$ [11], plotted in Fig. 11
for resonant velocities during early quasilinear evolution and a later time
when evolution deviates from quasilinear diffusion. During the quasilinear
phase the autocorrelation signature is sinusoidal but during later evolution
the signature changes to that of a structure repeating at the plasma
frequency. The repeating signal corresponds to a temporal similarity in the
phase flow. That is, one can approximately trace the past and future behavior
of a phase space clump in a given potential well by looking at the neighboring
potential wells because clumps transit between wells in approximately one
plasma period. Figure 12b shows this effect in detail, tracing the approximate
future state of a clump as it transits a wavepacket with a series of small
circles. In the way visualized, resonant particle mixing becomes governed by
progressive eddy turnover through the distinct potential wells of a
wavepacket. The eddies visualized are not trapped in any particular potential
well.

Figure 11: Autocorrelation $\langle\delta f(x,v,t)\delta
f(x-v_{g}\tau,v,t-\tau)\rangle_{L}$ of the fluctuating distribution function
in the Vlasov-Poisson simulation as a function of time-delay $\tau$ plotted
for: (a) the quasilinear phase of instability evolution, and (b) the
development of anomalously high diffusivity. A sinusoidal signal is indicative
of an underlying oscillatory pattern, while a positively-peaked signal
indicates a repeating pattern. As the instability evolves this pattern grows
increasingly peaked. By saturation there is only recurrence with no
oscillatory signal in resonant velocities. The oscillatory signal occurs
during times in agreement with quasilinear theory, while recurrence coincides
with the increase of turbulent flux beyond the quasilinear value. The resonant
fluctuation recurs at the plasma frequency, decorrelating into the past and
future over about six plasma periods. The decay in autocorrelation peaks
corresponds to the approximate coherence time of phase space clumps.

The phase fluid patterns have a long coherence time because the degree of
phase mixing in a sinusoidal potential depends on the time spent participating
in the eddy motion. Particle trajectories in a given clump slowly decorrelate,
with a decorrelation time given by the decay of peaks in Fig. 11. The change
from a sinusoidal signal to that of a repeating structure in the
autocorrelation function at resonant velocities is concurrent with the time of
disagreement of the turbulent flux from the quasilinear prediction. This is
the time for which the particle bounce frequency (or equivalently, phase space
eddy turnover time) becomes smaller than but comparable to the wavepacket
autocorrelation and resonant particle transit times.

Figure 12: Detailed phase-space views of the fluctuation $\delta f\equiv
f-\langle f\rangle_{L}$ in the Vlasov-Poisson simulation: a) outside a
wavepacket; and b) within a wavepacket. Phase space is shown at
$t\omega_{p}=120$ for the $1000\lambda_{D}$ domain. Null-points between
wavepackets can be determined from the phase-space view as the regions where
the coherent plasma oscillation sustained by non-resonant electrons goes to
zero amplitude. The mixing mechanism depends on the local potential structure.
Outside wavepackets filaments mainly free-stream while within wavepackets
eddies mix. The long decorrelation time of a clump is illustrated in (b) by a
circle tracking a recurrent clump, demonstrating that an eddy can complete a
turn as it propagates from trough-to-trough, even though the rotating clump is
experiencing different potential amplitudes in each trough.

####  V.3.2 Wavepacket autocorrelation and quasilinear applicability

Besides its oscillation frequency a wavepacket is characterized by its
autocorrelation time $\tau_{\text{ac}}$, or the lifetime of the potential
profile in the wavepacket, and the lifetime of the wavepacket
$\tau_{\text{life}}$ as an amplitude envelope. The autocorrelation time,
similar to the particle potential transit time, is estimated as
$\tau_{ac}=2\pi|(\zeta-v_{g})\Delta k|^{-1}$ as in Section III.2. The quantity
$\Delta k$ is the spectral width of the energy-containing oscillations, and
can be estimated from the saturated spectrum, shown in Fig. 10(a), as $\Delta
k\approx 0.1$. The phase velocities range from $\zeta\approx 3$–$5v_{t}$, so
with $v_{g}\approx 1v_{t}$ the autocorrelation frequency is approximately
$\omega_{\text{ac}}=0.4$.

The importance of wavepacket lifetime can be estimated using the Hilbert
transform to obtain the wave envelope $\psi$, and computing the group
velocity-shifted autocorrelation function
$\langle\psi(x,t)\psi(x-v_{g}\tau,t-\tau)\rangle_{L}$. When this shifted
autocorrelation is compared against the unshifted autocorrelation of the
envelope during the saturated state of the instability a stationary signal is
observed. The slow evolution of the wavepacket envelope points to the
robustness of the saturated state of oscillation.

Returning to the question of quasilinear applicability, the relaxation time is
first estimated from finite differences of the Vlasov-Poisson averaged
distribution $\langle f\rangle_{L}$ as
$\tau_{r}^{-1}\approx\mathcal{O}(10^{-3})\omega_{p}$. The maximum growth rate
is initially $\omega_{i}\approx 0.03\omega_{p}$. Thus the inequality
$\tau_{r}^{-1}<\omega_{i}<\tau_{ac}^{-1}$ is satisfied as
$10^{-3}<10^{-2}<10^{-1}$. The rms bounce frequency of the perturbation is
$\omega_{b}\approx 0.05\omega_{p}$, but by saturation the bounce frequency has
grown to $\omega_{b}\approx 0.1\omega_{p}$ with maximum values of up to
$\omega_{b}=0.4\omega_{p}$ in the largest amplitude wavepackets. While the
instability saturates at only $0.2\%$ of the total thermal energy of the
system this is enough to make the bounce frequency comparable to the
autocorrelation frequency and the quasilinear flux is found to be modified
[14] as discussed previously.

##  VI Summary of findings

This study revisited the bump-on-tail instability with a conservative high-
order discontinuous Galerkin method applied to its quasilinear theory as well
as to a direct simulation of the Vlasov-Poisson equations. The initial
condition was chosen in the regime of beam parameters for which quasilinear
theory should be applicable according to Ref.[28]. It was found that
quasilinear theory is initially in good agreement with the system evolution
but underestimates the rate of turbulent momentum flux beginning at field
energies halfway to instability saturation. The result of the correction
beyond the quasilinear approximation is a quenching of the instability in a
finite time and the replacement of the unstable initial condition by a robust
state of oscillation. This correction occurs as the phase space eddy turnover
time in the largest amplitude wavepackets becomes comparable to the phase
fluid transit time through the wavepacket potentials. Progressive eddy
turnover transitions diffusive scattering of resonant electrons by randomly-
phased waves to their nonlinear phase mixing in potential wells. That the
wavepackets, or localized concentrations of fluctuation-to-mean energy
transfer, are responsible for the deviation from quasilinear diffusion can
also be understood as the breakdown of the quasilinear closure assumption by
the development of significant spatial variance to the turbulent momentum
flux. The findings of this article support the theory of phase space density
granulations [11, 48, 49] and shed light on the distribution of fluctuations
with velocity in collisionless plasma turbulence.

Limitations of quasilinear theory have been known for some time. Theories to
account for the discrepancies discussed in this article include resonance
broadening theory and expansion to higher-order in the field amplitude [3, 12,
11]. These approaches can considerably complicate the analysis without a
necessarily rewarding conclusion, as energy conservation in resonance
broadening theory is still actively researched [3] and higher-order expansions
may lead to misleading analysis for nonlinear Landau damping processes [52].
The findings of this study suggest that using quasilinear theory to predict
plasma processes for all but the weakest instabilities should be considered as
only order-of-magnitude accurate.

##  VII Acknowledgments

The authors would like to thank I.A.M. Datta, W.H. Thomas, A. Ho, E.T. Meier,
A.D. Stepanov, and G.V. Vogman for helpful discussions. The information, data,
or work presented herein was funded in part by the Air Force Office of
Scientific Research under award No. FA9550–15–1–0271. This material is also
based upon work supported by the National Science Foundation under Grant No.
PHY-2108419.

##  VIII Author declarations

The authors have no conflicts to disclose.

##  IX Data availability statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Drummond and Pines [1962] W. E. Drummond and D. Pines, Nuclear Fusion Supplement 3, 1049 (1962).
  * Vedenov, Velikhov, and Sagdeev [1962] V. V. Vedenov, E. P. Velikhov, and R. Z. Sagdeev, Nuclear Fusion Supplement 2, 465 (1962).
  * Krommes [2015] J. A. Krommes, Journal of Plasma Physics 81 (2015).
  * Balescu [2005] R. Balescu,  _Aspects of anomalous transport in plasmas_ (CRC Press, 2005).
  * Tobias and Marston [2017] S. M. Tobias and J. B. Marston, Journal of Fluid Mechanics 810, 412–428 (2017).
  * Davidson [1978] R. C. Davidson, The Physics of Fluids 21, 1375 (1978), https://aip.scitation.org/doi/pdf/10.1063/1.862379 .
  * Seough, Yoon, and Hwang [2015] J. Seough, P. H. Yoon, and J. Hwang, Physics of Plasmas 22, 012303 (2015), https://doi.org/10.1063/1.4905230 .
  * Allanson  _et al._ [2021] O. Allanson, C. E. J. Watt, H. J. Allison, and H. Ratcliffe, Journal of Geophysical Research: Space Physics 126, e2020JA028793 (2021).
  * Albert [2010] J. M. Albert, Journal of Geophysical Research: Space Physics 115 (2010), https://doi.org/10.1029/2009JA014732.
  * Cardy, Falkovich, and Gawedzki [2008] J. Cardy, G. Falkovich, and K. Gawedzki,  _Non-equilibrium statistical mechanics and turbulence_ (London Mathematical Society, 2008).
  * Diamond, Itoh, and Itoh [2010] P. Diamond, S. Itoh, and K. Itoh,  _Modern plasma physics_ , Vol. 1 (Cambridge University Press, 2010).
  * Davidson [1972] R. C. Davidson,  _Methods in nonlinear plasma theory_ (Academic Press, 1972).
  * Henri  _et al._ [2011] P. Henri, F. Califano, C. Briand, and A. Mangeney, EPL (Europhysics Letters) 96, 55004 (2011).
  * Vlad  _et al._ [2004] M. Vlad, F. Spineanu, J. Misguich, J. D. Reuss, R. Balescu, K. Itoh, and S. Itoh, Plasma physics and controlled fusion 46, 1051 (2004).
  * Seough, Yoon, and Hwang [2014] J. Seough, P. H. Yoon, and J. Hwang, Physics of Plasmas 21, 062118 (2014), https://doi.org/10.1063/1.4885359 .
  * Lafleur  _et al._ [2018] T. Lafleur, R. Martorelli, P. Chabert, and A. Bourdon, Physics of Plasmas 25, 061202 (2018), https://doi.org/10.1063/1.5017626 .
  * Rajawat and Sengupta [2017] R. S. Rajawat and S. Sengupta, Physics of Plasmas 24, 122103 (2017), https://doi.org/10.1063/1.5006463 .
  * Cohen  _et al._ [1989] B. I. Cohen, A. B. Langdon, D. W. Hewett, and R. J. Procassini, Journal of Computational Physics 81, 151 (1989).
  * Brackbill [2016] J. U. Brackbill, Journal of Computational Physics 317, 405 (2016).
  * Horkỳ, Miloch, and Delong [2017] M. Horkỳ, W. J. Miloch, and V. A. Delong, Physical Review E 95, 043302 (2017).
  * Heath  _et al._ [2012] R. E. Heath, I. M. Gamba, P. J. Morrison, and C. Michler, Journal of Computational Physics 231, 1140 (2012).
  * Cheng, Gamba, and Morrison [2013] Y. Cheng, I. M. Gamba, and P. J. Morrison, Journal of Scientific Computing 56, 319 (2013).
  * Juno  _et al._ [2018] J. Juno, A. Hakim, J. TenBarge, E. Shi, and W. Dorland, Journal of Computational Physics 353, 110 (2018).
  * Goldston and Rutherford [1995] R. J. Goldston and P. H. Rutherford,  _Introduction to plasma physics_ (Taylor and Francis, 1995).
  * Landau [1946] L. D. Landau, J. Phys. USSR 10 (1946).
  * Case [1959] K. M. Case, Annals of physics 7, 349 (1959).
  * van Kampen [1955] N. G. van Kampen, Physica 21, 949 (1955).
  * Thorne and Blandford [2017] K. S. Thorne and R. D. Blandford,  _Modern classical physics: optics, fluids, plasmas, elasticity, relativity, and statistical physics_ (Princeton University Press, 2017).
  * Frisch [1995] U. Frisch,  _Turbulence: the legacy of A.N. Kolmogorov_ (Cambridge university press, 1995).
  * van Kampen [2007] N. G. van Kampen,  _Stochastic processes in physics and chemistry_ , 3rd ed. (Elsevier, 2007).
  * Howes, Klein, and Li [2017] G. G. Howes, K. G. Klein, and T. C. Li, Journal of Plasma Physics 83 (2017), 10.1017/S0022377816001197.
  * Li  _et al._ [2019] T. C. Li, G. G. Howes, K. G. Klein, Y. H. Liu, and J. M. TenBarge, Journal of Plasma Physics 85 (2019), 10.1017/S0022377819000515.
  * Klein  _et al._ [2020] K. G. Klein, G. G. Howes, J. M. TenBarge, and F. Valentini, Journal of Plasma Physics 86 (2020), 10.1017/S0022377820000689.
  * Juno  _et al._ [2021] J. Juno, G. G. Howes, J. M. TenBarge, L. B. Wilson, A. Spitkovsky, D. Caprioli, K. G. Klein, and A. Hakim, Journal of Plasma Physics 87 (2021), 10.1017/S0022377821000623.
  * Schmitt [2007] F. G. Schmitt, Comptes Rendus Mécanique 335, 617 (2007).
  * Kadomtsev [1965] B. B. Kadomtsev,  _Plasma turbulence_ (Academic Press, 1965).
  * Dewar [1973] R. L. Dewar, The Physics of Fluids 16, 1102 (1973).
  * Giga, Giga, and Kobayashi [1999] M. H. Giga, Y. Giga, and R. Kobayashi, Hokkaido University Preprint Series in Mathematics 461, 1 (1999).
  * Escande and Sattin [2007] D. F. Escande and F. Sattin, Physical Review Letters 99 (2007).
  * Besse  _et al._ [2011] N. Besse, Y. Elskens, D. F. Escande, and P. Bertrand, Plasma Physics and Controlled Fusion 53 (2011).
  * Elskens and Escande [2019] Y. Elskens and D. Escande,  _Microscopic dynamics of plasmas and chaos_ (CRC Press, 2019).
  * Hesthaven and Warburton [2007] J. S. Hesthaven and T. Warburton,  _Nodal discontinuous Galerkin methods: algorithms, analysis, and applications_ (Springer Science & Business Media, 2007).
  * Shumlak  _et al._ [2011] U. Shumlak, R. Lilly, N. Reddell, E. Sousa, and B. Srinivasan, Computer Physics Communications 182, 1767 (2011).
  * Ho, Datta, and Shumlak [2018] A. Ho, I. A. M. Datta, and U. Shumlak, Frontiers in Physics 6, 105 (2018).
  * Datta, Crews, and Shumlak [2021] I. A. M. Datta, D. W. Crews, and U. Shumlak, Physics of Plasmas 28, 072112 (2021).
  * Boyd [2000] J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ , Second ed. (Dover, 2000).
  * LeVeque  _et al._ [2002] R. J. LeVeque _et al._ ,  _Finite volume methods for hyperbolic problems_ , Vol. 31 (Cambridge University Press, 2002).
  * Dupree [1970] T. H. Dupree, Physical Review Letters 25, 789 (1970).
  * Dupree [1972] T. H. Dupree, The Physics of Fluids 15, 334 (1972).
  * Boutros-Ghali and Dupree [1981] T. Boutros-Ghali and T. H. Dupree, The Physics of Fluids 24, 1839 (1981).
  * Cary and Kaufman [1981] J. R. Cary and A. N. Kaufman, The Physics of Fluids 24, 1238 (1981).
  * Mouhot and Villani [2011] C. Mouhot and C. Villani, Acta Mathematica 207, 29 (2011).