2.1 Magnetic coordinates

We can express the magnetic field using the Clebsch representation

(2.20)\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \psi , \end{equation}

where $\psi$ and $\alpha$ are the flux surface and field line labels, respectively. We choose to work in field-aligned coordinates (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995) denoted by $(x,y,z)$, with $z$ the position along the magnetic field line and $(x,y)$ the position in the plane perpendicular to $\hat {\boldsymbol {b}}$. Together, the unit vectors $\{ \hat {\boldsymbol {x}}, \hat {\boldsymbol {y}}, \hat {\boldsymbol {b}}\}$ form a left-handed, orthonormal basis. We relate the coordinates $(x,y)$ to $(\psi, \alpha )$ using

(2.21)\begin{equation} \left. \begin{aligned} x & = \frac{{\rm d} x}{{\rm d} \psi} (\psi - \psi_0) , \\ y & = \frac{{\rm d} y}{{\rm d} \alpha} (\alpha - \alpha_0) , \end{aligned} \right\} \end{equation}

with $(\psi _0, \alpha _0)$ denoting the values at the centre of the domain.

We take the discrete Fourier transform in $\{x,y\}$ whilst retaining the real-space coordinate for the parallel direction, and evaluate nonlinear terms pseudo-spectrally in order to retain spectral accuracy in spatial derivatives when implementing the equations numerically.Footnote ⁴ The discrete Fourier transform in $\{x,y\}$

(2.22)\begin{equation} g_{\nu} (x,y,z, v_{{\parallel}}, \mu_{\nu},t) = \sum_{k_x, k_y} \hat{g}_{\boldsymbol{k},\nu} (z,v_{{\parallel}}, \mu_\nu,t) e^{{({\rm i}k_x x+ {\rm i} k_y y)}}, \end{equation}

is justified provided the condition $k_x \sim k_y \gg 1/L$ is satisfied. This translates to a requirement that the turbulent fluctuations at the edges of our domain are decorrelated, such that they may be considered statistically identical, and periodic boundary conditions can be enforced in $(x,y)$ (Beer et al. Reference Beer, Cowley and Hammett1995).

2.2 System equations

We next transform our evolution equations (2.14), and (2.17)–(2.19). The definition of $\chi$ is used to find the Fourier transform of the quantity $\langle \chi \rangle _{\boldsymbol {R}_{\nu }} = \langle \phi \rangle _{\boldsymbol {R}_{\nu }} - v_{\parallel } \langle A_{\parallel } \rangle _{\boldsymbol {R}_{\nu }}/c - \langle \boldsymbol {v}_{\perp } \boldsymbol {\cdot } \boldsymbol {A}_{\perp } \rangle _{\boldsymbol {R}_{\nu }}/c$

(2.23)\begin{equation} \mathcal{F}_{\boldsymbol{k}} \left[ \left\langle \chi \right\rangle_{\boldsymbol{R}_{\nu}} \right] \doteq \hat{\chi}_{\boldsymbol{k}, \nu} = \left[ {\rm J}_{0,\boldsymbol{k}, \nu} \hat{\phi}_{\boldsymbol{k}} - \frac{v_{{\parallel}}}{c} {\rm J}_{0,\boldsymbol{k}, \nu} \hat{A}_{{\parallel}, \boldsymbol{k}} + 2 \frac{\mu_{\nu}}{Z_{\nu} e} \frac{{\rm J}_{1,\boldsymbol{k},\nu} }{a_{\boldsymbol{k},\nu }} \delta \hat{B}_{{\parallel}, \boldsymbol{k}} \right] , \end{equation}

where we define $\mathcal {F}_{\boldsymbol {k}} [ \langle \chi \rangle _{\boldsymbol {R}_{\nu }} ] \doteq \hat {\chi }_{\boldsymbol {k}, \nu }$ to be the Fourier components of $\langle \chi \rangle _{\boldsymbol {R}_{\nu }}$. The variables $\textrm {J}_{n, \boldsymbol {k}, \nu }$ above are the $n$th-order Bessel functions of the first kind for species $\nu$. The Bessel functions arise naturally as a result of the gyroaverages that appear in (2.14), and have argument $a_{\boldsymbol {k}, \nu } = k_{\perp } v_{\perp } / \varOmega _{\nu }$. Utilising this, we can write each of our functional operators in terms of $\{ \hat {g}_{\boldsymbol {k},\nu }, \hat {\phi }_{\boldsymbol {k}}, \hat {A}_{\parallel, \boldsymbol {k}}, \delta \hat {B}_{\parallel, \boldsymbol {k}} \}$. The resulting gyrokinetic equation is

(2.24)\begin{align} \hat{G}_{\boldsymbol{k},\nu} & = \frac{\partial{\hat{g}_{\boldsymbol{k},\nu}}}{\partial{t}} + v_{{\parallel}} \hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} z \left[\frac{\partial{\hat{g}_{\boldsymbol{k},\nu}}}{\partial{z}} + \frac{Z_{\nu} e}{T_{\nu}} \frac{\partial{\hat{\chi}_{\boldsymbol{k}, \nu}}}{\partial{z}} F_{0,\nu} \right] \nonumber\\ & \quad + {\rm i} \omega_{*,\boldsymbol{k},\nu} F_{0,\nu} \hat{\chi}_{\boldsymbol{k}, \nu} + {\rm i} \omega_{d,\boldsymbol{k},\nu} \left[ \hat{g}_{\boldsymbol{k},\nu} + \frac{Z_{\nu} e}{T_{\nu}} \hat{\chi}_{\boldsymbol{k}, \nu} F_{0,\nu} \right] \nonumber\\ & \quad - \frac{\mu_{\nu}}{m_{\nu}} (\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B_0) \frac{\partial{\hat{g}_{\boldsymbol{k},\nu}}}{\partial{v_{{\parallel}}}} + \frac{Z_{\nu} e}{T_{\nu}} \frac{\mu_{\nu}}{m_{\nu} c } (\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B_0) F_{0,\nu} {\rm J}_{0,\boldsymbol{k}, \nu} \hat{A}_{{\parallel}, \boldsymbol{k}} - \hat{C}_{\boldsymbol{k},\nu} [\hat{g}_{\boldsymbol{k},\nu}] \nonumber\\ & = 0. \end{align}

The drift frequencies, $\omega _{d,\nu }$ and $\omega _{*,\nu }$, correspond to the magnetic drift frequencies resulting from the gradient and curvature of the magnetic field and the diamagnetic drift respectively, and are given by

(2.25)\begin{gather} \omega_{d, \boldsymbol{k}, \nu} = \frac{1}{\varOmega_{\nu}} (v_{{\parallel}}^2 \boldsymbol{v}_{\kappa} + \mu_{\nu} \boldsymbol{v}_{\boldsymbol{\nabla} B} ) \boldsymbol{\cdot} ( k_x \boldsymbol{\nabla} x + k_y \boldsymbol{\nabla} y ) , \end{gather}

(2.26)\begin{gather}\omega_{*,\boldsymbol{k}, \nu} = \frac{c k_y } {B_{0} } \frac{{\rm d} y}{{\rm d} \alpha} \frac{{\rm d} \ln F_{0,\nu}}{{\rm d} \psi} , \end{gather}

with $\boldsymbol {v}_{\kappa } = \hat {\boldsymbol {b}} \times ( \hat {\boldsymbol {b}} \boldsymbol {\cdot } \boldsymbol {\nabla } \hat {\boldsymbol {b}} )$, and $\boldsymbol {v}_{\boldsymbol {\nabla } B} = \hat {\boldsymbol {b}} \times \boldsymbol {\nabla } B$.

The distribution function and field quantities are time-dependent functions, and their evolution equations admit normal mode solutions. For a generic set of initial conditions, the fastest growing (or slowest decaying) mode is likely to have an initial amplitude that is comparable to or smaller than other normal modes in the system. This results in an initial transient period due to the competition between these different normal modes, followed by exponential growth or decay once the fastest growing (or slowest decaying) mode has a much larger amplitude than all others. Therefore, we decompose $\hat {f}_{\boldsymbol {k}} = \sum _{j} \bar {f}_{\boldsymbol {k},j} \textrm {e}^{\gamma _{\boldsymbol {k},j} t}$Footnote ⁵ with $\gamma _{\boldsymbol {k},j} \in \mathbb {C}$ the complex frequencies, and $\bar {f}_{\boldsymbol {k},j}$ the amplitude for the $j$th normal mode.Footnote ⁶ Quasineutrality and Ampère's law, along with the linearity of the gyrokinetic equation, ensure that all fluctuating quantities for an equivalent normal mode share the same complex frequency and will thus exhibit the same time-dependent behaviour during the period of exponential growth or decay. After a sufficiently long period of time (the exact length of which will depend on both the relative growth rates and starting amplitudes of each mode) the fastest growing mode will dominate, meaning we can approximate the time dependent behaviour after large times using a single temporal mode; $\hat {f}_{\boldsymbol {k}} \approx \bar {f}_{\boldsymbol {k},0} \textrm {e}^{\gamma _{\boldsymbol {k}, 0} \tilde {t}}$, with $\mathrm {Re}(\gamma _{\boldsymbol {k},0}) > \mathrm {Re}(\gamma _{\boldsymbol {k},j})$, $\forall j$, $j \neq 0$. Using this normal mode decomposition in the gyrokinetic equation and considering behaviour beyond the transient period gives

(2.27)\begin{align} \hat{G}_{\boldsymbol{k},\nu} & = \gamma_{\boldsymbol{k}, 0} \bar{g}_{\boldsymbol{k} , \nu, 0}+ v_{{\parallel}} \hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} z \left[\frac{\partial{\bar{g}_{\boldsymbol{k} , \nu, 0}}}{\partial{z}} + \frac{Z_{\nu} e}{T_{\nu}} \frac{\partial{\bar{\chi}_{\boldsymbol{k}, 0}}}{\partial{z}} F_{0,\nu} \right] \nonumber\\ & \quad + {\rm i} \omega_{*,\boldsymbol{k},\nu} F_{0,\nu} \bar{\chi}_{\boldsymbol{k}, 0} + {\rm i} \omega_{d,\boldsymbol{k},\nu} \left[ \bar{g}_{\boldsymbol{k} , \nu, 0} + \frac{Z_{\nu} e}{T_{\nu}} \bar{\chi}_{\boldsymbol{k}, 0} F_{0,\nu} \right] \nonumber\\ & \quad - \mu_{\nu} \hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B_0 \frac{\partial{\bar{g}_{\boldsymbol{k} , \nu, 0}}}{\partial{v_{{\parallel}}}} + \frac{Z_{\nu} e}{T_{\nu}} \frac{\mu_{\nu}}{m_{\nu} c } (\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} B_0) F_{0,\nu} {\rm J}_{0,\boldsymbol{k}, \nu} \bar{A}_{{\parallel}, \boldsymbol{k}, 0} - \hat{C}_{\boldsymbol{k},\nu} [\bar{g}_{\boldsymbol{k} , \nu, 0}]. \end{align}

