Pad\'e-based arbitrary wavelength polarization closures for full-F gyro-kinetic and -fluid models

TL;DR

This paper introduces Padé-based polarization closures for full-F gyro-kinetic and gyro-fluid models, improving accuracy and energy conservation across arbitrary wavelengths, addressing a key limitation in existing models.

Contribution

It presents a novel quadratic form and Padé approximations for polarization closures that enhance accuracy and energy conservation in gyro-fluid models.

Findings

Achieves polarization charge density with desired accuracy

Retains linear polarization effects at arbitrary wavelengths

Ensures proper energy conservation and correct long wavelength limits

Abstract

We propose a solution to the long-standing short wavelength polarization closure shortfall of full-F gyro-fluid models. This is achieved by first finding an appropriate quadratic form of the gyro-fluid moment over the polarization part of the gyro-center Hamiltonian. Secondly, we deduce Pad\'e-based approximations to the latter expression that produce a polarization charge density with the desired order of accuracy and retain linear polarization effects for arbitrary wavelengths. The proposed closures feature proper energy conservation and the anticipated Oberbeck-Boussinesq and long wavelength limits.

Tables

Table 1. Table 1: Padé-approximations for FLR operators

$q$	${}_q{}^1\check{\Gamma }_1^{}$	${}_q{}^1\check{\Gamma }_1^{†}$	${}_q{}^1\check{\Gamma }_2^{}$	${}_q{}^1\check{\Gamma }_3^{}$
$2$	$\frac{1}{\sqrt{1-{\rho }^2{\Delta }_{⟂}}}$	$\frac{1}{\sqrt{1-{\Delta }_{⟂}{\rho }^2}}$	$\frac{{\rho }^2}{2}\frac{1}{\sqrt{{\sum }_{i=0}^3\left(\frac{3}{i}\right){\left(-{\rho }^2{\Delta }_{⟂}\right)}^i}}{\Delta }_{⟂}$	${\rho }^2\frac{1}{\sqrt{{\sum }_{i=0}^5\left(\frac{5}{i}\right){\left(-{\rho }^2{\Delta }_{⟂}\right)}^i}}{\Delta }_{⟂}-\frac{{\rho }^4}{4}\frac{1}{\sqrt{{\sum }_{i=0}^5\left(\frac{5}{i}\right){\left(-{\rho }^2{\Delta }_{⟂}\right)}^i}}{\Delta }_{⟂}^2$
$4$	$\frac{1}{1-{\rho }^2/2{\Delta }_{⟂}}$	$\frac{1}{1-{\Delta }_{⟂}{\rho }^2/2}$	$\frac{{\rho }^2}{2}\frac{1}{{\sum }_{i=0}^2\left(\frac{2}{i}\right){\left(-{\rho }^2/2{\Delta }_{⟂}\right)}^i}{\Delta }_{⟂}$	${\rho }^2\frac{1}{{\sum }_{i=0}^3\left(\frac{3}{i}\right){\left(-{\rho }^2/2{\Delta }_{⟂}\right)}^i}{\Delta }_{⟂}$

Equations

$$S:=$$

$$\gamma$$

$$H$$

$${\Psi }_1$$

$$⟨f⟩(\mathbf{X},v_{\parallel },\mu ,t)$$

$${\text{vol}}_{\mathbf{Z}}:=\sqrt{\text{det}\left(\omega \right)}{d}^{6}Z=m^2|B_{\parallel }^{\ast }|\sqrt{g}{d}^{3}Xd{v}_{\parallel }d\mu d\theta ,$$

$${\text{vol}}_{\mathbf{X}}$$

$$\dot{\mathbf{X}}$$

$${\dot{v}}_{\parallel }$$

$$\dot{\mu }$$

$$\dot{\theta }$$

$$\frac{\partial }{\partial t}\left(BF\right)+\nabla \cdot \left(\dot{\mathbf{X}}BF\right)+\frac{\partial }{\partial v_{\parallel }}\left({\dot{v}}_{\parallel }BF\right)$$

$$\sum _{\alpha }\frac{\delta }{\delta \varphi (\mathbf{x},t)}\int {\text{vol}}_{\mathbf{Z}}Fq\Psi$$

$$\sum _{\alpha }\int \frac{{\text{vol}}_{\mathbf{p}}}{B}\left[q⟨BF{⟩}^{†}+q^2\left(\varphi ⟨{\partial }_{\mu }F{⟩}^{†}-⟨⟨\varphi ⟩{\partial }_{\mu }F{⟩}^{†}\right)\right]$$

$$⟨f{⟩}^{†}(\mathbf{x},v_{\parallel },t)$$

$$\displaystyle\big{\|}\zeta\big{\|}$$

$$N$$

$$P_{\perp }$$

$$Q_{\perp }$$

$$\displaystyle\frac{\partial}{\partial t}\big{\|}\zeta\big{\|}+{\boldsymbol{{\nabla}}}\cdot\big{\|}\zeta\dot{{{\boldsymbol{X}}}}\big{\|}$$

$$\frac{\partial }{\partial t}$$

$$\displaystyle+\frac{B{\boldsymbol{\mathcal{K}}}_{{\boldsymbol{{\nabla}}}B}}{q}\big{\|}{\mu}^{k+1}{v}_{\parallel}^{l}\big{\|}+\frac{\hat{{\boldsymbol{b}}}}{B}\times\big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|}\bigg{)}$$

$$=$$

$$\displaystyle+{\boldsymbol{\mathcal{K}}}_{{\boldsymbol{\kappa}}}\cdot\big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|}+\frac{{\boldsymbol{\mathcal{K}}}_{{\boldsymbol{\kappa}}}\cdot{\boldsymbol{{\nabla}}}B}{q}\big{\|}{\mu}^{k+1}{v}_{\parallel}^{l}\big{\|}\bigg{)}.$$

$$R_{\perp }$$

$$R_{\parallel }$$

$$S_{\perp }$$

$$\sum _{\alpha }qN-\nabla \cdot \mathbf{P}$$

$$\displaystyle{\boldsymbol{P}}_{1}:=\sum_{\alpha}{\boldsymbol{\nabla}}_{\perp}\int\text{vol}_{{{\boldsymbol{X}}}}q\frac{\delta\big{\|}\phi-\Psi_{1}\big{\|}}{\delta(\Delta_{\perp}\phi({\boldsymbol{x}},t))}$$

$${\mathbf{P}}_2$$

$$\displaystyle\big{\|}\Psi_{1}\big{\|}=$$

$$\displaystyle\big{\|}\Psi_{2}\big{\|}=$$

$$\displaystyle+4{\boldsymbol{{\nabla}}}_{\perp}\phi\cdot{\Delta}_{\perp}{\boldsymbol{{\nabla}}}_{\perp}\phi-\left({\Delta}_{\perp}\phi\right)^{2}\Big{]}$$

$$\displaystyle-\frac{qR_{\perp}}{48m^{2}\Omega^{4}}\Big{[}2\left|{\boldsymbol{{\nabla}}}_{\perp}^{3}\phi\right|^{2}+3{\boldsymbol{{\nabla}}}_{\perp}\phi\cdot{\Delta}_{\perp}^{2}{\boldsymbol{{\nabla}}}_{\perp}\phi$$

$$+3{\left|{\Delta }_{\perp }{\nabla }_{\perp }\varphi \right|}^2+6{\nabla }_{\perp }^2\varphi :{\Delta }_{\perp }{\nabla }_{\perp }^2\varphi$$

$$\displaystyle-3{\Delta}_{\perp}\phi{\Delta}_{\perp}^{2}\phi\Big{]}+\dots$$

