Beyond the Oberbeck-Boussinesq and long wavelength approximation

We present the first simulations of a reduced magnetized plasma model that incorporates both arbitrary wavelength polarization and non-Oberbeck-Boussinesq effects. Significant influence of these two effects on the density, electric potential and ExB vorticity and non-linear dynamics of interchange blobs are reported. Arbitrary wavelength polarization implicates so-called gyro-amplification that compared to a long wavelength approximation leads to highly amplified small-scale ExB vorticity fluctuations. These strongly increase the coherence and lifetime of blobs and alter the motion of the blobs through a slower blob-disintegration. Non-Oberbeck-Boussinesq effects incorporate plasma inertia, which substantially decreases the growth rate and linear acceleration of high amplitude blobs, while the maximum blob velocity is not affected. Finally, we generalize and numerically verify unified scaling laws for blob velocity, acceleration and growth rate that include both ion temperature and arbitrary blob amplitude dependence.


I. INTRODUCTION
The Oberbeck-Boussinesq and long wavelength approximation are two widely adopted simplifications in plasma theory with the aim to reduce the model's algebraic complexity and computational burden.
The Oberbeck-Boussinesq [1,2] approximation origins from neutral fluid dynamics and assumes that density variations only fully enter the buoyancy force and are everywhere else, i.e. the inertia terms, assumed constant or linearly dependent on the temperature difference. This could be interpreted as if only a thin layer of a fluid is considered where no large inhomogeneities occur, resulting in the alternative term thin layer approximation. In reduced plasma fluid theories the Oberbeck-Boussinesq approximation enters through the polarization charge (or polarization current) that represents inertia. The Oberbeck-Boussinesq approximation is naturally inherent to δF models. However, the validity of this assumption is questionable if large inhomogeneities appear. Such a large inhomogeneity occurs if large relative fluctuation amplitudes emerge or if their size compares to the perpendicular e-folding length L ⊥ so that k ⊥ L ⊥ ∼ 1. Large inhomogeneities can arise in a variety of magnetized plasma. They are typically observed in the edge and scrape-off layer of magnetically confined fusion plasmas [3][4][5][6][7][8][9][10][11][12][13]. The physics of the non-Oberbeck-Boussinesq regime is rich, enabling radial zonal flow advection [14], radially inhomogeneous zonal flows [15], blob-hole symmetry breaking [16] and decreased blob or increased hole acceleration [17].
The long wavelength approximation [18] assumes that scales much larger than the gyro-radius k ⊥ ρ 1 dom- * E-mail: markus.held@uit.no † E-mail: mattwi@fysik.dtu.dk inate the plasma dynamics. It typically neglects gyroradius effects above powers (k ⊥ ρ) 2 . This long wavelength ordering is inherent to reduced plasma descriptions, such as the drift-kinetic (or -fluid) description. However, it is also virtually always enforced for polarization effects in full-F gyro-kinetic or (-fluid) models, that are by construction k ⊥ ρ ∼ 1. For full-F gyro-kinetic models this simplification is reasoned in the computational cost, while for full-F gyro-fluid models theory and closures have been only extended recently [19]. However, turbulence in magnetized plasmas, e.g. driven by an ion temperature gradient, is active at k ⊥ ρ ∼ 1 [20,21] and the interplay between small and large scales determines the overall plasma dynamics and transport. In particular, electron temperature gradient driven modes that are active on the much smaller electron gyroradius interact with ion temperature gradient modes [22]. Further, order unity ion to electron temperature ratios appear [23][24][25], such that the ion gyro-radius is of the order of the drift scale. Thus k ⊥ ρ ∼ 1 turbulence influences e.g. also drift-wave and interchange turbulence [26]. This evidence suggests that the long wavelength approximation is inadmissible in turbulence models of magnetized plasmas. The application of the long wavelength approximation significantly reduces the lifetime and compactness of interchange blobs within the Oberbeck-Boussinesq regime [27]. Further, the long wavelength approximation crucially affects the normal modes of ion temperature gradient or trapped electron mode instabilities [28][29][30][31].
So far, the combined arbitrary wavelength and non-Oberbeck-Boussinesq regime is terra nova and it is unknown whether a synergy or superposition of known effects in the individual regimes exists. In the following, we present the first numerical investigation of a reduced magnetized plasma model that captures both regimes consistently. To this end, we utilize a recently developed full-F gyro-fluid model [19] that incorporates non-Oberbeck-Boussinesq and arbitrary wave-length polarization. Previous studies using gyro-fluid models were either limited to Oberbeck-Boussinesq [32] or long wavelength approximations [14,16,27,33]. Gyrokinetic studies relax the Oberbeck-Boussinesq approximation rarely consistently and then only in the long wavelength limit [34,35]. Attempts to include also arbitrary wavelength effects are based on an unphysical ad-hoc Padé-approximation that lacks a positive definite (e.g. quadratic) form of the kinetic energy [30,31,[36][37][38]. Studies based on drift-kinetic or -fluid models are inherently derived under the assumption k ⊥ ρ i 1 (see for example [39,40]) and are thus always in the longwavelength limit. We limit our study to interchange blob dynamics since they are a well known and understood benchmark test case and the latter combined regimes are of particular relevance for filamentary transport in the scrape-off layer of magnetically confined fusion devices. By means of theoretical estimates and numerical experiments that spawn a large parameter space we show that arbitrary wavelength polarization overlays with non-Oberbeck-Boussinesq effects in their contribution to the nonlinear blob dynamics. On the one hand, arbitrary wavelength effects introduce gyro-amplification of smallscale E×B vorticity. This manifests in strong small-scale E × B shear flows at the blob edge that result in highly coherent and long-living blobs. On the other hand, non-Oberbeck-Boussinesq effects are accompanied by plasma inertia that influences the acceleration and growth rate of the blobs with increasing blob amplitude.
The remainder of the manuscript is organized as follows. The full-F gyro-fluid model is presented in Section II, where we elaborate on the accuracy of the arbitrary wavelength polarization closure and for comparison also introduce and discuss the long-wavelength and Oberbeck-Boussinesq approximated models. Further general relative error estimates for the Oberbeck-Boussinesq and long wavelength limit are derived that allow to estimate the validity of the latter approximations. On top of that, so-called gyro-amplification is introduced, which boosts small-scale E × B vorticity fluctuations when arbitrary wavelength polarization is retained. At the end of this section the invariants of our model are introduced and the blob initial condition is stated. In Section III unified scaling laws for the blob center of mass velocity, acceleration and interchange growth rates are deduced from the full-F gyro-fluid model, its invariants and initial condition. Our numerical experiments in Section IV reinforce the fundamental difference in the nonlinear dynamics of cold and hot blobs. We there also expand in detail on the influence of arbitrary wavelength and non-Oberbeck-Boussinesq effects on the blob dynamics, shape, pattern and compactness. Further, the previously derived unified scaling laws are verified by our numerical experiments. Finally, we summarize our main findings in Section V.

II. FULL-F GYRO-FLUID MODEL
We consider a full-F gyro-fluid model, that is founded on the standard gyro-kinetic ordering [41] and fully encompasses effects down to the gyro-radius scale [19]. The gyro-fluid moment hierarchy is closed by truncation, in particular an isothermal assumption. For the sake of simplicity, we neglect parallel dynamics and electromagnetic effects. The full-F gyro-center continuity equations for a gyro-center species density N 1 is given by where the gyro-center E × B and ∇B drifts are with the gyro-average and polarization part of the gyrofluid potential The magnetic field points into z and its magnitude is Here we defined the reference radius R 0 , the initial position x 0 and the reference magnetic field magnitude B 0 . The magnetic field unit vector is defined byb := B/B = e z . It is utilized for the perpendicular projection of a vector h according to h ⊥ := −b × (b × h) = P · h, where P := g −bb and g is the projection and metric tensor, respectively. We also introduced the particle charge q, mass m, gyro-frequency Ω := qB/m and perpendicular gyro-center temperature T ⊥ . The finite Larmor radius (FLR) and polarization operators are included through Padé approximations [19] and are taken in this work in the ρ = const. limit. As a consequence these operators are self-adjoint. Here, we introduce the thermal gyro-radius and also define the the drift scale respectively. The utilized Padé approximations in the FLR and polarization operators (Eq. (4a) and (4b)) are excellent approximations to the exact operators where I 0 is the zeroth order modified Bessel function [19].
In particular, the Padé approximations are fully accurate up to O(k 2 ⊥ ρ 2 ) and closely mimic the behaviour at arbitrary wavelengths [19]. We emphasize that in the gyro-center continuity Eq. (1) only the compression of the E × B drift, ∇ · U E = − 1 B0R0 ∂ y (ψ 1 + ψ 2 ), is non-vanishing. This is important for the blob dynamics at small fluctuation amplitudes [17,42]. The full-F gyro-fluid model is completed by the Poisson equation with the polarization densities [19] Note that the use of Padé approximations allows us to write the polarization charges as simple divergence relations of the polarization densities (e.g. −∇ · P 1 = s q(Γ 1 − 1)N ≡ ∇ · s q∇ ⊥ Γ 1 ρ 2 N/2). In Eq. (10b) and (2a) we adopt the second order accurate Padé approximation of Ref. [19] that gives rise to the square root of the Padé approximated polarization operator √ Γ 0 . A square root operation could be avoided by using the fourth order accurate Padé approximation of Ref. [19]. However, this results in a fourth order elliptic equation instead of a second order elliptic equation in Eq. 9 (cf. Ref. [19]). Further, it does not resemble the commonly used second order accurate Padé approximation of Eq. (4b) of Oberbeck-Boussinesq approximated gyro-fluid models in the Oberbeck-Boussinesq limit. The herein utilized gyro-fluid model resorts to a constant thermal gyro-radius (ρ = const.) and a magnetic field aligned coordinate system (Cartesian coordinates andb =ê z ). These simplifications ease the analytical and numerical treatment of the arbitrary wavelength polarization closure. The ρ = const. approximation in the FLR and polarization operators (Eqs. (4a) and (4b)) permits to commute the gyro-radius with the perpendicular Laplacian, so that these operators are self-adjoint. The additional choice of a magnetic field aligned coordinate system allows us to commute the polarization operators through the divergence ∇ · √ Γ 0 f = √ Γ 0 ∇ · f and spatial derivatives In general this is permitted due a spatially dependence in the thermal gyro-radius ρ and in the generally non-orthogonal projection tensor P. However, we remark that cylindrical coordinates and the straight field line approach are commonly exploited for realistic three dimensional computations with toroidal magnetic fields, such as Tokamaks. For this approach the simplifications above can be readily exploited within the ρ = const. approximation. We emphasize that the arbitrary wavelength polarization closure is by construction only fully accurate up to O(k 2 ⊥ ρ 2 ) but resembles the Padé approximated Oberbeck-Boussinesq limit (Eq. (4b)) at arbitrary k ⊥ ρ. This enables us to bridge the gap between long wavelength approximated full-F and arbitrary wavelength Oberbeck-Boussinesq approximated full-F (or δF ) gyrofluid models. The numerical implementation and study of a polarization closure that is fully accurate or achieves better accuracy to arbitrary k ⊥ ρ is shifted to future work. In Sec. II A we elaborate more on the accuracy of the chosen arbitrary wavelength polarization closure. We remark that we do not consider the ad-hoc approximation ∇ · P 2 ≈ − s [30,31,[36][37][38] since it can not be derived from a particular choice of ψ 2 by means of field theory.
A. Accuracy of the arbitrary wavelength polarization closure

Exact arbitrary wavelength polarization
The exact expressions for the arbitrary wavelength polarization treatment with a near Maxwellian gyrocenter distribution function have been derived recently in Ref. [19], which provide the gyro-moment expressions in Fourier space as well as its Taylor series expansion in configuration space. In the following we list the arbitrary wavelength expressions of the polarization part of the basic gyro-fluid potential and its associated polarization charge density which results from variational calculus [19]. Without loss of generality in the following discussion we adhere to ρ = const.. Here, we defined the wave-vector K := k+k and the non-Oberbeck-Boussinesq form of the exact polarization operator First, we notice that Eqs. (13), i.e. Γ ex 0 (b, b ) − 1, is a convolution kernel of infinite rank due to the bracket e b ⊥ ·b ⊥ − I 0 (b ⊥ b ⊥ ) and thus is not separable. Second, solving the exact polarization part of the basic gyro-fluid potential (Eq. (11)) and the gyro-fluid Poisson Eq. (12) necessitates an efficient numerical solution of a convolution and deconvolution, respectively. By contrast a separable approximation to the exact convolution kernel allows to bypass these cumbersome convolutions within a treatment in configuration space. In Fig. 1 (top row) we show the exact polarization kernel for three meaningful phase angles θ. Here, the phase angle θ is the angle between The exact polarization kernel features a sharp finite tail for all presented phase angles around roughly b ⊥ ∼ b ⊥ . This sharp tail appears broadly around b 1 and extends narrowly to b → ∞. The polarization operator of Eq. (13) reduces to the conventional and exact polarization operator (Eq. (8)) in the Oberbeck-Boussinesq limit. There, Eq. (11) vanishes and we neglect the spatial dependence in the gyrofluid moment variables and the magnetic field magnitude in the polarization charge density, so that Γ ex The long wavelength limit expressions are derived by Taylor expanding the exact polariza-

Non-separable Padé approximation
Exploiting a symmetric bivariate Padé approximation of order (1,2) to each term of the exact operator (by expanding Γ ex 0 (tb ⊥ , tb ⊥ , θ) around t = 0 and setting t = 1 afterwards) yields Note that the bivariate Padé approximation is not separable, fully accurate to O(k 2 ⊥ ρ 2 ) and resembles Eq. (4b) in the Oberbeck-Boussinesq limit. In Fig. 1 (center row) it is shown that the sharp tail of the exact polarization kernel is well approximated to arbitrary wavelengths for parallel and anti-parallel wave-vectors (θ = 0 and θ = π). However for orthogonal wave-vectors the kernel vanishes and does not resemble the finite peak of the exact polarization kernel. Note that a bivariate Padé approximation of order (1,2) to the full expression yields Γ ex which is not separable and also less accurate. Further, a nested Padé approximation of order (1, 2) also produces a non-separable approximation, which is too complex to implement while being less accurate than Eq. (14).

Separable Padé approximation
The Padé approximation proposed in Ref. [19] is by construction separable, exactly accurate to order O(k 2n ⊥ ρ 2n ) and reduces in the Oberbeck-Boussinesq limit to the Padé approximated polarization operator. The most simple yet O(k 2 ⊥ ρ 2 ) accurate form of the Padé approximation is [19] Observe that the latter separable operator resembles the long wavelength limit and the Oberbeck-Boussinesq limit of Γ 0 . Further, we find the polarization relevant terms in configuration space (Eqs. (3b) and (10b) with Eq. (4b)) from an inverse Fourier transform of Eqs. (11) and (12) with Eq. (15). The separable Padé approximated polarization kernel is depicted in Fig. 1 (bottom row). We notice that the sharp region of the exact polarization kernel, i.e for b ⊥ ∼ b ⊥ 1, is not well captured by the separable Padé approximation. Similarly to the non-separable Padé approximated polarization kernel its separable equivalent vanishes for orthogonal wave-vectors.

B. Oberbeck-Boussinesq limit
The Oberbeck-Boussinesq approximated full-F gyrofluid model arises by assuming (i) small relative fluctuation amplitudes in the gyro-center density and perpendicular pressure, |δN | 1 and |δP ⊥ | 1, (ii) large density and perpendicular pressure gradient lengths, k ⊥ N/|∇N | 1 and k ⊥ P ⊥ /|∇P ⊥ | 1, (iii) small spatial variations in the magnetic field magnitude 2 and (iv) a vanishing second order polarization term (ψ 2 = 0) in the gyro-fluid E × B drift of Eq. (2a) to restore energetic consistency. Here, we introduced the relative fluctuation δQ := Q/Q 0 − 1 of a quantity Q with respect to its stationary part Q 0 . These assumptions result in a Oberbeck-Boussinesq approximated full-F gyrofluid model that is similar to the δF gyro-fluid model. We emphasize that the Oberbeck-Boussinesq approximation necessitates the splitting of the gyro-center distribution function into a stationary and relative fluctuation part [19]. This stationary part of the gyro-center distribution function then appears as a gyro-moment quantity in the polarization density. In general this stationary state is not known a priori or there is no stationary state at all. For the first case the stationary state must be determined by a (non-Oberbeck-Boussinesq or full-F ) model that avoids such a separation. For the second case an approximate stationary state must be assumed.  14) and (15), respectively) are shown for three different phase angles θ = (0, π/2, π). In contrast to the exact polarization kernel its Padé approximations vanish for orthogonal wave-vectors θ = π/2. Further the separable Padé approximation is not able to capture the sharpening of the finite tail of the exact and non-separable Padé approximation towards arbitrary wavelengths for parallel (θ = 0) and anti-parallel (θ = π) wave-vectors.
The Oberbeck-Boussinesq approximated full-F gyro-fluid model consists of the gyro-center continuity equations Note that we do not drop the E ×B compression term in Eq. (16), in contrast to the original Oberbeck-Boussinesq approximation.

