Saturation of magnetised plasma turbulence by propagating zonal flows
Richard Nies Felix Parra Michael Barnes Noah Mandell William Dorland


APS/123-QED










 rnies@pppl.gov

  ^1Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA

  ^2Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA

  ^3Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3NP, United Kingdom

  ^4Department of Physics, University of Maryland, College Park, MD 20742, USA


















             

§ ABSTRACT 
Strongly driven ion-scale turbulence in tokamak plasmas is shown to be regulated by a new propagating zonal flow mode, the toroidal secondary, which is nonlinearly supported by the turbulence. The mode grows and propagates due to the combined effects of zonal flow shearing and advection by the magnetic drift. Above a threshold in the turbulence level, small-scale toroidal secondary modes become unstable and shear apart turbulent eddies, forcing the turbulence level to remain near the threshold. By including the new zonal flow physics into a theory of turbulence saturation based on the critical balance conjecture, scaling laws for the turbulent heat flux, fluctuation spectra, and zonal flow amplitude are derived and shown to be satisfied in nonlinear gyrokinetic simulations.











                              
Saturation of magnetised plasma turbulence by propagating zonal flows
    W. Dorland^4
    September 9, 2024
=====================================================================



Introduction.– Efforts to achieve the high temperatures required for thermonuclear fusion in magnetically confined plasmas are stymied by the large heat flux caused by turbulent mixing. The turbulent fluctuations are driven by micro-instabilities, most notably the ion temperature gradient (ITG) instability <cit.>. The saturation level of such ion-scale modes crucially depends on zonal flows (ZFs) <cit.>, flow bands that are nonlinearly generated by the turbulence <cit.> and can shear apart turbulent eddies. In toroidal plasmas, the linear ZF physics allows for stationary ZFs <cit.> and fast oscillating geodesic acoustic modes (GAMs) <cit.>.

In this Letter, we show that strongly driven tokamak ITG turbulence exhibits a new small-scale propagating ZF mode, the toroidal secondary. We elucidate the physics of the new mode, and demonstrate its relevance to turbulence saturation. We then show how the application of the critical balance conjecture to ITG turbulence <cit.> must be altered to account for the role of the toroidal secondary. The revised theory leads to new scalings of the fluctuation amplitude and length scales, which are found to be in agreement with nonlinear gyrokinetic simulations and with past experimental observations <cit.>. 

