Gyro-fluid equations are velocity space moments of the gyrokinetic equations. Special gyro-Landau-fluid closures have been developed that include the damping due to kinetic resonances by fitting to the collisionless local plasma response functions. This damping allows for accurate linear eigenmodes to be computed with a relatively low number of velocity space moments compared to the number of velocity quadrature points in gyrokinetic codes. However, none of the published gyro-Landau-fluid closure schemes considers the Onsager symmetries of the resulting quasi-linear fluxes as a constraint. Onsager symmetry guarantees that the matrix of diffusivities is positive definite, an important property for the numerical stability of a transport solver. A two parameter real closure for improving the accuracy of low resolution gyro-fluid equations, which preserves the Onsager symmetry and allows higher velocity space moments, is presented in this paper. The new linear gyro-fluid system (GFS) is used to extend the TGLF quasi-linear transport model so that it can compute the energy and momentum fluxes due to parallel magnetic fluctuations, completing the transport matrix. The GFS equations do not use a bounce average approximation. The GFS equations are fully electromagnetic with general flux surface magnetic geometry, pitch angle scattering for electron collisions, and subsonic equilibrium toroidal rotation. Using GFS eigenmodes in the quasi-linear TGLF model will be shown to yield a more accurate match to fluxes computed by CGYRO turbulence simulations. Prospects for future applications of a quasi-linear theory to new plasma transport regimes and magnetic confinement devices in addition to tokamaks are opened by the flexibility of the GFS eigensolver.

## I. FLUID MODELS OF THE GYROKINETIC EQUATION

The gyrokinetic equations^{1} are a reduction of the full 6D Vlasov–Fokker–Plank equations valid for strongly magnetized plasmas. The fast gyro-orbit motion of the particles is averaged, and an expansion to first order in the ratio of the ion gyro-radius to the equilibrium pressure gradient scale length is made. The resulting gyrokinetic system is valid for wave frequencies small compared to the gyro-frequency. There is also an approximation that the wavenumber parallel to the magnetic field is smaller than the wavenumber perpendicular to the magnetic field. The gyrokinetic equation is usually written with the energy and magnetic moment of the plasma particles as the velocity space variables since these are conserved by the gyro-averaged guiding center motion.

When computers were not powerful enough to simulate turbulence with the non-linear continuum gyrokinetic equations, reduced systems of velocity space fluid moment equations were developed. It was shown by Hammett and Perkins^{2} that the kinetic effect of Landau damping could be very well captured by a low order system of fluid moment equations. The idea was to fit the analytic density kinetic response integral by choosing a closure of the parallel velocity moment fluid system. The highest truncated moment is replaced with a sum of the next two lower moments with fitting coefficients. These closure coefficients could be evaluated analytically for the pole of the parallel advection term $ v | | k | |$ in a simplified kinetic model without magnetic drift terms. This Landau-fluid model was shown to be an accurate fit to the analytic density kinetic response, even for only four parallel velocity moments.

The inclusion of the magnetic curvature drifts introduces additional resonances in the velocity space. A closure for the magnetic drift resonances was found by the same method of taking the highest truncated velocity moment to be a linear sum of lower moment terms.^{3,4} The closure coefficients could not be fit analytically to the kinetic density response, like the Landau pole, and had to be carefully fit with a minimization method.^{5} The fluid moment model with Landau and magnetic drift resonances is called the gyro-Landau-fluid (GLF) model. The lowest number of velocity space moments that gave a reasonable accuracy (GLF23) was employed for state of the art numerical turbulence simulations.^{3} The GLF23 model was later used for a pioneering quasi-linear reduced model of turbulent transport.^{6} These GLF models neglected trapped ions, giving a Landau damping rate to all of the ions without the bounce averaging effect. Bounce average electron models were included,^{3,7} but this introduced the problem of having to change the equations somewhere between ion and electrons scales in order to model the electron temperature gradient instability due to circulating electrons.

The trapped gyro-Landau fluid (TGLF) equations^{8} are designed to include both trapped and passing particles for all species in a single system bridging the gap between ion and electron wavenumber scales. The TGLF equations generalize the GLF closures of Beer^{4} to be dependent on the trapped-passing boundary in velocity space and introduced velocity moment equations covering only the trapped or passing regions. This continuity of scales comes at the expense of complexity of the 20 magnetic drift closure coeffcients, which had to be fit for 21 different trapped fractions by a painstaking minimization process. The velocity moments of the gyro-averaging operator (Bessel functions) in the gyrokinetic equations also become dependent on the trapped-passing boundary. The more accurate linear system TGLF does result in a more accurate quasi-linear representation of gyrokinetic turbulent fluxes. The prediction of the temperature, density, and rotation profiles in the plasmas were also found to be more accurate with TGLF^{9} validating the gyro-kinetic turbulent transport theory.

