Numerical calculation of the neoclassical electron distribution function in an axisymmetric torus

We solve for a stationary, axisymmetric electron distribution function ($f_e$) in a torus using a drift-kinetic equation (DKE) with complete Landau collision operator. All terms are kept to gyroradius and collisionality orders relevant to high- temperature tokamaks (i.e., the neoclassical banana regime for electrons). A solubility condition on the DKE determines the non-Maxwellian pieces of $f_e$ (called $f_{NMe}$) to all relevant orders. We work in a 4D phase space $(\psi, \theta, v, \lambda)$, where $\psi$ defines a flux surface, $\theta$ is the poloidal angle, $v$ is the total velocity, and $\lambda$ is the pitch angle parameter. We expand $f_{NMe}$ in finite elements in both $v$ and $\lambda$. The Rosenbluth potentials, $\Phi$ and $\Psi$, which define the collision operator, are expanded in Legendre series in $\cos \chi$, where $\chi$ is the pitch angle, Fourier series in $\cos \theta$, and finite elements in $v$. At each $\psi$, we solve a block tridiagonal system for $f_{NMe}$, $\Phi$, and $\Psi$ simultaneously, resulting in a neoclassical $f_e$ for the entire torus. Our goal is to demonstrate that such a formulation can be accurately and efficiently solved numerically. Results will be compared to other codes (e.g., NCLASS, NEO) and could be used as a kinetic closure for an MHD code (e.g., M3D-C1).