Analytical solutions for nonlinear plasma waves with time-varying complex frequency

TL;DR

This paper develops analytical methods to study the evolution of nonlinear plasma waves with time-varying complex frequencies, extending classical BGK solutions to include frequency sweeping modes and their nonlinear dynamics.

Contribution

It introduces a formalism for analyzing plasma waves with time-dependent complex frequencies using Taylor expansions and basis decompositions, providing new analytical solutions and insights.

Findings

Derived solutions for time-varying complex frequencies in plasma waves.

Presented approximate solutions for constant f contours in nonlinear regimes.

Identified highly nonlinear orbits with larger velocity widths than BGK islands.

Abstract

Bernstein-Kruskal-Greene (or BGK) modes are ubiquitous nonlinear solutions for the 1D electrostatic Vlasov equation, with the particle distribution function $f$ given as a function of the particle energy. Here, we consider other solutions $f = f[\epsilon]$ where the particle energy is equal to the second-order velocity space Taylor expansion of the function $\epsilon(x,v,t)$ near the wave-particle resonance. This formalism allows us to analytically examine the time evolution of plasma waves with time-varying complex frequency $\omega(t) + i \gamma(t)$ in the linear and nonlinear phases. Using a Laplace-like decomposition of the electric potential, we give allowed solutions for the time-varying complex frequencies. Then, we show that $f$ can be represented analytically via a family of basis decompositions in such a system. Using a Gaussian decomposition, we give approximate solutions for contours of constant $f$ for a single stationary frequency mode, and derive the evolution equation for the nonlinear growth of a frequency sweeping mode. For this family of modes, highly nonlinear orbits are found with the effective width in velocity of the island roughly a factor of $\sqrt{2}$ larger than the width of a BGK island.

Figures

Equations

$$\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}-\frac{q}{m}\frac{\partial \varphi }{\partial x}\frac{\partial f}{\partial v}=\hat{C}f$$

$$\underset{v\to \pm \infty }{\lim }f(x,v,t)\equiv 0$$

$$\left|\underset{C}{\int }f(x,v,t)\mathrm{d}\ell \right|<\infty$$

$$M_l[S]\equiv \underset{-\infty }{\overset{\infty }{\int }}{v}^{l}S(v)\mathrm{d}v$$

$$|M_l[f]|<\infty \forall l\ge 0$$

$$M_0[f]\equiv n(x,t)$$

$$M_1[f]\equiv \frac{1}{q}J(x,t)$$

$$M_2[f]\equiv \frac{2}{m}T(x,t)$$

$$0={\mu }_0\underset{-\infty }{\overset{\infty }{\int }}qvf\mathrm{d}v-{\mu }_0{ϵ}_0\frac{\partial }{\partial t}\frac{\partial \varphi }{\partial x}$$

$$\underset{-\infty }{\overset{\infty }{\int }}\left[\frac{1}{2}m{v}^2\frac{\partial f}{\partial t}+qvf\frac{\partial \varphi }{\partial x}\right]\mathrm{d}v\equiv {\left(\frac{\partial T}{\partial t}\right)}_{\text{C,T}}$$

$${\left(\frac{\partial T}{\partial t}\right)}_{\text{C,T}}\equiv \frac{1}{2}m\underset{-\infty }{\overset{\infty }{\int }}{v}^2({\hat{C}}^{\prime }f)\mathrm{d}v$$

$$\underset{-\infty }{\overset{\infty }{\int }}\left[\frac{1}{2}m{v}^2\frac{\partial f}{\partial t}\right]\mathrm{d}v+\frac{{ϵ}_0}{2}\frac{\partial }{\partial t}{\left(\frac{\partial \varphi }{\partial x}\right)}^2={\left(\frac{\partial T}{\partial t}\right)}_{\text{C,T}}$$

$${\left(\frac{\partial T}{\partial t}\right)}_{\text{C,T}}=-{\gamma }_D\frac{{ϵ}_0}{2}{\left(\frac{\partial \varphi }{\partial x}\right)}^2$$

$$\frac{\partial }{\partial t}({\hat{C}}^{\prime }f)=-{\gamma }_D\left[({\hat{C}}^{\prime }f)-\frac{\partial f}{\partial t}\right]$$

$${\hat{C}}^{\prime }f=e^{-{\gamma }_{D}t}\left[{\frac{\partial f}{\partial t}|}_{t=0}+{\gamma }_D\underset{0}{\overset{t}{\int }}{e}^{{\gamma }_D\tau }\frac{\partial f}{\partial \tau }\mathrm{d}\tau \right]$$

$$\hat{C}f\approx {\frac{\partial f}{\partial t}|}_{t=0}{e}^{-{\gamma }_{D}t}+{\gamma }_D(f-F)+\mathcal{O}({\gamma }_D^2)$$

$$f(x,v,t)=\frac{1}{2}\sum _j\left[f_j(v,t)e^{\text{i}{k}_{j}x}+\text{c.c.}\right]$$

$$p{|}_{t=0}\approx \frac{q^2}{m{ϵ}_0}\underset{\Omega }{\int }{\frac{\partial f_0}{\partial v}|}_{t=0}\frac{v\mathrm{d}v}{p{|}_{t=0}-{\gamma }_D+\text{i}{k}_{j}v}$$

