Padé approximation of the adiabatic electron contribution to the gyrokinetic quasi-neutrality equation in the ORB5 code

This work aims at completing the implementation of a solver for the quasineutrality equation using a Padé approximation in the global gyrokinetic code ORB5. Initially [Dominski, Ph.D. thesis, 2016], the Pade approximation was only implemented for the kinetic electron model. To enable runs with adiabatic or hybrid electron models while using a Pade approximation to the polarization response, the adiabatic response term of the quasi-neutrality equation must be consistently modified. It is shown that the Pade solver is in good agreement with the arbitrary wavelength solver of ORB5 [Dominski, Ph.D. thesis, 2016]. To perform this verification, the linear dispersion relation of an ITG-TEM transition is computed for both solvers and the linear growth rates and frequencies are compared.


Introduction
It is well known that so-called anomalous transport due to micro-turbulence plays a dominant role in the particle and heat transport across the magnetic field lines of a magnetically confined plasma. This turbulence is mainly caused by small-scale instabilities such as the Ion Temperature Gradient (ITG) mode, the Trapped Electron Mode (TEM), and the Electron Temperature Gradient (ETG) mode that are driven by the plasma density and temperature gradients.
To describe efficiently these instabilities that typically evolve on slow time scales as compared to the ion cyclotron frequency, the gyrokinetic framework [1] is used. In this model, the fast gyro motion is averaged out, thus reducing the dimensionality of the phase space from 6D to 5D and eliminating fast time scales, resulting in a problem of decreased complexity. The distribution functions of each plasma species are evolved with the reduced set of gyrokinetic variables and the electrostatic field is solved consistently with the gyrokinetic Quasi-Neutrality Equation (QNE).
In the global gyrokinetic code ORB5, the field solver was originally using the Long Wavelength Approximation (LWA) for the polarization density [2,3]. Recently, an Arbitrary Wavelength Solver (AWS) has been implemented based on previous works, see references [4,5]. A preliminary version of a Padé solver has also been implemented which was only valid with the kinetic electron model [6,7].
The Padé scheme involves approximating the polarization response with the inverse of a factor that includes the Laplacian; to make the resulting quasi-neutrality equation explicit, the whole equation is multiplied by this prefactor. The term in the quasi-neutrality equation involving adiabatic electron response must therefore also be consistently modified in this framework.
In this work, we have completed the former Padé solver by adding correction terms to the adiabatic electron response, allowing us to simulate perturbations down to small scales using the adiabatic electron model or a "hybrid" model in which trapped electrons are kinetic while passing electrons are assumed adiabatic. The complete Padé solver is of practical interest because it produces smaller matrices in terms of non-zero elements than the AWS which are, in turn, faster to assemble. As an example, for the runs performed during this work, the Padé solver required 1.5 s to assemble the matrix whereas the AWS took 145 s. Furthermore, the Padé matrix needed ∼ 9 times less memory to be stored than the AWS matrix. In the present version of the ORB5 code, the QNE is linearized, and the matrix needs to be computed only once. However, in view of possible extensions of the code, this performance gain can become crucial. For example, as we approach the tokamak edge region, the perturbation amplitude becomes comparable to the background equilibrium and the non-linear terms of the QNE need to be computed. This would require assembling the matrix at every time step of the simulation which is not practical with the AWS due to its long assembly time.
This paper is organized as follows. Section 2 describes the gyrokinetic QNE, as well as its long wavelength and Padé approximations. In Section 3, we present the discretization of the QNE and its implementation in the code. Finally, the code verification in the linear regime is presented in Section 4.

Computational Model
The linearized QNE for a two-species, electrostatic plasma can be written and δn e = δn gyr e + δn ad e . Here, n σ = n 0σ + δn σ with n 0σ being the background density, δn σ the perturbed density. The background densities are assumed to satisfy the QNE, i.e. q i n i0 + q e n e0 = 0.
In Eq. (1), δn ad e is the adiabatic contribution to the perturbed electron density given by where x = (s, θ , ϕ) with s, θ , and ϕ being respectively the radial, poloidal and toroidal coordinates, T e,0 is the equilibrium electron temperature, δφ is the flux-surface averaged perturbed potential, and α ad is the fraction of adiabatic electrons (α ad = 1 for the adiabatic electron model, α ad = 0 for the fully kinetic electron model, and α ad is equal to the fraction of passing electrons for the hybrid model). Both δn gyr i and δn gyr e are the ion and electron gyro densities: where δf σ is the perturbed distribution function, J the Jacobian relative to the transformation from particle to guiding-center variables, R the guiding-center position, ρ i the ion Larmor radius, and finally d Z the infinitesimal phase-space volume. Note that, in the last integral of Eq. (3), only the phase-space region of non-adiabatic electrons is considered. Finally, the gyrokinetic ion polarization density is given by where B 0 is the equilibrium magnetic field amplitude, and δφ = δφ − δφ α with δφ α being the gyro-averaged electrostatic potential. In ORB5, the quasi-neutrality equation (1) with ion polarization density given by (4) is solved in its integral form by the AWS recently implemented, see references [6,7]. This upgrade was added in order to overcome the original LWA made in the code, in which the ion polarization density contribution is given by where Ω i,0 = q i B 0 /m i and m i is the ion mass. In the Padé approximation [4,8], one considers the improved approximation At this point, all the QNE is multiplied by In Eq. (7), the first term is the LWA of the ion polarization density, the second term is the adiabatic electron density contribution and the third term is the Padé correction to the adiabatic electron density contribution. In Eq. (8), the first and third terms are respectively the ion and electron gyro densities, and the second and third terms are the Padé corrections to these gyro densities.
The first version of the Padé solver, already implemented in ORB5, was actually only valid for the kinetic electron model [6]. In the present work, we are thus interested in implementing the Padé correction terms related to the adiabatic electron contribution to the QNE. These corrections terms are necessary for both the adiabatic and the hybrid electron models.