Recently, there has been a renewed interest in higher velocity space resolution gyro-fluid models.^{10,11} These papers have in common with the present paper the use of parallel velocity (Hermite) and perpendicular energy (Laguerre) polynomials extending the work of Beer.^{4} The aim of this paper is to have a reduced model of sufficient speed and accuracy for quasi-linear flux evaluation rather than for non-linear time-dependent gyrokinetic turbulence modeling. The eigenvalue matrix solution of the linear gyro-fluid system does not require adding dissipation that is typically required for stability of a time dependent solution. It is found that a purely real closure scheme can be used without explicit added damping. The pseudo-spectral approach is not employed in this paper giving more accuracy of the analytically evaluated Bessel function integrals for a small number of perpendicular energy moments. The careful preservation of Onsager symmetry by the numerical solution scheme is also a new contribution.

Extensive use of TGLF to predict plasma profiles in tokamaks has uncovered three problems that will be resolved with the new flexible gyro-fluid system (GFS) eigensolver presented in this paper.

First, the bounce averaging model is not absolute. The ability of trapped particles to average away Landau damping depends on the mode frequency relative to the trapped particle bounce frequency. The original TGLF equations have a model for the loss of bounce averaging when the parallel wavenumber is increased. Recently, this model was amended to include a fit correction for high diamagnetic drift frequencies.^{12} The largest errors between the gyrokinetic linear growthrates of CGYRO^{13} and TGLF occur because of the loss of bounce averaging effect. The bounce averaging approximation will be eliminated in the new GFS equations greatly simplifying the closure scheme and making it possible to vary the velocity space resolution.

Second, the parallel magnetic field fluctuations have been shown to be important for stability at high beta.^{14} Although TGLF is fully electromagnetic, the linear eigenvalues do not accurately capture the parallel magnetic flutter impact when it is important. There are also not enough perpendicular energy moments in TGLF to compute the energy and momentum flux due to parallel magnetic fluctuations. The new flexible GFS in this paper is able to include more perpendicular energy moments and compute these fluxes.

Third, it would be beneficial to the numerical convergence of the transport code solution with TGLF if the effective diffusion and convection coefficients could be output. Newton-type transport solvers take the derivative of the turbulent fluxes with respect to the plasma gradients to find the Jacobian that determines the direction to advance the gradients to match the target fluxes. If the Jacobian is not positive definite the Newton direction will fail to find the solution. Negative Jacobian eigenvalues are not uncommon for TGLF fluxes. It has been shown by Sugama and Horton (SH) that the turbulent fluxes satisfy Onsager symmetries^{15} like a neoclassical transport theory. The Onsager symmetries ensure that the effective diffusion coefficients in the flux-gradient relations form a positive definite, symmetric matrix, giving positive entropy production. This property was shown for the quasi-linear continuum gyrokinetic theory. It will be shown below that the GLF closures break these symmetries for the moment equations. The GFS equations of this paper will carefully preserve the Onsager symmetries so that the effective diffusion and convection coefficients can be output by the TGLF quasi-linear model.

These are the three main objectives in building the new flexible GFS moment equations of this paper. Many new applications of these equations present themselves for the future as discussed in Sec. VI.

## II. MANIFESTLY SYMMETRIC FORM OF THE GYROKINETIC EQUATIONS

The electromagnetic gyrokinetic equations with an equilibrium flow common to all species have been transformed into a manifestly Onsager symmetric form by Sugama and Horton (SH).^{15} Their work extends the well-known Onsager symmetry^{16} of the collisional transport system to include the quasi-linear fluxes driven by gyrokinetic turbulence. In the quasi-linear approximation, the formula for the non-linear turbulent fluxes is evaluated with linear gyrokinetic eigenmodes and a model for the saturated fluctuation intensity. The saturation model is constructed from the linear eigenmode spectrum of wavenumbers and growthrates. Hence, the quasi-linear matrix of coefficients of the thermodynamic forces has an implicit non-linear dependence on the local plasma gradients. In this paper, only the special case of axisymmetric toroidal geometry with a toroidal rotation common to all species (ExB) is considered here. The axisymmetric SH equations are used in the CGYRO gyrokinetic code^{17} but with energy and parallel velocity as independent variables and keeping the centrifugal effects at high Mach number. Only sub-thermal velocity toroidal rotation is considered in this paper so the equilibrium density and electric potential are flux functions and there are no centrifugal effects. The SH continuum gyrokinetic equations, with energy and magnetic moment as independent variables, will be put into the form that will be used for generating fluid matrix equations from the velocity moments using the parallel velocity and perpendicular energy as independent variables in Sec. III.