C. Long wavelength limit
In the long wavelength limit we assume in the polarization terms of Eq. (10b) and (3b), which yields We introduced the E × B vorticity where u E = (b × ∇φ)/B. The long wavelength approximation of Eq. (18) is in line with previous full-F gyro-fluid studies that equally adopt the long wavelength limit only in the 2nd order polarization terms [14-16, 27, 33, 43]. By contrast, a full-F drift-fluid formalism is recovered if the long wavelength approximation is additionally applied to the 1st order polarization (also known as diamagnetic) terms of Eq. (10a) and (3a) together with replacing the gyro-center density N by a particle density n in Eq. (18).

D. Long wavelength and Oberbeck-Boussinesq limit
For comparison also the long wavelength limit of the Oberbeck-Boussinesq approximated full-F gyro-fluid model is considered, where the polarization terms of Eq. (10b) and (3b) are approximated according For the sake of clarity, we summarize the various treatments of the polarization terms that appear in the original full-F gyro-fluid model, as well as its three approximations in Table I.

E. Gyro-amplification
Gyro-amplification refers to the increase of E × B vorticity when retaining arbitrary perpendicular wavelength polarization. This can be understood from solving the Oberbeck-Boussinesq approximated Poisson equation with Eq. (17) for the E × B vorticity Ω E ∼ ∆ ⊥ φ. This results in an additional factor Ω E,k ∼ (1+ρ 2 k 2 ⊥ ) for arbitrary wavelength polarization, stemming from the use of the polarization operator Γ 0 (4b). Thus, E × B vorticity structures at and below the drift scale ρ s0 are greatly enhanced for typical values of T i⊥ /T e⊥ = 1, which is depicted in Fig. 2. We will study this behaviour further in Section IV, where we will find that the amplification is typically around a factor 10 at the dominant scale in comparison to the long wavelength treatment. In order to quantify the applicability of the long wavelength and Oberbeck-Boussinesq approximation rigorously we introduce relative errors for the respective regime. These relative errors are deduced from the particular approximation in the Poisson Eq. (9) by expanding up to O(ρ 2 k 2 ⊥ ). Additionally, spatial variations in the magnetic field magnitude are assumed to be small. According to this principle, the relative errors of the long wavelength and Oberbeck-Boussinesq approximation are derived to where f 2 denotes the 2-norm. The relative error LW L depends on the inverse perpendicular temperature gradi-ent length, the thermal gyro-radius and the spatial scale of the electric potential (and thus the E × B vorticity). The relative error LW L is large if In the Oberbeck-Boussinesq limit LW L of Eq. (21a) recovers LW L|OB of Eq. (22a). An error of order unity is approached if one of the following conditions is fulfilled: In the long wavelength limit OB of Eq. (21b) recovers OB|LW L of Eq. (22b). Analogously, we define the relative error for long wavelength approximation under the condition of the Oberbeck-Boussinesq approximation and vice versa by Note that LW L|OB depends on the thermal gyro-radius and the spatial scale of the electric potential (and consequently the E × B vorticity). More specifically, the long wavelength approximation is questionable if The relative error OB|LW L depends on both the inverse density gradient length and the relative density fluctuation amplitude. The Oberbeck-Boussinesq approxima-