Tokamak turbulence.– To avoid fast parallel losses, the magnetic field of tokamaks is made to lie on nested toroidal `flux surfaces'. The plasma quickly reaches local thermodynamic equilibrium and is thus approximately Maxwellian, with the density n_s and temperature T_s of each species s uniform on flux surfaces. The resulting radial gradients between the hot dense core and the cold dilute edge drive turbulent fluctuations which cause transport across flux surfaces. These fluctuations are slow compared to the Larmor gyration frequency Ω_s = e_s B /m_s, with the particle mass m_s and charge e_s, and B the magnetic field strength. The fast gyratory motion may thus be averaged over, and the turbulence may be modelled using gyrokinetics <cit.>, which describes rings of charge centered on the gyrocentre position R, shifted from the particle position r by the gyroradius vector.

Moreover, in the tokamak core, the fluctuations in the distribution function f_s are generally small compared to the Maxwellian background F_Ms, i.e., δ f_s = f_s - F_Ms≪ F_Ms, and they are highly anisotropic. Indeed, the typical turbulence length scales parallel to the magnetic field are of the order of the tokamak size, as measured e.g. by the major radius R, while in the perpendicular direction they are on the much smaller gyroradius scale ρ_s = v_Ts/Ω_s ≪ R, with the thermal speed v_Ts=(2 T_s/m_s)^1/2. The fluctuations may thus be described locally in a flux-tube domain, as shown in Fig. <ref>, taking the magnetic field geometry to only vary along the field line on a chosen flux surface denoted by r_0. Here, we assume the tokamak flux surfaces to have circular cross-section, and we choose the radial coordinate to be x=r-r_0, with r the radial distance from the magnetic axis. As x labels flux surfaces, its gradient is parallel to the density and temperature gradients. The binormal coordinate y and the parallel coordinate θ determine the position within a flux-surface: y, like x, is perpendicular to B, whereas θ gives the location along B. Then, the flux-surface average ⟨ A ⟩_yθ of a given quantity A is given by its average in y and θ, the latter including the Jacobian 1/B·∇θ. Finally, we also define the zonal and nonzonal components of A through A^Z≡⟨ A ⟩_y = A - A^NZ.

In this Letter, we consider electrostatic and collisionless ion-scale fluctuations, for a single ion species with charge e_i = Z_i e. The gyrokinetic equation may be written in terms of the non-adiabatic part of the perturbed ion distribution function h_i = δ f_i + F_Mi e_i φ/T_i, using the particle energy and magnetic moment as velocity-space coordinates,

    ∂_t ( h_i-  F_Mi e_i φ /T_i  ) + (v_∥ + ṽ_M) ·∇ h_i
         + v_E·∇( F_Mi + h_i) = 0,

where the fluctuating electrostatic potential φ is self-consistently determined by quasineutrality

    T_i/e_i n_i∫d^3 v h_i  = τ( φ - ⟨φ⟩_yθ) + φ,

with τ = T_i/Z_i T_e and = B/B. The overlines indicate gyro-averages, with the φ average taken at fixed gyrocentre position R, and h_i at fixed real space position r. The electrons were assumed in (<ref>) to respond adiabatically only to potential variations within flux surfaces φ-⟨φ⟩_yθ, as they cannot move radially due to their small gyroradii. The gyrokinetic equation includes advection by the magnetic drift ṽ_M and by the gyro-averaged E×B-drift v_E = ×∇φ/B, with the latter determining the time-and-volume-averaged radial heat flux ⟨ Q ⟩_txyθ, with

    Q =  1/⟨ ∇ x⟩_θ∫d^3 v  m_i v^2/2h_iv_E ·∇ x ,

which is of primary interest to ascertain the level of transport caused by turbulent fluctuations.




Propagating zonal flows.– In gyrokinetic simulations of strongly driven ITG turbulence, ZF activity is apparent at multiple scales, see Fig.<ref>, showing large-scale stationary ZFs and small-scale propagating ZFs. The corresponding Fourier spectrum in Fig.<ref> exhibits stationary ZFs and GAMs at large radial scales k_x ρ_i ∼ 0.1, and a new propagating ZF mode, the toroidal secondary. It is predominant at smaller scales k_x ρ_i ∼ 0.5, and has a frequency ω∼ k_x v_Mx set by the radial magnetic drift velocity v_Mx≡ρ_i v_Ti/R, defined such that

    ṽ_Mx = v_Mxsinθv_∥^2 + v_⊥^2/2/v_Ti^2

in a large aspect-ratio tokamak of circular cross-section.

Unless specified otherwise, the gyrokinetic simulations presented in this Letter are performed using  <cit.>. The simulations model a flux surface at half-radius r_0=a/2 of the Cyclone Base Case, a tokamak with flux surfaces of circular cross-section and inverse aspect ratio R/a=2.8, where a is the tokamak minor radius. The safety factor q gives the magnetic field line pitch and sets the connection length L_∥ = q R, with 2π L_∥ the distance along the field line corresponding to one poloidal turn. The safety factor is varied between simulations, as is the ion temperature gradient R/L_T≡ -R  dln T_i/dr, while the density gradient and other quantities are held fixed at the usual Cyclone Base Case values. These are listed in the Supplementary Material, alongside the numerical parameters used in all simulations.

Compared to the small-scale, fast-oscillating toroidal secondary, the background turbulence has long radial wavelengths and evolves slowly in time, as will be shown explicitly below. The physics of the toroidal secondary may therefore be understood by considering a secondary mode growing and propagating over a background stationary (∂_t=0) streamer (k_x=0) mode, referred to as the primary mode and taken to be representative of the turbulence. This secondary model was originally considered by <cit.>, where the magnetic drift and parallel streaming in (<ref>) were ignored, leading to a purely growing mode. As shown in Fig. <ref>, gyrokinetic simulations of the secondary additionally exhibit the small-scale oscillating ZFs observed in the fully nonlinear gyrokinetic simulations, and a purely growing mode is observed only at large primary drive or long ZF wavelengths. We note that dispersion relations for the various modes may be derived analytically using the secondary formalism. These will be presented in future work, as we here focus on the phenomenology and physical mechanism of the toroidal secondary, and its effect on the turbulence.



To understand the mechanism for ZF drive and propagation, let us consider the vorticity equation at long perpendicular wavelengths,

    n_i ∂_t ⟨ρ_i ∇ x∂_x v_E^Z⟩_θ = 
       v_Ti∂_x ⟨∫d^3 v (  (v_⊥^2 /2Ω_i^2) ∇ x ·∇v_E ·∇ h_i + ṽ_Mxh_i )⟩_yθ,

which is derived from the time-derivative and flux-surface average of (<ref>). The vorticity equation describes the evolution of the ZF velocity v_E^Z=∇ x∂_x ⟨φ⟩_y/B due to nonlinear Reynolds and diamagnetic stresses, and due to the ṽ_Mx-induced Stringer-Winsor (SW) force <cit.>, which is dominant for the toroidal secondary. The SW force results from up-down asymmetries in (P+e_i n_i φ)^Z, with the fluctuating pressure

    P ≡∫d^3 v m_i/2( v_∥^2 + v_⊥^2/2) δ f_i,

as the magnetic drift points radially outward above the tokamak midplane and inward below it, reflected in the sinθ of (<ref>) and sketched in Fig. <ref>.

As illustrated in Fig. <ref>, the (P+e_i n_i φ)^Z asymmetry can be generated by the ZF itself, through the compression induced by the varying ZF velocity v_E^Z∝ 1/B between the tokamak outboard and inboard sides. In the case of the toroidal secondary, up-down asymmetry in P^Z is also generated nonlinearly, due to the radial magnetic drift (<ref>) advecting eddies sheared by the zonal flow, as explained in Fig. <ref>. The (P+e_i n_i φ)^Z asymmetry is generated on flux surfaces neigbouring the ZF perturbation, causing the mode to propagate and also grow provided the up-down asymmetry is sufficiently large. At long radial wavelengths, the ZF inertia in (<ref>) (the term with ∂_t) becomes negligible, such that the two mechanisms of (P+e_i n_i φ)^Z asymmetry generation described in Fig. <ref> must cancel each other. Toroidal secondary modes therefore oscillate more slowly than GAMs (see Fig. <ref>), for which the compressibility-induced SW force is balanced by the ZF inertia. We note that the toroidal secondary mechanism is distinct from that of <cit.>, who consider the pressure asymmetry generation due to a combination of ZF shearing and the effect of magnetic shear ŝ.

The physical picture outlined above did not include the effects of parallel streaming, which acts to short-circuit the up-down (P+e_i n_i φ)^Z asymmetries. Therefore, the toroidal secondary is only found at short radial wavelengths, where ω∼ k_x v_Mx≳ v_∥·∇∼ v_Ti/qR. At longer radial wavelengths, the ZFs do not oscillate, as seen in the gyrokinetic simulations of the secondary (Fig. <ref>) or of the nonlinear turbulence (Fig. <ref>). Furthermore, to avoid kinetic damping by the magnetic drift, the toroidal secondary mode requires a sufficiently large primary drive v_Ex≳ v_Mx to become unstable, as shown in Fig. <ref>. Finally, the toroidal secondary is stabilised by finite Larmor radius effects, such that the growth rate peaks at k_x ρ_i ∼ 0.5 in Fig. <ref>, explaining the prominent ZFs at this scale in fully nonlinear simulations (Fig. <ref>).




Zonal flows and turbulence saturation.– To ascertain the role of the various zonal flow modes in setting the turbulence saturation level, we consider the energy transfer in k_x caused by ZFs at different scales. Let us define the low-pass-filtered distribution

    h_i, K_x≡∫_-K_x^K_xdk_x ĥ_i(k_x) e^i k_x X,

where the Fourier modes ĥ_i(k_x) are normalised such that h_i = lim_K_x →∞ h_i, K_x. Then, the low-pass filtered gyrokinetic free-energy

    ℰ_K_x≡∑_s T_s ⟨∫d^3 v ( δ f_s,K_x)^2/2 F_Ms⟩_xyθ

may be shown to evolve as 

    ∂_t ℰ_K_x + 𝒯_K_x = ℐ_K_x - 𝒟_K_x,

with the energy transfer rate

    𝒯_K_x = ⟨ T_i ∫d^3 v h_i, K_x/F_Miv_E ·∇ h_i ⟩_xyθ,

the injection rate

    ℐ_K_x = L_T^-1⟨∫d^3 v  m_i v^2/2h_i, K_xv_E ·∇ x ⟩_xyθ,

and a dissipation rate 𝒟_K_x, e.g. from collisions or numerical dissipation.


In Fig. <ref>, the small-scale ZFs are shown to play a crucial role in transferring energy from large to small radial scales, as they are the largest contributors to 𝒯_K_x at radial wavenumbers larger than the toroidal secondary scale K_x ρ_i ∼ 0.5 - 0.7; for larger K_x the nonzonal contribution is dominant. The contribution to 𝒯_K_x from the large-scale ZFs is subdominant, though we note that numerical experiments artificially removing the ZFs at small k_x were found to alter the turbulence saturation level.




Turbulence scaling laws.– We now consider how the toroidal secondary modes may be incorporated in a theory of turbulence saturation. Strongly driven ITG turbulence in tokamaks has previously been modelled <cit.> using the critical balance conjecture <cit.>, which posits the nonlinear transfer rate ω_NL and parallel propagation rate ω_∥ of the saturated turbulence to be of the same order: at every scale, ω_NL∼ v_Ex/l_x ∼ω_∥∼ v_Ti / l_∥, with the rates estimated from (<ref>). Here, l_x, l_y, l_∥ are the radial, binormal, and parallel length scales of the turbulence. At the turbulence outer scale, denoted by `o' superscripts, ω_NL must also balance the energy injection rate ω_⋆^T, such that ω_NL^o  ∼ω_⋆^T,o∼ (v_Ti/L_T) ρ_i/l_y^o. Furthermore, the parallel outer length scale in a tokamak is assumed to be the connection length l_∥^o ∼ L_∥ = q R.

