Gyrokinetic simulations of neoclassical electron transport and bootstrap current generation in tokamak plasmas in the TRIMEG code

TL;DR

This paper details the implementation of neoclassical physics, including a simplified collision operator, in the TRIMEG gyrokinetic code to simulate electron transport and bootstrap current in tokamak plasmas, enabling more comprehensive fusion plasma modeling.

Contribution

It introduces a novel implementation of neoclassical physics in the unstructured mesh gyrokinetic code TRIMEG, including flux surface averaging without poloidal coordinates.

Findings

Good agreement with neoclassical theory at large aspect ratio

Discrepancies observed at moderate aspect ratio and realistic geometry

Demonstrates capability for self-consistent neoclassical effects in gyrokinetic simulations

Abstract

For magnetic confinement fusion in tokamak plasmas, some of the limitations to the particle and energy confinement times are caused by turbulence and collisions between particles in toroidal geometry, which determine the "anomalous" and the neoclassical transport, respectively. In this work, we focus on the implementation of neoclassical physics in the gyrokinetic code TRIMEG, which is a TRIangular MEsh-based Gyrokinetic code that can handle both the closed and open field line geometries of a divertor tokamak. We report on the implementation of a simplified Lorentz collision operator in TRIMEG. Since the code uses an unstructured mesh, a procedure for calculating the flux surface averages of particle and energy fluxes and the bootstrap current is derived without relying on the poloidal coordinate, which is useful also for other simulations in unstructured meshes. With the newly implemented collision operator, we study electron transport and bootstrap current generation for various simplified and realistic geometries. In comparison to neoclassical theory, good agreement is obtained for the large aspect ratio case regarding the particle and energy fluxes as well as the bootstrap current. However, some discrepancies are observed at moderate aspect ratio and for a case with the realistic geometry of the ASDEX Upgrade tokamak. These deviations can be explained by different treatments and approximations in theory and simulation. In this paper, we demonstrate the capability to calculate the electron transport and bootstrap current generation in TRIMEG, which will allow for the self-consistent inclusion of neoclassical effects in gyrokinetic simulations in the future.

Tables

Table 1. Table 1: Normalization used in TRIMEG.

TRIMEG Normalization units
$m_N$	$m_i$
$R_N$	$1$ m
$T_N$	$T_{i,ref}$
$v_N$	$v_{th,i}=\sqrt{\frac{2T_N}{m_N}}$
$t_N$	$\frac{R_N}{v_N}$
$B_N$	$1$ Tesla

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\delta f(t)=-\frac{\mathrm{d}}{\mathrm{d}t}{f}_0+C(f,f),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}=\frac{\partial }{\partial t}+\dot{\mathbf{R}}\cdot \nabla +{\dot{v}}_{\Vert }\frac{\partial }{\partial v_{\Vert }},$$

$$f_{\mathrm{M}}=\frac{n_0}{(2T/m{)}^3{\pi }^{3/2}}\exp \left\{-\frac{m{v}_{\Vert }^2}{2T}-\frac{m\mu B}{T}\right\},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\ln f_{\mathrm{M}}$$

$$f_0(z,t)=\frac{N_{ph}}{N}P(z,t)g(z,t)\approx \frac{N_{ph}}{N}\sum _{i=1}^N{p}_i\frac{\delta (z-z_i(t))}{J(z)},$$

$$\delta f(z,t)=\frac{N_{ph}}{N}W(z,t)g(z,t)\approx \frac{N_{ph}}{N}\sum _{i=1}^N{w}_i\frac{\delta (z-z_i(t))}{J(z)},$$

$$\frac{d}{dt}{w}_i(t)=-p_i\frac{d}{dt}\ln f_0(z_i(t)),$$

$$\frac{d}{dt}{p}_i(t)=p_i\frac{d}{dt}\ln f_0(z_i(t)),$$

$$\frac{d}{dt}(w_i(t)+p_i(t))=0.$$

$$\frac{dw}{dt}=-(1-w){\vec{v}}_d\cdot \nabla \ln f_0.$$

$$C_{ei}(\delta f_e)={\nu }_0\frac{1}{2}\frac{\partial }{\partial \xi }(1-{\xi }^2)\frac{\partial }{\partial \xi }.$$

$$\xi=\xi_{0}(1-\nu\Delta t)+({\color[rgb]{0,0,0}\mathcal{R}}-0.5)[12(1-\xi_{0}^{2})\nu\Delta t]^{1/2},$$

$$\frac{dw}{dt}=(1-w)({\vec{v}}_d\cdot \vec{\kappa })=(1-w)(v_{d,R}{\kappa }_R+v_{d,Z}{\kappa }_Z),$$

$$\begin{array}{rl}\vec{\kappa }& =-\frac{\partial \ln f_0}{\partial \vec{x}}=-\left(\frac{\partial \ln f_0}{\partial n}\frac{\partial n}{\partial R}+\frac{\partial \ln f_0}{\partial T}\frac{\partial T}{\partial R}\right)\hat{R}-\left(\frac{\partial \ln f_0}{\partial n}\frac{\partial n}{\partial Z}+\frac{\partial \ln f_0}{\partial T}\frac{\partial T}{\partial Z}\right)\hat{Z}\\ & =-\left[\frac{\partial \ln n}{\partial R}+\left(\frac{m{v}^2}{2T}-\frac{3}{2}\right)\frac{\partial \ln T}{\partial R}\right]\hat{R}-\left[\frac{\partial \ln n}{\partial Z}+\left(\frac{m{v}^2}{2T}-\frac{3}{2}\right)\frac{\partial \ln T}{\partial Z}\right]\hat{Z}.\end{array}$$