G. Conserved quantities
In the absence of dissipation and surface integral terms the presented gyro-fluid models possess a number of conserved quantities. These are the gyro-center total particle number M , total polarization charge Q, Helmholtz free energy F and mechanical energy Z where here and in the following we integrate over the entire domain. Here, we defined the kinetic energy E k , 3 More specifically, |δn/(1 + δn)| ∼ 1/2 so that e.g. for blobs and holes δn ∼ 1 and δn ∼ −1/3, respectively.
the entropy S and the potential energy H The kinetic energy and the entropy obey E k (t) ≥ 0 and S(t) ≤ 0, respectively. Note that in the Oberbeck-Boussinesq limit Eqs. (24b) and (24c) reduce to by Taylor expanding for small fluctuation amplitudes or small magnetic field variations, respectively.

H. Initialization
We will consider two initial conditions that differ in the initial ion gyro-center density N i (x, 0) and as a consequence the initial electric potential φ(x, 0). However, the electron density field n e (x, 0) is always initialized by a Gaussian with size σ, amplitude ∆n e and initial position x 0 on top of a constant background n e0 . The initial electron density of Eq. (26) determines the initial total particle number M e (0) = 2πσ 2 ∆n e and the initial electron entropy S e (0) = −2πσ 2 f (∆n e ). Here, we defined the function f (x) : where Li n (x) is the polylogarithm function. For small relative fluctuation amplitudes the initial electron entropy recovers its δF limit S e (0) ≈ −π/2σ 2 ∆n 2 e . The initial electron potential energies is H e (0) ≈ 0, where we used again the Taylor expanded integral.