In the original theory of <cit.>, the turbulent eddies were assumed to be isotropic in the directions perpendicular to the magnetic field, i.e., l_x ∼ l_y. In fact, the turbulence is anisotropic at the outer scale due to the small-scale ZFs, with the turbulence level regulated to be near the marginal stability threshold of the toroidal secondary, v_Ex∼ v_Mx (see Fig. <ref>). This corresponds to a grand critical balance at the outer scale, ω_NL^o ∼ω_∥^o ∼ω_⋆^T,o∼ω_Mx^o with ω_Mx∼ v_Mx/l_x, in agreement with previously unexplained experimental observations <cit.>. 

Grand critical balance leads to the following scalings for the turbulence length scales, amplitude, and heat flux at the outer scale:

    l_y^o/ρ_i∼q R/L_T, l_x^o/ρ_i∼ q, e_i φ^o/T_iR/ρ_i∼q R/L_T,  ⟨ Q ⟩_txyθ/Q_gB∼q R/L_T,

with the gyro-Bohm heat flux Q_gB = n_i T_i v_Ti (ρ_i/R)^2. Notably, the heat flux is predicted to scale only linearly with the background temperature gradient, instead of the cubic scaling predicted for perpendicularly isotropic eddies <cit.>[The study by <cit.> also compared their scalings with gyrokinetic simulations of the Cyclone Base Case. However, due to computational limitations, the range of temperature gradients was limited to R/L_T ≤ 17.5, where a cubic fit is a good approximation to an offset linear dependence, see Fig. <ref>.]. This revised heat flux scaling is satisfied in gyrokinetic simulations, as shown in Fig. <ref>. Furthermore, we verify the grand critical balance scalings for the turbulence amplitude and outer scale by considering the fluctuation spectra, as shown in Figs. <ref> and <ref>. 



