Global gyrokinetic simulations of electromagnetic turbulence in stellarator plasmas

Abstract

Global electromagnetic turbulence is simulated in stellarator geometry using the gyrokinetic particle-in-cell code EUTERPE. The evolution of the turbulent electromagnetic field and the plasma profiles is considered at different values of the plasma beta and for different magnetic configurations. It is found that turbulence is linearly driven at relatively high toroidal mode numbers. In the nonlinear regime, lower toroidal mode numbers, including zonal flows, are excited resulting in a quench of the linear instability drive. The turbulent heat flux is outward and leads to the nonlinear 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.

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. 2021). Confinement has been studied in Wendelstein 7-X (W7-X) experiments (Klinger et al. 2016; Dinklage et al. 2018; Wolf et al. 2019; Carralero et al. 2021), 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. 2020; Zocco et al. 2022). 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. (2020) and Sánchez et al. (2020). 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. 2021; Wilms et al. 2021; Mishchenko et al. 2022), 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. (2018) 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. 2013).

A fully electromagnetic extension of the global nonlinear gyrokinetic code EUTERPE (Kornilov et al. 2004) 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 2022). The parallel component of the magnetic field perturbation, known to be important for KBMs, especially for the large-scale modes (Aleynikova et al. 2018), 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. 2019) in the plasma core as expected from the neoclassical theory (Beidler et al. 2018). 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. (2014) for details on the derivation and notations.

The equations include the gyrokinetic Vlasov equation,