$$F=F_M\sum _{i=0}^{\infty }\sum _{j=0}^{\infty }{\mathcal{N}}_{ij}{L}_i(c_{\perp })H_j(c_{\parallel }),$$

$$F_M$$

Full text

I Introduction:

Gyro-fluid models are an extremely useful approach to provide insights into the behavior of magnetized plasmas. They are widely applied to study turbulent transport in the tokamak core Waltz et al. (1997); Staebler, Kinsey, and Waltz (2005), edge Scott (2000); Held, Wiesenberger, and Kendl (2019) and scrape-off layer Madsen et al. (2011); Wiesenberger, Madsen, and Kendl (2014); Kendl (2015); Held et al. (2016); Wiesenberger et al. (2017) and phenomena like collisionless reconnection Comisso et al. (2013), zonal flows Hahm et al. (1999); Held et al. (2018) and edge localized modes Kendl, Scott, and Ribeiro (2010); Xu et al. (2013). Gyro-fluid models origin from gyro-kinetic theory Frieman and Chen (1982); Dubin et al. (1983); Lee (1983) and rely on precise fluid closures to incorporate kinetic effects. They excel due to their vastly reduced computational cost in comparison to gyro-kinetic models and algebraic simplicity in comparison to drift-fluid models since the gyro-viscous cancellations emerge automatically.

Gyro-fluid models are particularly characterized by closures, which include finite Larmor radius (FLR) and linear polarization density effects down to the gyro-radius scale Knorr et al. (1988); Brizard (1992); Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993). The latter represent the major hallmarks of gyro-kinetic theory. Additionally, gyro-fluid closures can also encompass kinetic collisionless dissipation channels like Landau damping or FLR phase mixing Hammett and Perkins (1990); Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993); Hunana et al. (2018).

Full-F gyro-fluid models Strintzi and Scott (2004); Strintzi, Scott, and Brizard (2005); Madsen (2013), as opposed to their $\delta\text{F}$ counterpart Knorr et al. (1988); Brizard (1992); Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993); Waltz, Kerbel, and Milovich (1994); Beer and Hammett (1996); Scott (2000); Snyder and Hammett (2001); Scott (2010), avoid the separation of scales and the concomitant Oberbeck-Boussinesq approximation Oberbeck (1879); Boussinesq (1903). The resulting highly non-linear nature of full-F gyro-fluid models complicates the development of fluid closures for kinetic effects. For this reason, polarization in full-F gyro-fluid models is treated within a long perpendicular wavelength approximation Strintzi and Scott (2004); Strintzi, Scott, and Brizard (2005); Madsen (2013), which dates back to the beginning of gyro-kinetic theory and neglects polarization at short perpendicular wavelengths Dubin et al. (1983). This short perpendicular wavelength polarization shortfall of current full-F gyro-fluid models is first of all a fundamental theoretical issue since the $\delta\text{F}$ polarization charge density is not recovered in the Oberbeck-Boussinesq limit but only in the additional long perpendicular wavelength limit. Secondly, the accomplishment of this shortfall is crucial to accurately predict the stability and transport of single- or multi-scale drift wave modes for arbitrary perpendicular wavelengths. In particular, short perpendicular wavelength structures, which emerge e.g. from ion temperature gradient, trapped electron or interchange modes, can significantly affect the turbulent transport Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993); Dominski et al. (2017); Mishchenko et al. (2019); Wiesenberger, Madsen, and Kendl (2014). Moreover, non-Oberbeck-Boussinesq polarization effects strongly alter the dynamics of zonal flows and filaments Held et al. (2018); Kendl (2015); Wiesenberger et al. (2017).

In this contribution, we overcome the short perpendicular wavelength polarization shortfall and present novel fluid closures for polarization effects for arbitrary perpendicular wavelengths. We find that the exact polarization closure, which rests upon a near Maxwellian distribution function, is not suitable for numerical computations. Consequently, we refine the latter closure to a Padé-based approximation of desired accuracy, which accurately mimics arbitrary perpendicular wavelength polarization effects in the Oberbeck-Boussinesq limit and retains the original non-linear structure. We prove that the latter is pivotal for energetic consistency - a missing feature of previous ad-hoc approximations Idomura, Tokuda, and Kishimoto (2003); Bottino et al. (2004); Ku, Chang, and Diamond (2009); Dominski et al. (2017); Mishchenko et al. (2019).

The remainder of the manuscript is organized as follows. We state the full-F gyro-kinetic framework in Sec. II from which we derive the full-F gyro-fluid moment hierarchy and Poisson equation in Sec. III. The gyro-fluid closure of the full-F gyro-fluid moment hierarchy and Poisson equation is discussed in Sec. IV, which contains the principal results of the manuscript. In particular, it includes explicit closures for gyro-average effects in Sec. IV.1, which enter also the novel closures for polarization effects in Sec. IV.2. The presented closures exhibit the correct Oberbeck-Boussinesq and long perpendicular wavelength limit and retain energetic consistency, which is discussed in Sec. IV.2.4 and Sec. IV.3, respectively. Finally, we conclude our findings in Sec. V.

II Gyro-kinetic fundament

Our discussion is based on the nonlinear electrostatic collisionless gyro-kinetic Vlasov-Poisson system Frieman and Chen (1982); Dubin et al. (1983); Lee (1983), which we consistently rederive in this section via field theoretical methods Sugama (2000); Brizard (2000); Brizard and Hahm (2007); Scott (2010). We start right after the transformation of the Lagrangian from particle phase-space coordinates ${\boldsymbol{z}}$ to the gyro-center phase-space corodinates ${{\boldsymbol{Z}}}$, which can be accomplished via Lie-transformation methods. Following the approach developed in Ref. Sugama (2000) we define the gyro-center action

[TABLE]

which consists of the gyro-center particle Lagrangian $L_{p}:=\gamma_{i}\dot{{Z}}^{i}-H$, a vanishing gyro-center field Lagrangian and the gyro-center distribution function ${F}({{\boldsymbol{Z}}},t)$. Here, the [math]-subscript denotes the initial values $({{\boldsymbol{Z}}}_{0},t_{0})$ and we introduced the gyro-center phase-space coordinates ${{\boldsymbol{Z}}}:=\left({{\boldsymbol{X}}},{{\boldsymbol{p}}}\right)\equiv\left({{\boldsymbol{X}}},{\mu},{v}_{\parallel},{\theta}\right)$, which encompass the gyro-center position ${{\boldsymbol{X}}}$, magnetic moment ${\mu}$, parallel velocity ${v}_{\parallel}$ and gyro-angle ${\theta}$. For convenience we omit the species subscript $\alpha$. In the gyro-center particle Lagrangian the fundamental gyro-center Poincaré 1-form and the Hamiltonian are given by

[TABLE]

where we introduced the particle charge $q$ and mass $m$ and the magnitude $B:=|{\boldsymbol{B}}|$ of the magnetic field ${\boldsymbol{B}}$. In the field part the gyro-center Hamiltonian of Eq. (3) appears the central quantity of our work - the gyro-center potential $\Psi:=\Psi_{1}+\Psi_{2}$. It is the superposition of the gyro-average and polarization contribution of the electric potential $\phi$,

[TABLE]

