1. Introduction
Gyrokinetic modelling is the most advanced and first principles-based method for predicting turbulent transport driven by microinstabilities in magnetic fusion devices (Brizard & Hahm Reference Brizard and Hahm2007; Cary & Brizard Reference Cary and Brizard2009). Although nonlinear simulations can accurately predict particle, momentum and heat transport (Bourdelle et al. Reference Bourdelle, Citrin, Baiocchi, Casati, Cottier, Garbet, Imbeaux and Contributors2015), they are often too computationally expensive for integrated modelling purposes. The cost of any given nonlinear simulation is typically of the order of
$10^4$
CPU hours to
$10^5$
CPU hours at a single radial point (Citrin et al. Reference Citrin2017). Integrated modelling frameworks, meanwhile, require thousands of flux calculations for every second of a plasma discharge in a large magnetic confinement fusion device. Therefore, it is computationally advantageous to approximate the turbulent transport from local linear simulations and make use of a mixing length estimate (Casati et al. Reference Casati2009; Bourdelle Reference Bourdelle2015; Staebler et al. Reference Staebler, Bourdelle, Citrin and Waltz2024). Quasilinear models, such as TGLF and QuaLiKiz, heavily approximate the linearized gyrokinetic equations to further reduce the computational cost of any given simulation (Waltz et al. Reference Waltz, Staebler, Dorland, Hammett, Kotschenreuther and Konings1997; Staebler, Kinsey & Waltz Reference Staebler, Kinsey and Waltz2007; Citrin et al. Reference Citrin2017; Stephens et al. Reference Stephens, Garbet, Citrin, Bourdelle, van de Plassche and Jenko2021). These models typically formulate the gyrokinetic equations as an eigenvalue problem instead of an initial value problem. This approach allows for systematic reductions of the resulting equations. Moreover, eigenvalue codes are often necessary since accurate quasilinear models can require the solutions of multiple eigenmodes rather than just that of the most dominant instability (Pueschel et al. Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016). In contrast, initial value codes can only obtain the dominant instability.
One popular way of formulating the gyrokinetic eigenvalue problem is with the use of a Ritz variational principle, which shares some similarities with the action principle. In previous work, the variational principle has also been used to analyse other modes in plasma physics, such as drift-tearing modes and drift waves (Hazeltine & Ross Reference Hazeltine and Ross1978; Ross & Mahajan Reference Ross and Mahajan1978; Mahajan et al. Reference Mahajan, Hazeltine, Strauss and Ross1979). Some example applications of this principle to gyrokinetics can be found in Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990), Bourdelle et al. (Reference Bourdelle, Garbet, Hoang, Ongena and Budny2002), Citrin et al. (Reference Citrin2017), Hamed et al. (Reference Hamed, Muraglia, Camenen, Garbet and Agullo2019), Kotschenreuther et al. (Reference Kotschenreuther, Liu, Mahajan, Hatch and Merlo2024) and Morren et al. (Reference Morren, Proll, van Dijk and Pueschel2024). The idea is as follows: suppose we have an eigenvalue problem of the form
$\mathcal{L}(\omega ) \chi = 0$
, where
$\omega = \omega _r + i \gamma$
is a complex eigenvalue where
$\omega _r$
and
$\gamma$
are real,
$\mathcal{L}$
is an integrodifferential operator and
$\chi$
is a complex-valued field that represents the eigenmode. One then left-multiplies by
$\chi ^{\ast }$
and integrates over the entire domain with a suitable weighting function to obtain a weak form of the problem. For gyrokinetics, one usually solves for the distribution function in terms of the perturbed electromagnetic potentials and fields and then substitutes it into the gyrokinetic field equations. In practice,
$\chi$
is either guessed from other methods or taken from a gyrokinetic simulation, though in principle one can solve for
$\chi$
self-consistently with the use of trial functions. The mode is considered unstable if the growth rate
$\gamma$
is positive; otherwise it is considered stable.
Unfortunately, for the gyrokinetic field equations, the above method is mathematically incorrect. Although the above specific formulation of the principle was likely taken either from Ritz methods in quantum mechanics or the action principle in field theory, a closer examination of nonlinear eigenvalue problems reveals that it is simply incorrect in gyrokinetics. The fatal flaw is that the gyrokinetic field equations are a nonlinear eigenvalue problem (Voss Reference Voss2013), while the above method was formulated for linear Hermitian eigenvalue problems. The gyrokinetic field equations are not linear in the eigenvalue
$\omega$
due to the Landau resonance. In addition, gyrokinetic field equations are not Hermitian because one needs to analytically continue the integral equations with the use of a Landau contour for damped eigenmodes, i.e. modes where the eigenvalues have negative imaginary part. (Note that in nonlinear eigenvalue problems, we take ‘Hermitian’ to mean that
$\hat {L}(\omega )^{\dagger } = \hat {L}(\omega ^{\ast })$
, where
$\dagger$
denotes the conjugate transpose and the asterisk denotes the scalar complex conjugate.)
This work has two key goals. The first is to derive and state a mathematically rigorous variational principle for the collisionless local linear gyrokinetic field equations. In general, the rigorous principle for nonlinear non-Hermitian eigenvalue matrix problems requires the use of a left-eigenfunction that solves the adjoint problem (Voss Reference Voss2013). Fortunately, although the gyrokinetic field equations are not Hermitian, they are indeed complex symmetric, meaning that
$\hat {L}(\omega ) = \hat {L}^T (\omega )$
, where
$T$
denotes the ordinary transpose. Remarkably, we can then rescue the above incorrect variational principle with a slight modification. Instead of left-multiplying by
$\chi ^{\ast }$
, one merely needs to left-multiply by
$\chi$
without taking the complex conjugate and then proceed as usual. We present a heuristic proof of this complex symmetry for the collisionless gyrokinetic field equations, as well as a brief overview of nonlinear eigenvalue problems in § 2.
The second goal is to derive the specific equations to be used in the variational principle. We rigorously prove the complex symmetry by inspection and compactly represent a weak form for the linearized field equations in general geometry and with electromagnetic effects in the local collisionless limit. Not only is the variational principle now formulated correctly, but it also allows for generic trial functions with linear and nonlinear undetermined coefficients. One can then use a finite element method with adaptive mesh refinement, for example, to solve the gyrokinetic field equations, although we do not do so here. In deriving the equations, we also obtain physical insight into the system via the connection between the gyrokinetic field equations and the use of action angle variables in guiding centre orbit theory (Stephens, Garbet & Jenko Reference Stephens, Garbet and Jenko2020). A weak form of the field equations is rigorously derived in § 4, and we also present a brief analysis of guiding centre orbits in Appendix A with special attention paid to tokamak and stellarator geometries. Also, although we include the parallel perturbed magnetic field
$\delta B_{\parallel }$
in § 2, for the sake of brevity we do not include it in § 4. Extending the derived equations to include
$\delta B_{\parallel }$
is straightforward.
We note that the general methods and integral transforms used to derive the equations are not new. Many of the techniques can be found in Rewoldt, Tang & Chance (Reference Rewoldt, Tang and Chance1982) and Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990). We consider our work an extension of this previous work. First, we have generalized their approach to account for (periodized) stellarator geometry, when their methods were originally designed for tokamak geometry. Second, the integral transforms in Rewoldt et al. (Reference Rewoldt, Tang and Chance1982) have been modified such that the dispersion relation resembles the form presented in Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990). In particular, the drift resonance for passing particles is now present in the resonant denominator with a continuum integral, which should be more convenient for numerical applications. Last, while both Rewoldt et al. (Reference Rewoldt, Tang and Chance1982) and Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990) only allowed for the variation of linear coefficients, our formulation allows for both linear and nonlinear varying parameters with any suitable set of trial functions.
After we discuss the results in § 5, we also present a similar analysis for the linearized gyrokinetic equation proper in Appendix B using the perturbed distribution function. The adjoint state method has been recently applied to the linearized gyrokinetic equations (Acton et al. Reference Acton, Barnes, Newton and Thienpondt2024), where an adjoint equation is derived and a corresponding adjoint solution is obtained via simulation. We briefly derive a connection between the right eigenfunction and the left eigenfunction for the linearized gyrokinetic equation; this may be useful when applying adjoint state methods or if one wishes to use the distribution function to formulate a variational principle.
2. The variational principle for nonlinear eigenvalue problems
2.1. Linear Hermitian eigenvalue problems
To begin, we review the variational principle for linear Hermitian eigenvalue problems. Since this method is often used in introductory quantum mechanics textbooks, we shall use it to build intuition (Sakurai & Napolitano Reference Sakurai and Napolitano2011; Griffiths & Schroeter Reference Griffiths and Schroeter2018).
Let
$H$
be an
$n\times n$
Hermitian matrix (where
$H^\dagger = H$
). The linear Hermitian eigenvalue problem is to find
$\lambda$
such that
\begin{equation} H \boldsymbol{x} = \lambda \boldsymbol{x}, \end{equation}
has a non-trivial solution
$\boldsymbol{x}$
. This problem is linear because it depends only linearly on
$\lambda$
. The adjoint problem is
\begin{equation} \boldsymbol{y}^{\dagger } H = \lambda \boldsymbol{y}^{\dagger }. \end{equation}
A non-trivial solution
$\hat {\boldsymbol{x}}$
of (2.1) is a right eigenvector of
$H$
, while a non-trivial solution
$\hat {\boldsymbol{y}}$
of (2.2) is a left eigenvector of
$H$
. We call a solution pair
$(\hat {\lambda }, \hat {\boldsymbol{x}})$
of (2.1) an eigenpair. Meanwhile, if
$\hat {\boldsymbol{y}}$
and
$\hat {\boldsymbol{x}}$
share an eigenvalue
$\hat {\lambda }$
, then
$(\hat {\boldsymbol{y}}, \hat {\lambda }, \hat {\boldsymbol{x}})$
is an eigentriplet. Because
$H$
is Hermitian, we know that the eigenvalues are real and that the left and right eigenvectors corresponding to the same eigenvalue coincide. Therefore, it suffices only to consider (2.1).
We now describe the variational method for Hermitian linear eigenvalue problems. First, we define the Rayleigh–Ritz quotient
\begin{equation} R(\boldsymbol{x}) = \frac {\boldsymbol{x}^\dagger H \boldsymbol{x}}{\boldsymbol{x}^\dagger \boldsymbol{x}}. \end{equation}
The variational theorem tells us that
\begin{equation} \lambda _{\text{min}} \leqslant R(\boldsymbol{x}) \leqslant \lambda _{\text{max}}. \end{equation}
Moreover, if
$(\hat {\lambda }, \hat {\boldsymbol{x}})$
is an eigenpair, then
$R(\boldsymbol{x})$
is stationary at
$\boldsymbol{x} = \hat {\boldsymbol{x}}$
.
The variational principle is commonly used in quantum mechanics, where the variational principle can be extended for infinite-dimensional systems. Assume we are given a one-dimensional Hamiltonian operator
$\hat {H}$
and complex-valued trial function
$\psi (x; \alpha )$
that is a function of position
$x$
and parameterized by complex numbers
$\alpha$
. We can then estimate the ground state by minimizing
\begin{equation} R(\alpha ) = \frac {\int _{-\infty }^{\infty } \mathrm{d}{x} \psi (x; \alpha )^\ast \hat {H} \psi (x; \alpha )}{\int _{-\infty }^{\infty } \mathrm{d}{x} \psi (x; \alpha )^\ast \psi (x; \alpha )}. \end{equation}
In this case, if we find a solution
$\hat {\alpha }$
such that
\begin{equation} \frac {\partial R}{\partial \alpha } = 0, \end{equation}
and such that
$R(\alpha )$
is minimized, then the estimates for the ground state eigenfunction
$\psi _0(x)$
and eigenvalue
$\lambda _0$
are
\begin{align} \psi _0(x) &\approx \psi (x, \hat {\alpha }),\\[-10pt]\nonumber \end{align}
\begin{align} \lambda _0(x) &\approx R(\hat {\alpha }). \end{align}
Importantly, excited states will also satisfy
$\partial _{\alpha } R = 0$
. Therefore, additional eigenvalues and eigenvectors can be obtained by finding other values of
$\alpha$
that locally extremize
$R$
. Recalling that
$\alpha$
is in general a tuple of complex numbers, one rudimentary way of doing this would be to expand the wavefunction in an orthonormal basis,
\begin{align} \psi (x; \alpha ) & = \sum _{n=1}^{N} \alpha _n \phi _n(x), \end{align}
\begin{align} \int _{-\infty }^{\infty } \mathrm{d}{x} \phi _m(x)^\ast \phi _{n}(x) & = \delta _{mn}, \end{align}
where
$\delta$
is the Kronecker delta. Then, assuming the Hamiltonian is diagonalizable in this basis, one can obtain
$n$
eigenvalues and eigenfunctions by solving
$\partial _{\alpha } R = 0$
.
2.2. Nonlinear eigenvalue problems
We now extend the discussion to nonlinear eigenvalue problems. A more thorough introduction to the subject can be found in Voss (Reference Voss2013). Let
$T(\lambda )$
be an analytic
$n\times n$
matrix function that depends on the complex parameter
$\lambda$
(in general, nonlinearly). The nonlinear eigenvalue problem is to find a non-trivial eigenpair
$(\hat {\lambda }, \hat {\boldsymbol{x}})$
that solves the equation
\begin{equation} T(\lambda ) \boldsymbol{x} = \boldsymbol{0}. \end{equation}
Likewise, the adjoint problem is given by
\begin{equation} \boldsymbol{y}^{\dagger } T(\lambda ) = \boldsymbol{0}. \end{equation}
We note that the linear Hermitian eigenvalue problem is simply a special case of the nonlinear eigenvalue problem with
$T(\lambda ) = H - \lambda I$
. However, nonlinear eigenvalue problems have several peculiar properties in comparison with linear eigenvalue problems. In particular:
-
(i) Eigenvectors corresponding to distinct eigenvalues are not necessarily linearly independent (they are in linear eigenvalue problems);
-
(ii) Left and right eigenvectors are not necessarily orthogonal to one another (they are in linear eigenvalue problems);
-
(iii) Even in finite dimensions, there may exist an infinite number of distinct eigenvalues (there exist only a finite number of distinct eigenvalues in finite-dimensional linear eigenvalue problems).
Moreover, a variational characterization of the nonlinear eigenvalue problem is more complicated than that of the linear Hermitian problem. Because the eigenvalues are in general complex, a min–max characterization may not exist. We can, however, still construct a variational principle. Consider the Rayleigh quotient defined by
\begin{equation} R(\boldsymbol{y}, \lambda , \boldsymbol{x}) = \lambda - \frac {\boldsymbol{y}^{\dagger } T(\lambda ) \boldsymbol{x}}{\boldsymbol{y}^\dagger T'(\lambda ) \boldsymbol{x}}, \end{equation}
where
$T'(\lambda )$
is the derivative of
$T(\lambda )$
with respect to
$\lambda$
. If
$(\hat {\boldsymbol{y}}, \hat {\lambda }, \hat {\boldsymbol{x}})$
is an eigentriplet and
$\hat {\boldsymbol{y}}^\dagger T'(\lambda ) \hat {\boldsymbol{x}} \neq 0$
, then
$R$
is stationary at
$(\hat {\boldsymbol{y}}, \hat {\lambda }, \hat {\boldsymbol{x}})$
.
Unlike the linear Hermitian eigenvalue problem, the nonlinear eigenvalue problem uses a Rayleigh quotient that involves both the left and the right eigenvector. This is because, in general, the left and right eigenvectors no longer coincide. Fortunately, a simplification is possible for specific matrices. First, we say
$T(\lambda )$
is complex symmetric if
\begin{equation} T(\lambda ) = T(\lambda )^T \end{equation}
for every
$\lambda$
in our domain of interest. It then follows that if
$(\hat {\lambda }, \hat {\boldsymbol{x}})$
is an eigenpair, then
$(\hat {\bar {\boldsymbol{x}}}, \hat {\lambda }, \hat {\boldsymbol{x}})$
is an eigentriplet, where
$\bar {\boldsymbol{x}}$
denotes the complex conjugate of
$\boldsymbol{x}$
without taking the transpose. In this case, we can simplify the Rayleigh quotient to
\begin{equation} R(\lambda , \boldsymbol{x}) = \lambda - \frac {\boldsymbol{x}^T T(\lambda ) \boldsymbol{x}}{\boldsymbol{x}^T T'(\lambda ) \boldsymbol{x}}. \end{equation}
Unlike the Hermitian linear problem, we now consider only the transpose of
$\boldsymbol{x}$
and not the conjugate transpose. We are therefore making use of a null product since
$\boldsymbol{x}^T \boldsymbol{x} = 0$
has non-trivial solutions. We note that this is commonly used in non-Hermitian quantum mechanics, where it is called the c-product (Moiseyev Reference Moiseyev2011).
It is important to note that although have formulated a variational principle for finite-dimensional matrix problems, the gyrokinetic equation is a partial differential equation and thus infinite-dimensional. Variational characterizations and min–max theorems of the real eigenvalues of the nonlinear eigenvalue problem for self-adjoint operators have been rigorously proven in the infinite-dimensional setting (Hadeler Reference Hadeler1968; Voss, Werner & Hadeler Reference Voss, Werner and Hadeler1982). However, the present authors are unaware of a rigorous mathematical proof that the variational principle holds for more general infinite-dimensional operators. After discretization, though, the nonlinear eigenvalue problem will be rendered finite-dimensional and the variational principle will then hold. We therefore extend the principle and the definition of the Rayleigh quotient to the infinite-dimensional problem without proof.
To treat infinite-dimensional problems on the real line, we first define an inner product. Given functions
$f(x)$
and
$g(x)$
, we define the inner product
$\left \langle g, f \right \rangle$
as
\begin{equation} \left \langle g, f \right \rangle = \int _{-\infty }^{\infty } {\rm d}{x} w(x) g(x)^{\ast } f(x), \end{equation}
where
$w(x)$
is a positive definite weighting function. We then say that the operator
$\hat {T}(\lambda )$
is complex symmetric if
\begin{equation} \big\langle g, \hat {T}(\lambda ) f\big\rangle = \big\langle f^\ast , \hat {T}(\lambda ) g^{\ast }\big\rangle \end{equation}
for all complex
$\lambda$
in our domain of interest. Rather than computing the full Rayleigh quotient, the variational principle can be stated more simply. Let
$\psi (x; \alpha )$
be a complex-valued function parameterized by complex numbers
$\alpha$
. We define the functional
$p$
as
\begin{equation} p(\lambda , \alpha ) = \big\langle \psi ^\ast , \hat {T}(\lambda ) \psi \big\rangle = \int _{-\infty }^{\infty } {\rm d}{x} w(x) \psi (x; \alpha ) \hat {T}(\lambda ) \psi (x; \alpha ). \end{equation}
Note the absence of an explicit complex conjugate inside the integral. We can obtain estimate right eigenfunctions
$\psi$
, left eigenfunctions
$\psi ^\ast$
and eigenvalues
$\lambda$
by determining values of
$\alpha$
and
$\lambda$
that simultaneously satisfy the following equations:
\begin{align} p &= 0,\\[-8pt]\nonumber \end{align}
\begin{align} \frac {\partial p}{\partial \alpha } & = 0 ,\\[-8pt]\nonumber \end{align}
\begin{align} \frac {\partial p}{\partial \lambda } & \neq 0. \end{align}
Multiple eigenmodes can be obtained by finding multiple value of
$\alpha$
and
$\lambda$
that satisfy the above equations. These are the relevant general equations to be solved. We now seek to prove that the gyrokinetic field equations are indeed complex symmetric.
3. Complex symmetry of the gyrokinetic field equations
When writing down the gyrokinetic equations, we work in the ballooning representation where the perturbed fields
$(\delta \phi , \delta A_{\parallel }, \delta B_{\parallel })$
are functions of the distance
$l$
along the field line. The fields correspond to the electrostatic potential, the parallel vector potential and the parallel magnetic field, respectively. The field equations can be obtained by substituting the perturbed distribution function from the linearized gyrokinetic equation into the quasineutrality equation and Ampere’s law. We shall do this explicitly in § 4, but for now, we simply write the result
\begin{align} \sum _s \frac {Z_{s}^2 e^2 n_{0 s}}{T_{0 s}} \delta \phi & - \int {\rm d}^3{v} \frac { Z_{s}^2 e^2 F_{0 s}}{T_{0 s}} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}}\nonumber\\&\quad \times \left (J_{0 s} \delta \phi - J_{0 s} v_{\parallel } \delta A_{\parallel } + \frac {2 J_{1s}}{k_{\perp } \rho _s} \frac {\mu }{Z_{s} e } \delta B_{\parallel } \right ) = 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {k_{\perp }^2}{\mu _0} \delta A_{\parallel } & - \sum _s \int {\rm d}^3{v} \frac {Z_{s}^2 e^2 F_{0 s}}{T_{0 s}} v_{\parallel } J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}}\nonumber\\&\quad \times \left (J_{0 s} \delta \phi - J_{0 s} v_{\parallel } \delta A_{\parallel } + \frac {2 J_{1s} }{k_{\perp } \rho _s} \frac {\mu }{Z_{s} e} \delta B_{\parallel } \right ) = 0,\\[-10pt]\nonumber \end{align}
\begin{align} & - \frac {\delta B_{\parallel } }{\mu _0} - \sum _s \int {\rm d}^3{v} \frac {Z_{s} e F_{0 s}}{T_{0 s}} \frac {2 J_{1s} }{k_{\perp } \rho _s} \mu \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}}\nonumber\\&\quad \times \left (J_{0 s} \delta \phi - J_{0 s} v_{\parallel } \delta A_{\parallel } + \frac {2 J_{1s}}{k_{\perp } \rho _s} \frac {\mu }{Z_{s} e} \delta B_{\parallel } \right ) = 0. \end{align}
All species-specific quantities are labelled with the subscript
$s$
. The background quantities are the density
$n_0$
, the temperature
$T_{0}$
, the Maxwellian
$F$
and the background magnetic field
$\boldsymbol{B}$
. The constants
$Z$
and
$m$
represent the charge number and mass of each species, where
$e$
is the elementary charge. Next,
$v_\parallel$
and
$\mu = m v_{\perp }^2 / 2 B$
are the velocity parallel to
$\boldsymbol{B}$
and the magnetic moment, respectively, where
$v_{\perp }$
is the speed perpendicular to the magnetic field. The frequencies present are the eigenvalue
$\omega$
, the diamagnetic drift frequency
$\omega _{\ast }$
and the drift frequency
$\omega _{d}$
; meanwhile,
$n$
is the toroidal mode number. The wavenumber has been separated into a perpendicular part
$k_{\perp }$
and a parallel part
$k_{\parallel }$
. The first is simply a function of
$l$
in ballooning space, whereas the second is a differential operator. The Bessel functions
$J_0$
and
$J_1$
represent the effect of gyroaveraging and are evaluated at
$k_{\perp } \rho$
, where
$\rho = m v_{\perp } / Z e B$
is the gyroradius. Lastly,
$\mu _0$
is the magnetic constant. It is to be understood that the inverse of a differential operator is an integral operator; this operator is written explicitly in § 4.
To simplify notation, represent the three fields with the normalized column vector
$\boldsymbol{\chi } = (e \delta \phi / T_r, \delta A_{\parallel } / \rho _r B_r, \delta B_{\parallel } / B_r)$
. The fields are subject to the boundary condition that
$\boldsymbol{\chi } \to \boldsymbol{0}$
as
$\left |l\right | \to \infty$
. Here,
$n_r, T_r, \rho _r, B_r$
are arbitrary reference density, temperature, gyroradius and magnetic field, respectively. It is also useful to define a reference sound speed
$c_r = \sqrt {m_r / T_r}$
and plasma beta
$\beta _r = 2 \mu _0 n_r T_r / B_r^2$
, where
$m_r$
is a reference mass. We define our inner product with the field line integral,
\begin{equation} \left \langle \boldsymbol{\chi }_2, \boldsymbol{\chi }_1 \right \rangle = \int _{-\infty }^{\infty } \frac {{\rm d}{l}}{B} \left (\frac {e^2}{T_r^2}\delta \phi _2^\ast \delta \phi _1 + \frac {1}{\rho _r^2 B_r^2} \delta A_{\parallel 2}^\ast \delta A_{\parallel 1} + \frac {1}{B_r^2} \delta B_{\parallel 2}^\ast \delta B_{\parallel 1}\right ). \end{equation}
The field equations can then be expressed in a more convenient matrix form,
\begin{equation} \hat {L} \boldsymbol{\chi } = \boldsymbol{0}, \end{equation}
where
\begin{equation} \hat {L} = \begin{pmatrix} \hat {L}_{\phi \phi } & \quad \hat {L}_{\phi A} & \quad \hat {L}_{\phi B} \\[3pt] \hat {L}_{A \phi } & \quad \hat {L}_{A A} & \quad \hat {L}_{A B} \\[3pt] \hat {L}_{B \phi } & \quad \hat {L}_{B A} & \quad \hat {L}_{B B} \end{pmatrix} \end{equation}
and
\begin{align} \hat {L}_{\phi \phi } &= \sum _s \frac {Z_s^2 n_{0 s} T_{r}}{n_{r} T_{0 s}} - \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s}} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi A} & = \sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s} } J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s} \frac {v_{\parallel }}{c_r},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi B } & = - \sum _s \int {\rm d}^3{v} \frac {Z_s F_{0 s} T_r}{n_r T_{0 s} } J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} \frac {2 J_{1s}}{k_\perp \rho _s} \frac {\mu B_r}{T_r},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{A \phi } & = - \sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s} }\frac {v_{\parallel }}{c_r} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{A A} & = \frac {2 k_{\perp }^2 \rho _{r}^2}{\beta _r} + \sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s}} \frac {v_{\parallel }}{c_r} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s} \frac {v_{\parallel }}{c_r} ,\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{A B } & = - \sum _s \int {\rm d}^3{v} \frac {Z_s F_{0 s} T_r}{n_r T_{0 s}} \frac {v_{\parallel }}{c_r} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} \frac {2 J_{1s}}{k_{\perp } \rho _s} \frac {\mu B_r}{T_r},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{B \phi } & =- \sum _s \int {\rm d}^3{v} \frac {Z_s F_{0 s} T_r}{n_r T_{0 s} } \frac {2 J_{1s}}{k_\perp \rho _s} \frac {\mu B_r}{T_r} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{B A} & = \sum _s \int {\rm d}^3{v} \frac {Z_s F_{0 s} T_r}{n_r T_{0 s}} \frac {2 J_{1s}}{k_{\perp } \rho _s} \frac {\mu B_r}{T_r} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} \frac {v_{\parallel }}{c_r} J_{0 s} ,\\[-12pt]\nonumber \end{align}
\begin{align} \hat {L}_{B B } & = - \frac {2}{\beta _r} - \sum _s \int {\rm d}^3{v} \frac {F_{0 s} T_r} {n_r T_{0 s}} \frac {2 J_{1s}}{k_\perp \rho _s} \frac {\mu B_r}{T_r} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} \frac {2 J_{1s}}{k_\perp \rho _s} \frac {\mu B_r}{T_r}.\\[-15pt]\nonumber \end{align}
This block matrix form will be convenient for proving the complex symmetry of the gyrokinetic field equations.
To prove complex symmetry, we must first examine the operator
$v_{\parallel } k_{\parallel }$
. The operator
$k_\parallel$
is defined such the energy and magnetic moment are held constant when taking the partial derivative with respect to the field line variable. Therefore,
$v_{\parallel }$
and
$k_{\parallel }$
do not commute, so the whole operator must be examined together. However, the operator
$v_{\parallel } k_{\parallel }$
is effectively antisymmetric when the integrals over velocity space and along the field line are taken into account. To see this, first take
$\sigma$
to represent the sign of
$v_{\parallel }$
and define the energy
$\varepsilon$
as
\begin{equation} \varepsilon = \frac {1}{2} m v_{\parallel }^2 + \mu B. \end{equation}
Recall that
\begin{equation} {\rm d}^3{v} = \sum _{\sigma = \pm 1} \frac {2 \pi B}{m_s^2 \left |v_{\parallel }\right |} {\rm d}{\varepsilon } {\rm d}{\mu }. \end{equation}
Let
$f$
and
$g$
be arbitrary functions of
$(l, \varepsilon , \mu , \sigma )$
that decay away at infinity. Then,
\begin{align} \sum _{\sigma = \pm 1} \int \frac {{\rm d}{l}}{B} \frac {2 \pi B}{m_s^2 \left |v_{\parallel }\right |} {\rm d}{\varepsilon } {\rm d}{\mu } g^\ast v_{\parallel } k_{\parallel } f & = \sum _{\sigma = \pm 1} \int {\rm d}{l} \frac {2 \pi }{m_s^2} {\rm d}{\varepsilon } {\rm d}{\mu } g^\ast \sigma k_{\parallel } f\nonumber \\ & = \sum _{\sigma = \pm 1} \int {\rm d}{l} \frac {2 \pi }{m_s^2} {\rm d}{\varepsilon } {\rm d}{\mu } f (-\sigma k_{\parallel }) g^\ast\nonumber \\ & = \sum _{\sigma = \pm 1} \int \frac {{\rm d}{l}}{B} \frac {2 \pi B}{m_s^2 \left |v_{\parallel }\right |} {\rm d}{\varepsilon } {\rm d}{\mu } f (-v_{\parallel } k_{\parallel }) g^{\ast }, \end{align}
where we integrated by parts and set the boundary terms to zero. This proves that
$v_{\parallel } k_{\parallel }$
is antisymmetric. Note that one can bring
$v_{\parallel } k_{\parallel }$
under the complex conjugate to show that it is also Hermitian, as expected.
To show that
$\hat {L}$
is complex symmetric under our inner product, we can take advantage of the block matrix form,
\begin{equation} \hat {L}^T = \left(\begin{array}{c@{\quad}c@{\quad}c} \hat {L}_{\phi \phi }^T & \hat {L}_{A \phi }^T & \hat {L}_{B \phi }^T \\[5pt] \hat {L}_{\phi A}^T & \hat {L}_{A A}^T & \hat {L}_{B A}^T \\[5pt] \hat {L}_{\phi B}^T & \hat {L}_{A B }^T & \hat {L}_{B B}^T \end{array}\right). \end{equation}
We can check this term by term as follows. First, we note that when taking the transpose, we must reverse the order of all operators,
\begin{equation} \left (\int {\rm d}^3{v} A B C \right )^T = \int {\rm d}^3{v} C^T B^T A^T. \end{equation}
We also note that any function that is independent
$l$
commutes with
$k_{\parallel }$
and all other functions independent of
$l$
. Next, the transpose of the resonant denominator can be computed by simply substituting
$v_{\parallel } k_{\parallel } \to - v_{\parallel } k_{\parallel }$
, since this operator is complex antisymmetric as shown above in (3.19). Lastly, since we are integrating over all of velocity space, we can also simply substitute
$v_{\parallel } \to - v_{\parallel }$
everywhere and the value of the integrals will not change.
We demonstrate this explicitly with two examples. For
$\hat {L}_{\phi \phi }$
, we have
\begin{align} \hat {L}_{\phi \phi }^T &= \left (\sum _s \frac {Z_s^2 n_{0 s} T_{r}}{n_{r} T_{0 s}} - \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s}} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s}\right )^T\nonumber\\[3pt] & = \sum _s \frac {Z_s^2 n_{0 s} T_{r}}{n_{r} T_{0 s}} - \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s}} \left (J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s}\right )^T\nonumber \\[3pt] & = \sum _s \frac {Z_s^2 n_{0 s} T_{r}}{n_{r} T_{0 s}} - \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s}} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega + v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s}\nonumber \\[3pt] & = \sum _s \frac {Z_s^2 n_{0 s} T_{r}}{n_{r} T_{0 s}} - \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s}} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s}\nonumber\\[3pt] & = \hat {L}_{\phi \phi }. \end{align}
For
$\hat {L}_{\phi A}$
, we have
\begin{align} \hat {L}_{\phi A}^T & = \left (\sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s} } J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s} \frac {v_{\parallel }}{c_r}\right )^T \nonumber\\ & = \sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s} } \left (J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s} \frac {v_{\parallel }}{c_r}\right )^T \nonumber\\ & = \sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s} } \frac {v_{\parallel }}{c_r} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega + v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s} \nonumber\\ & = - \sum _s \int {\rm d}^3{v} \frac {Z_s^2 F_{0 s} T_r}{n_r T_{0 s} } \frac {v_{\parallel }}{c_r} J_{0 s} \frac {\omega - n \omega _{\ast s}}{\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s}} J_{0 s} \nonumber\\ & = \hat {L}_{A \phi }. \end{align}
It should now be clear by inspection that
$\hat {L}$
is indeed complex symmetric by simply reversing the order of all functions and operators of
$l$
and taking
$v_{\parallel } \to - v_{\parallel }$
.
We now state the variational principle. Let
$\boldsymbol{\chi }$
be parameterized by complex numbers
$\alpha$
. We define the functional
$p(\omega , \alpha )$
as
\begin{equation} p(\omega , \alpha ) = \int _{-\infty }^{\infty } \frac {{\rm d}{l}}{B} \boldsymbol{\chi }^T \hat {L}(\omega ) \boldsymbol{\chi }. \end{equation}
We then seek to solve the following equations:
\begin{align} p &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {\partial p}{\partial \alpha } &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {\partial p}{\partial \omega } &\neq 0. \end{align}
Doing so will give us the right eigenfunctions
$\boldsymbol{\chi }$
, the left eigenfunctions
$\bar {\boldsymbol{\chi }}$
and the eigenvalue
$\omega$
(here, the overbar represents complex conjugation without taking the transpose).
Note that the only assumptions we have made are (i) the background is Maxwellian, (ii) there is no rotation or background electric field
$\boldsymbol{E}$
, (iii) there are no collisions and (iv) the ballooning representation correctly describes the local instability. Impurities, arbitrary geometry and electromagnetic effects are included in the above general formulation. Including a non-trivial collision operator (e.g. pitch-angle scattering) would require expanding the distribution function in terms of basis functions and coupling them together (Hamed et al. Reference Hamed, Muraglia, Camenen and Garbet2018). On the other hand, the inclusion of
$E \times B$
shear through the frequency
$\gamma _E$
is problematic. The correct analysis of the problem requires time-dependent
$k_\perp$
or magnetic field strengths (depending on the specific representation used), so it is no longer a simple eigenvalue problem (Maeyama et al. Reference Maeyama, Watanabe, Nakata, Nunami, Asahi and Ishizawa2025). Even in the small
$\gamma _E$
limit where one just wants an ‘instantaneous’ growth rate, the formulation would still break down because the extra differential operators proportional to
$\gamma _E$
break the complex symmetry of the gyrokinetic field equations.
Note that we have not explicitly written down the gyrokinetic field equations, but rather neatly represented them to prove complex symmetry. Actually writing down the solutions requires some background knowledge in guiding centre dynamics and a detailed analysis of bounce-transit orbits.
4. A weak form of the gyrokinetic field equations
4.1. Preliminaries
All the guiding centre orbit analysis necessary to write down explicit solutions to the gyrokinetic equation is presented in Appendix A. The methods used here are inspired by Tang, Connor & Hastie (Reference Tang, Connor and Hastie1980), Rewoldt et al. (Reference Rewoldt, Tang and Chance1982), Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990) and Chandran & Schekochihin (Reference Chandran and Schekochihin2024). The background Maxwellian is
\begin{equation} F_{0 s} = \frac {n_{0 s}(\psi )}{\pi ^{3/2} v_{T_{s}}^3} \exp (- \varepsilon / T_{0s}(\psi )), \end{equation}
where
$v_{T_{s}} = \sqrt {2 T_{0s} / m_s}$
. We also work in non-orthogonal Clebsch coordinates where the magnetic field
$\boldsymbol{B}$
can be written as
\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \psi , \end{equation}
where
$\alpha$
is a field line label and
$\psi$
is the poloidal flux normalized by
$2 \pi$
. We also use a field line following coordinate
$\theta$
such that all physical quantities obey
$f(\psi , \theta + 2 \pi , \alpha ) = f(\psi , \theta , \alpha )$
and
$\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta = J_{\psi }^{-1} \gt 0$
, where
$J_{\psi }$
is the Jacobian of the coordinate system. All perturbed fields are subject to the ballooning transformation,
\begin{equation} u(\psi , \theta , \varphi ) = \sum _{p} \hat {u}(\psi , \theta + 2 \pi p) e^{i n S(\psi , \theta + 2 \pi p, \varphi )}, \end{equation}
where
$\hat {u}$
is the slowly varying envelope,
$S$
is the eikonal such that
$\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} S = 0$
and
$n$
is the toroidal mode number. This guarantees that
$u$
is periodic in
$\theta$
. We will work exclusively in the ballooning representation, so we shall drop the hat for convenience. By convention,
$S$
is defined such that
\begin{equation} \boldsymbol{\nabla} S = \boldsymbol{\nabla} \alpha + \theta _0 q' \boldsymbol{\nabla} \psi , \end{equation}
where
$\theta _0$
is the ballooning angle and
$q'$
is the derivative of the saftey factor
$q$
with respect to
$\psi$
. We define the perpendicular wavevector as
\begin{equation} \boldsymbol{k}_{\perp } = n \boldsymbol{\nabla} S. \end{equation}
We also define the pitch angle parameter
$\lambda$
as
\begin{equation} \lambda = \frac {\mu }{\varepsilon } .\end{equation}
Using the coordinates
$(\theta , \varepsilon , \lambda )$
, which is equivalent to using
$(\theta , \varepsilon , \mu )$
, the linearized gyrokinetic equation can then be written as
\begin{equation} \frac {v_{\parallel }}{J_{\psi } B} \frac {\partial h_s}{\partial \theta } - i \left (\omega - n \omega _{\mathrm{d} s} \right ) h_s = - i \frac {Z_s e}{T_{0s}} F_{0 s} \left (\omega - \omega _{\ast s} \right ) J_0 (k_{\perp } \rho _s) \left (\phi - v_{\parallel } A_{\parallel } \right ), \end{equation}
where
\begin{align} J_{\psi } &= \left (\boldsymbol{\nabla} \theta \boldsymbol{\cdot} (\boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \psi ) \right )^{-1}\!,\\[-10pt]\nonumber \end{align}
\begin{align} n \omega _{\mathrm{d} s} & = \boldsymbol{k}_{\perp } \boldsymbol{\cdot} \boldsymbol{v}_{\mathrm{d} s},\\[-10pt]\nonumber \end{align}
\begin{align} n \omega _{\ast s} &= \frac {n T_{s}}{Z_{s} e} \frac {1}{n_{0 s}} \frac {\partial n_{0 s}}{\partial \psi } \left (1 + \eta _s \left (\frac {\varepsilon }{T_{0s}} - \frac {3}{2} \right )\right )\!,\\[-10pt]\nonumber \end{align}
\begin{align} \eta_s &= \left(\frac{1}{T_{0s}}\frac{\partial T_{0s}}{\partial \psi}\right)\!\bigg/\!\!\left(\frac{1}{n_{0s}}\frac{\partial n_{0s}}{\partial \psi}\right)\!,\\[-10pt]\nonumber \end{align}
\begin{align} \rho _s & = \frac {m_s v_{\perp }}{e Z_s B} = \frac {m_s}{e Z_s B} \sqrt {\frac {2 \varepsilon }{m_s}} \sqrt {\lambda B}, \end{align}
where the guiding centre velocity
$\boldsymbol{v}_{\mathrm{d} s}$
is given by
\begin{equation} \boldsymbol{v}_{\mathrm{d} s} = \frac {m_{s}}{Z_{s} e B} \frac {v_{\parallel }^2 + \mu B / m_{s}}{B^2} \boldsymbol{B} \times \boldsymbol{\nabla} B + \frac {m_{s}}{Z_{s} e B} \frac {\mu _0 v_{\parallel }^2}{B^3} \boldsymbol{B} \times \boldsymbol{\nabla} p, \end{equation}
where
$p$
is the total plasma pressure. For convenience, we have dropped
$\delta$
from any perturbed fields. Note that in this coordinate system, the parallel velocity is a function defined by
\begin{equation} v_{\parallel } = \sqrt {\frac {2 \varepsilon }{m}} \sqrt {1 - \lambda B}. \end{equation}
4.2. Passing particles
4.2.1. Passing solution of the gyrokinetic equation
For the passing solution, we use the boundary condition that
$h_s \to 0$
as
$|\theta | \to \infty$
. We write the solution assuming
$\mathrm{Im} (\omega) \gt 0$
(otherwise we need to analytically continue the solution) and assuming
$J_{\psi } \gt 0$
. We also need to consider both signs of
$v_{\parallel }$
separately. Defining
$\theta _{-} = \theta - \theta '$
, the solution is
\begin{equation} h_{s, p, \pm } = \mp i \frac {Z_s e F_{0 s} }{T_{0s}} \left (\omega - n \omega _{\ast s} \right ) \int _{\mp \infty }^{\theta } {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (\frac {\phi }{|v_{\parallel }|} \mp A_{\parallel } \right ) e^{\pm i I^{\theta }_{\theta '}}, \end{equation}
where
$I$
will be defined shortly. Then,
\begin{equation} h_{s, p, +} + h_{s, p, -} = -i \frac {Z_s e F_{0 s}}{T_{0s}} \left (\omega - n \omega _{\ast s} \right ) \int _{-\infty }^{\infty } {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (\frac {\phi }{|v_{\parallel }|} - \textit{sgn}(\theta _{-}) A_{\parallel } \right ) e^{i \textit{sgn}(\theta _{-}) I^{\theta }_{\theta '}}. \end{equation}
Meanwhile,
\begin{equation} h_{s, p, +} - h_{s, p, -} = -i \frac {Z_s e F_{0 s}}{T_{0s}} \left (\omega - n \omega _{\ast s} \right ) \int _{-\infty }^{\infty } {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (\textit{sgn}(\theta _{-}) \frac {\phi }{|v_{\parallel }|} - A_{\parallel } \right ) e^{i \textit{sgn}(\theta _{-}) I^{\theta }_{\theta '}}. \end{equation}
The quantity
$I$
is an integral defined as
\begin{equation} I_{a}^{b} = \int _{a}^{b} {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{\mathrm{d} s}\right ). \end{equation}
4.2.2. Quasineutrality
For quasineutrality, we will need to integrate
\begin{equation} \mathcal{L}_{\phi , s, p} = Z_s e \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \int {\rm d}^3{v} \phi J_{0 s} \left (h_{s, p, +} + h_{s, p, -}\right ), \end{equation}
where
\begin{equation} {\rm d}^3{v} = 2 \pi {\rm d}{\varepsilon } {\rm d}{\lambda } \frac {B}{m_s^2} \frac {\varepsilon }{|v_{\parallel }|} \end{equation}
and the integration limits for passing particles are
\begin{align} 0 &\leqslant \lambda \leqslant \frac {1}{B_{\text{max}}},\\[-10pt]\nonumber \end{align}
\begin{align} 0 &\leqslant \varepsilon \leqslant \infty . \end{align}
For notational convenience, we have defined
${\rm d}^3{v}$
such that we only integrate over half of the velocity space without an explicit sum over the sign of the parallel velocity. We will sum over both signs of the parallel velocity separately. It is useful to perform a change of variables using the action angle
$\alpha _t$
,
\begin{equation} \frac {{\rm d}{\theta } J_{\psi } B}{|v_{\parallel }|} = \frac {{\rm d}{\alpha _t}}{\Omega _t}, \end{equation}
where
$\Omega _t$
is the transit frequency. We also note the following:
\begin{equation} I_{0}^{\theta } = \left (\frac {\omega - n \Omega _{d}}{\Omega _t }\right ) \alpha _t - n \Lambda _{\mathrm{d} s}(\alpha _t), \end{equation}
where
$\Omega _{d}$
is the transit averaged drift frequency and
$\Lambda _{d}$
contains the drift deviations, defined in (A.57) and (A.60), respectively. We then define the following Fourier transform pairs:
\begin{align} \phi ^{+}_{s}(k) &= \int _{-\infty }^{\infty } \frac {{\rm d}{\alpha _t}}{2 \pi } J_{0 s} \phi e^{-i k \alpha _t + i n \Lambda _{\mathrm{d} s}(\alpha _t)}, \end{align}
\begin{align} \phi ^{-}_{s}(k) &= \int _{-\infty }^{\infty } \frac {{\rm d}{\alpha _t}}{2 \pi } J_{0 s} \phi e^{i k \alpha _t - i n \Lambda _{\mathrm{d} s}(\alpha _t)}, \end{align}
\begin{align} A^{+}_{s}(k) &= \int _{-\infty }^{\infty } \frac {{\rm d}{\alpha _t}}{2 \pi } J_{0 s} |v_{\parallel }| A_{\parallel } e^{-i k \alpha _t + i n \Lambda _{\mathrm{d} s}(\alpha _t)}, \end{align}
\begin{align} A^{-}_{s}(k) &= -\int _{-\infty }^{\infty } \frac {{\rm d}{\alpha _t}}{2 \pi } J_{0 s} |v_{\parallel }| A_{\parallel } e^{i k \alpha _t - i n \Lambda _{\mathrm{d} s}(\alpha _t)}, \end{align}
\begin{align} J_{0 s} \phi e^{i n \Lambda _{\mathrm{d} s}(\alpha _t)} & = \int _{-\infty }^{\infty } {\rm d}{k} e^{i k \alpha _t} \phi ^{+}_{s}(k), \end{align}
\begin{align} J_{0 s} \phi e^{-i n \Lambda _{\mathrm{d} s}(\alpha _t)} & = \int _{-\infty }^{\infty } {\rm d}{k} e^{-i k \alpha _t} \phi ^{-}_{s}(k), \end{align}
\begin{align} J_{0 s} |v_{\parallel }| A_{\parallel } e^{i n \Lambda _{\mathrm{d} s}(\alpha _t)} & = \int _{-\infty }^{\infty } {\rm d}{k} e^{i k \alpha _t} A^{+}_{s}(k), \end{align}
\begin{align} J_{0 s} |v_{\parallel }| A_{\parallel } e^{-i n \Lambda _{\mathrm{d} s}(\alpha _t)} & = -\int _{-\infty }^{\infty } {\rm d}{k} e^{-i k \alpha _t} A^{+}_{s}(k). \end{align}
We then obtain the following identities:
\begin{equation} \begin{aligned} &\int _{-\infty }^{\infty } {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} \int _{-\infty }^{\infty } {\rm d}{\theta '} \frac {J_{\psi }' B'}{|v_{\parallel }|'} J_{0 s} \phi J_{0 s}' \phi ' e^{i \textit{sgn}(\theta _{-}) I_{\theta '}^{\theta }} \\[6pt] &\qquad\qquad\qquad\ = \frac {4 \pi i}{\Omega _{t s}} \int _{-\infty }^{\infty } {\rm d}{k} \frac {\phi ^{+}_{s}(k) \phi ^{-}_{s}(k)}{\omega - k \Omega _{t s} - n \Omega _{\mathrm{d} s}} \end{aligned} \end{equation}
and
\begin{equation} \begin{aligned} &\int _{-\infty }^{\infty } {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} \int _{-\infty }^{\infty } {\rm d}{\theta '} J_{\psi }'B' \textit{sgn}(\theta _{-}) J_{0 s} \phi J_{0 s}' A_{\parallel }' e^{i \textit{sgn}(\theta _{-}) I_{\theta '}^{\theta }} \\[6pt] &\qquad\qquad\qquad\ = \frac {2 \pi i}{\Omega _{t s}} \int _{-\infty }^{\infty } {\rm d}{k} \frac {A^{+}_{s}(k) \phi ^{-}_{s}(k) + A^{-}_{s}(k) \phi ^{+}_{s}(k)}{\omega - k \Omega _{t s} - n \Omega _{\mathrm{d} s}}. \end{aligned} \end{equation}
Substituting this into the passing solution, we find
\begin{align} \mathcal{L}_{\phi , p, s} &= \frac {2 \pi ^2 Z_s^2 e^2}{m_s^2 T_{0s}} \int _{0}^{\infty }\! {\rm d}{\varepsilon } \varepsilon F_{0 s} (\omega - n \omega _{\ast s}) \int _0^{1 / B_{\text{max}}} \frac {{\rm d}{\lambda }}{\Omega _{t s}} \int _{-\infty }^{\infty } {\rm d}{k} \frac {2 \phi ^{+}_{s} \phi ^{-}_{s} - A^{+}_{s} \phi ^{-}_{s} - A^{-}_{s} \phi ^{+}_{s}}{\omega - k \Omega _{t s} - n \Omega _{\mathrm{d} s}}. \end{align}
This compactly takes into account both signs of the parallel velocity. As noted in Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990), this form is more numerically advantageous than the original form. In the low
$\varepsilon$
limit, the exponentials in (4.14) are highly oscillatory due to the term proportional to
$\omega / |v_{\parallel }|$
. We have moved this term into the Landau resonance. Meanwhile, the remaining eikonals in the Fourier transform pairs are well behaved as
$\varepsilon \to 0$
, since
$\omega _{d} / |v_{\parallel }| \propto \sqrt {\varepsilon }.$
4.2.3. Ampere’s law
For Ampere’s law, we will need to integrate
\begin{equation} \mathcal{L}_{A, p, s} = Z_s e \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \int {\rm d}^3{v} A_{\parallel } |v_{\parallel }| J_{0 s} \left (h_{s, p, +} - h_{s, p, -}\right ). \end{equation}
The steps are similar to that of quasineutrality, so we just write the final result,
\begin{align} \mathcal{L}_{A, p, s} &= \frac {2 \pi ^2 Z_s^2 e^2}{m_s^2 T_{0s}} \!\int _{0}^{\infty } \!{\rm d}{\varepsilon } \varepsilon F_{0 s} (\omega - n \omega _{\ast s}) \!\int _0^{1 / B_{\text{max}}} \frac {{\rm d}{\lambda }}{\Omega _{t s}} \int _{-\infty }^{\infty } {\rm d}{k} \frac {2 A^{+}_{s} A^{-}_{s} - A^{+}_{s} \phi ^{-}_{s} - A^{-}_{s} \phi ^{+}_{s}}{\omega - k \Omega _{t s} - n \Omega _{\mathrm{d} s}}.\nonumber\\[8pt] \end{align}
4.3. Trapped particles
4.3.1. Trapped solution of the gyrokinetic equation
For the trapped solution, we instead need to use the boundary conditions
\begin{align} h_{s, tr, +} (\theta _1) &= h_{s, tr, -} (\theta _1),\\[-10pt]\nonumber \end{align}
\begin{align} h_{s, tr, +} (\theta _2) &= h_{s, tr, -} (\theta _2), \end{align}
where
$\theta _1$
and
$\theta _2$
are bounce points such that
$\theta _1 \lt \theta \lt \theta _2$
. This leads to
\begin{align} &h_{s, tr, +} + h_{s, tr, -} \nonumber\\[6pt] &\quad = \frac {2 Z_s e F_{0 s}}{T_{0s}} \frac {\omega - n \omega _{\ast s}}{\sin \big(I_{\theta _1}^{\theta _2}\big) } \Bigg [ \int _{\theta _1}^{\theta } {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (\frac {\phi }{|v_{\parallel }|} \cos \big(I_{\theta _1}^{\theta '}\big) \cos \big(I_{\theta }^{\theta _2}\big)\right.\nonumber\\& \qquad\qquad\qquad\qquad\qquad + i A_{\parallel } \sin \big(I_{\theta _1}^{\theta '}\big) \cos \big(I_{\theta }^{\theta _2}\big) \bigg)\nonumber \\&\qquad + \int _{\theta }^{\theta _2} {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (\frac {\phi }{|v_{\parallel }|} \cos \big(I_{\theta _1}^{\theta }\big) \cos \big(I_{\theta '}^{\theta _2}\big) - i A_{\parallel } \cos \big(I_{\theta _1}^{\theta }\big) \sin \big(I_{\theta '}^{\theta _2}\big) \right )\Bigg ], \end{align}
and
\begin{align} &h_{s, tr, +} - h_{s, tr, -} \nonumber\\[6pt] &\quad = \frac {2 Z_s e F_{0 s}}{T_{0s}} \frac {\omega - n \omega _{\ast s}}{\sin \big(I_{\theta _1}^{\theta _2}\big) } \Bigg [ \int _{\theta _1}^{\theta } {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (- \frac {\phi }{|v_{\parallel }|}\cos \big(I_{\theta _1}^{\theta '}\big) \sin \big(I_{\theta }^{\theta _2}\big)\right.\nonumber\\& \qquad\qquad\qquad\qquad\qquad + i A_{\parallel } \sin \big(I_{\theta _1}^{\theta '}\big) \sin \big(I_{\theta }^{\theta _2}\big) \bigg)\nonumber \\&\qquad + \int _{\theta }^{\theta _2} {\rm d}{\theta '} J_{\psi } B J_{0 s} \left (\frac {\phi }{|v_{\parallel }|} \sin \big(I_{\theta _1}^{\theta }\big) \cos \big(I_{\theta '}^{\theta _2}\big) + i A_{\parallel } \sin \big(I_{\theta '}^{\theta _2}\big) \sin \big(I_{\theta }^{\theta _1}\big) \right )\Bigg ]. \end{align}
Note that we need to pick
$\theta _1, \theta _2$
such that
$\theta _1 \leqslant \theta \leqslant \theta _2$
.
4.3.2. Quasineutrality
Integrating over the trapped part of velocity space requires much more work than for passing particles. We define
\begin{equation} {\rm d}^3{v} = 2 \pi {\rm d}{\varepsilon } {\rm d}{\lambda } \frac {B}{m_s^2} \frac {\varepsilon }{|v_{\parallel }|}. \end{equation}
The limits of integration for trapped particles, however, are
\begin{align} \frac {1}{B_{\text{max}}} &\leqslant \lambda \leqslant \frac {1}{B(\theta )},\\[-5pt]\nonumber \end{align}
\begin{align} 0 & \leqslant \varepsilon \leqslant \infty . \end{align}
Thus, the upper limit of the pitch angle parameter depends on the current position along the field line. We can, however, use a trick to make the integration limits more manageable. Let
$\theta _{\text{min}} \leqslant \theta _1 \leqslant \theta _2 \leqslant \theta _{\text{max}}$
, where
$\theta _{\text{min}}$
and
$\theta _{\text{max}}$
correspond to the extent of the bounce well. Then,
\begin{equation} \int _{\theta _{\text{min}}}^{\theta _{\text{max}}} {\rm d}{\theta } \int _{1/ B_{\text{max}}}^{1/ B(\theta )} {\rm d}{\lambda } = \int _{1/B_{\text{max}}}^{1/B_{\text{min}}} {\rm d}{\lambda } \int _{\theta _1}^{\theta _2} {\rm d}{\theta }. \end{equation}
We also note the following identity:
\begin{align} &\int _{\theta _1}^{\theta _2} {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} J_{0 s} \phi \int _{\theta _1}^{\theta } {\rm d}{\theta '} \frac {J_{\psi }' B'}{|v_{\parallel }|'} J_{0 s}' \phi '\cos \big(I_{\theta _1}^{\theta '}\big) \cos \big(I_{\theta }^{\theta _2}\big) \nonumber\\[5pt] &\qquad{} + \int _{\theta _1}^{\theta _2} {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} J_{0 s} \phi \int _{\theta }^{\theta _2} {\rm d}{\theta '} \frac {J_{\psi }' B'}{|v_{\parallel }|'} J_{0 s}' \phi ' \cos \big(I_{\theta _1}^{\theta }\big) \cos \big(I_{\theta '}^{\theta _2}\big) \nonumber\\[9pt] &\quad{} = \frac {\pi \sin \!\big(I_{\theta _1}^{\theta _2}\big) }{\Omega _{b s}} \sum _{n_2 = -\infty }^{\infty } \frac {\left (\phi ^{n_2}_{s}\right )^2}{\omega - n_2 \Omega _{b s}- n \Omega _{\mathrm{d} s}}, \end{align}
where
$\phi ^{n_2}_{s}$
is defined as
\begin{equation} \phi ^{n_2}_{s} = \left \langle J_{0 s} \phi \cos (n \tilde {w}_{\mathrm{d} s} - n_2 \alpha _b)\right \rangle _b. \end{equation}
Here,
$\Omega _{d}$
and
$\tilde {w}_{d}$
are the bounce averaged drift frequency and drift deviations, defined in (A.96)–(A.97). We are essentially computing the bounce harmonics of the electrostatic potential modulated by finite Larmor radius effects via the Bessel function and also finite banana width effects via
$\tilde {w}_{\mathrm{d} s}$
. Dropping the
$s$
subscript momentarily, one can prove this by writing
\begin{align} J_{0} \phi \cos (n \tilde {w}_{d}) &= \phi _0^c + \sum _{n_2=1}^{\infty } 2 \phi _{n_2}^{c}\cos (n_2 \alpha _b) , \end{align}
\begin{align} J_{0} \phi \sin (n \tilde {w}_{d}) &= \sum _{n_2=1}^{\infty } 2 \phi _{n_2}^{s} \sin (n_2 \alpha _b), \end{align}
as well as
\begin{align} \phi ^{n_2} &= \phi _{n_2}^{c} + \phi _{n_2}^{s},\\[-8pt]\nonumber \end{align}
\begin{align} \phi ^{-n_2} &= \phi _{n_2}^{c} - \phi _{n_2}^{s}, \end{align}
and computing the integrals manually. Meanwhile, we define a similar bounce average for the electromagnetic piece,
\begin{equation} A^{n_2}_{s} = \left \langle i J_{0 s} v_{\parallel } A_{\parallel } \sin (n \tilde {w}_{\mathrm{d} s} - n_2 \alpha _b )\right \rangle _b. \end{equation}
(Note the factor of
$i$
.) One can then write (again dropping the
$s$
subscript for clarity)
\begin{align} i J_{0} |v_{\parallel }| A_{\parallel } \sin (n \tilde {w}_{d}) & = A_0^c + \sum _{n_2=1}^{\infty } 2 A_{n_2}^c \cos (n_2 \alpha _b), \end{align}
\begin{align} i J_{0} |v_{\parallel }| A_{\parallel } \cos (n \tilde {w}_{d} ) & = \sum _{n_2 = 1}^{\infty } 2 A_{n_2}^s \sin (n_2 \alpha _b), \end{align}
as well as
\begin{align} A^{n_2} &= A_{n_2}^c - A_{n_2}^s,\\[-10pt]\nonumber \end{align}
\begin{align} A^{-n_2} &= A_{n_2}^c + A_{n_2}^s. \end{align}
This leads to the following identity:
\begin{align} &\int _{\theta _1}^{\theta _2} {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} J_{0 s} \phi \int _{\theta _1}^{\theta } {\rm d}{\theta '} J_{\psi }' B' J_{0 s}' i A_{\parallel }' \sin \big(I_{\theta _1}^{\theta '}\big) \cos \big(I_{\theta }^{\theta _2}\big) \nonumber\\[5pt] &\qquad{} - \int _{\theta _1}^{\theta _2} {\rm d}{\theta } \frac {J_{\psi } B}{|v_{\parallel }|} J_{0 s} \phi \int _{\theta }^{\theta _2} {\rm d}{\theta '} J_{\psi }' B' J_{0 s}' i A_{\parallel }' \cos \big(I_{\theta _1}^{\theta }\big) \sin \big(I_{\theta '}^{\theta _2}\big) \nonumber\\[5pt] &\quad{} =- \frac {\pi \sin \big(I_{\theta _1}^{\theta _2}\big) }{ \Omega _{b s}} \sum _{n_2 = -\infty }^{\infty } \frac {\phi ^{n_2}_{s} A^{n_2}_{s}}{\omega - n_2 \Omega _{b s}- n \Omega _{\mathrm{d} s}}. \end{align}
We can finally write down the trapped part of the integral we need to compute. First, define
\begin{equation} \mathcal{L}_{\phi , tr, s} = Z_s e \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \int {\rm d}^3{v} \phi J_{0 s} \left (h_{s, tr, +} + h_{s, p, -}\right ). \end{equation}
Then, we find
\begin{align} \mathcal{L}_{\phi , tr, s} &= \sum _p \sum _{n_2} \frac {4 \pi ^2 Z_s^2 e^2}{m_s^2 T_{0s}} \int _{0}^{\infty } {\rm d}{\varepsilon } \varepsilon F_{0 s} (\omega - n \omega _{\ast s} ) \int _{1/B_{\text{max}}}^{1/B_{\text{min}}} \frac {{\rm d}{\lambda }}{\Omega _{b s}} \frac {\phi ^{n_2}_{s} (\phi ^{n_2}_{s} - A^{n_2}_{s})}{\omega - n_2 \Omega _{b s} - n \Omega _{\mathrm{d} s}}, \end{align}
where the sum over
$p$
denotes a sum over all bounce wells along the field line. For stellarator geometry, this includes trapped particles that traverse multiple bounce wells before bouncing back. The sum over
$p$
for stellarators, then, is effectively a sum over all possible pairs of local magnetic field strength minima and maxima that lead to a trapped particle orbit. Note that we are also summing over all bounce harmonics via
$n_2$
. As in the passing case, the exponentials (4.39)–(4.40) are highly oscillatory due to the term proportional to
$\omega / |v_{\parallel }|$
. We have moved this term into the Landau resonance. Meanwhile, the remaining eikonals in the Fourier transform pairs are well behaved as
$\varepsilon \to 0$
, since
$\omega _{d} / |v_{\parallel }| \propto \sqrt {\varepsilon }.$
4.3.3. Ampere’s law
The process for calculating
\begin{equation} \mathcal{L}_{A, tr, s} = Z_s e \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \int {\rm d}^3{v} |v_{\parallel }| A_{\parallel } J_{0 s} \left (h_{s, tr, +} - h_{s, tr, -}\right ) \end{equation}
is similar, so we quote the final result,
\begin{equation} \mathcal{L}_{A, tr, s} = \sum _p \sum _{n_2} \frac {4 \pi ^2 Z_s^2 e^2}{m_s^2 T_{0s}} \int _{0}^{\infty } {\rm d}{\varepsilon } \varepsilon F_{0 s} (\omega - n \omega _{\ast s} )\int _{1/B_{\text{max}}}^{1/B_{\text{min}}} \frac {{\rm d}{\lambda }}{\Omega _{b s}} \frac {A^{n_2}_{s} \big(A^{n_2}_{s} - \phi ^{n_2}_{s}\big)}{\omega - n_2 \Omega _{b s} - n \Omega _{\mathrm{d} s}}. \end{equation}
4.4. Full solution and weak form
The quasineutrality equation and parallel Ampere’s law read
\begin{align} \sum _s \frac {n_{0 s}^2 Z_s^2 e^2 }{T_{0s}} \phi - \sum _s Z_s e \int {\rm d}^3{v} J_{0 s} (h_{s, +} + h_{s, -} )&= 0, \end{align}
\begin{align} \frac {k_{\perp }^2}{\mu _0} A_{\parallel } - \sum _s Z_s e \int {\rm d}^3{v} J_{0 s} |v_{\parallel }| (h_{s, +} - h_{s, -} ) & = 0. \end{align}
The proper variational principle is to define the integral
\begin{equation} \begin{split} \mathcal{L} = \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \phi \left (\sum _s \frac {n_{0 s}^2 Z_s^2 e^2 }{T_{0s}} \phi - \sum _s Z_s e \int {\rm d}^3{v} J_{0 s} (h_{s, +} + h_{s, -} ) \right ) \\ + \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } A_{\parallel } \left ( \frac {k_{\perp }^2}{\mu _0} A_{\parallel } - \sum _s Z_s e \int {\rm d}^3{v} J_{0 s} |v_{\parallel }| (h_{s, +} - h_{s, -} )\right ). \end{split} \end{equation}
We have already computed the resonant components. We define the non-resonant pieces as
\begin{equation} \mathcal{L}_{0} = \sum _s \frac {n_{0 s}^2 Z_s^2 e^2 }{T_{0s}} \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \phi ^2 + \int _{-\infty }^{\infty } {\rm d}{\theta } J_{\psi } \frac {k_{\perp }^2}{\mu _0} A_{\parallel }^2. \end{equation}
Then we simply write
\begin{equation} \mathcal{L} = \mathcal{L}_{0} - \sum _s \left (\mathcal{L}_{p, s} + \mathcal{L}_{tr, s}\right ), \end{equation}
where
\begin{equation} \begin{split} \mathcal{L}_{p, s} & = \mathcal{L}_{\phi , p, s} + \mathcal{L}_{A, p, s} \\[5pt] & = \frac {4\pi ^2 Z_s^2 e^2}{m_s^2 T_{0s}} \int _{0}^{\infty } {\rm d}{\varepsilon } \varepsilon F_{0 s} \int _0^{1 / B_{\text{max}}} \frac {{\rm d}{\lambda }}{\Omega _{t s}} \int _{-\infty }^{\infty } {\rm d}{k} \frac {\left (\omega - n \omega _{\ast s} \right ) \left (\phi ^{+}_{s} - A^{+}_{s} \right ) \left (\phi ^{-}_{s} - A^{-}_{s} \right ) }{\omega - k \Omega _{t s} - n \Omega _{\mathrm{d} s}} \end{split} \end{equation}
and
\begin{equation} \begin{split} \mathcal{L}_{tr, s} & = \mathcal{L}_{\phi , tr, s} + \mathcal{L}_{A, tr, s} \\[5pt] & = \sum _{p} \sum _{n_2} \frac {4 \pi ^2 Z_s^2 e^2}{m_s^2 T_{0s}} \int _{0}^{\infty } {\rm d}{\varepsilon } \varepsilon F_{0 s} \int _{1/B_{\text{max}}}^{1/B_{\text{min}}} \frac {{\rm d}{\lambda }}{\Omega _{b s}} \frac {\left (\omega - n \omega _{\ast s} \right )\left (\phi ^{n_2}_{s} - A^{n_2}_{s}\right )^2 }{\omega - n_2 \Omega _{b s} - n \Omega _{\mathrm{d} s}}. \end{split} \end{equation}
In effect, we are just computing transit and bounce averages of the potentials along the field line, taking into account finite Larmor radius effects and drift deviations. The form for passing particles and trapped particles is essentially the same, except that for trapped particles we discretely sum over all bounce harmonics and bounce wells. In contrast, for passing particles we must compute a continuum integral in Fourier space.
If we parameterize
$(\phi , A_{\parallel })$
with tuples of complex numbers
$\alpha$
, then the variational principle is
\begin{align} \mathcal{L}(\omega , \alpha ) &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {\partial \mathcal{L}}{\partial a} (\omega , \alpha ) & = 0. \end{align}
If we assume local analyticity of
$\omega$
in terms of
$\alpha$
, this essentially implies the following stationary principle:
\begin{equation} \frac {\partial \omega }{\partial \alpha } = 0. \end{equation}
One can see this and naively assume a min–max principle applies as in quantum mechanics. For example, one could think that the growth rate
$\gamma$
might be maximized at the solution. However, we note that by the maximum modulus principle, the appropriate solution will generally correspond to
$\gamma (\alpha )$
being a saddle point. Essentially, because the eigenvalues are complex and we have a non-Hermitian system, we cannot guarantee a min–max principle for the growth rate.
Finally, we point out that we can extend the analysis to include
$\delta B_{\parallel }$
by considering the parallel part of Ampere’s law and including the appropriate terms as seen in § 2.
5. Conclusions
In this work, we have corrected the variational principle for local linear collisionless gyrokinetics. In doing so, we have also derived a rigorous weak form for the field equations in general geometry and with electromagnetic effects with no other approximations. While the methods in Rewoldt et al. (Reference Rewoldt, Tang and Chance1982) and Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990) were originally designed for tokamak geometry, we have extended them to account for (periodized) stellarator geometry. Moreover, we have modified the methods and integral transforms in the preceding works to produce more convenient expressions for passing particles and also allow for the variation of nonlinear coefficients for any set of suitable trial functions. For example, rather than setting a constant Gaussian width as in Rewoldt et al. (Reference Rewoldt, Tang and Chance1982), one can now directly vary the width itself per the variational principle (as one often does in introductory quantum mechanics).
The equations as written can be further reduced and modified with various approximations, as done in QuaLiKiz (Stephens et al. Reference Stephens, Garbet, Citrin, Bourdelle, van de Plassche and Jenko2021). Moreover, there exist two obvious theoretical extensions. The first would be to include a non-trivial collision operator, such as pitch angle scattering. Such a method would require that the distribution function be expanded in terms of basis functions, such as eigenfunctions of the collision operator. One would then obtain coupled linear differential equations for components of the distribution function. Then, one would either need to prove complex symmetry of the resulting field equations or else derive another variational principle. The second extension would be to relax the local assumption and examine global eigenmodes with a finite radial width. It would be interesting to determine if a similar variational principle exists for global eigenmodes.
The derived equations can also be used to develop codes to solve the eigenvalue problem with fewer assumptions. For example, QuaLiKiz assumes a shifted circular geometry and does not include electromagnetic effects. A variational code that included both electromagnetic effects and general geometry would prove advantageous in simulating microinstabilities in stellarator geometries. Moreover, the variational approach allows for a self-consistent determination of eigenvalues and eigenfunctions without requiring a completely predetermined trial function.
Lastly, one could also develop codes that utilize state-of-the-art finite element techniques to obtain fast but highly accurate eigenvalues and eigenmodes. A conventional approach to solving the eigenvalue equation (such as in the gyrokinetic code GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000)) is to employ a fixed finite difference discretization in phase space (including velocity space). In contrast, the above variational approach discretizes velocity space via numerical integration, which generically can be adaptive. Moreover, when discretizing the field line following coordinate via a polynomial finite element method, one can adaptively refine the mesh with a combination of h-refinement (subdividing the elements) and p-refinement (using higher-order polynomials), called hp refinement (Szabó & Babuška Reference Szabó and Babuška2021). Moreover, since we only need to consider finite elements along a single dimension (the field line following coordinate), the number of refinement candidates would be relatively modest.
Acknowledgments
The authors thank J. Burby, N. Cao, W. Barham and M. Ruth for insightful discussions. We also thank D. Hatch and S. Mahajan for providing useful references.
Editor Peter Catto thanks the referees for their advice in evaluating this article.
Funding
This work was supported by the U.S. Department of Energy under grant no. DE-FG02-04ER54742 (IFS).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Guiding centre orbit analysis
A.1. Preliminaries
To write down an explicit solution to the gyrokinetic equation in variational form, it is most convenient to take advantage of various integral transforms. These integral transforms require solving the guiding centre equations to zeroth and first order. These ordinary differential equations can be solved ahead of time for a given geometry with a numerical integrator, and their solutions define the integral transforms used in § 4. Using the Clebsch coordinates defined in § 4, the guiding centre equations of motion are then
\begin{align} \frac {\mathrm{d}{\psi }}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi ,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\alpha }}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha ,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\theta }}{{\rm d}{t}} & = v_{\parallel } \boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta + \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta ,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {\mathrm{d}{v_{\parallel }}}{{\rm d}{t}} & = - \frac {\mu }{m} \boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B, \end{align}
where
$t$
represents time. The magnetic drift velocity
$\boldsymbol{v}_{d}$
is given by
\begin{equation} \boldsymbol{v}_d = \frac {m}{Z e B} \frac {v_{\parallel }^2 + \mu B / m}{B^2} \boldsymbol{B} \times \boldsymbol{\nabla} B + \frac {m}{Z e B} \frac {\mu _0 v_{\parallel }^2}{B^3} \boldsymbol{B} \times \boldsymbol{\nabla} p. \end{equation}
Both the energy
$\varepsilon$
and the magnetic moment
$\mu$
are constants of motion. We note that
$\psi$
and
$\alpha$
only change due to the drift motion, whereas
$\theta$
changes both due to the parallel motion and the drift motion. The parallel velocity in turn changes due to the mirror force. Also, note that
\begin{equation} \boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} f = \frac {1}{J_{\psi } B} \frac {\partial f}{\partial \theta }. \end{equation}
It will be convenient to define the pitch angle parameter,
\begin{equation} \lambda = \frac {\mu }{\varepsilon }, \end{equation}
which is also a constant of motion. This allows us to write the parallel velocity as
\begin{equation} v_{\parallel }^2 = \frac {2 \varepsilon }{m} \left (1 - \lambda B \right ). \end{equation}
Rather than solve these equations of motion self-consistently, it is advantageous to separate the system into fast time scales and slow time scales. This is much like how the fast cyclotron motion was averaged out to obtain the guiding centre equations of motion. Therefore, we solve the equations perturbatively: the fast field-line motion is determined by
\begin{align} \frac {\mathrm{d}{\psi }}{{\rm d}{t}} & = 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\alpha }}{{\rm d}{t}} & = 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\theta }}{{\rm d}{t}} & = v_{\parallel } \boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta ,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{v_{\parallel }}}{{\rm d}{t}} & = - \frac {\mu }{m} \boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B. \end{align}
These are the zeroth-order guiding centre equations of motion. The first two equations can be trivially integrated out as
$\psi = \psi _0$
and
$\alpha = \alpha _0$
, which determine the reference flux surface and the reference field line. To lowest order, the particle is bound to the field line according to the mirror force. All magnetic field quantities are evaluated at
$\psi _0, \alpha _0$
, meaning that they are purely functions of
$\theta .$
This reduced system is integrable and Hamiltonian, where the Hamiltonian
$H$
is
\begin{equation} H = \frac {p_{\theta }^2}{2 m J_{\psi }^2 B^2} + \mu B. \end{equation}
Notice that the kinetic term depends on
$\theta$
. This Hamiltonian can also be put into a separable form (meaning
$H = T(p) + V(q)$
), which is convenient for certain ordinary differential equation integrators. Define the arclength
$l$
along the field line as
\begin{equation} l(\theta ) = \int _0^{\theta } {\rm d}{\theta '} J_{\psi } B, \end{equation}
meaning
$\boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} l = 1$
and
$v_{\parallel } = \dot {l}$
. The Hamiltonian can then be rewritten as
\begin{equation} H' = \frac {p^2}{2 m} + \mu B, \end{equation}
where
$B$
is now a function of
$l$
. The Hamiltonian is now separable.
When analysing the field-line motion, there are two types of particles: trapped particles and passing (circulating) particles. A trapped particle is trapped in a magnetic well and thus periodically bounces between two mirror points, akin to a pendulum; this is termed a bounce orbit. At these mirror points,
$v_{\parallel } = 0$
and changes sign. A passing particle, instead, has just the right pitch angle that
$|v_{\parallel }| \gt 0$
always, meaning it never changes direction; this is termed a transit orbit. It is similar to giving a pendulum enough kinetic energy such that it goes round and round over the top. Therefore, it is convenient to analyse trapped particles and passing particles separately.
To first order, we calculate the magnetic drift velocity and evaluate it at
$(\psi _0, \theta (t), \alpha _0)$
. We then solve
\begin{align} \frac {{\rm d}{\tilde {\psi }}} {{\rm d}{t}} &= \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \Big |_{\psi _0, \theta (t), \alpha _0},\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\tilde {\alpha }}}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha \Big |_{\psi _0, \theta (t), \alpha _0}. \end{align}
Here,
$\tilde {\psi }$
and
$\tilde {\alpha }$
are the deviations from the flux surface and field line, respectively. These are first-order equations of motion where the right-hand side is purely a function of
$t$
, so they can be directly integrated. For omnigenous fields, the radial deviations will tend to average to zero over short time scales (Helander Reference Helander2014), whereas the field line deviations will not. This means that the particle will slowly precess around the flux surface. Over very long time scales, the particle may end up on a different field line on the flux surface, in turn affecting the bounce-transit dynamics. So, we expect this analysis to only hold on time scales shorter than the drift time. Moreover, it only holds if
$\rho ^\ast \ll 1$
, where
$\rho ^\ast$
is the gyroradius divided by some macroscopic length scale, such as the machine size; this is because
$v_{\parallel }/v_{d} \sim \rho ^\ast$
.
For certain geometries, the solutions to
$\tilde {\psi }$
and
$\tilde {\alpha }$
will consist of a secular piece and a periodic piece. We will denote the periodic piece of each as
$\hat {\psi }$
and
$\hat {\alpha }$
, respectively, when appropriate.
A.2. Passing particle orbits
A.2.1. Zeroth-order transit motion
Passing particles have to be treated slightly differently in stellarators versus tokamaks; the periodicity of the magnetic field in a tokamak makes the motion along the file line essentially periodic (modulo poloidal turns), which is not the case for a general stellarator geometry. Therefore, certain details apply only to tokamak geometries.
A passing particle has non-vanishing parallel velocity, meaning that along the field
$v_{\parallel }$
never changes sign,
\begin{equation} v_{\parallel } = \pm \sqrt {\frac {2 \varepsilon }{m}} \sqrt {1 - \lambda B}. \end{equation}
Let
$B_{\text{max}}$
denote the maximum magnetic field strength along the field line. The particle is then passing if
\begin{equation} 0 \leqslant \lambda \leqslant \frac {1}{B_{\text{max}}}. \end{equation}
Thus,
$\lambda = 0$
corresponds to a deeply passing particle that does not feel the mirror force whatsoever, and
$\lambda = 1/B_{\text{max}}$
corresponds to a barely passing particle that rests in unstable equilibrium at the top of the global magnetic field well.
Because the magnetic field is periodic in a tokamak, the transit motion is essentially periodic. The transit period is defined as
\begin{equation} \tau _t = \oint {\rm d}{\theta } \left (\frac {{\rm d}{\theta }}{{\rm d}{t}}\right )^{-1} = \sqrt {\frac {m}{2 \varepsilon }} \int _{-\pi }^{\pi } \frac {J_{\psi } B {\rm d}{\theta }}{\sqrt {1 - \lambda B}}. \end{equation}
Since all quantities in the integrand are
$2 \pi$
periodic, the limits of the integral do not matter as long as the upper limit is
$2 \pi$
greater than the lower limit. This is analogous to finding the period of a pendulum where the bob can swing through the top. It is convenient to define a dimensionless transit time,
\begin{equation} \hat {\tau }_t = \frac {1}{q R_0} \int _{-\pi }^{\pi } \frac {J_{\psi } B {\rm d}{\theta }}{\sqrt {1 - \lambda B}}. \end{equation}
The transit frequencies (both with and without dimensions) are defined as
\begin{align} \Omega _t &= \frac {2 \pi }{\tau _t},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {\Omega }_t &= \frac {2 \pi }{\hat {\tau }_t}. \end{align}
Solving for the transit motion as a function of time introduces redundancies. If one varies the energy without varying
$\lambda$
, then one obtains the same motion, just on a different time scale. (This is similar to adjusting the length of a pendulum and adjusting the kinetic energy to compensate.) We therefore introduce the action angle variable
$\alpha _t$
, defined as
\begin{equation} \frac {{\rm d}{\alpha }_t}{{\rm d}{t}} = \Omega _t = \frac {1}{q R_0} \sqrt {\frac {2 \varepsilon }{m}} \hat {\Omega }_t. \end{equation}
We can then rewrite the second-order equation of motion as
\begin{equation} \frac {{\rm d}}{{\rm d}{\alpha _t}} \left (\frac {J_\psi B}{q R_0} \frac {{\rm d}{\theta }}{{\rm d}{\alpha _t}}\right ) = - \frac {q R_0}{J_\psi B} \frac {\lambda }{2 \hat {\Omega }_t^2} \frac {\partial B}{\partial \theta }. \end{equation}
The modified Hamiltonian is then
\begin{equation} \hat {H}(\theta , \hat {p}_{\theta }) = \frac {1}{2} \left (\frac {\hat {p}_{\theta } q R_0}{J_{\psi } B}\right )^2 + \frac {\lambda B}{2 \hat {\Omega }_2^2}. \end{equation}
Then,
$\theta (\alpha _t)$
will be
\begin{equation} \theta (\alpha _t) = \alpha _t + \tilde {\theta }(\alpha _t), \end{equation}
where
$\tilde {\theta }$
is
$2 \pi$
periodic. So the solution is the sum of a linear secular term and a periodic term. At this point, we can solve for
$\alpha _t$
via any ordinary differential equation integrator. Note that if we wish to make the Hamiltonian separable, that is of the form
\begin{equation} H(l, p_l) = \frac {p_l^2}{2} + \frac {f(l)}{2}, \end{equation}
we can transform to the normalized arclength
$l$
via
\begin{equation} l(\theta ) =\frac {1}{q R_0} \int _0^{\theta } {\rm d}{\theta }' J_{\psi } B, \end{equation}
in which case
\begin{equation} \hat {H} = \frac {p_l^2}{2} + \frac {\lambda B}{2 \hat {\Omega }_2^2}. \end{equation}
The initial conditions turn out to be arbitrary. For simplicity, we set
\begin{align} \theta (\alpha _t = 0) &= 0,\end{align}
\begin{align} \frac {{\rm d}{\theta }}{{\rm d}{\alpha _t}}(\alpha _t = 0) &= \left (\frac {q R_0}{J_{\psi } B} \frac {\sqrt {1 - \lambda B}}{\hat {\Omega }_t} \right ) \Bigg |_{\theta = 0}.\end{align}
The second initial condition can be obtained by considering we want
\begin{equation} \frac {{\rm d}{\theta }}{{{\rm d}{t}}}(t=0) = \left ( \sqrt {\frac {2 \varepsilon }{m}} \frac {\sqrt {1 - \lambda B}}{J_{\psi } B} \right ) \Bigg |_{\theta = 0}. \end{equation}
Essentially, we are setting
$v_{\parallel }$
at
$\theta = 0$
to be consistent with our choice in
$\varepsilon$
and
$\lambda$
. Due to time reversibility and periodicity in tokamak geometry, we do not need to separately consider passing particles travelling backwards.
Passing particle orbits in stellarators are complicated by the fact that the magnetic field will in general not be periodic. For a stellarator, we in principle need to calculate the orbit along the entire field line in the domain under consideration. Stellarators also allow for one-sided bounce orbits, which are analogous to hyperbolic orbits in celestial mechanics; we do not consider those here.
For stellarators, it still can pay to normalize the energy out of the problem and consider characteristic length and time scales. We define the characteristic transit time as
\begin{equation} \tau _t = \frac {1}{2 N + 1} \sqrt {\frac {m}{2 \varepsilon }} \int _{-(2 N + 1) \pi }^{(2 N + 1) \pi } \frac {J_{\psi } B {\rm d}{\theta }}{ \sqrt {1 - \lambda B}}, \end{equation}
where
$N$
is a non-negative integer. This is the amount of time it takes to traverse a full poloidal turn averaged over
$2 N + 1$
turns. We can also define a characteristic length as
\begin{equation} L = \frac {1}{2 N + 1} \int _{-(2 N + 1) \pi }^{(2 N + 1) \pi } J_{\psi } B {\rm d}{\theta }. \end{equation}
The normalized transit time is then
\begin{equation} \hat {\tau }_t = \frac {1}{(2 N + 1) L} \int _{-(2 N + 1) \pi }^{(2 N + 1) \pi } \frac {J_{\psi } B {\rm d}{\theta }}{\sqrt {1 - \lambda B}}. \end{equation}
The transit frequency can then be defined similarly as in the tokamak case. We then rescale the time coordinate such that
\begin{equation} \frac {{\rm d}{\alpha _t}}{{\rm d}{t}} = \Omega _t = \frac {1}{L} \sqrt {\frac {2 \varepsilon }{m}} \hat {\Omega }_t. \end{equation}
Then, one solves the equations of motion as in the tokamak case with one’s preferred integrator. Note that, unlike in the tokamak case, one also needs to consider both signs of the parallel velocity separately. One can then obtain
$\theta (\alpha _t)$
.
A.2.2. First-order precession motion
We now examine the first-order precession motion. The first-order guiding centre equations of motion satisfy
\begin{align} \frac {{\rm d}{\tilde {\psi }}}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi ,\\[-10pt]\nonumber\end{align}
\begin{align} \frac {{\rm d}{\tilde {\alpha }}}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha . \end{align}
Although these equations can be simply integrated over time, some physical intuition is in order. As above, periodicity arguments can be exploited in tokamak geometry.
When analysing the drift motion, it is useful to consider how it varies over the transit time. In tokamaks, we thus define the transit average: for any function of
$\theta$
, we define the transit average as
\begin{equation} \left \langle f(\theta )\right \rangle _t = \frac {1}{\tau _t} \int _{0}^{2 \pi } \frac {{\rm d}{\theta } J_{\psi } B}{|v_{\parallel }|} f(\theta ). \end{equation}
Equivalently, if we instead know
$f$
as a function of
$\alpha _t$
, then
\begin{equation} \left \langle f(\alpha _t) \right \rangle _t = \int _0^{2 \pi } \frac {{\rm d}{\alpha _t}}{2 \pi } f(\alpha _t). \end{equation}
First, we note that in tokamaks
\begin{equation} \Omega _{d \psi } = \left \langle \boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \right \rangle _t = 0, \end{equation}
meaning that to first order, the particle does not drift away from the flux surface. As for trapped particles, this would be exactly zero if we integrated the guiding centre equations exactly in a tokamak. This means that
$\tilde {\psi }$
is a periodic function of
$\alpha _t$
. In stellarators, this is only true if the geometry is omnigenous, but we can still split
$\tilde {\psi }$
into a secular and periodic component regardless. For simplicity, we just set
$\tilde {\psi }(\alpha _t = 0) = 0$
as an initial condition. The solution can be written as
\begin{equation} \tilde {\psi }(\alpha _t) = \frac {\Omega _{d \psi }}{\Omega _{t}} \alpha _{t} + \hat {\psi }(\alpha _t), \end{equation}
where
$\hat {\psi }$
is a periodic function of
$\alpha _t$
and can be computed as
\begin{equation} \hat {\psi }(\alpha _t) = \int _0^{\alpha _t} \frac {{\rm d}{\alpha '_t}}{\Omega _t} \left (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - \Omega _{\mathrm{d} \psi } \right ). \end{equation}
Meanwhile, the transit average of
$\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha$
will generally not be zero. Moreover, this quantity is actually not particularly useful to use in tokamaks due to the secular term present in the definition of
$\alpha$
. We define a different drift frequency as
\begin{equation} \Omega _{\mathrm{d} \alpha } = \big\langle \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha + q' \Omega _t \hat {\psi } + q' \alpha _t \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \big\rangle _t. \end{equation}
Then, one can show
\begin{equation} \tilde {\alpha } = \left (\frac {\Omega _{\mathrm{d} \alpha }}{\Omega _t} - q' \hat {\psi }(\alpha _t) - q' \frac {\Omega _{d \psi }}{2 \Omega _{t}} \alpha _{t} \right ) \alpha _t + \hat {\alpha }(\alpha _t), \end{equation}
where
$\hat {\alpha }$
is
$2 \pi$
periodic. Using the initial condition that
$\tilde {\alpha }(\alpha _t = 0) = \hat {\alpha }(\alpha _t = 0)$
, we can write the solution as
\begin{equation} \hat {\alpha }(\alpha _t) = \int _0^{\alpha _t} \frac {{\rm d}{\alpha '_{t}}}{\Omega _t} \left (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha + q' \Omega _t \hat {\psi } + q' \alpha _t \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - \Omega _{\mathrm{d} \alpha } \right ). \end{equation}
To prove that this is indeed the solution, we need to explicitly write
$\alpha = \varphi - q (\theta + \lambda )$
, where
$\lambda$
is a periodic function that makes
$\chi = \theta + \lambda$
straight. Then,
\begin{equation} \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi - q \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} (\theta + \lambda ) - q' (\tilde {\theta } + \lambda ) \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - q' \alpha _t \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi , \end{equation}
where the last term on the right-hand side is the secular piece and everything else is periodic. We note that
\begin{equation} \frac {\alpha _t}{\Omega _t} \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi = \frac {\alpha _{t}}{\Omega _{t}} \Omega _{d \psi } + \frac {{\rm d}}{{\rm d}{\alpha _t}} (\alpha _t \hat {\psi } ) - \hat {\psi }. \end{equation}
Therefore, we can integrate by parts and find
\begin{align} \int _0^{\alpha _t} \frac {{\rm d}{\alpha _t'}}{\Omega _t} \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha &= -q' \alpha _{t}^2 \frac {\Omega _{d \psi }}{\Omega _{t}} - q' \alpha _t \hat {\psi }(\alpha _t)\nonumber\\[5pt] &\quad{}+ \int _0^{\alpha _t} \frac {{\rm d}{\alpha _t'}}{\Omega _t} (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha + q' \Omega _t \hat {\psi } + q' \alpha _t \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi ). \end{align}
The terms under the integral on the right-hand side are all periodic, so we can split it up into a monotonically increasing piece and a periodic piece, which is precisely what we did above,
\begin{equation} \int _0^{\alpha _t} \frac {{\rm d}{\alpha _t'}}{\Omega _t} \big(\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha + q' \Omega _t \hat {\psi } + q' \alpha _t \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \big) = \frac {\Omega _{\mathrm{d} \alpha }}{\Omega _t} \alpha _t + \hat {\alpha }(\alpha _t). \end{equation}
This completes the proof.
In stellarators, there is strictly speaking no periodic motion to exploit for passing particle trajectories. However, the entire field line is not typically simulated in flux-tube simulations. Instead, one includes only a finite number of poloidal turns (Sánchez et al. Reference Sánchez2021). In the ballooning representation, one can take advantage of this by periodizing the magnetic field after some number of poloidal turns. Since the magnetic field is then periodic along the field line, the analysis presented in the tokamak-focused sections can be extended to stellarator geometry. In general, it is also necessary to consider a transit average of the radial drift frequency
$\Omega _{\mathrm{d} \psi }$
in stellarator geometries.
A.2.3. Useful integral identities
Now that we better understand transit dynamics, we briefly derive some useful expressions. For stellarator geometries, we assume that the geometry has been periodized. We encounter two types of integrals in the main text when integrating over the field line. The first is
\begin{equation} \int _{0}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} f(\theta ') = \tau _t \int _0^{\alpha _t} \frac {{\rm d}{\alpha _t'}}{2 \pi } f\left (\theta (\alpha _t')\right ). \end{equation}
The second type directly involves the magnetic drift resonance,
\begin{equation} \int _{0}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{d} \right ), \end{equation}
where
$\omega$
is the mode frequency,
$n$
is the toroidal mode number and
$\omega _{d}$
is defined as
\begin{equation} \omega _{d} = \boldsymbol{v}_{d} \boldsymbol{\cdot} \left (\boldsymbol{\nabla} \alpha + \theta _0 q' \boldsymbol{\nabla} \psi \right ), \end{equation}
where
$\theta _0$
is the ballooning angle. This can be rewritten as
\begin{equation} \int _{0}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{d} \right ) = \frac {\omega }{\Omega _t} \alpha _t - n \tilde {\alpha }(\alpha _t) - n \theta _0 q' \tilde {\psi }(\alpha _t). \end{equation}
This simplifies to
\begin{equation} \int _{0}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{d} \right ) = \frac {\alpha _t}{\Omega _t} \left ( \omega - n \Omega _{d} + \frac {1}{2} n q' \Omega _{d \psi } \alpha _{t} + n q' \Omega _{t} \hat {\psi }(\alpha _t) \right ) - n \tilde {w}_{d}, \end{equation}
where
\begin{align} \Omega _{d} &= \Omega _{d \alpha } - \theta _0 q' \Omega _{d \psi },\\[-10pt]\nonumber \end{align}
\begin{align} \tilde {w}_{d} & = \hat {\alpha }(\alpha _t) + \theta _0 q' \hat {\psi }(\alpha _t). \end{align}
Using action angle variables, we have successfully split these integrals into a secular term and a periodic term. In the main text, it is useful to write
\begin{equation} \int _{0}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{d} \right ) = \frac {\alpha _t}{\Omega _t} \left ( \omega - n \Omega _{d} \right ) - n \Lambda _d, \end{equation}
where
$\Lambda _d$
contains all drift deviation terms,
\begin{equation} \Lambda _d = \tilde {w}_{d} - q' \frac {\Omega _{d \psi }}{2 \Omega _{t}} \alpha _{t}^2 - q' \hat {\psi } \alpha _{t}. \end{equation}
A.3. Trapped particle orbits
A.3.1. Zeroth-order bounce motion
A trapped particle has parallel velocity equal to zero at two distinct bounce points and to lowest order undergoes periodic motion between those bounce points. The parallel velocity is
\begin{equation} v_{\parallel } = \pm \sqrt {\frac {2 \varepsilon }{m}} \sqrt {1 - \lambda B}. \end{equation}
For a given
$\lambda$
, one can determine the bounce points by solving
\begin{equation} 1 - \lambda B = 0, \end{equation}
where
$B$
is purely a function of the field line following coordinate
$\theta$
.
Note that in tokamak geometry, we only need to calculate two bounce points; since the magnetic field is periodic, all other pairs of bounce points are equal modulo the period. In a stellarator, the magnetic field is no longer periodic and there will be multiple distinct magnetic wells. Moreover, a trapped particle can traverse over multiple bounce wells along an orbit.
It is often useful to use the trapped parameter
$\kappa$
instead of
$\lambda$
. One common definition is
\begin{equation} \kappa ^2 = \frac {1 - \lambda B_{\text{min}}}{ 2 \epsilon _B} \end{equation}
where
\begin{equation} 2 \epsilon _B = 1 - \frac {B_{\text{min}}}{B_{\text{max}}}. \end{equation}
Here,
$B_{\text{min}}$
and
$B_{\text{max}}$
are, respectively, the minimum and maximum of the magnetic well under consideration; this could be a large ‘well’ that is composed of smaller local bounce wells, thus accounting for trapped particles that traverse multiple bounce wells along an orbit. Let
$\theta _{\text{min}}$
be the location of the minimum and
$\theta _{\text{max}}$
the location of the smallest nearest maximum. (For a tokamak,
$\theta _{\text{max}}$
always comes in pairs that correspond to the same maximum. For a stellarator, there can be two distinct nearest maxima; we pick the smaller one.) The above definition then guarantees that
\begin{align} v_{\parallel }(\theta = \theta _{\text{min}}) &= 0 \iff \kappa = 0,\\[-10pt]\nonumber \end{align}
\begin{align} v_{\parallel }(\theta = \theta _{\text{max}}) &= 0 \iff \kappa = 1. \end{align}
Therefore, at
$\kappa = 0$
the particle is deeply trapped and confined to the bottom of the well. At
$\kappa = 1$
, the particle is barely trapped in unstable equilibrium at the top of the local well. Moreover,
\begin{align} v_{\parallel }^2 &= \frac {2 \varepsilon }{m} \left ( 1 - \left (1 - 2 \epsilon _B \kappa ^2 \right ) \frac {B}{B_{\text{min}}}\right ),\\[-12pt]\nonumber \end{align}
\begin{align} v_{\perp }^2 & = \frac {2 \varepsilon }{m} \left (1 - 2 \epsilon _B \kappa ^2 \right ) \frac {B}{B_{\text{min}}}. \end{align}
Therefore, all trapped particles satisfy
\begin{equation} 0 \leqslant \kappa \leqslant 1. \end{equation}
The zeroth-order bounce motion is periodic, with a bounce period equal to
\begin{equation} \tau _b = \oint {\rm d}{\theta } \left (\frac {{\rm d}{\theta }}{{\rm d}{t}}\right )^{-1} = 2 \sqrt {\frac {m}{2 \varepsilon }} \int _{\theta _1}^{\theta _2} \frac {J_{\psi } B {\rm d}{\theta }}{\sqrt {1 - \lambda B}}, \end{equation}
where
$\theta _1 \lt \theta _{\text{min}} \lt \theta _2$
are the two bounce points. Since we are considering the full forward and then backwards motion, a factor of 2 is included. For tokamaks especially, it is convenient to define a dimensionless bounce time,
\begin{equation} \hat {\tau }_b = \frac {2}{q R_0} \int _{\theta _1}^{\theta _2} \frac {J_{\psi } B {\rm d}{\theta }}{\sqrt {1 - \lambda B}}. \end{equation}
The bounce frequencies (both with and without dimensions) are defined as
\begin{align} \Omega _b & = \frac {2 \pi }{\tau _b},\\[-12pt]\nonumber \end{align}
\begin{align} \hat {\Omega }_b & = \frac {2 \pi }{\hat {\tau }_b}. \end{align}
Solving for the bounce motion as a function of time introduces redundancies. If one varies the energy without varying
$\lambda$
, then one obtains the same motion, just on a different time scale. (This is similar to adjusting the length of a pendulum but keeping the bounce angles the same.) We therefore introduce the action angle variable
$\alpha _b$
, defined as
\begin{equation} \frac {{\rm d}{\alpha }_b}{{\rm d}{t}} = \Omega _b = \frac {1}{q R_0} \sqrt {\frac {2 \varepsilon }{m}} \hat {\Omega }_b. \end{equation}
We can then rewrite the second-order equation of motion as
\begin{equation} \frac {{\rm d}}{{\rm d}{\alpha _b}} \left (\frac {J_\psi B}{q R_0} \frac {{\rm d}{\theta }}{{\rm d}{\alpha _b}}\right ) = - \frac {q R_0}{J_\psi B} \frac {\lambda }{2 \hat {\Omega }_b^2} \frac {\partial B}{\partial \theta }. \end{equation}
The modified Hamiltonian is then
\begin{equation} \hat {H}(\theta , \hat {p}_{\theta }) = \frac {1}{2} \left (\frac {\hat {p}_{\theta } q R_0}{J_{\psi } B}\right )^2 + \frac {\lambda B}{2 \hat {\Omega }_b^2}. \end{equation}
Then,
$\theta (\alpha _b)$
is
$2 \pi$
periodic. At this point, we can solve for
$\alpha _b$
via any ordinary differential equation integrator. Note that if we wish to make the Hamiltonian separable, that is of the form
\begin{equation} H(l, p_l) = \frac {p_l^2}{2} + \frac {f(l)}{2}, \end{equation}
we can transform to the normalized arclength
$l$
via
\begin{equation} l(\theta ) =\frac {1}{q R_0} \int _0^{\theta } {\rm d}{\theta }' J_{\psi } B, \end{equation}
in which case
\begin{equation} \hat {H} = \frac {p_l^2}{2} + \frac {\lambda B}{2 \hat {\Omega }_b^2}. \end{equation}
In terms of initial conditions, it is useful to start at the first bounce point,
\begin{align} \theta (\alpha _b = 0) &= \theta _1,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\theta }}{{\rm d}{\alpha _b}}(\alpha _b = 0) &= 0. \end{align}
This guarantees
$\theta (\alpha _b = \pi ) = \theta _2$
. Due to time-reversibility,
$\theta$
is also an even function of
$\alpha _2$
. (In fact,
$\theta (n \pi + \alpha _b) = \theta (n \pi - \alpha _b)$
for all integers
$n$
.)
A.3.2. First-order drift motion
The first-order guiding centre equations of motion satisfy
\begin{align} \frac {{\rm d}{\tilde {\psi }}}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi ,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{\rm d}{\tilde {\alpha }}}{{\rm d}{t}} & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha , \end{align}
where the magnetic drift velocity
$\boldsymbol{v}_{d}$
and the quantities
$\boldsymbol{\nabla} \theta , \boldsymbol{\nabla} \psi $
and
$\boldsymbol{\nabla} \alpha$
are all evaluated at
$\theta (t), \psi _0, \alpha _0$
. Although the second set of equations can be simply integrated over time, some physical intuition is in order.
When analysing the drift motion, it is useful to consider how it varies over the bounce time. We thus define the bounce average: for any function of
$\theta$
, we define the bounce average as
\begin{equation} \left \langle f(\theta )\right \rangle _b = \frac {1}{\tau _b} \oint \frac {{\rm d}{\theta } J_{\psi } B}{|v_{\parallel }|} f(\theta ), \end{equation}
where the closed integral signifies we need to consider forwards and backwards motion. (So for example,
$\left \langle v_{\parallel }\right \rangle _b = 0$
). Equivalently, if we instead know the function as a function of
$\alpha _b$
, then
\begin{equation} \left \langle f(\alpha _b) \right \rangle _b = \int _0^{2 \pi } \frac {{\rm d}{\alpha _b}}{2 \pi } f(\alpha _b). \end{equation}
In tokamaks and omnigenous stellarators, we expect
\begin{equation} \Omega _{d \psi } = \left \langle \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \right \rangle _b = 0, \end{equation}
meaning that to first order, the particle does not drift away from the flux surface (Helander Reference Helander2014). Note that if we solved the full guiding centre equations of motion in a tokamak or a quasisymmetric stellarator, this would be exact due to the existence of a third invariant (in a tokamak, this would be the canonical toroidal momentum).
If we set the initial condition to be
$\tilde {\psi }(\alpha _b = 0) = 0$
, then
\begin{equation} \tilde {\psi }(\alpha _b) = \frac {\Omega _{d \psi }}{\Omega _b} \alpha _b + \hat {\psi }(\alpha _b), \end{equation}
where,
$\hat {\psi }$
is
$2 \pi$
periodic and an odd function of
$\alpha _b$
. (In fact,
$\hat {\psi }(n \pi + \alpha _b) = - \hat {\psi }(n \pi - \alpha _b)$
for all integers
$n$
.) The solution can be written as
\begin{equation} \hat {\psi }(\alpha _b) = \int _0^{\alpha _b} \frac {{\rm d}{\alpha _b'}}{\Omega _b} \left (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - \Omega _{\mathrm{d} \psi }\right ). \end{equation}
When calculating
$\tilde {\alpha }$
, we cannot guarantee that the bounce average of
$\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha$
will be zero. This quantity is defined as the bounce-averaged drift frequency and characterizes the drift from field line to field line,
\begin{equation} \Omega _{d \alpha } = \left \langle \boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha \right \rangle _b. \end{equation}
Therefore, the appropriate solution is
\begin{equation} \tilde {\alpha }(\alpha _b) = \frac {\Omega _{d \alpha }}{\Omega _b} \alpha _b + \hat {\alpha }(\alpha _b), \end{equation}
where
\begin{equation} \hat {\alpha }(\alpha _b) = \int _0^{\alpha _b} \frac {{\rm d}{\alpha _b'}}{\Omega _b} \left (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha - \Omega _{\mathrm{d} \alpha } \right ). \end{equation}
As before, we have chosen
$\tilde {\alpha }(\alpha _b = 0)$
as our initial condition, and
$\hat {\alpha }$
is a periodic odd function of
$\alpha _b$
such that
$\hat {\alpha }(n \pi + \alpha _b) = - \hat {\alpha }(n \pi - \alpha _b)$
for all integers
$n$
.
A.3.3. Useful integral identities
Now that we better understand bounce dynamics, we briefly derive some useful expressions. We encounter two types of integrals in the main text when integrating over the field line. The first is
\begin{equation} 2 \int _{\theta _1}^{\theta _2} \frac {{\rm d}{\theta } J_{\psi } B}{|v_{\parallel }|} f(\theta ) = \tau _b \left \langle f(\theta ) \right \rangle _b = \tau _b \int _0^{2 \pi } \frac {{\rm d}{\alpha _b}}{ 2 \pi } f\left (\theta (\alpha _b)\right ). \end{equation}
The second type directly involves the magnetic drift resonance,
\begin{equation} \int _{\theta _1}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{d} \right ), \end{equation}
where
$\omega$
is the mode frequency,
$n$
is the toroidal mode number and
$\omega _{d}$
is defined as
\begin{equation} \omega _{d} = \boldsymbol{v}_{d} \boldsymbol{\cdot} \left (\boldsymbol{\nabla} \alpha + \theta _0 q' \boldsymbol{\nabla} \psi \right ), \end{equation}
where
$\theta _0$
is the ballooning angle. This integral can be rewritten as
\begin{equation} \int _{\theta _1}^{\theta } \frac {{\rm d}{\theta }' J_{\psi } B}{|v_{\parallel }|} \left (\omega - n \omega _{d} \right ) = \frac {\alpha _b}{\Omega _b} \left (\omega - n \Omega _{d} \right ) - n \tilde {w}_{d}, \end{equation}
where
\begin{align} \Omega _{d} &= \Omega _{\mathrm{d} \alpha } - \theta _0 q' \Omega _{\mathrm{d} \psi },\\[-8pt]\nonumber \end{align}
\begin{align} \tilde {w}_d &= \int _0^{\alpha _b} {\rm d}{\alpha _b'} \frac {\omega _{d} - \Omega _{d \alpha } + \theta _0 q' \Omega _{\mathrm{d} \psi }}{\Omega _b} = \hat {\alpha }(\alpha _b) + \theta _0 q' \hat {\psi }(\alpha _b). \end{align}
Therefore, when computing the resonant contribution to the gyrokinetic field equations, we are essentially computing the bounce averaged drift frequency as well as the drift deviations along the guiding centre orbit. Notice that the equations greatly simplify when working in
$\alpha _b$
instead of
$\theta$
; it pays to use action angle variables!
Lastly, we note a useful duplication formula for tokamak geometry. Because all bounce wells are identical and spaced
$2 \pi$
apart along the field line, we can handle all bounce wells simultaneously. Denote the orbits in the well
$2 \pi p$
away from the reference centre well with the subscript
$p$
, where
$p$
is any integer. Quantities without a
$p$
subscript correspond to the orbits of the reference centre well. Then
\begin{align} \Omega _{b, p} & = \Omega _{b},\\[-10pt]\nonumber \end{align}
\begin{align} \Omega _{d \psi , p} & = 0,\\[-10pt]\nonumber \end{align}
\begin{align} \Omega _{d \alpha , p} & = \Omega _{d \alpha },\\[-10pt]\nonumber \end{align}
\begin{align} \theta _p(\alpha _b) &= 2 \pi p + \theta (\alpha _b),\\[-10pt]\nonumber \end{align}
\begin{align} \tilde {\psi }_{p}(\alpha _b) & = \hat {\psi }(\alpha _b),\\[-10pt]\nonumber \end{align}
\begin{align} \tilde {\alpha }_{p} ( \alpha _b) & = \frac {\Omega _{d \alpha }}{\Omega _b} \alpha _b + \hat {\alpha } (\alpha _b) - 2 \pi p q' \hat {\psi }(\alpha _b). \end{align}
This is because
\begin{align} (v_{\parallel })_p & = v_{\parallel },\\[-10pt]\nonumber \end{align}
\begin{align} (\boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B)_p & = \boldsymbol{\hat {b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B,\\[-10pt]\nonumber \end{align}
\begin{align} (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi )_p & = \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi ,\\[-10pt]\nonumber \end{align}
\begin{align} (\boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha )_p &= \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha - 2 \pi p q' \boldsymbol{v}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi . \end{align}
Therefore, we only need to compute the orbits for the
$p = 0$
bounce well in tokamak geometry.
Appendix B. Left eigenvector for the gyrokinetic equation
The electrostatic gyrokinetic equation for each species reads
\begin{equation} \left (\omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} s} \right ) h_{s} - \left (\omega - n \omega _{\ast s}\right ) \frac {Z_{s} e F_{0 s}}{T_{0s}} J_{0 s} \phi = 0. \end{equation}
For quasineutrality, we have
\begin{equation} \sum _{s} \frac {Z_{s}^2 e^2 n_{0 s}}{T_{0 s}} \phi - \int {\rm d}^3{v} Z_{s} e J_{0 s} h_{s} = 0. \end{equation}
For simplicity, we work with two species, ions and electrons. We also set
$Z_i = 1$
,
$n_0 = n_{0i} = n_{0e}$
,
$T_{0} = T_{0i} = T_{0e}$
, and normalize all quantities. Let
$\boldsymbol{\chi }$
be the normalized column vector
$ \left (h_i, h_e, \phi \right )^T$
. We define our inner product as
\begin{equation} \langle \boldsymbol{\chi }_2, \boldsymbol{\chi }_1 \rangle = \int _{-\infty }^{\infty } \frac {{\rm d}{l}}{B} \left (\int {\rm d}^3{v} \frac { h_{i2}^\ast h_{i1}}{n_0 F_{0i}} + \frac {h_{e2}^\ast {h_{e1}}}{n_0 F_{0e}} \right ) + \int \frac {{\rm d}{l}}{B} \frac {e^2}{T_{0}^2} \phi _2^\ast \phi _1 .\end{equation}
The division by the Maxwellian for the distribution functions ensures that the product converges. The above equations can then be written in block form,
\begin{equation} \hat {L} \boldsymbol{\chi } = \begin{pmatrix} \hat {L}_{ii} & \quad \hat {L}_{ie} & \quad \hat {L}_{i\phi } \\[4pt] \hat {L}_{ei} & \quad \hat {L}_{ee} & \quad \hat {L}_{e\phi } \\[4pt] \hat {L}_{\phi i} & \quad \hat {L}_{\phi e} & \quad \hat {L}_{\phi \phi } \end{pmatrix} \begin{pmatrix} h_i \\[4pt] h_e \\[4pt] \phi \end{pmatrix} = \boldsymbol{0}. \end{equation}
Here
\begin{align} \hat {L}_{ii} &= \omega - v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} i} ,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{ie} &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{i \phi } &= - (\omega - n \omega _{\ast i}) F_{0i} J_{0i},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{ei} &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{ee} &= \omega - v_{\parallel } k_{\parallel } - n \omega _{d e},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{e\phi } &= (\omega - n \omega _{\ast e}) F_{0 e} J_{0 e},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi i} &= - \int {\rm d}^3{v} J_{0 i} ,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi e} &= \int {\rm d}^3{v} J_{0 e} ,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi \phi } &= 2. \end{align}
The transpose is
\begin{align} \hat {L}_{ii}^T & = \omega + v_{\parallel } k_{\parallel } - n \omega _{\mathrm{d} i} ,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{ie}^T &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{i \phi }^T &= -\int {\rm d}^3{v} (\omega - n \omega _{\ast i}) J_{0 i},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{ei}^T &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{ee}^T &= \omega + v_{\parallel } k_{\parallel } - n \omega _{d e},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{e\phi }^T &= \int {\rm d}^3{v} (\omega - n \omega _{\ast e}) J_{0 e},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi i}^T &= - J_{0 i} F_{0 i},\\[-10pt]\nonumber \end{align}
\begin{align} \hat {L}_{\phi e}^T &= J_{0 e} F_{0 e}, \\[-10pt]\nonumber\end{align}
\begin{align} \hat {L}_{\phi \phi }^T &= 2. \\[-10pt]\nonumber\end{align}
Note the sign change in
$v_{\parallel }$
. Using adjoint fields
$(g_i, g_e, \psi )$
, the transposed equation is then
\begin{align} \left (\omega + v_{\parallel } k_{\parallel } - n \omega _{d i} \right ) g_i - J_{0 i} F_{0 i} \psi &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \left (\omega + v_{\parallel } k_{\parallel } - n \omega _{d e} \right ) g_e - J_{0 e} F_{0 e} \psi &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} 2 \psi - \int {\rm d}^3{v} \left (\omega - n \omega _{\ast i} \right ) J_{0 i} g_i + \int {\rm d}^3{v} \left (\omega - n \omega _{\ast e} \right ) J_{0 e} g_e & = 0. \end{align}
Notice that the equations are markedly different, so the full system is not symmetric. We can, however, make a connection between the old system and the new system. Assume that the fields
$(h_i, h_e, \phi )$
solve the old system, and let
\begin{align} g_i(l, v_{\parallel }, \mu ) &= \frac {1}{\omega - n \omega _{\ast i}} h_i(l, -v_{\parallel }, \mu ),\\[-10pt]\nonumber \end{align}
\begin{align} g_e(l, v_{\parallel }, \mu ) &= \frac {1}{\omega - n \omega _{\ast e}} h_e(l, -v_{\parallel }, \mu ),\\[-10pt]\nonumber \end{align}
\begin{align} \psi (l) & = \phi (l). \end{align}
Plugging these in (and taking
$v_{\parallel } \to - v_{\parallel }$
), we obtain
\begin{align} \left (\omega + v_{\parallel } k_{\parallel } - n \omega _{d i} \right ) g_i - J_{0i} F_{0i} \psi & \to \left (\omega - v_{\parallel } k_{\parallel } - n \omega _{d i}\right ) \frac {h_i}{\omega - n \omega _{\ast i}} - J_{0 i} F_{0 i} \phi ,\\[-10pt]\nonumber \end{align}
\begin{align} \left (\omega + v_{\parallel } k_{\parallel } - n \omega _{d e} \right ) g_e - J_{0 e} F_{0 e} \psi & \to \left (\omega - v_{\parallel } k_{\parallel } - n \omega _{d e}\right ) \frac {h_e}{\omega - n \omega _{\ast e}} - J_{0 e} F_{0 e} \phi \end{align}
and
\begin{align} 2 \psi &- \int {\rm d}^3{v} \left (\omega - n \omega _{\ast i} \right ) J_{0i} g_i + \int {\rm d}^3{v} \left (\omega - n \omega _{\ast e} \right ) J_{0e} g_e \to 2 \phi \nonumber\\[5pt] &- \int {\rm d}^3{v} J_{0 i} h_i + \int {\rm d}^3{v} J_{0 e} h_e. \end{align}
Note that
$k_{\parallel }$
commutes with
$\omega - n \omega _{\ast }$
. Given that
$(g_i, g_e, \phi )$
solve the old system, each of the above equations is zero. (The quasineutrality equation is obvious, for the other two just multiply by
$\omega - n \omega _{\ast s}$
.) Therefore, the above guess solves the transposed equations and thus corresponds to the left eigenvector of the system.
Open access
