1. Introduction
Turbulence plays an important role in stellarators, particularly if the neoclassical losses have been reduced by optimization of the magnetic field (Beidler et al. Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021). Confinement has been studied in Wendelstein 7-X (W7-X) experiments (Klinger et al. Reference Klinger, Alonso, Bozhenkov, Burhenn, Dinklage, Fuchert, Geiger, Grulke, Langenberg and Hirsch2016; Dinklage et al. Reference Dinklage, Beidler, Helander, Fuchert, Maaßberg, Rahbarnia, Sunn Pedersen, Turkin, Wolf and Alonso2018; Wolf et al. Reference Wolf, Alonso, Äkäslompolo, Baldzuhn, Beurskens, Beidler, Biedermann, Bosch, Bozhenkov and Brakel2019; Carralero et al. Reference Carralero, Estrada, Maragkoudakis, Windisch, Alonso, Beurskens, Bozhenkov, Calvo, Damm and Ford2021), and these have been accompanied by nonlinear gyrokinetic simulations, but most of the latter have been performed in a flux tube (Alcusón et al. Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020; Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). However, in contrast to an ideally axisymmetric tokamak, stellarator flux tubes are not equivalent to each other and there is no good way to incorporate the ambient radial electric field. Global electrostatic turbulence has been addressed in stellarator plasmas by Navarro et al. (Reference Navarro, Merlo, Plunk, Xanthopoulos, von Stechow, Siena, Maurer, Hindenlang, Wilms and Jenko2020) and Sánchez et al. (Reference Sánchez, Mishchenko, García-Regaña, Kleiber, Bottino and Villard2020). Electromagnetic turbulence, involving fluctuations in both the electric and magnetic fields, arises at sufficiently high plasma pressure and is therefore expected to be important in well-confined fusion plasmas. For stellarators, numerical treatment of such turbulence in a global domain is still in its infancy (Mishchenko et al. Reference Mishchenko, Biancalani, Bottino, Hayward-Schneider, Lauber, Lanti, Villard, Kleiber, Könies and Borchardt2021; Wilms et al. Reference Wilms, Navarro, Merlo, Leppin, Görler, Dannert, Hindenlang and Jenko2021; Mishchenko et al. Reference Mishchenko, Bottino, Hayward-Schneider, Poli, Wang, Kleiber, Borchardt, Nührenberg, Biancalani and Könies2022), due in part to computational complexity and the technical challenge of including magnetic fluctuations. One problem requiring such a treatment is the transition between ion temperature gradient (ITG) and kinetic mode (KBM)-driven turbulence. It has been shown by Aleynikova et al. (Reference Aleynikova, Zocco, Xanthopoulos, Helander and Nührenberg2018) that the linear growth rate of the KBMs has its maximum at large-scale modes in the W7-X. Such modes cannot be properly treated by local gyrokinetic codes, resulting in a failure of the turbulence to saturate at finite amplitude (Pueschel et al. Reference Pueschel, Hatch, Görler, Nevins, Jenko, Terry and Told2013).
A fully electromagnetic extension of the global nonlinear gyrokinetic code EUTERPE (Kornilov et al. Reference Kornilov, Kleiber, Hatzky, Villard and Jost2004) has been developed, which includes all components of the electromagnetic field perturbations. In this paper, results have been obtained for the first time showing saturated KBM turbulence in stellarator plasmas, even for plasma equilibria that have large-scale magnetohydrodynamic instability (Roberg-Clark, Plunk & Xanthopoulos Reference Roberg-Clark, Plunk and Xanthopoulos2022). The parallel component of the magnetic field perturbation, known to be important for KBMs, especially for the large-scale modes (Aleynikova et al. Reference Aleynikova, Zocco, Xanthopoulos, Helander and Nührenberg2018), has been included in the long-wavelength approximation. The nonlinear saturation of the KBM turbulence via a combination of the zonal-flow excitation and profile relaxation has been observed with the heat flux pointing outwards and the particle flux inwards. The latter so-called ‘particle pinch’ may explain the cases where the hollow density profiles have not been experimentally observed (Wolf et al. Reference Wolf, Alonso, Äkäslompolo, Baldzuhn, Beurskens, Beidler, Biedermann, Bosch, Bozhenkov and Brakel2019) in the plasma core as expected from the neoclassical theory (Beidler et al. Reference Beidler, Feng, Geiger, Köchl, Maßberg, Marushchenko, Nührenberg, Smith and Turkin2018). The effect of the parallel perturbation of the magnetic field on the stellarator turbulence is addressed.
The paper is organized as follows. In § 2, we briefly describe the equations solved and the discretization employed. In § 3, simulation results are presented and discussed. In § 4, we draw some conclusions.
2. Equations and discretization
In this section, we present for completeness the equations solved by the code. We use the mixed-variable formulation of gyrokinetic theory. This formulation is based on the splitting of the magnetic potential into the two parts $A_{\|} = A_{\|}^{{(h)}} + A_{\|}^{{(s)}}$, see Mishchenko et al. (Reference Mishchenko, Koenies, Kleiber and Cole2014) for details on the derivation and notations.
The equations include the gyrokinetic Vlasov equation,
Here, $f_{1s}$ is the perturbed part of the distribution function with the species index $s=i,e$, and $F_{0s}$ is the ambient distribution function (the control variate), a Maxwellian in this paper. The gyrocentre phase-space coordinates are used with the gyrocentre ${\boldsymbol {R}}$, the parallel velocity $v_{\|}$, the magnetic moment $\mu = m_s v_{\perp }^2/(2 B)$. Here, $v_{\perp }$ denotes the perpendicular velocity and $m_s$ is particle's mass. The equations for the gyrocentre orbits are
Here, the superscript $(0)$ denotes unperturbed orbits and $(1)$ corresponds to the orbit perturbations by the self-consistent electromagnetic field which has three physical components: the electrostatic potential $\phi$; the parallel component of the magnetic vector potential $A_{\|} = A_{\|}^{{(h)}} + A_{\|}^{{(s)}}$; and the parallel component of the magnetic field perturbation $\delta B_{\|}$. The angular brackets $\langle \ldots \rangle$ denote the average over the gyrophase $\alpha$, $B$ is the absolute value of the ambient magnetic field ${\boldsymbol {B}}$, and ${\boldsymbol {b}} = {\boldsymbol {B}}/B$.
The gyrokinetic system of equations is closed with the field equations:
Here, $\beta _s = \mu _0 n_{0s} T_{0s}/B^2$, $\rho _s = \sqrt {m_s T_s}/(q_s B)$, $\omega _{{\rm ci}} = q_{i} B/m_{i}$ with $q_s$ the particle charge of the species $s=i,e$ and $m_s$ the particle mass, $\delta ^3_{{\rm gy}} = \delta ({\boldsymbol {R}} + \boldsymbol {\rho }-{\boldsymbol {x}})$, ${\boldsymbol {\rho }} = v_{\perp } \hat {\boldsymbol \zeta }/\omega _{{c}s}$, $\hat {\boldsymbol \zeta }$ is the unit vector following the gyromotion of a particle, ${\rm d}^6 Z = B_{\|}^* \, {\rm d} {\boldsymbol {R}} \, {\rm d} v_{\|} \, {\rm d}\mu \, {\rm d}\alpha$, and $\alpha$ is the gyrophase. Further, an equation for $\partial A_{\|}^{(s)} / \partial t$,
and the perpendicular pressure balance condition for $\delta B_{\|}$,
are needed. Note that the terms containing $\delta B_{\|}$ are included to the lowest order in $k_{\perp }\rho _{i}$, see appendix A for details. For the cancellation problem, the pullback mitigation (Mishchenko et al. Reference Mishchenko, Koenies, Kleiber and Cole2014) is applied.
3. Simulations
3.1. Electromagnetic turbulence in W7-X
We consider a high-mirror (KJM) configuration (Geiger et al. Reference Geiger, Beidler, Feng, Maaßberg, Marushchenko and Turkin2014; Aleynikova, Zocco & Geiger Reference Aleynikova, Zocco and Geiger2022) of the W7-X down-scaled to $L_x = 400$ where $L_x$ is the radial size of the simulation box measured in the sound gyroradii. The simulation box is chosen in such a way that the whole cross-section of the W7-X fits into it. The reduced mass ratio $m_{i}/m_{e} = 200$ is used throughout the paper. We choose a flat density and the temperature profile following Riemann, Kleiber & Borchardt (Reference Riemann, Kleiber and Borchardt2016),
Here, $s$ is the toroidal flux normalized to its value on the plasma edge. Throughout the paper, we use the following numerical parameters: $0.75\times 10^9$ of ion markers; $10^9$ of electron markers; 128 radial grid points; 512 grid points in the poloidal direction; and 128 points covering one of the five ambient magnetic field periods toroidally. Three-dimensional tensor products of second-order B splines are used as the finite elements, the diagonal poloidal Fourier filter following the rotation transform profile is applied with the filter width of $35$ modes centred at the poloidal mode number locally satisfying the condition $k_{\|} = 0$. In the toroidal direction, Fourier modes $-30 \le n \le 30$ are resolved in a single period of W7-X (which has five periods in total). The time step used is $\Delta t = 0.5 \omega _{ci}^{-1}$. The simulation domain is restricted to the flux surfaces $0.25 \le s \le 0.75$ with the Dirichlet boundary conditions applied on both boundaries. For the noise control, we use the Krook operator with $\gamma _{{\rm Krook}} = 10^{-4}\omega _{ci}$ for both species. No other sources are applied in the simulations described in this paper.
In figure 1(a), the value of the plasma density is scanned which is equivalent to a scan in the plasma beta. One can see a characteristic shape of the dependence of the growth rate on beta, indicating the electromagnetic stabilization of the ITG branch and a destabilization of the KBM branch with increasing beta. Here, beta is defined as $\beta = 4 \mu _{0} n_* T_{*} / B_*^{2}$ with
with the toroidal angle $\varphi$, the total number of the ions $N_{{\rm ph},i}$ and the plasma volume $V$. The electron heat flux at $\beta = 2.8\,\%$ (KBM regime) and at $\beta = 0.28\,\%$ (electromagnetic ITG regime) is shown in figure 1(b). One sees that the simulations clearly enter the nonlinear phase, and for the KBM case ($\beta =2.8\,\%$), a comparison of the simulation including $\delta B_{\|}$ with a simulation neglecting $\delta B_{\|}$ is also shown. One sees that the quantitative difference between these cases is around $10\,\%$.
In figure 2(a), one sees that the electrostatic potential evolves nonlinearly and saturates after the initial linear growth. Here, the mode evolution is shown for all toroidal harmonics included in the simulations at the flux surface $s=0.5$ and the poloidal angle $\theta = 0$. The radial dependency of these toroidal harmonics is plotted at one time step in figure 2(b). One sees that all the modes are well resolved, smooth and confined inside the radial simulation domain which indicates robustness of our numerics.
In figure 3, the time evolution of the electrostatic and electromagnetic energies is shown for different $\beta$ values. One sees how the magnetic component of the perturbed-field energy increases with $\beta$. The electrostatic and the electromagnetic contributions into the perturbed-field energy become comparable near the transition point $\beta = 1.68\,\%$ where the KBM branch becomes dominant as can be seen in figure 1(a).
In figure 4, the evolution of the electrostatic potential $\phi _n(s)$ (contour plot) and the second radial derivative of its zonal component $\partial ^2\phi _{n=0}/\partial s^2$ (white curve) are shown. The ambient temperature profile is plotted, too (magenta curve). The electrostatic potential is measured at the poloidal angle $\theta = 0$ and shown as a function of the flux surface coordinate and the toroidal mode number. The second radial derivative of the zonal potential is related to the shearing rate of the zonal flow. One can see how the simulation starts with the linear instability with a toroidal mode number $n \approx 15$ developing around the position of the maximal ambient temperature gradient. Only a single period of the five-period W7-X has been simulated here. It implies that for the whole torus, the toroidal mode number must be multiplied with the total number of the periods $N_{\rm per} = 5$ resulting in the whole-torus toroidal mode number $\bar {n}_{\rm tot} = 75$. The linear instability evolves, spreading radially and generating low-mode-number components including the zonal flow at the later times. The second radial derivative of the zonal potential, initially located at the position of the maximum temperature gradient, also spreads radially following the turbulence evolution and the resulting zonal-flow generation.
In figure 5(a), the total temperature radial profile is plotted for both species together with the initial temperature. One sees that the temperature profile evolves in response to the turbulence. The underlying colour plot shows the turbulent electrostatic potential $\phi _n(s)$ as a function of the radius and the toroidal mode number at a given time chosen to be in the saturation phase, see figure 2(a). The same turbulent electrostatic potential is shown as a background colour plot also in figures 5(b), 6(a) and 6(b) where the colour map is switched to black–white to simplify the graph. The energy flux plotted for both species in figure 5(b) has a positive sign, which means that it is directed outward. Its values in gyro-Bohm units are indicated on the right-hand vertical axis of figure 5(b).
In figure 6(a), the density profile is shown for the same time. The particle flux is plotted in figure 6(b) with its values in gyro-Bohm units shown on the right-hand vertical axis. One sees that the particle flux is negative, in other words it is directed inwardly, in contrast to the energy flux. As a consequence of the inward turbulent particle flux, the density profile, figure 6(a), is modified with a ‘hump’ appearing in the inner part of the simulation domain and a ‘hole’ which the turbulence ‘digs’ in the outer part. Note that these results qualitatively resemble observations from W7-X and Large Helical Device where hollow density profiles, predicted by neoclassical theory (Maaßberg, Beidler & Simmet Reference Maaßberg, Beidler and Simmet1999; Yamagishi et al. Reference Yamagishi, Yokoyama, Nakajima and Tanaka2007; Beidler et al. Reference Beidler, Feng, Geiger, Köchl, Maßberg, Marushchenko, Nührenberg, Smith and Turkin2018), are not seen experimentally (Tanaka et al. Reference Tanaka, Kawahata, Tokuzawa, Akiyama, Yokoyama, Shoji, Michael, Vyacheslavov, Murakami and Wakasa2010; Wolf et al. Reference Wolf, Alonso, Äkäslompolo, Baldzuhn, Beurskens, Beidler, Biedermann, Bosch, Bozhenkov and Brakel2019). A possible explanation for these results is an inward turbulent particle pinch such as seen here numerically (also observed recently by Thienpondt et al. (Reference Thienpondt, Garcia-Regana, Calvo, Alonso, Velasco, Gonzalez-Jerez, Barnes, Brunner, Ford and Fuchert2022)). We observe the inward direction of the particle transport at all beta values including those at which W7-X has already been operated (Klinger et al. Reference Klinger, Andreeva, Bozhenkov, Brandt, Burhenn, Buttenschön, Fuchert, Geiger, Grulke and Laqua2019). Note that no particle source has been implemented in our simulations which may also affect the density profile evolution. The turbulent particle transport and the inward turbulent particle pinch in stellarator plasmas will be a subject of further studies in future. It seems probable that the numerical results can be understood in terms of the quasilinear theory (Helander & Zocco Reference Helander and Zocco2018). Another interesting application of our computations will be to consider the impurity minority transport by electromagnetic turbulence. This task is technically straightforward and may help understanding the absence of the impurity accumulation in W7-X (Langenberg et al. Reference Langenberg, Wegner, Marchuk, Garcıa-Regaña, Pablant, Fuchert, Bozhenkov, Damm, Pasch and Brunner2021) although it is predicted by the neoclassical theory. A possible explanation is the outward impurity transport by the turbulence. We will report on this issue in future publications.
3.2. Electromagnetic turbulence in a compact stellarator (W7-K)
As the next step, we consider a compact stellarator configuration, W7-K, which is shown in figure 7. This configuration has been proposed by Roberg-Clark et al. (Reference Roberg-Clark, Plunk and Xanthopoulos2022) in the course of the stellarator optimization effort minimizing ITG turbulence. It has a particularly high threshold for the linear destabilization of the ITGs and a strong drive of the zonal flow in the nonlinear regime. In contrast to the KJM configuration of W7-X, considered in the previous section, W7-K is Mercier unstable already at a relatively small beta. Therefore, the electromagnetic turbulence may be a more important issue for configurations of this type than the electrostatic ITG. In this paper, we show the first simulations of electromagnetic turbulence in the W7-K configuration.
We consider a W7-K configuration with the machine size $L_x = 488.5$ (measured in the sound gyroradii) and employ the same temperature profile (3.1), and the flat density as in the previous section. In figure 8, the dependency of the growth rate on the plasma beta is shown. One can see a characteristic transition of ITG instability to KBM instability. In figure 8(a), simulations including $\delta B_{\|}$ are compared with runs using the same parameters but without this field. One can see that there is a clear difference, in particular at higher values of beta, in agreement with the conventional wisdom that the inclusion of $\delta B_{\|}$ decreases the growth rate (Belli & Candy Reference Belli and Candy2010). The equilibrium used in figure 8(a) was fixed to the vacuum W7-K. This has been corrected in figure 8(b) where the results of the simulations applying consistent equilibria are shown, i.e. the same plasma profiles were used in the Variational Moments Equilibrium Code (VMEC) equilibrium computation as in the gyrokinetic code. One sees that this also leads to a different result. Unsurprisingly, the discrepancy between the simulations using the vacuum and the consistent equilibria increases with beta, reaching a 20 % difference (or even slightly more) at the highest beta values. These results are in agreement with the findings of Aleynikova et al. (Reference Aleynikova, Zocco, Xanthopoulos, Helander and Nührenberg2018) which have stressed importance of using consistent equilibria and including the parallel component of the magnetic field perturbation.
In figure 9, the temporal evolution of the averaged radial heat flux is shown in the electromagnetic ITG and in the KBM regimes. Again, one sees that the simulations reach the nonlinear stage. Note that the heat flux in W7-K, figure 9(a), is plotted at $\beta = 1.28\,\%$, corresponding to the electromagnetic ITG regime, as can be concluded from figure 8(a). One sees that for this beta value, the heat flux in W7-K is comparable to the heat flux in W7-X shown in figure 1(b) for $\beta = 0.28\,\%$ (electromagnetic ITG) and $\beta = 2.8\,\%$ (KBM). For W7-K, the heat flux becomes considerably larger in the KBM regime which requires a higher beta because of the higher value of the transition threshold in W7-K for the parameters considered. Higher values of the KBM transition threshold in W7-K compared with W7-X is an interesting feature. We hypothesize that this might be related to the specific optimization (Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2022) of W7-K seeking to increase the marginal stability of the ITGs compared with the W7-X values. It appears that this optimization may have a similar effect on the KBM marginal stability. This aspect needs a more thorough study in future.
In figure 10, the evolution of the electrostatic potential $\phi _n(s)$ is shown for the poloidal angle $\theta =0$ as a function of the flux surface and the toroidal mode number. Similarly to the W7-X simulations, described in the previous section, the turbulence starts in the linear phase with a coherent mode at the toroidal mode numbers around $n \approx 17$ computed for a single period. For the three-period W7-K geometry, this corresponds to the toroidal mode numbers around $\bar {n}_{{\rm tot}}=51$ when calculated for the whole torus. The linear instability evolves, spreading radially and generating zonal flows and other components with the lower toroidal mode numbers. This process leads to a saturation of the electromagnetic turbulence as in W7-X, see figure 4.
Finally, in figure 11, the evolution of the electrostatic potential is shown in the poloidal cross-section. One sees that the turbulence spreads both poloidally and radially in the course of its nonlinear evolution although it remains quite localized in the real space for the plasma profiles considered in this paper.
In future work, we will extend our simulations to other configurations featuring a very high ITG destabilization threshold (Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2022) with a possibly less benign Mercier stability. Characterizing transport in such stellarator geometries in comparison with the transport in the W7-X will provide valuable information for the choice of the magnetic geometry in optimized stellarator design.
4. Conclusions
First global fully electromagnetic simulations of the KBM turbulence in stellarators have been performed which include all components of the magnetic field perturbations. This extends earlier work (Mishchenko et al. Reference Mishchenko, Biancalani, Bottino, Hayward-Schneider, Lauber, Lanti, Villard, Kleiber, Könies and Borchardt2021; Wilms et al. Reference Wilms, Navarro, Merlo, Leppin, Görler, Dannert, Hindenlang and Jenko2021; Mishchenko et al. Reference Mishchenko, Bottino, Hayward-Schneider, Poli, Wang, Kleiber, Borchardt, Nührenberg, Biancalani and Könies2022) where the global electromagnetic ITG turbulence was simulated in W7-X neglecting the parallel component of the magnetic field perturbation.
Turbulence plays an important role in stellarators, particularly if the neoclassical losses have been reduced by optimization of the magnetic field. In this paper, we have studied the global electromagnetic turbulence using the gyrokinetic particle-in-cell code EUTERPE. The evolution of the turbulent electromagnetic field and the plasma profiles have been considered at varying plasma beta and for different stellarator magnetic configurations. It is found that the turbulence is linearly driven at relatively high toroidal mode numbers where it is excited through the choice of the initial conditions with the corresponding mode structure. In the later phase, lower mode numbers are generated in the perturbed field which may result from the nonlinear couplings (e.g. the zonal flows) but could also be a consequence of the linear stability properties in case of the KBMs which favour the lower part of the mode-number spectrum (Aleynikova et al. Reference Aleynikova, Zocco, Xanthopoulos, Helander and Nührenberg2018). The turbulent heat flux is outward and leads to a relaxation of the plasma temperature profile. The particle flux is inward for the parameters considered. The effect of the parallel perturbation of the magnetic field on the stellarator turbulence is addressed. This perturbation component can change the growth rate and the saturation level, making its inclusion potentially relevant for an accurate description of the KBM dynamics.
In future, we will study the electromagnetic turbulence in novel ITG-optimized stellarator configurations (Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2022) and consider the turbulent particle transport including the inward particle pinch (Tanaka et al. Reference Tanaka, Kawahata, Tokuzawa, Akiyama, Yokoyama, Shoji, Michael, Vyacheslavov, Murakami and Wakasa2010) and the impurity accumulation problem (Langenberg et al. Reference Langenberg, Wegner, Marchuk, Garcıa-Regaña, Pablant, Fuchert, Bozhenkov, Damm, Pasch and Brunner2021). Also, we plan to compare the role played by the parallel component of the magnetic perturbations in local and global KBM turbulence, and the relative importance of radially localized zonal-flow drive and non-local radial profile relaxation for the turbulence saturation.
The computational progress reported here has been made possible by a substantial increase in the computing power employed in the simulations. This increase in the resources applied resulted from EU's PRACE project EMGKPIC and EUROfusion's E-TASC TSVV projects. Typically, several thousands of central processing unit cores have been used for many days. Note that the future exascale computing will make it possible to perform simulations employing several thousands of computing nodes resulting in two orders of magnitude increase in the computing power. This highlights the importance of making the EUTERPE code fit for the exascale challenge which includes its graphics-processing-unit-enabling but also issues related to the massive parallel output and memory requirements.
Acknowledgements
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (grant agreement no. 101052200 – EUROfusion). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. Simulations presented in this work were performed on the MARCONI FUSION HPC system at CINECA. We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France.
Editor A. Schekochihin thanks the referees for their advice in evaluating this article.
Funding
This work was supported by the European Union via the Euratom Research and Training Programme (grant agreement no. 101052200 – EUROfusion).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available on request from the corresponding author.
Appendix A. Derivation of $\delta B_{\|}$ terms
Consider the gyrokinetic Ampere's law,
Following Brizard & Hahm (Reference Brizard and Hahm2007), we can write the perturbed current in the form
In our model, we neglect the higher-order terms, approximating
This expression can further be simplified expanding it around the gyrocentre position and assuming the gyroradius to be small (we keep the terms $\sim k_{\perp }^2 \rho ^2$),
Here, we have introduced the moments
Neglecting the pressure anisotropy $P_{1s\|} - P_{1s\perp } \approx 0$ and employing the approximation $\boldsymbol {\nabla } \times (\delta B_{\|} {\boldsymbol {b}}) \approx \boldsymbol {\nabla } \delta B_{\|} \times {\boldsymbol {b}}$, we obtain $\mu _0 \boldsymbol {\nabla }_{\perp } \sum _s P_{1s\perp } + B \boldsymbol {\nabla }_{\perp } \delta B_{\|} = 0$. Neglecting $\boldsymbol {\nabla }_{\perp }$ in this equation results in (2.11a,b) with $\delta p_{\perp } = \sum _s P_{1s\perp }$.
Consider now the perturbed equations of motion,
Here, $\delta B_{\|}$ enters through the contributions proportional to $\left \langle {\boldsymbol {v}}_{\perp }\boldsymbol {\cdot } {\boldsymbol {A}}_{\perp } \right \rangle$. Expanding this expression for the small gyroradius, we obtain $e \left \langle {\boldsymbol {v}}_{\perp } \boldsymbol {\cdot } {\boldsymbol {A}}_{\perp } \right \rangle \approx - \mu \delta B_{\|}({\boldsymbol {R}})$ at the same level of accuracy (${\sim }k_{\perp }^2 \rho ^2$) which was used to derive (2.11a,b). Substituting this result into (A7) and (A8), we obtain (2.4) and (2.5).