The corresponding transformed, normalised field equations are $\hat {Q}_{\boldsymbol {k}} = \hat {M}_{\boldsymbol {k}} = \hat {N}_{\boldsymbol {k}} = 0$, where

(2.28)\begin{align} \hat{Q}_{\boldsymbol{k}} & = \sum_{\nu} Z_{\nu} e \left\{ \frac{2 {\rm \pi}B_0}{m_{\nu} } \int{\mathrm{d}^2 v \, {\rm J}_{0,\boldsymbol{k}, \nu} \bar{g}_{\boldsymbol{k} , \nu, 0} } + \frac{Z_{\nu} e n_{\nu} }{T_{\nu}} \left( \varGamma_{0, \boldsymbol{k}, \nu} - 1 \right) \bar{\phi}_{\boldsymbol{k}, 0} + \frac{n_{\nu} }{B_0} \varGamma_{1, \boldsymbol{k}, \nu} \delta \bar{B}_{{\parallel}, \boldsymbol{k}, 0} \right\} , \end{align}

(2.29)\begin{align} \hat{M}_{\boldsymbol{k}} & =- \frac{4 {\rm \pi}}{k_{{\perp}}^2 c} \sum_{\nu} Z_{\nu} e \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v v_{{\parallel}} {\rm J}_{0,\boldsymbol{k}, \nu} \bar{g}_{\boldsymbol{k} , \nu, 0} \nonumber\\ & \quad + \left[1 + \frac{4 {\rm \pi}}{ k_{{\perp}}^2 c^2} \sum_{\nu} \frac{ (Z_{\nu} e)^2 n_{\nu}}{m_{\nu}} \varGamma_{0,\boldsymbol{k}, \nu} \right] \bar{A}_{{\parallel}, \boldsymbol{k}, 0} , \end{align}

and

(2.30) \begin{align} &\hat{N}_{\boldsymbol{k}} = 8 {\rm \pi}\sum_{\nu} \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v \frac{{\rm J}_{1, \boldsymbol{k}, \nu} }{a_{\boldsymbol{k}, \nu}} \mu_{\nu} \bar{g}_{\boldsymbol{k} , \nu, 0} + \left[ 4 {\rm \pi} \sum_{\nu} \frac{Z_{\nu} e n_{\nu}}{B_0} \varGamma_{1,\boldsymbol{k}, \nu} \right] \bar{\phi}_{\boldsymbol{k}, 0} \nonumber\\ & + \left[1+ 16 {\rm \pi}\sum_{\nu} \frac{ n_{\nu} T_{\nu}}{B_0^2} \varGamma_{2,\boldsymbol{k}, \nu} \right] \delta \bar{B}_{{\parallel}, \boldsymbol{k}, 0} , \end{align}

with $\int \,\mathrm {d}^2 v \doteq \int \,\mathrm {d} \mu \int \,\mathrm {d} v_{\parallel }$. In the above we have defined

(2.31)\begin{gather} \varGamma_{0,\boldsymbol{k}, \nu} = {\rm I}_0 (\alpha_{\boldsymbol{k},\nu} ) {\rm e}^{- \alpha_{\boldsymbol{k},\nu} }, \end{gather}

(2.32)\begin{gather}\varGamma_{1,\boldsymbol{k},\nu} = [{\rm I}_0(\alpha_{\boldsymbol{k},\nu} ) - {\rm I}_1(\alpha_{\boldsymbol{k},\nu} ) ] {\rm e}^{- \alpha_{\boldsymbol{k},\nu} }, \end{gather}

(2.33)\begin{gather}\varGamma_{2,\boldsymbol{k},\nu} = {\rm I}_1 (\alpha_{\boldsymbol{k},\nu}) {\rm e}^{- \alpha_{\boldsymbol{k},\nu} }, \end{gather}

where $\textrm {I}_0$ and $\textrm {I}_1$ are modified Bessel functions of the first kind, and $\alpha _{\boldsymbol {k}, \nu } = k_{\perp }^2 \rho _{\nu }^2 /2$.

3.1 General formalism

Since our equations are limited to the linear regime, there is no coupling of different Fourier modes and it is thus possible to consider each perpendicular wavenumber individually. Given that we are also considering the post-transient limit with only one dominant growth rate, it will henceforth be assumed that only a single perpendicular wavenumber is being considered, and we will thus drop the $\boldsymbol {k}$ subscript, along with the subscript that denotes the dominant mode for the distribution function and field quantities. For notational brevity we shall also drop the overbars that appear on the distribution function and fields that denoted normal mode decomposition. Hence, everywhere $g_{\nu }$, $\phi$, $A_{\parallel }$ and $\delta B_{\parallel }$ appear it shall be assumed that they contain suppressed spatial and temporal Fourier subscripts.

Consider now a set of variables that influences our system and which exists within a parameter space spanned by all possible $\boldsymbol {p}$. We take the set $\{ p_i\}$, $i \in [1, \mathcal {N}_p ]$, to be linearly independent, with no time variation, and explore how variations within this space affect the linear growth rate. At present we need not specify which variables are denoted by $\boldsymbol {p}$, and thus derive a general set of adjoint equations for the gyrokinetic-Maxwell system above, (2.27)–(2.30).

Our objective functional, $\hat {G}_{\nu }$, is a linear functional of $\{g_{\nu }, \phi, A_{\parallel }, \delta B_{\parallel } \}$, which are coupled to the fields through the field equations, (2.28)–(2.30). When taking derivatives of $G_{\nu }$ we invariably end up with derivatives acting on all four of these variables. This is undesirable because it requires we calculate the gradients of $\hat {g}_{\nu }$ and the fields, which in turn requires us to solve the gyrokinetic system $\mathcal {N}_p + 1$ times, as discussed above. In order to eliminate these four derivatives we use something akin to the method of Lagrange multipliers and introduce four adjoint variables which multiply the corresponding constraint equations. The optimisation Lagrangian is thus

(3.1)\begin{equation} \mathcal{L} := \left\langle \hat{G}_{\nu}, \lambda_{\nu}\right\rangle_{z, v, {\nu}} + \left\langle \hat{Q}, \xi\right\rangle_{z}+ \left\langle \hat{M}, \zeta \right\rangle_{z}+ \left\langle \hat{N}, \sigma\right\rangle_{z} , \end{equation}

with the angle brackets representing inner products defined throughFootnote ⁷

(3.2)\begin{equation} \left. \begin{aligned} \left\langle a, b\right\rangle_{z} & = \int{ \frac{\mathrm{d} z }{B_0 \,\hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z} a b^* } , \quad \left\langle a, b\right\rangle_{v, {\nu}} = \sum_{\nu} \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v a b^* ,\\ \left\langle a, b\right\rangle_{z, v, {\nu}} & = \sum_{\nu} \int \frac{\mathrm{d} z}{B_0 \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z } \frac{2{\rm \pi} B_0}{m_{\nu} } \int \,\mathrm{d}^2 v a b^* . \end{aligned} \right\} \end{equation}

We have introduced the set of adjoint variables $\lambda _{\nu }$, $\xi$, $\zeta$ and $\sigma$, whose forms are to be determined.Footnote ⁸ We identify $\lambda _{\nu }$ as the adjoint variable to the distribution function, $g_{\nu }$, whereas $\xi$, $\zeta$ and $\sigma$ are adjoint to the field variables (referred to henceforth as adjoint fields). The quantities $\hat {G}_{\nu }$, $\hat {Q}$, $\hat {M}$ and $\hat {N}$ are our objective functions, and we remark that since $\hat {G}_{\nu }= \hat {Q} = \hat {M} = \hat {N} = 0$ for a consistent set of $\{ \boldsymbol {p}_0, g_{\nu }(\boldsymbol {p}_0), \phi (\boldsymbol {p}_0), A_{\parallel } (\boldsymbol {p}_0), \delta B_{\parallel } (\boldsymbol {p}_0)\}$ we have $\mathcal {L}|_{\boldsymbol {p}_0} = 0$ for all choices of adjoint variables. For later convenience we decompose the functional operators into their components that act on $g_{\nu }$, $\phi$, $A_{\parallel }$ and $\delta B_{\parallel }$ separately

(3.3)\begin{align} \hat{G}_{\nu} [\boldsymbol{p}; g_{\nu}, \phi, A_{{\parallel}}, \delta B_{{\parallel}}] & = \hat{G}_{g, \nu} [\boldsymbol{p} ; g_{\nu}] + \hat{G}_{\phi, \nu} [\boldsymbol{p} ; \phi] + \hat{G}_{A_{{\parallel}}, \nu} [\boldsymbol{p} ; A_{{\parallel}}] + \hat{G}_{B_{{\parallel}}, \nu} [\boldsymbol{p} ;\delta B_{{\parallel}}]\nonumber\\ & \quad - \hat{C}_{\nu}[\boldsymbol{p};g_{\nu}], \end{align}

(3.4)\begin{align}\hat{Q} [\boldsymbol{p} ; g_{\nu}, \phi, \delta B_{{\parallel}}] & = \left\langle \hat{Q}_{g, \nu} [\boldsymbol{p} ; g_{\nu}], \mathbb{I}\right\rangle_{v, {\nu}} + \hat{Q}_{\phi} [\boldsymbol{p} ; \phi ] + \hat{Q}_{B_{{\parallel}}} [\boldsymbol{p} ; \delta B_{{\parallel}}], \end{align}

(3.5)\begin{align}\hat{M} [\boldsymbol{p} ; g_{\nu}, A_{{\parallel}}] & = \left\langle \hat{M}_{g,\nu} [\boldsymbol{p} ; g_{\nu}] , \mathbb{I} \right\rangle_{v, {\nu}} + \hat{M}_{A_{{\parallel}}} [\boldsymbol{p} ; A_{{\parallel}} ], \end{align}