$$\begin{split}\frac{dw}{dt}=(1-w)\left\{-\bar{v}_{d,R}\left[\frac{\partial\ln{\color[rgb]{0,0,0}\bar{n}}}{\partial{\color[rgb]{0,0,0}\bar{R}}}+\left(\frac{\bar{m}\bar{v}_{\parallel}^{2}}{\bar{T}}+\frac{2\bar{\mu}\bar{m}\bar{B}}{\bar{T}}-\frac{3}{2}\right)\frac{\partial\ln{\color[rgb]{0,0,0}\bar{T}}}{\partial{\color[rgb]{0,0,0}\bar{R}}}\right]\right.\\ \left.-\bar{v}_{d,Z}\left[\frac{\partial\ln{\color[rgb]{0,0,0}\bar{n}}}{\partial{\color[rgb]{0,0,0}\bar{Z}}}+\left(\frac{\bar{m}\bar{v}_{\parallel}^{2}}{\bar{T}}+\frac{2\bar{\mu}\bar{m}\bar{B}}{\bar{T}}-\frac{3}{2}\right)\frac{\partial\ln{\color[rgb]{0,0,0}\bar{T}}}{\partial{\color[rgb]{0,0,0}\bar{Z}}}\right]\right\}.\end{split}$$

$$\xi=\xi_{0}(1-\bar{\nu}\Delta{\color[rgb]{0,0,0}\bar{t}})+({\color[rgb]{0,0,0}\mathcal{R}}-0.5)[12(1-\xi_{0}^{2})\bar{\nu}\Delta{\color[rgb]{0,0,0}\bar{t}}]^{1/2},$$

$$⟨F⟩=\frac{{\int }_{\Delta V}{d}^3\vec{R}F}{{\int }_{\Delta V}{d}^3\vec{R}}=\frac{d\psi }{dV}\int \frac{dS}{|\nabla \psi |}F,$$

$$⟨F⟩=\frac{1}{\Delta V}{\int }_{\Delta V}{d}^3{\vec{R}}^{\prime }\frac{1}{\Delta \hat{V}}{\int }_{\Delta \hat{V}}{d}^3\vec{R}\frac{1}{J_{\vec{R}}}\frac{N_{ph}}{N_m}\sum _{p=1}^{N_m}{w}_p\delta (\vec{R}-{\vec{R}}_p)\delta (\vec{v}-{\vec{v}}_p)X.$$

$$⟨F⟩=\frac{1}{\Delta V}\frac{N_{ph}}{N_m}\sum _{p\in \Delta V}{w}_p{X}_p.$$

$$F(\psi ,\theta ,\varphi )=\frac{{\int }_{\Delta \tilde{V}}{d}^3\vec{R}\delta fX}{{\int }_{\Delta \tilde{V}}{d}^3\vec{R}},$$

$$⟨F{⟩}_{\varphi }=\frac{{\int }_0^{2\pi }d\varphi F}{{\int }_0^{2\pi }d\varphi }.$$

$$⟨F{⟩}_{\varphi }=\frac{1}{{\int }_0^{2\pi }d\varphi }\frac{1}{\Delta \tilde{V}}\frac{N_{ph}}{N_m}{\int }_{\Delta \tilde{V}}{d}^3\vec{R}\int d\varphi \sum _{p=1}^{N_m}{w}_p\delta (\vec{R}-{\vec{R}}_p)\delta (\vec{v}-{\vec{v}}_p)X,$$

$$⟨F{⟩}_{\varphi }=\frac{1}{2\pi \Delta \tilde{S}}\Delta \varphi \frac{N_{ph}}{N_m}\sum _{p\in \Delta \tilde{S}}\frac{1}{R_p}{w}_p{X}_p,$$

$$\frac{⟨F{⟩}_{\varphi }}{⟨n{⟩}_V}=\frac{1}{2\pi \Delta \tilde{S}}\frac{V_{tot}}{N_m}\sum _{p\in \Delta \tilde{S}}\frac{1}{R_p}{w}_p{X}_p,$$

$$A_{1e}^{\prime }$$

$${\Gamma }_e$$

$$q_e$$

$$⟨(J_{\parallel }-J_{\parallel s})/h⟩$$

$$K_{11}=1.04\left(\frac{1}{1+2.01{\nu }_{\ast e}^{1/2}+1.53{\nu }_{\ast e}}+\frac{{ϵ}^{3/2}(0.89^2/1.53){\nu }_{\ast e}{ϵ}^{3/2}}{1+0.89{\nu }_{\ast e}{ϵ}^{3/2}}\right),$$

$$K_{12}=1.20\left(\frac{1}{1+0.76{\nu }_{\ast e}^{1/2}+0.67{\nu }_{\ast e}}+\frac{{ϵ}^{3/2}(0.56^2/0.67){\nu }_{\ast e}{ϵ}^{3/2}}{1+0.56{\nu }_{\ast e}{ϵ}^{3/2}}\right),$$

$$K_{22}=2.55\left(\frac{1}{1+0.45{\nu }_{\ast e}^{1/2}+0.43{\nu }_{\ast e}}+\frac{{ϵ}^{3/2}(0.43^2/0.43){\nu }_{\ast e}{ϵ}^{3/2}}{1+0.43{\nu }_{\ast e}{ϵ}^{3/2}}\right),$$

$$K_{13}=2.3[1+1.02{\nu }_{\ast e}^{1/2}+1.07{\nu }_{\ast e}{]}^{-1}[1+1.07{\nu }_{\ast e}{ϵ}^{3/2}{]}^{-1},$$

$$K_{23}=4.19[1+0.57{\nu }_{\ast e}^{1/2}+0.61{\nu }_{\ast e}{]}^{-1}[1+0.61{\nu }_{\ast e}{ϵ}^{3/2}{]}^{-1},$$

$$K_{33}=1.83[1+0.68{\nu }_{\ast e}^{1/2}+0.32{\nu }_{\ast e}{]}^{-1}[1+0.66{\nu }_{\ast e}{ϵ}^{3/2}{]}^{-1}.$$