Non-rotating Gaussian
For the first initial condition the initial electric field vanishes (∇ ⊥ φ(x, 0) = 0) so that the initial ion gyrocenter densities fulfill the condition n e (x, 0) = Γ 1,i N i (x, 0).
The initial ion entropy and potential energy is S i (0) ≈ S e (0) and H i (0) ≈ 0, respectively. Further, the initial ion kinetic energy is vanishing E k,i (0) = 0.

Rotating Gaussian
For the second initial condition the total polarization density is constant P 1 (x, 0) + P 2 (x, 0) = const., which gives rise to an initial electric field for finite ion temperature. The constant polarization density condition is fullfilled by In the long wavelength limit we obtain the initial E × B rotation is ∇ ⊥ φ ≈ T i⊥ /(2q i )∇ ⊥ ln n + const.. Thus the sum of the E × B and diamagnetic vorticity density is non-vanishing and no force balance is established initially [33,44]. From Eq. (28) follows that the initial ion entropy and potential energy is given by S i (0) = S e (0) and H i (0) = H e (0), respectively. Further, the initial ion kinetic en- πg(∆n e ), where we used the long wavelength limit expression, ∇ ⊥ φ = T i⊥ /(2q i )∇ ⊥ ln N i , assumed small variations in the magnetic field magnitude B −2 ≈ B −2 0 and defined g(x) := (6x − π 2 + 3 ln(1 + x) ln((1 + x)/x 2 ) + 6Li 2 (1/(1 + x)))/3. For small fluctuation amplitudes the initial ion kinetic energy quadratically depends on the relative fluctuation Note that the initial entropies are much larger than the initial ion kinetic energy, s T ⊥ S(0) −E k,i (0), as long as . This condition is fulfilled for all parameters that we consider. Consequently, we neglect the contribution of the initial ion kinetic energy in our further analysis.