$$\frac{\partial }{\partial t}\frac{\partial {\varphi }_{\text{alt.}}}{\partial x}=\frac{q}{{ϵ}_0}\underset{-\infty }{\overset{\infty }{\int }}vf\mathrm{d}v-{\gamma }_D\frac{\partial {\varphi }_{\text{alt.}}}{\partial x}$$

$$p_{\text{alt.}}{|}_{t=0}+{\gamma }_D=\frac{q^2}{m{ϵ}_0}\underset{\Omega }{\int }{\frac{\partial f_0}{\partial v}|}_{t=0}\frac{v\mathrm{d}v}{p_{\text{alt.}}{|}_{t=0}+\text{i}{k}_{j}v}$$

$$\frac{{\partial }^2}{\partial t^2}\left(\frac{\partial \varphi }{\partial x}\right)=\frac{q}{{ϵ}_0}\underset{-\infty }{\overset{\infty }{\int }}v\frac{\partial f}{\partial t}\mathrm{d}v$$

$$\begin{array}[]{r l}\dfrac{\partial^{2}}{\partial t^{2}}\left(\dfrac{\partial\phi}{\partial x}\right)&=-\overbrace{\dfrac{q^{2}}{m\epsilon_{0}}n}^{\textrm{oscillation}}\left(\dfrac{\partial\phi}{\partial x}\right)\\ &\hskip 30.0pt+\underbrace{\dfrac{q}{\epsilon_{0}}\int\limits_{-\infty}^{\infty}v(\hat{C}^{\prime}f)\mathop{}\!\mathrm{d}v}_{\textrm{drive}}\end{array}$$

$$\varphi (x,t)\equiv \sum _l{\varphi }_l(t)$$

$$\begin{array}[]{r l}\phi_{l}(x,t)&\equiv\dfrac{1}{2}\displaystyle\sum\limits_{s}\Bigg{\{}\phi_{sl}\exp\left[ik_{l}x+\int\limits_{0}^{t}p_{sl}(\tau)\mathop{}\!\mathrm{d}\tau\right]\\ &\hskip 50.0pt+\textrm{ c.c.}\Bigg{\}}\end{array}$$

$$p_{sl}(t)={\gamma }_{sl}(t)-i{\omega }_{sl}(t)$$

$$\frac{\mathrm{d}{p}_{sl}}{\mathrm{d}t}=-p_{sl}^2+{\Gamma }_{sl}^2$$

$$\begin{array}[]{r l}\Gamma_{sl}^{2}&=-\omega_{\textrm{pl.}}^{2}+\dfrac{\omega_{\textrm{pl.}}^{2}}{\bar{n}_{0}q(\partial_{x}\phi)_{sl}}\bigg{\{}\displaystyle\int\limits_{-\infty}^{\infty}\bigg{[}(f-F)\partial_{x}\phi\\ &\hskip 20.0pt-mv[\hat{C}^{\prime}f]_{sl}\bigg{]}\mathop{}\!\mathrm{d}v\bigg{\}}_{sl}\end{array}$$

$$p_{sl}(t)={\Gamma }_{sl}\tanh \left[\phi +{\Gamma }_{sl}t\right]$$

$$\forall {\Gamma }_{sl}\notin \{0,p_{sl}\}:\underset{t\to \infty }{\lim }{p}_{sl}(t)=|\Re ({\Gamma }_{sl})|\pm i\Im ({\Gamma }_{sl})$$

$${\varphi }_l(x,t)=\sum _s\left\{\frac{{\varphi }_{sl}}{2}{\text{e}}^{i{k}_{l}x}\cosh \left[\phi +{\Gamma }_{sl}t\right]+\text{c.c.}\right\}$$

$$\underset{\phi \to \infty }{\lim }\cosh \left[\phi +{\Gamma }_{sl}t\right]=:e^{p_{sl}(0)t}$$

Full text

I Introduction

Kinetic plasmas can be described by a Boltzmann equation wesson2011tokamaks that describes the evolution of the particle distribution function, $f$. The force exerted on each of the particles in the system manifests as a term acting as $\mathbf{F}\cdot\nabla_{v}$ on the distribution function, where $\nabla_{v}$ denotes a velocity-space gradient, and $\mathbf{F}$ is typically taken to be the Lorentz force. This leads to nonlinear coupling between the electromagnetic field and the particle distribution function. As a result of this nonlinearity, it is difficult to solve systems with kinetic plasmas analytically.

The method of linearisation allows one to derive Landau damping or inverse Landau damping of electrostatic plasma waves in the system under the assumption that the nonlinear coupling is sufficiently small landau1946vibrations . However, this coupling leads to wave-wave instabilities which affect the nonlinear stability of the system. This can lead to destabilisation rosenbluth1997theory ; lesur2016nonlinear as well as stabilisation lin1998turbulent ; anderson2017statistical . These nonlinear scenarios have been partially investigated analytically, however these systems are typically examined using computational means vann2003fully ; lang2010nonlinear ; gorelenkov2018resonance ; woods2018stochastic .

The phenomenon of solitons is widely documented korteweg1895change ; zabusky1965interaction ; emplit1987picosecond , where nonlinear structures form due to spatial coupling between waves in the system. In kinetic plasmas, nonlinear coupling in momentum space allows for structures to form on the distribution function known as ‘holes’ and ‘clumps’, referring to a local decrease or increase respectively on the spatially averaged distribution function berk1997spontaneous . These structures propagate in a non-dispersive manner similar to solitons. However, as they move through momentum space they draw free energy from the system even when at constant amplitude.

