Semianalytical calculation of the zonal-flow oscillation frequency in stellarators

TL;DR

This paper introduces a semianalytical numerical method to efficiently and accurately calculate the zonal-flow oscillation frequency in stellarators, aiding the understanding of turbulence regulation in fusion devices.

Contribution

A novel semianalytical approach for calculating zonal-flow oscillation frequency in stellarators, improving speed and accuracy over existing gyrokinetic simulation methods.

Findings

The semianalytical method matches simulation results accurately.

It is faster than traditional gyrokinetic simulation-based calculations.

The approach provides physical insights into zonal-flow dynamics.

Abstract

Due to their capability to reduce turbulent transport in magnetized plasmas, understanding the dynamics of zonal flows is an important problem in the fusion programme. Since the pioneering work by Rosenbluth and Hinton in axisymmetric tokamaks, it is known that studying the linear and collisionless relaxation of zonal flow perturbations gives valuable information and physical insight. Recently, the problem has been investigated in stellarators and it has been found that in these devices the relaxation process exhibits a characteristic feature: a damped oscillation. The frequency of this oscillation might be a relevant parameter in the regulation of turbulent transport, and therefore its efficient and accurate calculation is important. Although an analytical expression can be derived for the frequency, its numerical evaluation is not simple and has not been exploited systematically so far. Here, a numerical method for its evaluation is considered, and the results are compared with those obtained by calculating the frequency from gyrokinetic simulations. This "semianalytical" approach for the determination of the zonal-flow frequency reveals accurate and faster than the one based on gyrokinetic simulations.

Tables

Table 1. Table 1: CPU time (total core hours) required to obtain the zonal-flow frequency with CAS3D-K, EUTERPE and GENE in the magnetic configurations considered in the paper. Note that the CAS3D-K and GENE calculations are radially local and the values in the table correspond to the time needed for the calculation at a single radial position. EUTERPE is radially global and the times in the table correspond to the full radius calculation.

	CAS3D-K	EUTERPE	GENE
LART (large ripple)	$3-4$	$\sim 2000$	—
LART (small ripple)	$2-3$	$\sim 10000$	—
W7-X SD	$2-3$	$\sim 2000$	$\sim 180$
W7-X HM	$2-3$	$\sim 2000$	—
TJ-II SD	$3-4$	$\sim 1000$	—
LHD SD, OS & IS	$1-3$	$\sim 3000$	—

Equations

$${\Omega }_{\mathrm{ZF}}\sim \frac{\overline{{\mathbf{v}}_{di}\cdot \nabla \psi }}{{\mathbf{v}}_{di}\cdot \nabla \psi }\frac{v_{ti}}{L}.$$

$$\frac{1}{L}\ll k_{\perp }\ll \frac{1}{{\rho }_{ti}},$$

$$\mathbf{B}=-{\Psi }_p^{\prime }(\psi )\nabla \psi \times \nabla (\zeta -q(\psi )\theta ),$$

$$v_{\parallel }(\psi ,\theta ,\alpha ,v,\lambda ,\sigma )=\sigma v\sqrt{1-\lambda B(\psi ,\theta ,\alpha )}.$$

$$\phi ={\phi }_k(\psi ,t){\mathrm{e}}^{\mathrm{i}{k}_{\psi }\psi }.$$

$$F_s(\mathbf{x},v,\lambda ,\sigma ,t)=F_{s0}(\mathbf{x},v)+F_{s1}(\mathbf{x},v,\lambda ,\sigma ,t)+O({\rho }_{ts\star }^2{F}_{s0}),$$

$$F_{s0}(\mathbf{x},v):=\frac{n_s(\psi (\mathbf{x}))}{(\sqrt{2\pi }{v}_{ts}(\psi (\mathbf{x})){)}^3}\exp \left(-\frac{v^2}{2v_{ts}^2(\psi (\mathbf{x}))}\right),$$

$$F_{s1}(\mathbf{x},v,\lambda ,\sigma ,t)=f_s(\mathbf{x},v,\lambda ,\sigma ,t){\mathrm{e}}^{\mathrm{i}{k}_{\psi }\psi (\mathbf{x})}.$$

$$h_s=f_s+\frac{Z_{s}e}{T_s}{\phi }_k{J}_0(k_{\perp }{\rho }_s)F_{s0},$$

$$\left({\partial }_t+v_{\parallel }\hat{\mathbf{b}}\cdot \nabla +\mathrm{i}{k}_{\psi }{\omega }_s\right)h_s=\frac{Z_{s}e}{T_s}{J}_0(k_{\perp }{\rho }_s)F_{s0}{\partial }_t{\phi }_k.$$

$${\mathbf{v}}_{ds}=\frac{v^2}{{\Omega }_s}\hat{\mathbf{b}}\times \left[(1-\lambda B)\hat{\mathbf{b}}\cdot \nabla \hat{\mathbf{b}}+\frac{\lambda }{2}\nabla B\right]$$

$$\sum _s\frac{Z_s^{2}e}{T_s}{n}_s{\phi }_k={⟨\sum _s{Z}_s\int J_{0s}{h}_s{\mathrm{d}}^{3}v⟩}_{\psi }.$$

$$⟨G{⟩}_{\psi }=V^{\prime }(\psi {)}^{-1}{\int }_0^1\mathrm{d}\theta {\int }_0^1\mathrm{d}\zeta \sqrt{g}G(\psi ,\theta ,\zeta ),$$

$$\frac{1}{L}\ll k_{\perp }\lesssim \frac{1}{{\rho }_{ts}}.$$

$$\frac{L{\mathbf{v}}_E\cdot \nabla \alpha }{v_{ts}}<k_{\perp }{\rho }_{ts}.$$

$${\hat{\phi }}_k(p)=\frac{{\sum }_s{Z}_s{\left\{\frac{1}{p+\mathrm{i}{k}_{\psi }\overline{{\omega }_s}}{\mathrm{e}}^{-\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}\overline{{\mathrm{e}}^{\mathrm{i}{k}_{\psi }{\delta }_s}{f}_s(0)}/F_{s0}\right\}}_s}{{\sum }_s\frac{Z_s^{2}e}{T_s}{\left\{1-\frac{p}{p+\mathrm{i}{k}_{\psi }\overline{{\omega }_s}}{\mathrm{e}}^{-\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}\overline{{\mathrm{e}}^{\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}}\right\}}_s},$$

$${\left\{Q\right\}}_s:={⟨\sum _{\sigma =-1}^1{\int }_0^{\infty }\mathrm{d}v{\int }_0^{1/B}\mathrm{d}\lambda \frac{\pi v^{2}B}{\sqrt{1-\lambda B}}Q(\psi ,\theta ,\alpha ,v,\lambda ,\sigma )F_{s0}⟩}_{\psi }$$