III. UNIFIED SCALING LAWS FOR INTERCHANGE BLOB DYNAMICS
We define three measures for the interchange blob dynamics. These are the center of mass position, velocity and acceleration The x-component of the center of mass position is related to the potential energy of the electrons via H e (t) ≈ −M e T e (X − x 0 )/R 0 . Thus, the potential energy is decreasing in time s H(t) ≤ s H(0) for blobs (and holes). This allows us to show that the kinetic energy as well as the entropy is increasing in time s E k (t) ≥ s E k (0) and s T ⊥ S(t) ≥ s T ⊥ S(0), respectively. Thus, the second law of thermodynamics is fulfilled. In the following we will deduce scaling laws for the xcomponent of the center of mass acceleration A x and velocity V x as well as for the interchange growth rate γ. We base these scaling laws on the following estimate for the squared center of mass momentum [17] ( where we use the Cauchy-Schwartz inequality and the estimate dAn e (B −1 ∂ y φ) 2 ≤ 2E k,i (t)/m i . Further, we introduce the integral D(t) := dA(n e − n e0 ) 2 /n e . For blobs (n e ≥ n e0 ) there exist two close upper bounds for the expression inside the latter integral, which are depicted in Fig. 3. These are the (scaled) electron entropy density 2 [n e ln(n e /n e0 ) − (n e − n e0 )] for n e0 ≤ n e < n e,c and the relative density n e − n e0 for n e ≥ n e,c . Here, we defined the critical density n e,c := {1 + exp [3/2 + W (−3/2/ exp (3/2))]} n e0 ≈ 2.397n e0 with the product logarithm W (x), which follows from equating the scaled electron entropy and relative density. Thus, we can further assess the integral D(t) to where we used the identities S e (t) ≥ S e (0) and M e (t) = M e (0).

A. Acceleration
with the ion acoustic speed c s := (T e⊥ + T i⊥ )/m i and its cold ion limit equivalent c s0 := T e⊥ /m i . Assuming now a linear acceleration V x = A x t and X(t) = A x t 2 /2 in Eq. (32) results into an estimate for the acceleration [17,42] Here, we introduced the parameter Q ∈ (0, 1] and used the (1,1) Padé approximation for Se(0) Me ≈ − 1 4 ∆ne ne0+2/9∆ne . Analogously, we obtain with the second estimate of Eq (31) the upper bound for the acceleration A x ≤ c 2 s /R 0 . This upper bound can be interpreted as an effective gravity g eff := c 2 s /R 0 that would arise when replacing the ∇B drift with a gravitational drift.
B. Velocity

Linear velocity scaling
From the conservation of energy (Eq. (23c)) results that the kinetic energy of the ions is bounded by . This together with the first estimate of Eq (31) in Eq. (30) yields the linear velocity scaling [17,42] where the Padé approximation is Taylor expanded for small relative density fluctuations.

Square root velocity scaling
The square root scaling can be deduced by taking the time derivative of the full-F Poisson Eq. (9) and applying the long wavelength approximation [33,45] where we introduced the vorticity density W :=b · ∇ × j = ∇ · ω and the vector ω := −b × j with the sum of the E × B and ion diamagnetic current density j := n e (u E + u D,i ) and the ion diamagnetic drift . A scale analysis of Eq. (35) yields the square root velocity scaling [32,33,46,47] with parameter R ∈ (0, 1]. Interestingly, a square root scaling results also from similar methods as we used for the linear velocity scaling law (Eq. (34)). To this end, we use the same bound E k, together with the second estimate of Eq (31) in Eq. (30). This yields max is estimated through the scale analysis expression of Eq. (36). This scaling law agrees in the Oberbeck-Boussinesq limit with Eq. (36). However, in the cold ion limit the numerical data clearly supports the square root velocity scaling law of Eq. (36) and does not justify the inclusion of the non-Oberbeck-Boussinesq factor.

Unified velocity scaling
Following similar methods as in Ref. [17,48] we can combine the linear and square root velocity scaling law into a unified scaling law that accounts for finite ion temperatures The latter resembles the small and large relative density fluctuation limits, Eq. (34) respectively Eq. (36), correctly. Note that the only difference to the proposed unified scaling law in Ref. [17] is the factor 1 + T i⊥ /T e⊥ that is inherent in the ion acoustic speed c s .

C. Interchange growth rate
The time at which the maximum velocity occurs is estimated by the inverse interchange growth rate t max(Vx) ∼ γ −1 , which we define as

IV. NUMERICAL EXPERIMENTS
The full-F gyro-fluid model and its three approximations (long wavelength, Oberbeck-Boussinesq, long wavelength + Oberbeck-Boussinesq) is numerically solved by using the open source library Feltor [49,50]. We choose a discontinuous Galerkin discretization on a rectangular grid in space and use an explicit adaptive timestepper based on the Bogacki-Shampine embedded Runge-Kutta method in time. The convective terms in Eq. (1) are thus discretized using a discontinuous Galerkin upwind scheme [51], while the elliptic polarization equation (9) is discretized using local discontinuous Galerkin methods [52]. The advantages of these methods are their high order, low numerical diffusion and ease of parallelization. We thus take advantage of Nvidia's V100 GPUs on the Marconi M100 supercomputer for all our simulations.
We remark that the numerical application of the √ Γ 0 operation demands to efficiently solve a matrix function equation of the form √ Ax = b, where A is the discrete form of 1 − ρ 2 ∆ ⊥ in configuration space and the matrix function is the matrix square root. Our approach is based on a Krylov approximation in order to drastically reduce the original matrix size for the matrix function (i.e. square root function) computation by projecting onto a Krylov subspace. More specifically, we use the symmetric Lanczos algorithm to obtain a matrix decomposition of A in tridiagonal form. The matrix function is then applied to the small tridiagonal matrix in terms of either an eigenvalue decomposition or a Cauchy-contour integral [53]. This is by far more efficient than for example an eigenvalue decomposition of the original large sparse matrix A. The remaining problem is the construction of a judicious stopping criterion for the Lanczos algorithm. For this the recently proposed solution of Ref. [54] is utilized. We remark that all of our implementations are freely available [50] and well-documented. We employ unit and integration tests and we observe the expected order of convergence in all manufactured test problems. An in depth discussion of the implementation and application of matrix-function computation will be given in a separate manuscript. The supplemental dataset to this contribution ensures bitwise reproducibility of the herein reported results [55].
We resolve all herein presented simulations with n = 5 polynomial coefficients and N x = N y = 300 grid cells with a box size of L x = L y = 40σ. In total this amounts to approximately 2 million grid points. This high order and fine resolution is necessary to resolve very small scale E × B vorticity structures in the numerical solutions as will become evident in this section. Boundary conditions are Dirichlet in x-direction and periodic in y-direction. For the sake of regularizing the resulting 5-th order upwind discretization we add a very small hyperdiffusion term of 2nd order to the continuity equations. The mass diffusion coefficient ν is determined through the Rayleigh number Ra := g eff σ 3 ∆ne ν 2 ne0 , which is fixed to Ra = 2 × 10 9 . The Schmidt number is set to unity so that the viscosity coefficient equals the mass diffusion coefficient. As a consequence, the chosen parameters span the turbulent regime [47,56]. The chosen physical parameter space lies in a typical experimental regime, but also in a regime where the ap-proximations in question are expected to fail. This is for small blobs with blob size σ = 5ρ s0 , hot ion temperature T i⊥0 /T e⊥0 = 4 and density perturbations that encompass values in the range ∆n e /n e0 ∈ 10 −3 , 10 . The parameter for the magnetic field gradient is chosen to ρs0 R0 = 0.00015 to allow direct comparison to previous studies [33].

A. Fundamental dynamics of cold and hot blobs
It is well known that the finite ion temperature effects lead to fundamentally different blob dynamics in comparison to cold ion temperature [27,32,33], which are visualized in Fig. 4. The respective movies ("Full-F, amp=1, tau=0", "Full-F, amp=1, tau=4", "Full-F, amp=1, tau=4, 0pol") to this figure can be found in the supplementary material. In the cold ion limit a nearly up-down symmetric mushroom like blob develops that is advected purely into the x-direction [47]. The mushroom like structure consists of a steepening blob front and lobes that roll itself up. This structure quickly disintegrates after the initial linear acceleration phase. The x-directed E × B propagation is a consequence of a dipole in the electric field that emerges due to the interplay of polarization and the ∇B drift. By contrast, for finite ion temperature the blob strongly retains its initial shape both for the non-rotating and rotating Gaussian initial condition. The tendency to maintain its initial shape is attributed to small scale E × B shear flows 4 that suppress the removal of mass via small eddy-satellites. These shear flows can be attributed to the gyro-amplification discussed in Section II E. The initial dipole in the vorticity rolls up into a spiral of strongly sheared flows. These are strongest at the edge of the blob and decrease in magnitude towards the blob center. The small scale blob E ×B shear flows share features of zonal flows of magnetized plasma turbulence [14]. The emerging blob shape is most reminiscent of a jellyfish, in particular for the rotating Gaussian initial condition [57]. The additional blob motion in y-direction strongly depends on the initial condition. In particular the blob is propagating also into the negative y-direction for the non-rotating initial condition, while the y-motion mostly disappears for an initially rotating Gaussian.

B. Arbitrary wavelength polarization and non-Oberbeck-Boussinesq effects on nonlinear blob dynamics
We now turn our attention to the differences between the considered models (cf. Table I) for various initial FIG. 4. The evolution of the electron density ne and the E × B vorticity is shown for two different ion to electron temperature ratios T i⊥ /T e⊥ = (0.0, 4.0) and two different initial conditions for the electric field. The initial blob amplitude and blob size are ∆ne/ne0 = 1.0 and σ = 5ρs0, respectively. The respective movies ("Full-F, amp=1, tau=0", "Full-F, amp=1, tau=4", "Full-F, amp=1, tau=4, 0pol") to this figure can be found in the supplementary material. blob amplitudes and for the non-rotating Gaussian initial condition (Eqs (26) and (27)). In Figures 5, 6 7 the temporal blob evolution is shown for three typical blob amplitudes ∆n e /n e0 = (5.0, 1.0, 0.1). The respective movies ("Full-F, amp=5, tau=4", "OB, amp=5, tau=4", "LWL, amp=5, tau=4", "OB+LWL, amp=5, tau=4", "Full-F, amp=1, tau=4", "OB, amp=1, tau=4", "LWL, amp=1, tau=4", "OB+LWL, amp=1, tau=4", "Full-F, amp=0.1, tau=4", "OB, amp=0.1, tau=4", "LWL, amp=0.1, tau=4", "OB+LWL, amp=0.1, tau=4") that highlight the differences between the considered models can be found in the supplementary material. For all models we observe the initial rolling up or spiraling in the E × B vorticity Ω E :=b · ∇ × u E and the consequent x-and y-directed blob motion [27,32,33]. The amplitude of the E × B shear flows decreases by roughly an order of magnitude if the long wavelength approximation is applied. Further, for the long wavelength approximated models the E × B shear flows are at much larger scales and no longer show up a clear gradient in magnitude from the blob center towards the blob edge. As a consequence, the long wavelength approximation leads to less coherent blobs that disintegrate more quickly. The reduction in blob coherence is most pronounced for small amplitudes and weakens with increasing initial blob am-plitude. Apart from the effect of blob shape and structure the blob propagation is also significantly affected by the considered approximations. The Oberbeck-Boussinesq approximation leads to an increased blob displacement for very large blob amplitudes (∆n e /n e0 = (5.0, 1.0) in Fig 5  and 6). This is best visible at the end of the linear acceleration phase (roughly coinciding with the second time point in Fig 5 and 6). After the linear acceleration phase the long wavelength approximation results in different blob positions due to the dissociation of smaller eddy-satellites that can change the movement of the blob or due to complete disintegration of the blob. This is most clearly recognizable for the small amplitude blob (∆n e /n e0 = 0.1 in Fig 7). Finally, we remark that the rotating Gaussian initial condition is not included in this analysis, since in this case one of the physical initial fields, specifically the electric potential φ(x, 0), changes if the long wavelength or Oberbeck-Boussinesq approximation is applied. However, we observed that the long wavelength approximation results in less coherent blobs in comparison to the non-rotating Gaussian initial condition. This underpins the results that are obtained for the non-rotating Gaussian intial condition.  Table I). Note the change in vorticity pattern and magnitude. The respective movies ("Full-F, amp=5, tau=4", "OB, amp=5, tau=4", "LWL, amp=5, tau=4", "OB+LWL, amp=5, tau=4") to this figure can be found in the supplementary material.