(2.1)\begin{equation} \frac{\partial f_{1s}}{\partial t} + \dot{{\boldsymbol{R}}} \boldsymbol{\cdot} \frac{\partial f_{1s}}{\partial {\boldsymbol{R}}} + \dot{v}_{\|} \frac{\partial f_{1s}}{\partial v_{\|}} = {}- \dot{{\boldsymbol{R}}}^{(1)} \boldsymbol{\cdot} \frac{\partial F_{0s}}{\partial {\boldsymbol{R}}} -\dot{v}_{\|}^{(1)} \frac{\partial F_{0s}}{\partial v_{\|}}. \end{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

(2.2a,b)\begin{gather} \dot{{\boldsymbol{R}}} = \dot{{\boldsymbol{R}}}^{(0)} + \dot{{\boldsymbol{R}}}^{(1)},\quad \dot{v}_{\|} = \dot{v}_{\|}^{(0)} + \dot{v}_{\|}^{(1)}, \end{gather}

(2.3a,b)\begin{gather}\dot{{\boldsymbol{R}}}^{(0)} = v_{\|} {\boldsymbol{b}}^*_0 + \frac{1}{q B_{\|}^*}{\boldsymbol{b}}\times \mu \boldsymbol{\nabla} B,\quad \dot{v}_{\|}^{(0)} ={-}\frac{\mu}{m}{\boldsymbol{b}}^*_0 \boldsymbol{\cdot} \boldsymbol{\nabla} B, \end{gather}

(2.4)\begin{gather}\dot{{\boldsymbol{R}}}^{(1)} = \frac{{\boldsymbol{b}}}{B_{\|}^*} \times \boldsymbol{\nabla} \left\langle \phi - v_{\|} A_{\|}^{(s)} - v_{\|} A_{\|}^{(h)} \right\rangle - \frac{q}{m} \langle A^{(h)}_{\|} \rangle {\boldsymbol{b}}^* + \frac{1}{q B_{\|}^*} {\boldsymbol{b}}\times \mu \boldsymbol{\nabla} \delta B_{\|}, \end{gather}

(2.5)\begin{align} \dot{v}_{\|}^{(1)} & ={-}\frac{q}{m} \left[{\boldsymbol{b}}^* \boldsymbol{\cdot}\boldsymbol{\nabla} \left\langle \phi - v_{\|} A_{\|}^{(h)} \right\rangle + \frac{\partial }{\partial t}\left\langle A_{\|}^{(s)} \right\rangle \right] \nonumber\\ & \quad - \frac{\mu}{m} \left( \frac{{\boldsymbol{b}} \times \boldsymbol{\nabla} B}{B_{\|}^*} \boldsymbol{\cdot}\boldsymbol{\nabla}\left\langle A_{\|}^{(s)} \right\rangle + {\boldsymbol{b}}^*_0 \boldsymbol{\cdot}\boldsymbol{\nabla}\delta B_{\|} \right), \end{align}

(2.6a,b)\begin{gather} {\boldsymbol{b}}^* = {\boldsymbol{b}}^*_0 + \frac{\boldsymbol{\nabla}\langle A_{\|}^{{(s)}} \rangle\times {\boldsymbol{b}}}{B_{\|}^*},\quad {\boldsymbol{b}}^*_0 = {\boldsymbol{b}} + \frac{m v_{\|}}{q B_{\|}^*} \boldsymbol{\nabla}\times{\boldsymbol{b}}, \end{gather}

(2.7)\begin{gather}B_{\|}^* = B + \frac{m v_{\|}}{q} {\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\times{\boldsymbol{b}}. \end{gather}

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:

(2.8a,b)\begin{gather} -\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{n_0}{B \omega_{{\rm ci}}} \boldsymbol{\nabla}_{{\perp}}\phi\right) = \bar{n}_{1{\rm i}} - \bar{n}_{1e},\quad \bar{n}_{1s} = \int f_{1s} \delta^3_{{\rm gy}} {\rm d}^6 Z, \end{gather}

(2.9a,b)\begin{gather}\left(\sum_{s = i,e}\frac{\beta_s}{\rho_{s}^{2}} - \nabla_{{\perp}}^{2} \right) A_{\|}^{(h)} = \mu_{0} \sum_{s = i,e} \bar{j}_{\|1s} + \nabla_{{\perp}}^2 A_{\|}^{(s)},\quad \bar{j}_{\|1s} = \int v_{\|} f_{1s} \delta^3_{{\rm gy}} {\rm d}^6 Z. \end{gather}

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$,

(2.10)\begin{equation} \frac{\partial }{\partial t}A_{\|}^{(s)} + {\boldsymbol{b}} \boldsymbol{\cdot}\boldsymbol{\nabla} \phi = 0, \end{equation}

and the perpendicular pressure balance condition for $\delta B_{\|}$,

(2.11a,b)\begin{equation} \delta B_{\|} \approx{-} \mu_0 \delta p_{{\perp}} / B,\quad \delta p_{{\perp}} = \sum_{s = i,e} \int \mu B f_{1s} \delta({\boldsymbol{R}} -{\boldsymbol{x}}) {\rm d}^6 Z, \end{equation}

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. 2014) is applied.

3. Simulations

3.1. Electromagnetic turbulence in W7-X

We consider a high-mirror (KJM) configuration (Geiger et al. 2014; Aleynikova, Zocco & Geiger 2022) 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 (2016),

(3.1)\begin{equation} \frac{{\rm d} \ln T_{{\rm i,e}}}{{\rm d} s} = \begin{cases} - 20 \left(\left.\dfrac{1}{4} - \right|s - \left.\dfrac{1}{2}\right|\right), & 0.25 \le s \le 0.75,\\ 0, & {\rm otherwise.} \end{cases} \end{equation}

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

(3.2a–c)\begin{equation} T_*=T_{e}(s = 0.5),\quad n_*=\frac{N_{\rm ph},i}{V},\quad B_*= \frac{1}{2{\rm \pi}} \int_{0}^{2{\rm \pi}} B(s=0, \varphi)\,{\rm d} \varphi \end{equation}

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\,\%$.