*e*, and magnetic moment

_{is}*μ*through: $ e i s = m i s v \u2032 2 2 T i s = u | | 2 + e i s \u22a5 , \u2009 \mu i s = m i s v \u22a5 \u2032 2 2 T i s B = e i s \u22a5 B$.

_{is}*n*, temperature

_{is}*T*, charge

_{is}*Z*, and mass

_{is}*m*have been introduced.

_{is}## III. GYRO-FLUID SYSTEM (GFS) OF LINEAR GYROKINETIC EQUATIONS

The analytic and numerical evaluation of the transformation of the gyrokinetic equation into a gyro-fluid system of matrix equations in this paper was performed using the Mathematica software.^{18} Orthonormal polynomial sets are generated with exact coefficients and fluid moments of the gyrokinietic equations evaluated analytically. Only the collision coefficients were reduced to numerical values and written into a Fortran data statement form. Taking perpendicular velocity moments of the Bessel functions analytically gives exact functions of the spatially dependent coefficient of the perpendicular energy in the argument of the Bessel functions. This analytic evaluation allows the use of a far fewer perpendicular moments than would be required using a pseudo-spectral Gauss–Laguerre quadrature method.

The Onsager symmetry of the fluxes is directly related to the symmetry of the matrix representation of the various terms. In general, if the continuum term is real then the matrix should be Hermitian. This is particularly important for the mirror force terms. The mirror force should modify the mode frequency not directly drive or damp the equations. If the continuum term is imaginary then the matrix form should be anti-hermitian. The collision operator is the only natural anti-hermitian matrix in the equations. You can verify these properties by examining the symmetry of the inverse of the matrix $ [ G i s ]$.

A Mathematica notebook version of the GFS equations, for simplified s-alpha geometry, was first used to test and refine the organization of the program. This notebook was then used to verify the final Fortran code. In this section, the truncated fluid moment equations will be presented. The closure scheme is given in Sec. IV.

*θ*is replaced with a new coordinate

*z*defined by the integral transformation

*A*are written as

_{is}*is*=

*js*) will be included for the collision operator. The electron collisions are given by

^{19}

*ν*is a function of the electron speed $ u e = e 1 s$,

_{e}The moments of the collision operator are evaluated numerically with Mathematica by substituting $ u 1 s | | = \xi e 1 s$ and $ e 1 s \u22a5 = ( 1 \u2212 \xi 2 ) e 1 s$ into the polynomial expansion of $ H \u0303 1 s$. The resulting matrix is written to a file as a Fortran data statement. The collision operator matrix is symmetric in the velocity space so its contribution to $ G 1 s$ is anti-Hermitian. The eigenvalues of the pitch angle scattering operator are all positive so it is a pure damping term in $ G 1 s$. The like species ion–ion collisions give the same matrix as the electrons without the *Z _{eff}* term and with $ \nu i s , i s = \nu e e m e / m i s ( T e / T i s ) 1.5$.

*ω*being the eigenvalue to be determined. The fluxes and coefficients are evaluated for each of the unstable eigenvlaues and divided by the intensity of the electric potential fluctuation to obtain quasi-linear weights. These are then multiplied by the SAT2 model for the saturated turbulence intensity specturm in TGLF to compute the quasi-linear fluxes that will be compared to non-linear CGYRO results to verify the GFS eigenmodes below.

## IV. RESOLUTION REQUIREMENTS FOR THE TRUNCATED GFS MOMENT SYSTEM

Without any closure, the only dissipation in the truncated GFS moment equations is due to collisions. All of the other matrices contributing to $ [ G i s ] i z , i u , i e , j z , j u , j e$ are Hermitian, so they have real eigenvalues when diagonalized, but the collisions contribute imaginary eigenvalues that are always of a sign to produce damping. It is instructive to compare the most unstable linear eigenmodes from the truncated GFS equations to the initial value eigenmodes of CGYRO. The ion-scale poloidal wavenumber *k _{y}* spectrum for the GASTD case ( $ a / L n = 1 , a / L T = 3 , q = 2 , s \u0302 = 1 , \kappa = 1 , r / a = 0.5 , R / a = 3 , \nu e e = 0.1$) is shown in Fig. 1. All of these cases were chosen to have the same number of parallel energy moments as perpendicular energy moments (ne) so that the total energy conservation of the mirror force and pitch angle collisions is preserved for $ n e tot = n e$ energy moments. This case has a branch jump from the ion temperature gradient (ITG) mode at low