$${\bar{A}}_{1e}^{\prime }=R_N{A}_{1e}^{\prime }=\frac{\partial A_{1e}}{\partial \bar{r}},$$

Full text

I Introduction

Neoclassical transport determines the minimum level of transport for plasma confinement in tokamaks or stellarators and is one of the key issues in the optimization of experimental devices Beidler et al. (2021). The bifurcation of neoclassical transport for steep equilibrium profiles has also been proposed as a possible mechanism for the transition to the high confinement regime in tokamak plasmas Helander (1998). The neoclassical theory has been developed for the comparison and the prediction of the confinement properties in toroidally confined plasmas Hinton and Hazeltine (1976); Hirshman and Sigmar (1981). Standard neoclassical theory, while illustrating the dominant physics mechanism, assumes a static magnetic equilibrium and a steady-state solution and solves the gyrokinetic/drift kinetic equation asymptotically using the poloidal Larmor radius as a small parameter. It is demonstrated from the gyrokinetic simulation that the large orbit effects can lead to significant discrepancy compared with traditional neoclassical theory Lin, Tang, and Lee (1997). It is also found that the nonlocal effect of neoclassical transport is important in matching experimental results Wang et al. (2006a). Recently, neoclassical transport in the plasma edge with open field lines has attracted significant attention Chang, Ku, and Weitzner (2004). In edge transport studies, a fully nonlinear collision operator has been developed to handle the non-Maxwellian electron distribution function in the plasma edge Zhao, Chankin, and Coster (2019). In all these aspects, the traditional neoclassical theory is not applicable. Furthermore, the self-consistent interaction of neoclassical processes and turbulence is eventually needed to capture the full picture of dynamical processes in fusion devices and provide reliable predictions.

Particle simulations constitute a powerful tool for the studies of neoclassical physics Lee (1983); Lin, Tang, and Lee (1995); Wang et al. (2006a); Bergmann, Peeters, and Pinches (2001). The neoclassical simulation can also provide a more consistent initial condition for turbulence simulations with a neoclassical particle distribution Wang et al. (2006b); Vernay et al. (2010). In this work, we focus on the implementation of the particle scheme in an unstructured mesh, following our previous work of the TRIMEG (TRIangular MEsh based Gyrokinetic) code Lu et al. (2019). The intention of this paper is to demonstrate the capabilities of the TRIMEG code and benchmark the implementation of the collision operator. The Monte-Carlo integration for the calculation of the particle/energy fluxes and bootstrap current is developed in unstructured meshes. In addition, the traditional local neoclassical theory is compared with the global particle simulations. While the eventual goal is the full treatment of multiple species with the neoclassical electric field, in this work, we focus on electron transport without coupling to the neoclassical electric field and ion dynamics.

The remainder of this paper is organized as follows. In Sec. II, the basic equations are presented. We start from general equations and their reduction to those that are implemented in the TRIMEG code. The gyrokinetic equations and the collision operator are listed. Specific issues such as the flux surface average and toroidal average in an unstructured mesh are discussed. In Sec. III, three typical cases with increasing complexity are studied. The local theoretical results of particle/energy fluxes and bootstrap current are calculated and compared to those obtained from the gyrokinetic simulations. The agreement and the discrepancy between the local theory and the global particle simulations, as well as the reasons or the possible limitation of the local theory, are discussed. Conclusions as well as an outlook to future work are given in Sec. IV.

II Models and Equations

II.1 General equations in $\delta f$ simulations

Using the $\delta f$ scheme, the distribution function $f$ is split into a fixed part $f_{0}$ and a perturbation $\delta\!f$, namely, $f=f_{0}+\delta f$. $f_{0}$ is a known distribution and $\delta f$ is obtained from

[TABLE]

where $C(f,f)$ is the collision operator, ${\bf R},v_{\|},\mu$ are the guiding center coordinates, the magnetic moment devided by the particle mass is adopted as one of the phase space coordinates defined as $\mu=v_{\perp}^{2}/(2B)$, $v_{\|}$ and $v_{\perp}$ are the velocities along and perpendicular to the magnetic field, respectively.

In this work, the Maxwell distribution is chosen ($f_{0}=f_{\rm{M}}$),

[TABLE]

where $T$ and $n$ are functions of the coordinates $(R,Z)$ in the poloidal cross-section, and thus

[TABLE]

where $\vec{\kappa}_{n,T,B}\equiv\nabla\ln\{n,T,B\}$, ${\bf\dot{R}}_{d}$ is the magnetic drift velocity, $\delta{\bf\dot{R}}$ is the perturbed velocity due to the wave field and the neoclassical electric field, $\dot{v}_{\parallel,0}$ and $\delta\dot{v}_{\|}$ are due to the mirror force and the wave/neoclassical fields, respectively. Furthermore, the perturbed parts $\delta\dot{\mathbf{R}}$ and $\delta\dot{v_{\parallel}}$ contain contributions from the turbulence as well as the neoclassical electric field and can be represented as $\delta\dot{\mathbf{R}}=\delta\dot{\mathbf{R}}_{turb}+\delta\dot{\mathbf{R}}_{E_{ZF}}+\delta\dot{\mathbf{R}}_{E_{nc}}$, where $\delta\dot{\mathbf{R}}_{turb}$ is the non axi-symmetric turbulence component while $\delta\dot{\mathbf{R}}_{E_{ZF}}$ and $\delta\dot{\mathbf{R}}_{E_{nc}}$ are related the $n=0$ axi-symmetric field due to turbulence and collisions between particles, respectively. For models with \Phi{\color[rgb]{0,0,0}{}_{nc}}=\Phi{\color[rgb]{0,0,0}{}_{nc}}(r), where $\Phi$ is the scalar potential and the radial-like coordinate $r$ is a function of the magnetic flux, we readily have $\delta\dot{v}_{\parallel}=0$, which is the assumption used in some neoclassical studies Hinton and Hazeltine (1976). More generally, the neoclassical scalar potential can vary in the poloidal direction and it is denoted as $\Phi_{nc}=\Phi_{nc}(R,Z)$. In the previous TRIMEG work Lu et al. (2019, 2021, 2023), neoclassical effects were ignored. For the neoclassical studies in this work, we ignore the $n\neq 0$ modes and we have