C. Blob compactness
The blobs ability to retain its initial (Gaussian) shape is quantified by the blob compactness Here, we introduced the Heaviside function  Table I). Note the change in vorticity pattern and magnitude. The respective movies ("Full-F, amp=0.1, tau=4", "OB, amp=0.1, tau=4", "LWL, amp=0.1, tau=4", "OB+LWL, amp=0.1, tau=4") to this figure can be found in the supplementary material. and the position of the maximum electron density X max (t).
In Figure 8 the time evolution of the blob compactness I c (t) is shown for three different initial blob amplitudes ∆n e /n e0 = (5.0, 1.0, 0.1). Large initial amplitude blobs decrease their compactness slower than their small initial amplitude counterparts. Note that a slower decrease in blob compactness results in an increase in the bloblifetime and likewise blob-coherence. The blob compactness is significantly reduced if the long wavelength approximation is applied, while the Oberbeck-Boussinesq approximation only slightly increases the compactness for the large initial amplitude. In Figure 9 the blob compactness I c (t) is shown at time t = 3/γ as a function of the initial blob amplitude ∆n e /n e0 for the complete studied parameter space. A transition from low to high compactness takes place between roughly ∆n e /n e0 ≈ 0.1. Clearly, large initial blob amplitudes lead to more coherent blobs than small initial blob amplitudes. The long wavelength approximation results in significant reduction in compactness within the transition region. For very small blob amplitudes no significant deviations in the compactness appear if the long wavelength approximation is made, since in both cases the blobs have already largely disintegrated at t = 3/γ. On the other hand, for large initial blob amplitudes the disintegration of the long wavelength approximated has not been initiated at t = 3γ. Thus, in the large initial amplitude limit the long wavelength approximation introduces only a slight decrease in compactness at that time point. The increase in coherency is attributed to  (Table I) and the nonrotating Gaussian initial conditions. the increase in E × B vorticity when retaining arbitrary perpendicular wavelength polarization, which we introduced as gyro-amplification in Section II E.