$$\overline{Q}:=\left\{\begin{array}[]{lcr}\langle B\,Q/|v_{\parallel}|\rangle_{\psi}/\langle B/|v_{\parallel}|\rangle_{\psi}&&\textrm{passing particles,\,\,}\vspace{0.2cm}\\ \omega_{b}\oint\mathrm{d}\theta\,Q/(v_{\parallel}\,\hat{\mathbf{b}}\cdot\nabla\theta)&&\textrm{trapped particles,}\end{array}\right.$$

$$\overline{v_{\parallel }\hat{\mathbf{b}}\cdot \nabla Q}=0$$

$${\omega }_s=\overline{{\omega }_s}+v_{\parallel }\hat{\mathbf{b}}\cdot \nabla {\delta }_s,$$

$$\mathbf{B}=I_t\nabla \theta -I_p\nabla \zeta +\tilde{\beta }\left(\psi ,\theta ,\zeta \right)\nabla \psi ,$$

$$\sqrt{g}=\frac{I_t{\Psi }_p^{\prime }-I_p{\Psi }_t^{\prime }}{B^2}.$$

$${\delta }_s=-\frac{I_p}{{\Psi }_p^{\prime }}{\rho }_{\parallel s}+\widetilde{{\delta }_s}.$$

$$\widetilde{{\delta }_s}=\left(\frac{I_p}{{\Psi }_p^{\prime }}-\frac{I_t}{{\Psi }_t^{\prime }}\right)\sum _{m,n\ne 0}\left(\frac{qn}{m+qn}\right){({\rho }_{\parallel s})}_{mn}{\mathrm{e}}^{2\pi \mathrm{i}(m\theta +n\zeta )},$$

$${\rho }_{\parallel s}=\sum _{m,n}({\rho }_{\parallel s}{)}_{mn}{\mathrm{e}}^{2\pi \mathrm{i}(m\theta +n\zeta )}.$$

$$\widetilde{{\delta }_s}={\int }_0^{\tau }{\tilde{\omega }}_{s\alpha }\mathrm{d}{\tau }^{\prime },$$

$${\tilde{\omega }}_{s\alpha }=\frac{q{I}_p-I_t}{{\Psi }_p^{\prime }}\left[{\tau }_b^{-1}{\partial }_{\alpha }{\rho }_{\parallel s}-\overline{{\tau }_b^{-1}{\partial }_{\alpha }{\rho }_{\parallel s}}\right].$$

$$\tau=\left\{\begin{array}[]{lcr}\int_{\theta_{b1}}^{\theta}|\tau_{b}|\,\mathrm{d}\theta^{\prime}&\mathrm{when}\quad\sigma>0\vspace{0.2cm}\\ \int_{\theta_{b1}}^{\theta_{b2}}|\tau_{b}|\,\mathrm{d}\theta-\int_{\theta_{b2}}^{\theta}|\tau_{b}|\,\mathrm{d}\theta^{\prime}&\mathrm{when}\quad\sigma<0,\end{array}\right.$$

$$f_s(0)=\frac{Z_{s}e}{T_s}\frac{{⟨1-{\Gamma }_{0s}⟩}_{\psi }}{{\Gamma }_{0s}}{J}_{0s}{F}_{s0}{\phi }_k(0)$$

$${\hat{\phi }}_k(p)=\frac{{\sum }_s\frac{Z_s^2}{T_s}{\left\{\frac{1}{p+\mathrm{i}{k}_{\psi }\overline{{\omega }_s}}{\mathrm{e}}^{-\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}\overline{{\mathrm{e}}^{\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}{⟨1-{\Gamma }_{0s}⟩}_{\psi }/{\Gamma }_{0s}}\right\}}_s}{{\sum }_s\frac{Z_s^2}{T_s}{\left\{1-\frac{p}{p+\mathrm{i}{k}_{\psi }\overline{{\omega }_s}}{\mathrm{e}}^{-\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}\overline{{\mathrm{e}}^{\mathrm{i}{k}_{\psi }{\delta }_s}{J}_{0s}}\right\}}_s}{\phi }_k(0).$$

Full text

1 Introduction

The reduction of turbulent transport by shearing of zonal flows in tokamaks and stellarators has been studied for decades now [1, 2], and their presence has been observed experimentally in several fusion devices [3]. These flows are associated to electrostatic potential perturbations constant on flux surfaces, which are generated from drift-waves by inverse cascade of energy in a turbulent plasma. Although the generation of zonal flows by plasma turbulence is a non-linear process, the study of its linear relaxation provides insight into the problem at an affordable cost. Furthermore, the linear evolution of zonal flows could be relevant for the saturation and regulation of the turbulence in some situations; e.g. under marginal stability conditions in which the turbulence drive, and consequently the zonal-flow drive, is small.

The linear evolution of zonal flows is greatly influenced by the magnetic geometry. Therefore, understanding its behavior in different geometries could help in the design of future devices optimized for enhanced zonal flows and reduced turbulent transport. Rosenbluth and Hinton [4] studied the evolution of an initial zonal potential perturbation with small wavenumber in tokamak geometry, showing that the perturbation is not completely damped by collisionless processes but reaches a finite value at long times, the so-called zonal-flow residual level. Since then, the long-time collisionless evolution of zonal potential perturbations111That is, electrostatic potential perturbations that are constant on flux surfaces, while having a finite radial scale. has attracted the attention of different authors that continued the analysis of the problem in axisymmetric configurations [5, 6, 7] and more recently in stellarators [8, 9, 10, 11, 12, 13, 14]. In [14] a method for the fast numerical calculation of the residual level in tokamaks and stellarators for arbitrary wavenumbers of the perturbation was presented.

In general, the collisionless evolution of an initial zonal-flow perturbation in a stellarator involves two oscillations with different typical values of the frequency and also different physical origin. On the one hand, a decaying geodesic acoustic mode (GAM) oscillation [15], mainly caused by the dynamics of passing ions (at least in large aspect ratio devices). The collisionless damping of the GAM strongly depends on the safety factor [8, 10, 16] and the frequency, $\Omega_{\rm GAM}$, is on the order of the characteristic frequencies of the turbulence; that is, $\Omega_{\rm GAM}\sim v_{ti}/L$. Here, $v_{ti}$ is the ion thermal speed and $L$ is a characteristic macroscopic scale of the system222Typically, $L\sim R$, the major radius of the torus.. On the other hand, in [11] it was found that in stellarator geometry an additional oscillation takes place. Its frequency, that we denote by $\Omega_{\rm ZF}$, is significantly smaller than $\Omega_{\rm GAM}$. Unlike the GAM oscillation, that happens in tokamaks and stellarators, this slower oscillation is characteristic of stellarator geometries333It is more precise to say that the oscillation with frequency $\Omega_{\rm ZF}$ exists only in non-omnigeneous devices. In particular, it can also manifest in rippled tokamaks. This will be explained below. and has been experimentally observed in [17]. If we denote by $\mathbf{v}_{di}$ the ion magnetic drift and by $\psi$ a radial coordinate,

[TABLE]

Here, $\overline{\mathbf{v}_{di}\cdot\nabla\psi}$ stands for a typical value of the radial magnetic drift averaged over a trapped trajectory. Since, in general, $\overline{\mathbf{v}_{di}\cdot\nabla\psi}/\mathbf{v}_{di}\cdot\nabla\psi\ll 1$, we have $\Omega_{\rm ZF}\ll\Omega_{\rm GAM}$. From now on, when we refer to a zonal flow oscillation (or simply to an oscillation), and if not stated otherwise, we will understand that we are talking about the low frequency one, which is the main subject of this paper.

The interest in the computation of $\Omega_{\rm ZF}$ resides in the potential role of the oscillation for the regulation of turbulent transport, pointed out in [13]. Since the residual value of the zonal flow is reached at times longer than the typical saturation time of the turbulence, the oscillatory phase of the zonal flow relaxation is likely to be physically more important than the value achieved when $t\to\infty$.

An analytical expression for $\Omega_{\rm ZF}$ was first derived in [11]. The oscillation was further studied, including its damping, in [12]. As we will see in Section 2, the derivation of the expression for $\Omega_{\rm ZF}$ employs a local gyrokinetic equation. The expansions used to arrive at the expression for $\Omega_{\rm ZF}$ require

[TABLE]

where $k_{\perp}$ is the wavenumber of the zonal potential perturbation, $\rho_{ti}=v_{ti}/\Omega_{i}$ is the thermal ion gyroradius, $\Omega_{i}=Z_{i}eB/m_{i}$ is the ion gyrofrequency, $Z_{i}e$ is the ion charge, $e$ is the proton charge, $B$ is the magnitude of the magnetic field $\mathbf{B}$ and $m_{i}$ is the ion mass. In reference [12], it was pointed out that it is unclear whether or not the analytical expression for the frequency is quantitatively accurate in actual devices. Anyway, the accuracy of the analytical expression has not been systematically checked so far, mainly because its evaluation is non-trivial: it involves phase-space averages that cannot be computed analytically in stellarator geometry, and therefore the evaluation must be carried out numerically. Here, we investigate this “semianalytical” method to determine the value of $\Omega_{\rm ZF}$ by using an extension [14] of CAS3D-K [18, 19]. The results of the semianalytical calculation are compared to those obtained from gyrokinetic simulations with the radially global code EUTERPE [20, 21] and the radially local code GENE [7, 22, 23, 24] in the W7-X, TJ-II and LHD stellarators and also in a series of rippled tokamak equilibria.

The paper is organized as follows. In Section 2 we derive the expression for the frequency of the zonal flow oscillation, $\Omega_{\rm ZF}$. In Section 3, we describe the numerical methods to evaluate this expression with the code CAS3D-K and we comment on how to obtain $\Omega_{\rm ZF}$ from gyrokinetic simulations with the global gyrokinetic code EUTERPE and the full flux-surface version of GENE. In Section 4 we compute $\Omega_{\rm ZF}$, with both methods (that is, the semianalytical approach and the one based on direct gyrokinetic simulations), in a tokamak configuration with different ripple values, which we use as a test to gain insight into the problem and check the limitations of the numerical methods. As explained in Section 4, we will use an ‘academic’ tokamak configuration, with unrealistic safety factor profile, so that in the simulations the zonal flow oscillation is neatly observed and we can clearly illustrate how $\Omega_{\rm ZF}$ varies as a function of the ripple size. In Section 5 we calculate the zonal-flow frequency in the W7-X, TJ-II and LHD stellarators. The computational requirements of the two approaches to the calculation of $\Omega_{\rm ZF}$ are compared in Section 6. Finally, the conclusions are given in Section 7.

2 Derivation of the analytical expression for the zonal-flow oscillation frequency

In this section we briefly explain how to calculate the long-time, linear and collisionless evolution of a zonal potential perturbation. More details can be found in [14]. The focus here will be on the derivation of an expression for $\Omega_{\rm ZF}$ and on the approximations made to obtain it. In later sections the numerical evaluation of this expression is compared to gyrokinetic simulations for several toroidal devices.

2.1 Long-time collisionless relaxation of zonal electrostatic perturbations

We use straight-field-line coordinates $\{\psi,\theta,\zeta\}$ to locate a point in space, $\mathbf{x}(\psi,\theta,\zeta)$. Here, $\psi=\Psi_{t}/\Psi_{t}^{\mathrm{edge}}$ is the radial coordinate, defined as the toroidal flux $\Psi_{t}$ normalized by its value at the last closed flux surface $\Psi_{t}^{\mathrm{edge}}$, and $\theta$ and $\zeta$ are, respectively, poloidal and toroidal angles normalized such that $\theta,\zeta\in[0,1)$. In these coordinates, the contravariant form of the magnetic field reads

[TABLE]

where $q(\psi)=\Psi^{\prime}_{t}(\psi)/\Psi^{\prime}_{p}(\psi)$ is the safety factor and $\Psi_{p}^{\prime}(\psi)$ and $\Psi_{t}^{\prime}(\psi)$ are the derivatives of the poloidal and toroidal fluxes with respect to $\psi$. It is convenient to define the coordinate $\alpha=\zeta-q(\psi)\theta$, that labels magnetic field lines on a surface once $\psi$ has been fixed. Unless otherwise specified, below we use $\{\psi,\theta,\alpha\}$ as independent spatial coordinates. Observe that, given $\psi$ and $\alpha$, the coordinate $\theta$ locates the point along the field line. The independent velocity coordinates (all functions in this paper are independent of the gyrophase) are $\{v,\lambda,\sigma\}$, where $v$ is the magnitude of the velocity, $\lambda=B^{-1}v_{\perp}^{2}/v^{2}$ is the pitch-angle coordinate, $v_{\perp}$ is the component of the velocity perpendicular to the magnetic field and $\sigma=v_{\parallel}/|v_{\parallel}|$ is the sign of the parallel velocity,

[TABLE]

We are interested in studying the linear and collisionless evolution of an electrostatic potential perturbation of the form

[TABLE]

Let us denote the phase-space distribution function of species $s$ by $F_{s}=F_{s}(\mathbf{x},v,\lambda,\sigma,t)$ and the thermal gyroradius by $\rho_{ts}=v_{ts}/\Omega_{s}$, where $v_{ts}=\sqrt{T_{s}/m_{s}}$, $T_{s}$ and $m_{s}$ are the thermal speed, temperature and mass of species $s$, and $\Omega_{s}=Z_{s}eB/m_{s}$ and $Z_{s}e$ are the gyrofrequency and charge of species $s$. We expand $F_{s}$ in $\rho_{ts\star}=\rho_{ts}/L\ll 1$,

[TABLE]

where $F_{s1}\sim O(\rho_{ts\star}F_{s0})$ and $F_{s0}$ is a Maxwellian distribution with density $n_{s}(\psi)$ and temperature $T_{s}(\psi)$ constant on flux surfaces. That is,

[TABLE]

with $\sum_{s}Z_{s}n_{s}=0$ due to quasineutrality. For the sake of consistence with (5), we also assume that the perturbation to the Maxwellian can be expressed as

[TABLE]

Then, the collisionless gyrokinetic equation, written in terms of the non-adiabatic component

[TABLE]

reads [14]

[TABLE]

Here, $k_{\perp}:=k_{\psi}|\nabla\psi|$, $\hat{\mathbf{b}}$ is the unit vector along the magnetic field, $J_{0}$ is the zeroth-order Bessel function of the first kind and $\rho_{s}=v_{\perp}/\Omega_{s}$ is the gyroradius of species $s$. Below, we often use the short-hand notation $J_{0s}\equiv J_{0}(k_{\perp}\rho_{s})$. Finally, $\omega_{s}:=\mathbf{v}_{ds}\cdot\nabla\psi$ is the radial magnetic drift frequency of species $s$, where

[TABLE]

is the magnetic drift velocity.

Equation (10) is solved along with the flux-surface averaged quasineutrality equation

[TABLE]

The flux-surface average is defined, for any given function $G=G(\psi,\theta,\zeta)$, as

[TABLE]

where $\sqrt{g}=[\nabla\psi\cdot(\nabla\theta\times\nabla\zeta)]^{-1}$ is the square root of the metric determinant and $V^{\prime}(\psi)=\int_{0}^{1}\mathrm{d}\theta\int_{0}^{1}\mathrm{d}\zeta\sqrt{g}$ is the derivative of the volume enclosed by the flux surface labeled by $\psi$.

Note that equation (10) is based on an eikonal representation of the fields (see (5) and (8)), where the variation on small scales is contained in $\exp(\mathrm{i}k_{\psi}\psi)$, and $\varphi_{k}$, $f_{s}$ and $h_{s}$ vary on the macroscopic scale $L$. Equation (10) is valid as long as [14]

[TABLE]

Before proceeding to give the solution of equations (10) and (12) at long times, it seems timely to discuss why the role of the background radial electric field, $E_{\psi}$, has not been considered. The background radial electric field enters the collisionless gyrokinetic equation via terms like $\mathbf{v}_{E}\cdot\nabla\alpha\partial_{\alpha}h_{s}$, where $\mathbf{v}_{E}\cdot\nabla\alpha=E_{\psi}/\Psi^{\prime}_{p}$ is the component of the $E\times B$ drift $\mathbf{v}_{E}$ in the direction $\alpha$. This term is smaller than the last term on the left side of (10) as long as

[TABLE]

Hence, when $k_{\perp}\rho_{ts}\sim 1$ the effect of $E_{\psi}$ can be safely neglected444This is true in the strictly collisionless limit and for the time scales considered in this paper. It is known that the background radial electric field is essential to calculate the neoclassical equilibrium in the $\sqrt{\nu}$ collisionality regime, for example [26].. Of course, when $k_{\perp}$ approaches the value $1/L$ (but recall that (14) has to be satisfied), the background radial electric field can introduce corrections to the results that we will present in this article. However, a rigorous treatment of these corrections exceeds the scope of this work. In particular, if $E_{\psi}$ is included, the partial differential equation to be solved has higher dimensionality because differential terms in $\alpha$ appear. Some analytical progress in this direction is provided in [25], using a model for stellarator geometry and making some additional simplifying assumptions. In what follows, we restrict ourselves to the set of equations consisting of (10) and (12).

In order to solve equations (10) and (12) it is convenient to use the Laplace transform, defined for any function $Q(t)$ as $\widehat{Q}(p)=\int_{0}^{\infty}Q(t)\mathrm{e}^{-pt}\mathrm{d}t$, where $p$ denotes the variable in Laplace space. Solving for $p/(v_{ts}L^{-1})\ll 1$, we find that [14]

[TABLE]

where $f_{s}(0)\equiv f_{s}(\mathbf{x},0)$ is the initial condition for $f_{s}$ and we have simplified the notation by defining the phase-space integration

[TABLE]

for any phase-space function $Q$.

The derivation of (16), and the expression (16) itself, makes use of the orbit average, defined for a function $Q(\psi,\theta,\alpha,v,\lambda,\sigma,t)$ by

[TABLE]

where $\omega_{b}:=[\oint\mathrm{d}\theta/(v_{\parallel}\,\hat{\mathbf{b}}\cdot\nabla\theta)]^{-1}$ is the bounce frequency. The orbit average removes frequencies that are $O(v_{ts}/L)$ or higher, and in particular it removes the GAM oscillation. Here, the symbol $\oint$ stands for integration (at fixed $\psi$ and $\alpha$) over the trapped trajectory.

The quantity $\delta_{s}$ entering (16) is defined as follows. The orbit average acting on the parallel streaming operator has the useful property

[TABLE]

for any single-valued function $Q$. Using this property, we can write the radial magnetic drift frequency as

[TABLE]

where $\delta_{s}=\delta_{s}(\psi,\theta,\alpha,v,\lambda,\sigma)$, that we choose to be odd in $\sigma$, is the radial displacement of the particle’s gyrocenter from its mean flux surface.

Note that in axisymmetric tokamaks and in omnigenous stellarators $\overline{\omega_{s}}=0$ holds for all trajectories; in a tokamak with ripple and in non-omnigenous stellarators, one can only guarantee $\overline{\omega_{s}}=0$ for passing particles. In subsection 2.2 we explain how to calculate $\delta_{s}$ from (20).

2.2 Calculation of $\delta_{s}$

In this work, we have employed and implemented in CAS3D-K a faster method to compute $\delta_{s}$ for trapped trajectories with respect to that used in [14]. We point out this below.

If we take $\theta$ and $\zeta$ to be Boozer angles, the contravariant form of $\mathbf{B}$ is given by (3) and its covariant form is given by

[TABLE]

with $I_{t}=I_{t}(\psi)$ and $I_{p}=I_{p}(\psi)$ the toroidal and poloidal currents, respectively. It is straightforward to realize that the metric determinant $\sqrt{g}=[\nabla\psi\cdot(\nabla\theta\times\nabla\zeta)]^{-1}$ can be written as

[TABLE]

We split the solution to (20) as

[TABLE]

Following [14], we have that for passing particles

[TABLE]

where

[TABLE]

For trapped particles

[TABLE]

with

[TABLE]

Here, we have used the definition $\tau_{b}:=B\sqrt{g}/(v_{\parallel}\Psi^{\prime}_{p})$ and the parameter $\tau$ given by

[TABLE]

which is monotonic over the orbit. Finally, $\theta_{b1}$ and $\theta_{b2}$ are the bounce points of the orbit; that is, the solutions for $\theta$ of the equation $1-\lambda B(\psi,\theta,\alpha)=0$. For the calculations of this paper we have implemented in CAS3D-K the expression (26), instead of using an expansion in bounce harmonics as in [14].

2.3 Expression for the zonal flow oscillation frequency

Using the initial condition555The only relevant feature of this initial condition, as far as the calculation of the zonal-flow frequency is concerned, is that $f_{s}(0)\approx(Z_{s}e/T_{s})k_{\perp}^{2}\rho_{ts}^{2}F_{s0}\varphi_{k}(0)$ for $k_{\perp}^{2}\rho_{ts}^{2}\ll 1$.

[TABLE]

in (16), where $\Gamma_{0}\left(k_{\perp}^{2}\rho_{ts}^{2}\right):=\exp{\left(-k_{\perp}^{2}\rho_{ts}^{2}\right)}\,I_{0}\left(k_{\perp}^{2}\rho_{ts}^{2}\right)$ and $I_{0}$ is the zeroth order modified Bessel function, equation (16) reads

[TABLE]

The residual value of the electrostatic potential is obtained [14] from (30) by recalling the property of the Laplace transform $\lim_{t\to\infty}\varphi_{k}(t)=\lim_{p\to 0}p\widehat{\varphi}_{k}(p)$. Here, our goal is the derivation of an expression for $\Omega_{\rm ZF}$. This is found by expanding equation (30) in powers of $k_{\psi}\overline{\omega_{s}}/p\ll 1$, and keeping the lowest order terms in $k_{\perp}\rho_{ts}\sim k_{\psi}\delta_{s}\ll 1$. A direct check shows that, with these approximations,

[TABLE]

where

[TABLE]

and

[TABLE]

Taking the inverse Laplace transform of equation (31), we obtain the same result as in references [11, 12], namely

[TABLE]

where the zonal-flow frequency is given by

[TABLE]

and the amplitude of the zonal-flow oscillation can be written as

[TABLE]

Note that the expressions for the zonal-flow frequency and its amplitude do not depend on the radial wavenumber $k_{\psi}$, and that $A_{\mathrm{ZF}}$ is independent of the temperatures and the masses of the species.

From equation (36), the estimate (1) is now understood. Clearly, a reduction of the bounce-averaged magnetic drift frequency, $\overline{\omega_{s}}$, leads to a decrease of the zonal-flow frequency, $\Omega_{\mathrm{ZF}}$. This allows to have larger values of the zonal flow in time scales comparable to those of the turbulence, as pointed out in [13]. Therefore, reducing $\overline{\omega_{s}}$ not only leads to a reduction of neoclassical transport but it might also be beneficial for the reduction of turbulent transport via zonal flows. This is one of the reasons that make the calculation of $\Omega_{\mathrm{ZF}}$ an interesting physical problem. The simultaneous optimization of neoclassical and turbulent transport has been addressed as well in [27, 8, 28, 29].

Whereas the frequency is a quantity that can be meaningfully compared with gyrokinetic simulations, the situation is different regarding the amplitude of the oscillation. The problem is that in practice the amplitude is damped [12], and in the analytical calculation above we have not considered the damping. For this reason, and due to the fact that $\Omega_{\rm ZF}$ seems to be more physically relevant, we focus on it in what follows.

We turn to the validity of the expression derived for $\Omega_{\rm ZF}$. The result (36) is consistent with the orderings enumerated just before (31) because

[TABLE]

If we recall (14), we deduce that (36) is correct as long as

[TABLE]

as advanced in the Introduction.

The zonal-flow frequency (36) and the amplitude (37) were first derived in [11] and the problem was analyzed in more detail in [12], where it was pointed out that (39) (or, for what matters, (2)) might be difficult to satisfy in actual devices, and therefore the usefulness of (36) was unclear in quantitative terms. In this work, we address a systematic comparison of the right side of (36) with the value of $\Omega_{\rm ZF}$ obtained through gyrokinetic simulations. It turns out that (36) is accurate (although there are a number of nuances about the assessment of the accuracy, depending on the device and associated to both the analytical calculation and the fitting of the frequency from gyrokinetic simulations, that are pointed out in the following sections) and that the numerical methods that we present to evaluate the right side of (36) are faster than the determination of the frequency from gyrokinetic codes.

In the gyrokinetic simulations of this paper, plasmas consisting of one ion species and adiabatic electrons are employed. Then, in order to compare the semianalytical calculations with the results from gyrokinetic codes, the sum over species in equations (32), (33) and (34) is then limited to one term, $s=i$.

3 Two kinds of numerical tools for the calculation of the zonal-flow oscillation frequency

The calculation of the zonal-flow frequency through (36) involves the evaluation of the right sides of (32), (33) and (34). In general, these expressions cannot be easily evaluated analytically. We use the code CAS3D-K [18, 19, 14] to evaluate them numerically. In order to prove the accuracy of (36) to determine the oscillation frequency in actual stellarator configurations, we compare its evaluation with CAS3D-K against the frequency obtained from simulations with the radially global gyrokinetic code EUTERPE [20, 21] and the radially local, full flux-surface version of GENE [7, 22, 23, 24].

The code CAS3D-K allows to perform averages of phase-space functions over the lowest order trajectories. As these trajectories lie on flux surfaces, the averages involve basically two-dimensional integrations over the spatial coordinates. We use $\{\psi,\theta,\zeta\}$ as independent spatial coordinates for passing particles and $\{\psi,\theta,\alpha\}$ for trapped particles, and $\{v,\lambda,\sigma\}$ as independent velocity space coordinates. In particular, the averaging operations implemented in CAS3D-K are the orbit-average, described in equation (18), and the phase-space average operation which is given in (17). A more detailed description of how these averages are performed in CAS3D-K can be found in reference [14]. In this work, the calculation of $\delta_{s}$ for trapped particles has been significantly simplified in CAS3D-K by using the direct integration given by equation (26) instead of the Fourier method originally implemented in CAS3D-K and described in [14]. The values of $\delta_{s}$ obtained with both methods are the same while the new one is much simpler to implement and faster.

EUTERPE is a global $\delta f$ gyrokinetic code in 3D geometry with a Lagrangian Particle In Cell (PIC) scheme. The gyrokinetic simulations carried out with EUTERPE in this work are linear and collisionless with a plasma that consists of singly charged ions and adiabatic electrons with equal temperatures. The simulations are initiated with a zonal potential perturbation, which is produced by taking a perturbed distribution function for the ions

[TABLE]

Here, $F_{i0}$ is a Maxwellian distribution and $\epsilon$ is a small factor (on the order of $10^{-3}$) that makes the perturbation to the distribution function, $F_{i1}$, much smaller than the equilibrium distribution, $F_{i0}$.

A long-wavelength approximation is used in the quasineutrality equation in this work. Namely, the $\Gamma_{0}(x)=\mathrm{e}^{-x}I_{0}(x)$ function is approximated as $\Gamma_{0}(x)\approx 1-x$, with $x=k_{\perp}^{2}\rho_{ti}^{2}$ and $I_{0}$ being the modified Bessel function, which is valid for $k_{\perp}\rho_{ti}\ll 1$. The initial condition (40) is a good approximation of (29) under these conditions. Note also that the simulation is carried out with a fixed value of $k_{\psi}$, so that $k_{\perp}\rho_{ti}$ has a radial dependence due to the radial variation of the magnetic field. The radial variation of $k_{\perp}\rho_{ti}$ is not relevant, however, because the oscillation frequency does not depend on the radial scale (see Section 2). This was confirmed by carrying out several simulations for different values of $k_{\psi}$, from $0.5\pi$ to $4.5\pi$, with the same configuration and same numerical parameters.

With the initial condition (40) for ions, the collisionless simulation is linearly evolved, retaining just a few ($5$ to $10$) toroidal and poloidal Fourier modes of the potential. The $m=0,n=0$ mode, which we can identify with $\varphi_{k}$ in equation (5), is the dominant component of the potential spectrum during the simulation, thus proving that the assumption of an eikonal form in equation (5) is appropriate. Here, $m$ and $n$ label poloidal and toroidal modes, respectively. The time evolution of $\varphi_{k}$ at a number of radial positions is tracked.

To obtain the frequency of the zonal-flow oscillation, the time trace of the zonal potential, or its radial derivative666Assuming the eikonal form of the potential in equation (5), both the normalized potential and its radial derivative have a similar time evolution, with the same oscillation frequency. (normalized to its initial value), is fitted to a model function with the form

[TABLE]

where $\Omega_{\mathrm{ZF}}$ is the zonal-flow oscillation frequency, $A_{\mathrm{ZF}}$ is the amplitude of the oscillation, $\gamma_{\mathrm{ZF}}$ is the damping rate and $R_{\mathrm{ZF}}$ is the residual level, and ′ means derivative with respect to the radial coordinate (normalized toroidal flux in EUTERPE). The last term in equation (41) accounts for the decay to the zonal-flow residual level, to which the oscillations are superimposed. This term is important for cases in which the decay to this residual level is slow; e.g. in tokamak configurations with a small ripple added (see Section 4). The fit is performed with the non-linear fitting routine fit of the MATLAB software package. Note that, in principle, a similar value of the zonal-flow frequency could be obtained from a FFT of the time signal, but this method shows to be less precise in practice.

There are sources of error and uncertainty in the oscillation frequency obtained by this procedure, related to the simulation itself and the approximations in the gyrokinetic code, and also associated to the fitting process. Several dynamics are mixed in the zonal flow relaxation, which makes the fitting process non trivial. First, there is a decay to the residual level and, superimposed to it, several oscillations appear with different amplitudes and characteristic times (see figures 4 and 7). Since we are interested in the low frequency oscillation, the faster GAM oscillation, usually appearing at the beginning of the relaxation, is excluded in the fitting process. A very similar value is obtained if instead of using the $m=0,n=0$ component of the potential its flux-surface average is used777Note that the $m=0,n=0$ component of the potential is not strictly equal to the flux-surface average of the potential, in general..

Although the comparison of the frequency calculations is focused on the codes CAS3D-K and EUTERPE, calculations for the W7-X standard configuration have also been carried out with the code GENE. This is an Eulerian gyrokinetic $\delta f$ code which can be run in radially global, full flux surface or flux tube simulation domains. In the GENE simulations shown in this work, we use the full-flux surface version and adiabatic electrons. In this version the GENE code is spectral in the radial coordinate, while a finite difference scheme is used for the coordinates along the flux surface. The distribution function is initialized according to (29). From the GENE simulations, as in the case of EUTERPE, the frequency is obtained by fitting the normalized zonal potential time trace to a model like (41).

In all cases we start from a magnetic equilibrium calculated with the 3D MHD equilibrium code VMEC [30]. Then, the relevant equilibrium quantities are mapped from VMEC to PEST coordinates used in EUTERPE by means of an intermediate code; from VMEC to GENE coordinates by means of the GIST [24] package; and from VMEC to Boozer coordinates used in CAS3D-K by employing the MC3D code which is part of the CAS3D [31] code package, and also part of CAS3D-K.

For convenience, we use flat density and temperature profiles and assume a plasma with adiabatic electrons and singly charged ions in all cases. Although the normalized toroidal flux, $\psi$, is used as a radial coordinate in the codes, in what follows we show the results in the more natural coordinate $r/a:=\sqrt{\psi}$, where $a$ is the minor radius of the device.

4 Zonal-flow oscillation frequency in tokamaks with ripple

In order to understand the influence of the magnetic configuration on the zonal-flow oscillation, we take an axisymmetric equilibrium and modify it by adding different amounts of ripple. The resulting set of configurations allows us to study the dependence of the frequency on the ripple size and calibrate the calculation methods in a controlled path continuously departing from the axisymmetric case. As we explain below, the axisymmetric tokamak configuration employed here is completely academic because we use a safety profile that is not realistic. The reason for this is that, in this section, we are interested in showing in a pedagogical way how $\Omega_{\rm ZF}$ increases as the magnitude of the axisymmetry-breaking terms grows, both using the semianalytical method and gyrokinetic simulations. Regarding the latter, the safety profile chosen helps remove the GAM oscillation, and therefore allows to clearly observe the zonal flow oscillation and to measure its frequency.

We start with the VMEC input for the equilibrium of an axisymmetric large aspect ratio tokamak (LART) with major radius $R=5$ m and minor radius $a=0.5$ m, and modify the plasma boundary to include ripple. The ripple is generated by adding a poloidally symmetric perturbation to the boundary shape through non-zero coefficients $(R_{\rm BC})_{0,2}$ and $(Z_{\rm BS})_{0,2}$. We always take $(R_{\rm BC})_{0,2}=(Z_{\rm BS})_{0,2}$. Here, $(R_{\rm BC})_{m,n}$ and $(Z_{\rm BS})_{m,n}$ are the cosine and sine components of the coordinates $R_{\rm BC}$ and $Z_{\rm BS}$ at the boundary, respectively, and $m$ and $n$ are poloidal and toroidal mode numbers (see [32]). We calculate the VMEC equilibrium employing 36 Fourier cosine components888$|B|$ can be written as a Fourier cosine series with coefficients $B_{m,n}$ as $|B|=\sum_{m=0}^{3}\sum_{n=-4}^{4}B_{m,n}\cos(m\theta-n\zeta)$. for $B$; specifically, $0\leq|n|\leq 4$ and $0\leq m\leq 3$. After including the $(Z_{\rm BS})_{0,2}\neq 0$ perturbation to the boundary, several modes $B_{0,n}$, $n\neq 0$, are generated, whose size increases with the size of $(Z_{\rm BS})_{0,2}$. The non-axisymmetric perturbation also modifies, although only slightly, the axisymmetric components $B_{00}$, $B_{10}$, $B_{20}$ and $B_{30}$. We show in figure 1 some of the Fourier cosine coefficients of $B$ at the radial position $r/a=0.8$ for different values of the perturbation $(Z_{\rm BS})_{0,2}$. In the figure, we give the axisymmetric components, the largest non-axisymmetric coefficients, $B_{12}$ and $B_{14}$, and also $\sqrt{\sum_{n\neq 0}B_{mn}^{2}}$, where the sum includes all non-axisymmetric components. For the cases under study, we define the ripple by

[TABLE]

This parameter will be used as a measure of the perturbation to the axisymmetric equilibrium.

The average magnetic field strength at the axis is $B_{0}\simeq 2.48$ T for the non-perturbed case and decreases slightly as the ripple amplitude is increased in the perturbed equilibria. The $q$ profile, which is the same for all these configurations, is shown in figure 2. This $q$ profile is not standard for a tokamak configuration. Actually, it corresponds to a TJ-II configuration with high rotational transform. The collisionless Landau damping of the GAM oscillation strongly depends on the safety factor [33]. Therefore, as advanced above, we use unusually small safety factor values to reduce the amplitude of the GAM oscillation, which allows to observe the low frequency oscillation. For a $q$ profile typical from tokamaks, the GAM is weakly damped even at the center where $q$ is small, and the low frequency oscillation in which we are interested is undetectable even for large values of ripple added.

In this set of configurations we calculate the zonal flow oscillation frequency with CAS3D-K and also by means of EUTERPE simulations. We use flat density and temperature profiles with $T_{e}=T_{i}=5$ keV. In EUTERPE we use an initial perturbation with $k_{\psi}=0.5\pi$, so that the normalized radial scale of the perturbation is $\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}<0.026$.

The radial dependence of the zonal-flow frequency calculated with CAS3D-K and EUTERPE for the tokamaks described above is shown in figure 3 for different ripple values. The oscillation frequency shows rather flat radial profiles with a slight increase near the magnetic axis. The calculations with both codes show good agreement and this is better in configurations with larger ripple amplitudes; the reason is given next. The time traces of the normalized zonal electric field corresponding to the cases shown in figure 3 are provided in figure 4 at a radial position of $r/a=0.5$. As can be seen in this figure, the initial time steps are dominated by a fast oscillation followed by a decaying smaller-frequency oscillation which has larger amplitude in the tokamak configurations with larger ripple values. It is then clear from figure 4 that the estimation of the frequency with the fitting method is less precise in configurations with smaller ripple values because the amplitude of the zonal-flow oscillation in those configurations is actually small and it is superimposed to a slow decay to the residual level. In particular, the fitting method fails for sufficiently small $\xi$. Hence, the slight differences in the calculation of $\Omega_{\rm ZF}$ from the semianalytical approach and from the gyrokinetic simulations seem to be simply due to the uncertainties in the fitting method.

With CAS3D-K we can compute $\Omega_{\rm ZF}$ for arbitrarily small values of the ripple. Theoretically, a linear dependence of $\Omega_{\rm ZF}$ on $\xi$ is expected for $\xi\ll 1$. Note that the results of [34] are applicable here and imply $\overline{\omega_{i}}\sim\xi\rho_{ti}v_{ti}/L$ for $\xi\ll 1$. Since $A_{0}$ and $A_{1}$ (see (32) and (33)) are dominated by the axisymmetric terms, the expansion of (36) for $\xi\ll 1$ yields

[TABLE]

This is confirmed by figure 5.

The series of perturbed tokamak equilibria used in this section have allowed us to clarify the dependence of the oscillation frequency with the ripple and understand some difficulties that appear in the process of comparing the two approaches that are the subject of this paper. In particular, the calculation of the oscillation frequency by means of gyrokinetic simulations showed to have problems when the ripple size is too small because the fit of the potential time traces becomes complicated. The calculation with CAS3D-K, however, does not have this limitation and proved to be robust even with very small ripple amplitudes.

5 Zonal-flow oscillation frequency in stellarators

In this section, we calculate the zonal-flow frequency in stellarator geometries with CAS3D-K and compare the results with those obtained from gyrokinetic simulations using the codes EUTERPE and GENE. Specifically, we work out the frequency in several magnetic configurations of the W7-X, TJ-II and LHD stellarators.

As discussed in previous sections, the time evolution of an initial zonal-flow perturbation shows, in general, an initial damped GAM oscillation followed by a decaying low frequency oscillation. The safety factor profiles for all the stellarator configurations studied in this section are shown in figure 6. An example of the time traces of the zonal electric field obtained with EUTERPE is given in figure 7 for some of the considered stellarator configurations. As can be seen in this figure, the evolution of the initial zonal-flow perturbation is qualitatively different in each device. The standard and high mirror configurations of W7-X show a very small amplitude GAM oscillation before $t=5R/v_{ti}$ and a much slower and higher-amplitude oscillation afterwards. The curve corresponding to LHD shows up to six clear cycles of the GAM oscillation, while the low frequency one is almost imperceptible. TJ-II represents an intermediate situation: the amplitude of the GAM oscillation is larger than in W7-X but smaller than in LHD. As for the low frequency oscillation, two oscillation cycles are perfectly observable. Therefore, obtaining an accurate value of the zonal-flow frequency from the fitting method requires a careful analysis case by case.

5.1 W7-X stellarator: standard configuration

In this subsection we calculate the oscillation frequency in the standard configuration of the W7-X stellarator (W7-X SD). In the W7-X SD configuration the average magnetic field strength at the magnetic axis is $B_{0}=2.42$ T. The $q$ profile of this configuration is given in figure 6. In the calculations presented here we use flat density and temperature profiles with $T_{e}=T_{i}=5$ keV.

The calculations of the zonal-flow frequency with CAS3D-K, EUTERPE and GENE in the W7-X SD configuration are shown in figure 8. The oscillation frequency in all cases shows a rather flat radial profile and the agreement between them is very good. In the simulations carried out with EUTERPE we use an initial perturbation with $k_{\psi}=0.5\pi$, so that the normalized radial scale of the perturbation varies radially, with $\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}<0.032$. Very similar values of the oscillation frequency were obtained in simulations with different (small) values of $k_{\psi}$ (not shown here). Analogously, several simulations were carried out with GENE for different radial scales of the perturbation with $0.01\leq\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}\leq 0.05$ obtaining very similar values of the oscillation frequency (see figure 8). These results are in agreement with equation (36), which is independent of $k_{\perp}$ in the long-wavelength limit.