Understanding how these structures are formed may prove crucial for mitigating fast ion loss in tokamaks, as well as other types of kinetic instability. In tokamaks, wave-particle resonance allows for free energy exchange between particles and the electromagnetic field, depending on momentum space gradients. As holes and clumps form, they perturb the local distribution function, allowing for destabilisation of other waves in the system. This cascading effect may play a role in the observation of abrupt large events or mode avalanching in tokamaks sharapov2002alfven ; gryaznevich2006perturbative ; podesta2009experimental , where large spikes in the intensity of magnetic perturbations over a broadband of frequencies can suddenly occur. This is correlated with large amounts of fast ion loss, sharply degrading operational performance of the fusion device.

The Vlasov equation vlasov1938vibration has been solved in previous work for the case where the particles follow orbits on which their energy is conserved. bernstein1957exact . These types of plasma wave, known as Bernstein-Green-Kruskal (or BGK) modes, result in periodic orbits in phase-space ‘islands’ where trapped particles resonate with the wave. Further work by Berk, Breizman and others berk1995numerical ; berk1996nonlinear ; berk1997spontaneous has explored the movement of holes and clumps through phase-space by assuming that the hole and clump evolve via a temporal continuum of BGK modes. While this work has proven very successful for predicting the frequency sweep associated with these structures, the work does not fully explain how these structures form, and assume a constant slope on the distribution function during structure formation. In addition, this theory is limited to systems near marginal stability; rapidly growing waves and rapidly sweeping waves break the required adiabaticity conditions.

In this vein, work by Lilley et al. has shown computationally how holes and clumps form via phase mixing of the wave within a phase-space island lilley2010effect . Work by Wang et al. has recently examined frequency chirping for the late-time evolution of plasma waves using a mixture of analytical theory and the CHIRP code wang2018frequency . However, the threshold for hole-and-clump formation and the formation of holes and clumps on distribution functions with large curvature are both still relatively unexplored. As a result, understanding of hole-and-clump destabilisation is limited, and associated stability thresholds have not been fully investigated.

In Section II, we review basic properties of the distribution function and its moments. In Section III, we derive a Krook-like collision operator, and link this exact form to the complex linear dispersion relation for waves near marginal stability in the literature berk1997spontaneous . We show that a similar dispersion relation for modes undergoing Landau resonance can be captured in a heuristically similar collisionless model.

In Section IV, we formalise a nonlinear spectral decomposition which allows for nonlinear damped dynamics for the electrostatic potential. Then, in Section V we show a general family of exact solutions for the Vlasov equation for waves with time-varying complex frequency, defined by a non-holonomic constraint on the system.

In Section VI, we show that it is possible to restrict the family solutions further such that $f=f[\epsilon]$, where the particle energy is equal to the second-order velocity space Taylor expansion of the function $\epsilon(x,v,t)$ near a wave-particle resonance at $v=u$. Here, the single particle orbits are such that the effective potential is given by $\phi$ plus a velocity dependent term (which is identically zero in BGK models), allowing for velocity space ‘screening’ of the resonance.

Later in Section VI, we represent $f$ as a function of the energy-like function $\epsilon$, then show via a suitable basis expansion that one can easily yield approximate solutions for static frequency waves and chirping waves, with an ‘effective’ width in $v$ of the island a factor of $\sqrt{2}$ larger than the BGK island width. Finally, we show that it is possible to construct a ‘first-order’ distribution function that accurately reproduces every moment of the distribution function.

II The distribution function

The evolution of the single species particle distribution function is given by the 1D electrostatic Boltzmann equation with the force field given by the Lorentz force wesson2011tokamaks :

[TABLE]

where $\hat{C}$ denotes a differential operator acting on the particle distribution function $f(x,v,t)$, $x\in\mathbb{R}$ is position, $v\in\mathbb{R}$ is velocity, $t\in\mathbb{R}_{\geq 0}$ is time, $q$ is the particle charge, $m$ is the particle mass, and $\phi(x,t)$ is the electric potential. The $\hat{C}$ operator encapsulates particle collisions, and other terms which allow for conservation of the rate of change of every moment of the distribution function.

II.1 Properties of the distribution function

II.1.1 Vanishing at infinity

The scalar field $f(x,v,t)$ is required to vanish as $|v|\to\infty$:

[TABLE]

II.1.2 Integrable

The scalar field $f(x,v,t)$ is integrable on its submanifolds $x\in\mathbb{R}$ and $v\in\mathbb{R}$:

[TABLE]

where $C$ is a contour that spans two points on $\mathbb{R}\times\mathbb{R}$.

II.1.3 Moments

The $l^{\textrm{th}}$ moment of the scalar field $S(v)$ is defined as:

[TABLE]

Each and every moment $l\geq 0$ of $f(x,v,t)$ is finite, that is to say:

[TABLE]

We list the following useful moments of the distribution function:

[TABLE]

where $n(x,t)$ is the particle density, $J(x,t)$ is the current density, and $T(x,t)$ is the kinetic energy density.

III Evolution equations

The evolution of the electric potential is given here by the Maxwell-Ampère law:

[TABLE]

where we have taken ${B}=0$.

III.1 Energy balance

By taking the second moment of equation (1):

[TABLE]

where $(\partial T/\partial t)_{\textrm{C,T}}$ allows for sources and sinks of energy density:

[TABLE]

where $\hat{C}^{\prime}f\equiv\hat{C}f-v\partial f/\partial x$ encapsulates collisions as well as spatial advection. Use of equation (7) allows one to find an energy conservation equation:

[TABLE]

where the first term denotes the rate of change of energy density in the distribution function, the second term denotes the rate of change of electric field energy density, and the third term denotes sources and sinks of energy density.

III.2 Berk-Breizman sink

Berk and Breizmanberk1995numerical use an energy sink which dissipates plasma waves at a rate $\gamma_{D}$. This allows for a region of nonlinear stability for the system in the presence of instabilities.

One can represent this sink as the following:

[TABLE]

Therefore, by using equations (9), (10) and (11) one finds that for a vanishing integrand under $v$:

[TABLE]

This is an inhomogeneous decay equation of the form $y(t)^{\prime}=-ay+b(t)$. The solution, given by integrating factor, is:

[TABLE]

If one assumes the contribution from heat diffusion to be negligible in comparison to the contribution from collisions and other sources and sinks, then a Maclaurin expansion in $\gamma_{D}$ yields:

[TABLE]

where $F\equiv f(x,v,0)$, such that a Krook-like collision operator safely approximates the behaviour to first order in $\gamma_{D}$. The higher order terms also vanish for early times where $\gamma_{D}t\ll 1$.

III.3 Complex linear dispersion relation

We spectrally decompose the distribution function and other spatially dependent quantities as follows:

[TABLE]

where $\{k_{j}\}$ is a set of wavenumbers, and c.c. denotes the complex conjugate. Using the Berk-Breizman sink, it is possible to linearise equations (1) and (7) to produce the following complex linear dispersion relation for a single mode, under the Krook-like approximation:

[TABLE]

where $\Omega$ is the suitable Landau contour for the problem, and $p=\gamma-\textrm{i}\omega$. The linear growth rate is given by $\gamma_{L}=\gamma(t=0)$, and $\omega(t)$ is the frequency of the mode.

Commonly in the literature, the ‘modified’ version of the Maxwell-Ampère law is used instead to generate the same energy balance equation as equation (10):

[TABLE]

where $\phi_{\textrm{alt.}}$ is a modified electric potential. For these models to be consistent, Maxwell’s equations must retain the same canonical form; previous work shows that transformation of the system Lagrangian that preserves the system Hamiltonian is indeed possible under certain conditions. However, for all counterexamples, these models do not form a fully consistent set of equations. woods2018stochastic

Using this alternative model, one obtains a very similar complex linear dispersion relation woods2018stochastic :

[TABLE]

This means that while the model is inconsistent (energy is typically not conserved), it produces similar dynamics to a model with Krook-like collisions; the primary difference is that the location of the Landau resonance is shifted from $p=\gamma_{D}-\textrm{i}k_{j}v$ to $p=-\textrm{i}k_{j}v$.

In the nonlinear phase of the systems evolution, it is a well known result that the original resonance undergoes broadening. Therefore, provided that $\gamma_{D}/k_{j}$ is much smaller than the resonance width, these models perform well except for energy conservation; as we have already utilised a small $\gamma_{D}$ approximation to assume Krook-like dissipation, we are safe to make this assumption.

For this reason, in the rest of the paper we will explore $\hat{C}f=0$, such that the original and ‘modified’ versions of the Maxwell-Ampère law are the same. We expect that an extension to this work featuring a Berk-Breizman sink such that $\hat{C}f\neq 0$ could be investigated approximately by still considering a collisionless kinetic equation, and considering a model using the ‘modified’ Maxwell-Ampére law. This would allow one to retain the formalism derived here for the collisionless case, provided that the damping is small.

IV Nonlinear damped dynamics

IV.1 Complex frequency evolution

By coupling the Maxwell-Ampère law to the Boltzmann equation, we obtain a system of differential equations which features strong nonlinearities. By taking the time derivative of equation (7):

[TABLE]

From here, one can show that substituting $\partial_{t}f$ for equation (1) gives an equation for $(\partial_{x}\phi)$ similar to driven simple harmonic motion:

[TABLE]

However, one should note that $f=f(x,v,t)$, and therefore this is only part of a coupled system of equations. We shall solve equation (20) by representing the electric potential $\phi$ via a set of wavepackets $\{\phi_{l}\}$:

[TABLE]

Then, we employ the following solution for each wavepacket:

[TABLE]

where $\{k_{l}\}$ is a set of wavenumbers, $\{\phi_{sl}\}$ is a set of wave amplitudes, and $\{p_{sl}\}$ is a set of complex frequencies:

[TABLE]

where $\{\gamma_{sl}(t)\}$ is a set of nonlinear growth rates, and $\{\omega_{sl}(t)\}$ is a set of constituent frequencies. Using this solution, one finds the following equation:

[TABLE]

where $\{\Gamma_{sl}\}$ is a set of constants, corresponding to rates of change that determine the dynamics of the system:

[TABLE]

where $\omega_{\textrm{pl.}}=\sqrt{\bar{n}_{0}q^{2}/m\epsilon_{0}}$ is the plasma frequency, $F(v)=f(x,v,t=0)$ is the distribution function at a time $t=0$, and $\bar{n}_{0}$ is the spatial average of the particle density at a time $t=0$.