respectively. For the definition of the gyro-average $\langle f\rangle({{\boldsymbol{X}}},{v}_{\parallel},{\mu},t):=\frac{1}{2\pi}\int_{0}^{2\pi}d{\theta}f({\boldsymbol{x}},{v}_{\parallel},t)$ and all therefrom deduced operators a Taylor series representation in configuration space is adopted, since in Fourier space convolution integrals appear, which complicate a simple and clear presentation. Thus, the Taylor series expansion of the gyro-average reads

[TABLE]

where the coefficients belong to the zeroth Bessel function $J_{0}(ix)$ around $x=0$ Frieman and Chen (1982); Lee (1983); Dubin et al. (1983). For the sake of clarity, we define the particle position ${\boldsymbol{x}}={{\boldsymbol{X}}}+{\boldsymbol{\rho}}_{0}({{\boldsymbol{X}}},{\mu},{\theta})$, gyro-radius $\varrho:=\sqrt{2B{\mu}/m}/\Omega$, gyro-frequency $\Omega:=qB/m$, gyro-arm ${\boldsymbol{\varrho}}:=\varrho\hat{{\boldsymbol{a}}}$ and magnetic field unit vector $\hat{{\boldsymbol{b}}}:={\boldsymbol{B}}/B$. The latter forms with $\hat{{\boldsymbol{a}}}$ and $\hat{{\boldsymbol{c}}}$ an orthonormal coordinate system $(\hat{{\boldsymbol{a}}},\hat{{\boldsymbol{b}}},\hat{{\boldsymbol{c}}})$. Further, we introduce the abbreviation for the parallel component of a vector ${\boldsymbol{h}}$ with $h_{\parallel}={\boldsymbol{b}}\cdot{\boldsymbol{h}}$ and its perpendicular projection with ${\boldsymbol{h}}_{\perp}:=-\hat{{\boldsymbol{b}}}\times(\hat{{\boldsymbol{b}}}\times{\boldsymbol{h}})$. This defines the perpendicular Laplacian to $\Delta_{\perp}:={\boldsymbol{\nabla}}\cdot{\boldsymbol{\nabla}}_{\perp}$.

In the gyro-center Poincaré 1-form of Eq. (2) the gyro-center vector potential ${\boldsymbol{A}}^{*}:={\boldsymbol{A}}+\frac{m{v}_{\parallel}}{q}\hat{{\boldsymbol{b}}}$ appears, which determines the associated gyro-center magnetic field ${\boldsymbol{B}}^{*}:={\boldsymbol{B}}+\frac{m{v}_{\parallel}}{q}{\boldsymbol{{\nabla}}}\times\hat{{\boldsymbol{b}}}$. The parallel component of ${\boldsymbol{B}}^{*}$ enters the gyro-center phase-space volume

[TABLE]

where $\omega:=d\gamma$ is the symplectic 2-form and $g$ is the determinant of the metric tensor. In order to simplify the notation we also define the spatial and the conjugate momentum space volume

[TABLE]

In this work we assume ${\boldsymbol{B}}^{*}\approx{\boldsymbol{B}}+\frac{m{v}_{\parallel}}{q}({\boldsymbol{{\nabla}}}\times\hat{{\boldsymbol{b}}})_{\perp}$ so that $B^{*}_{\parallel}\approx B$. This simplifies computations and avoids the coordinate singularity at {v}_{\parallel}=-qB/\big{[}m({\boldsymbol{{\nabla}}}\times\hat{{\boldsymbol{b}}})_{\parallel}\big{]}, where without the approximation $B^{*}_{\parallel}=\text{vol}_{{{\boldsymbol{Z}}}}=0$.

We now specified all terms induced by the gyro-center action, which enables us to compute the explicit expression of the gyro-center Vlasov-Poisson system by means of field theory. The gyro-center equations of motion, which enter the gyro-center Vlasov equation, can be derived by the variation of the gyro-center action of Eq. (1) with respect to ${{\boldsymbol{Z}}}$ Sugama (2000). Constraining this action to be vanishing yields the Euler-Lagrange equation $\frac{\partial L_{p}}{\partial{Z}_{i}}-\frac{d}{dt}\frac{\partial L_{p}}{\partial\dot{{Z}}_{i}}=0$. The latter can be manipulated to Hamilton’s equation $\frac{d{Z}^{i}}{dt}=J^{ij}\frac{\partial H}{\partial{Z}^{j}}$, where the components $J^{ij}$ of the Poisson matrix $J:=\omega^{-1}$ emerge. The evaluation of the components of Hamilton’s equation provides the gyro-center equations of motion

[TABLE]

The gyro-center Vlasov equation expresses the conservation of the gyro-center particle distribution function along the particle trajectories $\frac{d}{dt}{F}=0$. We rearrange this to the conservative form of the gyro-center Vlasov equation

[TABLE]

with the gyro-center equations of motion given by Eqs. (8).

In gyro-kinetic field theory the gyro-center Poisson equation follows from the vanishing variation of the gyro-center action of Eq. (1) with respect to $\phi({\boldsymbol{x}},t)$, which produces

[TABLE]

Using Eqs. (4) the variation along the electric potential of Eq. (10) is explicitly evaluated to

[TABLE]

Here, we defined the adjoint of the gyro-average $\langle f\rangle^{\dagger}({\boldsymbol{x}},{v}_{\parallel},t):=\int\text{vol}_{{{\boldsymbol{X}}}}\langle f({{\boldsymbol{X}}},{\mu},{v}_{\parallel},t)\delta\left({\boldsymbol{x}}-{{\boldsymbol{X}}}-{\boldsymbol{\rho}}_{0}\right)\rangle$, since $\int\text{vol}_{{{\boldsymbol{X}}}}\langle f\rangle h({{\boldsymbol{X}}})=\int\text{vol}_{{\boldsymbol{x}}}\langle h\rangle^{\dagger}f({\boldsymbol{x}})$. Note that analogously to the gyro-average of Eq. (5) we can expand the adjoint gyro-average in a Taylor series

[TABLE]

which reveals that the spatial operators act in this case also on the gyro-radius.

III Gyro-fluid moment hierarchy and Poisson equation

Gyro-fluid models are derived by the gyro-center momentum space integrals over the gyro-center distribution function ${F}({{\boldsymbol{Z}}},t)$ times an arbitrary gyro-center phase-space function $\zeta({{\boldsymbol{Z}}},t)$, which defines the gyro-fluid moment Brizard (1992)

[TABLE]

The latter gyro-fluid moment is exploited for the basic gyro-fluid moment quantities

[TABLE]

the gyro-center density, parallel velocity and the perpendicular and parallel pressure and the perpendicular and parallel contributions of the parallel heat flux, respectively. Here, we defined $\tilde{{v}}_{\parallel}:={v}_{\parallel}-U_{\parallel}$. Additionally we relate the perpendicular and parallel pressure to the corresponding temperature via the ideal gas law $P_{\perp}=NT_{\perp}$ and $P_{\parallel}=NT_{\parallel}$. Further, the gyro-fluid moment over the Vlasov-equation (9) multiplied by $\zeta$ yields the general expression for the time evolution of the gyro-fluid moments Brizard (1992)

[TABLE]

with $\frac{d}{dt}\zeta=\left(\frac{\partial}{\partial t}+\dot{{{\boldsymbol{X}}}}\cdot{\boldsymbol{{\nabla}}}+\dot{{v}}_{\parallel}\frac{\partial}{\partial{v}_{\parallel}}\right)\zeta$. Inserting the gyro-center equations of motions of Eqs. (8) into the gyro-fluid moment evolution Eqs. (14) and restricting the phase space function to $\zeta={\mu}^{k}{v}_{\parallel}^{l}$ yields the gyro-fluid moment hierarchy evolution equation