[TABLE]

where $\delta\dot{\mathrm{R}}_{E_{nc}}$ and $\delta\dot{v}_{\|E_{nc}}$ are the perturbed velocity and acceleration due to the neoclassical electric field, respectively. In this work, since we focus on electron transport, $E_{nc}$ is not taken into account, as adopted by previous work Lin (1996). In addition, since the electron Larmor radius is negligible compared to the scale length of the equilibrium gradient, the drift kinetic equations are solved. Note that Eqs. 4–5 are written in $({\bf R},v_{\|},\mu)$ coordinates, for demonstrating the general form of the models without and with neoclassical physics. The right-hand side can be also written in $({\bf R},E,\mu)$ coordinates, in order to use the constants of motion $(E,\mu)$, where $E$ is the energy. Further simplifications, e.g. omission of the electric field, is adopted to get the equation implemented in TRIMEG in this work, as shown in Eq. 14 in the next section.

II.2 Discretization of the distribution function

When representing the distribution function in simulations, we need to discretize it. In the $\delta f$ scheme, $f_{0}$ and $\delta f$ can be represented as follows with respective weight fields $P(z,t)$ and $W(z,t)$ Lanti (2019) :

[TABLE]

where $N_{ph}$ and $N$ are the physics particle number and marker number, respectively, $z$ denotes the phase space coordinate $({\bf R},v_{\|},\mu)$, whose Jacobian is $J(z)$. Furthermore, the time evolution for the weights in the collisionless case can be written as

[TABLE]

Assuming that the background distribution is of zeroth order in $\rho_{p}/R_{0}$, where $\rho_{p}$ is the poloidal gyro-radius, the time evolution for the weight equation can be written as Lin (1996)

[TABLE]

To take into account the collision effects, at each time step, the collision operators can be applied after the collisionless dynamics Lanti (2019). Theoretically, the collision operator is taken into account in particle simulations by solving the Langevin equation, as we will discuss in the following section.

II.3 Collision operator

Realistic non-liniear collision operators have been implemented in the XGC code recently Hager and Chang (2016). However, in our studies for simplicity we use the linearized collision operator. By treating the heavy species as infinitely massive with negligible thermal velocity, the collision operator becomes much simpler and the Lorentz collision operator is readily obtained Hazeltine and Waelbroeck (2019)

[TABLE]

In this work, only the electron-ion Lorentz collision operator is used, given by Eq. 12, which is written using the Monte-Carlo method as Lin (1996)

[TABLE]

where $\xi=v_{\parallel}/v$, $\nu$ is the collision frequency, and $\mathcal{R}$ is a uniform random number between 0 and 1 Lin (1996). Since the Lorentz model is valid in the limit of large ion-to-electron mass ratio, the focus of this work is the benchmark of the code with the theoretical results in the simplified case while more realistic collision operator needs to be implemented in future in order to compare with experiments.

II.4 Implemented Equations and normalizations

For our neoclassical studies, the TRIMEG code was modified as described in this section. In addition to the discussions in II.1, the $\delta f$ model was adopted and a weight equation was implemented as given by Eq. 11 in $(R,Z,\phi)$ coordinates

[TABLE]

where $\kappa$ is the negative gradient of the Logarithm of the equilibrium distribution function given by

[TABLE]

Taking into account the normalization given in table 1, the normalized equation implemented in the code is

[TABLE]

To take into account the collisional effects, the discretized Lorentz collision operator given by Eq. 13 is implemented in the code,

[TABLE]

where $\xi$ is the particle pitch as defined earlier, $\bar{\nu}=\nu\,t_{N}$ is the normalized collision frequency and $\bar{t}$ is the normalized time. Furthermore, regardless of the choice of {\color[rgb]{0,0,0}\nu}\,\Delta t, due to the choice of $\mathcal{R}$ being random, $|\xi|$ can become greater than 1, as shown in Fig. 1. This would cause a nonphysical solution with $|v_{\|}/v|>1$, which was fixed by re-setting $|\xi|$ equal to exactly 1 in those cases. Frequent re-setting of the pitch would break the uniformity of the random distribution due to over-sampling of $1,-1$. Hence, we want to avoid it by limiting the number of times the pitch becomes greater than 1. For large values of {\color[rgb]{0,0,0}\nu}\,\Delta t(>0.1{{\color[rgb]{0,0,0}5}}), more markers end up with $|\xi|>1$ after the collision, and the re-setting operation is needed more frequently. For small values of {\color[rgb]{0,0,0}\nu}\,\Delta t(<0.1{\color[rgb]{0,0,0}5}), only a small portion of markers enter the $|\xi|>1$ zone, and the re-setting operation is needed less frequently. In the simulations, the value for {\color[rgb]{0,0,0}\nu}\,\Delta t is chosen to be be small enough ({\color[rgb]{0,0,0}\nu}\,\Delta t{\color[rgb]{0,0,0}<0.15}) to avoid frequent use of re-setting. While in Fig. 1, the random number is close to $1$ ($\mathcal{R}=0.9$), the total portion of the $|\xi|>1$ particles is much smaller than $50\%$ for $\nu\Delta t=0.15$ since $\mathcal{R}\in[0,1)$.

II.5 Diagnosis for axisymmetric components in TRIMEG