Discretization of the QNE equation
In ORB5, the QNE is discretized using a finite-element representation of the electrostatic field where Λ ν (s, θ , ϕ) = Λ i (s)Λ j (θ )Λ l (ϕ) are a tensor product of one-dimensional B-splines basis functions [2,9]. With this representation the QNE is discretized into a linear system of equations following the Ritz-Galerkin method: The last term of this equation, depending on the equilibrium profile gradient, is not considered in this work. It can be shown that it is ε = ρ i /L T,n smaller than the first term, and can therefore be neglected. Because of the toroidal and poloidal periodicity, the discretized QNE can be Fourier transformed using a DFT. The advantage of this procedure is multiple. First, as the system is toroidally axisymmetric, the toroidal mode numbers n are decoupled [2]: Furthermore, based on a work presented in Ref. [10], only the relevant poloidal modes are assembled. Indeed, since the turbulence of interest is principally aligned with the magnetic field, then for each toroidal mode number n only a few poloidal mode numbers m = nq ± ∆m need to be retained (for most simulations ∆m = 5 is enough to converge the results), where q is the safety factor. This can be used during the assembly in order the reduce the size of the matrix and thus, speeding up the computations.

Code Verification in the linear regime
The new version of the Padé solver is verified in the linear regime against the AWS which was benchmarked [6] against the GENE code [11]. To this end, a modified CYCLONE test case [12] is whereρ = ρ/a is a normalized radial coordinate,ρ 0σ = 0.5, ∆ Tσ = ∆ n = 0.30, T i = T e , L n /L T i = 3.14 is the density to temperature characteristic length ratio, and R 0 /L T i = 6.92. The numerical parameters are the following: N M = 2 · 10 6 is the number of numerical markers, ∆t = 16Ω −1 i is the time step, n s = 120, n θ = 512 and n ϕ = 128 are the number of grid intervals in the s, θ and ϕ directions, the simulation domain is annular from s = 0.3 to s = 1 and the hybrid electron model is used.
The simulation results are plotted in Figure 1, where the frequency and growth rate are shown as a function of k θ ρ i , where k θ is the poloidal mode wave number.
The left plot shows the linear growth rates of the most unstable modes. The Padé solver (solid green line) is in very good agreement with the AWS (solid blue line). For values k θ ρ i > 0.6, the Padé solver slightly overestimates the growth rate but the difference is less than 3.5%. To emphasize the importance of the Padé correction to the adiabatic electron contribution, a scan is made neglecting L ad,corr νν , i.e. only the drift-kinetic electrons have the Padé correction. It is seen that the growth rate is always overestimated, especially for k θ ρ i > 0.4. The right plot of Figure  1 shows the linear frequencies of the modes for a similar scan in k θ ρ i . The Padé solver is also in agreement with the AWS. More surprisingly, the Padé correction to the adiabatic electron contribution seems to be negligible for the estimation of the frequency.
In Figure 2, we represent a poloidal cut for the three solvers: Long Wavelength Solver (LWS), AWS and Padé. For all cases, the cut shows the electrostatic potential and is taken at t = 8000 Ω −1 ci and for k θ ρ i = 0.62. The perturbation computed by the Padé solver (right) has a slightly higher amplitude due to the growth rate that is slightly overestimated by Padé as compared to the AWS (center) (γ P ade = 0.27 c s /R 0 and γ AW S = 0.26 c s /R 0 ). On the top-left of each figure is shown a zoom of the poloidal mode structure. Both Padé solver and AWS resolve correctly the fine structures (radial ripples in the zoomed pictures). On the other hand, the LWS (left) fails at solving the fine structures and the amplitude of the perturbation is around ten times smaller than the AWS because the LWS largely underestimates the growth rate: we are beyond the validity limit of the long wavelength approximation.
Thanks to this new adiabatic electron correction term, the Padé solver allows one to solve for short-scales instabilities that arise around mode rational flux surfaces when using both the adiabatic and hybrid electron models.  Figure 2. Poloidal cuts for t = 8000 Ω −1 ci and k θ ρ i = 0.62 for all three cases. The LWS is on the left, the AWS in the center and the Padé version of the solver on the right.

Conclusions
The new adiabatic electron correction term for the Padé version of the solver has been successfully implemented in ORB5 for n = 0 mode numbers and allows one to accurately simulate the ITG-TEM transition. Developments related to the n = 0 modes are under progress, which are necessary for the simulation of turbulence with the adiabatic electron approximation when using the Padé version of the solver. Other tests (e.g. Rosenbluth-Hinton [13]) are being carried out to finalize the verification of the Padé solver in the linear regime.
Linear dispersion relation (growth rate and frequency) shows that Padé solver has a good agreement with the AWS. For k θ ρ i > 0.6, Padé solver slightly overestimates the growth rate as compared to AWS. Nevertheless, the difference is always smaller than 3.5%. The contribution of the Padé correction to the adiabatic electron response is shown to be non-negligible for the growth rate estimation. The mode frequencies obtained with both Padé solver and AWS are in good agreement even without the Padé correction to the adiabatic response.
Finally, we have shown that the Padé solver, as for the AWS, is capable of solving the fine structures that arise around mode rational surfaces.