[TABLE]

Here, we defined the curvature ${\boldsymbol{\kappa}}:=\hat{{\boldsymbol{b}}}\cdot{\boldsymbol{{\nabla}}}\hat{{\boldsymbol{b}}}$ and the ${\boldsymbol{{\nabla}}}B$ drift related terms ${\boldsymbol{\mathcal{K}}}_{{\boldsymbol{\kappa}}}:=\frac{1}{B}\hat{{\boldsymbol{b}}}\times{\boldsymbol{\kappa}}$ and ${\boldsymbol{\mathcal{K}}}_{{\boldsymbol{{\nabla}}}B}:=\frac{1}{B}\hat{{\boldsymbol{b}}}\times{\boldsymbol{{\nabla}}}\ln{(B)}$. The moment hierarchy evolution Eqs. (III) encompass in general an infinite set of evolution equations for all non-negative integer powers of \big{\|}{\mu}^{k}{v}_{\parallel}^{l}\big{\|}. We target the formulation of a six moment gyro-fluid model for the basic gyro-fluid moment quantities of Eq. (13) so that $k\leq 1$ and $l\leq 3$. As a consequence we need to close in Eq. (III) not only the higher gyro-fluid moment quantities

[TABLE]

but also $\Psi$ related terms of the form \big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|}.

The remaining part of the gyro-fluid model is the gyro-fluid Poisson equation. A crucial move in the deduction of gyro-fluid closures for gyro-average and polarization effects is to express the gyro-center Poisson Eq. (10) in terms of the gyro-fluid moment \big{\|}\Psi\big{\|} before the variation along $\phi$ is explictly evaluated. Accordingly, we express the gyro-center Poisson equation in terms of the gyro-fluid moment \big{\|}\Psi\big{\|} and split it into a gyro-center charge density $\sum_{\alpha}qN$ and a polarization charge density $-{\boldsymbol{\nabla}}\cdot{\boldsymbol{P}}$ contribution

[TABLE]

The polarization density ${\boldsymbol{P}}$ can be decomposed into first and second order polarization contributions ${\boldsymbol{P}}={\boldsymbol{P}}_{1}+{\boldsymbol{P}}_{2}$, which are associated with $\Psi_{1}$ and $\Psi_{2}$, respectively. The first order polarization density

[TABLE]

contains FLR effects and is independent of the electric potential $\phi$. The second order polarization density

[TABLE]

can be rewritten into a linear polarization density ${\boldsymbol{P}}_{2}=\epsilon_{0}\mathbb{X}\cdot{\boldsymbol{\nabla}}_{\perp}\phi$ and consequently contributes only for a non-vanishing electric field Krommes (1993). Here, $\epsilon_{0}$ is the vacuum permittivity and $\mathbb{X}$ is the rank-2 electric susceptibility tensor, which contains 2nd order polarization effects for arbitrary perpendicular wavelengths. From now on we associate first and second order polarization contributions with gyro-average (or FLR) and polarization effects, respectively.

The gyro-fluid moment hierarchy evolution Eqs. (III) and the Poisson equation (17) demand a closure for terms of the form \big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|} and \big{\|}\Psi\big{\|}. The accurate evaluation of the closure terms \big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|} and \big{\|}\Psi\big{\|} needs an infinite set of gyro-fluid moment quantities, even if a closed set of evolution equations, as for instance for the basic gyro-fluid moment quantities of Eq. (13) with $k\leq 1$ and $l\leq 3$, is calculated. This peculiarity of gyro-fluid models is reasoned in the gyro-average of Eq. (5), which entails all non-negative integer powers of the gyro-center magnetic moment ${\mu}$ Dorland and Hammett (1993). In particular, the gyro-fluid moment over the gyro-average and polarization part of the gyro-center potential $\Psi$ can be Taylor expanded according to

[TABLE]

Here, the norm of a rank-n tensor $\mathbb{T}$ is defined by $\left|\mathbb{T}\right|:=\sqrt{\mathbb{T}_{i_{1},\dots,i_{n}}\mathbb{T}^{i_{1},\dots,i_{n}}}$ and ${\boldsymbol{{\nabla}}}_{\perp}^{i}$ represents $(i-1)$-times a tensor product, e.g. for $i=2$ we get ${\boldsymbol{{\nabla}}}_{\perp}^{2}={\boldsymbol{{\nabla}}}_{\perp}{\boldsymbol{{\nabla}}}_{\perp}$.

IV Gyro-fluid closure

In full-F theory so far two different strategies are pursued to resolve this infinite hierarchy, which is triggered by the \big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|} and \big{\|}\Psi\big{\|} terms. The first most trivial approach truncates the Taylor series expansion at a specific order and opts usually towards a long perpendicular wavelength approximation, so that \big{\|}\Psi_{1}\big{\|}\approx\left(N+\frac{P_{\perp}{\Delta}_{\perp}}{2m\Omega^{2}}\right)\phi and \big{\|}\Psi_{2}\big{\|}\approx-\frac{qN}{2m\Omega^{2}}\left|{\boldsymbol{{\nabla}}}_{\perp}\phi\right|^{2}. However, such approximated gyro-fluid models are not advantageous over drift-fluid models, since they treat the FLR and polarization terms with the same level of detail. As opposed to this, the second approach retains arbitrary perpendicular wavelength effects by restricting the gyro-center distribution function ${F}$ to a specific form, e.g. a Maxwellian as was originally used by Refs. Knorr et al. (1988); Brizard (1992). In general, this approach truncates the Laguerre-Hermite expansion of the gyro-center distribution function ${F}$ at a particular polynomial order. As a consequence, FLR and polarization effects are kept for full-F to a specific order of the Taylor series expansion and for truncated ${F}$ above this order.

In $\delta\text{F}$ theory also other closure strategies can be utilized, which rest upon the linear gyro-center solution of the Vlasov equation Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993). However, such methods are inapplicable in full-F theory, since full-F closures must avoid both the splitting of the gyro-center distribution function ${F}={F_{0}}(1+\delta F)$ into a stationary Maxwellian ${F_{0}}$ and small fluctuating non-Maxwellian $\delta F$ part and the Oberbeck-Boussinesq approximation.

We follow now the second approach and evaluate the gyro-fluid moments over the various closure terms \big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|} and \big{\|}\Psi\big{\|} to arbitrary perpendicular wavelengths by specifing the gyro-center distribution function ${F}$. First we expand the gyro-center distribution function in a Laguerre-Hermite polynomial in $({\mu},\tilde{{v}}_{\parallel})$ space  Jorge, Ricci, and Loureiro (2017); Mandell, Dorland, and Landreman (2018); Frei, Jorge, and Ricci (2019)

[TABLE]

with Maxwellian

[TABLE]

expansion coefficients $\mathcal{N}_{ij}$, ${c}_{\parallel}:=\tilde{{v}}_{\parallel}\sqrt{m/(2T_{\parallel})}$ and ${c}_{\perp}:={\mu}B/T_{\perp}$. Here, we introduced the Laguerre polynomial $L_{i}(x):=\frac{\exp{\left({x}\right)}}{i!}\frac{d^{i}}{dx^{i}}\left[\exp{\left({-x}\right)}x^{i}\right]$ and the physicist’s Hermite polynomial $H_{j}(x):=(-1)^{j}\exp{\left({x^{2}}\right)}\frac{d^{j}}{dx^{j}}\exp{\left({-x^{2}}\right)}$. Secondly, we approximate the gyro-center distribution function by truncating the Laguerre-Hermite expansion at order $\left(I,J\right)$ so that

