The collisional drift wave instability in steep density gradient regimes

The collisional drift wave instability in a straight magnetic field configuration is studied within a full-F gyro-fluid model, which relaxes the Oberbeck-Boussinesq (OB) approximation. Accordingly, we focus our study on steep background density gradients. In this regime we report on corrections by factors of order one to the eigenvalue analysis of former OB approximated approaches as well as on spatially localised eigenfunctions, that contrast strongly with their OB approximated equivalent. Remarkably, non-modal phenomena arise for large density inhomogeneities and for all collisionalities. As a result, we find initial decay and non-modal growth of the free energy and radially localised and sheared growth patterns. The latter non-modal effect sustains even in the nonlinear regime in the form of radially localised turbulence or zonal flow amplitudes.


I. INTRODUCTION
In a series of seminal theoretical works in the early 1960s it has been established that low-frequency small scale instabilites are naturally immanent to magnetically confined plasmas [1][2][3][4][5][6].This is due to the inherent plasma pressure gradients, which nurture these so called drift wave (DW) instabilities.These DW instabilities drive the turbulent cross-field transport of particles and heat, which exceeds predictions from classical and neo-classical theory and remains a serious barrier for sufficient plasma confinement in laboratory plasmas.
After nearly 60 years density gradient driven DW instabilities remain an active theoretical research field.The latest efforts focus on a unified description of the collisional and collisionless DW instability [28,29] or proof instability for collisionless DWs in sheared magnetic fields after decades of misconception [19].
Another outstanding regime of interest is that of large inhomogeneities, in particular if the background density varies over more than one order of magnitude.This for example prevails in tokamak fusion plasmas, where steep background density gradients can emerge due to the for- * E-mail: markus.held@uibk.ac.at mation of internal or edge transport barriers [30,31].Under these circumstances the Oberbeck-Boussinesq (OB) approximation [32,33] breaks down, which is throughout applied in former studies of DW instabilites.Thus, a rigorous stability analysis for large inhomogeneities must relax the latter assumption.
In spite of that, recent OB approximated analysis of large inhomogeneity effects on trapped electron or ion temperature gradient driven DW instabilities indicate their relevance for transport bifurcation in the edge of tokamaks [34][35][36].Initial investigations of large density inhomogeneity effects on the stability and dynamics of collisional DWs rest partly upon the OB approximation and exploit further approximations in their linear analysis [37,38].
In this contribution we investigate the linear dynamics of the collisional DW instability for steep background density gradients.This is achieved by consistently linearising a non-Oberbeck-Boussinesq (NOB) approximated full-F gyro-fluid model [39], which accurately accounts for collisional friction between electrons and ions along the magnetic field.We numerically solve the generalised eigenvalue problem to show that the growth rate and real frequency deviate by factors of order one from the OB approximated case in the steep background density gradient regime.In this regime, the NOB approximated model is strongly non-normal.As a consequence the eigenfunctions are spatially localised, as opposed to the OB limit.Moreover, non-modal effects arise, resulting in initial decay and transient non-modal growth of the free energy and radially localised and sheared growth of an initially unstable random perturbation.The latter NOB signature lingers into the nonlinear regime, where radially localised turbulence amplitudes appear.

II. GYRO-FLUID MODEL
Our analysis is based on an energetically consistent full-F gyro-fluid model [39], which is derived by taking the gyro-fluid moments over the gyro-kinetic Vlasov-Maxwell equations [40].In order to ease the following discussion we assume constant temperatures, cold ions and a constant magnetic field B = B 0 with straight and unsheared unit vector b := B/B = êz .The resulting set of gyro-fluid equations consists of continuity equations for electron density n and ion gyro-center density N and the quasi-neutrality constraint where φ is the electric potential, Ω 0 := eB/m i is the ion gyro-frequency and ∇ ⊥ := − b × ( b × ∇) and ∇ := b • ∇ are the perpendicular and parallel gradient, respectively.The E × B drift velocity is defined by u E := b × ∇φ/B 0 .As opposed to this, the gyro-center E × B drift velocity ). Parallel collisional friction between electrons and ions is introduced on the right hand side of Eq. (1a) with parallel Spitzer resistivity η := 0.51m e ν e /(ne 2 ) [41,42].These closures of the Hasegawa-Wakatani (HW) type are obtained from the evolution equation for the parallel electron velocity.In the electron collision frequency ν e := ne 4 ln(Λ)/(3 (2π) 3 4  0 m e T 3 e ) the Coulomb logarithm ln(Λ) is treated as a constant [43].This is a reasonable approximation even if the density profile varies over several orders of magnitude.Consequently, the parallel Spitzer resistivity η has no explicit dependence on the electron density n, since we only retain the electron density n proportionality in the electron collision frequency ν e .
The two dimensional form of the presented full-F extension of the ordinary HW (OHW) model is obtained by rewriting ∇ 2 ln (n) to ∇ 2 ln (n/ n ) and by replacing the parallel derivative with the characteristic parallel wave-number, so that ∇ 2 = −k 2 .Here, the average over the "poloidal" y coordinate h := L −1 y Ly 0 dy h is introduced, which is the 2D equivalent of a flux surface average.With these manipulations Eq. (1a) reduces to [44] where the adiabaticity parameter is The nonlinear gyro-fluid model of Eqs. ( 2), (1b) and (1c) allows to study NOB effects on collisional DWs, since in this regime of large collisionality ν e ω and small Knudsen-number Kn := k λ mf p 1 the presented fluid approach is valid.