In the studies of neoclassical transport, the particle and energy fluxes and the bootstrap current are all axisymmetric variables and need to be calculated numerically. An important routine implemented in the code during our studies is the flux surface average calculation. For any function $F$, the flux surface average is Hinton and Hazeltine (1976)

[TABLE]

where $\Delta V$ is the small volume between two adjacent flux surfaces, and $dS$ is the area element on the flux surface. For general geometry, we can substitute $\delta f$ with Eq. 7 and re-write the flux surface average given by Eq. 18 as

[TABLE]

After taking the integrals this is simplified to

[TABLE]

We can also readily calculate the flux surface average for 3D variables, namely,

[TABLE]

where $\psi$ is the poloidal magnetic flux, $\theta$ is a poloidal-like coordinate, and $\Delta\tilde{V}$ is an infinitesimal volume. The flux surface average of $F$ yields the same result in Eq. 20.

In some cases, the toroidal average is used for the analysis instead of the flux surface average (which is relevant for future applications with asymmetric poloidal structures):

[TABLE]

When plugging in the definitions for $F$, the average can be written as

[TABLE]

where if integrals over the toroidal direction and over an arbitrarily small volume are taken, this equation can be rewritten as

[TABLE]

where $\Delta\tilde{V}=\tilde{R}\Delta\tilde{S}\Delta\phi$ has been adopted.

In our analysis, we deal with the toroidal average or the flux surface average normalized by average density. Hence, if we take into account that $N_{ph}=V_{tot}\langle n\rangle_{V}$, where $V_{tot}$ is the total volume and $\langle n\rangle_{V}$ is the average density in this volume, the toroidal average equation can be rewritten as

[TABLE]

where $\Delta\tilde{V}=\Delta\tilde{S}\tilde{R}\Delta\phi$ was used. The flux surface average Eq. 20 can be also expressed in the same way. For the calculation of the annulus area in shaped tokamak geometry, the Monte-Carlo integration method is used and the annulus area is calculated at the beginning of the simulation when the markers are distributed uniformly.

II.6 Diagnosis and benchmark using the local neoclassical transport theory

The benchmark is done by comparing the simulation results to the theoretical local electron transport model Hinton and Hazeltine (1976); Lin (1996). We can also write simplified equations for electrons in a circular geometry, where ion charge is equal to one, and assuming constant pressure and temperature for the ion species. From the discussions in previous work Hinton and Hazeltine (1976), we can summarize these equations as follows

[TABLE]

where the dimensionless $K$ coefficients are fitted analytically using values for the numerical coefficientsHinton and Hazeltine (1976) for charge number equal to 1 with $\nu_{\ast}=\sqrt{2}\nu qR/(\epsilon^{3/2}v_{\rm th})=3\sqrt{\pi}\nu_{*e}/4$ :

[TABLE]

The normalized equations used in diagnostics are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the normalized quantities are denoted by a bar over the variables, and $\left<(J_{\parallel}-J_{\parallel s})/h\right>\equiv j_{b}$. In the asymptotic limits of collisionality, these formulae can also be rewritten as described in the paper by Lin Lin (1996). For the banana regime, in the limit of $\nu\rightarrow 0$, the following analytical equations for particle flux $\Gamma$, energy flux $Q$, and the bootstrap current $j_{b}$ can be obtained Lin (1996):

[TABLE]

where $\langle...\rangle\equiv\int_{0}^{2\pi}h{d\theta}/{(2\pi)}$ represents the flux surface average. Additionally, $\epsilon=r/R_{0}$, $h\equiv 1+\epsilon\cos\theta$, $\rho=mv_{th,e}c/(eB_{0})$, $q=rB_{0}/(R_{0}B_{p0})$, $p=(3/2)nT$, and to the lowest order in $\epsilon$,

[TABLE]

For the collisional limit $\nu qR/v_{\rm th}\gg 1$, the following analytical solutions are obtained Lin (1996)

[TABLE]

The normalized form of the fluxes and bootstrap current are as follows, for the banana regime,

[TABLE]

[TABLE]

[TABLE]

where $\bar{\kappa}_{(n,T)}=R_{N}\kappa_{(n,T)}$. For the collisional regime,

[TABLE]

[TABLE]

Comparisons to the simulation results are done by calculating values using the flux surface average described in section II.5 during the simulation, where the values for $X$ are chosen as $X=v_{r}mv^{2}/2$ for the energy flux, and $X=v_{r}$ for the particle flux, where $v_{r}=\vec{v}\cdot\nabla\psi/|\nabla\psi|$ is the radial drift velocity. The analytical formulas described in this chapter are implemented in the Fortran code, where the values for $\kappa_{n}$ and $q$ are interpolated from the density profile and the magnetic equilibrium. These results are described in sections III.

III Simulation setup and results

III.1 Electron transport results for the larger aspect ratio case

In this section, we will consider the International Tokamak Physics Activity (ITPA) case, which has been defined in the benchmark of the Toroidal Alfvén Eigenmode (TAE) driven by energetic particles (EPs) Könies et al. (2018). This is a Tokamak plasma with large aspect ratio ($a/R_{0}=0.1)$ and concentric circular magnetic surfaces. The on-axis magnetic field is $3$ Tesla. The major radius $R_{0}=10$ m. The nominal safety factor is $q=1.71+0.16(r/a)^{2}$ with low magnetic shear. In generating the EQDSK file for the ad-hoc equilibrium, the analytical form

[TABLE]

is adopted, where $\bar{q}=q\sqrt{1-r^{2}/R_{0}^{2}}$. A simplified match to the $q$ profile is used by letting

[TABLE]

which is a good approximation for moderate or large aspect ratios. In our study of neoclassical transport and bootstrap current generation, we only keep the electrons as the kinetic species and assume a uniform electron temperature. The electron density gradient is nonuniform and the following profile is used for solving the weight equation Lu et al. (2021)

[TABLE]