[TABLE]

This is known as closure by truncation, since all expansion coefficients above $\left(I,J\right)$ are set to zero. However, this approach does not rule out polynomial or asymptotic closures for smaller $\left(I,J\right)$. For the sake of clarity we define also the truncated gyro-fluid moment

[TABLE]

which rests upon the truncated gyro-center distribution function of Eq. (23).

For the remainder of the manuscript we truncate the gyro-center distribution function of Eq. (23) at $\left(I,J\right)=\left(1,3\right)$ and fix the only non-vanishing expansion coefficients to $\mathcal{N}_{00}:=1$, $\mathcal{N}_{11}:=-\frac{Q_{\perp}}{P_{\perp}}\sqrt{\frac{m}{2T_{\parallel}}}$ and $\mathcal{N}_{03}:=\frac{Q_{\parallel}}{12P_{\parallel}}\sqrt{\frac{m}{2T_{\parallel}}}$. This results in a near Maxwellian distribution function ${F}_{13}$, which could be written as in Ref. Madsen (2013)

[TABLE]

with $\xi:=\sqrt{\frac{2m}{T_{\parallel}}}{c}_{\parallel}\left[\frac{Q_{\perp}}{P_{\perp}}\left({c}_{\perp}-1\right)+\frac{Q_{\parallel}}{6P_{\parallel}}\left(2{c}_{\parallel}^{2}-3\right)\right]$. As a consequence of the near Maxwellian assumption of Eq. (25) the gyro-fluid hierarchy is closed by truncation, because the higher gyro-fluid moment quantities of Eq. (16) are expressed through the basic gyro-fluid moments

[TABLE]

In certain cases a different closure, such as a polynomial or asymptotic (e.g. collisional) closure, is intended instead of the closure by truncation of Eq. (26). Then one must either resort to a four moment model with $(N,U_{\parallel},P_{\perp},P_{\parallel})$, where the chosen closure is determined by $(Q_{\perp},Q_{\parallel})$, or adopt the higher truncation $\left(I,J\right)=\left(2,5\right)$, so that $(R_{\perp},R_{\parallel,\perp},R_{\parallel},S_{\parallel,\perp},S_{\perp},S_{\parallel})$ retake the role of unspecified closure variables. However, the gyro-average and polarization closures are increasingly complex to implement numerically for higher truncations, such as $\left(I,J\right)=\left(2,5\right)$, than those at $\left(I,J\right)=\left(1,3\right)$. Thus, from the practical point of view it is reasonable to close the higher gyro-fluid moment quantities in the gyro-average and polarization terms differently than in the remaining gyro-fluid model, even though full consistency within the gyro-fluid model is lost. The derivation of consistent full-F gyro-fluid closures for e.g. $(Q_{\perp},Q_{\parallel})$ at arbitrary collisionalities is an ongoing effort, which will be reported in a future work.

The near Maxwellian assumption of Eq. (25) allows us to evaluate the various gyro-fluid moment closure terms \big{\|}{\mu}^{k}{v}_{\parallel}^{l}{\boldsymbol{{\nabla}}}\Psi\big{\|} and \big{\|}\Psi\big{\|}, which generate a hierarchy of gyro-fluid potentials. These gyro-fluid moment closure terms are evaluated for a six moment gyro-fluid model to

[TABLE]

with $\eta:=\ln{(B/T_{\perp})}$. The gyro-fluid moment closure terms of Eqs (27) produce the basic $\psi:=\psi_{1}+\psi_{2}$ and higher gyro-fluid potentials $\chi:=\chi_{1}+\chi_{2}$ and $\iota:=\iota_{1}+\iota_{2}$. The gyro-fluid potentials $(\psi,\chi,\iota)$ arise from the gyro-fluid moment over the gyro-center potential $\Psi$. The respective gyro-average and polarization parts of the gyro-fluid potentials are defined by

[TABLE]

with $i\in\left\{1,2\right\}$. This hierarchy retains gyro-average ($\Psi_{1}$) as well as polarization ($\Psi_{2}$) effects through $\left(\psi_{1},\chi_{1},\iota_{1}\right)$ and $\left(\psi_{2},\chi_{2},\iota_{2}\right)$, respectively.

The gyro-average of Eq. (5) is inherent to both the closure for gyro-average and polarization effects (cf. (4)). As a consequence, the higher gyro-fluid potentials $\left(\chi,\iota\right)$ of not only the gyro-average $\left(\chi_{1},\iota_{1}\right)$ but also of the polarization $\left(\chi_{2},\iota_{2}\right)$ can be related to the basic gyro-fluid potential $\psi$ by the simple and exact recursive identities Dorland and Hammett (1993)

[TABLE]

Here, we defined the thermal gyro-radius $\rho:=\sqrt{T_{\perp}m}/(qB)$. These formulas provide consistent closures for the gyro-average and polarization contributions of Eqs. (28b) and (28c), even if an approximated closure of the lowest moment $\psi$ is utilized. This result extends the recursive identities, originally derived for the gyro-average contributions Dorland and Hammett (1993), to the polarization contributions.

IV.1 Gyro-average closures

The gyro-average contributions in the gyro-fluid potentials $\psi_{1}$ and $\chi_{1}$, $\iota_{1}$ give rise to the basic and higher FLR operators, $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}$ and $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{2}$, $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{3}$, respectively. These FLR operators are defined via the gyro-fluid potentials according to

[TABLE]

The near Maxwellian distribution function of Eq. (25) allows us to evaluate these operators for arbitrary perpendicular wavelengths. The consequent basic FLR operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}$ is given in configuration space by a Taylor series

[TABLE]

with the familiar Taylor series coefficients of the exponential $\exp{(x^{2}/2)}$ at $x=0$ Knorr et al. (1988); Brizard (1992). Note that Eq. (31) renders the Taylor expansion for full-F of Eq. (20a) up to $\mathcal{O}(b^{2})$ with $b:=\rho k_{\perp}$. As a result it features full-F FLR effects up to $\mathcal{O}(b^{2})$ and near Maxwellian FLR effects for all higher even orders of $b$.

As soon as the basic FLR operator is determined either by the near Maxwellian assumption (Eq. (31)) or by an additional approximation (discussed later), the higher FLR operators and the first order polarization (FLR) contribution in the Poisson equation follow immediately. For the higher FLR operators $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{2}$ and $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{3}$ the recursive identities of Eqs. (29) are utilized in combination with Eq. (31). This yields $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{2}=\sum_{i=0}^{\infty}\frac{(\rho^{2}\Delta_{\perp})^{i}}{2^{i}(i-1)!}$ and $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{3}=\sum_{i=0}^{\infty}\frac{(1+i)(\rho^{2}\Delta_{\perp})^{i}}{2^{i}(i-1)!}$. The first order polarization charge density is produced by the adjoint $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}^{\dagger}$ of the basic FLR operator Strintzi and Scott (2004); Madsen (2013)

[TABLE]

and emerges from the variationial of Eq. (18).

