Quasilinear modelling of collisional trapped electron modes

Abstract

The simplification of collision operators is necessary for quasilinear turbulence modelling used with integrated modelling frameworks, such as the gyrokinetic code QuaLiKiz. The treatment of collisions greatly impacts the accuracy of trapped electron mode (TEM) modelling, which is necessary to predict the electron heat flux and the balance between inward and outward particle fluxes. In particular, accurate particle flux predictions are necessary to successfully model density peaking in the tokamak core. We explored two ways of improving collisional TEM model reduction for tokamak core plasmas. First, we carried out linear GENE simulations to study the complex interplay between collisions and trapped electrons. We then used these simulations to define an effective trapped fraction to characterize the collisional TEM based on two key parameters, the local inverse-aspect ratio $\epsilon$ and the collisionality $\nu ^\ast$. One aspect missing from analytical TEM research is that the collisional frequency and the bounce-transit frequency are both velocity dependent; this effective trapped fraction takes both into account. In doing so, we determined that two parameters are not enough to model the collisional TEM, as an additional third free parameter was necessary. We determined that this model, as currently formulated, is not suitable for integrated modelling purposes. Second, we directly improved QuaLiKiz's Krook operator, which relies on two free parameters. We determined that these parameters required adjustments against higher-fidelity collisional models. In order to improve density profile predictions when paired with integrated models, we refined the Krook operator by using GENE simulations as a higher-fidelity point of comparison. We then demonstrate strong improvement of density peaking predictions of QuaLiKiz within the integrated modelling framework JETTO.

1 Introduction

In tokamaks, electrostatic drift instabilities constitute an important component of turbulence. Specifically, microinstabilities such as trapped electron modes (TEMs) and ion temperature gradient (ITG) driven modes are responsible for driving turbulence, which dominates tokamak transport, and thus are key to modelling confinement physics (Liewer et al. 1986; Kadomtsev & Pogutse 1995). Many effects such as $E\times B$ shearing ($E$ and $B$ being the electromagnetic fields), collisions and high magnetic shear play a role in setting the critical temperature gradient thresholds of the transport-driving instabilities. In particular, collision detrapping of electrons can stabilize TEM turbulence. Experimentally, it has been shown that turbulent convection is responsible for the peaking of the density profile in the tokamak core (Hoang et al. 2003). The density peaking can also be impacted by the balance between ITG and TEM instabilities, since the particle flux's specific dependence on plasma parameters can depend strongly on the turbulent regime (Hoang et al. 2003; Angioni et al. 2009; Fable, Angioni & Sauter 2010). Consequently, in the context of integrated modelling, an inaccurate treatment of collisionality can lead to severe discrepancies in density profiles in the core. While collisions can in principle be accounted for through the use of sufficiently complex collision operators when simulating microinstabilities, the use of such operators incurs a large computational cost, restricting their usefulness when considering their inclusion in reduced transport models. Simplified models of collisions are therefore necessary when computational speed is a priority.

Motivated by the aforementioned potential discrepancies in integrated modelling, our study of TEMs is the primary focus of this work. The TEM instabilities are driven by the presence of a trapped electron population in the tokamak whereby the electrons perform a bounce motion along the magnetic field lines (Kadomtsev & Pogutse 1970). Collisions serve to detrap particles by scattering them into the passing part of velocity space; as a result, TEMs are typically divided into collisionless TEMs (CTEMs) and dissipative TEMs (DTEMs) (Kadomtsev & Pogutse 1971). In high-collisionality regimes (where the collisionality parameter $\nu ^\ast$, defined in (3.4) as the ratio of the detrapping collision frequency to the bounce frequency, is such that $\nu ^\ast \not \ll 1$), the particle undergoes many collisions before completing a full bounce motion, thus leading to detrapping of the particle and thereby producing a stabilizing effect. However, it is also important to recognize that passing electrons play a role in TEM instabilities, thus creating a complicating interaction between the trapped and passing electron populations that leads to collisions becoming a destabilizing effect in certain low-collisionality regimes (Connor, Hastie & Helander 2006). Moreover, the trapped–passing boundary in velocity space constitutes a region where collisional effects are amplified.