As there are multiple frequencies for $\phi_{l}(t)$, it is instructive to note that the complex frequency of the $l^{\textrm{th}}$ mode is made up of contributions from all of the branches.

Theorem 1**.**

There is no solution where $\lim_{t\to\infty}\gamma_{sl}\leq 0$ except for the case where $\gamma_{sl}$ is constant.

Proof.

Equation (24) yields:

[TABLE]

where the complex angle $\varphi\equiv\textrm{arctanh}(p_{sl}(0)/\Gamma_{sl})$. We show plots of $p_{sl}$ in Figure 2. If $\Gamma_{sl}$ is equal to the initial complex frequency, then $p_{sl}$ does not change value. For the case where $\Gamma_{sl}=0$, and $p_{sl}\neq 0$, there is no stable solution. For all other cases:

[TABLE]

where $\Re(\Gamma_{sl})$ denotes the real part of $\Gamma_{sl}$, $\Im(\Gamma_{sl})$ denotes the imaginary part of $\Gamma_{sl}$, and $\pm$ corresponds to the sign of $\Re(\Gamma_{sl})$. Therefore, no solution exists where $\lim_{t\to\infty}\gamma_{sl}\leq 0$ except for the case where $\gamma_{sl}$ is temporally static. ∎

By integrating equation (26), one can rewrite equation (22) as the following:

[TABLE]

We note that for the case where $\varphi\to\infty$, $p_{sl}(t)\to p_{sl}(0)$ and therefore:

[TABLE]

We also impose the following physical limit:

[TABLE]

such that $\lim_{t\to\infty}|\phi|<\infty$. However, finding solutions that satisfy this limit is not trivial. One can decompose equation (24) into two real equations:

[TABLE]

These equations can be linearised from a time $t=0$ where $\mathop{}\!\mathrm{d}\gamma_{sl}/\mathop{}\!\mathrm{d}t=0$. From equation (30a):

[TABLE]

where we have implicitly discarded some nonlinear terms $\mathcal{O}[(f-F)^{2}]$, and where the dimensionless constant $\lambda^{2}_{sl}$ is given by:

[TABLE]

As such, one can see that at $t=0$, for $\gamma_{sl}\ll\omega_{\textrm{pl.}}$ the frequency is approximately equal to the plasma frequency.

V Basis expansion of the Vlasov equation

The Vlasov equation is given byvlasov1938vibration :

[TABLE]

We begin by stating:

[TABLE]

This new variable $\epsilon$ is defined everywhere. Let us further impose that $f(\epsilon)$ is analytic for $\epsilon\in\mathbb{R}$. Then:

[TABLE]

Ignoring the trivial solution $\mathop{}\!\mathrm{d}f/\mathop{}\!\mathrm{d}\epsilon=0$, we find that $\epsilon$ also satisfies the Vlasov equation:

[TABLE]

By using this approach, one can represent $f$ as a functional of $\epsilon$, where $\epsilon$ can be chosen to have some physical significance. We use a basis decomposition in $v$ of $\epsilon$:

[TABLE]

This form of the basis decomposition is key, and corresponds to a different basis decomposition at every point in $x$ and $t$. By allowing each $(x,t)$ to permit a different basis decomposition, we allow for nonlinear coupling which is typically not allowed using a variables separable method of the form $f=XT(x,t)V(v)$.

To satisfy equations (2) and (3), we desire that $\{g_{j}(x,v,t)\}$ are square-integrable functions on $v\in\mathbb{R}$, and to do so we make the constraint that $\{g_{j}(x,v,t)\}$ are compactly supported. Accordingly, $\{g_{j}(x,v,t)\}$ also vanish at infinity:

[TABLE]

Then, if we take the $0^{\textrm{th}}$ moment of the reduced Vlasov equation:

[TABLE]

Using our basis expansion:

[TABLE]

Next, we perform a co-moving Galilean transform to $g(x,v,t)=g_{j}(x,v-u_{j}(t),t)$, where $u_{j}$ is at this point an undetermined function. Then:

[TABLE]

If $g_{j}$ is symmetric about $v-u_{j}(t)=0$, then the second integral vanishes due to antisymmetry. Then:

[TABLE]

While it appears that this is a restriction on the solutions that are allowed, in fact any $\epsilon(x,v,t)$ can be represented via such a decomposition. The form of the basis functions has been restricted, but the function it represents has not (c.f. Fourier expansion with phase angles versus Fourier expansion with sines and cosines only).

Therefore from equation (37), one finds:

[TABLE]

This takes a similar form to a fluid advection equation with a source term.

At this point, we define an integral mapping $x,t\longmapsto\chi$, where $\chi$ is a matrix of coordinates in instantaneous comoving frames:

[TABLE]

where $\chi_{jl}/k_{l}$ represents a Galilean transform to the comoving frame with variable wave speed $u_{j}\equiv\omega_{jl}(t)/k_{l}$. In this sense, $\chi$ contains all of the transformations possible with the available frequencies $\{\omega_{j}\}$ and available wavenumbers $\{k_{l}\}$ in the system.

Our basis decomposition is therefore in such a form that the symmetry point in $v$ of the even functions $\{g_{j}\}$ is the variable wave speed of a wave with frequency $\omega_{jl}$ and wavenumber $k_{l}$. Therefore, one can intuitively posit that such a decomposition yields phase-space structures which move with the frequency of some form of wave in the system. If these waves are chosen to be electrostatic waves in the system, then one expects that the migration of phase-space structures corresponds to nonlinear frequency sweeping moving the wave-particle resonances.

