Gradient- and flux-driven global gyrokinetic simulations of ITG and TEM turbulence with an improved hybrid kinetic electron model

ORB5 is a global gyrokinetic code being developed by many scientists over the last 20 years. It allows performing so-called gradient-driven simulations, in which initial profiles such as density and temperature are specified and an adaptative source is used to relax the zonal component of the perturbed distribution function so that the total flux-surface-averaged distribution function f¯ does not deviate too much from the initial one, f¯0 . Thus, f0 can be used as a control variate to reduce the sampling noise. Even though the gradient-driven model is useful due to its relatively low numerical cost (there is no need to perform simulations over transport time scales) it does not reflect the experimental reality where fluxes are prescribed by the actual heat sources. Furthermore, comparison between experimental and gradient-driven numerical predictions of the heat fluxes is a challenging task because of the profile stiffness: a small change, within experimental error bars, in the imposed temperature gradients leads to a large variation of the computed heat fluxes. To address the stiffness problem, the ORB5 code has been adapted to be able to run flux-driven simulations. In this model, fluxes are prescribed by a radially localized heat source and sink and the background profiles slowly evolve toward a quasi-steady state. As a first step toward complete flux-driven simulations, a mixed approach is used, in which the run is started in the gradient-driven mode providing estimates for the heat sources/sinks, which are then used for pursuing the simulation in flux-driven mode. This allows us to keep the numerical cost relatively low as a quasi-steady state is more quickly reached. In this work, we present the heat source implemented in ORB5 to perform the flux-driven simulations as well as the numerical technique used to ensure that no other moment is injected by the source up to machine precision. Furthermore, simulation results illustrating the mode switching between gradient- and flux-driven will be shown as well as preliminary results of TEM dominated cases done using an upgraded hybrid electron model.


Introduction
Many experimental and theoretical investigations have shown that the plasma density, temperature, and pressure profiles exhibit a stiff behaviour in most of the plasma volume. Stiffness can be defined as the high sensitivity, above a critical gradient, of the heat/particle flux to a small change in the temperature/density gradient, thus limiting e.g. the effect of input power on the profile shape. This is especially true for L-mode discharges made in the TCV tokamak in a Trapped Electron Mode (TEM) regime [1] where the plasma was found to have a stiff region ranging from the sawtooth inversion radius to the near edge (ρ V ∼ 0.8) and a nonstiff region in the edge (0.8 ≤ ρ V ≤ 1), where ρ V = V /V edge is a radial coordinate. In the stiff region, the plasma density and temperature profiles have a constant logarithmic gradient which is little influenced by the value of the input power injected in the system. On the other hand, in the non-stiff region where the gradients are constant, a substantial increase of the "pedestal" temperature with input power was observed. Note that even in the particular context of the Lmode experiments from [1] a "pedestal" location was defined where the transport characteristics changed from constant logarithmic gradient to constant gradient.
Based on this experimental observation, different numerical gyrokinetic studies have been made in both the flux-tube limit [2] as well as using global simulations [3] to better understand the implication of profile stiffness on tokamak performance. In the latter study, simulations in the Ion Temperature Gradient (ITG) regime with adiabatic electrons have been carried out using the global Particle-In-Cell (PIC) code ORB5 to verify if such density and temperature profiles are compatible with the computed transport.
Both numerical studies cited above were done using a gradient-driven approach which makes the profile stiffness investigation rather difficult because a small variation in the profile gradients -within the experimental error bars -leads to drastically different heat and particle fluxes. To address this issue, the ORB5 code has been adapted to run flux-driven simulations in which the fluxes are imposed and not the gradients.
The aim of this work is twofold. First, the flux-driven approach used in ORB5 as well as the different heat sources will be presented together with illustrating simulations. Second, in order to be as close as possible to the experiment, we will present an upgraded hybrid electron model allowing ORB5 to run nonlinear TEM simulations. The original hybrid model was not consistent for nonlinear runs as it did not ensure the ambipolarity condition and did not indeed conserve the toroidal momentum as already pointed out by Idomura [4].
The remainder of this paper is organized as follows. In Section 2 the ORB5 code is described with an emphasis on its hybrid electron model, its heat source and Krook noise control. In Section 3, flux-driven simulations of an ITG dominant case as well as gradient-driven TEM simulations using the new upgraded hybrid electrons model are presented. Finally, conclusions are made in Section 4.