In the EUTERPE calculations, the error when obtaining the frequency is larger in the outer region of the plasma because the fit to the model function gives a less precise value. The treatment of lost particles at the outer boundary can also introduce some bias in the oscillation frequency. In the EUTERPE calculations an error on the order of a $20\%$ of its value can be assumed in these radial positions. It is in these radial positions where the CAS3D-K calculation shows a significant increase in the frequency and the differences with the EUTERPE calculations are on the order of the error.

Finally, it might be useful to provide values for the frequency in physical units. Going back to figure 8, recalling that in this case $T_{i}=5$ keV and employing that the major radius of W7-X is $R=5.5$ m, we find $\Omega_{{\rm ZF}}/(2\pi)\sim 2.5$ kHz.

5.2 TJ-II stellarator

We have calculated the zonal flow oscillation frequency in the standard configuration of the TJ-II stellarator with CAS3D-K and EUTERPE. The $q$ profile of this configuration is given in figure 6. We also use flat density and temperature profiles, in this case with temperatures $T_{e}=T_{i}=100$ eV, close to typical ion temperature values in ECRH plasmas of TJ-II. In EUTERPE we use, as in the previous case, an initial perturbation with $k_{\psi}=0.5\pi$, so that the normalized radial scale of the perturbation is $\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}<0.023$.