Dimensionless quantities that attempt to characterize the collisional regime are typically constructed by calculating these frequencies for thermal particles. While these quantities are useful to describe the overall effect of collisions, it is important to realize that collisions do not impact trapped particles in a uniform manner. Thus, we expect reasonably realistic collision operators to incorporate both energy and pitch-angle dependence in order to treat collisions properly for TEMs. An overview of collision operators used in this work is given in § 2.

There are two primary goals for this work. First, we attempt to construct a model to characterize DTEM growth rates in the regime where collisions are stabilizing. The intended goal of this model was to simplify the treatment of trapped electron collisions in an integrated modelling context. We do so with the aid of the gyrokinetic electromagnetic numerical experiment (GENE) code (Jenko et al. 2000) by simulating a number of linear TEMs. The GENE code solves the gyrokinetic equation in field-line coordinates by using the well-known $\delta f$ approximation, where $\delta f$ refers to the perturbed distribution function. We use the $s$–$\alpha$ model to conduct linear simulations using a linearized Landau–Boltzmann collision operator, where s and $\alpha$ are the magnetic shear and Shafranov parameter respectively. The core idea is to construct an effective trapped electron fraction that takes into account the fact that not all trapped electrons contribute to TEM destabilization at high collisionality. By scanning over the conventional trapped electron fraction as well as the collisionality for numerous base cases, we show that collisions and the conventional trapped electron fraction do not independently determine the growth rate of a DTEM. Instead, the effective trapped fraction, a function of both the conventional trapped electron fraction and the collisionality, suitably determines the growth rate of the instability (along with other independent parameters such as the density gradient, temperature gradient and so on). We determined that this model requires a free parameter, which we term $a_c$, to accurately model the effect of collisions; this parameter depends on other plasma parameters such as gradients and geometry. Because this free parameter varies over an order of magnitude for the examined cases with no clear pattern, we determined that this model is not suitable for integrated modelling. Nevertheless, we present the results of this study, as it still contained valuable physics insights and could present opportunities for improvement in the future.

The second goal is address the failure of the quasilinear gyrokinetic code QuaLiKiz (Bourdelle et al. 2007; Bourdelle 2015; Citrin et al. 2017) to predict density peaking in the tokamak core in high-collisionality regimes when coupled with integrated modelling suites. Originally based on the eigenvalue code Kinezero (Bourdelle et al. 2002), QuaLiKiz calculates the anomalous quasilinear transport arising from ITG, TEM and electron temperature gradient microinstabilities. Because calculation speed is key for integrated modelling, QuaLiKiz is designed to be orders of magnitude faster than first-principles-based linear gyrokinetic codes such as GENE via numerous approximations. The kinetic dispersion relation is separated species by species and also separated between trapped particles and passing particles. Approximations such as using the $s$–$\alpha$ equilibrium, considering only electrostatic fluctuations and assuming Gaussian eigenfunctions allows the code to perform full computations in ${\sim }1$ CPU per wavenumber. The full derivation can be found in Stephens et al. (2021). The simplified collision operator in QuaLiKiz is Krook-like in that it contains a simple velocity dependence and no differentiation of the perturbed distribution function. In addition, the Krook-like operator retains pitch-angle dependence to account for particle detrapping near the trapped–passing boundary; the operator also depends on two free parameters. In investigations of high collisionality in Joint European Torus (JET) and the tungsten steady-state tokamak (WEST), it has been found that QuaLiKiz quasilinear flux values lead to inaccurate density profile flattening predictions when used with integrated modelling suites (Tala et al. 2019). Predicting the correct density peaking is essential for performance in terms of pulse length and energy content, as the density will impact the current drive and the fusion power scales quadratically with the fuel density. We also note that electromagnetic effects are also important to include for reactor-like plasmas (Hein et al. 2010; Fable et al. 2019), although these effects are not included in QuaLiKiz. Given previous collisionality studies as the numerous approximations applied to TEMs in QuaLiKiz, we have high confidence that the improper treatment of DTEMs is the primary culprit for the mismatch. To address the inaccuracy in density profile predictions, we seek to improve the Krook-like collision operator in QuaLiKiz to properly treat DTEMs. To tune QuaLiKiz's collision operator, we use GENE simulations as a point of comparison. We mandate that DTEMs simulated by QuaLiKiz exhibit a collisionality dependence in the DTEM growth rates matching that of GENE. To verify the improvement, we use the newly improved version of QuaLiKiz to simulate heat and particle transport of JET H-mode and L-mode collisionality scans within the core transport solver JETTO as part of the JET integrated transport code (JINTRAC) suite (Cenacchi & Taroni 1998; Romanelli et al. 2014; Citrin et al. 2017), many of which were investigated in Tala et al. (2019). Thanks to the improved Krook-like operator obtained from comparisons with GENE simulations, we attain proper predictions of density peaking in the core when using QuaLiKiz in predicting JET profiles.

