2.1 The GENE and QuaLiKiz simulations

The simulations conducted in this work are based on the gyrokinetic code GENE and the quasilinear gyrokinetic code QuaLiKiz. We use GENE to solve the linearized gyrokinetic equations in the electrostatic limit in a flux tube; this limit is justified for the low-$\beta$ tokamak cases analysed in this work, where $\beta$ is the ratio of plasma pressure to magnetic pressure. For a given species $s$, the distribution $f_s$ is split into a background Maxwellian $f_{0s}$ and a fluctuating part $f_{1s}$ such that $f_s = f_{0s} + f_{1s}$. The electrostatic potential is self-consistently solved via the quasineutrality equation. In this work, we use both the initial value solver and the eigenvalue solver versions of GENE. More details can be found in Jenko et al. (Reference Jenko, Dorland, Kotschenreuther and Rogers2000).

In contrast, QuaLiKiz utilizes a variational approach and is based on the action-angle formalism. Using this formalism, the kinetic equation is easily linearized and then coupled to the quasineutrality equation. Therefore, instead of working with the distribution function directly, we multiply the quasineutrality equation by the complex conjugate of the electrostatic potential and integrate over the domain, thus obtaining a zero-dimensional dispersion relation instead of a differential equation.

Note that, in this work, we have neglected toroidal rotation and radial electric shear. In Tala et al. (Reference Tala2019), it was found that finite rotation only has a minor impact on the density peaking, so we neglect these effects to simplify the analysis. Once the mode frequency is found, we use a quasilinear saturation rule to compute the turbulent particle, angular momentum and heat fluxes. We take the functional form of the eigenfunction to be a Gaussian in ballooning space, where the parameters of the Gaussian are determined via a high-frequency fluid approximation. As described in Cottier et al. (Reference Cottier, Bourdelle, Camenen, Gürcan, Casson, Garbet, Hennequin and Tala2014), the Gaussian ansatz agrees well with linear gyrokinetic simulations near the peak of the eigenfunction, but TEM eigenfunctions tend to be more extended in ballooning space compared with a Gaussian. Assuming a shifted-circular geometry, we then obtain tractable two-dimensional integrals that can be performed numerically, where passing and trapped particles are treated separately. More details as to the particularities of the derivation can be found in Stephens et al. (Reference Stephens, Garbet, Citrin, Bourdelle, van de Plassche and Jenko2021).

2.2 Collision operators

In this work, we make use of two collision operators: the Landau–Boltzmann collision operator (Hazeltine & Waelbroeck Reference Hazeltine and Waelbroeck2004) in GENE and a Krook-like operator (Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995) in QuaLiKiz. In general, a collision operator enters the Boltzmann equation as

(2.1)\begin{equation} \frac{{\rm d} f_s}{{\rm d} t} = \left(\frac{\partial f_s}{\partial t}\right)_{\text{coll}} = \sum_{s'} C_{ss'}, \end{equation}

where $f_s$ is the distribution function of species $s$ and the collision operator $C_{ss'}$ takes into account collisions between different species. The full Landau–Boltzmann operator is written as

(2.2)\begin{equation} C_{ss'} = C_{ss'}(f_s, f_s') = \frac{\gamma_{ss'}}{m_s} \frac{\partial}{\partial \boldsymbol{v}} \boldsymbol{\cdot} \int {\rm d}^{3} {v'} \left(\frac{u^2 {\boldsymbol{\mathsf{I}}} - \boldsymbol{u}\boldsymbol{u}}{u^3}\right) \boldsymbol{\cdot} \left(\frac{f_{s'}'}{m_s} \frac{\partial f_s}{\partial \boldsymbol{v}} - \frac{f_s}{m_{s'}} \frac{\partial f_{s'}'} {\partial \boldsymbol{v}'}\right), \end{equation}

where we use notation such that $f = f(\boldsymbol {v})$ while $f' = f(\boldsymbol {v}')$ and $\boldsymbol {v}$ is the velocity vector (Hazeltine & Waelbroeck Reference Hazeltine and Waelbroeck2004). We also use standard dyadic notation with $\boldsymbol *{I}$ being the three-dimensional identity matrix and define the relative velocity $\boldsymbol {u} = \boldsymbol {v} - \boldsymbol {v}'$. Meanwhile, the factor $\gamma _{ss'}$ is written as