*k*(negative frequency) to the trapped electron mode at higher

_{y}*k*(positive frequency). This jump occurs off-scale for $ k y = 1.1$ for CGYRO. The wavenumber where this jump occurs is sensitive to the velocity space resolution. Notice that the cases with an even number of parallel velocity moments (lower 2 panels) only find the TEM mode at higher velocity resolution compared to the cases with an odd number of parallel moments (upper 2 panels). For the highest resolution (ne = 6), there is little difference between even and odd number of parallel moments and both agree well with the CGYRO eigenvalues. The matrix $ [ u i s | | ] i u , j u$ is symmetric with pairs of positive and negative eigenvalues. When nu is odd there is a zero eigenvalue but not when nu is even. Hence, the odd moment matrix has better coverage of the deeply trapped particles with small parallel velocity. This explains why the odd nu eigenvalues have stronger trapped particle effects (higher growthrate at higher

_{y}*k*and a TEM branch jump) at low resolution than the even nu eigenvalues with the same ne. Note that the growth rates at low

_{y}*k*agree well with CGYRO even for low velocity space resolution ( $ ne \u2265 3$).

_{y}The convergence of the quasi-linear fluxes with perpendicular energy resolution (ne) is shown in Fig. 2. The quasi-linear fluxes were computed using the SAT2 model for the saturated fluctuation intensity^{12,20} that was fit to a database of CGYRO turbulence fluxes using the linear CGYRO eigenmodes. At low resolution, the even nu case (dashed) has a stronger ion energy flux, weaker electron energy flux and more positive electron particle flux than the odd nu case (black). This is consistent with the stronger trapped electron response for the odd nu cases. The odd nu cases for the parallel velocity resolution converge to almost the CGYRO non-linear result (Gray lines). The CGYRO runs used 8 energy, 16 pitch angles and 28 parallel grid points for this case. From this study it is concluded that the truncated GFS requires at least 6 energy, 11 parallel velocity moments, and 12 Hermite moments along the magnetic field. It will be shown that the velocity space resolution requirements for GFS can be reduced with a real closure that preserves the Onsager symmetry of the fluxes in Sec. V.

## V. SYMMETRIC CLOSURE

^{2}a reduced kinetic model with drifts driven by electric potential fluctuations given below is employed,

*nu*= 3, this gives

*nu*> 3, the two closure coefficients are determined by fitting the matrix equation version of the plasma dispersion function and its derivative with respect to

*ζ*to the analytic one Eq. (60) at zero frequency.

*nu*row of the response matrix is not symmetric, in general,

^{2}for

*nu*= 4. Note that the number of parallel moments contained in the background distribution function is typically small (3) so the

*nu*= 4 closure does not violate the Onsager symmetry for the three fluxes in this case. However, the symmetry of the matrix system of equations has been broken, and this could have consequences for the linear or non-linear solutions. Of particular concern would be the non-linear gyrokinetic coupling that contains all of the velocity moments of the fluctuating distribution function.

*nu*= 3, the modified parallel velocity matrix looks like

Attempts to use this symmetric form of the Hammett–Perkins closure with the mirror force contribution (without bounce averaging) yield damping of the linear eigenmodes that is too strong. The delicate interaction of the HP closure with the magnetic drift closures is made even more difficult by the mirror force. The too strong damping of this closure is because all of the particles are assumed to be circulating particles that contribute to the Landau damping of the density response function. The TGLF equations include trapped particles by splitting the velocity moments into trapped and passing regions and assuming that the trapped particles average away the parallel advection term and, hence, are not Landau damped. This way of including trapping is at the cost of a very complex closure once the magnetic drifts are included, and the need to model the loss of bounce averaging at higher frequency. TGLF achieves higher accuracy than previous GLF models but is very difficult to extend to higher velocity space resolution. TGLF does not preserve the generalized Onsager symmetry and does not have a tensor product structure since it has fewer perpendicular energy moments for the higher parallel velocity moment equations.

^{4,21}by introducing complex coefficients to the highest row diagonal and off diagonal terms. It has already been shown that the GFS model will converge to an accurate model of the gyrokinetic theory without any closure (truncation) for high enough resolution, so evidently the complex closure is only needed in order to reduce the resolution required. There is a truncation error for the highest parallel velocity moment that makes the square of the $ [ u i s | | ]$ matrix not agree with the matrix of the velocity squared