While the latter FLR operators accurately capture gyro-averaging effects for the near Maxwellian distribution function, they fail to match the linear gyro-kinetic solution and consequently the ion temperature gradient marginal stability relation for a finite set of gyro-moments in a slab and constant magnetic field Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993); Mandell, Dorland, and Landreman (2018). Therefore, approximations to Eq. (31) are needed, which are well behaved and capture both gyro-averaging effects and the linear response and consequently the ITG instability properly.

The $\sqrt{\Gamma_{0}}$ approximation overcomes this drawback Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993) and replaces the basic FLR operator by

[TABLE]

Here, the linear polarization operator is defined by

[TABLE]

where the Gamma function $\Gamma$ is not to be mistaken with the gyro-average or polarization operator. It originates from the Oberbeck-Boussinesq limit of the Poisson equation \prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}(\phi)=\big{\|}\langle\langle\phi\rangle\rangle^{\dagger}\big{\|}_{13}/N, which we relax after integration for the thermal gyro-radius $\rho$. The Taylor series coefficients of the linear polarization operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$ correspond to the well-known function $I_{0}(-x^{2})\exp{(x^{2})}$ at $x=0$ Dubin et al. (1983); Lee (1983), where $I_{n}$ is the modified Bessel function.

From the numerical point of view neither the near Maxwellian basic FLR operator (Eq. (31)) nor its $\sqrt{\Gamma_{0}}$ approximation (Eq. (33)) are practical in full-F gyro-fluid models. This is because in configuration space accuracy to arbitrary perpendicular wavelengths is lost due to a truncation of the Taylor series and in Fourier space computationally demanding convolution integrals emerge. Padé-approximations of order $(P,Q)$ in configuration space offer a way out of this difficulty. The $(P,Q)$ Padé-approximation is a uniquely determined rational approximation to a function or in general an operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Upsilon}$, which we define as

[TABLE]

Here, the series coefficients $\alpha_{i}$ and $\beta_{i}$ are determined from the condition $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Upsilon}^{(P+Q)}(0)=\left(\prescript{\scriptscriptstyle P}{\scriptscriptstyle Q}{\Upsilon}\right)^{(P+Q)}(0)$. In this work these coefficients are deduced by the PadeApproximant routine of Wolfram Mathematica Wolfram Research, Inc., . Proper approximations to the basic FLR operator $\Gamma_{1}$ are based on a suitable $(P,Q)$ Padé-approximation, abbreviated by $\prescript{\scriptscriptstyle P}{\scriptscriptstyle Q}{\Upsilon}$, to the chosen operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Upsilon}$. For the gyro-average part of the basic gyro-fluid potential we choose for $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Upsilon}$ the polarization operator $\Gamma_{0}$ (or its square root) whereas for the polarization part of the basic gyro-fluid potential, as we discuss later, we pick the basic FLR operator $\Gamma_{1}$ (or its square). Two $\mathcal{O}(b^{2})$ accurate and at $b\rightarrow\infty$ vanishing Padé-approximations for the basic FLR operator $\prescript{\scriptscriptstyle P}{\scriptscriptstyle Q}{\check{\Gamma}}_{1}$ emerge at order $(1,2)$ Hammett, Dorland, and Perkins (1992)

[TABLE]

and $(1,4)$ Dorland and Hammett (1993); Strintzi and Scott (2004); Madsen (2013)

[TABLE]

Note that for the latter $(1,2)$ and $(1,4)$ Padé-approximation the $\sqrt{\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}}$ relationship of Eq. (33) is exactly fulfilled, so that $\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\left(\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}^{2}\right)}=\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\Gamma}_{0}$ and $\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\Gamma}_{1}=\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\sqrt{\Gamma}_{0}}$, respectively. The higher Padé approximated FLR operators follow from the recursive identities of Eq. (29) and are summarized together with the basic Padé approximated FLR operator and its adjoint in Table 1 Held et al. (2016).

The Padé approximated FLR operators of Table 1 are also utilized for the polarization closures, which are discussed in the next Sec. IV.2.

Finally, we depict the near Maxwellian FLR operators and their proposed approximations in Fig. 1. Here, the Padé-approximations qualitatively agree with the square root approximated FLR operators.

IV.2 Polarization closures

Analogous to the gyro-average closures of Sec. IV.1 the polarization closures rely on the near Maxwellian assumption of Eq. (25). The elementary quantity for the polarization closures is the polarization part of the basic gyro-fluid potential $\psi_{2}$ (Eq. (28a). The derivation of the near Maxwellian polarization closure follows the same methodology than that of the near Maxwellian gyro-average closure, but is more involved. We calculate $\psi_{2}$ by (i) transforming only the gyro-average terms in Eq. (4) to Fourier space, (ii) evaluating the gyro-fluid moment with the help of the integral identities \big{\|}J_{0}(\varrho K_{\perp})\big{\|}_{13}=e^{-\rho^{2}K_{\perp}^{2}/2}=e^{-\rho^{2}(k_{\perp}^{2}+k_{\perp}^{\prime 2})/2}e^{-\rho^{2}{\boldsymbol{k}}_{\perp}\cdot{\boldsymbol{k}}_{\perp}^{\prime}} and \big{\|}J_{0}(\varrho k_{\perp})J_{0}(\varrho k^{\prime}_{\perp})\big{\|}_{13}=e^{-\rho^{2}(k_{\perp}^{2}+k_{\perp}^{\prime 2})/2}I_{0}(\rho^{2}k_{\perp}k_{\perp}^{\prime}) Brizard (1992) where ${\boldsymbol{K}}_{\perp}={\boldsymbol{k}}_{\perp}+{\boldsymbol{k}}_{\perp}^{\prime}$, (iii) Taylor expanding the function $e^{-\rho^{2}{\boldsymbol{k}}_{\perp}\cdot{\boldsymbol{k}}_{\perp}^{\prime}}-I_{0}(\rho^{2}k_{\perp}k_{\perp}^{\prime})=2\sum_{i=1}^{\infty}\rho^{2i}\left(a_{i}\mathrm{i}^{2i}\langle{\boldsymbol{k}}_{\perp}^{i},{\boldsymbol{k}}_{\perp}^{\prime i}\rangle-b_{i}k_{\perp}^{i}k_{\perp}^{\prime i}\right)$, where $\langle\mathbb{A},\mathbb{B}\rangle:=\mathbb{A}_{i_{1},\dots,i_{n}}\mathbb{B}^{i_{1},\dots,i_{n}}$ fulfills the identity $\langle\mathbb{A},\mathbb{A}\rangle=\left|\mathbb{A}\right|^{2}$ and finally (iv) performing the inverse Fourier transformation. This yields the polarization part of the basic near Maxwellian gyro-fluid potential $\psi_{2}$, which includes only the square of linear differential operators:

[TABLE]

Here, we introduced the Taylor series coefficients

[TABLE]

of the function $\left[\exp(x^{2})-1\right]/2$ and $\left[I_{0}(x^{2})-1\right]/2$ around $x=0$, respectively. The first few coefficients are explicitly given by $a_{i}=\left\{1/2,1/4,1/12,1/48,\dots\right\}$ and $b_{i}=\left\{0,1/8,0,1/128,\dots\right\}$. Analogously, to the near Maxwellian gyro-average closure of Eq. (31) the near Maxwellian polarization closure of Eq. (38) reproduces the full-F expression of Eq. (20) to $\mathcal{O}(b^{4})$. Consequently, it is hallmarked by full-F polarization effects up to $\mathcal{O}(b^{4})$ and near Maxwellian polarization effects for all higher even orders of $b$.

From the variational of Eq. (19) with the near Maxwellian polarization closure of Eq. (38) we obtain the near Maxwellian polarization density

[TABLE]

which agrees with the general Laguerre-Hermite expanded expression of Ref. Frei, Jorge, and Ricci (2019) in the near Maxwellian and $B^{*}_{\parallel}\approx B$ limit. Here, we introduced the notation $({\boldsymbol{\nabla}}\cdot)^{i}$ for $i$-times a dot product. For example, for $i=2$ we obtain $\left({\boldsymbol{\nabla}}\cdot\right)^{2}={\boldsymbol{\nabla}}\cdot{\boldsymbol{\nabla}}\cdot$.

IV.2.1 Arbitrary order approximation

The near Maxwellian polarization closure of Eq. (38) together with the higher polarization closures, based on Eqs. (29), and the polarization density of Eq. (IV.2) are unpractical for numerical computations due to an infinite set of differential operators. In Fourier space this statement converts to an infinite number of convolutions, which pertains also for the Hermite-Laguerre expanded formulation Frei, Jorge, and Ricci (2019). Thus, the remaining task is to find truncations of the infinite series expression of Eq. (38), which (i) retain the basic quadratic structure of Eq. (38) for energetic consistency and (ii) feature a polarization density with the proper Oberbeck-Boussinesq and long perpendicular wavelength limit. These requirements are fulfilled if we truncate at $K\in\left\{2\mathbb{N}+1\right\}\cup\left\{2\right\}$ and utilize the Padé-approximations

[TABLE]

instead of the basic FLR operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}$. This yields the $\mathcal{O}\left(b^{2K}\right)$ accurate rational approximation to the polarization part of the basic gyro-fluid potential