Finally, we note that, as q is increased, there is a clear scale separation between the turbulence outer scale k_x^o ρ_i ∼ 1/q, and the toroidal secondary scale k_x ρ_i ∼ 0.5, where the turbulence spectrum exhibits only a small local peak, see Fig. <ref>. Therefore, the grand critical balance scalings (<ref>) also justify the secondary model considered in Figs. <ref> and <ref> to describe the toroidal secondary mode, as the turbulence evolves over long spatio-temporal scales k_x^o ρ_i ∼ 1/q and ω_NL^o ∼ v_Ti/qR compared to the toroidal secondary's high wavenumber k_x ρ_i ∼ 1 and frequency ω^ZF∼ v_Ti/R. Due to the fast ZF oscillation, the shearing of turbulent eddies is not coherent, and must be modelled diffusively instead. To affect turbulence saturation, the ZFs must cause diffusion over an eddy length scale l_y^o in a nonlinear time, i.e., (l_y^o)^2 ∼ D/ω_NL^o, where D∼ (v_E^Z)^2 / ω^ZF is the diffusivity due to ZFs. Using (<ref>), the zonal flow energy 

    E^ZF≡e_i^2/2n_i T_i⟨∫d^3 v φ^Z( φ^Z - φ^Z) F_M ⟩_txθ