where $c_{0}=0.49123$, $c_{1}=0.298228$, $c_{2}=0.198739$, and $c_{3}=0.521298$. In the original benchmark study, this density profile is used as the EP density, while in our work we use it as the electron density. In gyrokinetic simulations, the temperature and density are given in terms of the Larmor radius and $\beta$ (the ratio of the plasma pressure to the magnetic pressure). While the nominal values are $\beta_{e}=8\pi n_{e}T_{e}/B_{0}^{2}=9\times 10^{-4}$, $\rho_{N}=$0.001,52\text{,}\mathrm{m} for the electromagnetic simulations, in this work, we adopt an enhanced density and temperature case with $\beta_{e}=0.03$, $\rho_{N}=7.8767\times 10^{-3}$\text{\,}\mathrm{m}. Note that the value of $\beta_{e}$, however, does not change the values of the flux and bootstrap current normalized by the equilibrium density for the present model.

We simulate cases with high and low collision frequencies and compare the results with the local theory. The time step sizes were chosen such that $\nu\,\Delta t$ ranges from $6.7\cdot 10^{-6}$ to $0.034$. During the simulation, the density profile changes gradually because of the particle transport due to collisions. At the end of the simulation, the density profiles are analyzed, as can be seen in Fig. 2. For the high collision case, the density change is of the order of $\sim 3\%$, while for the low collision case, the change is $\sim 2.5\%$ near the axis and $|\delta n|<1\%$ in the rest of the volume. This density variation has negligible effects on the density gradient as given in Eqs. 50 with the chosen coefficients and thus the particle and energy fluxes and the bootstrap current stay at the same level after a ramp-up phase until the end of the simulations. In addition, as our simulations are not in the long time scales we do not consider the known issues of unknown marker distribution in the Monte-Carlo delta-$f$ collision operators Chen, Cheng, and Parker (2022) Wang et al. (2004). However, for longer time scale simulations the deviation of the initial marker distribution for the initial one and the noise level also need to be systematically analysed in the future.

The radial profiles of particle flux, energy flux, and bootstrap current are also analyzed, as shown in Figs. 3-4. For the low collision frequency limit, we observe good agreement between the theoretical calculation and the simulation, as seen in Fig. 3. Blue dashed lines indicate the results from particle simulations. The green lines for the particle flux and energy flux indicate solutions in the low collisionality limit given by Eqs. 43-44, while the red lines are in the high collisionality limit given by Eqs. 46-47. Note that although the value of the collision is set to be as close as possible to the collisionless or collisional limit in the simulations, it is still not the limit of $\nu\rightarrow 0$ and $\nu\rightarrow\infty$, as adopted in the theoretical formulae, which can lead to minor or moderate discrepancies between the simulation and theoretical results, as shown in Figs. 3-4 and other cases in the Cyclone case and the ASDEX-Upgrade case. As for the bootstrap current, the green line indicates the low collision limit given by Eq. 45, and the orange line is given by the analytical formula by Hinton which takes into account finite collision frequency given in Eqs. 39 with coefficients in Eqs. 30 - 35. We observe a big difference in the energy and particle flux near the axis for low collisionality. The discrepancy between the standard neoclassical theory and the simulation is also observed in the study of ion transport in previous work Lin, Tang, and Lee (1997). Better agreement can be obtained in principle by considering the theoretical formula closer to the axis following the derivation in the neoclassical theory Hinton and Hazeltine (1976); Chang and Hinton (1982); Sauter, Angioni, and Lin-Liu (1999) or by comparing with other codes Wang et al. (2004); Lin, Tang, and Lee (1997); Belli and Candy (2008) but it is beyond the scope of this work.

As for the high collisionality case, the results are given in Fig. 4. For the particle flux and the energy flux, a very good agreement with the analytical solutions is observed. The bootstrap current is much lower than that in the banana limit and thus more markers ($10^{7}$) are used in this simulation to enhance the signal-to-noise ratio. The bootstap current from the simulation is of the same order of magnitude as but higher than the analytical result.

We also investigate the fluxes and bootstrap current as a function of the collision frequency as shown in Fig. 5. The reference radial location was chosen arbitrarily to be $r\approx$0.23\text{,}\mathrm{m}$$, and the fluxes and the current are compared with the values from the theoretical radial profile at this reference location. We observe that the bootstrap current from simulations decreases much faster as collision frequency increases than expected from the analytical interpolation solution in the plateau and the collisional regimes. However, we observe good agreement for the particle and energy fluxes with the analytical solutions, in the low ($\nu^{*}<1$) and high ($\nu^{*}>10$) collisionality limits. Furthermore, the plateau ($1<\nu^{*}<10$) regime is also visible where the fluxes stay almost constant for different collision frequencies, and the overall behavior is qualitatively consistent with previous results Hinton and Hazeltine (1976); Lin (1996).

III.2 Electron transport results for the moderate aspect ratio case

In this section, we study the electron transport for the moderate aspect ratio case featured by $a/R_{0}=0.36$. The details of the magnetic equilibrium and the density profile are listed in the reference Lu et al. (2019) and the main parameters are briefly summarized as follows. The major radius and the minor radius are $R_{0}=1.67$ m and $a=0.6012$ m, respectively. In generating the EQDSK file as the input of the magnetic equilibrium, the $\bar{q}$ profile in Eq. 48 is adopted with

[TABLE]

where $r_{c}=0.5a$. The radial profile of the density and its gradient are given analytically as follows Lu et al. (2019),

[TABLE]