[TABLE]

which mimics arbitrary perpendicular wavelength polarization effects through the $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}\approx\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2K}{\check{\Gamma}}_{1}$ approximation. Analogously to Eq. (IV.2) the corresponding $\mathcal{O}\left(b^{2K}\right)$ approximated polarization density to Eq. (46) is derived to

[TABLE]

Note that the $\mathcal{O}\left(b^{2K}\right)$ polarization closures retain only $\mathcal{O}\left(b^{4}\right)$ full-F polarization effects due to the near Maxwellian assumption. Thus, we derive in the following explicit and numerically feasible expressions for the truncated polarization closures at $K\in\left\{1,2\right\}$.

IV.2.2 Second order approximation

In the $K=1$ case Eq. (46) reduces to an $\mathcal{O}(b^{2})$ accurate approximation

[TABLE]

which features arbitrary wavelength polarization effects through the $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}\approx\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\check{\Gamma}}_{1}$ approximation. The remaining polarization parts of the higher gyro-fluid potentials $\chi_{2},\iota_{2}$ are consistently derived from the recursive closure formulas of Eqs. (29) together with the polarization part of the basic gyro-fluid potential $\psi_{2}$ of Eq. (48):

[TABLE]

The associated $\mathcal{O}(b^{2})$ accurate polarization density

[TABLE]

follows from Eq. (IV.2.1). The Padé-approximated FLR operators appearing in Eqs. (48)-(50) are to be found in the $q=2$ row of Table 1.

We stress that the herein proposed second order polarization charge density $-{\boldsymbol{\nabla}}\cdot{\boldsymbol{P}}_{2}\approx\sum_{\alpha}\frac{q^{2}}{m}{\boldsymbol{\nabla}}\cdot\left({\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\check{\Gamma}}_{1}^{\dagger}}\frac{N}{\Omega^{2}}\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\check{\Gamma}}_{1}{\boldsymbol{\nabla}}_{\perp}\phi\right)$ does not agree with the widely used ad-hoc second order Padé-approximation in gyro-kinetic models $-{\boldsymbol{\nabla}}\cdot{\boldsymbol{P}}_{2}\approx\sum_{\alpha}\frac{q^{2}}{m}\left[1-{\boldsymbol{\nabla}}\cdot(\rho^{2}{\boldsymbol{\nabla}}_{\perp})\right]^{-1}{\boldsymbol{\nabla}}\cdot\left(\frac{N}{\Omega^{2}}{\boldsymbol{\nabla}}_{\perp}\phi\right)$ Idomura, Tokuda, and Kishimoto (2003); Bottino et al. (2004); Ku, Chang, and Diamond (2009); Dominski et al. (2017); Mishchenko et al. (2019), which is not producing a quadratic kinetic ${\boldsymbol{E}}\times{\boldsymbol{B}}$ energy. Non-linear gyro-kinetic simulations with this ad-hoc Padé-approximation and a full-F arbitrary wavelength treatment of the polarization charge density show only slight differences for typical turbulent transport observables, such as the heat diffusivity Dominski et al. (2017).

IV.2.3 Fourth order approximation

In the $K=2$ case the closure derived from Eq. (46) reads

[TABLE]

has $\mathcal{O}\left(b^{4}\right)$ accuracy and keeps arbitrary wavelength effects due to the $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}\approx\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 4}{\check{\Gamma}}_{1}$ approximation. Analogous to the second order approximation we derive from the fourth order approximation of Eq. (IV.2.3) the polarization parts of the higher gyro-fluid potentials

[TABLE]

and the polarization density

[TABLE]

Again we refer the reader for the various Padé-approximated FLR operators in Eqs. (IV.2.3)-(IV.2.3) to the $q=4$ row of Table 1.

The fourth order approximation ($K=2$) is more complex to implement numerically than the second order approximation ($K=1$) due to higher order spatial derivatives, but it does not require to compute the square root of an operator (cf. Table 1).

In summary, the proposed approximated polarization closures of Eqs. (46)-(IV.2.1) or more explicitly of Eqs. (48)-(50) and Eqs. (IV.2.3)-(IV.2.3) are accurate up to $\mathcal{O}(b^{2K})$, $\mathcal{O}(b^{2})$ and $\mathcal{O}(b^{4})$, respectively. Further they imply accurate polarization effects at arbitrary perpendicular wavelengths in the Oberbeck-Boussinesq limit, the correct long perpendicular wavelength limit and produce an appropriate energy conservation law, which is demonstrated in the following.

IV.2.4 Oberbeck-Boussinesq and long perpendicular wavelength limit

In the Oberbeck-Boussinesq limit the spatial dependence of gyro-fluid moment variables ($N,P_{\perp},\dots$) and the magnetic field magnitude $B$ is neglected in the polarization closure terms. Further, only the stationary contributions of the gyro-fluid moment variables ($N_{0},P_{\perp 0},\dots$), which result from the various gyro-fluid moments with the stationary gyro-center Maxwellian ${F_{0}}$, are taken into account. As a consequence, the polarization and FLR operators are self-adjoint (e.g. $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}=\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{1}^{\dagger}$), commute with spatial derivatives and are to be understood to contain only contributions from the stationary thermal gyro-radius $\rho_{0}:=\sqrt{T_{\perp 0}m}/(qB)$. Thus, the basic near Maxwellian second order polarization charge density resulting from Eq. (IV.2) reduces in the Oberbeck-Boussinesq limit to

[TABLE]

where we used $({\boldsymbol{\nabla}}\cdot)^{i}{\boldsymbol{\nabla}}_{\perp}^{i}=\Delta_{\perp}^{i}$. Note now that the polarization charge density of Eq. (54) is similar to the one directly obtained from the gyro-center Poisson Eq. (10) with ${F}\approx{F_{0}}$.