Then, from equation (39):

[TABLE]

This forms a non-holonomic constraint on the system. This result is particularly interesting as it implies that it is not possible for one to set up a suitable Lagrangian for the system with $c_{j}M_{0}[g_{j}]$ as generalised fields. It is therefore our belief that this formalism cannot be derived from an action-angle formulation. In such a formulation, Hamilton’s principle assumes only holonomic constraints on the system. In contrast, our formalism yields an analog of d’Alembert’s principle, allowing for more exotic scenarios involving non-conservative forces.

We have shown that $\epsilon$ solves the $0^{\textrm{th}}$ moment of the Vlasov equation; by re-examining the Vlasov equation, we can now explicitly represent the mapping between the electric potential $\phi$ and $\epsilon$.

VI Effective potential

VI.1 Motivation from BGK theory

We expect to recover the BGK solution bernstein1957exact near a given resonance $v=u$:

[TABLE]

where $\phi(x,t)$ is given by:

[TABLE]

and where $u\equiv\omega_{l}/k_{l}$ for all frequencies $\{\omega_{l}\}$ and all wavenumbers $\{k_{l}\}$. The BGK solution is only for a single phase velocity. Therefore, we Taylor expand $\epsilon$ around $v=u$:

[TABLE]

By inspection, we find the possible approximation:

[TABLE]

With a form satisfying this being a Gaussian:

[TABLE]

This is not a particularly robust solution: the function is not analytic in $\chi$ at $\phi(x,t)=0$, Gaussian basis functions are not orthonormal, and $g_{j}$ does not have the correct growth rate. However, every moment of this function is analytic in $x$ if $\phi(x)\leq 0$:

[TABLE]

Furthermore, in this approximate form, particles enjoy an effective potential in phase-space which is dependent on their velocity. Particles far from the resonance are ‘screened’ in a manner similar to Debye shielding, however this phenomenon occurs in momentum space, not real space.

To illustrate this further, one can represent $\epsilon$ in the form:

[TABLE]

where $\epsilon_{0}=q\phi+\frac{1}{2}m(v-u)^{2}$ is the particle energy. As particles follow phase-space orbits given by contours of constant $f$, the orbit is given by an orbit where the potential energy have effectively been modified by terms of $\mathcal{O}(v-u)^{3}$.

These extra terms allow for nonlinear orbits where particles exchange differing amounts of energy with the wave depending on how close they are to the resonance; this is in better accordance with the Landau damping landau1946vibrations .

VI.2 Distribution function construction

Lemma 1**.**

The distribution function $f$ can be represented as a functional in the distributive form:

[TABLE]

where $\epsilon_{0}$ is an arbitrary function, and $\Delta\epsilon=\epsilon-\epsilon_{0}$.

Proof.

First, we state that the vector \vec{\mu}\equiv(\begin{array}[]{c c c}x&v&t\end{array}). We extend our statement that $f$ is analytic in $\epsilon$ to instead require that $f$ permits a functional Taylor expansion of $f$ about a function $\epsilon_{0}$: dreizler2013density

[TABLE]

where $\Omega\in\mathbb{R}^{D}$ is a D-dimensional closed volume in which $d\vec{\mu}_{j}$ is the infinitessimal volume element, $\Delta\epsilon_{j}\equiv\Delta\epsilon(\vec{\mu}_{j})$ and $(\delta/\delta\epsilon)$ denotes the functional derivative with respect to $\epsilon$, given by:

[TABLE]

where $\varepsilon\eta(\mu)$ is equivalent to $\Delta\epsilon(\mu)$. It is then possible for us to express equation (48) in the form:

[TABLE]

where $\mathcal{C}_{l}$ denotes a functional of many functions. Each function has a different domain, but an identical codomain. This codomain is the Banach space defined by $\epsilon_{0}$. Therefore by inspection, the nonlocal functional $\Delta f$ is given by:

[TABLE]

∎

Then, it is possible to show that our perturbation to the distribution function $\delta f$ can be represented using an entirely different basis than that which is employed with $\epsilon_{0}$.

Lemma 2**.**

If $\exists\epsilon_{0}:\partial f[\epsilon_{0}]/\partial t=0$, then the following Dirichlet boundary condition must be satisfied:

[TABLE]

Proof.

The functional form of $\epsilon_{0}$ is solely a function of $x$ and $v$. Therefore, it is possible to enforce the stronger constraint:

[TABLE]

Accordingly, as $f$ is analytic in $t$, it permits a Maclaurin expansion in $t$. Expanding $f$ as such gives:

[TABLE]

By setting $t=0$, the aforementioned Dirichlet boundary condition is found. ∎

It is therefore possible to construct a ‘first-order’ distribution function by using Lemmas 1 and 2. By splitting $\Delta f$ into a linear and nonlinear contribution:

[TABLE]

By discarding the terms which nonlinearly depend on $\Delta\epsilon$, one can approximate the distribution function. This generalises the technique used by Bernstein, Greene and Kruskal to obtain ‘first-order’ distribution functions for so-called BGK modes.bernstein1957exact

It is however also worth noting that under the approximation that $\Delta f$ is a linear functional, the functional is distributive under addition:

[TABLE]

VI.3 Gaussian expansion