(3.6)\begin{align}\hat{N} [\boldsymbol{p} ; g_{\nu}, \phi, \delta B_{{\parallel}} ] & = \left\langle \hat{N}_{g,\nu} [\boldsymbol{p} ; g_{\nu} ], \mathbb{I}\right\rangle_{v, {\nu}} + \hat{N}_{\phi} [\boldsymbol{p} ; \phi ] + \hat{N}_{B_{{\parallel}}} [\boldsymbol{p}; \delta B_{{\parallel}}] , \end{align}

with $\mathbb {I}$ simply equal to one. Explicit expressions for these operators are given in Appendix A.

We next consider taking the gradient of (3.1) with respect to the variables in our parameter space, $\boldsymbol {p}$. We now isolate all terms multiplying derivatives of $\{g_{\nu }, \phi, A_{\parallel }, \delta B_{\parallel } \}$ so that each of their coefficients can be set to zero. We are at liberty to do this because of the freedom that exists in choosing the adjoint variables introduced in (3.1). The gradient of (3.1) is expanded using (3.3)–(3.6)

(3.7)\begin{align} \boldsymbol{\nabla}_{\boldsymbol{p}} \mathcal{L} & = \partial_{\boldsymbol{p}} \mathcal{L} + \left\langle \boldsymbol{\nabla}_{\boldsymbol{p}} g_{\nu} , \hat{G}^{{\dagger}}_{g,\nu} [\boldsymbol{p}; \lambda_{\nu}] - \hat{C}^{{\dagger}}_{\nu} [\boldsymbol{p} ; \lambda_{\nu}] + \hat{Q}^{{\dagger}}_{g,\nu} [\boldsymbol{p}; \xi] + \hat{M}^{{\dagger}}_{g,\nu} [\boldsymbol{p}; \zeta] + \hat{N}^{{\dagger}}_{g,\nu} [\boldsymbol{p}; \sigma]\right\rangle_{z, v, {\nu}} \nonumber\\ & \quad + \left\langle \boldsymbol{\nabla}_{\boldsymbol{p}} \phi , \hat{G}^{{\dagger}}_{\phi,\nu} [\boldsymbol{p} ; \lambda_{\nu}] + \hat{Q}^{{\dagger}}_{\phi} [\boldsymbol{p} ; \xi] + \hat{N}^{{\dagger}}_{\phi} [\boldsymbol{p} ; \sigma] \right\rangle_{z} + \left\langle \boldsymbol{\nabla}_{\boldsymbol{p}} A_{{\parallel}} , \hat{G}^{{\dagger}}_{A_{{\parallel}},\nu} [\boldsymbol{p} ; \lambda_{\nu}] + \hat{M}^{{\dagger}}_{A_{{\parallel}}} [\boldsymbol{p} ; \zeta]\right\rangle_{z} \nonumber\\ & \quad + \left\langle \boldsymbol{\nabla}_{\boldsymbol{p}} \delta B_{{\parallel}} , \hat{G}^{{\dagger}}_{B_{{\parallel}},\nu} [\boldsymbol{p} ; \lambda_{\nu}] + \hat{Q}^{{\dagger}}_{B_{{\parallel}}} [\boldsymbol{p} ; \xi] + \hat{N}^{{\dagger}}_{B_{{\parallel}}} [\boldsymbol{p} ; \sigma] \right\rangle_{z} + \mathcal{B} . \end{align}

Here, the daggers that appear on the functional operators denote the adjoint of those operators with respect to the relevant inner product. The partial derivative, $\partial _{\boldsymbol {p}}$, is taken at fixed $g_{\nu }$, $\phi$, $A_{\parallel }$ and $\delta B_{\parallel }$. Note that it acts on the inner product itself in addition to the terms within it to account for the $\boldsymbol {p}$-dependence of the JacobiansFootnote ⁹ and functional operators. The term ‘$\mathcal {B}$’ accounts for boundary terms that arise when we integrate by parts to invert operators such as $\left \langle \boldsymbol {\nabla }_{\boldsymbol {p}} \partial _z g_{\nu } , \lambda _{\nu }\right \rangle _{z}$ onto the adjoint variables, $\left \langle \boldsymbol {\nabla }_{\boldsymbol {p}} g_{\nu }, \partial _z \lambda _{\nu }\right \rangle _{z}$, and is given by

(3.8)\begin{align} \mathcal{B} & = \sum_{\nu} \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v v_{{\parallel}} \lambda_{\nu}^* \left[ \boldsymbol{\nabla}_{\boldsymbol{p}} g_{\nu} + \frac{Z_{\nu} e}{T_{\nu} } {\rm J}_{0,\nu} F_{0,\nu} \boldsymbol{\nabla}_{\boldsymbol{p}} \phi \right. \nonumber\\ & \quad \left. \left.- 2 \frac{Z_{\nu} e}{T_{\nu} } v_{{\parallel}} {\rm J}_{0,\nu} F_{0,\nu} \boldsymbol{\nabla}_{\boldsymbol{p}} A_{{\parallel}} + 4 \mu_{\nu} \frac{{\rm J}_{1,\nu}}{a_{\nu}} F_{0,\nu} \boldsymbol{\nabla}_{\boldsymbol{p}} \delta B_{{\parallel}} \right] \right|_{z =-\infty}^{z = \infty} \nonumber\\ & \quad \left. - \sum_{\nu} \frac{2{\rm \pi} B_0}{m_{\nu} } \int \,\mathrm{d} z \int \,\mathrm{d} \mu\, \mu_{\nu} \frac{\partial{B_0}}{\partial{z}} \lambda_{\nu}^* \boldsymbol{\nabla}_{\boldsymbol{p}} g_{\nu} \right|_{v_{{\parallel}} =- \infty}^{v_{{\parallel}} = \infty} . \end{align}

We can conveniently set these terms to zero by applying appropriate restrictions on our adjoint variables. These terms consequently define the boundary conditions we apply to our adjoint variables. The incoming boundary conditions along the magnetic field on $g_{\nu }$ are taken to be $g_{\nu }(z\rightarrow -\infty, v_{\parallel } > 0, \mu _{\nu }), g_{\nu }(z \rightarrow \infty, v_{\parallel } <0, \mu _{\nu }) \rightarrow 0$, independently of $\boldsymbol {p}$, such that $d_{\boldsymbol {p}} g_{\nu } = 0$ at these limits. We impose the boundary condition on $\lambda _{\nu }^*$ to be $\lambda _{\nu }^* (z \rightarrow -\infty, v_{\parallel } <0, \mu _{\nu }), \lambda _{\nu }^* (z \rightarrow \infty,v_{\parallel }>0, \mu _{\nu }) \rightarrow 0$ in order to eliminate the boundary term arising from the $z$ integration by parts, and hence remove the need to calculate $\boldsymbol {\nabla }_{\boldsymbol {p}} \{g_{\nu },\phi,A_{\parallel },\delta B_{\parallel } \}$ at the boundaries in $z$. The boundary term arising from integration by parts in $v_{\parallel }$ is automatically satisfied as it is assumed that $g_{\nu }(z, v_{\parallel } \rightarrow \pm \infty, \mu _{\nu }) = 0$, $\forall \{z, \mu _{\nu } \}$ independently of $\boldsymbol {p}$. However, it is convenient to impose that $\lambda _{\nu }^* (z, v_{\parallel } \rightarrow \pm \infty, \mu _{\nu }) = 0$ such that $\lambda _{\nu }^*$ and $g_{\nu }$ satisfy similar boundary conditions, whilst also ensuring $\lambda _{\nu }$ is sensibly defined and normalisable.Footnote ¹⁰ The substitution $\lambda ^{\leftrightarrow, *}_{\nu } = \lambda _{\nu }^*(z,-v_\parallel,\mu _{\nu })$ is made such that the adjoint equations more closely resemble those in the original gyrokinetic system. This redefines the $z$-boundary condition on the adjoint variable $\lambda ^{\leftrightarrow, *}_{\nu } (z \rightarrow - \infty, v_{\parallel } > 0, \mu _{\nu }), \lambda ^{\leftrightarrow, *}_{\nu } (z \rightarrow \infty, v_{\parallel } < 0, \mu _{\nu }) \rightarrow 0$, which now mirrors those satisfied by $g_{\nu }$.

Setting the remaining coefficients of $\boldsymbol {\nabla }_{\boldsymbol {p}} g_{\nu }$, $\boldsymbol {\nabla }_{\boldsymbol {p}} \phi$, $\boldsymbol {\nabla }_{\boldsymbol {p}} A_{\parallel }$ and $\boldsymbol {\nabla }_{\boldsymbol {p}} \delta B_{\parallel }$ in (3.7) equal to zero yields the constraint equations for the adjoint variables

(3.9)\begin{gather} \hat{G}^{{\dagger}}_{g,\nu} [ \boldsymbol{p} ; \lambda^{\leftrightarrow}_{\nu}] + \hat{Q}^{{\dagger}}_{g,\nu} [\boldsymbol{p} ; \xi] + \hat{M}^{{\dagger}}_{g,\nu} [\boldsymbol{p}; \zeta ] + \hat{N}^{{\dagger}}_{g,\nu} [\boldsymbol{p} ; \sigma] - \hat{C}^{{\dagger}}_{\nu} [\boldsymbol{p} ; \lambda^{\leftrightarrow}_{\nu}] = 0 , \end{gather}

(3.10)\begin{gather}\langle \hat{G}^{{\dagger}}_{\phi,\nu} [\boldsymbol{p} ; \lambda^{\leftrightarrow}_{\nu} ] \rangle_{v, {\nu}} + \hat{Q}^{{\dagger}}_{\phi} [ \boldsymbol{p} ; \xi ] + \hat{N}^{{\dagger}}_{\phi} [\boldsymbol{p} ; \sigma] = 0 , \end{gather}

(3.11)\begin{gather}\langle \hat{G}^{{\dagger}}_{A_{{\parallel}},\nu} [ \boldsymbol{p} ; \lambda^{\leftrightarrow}_{\nu}] \rangle_{v, {\nu}} + \hat{M}^{{\dagger}}_{A_{{\parallel}}} [\boldsymbol{p} ; \zeta] = 0 , \end{gather}

(3.12)\begin{gather}\langle \hat{G}^{{\dagger}}_{B_{{\parallel}},\nu} [ \boldsymbol{p} ; \lambda^{\leftrightarrow}_{\nu}] \rangle_{v, {\nu}} + \hat{Q}^{{\dagger}}_{B_{{\parallel}}} [\boldsymbol{p} ; \xi] + \hat{N}^{{\dagger}}_{B_{{\parallel}}} [ \boldsymbol{p} ; \sigma] = 0 . \end{gather}

