Comparison of local and global gyrokinetic calculations of collisionless zonal flow damping in quasi-symmetric stellarators

The linear collisionless damping of zonal flows is calculated for quasi-symmetric stellarator equilibria in flux-tube, flux-surface, and full-volume geometry. Equilibria are studied from the quasi-helical symmetry configuration of the Helically Symmetric eXperiment (HSX), a broken symmetry configuration of HSX, and the quasi-axial symmetry geometry of the National Compact Stellarator eXperiment (NCSX). Zonal flow oscillations and long-time damping affect the zonal flow evolution, and the zonal flow residual goes to zero for small radial wavenumber. The oscillation frequency and damping rate depend on the bounce-averaged radial particle drift in accordance with theory. While each flux tube on a flux surface is unique, several different flux tubes in HSX or NCSX can reproduce the zonal flow damping from a flux-surface calculation given an adequate parallel extent. The flux-surface or flux-tube calculations can accurately reproduce the full-volume long-time residual for moderate $k_x$, but the oscillation and damping time scales are longer in local representations, particularly for small $k_x$ approaching the system size.


I. INTRODUCTION
Control of turbulent transport is a crucial step in the development of fusion energy.
In many cases, zonal flows can play a role in the regulation and reduction of turbulent transport. A zonal flow is a toroidally and poloidally symmetric E × B flow that can be driven by electric fields that develop from fluctuations of the plasma potential, such as is the case for most drift wave turbulence [1]. The zonal flow does not drive transport itself, but by facilitating transfer of energy between radial wavenumbers it can regulate the linear instability and affect turbulence saturation [2,3]. Strong zonal flows have been found to be important in configurations such as tokamaks [4] or the reversed-field pinch [5][6][7], and they can affect turbulence saturation in stellarators as well. [8][9][10][11] Linear zonal flow damping is often examined as a proxy for the full zonal flow evolution [12] and is used in models to predict turbulent transport [13,14]. The Rosenbluth-Hinton model [15] provides a zonal flow residual that describes the undamped part of the poloidal flow in a large-aspect-ratio tokamak. This undamped flow acts to saturate drift wave turbulence, and the residual is used as a measure of the amplitude that zonal flows achieve in nonlinear simulations. In axisymmetric systems, this is commonly the case, and the residual is sometimes used as a proxy for the resulting turbulence saturation [16]. However, this is unlikely to be true if the collisionless damping to the residual is slow compared to the rate at which turbulence injects energy into the zonal flow. In non-axisymmetric devices, the radial drift of trapped particles can drive long-time damping and oscillations of the zonal flow [12,[17][18][19][20], as will be discussed in Sec. II. These features can disassociate the zonal flow residual from saturated turbulence. Calculations in this paper are linear and do not address the transfer of energy between modes, but can examine changes in the collisionless damping of the zonal flow.
Zonal flow damping and the driving turbulence both depend on aspects of the magnetic geometry, such as the rotational transform and trapped particle regions. Due to the large number of parameters that can describe the plasma boundary, the 3D shaping of stellarators offers a large parameter space to search for configurations that can benefit specific turbulence or zonal flow properties. Particularly in helical systems optimized to reduce neoclassical transport, stronger zonal flows may reduce turbulent transport [21]. However, nonlinear simulations are expensive to include in an iterative optimization loop, and an efficient, general, linear proxy for turbulent transport would be a powerful tool. In order to obtain such a proxy, a thorough understanding of zonal flow dynamics in stellarators is required.
Zonal flows have been studied numerically in flux-tube geometry for the Large Helical Device [22] and Wendelstein 7-X [12], and in full-volume geometry for TJ-II [19,23], the Large Helical Device [20], and Wendelstein 7-X [20,24]. As part of benchmarking gyrokinetic codes, full-volume linear calculations of zonal flow damping have been compared to local analytic theory [20,[24][25][26]. However, quasi-symmetric configurations are absent from previous studies, despite the expectation that a perfectly quasi-symmetric configuration will support an undamped zonal flow similar to a tokamak. A quasi-symmetric stellarator has a symmetry in the magnitude of the magnetic field |B|, and the magnetic spectrum is dominated by a single mode B mn , so that B ≈ B mn cos(nφ − mθ). (1) Here, n and m are the toroidal and poloidal mode numbers, and φ and θ are the toroidal and poloidal coordinates, respectively. When the magnetic field is described by a single mode, the collisionless bounce-averaged drift of trapped particles from a flux surface goes to zero, reducing neoclassical transport and flow damping. Different quasi-symmetries are defined by the choice of dominant mode in the magnetic spectrum. Quasi-polidal symmetry has a dominant m = 0 mode, and is not included here. Quasi-axial symmetry has a dominant n = 0 mode, similar to a tokamak, as seen in Fig. 1. Quasi-helical symmetry (QHS) uses a single n = 0, m = 0 mode, creating the helical shape of the |B| contours in Fig. 2.
In this paper, the zonal flow damping is numerically calculated in flux-tube, flux-surface, and full-volume geometry representations for quasi-symmetric configurations. We look to understand how much geometry information is required for an accurate determination of the zonal flow time evolution. Although neoclassical transport and flow damping in quasisymmetric stellarators is more similar to tokamaks than to classical stellarators, we show that the linear zonal flow response for a realistic but almost quasi-symmetric geometry still resembles a classical stellarator.
The paper is structured as follows. Sec. II reviews collisionless zonal flow damping in non-axisymmetric equilibria and introduces the geometries and numerical tools used in this work. Sec. III identifies the differences in calculations of zonal flow damping in full-volume, flux-surface, and flux-tube frameworks. In Sec. IV, calculations of zonal flow damping in flux tubes are shown to reproduce the zonal flow residual from flux-surface calculations, but only for sufficiently long flux tubes. Sec. V presents results from the quasi-symmetric and broken-symmetry configurations of HSX and compares them to the NCSX zonal flow evolution.

II. COLLISIONLESS ZONAL FLOW DAMPING
The Rosenbluth-Hinton model [15] quantifies the long-time linear response of the zonal flow to a large-radial-scale potential perturbation in a collisionless, axisymmetric system.
The initial amplitude of the perturbation is reduced by plasma polarization and undergoes geodesic acoustic mode (GAM) oscillations before relaxing to a steady-state residual. The long-time residual zonal flow is defined as the ratio of the zonal potential in the long-time limit to the initial zonal potential. In a large-aspect-ratio tokamak, the residual amplitude depends on the geometry as [15], and can be interpreted as a measure of how strongly collisionless processes modify the zonal flow. Here, ϕ is the zonal potential, ϕ 0 is its initial amplitude, q is the safety factor, and ǫ t is the inverse aspect ratio of the flux surface of interest. The term 1.6q 2 /ǫ 1/2 t results from the neoclassical polarization due to toroidally trapped ions. When the Rosenbluth-Hinton residual is high, collisionless zonal flow damping is small and the system can support strong zonal flows. When the residual is small, damping is significant, and the existence of strong zonal flows will depend on strong pumping from the turbulence.
The zonal flow response in non-axisymmetric systems is significantly modified by neoclassical effects. The zonal flow amplitude after the polarization decay is no longer the Rosenbluth-Hinton residual. Instead, helical systems exhibit decay described by a timescale τ c ∼ 1/|k xvdr | to reach a residual in a long-time limit [27,28]. Here, k x is a radial wavenumber, andv dr is the bounce-averaged radial drift velocity. In an unoptimized device,v dr is large and the zonal flow will decay quickly to a residual, whereas a well-optimized device will have very long decay times. In a perfectly symmetric device, no long-time decay is observed, corresponding to the limit of infinitely long decay times. In this case, any geometry with finite radial particle drifts will decay to zero residual as k x → 0. Defining the residual as the zonal potential at some time shorter than t → ∞ necessarily involves neoclassical effects, as discussed in Sec. V. The long-time decay in helical systems could prevent any connection between the zonal flow residual and saturated turbulence. If a system takes a long time to decay, nonlinear energy transfer will become important before the decay has dissipated energy in the zonal mode, and the residual no longer relates to the zonal flow amplitude in the quasi-stationary state.
Furthermore, an oscillation in the zonal flow is caused by neoclassical effects [17,18].
This oscillation is characterized by the radial drift of trapped particles, as opposed to the passing particle dependence of the GAM. Drifting trapped particles cause a radial current that interacts with the zonal potential perturbation to cause zonal flow oscillations. This oscillation is damped by Landau damping on trapped particles. The oscillation damping and frequency both increase with the radial particle drift, or equivalently, neoclassical transport.
For an unoptimized device, the zonal flow oscillation is of higher frequency but damped more quickly. In a well-optimized device, the zonal flow oscillation is prominent due to the small damping, but the oscillation frequency is small compared to characteristic inverse time scales in fully developed turbulence. In a perfectly symmetric device, the zonal flow oscillation vanishes.
For both stellarators and tokamaks, the zonal flow residual depends on the radial wavenumber of the zonal perturbation, but this dependency is stronger in stellarators than in tokamaks [24,[29][30][31]. In this paper, radial wavenumbers are normalized as k x ρ s , where ρ s is the ion sound gyroradius. The numerical calculation of the zonal flow residual in a tokamak matches the Rosenbluth-Hinton residual as k x approaches zero. However, any geometry with finite radial particle drifts causes the zonal flow residual to vanish as k x → 0, although the long-time decay to that residual can be very slow in a well-optimized device.
Zonal flow oscillations, zonal flow damping, and even the zonal flow residual are further modified by the inclusion of a radial electric field [32][33][34]. The radial electric field drives coupling across field lines in the poloidal direction, and it is likely this would be visible in the difference between flux-tube and flux-surface calculations. Radial electric fields are not included here, and their effect on calculations in quasi-symmetric devices, or in local and global geometry representations, is left for future work.

A. Simulations in local and global geometry representations
A zonal flow is a toroidally and poloidally symmetric potential perturbation, and the local geometry anywhere on the surface can potentially be important to determine its response. In an axisymmetric geometry, a field line followed for one poloidal transit samples all unique magnetic geometry on a surface, as would any other field line on the same surface.
However, different field lines in a stellarator do not generally sample the same geometry.
Local geometry variations that may be important for the zonal flow may not be sampled by a is the minor radius, and c s is the ion sound speed.

Full-volume geometry
Full-volume calculations of zonal flow damping provide the most complete representation of geometry effects on the zonal flow. In this work, these calculations are carried out with the δf gyrokinetic particle-in-cell code EUTERPE [35,36]. The details of the zonal flow calculation are discussed in Refs. 20 and 24. The full-volume geometry of the fields is represented in real space using the PEST magnetic coordinates [37], where φ is the toroidal angle, s is the toroidal normalized flux, and θ * is the poloidal angle defined such that field x ρ 2 s and I 0 is the modified Bessel function [20]. Then the quasi-neutrality equation for adiabatic electrons reads: where n 0 is the equilibrium density, B is the magnetic field, n i is the gyroaveraged ion The flux tube is a reduced-geometry representation for toroidal magnetic geometries, and is constructed from a sheared box around a single field line identified by a field line label in PEST [37] coordinates as α = √ s 0 /q 0 (qθ * − φ). When the flux tube is centered on the outboard midplane, α is also a toroidal coordinate of the center point of the flux tube. The box uses non-orthogonal coordinates x in the radial direction, y in the flux surface, and z along the field line. In a Gene flux tube, a spectral representation is used for the x and y directions, while the z direction is in real space. For k y = 0 zonal modes, boundary conditions in x, y, and z are periodic. In axisymmetric systems, a flux tube of one poloidal turn samples all unique geometry on the flux surface. In a non-axisymmetric system, a flux tube does not necessarily close upon itself. The standard approach to using flux tubes in stellarators does not require true geometric periodicity of the flux tube [40]. A stellarator-symmetric flux tube, or one that is symmetric about the midpoint z = 0, provides continuous, but not necessarily smooth, geometry at the flux tube endpoints. However, k z = 0 modes, such as zonal flows, may be sensitive to the geometry at this boundary. True geometric periodicity requires that qn pol N is an integer, where N is the toroidal periodicity of the geometry. We The zonal flow damping, and resulting residual and oscillations, is calculated by initializing a flux-surface-symmetric impulse to the distribution function at a single radial mode and allowing the amplitude to collisionlessly decay due to classical and neoclassical polarization without further energy input. In Gene, the zonal perturbation is implemented by initializing only one k x = 0, k y = 0 mode. In a Fourier representation, the wavenumber is explicit, and a long-wavelength approximation is not used. The perturbation is introduced in the density with Maxwellian velocity space, which produces an equivalent potential perturbation to that used in EUTERPE. Note that for the present case of adiabatic electrons, the two initial conditions discussed in Ref. [41] are identical.
Numerical calculations for linearly stable systems without dissipation may have to contend with numerical recurrence phenomena. Such recurrence, which results from a reestablishment of phases from the initial condition and concomitant unphysical temporary increase in amplitudes, can be eliminated by including numerical spatial or velocity hyper-diffusion [42]. However, numerical diffusion is not an appropriate solution in nearly quasi-symmetric stellarators, as calculation times are very long and even a small amount of diffusion will cause significant damping of the zonal-flow residual. To solve the problem, the parallel velocity space grid spacing ∆v can be decreased sufficiently that the recurrence time, τ rec = 2π/(k z ∆v ), exceeds the duration of the simulation [43]. In the present work, most flux- leading to ∆v = 0.015 v T i . Here, v T i = (2T i /m i ) 1/2 is the ion thermal velocity. We take k z to be the wavenumber of the periodicity of the magnetic structure, as seen in

B. The HSX and NCSX geometries
The zonal flow response is studied by means of flux-tube, flux-surface, and full-volume gyrokinetic simulations in the Helically Symmetric eXperiment (HSX) [44] and the National Compact Stellarator eXperiment (NCSX) [45]. VMEC [46] is used to calculate the HSX and NCSX equilibria.
HSX is a four-field-period optimized stellarator, designed to improve single-particle confinement through quasi-helical symmetry. The main coils produce the helically symmetric field, and a set of auxiliary coils can be energized to modify the magnetic spectrum. The symmetry is broken by adding mirror terms to the spectrum, in which case transport is similar to a classical stellarator with large neoclassical transport in the low collisionality regime. HSX has demonstrated reduced neoclassical flow damping [47] and transport [48] in the QHS configuration. It has also been hypothesized that neoclassical optimization may reduce anomalous transport through stronger zonal flows [21]. Zonal flows are clearly ob- Colors correspond to |B|, where blue is the minimum field strength. A flux tube of one poloidal turn is shown in red, and one of 4 poloidal turns in green.
-4 π -3 π -2 π -1 π 0 π 1 π 2 π 3 π 4 π θ 0 |B|, HSX flux tubes -8 π -6 π -4 π -2 π 0 π 2 π 4 π 6 π 8 π θ 0.8 The NCSX configuration was also optimized for neoclassical transport, but is a threeperiod device designed for quasi-axial symmetry. The equilibrium used here has total normalized plasma pressure β ≈ 4%. In the NCSX configuration, we use the flux tubes symmetric about the midpoint z = 0 (α = 0 and α = π/N, with N = 3 for NCSX), as well as one non-symmetric flux tube (α = π/2). The radial particle drift for the surface at s = 0.5 is between that of the QHS and Mirror configurations of HSX. The rotational transform in NCSX is roughly half that in HSX, as seen in Fig. 4. The difference in rotational transform means that the part of the surface sampled by a flux tube is very different. In Fig. 2, the multiple turns of a flux tube in HSX cluster together in a band around the device. In Fig. 1 ing. Following Monreal [20], curve fitting is used to extract the residual and parameters of the zonal flow oscillation. The time evolution during the post-GAM phase is described by, The zonal flow oscillation is parameterized by an amplitude A ZF , oscillation frequency Ω ZF , and damping rate γ ZF . The long-time decay follows an algebraic decay according to Ref. [18], but is well approximated by an exponential decay Ce −γ to avoid an abundance of fitting parameters. The zonal flow residual is R ZF . The evolution of normalized zonal potential and normalized zonal electric field are equivalent for a zonal potential with only a single k x [20].
However, in practice, the zonal electric field is preferred in global EUTERPE simulations to simplify the fitting.
We are primarily interested in the zonal flow oscillation here, not the GAM. The damping of the GAM increases with decreasing rotational transform [27]. In HSX, the rotational transform is about one and the GAM is damped on time scales of the order 10 a/c s . GAM oscillations are more apparent in NCSX calculations, as the rotational transform is about twice that in HSX, but are still damped quickly compared to the zonal flow oscillations.
Fitting starts after the GAM oscillations have damped away to avoid further complexity in curve fitting. The zonal flow oscillation is Landau-damped by trapped particles, and depends on the radial drift off of the flux surface [18]. Neoclassically optimized devices can have long-lived zonal flow oscillations as neoclassical transport is reduced, which also reduces the oscillation frequency to well below the GAM frequency. In the configurations examined here, fitting the zonal flow oscillations is important to accurately fit the zonal flow residual.

III. COMPARISON OF LOCAL AND GLOBAL CALCULATIONS
The NCSX configuration is quasi-axisymmetric, which, among three-dimensional geometries, most closely resembles a tokamak. As discussed in Sec. II, the zonal flow residual as k x → 0 is a key difference between symmetric and non-symmetric systems. The time traces for full-volume, flux-surface, and flux-tube simulations are fitted to extract the zonal flow residual plotted in Fig. 5. A discussion of different flux tubes in the NCSX configuration is provided in Sec. IV B. In NCSX, the zonal flow residual goes to zero for small k x , just as it does for classical stellarators. The limit as k x → 0, as well as a peak residual around k x ρ s ≈ 0.5, is reproduced in all three geometry representations. However the amplitude of the residual differs between the local and global representations particularly for very small k x . In the full-volume calculations, the peak of the residual is slightly lower, while the smallest k x support a larger residual than the flux-tube calculations. A long-wavelength approximation valid for k x ρ s < 1 is used in the global simulations, which may be approaching its limit of validity towards the peak. Without flux surface calculations at small k x , we cannot constrain the physical cause for differences between full-volume and flux-surface results. Coupling between surfaces may be occurring, but the same disagreement is not seen for HSX configurations in Fig. 17. More importantly, the flux-surface approximation will break down when scales are large enough to involve profile effects, and the smallest k x values examined approach the machine size. Thus the observed discrepancies are to be expected given the limitations of the frameworks. flux-tube calculations, and in Fig. 7 for the full-volume and flux-tube calculations. tube is centered at α = 0, while the QHS-t "triangle" flux tube is centered at α = π/4. In  The QHS and Mirror configurations of HSX have been designed specifically to study differences in neoclassical transport and flow damping. As discussed in Sec. II, the zonal flow long-time damping and oscillation frequency are related to neoclassical transport. According to theory [18,27], the more optimized QHS configuration should exhibit lower-frequency zonal flow oscillations as well as slower long-time decay to the residual level. These expectations are verified in Fig. 16. The zonal flow oscillation frequency is higher by a factor of 2.5 in Mirror than QHS at k x ρ s = 0.02, and the long-time damping is significantly faster at k x ρ s = 0.3.
As observed in Fig. 17, the zonal flow residual does not differ between QHS and Mirror.
The Rosenbluth-Hinton residual in Eq.
(2) depends primarily on the safety factor q, a con- sequence of the ratio of the banana-induced polarization relative to the gyro-orbit-induced polarization [27]. In a non-axisymmetric system, additional trapped particles further modify the zonal flow evolution through their polarization and radial drift. While the radial drift is important for the time evolution, the zonal flow residual primarily depends upon the polarization effects provided by the trapped particles [12]. The broken symmetry of the mirror configuration increases the trapped particle radial drift, as demonstrated by the zonal flow oscillation frequency and damping, but does not change the zonal flow residual. We conclude that the helically trapped particles dominate the polarization drift and set the zonal flow residual in both systems.
As compared to NCSX, the QHS configuration produces less trapped-particle radial drift  Calculation of the zonal flow decay in HSX captures the expected neoclassical effects on decay rate and oscillation frequency. However, the saturation of drift wave turbulence is a strong motivation for the study of zonal flows. Turbulence will transfer energy and reorganize the system within a correlation time, effectively resetting the zonal flow time evolution. In nonlinear simulations of trapped electron mode (TEM) turbulence in HSX [50], the correlation time is on the order of 10 a/c s . In Fig. 18, the short-time damping of the zonal flow is plotted, but again, there is no difference between the QHS and Mirror configurations.
Depending on driving gradients, the heat flux from nonlinear TEM turbulence simulations differs between these configurations, and is not explained by the linear growth of the most unstable mode [50]. If a difference in heat fluxes between configurations is due to the linear coils [51], by magneto-hydrodynamic activity [7], or by microturbulence itself [52], erode the zonal flow residual and lead to increased turbulent transport [41,53]. However, erosion time scales in these scenarios were on the order of the turbulent correlation time, giving further credence to the idea that the long-time decay present in the systems investigated here is unlikely to affect transport directly. Future work should include external radial electric fields, which can strongly modify the zonal flow decay and residual. The radial electric field in a stellarator is usually determined by an ambipolarity constraint on neoclassical transport, which can differ between configurations but requires knowledge of density and temperature profiles.

DATA AVAILABILITY
The data supporting the findings of this study is available from the authors upon reasonable request.