Non-Oberbeck–Boussinesq zonal flow generation

Novel mechanisms for zonal flow (ZF) generation for both large relative density fluctuations and background density gradients are presented. In this non-Oberbeck–Boussinesq (NOB) regime ZFs are driven by the Favre stress, the large fluctuation extension of the Reynolds stress, and by background density gradient and radial particle flux dominated terms. Simulations of a nonlinear full-F gyro-fluid model confirm the predicted mechanism for radial ZF propagation and show the significance of the NOB ZF terms for either large relative density fluctuation levels or steep background density gradients.


I. INTRODUCTION
Self-organization from turbulent to coherent states is a ubiquitous process in fluids.In particular, much interest and effort has been drawn to the formation of zonal flows (ZFs) [1][2][3].These coherent flows arise in atmospheres, in the form of banded cloud structures on Jupiter [4], Saturn's north-polar hexagon [5] or mid-latitude westerlies on earth and in the ocean as stationary jets [6].In magnetized fusion plasmas ZFs are key players for the reduction of the radial transport of particles and heat and for the transition to improved confinement regimes in tokamaks [7][8][9][10][11][12].
Reynolds stress is quintessential for ZF generation in all fluids [1][2][3][13][14][15][16][17], but in magnetized plasmas also other stresses like the Maxwell [18,19] or the diamagnetic stress [20,21] can become significant.Virtually all of the work on ZF theory so far rely on δf models [1,13,22], which invoke the so called Oberbeck-Boussinesq (or thin layer) approximation [23,24].However, the latter breaks down, if the background density varies over more than one order of magnitude or if the relative density fluctuations exceed roughly 10 percent.This for example prevails in the edge of tokamak fusion plasmas, where experimental measurements typically feature relative density fluctuation levels around the order 0.1 in the edge and up to unity at the last closed flux surface [25][26][27][28][29][30][31][32][33].Moreover, typical edge background density gradient (e-folding) lengths reach from 50ρ s0 in lowconfinement to 10ρ s0 in high-confinement tokamak plasmas [34,35].Here, ρ s0 := √ T e0 m i /(eB 0 ) is the drift scale with reference electron temperature T e0 , ion particle charge e, ion mass m i and reference magnetic field B 0 .
Non-Oberbeck-Boussinesq (NOB) effects on ZF generation are an unresolved issue.However, theoretical and experimental studies of poloidal ZFs in the edge of fusion plasmas indicate that unknown mechanisms beyond the Reynolds stress exist [36] and that steep back- * E-mail: markus.held@uibk.ac.at ground density gradients and large relative density fluctuations affect the poloidal ZF dynamics [15,37,38].Moreover, the importance of large relative density fluctuations for toroidal momentum transport, as suggested by theoretical estimates in the strong and weak turbulence regime [39,40] and experimental measurements in the TORPEX and PANTA device [41,42], point towards a similar significance for poloidal momentum transport.
In the following we generalize the theory of ZFs to NOB effects.To this end, we decompose the density and electric potential of a full-F gyro-fluid model of a magnetized plasma [43] with the help of a density weighted Favre average [44].This well known decomposition strategy in compressible fluid dynamics (see e.g.[45]) is here for the first time introduced to plasma physics and enables us to disentangle the density fluctuations from the ZF dynamics, while retaining the relevant physical effects.As a result, we identify novel agents in the poloidal ZF dynamics, which become significant for high relative density fluctuations or steep background density gradients.We confirm the herein proposed NOB mechanism for radial advection of ZFs with the help of numerical simulations of a fully nonlinear model for drift wave-ZF dynamics.The exploited model is based on the specified extension of the Hasegawa-Wakatani model to the full-F framework.Additionally, we show how the ZF dynamics is distributed among the proposed NOB actors and provide scalings with collisionality, reference background density gradient length and the maximum of the relative density fluctuation amplitude.