The expressions for these adjoint operators are stated in Appendix B. The $z$-derivatives which appear in the $\hat {G}_{\nu }$ operators, under the velocity integrals, in (3.10)–(3.12) pose a potential difficulty; to calculate the adjoint fields, information is required for $\lambda _{\nu }^*$ at all $z$, for both positive and negative velocities. Given the boundary conditions on $\lambda _{\nu }$ and the propagation of information by advection this information is not readily available. This is akin to the problem faced when solving the gyrokinetic system; an incoming boundary condition is imposed on the distribution function at $z \rightarrow \pm \infty$ when $v_{\parallel } <0$ and $v_{\parallel } > 0$, respectively. The outgoing boundary information is not known a priori but must instead be solved for. Anticipating the difficulties this will create for numerical solution, we take moments of (3.9) are taken to simplify the adjoint equations. This simplifies the adjoint equations, bringing them into a form more closely resembling the that of the gyrokinetic field equation. It should be emphasised that we retain a fully kinetic treatment in doing so.

A summary of this calculation can be found in Appendix C, and the result after algebraic manipulation is

(3.13)\begin{gather} \gamma^* \lambda^{\leftrightarrow}_{\nu} + v_{{\parallel}} \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z \frac{\partial{\lambda^{\leftrightarrow}_{\nu}}}{\partial{z}} - \frac{\mu_{\nu}}{m_{\nu}} \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z \frac{\partial{B_0}}{\partial{z}} \frac{\partial{\lambda^{\leftrightarrow}_{\nu}}}{\partial{v_{{\parallel}}}} - {\rm i} \omega_{d, \nu} \lambda^{\leftrightarrow}_{\nu} + Z_{\nu} e {\rm J}_{0,\nu} \xi\nonumber\\ - \frac{4 {\rm \pi}}{k_{{\perp}}^2} Z_{\nu} e {\rm J}_{0,\nu} \frac{v_{{\parallel}}}{c} \zeta + 8 {\rm \pi}\frac{{\rm J}_{1,\nu}}{a_{\nu}} \mu_{\nu} \sigma - \hat{C}_{\nu}[\lambda^{\leftrightarrow}_{\nu}] = 0 , \end{gather}

(3.14)\begin{gather}\xi + \frac{1}{\eta} \sum_{\nu} \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v \left[{\rm i} \omega_{*,\nu} + \frac{Z_{\nu}e}{T_{\nu}} \gamma^* \right] {\rm J}_{0,\nu} F_{0,\nu} \lambda^{\leftrightarrow}_{\nu} =0 , \end{gather}

(3.15)\begin{gather}\zeta - \frac{1}{k_{{\perp}}^2 } \sum_{\nu} \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v \frac{v_{{\parallel}}}{c} \left[ {\rm i} \omega_{*,\nu} + \frac{Z_{\nu}e}{T_{\nu}} \gamma^* \right] {\rm J}_{0,\nu} F_{0,\nu} \lambda^{\leftrightarrow}_{\nu} = 0 , \end{gather}

(3.16)\begin{gather}\sigma - \sum_{\nu} \frac{2 {\rm \pi}B_0}{m_{\nu} } \int \,\mathrm{d}^2 v \left( 2 \frac{\mu_{\nu} }{Z_{\nu} e}\frac{{\rm J}_{1,\nu}}{a_{\nu}} \right) \left[ {\rm i} \omega_{*,\nu} + \frac{Z_{\nu} e}{T_{\nu} } \gamma^* \right] F_{0,\nu} \lambda^{\leftrightarrow}_{\nu} = 0 , \end{gather}

with $\eta = \sum _{\nu } (Z_{\nu } e)^2 n_{\nu }/ T_{\nu }$. Noting that we can also rewrite $\hat {G}_{\nu } [\boldsymbol {p}; g_{\nu }, \phi, A_{\parallel }, \delta B_{\parallel }]= \gamma g_{\nu } + \hat {L}_{\nu } [\boldsymbol {p}; g_{\nu }, \phi, A_{\parallel }, \delta B_{\parallel }]$, and using $\boldsymbol {\nabla }_{\boldsymbol {p}} \mathcal {L} |_{\boldsymbol {p}_0} = 0$, we can rearrange equation (3.7) to obtain

(3.17)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{p}} \gamma \langle g_{\nu}, \lambda_{\nu} \rangle_{z,v, {\nu}} =- \left. \left[ \langle \partial_{\boldsymbol{p}} \hat{L}_{\nu} , \lambda_{\nu} \rangle_{z,v, {\nu}} + \langle \partial_{\boldsymbol{p}} \hat{Q}, \xi \rangle_{z} + \langle \partial_{\boldsymbol{p}} \hat{M} , \zeta \rangle_{z} + \langle \partial_{\boldsymbol{p}} \hat{N}, \sigma \rangle_{z} \right] \right|_{\boldsymbol{p}_0} , \end{equation}

where the partial derivatives have been pulled inside the inner products as the contribution arising from the Jacobian derivatives is zero by virtue of $\hat {G}_{\nu }(\boldsymbol {p}_0) = \hat {Q}(\boldsymbol {p}_0) = \hat {M}(\boldsymbol {p}_0) = \hat {N}(\boldsymbol {p}_0) = 0$.

To solve for the derivative of the linear growth rate, the following procedure is taken: the gyrokinetic equations, (2.27)–(2.30), are solved to obtain $g_{\nu }$ and $\gamma$, and then the equations (3.13)–(3.16) are used to solve for the adjoint variables. These quantities are then all fed into (3.17) to compute $\boldsymbol {\nabla }_{\boldsymbol {p}} \gamma$.

3.2 Electrostatic, collisionless limit

We consider now the electrostatic, collisionless limit, as we shall perform numerical tests in this regime. In the limit of small plasma $\beta$, meaning the ratio of plasma to magnetic pressure tends to zero, the magnetic field perturbations tend to zero. Hence, to extract the electrostatic, collisionless limit from these equations we set $A_{\parallel } = \delta B_{\parallel } = 0$, and $C_{\nu, \nu '} =0$. The linear, collisionless, electrostatic gyrokinetic equation is

(3.18)\begin{align} \hat{G}_{\nu} [\boldsymbol{p}; g_{\nu}, \phi] & = \gamma g_{\nu} + v_{{\parallel}} \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z \left[\frac{\partial{g_{\nu}}}{\partial{z}} + \frac{Z_{\nu} e}{T_{\nu} } \frac{\partial{{\rm J}_{0,\nu} \phi}}{\partial{z}} F_{0,\nu} \right] + {\rm i} \omega_{*,\nu} {\rm J}_{0,\nu} \phi F_{0,\nu}\nonumber\\ & \quad - \frac{\mu_{\nu}}{m_{\nu} } \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z \frac{\partial{B_0}}{\partial{z}} \frac{\partial{g_{\nu}}}{\partial{v_{{\parallel}}}} + {\rm i} \omega_{d,\nu} \left[ g_{\nu} + \frac{Z_{\nu} e}{T_{\nu} } {\rm J}_{0,\nu} \phi F_{0,\nu} \right], \end{align}

which is closed by the electrostatic limit of quasineutrality

(3.19)\begin{equation} \hat{Q} [\boldsymbol{p};g_{\nu}, \phi] = \sum_{\nu}Z_{\nu} e \left[ \frac{2{\rm \pi} B_0}{m_{\nu} } \int \,{\rm d}^2 v {\rm J}_{0,\nu} g_{\nu} + \frac{Z_{\nu} n_{\nu}}{T_{\nu}} (\varGamma_{0,\nu} -1) \phi \right] , \end{equation}

with $\hat {G}_{\nu }$ and $\hat {Q}$ identically zero, and $\hat {M}$, $\hat {N}$ providing no contribution in the electrostatic limit. In the above we are once again considering the long time behaviour of a single wavenumber, and have suppressed the associated subscripts.

As in § 3.1 we decompose the functional operators $\hat {G}_{\nu }[\boldsymbol {p}; g_{\nu }, \phi ]$ and $\hat {Q}[\boldsymbol {p}; g_{\nu }, \phi ]$ into components that act on $g_{\nu }$ and $\phi$ separately, with all other operators in (3.3)–(3.6) set to zero. The derivation in § 3.1 is unchanged, with the exception that some terms may now be omitted. The resulting derivative of the growth rate in the electrostatic, collisionless regime is

(3.20)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{p}} \gamma \langle g_{\nu}, \lambda_{\nu} \rangle_{z,v, {\nu}} =- \left. \left[ \langle \partial_{\boldsymbol{p}} \hat{L}_{\nu}, \lambda_{\nu} \rangle_{z,v, {\nu}} + \langle \partial_{\boldsymbol{p}} \hat{Q}, \xi \rangle_{z} \right] \right|_{\boldsymbol{p}_0}, \end{equation}

where $\hat {L}_{\nu }= \hat {G}_{\nu } - \gamma g_{\nu }$ is given by equation (3.18), and the adjoint equations are

(3.21)\begin{gather} \gamma^* \lambda^{\leftrightarrow}_{\nu} + v_{{\parallel}} \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z \frac{\partial{\lambda^{\leftrightarrow}_{\nu}}}{\partial{z}} - \frac{\mu_{\nu}}{m_{s}} \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} z \frac{\partial{B_0}}{\partial{z}} \frac{\partial{\lambda^{\leftrightarrow}_{\nu}}}{\partial{v_{{\parallel}}}} - {\rm i} \omega_{d, \nu} \lambda^{\leftrightarrow}_{\nu} + Z_{\nu} e {\rm J}_{0,\nu} \xi = 0 , \end{gather}

(3.22)\begin{gather}\xi + \frac{1}{\eta} \sum_{\nu} \frac{2 {\rm \pi}B_0} {m_{\nu} } \int \,\mathrm{d}^2 v \left[{\rm i} \omega_{*,\nu} + \frac{Z_{\nu} e}{T_{\nu} } \gamma^* \right] {\rm J}_{0,\nu} F_{0,\nu} \lambda^{\leftrightarrow}_{\nu} = 0 , \end{gather}

with $\lambda ^{\leftrightarrow }_{\nu } (v_{\parallel } ) = \lambda _{\nu }(-v_{\parallel } )$, and $\eta = \sum _{\nu } (Z_{\nu } e)^2 n_{\nu }/ T_{\nu }$ as before.