Figure 1. (a) Growth rate computed from the linear stage of the perturbed energy flux evolution. (b) Electron heat flux in gyro-Bohm units $Q_{{\rm GB}} = n_* T_* c_{s} (\rho _{s}/a)^2$ with the characteristic sound speed $c_{s}^2 = T_*/m_{i}$, the characteristic ion sound gyroradius $\rho _{s} = c_{s}/\omega _*$ and $\omega _* = q_{i} B_*/m_{i}$ the characteristic ion-cyclotron frequency. Here, $a$ is the horizontal (radial) half-extent of the simulation domain. For the definitions of $B_*$, $T_*$ and $n_*$, see the main text.

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.

Figure 2. (a) Temporal evolution and (b) radial structure of the electrostatic potential. Individual toroidal modes are plotted for $\beta = 2.8\,\%$. In (a), one can see that the turbulence saturates in time. In (b), it is shown that the mode structure is smooth and located well inside the simulation domain.

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).

Figure 3. (a) Energy of the perturbed electrostatic field and (b) the perturbed magnetic field for different $\beta$ values.

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.

Figure 4. The evolution of the electrostatic potential (contour plot) and the second radial derivative of its zonal component (white curve) plotted for $\beta = 2.8\,\%$. The ambient temperature profile is plotted, too (the 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 (computed for a single period of W7-X).

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).

Figure 5. (a) The instantaneous temperature profile plotted versus the electrostatic potential mode structure (colour plot). (b) The energy flux plotted at the same time ($\beta = 2.8\,\%$).

Figure 6. (a) Instantaneous density profile plotted versus the electrostatic potential mode structure (colour plot). (b) The particle flux plotted at the same time ($\beta = 2.8\,\%$).

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 1999; Yamagishi et al. 2007; Beidler et al. 2018), are not seen experimentally (Tanaka et al. 2010; Wolf et al. 2019). A possible explanation for these results is an inward turbulent particle pinch such as seen here numerically (also observed recently by Thienpondt et al. (2022)). 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. 2019). 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 2018). 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. 2021) 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. (2022) 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.

Figure 7. New compact stellarator configuration (W7-K) optimizing the plasma properties with respect to the electrostatic ITG turbulence. It has the major radius $R_0 = 1.98$ m, the minor radius $a = 0.49$ m, the volume $V = 9.6\,{\rm m}^3$ and $N_{\rm per} = 3$ periods. This configuration is Mercier-unstable at relatively low beta implying that the electromagnetic turbulence may become an issue.

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 2010). 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. (2018) which have stressed importance of using consistent equilibria and including the parallel component of the magnetic field perturbation.

Figure 8. (a) Growth rate in the W7-K configuration plotted as a function of beta. The case including $\delta B_{\|}$ is compared with the case without it. Vacuum W7-K equilibrium is used. (b) Growth rate in W7-K as a function of beta plotted for consistent equilibria. The growth rate is computed from the linear simulations with the dominant poloidal and toroidal mode numbers $m_0=85$ and $n_0=-51$.

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. 2022) 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.

Figure 9. The heat fluxes for the both species in the W7-K configuration for (a) $\beta = 1.28\,\%$ (electromagnetic ITG regime) and (b) $\beta = 6.4\,\%$ (KBM regime).

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.

Figure 10. Electrostatic potential evolution in W7-K plotted as a function of the flux surface and the toroidal mode number (at the poloidal angle $\theta =0$). In this simulation, $\beta = 1.28\,\%$.

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.

Figure 11. Electrostatic potential evolution in W7-K plotted in the real-space poloidal cross-section at the toroidal angle $\varphi =0$. Here, $\beta = 1.28\,\%$, vacuum equilibrium is used, and $\delta B_{\|}$ is included.

In future work, we will extend our simulations to other configurations featuring a very high ITG destabilization threshold (Roberg-Clark et al. 2022) 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. 2021; Wilms et al. 2021; Mishchenko et al. 2022) 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. 2018). 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. 2022) and consider the turbulent particle transport including the inward particle pinch (Tanaka et al. 2010) and the impurity accumulation problem (Langenberg et al. 2021). 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.