III. LINEARISED GYRO-FLUID MODEL
The linearised gyro-fluid model is obtained by expanding the nonlinear gyro-fluid model of Eqs. ( 2), (1b) and (1c) around a reference background density profile n G (x) according to The chosen exponential reference background density profile n G yields a constant density gradient (e-folding) length as in δf theory.Note that we found similar trends in our results for non-exponential reference background density profiles.Assuming that the relative fluctuation amplitudes δn ∝ δN 1 are small and the averaged and reference background density profiles coincide n ≈ n G yields the final form of the linearised gyro-fluid model As opposed to previously exploited linearised models [28,37,38,45,46] we do not apply the OB approximation in the linearised Eqs.(6a), (6b) and (6c) or in the further course of the calculation.Thus the derived set of linearised Eqs.(6a), (6b) and (6c) differs from linearised OB approximated models in two substantial aspects.First, we retain the background density on the right hand side of Eq. (6a) instead of assuming a constant reference density n G ≈ n 0 .Secondly, we preserve the term − 1 Ln ∂ ∂x φ on the left hand side of Eq. (6c), which originates from the nonlinear contribution of the polarization charge density.Both of these terms, but especially the NOB approximated resistive term, produce novel linear effects for collisional DWs in steep background density gradients as is shown in section IV.

IV. LINEAR EFFECTS
In the following we want to gain insight into the linear dynamics, in particular in stability and transient time behavior, of the linear model Eqs.(6a), (6b) and (6c).This is accomplished within a discrete approach, which utilises a Fourier transformation of the form e ikyy in the periodic poloidal coordinate y and a Galerkin approach with a sine basis in radial direction.This fulfills the chosen Dirichlet boundary conditions in radial direction.The resulting linear equation reads in matrix form with vector v := (δn ky , eφ ky /T e0 ) T and matrix with matrix The coefficients of the symmetric matrix E and D xx and of the skew-symmetric matrix D x are derived to 4(e Lx/Ln (−1) with radial box size L x .The radial wave-numbers are defined by k x := πm/L x and k x := πm /L x with modenumbers m ∈ N and m ∈ N, respectively.Note, that the matrix A of Eq. ( 8) is in general far from normal (AA † = A † A) with non-orthogonal eigenvectors but approaches a normal matrix (AA † = A † A) in the OB limit of very flat background density profiles [21]).This is best shown by characterizing the departure from normality by the the condition number of the matrix of eigenvectors V of A. Here, the matrix norm is induced by the free energy vector norm ||v|| := √ v † Mv with [21] M := 1 2 This normalises the condition number κ to unity for a normal system.In Fig. 1 we plot the condition number κ as a function of adiabaticity α and background density length L n , respectively.In contrast to the OB case, we universally obtain non-normality (κ > 1) for the NOB case for steep background density gradients (L n = 16ρ s0 ).In this regime the NOB condition number κ is at least a magnitude higher than its OB equivalent and approaches extremely large values for adiabaticities below α ≤ 0.01.For a fixed adiabaticity α = 0.005 the NOB dynamics are normal in the OB limit (L n ≥ 10 4 ), but are strongly non-normal for steep background density gradients.This behavior of the condition number κ suggests much larger non-modal effects for the NOB model than for the OB model.However, the condition number κ does not give insight into how the departure from normality affects the linear dynamics.Thus, we analyze in the following the modal and non-modal behavior of Eqs.(7).