4.1 Normalisation

4.1.1 Gyrokinetic normalisations

Here, we normalise the full electromagnetic gyrokinetic system and give the reduced electrostatic limit at the end of the section. We choose our reference quantities to be the same as in stella to aid in numerical implementation, and denote these with a subscript ‘$r$’. A reference length scale that characterises the parallel lengths in the simulation is introduced as $L=a$; for the Miller formalism stella takes $a$ to be the half diameter of the plasma volume (minor radius), at the height of the magnetic axis. The perpendicular reference length scale is taken as a reference gyroradius $\rho _r \doteq v_{\textrm {th},r}/\varOmega _r$, with $v_{\textrm {th},r} = \sqrt {2 T_r/m_r}$ the reference velocity, and $\varOmega _r = e B_r/m_{r}c$, with $T_r$, $n_r$, $B_r$ and $m_r$ being user specified quantities. The normalised variables are denoted with a tilde, and are provided in table 1.

Table 1.

List of normalised parameters and variables.

Multiplying the gyrokinetic equation, rewritten in normalised coordinates, by the factor $(a^2/\rho _r v_{\textrm {th},r}) \textrm {exp}(-v^2/v_{\textrm {th},\nu }^2)/F_{0,\nu }$, we obtain the normalised low-flow, electromagnetic gyrokinetic equation taken in the long time limit

(4.1)\begin{align} \hat{G}_{\boldsymbol{k},\nu} & = \tilde{\gamma} \tilde{g}_{\boldsymbol{k}, \nu} + \tilde{v}_{{\rm th},\nu} \tilde{v}_{{\parallel}} \, \hat{\boldsymbol{b}}\boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \left[\frac{\partial{\tilde{g}_{\boldsymbol{k}, \nu}}}{\partial{\tilde{z}}} + \frac{Z_{\nu}}{\tilde{T}_{\nu}} \frac{\partial{\left\langle \tilde{\chi} \right\rangle_{\boldsymbol{k}, \nu}}}{\partial{\tilde{z}}} {\rm e}^{-\tilde{v}_{\nu}^2} \right] \nonumber\\ & \quad + {\rm i} \tilde{\omega}_{*,\boldsymbol{k},\nu} {\rm e}^{-\tilde{v}_{\nu}^2} \left\langle \tilde{\chi} \right\rangle_{\boldsymbol{k}, \nu} + {\rm i} \tilde{\omega}_{d,\boldsymbol{k},\nu} \left[ \tilde{g}_{\boldsymbol{k}, \nu} + \frac{Z_{\nu}}{\tilde{T}_{\nu}} \left\langle \tilde{\chi} \right\rangle_{\boldsymbol{k}, \nu} {\rm e}^{-\tilde{v}_{\nu}^2} \right] \nonumber\\ & \quad - \tilde{v}_{\text{th},\nu} \tilde{\mu}_{\nu} \hat{\boldsymbol{b}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{B}_0 \frac{\partial{\tilde{g}_{\boldsymbol{k}, \nu}}}{\partial{\tilde{v}_{{\parallel}}}} + 2 \frac{Z_{\nu}}{\tilde{m}_{\nu}} \tilde{\mu}_{\nu} \hat{\boldsymbol{b}} \boldsymbol{\cdot} \tilde{B}_0 {\rm e}^{-\tilde{v}_{\nu}^2} {\rm J}_{0,\boldsymbol{k}, \nu} \tilde{A}_{{\parallel}, \boldsymbol{k}} - \hat{C}_{\boldsymbol{k},\nu} [\tilde{g}_{\boldsymbol{k}, \nu}]. \end{align}

The corresponding transformed, normalised field equations are given by

(4.2)\begin{align} \hat{Q}_{\boldsymbol{k}} & = \sum_{\nu} Z_{\nu} \tilde{n}_{\nu} \left\{ \frac{2\tilde{B}_0}{\sqrt{\rm \pi}} \int{\,\mathrm{d}^2 \tilde{v} {\rm J}_{0,\boldsymbol{k}, \nu} \tilde{g}_{\boldsymbol{k},\nu, 0} } + \frac{Z_{\nu}}{T_{\nu}} \left( \varGamma_{0, \boldsymbol{k}, \nu} - 1 \right) \tilde{\phi}_{\boldsymbol{k},0} + \frac{1}{\tilde{B}_0} \varGamma_{1, \boldsymbol{k}, \nu} \delta \tilde{B}_{{\parallel}, \boldsymbol{k}, 0} \right\} , \end{align}

(4.3)\begin{align} \hat{M}_{\boldsymbol{k}} & =- \frac{\beta_r}{\left(k_{{\perp}}\rho_r\right)^2} \sum_{\nu} Z_{\nu} \tilde{n}_{\nu} \tilde{v}_{\text{th},\nu} \frac{2\tilde{B}}{\sqrt{\rm \pi}} \int \,\mathrm{d}^2 \tilde{v} \tilde{v}_{{\parallel}} {\rm J}_{0,\boldsymbol{k}, \nu} \tilde{g}_{\boldsymbol{k},\nu, 0} \nonumber\\ & \quad + \left[1+ \frac{\beta_r}{ \left(k_{{\perp}}\rho_r\right)^2} \sum_{\nu} \frac{Z_{\nu} \tilde{n}_{\nu}}{\tilde{m}_{\nu}} \varGamma_{0,\boldsymbol{k}, \nu} \right] A_{{\parallel}, \boldsymbol{k},0 } , \end{align}

(4.4)\begin{align} \hat{N}_{\boldsymbol{k}} & = 2 \beta_r \sum_{\nu} \tilde{n}_{\nu} \tilde{T}_{\nu} \frac{ 2 \tilde{B}_0}{\sqrt{\rm \pi}} \int \,\mathrm{d}^2 \tilde{v} \tilde{\mu}_{\nu} \frac{{\rm J}_{0,\boldsymbol{k}, \nu}}{\tilde{a}_{\boldsymbol{k}, \nu}} \tilde{g}_{\boldsymbol{k},\nu, 0} + \left[\frac{\beta_r}{2\tilde{B}_0} \sum_{\nu} Z_{\nu}\tilde{n}_{\nu} \varGamma_{1,\boldsymbol{k}, \nu} \right] \tilde{\phi}_{\boldsymbol{k}, 0} \nonumber\\ & \quad + \left[1+ \frac{\beta_r}{2\tilde{B}_0} \sum_{\nu} Z_{\nu}\tilde{n}_{\nu} \tilde{T}_{\nu} \varGamma_{2,\boldsymbol{k}, \nu} \right] \delta \tilde{B}_{{\parallel}, \boldsymbol{k},0} , \end{align}

with the reference plasma beta, $\beta _r = 8{\rm \pi} n_r T_r/B_r^2$.

4.1.2 Adjoint normalisation

We choose the normalisation of the adjoint variables in such a way that our optimisation Lagrangian is dimensionless. In general, this is achieved by enforcing that the dimension of the adjoint variables in § 3 satisfy

(4.5)\begin{equation} \left. \begin{array}{ll@{}} {}[\lambda_{\nu}] = [\hat{G}_{\nu}]^{-1} & [\xi] = [\hat{Q}]^{-1}, \\ {}[\zeta] = [\hat{M}]^{-1} & [\sigma ] = [\hat{N}]^{-1}, \end{array} \right\} \end{equation}

with $[A]$ denoting the dimensionality of $A$. The normalised adjoint variables should then satisfy

(4.6)\begin{equation} \left. \begin{array}{ll@{}} \tilde{\lambda}_{\nu} = \dfrac{\lambda_{\nu}} {[\lambda_{\nu}] } & \tilde{\xi} = \dfrac{\xi}{[\xi] },\\ \tilde{\zeta} = \dfrac{\zeta}{[\zeta]} & \tilde{\sigma} = \dfrac{\sigma}{ [\sigma ] } . \end{array} \right\} \end{equation}

Analysis of (2.27), and (2.28)–(2.30) gives the normalisations consistent with those chosen for the gyrokinetic variables

(4.7)\begin{equation} \left. \begin{array}{ll@{}} \tilde{\lambda}_{\nu} = \lambda_{\nu} \dfrac{F_{0,\nu}} {{\rm e}^{-\tilde{v}_{\nu}^2} } \dfrac{\rho_r v_{{\rm th},r} }{a^2} & \tilde{\xi} = \xi \dfrac{n_r e \rho_r}{ a }, \\ \tilde{\zeta} = \zeta \dfrac{ B_r \rho_r^2}{a} & \tilde{\sigma} = \sigma \dfrac{B_r \rho_r }{a} . \end{array} \right\} \end{equation}

We then multiply the adjoint equation (3.13), written in terms of normalised coordinates, by a factor of $( \rho _r/a ) F_{0,\nu }/ \textrm {e}^{-\tilde {v}_{\nu }^2}$ to obtain the electromagnetic, collisional adjoint equations in normalised units

(4.8)\begin{gather} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}} + \tilde{\gamma}^* \tilde{\lambda}^{\leftrightarrow}_{\nu} + \tilde{v}_{{\rm th},\nu} \tilde{v}_{{\parallel}} \, \hat{\boldsymbol{b}}\boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{z}}} - \tilde{v}_{{\rm th},\nu} \, \tilde{\mu}_{\nu} \, \hat{\boldsymbol{b}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \frac{\partial{{\tilde{B}_0}}}{\partial{{\tilde{z}}}} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{v}_{{\parallel}}}} - {\rm i} \tilde{\omega}_{d, \nu} \tilde{\lambda}^{\leftrightarrow}_{\nu} \nonumber\\ + Z_{\nu} \tilde{n}_{\nu} {\rm J}_{0,\nu} \tilde{\xi} - \frac{\beta_r}{(k_{{\perp}} \rho_r)^2} Z_{\nu} \tilde{n}_{\nu} \tilde{v}_{\text{th},\nu}{\rm J}_{0,\nu} \tilde{v}_{{\parallel}} \tilde{\zeta} + 2\beta_r \tilde{T}_{\nu} \tilde{\mu}_{\nu} \frac{{\rm J}_{1,\nu}}{\tilde{a}_{\nu}} \tilde{\sigma} - \hat{C}_{\nu}[\lambda_{\nu}] = 0 , \end{gather}