This paper is organized as follows. First, we briefly describe the gyrokinetic modelling equations and the collision operators used in this work in § 2. We then attempt to construct an effective trapped fraction model in § 3. Next, in § 4 we discuss improvements made to QuaLiKiz and analyse the ramifications of these improvements via the integrated modelling framework. Finally, we present conclusions in § 5. The GENE simulations settings can be found in Appendix A. Due to the large number of simulations produced for this work, we present figures for isolated GENE simulations in Appendix B, figures involving QuaLiKiz growth rate predictions in Appendix C, figures involving QuaLiKiz flux predictions in Appendix D and figures produced from JETTO simulations in Appendix E.

2 Gyrokinetic modelling

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. (2000).

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. (2019), 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. (2014), 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. (2021).

2.2 Collision operators

In this work, we make use of two collision operators: the Landau–Boltzmann collision operator (Hazeltine & Waelbroeck 2004) in GENE and a Krook-like operator (Kotschenreuther, Rewoldt & Tang 1995) 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 2004). 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 (2019).

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 1954; DeLucia & Rewoldt 1981) 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 1986). 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 2007). 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 (1981) and Rewoldt, Tang & Hastie (1987). 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$.¹ 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 1972). 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 1981)

(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. 1986, 1987) and then modified in Kotschenreuther et al. (1995) 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. (1995) 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 Effective trapped electron fraction

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. (2021). 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. 2006; Zhao, Zhang & Xiao 2017), collisions can destabilize the TEM (Kadomtsev & Pogutse 1967, 1970). In other cases, collisions instead stabilize the mode (Gang, Diamond & Rosenbluth 1991; Romanelli et al. 2007; Manas et al. 2015). 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 (1970). 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 1971). 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 2012) 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 2020) 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. 2020). 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. 1997) and a WEST case based on an electron-heated long L-mode pulse (Yang et al. 2020). We also use three additional experimentally motivated cases based upon JET profiles where DTEMs play an important role in turbulent transport (Keilhacker et al. 2001; Tala et al. 2019). The JET pulse numbers are 73 342 (Baiocchi et al. 2015), 95 272 (Tala et al. 2022), and 94 875 (Citrin et al. 2022). 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 2010). 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. 1995; Manas et al. 2015). 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.

3.3 Discussion

Our simple parameter model successfully describes the collisionality dependence of the growth rate. This may come as a surprise, given that the Landau–Boltzmann operator contains physics such as pitch-angle scattering and energy diffusion. For TEMs, however, the dominant effect of collisions is due to the velocity dependence of the collision frequency and the scattering of trapped particles into the passing part of velocity space. Both effects are included in conventional trapped particle Krook operators, where the effective velocity-dependent collision frequency is obtained by simply dividing by the inverse-aspect ratio. These dominant collisional effects are also encoded in the effective trapped fraction model. Although we can successfully characterize the TEM growth rates with this approach, there are a few caveats. While $a_c\sim \mathcal {O}(1)\unicode{x2013}\mathcal {O}(10)$ for all cases, its value varies from 1.2 to 8.1 as seen in table 2. The wide variance indicates that we have no guarantee that $a_c$ will be of $\mathcal {O}(1)$ in all cases; parameter regimes with even larger values of $a_c$ may exist. We expect fundamental parameters such as the temperature gradients, density gradient, and the magnetic geometry to strongly affect the value of $a_c$. For instance, the best fit for the WEST-based case leads to $a_c = 1.2$, substantially smaller than the other cases analysed; distinguishing factors of the WEST-based case include a large electron temperature gradient and a non-unity ion-to-electron temperature ratio. In particular, (3.14) indicates that small values of $a_c$ are associated with weaker collisional damping.