is thus predicted to scale as

    E^ZF/T_i≈⟨(v_E^Z/v_Ti)^2 ⟩_txθ∼ q ( R/L_T)^2 ( ρ_i/R)^2.

Here, the quantity φ^Z in (<ref>) corresponds to two successive gyro-averages, first at fixed R and then at fixed r and were approximated in the long-wavelength limit k_⊥ρ_i ≪ 1 in (<ref>). As shown in Fig. <ref>, the predicted ZF scaling (<ref>) is well satisfied by the small-scale toroidal secondary modes in gyrokinetic simulations.

Discussion.– In this Letter, we have presented the toroidal secondary, a new mode which leads to small-scale propagating zonal flows in nonlinear gyrokinetic simulations of ITG turbulence (Fig. <ref>). The toroidal secondary regulates the turbulence saturation level to be near the threshold turbulence amplitude v_Ex∼ v_Mx for which the toroidal secondary mode becomes unstable (Fig. <ref>). As a consequence, the turbulence follows the grand critical balance previously observed experimentally <cit.>, leading to turbulence scaling laws (<ref>) satisfied in gyrokinetic simulations (Fig. <ref>). 

We note that toroidal secondary modes resulting from electron-scale turbulence could differ considerably in character from those presented in this Letter, as the ZF inertia will play a more important role. For ion scale modes, the ZF inertia does not contribute significantly to the toroidal secondary mechanism due to the modified adiabatic electron response in (<ref>), which leads to a small ZF inertia at long wavelengths, see (<ref>).

The toroidal secondary modes shown in Fig. <ref> are reminiscent of the avalanches previously observed in gyrokinetic simulations <cit.>, whose physical mechanism remains debated. Recent efforts <cit.> to model the avalanches have considered the coupling of zonal flows with both the turbulence and with toroidal geometric effects, though nonlinear fluxes were neglected, precluding the toroidal secondary modes presented here (see Fig. <ref>). From Fig. <ref>, the toroidal secondary is seen to have a typical propagation velocity ∼ 2 v_Tiρ_i/R, comparable to previously reported avalanche propagation speeds <cit.>. Furthermore, the turbulence amplitude threshold required to destabilise the toroidal secondary could explain the sandpile-like behaviour associated with avalanches.


This work was supported by U.S. DOE DE-AC02-09CH11466 and DE-FG02-93ER54197, by Scientific Discovery Through Advanced Computing (SciDAC) Grant No. UTA18000275, and by the Engineering and Physical Sciences Research Council (EPSRC) [EP/R034737/1]. The simulations presented in this article were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.












§ GYROKINETIC SIMULATION DETAILS


All simulations presented in this Letter modelled a flux-tube extending for a single poloidal turn in the Cyclone Base Case tokamak. In addition to the parameters given in the main text, a magnetic shear of ŝ = dln q/dln r = 0.8 was used, and the density gradient was held fixed at R/L_n ≡ -R  dln n_i/dr = 2.2. A hydrogen plasma with Z_i=1 was considered, and the electron to ion temperature ratio was set to τ=1.

The numerical parameters employed in the  <cit.> gyrokinetic simulations are provided in Table <ref>. The quantities L_x/ρ_i and L_y/ρ_i denote the box size in the radial and binormal directions, respectively, and may be related to the minimum wavenumber through k_min = 2π/L. Furthermore, N_θ gives the parallel resolution, while N_x and N_y denote the number of Fourier modes in the radial and binormal directions. The velocity resolution is given by N_μ and N_v_∥ for the magnetic moment and the parallel velocity, respectively, with v_∥/v_Ti∈ [-3.5,3.5] for all simulations. Finally, D_h denotes the hyperdissipation parameter.



As shown in Table <ref>, the simulations using the  code <cit.> could be run at fixed numerical parameters as κ and q were varied, as the computational cost is small owing to the use of GPUs. In the κ-scan, temperature gradient values of κ∈[6.9, 13.9, 20.9, 27.8, 34.7, 41.7, 48.6] were considered, at fixed safety factor q=1.4, while for the q-scan, κ=13.9 and q ∈ [1.4, 2.8, 4.2, 5.6, 7.0]. The velocity resolution is here given by the number of Hermite and Laguerre moments N_h and N_l, respectively.