D. Verification of unified scaling laws for hot blobs
In the following we verify the derived unified scaling laws of Section III. We do not attempt to verify any FIG. 9. The dependence of the blob compactness Ic at t = 3/γ on the initial density amplitude is shown for all considered models (Table I).
scaling laws for the y-directed (or total) center of mass dynamics [33], since this motion can depend strongly on e.g. the initial condition or on superimposed background flows of the plasma. The fitting constants (R, Q) of the unified blob scaling laws of Eqs. (37), (33) and (38) can be already determined in the cold ion limit [17]. Previously, these constants have been determined by the best fit of the acceleration and velocity scaling to (R, Q) = (0.85, 0.32), while neglecting the fit on the growth rate. We here present improved fitting constants (R, Q) = (0.95, 0.31) determined by a best fit to all three scaling laws (velocity, acceleration and growth rate). In Figure 10 we verify the velocity scaling of Eq. (37) by the measured maximum of the center of mass velocities for varying initial blob amplitudes and for all considered models. The unified velocity of Eq. (37) accurately  (Table I). The line represents the unified velocity scaling of Eq. (37). The shading marks a deviation by 25%.
captures the behaviour for all amplitudes within a rela-tive error of 25%. The highest accuracy is achieved for very small amplitudes. For increasing amplitudes the unified velocity scaling of Eq. (37) slightly overestimates the center of mass blob velocities and the relative error approaches roughly 25%. No significant deviations due to the long wavelength or Oberbeck-Boussinesq approximation appear in the maximum of the center of mass velocities. In Figure 11 we compare the acceleration scaling of Eq. (33) with the measured linear center of mass acceleration max(V x )t −1 max(Vx) for varying blob amplitudes and for all considered models. We emphasize here that we do not take the absolute maximum value of the measured center of mass acceleration since our scaling theory assumes a linear acceleration (cf. Sec. III). In agreement with previous studies we find that the linear acceleration is increased by up to a factor two for increasing amplitudes if the Oberbeck-Boussinesq approximation is utilized. Further, the acceleration scaling estimate of Eq. (33) matches the measured linear center of mass acceleration excellently for small amplitudes. For large amplitudes we find a qualitative agreement.
FIG. 11. The measured average radial center of mass acceleration as a function of the blob amplitude is depicted for all considered models (Table I) and both initial conditions. The solid and dashed line represents the acceleration scaling estimate of Eq. (33) and its Oberbeck-Boussinesq limit, respectively. The shadings indicate a relative error of 25%.
In Figure 12 we show the growth rate scaling of Eq. (38) and the measured growth rate, t −1 max(Vx) , for varying blob amplitudes and for all considered models. The growth rate estimate of Eq. (38) matches the data in the non-Oberbeck-Boussinesq regime by up to 25%. In the Oberbeck-Boussinesq limit the agreement is slightly above the 25% for large amplitudes, but resembles very well the respective limit of the growth rate scaling law. Note that we take always the first local maximum of the center of mass velocity and not the total maximum of the center of mass velocity, which could occur at later times for large blob amplitudes (cf. [33]). This is because our scaling theory assumes a linear acceleration phase that is only approximately fulfilled up to the first maximum FIG. 12. The measured average growth rate as a function of the blob amplitude is depicted for all considered models (Table I) and both initial conditions. The solid and dashed line represents the growth rate scaling estimate of Eq. (38) and its Oberbeck-Boussinesq limit, respectively. The shadings indicate a deviation by 25%.
of the center of mass velocity. However, if we take the total maximum of the center of mass velocity we obtain the same agreement of 25% for the velocity scaling, but weaker agreement for large amplitudes for the acceleration and growth rate scaling.