II. ZF THEORY A. δf formalism
We start our discussion with a short re-derivation of the conventional ZF equation and Reynolds stress from a cold ion δf gyro-fluid model, which couples small relative density fluctuations to the electric potential via E × B advection and linear polarization [46][47][48] Here, δn := n/n G − 1 is the relative electron density fluctuation, δN := N/n G − 1 is the relative ion gyrocenter density fluctuation, φ is the electric potential and Ω 0 := eB 0 /m i is the ion gyro-frequency.The reference background density n G (x) refers to a constant reference background gradient length L n := −1/∂ x ln (n G /n 0 ) with constant reference density n 0 .For the sake of simplicity the magnetic field B = B 0 is assumed constant and the unit vector in the magnetic field direction is b := B/B 0 = êz .The perpendicular gradient and the E × B drift velocity are defined by and u E := b×∇φ/B 0 , respectively.The term Λ δ denotes a closure for the parallel dynamics, which is discussed later in more detail.Taking the time derivative over the polarization equation (1c) yields the δf drift-fluid vorticity density equation with the linear E × B vorticity density ).Now we apply the average over the "poloidal" y coordinate h := L −1 y Ly 0 dy h to Eq. (2), which is the 2D equivalent of a flux surface average.Reynolds decomposition h = h + h and integration over the "radial" coordinate x result in the δf evolution equation for poloidal ZFs [13] Here, we introduced u x := −∂ y φ/B 0 , u y := ∂ x φ/B 0 and the anticipated Reynolds stress R := u x u y [49], where u x u y = u x u y + u x u y and u x = 0 was used.In passing we note that we assume that radial boundary conditions give rise to no additional terms in Eq. ( 3) and for the remainder of this letter.

B. Full-F formalism
In full-F theory the splitting of the gyro-fluid moment variables into fluctuating and background parts is avoided and the quasi-neutrality constraint for electrons and ions is rendered by the nonlinear polarization equation [43].The cold ion full-F gyro-fluid model [50,51] evolves the full electron density n and ion gyro-center density N .In the gyro-center E × B drift velocity by Both, the latter ponderomotive correction and the polarization charge nonlinearity on the left hand side of Eq. (4c) are crucial for energetic consistency and an exact momentum conservation law [52].We refer to the parallel closure term Λ later on.In the long wavelength limit we can again reformulate Eqs.(4b) and (4c) into a drift-fluid vorticity density equation where the nonlinear E × B vorticity density is given by As above, we obtain the averaged poloidal momentum equation [38,53,54] The divergence of the full-F stress drive terms of Eq. ( 6) is related to the averaged radial flux of vorticity density minus the ponderomotive correction via the full-F Taylor identity The interpretation of Eq. ( 6) is problematic since (i) absolute density fluctuations n arise instead of relative density fluctuations n/ n , (ii) the time evolution of the averaged poloidal momentum nu y is given in terms of the averaged poloidal velocity u y and (iii) background density gradient ∂ x ln n effects are not obvious.Despite these obstacles Eq. ( 6) has been recently used to show that the second term, occasionally misinterpreted as advective, and the cubic term can be comparable to the Reynolds stress related drive ∂ x ( n R) [15,37,38].Thus, we go a step further and utilize a density weighted Favre decomposition instead of the Reynolds decomposition according to h := [[h]] + h and [[h]] := nh / n [44].Note that the Favre decomposition reduces to the Reynolds decomposition if the density n is only a function of x.Now we combine the poloidal average of Eq. (4a) divided by n with Eq. ( 6) divided by n to obtain a ZF evolution equation for the Favre averaged poloidal velocity where we used ] can be rewritten into Consequently, the first term −∂ x F on the right hand side of Eq. ( 9) is the superposition of the conventional Reynolds stress drive ) and the triple fluctuation drive The novel second term on the right hand side of Eq. ( 9) represents radial advection of poloidal ZFs [[u y ]].Its direction depends on the sign of the averaged radial particle flux , which is typically positive, so that T 4 describes an outward pinch of ZFs.The novel third term T 5 := F/L n on the right hand side of Eq. ( 9) is proportional to the inverse of the background density gradient length 1/L n := −∂ x ln n .This term is large for small reference background density gradient lengths L n , or has large radially localized values if the density profile n develops into a staircase like pattern [55].In contrast to the Favre stress drive −∂ x F, the background density gradient drive T 5 contributes to the ZF generation even if the Favre stress is radially homogeneous ∂ x F = 0. Remarkably, the background density gradient drive remains finite in the small relative density fluctuation limit, where the density n is only a function of x and the Favre stress F resembles the conventional Reynolds stress R.
In order to interpret the dynamics of the background density gradient drive T 5 let us assume for a moment that the turbulent viscosity hypothesis [13,56].In this case Eq. ( 9) reduces to a simple advection-diffusion equation for ZFs where the background density gradient pinch velocity V := ν T /L n appears now in addition to the radial outward pinch velocity [[u x ]].The direction of the additional pinch depends on the sign of the turbulent viscosity ν T .Finally, we extend the theory for energy transfer inside the kinetic E × B energy E(t) := m i dA nu 2 E /2 to the full-F formalism.Here, the Favre decomposition E = E 0 + E 1 is pivotal to derive the conservation laws for the zonal (or mean) 2 /2 and turbulent part E 1 (t) := m i dA n u 2 E /2 of the kinetic E×B energy and supersedes the Reynolds decomposition in the δf formalism [19,21].With the help of Eqs. ( 8) and ( 9) we obtain the conservation laws for the zonal and turbulent kinetic E × B energy This unveils that the Favre stress term n F ∂ ∂x [[u y ]] is the central mechanism for energy transfer between the zonal and turbulent kinetic E × B energy.As a consequence, density fluctuations (cf.Eq. ( 10)) manifest as an additional transfer channel in the full-F formalism.