The results of the zonal-flow frequency calculated with CAS3D-K and EUTERPE for this configuration are shown in figure 9. As can be seen in this figure, the zonal-flow frequency in the TJ-II stellarator increases with the radial coordinate $r/a$ and the calculations with the different numerical methods show a good agreement.

It is important to point out that in this configuration the calculation of $\Omega_{\rm ZF}$ from the fit to the model function (41) is not as accurate as in the W7-X SD case because the damping of the low frequency oscillation is larger in TJ-II and fewer oscillation cycles are observed (see figure 7). However it can be calculated with reasonable accuracy from the EUTERPE electric field time traces for all radial positions because the initial GAM oscillation is quickly damped leaving a couple of cycles of neat zonal-flow low frequency oscillation. The first part of the time trace, in which the GAM oscillation is present, is not included in the fit.

From figure 9 we can deduce that in TJ-II, whose major radius is $R=1.5$ m, the frequency varies radially in the range $\Omega_{{\rm ZF}}/(2\pi)\sim 5-8$ kHz for $T_{i}=100$ eV.

5.3 LHD stellarator

In the LHD stellarator, we calculate the zonal-flow frequency in three magnetic configurations with different positions of the magnetic axis, $R_{\mathrm{ax}}$, and different values of the magnetic field strength at the magnetic axis, $B_{0}$. These are: the standard configuration (LHD SD) with $R_{\mathrm{ax}}=3.74$ m and $B_{0}=2.53$ T; an outward-shifted configuration (LHD OS) with $R_{\mathrm{ax}}=3.91$ m and $B_{0}=1.48$ T; and an inward-shifted configuration (LHD IS) with $R_{\mathrm{ax}}=3.57$ m and $B_{0}=1.57$ T. The $q$ profiles of these configurations are given in figure 6 and we take flat density and temperature profiles with $T_{i}=T_{e}=5$ keV. In EUTERPE simulations we use an initial perturbation with $k_{\psi}=0.5\pi$, so that the normalized radial scale of the perturbation is $\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}<0.018$.