(4.9)\begin{gather}\tilde{\xi} + \frac{1}{\tilde{\eta}} \sum_{\nu} \frac{2\tilde{B}}{\sqrt{\rm \pi}} \int \,\mathrm{d}^2 \tilde{v} \left[{\rm i} \tilde{\omega}_{*,\nu} + \frac{Z_{\nu}}{\tilde{T}_{\nu}} \tilde{\gamma}^* \right] {\rm J}_{0,\nu} {\rm e}^{-\tilde{v}_{\nu}^2} \tilde{\lambda}^{\leftrightarrow}_{\nu} =0 , \end{gather}

(4.10)\begin{gather}\tilde{\zeta} - \sum_{\nu} \frac{2\tilde{B}}{\sqrt{\rm \pi}} \int \,\mathrm{d}^2 \tilde{v} (2 \tilde{v}_{\text{th},\nu} \tilde{v}_{{\parallel}}) \left[ {\rm i} \tilde{\omega}_{*,\nu} + \frac{Z_{\nu}}{\tilde{T}_{\nu}} \tilde{\gamma}^* \right] {\rm J}_{0,\nu} {\rm e}^{-\tilde{v}_{\nu}^2} \tilde{\lambda}^{\leftrightarrow}_{\nu} = 0 , \end{gather}

(4.11)\begin{gather}\tilde{\sigma} - \sum_{\nu} \frac{2\tilde{B}}{\sqrt{\rm \pi}} \int \,\mathrm{d}^2 \tilde{v} \left( 4\tilde{\mu}_{\nu} \frac{\tilde{T}_{\nu}}{Z_{\nu}} \frac{{\rm J}_{1,\nu}}{\tilde{a}_{\nu}} \right) \left[ {\rm i} \tilde{\omega}_{*,\nu} + \frac{Z_{\nu}}{\tilde{T}_{\nu}} \tilde{\gamma}^* \right] {\rm e}^{-\tilde{v}_{\nu}^2} \tilde{\lambda}^{\leftrightarrow}_{\nu} = 0 , \end{gather}

with $\tilde {\lambda }^{\leftrightarrow }_{\nu } = \tilde {\lambda }_{\nu }(\tilde {z},-\tilde {v}_{\parallel },\tilde {\mu }_{\nu })$, and $\tilde {\eta } = \sum _{\nu } Z_{\nu }^2 \tilde {n}_{\nu }/ \tilde {T}_{\nu }$. This change of variables has been made such that the adjoint equations more closely resemble the gyrokinetic equations, with the boundary conditions on $\tilde {\lambda }^{\leftrightarrow }_{\nu }$ mimicking those on $\tilde {g}_{\boldsymbol {k},\nu, 0}$. This is convenient as the numerical machinery from existing gyrokinetic codes can be re-used. An artificial time dependence has also been introduced in (4.8) in order to make computation easier, and we solve for the steady-state solution for $\tilde {\lambda }_{\nu }$, when $\partial _{\tilde {t}} \tilde {\lambda }_{\nu } = 0$.Footnote ¹¹

As before, we can use ${\hat {G}_{\nu } [\boldsymbol {p}; \tilde {g}_{\nu }, \tilde {\phi }, \tilde {A}_{\parallel }, \delta \tilde {B}_{\parallel }]= \tilde {\gamma } \tilde {g}_{\nu } + \hat {L}_{\nu } [\boldsymbol {p}; \tilde {g}_{\nu }, \tilde {\phi }, \tilde {A}_{\parallel }, \delta \tilde {B}_{\parallel }]}$, with the added constraint of ${ \boldsymbol {\nabla }_{\boldsymbol {p}} \mathcal {L} |_{\boldsymbol {p}_0} = 0}$ to rearrange the above and obtain

(4.12)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{p}} \tilde{\gamma} \langle \tilde{g}_{\nu}, \tilde{\lambda}_{\nu} \rangle_{\tilde{z},\tilde{v}, {\nu}} =- \left. \left[ \langle \partial_{\boldsymbol{p}} \hat{L}_{\nu}, \tilde{\lambda}_{\nu} \rangle_{\tilde{z},\tilde{v}, {\nu}} + \langle \partial_{\boldsymbol{p}} \hat{Q}, \tilde{\xi} \rangle_{\tilde{z}} + \langle \partial_{\boldsymbol{p}} \hat{M} , \tilde{\zeta} \rangle_{\tilde{z}} + \langle \partial_{\boldsymbol{p}} \hat{N} , \tilde{\sigma} \rangle_{\tilde{z}} \right] \right|_{\boldsymbol{p}_0} , \end{equation}

with closure equations provided by (4.8)–(4.11).

4.1.3 Electrostatic collisionless limit

Finally, in the electrostatic, collisionless regime, the system of equations to solve in the normalised stella coordinates is given by

(4.13)\begin{gather} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}} + \tilde{\gamma}^* \tilde{\lambda}^{\leftrightarrow}_{\nu} + \tilde{v}_{{\rm th},\nu} \tilde{v}_{{\parallel}} \, \hat{\boldsymbol{b}}\boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{z}}} - \tilde{v}_{{\rm th},\nu} \, \tilde{\mu}_{\nu} \, \hat{\boldsymbol{b}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \frac{\partial{{\tilde{B}_0}}}{\partial{{\tilde{z}}}} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{v_{{\parallel}}}} - {\rm i} \tilde{\omega}_{d, \nu} \tilde{\lambda}^{\leftrightarrow}_{\nu} \nonumber\\ + Z_{\nu} \tilde{n}_{\nu} {\rm J}_{0,\nu} {\rm e}^{-\tilde{v}_{\nu}^2} \xi = 0 , \end{gather}

(4.14)\begin{gather}\xi + \frac{1}{\tilde{\eta}} \sum_{\nu} \frac{2\tilde{B}}{\sqrt{\rm \pi}} \int \,\mathrm{d}^2 \hat{v} \left[{\rm i} \tilde{\omega}_{*,\nu} + \frac{Z_{\nu}}{\tilde{T}_{\nu}} \tilde{\gamma}^* \right] {\rm J}_{0,\nu} {\rm e}^{-\tilde{v}_{\nu}^2} \tilde{\lambda}^{\leftrightarrow}_{\nu} = 0 , \end{gather}

and

(4.15)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{p}} \tilde{\gamma} \langle \tilde{g}_{\nu}, \tilde{\lambda}_{\nu} \rangle_{\tilde{z},\tilde{v}, {\nu}} =- \left. \left[ \langle \partial_{\boldsymbol{p}} \hat{L}_{\nu}, \tilde{\lambda}_{\nu} \rangle_{\tilde{z},\tilde{v}, {\nu}} + \langle \partial_{\boldsymbol{p}} \hat{Q}, \tilde{\xi} \rangle_{\tilde{z}} \right] \right|_{\boldsymbol{p}_0} . \end{equation}

4.2 Magnetic geometry

The coefficients in equations (4.1)–(4.4), and thus the associated linear growth rates, are implicitly dependent on the magnetic geometry. For the remainder of the paper we will take $\boldsymbol {p}$ to be an appropriate set of parameters that specifies the local magnetic geometry. In particular, we will use the Miller formalism (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) to parameterise the magnetic field on the flux surface of interest. The Miller approach ensures that the Grad–Shafranov (Shafranov Reference Shafranov1966) equation is always satisfied locally by using a set of independent parameters to define a single flux surface in an axisymmetric device. The model equations describing the shape of the flux surface with flux label $r$ are

(4.16)\begin{gather} R(r,\theta) = R_0 (r) + r \cos \left\{ \theta + \sin(\theta) \delta(r) \right\}, \end{gather}

(4.17)\begin{gather}Z(r,\theta) = r \kappa (r) \sin (\theta). \end{gather}

Here, $R(r,\theta )$ and $Z(r,\theta )$ define the major radial and vertical locations for a given poloidal location, $\theta$, which is related to the cylindrical angle, and $\delta$ and $\kappa$ indicate the triangularity and elongation of the surface, respectively. Specifying the full set of Miller parameters provides all of the information required to compute the geometric coefficients in the gyrokinetic-adjoint system consistent with the local magnetohydrodynamic (MHD) equilibrium (though consistency with a global MHD equilibrium is not guaranteed).

The user specified input parameters used in the local version of stella are $\{ r_{\psi _0}, R_{\psi _0}, \varDelta _{\psi _0}, q_{\psi _0}, \hat {s}_{\psi _0}, \kappa _{\psi _0}, \kappa '_{\psi _0}, I_{\psi _0}, \delta _{\psi _0}, \delta '_{\psi _0}, \beta '_{\psi _0} \}$, corresponding to markers for the minor and major radii, horizontal Shafranov shift ($\varDelta _{\psi _0} = R_{\psi _0}'$), safety factor, magnetic shear ($\hat {s} \doteq (r/q) q'$), elongation and its radial derivative, the external axial discharge current (which is used as a proxy for the reference magnetic field), triangularity and its radial derivative and the radial pressure derivative, respectively, with each being specified at the flux surface of interest. Further information about how the Miller geometry is treated in stella is given in Appendix E. We henceforth set $\boldsymbol {p} := \{ r_{\psi _0}, R_{\psi _0}, \varDelta _{\psi _0}, q_{\psi _0}, \hat {s}_{\psi _0}, \kappa _{\psi _0}, \kappa '_{\psi _0}, I_{\psi _0}, \delta _{\psi _0}, \delta '_{\psi _0}, \beta '_{\psi _0} \}$.

5.1 Initial simulation

Solving, at an initial set of $\boldsymbol {p}_0$, for $\tilde {\gamma }, \tilde {g}_{\nu }$ and $\tilde {\phi }$ is achieved by allowing stella to run for a sufficiently long time, such that the solution is dominated by a single normal mode. To determine the time at which this is satisfied, we employ a convergence test: the growth rate is calculated at each time step and if the value of this is constant in time (within a specified tolerance) then the system is deemed to be converged. There are two components of the convergence test. The first is to check that the growth rates calculated at adjacent time steps are within a given tolerance of each other. The second is to perform a windowed average to check that the growth rate remains consistent over a defined number of time steps. Two windowed averages are done; one over $N_t$ time steps and one over $\mathbb {Z}(N_t/2)$ time steps. When these two window-averaged growth rates agree within a set tolerance, $\tilde {g}_{\nu }$ and $\tilde {\phi }$ are taken to have converged. The corresponding growth rate is then calculated from the windowed average.Footnote ¹²

5.2 Adjoint simulation

As previously introduced, an artificial time dependence is added to the adjoint equations to facilitate computation. The solution is found in the steady-state limit, in which the time derivative appearing in the adjoint equations goes to zero.