(2.3)\begin{equation} \gamma_{ss'} = \frac{e^4 Z_{s}^2 Z_{s'}^2 \lambda_{ss'}}{8 {\rm \pi}\epsilon_0^2}, \end{equation}

where $m$ and $e$ are the mass of the indicated species respectively, $\epsilon _0$ is the vacuum permittivity and $\lambda = \ln (\varLambda )$ is the Coulomb logarithm. Because the simulations in GENE are linear, the collision operator is likewise linearized via the $\delta f$ approximation, as described in Crandall (Reference Crandall2019).

In QuaLiKiz, a different approach is used. Since any collision operator that uses derivatives would slow down QuaLiKiz considerably, we require a simpler approach to obtain a tractable dispersion relation. Energy-dependent collision operators, currently used in QuaLiKiz, already increase the computational complexity since Fried and Conte integrals can no longer be used to analytically simplify calculations of collisional species. Thus, any additional complexity in the operator must be carefully considered to avoid a computationally intractable code. We therefore make use of a Krook-style operator (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954; DeLucia & Rewoldt Reference DeLucia and Rewoldt1981) constructed to mimic the qualities of a pitch-angle operator. First, we neglect ion–ion collisions completely since they play little role for TEMs. Electron–electron collisions are also ignored as they provide only a small correction to the collision operator. Finally, this collision operator is only used for trapped particles, meaning that particle number, momentum and energy are not conserved. Although an artificial collision term could in principle be added to passing particles to ensure conservation properties, it has been shown that these terms are negligible corrections in the electrostatic limit (Rewoldt, Tang & Hastie Reference Rewoldt, Tang and Hastie1986). Since QuaLiKiz is an electrostatic code, these conservation properties can be safely ignored. To verify this, a brief implementation of passing electron collisions did not alter calculated growth and was thus discarded. Consequently, marginally passing particles deflecting into the trapped part of velocity space tend to have a weak effect on TEM growth rates for the considered parameter regimes. We then write for electrons

(2.4)\begin{equation} \left(\frac{\partial f_e}{\partial t}\right)_{\text{coll}} = \left\{\begin{array}{@{}ll@{}} - \nu_e \left(\delta f_e - \dfrac{e \phi}{T_e} f_{0e}\right) & \text{if trapped}, \\ 0 & \text{if passing}, \end{array}\right. \end{equation}

where $\phi$ is the perturbed electrostatic field and $T$ is the temperature (Romanelli, Regnoli & Bourdelle Reference Romanelli, Regnoli and Bourdelle2007). Meanwhile, because we have no collision operator for ions, we write

(2.5)\begin{equation} \left(\frac{\partial f_i}{\partial t}\right)_{\text{coll}} = 0. \end{equation}

Collisions thus drive the $\delta f_e$ to $e \phi / T_e$, which corresponds to the electron adiabatic response created by the electrostatic perturbation.

We proceed with the construction of the Krook-style collision operator, following the analysis presented by DeLucia & Rewoldt (Reference DeLucia and Rewoldt1981) and Rewoldt, Tang & Hastie (Reference Rewoldt, Tang and Hastie1987). For simplicity, we consider only one ion species; the final result can be generalized for multiple ion species by replacing $Z_i$ with $Z_\text {eff}$ for electron collisions and a sum over all ions for ion collisions. The key is to construct two separate preliminary collision operators for electrons and ions with collision frequencies

(2.6)\begin{gather} \nu_e = a_e \epsilon \frac{\nu_{ei}}{(\hat{v}_e)^3} \frac{Z_i + H(v / v_{T_{e}})}{|1 - 2 \epsilon - \lambda|^2}, \end{gather}

(2.7)\begin{gather}\nu_i = a_i \epsilon \frac{\nu_{ii}}{(\hat{v}_i)^3} \frac{H(\hat{v}_i)}{|1 - 2 \epsilon - \lambda|^2}. \end{gather}

Here, $\epsilon = r /R_0$ is the inverse-aspect ratio where $r$ is the local minor radius and $R_0$ is the major radius, $\lambda$ is the pitch angle parameter, $\nu _{ss'}$ denotes the characteristic collision frequency between species $s$ and $s'$, $\hat {v}_s = v / v_{\text {th}, s}$ is the speed normalized to the thermal velocity $v_{\text {th}, s} = \sqrt {2 T_s/m_s}$ and $a_s$ is a constant. Meanwhile the function $H$ is defined as

(2.8)\begin{equation} H(x) = \frac{{\rm e}^{{-}x^2}}{\sqrt{\rm \pi} x} + \left(1 - \frac{1}{2x^2} \right) \text{erf}(x). \end{equation}

In QuaLiKiz, we use the following definition for the pitch-angle parameter $\lambda$:

(2.9)\begin{equation} \lambda = \frac{\mu B_{\text{min}}}{E}, \end{equation}

where $\mu$ is the magnetic moment, $B_{\text {min}}$ is the magnetic field strength evaluated at its minimum on a given flux surface and $E$ is the kinetic energy. Because QuaLiKiz assumes circular flux surfaces, for small $\epsilon$, the magnetic field strength can be approximated as $B \approx B_0 (1 - \epsilon \cos (\theta ))$, where $\theta$ is the poloidal angle and $B_0$ is a constant characteristic field strength. Thus, $B_{\text {min}} = B_0 (1 - \epsilon )$. The pitch-angle parameter determines whether a given particle is either trapped ($\lambda > 1 - 2 \epsilon$) or passing ($\lambda < 1 - 2 \epsilon$). The trapped–passing boundary is given at $\lambda = 1- 2\epsilon$.Footnote ¹ The pitch-angle dependence in the collision operator mimics the effects of pitch-angle scattering in a crude way without the use of differential equations. Meanwhile, the collision frequency is constructed to diverge at the trapped–passing boundary. This accounts for enhancement of collisions at this boundary in velocity space. It is important for any implementation of this method to correctly captures this divergence. We also note for deeply trapped, thermal particles that

(2.10)\begin{gather} \nu_{e} \sim \frac{Z \nu_{ei}}{\epsilon}, \end{gather}

(2.11)\begin{gather}\nu_{i} \sim \frac{\nu_{ii}}{\epsilon}. \end{gather}

The collision frequencies are constructed such that they mimic the effective collision frequency for a simple Krook operator ($\nu _{\text {eff}} = \nu / \epsilon$). This takes into account the fact that the typical scattering process we are concerned with is not $90\,^{\circ }$ scattering but rather diffusion from the trapped part of velocity space to the passing part of velocity space.

The coefficients $a_e$ and $a_i$ are determined by comparing the collisional models with a dispersion relation obtained using a Lorentz collision operation by means of the variational principle (Rosenbluth, Ross & Kostomarov Reference Rosenbluth, Ross and Kostomarov1972). It is important to note that this is calculated in the regime where

(2.12)\begin{equation} \frac{\nu_{ii}}{\epsilon} \ll |\omega| \ll \frac{\nu_{ei}}{\epsilon}, \end{equation}

where $\omega$ is the complex mode frequency. DeLucia then writes (DeLucia & Rewoldt Reference DeLucia and Rewoldt1981)

(2.13)\begin{gather} a_e = 1.31, \end{gather}

(2.14)\begin{gather}a_i = 9.42 \times 10^{{-}3}. \end{gather}

These coefficients are chosen such that they reproduce the collisional scalings obtained from the more accurate Lorentz model. The next step is to construct an operator that is valid outside this limit and bridges the regime of low collisionality and high collisionality (relative to the mode frequency). The result is

(2.15)\begin{gather} \nu_e = \epsilon \frac{\nu_{ei}}{\hat{v}_e^3} \frac{Z_i + H(\hat{v}_e)}{|1 - 2 \epsilon - \lambda|^2} \frac{0.111 \delta + 1.31}{11.79 \delta + 1}, \end{gather}

(2.16)\begin{gather}\nu_i = \epsilon \frac{\nu_{ii}}{\hat{v}_i^3} \frac{H(\hat{v}_i)}{|1 - 2 \epsilon - \lambda|^2} \frac{0.111 \delta + 1.31}{11.79 \delta + 1}, \end{gather}

where $\delta$ will be defined momentarily. The numerical factors are constructed such that

(2.17)\begin{gather} \lim_{\delta \to 0} \frac{0.111 \delta + 1.31}{11.79 \delta + 1} \approx a_e = 1.31, \end{gather}

(2.18)\begin{gather}\lim_{\delta \to \infty} \frac{0.111 \delta + 1.31}{11.79 \delta + 1} \approx a_i = 9.42 \times 10^{{-}3}. \end{gather}

The $\delta$-dependent ratio ensures that electron and ion collision frequencies have similar limiting behaviour for arbitrary frequency $\omega$ while still independently abiding by the limiting behaviour found when enforcing the ordering in (2.12). We then require the numerical factor $\delta$ to depend on the collisional frequency and mode frequency such that

(2.19)\begin{gather} \lim_{|\omega| / \nu_{ss'}\to 0} \delta = 0, \end{gather}

(2.20)\begin{gather}\lim_{|\omega| / \nu_{ss'} \to \infty} \delta^{{-}1} = 0. \end{gather}

There is a wide degree of freedom in constructing a specific definition of $\delta$ such that one cannot prescribe a formula from analytical analysis alone. It is necessary to make reference to direct numerical simulations using more exact operators and perform a comparison with define $\delta$. The original form of this operator was applied to codes used by Rewoldt (Rewoldt et al. Reference Rewoldt, Tang and Hastie1986, Reference Rewoldt, Tang and Hastie1987) and then modified in Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995) based on comparisons with the more exact Lorentz collision operator.

We consider two functional forms of $\delta$. The first one is taken from Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995) and is defined as

(2.21)\begin{equation} \delta_K = \left(\frac{|\omega| \epsilon}{37.2 \nu_{ei} Z_{\text{eff}}}\right)^{1/3}, \end{equation}

whereas the new one implemented in QuaLiKiz through detailed comparisons with GENE is defined as

(2.22)\begin{equation} \delta_Q = 12.0 \left(\frac{|\omega| \epsilon}{\nu_{ei} Z_{\text{eff}}}\right)^{3/2}. \end{equation}

A detailed comparison of these two operators will be made in § 4. We also note that the collision operator as currently implemented in QuaLiKiz uses a slightly different definition for the electron–ion collision frequency. The definition QuaLiKiz uses is

(2.23)\begin{equation} \nu_{ei,Q} = \frac{4}{3 \sqrt{\rm \pi}} \nu_{ei}, \end{equation}

where in other works $\nu _{ei}$ is commonly written as

(2.24)\begin{equation} \nu_{ei} = \frac{n_e e^4 \lambda_e}{4 {\rm \pi}\epsilon_0^2 m_e^2 v_{\text{th}, e}^3}, \end{equation}

which is based on energy transfer times in collision. Here, $\lambda _e$ is the Coulomb logarithm for collisions involving electrons.

3.1 Derivation of model

The intended purpose of this model is to simplify the treatment of collisions for gyrokinetic calculations in an integrated modelling framework. In reduced kinetic modelling, it is typical to bounce average the trapped electron kinetic response, much like in the bounce-averaged drift kinetic equation. The bounce averaging approach QuaLiKiz can be found in Stephens et al. (Reference Stephens, Garbet, Citrin, Bourdelle, van de Plassche and Jenko2021). To summarize, the trapped electron part of the dispersion relation requires an integration over the trapped part of velocity space; the natural variables to integrate over are the energy and the pitch angle. In the collisionless limit, the energy integration can be performed analytically using the plasma dispersion function. The goal, then, is to devise a model that includes the effects of collisions while retaining the numerical advantage of working in the collisionless limit.

We aim to construct a reduced model for DTEMs that can accurately predict and characterize the growth rate spectrum. We first describe basic features of trapped particles and TEMs. Specifically, we want to characterize the how collisionality stabilizes DTEM growth rates by considering the effect of collisional detrapping on trapped particle bounce orbits. We then distinguish between deeply trapped particles, which are difficult to detrap via collisions, and marginally trapped particles, which can easily be deflected into the passing part of velocity space. We define the effective trapped fraction by quantifying the proportion of the trapped particle population that is deeply trapped, as opposed to marginally trapped. The larger the collisionality, the lower the effective trapped fraction. Since the trapped fraction drives TEM instabilities, we hypothesize that the effective trapped fraction should drive DTEM instabilities.

The driving idea is that the deeply trapped particles should be the subgroup that drives the DTEM. Our plan of attack, then, is to construct an effective trapped fraction by first defining a parameter, which we term $a_c$, that delineates between marginally trapped and deeply trapped particles by considering the effect of collisions. The effective trapped fraction would simply be the subgroup of trapped particles that are deeply trapped. The model parameter $a_c$ at this stage is undetermined. Our goal is to determine the value of $a_c$ with the hope that it can be held constant across core tokamak parameters, namely inverse minor aspect ratio, collisionality, temperature and density gradients and magnetic geometry. For the cases considered, we are able to determine a value of $a_c$ such that the effective trapped fraction could take into account variations in inverse-aspect ratio and collisionality. Each case differs in temperature and density gradients as well as magnetic geometry. Unfortunately, we determine that the parameter $a_c$ varied substantially from case to case. In the present study, we are unable to predict the appropriate value of $a_c$ given a set of gradient and geometry information. Therefore, to improve the collision operator in QuaLiKiz, we take a different (and successful) approach as detailed in § 4.

At a given poloidal angle $\theta$, we define the trapped fraction as the proportion of particles in velocity space that are trapped. It is given by

(3.1)\begin{equation} f_t (\theta) = \frac{\displaystyle \int_{\text{trapped}} {\rm d}^{3} {v} f (\boldsymbol{v})}{\displaystyle \int {\rm d}^{3} {v} f (\boldsymbol{v})} = \frac{1}{n} \int_{\text{trapped}} {\rm d}^{3} {v} f (\boldsymbol{v}). \end{equation}

The boundary in velocity space between trapped particles and passing particles can be characterized by the pitch angle. Since a curve of constant pitch angle forms a cone in velocity space, the trapped fraction can be calculated for an isotropic velocity distribution by simply taking the ratio of a spherical cone's volume to a sphere's volume. For an isotropic distribution in the small inverse-aspect-ratio limit, the trapped particle fraction at specific poloidal angle is

(3.2)\begin{equation} f_t (\theta) = \sqrt{\epsilon (1 + \cos(\theta))}. \end{equation}

However, we are typically interested in simulations that take into account the entire flux surface, not just a single position on the flux surface. By taking advantage of the fact that the density of particles is approximately a flux function, we can compute the flux surface average of the trapped particle fraction to obtain

(3.3)\begin{equation} \langle f_t\rangle = \frac{1}{4 {\rm \pi}^2} \int {\rm d}{S} f_t (\theta) = \frac{1}{2{\rm \pi}} \int_{0}^{2 {\rm \pi}} {\rm d}{\theta} (1 + \epsilon \cos(\theta)) f_t (\theta) = \frac{2 \sqrt{2 \epsilon}}{\rm \pi} + \frac{2 \sqrt{2} \epsilon^{3/2}}{3 {\rm \pi}}. \end{equation}

The trapped fraction drives both CTEM and DTEM instabilities. The specific effect collisions have on the instability can vary between different regimes; in some cases, especially those with particularly large density gradients (Connor et al. Reference Connor, Hastie and Helander2006; Zhao, Zhang & Xiao Reference Zhao, Zhang and Xiao2017), collisions can destabilize the TEM (Kadomtsev & Pogutse Reference Kadomtsev and Pogutse1967, Reference Kadomtsev and Pogutse1970). In other cases, collisions instead stabilize the mode (Gang, Diamond & Rosenbluth Reference Gang, Diamond and Rosenbluth1991; Romanelli et al. Reference Romanelli, Regnoli and Bourdelle2007; Manas et al. Reference Manas, Camenen, Benkadda, Hornsby and Peeters2015). The literature typically refers to a simple picture where collisions can stabilize the TEM by detrapping electrons and expelling them into the passing part of velocity space. Collisional effects are in actuality more complicated given that their net effect on the growth rate depends on the collisional regime; an elementary treatment that demonstrates both destabilizing and stabilizing effects can be found in Kadomtsev & Pogutse (Reference Kadomtsev and Pogutse1970). It is also important to keep in mind that, for extremely high collisionality, one cannot adequately speak of trapped electrons due to the lack of any characteristic bounce motion (Kadomtsev & Pogutse Reference Kadomtsev and Pogutse1971). We note that earlier studies often made use of Krook collision operators with no velocity dependence for analytical simplicity, while modern supercomputers allow for the use of more sophisticated collision operators such as the Lorentz or Landau–Boltzmann operator. For the simulations conducted with typical tokamak core parameters and more realistic collision operators, the instability generally stabilizes with increasing collisionality. We introduce the dimensionless quantity $\nu ^\ast$ to characterize the collisional regime and define it as

(3.4)\begin{equation} \nu^\ast{=} \frac{4}{3 \sqrt{\rm \pi}} \frac{\nu_{ei}}{\omega_{b0} \epsilon}. \end{equation}

This is the classic collisionality parameter encountered in neoclassical theory. The factor of $\epsilon ^{-1}$ takes into account the effective collisional frequency relevant for trapped electrons in velocity space as discussed in § 2. Meanwhile, the numerical prefactor $4/(3 \sqrt {{\rm \pi} })$ is included by convention, since the characteristic collision time (Wesson Reference Wesson2012) associated with energy transfer between electrons and ions is

(3.5)\begin{equation} \tau_e = 3 (2 {\rm \pi})^{3/2} \frac{\epsilon_0^2 m_e^{1/2} T_e^{3/2}}{n_e e^2 \lambda_e} = \frac{3 \sqrt{\rm \pi}}{4 \nu_{ei}}. \end{equation}

Meanwhile, the term $\omega _{b0}$ is the characteristic electron bounce frequency (Stephens, Garbet & Jenko Reference Stephens, Garbet and Jenko2020) and defined to be

(3.6)\begin{equation} \omega_{b0} = \frac{\sqrt{T_e / m_e} \sqrt{\epsilon}}{q R_0}, \end{equation}

where $q$ is the safety factor. Small values of $\nu ^\ast$ corresponds to the scenario where trapped electrons undergo many bounce motions before undergoing a significant collision, whereas large values of $\nu ^\ast$ imply that the trapped electrons undergo many collisions before completing a single bounce motion. However, this characterization is only valid for deeply trapped thermal electrons. For any given distribution of electrons, there will be a population of low-energy trapped electrons that undergo many collisions before completing a bounce motion. We therefore consider the velocity-dependent collisionality parameter $\hat {\nu }$ defined as

(3.7)\begin{equation} \hat{\nu}(v)= \frac{Z_{\text{eff}} \nu_{ei}}{\omega_b \epsilon} \frac{1}{\hat{v}_e^3}. \end{equation}

We have included $\hat {v}_e^{-3}$ to take into account the velocity dependence of the collisional frequency. Meanwhile, $\omega _b$ corresponds to the velocity dependent bounce frequency, defined for small inverse-aspect ratio as

(3.8)\begin{equation} \omega_b = \frac{\sqrt{E/m_e} \sqrt{\epsilon}}{q R_0} \frac{\rm \pi}{2 K(\kappa)} = \omega_{b0} \hat{v}_e \frac{\rm \pi}{2 K(\kappa)}. \end{equation}

Here, $K$ is the complete elliptic integral of the first kind. We also define the trapping parameter $\kappa$ such that

(3.9)\begin{equation} \lambda = 1 - 2 \epsilon \kappa^2. \end{equation}

The trapping parameter is related to the bounce angle $\theta _b$ via

(3.10)\begin{equation} \kappa = \sin(\theta_b/2), \end{equation}

where $\theta _b$ is the poloidal angle such that $v_{\parallel } = 0$ (e.g. the banana tip). Trapped particles at a specific poloidal angle $\theta$ have the property that

(3.11)\begin{equation} |\sin(\theta/2)| \le \kappa \le 1. \end{equation}

Since all of our simulations have sufficiently small $\epsilon$, this formula is accurate for our purposes This formula is well suited for our purposes, since it has been shown to be accurate for $\epsilon \lesssim 0.3$ in accordance with our simulations (Stephens et al. Reference Stephens, Garbet and Jenko2020). We thus use the approximate formula in (3.8).

Thus, while the bounce frequency is largely dependent on the velocity, it is also dependent on the pitch angle. All else being equal, particles that are close to the trapped–passing boundary have lower bounce frequency and thus correspond to larger values of $\hat {\nu }$. Intuitively, this corresponds to marginally trapped particles being more easily detrapped. The full expression for $\hat {\nu }$ is then

(3.12)\begin{equation} \hat{\nu} = \frac{\nu_{ei}}{\omega_{b0} \epsilon} \frac{2 K(\kappa)}{{\rm \pi} \hat{v}_e^4}= \nu^{{\ast}} \frac{2 K(\kappa)}{{\rm \pi} \hat{v}_e^4}. \end{equation}

We now use the notion of the trapped particle fraction as inspiration to define the marginally trapped fraction of particles. We call particles with low values of $\hat {\nu }$ marginally trapped as defined by the condition

(3.13)\begin{equation} \hat{\nu} a_c \ge 1, \end{equation}

where $a_c$ is a constant that is ${\sim }\mathcal {O}(1)$. Equivalently, we can write this as

(3.14)\begin{equation} \hat{v}_e \le \left(a_c \nu^{{\ast}} \frac{2 K (\kappa)}{\rm \pi}\right)^{1/4}. \end{equation}

We then define the marginal trapped fraction to be the fraction of total particles that meet the above condition, leading to

(3.15)\begin{equation} f_m = \frac{\displaystyle \int_{\text{marginal}} {\rm d}^{3} {v} f_{0e}}{\displaystyle \int {\rm d}^{3} {v} f_{0e}}, \end{equation}

where $f_{0e}$ is a Maxwellian given by

(3.16)\begin{equation} f_{0e} = \frac{n_{e}}{(\sqrt{\rm \pi} v_{\text{th,e}})^3} \,{\rm e}^{{-}v^2 / v_{\text{th,e}}^2} = \frac{n_{e}}{(\sqrt{\rm \pi} v_{\text{th,e}})^3} \,{\rm e}^{-\hat{v}_e^2}. \end{equation}

Essentially, we integrate over the speed $v$ up until the marginal condition is met and then integrate over $\kappa$. To lowest order in $\epsilon$, the integral simplifies to

(3.17)\begin{equation} f_m = \int_{|\sin(\theta / 2)|}^{1} \frac{{\rm d}{\kappa} \sqrt{2 \epsilon} \kappa }{\sqrt{\kappa^2 - \sin^2(\theta / 2)}} F\left(a_c \nu^{{\ast}} \frac{2 K (\kappa)}{\rm \pi}\right), \end{equation}

where the function $F$ is defined as

(3.18)\begin{equation} F (x) = {\rm erf}(x^{1/4}) - \frac{2}{\sqrt{\rm \pi}} x^{1/4} \,{\rm e}^{{-}x^{1/2}}. \end{equation}

The constant $a_c$ therefore sets the boundary at which we consider particles to be marginally trapped. We note that

(3.19)\begin{gather} \lim_{a_c \to 0} f_m = f_t, \end{gather}

(3.20)\begin{gather}\lim_{a_c \to \infty} f_m = 0. \end{gather}

We then define the effective trapped electron fraction that contributes to the trapped electron drive to be

(3.21)\begin{equation} f_t^{{\ast}} = f_t - f_m. \end{equation}

We next take flux-surface average of $f_m$ to remove the poloidal dependence, leading to

(3.22)\begin{align} \langle f_m\rangle & \approx \frac{1}{2 {\rm \pi}} \int_{0}^{2 {\rm \pi}} {\rm d}{\theta} \int_{|\sin(\theta / 2)|}^{1} \frac{{\rm d}{\kappa} \sqrt{2 \epsilon} \kappa }{\sqrt{\kappa^2 - \sin^2(\theta / 2)}} F\left(a_c \nu^{{\ast}} \frac{2 K (\kappa)}{\rm \pi}\right), \nonumber\\ & =\frac{2\sqrt{2 \epsilon}}{\rm \pi} \int_{0}^{1} {\rm d}{\kappa} K (\kappa) \kappa F\left(a_c \nu^{{\ast}} \frac{2 K (\kappa)}{\rm \pi}\right), \end{align}

where the simplified expression is obtained by changing the order of integration. Likewise, the flux surface-averaged effective trapped electron fraction is

(3.23)\begin{equation} \langle f_t^{{\ast}}\rangle = \langle f_t\rangle - \langle f_m\rangle = \langle f_t\rangle \left(1 - \int_{0}^{1} {\rm d}{\kappa} K (\kappa) \kappa F\left(a_c \nu^{{\ast}} \frac{2 K (\kappa)}{\rm \pi}\right)\right). \end{equation}

In the regime where collisions contribute to stabilization of DTEMs, we hypothesize that the effective trapped electron fraction drives DTEM instabilities in a way that is analogous to the trapped electron fraction drive in CTEMs. For simplicity, we fix other important parameters such as the normalized temperature and density gradients and consider the growth rate $\gamma$ on a case-by-case basis while allowing $\epsilon$ and $\nu ^{\ast }$ to vary. For the purpose of this work, we vary $\nu _{ei}$ for any given case by scanning over the reference temperature of the simulation. In general, the growth rate would of course depend on both $\epsilon$ and $\nu ^\ast$ in a non-trivial way. The growth rate for CTEMs can then be written as

(3.24)\begin{equation} \gamma = \gamma(\langle f_t\rangle), \end{equation}

where gradients, magnetic geometry parameters and so on are held constant. For DTEMs, we hypothesize that the growth rate can instead be written as

(3.25)\begin{equation} \gamma = \gamma(\langle f_t^\ast\rangle). \end{equation}

That is, by holding everything else fixed, the collisionality and $\epsilon$ dependence can be summarized by the effective trapped electron fraction. Intuitively, we claim collisions stabilize the mode via the detrapping effect where electrons that are close to the trapped–passing boundary or of particularly low energy are prone to detrapping. The dimensionless quantity $a_c$ determines the exact strength of this detrapping effect. We argue that in suitable cases, the collisionality and $\epsilon$ dependence can be reduced to a one-parameter model where the parameter $a_c$ is determined by the parameters of the case (e.g. the temperature and density gradients). Equivalently, we can determine the effective inverse-aspect ratio $\epsilon ^\ast$ by calculating

(3.26)\begin{equation} \epsilon^{{\ast}} = \left(\frac{\rm \pi}{2\sqrt{2}} \langle f_t^\ast\rangle\right)^2. \end{equation}

3.2 Model testing

To test the efficacy of this model, we conduct linear gyrokinetic simulations based off of five parameter sets using GENE. Not only do we test TEM-dominated regimes, but we also test ITG-dominated regimes where there exists a subdominant TEM. In ITG-dominated regimes, an eigenvalue solver is required to analyse any subdominant TEMs; otherwise an initial value solver suffices when the TEM is dominant. We aim to estimate the model parameter $a_c$ in different cases to determine whether it varies significantly from case to case. The parameter sets are summarized in table 1. We use two basic cases: the general atomics (GA) standard case (Waltz et al. Reference Waltz, Staebler, Dorland, Hammett, Kotschenreuther and Konings1997) and a WEST case based on an electron-heated long L-mode pulse (Yang et al. Reference Yang2020). We also use three additional experimentally motivated cases based upon JET profiles where DTEMs play an important role in turbulent transport (Keilhacker et al. Reference Keilhacker, Gibson, Gormezano and Rebut2001; Tala et al. Reference Tala2019). The JET pulse numbers are 73 342 (Baiocchi et al. Reference Baiocchi, Garcia, Beurskens, Bourdelle, Crisanti, Giroud, Hobirk, Imbeaux and Nunes2015), 95 272 (Tala et al. Reference Tala2022), and 94 875 (Citrin et al. Reference Citrin, Maeyama, Angioni, Bonanomi, Bourdelle, Casson, Fable, Görler, Mantica, Mariani, Sertoli, Staebler and Watanabe2022). In the table, $R/a$ is the aspect ratio of the machine where $a$ is the minor radius of the machine while $R$ is the major radius. We also define the parameters

(3.27)\begin{gather} \tau = \frac{T_i}{T_e}, \end{gather}

(3.28)\begin{gather}\frac{1}{L_{T_{s}}} ={-} \frac{1}{T_{s}} \frac{{\rm d} T_{s}}{{\rm d} r}, \end{gather}

(3.29)\begin{gather}\frac{1}{L_{n}} ={-} \frac{1}{n} \frac{{\rm d} n}{{\rm d} r}. \end{gather}

For simplicity, we include deuterium as the ion species with charge $Z_i = 1$ and nucleon number $A_i = 2$. We include the effects of multiple ion species via $Z_{\text {eff}}$. Assuming flat $Z_\text {eff}$, quasineutrality guarantees that $L_{n_{e}} = L_{n_{i}} = L_n$. The simulations also assume a circular equilibrium defined by the safety factor $q$ and the magnetic shear $s$. For JET pulse 73 342 and JET pulse 94 875 parameters we use the $s$–$\alpha$ equilibrium model, with $\alpha$ denoting the magnetohydrodynamic parameter that characterizes the pressure gradient. In all these linear runs, we take the electrostatic limit and thus neglect electromagnetic effects. For all cases, we use the Landau–Boltzmann collision operator. Various GENE settings can be found in Appendix A.

Table 1.

Summary of cases being simulated. For each case, we list the aspect ratio, the effective charge number, gradient length scales, parameters characterizing the magnetic geometry, the ion-to-electron temperature ratio and the inverse-aspect ratio at the reference radial location. Note that we scan over the collisionality by self-consistently varying the temperature.

For each case we perform a scan over $k_{\theta } \rho _s$, where $k_{\theta }$ is the poloidal component of the wavenumber of the mode and $\rho _s = \sqrt {T_e} / \varOmega _i$, where $\varOmega _i$ is the reference ion cyclotron frequency. For each case, the interval is $0.2 \le k_{\theta } \rho _s \le 0.5$ with increment $0.1$. We also scan over $\epsilon$ for each case, and for each value of $\epsilon$ we scan over $\nu ^\ast$ by varying the reference temperature. We perform convergence tests for each individual simulation by varying the resolution of the simulation until the growth rate has sufficiently converged. For the simulations where the TEM is dominant, we use initial value simulations to determine the growth rate. For runs where the TEM is subdominant to an ITG mode, we use the eigenvalue solver to solve for the TEM growth rate. We use the sign of the real frequency as well as the behaviour of the mode with respect to collisionality to determine whether it is a TEM or not for our purposes. It is important, however, to keep in mind that instabilities are not always easily separable from each other, so some caution must be taken when labelling a mode as TEM or ITG (Merz & Jenko Reference Merz and Jenko2010). This is most easily seen in the WEST electron-heated case where the real frequency of the dominant mode continuously changes with collisionality. Here, we call a mode ITG dominated if the sign of the real mode frequency aligns with the ion drift direction and TEM dominated if the sign of the real mode frequency is associated with the electron drift direction. The other four cases are typically ITG dominated with a subdominant TEM at reference parameters. For each mode we obtain the frequency $\omega = \omega _r + i \gamma$, where $\omega _r$ is the real frequency and $\gamma$ is the TEM-dominated growth rate.

First, we plot the growth rates for all five cases as a function of $k_{\theta } \rho _{s}$ as well as $\epsilon$ in figures 1 and 2 (shown here) and figures 22–24 (shown in Appendix B). We show representative examples of the five cases in the main body of the paper while plots for the remaining cases are shown in the appendices; these additional plots display similar trends as the representative examples. Note that some of these curves do not possess the same number of points or collisionality range. This is because for all collisionality scans, we use the same values of the temperature and density for convenience. Moreover, when $\epsilon$ is small, then smaller values of $\nu ^\ast$ are required to stabilize the mode; mismatches in the domain occur when the mode is stabilized. The trapped electron fraction strongly drives the mode, which can be seen by comparing the growth rate spectrum for different values of $\epsilon$. Next, we plot the growth rate as a function of collisionality for $k_{\theta } \rho _s = 0.3$ in figures 3 and 4 (shown here) and figures 25–27 (shown in Appendix B). It can be seen that collisions stabilize the mode for moderately high values of $\nu ^\ast$ and that, with exception of a small number of cases, the growth rate decreases monotonically with collisionality. This is in congruence with previously acquired results for typical tokamak core parameters (Kotschenreuther et al. Reference Kotschenreuther, Rewoldt and Tang1995; Manas et al. Reference Manas, Camenen, Benkadda, Hornsby and Peeters2015). In general, however, the instability's dependence on collisionality is non-trivial for all of these cases. To test the model, we scan over different values of $a_c$ for any given case, where $a_c$ does not vary within a case for different values of $\epsilon$, $\nu ^\ast$, or $k_{\theta } \rho _{s}$. We then obtain a best fit such that the growth rate dependence on $\epsilon ^\ast$ is nearly independent of $\epsilon$ for the collisional cases and matches the collisionless growth rate. Figures 5 and 6 (shown here) and figures 28–30 (in Appendix B) show the result for $k_{\theta } \rho _{s} = 0.3$. The only significant deviations occur for the GA and WEST based cases at fairly high values of $\epsilon$ as well as when the growth rate saturates with increasing collisionality. An example of this saturation can be found in the blue, orange and yellow curves in figure 26. Otherwise, the deviation from the collisionless calculation is less than $5\,\%$. The determined values for $a_c$ are shown in table 2.

Figure 1.

Collisionless TEM growth rates calculated by GENE for different values of $k_{\theta } \rho _s$ plotted against $\epsilon$ using GA standard parameters. Note that the growth rate increases monotonically with $\epsilon$.

Figure 2.

Collisionless TEM growth rates calculated by GENE for different values of $k_{\theta } \rho _s$ plotted against $\epsilon$ using WEST pulse 54 178 parameters. Note that the growth rate increases monotonically with $\epsilon$.

Figure 3.

Collisional TEM growth rates calculated by GENE for different values of $\epsilon$ plotted against $\nu ^{\ast }$ using GA standard parameters where $k_{\theta } \rho _s = 0.3$. Note that the growth rate decreases nearly monotonically with $\nu ^{\ast }$.

Figure 4.

Collisional TEM growth rates calculated by GENE for different values of $\epsilon$ plotted against $\nu ^{\ast }$ using WEST pulse 54178 parameters where $k_{\theta } \rho _s = 0.3$. Note that the growth rate decreases nearly monotonically with $\nu ^{\ast }$.

Figure 5.

The TEM growth rates calculated by GENE for different values of nominal $\epsilon$ plotted against effective $\epsilon ^{\ast } = ( {\rm \pi}\langle f_t^\ast \rangle / 2 \sqrt {2})^2$ using GA Standard parameters where $k_{\theta } \rho _s = 0.3$ and $\langle f_t^\ast \rangle$ is a function of $\nu ^\ast$.

Figure 6.

The TEM growth rates calculated by GENE for different values of nominal $\epsilon$ plotted against effective $\epsilon ^{\ast } = ( {\rm \pi}\langle f_t^\ast \rangle / 2 \sqrt {2})^2$ using WEST pulse 54 178 parameters where $k_{\theta } \rho _s = 0.3$ and $\langle f_t^\ast \rangle$ is a function of $\nu ^\ast$.

Table 2.

Value of $a_c$ for each simulated case. We estimate this parameter by performing a least-squares fit between the model data (derived from collisionless simulations) and finite-collisionality simulations. Error is calculated by computing the largest pointwise difference in growth rates between the collisionless curve and any of the collisional curves.

4.1 Context

Next, we discuss improvements to the collisional operator in QuaLiKiz. Recall that the trapped electron collision operator contained a term, $\delta$, which depends on free parameters. In QuaLiKiz, it enters into the collision operator as

(4.1)\begin{equation} \nu_e = \epsilon \frac{\nu_{ei}}{\hat{v}_e^3} \frac{Z_\text{eff}}{|1 - 2 \epsilon - \lambda|^2} \frac{0.111 \delta + 1.31}{11.79 \delta + 1}. \end{equation}

The strategy is to change the definition of $\delta$ such that the TEM collisionality dependence in QuaLiKiz matches the aforementioned GENE simulations that use the Landau–Boltzmann collision operator; we use the same exact simulations done in § 3 to perform the QuaLiKiz-GENE comparisons. To appropriately change $\delta$, we parameterize it such that

(4.2)\begin{equation} \delta = a_{\delta} \left(\frac{|\omega| \epsilon}{\nu_{ei,Q} Z_{\text{eff}}}\right)^{b_{\delta}}, \end{equation}

where $a_{\delta }$ and $b_{\delta }$ are tuneable constants. Since our goal is to improve the treatment of highly collisional DTEMs in QuaLiKiz, we must ensure that the behaviour for DTEMs is preserved for low collisionality as well. From this parametrization we see that the old definition of the Krook operator formulated in Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995) uses

(4.3)\begin{gather} a_{\delta,K} = 0.30, \end{gather}

(4.4)\begin{gather}b_{\delta,K} = \tfrac{1}{3}, \end{gather}

whereas the new definition of the Krook operator uses

(4.5)\begin{gather} a_{\delta,Q} = 12.0, \end{gather}

(4.6)\begin{gather}b_{\delta, Q} = \tfrac{3}{2}. \end{gather}

Notably, the values of these parameters found in Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995) are also different from the original values detailed in DeLucia & Rewoldt (Reference DeLucia and Rewoldt1981); these free parameters were tuned to match the Lorentz operator used in Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995) to better predict the growth rate for TEMs. In arriving at the new tuning, we scan over values of $a_{\delta }$ and $b_{\delta }$ and fit the resulting collisionality dependence of the TEM growth rates to linear GENE simulations. Essentially, we use the GENE simulations as a reference to compare the old and new versions of the Krook operator used in QuaLiKiz. Because the values of $a_{\delta }$ and $b_{\delta }$ cannot be derived analytically, we expect their optimal values to depend on the specific model using the Krook operator as every kinetic model will have different underlying assumptions. Model differences constitute one potential reason for the difference in the tuning parameters. In Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995), the collision operator was used in an eigenvalue code that treated electrons drift kinetically. Moreover, that eigenvalue coded computed the eigenfunctions explicitly by solving integrodifferential equation in a Ritz basis. In contrast, in QuaLiKiz trapped electrons are treated in a bounce-averaged way with significant approximations and the Gaussian eigenfunction is prescribed in advance. Moreover, their linear initial value solver code used a Lorentzian collision operator, whereas our GENE simulations use the Landau–Boltzmann collision operator. Our results indicate that both parameters should be increased for implementation in QuaLiKiz for the examined scenarios. We observe in (4.2) that $a_\delta$ and $b_\delta$ have competing effects on the collisionality dependence; large values of $a_\delta$ decrease the effect of collisionality, whereas large values of $b_\delta$ increase the effect of collisionality. Therefore, large departures in both parameters can lead to modest changes in the computed eigenvalues.

It is important to note that we only perform linear GENE simulations to benchmark QuaLiKiz. Nonlinear validations of QuaLiKiz's turbulent transport model have been conducted in the past such as in Casati et al. (Reference Casati, Bourdelle, Garbet, Imbeaux, Candy, Clairet, Dif-Pradalier, Falchetto, Gerbaud, Grandgirard, Gürcan, Hennequin, Kinsey, Ottaviani, Sabot, Sarazin, Vermare and Waltz2009), Citrin et al. (Reference Citrin, Bourdelle, Cottier, Escande, Gürcan, Hatch, Hogeweij, Jenko and Pueschel2012) and Bourdelle et al. (Reference Bourdelle, Citrin, Baiocchi, Casati, Cottier, Garbet and Imbeaux2015). To compensate for the lack of nonlinear gyrokinetic simulations, we validate the new treatment of collisions against experimental profiles using integrated modelling, detailed in § 4.4.

4.2 The TEM growth rates

As a starting reference, we compare the calculated CTEM growth rates between QuaLiKiz and GENE. Figures for $k_{\theta } \rho _{s} = 0.3$ are shown in figures 8 and 9 (shown here) and figures 31–33 (shown in Appendix C) where the CTEM growth rates are plotted against $\epsilon$. We see here that, even in the collisionless case, there are some discrepancies between the growth rates computed in QuaLiKiz and GENE; this is to be expected since the precise growth rate is sensitive to the input parameters even when confined to similar models. Some discrepancy likely due to approximations of the bounce-averaged electrostatic potential in QuaLiKiz and the particularities of the Gaussian eigenfunction ansatz, the improvement of which is the subject of current work The deviation tends to grow in both the limit of large $\epsilon$ and the limit of small $\epsilon$. Moreover, the minimum value of $\epsilon$ that destabilizes the TEM is slightly different between GENE and QuaLiKiz simulations. For values of $\epsilon$ corresponding to mid-radius positions (e.g. $\epsilon \lesssim 0.25$), however, the agreement is satisfactory. Some discrepancy is expected given the approximations used in QuaLiKiz. In particular, QuaLiKiz assumes that the trapped particles are deeply trapped when calculating the trapped part of the dispersion relation (Stephens et al. Reference Stephens, Garbet, Citrin, Bourdelle, van de Plassche and Jenko2021). Moreover, the eigenfunctions used by QuaLiKiz make use of a Gaussian ansatz and are calculated in the fluid limit; it is known that this approximation is rather crude for TEMs (Cottier et al. Reference Cottier, Bourdelle, Camenen, Gürcan, Casson, Garbet, Hennequin and Tala2014). We take these discrepancies into account by normalizing the collisional growth rate to the CTEM growth rate when comparing collisionality scans in QuaLiKiz and in GENE; our goal is to match the collisionality dependence using the collisionless growth rate as a given. Thus we compute $\gamma _{\text {ref}}$ to be the collisionless growth rate for a given case and value of $\epsilon$ for GENE and QuaLiKiz separately.

Figure 8.

Collisionless TEM growth rates calculated by GENE and QuaLiKiz using GA standard parameters plotted against $\epsilon$, where $k_{\theta } \rho _s = 0.3$. Since this is calculated with no collisions, both old Krook and new Krook in QuaLiKiz would agree perfectly.

Figure 9.

Collisionless TEM growth rates calculated by GENE and QuaLiKiz using WEST pulse 54178 parameters plotted against $\epsilon$, where $k_{\theta } \rho _s = 0.3$. Since this is calculated with no collisions, both old Krook and new Krook in QuaLiKiz would agree perfectly.

In these figures, we also compare with an older ‘bugged’ version of QuaLiKiz. QuaLiKiz's definition of the pitch-angle parameter $\lambda$ uses the minimum value of the magnetic field as the reference magnetic field. Therefore, the trapped–passing boundary occurs at $\lambda = 1 - 2 \epsilon$, hence the appearance of $|1 - 2 \epsilon - \lambda |^2$ in the denominator of the collision operator. Essentially, the singularity at the trapped–passing boundary in velocity space is important to mimic the effects of Lorentzian collision operator. However, an older version of QuaLiKiz incorrectly used $|1 - \epsilon - \lambda |^2$ in the denominator, thus misplacing the trapped–passing boundary. This bug was fixed in the early stages of this work, and comparisons with this bugged version of QuaLiKiz are included for completeness.

Taking inspiration from Kotschenreuther et al. (Reference Kotschenreuther, Rewoldt and Tang1995), we scanned over different values for the free parameters $a_{\delta, Q}$ and $b_{\delta, Q}$, where each parameter was initially varied separately; for $b_{\delta, Q}$, we varied the exponent slightly over rational numbers between $1/4$ and $2$. Meanwhile, we varied $a_{\delta, Q}$ over two orders of magnitude. After converging on a trial set of values, we scanned over a small region in parameter space to determine the final tuning. The final values are $a_{\delta,Q} = 12.0$, $b_{\delta, Q} = 3/2$. For completeness, we illustrate the new trend for the GA standard case in figure 10 without normalizing to a reference growth rate. In addition, we plot results for $k_{\theta } \rho _{s} = 0.3$ and nominal values of $\epsilon$ in figures 11 and 12 (shown here) and in figures 34–36 (shown in Appendix C) with the normalization to reference growth rates taken into account. The nominal values of $\epsilon$ are listed in table 1. We see marked improvement for all DTEM growth rates and identify the probable culprit behind the issues in integrated modelling when comparing GENE with the Landau–Boltzmann operator, QuaLiKiz with the old Krook operator ((4.3) and (4.4)), and QuaLiKiz with the new Krook operator ((4.5) and (4.6)). Essentially, the previous version of the collision operator produced damping that was too strong.

Figure 10.

Collisional TEM growth rates calculated by GENE and QuaLiKiz. Here, we use GA standard parameters where $k_{\theta } \rho _s = 0.3$ and we plot against $\nu ^{\ast }$ for $\epsilon = 0.1667$.

Figure 11.

Collisional TEM growth rates calculated by GENE and QuaLiKiz (relative to their reference growth rates). Here, we use GA standard parameters where $k_{\theta } \rho _s = 0.3$ and we plot against $\nu ^{\ast }$ for $\epsilon = 0.1667$.

Figure 12.

Collisional TEM growth rates calculated by GENE and QuaLiKiz (relative to their reference growth rates). Here, we use WEST pulse 54 178 parameters where $k_{\theta } \rho _s = 0.3$ and we plot against $\nu ^{\ast }$ for $\epsilon = 0.10$.

4.3 Particle and Heat Fluxes

Next, we validate these improvements by computing the complete particle, ion heat and electron heat fluxes in QuaLiKiz over the range $0.1 \le k_{\theta } \rho _{s} \le 2.0$ for the three JET-based cases (due to quasineutrality, the particle flux is ambipolar). This will allow us to better interpret our integrated modelling results. Because we use JETTO-QuaLiKiz for the integrated modelling, we focus on the JET-based cases in this section. We see in figures 13–16 (shown here) and in figures 37–44 (shown in Appendix D) that the important properties of the flux structure are shifted rightwards in collisionality with the new Krook operator, indicating that the effective impact of collisions in the new operator is reduced. We can see in figures 13 and 37 that the outward particle flux at the nominal value of $\nu ^{\ast } = 0.5$ for these simulations is reduced; an erroneously large outward particle flux would lead to a less peaked density profile in the core. We note that, in other work, the net impact of the new QuaLiKiz Krook operator on low-collisionality JET hybrid H-mode density peaking was indeed shown to be low (Citrin et al. Reference Citrin, Maeyama, Angioni, Bonanomi, Bourdelle, Casson, Fable, Görler, Mantica, Mariani, Sertoli, Staebler and Watanabe2022); it is evident from figure 36 that the collisionality dependence of the old Krook operator was already quite close to that of GENE. Indeed, it is important that we do not worsen the collision operator in the cases where collisionality dependence is correct. Since QuaLiKiz's collision operator functioned quite well in the low-collisionality limit, a crude approach such as artificially reducing $\nu ^{\ast }$ by a factor of $10$ would be inappropriate.

Figure 13.

Total integrated particle flux calculated in QuaLiKiz for JET pulse 73 342 parameters plotted against $\nu ^{\ast }$ where $\epsilon = 0.24$. Note that nominal $\nu ^{\ast } = 0.50$.

Figure 14.

Total integrated ion heat flux calculated in QuaLiKiz for JET pulse 73 342 parameters plotted against $\nu ^{\ast }$ where $\epsilon = 0.24$. Note that nominal $\nu ^{\ast } = 0.50$.

Figure 15.

Total integrated electron heat flux calculated in QuaLiKiz for JET pulse 73 342 parameters plotted against $\nu ^{\ast }$ where $\epsilon = 0.24$. Note that nominal $\nu ^{\ast } = 0.50$.

Figure 16.

Ratio between total integrated particle flux and total integrated electron heat flux calculated in QuaLiKiz for JET pulse 73 342 parameters plotted against $\nu ^{\ast }$ where $\epsilon = 0.24$. Note that nominal $\nu ^{\ast } = 0.50$.

4.4 Integrated Modelling

Finally, we implement the new collision operator by running JETTO-QuaLiKiz simulations for H-mode and L-mode $\nu ^{\ast }$ collisionality scans reported in Tala et al. (Reference Tala2019) as shown in figures 17–19 (shown here) and figures 45–47 (shown in Appendix E), as well as two high-collisionality H-modes (figures 48 and 49 in Appendix E), a very high-collisionality ohmic L-mode (figure 50 in Appendix E) and a mid-collisionality heated L-mode (figure 20). In these simulations, we include an neutral beam injection-driven particle source. Moreover, the profiles for $\rho \ge 0.8$ are treated as boundary conditions, and only the radial domain $\rho < 0.8$ is simulated in the integrated model. Note that finite rotation effects are not considered in this work. We see that the correct behaviour is preserved in low-collisionality cases. We also observe improvement in the density profile predictions for mid- to high-collisionality cases, L-mode cases without much change in the temperature profile prediction. For deeper insight as to why the new Krook operator improves the density peaking, we compute the instability spectrum for a high-collisionality L-mode case and plot the results in figure 21. We see here that for $k_{\theta } \rho _{s} > 0.2$ the mode switches to a proper TEM-dominated regime for the new Krook operator, while the ITG mode is dominant for low $k_{\theta } \rho _{s}$. In the old Krook operator, the TEM never becomes dominant for this example. It is known that the interplay between the two modes in the mixed ITG–TEM regime is important for density peaking (Fable et al. Reference Fable, Angioni and Sauter2010). We strongly suspect that the improper treatment of DTEMs in QuaLiKiz is responsible for a large part in the incorrect density profiles.

Figure 17.

Integrated modelling of a JET H-mode case as part of collisionality scan in Tala et al. (Reference Tala2019), with lower collisionality.

Figure 18.

Integrated modelling of a JET H-mode case as part of collisionality scan in Tala et al. (Reference Tala2019), with higher collisionality.

Figure 19.

Integrated modelling of a JET H-mode case as part of collisionality scan in Tala et al. (Reference Tala2019).

Figure 20.

Integrated modelling of a medium-collisionality heated L-mode (Weisen et al. Reference Weisen, Maggi, Oberparleiter, Casson, Camenen, Menmuir, Horvath, Auriemma, Bache and Bonanomi2020).The version of QuaLiKiz with the incorrect implementation of the collision operator is in red.

Figure 21.

Instability spectrum plotted against $k_{\theta } \rho _s$ for JET L-mode at $\rho _{\text {norm}} = 0.75$. Note that with new Krook we retain the correct TEM instability, whereas for old Krook the TEM is stable.

Note that we also show particularly poor cases which exhibit noticeable hollowing in the density profile between the core and the edge region; these cases are associated with the highest-collisionality regimes among all the discharges studied. We note that there is some improvement in the density profile behaviour, although the incorrect density hollowing is still present.

As mentioned previously, in older versions of QuaLiKiz, the Krook operator was implemented incorrectly. Essentially, the pitch-angle dependence of the old Krook operator was coded in QuaLiKiz such that the singularity in the Krook operator was not located at the trapped–passing boundary in velocity space, but some other region of velocity space. For completeness, we also include a comparison between this bugged version of the old Krook operator, the bug-fixed version of the old Krook operator and the new Krook operator in figure 20. In this case, the red curve indicates the version of QuaLiKiz that used the bug-fixed version of the old Krook operator. The bug fix alone resulted in marked improvement in this particular case. We also plot the growth rates predicted by the old bugged Krook operator in figures 11 and 12 (shown here) and in figures 34–36 (shown in Appendix C).

Lastly, we note that the purpose of this study was to improve the density peaking predictions in the core. Unfortunately, one drawback of the new changes is that the temperature profile prediction has degraded in some cases. This is most likely due to large number of approximations made in QuaLiKiz; for instance, even in the collisionless regime, QuaLiKiz's predicted TEM growth rates differ from GENE. Further improvements may require a more exact gyrokinetic code that better matches linear GENE simulations.