where $L_{ref}=R_{0}$, $W_{n}=0.3$, $\kappa_{n}=2.23$. $\beta_{e}=0.03$, $\rho_{N}=7.8767\times 10^{-3}$ m. For the sake of simplicity, a uniform temperature profile is assumed. Compared with the ITPA case ($R_{0}/a=10$, $q\approx 1.75$ at $r/a=0.5$) studied in Chapter III.1, the aspect ratio of the Cyclone case is smaller ($a/R_{0}=0.36$). In addition, the safety factor is lower than that of the ITPA case in the inner radial region ($r<0.5a$) but is larger near the edge. In the theoretical derivation, the small parameter $\rho_{p}=(qR_{0}/r)m_{s}v_{\perp}/(Z_{s}eB)$ is used as the expansion parameter, where $Z_{s}$ is the charge number, and the subscript ‘s’ indicates species ‘s’. As a result, for different values of $q$ and $r/R_{0}$, the accuracy of the theoretical result can be different, which is more relevant for ion transport. As $\rho_{p}$ is close to or even larger than the characteristic length of the equilibrium or the density/temperature profiles, the traditional neoclassical formulae are not valid and corrections are needed as shown in previous studies of ion transport Chang and Hinton (1982); Helander (1998). For electron transport, $\rho_{p}\ll 1$ is usually well satisfied except if it is very close to the magnetic axis where $R_{0}/r\rightarrow\infty$ and the traditional theoretical formulae can also break down. The time step sizes were chosen such that $\nu\,\Delta t$ ranges from $1.2\cdot 10^{-4}$ to $0.12$, which is the highest value among the three cases, hence in the highest collisionality case non-physical effects from the implementation of the collision operator need to be considered. The verification of the convergence of this $\nu\,\Delta t$ value is shown on the right-hand side of Fig.15.

In our simulation, we observe changes in the density of the order of $0.01\%$ in the low collision frequency case, and about $2\%$ in the high collisional case, as shown in Fig. 6. We also observe large density changes on the axis, which can be due to the discontinuity of the density profile near the axis (the radial gradient of the density profile is not exactly zero at the axis according to Eq. 54 and can cause non-physical $\partial\ln n/\partial R$ and $\partial\ln n/\partial Z$ values in the 2D interpolation). The positive $\delta n$ near the axis for the low collisionality can be due to the different parameter conditions near the axis such as $r/R_{0}$ and the consequent physical collisionality. Nevertheless, the on-axis physics is not the focus of this work and it does not affect the physics near the middle radius due to the small orbit width of electrons.

We again start by analyzing the radial particle flux, energy flux, and bootstrap current profiles of high and low collision frequency cases, as shown in Figs. 7–8. For the low collisionality case, the agreement is not as good as that in the ITPA large aspect ratio case. The reason can be the large aspect ratio approximation adopted in the analytical formulas and the different treatments in our code for the poloidal-angle dependent values such as $q$ in calculating the theoretical fluxes and the current. It merits more effort in the future to adopt analytical formulas for comparision which are valid for broad regimes of aspect ratiosChang and Hinton (1986) and for general axisymmetric equilibria and all collisionality regimes Sauter, Angioni, and Lin-Liu (1999). The local theory for the comparison in this work is derived in the large aspect ratio limit and the good agreement is observed in this parameter regime in the ITPA case. In addition, when calculating the theoretical values of the fluxes and current, we calculate the theoretical values on numerous points in one annulus; then the averaged value is calculated with proper weights and thus the flux surface average value is obtained using the Monte-Carlo integration. This method of calculating the flux-surface-averaged fluxes and current gives us a convenient and practical way of calculating the fluxes and current for shaped tokamak geometry, as we also adopted for the ASDEX Upgrade case in the next section. For the high collision frequency case, the agreement is better closer to the axis, which is expected as the local aspect ratio is larger. The bootstrap current is much larger than the analytical solution. The reason can be the approximation in the interpolation formula Eq. 39, which is derived to match the results at the low and high collision regime. Indeed, the discrepancies between the transition formula and the formulae at the low and high collisions are also observed theoretically Hinton and Hazeltine (1976). In the theoretical solutions, the interpolation formula is obtained by fitting the analytical results in the banana-plateau and plateau-collisional regimes, hence our simulation result is expected to be more exact than the analytical solution close to the plateau regime. However, there are analytical theories

Analyzing how fluxes change with different collision frequencies, we take the radial location close to the maximum value of the fluxes ($r=$0.21\text{,}\mathrm{m}$$), which in this case is also a larger aspect ratio case. We observe good agreement with the analytical results for the particle and energy fluxes, as shown in Fig. 9. The bootstrap current from the simulation follows the interpolation formula (red line) well. For $\nu^{*}>10^{2}$, the discrepancy between the simulation and the interpolation formula is larger and it can be due to the inaccuracy of the interpolation formula. For the bootstrap current, the lowest collision frequency case does not produce the highest current as would be expected. This is due to the bootstrap current still increasing slightly when the simulation is finished. We keep the lowest collision case as the indicator of the most expensive case that we can study with our present capability. Note that in this study of the dependency on collisionality, the collision frequency varies by a factor of $10^{5}$ and the needed minimum simulation time also varies significantly since several ($\sim 10$) collisional periods are needed to reach the saturated state of the fluxes and the current. The lowest collisional case is the most expensive one in order to observe reasonable results. The current simulation consumes 14 hours of computation time on 4 nodes with each node containing two Intel(R) Xeon(R) Gold 6130 processors (16 cores per processor, 2.10GHz, 22MB Cache). Both the simulation results and the theoretical interpolation formula results are lower than the collisionless approximation and have the expected behavior of decreasing when the collision frequency is increasing.

III.3 Electron transport results for the ASDEX Upgrade case