A. Modal analysis
First, we address the eigenvalues ω(A) and the radial eigenfunctions.The numerically calculated growth rate γ(A) := Re(ω(A)) and real frequency ω R (A) := −Im(ω(A)) of the NOB case are compared to the maximum of its OB counterpart.In the OB limit the dispersion relation can be simply derived analytically [45,46] where we defined and the perpendicular wave number k ⊥ := k 2 x + k 2 y .Consequently, the real frequency ω R,OB and growth rate γ OB are derived to with real part Re(z) := −B 2 , imaginary part Im(z) := 4Bω * /Ω 0 and argument θ := arg(z) of the complex number z.
In Fig. 2 we show the normalised growth rate γ/ max kx (γ OB ) and real frequency ω R / max kx (ω R,OB ) as a function of the radial wave number k x for various background density gradient lengths L n and adiabaticities α.
Here, the NOB growth rates γ and real frequencies ω R exhibit significant deviations from the OB limit in particular for steep background gradients and for a range of typical adiabaticity parameters.In particular, the magnitude of the fastest growing mode differs by up to roughly a factor five and the radial wave-number of the fastest growing mode is also different for certain parameters.However, the NOB eigenvalues resemble the OB limit for very flat background gradients.
The radial eigenfunctions for the relative density fluctuation δn ky of the fastest growing k x mode are depicted in Fig. 3 for two different density gradient lengths L n .Remarkably, for steep background density gradients these eigenfunctions are spatially localised and do no longer coincide with the ordinary sine like eigenfunction of the OB     model.Moreover, the phase shift between the real and imaginary parts of the eigenfunctions leads to shearing in the x-y plane as we illustrate in section IV B. Again, for flat background density profiles the eigenfunctions transition into the OB approximated equivalent.

B. Non-modal analysis
We now face the question how non-modal effects manifest in numerical simulations of the nonlinear full-F ordinary HW model, given by Eqs. ( 2), (1b) and (1c).In particular, we study if initial or transient dynamics are pronounced during the linear phase and if these non-modal effects survive into the nonlinear regime.
It is well known that for a normal matrix the time evolution of the norm of the linear Eq. ( 7) is bounded by the spectral abscissa for t ≥ 0 since the matrix exponential reduces to ||e At || = e β(A)t .However, for a non-normal matrix the maximum growth estimate of the modal analysis of section IV A, determined by the spectral abscissa β, only holds for t → ∞ and the matrix exponential reduces to a loose upper bound ||e At || = κ(V)e β(A)t for t ≥ 0. As a consequence, a non-normal system may exhibit pronounced initial or transient phenomena, for which also estimates and bounds exist.In particular, the numerical abscissa η(A) := max k ω (A † + A)/2 1 represents an upper bound for t = 0 and the so called -pseudospectral abscissa α (A) yields estimates for transient phenomena [47].Although, the latter two approaches are useful to detect or quantify non-modal effects, we do not make use of them in the following discussion.Instead, we present a direct numerical approach to the initial value problem.
The numerical implementation of the latter full-F gyrofluid model utilises the open source library Feltor [48].This initial value code relies on a discontinuous Galerkin discretization, which is also used for verification of the herein presented Galerkin approach for the modal analysis of section IV A. We limit our study to a single exemplary initial condition, but note that in general the initial condition can be optimised to produce maximum growth at small, intermediate or large time [21].The chosen initial conditions mimics a random perturbation δn(x, 0) = δN (x, 0) = af bath (x) of amplitude a with vanishing electric potential φ(x, 0) = 0.
In Fig. 4 we show the temporal behavior of the square root of the normalised free energy norm ||v(t)||/||v(0)|| for various adiabaticities α in the steep gradient regime (L n = 16ρ s0 ).Here, the random bath initial condition with the small amplitude a = 10 −5 limits us to linear effects only.Interestingly, two clear footprints of nonmodal behavior emerge in the linear dynamics.First, initial decay of the square root of the normalised free energy norm appears for both the OB and NOB case, despite the fact that all eigenvalues are unstable (cf.Fig. 2).Secondly, the transient exponential growth at later times either surpasses (NOB) or falls below (OB and NOB) the spectral abscissa β.Both of these effects are due to the shrinking of non-orthogonal eigenvectors, which is intrinsic to non-normal systems.We refer the interested reader to [49] for an illustrative sketch of this phenomenon.The observed initial decay of the collisional DW instability is similar to that of the collisionless DW instability [27].However, for the latter instability transient amplification can trigger subscritical turbulence in the absence of linear instability.
During this linear growth phase non-modal features may appear in the spatial structure of the relative density fluctuation δn.This is depicted in Fig. 5 for the turbulent bath initial condition with a = 0.01 and for a steep background density profile (L n = 16ρ s0 ) and typical adiabaticity (α = 0.005).In contrast to the OB approxi- mated model, we observe sheared and localised growth of the initial perturbation in the steep background density regime.This is in qualitativ agreement with the previously reported NOB shearing effect of DWs of [50].The radial location of the strongest growth coincides approximately with the maximum of the absolute background density gradient |∂ x n G |.These NOB effects are again reasoned in the strong non-orthogonality of the eigenvector Matrix V, which is pronounced for a system with strong non-normality and by implication high condition number κ (cf.Fig. 1) Finally, we study how non-modality affects the turbulence intensity in collisional DW turbulence without and with zonal flows.Here, the underlying models are the NOB-extended OHW (equations ( 2), (1b) and (1c)) and modified HW (MHW) model [44], respectively.In Fig. 6 we show the spatial structure of the relative density fluctuation δn for both models and the latter initial condition but during the nonlinear phase.Without zonal flows (OHW) we observe a radial peaking of the maximum of the relative density fluctuation amplitude.Analogously, with zonal flows (MHW) a similar radial peaking appears for the zonal flow amplitude.Note that the radial localization of the turbulence intensity continues to exist at late turbulence saturation times if the radial particle transport is weak.This occurs for large adiabaticity (small collisionality).For both model cases, this constitutes a clear non-modal footprint in saturated collisional DW turbulence and shows that indeed these non-modal effects can survive into the nonlinear regime.As opposed to this, the radial peaking of the turbulence or zonal flow intensity is again absent in the OB approximated OHW or MHW model (cf.[45,46,51]).