III. PARALLEL CLOSURES
Self-sustained drift wave turbulence is maintained by the non-adiabatic parallel coupling of the relative density fluctuations and the electric potential, which can arise due to various mechanisms.Here, we exemplarily consider resistive drift wave turbulence, which arises due to resistive friction between electrons and ions along the magnetic field line.This mechanism enters the 2D gyro-fluid models via the parallel closure terms (Λ δ or Λ) of the Hasegawa-Wakatani (HW) type as summarized in Table I.Here, we introduced the full-F adia- baticity parameter α := T e0 k 2 /(η e 2 n 0 Ω 0 ) with parallel wavenumber k and parallel Spitzer resistivity η := 0.51m e ν e /(ne 2 ) [61,62].In the electron collision frequency ν e the Coulomb logarithm is treated as a constant so that η has no explicit dependence on n.As opposed to this in δf models the density dependence in the collision frequency ν e (n) ≈ ν e0 is completely neglected so that α δ := T e0 k 2 /(0.51m e ν e0 Ω 0 ) reduces to a parameter.
Only then, the poloidal variations of the adiabaticity parameters vanish ( α δ = α = 0), and the full-F and δf closures coincide in the limit of n ≈ n G and δn 1.

IV. SIMULATIONS
We use the open source library Feltor [63] to numerically solve the full-F gyro-fluid Eqs. ( 4 NOB effects on drift wave-ZF dynamics, as it is described by Eqs. ( 8) and ( 9) with Λ = 0, are in this setup studied by varying the adiabaticity parameter (or inverse collisionality) α and the reference background gradient length L n .In Fig. 1 we show that L n crucially determines the time evolution of ZFs in the high collisionality regime with α = 0.0005.While stationary ZFs emerge for L n = 128ρ s0 , a radial outward pinch of ZFs occurs for a four times smaller reference background density gradient length L n = 32ρ s0 .In this steep gradient and high collisionality regime the ZF signature is no longer solely determined by the conventional Reynolds stress drive, which is illustrated in Fig. 2. The Reynolds stress drive T 1 is here comparable to the radial advection term T 4 , which explains the observed radial outward propagation of ZFs in Fig. 1.
The radial profile of the terms of the right hand side of Eq. ( 9) for α = 0.0005 and Ln = 32ρs0.The ZF signature of the radial advection term T4 is comparable to the Reynolds stress T1.
In the following the parametric dependence of each term T i on the right hand side of Eq. ( 9) is investigated.To this end, the contribution of each term T i on ZF evolution is measured by taking the L 2 norm, denoted by h 2 , of the time integrated contribution.Following this, we propose a measure of the relative ZF contribution In Fig. 3a we show that the relative contribution M i of the NOB ZF terms (T 2 , . . ., T 5 ) decreases with the reference background density gradient length L n in the high collisionality regime (α = 0.0005).The summed up relative contribution of the NOB ZF terms exceeds the one of the conventional Reynolds stress for L n = 32ρ s0 .For steep reference background density gradients the radial advection term T 4 exhibits the largest relative contribution to the ZF dynamics of all the NOB terms.
In Fig. 3b the dependence of the relative importance of each term on the adiabaticity parameter α is depicted for a fixed reference background density gradient length L n = 32ρ s0 .While the conventional Reynolds stress term is again the dominating ZF contributor in particular for small collisionalities, all the NOB terms, except the background density gradient drive T 5 , gain in importance for higher collisionalities.Interestingly, in the small collisionality regime (α ≥ 0.01) the background density gradient drive T 5 exceeds all the remaining NOB actors.The quadruple fluctuation drive T 2 is for all studied parameters the smallest contributor to the ZF dynamics.
The dependence of the ZF terms on the time averaged maximum of the relative density fluctuation level n/ n ∞ t is shown in Fig. 4. Here, we denote the The relative contributions of the NOB ZF terms increase with the relative density fluctuation amplitude.In particular they can amount to roughly two thirds of the ZF dynamics.

V. CONCLUSION
We have generalized the ZF equation (3) to account for NOB effects in Eq. ( 9).Most importantly, the former Reynolds stress R is replaced by the Favre stress F, which adds to its predecessor in case of high relative density fluctuations.The latter is accompanied by two new agents in the NOB ZF Eq. ( 9).The first of these radially advects ZFs by the Favre averaged radial drift velocity, which is proportional to the averaged radial particle flux.The second term scales inversely with the background density gradient length and affects the ZF dynamics even if the relative density fluctuations are small or if the Favre stress is radially homogeneous.Thus, this term may be of significance in or during the formation of radial transport barriers, where steep density profiles form with strongly reduced radial particle transport.
Additionally we extended the ordinary and modified HW model to the full-F theory.We simulated the full-F gyro-fluid model with the modified HW closure to numerically corroborate our theoretical results.The simulations successfully reproduced the predicted radial advection of ZFs, which appeared for small reference background density gradient lengths and large averaged radial particle flux.Moreover, our numerical parameter study showed that the NOB ZF drives can be comparable to the Reynolds stress drive in the herein scanned parameter range.In particular the deviation between the Reynolds and Favre stress drive increases with the relative density fluctuation amplitude, collisionality and inversely with the reference background density gradient length.This deviation is mainly reasoned in the triple fluctuation drive.Its importance in steep background density gradient regimes is in qualitative agreement with the theoretical estimate in the strong turbulence regime [38].A similar dependence as for the Favre stress drive is found for the radial ZF advection mechanism.For the background density gradient drive only a dependence on the reference background density gradient is observed.
The presented results strongly argue in favor of the development and application of full-F gyro-fluid or gyrokinetic models for simulation of fusion edge plasma turbulence, and in general demonstrate exemplarily the relevance of NOB effects for ZF formation in fluids and plasmas with large fluctuations and inhomogeneities.The latter conditions prevail e.g. during the low-to highconfinement mode transition.Thus, a consistent full-F simulation approach of this phenomenon is crucial to allow for the herein presented NOB ZF mechanisms.
Finally, we emphasize that the relative error between the Favre and Reynolds average of the poloidal velocity, derived to | [[ u y ]] / u y |, is typically below a few percent.Thus, our proposed NOB ZF theory is also applicable to experimental measurements of the Reynolds averaged poloidal velocity u y .
) with the modified HW parallel closure of Table I.Numerical stability is ensured by adding hyperdiffusive terms of second order −ν∇ 4 ⊥ n and −ν∇ 4 ⊥ N to the right hand side of Eqs.(4a) and (4b).Moreover, we append the right hand side of Eqs (4a) and (4b) by a density source of the form ω S z Θ (z) with z := g(x) (n G − n ) to maintain the initial profile in a small region x ∈ [0, x b ].Here, we defined the Heaviside function Θ(z) and g(x) := [1 − tanh (x − x b )/σ b ] /2.The corresponding parameters are fixed to ν = 5 × 10 −4 c s0 ρ 3 s0 , ω S = 0.1Ω 0 , x b = 0.1L x and σ b = 0.5ρ s0 with cold ion sound speed c s0 := ρ s0 Ω 0 .The box with size L x = L y = 128ρ s0 is resolved by a discontinuous Galerkin discretization with P = 3 polynomial coefficients and at least N x = N y = 256 equidistant grid cells.The initial (gyro-center) density fields n(x, 0) = N (x, 0) = n G (x) (1 + δn 0 (x)) consist of the reference background density profile n G , which is perturbed by a turbulent bath δn 0 (x).

FIG. 3 .
FIG. 3. (a)The NOB ZF terms decrease with the reference background gradient length Ln in the high collisionality regime (α = 0.0005).(b) For a fixed Ln = 32ρs0 all the NOB ZF terms significantly contribute to the ZF dynamics in the high collisionality regime.As opposed to this, only the background density gradient drive T5 remains alongside the Reynolds stress drive T1 in the small collisionality regime.

TABLE I .
HW closures for δf and full-F models.