The frequencies calculated with CAS3D-K and EUTERPE in the LHD SD configuration are shown in figure 10. Both calculations show a clear radial increase of the frequency. In this case, the agreement between calculations with CAS3D-K and EUTERPE cannot be expected to be optimal, as can be easily understood from a simple inspection of figure 7. This figure shows that the GAM oscillation is weakly damped and it is almost impossible to distinguish the low frequency oscillation. Therefore, in this configuration the fitting process turns out to be very complicated. In this sense, the agreement shown in figure 10 is actually remarkable.

In figure 11 we show the CAS3D-K calculations of the zonal-flow frequency in the LHD SD, LHD OS and LHD IS configurations. They also exhibit a radial increase of the frequency in all configurations, except for the LHD IS configuration in external radial positions, for $r/a>0.9$. The inward-shifted configuration is better optimized for neoclassical transport, which has lower values of the averaged magnetic drift frequency, $\overline{\omega_{s}}$. This, at the same time, shows a larger value of the zonal-flow residual level (see reference [29]). As can be seen in figure 11, the frequency of the zonal-flow oscillation in the LHD IS configuration is smaller than in the LHD SD or LHD OS configurations, which is considered also a direct consequence of the reduction in $\overline{\omega_{s}}$.

As we did for the stellarator configurations discussed in previous subsections, we give values of the oscillation frequency in physical units for the LHD configurations studied here. Recall that $T_{i}=5$ keV. In the standard and outward shifted configurations the frequency increases radially in the range $\Omega_{\rm ZF}/(2\pi)\sim 1.85-18.5$ kHz, while in the inward shifted configuration it lies in the range $\Omega_{\rm ZF}/(2\pi)\sim 1.9-8.5$ kHz according to CAS3D-K calculations.