The resulting adjoint equations are treated in a similar way to the treatment of the usual gyrokinetic system of equations in stella. The main aim is to ensure that the parallel streaming term may be treated separately from the rest of the dynamics through operator splitting. This is done by discretising in time, and splitting the time derivative into a series of three steps

(5.1)\begin{equation} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}} = \left(\frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}}\right)_1 + \left(\frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}}\right)_2 + \left(\frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}} \right)_3 , \end{equation}

where

(5.2a)\begin{gather} \left(\frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}}\right)_1 =- \tilde{\gamma}^* \tilde{\lambda}^{\leftrightarrow}_{\nu} + {\rm i} \tilde{\omega}_{d, \nu} \tilde{\lambda}^{\leftrightarrow}_{\nu} - Z_{\nu} \tilde{n}_{\nu} \tilde{\xi} , \end{gather}

(5.2b)\begin{gather}\left(\frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}}\right)_2 = \tilde{v}_{{\rm th},\nu} \, \tilde{\mu}_{\nu} \, \hat{\boldsymbol{b}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \frac{\partial{{\tilde{B}_0}}}{\partial{{\tilde{z}}}} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{v}_{{\parallel}}}}, \end{gather}

(5.2c)\begin{gather}\left(\frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{t}}} \right)_3 =- \tilde{v}_{{\rm th},\nu} \tilde{v}_{{\parallel}} \, \hat{\boldsymbol{b}}\boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{z} \frac{\partial{\tilde{\lambda}^{\leftrightarrow}_{\nu}}}{\partial{\tilde{z}}} . \end{gather}

Analogous to stella, the terms in (5.2a) are treated explicitly using a strong-stability-preserving, third-order Runge–Kutta method (Gottlieb, Shu & Tadmor Reference Gottlieb, Shu and Tadmor2001). In combination with the operator splitting employed within stella the overall algorithm is second-order accurate in time step size, $\varDelta \tilde {t}$.

The parallel streaming and mirror terms, given by (5.2c) and (5.2b), respectively, are treated separately, due to the presence of the prefactor $\tilde {v}_{\textrm {th},\nu }$, which increases the relative amplitude of these terms when considering electron dynamics. As a result these terms have the potential to exert a stringent Courant–Friedrichs–Lewy condition on the simulation, and require a small time step to be taken in order to retain accuracy. Thus these terms are treated implicitly in time to relax this constraint.

The mirror term, (5.2b), is a simple advection equation of $\tilde {\lambda }^{\leftrightarrow }_{\nu }$ in $\tilde {v}_{\parallel }$, which is treated using a semi-Lagrange method, akin to the algorithm used to advect the distribution function in $\tilde {v}_{\parallel }$ within stella. The mirror coefficients are independent of both time and $\tilde {v}_{\parallel }$ and hence the exact characteristics of this equation are known. The interpolation in $\tilde {v}_{\parallel }$ is forth order accurate in $\tilde {v}_{\parallel }$ grid spacing, $\varDelta \tilde {v}_{\parallel }$(Barnes et al. Reference Barnes, Parra and Landreman2019).

The streaming term, (5.2c), is also an advection equation in $\tilde {z}$, which is treated using the Thompson algorithm for tri-diagonal solve. When fully centred in time the discretisation reduces to the Crank–Nicolson method, which is second-order accurate in time, and $\tilde {z}$ grid spacings, $\varDelta \tilde {t}$ and $\varDelta \tilde {z}$, respectively.

We are seeking a steady-state solution to the adjoint equations. The same convergence test is performed on $\tilde {\lambda }^{\leftrightarrow }_{\nu }$ as that performed on the distribution function in order to check that the complex growth rate has converged to zero within a given tolerance. When this is satisfied, the resulting $\tilde {\lambda }^{\leftrightarrow }_{\nu }$ is used to solve for $\tilde {\lambda }_{\nu }$. This is then stored for use in the remainder of the calculation.

5.3 Optimisation loop

Once the gradients $\boldsymbol {\nabla }_{\boldsymbol {p}} \gamma$ are obtained we use them inside an optimisation loop to find the $\boldsymbol {p}$ that minimises $\gamma$. We employ the Levenberg–Marquardt (LM) algorithm, (Transtrum & Sethna Reference Transtrum and Sethna2012), to find the local minimum. The method adopts a steepest descent behaviour when the location in parameter space is considered to be far from the minimum, and progresses towards Gauss–Newton behaviour as one approaches the minimum. This is achieved by introducing a damping factor, $\varGamma$, which is updated with each iteration. The algorithm iteratively solves the following:

(5.3)\begin{equation} [\boldsymbol{H} + \varGamma \text{diag}(\boldsymbol{H}) ] \,\mathrm{d} \boldsymbol{p} =- \boldsymbol{\nabla}_{\boldsymbol{p}} \gamma , \end{equation}

with $\boldsymbol {H} = \boldsymbol {\nabla }_{\boldsymbol {p}}^2 \gamma \approx (\boldsymbol {\nabla }_{\boldsymbol {p}} \gamma )^{{\dagger} } \boldsymbol {\nabla }_{\boldsymbol {p}} \gamma$ the Hessian matrix. When $\varGamma$ is large we have $\boldsymbol {p}_{\textrm {new}} \approx \boldsymbol {p}_{\textrm {old}} - \boldsymbol {\alpha }\boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {p}} \gamma$, with $\boldsymbol {\alpha } \doteq [\varGamma \textrm {diag} (\boldsymbol {H})]^{-1}$, which mimics the gradient-descent algorithm. However, when $\varGamma$ is small, (5.3) reduces to $\boldsymbol {p}_{\textrm {new}} \approx \boldsymbol {p}_{\textrm {old}} - \boldsymbol {H}^{{\dagger} } \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {p}} \gamma$, which matches with the Gauss–Newton algorithm.

The LM formalism is derived using the Taylor expansion, and as such a trust region is included within the optimisation loop to ensure that the updated value of $\boldsymbol {p}$ is close enough to the previous, such that the Taylor approximation is valid within the limits for which the algorithm is applied. The trust region for $\boldsymbol {p}$ is defined via

(5.4)\begin{equation} \rho = \frac{0.5 \,\mathrm{d} \boldsymbol{p}^{{\dagger}} \boldsymbol{\cdot} \boldsymbol{H} \boldsymbol{\cdot} \mathrm{d} \boldsymbol{p}}{\mathrm{d} \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{p}} \gamma} < \bar{\epsilon} , \end{equation}

where $\bar {\epsilon }$ is a chosen tolerance. If $\rho > \bar {\epsilon }$ for a given $\mathrm {d} \boldsymbol {p}$ the algorithm rejects the output $p_{\textrm {new}}$ and increases the weight $\varGamma$ in an attempt to improve the accuracy of the approximation. This helps ensure that the updated value of $\boldsymbol {p}$ is a reasonable one.

It is worth noting that the LM algorithm is designed to find local minima, so there is no guarantee that the minimum obtained is the global minimum of the system.

We also emphasise that the gradient-based optimisation algorithm is independent of the adjoint method that has been developed for gyrokinetic microstability. The optimisation loop may be itself optimised to efficiently search for regions of stability given a gradient input. Different algorithms, and indeed different parameter choices within each algorithm, will yield different efficiencies in finding stable solutions. An illustrative example of this is given later in figure 1, where the step size for the optimisation loop is varied to yield two distinct paths through the parameter space using the adjoint gradient. However, this is not the focus of this paper and we will not labour on enhancing the gradient-based optimisation loop.

Figure 1.

Two-dimensional parameter scan over elongation and triangularity, with the colour map indicating the amplitude of the linear growth rate. Here, $k_y = 0.68$, $k_x = 0.0$, $m_{e}/m_{i} = 2.7 \times 10^{-4}$, $T_i = T_e = 1$, $n_i = n_e = 1$, $a/L_{T_i} = a/L_{T_e} = 2.42$, $a/L_{n_i} = a/L_{n_e} = 0.81$, with $a$ the minor radius of the last closed flux surface. The path taken by the optimisation algorithm is indicated in white, with the initial point $\kappa = 1.5$, and $\delta = 0.14$. A second path, drawn in grey, is shown indicating the adjoint optimisation with a different step size for the optimisation loop.

6.1 Initial benchmark

The first numerical check we perform is to ensure that the values of $\mathrm {d}_{\boldsymbol {p}} \gamma _0$ obtained from our adjoint method agree with those obtained using a finite-difference approach. Following this, we perform a more extensive benchmark by conducting a parameter scan in the growth rate using stella for different values of the Miller parameters. We then choose an initial set of parameters and perform the adjoint-optimisation scheme described above, and overlay the results of this with the parameter scan to see how these compare, whilst checking the growth rates at each point considered by the adjoint-optimisation scheme against those obtained using a finite-difference approach. As a proof of principle we chose to vary two parameters: triangularity, $\delta$, and elongation, $\kappa$, whilst holding the other Miller variables fixed. Given that the Miller parametrisation is local to a given flux surface, this variation is not necessarily consistent with a global solution to the Grad–Shafranov equation. The choice to vary these two parameters in isolation is driven by two considerations; first, it is only intended as a proof-of-principle check of the adjoint approach so simplicity is desirable, and second, previous research has shown that maximal shaping, with large elongation and triangularity, minimises the linear ITG instability, whereas parameters such as $\kappa '$ and $\delta '$ have an order $a/R_0 \ll 1$ is the inverse aspect ratio, with $a$ and $R_0$ the scales associated with the minor and major radii, respectively (Highcock et al. Reference Highcock, Mandell, Barnes and Dorland2018). Table 2 lists the values of input equilibrium variables in the Miller geometry. These have been chosen to coincide with values used in Beeke et al. (Reference Beeke, Barnes, Romanelli, Nakata and Yoshida2020) in order to verify the qualitative behaviour found.

Table 2.

List of Miller Parameters.

Equilibrium Miller parameter values used in the initial benchmark simulations.

Given these initial values of $\delta$ and $\kappa$, we determine the linear growth rates for a grid of perpendicular wavenumbers within $k_x<2.0$, $k_y<2.0$, which reveals that the most unstable mode is found at $\{k_y, k_x\} = \{0.68, 0.0\}$ for mass ratio $m_{e}/m_{i} = 2.7 \times 10^{-4}$, normalised species temperatures and densities of $T_i = T_e = 1$, $n_i = n_e = 1$ and normalised species temperature and density gradients of $a/L_{T_i} = a/L_{T_e} = 2.42$, $a/L_{n_i} = a/L_{n_e} = 0.81$.