V. DISUSSION AND CONCLUSION
We studied the collisional DW instability in a straight and unsheared magnetic field within a full-F gyro-fluid model, which relaxes the OB approximation.In the regime of steep background density gradients both the eigenvalues and eigenfunctions fundamentally deviated from former OB approximated investigations.In particular, our modal analysis demonstrated NOB corrections by factors of order one to the eigenvalues, highly non-orthogonal eigenvectors and spatially localised eigenfunctions for typical plasma parameters.Our non-modal analysis revealed initial damping and transient nonmodal growth of the free energy of an initially unstable random perturbation.Remarkably, this linear growth is radially localised and sheared.It was numerically shown that this NOB signature subsists into the nonlinear regime, where radially localised turbulence or zonal flow amplitudes emerge.
The herein presented results emphasise the need for NOB approximated models to consistently capture the (linear) dynamics of the collisional DW instability for large density inhomogeneities.For instance, this may prove necessary for the accurate calculation of transport levels in high-confinement tokamak plasmas within quasilinear gyro-kinetic or gyro-fluid models (e.g.: [52,53]).Finally, we conclude that the study of linear effects, that is solely based on a modal approach, may give misleading predictions since this approach overlooks non-modal features like initial and transient phenomena.
FIG. 4. The normalised energy norm as a function of normalised time is plotted for different adiabaticities α = {0.05,0.005, 0.0005} in the steep gradient regime Ln = 16ρs0.The initial amplitude is a = 10 −5 .

FIG. 5 .
FIG. 5.The spatial pattern of the relative density fluctuation are shown during the linear phase at time t = 750/Ω0 (NOB) and t = 500/Ω0 (OB).Localised growth and shearing of the initial relative density perturbation appear in the NOB approximated model.The background density gradient length and the adiabaticity is Ln = 16ρs0 and α = 0.005, respectively.

FIG. 6 .
FIG.6.The spatial pattern of the relative density fluctuation δn is shown during the nonlinear phase without (OHW) and with zonal flows (MHW), at time t = 4500/Ω0 and t = 9500/Ω0, respectively.Radial peaking of the turbulence intensity in the nonlinear regime proves to be a NOB effect.The background density gradient length and the adiabaticity is Ln = 16ρs0 and α = 0.005, respectively.