We do not show the calculation of $\Omega_{\rm ZF}$ with EUTERPE in the IS and OS configurations of LHD because a fit to a model like (41) is not reliable. Basically, under the conditions chosen, the low-frequency oscillation is undetectable in the gyrokinetic simulations in these two configurations. The usefulness of calculating the zonal-flow oscillation frequency with CAS3D-K in the IS and OS configurations of LHD could be questioned in view that it is difficult to observe this oscillation in practice. However, this low frequency oscillation constitutes a natural mode of oscillation that could manifest under particular conditions that enhance it and/or diminish the GAM oscillation. Oscillations in this frequency range could, in principle, be excited by some forcing mechanism, such as fast ions or electrons, and might be measured. With this idea in mind, the calculations shown in figure 11 could be useful even in this kind of configurations.

5.4 W7-X stellarator: high-mirror configuration

In figure 12 we show the values of $\Omega_{\rm ZF}$ for the high-mirror configuration of W7-X (W7-X HM) obtained with CAS3D-K (red line) and EUTERPE (blue line). In the calculations presented here we use, as in the W7-X SD case, flat density and temperature profiles with $T_{e}=T_{i}=5$ keV. In the W7-X HM configuration the average magnetic field strength at the magnetic axis is $B_{0}=2.35$ T. The $q$ profile is given in figure 6. In the EUTERPE simulations we use an initial perturbation with $k_{\psi}=0.5\pi$, so that the normalized radial scale of the perturbation varies radially, with $\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}<0.032$. Note that in the W7-X high-mirror configuration, and for the indicated temperatures, the values of the frequency in physical units are on the order of $\Omega_{\rm ZF}/(2\pi)\sim 3$ kHz, which are slightly larger than in the standard configuration.