The Oberbeck-Boussinesq limit of the proposed $\mathcal{O}\left(b^{2K}\right)$ approximation of Eq. (IV.2.1) produces a polarization charge density equal to Eq. (54) except that the near Maxwellian polarization operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$ is replaced by the truncated polarization operator ${\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0,K}}:=1+2({\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2K}{\check{\Gamma}}_{1}})^{2}\sum_{i=1}^{K}\left(b_{i}-\mathrm{i}^{2i}a_{i}\right)\left(\rho_{0}^{2}\Delta_{\perp}\right)^{i}$, which converges to $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$ for $K\rightarrow\infty$. Strikingly, for the second and fourth order approximations (Eq. (50) and Eq. (IV.2.3)) the truncated polarization operator resembles the $(1,2)$ or $(1,4)$ Padé-approximation of the polarization operator, so that $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0,1}=\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\Gamma}_{0}$ or $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0,2}=\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 4}{\Gamma}_{0}$. The latter Padé-approximations retain high accuracy to the near Maxwellian polarization operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$, which is depicted in Fig. 2. The relative error of the $(1,2)$ and $(1,4)$ Padé-approximation to $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$ are of comparable magnitude and are roughly 7% and 11%, respectively. Additionally, the $\mathcal{O}\left(b^{10}\right)$ and $\mathcal{O}\left(b^{18}\right)$ accurate approximations are shown.

In the long perpendicular wavelength limit the presented near Maxwellian and approximated polarization closures reduce to the long perpendicular wavelength polarization expressions $\psi_{2}\approx-\frac{q}{2m\Omega^{2}}\left|{\boldsymbol{{\nabla}}}_{\perp}\phi\right|^{2}$, ${\boldsymbol{P}}_{2}\approx-\sum_{\alpha}\frac{q^{2}N}{m\Omega^{2}}{\boldsymbol{\nabla}}_{\perp}\phi$ and $\chi_{2}\approx\iota_{2}\approx 0$ Madsen (2013).

IV.3 Gyro-fluid energy conservation

The proposed arbitrary perpendicular wavelength polarization closure extension manifests itself only in the kinetic ${\boldsymbol{E}}\times{\boldsymbol{B}}$ energy part of the total gyro-fluid energy conservation law. The conservation of the total energy $\frac{d}{dt}E=0$ is obtained from the volume integral and species sum over the gyro-fluid moment evolution Eq. (14), where we replace the phase-space function $\zeta$ by the Hamiltonian $H$ of Eq. (3) and the gyro-center distribution function ${F}$ by the near Maxwellian ${F}_{13}$ of Eq. (25). This yields the total energy

[TABLE]

which is the superposition of the internal perpendicular and parallel energy, the parallel kinetic energy and the positive definite kinetic ${\boldsymbol{E}}\times{\boldsymbol{B}}$ energy

[TABLE]

In the course of the derivation of energy conservation law we made use of the time evolution equation of the kinetic ${\boldsymbol{E}}\times{\boldsymbol{B}}$ energy, which is explicitly computed to $\frac{dE_{k}}{dt}=\sum_{\alpha}q\int\text{vol}_{{{\boldsymbol{X}}}}\left[\left(\psi-\chi\right)\frac{\partial N}{\partial t}+\frac{\chi}{T_{\perp}}\frac{\partial P_{\perp}}{\partial t}\right]$, the gyro-fluid Poisson equation of Eq. (17) and assumed vanishing surface integral contributions. We stress that two properties must hold for exact conservation of the total energy of Eq. (55). First, $\psi_{2}$ must be made of squares of linear differential operators, as given by both the near Maxwellian or truncated expression. Second, both $\chi_{1}$ and $\chi_{2}$ must be derived by the recursive identity of Eq. (29). This enables us to retain energy conservation even if simpler approximations for the remaining higher FLR operators are utilized, for instance the long perpendicular wavelength fit $\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\check{\Gamma}}_{3}\approx 2\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 2}{\check{\Gamma}}_{2}$ or $\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 4}{\check{\Gamma}}_{3}\approx 2\prescript{\scriptscriptstyle 1}{\scriptscriptstyle 4}{\check{\Gamma}}_{2}$ Dorland and Hammett (1993). Finally, we note that total energy of Eq. (55) agrees in the long perpendicular wavelength limit with Ref. Madsen (2013).

V Conclusions

In this work novel full-F gyro-fluid closures are derived for polarization effects for arbitrary perpendicular wavelengths, which overcome limitations of former truncated or ad-hoc polarization closures. Based on a near Maxwellian assumption explicit expressions for the polarization part of the basic gyro-fluid potential and the polarization density are derived in Eqs. (38)-(IV.2). However, these expressions contain an infinite set of spatial operators and require suitable approximations for numerical computations.

Thus, a general approximation is deduced from a series truncation of the polarization part of the near Maxwellian basic gyro-fluid potential of Eq. (38), which retains the original quadratic structure of the expression for energetic consistency and replaces the inherent FLR operators by appropriate Padé-approximations. The resulting Padé-based closure for the polarization part of the basic gyro-fluid potential of Eq. (46) is $\mathcal{O}\left(b^{2K}\right)$ accurate. Notably, the associated polarization density of Eq. (IV.2.1) comprises polarization effects to arbitrary perpendicular wavelengths since its Oberbeck-Boussinesq limit (Eq. (54)) yields an $\mathcal{O}\left(b^{2K}\right)$ accurate rational approximation to the polarization operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$. The truncated polarization closures are specified to the tractable limit of second and fourth order accuracy in Eqs. (48)-(50) respective Eqs. (IV.2.3)-(IV.2.3), where the polarization parts of the higher gyro-fluid potentials of Eqs. (49) and Eqs. (52) are consistently closed by the recursive identities of Eq. (29). In this limit the rational approximation of the polarization operator $\prescript{\scriptscriptstyle}{\scriptscriptstyle}{\Gamma}_{0}$ is optimal so that it reduces to its $(1,2)$ and $(1,4)$ Padé-approximation. The truncated near Maxwellian polarization closures ensure energy conservation with a positive definite kinetic ${\boldsymbol{E}}\times{\boldsymbol{B}}$ energy as given by Eq. (56).

We emphasize that the proposed approximations for gyro-averaging (Table. 1) and polarization (Eqs. (48)-(50) or Eqs. (IV.2.3)-(IV.2.3)) are the full-F generalization of the widely used $\delta\text{F}$ Padé-model Hammett, Dorland, and Perkins (1992); Dorland and Hammett (1993); Beer and Hammett (1996); Snyder and Hammett (2001); Scott (2010). Finally, the proposed second and fourth order Padé-based approximations for arbitrary perpendicular wavelength polarization effects can be also utilized in gyro-kinetic models, when a (near) Maxwellian gyro-center distribution function is assumed for the polarization contributions in the gyro-center Vlasov-Poisson system. In particular, we recommend that the ad-hoc second order Padé-approximation for arbitrary perpendicular wavelength polarization effects in gyro-kinetic codes Idomura, Tokuda, and Kishimoto (2003); Bottino et al. (2004); Ku, Chang, and Diamond (2009); Dominski et al. (2017); Mishchenko et al. (2019) is replaced with the herein proposed second order Padé-approximation of Eqs. (48)-(50) in order to restore energetic consistency.

VI Acknowledgements

The authors acknowledge helpful discussion with B. J. Frei.