Analysis of the WEST-based case also reveals a shortcoming of the model. In particular, the model does not predict the correct growth rate trends at very high collisionality for $\epsilon < 0.10$. The underlying cause is made clear in figure 4, where the growth rates for $\epsilon = 0.06$ and $\epsilon = 0.08$ flatten with increasing collisionality. These modes have entered an unorthodox resistive regime at high collisionality, likely because the drive from the trapped electron fraction is quite low. This can be further confirmed in figure 7, which plots the real frequency of these modes over collisionality. At high collisionality, the real frequencies for $\epsilon \le 0.10$ suddenly become much smaller in magnitude, indicative of a different mode branch than the conventional DTEM. As such, we do not expect our model to succeed in this regime.

Figure 7. Collisional TEM real frequencies calculated by GENE for different values of $\epsilon$ plotted against $\nu ^{\ast }$ using WEST pulse 54 178 parameters where $k_{\theta } \rho _s = 0.3$. Note that at high collisionality, some of the real frequencies trend far closer to $0$, indicating an unorthodox TEM. Convergence tests were performed to confirm the accuracy of these results.

Because the two parameters $\epsilon$ and $\nu ^{\ast }$ did not fully characterize the DTEM growth rates, the derived model could not be easily extended into a quasilinear code. Unfortunately, the fitted value of $a_c$ varies too strongly in terms of plasma parameters. A brief implementation of this model was attempted in QuaLiKiz to fix the aforementioned issue with core density peaking for JET pulses 73 342 and 95272 with various values of $a_c$. We omit this specific analysis for brevity, as the implementation attempt did not produce the requisite density peaking.

4 Improvements to QuaLiKiz's collisional operator

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. (1995) 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. (1995) are also different from the original values detailed in DeLucia & Rewoldt (1981); these free parameters were tuned to match the Lorentz operator used in Kotschenreuther et al. (1995) 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. (1995), 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. (2009), Citrin et al. (2012) and Bourdelle et al. (2015). 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. 2021). 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. 2014). 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. (1995), 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. 2022); 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. (2019) 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. 2010). 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. (2019), with lower collisionality.

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

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

Figure 20. Integrated modelling of a medium-collisionality heated L-mode (Weisen et al. 2020).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.

5 Conclusions

The complex interplay between the trapped electron dynamics and collisions pose a difficult challenge in TEM modelling, especially because a large degree of model reduction is mandatory for profile predictions in integrated modelling. Nonetheless, we attempted to construct a reduced collisionality model with the aid of GENE simulations using the notion of an effective trapped electron fraction. Thanks to this investigation, we were able to gain insight as to how collisions impact TEM growth rates by conducting scans over $\epsilon$ and $\nu ^{\ast }$. The resulting model was successful in characterizing the growth rate spectrum for DTEMs, providing a straightforward parametric relationship between CTEM growth rates and DTEM growth rates. However, the above two parameters were not enough to completely characterize the growth rate; since a free parameter remained, we could not converge upon a universally applicable collisionality model.

Since the effective trapped fraction model could not be used for reduced models such as QuaLiKiz, we therefore focused on improving the Krook-like operator directly. In this work, we successfully improved the collisional model in QuaLiKiz by tuning the dimensionless parameters of the Krook-like operator to GENE simulations using the Landau–Boltzmann operator. The result was an improvement in moderately collisional L-mode cases while also preserving the correct treatment for low collisionality. The improved collision operator allowed for improved particle flux predictions in experimentally relevant regimes by lessening outward particle transport, thus leading to more peaked density profiles as seen in experiments at both low- and high-collisionality regimes. We note further improvements can be made to QuaLiKiz given that $\nu ^{\ast } > 1$ cases still exhibit noticeable and unphysical density profile hollowing in between the core and the edge. However, we are likely hitting the limitations with regards to improvements of reduced gyrokinetic models such as QuaLiKiz. We have already seen in this work that changes made to improve the density peaking predictions can also lead to worse temperature profile predictions. More gains could potentially be made via techniques such as neural network surrogate modelling of higher-fidelity gyrokinetic codes as well as GPU acceleration of hybrid codes.