In figure 1 we show a scan in the linear growth in elongation and triangularity obtained by running stella with the Miller parameters specified in table 2, and the values of $\kappa$ and $\delta$ adjusted accordingly. The contour colour indicates the magnitude of the growth rate, and the plot extends over a range of values that has been set by reasonable physical constraints on devices.

Figure 1 shows that increasing the elongation and triangularity of our flux surface reduces the linear growth rate, and that there exists a region of stability when the shaping is maximal, in agreement with previous work. The path taken using the adjoint method is indicated in white. At the chosen starting point, located in the region of instability, the gradient $\mathrm {d}_{\boldsymbol {p}}\gamma$ is calculated and the value of $\boldsymbol {p} = \{\delta, \kappa \}$ is updated using the previously mentioned LM method. The final point is found to be locally stable as the growth rate here is negative. The algorithm then checks a nearby point to determine if a small region of stability exists and, once this has been verified, outputs this as the final $\boldsymbol {p}$ value. In this particular case of a two-dimensional parameter scan, a finite-difference approach could utilise the same optimisation loop, requiring only three simulations, plus one additional simulation performed at the next iteration point to verify the growth rate. The adjoint method must simulate the gyrokinetic equation once, and also solve the adjoint equations – which are of similar computational cost as the gyrokinetic equation – in order to achieve the same result. In this particular low-dimensionality case, the advantage of the adjoint approach is comparatively small over using a finite-difference approach. However, the computation time required for the traditional method scales linearly with the number of parameters, so that for an $N$-dimensional parameter space, $N+1$ simulations are necessary at each iteration point, whilst the adjoint method is effectively independent of the number of parameters, and only the steady-state solutions to the gyrokinetic and adjoint equations are required. It can therefore be readily applied to a high-dimensional parameter space without any significant increase in computational cost.

A second path is plotted on figure 1 in grey. This path is taken using the same adjoint technique, but increasing the step size within the optimisation loop. When the step size is small, as with the white path, the LM algorithm more closely resembles a gradient-descent method, however, when the step size is larger, as with the grey path, the LM algorithm resembles Newton's method for gradient optimisation. The figure illustrates how the adjoint algorithm can be combined with an optimiser to quickly and efficiently converge to a stable region of parameter space.

6.2 Increasing the critical temperature gradient

We have seen that the adjoint method is a powerful technique for computing stable points in a large parameter space; figure 1 shows that it can efficiently be used inside an optimisation loop to find a minimum of the linear growth rate for a given temperature gradient. Once the adjoint-optimisation loop locates a region in the parameter space with negative or zero growth rate the local temperature gradient (or other plasma profile variable of interest) may be increased and the process repeated. Conversely, if a minimum positive (unstable) growth rate is found, the temperature gradient can be reduced to seek out the optimal shape that maximises the critical temperature gradient at which linear instability occurs to seek out microstability.

Figure 2 demonstrates an iterative use of the adjoint optimisation. Here, we have iterated the temperature gradient, and inside each temperature-gradient iteration the adjoint method is used to find a stable geometric configuration. However, this principle could be extended further by continuing to increase the temperature gradient until no region of stability is available, indicating a limiting temperature gradient that can be achieved through geometric considerations alone. Figure 2 is a demonstration that the adjoint method may be used to increase the temperature gradient, whilst retaining stability using geometry. Though we have opted to consider only two parameters in the above as a demonstration and for clarity have focused on ITG, it is possible to employ the adjoint method to optimise over a large number of geometric parameters simultaneously. Such an exploration would be expensive using traditional finite-difference methods.

Figure 2.

Plots showing a parameter scan in elongation and triangularity, with a temperature gradient of $a/L_{T_i} = 3.80$, increased from the $a/L_{T_i} = 2.42$ value used in figure 1. The geometry of the initial point, located in the unstable region, is provided by final point in figure 1 and is now unstable due to the increased temperature gradient. Here, $k_y = 0.68$, $k_x = 0.0$, $m_{e}/m_{i} = 2.7 \times 10^{-4}$, $T_i = T_e = 1$, $n_i = n_e = 1$, $a/L_{T_i} = a/L_{T_e} = 3.80$, $a/L_{n_i} = a/L_{n_e} = 0.81$. Note that the colour scales used in the figures above are different than that used in figure 1. The right-hand side plot is a zoomed in figure of the left.

6.3 Negative triangularity

To form a final example here, we note there has been previous evidence that negative triangularity can offer improvements for microstability, Pochelon et al. (Reference Pochelon, Goodman, Henderson, Angioni, Behn, Coda, Hofmann, Hogge, Kirneva and Martynov1999), so we repeated the benchmark scan using negative triangularity. All values of equilibrium Miller parameters are the same as in table 2, except we now set the initial value of triangularity to $\delta = -0.14$.

The path taken using the adjoint-optimisation loop is again shown in figure 3 by the white path, starting in the dark red region at $\{\delta, \kappa \} = \{-0.14, 1.52\}$.

Figure 3.

Growth rate contours for a parameter scan with negative triangularity for $k_y = 0.68$, $k_x = 0.0$ and equilibrium parameters $m_{e}/m_{i} = 2.7 \times 10^{-4}$, $T_i = T_e = 1$, $n_i = n_e = 1$, $a/L_{T_i} = a/L_{T_e} = 3.80$, $a/L_{n_i} = a/L_{n_e} = 0.81$. The white line indicates the path taken by the optimisation algorithm. The initial values of $\{\delta, \kappa \}$ are taken to be $\{-0.14, 1.52\}$.

This highlights a key feature of the gradient-based optimisation method: the solution is not unique, and the output can depend on the starting region within parameter space. See figure 4 for an illustrative example of this.

Figure 4.

Growth rate contours for a parameter scan with both positive and negative triangularity for $k_y = 0.68$, $k_x = 0.0$ and equilibrium parameters $m_{e}/m_{i} = 2.7 \times 10^{-4}$, $T_i = T_e = 1$, $n_i = n_e = 1$, $a/L_{T_i} = a/L_{T_e} = 3.80$, $a/L_{n_i} = a/L_{n_e} = 0.81$. The white line indicates the two paths taken by the optimisation algorithm starting in different regions in parameter space. The initial values of $\{\delta, \kappa \}$ are taken to be $\{0.14, 1.52\}$ and $\{-0.14, 1.52\}$.

Finally, we consider increasing the temperature gradient for the negative triangularity case, and we repeat the procedure to find a stable region of parameter space. We again take the input parameters $\{\delta, \kappa \} = \{-0.9965, 2.4488\}$ to be the outputs from the previously optimised case at a temperature gradient of $a/L_{T_i} = 2.42$. Then we use the adjoint algorithm coupled with the gradient optimiser to look for a stable region of parameter space at an increased temperature gradient of $a/L_{T_i} = 3.8$. The results of this are shown in figure 5.

Figure 5.

Growth rate contours for a parameter scan with negative triangularity at a temperature gradient of $a/L_{T_i}= 3.80$, increased from $a/L_{T_i}= 2.42$. The geometry of the initial point is taken as the final point in figure 3, and is now unstable due to the increased temperature gradient. The scan is performed at the same parameter values as those in figure 3 – $k_y = 0.68$, $k_x = 0.0$ $m_{e}/m_{i} = 2.7 \times 10^{-4}$, $T_i = T_e = 1$, $n_i = n_e = 1$, $a/L_{T_i} = a/L_{T_e} = 3.80$, $a/L_{n_i} = a/L_{n_e} = 0.81$. The white line indicates the path taken by the optimisation algorithm. The initial values of $\{\delta, \kappa \}$ are taken to be $\{-0.9965, 2.4488\}$.

7.1 Numerical demonstration of improved efficiency

For the case where $\mathcal {N} =2$ we are optimising with respect to two parameters, $\{ \delta, \kappa \}$. We find that to calculate a gradient at each point in our parameter space requires a total of $\sim 4.752$ CPU (central processing unit) hours on 4 nodes, each with 48 cores on the supercomputer Marconi. The same calculation as performed using the adjoint method requires a total of ${\sim }3.168$ CPU hours. Though this is only a modest improvement, we shall show that as $\mathcal {N}$ increases the advantage of using the adjoint method, over a finite-difference scheme, becomes increasingly apparent.

Increasing the number of parameters to $\mathcal {N} = 4$, optimising over the miller parameters $\{ \varDelta, \kappa, q,\delta \}$, we find that the finite-difference approach requires ${\sim } 7.920$ CPU hours to compute a gradient at each point in parameter space. However, when computing the same gradient using the adjoint method the CPU time, to the precision of the CPU clock, is ${\sim }3.168$ CPU hours.

If we increase the number of parameters further to $\mathcal {N} = 7$, optimising over the miller parameters $\{R_0, R_{\mathrm {geo}}, \varDelta, \kappa, q, \hat {s}, \delta \}$ we find that using the finite-difference approach requires $\sim 12.672$ CPU hours to compute each gradient at a point in parameter space. When using the adjoint method to compute the same gradient, to the precision of the CPU clock, is still $\sim 3.168$ CPU hours.

Hence, we conclude that the numerical speed up is significant with increasing $\mathcal {N}$. This allows for the potential of including multiple harmonics within our shaping optimisation for very little additional computing cost.

For the example including seven parameters, we iterate the process of stepping through parameter space to find a point of stability. We take our set of initial parameters to be those given in table 2, and we simulate using a temperature and density gradient of $a/L_{T_i} = a/L_{T_e} = 2.42$ and $a/L_{n_i} = a/L_{n_e} = 0.81$, respectively. We perturb the magnetic geometry by varying $\{R_0, R_{\mathrm {geo}}, \varDelta, \kappa, q, \hat {s}, \delta \}$. Following iterations of the coupled adjoint-LM system we find a configuration that is stable when $\{R_0, I, \varDelta, \kappa, q, \hat {s}, \delta \} = \{2.979, 2.846, 0.562, 1.656, 2.085, 0.167, 0.225\}$. A cross-section of the initial, unstable configuration and the final, stable configuration, with their corresponding neighbouring flux surfaces, are shown in figure 6.

Figure 6.

Plots of the flux surfaces in the poloidal cross-section. The orange surface is the flux surface at $\rho = 0.5$, and the blue and green surfaces are the two adjacent flux surfaces. The image (a) is the initial unstable configuration, before optimising. The image (b) is the stable, optimised configuration found using the adjoint-LM system.