Hermite functions form an orthonormal basis, and have recently been explored in the literature by J. M. Heninger et al. for G-transform based solutions of the linearised Vlasov-Poisson system. heninger2018integral Gaussian functions have a similar form to Hermite functions, but are easier to work with due to the lack of preceding Hermite polynomials in the function definition. In addition, the previous subsection motivates the use of Gaussian functions via the ability to reduce to BGK theory.

In this subsection, we will examine solutions of the form $f=f_{(0)}(x,v,t)+\Delta f[\Delta\epsilon]$, under the following Gaussian expansion:

[TABLE]

where $v_{N,l}(\chi,v,t)$ is a normalising function with units of velocity, and $W_{l}(\chi,t)$ is an enveloping function with units of energy.

By inserting $\epsilon$ into equation (33) one finds that to solve the Vlasov equation, $v_{N,l}$ is an approximate solution of the following differential equation under $x,v,t$ coordinates:

[TABLE]

For the case of a single frequency wave in the system, (53) is exact. Otherwise, we have assumed that the overlap integral of two Gaussians is approximately zero. As such, for resonant structures which are close to each other in phase-space, the validity of (53) breaks down.

It is possible to expand the left hand side of equation (53) using a Laurent series:

[TABLE]

where $\{L_{l}^{(j)}\}$ are Laurent coefficients. Accordingly, one finds that for the analytic part of equation (53), in the limit that $v\to u_{l}$:

[TABLE]

From this point onward, we will seek solutions where $L_{l}^{(0)}$ is the only term in the Laurent expansion.

Lemma 3**.**

If the Laurent series expansion of the left hand side of equation (53) only has one term, $L_{l}^{(0)}(v-u_{l})^{0}$, then, equation (53) is approximately solved with:

[TABLE]

Proof.

Using the expansion given by equation (54), equation (53) takes the form:

[TABLE]

By substituting equation (55) the lemma is proved if the trivial solution of $(v-u_{l})\neq 0$ is ignored. ∎

VI.4 Example solutions

Here, we will examine the possible solution:

[TABLE]

where $\phi=\sum\limits_{l}\phi_{l}(x,t)$. If each $\phi_{l}$ has the form (from equation (22)):

[TABLE]

where $\theta_{l}$ is some initial phase, then from Lemma 3:

[TABLE]

We desire $\Delta\epsilon$ to be smooth, finite, and real-valued everywhere, which in turn requires $0\leq v_{N,l}<\infty$. To enforce $|v_{N,j}|<\infty$ everywhere:

[TABLE]

where $n$ is an integer. This can be satisfied in two ways. Equation (57) can be satisfied if:

[TABLE]

However, one cannot guarantee that this is true for each and every $n$. Therefore, we use the alternative solution:

[TABLE]

such that $\tilde{W}_{l}$ is the amplitude of the wave with phase velocity $u_{l}$. Therefore, by examining the rate of change of $u_{l}$, one finds the equations:

[TABLE]

The system of equations (58) must be consistent with the solution for $\{p_{sl}(t)\}$ given by (26), determining the evolution of $\phi$.

VI.4.1 Stationary frequency

For the case of a single wave in the system, equation (58) yields a stationary frequency wave. This exact solution is consistent with equation (24) if the growth rate is also static:

[TABLE]

Therefore, we find:

[TABLE]

This solution therefore gives to lowest order in $(v-u_{l})$:

[TABLE]

This gives an estimate for the separatrix width a factor of $\sqrt{2}$ higher, as the particle acts as though it is trapped in a potential well twice as large as the real potential.

In contrast, if one chose to use $W_{l}\to-W_{l}$ instead, $v_{N,l}^{2}$ would be negative, allowing for us to recover the BGK solution. In such a case, for the full Gaussian representation $\Delta\epsilon$ does not vanish as $v\to\infty$ (as is also the case for BGK modes), and moments of the $\epsilon$ are not real valued (as $v_{N,l}$ is imaginary).

In Figure 1, we show the corresponding full solution for $\Delta\epsilon(x,v,t=0)$ under this approximation, and how it deviates from BGK theory away from resonance. Here, we use normalised values of $|\phi_{ll}|=q=m=1$, and $k=2\pi/L$.

VI.4.2 Sweeping frequency

For the case of two waves in the system with the same wavenumber:

[TABLE]

where $\tilde{\omega}\equiv\omega_{+}-\omega_{-}$, and $\tilde{\theta}\equiv\theta_{+}-\theta_{-}$. If one examines the case where:

[TABLE]

such that the frequency sweep is roughly symmetric, then $\tilde{W}_{+}\approx-\tilde{W}_{-}$. This type of mode therefore shares similarities with symmetrically sweeping modes as is observed in experiments on devices such as the Mega Amp Spherical Tokamak (MAST) gryaznevich2006perturbative . For either of the waves to be linearly stable, we require $\exists\,t_{0}:\gamma_{\pm}(t=t_{0})=0$. However, the waves can still be nonlinearly unstable. To assess the nonlinear stability, we combine equations (30a) and (30b):

[TABLE]

where $\dot{\omega}=\pm\mathop{}\!\mathrm{d}\tilde{\omega}/\mathop{}\!\mathrm{d}t$. For nonlinear stability, we require that $\mathop{}\!\mathrm{d}\gamma_{\pm}/\mathop{}\!\mathrm{d}t<0$. Therefore:

[TABLE]