In this section, a realistic geometry for Tokamak plasmas is used. The ASDEX Upgrade (AUG) case with shot number 34924 at 3.600 s is chosen as adopted by the previous work for the development of the TRIMEG code for the studies of the ion temperature gradient mode Lu et al. (2019). This is a typical discharge for the study of energetic particles and turbulence physics Lauber et al. (2018). The EQDSK file is obtained from experimental data. The major radius $R_{0}\approx$1.71\text{,}\mathrm{m}. The $q$ profile and the poloidal magnetic flux function are shown in Fig. [10](#S3.F10). In the simulation, we use the experimental equilibrium but uniform temperature and the analytical density profiles in Eq. [53](#S3.E53) with the radial coordinate replaced with $\rho_{pol}=\sqrt{(\psi-\psi_{0})/(\psi_{b}-\psi_{0})}$, where the subscript [math] and $b$ indicate the values at the magnetic axis and the last closed surface, respectively. $\beta_{e}=0.03$, $\rho_{N}=$0.01\text{\,}\mathrm{m}. In this work, we focus on testing the capability of the code in treating realistic geometry with minimum technical complexity. The time step sizes were chosen such that $\nu\,\Delta t$ ranges from $2.78\cdot 10^{-5}$ to $0.0278$.

As in previous chapters, we first look at the density changes due to low and high collision frequencies as shown in Fig. 11. The density changes are much smaller than in the previous cases. The density change due to high collision frequency is $0.6\%$, while in the low collision frequency case, it is almost negligible.

Looking at the radial profile for the low collision case given in Fig. 12, we observe good agreement with the theory. The bootstrap current and the particle flux agree very well, while the energy flux has a higher value at the center than predicted by the neoclassical theory. As for the high collision frequency case, given in Fig. 13, the discrepancy between theory and simulation is larger than that in the ITPA and Cyclone cases. In addition to the reason we discussed in the Cyclone case, for the AUG case, the magnetic flux surfaces are not circular but are strongly shaped. The theoretical formulae of fluxes and bootstrap current were derived for the circular magnetic flux surface originally Hinton and Hazeltine (1976) and our scheme of flux surface average of the fluxes and current is one possible way of an estimate, which is more reasonable for circular magnetic flux surfaces and is for the verification of the implementation of the code. More accurate theoretical/numerical solutions for shaped tokamak plasmas can be found elsewhere Sauter, Angioni, and Lin-Liu (1999); Belli and Candy (2008) and the comparison with our simulation results is possible, but is beyond the scope of this work.

Furthermore, to investigate the dependence on the collision frequency we picked the radial coordinate where the highest values of fluxes were observed ($r_{LFS}=$0.2\text{,}\mathrm{m}$$), and compared the analytical and simulation results, as shown in Fig. 14. The agreement with the theory can still be seen. However, the agreement for large collision frequencies is worse. Nevertheless, the trends of the fluxes and bootstrap current follow the theoretical results. More issues can be studied in the future for understanding the connections and the differences between the global gyro-kinetic simulation and the local theory, in order to identify the origin of the differences between the theoretical results and the simulation results as well as the limitation of the local theory. For future studies, comparing results to other gyrokinetic codes, such as ORB5 Lanti (2019) or GENE Görler et al. (2012) would be the next step for verification in specific cases. Aditionally, actual experimental electron tranport levels in the ASDEX Upgrade can also be compared Ryter et al. (2003); Meyer et al. (2019), however this merits more effort as contributions from both turbulence and neoclassical physics need to be considered.

IV Conclusion and outlook

We studied the electron transport and bootstrap current generation by adding a pitch angle scattering operator to the TRIMEG code, which uses an unstructured mesh and equations in $(R,Z)$ coordinates. This provides a robust tool in a broad collisionality range and for flexible parameters such as tokamak geometry. In this work, we only considered the electron species, set the temperature gradient to zero, but take into account the density gradients. We first compared the simulation results to analytical calculations in the large aspect ratio ($R_{0}/a=10$) approximation and found good agreement in those limits between theory and simulation. For moderate aspect ratio ($R_{0}/a=1/0.36$) cases the agreement was worse, as expected, due to the approximations in the theory such as the large aspect ratio approximation. In the case with AUG experimental geometry, in the collisional regime, the radial profiles of bootstrap current and fluxes are simulated, with significant discrepancy between the local theory and the simulation results. In the banana regime, the agreement between the theory and the simulation results are reasonably good. In the studies of collisionality scan, for high collision frequencies, the energy and particle fluxes decrease in magnitude faster than the analytical solutions as the collision frequency increased.

This work has verified the capability of the TRIMEG code to study the electron transport and the bootstrap current generation in tokamak plasmas for further simulations with open field lines. Future steps for further investigation would be to investigate the neoclassical physics in realistic experimental geometry with experimental density and temperature profiles. It also merits more efforts to add like-particle collision operators Lin, Tang, and Lee (1995); Wang et al. (2004) and fully non-linear collision operators Hager and Chang (2016) for more realistic studies. The particle-field coupling also needs to be added, for the study of the neoclassical radial electric field in future, which plays a key role in the instability stabilization and the optimization of the confinement performance.

Acknowledgements.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 – EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

Appendix A Results related to numerical convergences

The main concern in the numerical convergence is related to the marker number and the time step size. We select two typical cases to demonstrate the principles we adopted when choosing values of the time step size and the marker number. While it is not practical to perform convergence studies for each case, the bootstrap current, particle flux and energy flux are compared for different values of the relevant parameters. As shown in the left frame of Fig. 15, we calculate the particle flux at the middle minor radius by interpolation ($r/a=0.5$). For the ITPA case ($\bar{\nu}\approx 3$, $\nu\Delta t=0.00581$), results are compared for different marker numbers as shown in the left frame. The results start to converge as $N_{p}\geq 10^{6}$. As showin in the right frame of Fig. 15, the particle fluxes are compared with different time step sizes. For this ASDEX Upgrade case ($\bar{\nu}\approx 7\cdot 10^{1}$), the results start to converge as $\nu\Delta t\leq 0.11133$. In our simulations, we choose $\nu\Delta t\leq 0.055664$ for $\bar{\nu}\approx 7\cdot 10^{1}$ and smaller $dt$ for other cases of ASDEX Upgrade.