The CAS3D-K and EUTERPE calculations are in very good agreement for $r/a<0.6$, but not for $r/a>0.6$. In this case the difference between calculations of CAS3D-K and EUTERPE seem too large to be attributed to fitting errors or boundary condition effects in the gyrokinetic simulations. Let us try to understand the disagreement.

Note that in equation (10) the component of the magnetic drift tangent to the flux surface has been neglected. This is formally correct as long as condition (14) is satisfied. We will see next, however, that in the W7-X HM configuration the neglected term contributes to the calculation of $\Omega_{\rm ZF}$. We cannot easily extend the analytical calculation of $\Omega_{\rm ZF}$ to include the effect of the tangential component of the magnetic drift, but we can carry out the EUTERPE simulations removing this component, and therefore leaving only the radial component of the magnetic drift, which is a more faithful comparison to (10).

The black curve in figure 12 corresponds to the EUTERPE simulation retaining only the radial component of the magnetic drift. In these simulations the damping of the low frequency oscillations increases significantly with respect to the cases in which the full magnetic drift is kept. This makes the fit to a damped oscillation model like (41) more difficult, but still possible. We can safely conclude that the role of the tangential magnetic drift in the W7-X HM configuration is needed to explain the radial dependence of $\Omega_{\rm ZF}$ for $r/a>0.6$.

The relevance of the tangential component of the magnetic drift in the W7-X HM configuration, by comparison to its irrelevance in the W7-X SD one, can be understood by looking at figure 13, where the flux-surface-averaged absolute value of the normal and geodesic components of the field line curvatures is shown for both configurations. While the geodesic component of the field line curvature is very similar in both cases, the normal component (related to the tangential component of the magnetic drift) is larger in the W7-X HM configuration than in the W7-X SD one. Actually, the difference is significant for $r/a\gtrsim 0.5$ and gets larger at outer positions (compare figure 13 with the black and blue curves in figure 12).

At this point, one might wonder why we do not do the same test for the TJ-II and LHD SD configurations. That is, one might think that perhaps it is possible to prove that the agreement between CAS3D-K and EUTERPE simulations cannot be improved for TJ-II and LHD SD because the tangential component of the magnetic drift also counts in these configurations. However, this check turns out to be not feasible. In the simulations carried out neglecting the tangential component of the magnetic drift, the damping of the oscillations is remarkably higher with respect to the cases in which the full magnetic drift is retained. As we have already explained, in W7-X HM this effect is not enough to prevent a reasonable fit. On the contrary, in the rest of configurations studied (W7-X SD, TJ-II and LHD) the damping is so large when the tangential magnetic drift is removed that the fit is not reliable.

6 Simulation details and computation time

One of the most important features of the semianalytical method herein proposed for the evaluation of the zonal flow oscillation frequency is that it is faster than calculating it by means of gyrokinetic simulations with EUTERPE or GENE. We show in table 1 the total CPU time required to obtain the frequency with each method. The times shown in the table for CAS3D-K and GENE correspond to the time required for calculations at just one radial position, while for EUTERPE, as it is a global code, these values are given for the simulation of all radial positions at the same time. The values shown in table 1 are determined by the resolution used in the different cases with CAS3D-K and by the simulation time required to obtain a large enough number of oscillation cycles to make a fit in the case of EUTERPE and GENE. In the EUTERPE case, as they are global simulations, the minimum time required is determined by the most demanding radial position (the innermost radial locations, in general). These times represent the requirements for simulations or calculations converged with respect to the different numerical parameters and resolutions in each code. The numerical details in each case are listed below.