Additionally, for stability we require that $\lim_{t\to\infty}\dot{\omega}=0$. Therefore, for nonlinear stability the following inequality must be satisfied:

[TABLE]

where one can identify that $\Gamma_{\pm}^{2}\equiv\lim_{t\to\infty}p_{\pm}^{2}$. In Figure 3, we sketch the regions defined by (61) with the corresponding sign of $\mathop{}\!\mathrm{d}\gamma_{\pm}/\mathop{}\!\mathrm{d}t$. Intersections between the line $\gamma_{\pm}$ and $\lim_{t\to\infty}\gamma_{\pm}$ show stable values for the growth rate, corresponding to nonlinearly stable (metastable limit with $\gamma_{\pm}<0$) and nonlinearly unstable (stable limit with $\gamma_{\pm}>0$) states.

Everywhere else above the line defined by the negative root of (61), the wave alternates between increasing and decreasing growth rates. As the choice of initial conditions influences $\Gamma_{\pm}$, some initial conditions will not have a limit for $\gamma_{\pm}$. This allows for either repeated chirping as observed in simulations and experiments vann2003fully ; vanzeeland2019alfven , or nonlinear instability (if $\int_{t_{1}}^{t_{2}}\gamma_{\pm}\mathop{}\!\mathrm{d}t>0$ for $t_{2}>t_{1}$). This may be a candidate for the rapid frequency chirping observed during abrupt large events or mode avalanching in tokamaks sharapov2002alfven ; gryaznevich2006perturbative ; podesta2009experimental .

Between $\gamma_{\pm}=0$ and the line defined by the negative root of (61), the growth rate is negative, but decreasing in amplitude. Waves here are linearly stable, but may become nonlinearly unstable if they are able to cross the $\gamma_{\pm}=0$ line with finite amplitude. Below the line defined by the negative root of (61), the growth rate is negative and always decreasing, corresponding to both linear and nonlinear stability.

From equation (60), we find that the magnitude of the sweeping rate is directly proportion to $\tilde{A}$. Therefore, if there is a large difference in the wave amplitudes, the system here is nonlinearly unstable.

It is worth noting that this solution is only approximately valid for structures which have a very small overlap integral in phase-space. As such, one can consider these long-range, coupled structures in a approximate fashion by evolving a continuum of superposed, stationary frequency solutions.

VI.5 Separatrix

We define the separatrix as the largest closed contour in phase-space. Here, we examine only a single value of $l$ active in the system. Here, we examine $\Delta\epsilon(\chi,v,t)=c$, an unspecified constant. Points on such a contour only exist for:

[TABLE]

as $\exp(-z)\in(0,1]$ for positive $z$. As such, if $c<W_{l}$, no point exists on the contour for the corresponding values of $\chi,t$. Therefore, closed contours exist for $c>\textrm{min}(W_{l})$, and the separatrix is given by $c=0$. At this fixed value of $\Delta\epsilon(\chi,v,t)$:

[TABLE]

As a result, all the particles here are ‘trapped’. However in reality, weakly bound particles would scatter out of the potential via neoclassical transport in a model that considers collisions.

In light of this, we use the following fit for an effective separatrix to compare to BGK theory:

[TABLE]

which is given when the structure is at $1/\textrm{e}$ height. In comparison to BGK theory, our value of $W_{l}$ has twice the amplitude of $\phi_{l}$. Therefore, the width of the separatrix given by this theory is roughly a factor of $\sqrt{2}$ larger than that which is expected for a BGK island.

VII Conclusions

In conclusion, we have solved a variety of problems related to 1D kinetic plasmas, allowing for a detailed description of frequency chirping and basis expansion of the distribution function.

We have shown that a Berk-Breizman sink that linearly dissipates electric field energy from a 1D kinetic plasma can be approximated as a Krook-like collision operator. We have shown in Section IV that it is possible to represent the time-varying frequency of waves in the system via the corresponding growth rates. It was therefore shown in Section VI that it is possible to analytically produce models that account for chirping and nonlinear saturation.

We have shown that the Vlasov equation can be exactly solved for the general case of varying complex frequency waves, and show that a Gaussian basis expansion allows one to find an exact solution for a single wave in the system, and to recover approximate solutions.

We find that under a Gaussian expansion for $\Delta\epsilon$, the separatrix width in velocity for a single mode is a factor of $\sqrt{2}$ larger than the width expected for a BGK island. The family of modes we derive therefore exhibit highly nonlinear orbits, where passing particles for a BGK wave are here interacting with the wave. This implies that mode overlap occurs for modes which are further apart in phase velocity to one another than arises from BGK theory. However, the modes found here may not typically be excited, the island width observed in simulations is typically that which is expected from BGK theory berk1997spontaneous ; meng2018resonance .

While not shown here, a logical next step would be to use this formalism to investigate the full Laurent expansion in equation (54), allowing for a greater range of acceptable solutions. Alternatively, using a Hermite expansion for $\epsilon$ would allow for one to find an exact form of equation (53) which allows one to examine frequency bifurcation.

VIII Acknowledgements

The author was funded by the EPSRC Centre for Doctoral Training in Science and Technology of Fusion Energy grant EP/L01663X. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The author would like to thank V. N. Duarte, Z. Lin and H. Hezaveh for useful discussions. The author would also like to thank K. Imada and R. G. L. Vann for valuable help with improving this manuscript.