V. DISCUSSION
This paper explores for the first time the regime beyond both the long wavelength and Oberbeck-Boussinesq approximation. This previously unresolved regime is of particular importance to filamentary transport in the scrape-off layer of magnetically confined fusion devices, but also of general importance to turbulence and structure formation in magnetized plasmas. Our study is enabled by a full-F gyro-fluid model that exploits a recently developed arbitrary wavelength polarization closure [19].
Most importantly, we find that the inclusion of polarization down to the thermal gyro-radius scale leads to highly coherent blob structures in the presence of a substantial background ion temperature. The long wavelength approximation significantly reduces the coherence and lifetime of blobs due to the neglect of gyroamplification. As a consequence, it modifies the motion of the blobs due to a faster disintegration of the blobs. The Oberbeck-Boussinesq approximation affects the propagation of the blobs by increased linear acceleration and growth rate at large initial blob amplitudes.
The blobs center of mass motion in the linear acceleration phase is very well captured by unified velocity, acceleration and growth rate scaling laws (Eq. (37), (33) and (38)) that we generalized to finite ion background temperature and that hold for arbitrary initial blob amplitudes.
We hypothesize that similarly to Ref. [33] the characteristic footprint of finite ion gyro-radius effects on the blob dynamics is amplified if a density perturbation is accompanied by a temperature perturbation. The numerical implementation and study of an arbitrary wavelength polarization closure with more advanced full-F gyro-fluid moment hierarchies [19], e.g.
including temperature dynamics, parallel dynamics and electromagnetic effects, is ongoing. This effort aims to assess the validity of the long wavelength and Oberbeck-Boussinesq approximation also on turbulence in magnetized plasmas. However, the herein presented results on the permissibility of these approximation questions current efforts in edge and scrape-off layer modeling of fusion plasma that are based on at least one of this approximations.