^{12}For the resolution $ n e = 3 , n u = 7$, the best fit for this real closure is found for $ c = 0.651 , c 2 = 0.496$. The reduction of the error in the fluxes over the whole SAT2 database with increasing number of Hermite basis functions along the magnetic field (nz) is shown in Fig. 4. The energy flux error is reduced with just nz = 8 but the particle flux continues to improve with nz.

The TGLF-SAT2 with GFS for (ne = 3, nu = 7, nz = 12) and the CGYRO fluxes for all of the database cases are shown in Fig. 5. The largest fluxes are for the scan in r/R that extends into the spherical torus regime r/R = 0.5 (case 44). It is remarkable how well the TGLF-SAT2 fluxes match these very high trapped fraction cases illustrating the accuracy of the mirror force term. The net error for the fluxes with GFS is (Qe, Qi, Ge) = (15.9%, 15.2%, 21.6%). This is better than with the TGLF eigensolver (16.6%, 16.2%, 33.5%) but not as good as SAT2 using CGYRO liner eigenmodes (10.9%, 10.5%, 18.9%).^{12} Note that the TGLF eigensolver used *nz* = 6 but first determines the width of the Guassian z-width that maximizes the nz = 2 linear growthrate before increasing the nz resolution. For GFS. the Gaussian z-width is fixed at 1.74. and a single pass solution with nz = 12 was used.

The same real closure applied to lower parallel resolution $ n e = 3 , n u = 5$ gives a somewhat larger error of (17.6%, 17.1%, 22.2%) for GFS. Going to higher parallel resolution $ n e = 3 , n u = 9$ reduces the error to (14.1%, 15.8%, 20.3%) without adjusting the closure coefficients.

## VI. SUMMARY AND FUTURE APPLICATIONS

The new GFS equations achieved all three of the goals that were the motivation for this work. First, the bounce averaging approximation was not employed for GFS. Instead the mirror force Eq. (50) captures all of the trapping effects. Second, it has the flexibility to increase the velocity space resolution so it can compute the fluxes due to parallel magnetic field fluctuations. Third, the Onsager symmetry of the fluxes is preserved by the formulation and closure used for GFS. The very simple two parameter closure Eqs. (66), (67), and (69) scales with the velocity space resolution. The GFS eigensolver used in TGLF with the SAT2 saturation model was shown to achieve higher fidelity to CGYRO non-linear simulations for the energy and particle fluxes than the original TGLF eigensolver. If the velocity space resolution is (ne = 3, nu = 5), GFS has the same number of moment equations as the original TGLF (15) and the computational speed of the new eigensolver is similar to the original TGLF if the same number of z-basis functions is used. However, since the GFS eigensolver can output the response matrix coefficients of the fluxes, eliminating the need to take numerical derivatives to find the Newton direction, a transport solver can be faster using GFS compared to the original TGLF eigensolver. It is expected that a transport solver utilizing the response matrix output would also converge more reliably to the target fluxes due to the positive definite property of the response matrix.

The flexibility of GFS makes new applications of quasi-linear transport theory possible. The formulation of the GFS system will make it easy to extend to background distribution functions with different parallel and perpendicular temperatures such as fast particle distributions produced by external heating sources, or near material boundaries, where trapped particle orbit losses are possible. By introducing a parallel z-grid, as is used in CGYRO, rather than using Hermite polynomials, extended interchange-like modes could be found with GFS. This would make it easier to resolve micro-tearing and energetic particle driven modes with GFS. Due to the tensor product formulation of GFS, it may be possible to speed up the eigensolver with GPUs and pre-diagonalization of the subspace velocity matrices. Finally, since bounce averaging is not used, it may be possible to use GFS in closed, nested, flux surface geometries of stellarators, allowing for quasi-linear modeling of stellarator turbulent transport. With the new GFS linear eigensolver, the applications of TGLF turbulent transport modeling are greatly expanded.

## ACKNOWLEDGMENTS

Some of the SAT2 database of CGYRO runs^{12} were performed by Jon Kinsey, Nicola Bonanomi, Bhavin Patel, and Harry Dudding. This work benefited from insight by Ron Waltz, Clarrisse Bourdelle, and Jan Weiland. The authors would like to thank the referees for questions that led to improvements of the paper and closure. This work was supported by a grant from the U.S. Department of Energy via Nos. DE-FG02-95ER54309 and DE-AC05-00OR22725.

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**G. M. Staebler:** Conceptualization (lead); Formal analysis (lead); Methodology (lead); Software (equal); Writing – original draft (lead); Writing – review & editing (lead). **Emily Ann Belli:** Data curation (lead); Validation (equal); Writing – review & editing (equal). **Jeff Candy:** Data curation (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Mathematica, Version 13.2*