The ORB5 code
ORB5 [5,6,7] is a global, multi-species, electromagnetic gyrokinetic code using a Lagrangian PIC method with a finite-element Galerkin method and a B-splines representation of the electromagnetic fields. The ideal MHD magnetic equilibrium is obtained from the CHEASE code [8] which solves the Grad-Shafranov equation. The Vlasov-Maxwell system of equations is obtained from variational principles [9,10], allowing one to include all the necessary approximations in the gyrokinetic action before deriving the equations of motion. This ensures, from the theoretical point of view, exact conservation of quantities such as energy and toroidal momentum.

Poisson equation and hybrid electron model
In the ORB5 code, the Quasi-Neutrality Equation (QNE) is used to solve for the self-consistent perturbed electrostatic potential δφ. Limiting our discussion to a single ion species for the sake of simplicity, the original hybrid electron model can be written as where q i and q e are the ion and electron charge, and δn pol i and δn gyr i are respectively the ion polarization and gyro densities. The second term represents the passing electron adiabatic response with where α P is the local fraction of passing electrons, n e,0 is the background electron density, T e,0 is the background electron temperature, s = ψ/ψ e is the normalized radial variable where ψ and ψ e are respectively the poloidal flux and its value at the edge, and the overbar symbol stands for the flux-surface average defined as where J(s, θ ) is the Jacobian of the coordinates (s, θ , ϕ) with θ the straight-field-line poloidal angle and ϕ the toroidal angle. Finally, the last term is the trapped electron drift kinetic density given by where d Z is the infinitesimal phase-space volume, R is the electron guiding center and the integral is done over the portion of the electron phase-space representing trapped electrons. Note that the code can also run with a full adiabatic electron response (α P = 1 and δn drift e,T = 0): or a full drift-kinetic electron dynamics (α P = 0 and δn drift e,T computed over the whole phase space): The hybrid model, Eq. (1), is inconsistent especially for nonlinear runs as it does not ensure the ambipolarity condition and toroidal momentum conservation because (i) no kinetic response of the passing electrons is considered and (ii) trapping-detrapping processes due to parallel and E × B nonlinearities and collisions lead to a fraction of the particle density to be not well accounted for because from the Poisson equation point of view there are sinks/sources of particles. Indeed, recently trapped particles suddenly appear in the kinetic response while recently detrapped particles disappear.
To overcome this issue a new hybrid model was proposed by Idomura [4] in which part of the passing electron response is treated as kinetic, taking great care to avoid the so-called ω H mode [11] -the electrostatic limit of the kinetic Alfvén wave -that would require to drastically limit the simulation timestep to resolve its frequency, i.e ω H = (k /k ⊥ ) m i /m e Ω ci , where k and k ⊥ are respectively the parallel and perpendicular mode numbers, m i /m e is the ion to electron mass ratio and Ω ci is the ion cyclotron frequency. This new hybrid model proved to solve the ambipolarity condition and conserve the toroidal momentum. However, linear studies showed that the Geodesic Acoustic Mode (GAM) frequency and damping rate are significantly affected by this model.
In view of possible comparisons between simulations and TCV experiments it is needed to correctly compute the GAM properties. To this end, we developed an upgraded hybrid electron model which is a further modification to the one presented in [4]. The rationale of this upgraded model is to separate the zonal and non-zonal contributions of the passing electrons and account for them differently. Noting n and m the toroidal and poloidal mode numbers, the model can be summarized as follows: For the non-zonal contribution (n = 0), the approach is thus the same as the original hybrid model. For the n = 0 case, the electron contribution to the m = 0 modes is treated as adiabatic while the contribution to the m = 0 mode is treated as drift-kinetic. Note that all the components k = 0 of the kinetic passing electron dynamics are removed from the system, thus avoiding the appearance of the ω H mode. As compared to the model of Ref. [4], this upgraded model also satisfies the ambipolarity condition but the GAM frequency is correctly captured.