In CAS3D-K, the computation time to obtain converged results varies among devices and it depends on the phase-space resolution used in each case. The minimum resolution for the integration over the $\lambda$ coordinate in all the calculations is $n_{\lambda}=24$. The integration over the magnitude of the velocity is performed analytically. For the spatial integrals we use $n_{\zeta}=128$ trajectories and $n_{\theta}=512$ integration points per trajectory, in the case of passing particles. For trapped particles, these numbers are $n_{\zeta}=16$ trajectories per integration group and $n_{\theta}=512$ integration points per trajectory. The calculations in CAS3D-K are radially local. Therefore, the values in the table corresponding to CAS3D-K calculations are given per radial position and the ranges given in table 1 account for the maximum and minimum CPU times among different radial positions at a given configuration.

The EUTERPE simulations presented here are not extremely demanding from a computational point of view as only the zonal ($m=0$, $n=0$ Fourier component) and several smaller amplitude sidebands are resolved. The simulations were carried out with the following numerical parameters. In all the cases the resolution in poloidal and toroidal angles in PEST coordinates was $n_{\theta}=32$, $n_{\phi}=32$. The radial resolution was $n_{s}=24$ for the rippled tokamaks, TJ-II and LHD, while it was $n_{s}=32$ for the W7-X configurations. Simulations in W7-X were carried out with more radial resolution because several simulations with different radial scales of the perturbation, and maintaining the numerical parameters, were carried out for comparison. The number of markers used was $n_{M}=40\mathrm{M}-50\mathrm{M}$.

In general, the EUTERPE calculations require larger computational resources for devices with lower zonal-flow frequencies as more simulation time is required to resolve a number of oscillation cycles that is large enough to make a fit. For the rippled tokamaks discussed in Section 4, the required computational time with EUTERPE increases by a factor of 5 for the case with smaller ripple with respect to that with the largest one because the oscillation frequency decreases by this factor, while all physical and numerical parameters in the simulation are the same. In the case of LHD configurations there is a large amplitude GAM oscillation, which makes difficult the fit of the potential time trace to a model like (41) for the low frequency oscillation, which has much smaller amplitude. A longer simulation time was required as compared to other configurations to obtain a reasonable fit in the LHD standard configuration. Only in the standard configuration of LHD a reliable fit of the potential time traces was possible.

The GENE calculations were performed employing its full flux-surface version. Therefore, the calculations are 2D in the spatial coordinates $\{y,z\}$ and 2D in velocity space $\{v_{\parallel},\mu\}$. Here, $y$ is the coordinate along the binormal direction, $z$ is the coordinate along the field line, $v_{\parallel}$ is the parallel velocity and $\mu$ is the magnetic moment coordinate. The results of the CPU time per radial position obtained with GENE in the W7-X SD configuration are given in table 1. An analysis of convergence in velocity space resolution was carried out and the results are shown in figure 14. In that figure, several electrostatic potential time traces at $r/a=0.5$ are plotted, obtained from simulations with GENE in the W7-X SD configuration with $\left\langle k_{\perp}\rho_{ti}\right\rangle_{\psi}=0.05$ and using different resolutions in velocity space. The CPU time given for GENE in table 1 corresponds to the simulation of figure 14 with $n_{v_{\parallel}}=128$ and $n_{\mu}=8$. All the simulations were carried out with a very small value of hyperdiffusivity (two orders of magnitude smaller than typical values used in turbulence simulations), which is required to reduce the damping of the zonal flow oscillation and allow a reliable fit [35].

The calculations/simulations in this work were carried out in different supercomputers, with different capabilities, so that the computing time in the table 1 has to be considered as indicative. The EUTERPE simulations were carried out in two different supercomputers, EULER and MareNostrum III [36]. EULER is equipped with Intel Xeon 5450 quadcore processors at 3.0 GHz and Infiniband 4X DDR and Mare Nostrum III is equipped with Intel SandyBridge-EP processors at 2.6 GHz and Infiniband FDR10 interconnection. From 32 to 64 computing cores were used in the simulations. All the CAS3D-K calculations shown in this work were run also in the EULER supercomputer. The GENE calculations were carried out in the Uranus supercomputer, which is equipped with Intel Xeon E5-2630 processors at 2.4 GHz interconnected by Infiniband FDR. The computing time shown in the table is always the total CPU time (summed for all the computing cores used).

From the numbers shown in table 1, even corresponding to different computing facilities, it is clear that the calculations of the zonal-flow frequency with CAS3D-K are faster than those employing the gyrokinetic codes EUTERPE and GENE. This result strongly supports the usage of CAS3D-K in this type of computations.

7 Discussion and conclusions

In this work, we have proven the efficiency of a semianalytical method for calculating the zonal flow oscillation frequency in stellarators and rippled tokamaks. It is based on the numerical evaluation of expression (36), which was first derived in [11] but had not been compared with the frequency obtained from direct gyrokinetic simulations so far.

We have extended the code CAS3D-K for the evaluation of expression (36) in general rippled tokamak and stellarator configurations. In particular, we have implemented in CAS3D-K, for both passing and trapped particles, the correct solution for $\delta_{s}$, the displacement from the initial flux surface at each point of the particle orbit. The accuracy of this semianalytical approach using the code CAS3D-K was checked by comparing its results against the frequency obtained from gyrokinetic simulations with the global code EUTERPE and the radially local code GENE in a wide range of configurations. Specifically, we have calculated in a series of large aspect ratio tokamaks with different ripple values, as well as in the W7-X, TJ-II and LHD stellarators.

When using gyrokinetic codes, we obtain the zonal-flow frequency by fitting the time trace of the normalized zonal electrostatic potential (or electric field) to a model function, given in (41), including a damped oscillation and the decay to a residual value. The evolution of an initial zonal perturbation is initially dominated by a decaying GAM-like oscillation followed by a damped low frequency oscillation. Therefore, obtaining a precise value of the zonal-flow frequency from the fitting method is non-trivial, in general, and requires a case by case discussion.

The good agreement between the two methods (the semianalytical approach and the calculation of the frequency from gyrokinetic simulations) in real stellarator configurations supports the validity of the approximations in the derivation of (36) and the accuracy of the semianalytical approach via the extension of the CAS3D-K code. Only in the outer radial region of the high-mirror configuration of W7-X the accuracy of the semianalytical method is clearly insufficient. The reason is that, as we have explained, in this case, the tangential component of the magnetic drift, that is dropped in the analytical calculation leading to (36), must be kept.

The advantage of using CAS3D-K is that the computation time can be reduced up to two orders of magnitude with respect to the gyrokinetic calculations. This makes this method an option to be included in an stellarator optimization loop, in which CAS3D-K could provide fast calculations of zonal flow relaxation properties (oscillation frequency and residual level [14]) to be used as figures of merit of stellarator configurations.

Finally, it is worth emphasizing that expression (36) captures the influence of the magnetic geometry on the zonal flow oscillation frequency. However, as explained in subsection 2.1, this value of the frequency can be modified in the presence of a background radial electric field, $E_{\psi}$. The rigorous calculation of the corrections introduced by $E_{\psi}$ to expression (36) are beyond the scope of this work.

Acknowledgments

P. M. thanks Per Helander and Tobias Görler for helpful discussions. The authors thank Antonio López-Fraguas for his help with the usage of VMEC and acknowledge the computer resources, technical expertise and assistance provided by the Barcelona Supercomputing Center (BSC) and the Computing Center of CIEMAT. GENE runs have been carried out in Uranus, a supercomputer cluster located at Universidad Carlos III de Madrid (Spain) funded jointly by EU FEDER funds and by the Spanish Government via the National Projects UNC313-4E-2361, ENE2009-12213-C03-03, ENE2012-33219 and ENE2012-31753.

This research has been funded in part by grants ENE2012-30832 and ENE2015-70142-P, Ministerio de Economía y Competitividad (Spain) and by an FPI-CIEMAT PhD fellowship. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.