Krook noise control and heat source
The inherent noise of the PIC method is reduced by two means. First, the code employs a control variate technique (δf -PIC method) which consists in separating the full distribution function f into a known background part f 0 and a perturbed part δf and only evolving the latter using numerical particles. The main requirement of this technique to efficiently reduce noise is that δf f 0 . Second, various noise control methods are implemented in the code, e.g. a modified Krook operator, coarse graining procedure, and a quadtree smoothing routine. For this work, only the Krook operator was used. For the sake of completeness, we briefly describe this noise control approach but a more detailed discussion can be found in [12].
The modified Krook operator is a source term, S K , which is added on the right-hand side of the Vlasov equation, i.e. df dt and has the form where γ K is the relaxation rate, typically γ K ∼ γ max /10, with γ max the maximum linear growth rate. It is crucial to ensure that the Krook operator does not modify too much the transport properties of the system. To this end, a correction term, S K corr , is added in order for the Krook operator to conserve various transport related quantities such as the density n, the parallel momentum v , the kinetic energy E, and the Zonal Flows (ZF) playing an important role in the regulation of transport. The correction term is defined by where the density, kinetic energy, parallel velocity and ZF respectively, R, and μ are the gyrocenter position, and the magnetic moment of the particle. The coefficient g i are computed such that Inserting Eqs. (6) and (7) into Eq. (8) leads to the following linear system of equations which is solved at every timestep to find the coefficients g i : Note that not conserving E allows running temperature-gradient-driven simulations, while conserving n, v , and the ZF allows the density and flow profiles to freely evolve far from their initial ones.
To perform temperature flux-driven simulations we need to inject heat in a localized and prescribed way which is done via a heat source defined by where γ H (s) is a radial deposition rate profile and T is the temperature. As for the modified Krook operator, the heat source is added on the right-hand side of the Vlasov equation, Eq.
(5). The functional form of S H ensures that, for an infinite number of particles, only heat will be injected and not particles or momentum. From the numerical point of view, it is impossible to have an infinite number of markers. To make sure that the heat source conserves exactly the other moments (n, v , and zonal flows) up to machine precision, a correction term similar to the Krook operator is added. Note that for simulations with multi kinetic species both the Krook operator and the heat source are applied successively and independently to each species.

Numerical simulations
All the simulations presented in this paper are based on a similar TCV reference case reflecting the experimental conditions met in [1]. More precisely we reproduce the same functional dependence of the density and temperature profiles without necessarily having the same experimental powers. The background temperature profile is defined by where T 0 , T 1 , κ T , μ T , and ρ ped = 0.8 are input parameters, T ped = T 1 + μ T (1 − ρ ped ). Here, the ped subscript refers to the "pedestal" values. The density profile is defined in a similar way using the parameters N 0 , N 1 , κ N , and μ N . The magnetic configuration is computed by the CHEASE code based on a TCV equilibrium and is an axisymmetric ideal MHD equilibrium with an aspect ratio of 3.64, an elongation of 1.44, and a triangularity of 0.20. The safety factor profile is monotonic with q 0 = 0.78 and q a = 3.29. The reference magnetic surface for the normalization is s 0 = 1, the reference species is the ions, and ρ (s 0 ) = ρ i (s 0 )/a = 1/245.
All the runs are made with a grid resolution of N s × N θ × N ϕ = 256 × 512 × 256. The field-aligned Fourier filter is set to keep only modes such that n = 0 to n = 64 and satisfying m = nq(s) ± 5. Finally, the Krook rate is set to around 10% of the maximum linear growth rate, i.e. γ K = 1.4 · 10 −4 [Ω ci ].

Flux-driven ITG simulations
The ITG simulations are done using adiabatic electrons, Eq.
The electrons have the same temperature and density profiles as the ions. The numerical resolution for these runs is N p = 128 · 10 6 markers and Δt = 20 [Ω −1 ci ]. In general, flux-driven simulations tend to be very expensive as compared to gradient-driven simulations because the plasma has to relax toward an equilibrium state. Usually this relaxation process takes a time of the order of the transport time scale that increases with the machine size. Furthermore, the final state can be very different than the initial state, which is problematic for PIC codes using control variate without an adaptative background scheme, e.g. [13], since δf can then become comparable to f 0 .
To mitigate these two problems we used a mixed approach for our simulations. First, the simulation is initiated in gradient-driven mode with the Krook operator acting as noise control and heat source (conserving n, v , and ZF). As already mentioned, the temperature and density profiles have the same functional dependence as the experiment but the κ and μ values have to be adjusted so that we obtain localized heat source/sink profiles from the modified Krook operator. Even for a comparison with the experiment, different values of κ and μ within the measurement error bars would lead to drastically different heat powers due to the profile stiffness. The effective heat source profile associated to the Krook operator is then time-averaged and radially localized using two Gaussian shapes, one for the actual heat source and one for the negative heat sink, Figure 1. The heat sink is computed so that the total injected power is zero and placed close to the radial edge but not to close to avoid problems with finite sampling and boundary conditions. This source and sink localization using analytic shapes is done to obtain source/sink-free region where turbulence can develop freely.   With the so-obtained radial heat deposition profile we continue the gradient-driven run in flux-driven mode. This is done by setting the Krook operator to conserve n, v , E, and the ZFit acts as noise control only -and turning on the heat source using the radial deposition profile. Figure 2 shows a proof of principle of this mixed method. We plot the time evolution of the effective heat diffusivity: where q H is the heat flux and ρ is any radial coordinate. Furthermore, the effective heat diffusivity is normalized to the Gyro-Bohm diffusivity χ GB = ρ 2 s0 c s0 /a, where ρ s0 is the ion sound Larmor radius, c s0 is the ion sound speed, and a is the minor radius on the equatorial plane. In Fig. 2, the green line shows the result for the gradient-driven run used to build the radial deposition profile for the heat source. At the end of this run (dashed line) the run is resumed with the flux-driven approach showing a smooth transition and a convergence toward a similar quasi-steady state. To show the robustness of this method we redo a similar simulation but entirely with the flux-driven method and using the same radial deposition profile for the heat source. The first few overshoots are slightly different than for the gradient-driven simulation but it converges to a similar quasi-steady state.

TEM simulations using the upgraded hybrid electron model
The TEM simulations are done using the upgraded hybrid electron model, Eqs. (3)-(4). The considered ion and electron temperature and density profiles are respectively defined by T i 0 = 5, ]. First, we present a few properties of the upgraded model against the original hybrid model. A more detailed study will be exposed in a future work. Figure 3 shows a comparison between the original hybrid model and the upgraded version.
On the left plot we show the time evolution of the radially averaged effective heat diffusivity. With the original hybrid model, there is a clear distinction between ion (blue) and electron (orange) heat transport. The case being TEM dominant, the electron diffusivity is around 2.2 times higher than the ion diffusivity for late times (t > 3 · 10 4 Ω −1 ci ). With the upgraded hybrid model, the electron diffusivity (red) is not much changed and roughly follows the same trend as with the original model. On the other hand, the ion diffusivity (green) slowly grows during a transition period (t ∈ [∼ 0.5, 3] · 10 4 Ω −1 ci ) until reaching the electron level. In this case, the turbulence regime appears to be not a pure TEM anymore but a mixed ITG/TEM regime.
The right plot shows the ion and electron density profiles in the radial domain s ∈ [0.5, 1.0]. The original hybrid model did not ensure the ambipolarity condition as the ion (blue) and electron (orange) density profiles were totally different in this region. On the other hand, with the upgraded version both ion (green) and electron (red) density profiles match closely. Surprisingly, for the considered TEM simulation, it is the ion density profile that is the most affected by this corrected model.
On Figure 4, we show the effective heat diffusivities as a function of the effective density gradient R/L n for two simulations with respectively κ N = 2.3, μ N = 5.0 and κ N = 1.3, μ N = 6.1. The first is the TEM case discussed until now and the second represents a new run with an initial background density profile matching the relaxed density profile at the end of the first simulation. This was done to continue the initial simulation, for which δf was around 20% of f 0 , with optimal noise conditions. The two runs match perfectly for similar R/L n values, showing a certain consistency with this continuation method. In this case it is clear that the density gradient is directly related to the increase of the ion diffusivity. One first hypothesis is that with a lower density gradient at same temperature gradient the η i ratio (density gradient length to temperature gradient length ratio) is higher and thus the ITG is more linearly destabilized.
This hypothesis is clearly contradicted by the linear dispersion relation shown in Figure 5. Indeed, for the base case with the highest density gradient around 30% of the total growth rate comes from the ions (except for n = 4, 8, 12 that are mostly stable). On the other hand, with the smaller density gradient there is no significant contribution from the ions to the growth rate and the total linear growth rates comes mainly from the electrons.
Another hypothesis to explain the increase in ion heat diffusivity is that the Zonal Flows (ZF) amplitude reduces with the density gradient. Indeed, it is known that ZFs are responsible for quenching ITG turbulence and reducing the subsequent transport. For the TEM turbulence   the situation is not so clear and depends on various parameters [14]. In Figure 6 the dependence of E × B shearing rate on the effective heat diffusivity is shown. On the left plot we show the ion heat diffusivity. It is clear that as the density gradient decreases, the E × B shearing rate decreases also thus increasing the ion heat transport as expected. However, for the electrons (right plot) the heat diffusivity is not significantly affected by the decrease of E × B shearing rate.

Conclusions
We have tested a mixed procedure to run flux-driven simulations as a continuation of gradientdriven runs. The switching between gradient-and flux-driven phases appears to be smooth  and we have tested that the same quasi-steady states are reached in both cases. To show the robustness of the method, a flux-driven only simulation is performed showing again the same quasi-steady state.
To correct the original hybrid electron model that was implemented in ORB5 but not consistent from the point of view of ambipolarity condition and toroidal momentum conservation, an upgraded hybrid electron model has been developed to run nonlinear TEM simulations. We found an increased ion heat transport that is explained by the reduction of the E × B shearing rate due to a decrease in the density gradient.