Neoclassical transport in strong gradient regions of large aspect ratio tokamaks

Abstract

We present a new neoclassical transport model for large aspect ratio tokamaks where the gradient scale lengths are of the size of the ion poloidal gyroradius. Previous work on neoclassical transport across transport barriers assumed large density and potential gradients but a small temperature gradient, or neglected the gradient of the mean parallel flow. Using large aspect ratio and low collisionality expansions, we relax these restrictive assumptions. We define a new set of variables based on conserved quantities, which simplifies the drift kinetic equation whilst keeping strong gradients, and derive equations describing the transport of particles, parallel momentum and energy by ions in the banana regime. The poloidally varying parts of density and electric potential are included. Studying contributions from both passing and trapped particles, we show that the resulting transport is dominated by trapped particles. We find that a non-zero neoclassical particle flux requires parallel momentum input which could be provided through interaction with turbulence or impurities. We derive upper and lower bounds for the energy flux across a transport barrier in both temperature and density and present example profiles and fluxes.

1. Introduction

The pedestal, and transport barriers in general, play an important role in tokamak performance (Wagner et al. 1984; Greenfield et al. 1997) and thus it is useful to find a comprehensive transport model for these regions. In pedestals, for example, strong gradients of temperature, density and radial electric field of the order of the inverse ion poloidal gyroradius are observed (Viezzer et al. 2013). Moreover, it has been found that the ion energy transport in pedestals is close to the neoclassical level (Viezzer et al. 2018). Measurements of H-mode pedestals in Alcator C-Mod (Theiler et al. 2014; Churchill et al. 2015) and Asdex-Upgrade (Cruz-Zabala et al. 2022) have shown poloidal variations of density, electric field and ion temperature that cannot be explained using standard neoclassical theory. It is thus desirable to extend neoclassical theory for stronger gradients, and logical to choose the ion poloidal gyroradius as the characteristic scale length. Comparisons of experimental data with standard neoclassical theory (Hinton & Hazeltine 1976) such as the one by Viezzer et al. (2018) miss finite poloidal gyroradius effects.

Setting the scale length in transport barriers to be the poloidal gyroradius implies that the poloidal component of the $E \times B$-drift in large aspect ratio tokamaks becomes of the order of the poloidal component of the parallel velocity. As a result, a strong radial electric field shifts the trapped–passing boundary (Shaing, Hsu & Dominguez 1994a), and causes an exponential decrease proportional to the radial electric field in plasma viscosity (Shaing et al. 1994a) and radial heat flux (Kagan & Catto 2010; Shaing & Hsu 2012). The mean parallel flow is also affected by a strong radial electric field and can change direction (Kagan & Catto 2010). A strong shear in radial electric field causes orbit squeezing, which reduces the heat flux and increases the trapped particle fraction for increasing radial electric field shear (Shaing & Hazeltine 1992; Shaing, Hsu & Hazeltine 1994b).

Combining all these effects, Shaing & Hsu (2012) calculated the heat flux and mean parallel velocity but they neglected the strong mean parallel velocity gradient and the poloidal variation of the electric potential. Kagan & Catto (2008) and Catto et al. (2013) have likewise developed extensions to neoclassical theory to allow for stronger density gradients to calculate fluxes. In Kagan & Catto (2010) and Catto et al. (2011, 2013), the density gradient was taken to be steep but the temperature gradient scale length had to be much larger than the ion orbit width. Furthermore, they assumed a quadratic electric potential profile and also neglected the poloidal variation of the potential.

Comparisons between analytical solutions and simulations have been carried out by Landreman et al. (2014), which demonstrated the significance of source terms.

We will assume that the gradient length scale of potential, density and temperature is of the order of the poloidal gyroradius and we will retain the poloidal variations of density and potential. Assuming a large aspect ratio tokamak with circular flux surfaces in the banana regime and including unspecified sources of particles, parallel momentum and energy, we find equations for the ion distribution function, and a set of transport relations for ions.

In § 2, we justify our choice of orderings physically, and we motivate our choice of sources of particles, momentum and energy by considering the transition from the core into a transport barrier. A more detailed discussion of trapped and passing particles follows in § 3, where the shift of the trapped–passing boundary is derived and a new set of variables based on conserved quantities is introduced. In § 4 we calculate the ion distribution function in the trapped–barely passing and freely passing regions. We also calculate the poloidally varying part of density and potential. The solvability conditions for the equation containing the distribution function of the bulk ions are the density, parallel momentum and energy conservation equations, calculated in § 5. The ion transport equations are discussed further in § 6. We find that a non-zero parallel momentum input is required to sustain a neoclassical particle flux and consider the possibility of interaction with turbulence. For the energy flux, we derive upper and lower bounds and relate the gradient lengths of temperature and density to the growth of neoclassical energy flux as one moves into the transport barrier. We conclude by presenting some example profiles for the ‘high flow’ case and the ‘low flow’ case. A summary of our results is given in § 7.

2. Orderings and phase space outline

In this paper we consider the transition from regions with large turbulent transport into strong gradient regions. In a region of large turbulent transport, for example the core, neoclassical transport gives a minor contribution because turbulent transport carries most particles, momentum and energy. With the transition into a regime of low turbulence, like a transport barrier, the same total fluxes must be kept but as turbulence decreases, we anticipate that the turbulent transport goes down, too, and instead the fluxes must be picked up by neoclassical transport. Thus, we expect a rise in neoclassical fluxes at the transition from core to, for example, a pedestal (see figure 1). This argument is consistent with the observation that the energy flux in the pedestal is close to its neoclassical value (Viezzer et al. 2018). We will see, however, that this simple picture of the top of a transport barrier has limitations. In § 6.1 we find constraints that prevent the neoclassical fluxes from growing with radius.

Figure 1. The total flux must be kept constant across the core and pedestal. The neoclassical contribution increases in the pedestal whereas the turbulent fluxes decrease as turbulence quenches. There is the possibility of interaction between turbulent and neoclassical transport in the pedestal.

Turbulence and neoclassical transport could interact in the transport barrier and hence we need to include a source $\varSigma$ in the neoclassical picture. This source represents any possible input from turbulence as well as external injection of particles, momentum and energy. The source must balance the neoclassical fluxes $\varSigma /f \sim n^{-1} |\boldsymbol {\nabla } \psi |(\partial \varGamma /\partial \psi )\sim n^{-1}T^{-1}|\boldsymbol {\nabla } \psi |(\partial Q/\partial \psi )$, where f is the distribution function, $\varGamma$ is the neoclassical particle flux, $Q$ is the neoclassical energy flux, $n$ is the density, $T$ is the ion temperature and $\psi$ is the poloidal flux divided by $2{\rm \pi}$, which we use as a flux surface label. To estimate the size of $\varSigma$, we need the size of the neoclassical particle and energy fluxes. We consider trapped and passing particles separately.

We can estimate the contributions from trapped and passing particles to particle and energy transport by making random walk estimates. The diffusion coefficient $D$ for a random walk is $D \sim (\Delta x)^2/\Delta t$, with $\Delta x$ and $\Delta t$ the random walk size and time, respectively. The neoclassical particle flux is thus

(2.1)\begin{equation} \varGamma \sim \frac{(\Delta x)^2}{\Delta t}\frac{n}{L_n}, \end{equation}

where $L_n=|\boldsymbol {\nabla }\ln n|^{-1}$. In a large aspect ratio tokamak, where $r/R\sim \epsilon \ll 1$, $r$ is the minor radius and $R$ is the major radius, the poloidal gyroradius is much bigger than the gyroradius. For passing particles we will show that the orbit widths are $\Delta x\sim \epsilon \rho _p$, where $\rho _p=qR\rho /r$ is the ion poloidal gyroradius, $q$ is the safety factor and $\rho$ is the ion gyroradius. The time between collisions is $\Delta t\sim 1/\nu$, where $\nu$ is the collision frequency. The gradient of density is assumed to be of the order of the poloidal gyroradius and so the particle flux due to passing particles is

(2.2)\begin{equation} \varGamma_p \sim (\epsilon \rho_p)^2\nu \frac{n}{\rho_p}\sim \epsilon q \nu n \rho. \end{equation}

The orbit width for trapped particles will turn out to be $\Delta x \sim \sqrt {\epsilon } \rho _p$, the collisional time is $\Delta t \sim \epsilon /\nu$ and again the density gradient length is $L_n\sim \rho _p$. The fraction of trapped particles in phase space is only $\sim \sqrt {\epsilon }$, and with that we arrive at a neoclassical particle flux due to trapped particles of order

(2.3)\begin{equation} \varGamma_t \sim \sqrt{\epsilon}(\sqrt{\epsilon} \rho_p)^2\frac{\nu}{\epsilon} \frac{n}{\rho_p}\sim \frac{q}{\sqrt{\epsilon}}\nu n \rho. \end{equation}

A comparison of the transport contribution from passing and trapped particles shows that the particle flux due to trapped particles is much larger,

(2.4)\begin{equation} \frac{\varGamma_p}{\varGamma_t}\sim\epsilon^{3/2}\ll1. \end{equation}

The same estimate can be performed for the neoclassical energy flux when substituting the energy gradient $nT/L_T$ for the particle density gradient $n/L_n$, where $L_T\sim L_n\sim \rho _p$. In § 5 we find transport equations that are consistent with this estimate and show that transport is dominated by trapped particles.

Using the sizes of particle and energy flux above, we can now give an estimate for the source $\varSigma$ that we have to introduce in the kinetic equation to mimic turbulence, particle, momentum and energy sources. The gradient of the particle flux is

(2.5)\begin{equation} |\boldsymbol{\nabla} \psi| \frac{\partial{\varGamma}}{\partial{\psi}}\sim\frac{\varGamma_t}{\rho_p}\sim \sqrt{\epsilon} n \nu, \end{equation}

and hence we include a source

(2.6)\begin{equation} \varSigma \sim \sqrt{\epsilon} \nu f. \end{equation}

The random walk estimate of fluxes including source terms is accurate in the region of strong gradients but it should be noted that, for weak gradients, random walk arguments overestimate the neoclassical particle fluxes due to constrains imposed by intrinsic ambipolarity. Intrinsic ambipolarity (Sugama & Horton 1998; Parra & Catto 2009; Calvo & Parra 2012) is a property of neoclassical and turbulent particle fluxes in perfectly axisymmetric tokamaks: these particle fluxes give zero radial current to lowest order in an expansion in $\rho /r$ regardless of the value of the radial electric field. This property is only satisfied when the gradient length scales are much larger than the ion poloidal gyroradius. When the gradient length scales are of the order of the ion poloidal gyroradius and sources are included, the intrinsic ambipolarity constraint is relaxed as is found in this work and before in Landreman et al. (2014). We will find the ion neoclassical particle flux to be non-vanishing to lowest order in the presence of a parallel momentum source and discuss these effects in more detail in § 6.1.

3. Fixed-$\theta$ variables

To calculate the particle orbits, we introduce a new set of variables: the fixed-$\theta$ variables, which are based on the conserved quantities energy $\mathcal {E}$, canonical angular momentum $\psi _{\ast}$ and magnetic moment $\mu$,

(3.1a–c)\begin{equation} \mathcal{E} = \frac{1}{2}v^2+\frac{Ze\varPhi}{m},\quad \psi_{\ast}=\psi-\frac{I v_{\parallel}}{\varOmega},\quad \mu=\frac{v^2_\perp}{2B}. \end{equation}

Here, $v$ is the ion velocity, $m$ is the ion mass, $Ze$ is the charge, $\psi$ is the flux function, $\varOmega$ is the Larmor frequency and $B$ is the magnetic field strength. The electric potential is $\varPhi =\phi +\phi _\theta$. The piece $\phi$ is a flux function, $\phi =\phi (\psi )$, and its size is given by $e\phi /T \sim 1$, whereas $\phi _\theta$ is the small poloidally varying part of the electric potential, so $\phi _\theta =\phi _\theta (\psi,\theta )$ and $e\phi _\theta /T\sim \epsilon$. Here, $\theta$ is the poloidal angle. Throughout this work we will use that the electric potential is of the form

(3.2)\begin{equation} \phi_\theta(\psi,\theta)=\phi_c(\psi)\cos\theta, \end{equation}

which we will prove to be true in the banana regime for circular flux surfaces in § 4.4. Energy, canonical angular momentum and magnetic moment are constant in time, so following the trajectory of a single particle, we find

(3.3)\begin{equation} \frac{1}{2}v_{\parallel}^2+\mu B+\frac{Ze}{m}\varPhi(\psi, \theta)= \frac{1}{2}v_{{\parallel} f}^2+\mu B_f+\frac{Ze}{m}\varPhi(\psi_f, \theta_f), \end{equation}

and

(3.4)\begin{equation} \psi-\frac{Iv_{\parallel}}{\varOmega}=\psi_f-\frac{I v_{{\parallel} f}}{\varOmega_f}, \end{equation}

where the subscript $f$ indicates the values of the respective quantities at a fixed poloidal angle $\theta _f$, which represents a reference point in the orbit of the particle. It is important to note that $\psi _f$ and $v_{\parallel f}$ are constants for each particle. For example, following the trajectory of a passing particle, its velocity will deviate from $v_{\parallel f}$, but, having assumed the conservation laws above, the particle returns to its initial position $\psi _f$ with the velocity $v_{\parallel f}$ after one complete poloidal turn. Another particle on a different orbit will have a different $v_{\parallel f}$ and $\psi _f$. Hence, the fixed-$\theta$ quantities can be understood as labels of orbits and will be used as new phase space variables later on. The angle $\theta _f$ is left as a choice at this point, because choosing $\theta _f=0$ only captures particles that are trapped on the low-field side whereas setting $\theta _f={\rm \pi}$ captures particles trapped on the high-field side. We show in Appendix D.1 that it is important to take both sides into account when calculating trapped particle effects.

Using the standard large aspect ratio, circular-flux-surface tokamak, we can write the magnitude of the magnetic field as

(3.5)\begin{equation} B\simeq B_0\left(1-\frac{r}{R}\cos\theta\right) \end{equation}

to first order in the inverse aspect ratio $\epsilon$. Here, $B_0$ is the magnetic field on the magnetic axis. For $\theta _f=0$, the magnetic field is

(3.6)\begin{equation} B\simeq B_f\left[1+\frac{r}{R}(1-\cos\theta)\right], \end{equation}

with $B_f=B_0(1-r/R)$, whereas for $\theta _f={\rm \pi}$ the magnetic field can be written as

(3.7)\begin{equation} B\simeq B_f\left[1-\frac{r}{R}(1+\cos\theta)\right], \end{equation}

with $B_f=B_0(1+r/R)$. Changing $\theta _f$ from $\theta _f=0$ to $\theta _f={\rm \pi}$ causes a jump in $B_f$ of $O(\epsilon )$. It will be important in Appendix D that this difference is small.

In transport barriers, strong gradients in density, pressure and electric potential are observed. We will assume that $L_n\sim L_T\sim L_\varPhi \sim \rho _p$. Ordering the characteristic length of the transport barrier to be of the order of the poloidal gyroradius implies that the poloidal component of the $E\times B$-drift is of the same order as the poloidal component of the parallel velocity. The poloidal component of the $E\times B$-drift is

(3.8)\begin{equation} \frac{c}{B}(\boldsymbol{E}\times\boldsymbol{\hat{b}})\boldsymbol{\cdot}\boldsymbol{\nabla} \theta= \frac{cI}{B}\frac{\partial{\varPhi}}{\partial{\psi}}\boldsymbol{\hat{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta\equiv u \boldsymbol{\hat{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta+ \frac{cI}{B}\frac{\partial{\phi_\theta}}{\partial{\psi}}\boldsymbol{\hat{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta. \end{equation}

Here, $E=-\boldsymbol {\nabla }\varPhi$ is the electric field, $c$ is the speed of light and $\boldsymbol {\hat {b}}=\boldsymbol {B}/B$, where the magnetic field is $\boldsymbol {B}=I\boldsymbol {\nabla } \zeta +\boldsymbol {\nabla }\zeta \times \boldsymbol {\nabla }\psi$ and $\zeta$ is the toroidal angle. We have defined the velocity

(3.9)\begin{equation} u=\frac{cI}{B}\frac{\partial{\phi}}{\partial{\psi}}. \end{equation}

Note that we use $\varDelta \psi \sim Iv_t/\varOmega$ and thus $u\sim v_{t}$, where $v_t$ is the thermal speed. Due to our choice of ordering, $u$ and the parallel velocity $v_{\parallel}$ are of the same size. The poloidal velocity in this case is

(3.10)\begin{equation} \left(v_{\parallel} \boldsymbol{\hat{b}}+\frac{c}{B}\boldsymbol{E}\times\boldsymbol{\hat{b}}\right) \boldsymbol{\cdot}\boldsymbol{\nabla}\theta\simeq\left(v_{{\parallel}}+u\right) \hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta. \end{equation}

Particles are trapped on banana orbits if their poloidal velocity goes to zero at any point on their orbit. In the case of strong radial electric field this requires $v_{\parallel} +u=0$ instead of the usual trapping condition $v_{\parallel} =0$, as was first argued by Shaing et al. (1994a). It follows that particles with a parallel velocity close to $-u$, where $u$ is not necessarily small, are trapped. It has been previously shown that in this case the width of the trapped–barely passing region in velocity space is ${\sim }\sqrt {\epsilon } v_t$ (Shaing & Hazeltine 1992). We re-derive this result by calculating the deviations in radial position and velocity of particles on trapped orbits in Appendix A. Passing particles do not get reflected. One can divide the phase space into the freely passing region where $v_{\parallel} +u\sim v_t$ and the trapped–barely passing region $v_{\parallel} +u\sim \sqrt {\epsilon }v_t$.

For freely passing particles, we show in Appendix A.1 that $v_{\parallel} -v_{\parallel f}\sim \epsilon v_t$ and $\psi -\psi _f\sim \epsilon \rho _p R B_p$, where $B_p$ is the poloidal magnetic field. Thus, the deviations in parallel velocity and radial location are small in $\epsilon$. The deviations become large and diverge when $v_{\parallel} +u$ becomes small. This is the trapped–barely passing region. For trapped–barely passing particles, the differences are still small but larger by $\sqrt {\epsilon }$, so $v_{\parallel} -v_{\parallel f}\sim \sqrt {\epsilon } v_t$ and $\psi -\psi _f\sim \sqrt {\epsilon } \rho _p R B_p$ as can be found in Appendix A.2.

From (A13), which was first derived in this form by Shaing et al. (1994a) (see their (22)), we can deduce that particles are trapped for

(3.11)\begin{equation} \frac{(v_{{\parallel} f}+u_f)^2}{2}\leq\left.\left\lbrace S_f \left[\left(\mu B_f-v_{{\parallel} f}u_f\right)\left(\frac{B}{B_f}-1\right)+ \frac{Ze}{m}(\phi_\theta-\phi_{\theta f})\right]\right\rbrace\right\rvert_\text{max}. \end{equation}

The quantity $S$ is the squeezing factor as defined by (Hazeltine 1989)

(3.12)\begin{equation} S=1+\frac{cI^2}{B\varOmega}\frac{\partial^2\phi}{\partial \psi^2}. \end{equation}

Equation (3.11) implies that $v_{\parallel f}+u_f\sim \sqrt {|S_f|\epsilon }v_t$, which is consistent with Shaing & Hazeltine (1992). In our case, $S_f\sim 1$ and $\epsilon \ll 1$ and hence $v_{\parallel f }\simeq -u_f$ holds, to lowest order, in the trapped–barely passing region. We can rewrite (3.11) setting $v_{\parallel f}\simeq -u_f$

(3.13)\begin{equation} \frac{(v_{{\parallel} f}+u_f)^2}{2}\leq\left.\left\lbrace S_f\left[\left(\mu B_f+u_f^2\right) \left(\frac{B}{B_f}-1\right)+\frac{Ze}{m}(\phi_\theta-\phi_{\theta f})\right]\right\rbrace\right\rvert_\text{max}. \end{equation}

Now we see that the term on the right-hand side containing $u_f^2$ is the centrifugal force that pushes particles towards the outboard midplane and is small in low flow neoclassical theory. Here, both the magnetic mirror force and the centrifugal force can trap particles on the outboard side. For $\phi _c>0$, the electric potential can oppose the magnetic mirror and the centrifugal force and if the electrostatic force is strong enough, it can cause trapping of particles on the inboard side. This will become relevant in Appendix D.

Example orbits for trapped and passing particles for a circular-flux-surface tokamak are shown in figure 2. In the figure, we emphasise the difference between the width of trapped and passing particle orbits.

Figure 2. Orbits of passing (green) and trapped (red) particles which follow from (A4) and (A16a,b) are shown for $r/R=0.1$ and circular flux surfaces (blue). We chose $\theta _f=0$, $\phi _\theta =0$, $\mu B_f/v_t^2=1$, $\varOmega _f \psi _f/(Iv_t)=1$, $u_f/v_t=1.5$ and $S_f=1.5$. We use $v_{\parallel f}/v_t=-u_f/v_t+5$ for the example passing particle trajectory and $v_{\parallel f}/v_t=-u/v_t+0.2$ for the trapped particle trajectory. The spatial coordinates $X$ and $Y$ determine the position in the poloidal plane with respect to the magnetic axis. To make the orbits visible, we have chosen a flux surface with radius $r=\sqrt {X^2+Y^2}= \varOmega \psi _f/(I v_t)$, but note that we assume $r\ll \varOmega \psi _f/(I v_t)$ in the rest of the paper. The deviation from the flux surface are much larger for trapped particles than for passing particles.

4. Banana regime

The drift kinetic equation follows from an expansion of the Vlasov equation in $\rho /L$. In our case, this expansion is equivalent to an expansion in $\epsilon$ because $\rho /L\sim \rho /\rho _p\sim \epsilon$, where $\rho /R\ll \epsilon ^2$. Keeping only terms to O $(\epsilon ^3\varOmega f)$, the steady state drift kinetic equation for an ion distribution function $f(\psi,\theta,v_{\parallel},\mu )$ is

(4.1)\begin{align} & \left(v_{\parallel} \boldsymbol{\hat{b}}+\boldsymbol{v}_E\right) \boldsymbol{\cdot}\boldsymbol{\nabla}\theta\frac{\partial{f}}{\partial{\theta}}+\left(\boldsymbol{v}_E+\boldsymbol{v}_M\right) \boldsymbol{\cdot}\boldsymbol{\nabla}\psi\frac{\partial{f}}{\partial{\psi}} \nonumber\\ & \quad +\left[\boldsymbol{\hat{b}}+\frac{v_{\parallel}}{\varOmega}\boldsymbol{\hat{b}}\times \left(\boldsymbol{\hat{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\hat{b}}\right)\right]\boldsymbol{\cdot} \left(-\mu \boldsymbol{\nabla} B+\frac{Ze}{m}\boldsymbol{E}\right)\frac{\partial{f}}{\partial{v_{\parallel}}}=C[f,f]+\varSigma, \end{align}

where $\boldsymbol {v}_E$ is the $E\times B$-drift, $\boldsymbol {v}_M=\mu \boldsymbol {\hat {b}}\times \boldsymbol {\nabla } B/\varOmega +v_{\parallel} ^2 \boldsymbol {\hat {b}}\times (\boldsymbol {\hat {b}}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {\hat {b}})/\varOmega$ is the magnetic drift, $C[f,f]$ is the Fokker–Planck ion–ion collision operator and we include a source $\varSigma \sim \sqrt {\epsilon }\nu f$, which is consistent with our estimate in § 2. Note that we neglect terms small in $\epsilon$ and ion–electron collisions that are small in $\sqrt {m_e/m}$, where $m_e$ is the electron mass. It is convenient to make a change of variables from $(v_{\parallel},\psi )$ to the fixed-$\theta$ variables $(v_{\parallel f},\psi _f)$. The resulting drift kinetic equation is

(4.2)\begin{equation} \left.\dot{\theta}\frac{\partial{f}}{\partial{\theta}}\right\vert_{v_{{\parallel} f},\psi_f}=C[f,f]+\varSigma, \end{equation}

where $\dot {\theta }=(v_{\parallel} \boldsymbol {\hat {b}}+\boldsymbol {v}_E)\boldsymbol {\cdot }\boldsymbol {\nabla }\theta$, $f=f(\psi _f,\theta,v_{\parallel f},\mu )$ and the derivative in $\theta$ is holding $v_{\parallel f}$ and $\psi _f$ fixed. To lowest order in the inverse aspect ratio, one can approximate $\dot {\theta }\simeq (v_{\parallel} +u)/qR \gtrsim \epsilon ^{1/2} v_t/qR$. Assuming that the collisionality is in the banana regime $qR\nu /v_t \ll \epsilon ^{3/2}$, the system is, to lowest order in collision frequency, described by

(4.3)\begin{equation} \frac{v_{\parallel}+u}{qR}\frac{\partial{f}}{\partial{\theta}}=0 \end{equation}

and, hence, $f$ is to lowest order independent of $\theta$. Thus, any poloidal variations in density, mean flow velocity or temperature must be small.

To determine the dependence of $f$ on $\psi _f$, $v_{\parallel f}$ and $\mu$, we define the transit average, which is the average over one orbit of a particle. For passing particles, the transit average is

(4.4)\begin{equation} \langle \mathcal{F}\rangle_\tau=\frac{1}{\tau} \int_0^{2{\rm \pi}} \frac{\mathrm{d}\theta}{\dot{\theta}}\mathcal{F}, \end{equation}

where

(4.5)\begin{equation} \tau=\int_0^{2{\rm \pi}}\frac{\mathrm{d}\theta}{\dot{\theta}}. \end{equation}

Using the approximate form of $\dot {\theta }$, the transit average for trapped particles is

(4.6)\begin{equation} \langle \mathcal{F}\rangle_\tau=\frac{qR}{\tau} \int_{-\theta_b}^{\theta_b} \frac{\mathrm{d}\theta}{v_{\parallel}+u}\mathcal{F}(v_{\parallel}+u>0)+\frac{qR}{\tau} \int_{-\theta_b}^{\theta_b}\frac{\mathrm{d}\theta}{|v_{\parallel}+u|}\mathcal{F}(v_{\parallel}+u<0), \end{equation}

where

(4.7)\begin{equation} \tau=2qR \int_{-\theta_b}^{\theta_b}\frac{\mathrm{d}\theta}{|v_{\parallel}+u|} \end{equation}

and $\theta _b$ is the bounce angle, determined by $v_{\parallel} +u=0$. Transit averaging (4.2) gives

(4.8)\begin{equation} \langle C[f,f]\rangle_\tau=-\langle \varSigma \rangle_\tau. \end{equation}

To lowest order in $\epsilon$, the source $\langle \varSigma \rangle _\tau \sim \sqrt {\epsilon }\nu f$ is negligible, and the solution is a $\theta$-independent Maxwellian in fixed-$\theta$ variables,

(4.9)\begin{equation} f_{M_f}=n(\psi_f)\left(\frac{m}{2{\rm \pi} T(\psi_f)}\right)^{3/2}\exp \left(-\frac{m\left(v_{{\parallel} f }-V_{\parallel}(\psi_f)\right)^2}{2 T(\psi_f)}- \frac{m\mu B_f}{T(\psi_f)}\right). \end{equation}

Note that, unlike usual neoclassical theory, we keep the mean parallel velocity $V_{\parallel} \sim v_{\parallel}$. To zeroth order in $\epsilon$ particles do not leave their flux surface or experience a change in their parallel velocity going through one orbit, that is, $\psi \simeq \psi _f$ and $v_{\parallel} \simeq v_{\parallel f}$.

The dependence of $T$ on $\psi _f$ might be surprising because strong temperature gradients usually drive deviations away from a Maxwellian equilibrium. If the time scale associated with the ion energy flux $Q$, given by $nT/|\boldsymbol {\nabla } \psi |(\partial Q/\partial \psi )$ is longer than the ion–ion collision time, and the orbit widths are of the same order as the transport barrier, there is no temperature gradient because all particles have reached thermodynamic equilibrium and have been able to sample the entire volume. This is why the temperature gradient was assumed to be small in Kagan & Catto (2010) and Catto et al. (2013). However, by having introduced the large aspect ratio expansion, the gradient lengths can be of the same size as the poloidal gyroradius whilst still being much larger than the ion orbit width. In this way, we can get a Maxwellian to lowest order and a strong temperature gradient at the same time.

We define the next-order solution as

(4.10)\begin{equation} f=f_{M_f}+h(\psi_f, v_{{\parallel} f}, \mu)=f_M+g(\psi,\theta, v_{\parallel}, \mu), \end{equation}

where $f_M$ is the Maxwellian in (4.9) evaluated at the particle variables $\psi, v_{\parallel}$ and $\mu$,

(4.11)\begin{equation} f_{M}=n(\psi)\left(\frac{m}{2{\rm \pi} T(\psi)}\right)^{3/2}\exp \left(-\frac{m\left(v_{{\parallel} }-V_{\parallel}(\psi)\right)^2}{2 T(\psi)}- \frac{m\mu B(\psi,\theta)}{T(\psi)}\right), \end{equation}

and $h\sim g\sim \sqrt {\epsilon }f_M$ are the $O(\sqrt {\epsilon })$ corrections to the Maxwellian. One needs to be careful about the distinction between $h$ and $g$. Whilst $h$ is the distribution function in the fixed-$\theta$ variables and can be interpreted as the distribution of orbits, $g$ is a function of the variables $\psi$, $v_{\parallel}$ and $\mu$ and it is the distribution function of particles in the classic sense.

In the banana regime, the collision frequency satisfies $qR\nu /v_t\ll \epsilon ^{3/2}$. The collisionality is small enough that, in both the freely passing and the trapped–barely passing region, orbits can be completed before particles collide. Consequently, $h$ does not depend on $\theta$ to next order as $\dot {\theta }\partial h/\partial \theta \sim \epsilon ^{1/2}v_t h/qR$, while $C[h, f_M]+C[f_M,h]\sim \nu h/\epsilon$. Thus, following (4.3), $h$ does not depend on $\theta$. The large aspect ratio expansion is crucial from here on. We expand $g=h+f_{M_f}-f_M$ in orders of $\sqrt {\epsilon }$,

(4.12)\begin{equation} g=g_0+g_1+\cdots\quad \text{where}\ g_n\sim \epsilon^{{(n+1)}/{2}} f_M. \end{equation}

We will call the solution in the freely passing region, where $|v_{\parallel f}+u_f|\gg \sqrt {\epsilon }v_t$, the freely passing distribution function $g^p$, and the solution in the trapped–barely passing region, where $|v_{\parallel f}+u_f|\sim \sqrt {\epsilon }v_t$, the trapped–barely passing distribution function $g^t$. Note that, for convenience, we use the superscript $t$ for the trapped–barely passing region even though $g^t$ also includes the distribution of barely passing particles. The function $g^t$ only exists in a small region of phase space, where $|v_{\parallel f}+u_f| \sim \sqrt {\epsilon } v_t$. Thus, the contribution of $g^t$ can be interpreted as a discontinuity in $g^p$. We will find that it is sufficient to set $g\approx g^p$ in the entire phase space and determine from the solution for $g^t$ the jump and derivative discontinuity conditions at $v_{\parallel} =-u$ for $g^p$. A sketch of $g$ and how $g^t$ is reduced to a discontinuity is shown in figure 3.

Figure 3. (a) This is a sketch of the distribution function $g$. The region of trapped–barely passing particles (pink) is small whereas the passing region (white) covers most of velocity space. (b) The contribution coming from trapped–barely passing particles is approximated as a discontinuity of the passing particle distribution function and its derivatives in velocity space.

Within the trapped–barely passing region only – the region shaded in pink in figure 3(a) – we introduce the velocity variable $w\equiv v_{\parallel } + u\sim \sqrt {\epsilon } v_t$ which is defined such that, within the trapped–barely passing region, the region of overlap with the passing particle region maps to $w\rightarrow \pm \infty$, whereas from the point of view of the passing particle region, the region of overlap is still located at $v_{\parallel} +u\rightarrow 0$. The new variable $w$ effectively stretches out the trapped–barely passing region. We require that the outer limiting solutions for $g^t$ match the two inner limiting solutions of $g^p$, such that

(4.13a,b)\begin{equation} g^t(w\rightarrow\infty)=g^p(v_{\parallel}\rightarrow-u^+) \quad \text{and}\quad g^t(w\rightarrow-\infty)=g^p(v_{\parallel}\rightarrow-u^-), \end{equation}

as well as

(4.14a,b)\begin{equation} \left.\frac{\partial{g^t}}{\partial{w}}\right\rvert_{w\rightarrow\infty}= \left.\frac{\partial{g^p}}{\partial{v_{\parallel}}}\right\rvert_{v_{\parallel}\rightarrow-u^+}\quad \text{and}\quad \left.\frac{\partial{g^t}}{\partial{w}}\right\rvert_{w\rightarrow-\infty}= \left.\frac{\partial{g^p}}{\partial{v_{\parallel}}}\right\rvert_{v_{\parallel}\rightarrow-u^-}. \end{equation}

The jump condition at the trapped–passing boundary becomes

(4.15)\begin{equation} \Delta g^p=g^t_0(w\rightarrow \infty)-g^t_0(w\rightarrow - \infty). \end{equation}

The jump condition measures the difference between the co- and counter-moving barely passing particle distribution across the trapped–barely passing region.

In order for this jump to remain finite, the derivative of $g^t_0$ must tend to zero at $\pm \infty$. The discontinuity condition in the derivatives thus requires the next-order correction

(4.16)\begin{equation} \varDelta \left(\frac{\partial{g^p}}{\partial{v_{\parallel}}}\right)=\left.\frac{\partial{g^t_1}}{\partial{w}}\right\rvert_{w\rightarrow \infty}- \left.\frac{\partial{g^t_1}}{\partial{w}}\right\rvert_{w\rightarrow - \infty}. \end{equation}

The jump and derivative discontinuity conditions follow from the solution of (4.8), for which we need an expression for the ion–ion collision operator. The lowest-order solution is a Maxwellian, so we can linearise the collision operator around $f_M$ using (4.10),

(4.17)\begin{equation} C[f,f]\simeq C[f_M,g]+C[g,f_M]\equiv C^{(l)}[g] . \end{equation}

Here, we have used that the collision operator acting on the Maxwellians vanishes. We neglect the smaller, nonlinear contribution $C[g,g]$. The linearised collision operator is

(4.18)\begin{align} C^{(l)}[g] & = \lambda \boldsymbol{\nabla}_v\boldsymbol{\cdot}\left[\int\mathrm{d}^3v'f_M f'_M \boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega \omega \boldsymbol{\cdot}\left(\boldsymbol{\nabla}_v\left(\frac{g}{\,f_M}\right)-\boldsymbol{\nabla}_{v'} \left(\frac{g'}{\,f_M'}\right)\right)\right] \nonumber\\ & = \boldsymbol{\nabla}_v\boldsymbol{\cdot}\left[f_M {\boldsymbol{\mathsf{M}}} \boldsymbol{\cdot}\boldsymbol{\nabla}_v\left(\frac{g}{\,f_M}\right)-\lambda f_M \int_{V_\text{tbp}}\mathrm{d}^3v' f'_M \boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega \omega \boldsymbol{\cdot} \boldsymbol{\nabla}_{v'} \left(\frac{g^{t'}}{\,f'_M}\right)\right. \nonumber\\ & \quad \left.-\lambda f_M\int_{V_{p}}\,\mathrm{d}^3v' f'_M \boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega \omega \boldsymbol{\cdot} \boldsymbol{\nabla}_{v'} \left(\frac{g^{p'}}{\,f'_M}\right)\right], \end{align}

where $\lambda =2{\rm \pi} Z^4e^4 \log \varLambda /m^2$ and $\log \varLambda$ is the Coulomb logarithm. The integrals are over the trapped–barely passing region $V_\text {tbp}$ and the freely passing region $V_{p}$, respectively, and $\boldsymbol {\omega }=\boldsymbol {v}-\boldsymbol {v}'$. We have introduced the matrix

(4.19)\begin{align} {\boldsymbol{\mathsf{M}}} & = \lambda \int\mathrm{d}^3v'f'_M\boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega \omega \nonumber\\ & =\frac{\nu_\perp}{4}\left(|\boldsymbol{v}-V_{\parallel}\boldsymbol{\hat{b}}|^2 {\boldsymbol{\mathsf{I}}}-(\boldsymbol{v}-V_{\parallel}\boldsymbol{\hat{b}}) (\boldsymbol{v}-V_{\parallel}\boldsymbol{\hat{b}})\right)+ \frac{\nu_{\parallel}}{2}(\boldsymbol{v}-V_{\parallel}\boldsymbol{\hat{b}}) (\boldsymbol{v}-V_{\parallel}\boldsymbol{\hat{b}}), \end{align}

(4.20a–c)\begin{align} \nu_\perp & =3\sqrt{\frac{\rm \pi}{2}}\nu\frac{\varXi(x)-\varPsi(x)}{x^3},\quad \nu_{\parallel}=3\sqrt{\frac{\rm \pi}{2}}\nu\frac{\varPsi(x)}{x^3},\quad\text{and}\quad \nu= \frac{4\sqrt{\rm \pi} Z^4e^4n\log\varLambda}{3T^{3/2}m^{1/2}}, \end{align}

where $x=\sqrt {m/(2T)}|\boldsymbol {v}-V_{\parallel} \boldsymbol {\hat {b}}|$, $\varXi (x)=\text {erf}(x)=(2/\sqrt {{\rm \pi} })\int _0^x\exp (-y^2)\,\mathrm {d}y$, $\varPsi (x)= (\varXi - x\varXi ')/ (2x^2)$. The term proportional to $\nu _\perp$ describes pitch angle scattering and the term proportional to $\nu _{\parallel}$ represents energy diffusion.

We proceed to find the correction $g$. We expand (4.8) in orders of $\sqrt {\epsilon }$ and find to $O(\nu f_M/\sqrt {\epsilon })$ the jump condition $\Delta g^p$ in § 4.1 and to $O(\nu f_M)$ the derivative discontinuity condition $\varDelta (\partial g^p/\partial v_{\parallel} )$ in § 4.2. The distribution function $g^p$ as well as poloidal variations of density and potential enter at $O(\sqrt {\epsilon }\nu f_M)$ and are presented in §§ 4.3 and 4.4

4.1. Jump condition

The solution in the trapped–barely passing region gives the jump and derivative discontinuity conditions for $g^p$. We start by finding an expression for the jump condition (4.15) by collecting terms to $O(\nu f_M/\sqrt {\epsilon })$ in (4.8). The results of this subsection were already derived in a similar way by Shaing et al. (1994a). We reproduce the calculations to this order before presenting the higher-order calculations where we find significant differences with previous work.

The equation to solve for $g_0^t$ is

(4.21)\begin{equation} \langle C^{(l)}[g]\rangle_\tau=0. \end{equation}

Changing to the fixed-$\theta$ variables and keeping only terms of $O(\nu f_M/\sqrt {\epsilon })$ of the collision operator in (4.18) yields

(4.22)\begin{equation} C^{(l)}[g] \simeq \boldsymbol{\nabla}_v w_f\boldsymbol{\cdot} \frac{\partial{}}{\partial{w_f}} \left[f_M {\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v w_f\frac{\partial{(g^t/f_M)}}{\partial{w_f}} \right]. \end{equation}

Only the derivatives with respect to $w_f \equiv v_{\parallel f}+u_f$ are kept because they are larger than the other velocity derivatives by $1/\sqrt {\epsilon }$. This is because in the trapped–barely passing region $w_f\sim \sqrt {\epsilon }v_t$ and hence we assume $\partial g^t/\partial w_f\sim g^t/(\sqrt {\epsilon }v_t)$. Using fixed-$\theta$ variables is also convenient because the matching between the trapped–barely passing and freely passing region will hold for all $\theta$. It follows from (A17) that

(4.23)\begin{equation} \boldsymbol{\nabla}_v w_f=\frac{\left[w+Su\left(B_f/B-1\right)\right]\boldsymbol{\hat{b}}-S\left(B_f/B-1\right) \boldsymbol{v}_\perp}{w_f}\simeq\frac{w}{w_f}\boldsymbol{\hat{b}}. \end{equation}

Thus, the linear collision operator to lowest order is

(4.24)\begin{equation} C^{(l)}[g] \simeq \frac{w}{w_f}\frac{\partial{}}{\partial{w_f}}\left[\textsf{M}_{\parallel} \frac{w}{w_f}\frac{\partial{g^t_{0}}}{\partial{w_f}}\right], \end{equation}

where we have introduced the parallel component of $\boldsymbol{\mathsf{M}}$

(4.25)\begin{equation} \textsf{M}_{\parallel} \equiv \boldsymbol{\hat{b}}\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot} \boldsymbol{\hat{b}}\simeq\frac{\nu_\perp}{2}\mu B+\frac{\nu_{\parallel}}{2}(u+V_{{\parallel}})^2. \end{equation}

Here, we have used that $v_{\parallel} \simeq -u$ for trapped–barely passing particles. The collision frequencies $\nu _{\parallel}$ and $\nu _\perp$ are evaluated at $x\simeq \sqrt {m[(u+V_{\parallel} )^2+2\mu B]/(2T)}$.

To determine $g_0^t$, we use (4.10) and expand the lowest-order solution around a Maxwellian in the variables $(\psi,v_{\parallel},\mu )$

(4.26)\begin{align} f & = f_{M_f}+h(\psi_f,v_{{\parallel} f},\mu)\simeq f_M+(\psi_f-\psi) \left[\frac{\partial{}}{\partial{\psi}}\ln p+\frac{m(v_{\parallel}-V_{\parallel})}{T}\frac{\partial{V_{\parallel}}}{\partial{\psi}} \right. \nonumber\\ & \quad \left.+\left(\frac{m(v_{\parallel}-V_{\parallel})^2}{2T}+ \frac{m\mu B}{T}-\frac{5}{2}\right)\frac{\partial{}}{\partial{\psi}}\ln T\right]f_M \nonumber\\ & \quad -\frac{m(v_{\parallel}-V_{\parallel})}{T}(v_{{\parallel} f}-v_{\parallel})f_M+h(\psi,v_{\parallel},\mu). \end{align}

The radial derivative of the magnetic field is small and the term $m\mu /T(\partial B/\partial \psi )\sim \partial /\partial \psi \ln B\sim 1/(Ir)$ can be dropped. This result can be rewritten using the velocity variable $w=v_{\parallel} +u$, the relations (A16a,b), and the fact that $v_{\parallel} \simeq -u$ in the trapped–barely passing region

(4.27)\begin{equation} f\simeq f_M - \frac{I}{S\varOmega}\left(w-w_f\right)\mathcal{D}f_M(v_{\parallel}=-u)+h, \end{equation}

where we have defined

(4.28)\begin{align} \mathcal{D}=\frac{\partial{}}{\partial{\psi}}\ln{p} -\frac{m(u+ V_{\parallel})}{T} \left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}-\frac{\varOmega}{I}\right) + \left(\frac{m(u+V_{\parallel})^2}{2T}+\frac{m\mu B}{T}-\frac{5}{2}\right)\frac{\partial{}}{\partial{\psi}}\ln{T}. \end{align}

To avoid cluttering our notation, we will not distinguish between fixed-$\theta$ variables and $(\psi, v_{\parallel}, \mu )$ in most terms as they are almost the same. We will only keep the distinction between the two types of variables in places where they appear subtracted from each other, e.g. when we need $v_{\parallel} - v_{\parallel f}$ or $\psi - \psi _f$.

One can define the auxiliary function $\bar {h}$, which is a function of fixed-$\theta$ variables only, as

(4.29)\begin{equation} \bar{h}=h+\frac{I}{\varOmega}\frac{w_f}{S}\mathcal{D} f_{M}(v_{\parallel}=-u), \end{equation}

and with that we find

(4.30)\begin{equation} g^t_0=\bar{h}-\frac{I}{\varOmega}\frac{w}{S}\mathcal{D}f_{M}(v_{\parallel}=-u). \end{equation}

The trapped–barely passing region contains both barely passing particles and trapped particles and we need to distinguish between the two. The trapped–barely passing boundary for ions trapped on the low- (high-) field side for $S>0$ ($S<0$) and $\theta _f=0$ is

(4.31)\begin{equation} w_ \text{tpb}^2=4 S\left[\left(\mu B+u^2\right)\frac{r}{R}-\frac{Ze\phi_c }{m}\right]. \end{equation}

The trapped–barely passing boundary for ions trapped on the low- (high-) field side for $S>0$ ($S<0$) and $\theta _f={\rm \pi}$ is

(4.32)\begin{equation} w_\text{tpb}^2=4 S\left[\frac{Ze\phi_c }{m}-\left(\mu B+u^2\right)\frac{r}{R}\right]. \end{equation}

A more detailed discussion about the distinction between the two cases, is presented in Appendix D. For barely passing particles, for which $w_f^2\geq w_\text {tpb}^2$ holds, one can change from transit averages to flux surface averages by using that

(4.33)\begin{equation} \left\langle\frac{w}{w_f}(\cdots)\right\rangle_\tau= \frac{1}{\tau}\int\frac{\mathrm{d}\theta}{w}qR\frac{w}{w_f}(\cdots)= \frac{2{\rm \pi} qR}{\tau w_f}\langle\cdots\rangle_\psi \end{equation}

where $\langle \cdots \rangle _\psi =1/(2{\rm \pi} )\int \mathrm {d}\theta (\cdots )$ is the flux surface average. Then, using expression (4.30) and $(\partial w/\partial w_f)\simeq w_f/w$, the transit averaged collision operator becomes

(4.34)\begin{equation} \langle C^{(l)}[g]\rangle_\tau\simeq \frac{2{\rm \pi} qR}{\tau w_f}\frac{\partial{}}{\partial{w_f}} \left\lbrace \textsf{M}_{\parallel}\left[\frac{\langle w\rangle_\psi}{w_f}\frac{\partial{\bar{h}}}{\partial{w_f}}- \frac{I}{\varOmega S}\mathcal{D}f_{M}(v_{\parallel}=-u)\right]\right\rbrace. \end{equation}

For trapped particles, which obey $w_f^2\leq w_\text {tpb}^2$, the contribution $g^t_0-\bar {h}$ is odd in $w$ and hence it follows from (4.6) and (4.30) that

(4.35)\begin{align} \langle C^{(l)}[g-\bar{h}]\rangle_\tau & =-\left\langle\frac{w}{w_f}\frac{\partial{}}{\partial{w_f}} \left[\textsf{M}_{\parallel} \frac{I}{\varOmega S}\mathcal{D}f_{M}(v_{\parallel}=-u)\right]\right\rangle_\tau \nonumber\\ & =\frac{1}{\tau}\int_{-\theta_b}^{\theta_b}\frac{\mathrm{d}\theta}{w_f}qR\frac{\partial{}}{\partial{w_f}} \left[\textsf{M}_{\parallel}\frac{I}{\varOmega S}\mathcal{D}f_{M}(v_{\parallel}=-u)\right] \nonumber\\ & \quad -\frac{1}{\tau}\int_{-\theta_b}^{\theta_b}\frac{\mathrm{d}\theta}{w_f}qR\frac{\partial{}}{\partial{w_f}} \left[\textsf{M}_{\parallel} \frac{I}{\varOmega S}\mathcal{D}f_{M}(v_{\parallel}=-u)\right]=0. \end{align}

It then follows from (4.21) and (4.34) that

(4.36)\begin{equation} \textsf{M}_{\parallel} \frac{ \tau \langle w^2\rangle_\tau}{w_f}\frac{\partial{\bar{h}}}{\partial{w_f}}=K, \end{equation}

where $K$ is a constant. $\textsf{M}_{\parallel}$ is constant in $w_f$ and

(4.37)\begin{equation} \frac{\tau \langle w^2\rangle_\tau}{w_f}=qR\int^{\theta_b}_{-\theta_b}\mathrm{d}\theta \frac{w}{w_f}=qR\int^{\theta_b}_{-\theta_b}\mathrm{d}\theta\sqrt{1-\kappa^2\sin^2(\theta/2)}, \end{equation}

where $\kappa ^2$ is defined in (D4), such that for $w_f\rightarrow 0$, $\kappa ^2\rightarrow \infty$ as $\theta _b\rightarrow 0$. Hence, $\tau \langle w^2\rangle _\tau /w_f\rightarrow 0$ for $w_f\rightarrow 0$ and consequently $K=0$ and $\partial \bar {h}/\partial w_f=0$. For trapped particles, we find from (4.30) that

(4.38)\begin{equation} \frac{\partial{g^t_{0}}}{\partial{w_f}}=-\frac{I}{\varOmega S}\frac{w_f}{w}\mathcal{D}f_{M}(v_{\parallel}=-u). \end{equation}

The contribution $\langle C^{(l)}[g-\bar {h}]\rangle _\tau$ is not zero for barely passing particles because particles do not bounce, so there is no change in the sign of $w$ and thus the transit average of a function that is odd in $w$ does not vanish. Using (4.34) with the boundary condition $\partial g_0^t/\partial w_f\rightarrow 0$ for $w_f\rightarrow \infty$, we find that the derivative of the distribution function for barely passing particles is

(4.39)\begin{equation} \frac{\partial{g^t_{0}}}{\partial{w_f}}=\frac{I}{\varOmega S}\left(\frac{w_f }{\langle w \rangle_\psi}- \frac{w_f}{w}\right)\mathcal{D}f_{M}(v_{\parallel}=-u), \end{equation}

where we have used $\partial w/\partial w_f\simeq w_f/w$. For the jump condition (4.15) we need to integrate (4.38) and (4.39) over $w_f$. We will show in § 4.3 that in the freely passing particle region, the distribution function is independent of $\theta$ to lowest order and hence the jump condition must be independent of $\theta$ as well. Thus, the jump condition must satisfy

(4.40)\begin{equation} \Delta g^p=\int_{-\infty}^\infty\mathrm{d}w_f \frac{\partial{g_0^t}}{\partial{w_f}}= \left\langle\int_{-\infty}^\infty\mathrm{d}w_f \frac{\partial{g_0^t}}{\partial{w_f}}\right\rangle_\psi. \end{equation}

We calculate this integral in Appendix D using the potential $\phi _\theta =\phi _c\cos \theta$ (see § 4.4). The final result is

(4.41)\begin{equation} \Delta g^p=-2.758\frac{I}{ \varOmega S}\sqrt{\left|S\left[(\mu B+u^2) \frac{r}{R}-\frac{Ze}{m}\phi_c\right]\right|} \mathcal{D}f_M(v_{\parallel}=-u). \end{equation}

The distribution function $g^t$ in (B1) can be plotted using the integrals from Appendix D. The results for different values of $\theta$ are shown in figure 4. We find that the derivative is discontinuous at the trapped–passing boundary, and that the jump (4.41) is the same for any value of $\theta$.

Figure 4. The distribution function $g^t$ in the trapped and barely passing region is symmetric around $w=0$ and goes towards the same constants for any value of $\theta$ at $w\rightarrow \pm \infty$. Here, we chose $Iv_t\mathcal {D}f_M(v_{\parallel} =-u)/(\varOmega S )=1$, $g(\psi _f,w_f=0,\mu )=-1.2$, $w_\text {tpb}=\pm 1.5$ and $w_f$ is in units of thermal velocity. The jump is $\Delta g^p=-2.0685$.

4.2. Derivative discontinuity condition

We proceed to derive an expression for the discontinuity condition (4.16). For the jump condition, we have to consider terms of $O(\nu f_M/\sqrt {\epsilon })$. For the derivative discontinuity condition, we still consider the trapped–barely passing particles but need to go to higher order in $\sqrt {\epsilon }$ and collect terms of $O(\nu f_M)$. Going back to (4.21), we perform the change of variables in the collision operator (4.18) and only keep terms of $O(\nu f_M)$ or larger to get

(4.42)\begin{align} C^{(l)}[g] & \simeq \frac{1}{\mathcal{J}}\frac{\partial{}}{\partial{w_f}}\left[\mathcal{J} f_M \boldsymbol{\nabla}_v w_f\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v\left(\frac{g^t}{\,f_M}\right)\right] \nonumber\\ & \quad + \frac{1}{\mathcal{J}}\frac{\partial{}}{\partial{\mu}}\left[\mathcal{J} f_M \boldsymbol{\nabla}_v \mu \boldsymbol{\cdot} {\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v w_f\frac{\partial{(g_0^t/f_M)}}{\partial{w_f}}\right] \nonumber\\ & \quad +\frac{1}{\mathcal{J}}\frac{\partial{}}{\partial{\psi_f}}\left[\mathcal{J} f_M \boldsymbol{\nabla}_v \psi_f \boldsymbol{\cdot} {\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v w_f\frac{\partial{(g_0^t/f_M)}}{\partial{w_f}}\right], \end{align}

where

(4.43)\begin{equation} \mathcal{J}=\det\left(\frac{\partial{(\boldsymbol{r},\boldsymbol{v})}}{\partial{(\psi_f,\theta,\zeta,w_f,\mu,\varphi)}}\right)\simeq \frac{1}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \theta} \frac{1}{\boldsymbol{\nabla}_vw_f\boldsymbol{\cdot}(\boldsymbol{\nabla}_v \mu\times \boldsymbol{\nabla}_v\varphi)} \simeq\frac{qRw_f}{w} \end{equation}

is the Jacobian (note that we used (4.23) to obtain the last equality), $\varphi$ is the gyroangle with $\boldsymbol {\nabla }_v\varphi =\boldsymbol {\hat {b}}\times \boldsymbol {v}/v_\perp ^2$ and

(4.44a,b)\begin{equation} \boldsymbol{\nabla}_v\mu=\frac{\boldsymbol{v}_\perp}{B},\quad \boldsymbol{\nabla}_v\psi_f=\boldsymbol{\nabla}_v(\psi_f-\psi)\simeq \frac{I}{\varOmega S}\left(\frac{w}{w_f}-1\right)\boldsymbol{\hat{b}}, \end{equation}

for which we have used (A16a,b). The Maxwellians in the second and third terms of (4.42) are evaluated at $v_{\parallel }=-u$. Recall that the derivatives with respect to $w_f$ are bigger by $1/\sqrt {\epsilon }$ than the derivatives with respect to $\mu$ and $\psi _f$.

We argued in (4.16) that the parallel velocity derivative of $g^t_1$ is required for the derivative discontinuity condition. This derivative is of order $\sqrt {\epsilon } f_M$ and hence $g^t_1$ only appears in the first term of (4.42), where the second derivative in parallel velocity of $g^t_1$ produces a term of $O(\nu f_M)$. In all other terms that involve smaller derivatives with respect to $\mu$ and $\psi _f$, only $g^t_0$ enters to this order. We show in Appendix C that taking the transit average of the collision operator yields

(4.45)\begin{align} \langle C^{(l)}[g]\rangle_\tau & \simeq \frac{1}{w_f \tau}\frac{\partial{}}{\partial{w_f}} \left[f_Mw_f\tau\left\langle\boldsymbol{\nabla}_vw_f\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v \left(\frac{g^t}{\,f_M}\right)\right\rangle_\tau\right] \nonumber\\ & \quad +\frac{1}{w_f \tau}\frac{\partial{}}{\partial{\mu}}\left[2\mu f_Mw_f\tau \textsf{M}_{{\perp}} \left\langle\frac{w}{w_f}\frac{\partial{(g^t_0/f_M)}}{\partial{w_f}}\right\rangle_\tau\right] \nonumber\\ & \quad +\frac{1}{w_f\tau}\frac{\partial{}}{\partial{\psi_f}}\left[f_Mw_f\tau\frac{I}{\varOmega S}\textsf{M}_{\parallel} \left\langle\left(\frac{w}{w_f}-1\right)\frac{w}{w_f}\frac{\partial{(g^t_0/f_M)}}{\partial{w_f}}\right\rangle_\tau\right]=0. \end{align}

Here, we introduced the component of $\boldsymbol{\mathsf{M}}$

(4.46)\begin{equation} \textsf{M}_\perp \equiv \frac{\boldsymbol{v}_\perp}{|\boldsymbol{v}_\perp|^2}\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot} \boldsymbol{\hat{b}}\simeq(-u-V_{{\parallel}})\left(-\frac{\nu_\perp}{4}+\frac{\nu_{\parallel}}{2}\right), \end{equation}

and set $v_{\parallel} =-u$ in the arguments of $\nu _{\parallel}$ and $\nu _\perp$, which is a good approximation in the trapped–barely passing region.

The first term in (4.45) contains the derivative of $g^t_1$ that is needed for the discontinuity condition. The distribution function for trapped–barely passing particles, $g^t$, has to match with $g^p$ at the boundary between the trapped–barely passing region and the freely passing region, and thus

(4.47)\begin{equation} w_f\left\langle\boldsymbol{\nabla}_vw_f\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v\left(\frac{g^t}{\,f_M}\right)\right\rangle_\tau \simeq w\boldsymbol{\hat{b}}\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v\left(\frac{g^p}{\,f_M}\right) \end{equation}

for $w\rightarrow \pm \infty$. Hence, the solution for the discontinuity condition (4.16) in the banana regime takes the form

(4.48)\begin{align} & \varDelta\left[ w\tau\boldsymbol{\hat{b}}\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v \left(\frac{g^p}{\,f_M}\right) \right]f_M \nonumber\\ & \quad =-\frac{\partial{}}{\partial{\mu}}\left[f_M \int_{-\infty}^{\infty}\mathrm{d}w_f w_f\tau 2\mu \textsf{M}_{{\perp} } \left\langle \frac{w}{w_f}\frac{\partial{(g^t_0/f_M)}}{\partial{w_f}}\right\rangle_\tau\right] \nonumber\\ & \qquad -\frac{\partial{}}{\partial{\psi_f}}\left[f_M\frac{I}{\varOmega S} \int_{-\infty}^\infty\mathrm{d}w_f\:w_f\tau \textsf{M}_{{\parallel} } \left\langle\left(\frac{w}{w_f}-1\right)\frac{w}{w_f}\frac{\partial{(g^t_0/f_M)}}{\partial{w_f}}\right\rangle_\tau\right], \end{align}

where we have multiplied (4.45) by $w_f\tau$ and integrated over $w_f$. Note that on the left-hand side of the equation $w\tau \simeq 2{\rm \pi} qR$. Following the steps in Appendix D.2 and recalling (4.40), we arrive at

(4.49)\begin{equation} \varDelta \left[\boldsymbol{\hat{b}}\boldsymbol{\cdot}{\boldsymbol{\mathsf{M}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_v \left(\frac{g^p}{\,f_M}\right)\right]f_M=-\frac{\partial{}}{\partial{\mu}}\left(2\mu \textsf{M}_{{\perp}} \Delta g^p\right)+ \frac{\partial{}}{\partial{\psi_f}}\left(\frac{I}{\varOmega S}\textsf{M}_{{\parallel}}\Delta g^p\right), \end{equation}

where $\Delta g^p$ is given in (4.41).

We have found the jump and derivative discontinuity conditions. Next, an equation for the freely passing region is derived which completes an approximate description of the entire velocity space.

4.3. The freely passing region

The freely passing particle distribution function enters to $O(\sqrt {\epsilon }\nu f_M)$ in (4.8). The explicit expression of the collision operator in (4.18) is substituted into the simplified drift kinetic equation (4.8), which gives

(4.50)\begin{align} & \left\langle\boldsymbol{\nabla}_v\boldsymbol{\cdot}\left[f_M {\boldsymbol{\mathsf{M}}} \boldsymbol{\cdot}\boldsymbol{\nabla}_v \left(\frac{g^p}{\,f_M}\right)-\lambda f_M\int_{V_{p}}\mathrm{d}^3v' f'_M \boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega \omega \boldsymbol{\cdot} \boldsymbol{\nabla}_{v'} \left(\frac{g^{p'}}{\,f'_M}\right)\right]\right\rangle_\tau \nonumber\\ & \quad -\lambda\left\langle \boldsymbol{\nabla}_v \boldsymbol{\cdot}\left[ f_M\int_{V_\text{tbp}}\mathrm{d}^3v'f'_M \boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega\omega\boldsymbol{\cdot}\boldsymbol{\nabla}_{v'} \left(\frac{g^{t'}}{\,f'_M}\right)\right]\right\rangle_\tau=-\langle \varSigma \rangle_\tau. \end{align}

The distribution function $g\sim \sqrt {\epsilon } f_M$ and the gradient acting on $g^{t'}$ gives a factor of $1/\sqrt {\epsilon }v_t$. In the third term on the right-hand side $\boldsymbol {\nabla }_{v'}g^{t'}\sim f_M/v_t$ and $V_\text {tbp}\sim \sqrt {\epsilon }v_t^3$, so all three terms on the left-hand side are of $O(\sqrt {\epsilon }\nu f_M)$.

We combine the first two terms in (4.50) and define the linearised freely passing collision operator

(4.51)\begin{equation} C_p^{(l)}[g]\equiv\boldsymbol{\nabla}_v\boldsymbol{\cdot}\left[f_M {\boldsymbol{\mathsf{M}}} \boldsymbol{\cdot}\boldsymbol{\nabla}_v \left(\frac{g^p}{\,f_M}\right)-\lambda f_M\int_{V_{p}}\mathrm{d}^3v' f'_M \boldsymbol{\nabla}_\omega\boldsymbol{\nabla}_\omega \omega \boldsymbol{\cdot} \boldsymbol{\nabla}_{v'}\left(\frac{g^{p'}}{\,f'_M}\right)\right] \end{equation}

to write (4.50) as

(4.52)\begin{equation} \langle C_p^{(l)}[g] \rangle_\tau -\lambda \left\langle \boldsymbol{\nabla}_v \boldsymbol{\cdot}\left[ f_M\int_{V_\text{tbp}}\mathrm{d}^3v'f'_M\boldsymbol{\nabla}_\omega \boldsymbol{\nabla}_\omega\omega\boldsymbol{\cdot}\boldsymbol{\nabla}_{v'}\left(\frac{g^{t'}}{\,f'_M}\right)\right]\right\rangle_\tau=-\langle \varSigma \rangle_\tau. \end{equation}

This is the equation for the passing distribution function. Equation (4.52) has solvability conditions, which are the moment equations we calculate in § 5. To obtain the moment equations, the jump and derivative discontinuity conditions in (4.41) and (4.48) are needed.

We are interested in the poloidal variations of density, mean parallel flow velocity, temperature and electric potential, for which the $\theta$-dependent part of $g^p$, $g^p-\langle g^p \rangle _\psi$, is of interest. We argued that $h$ only depends on $\theta$ via the dependence of $\psi _f$ and $v_{\parallel f}$ on $\theta$. Since $g^p=h^p+f_{M_f}-f_M$ and $f_{M_f}-f_M\sim \epsilon f_M$, $g^p \simeq h^p$ to lowest order. The $\theta$-dependent part of $g^p$ is given by the next order

(4.53)\begin{align} & g^p-\langle g^p\rangle_\psi\simeq f_{M_f}-f_M-\langle \,f_{M_f}-f_M\rangle_\psi \nonumber\\ & \quad \simeq \left(\psi_f-\psi-\langle \psi_f-\psi\rangle_\psi\right)\frac{\partial{\,f_M}}{\partial{\psi}}+ \left(v_{{\parallel} f}-v_{\parallel}-\langle v_{{\parallel} f}-v_{\parallel}\rangle_\psi\right) \frac{\partial{\,f_M}}{\partial{v_{\parallel}}}+\frac{m \mu}{T}\left(B-\langle B\rangle_\psi\right)f_M \nonumber\\ & \quad =-\frac{Ir}{\varOmega R}\frac{\left(v_{\parallel}^2+\mu B\right) \cos{\theta}-Ze\phi_\theta R/(mr)}{v_{\parallel}+u}\left[\frac{\partial{}}{\partial{\psi}}\ln{p} + \frac{m(v_{\parallel}- V_{\parallel})}{T}\left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}-\frac{\varOmega}{I}\right) \right.\nonumber\\ & \qquad \left.+\left(\frac{m(v_{\parallel}-V_{\parallel})^2}{2T}+ \frac{m\mu B}{T}-\frac{5}{2}\right)\frac{\partial{}}{\partial{\psi}}\ln{T}\right]f_M- \frac{r}{R}\cos{\theta}\frac{m}{T}\left[v_{\parallel}(v_{\parallel}-V_{\parallel})+\mu B \right]f_M, \end{align}

where we have used the relations (A5) and (A6) as well as $B_f/B-\langle B_f/B\rangle _\psi =(r/R)\cos \theta$. The $\theta$-dependent part of the distribution function is of $O(\epsilon f_M)$ and consequently the $\theta$-independent part of $g^p$ is bigger than $g^p-\langle g^p \rangle _\psi$ by order $\sqrt {\epsilon }$. In Appendix B we show that the $\theta$ dependent part of the solution for $g^t_0$ matches with (4.53).

4.4. Poloidal variations and electric potential

In the tokamak core, trapped particles are located around $v_{\parallel} =0$, and for a Maxwellian with $V_{\parallel} =0$ the number of passing particles with $v_{\parallel} >0$ and $v_{\parallel} <0$ is the same to lowest order. The trapped–passing boundary in our ordering is shifted such that trapped particles are located around $v_{\parallel} =-u$. The lowest-order distribution function is still a Maxwellian, but it has a mean parallel velocity $V_{\parallel}$. For $V_{\parallel} \neq -u$, this implies that the number of passing particles with $v_{\parallel} +u>0$ and $v_{\parallel} +u<0$ is different. This discrepancy causes a poloidal variation in density, mean parallel velocity, temperature and poloidal potential.

If, for example, the magnetic drifts are pointing downwards, as shown in figure 5, particles with a positive (negative) poloidal velocity are being pushed inwards (outwards) with respect to their flux surface at $\theta =0$ and outwards (inwards) at $\theta ={\rm \pi}$. Let us assume a density gradient such that there is higher density inside a flux surface than there is outside. In this case, there are more particles with positive poloidal velocity at $\theta =0$ than there are particles with negative poloidal velocity (see figure 5a), because particles with positive poloidal velocity come from the high density region. At $\theta ={\rm \pi}$, the opposite is true, because the orbits of particles with positive poloidal velocity come from the low density region (see figure 5b). Thus, for a shifted trapped–passing boundary in the strong gradient case, the number of particles with positive and negative poloidal velocity are different to lowest order in $\epsilon$ and $\rho /r$ and density varies poloidally within a flux surface. For comparison, the same effect occurs in standard low flow neoclassical theory, but the number of particles with positive and negative poloidal velocity is the same to lowest order in $\rho /r$ and these effects cancel out. The asymmetry in the passing particle distribution function in the strong gradient case gives a poloidal density variation of $O(\epsilon )$, whereas, in standard low flow neoclassical theory, the poloidal density variation is much smaller. The same argument can be constructed for poloidal variation of temperature and mean parallel flow.

Figure 5. (a) At $\theta =0$, particles with a positive poloidal velocity (red) are pushed inwards, completing their orbits through a region of higher density, and particles with a negative poloidal velocity (blue) are pushed outwards, completing their orbits through a region of lower density. Hence, red particles are more numerous than blue particles. (b) At $\theta ={\rm \pi}$ on the same flux surface, the opposite is the case and there are fewer red particles than there are blue particles. If red particles are more numerous than blue particles and the density is higher at smaller radii, there will be a higher density at $\theta =0$ than at $\theta ={\rm \pi}$ and there is poloidal variation of density within a flux surface.

The small poloidal variation of density, $n_\theta$, is

(4.54)\begin{equation} n_\theta(\psi,\theta)=\int\mathrm{d}^3v g-\left\langle \int\mathrm{d}^3vg \right\rangle_\psi. \end{equation}

The integration is over the entire range of the parallel velocity and hence over both, the trapped–barely passing and freely passing regions. The freely passing region is the part of velocity space for which $v_{\parallel}$ is not close to $-u$. Importantly, the freely passing distribution function (4.53) diverges at $v_{\parallel} =-u$. This divergence is picked up by the trapped distribution function $g^t$. As a result, the integration over phase space is split up into an integration over $g^t$ in the trapped–barely passing region and a principal value integral over $g^p$ which captures the freely passing region while ignoring the divergence near $v_{\parallel} =-u$. Contribution from the divergence is accounted for by the integral of the distribution function $g^t$ in the trapped–barely passing region. For trapped particles, it follows directly from (B1) that

(4.55)\begin{equation} \int\mathrm{d}\mu\int_{\text{trapped}}\mathrm{d}w\,2{\rm \pi} B g^t-\left\langle \int\mathrm{d}\mu\int_{\text{trapped}}\mathrm{d}w\,2{\rm \pi} B g^t\right\rangle_\psi=0. \end{equation}

For barely passing and freely passing particles, the flux surface average of the density can be replaced by the integral over the flux surface averaged distribution function because the $\theta$-dependence of $B$ is small. Thus, (4.54) can be written as

(4.56)\begin{align} n_\theta=\int\mathrm{d}\mu\int_{\text{barely passing}}\mathrm{d}w\,2{\rm \pi} B (g^t-\langle g^t\rangle_\psi)+\int\mathrm{d}\mu \left[\text{PV}\!\int\mathrm{d}v_{\parallel} 2{\rm \pi} B(g^p-\langle g^p\rangle_\psi)\right], \end{align}

where the first term only contains the barely passing particles. However, this contribution vanishes to lowest order because $g^t_0-\langle g^t_0\rangle _\psi$ is odd in $w$, which follows from (B10). The integration of the second term in (4.56) is performed in Appendix E, where the $\theta$-dependent part of the distribution function is taken from (4.53). The result is

(4.57)\begin{align} & n_\theta=- n\frac{Ir}{\varOmega R}\left\lbrace \sqrt{\frac{2T}{m}}J \left[\left(\frac{mV_{\parallel}^2}{T}\cos\theta+\cos\theta -\frac{Ze\phi_\theta R}{Tr}\right) \left(\frac{\partial{}}{\partial{\psi}}\ln p-\frac{3}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right) \right.\right. \nonumber\\ & \quad \left.+\cos\theta\frac{\partial{}}{\partial{\psi}}\ln T \right]+\left[1-2\sqrt{\frac{m}{2T}}(V_{\parallel}+u)J\right] \left\lbrace(V_{\parallel}-u)\cos\theta\left(\frac{\partial{}}{\partial{\psi}}\ln p-\frac{3}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right) \vphantom{\sqrt{\frac{2T}{m}}}\right. \nonumber\\ & \quad -(V_{\parallel}+u)\left[\left(\frac{mV_{\parallel}^2}{2T}+\frac{1}{2}\right) \cos\theta-\frac{Ze\phi_\theta R}{2Tr}\right]\frac{\partial{}}{\partial{\psi}}\ln T+ \left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}-\frac{\varOmega}{I}\right) \left[\left(\frac{m u^2}{T}+1 \right.\vphantom{\sqrt{\frac{2T}{m}}}\right. \nonumber\\ & \quad \left.\vphantom{\sqrt{\frac{2T}{m}}} \left.\left.-\frac{m(V_{\parallel}+u)^2}{T}\right) \cos\theta-\frac{Ze\phi_\theta R}{Tr}\right]\right\rbrace+ \left[1+2\frac{m(V_{\parallel}+u)^2}{2T}-4\left(\frac{m}{2T}\right)^{3/2}(V_{\parallel}+u)^3J\right] \nonumber\\ & \quad \left.\vphantom{\sqrt{\frac{2T}{m}}} \times\cos\theta\left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}-\frac{\varOmega}{I}+ \frac{V_{\parallel}-u}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right)\right\rbrace-2n\frac{r}{R}\cos\theta, \end{align}

where we introduced the function

(4.58)\begin{equation} J=\frac{\sqrt{\rm \pi} }{2}\exp\left(-\frac{m(u+V_{\parallel})^2}{2T}\right) \operatorname{erfi}\left(\sqrt{\frac{m}{2T}}(u+V_{\parallel})\right), \end{equation}

which is plotted in figure 6 and $\operatorname {erfi}(x)=(2/\sqrt {{\rm \pi} })\int _0^x\exp (t^2)\,\mathrm {d}t$. The orbit width of passing particles is of order $\epsilon$ and hence the poloidal variation in density is of order $\epsilon$ as well.

Figure 6. The function $J$ defined in (4.58) as a function of $\bar {y}=\sqrt {m/(2T)}(u+V_{\parallel} )$.

The poloidal variation in density creates a poloidal variation in electric potential $\phi _\theta$ that is determined via quasineutrality. Assuming a Boltzmann response of the electrons, the quasineutrality condition yields

(4.59)\begin{equation} Z\int\mathrm{d}^3v\,g-\left\langle Z\int\mathrm{d}^3v\,g \right\rangle_\psi=\frac{en_e}{T_e}\phi_\theta. \end{equation}

Looking at (4.57) we find that the potential has the form $\phi _\theta =\phi _c \cos \theta$, and the quasineutrality condition (4.59) yields

(4.60)\begin{align} & \left\lbrace\frac{en_e}{T_e}-\frac{Z^2ne I}{T\varOmega} \left[\sqrt{\frac{2T}{m}}J\left(\frac{\partial{}}{\partial{\psi}}\ln p-\frac{3}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right)+ \left[1-2\sqrt{\frac{m}{2T}}(V_{\parallel}+u)J\right]\left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}- \frac{\varOmega}{I} \right.\right.\vphantom{\sqrt{\frac{2T}{m}}}\right.\nonumber\\ & \quad \left.\vphantom{\sqrt{\frac{2T}{m}}}\left.\vphantom{\sqrt{\frac{2T}{m}}}\left. -\frac{(V_{\parallel}+u)}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right)\right]\right\rbrace\phi_c=-Zn \frac{Ir}{\varOmega R}\left\lbrace \sqrt{\frac{2T}{m}}J \left[\left(\frac{mV_{\parallel}^2}{T}+1\right) \right.\right.\nonumber\\ & \quad \times\left.\vphantom{\sqrt{\frac{2T}{m}}} \left(\frac{\partial{}}{\partial{\psi}}\ln p-\frac{3}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right) +\frac{\partial{}}{\partial{\psi}}\ln T\right]+ \left[1-2\sqrt{\frac{m}{2T}}(V_{\parallel}+u)J\right] \nonumber\\ & \quad \times\left[(V_{\parallel}-u)\left(\frac{\partial{}}{\partial{\psi}}\ln p-\frac{3}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right) \right. \nonumber\\ & \quad \left.+\left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}-\frac{\varOmega}{I}\right) \left(\frac{mu^2}{T}+1-\frac{m(V_{\parallel}+u)^2 }{T}\right)-\frac{V_{\parallel}+u}{2} \left(\frac{mV_{\parallel}^2}{T}+1\right)\frac{\partial{}}{\partial{\psi}}\ln T\right] \nonumber\\ & \left.\vphantom{\sqrt{\frac{2T}{m}}} +\left[1+2\frac{m}{2T}(V_{\parallel}+u)^2-4\left(\frac{m}{2T}\right)^{3/2}(V_{\parallel}+u)^3J\right] \left(\frac{\partial{V_{\parallel}}}{\partial{\psi}}-\frac{\varOmega}{I}+ \frac{V_{\parallel}-u}{2}\frac{\partial{}}{\partial{\psi}}\ln T\right)\right\rbrace \nonumber\\ & \quad -2Zn\frac{r}{R}. \end{align}

For $\phi _c>0$, the maximum of the potential is on the low-field side of the plasma, so the potential can trap particles on the high-field side for $S>0$. For $\phi _c<0$ and $S>0$, the potential reaches its maximum on the high-field side and it can trap particles on the low-field side if electrostatic trapping dominates over magnetic trapping and centrifugal force.

Charge exchange recombination spectroscopy measurements in both Alcator C-Mod (Theiler et al. 2014; Churchill et al. 2015) and ASDEX-Upgrade (Cruz-Zabala et al. 2022) have observed poloidal variation in impurity density and temperatures in the pedestal of H-mode plasmas. These experiments also demonstrated that the main ion temperature and radial electric field cannot simultaneously be flux functions. This is consistent with our calculation and argumentation of poloidal variation in the electric potential and density.

We have found expressions for the distribution function in the passing region and the jump and derivative discontinuity condition given by the trapped–barely passing region, and we have found the form of the poloidally varying component of the electric potential. These expressions are needed to calculate the solvability conditions for (4.52).

5. Moment equations

In order to study the transport in the pedestal, we want to find particle, parallel momentum and energy fluxes and how they give rise to profiles of $n$, $T$, $u$, $V_{\parallel}$ and $\phi _c$. First, we integrate (4.52), for which the jump and derivative discontinuity conditions are required, and find the solvability conditions, which are the equations for particle, parallel momentum and energy conservation.

The full derivation is explained in Appendix F, where we show that the particle conservation equation

(5.1)\begin{equation} \frac{\partial{}}{\partial{\psi_f}}\left(-\frac{I }{\varOmega m}F_{{\parallel} }\right)=\int \mathrm{d}^3v_f\langle \varSigma\rangle_\tau, \end{equation}

is the result of integrating (4.52) over velocity space. Here,

(5.2)\begin{equation} F_{{\parallel} }=-\int\mathrm{d}\mu \frac{2{\rm \pi} m B}{S}\textsf{M}_{{\parallel} } \Delta g^p \end{equation}

is the parallel force due to the friction between passing and trapped particles, and $\Delta g^p$ is the jump condition given in (4.41). The integration over $\mathrm {d}^3v_f$ is an integration over velocity space in the fixed-$\theta$ variables, where $\mathrm {d}^3v_f=2{\rm \pi} B \,\mathrm {d}\mu \,\text {d}v_{\parallel f}$. The integration eliminated the contribution from the freely passing particle distribution $g^p$ to the particle transport and for this reason is a solvability condition: it must be satisfied regardless of the value of $g^p$. Trapped and barely passing particles dominate transport as we have estimated in § 2. We can compare (5.1) with a typical continuity equation

(5.3)\begin{equation} \frac{\partial{\varGamma}}{\partial{\psi_f}}=\left\langle \int\mathrm{d}^3v_f \,\varSigma\right\rangle_\psi, \end{equation}

where the term on the left-hand side is the divergence in $\psi _f$ of a particle flux $\varGamma$ and the term on the right-hand side is a source of particles. It follows directly from (5.1) and (5.3) that the neoclassical ion particle flux is

(5.4)\begin{equation} \varGamma=-\frac{I}{m\varOmega}F_{{\parallel}}. \end{equation}

The parallel force $F_{\parallel}$ can drive a radial particle flux via an effect similar to the one that gives the Ware pinch (Ware 1970).

The parallel momentum equation is the result of multiplying (4.52) by $mv_{\parallel f}$ and integrating over velocity space. The equation becomes

(5.5)\begin{equation} \frac{\partial{}}{\partial{\psi_f}}\left(\frac{I}{\varOmega}u F_{{\parallel} }\right)+F_{{\parallel} }=\gamma, \end{equation}

where $\gamma =\int \mathrm {d}^3v \,mv_{\parallel f} \langle \varSigma \rangle _\tau$ is the parallel momentum input per unit volume. The calculation that leads to (5.5) is presented in Appendix F.2. We can use the particle flux (5.4) in (5.5) and arrive at

(5.6)\begin{equation} \frac{\partial{}}{\partial{\psi_f}}\left(mu\varGamma\right)+\frac{m\varOmega}{I}\varGamma=-\gamma, \end{equation}

which is a relation purely between the particle flux, parallel momentum input and $u$. The first term on the left-hand side of (5.6) is the flux of parallel momentum carried by the trapped particles. The second term on the left is the force due to the friction between trapped and passing particles. The term on the right-hand side of the equation is a source of parallel momentum.

As for the particle and parallel momentum equations, one can find the energy equation by multiplying (4.52) by $mv_f^2/2$ and integrating over velocity space to arrive at

(5.7)\begin{equation} \frac{\partial{}}{\partial{\psi_f}}\left(\frac{IT }{m\varOmega}\varTheta\right)-u F_{{\parallel}}= \int \mathrm{d}^3v_f \frac{mv_f^2}{2} \langle\varSigma\rangle_\tau, \end{equation}

where

(5.8)\begin{equation} \varTheta=\int\mathrm{d}\mu \frac{2{\rm \pi} mB}{S}\left(\frac{m\mu B}{T}+ \frac{mu^2}{2T}\right)\textsf{M}_{{\parallel}} \Delta g^p. \end{equation}

The energy flux $Q$ is defined similarly to $\varGamma$ as

(5.9)\begin{equation} \frac{\partial{Q}}{\partial{\psi_f}}-Ze\varGamma \frac{\partial{\phi}}{\partial{\psi_f}}= \left\langle\int\mathrm{d}^3v_f \frac{m v_f^2}{2}\varSigma\right\rangle_\psi. \end{equation}

A comparison with (5.7) gives

(5.10)\begin{equation} Q=\frac{TI}{m\varOmega}\varTheta. \end{equation}

The flux of energy on the left-hand side of (5.7) contains both convective energy flux, which is the energy carried by the particle flux, and a conduction energy flux. The second term on the left of (5.9) is the work done by the radial electric field. The term on the right represents energy injection.

The same equations for particle, parallel momentum and energy (5.1), (5.5) and (5.7) can be found using moments of the original Fokker–Planck kinetic equation. At this point we can switch from fixed-$\theta$ variables to normal variables and drop the subscript $f$ because the difference is small in $\epsilon$.

We can substitute (4.25) for $\textsf{M}_{\parallel}$ and the jump condition (4.41) into (5.2) to find the particle flux from (5.4)

(5.11)\begin{align} \varGamma & =-2.758\frac{I^2{\rm \pi} B }{|S|^{3/2}\varOmega^2}n \left(\frac{m}{2{\rm \pi} T}\right)^{3/2} \sqrt{\frac{T}{m}}\int\mathrm{d}\mu \exp\left(-\frac{m(u+V_{{\parallel} })^2}{2T}-\frac{m\mu B}{T}\right) \nonumber\\ & \quad \times \sqrt{\left|\left(\frac{m\mu B}{T}+\frac{m u^2}{T}\right) \frac{r}{R}-\frac{Ze}{T}\phi_c\right|} \left(\nu_\perp\mu B+\nu_{\parallel}(u+V_{{\parallel} })^2\right)\mathcal{D}. \end{align}

Integration over $x=\sqrt {m\mu B/T+m(u+V_{\parallel })^2/(2T)}$ gives the final form of $\varGamma$

(5.12)\begin{align} \varGamma & =-1.102\sqrt{\frac{r}{R}}\frac{\nu I^2 p}{|S|^{3/2}m\varOmega^2} \left\lbrace \left[\frac{\partial{}}{\partial{\psi}}\ln p-\frac{m(u+V_{{\parallel}})}{T} \left(\frac{\partial{V_{{\parallel} }}}{\partial{\psi}}-\frac{\varOmega}{I}\right)\right]G_1(\bar{y},\bar{z}) \right. \nonumber\\ & \quad \left. -1.17\frac{1}{T}\frac{\partial{T}}{\partial{\psi}} G_2(\bar{y},\bar{z})\right\rbrace, \end{align}

where $\bar {y}=\sqrt {m/(2T)}(u+V_{\parallel })$, $\bar {z}=m{u}^2/T-Ze\phi _c R/(T r)$,

(5.13)\begin{gather} G_1(\bar{y},\bar{z})\equiv\frac{\displaystyle\int_{|\bar{y}|}^\infty\mathrm{d}x\,k(x,\bar{y},\bar{z})}{ \displaystyle\int_0^\infty\mathrm{d}x\,0.5x\,{\rm e}^{-x^2}\left[\varXi(x)-\varPsi(x)\right]}=7.51 \int_{|\bar{y}|}^\infty\mathrm{d}x\,k(x,\bar{y},\bar{z}), \end{gather}

(5.14)\begin{gather} G_2(\bar{y},\bar{z})\equiv\frac{\displaystyle\int_{|\bar{y}|}^\infty\mathrm{d}x \left(x^2-5/2\right)k(x,\bar{y},\bar{z})}{\displaystyle\int_0^\infty\mathrm{d}x\,0.5x \left(x^2-5/2\right) {\rm e}^{-x^2}\left[\varXi(x)-\varPsi(x)\right]} \nonumber\\ =-6.40\int_{|\bar{y}|}^\infty\mathrm{d}x \left(x^2-\frac{5}{2}\right)k(x,\bar{y},\bar{z}), \end{gather}

and

(5.15)\begin{equation} k(x,\bar{y},\bar{z})=\sqrt{|x^2+\bar{z}-\bar{y}^2|}{\rm e}^{-x^2} \left\lbrace\left(\frac{1}{2}-\frac{\bar{y}^2}{2x^2}\right) \left[\varXi(x)-\varPsi(x)\right]+\frac{\bar{y}^2}{x^2}\varPsi(x)\right\rbrace. \end{equation}

The functions $G_1$ and $G_2$ are normalised to recover the standard neoclassical results when $\bar {y}=0=\bar {z}$, $G_1(0,0)=1=G_2(0,0)$. The neoclassical ion particle flux in (5.12) depends on the radial electric field through $u$ (see (3.9)) and thus also through $\bar {y}$ and $\bar {z}$. We note that the term in (5.12) proportional to $[V_{\parallel} -(-u)]\varOmega /I$ is the particle flux due to the parallel friction between trapped particles located around $-u$ and the passing particles with a mean velocity $V_{\parallel}$. The term proportional to $(u+V_{\parallel} )\partial V_{\parallel} /(\partial \psi )$ is related to a shift in the Maxwellian and hence to the density gradient if the Maxwellian is not centred around the trapped region, i.e. if $V_{\parallel} +u$ is not small (see figure 7). The remaining terms include the pressure and temperature gradients that usually drive radial particle flux but here are modified by the integrals $G_1$ and $G_2$. Note that also the poloidal potential affects transport as it enters in $\bar {z}$.

Figure 7. (a) A small shift in $V_{\parallel}$ for $V_{\parallel}$ not close to $-u$ going from one surface (solid line) to another flux surface (dashed line) causes a strong change of the number of trapped particles (red area between curves) in the trapped–barely passing region (pink). (b) A small shift in $V_{\parallel}$ for $V_{\parallel}$ close to $-u$ gives only a small change in the number of trapped–barely passing particles (red areas between curves cancel) in the trapped–barely passing region.

Similarly, $Q$ is

(5.16)\begin{align} Q& =\frac{mu^2}{2}\varGamma-1.463\sqrt{\frac{r}{R}} \frac{\nu I^2 p T}{|S|^{3/2}m\varOmega^2} \left\lbrace\left[\frac{\partial{}}{\partial{\psi}}\ln p-\frac{m(u+V_{{\parallel} })}{T} \left(\frac{\partial{V_{{\parallel} }}}{\partial{\psi}}-\frac{\varOmega}{I}\right)\right]H_1(\bar{y},\bar{z})\right. \nonumber\\ & \quad \left.-0.25\frac{1}{T}\frac{\partial{T}}{\partial{\psi}} H_2(\bar{y},\bar{z})\right\rbrace, \end{align}

where

(5.17)\begin{equation} H_1(\bar{y},\bar{z})\equiv\frac{\displaystyle\int_{|\bar{y}|}^\infty\mathrm{d}x \left(x^2-\bar{y}^2\right)k(x,\bar{y},\bar{z})}{\displaystyle\int_0^\infty\mathrm{d}x\,0.5x^3 {\rm e}^{-x^2}[\varXi(x)-\varPsi(x)]}=5.66\int_{|\bar{y}|}^\infty\mathrm{d}x \left(x^2-\bar{y}^2\right)k(x,\bar{y},\bar{z}), \end{equation}

and

(5.18)\begin{align} H_2(\bar{y},\bar{z})& \equiv\frac{\displaystyle\int_{|\bar{y}|}^\infty\mathrm{d}x \left(x^2-\bar{y}^2\right)\left(x^2-5/2\right)k(x,\bar{y},\bar{z})}{\displaystyle\int_0^\infty\mathrm{d}x\,0.5x^3 \left(x^2-5/2\right){\rm e}^{-x^2}[\varXi(x)-\varPsi(x)]} \nonumber\\ & =-22.63 \int_{|\bar{y}|}^\infty\mathrm{d}x \left(x^2-\bar{y}^2\right)\left(x^2-\frac{5}{2}\right)k(x,\bar{y},\bar{z}). \end{align}

Again, we introduce a convenient normalisation such that $H_1(0,0)=1=H_2(0,0)$ in the standard neoclassical limit. The dependence of the neoclassical ion energy flux (5.16) on the radial electric field is hidden in $u$, $\bar {y}$ and $\bar {z}$.

We have found explicit expressions for particle (5.1), parallel momentum (5.5) and energy conservation (5.7). Next, we want to compare our results with previous work. First, we take the high flow and low flow neoclassical limit, and then we give a comparison of our results with those by Catto et al. (2013) and Shaing & Hsu (2012).

In the high flow regime of the usual neoclassical theory (Hinton & Wong 1985), $V_{\parallel} +u$ and all gradients as well as source terms are small. If we take this limit in (5.6) and assume that the source of parallel momentum $\gamma$ is small, we find that

(5.19)\begin{equation} \varGamma=0, \end{equation}

which is consistent with the usual result in the high flow regime (Hinton & Wong 1985; Catto, Bernstein & Tessarotto 1987). Using the particle flux equation (5.12), $\bar {\varGamma }=0$ gives

(5.20)\begin{equation} \frac{\partial{}}{\partial{\psi}}\ln p+\frac{m\varOmega(V_{\parallel}+u)}{IT}=1.17 \frac{G_2(\bar{y},\bar{z})}{G_1(\bar{y},\bar{z})}\frac{\partial{}}{\partial{\psi}}\ln T. \end{equation}

We can use this in (5.16) to get the high flow energy flux

(5.21)\begin{equation} Q=-1.71\sqrt{\frac{r}{R}}\frac{I^2\nu T p}{m\varOmega^2}\Delta \bar{Q}\frac{\partial{}}{\partial{\psi}}\ln T, \end{equation}

where

(5.22)\begin{equation} \Delta\bar{Q}\equiv\frac{H_1(\bar{y},\bar{z} )G_2(\bar{y},\bar{z})-0.21H_2 (\bar{y},\bar{z} )G_1(\bar{y},\bar{z} )}{G_1(\bar{y},\bar{z})}\geq0. \end{equation}

The quantity $\Delta \bar {Q}$ is positive, which follows from

(5.23)\begin{equation} \Delta \bar{Q}=-4.82\frac{\displaystyle\left(\int_{|\bar{y} |}^\infty\mathrm{d}x \,x^2 k(x,\bar{y},\bar{z})\right)^2- \int^\infty_{|\bar{y}|}\mathrm{d}x\,k(x,\bar{y},\bar{z}) \int_{|\bar{y}|}^\infty\mathrm{d}x\,x^4k(x,\bar{y},\bar{z})}{ \displaystyle\int_{|\bar{y}|}^\infty\mathrm{d}x\,k(x,\bar{y},\bar{z})} \end{equation}

and the Cauchy–Schwarz inequality

(5.24)\begin{equation} \left(\int_{|\bar{y}|}^\infty\mathrm{d}x\,x^2 k(x,\bar{y},\bar{z})\right)^2\leq \int^\infty_{|\bar{y}|}\mathrm{d}x\,k(x,\bar{y},\bar{z}) \int_{|\bar{y} |}^\infty\mathrm{d}x\,x^4k(x,\bar{y},\bar{z}). \end{equation}

Here, $k(x,\bar {y} ,\bar {z} )$ is given in (5.15). Note that $k>0$ because $\varXi -\varPsi >0$, $\varPsi >0$ and $x\geq |\bar {y}|$.

The quasineutrality condition (4.60) gives the poloidally varying electric potential in the high flow limit,

(5.25)\begin{equation} \left(\frac{en_e}{T_e}+\frac{Z^2n_i e}{T}\right)\phi_c=-Zn_i\frac{r}{R}\frac{mu^2}{T}. \end{equation}

The only contribution to the potential comes from the centrifugal force as all gradients and $m(V_{\parallel} +u)^2/T$ terms are small while $V_{\parallel} \simeq -u \sim v_t$.

The low flow neoclassical results can be retrieved by taking the limit of small radial electric field, $u/v_t \ll 1$, small mean parallel flow, $V_{\parallel} /v_t\ll 1$, and small gradients. It follows from (5.25) that the poloidal variation of the potential is small so that we can set $\bar {z}=0$ in the arguments of $G_1$, $G_2$, $H_1$ and $H_2$. Without a source of parallel momentum $\gamma =0$, (5.6) gives $\varGamma =0$, so the mean parallel flow follows directly from (5.20)

(5.26)\begin{equation} V_{{\parallel} }=-\frac{IT}{m\varOmega}\left(\frac{\partial{}}{\partial{\psi}}\ln p+ \frac{Ze}{T}\frac{\partial{\varPhi}}{\partial{\psi}}-1.17\frac{\partial{}}{\partial{\psi}}\ln T\right). \end{equation}

The neoclassical energy flux $Q$ then follows directly from (5.21) for $\Delta \bar {Q}=0.79$ and reads

(5.27)\begin{equation} Q=-1.35\sqrt{\frac{r}{R}}\frac{I^2\nu pT}{m \varOmega ^2}\frac{\partial{}}{\partial{\psi }}\ln T, \end{equation}

in agreement with Hinton & Wong (1985) and Catto et al. (1987).

We can compare our results with those of Catto et al. (2013) by taking the limit of small temperature gradient and small ${V}_{\parallel}$. We are able to retrieve the same energy flux if we set $\phi _\theta =0$, $\varGamma =0$ and correct an error in Catto et al. (2011) and pointed out by Shaing & Hsu (2012). The calculation is presented in detail in Appendix G.1.

The energy flux $Q$ in (5.16) is proportional to $|S|^{-3/2}$ and decays $\sim \exp (-\bar {y}^2)$ which is consistent with the results of strong radial electric field and radial electric field shear obtained by Shaing & Hazeltine (1992) and Shaing & Hsu (2012). We compare our results in the limit $(I/\varOmega )(\partial V_{\parallel} /\partial \psi ) \ll 1$ and $\varGamma =0$ with those of Shaing & Hsu (2012) in Appendix G.2. We find the same particle and energy equations if we account for a discrepancy in the function $k(x,\bar {y} ,\bar {z} )$.

Comparisons with numerical results can be made in certain limits. The global code PERFECT requires weak temperature gradients and could be checked against our results in the limit of small temperature gradient (Landreman et al. 2014). Other codes such as the axisymmetric versions of XGC (Chang, Ku & Weitzner 2004), Gkeyll (Hakim et al. 2020) and COGENT (Dorf et al. 2012) could be used to reproduce aspects of the strong gradient fluxes and poloidal variation. In the next section, we impose radial force balance (see (6.1)) which needs to be reconsidered carefully when comparing the following results with numerical evaluations of fluxes.

6. Transport equations and flux conditions

We work with (5.6), (5.12) and (5.16) to find relations between the particle flux $\varGamma$, the parallel momentum input $\gamma$, the energy flux $Q$ and the physical quantities $T$, $n$, $u$, $V_{\parallel}$ and $\phi _c$. Given $\varGamma$ and $Q$ as functions of $\psi$, and boundary conditions at the top or bottom of the transport barrier, we can integrate the equations to obtain the profiles of $T$, $n$, $u$, $V_{\parallel}$ and $\phi _c$.

So far, we have an equation for the particle flux (5.12), the parallel momentum equation (5.6), the energy flux (5.16) and quasineutrality (4.60). We are missing an equation for the radial electric field to be able to relate $\varGamma$, $\gamma$ and $Q$ with $T$, $n$, $u$, $V_{\parallel}$ and $\phi _c$. The equation for the radial electric field is provided by the conservation of toroidal angular momentum, but the necessary derivation is beyond the scope of this paper. For the purpose of the following calculations, we assume that for the ions, the pressure gradient is the dominant contribution in the radial force balance (Kagan & Catto 2008; McDermott et al. 2009; Viezzer et al. 2013). Hence, we impose

(6.1)\begin{equation} Z e n \frac{\partial{\varPhi}}{\partial{\psi}}+\frac{\partial{p}}{\partial{\psi}}=0, \end{equation}

which can be written as

(6.2)\begin{equation} \frac{\partial{}}{\partial{\psi}}\ln n=-\frac{\varOmega m u}{IT}-\frac{\partial{}}{\partial{\psi}}\ln T. \end{equation}

We introduce the new, dimensionless quantities

(6.3a–e)\begin{gather} \bar{u}=\sqrt{\frac{m}{2T_0}}u,\quad \bar{V}=\sqrt{\frac{m}{2T_0}}V_{{\parallel}},\quad \bar{T}=\frac{T}{T_{0}},\quad \bar{n}=\frac{n}{n_{0}},\quad \bar{\phi}_c=\frac{Ze\phi_cR}{T_0r} \end{gather}

(6.4a–c)\begin{gather}\frac{\partial{}}{\partial{\bar{\psi}}}=\frac{I}{\varOmega}\sqrt{\frac{2T_0}{m}}\frac{\partial{}}{\partial{\psi}},\quad \bar{z}=\bar{T}^{-1}\left(2\bar{u}^2-\bar{\phi}_c\right),\quad \bar{y}=\frac{\bar{u}+\bar{V}}{\sqrt{\bar{T}}}, \end{gather}

where $T_0$ is the ion temperature and $n_0$ the ion density at the boundary $\psi =0$. In the banana regime, the normalised fluxes are

(6.5a–c)\begin{equation} \bar{\varGamma}=\frac{\varGamma}{n_0I\sqrt{\dfrac{2T_0r}{mR}}\dfrac{\nu_0}{\varOmega}},\quad \bar{Q}=\frac{Q}{n_0I\sqrt{\dfrac{2T_0r}{mR}}T_0\dfrac{\nu_0}{\varOmega}},\quad \bar{\gamma}=\frac{\gamma}{\nu_0 n_0\sqrt{2mT_0r/R}}, \end{equation}

where $\nu _0$ is the collision frequency at the boundary. Changing to these dimensionless variables, we arrive at the following set of equations for the banana regime: The particle flux equation from (5.12) and (6.2) is

(6.6)\begin{align} \bar{\varGamma}=-0.55\frac{\bar{n}^2}{|S|^{3/2}\bar{T}^{3/2}} \left\lbrace\left[-2\bar{u}-2(\bar{u}+\bar{V}) \left(\frac{\partial{\bar{V}}}{\partial{\bar{\psi}}}-1\right)\right]G_1(\bar{y},\bar{z}) -1.17\frac{\partial{\bar{T}}}{\partial{\bar{\psi}}}G_2(\bar{y},\bar{z})\right\rbrace. \end{align}

The parallel momentum equation from (5.6) is

(6.7)\begin{equation} \bar{u} \frac{\partial{\bar{\varGamma}}}{\partial{\bar{\psi}}}+S\bar{\varGamma}=-\bar{\gamma}. \end{equation}

The energy flux equation from (5.12), (5.16) and (6.2) is

(6.8)\begin{align} \bar{Q}& = \bar{u}^2\bar{\varGamma}-0.73\frac{\bar{n}^2}{|S|^{3/2}\bar{T}^{1/2}} \left\lbrace\left[-2\bar{u}-2(\bar{u}+\bar{V})\left(\frac{\partial{\bar{V}}}{\partial{\bar{\psi}}}-1\right)\right] H_1(\bar{y},\bar{z}) \right. \nonumber\\ & \quad \left.-0.25\frac{\partial{\bar{T}}}{\partial{\bar{\psi}}}H_2(\bar{y},\bar{z})\right\rbrace. \end{align}

The pressure balance equation from (6.2) gives

(6.9)\begin{equation} \frac{\partial{}}{\partial{\bar{\psi}}}\ln \bar{n}=-2\frac{\bar{u}}{\bar{T}}-\frac{\partial{}}{\partial{\bar{\psi}}}\ln \bar{T}, \end{equation}

and the equation for the potential, which can be derived from (4.60), is

(6.10)\begin{align} & \left\lbrace 1-\frac{Z T_{e}}{T_{0}}\frac{1}{\bar{T}} \left[\sqrt{\bar{T} }J \left(-2\frac{\bar{u} }{\bar{T} }-\frac{3}{2}\frac{\partial{}}{\partial{\bar{\psi} }}\ln \bar{T}\right) \right.\right.\nonumber\\ & \qquad \left.\left.+\left(1-2\frac{\bar{V} +\bar{u}}{\sqrt{\bar{T}}}J\right) \left(\frac{\partial{\bar{V} }}{\partial{\bar{\psi} }}-1-\frac{\bar{V} +\bar{u} }{2} \frac{\partial{}}{\partial{\bar{\psi}}}\ln T\right)\right]\right\rbrace\bar{z} \nonumber\\ & \quad =\frac{Z T_{e}}{T_{0}}\frac{1}{\bar{T} }\left\lbrace \sqrt{\bar{T} }J \left[\left(2\frac{\bar{V} ^2-\bar{u} ^2}{\bar{T} }+1\right) \left(-2\frac{\bar{u} }{\bar{T} }-\frac{3}{2}\frac{\partial{}}{\partial{\bar{\psi} }} \ln \bar{T} \right)+ \frac{\partial{}}{\partial{\bar{\psi} }}\ln \bar{T}\right] \right.\nonumber\\ & \qquad +\left[1-2\frac{\bar{V} +\bar{u} }{\sqrt{\bar{T} }}J\right]\left[(\bar{V} -\bar{u} ) \left(-2\frac{\bar{u} }{\bar{T} }-\frac{3}{2}\frac{\partial{}}{\partial{\bar{\psi} }} \ln \bar{T} \right) \right.\nonumber\\ & \qquad \left.+\left(\frac{\partial{\bar{V} }}{\partial{\bar{\psi} }}-1\right) \left(1-2\frac{(\bar{V} +\bar{u} )^2}{\bar{T} }\right)- (\bar{V} +\bar{u} )\left(\frac{\bar{V} ^2-\bar{u} ^2}{\bar{T} }+ \frac{1}{2}\right)\frac{\partial{}}{\partial{\bar{\psi} }}\ln\bar{T}\right] \nonumber\\ & \qquad \left.+\left[1+2\frac{(\bar{V} +\bar{u} )^2}{\bar{T} }-4 \frac{(\bar{V} +\bar{u} )^3}{\bar{T} ^{3/2}}J\right] \left(\frac{\partial{\bar{V} }}{\partial{\bar{\psi} }}-1+\frac{\bar{V} -\bar{u} }{2}\frac{\partial{}}{\partial{\bar{\psi} }}\ln \bar{T} \right)+2+2 \frac{T_{0}}{ZT_e}\bar{u} ^2 \right\rbrace. \end{align}

The functions $J$, $G_1$, $G_2$, $H_1$ and $H_2$ are given in (4.58), (5.13), (5.14), (5.17) and (5.18). This set of equations is the most important result of our calculation and allow a discussion of the neoclassical transport of ions in strong gradient regions.

We can integrate equations (6.6)–(6.10) relating $\bar {n}$, $\bar {T}$, $\bar {u}$, $\bar {V}$ and $\bar {z}$ numerically by imposing boundary conditions at the top of the transport barrier and specifying particle, parallel momentum and energy sources to find profiles in the pedestal. We discuss the implications for particle (§ 6.1) and energy flux (§ 6.2) before presenting some example profiles (§ 6.3).

6.1. Particle flux and parallel momentum injection

In order to understand the appearance of a neoclassical particle flux, we analyse the parallel momentum equation (6.7). In edge transport barriers, measurements of the radial electric field have shown that in the pedestal $\partial \phi /\partial \psi >0$ and thus $\bar {u} >0$ (McDermott et al. 2009). We assumed that in the pedestal $\bar {u} \sim 1$. However, at the boundary to the large turbulent transport region, where our model connects to the usual neoclassical regime of small gradients in density and temperature, $\bar {u} \ll 1$. Thus, we are looking for solutions with a growing positive $\bar {u}$ as one moves into the transport barrier. Importantly, if there is no parallel momentum input, the particle flux must decay to ensure that $\bar {u}$ grows, because it follows from (6.7) that

(6.11)\begin{equation} \bar{\varGamma} \propto \exp\left(-\int\mathrm{d}\bar{\psi}\frac{S }{\bar{u}}\right). \end{equation}

For $\bar {u} >0$ and $S >0$, the neoclassical particle flux $\bar {\varGamma }$ decreases and is even smaller inside a transport barrier than outside when $\bar {\gamma }=0$.

We argued in § 2 that at the inner edge of a transport barrier there is a region of large turbulent transport and small collisional transport whereas in the transport barrier, we find a region of low turbulence. In order to keep up the same total flux, the neoclassical fluxes must increase and pick up the decreasing turbulent fluxes (see figure 1). However, this initial picture is too simple as it disagrees with our analysis of the decreasing particle flux. One option to solve the contradiction is that the particle flux is still carried by turbulence because the neoclassical fluxes never pick up the turbulent contribution. There must be enough turbulence in the transport barrier to carry the entire particle flux – recall that at this point we are only discussing the particle flux and not the energy flux. So even if the entire particle flux is carried by turbulence, the energy flux could still be neoclassical (see figure 8a). The second option is that the particle flux is truly neoclassical in the transport barrier, but turbulence or impurities supply the necessary parallel momentum source $\gamma$ so that (6.11) is not valid. Somehow, and we cannot specify at this point how exactly, turbulence or impurities interact with neoclassical transport and appear as a source of parallel momentum (see figure 8b). The difference between the two options is that in the first picture, the neoclassical particle flux is close to zero whereas in the second picture the particle flux is in large part neoclassical because turbulence or impurities produce $\bar {\gamma }$. This picture is consistent with previous results by Landreman & Ernst (2012) about the necessity of sources for non-zero steady state transport in the edge. Without the source, no ion neoclassical particle flux develops in the pedestal.

Figure 8. (a) The entire particle flux is carried by turbulence and the neoclassical particle flux stays negligible. (b) Turbulence interacts with neoclassical physics and supplies a parallel momentum source that allows a growing neoclassical particle flux.

The neoclassical ion particle flux is larger than the electron particle flux by $O(\sqrt {m/m_e})$. Unlike in the weak gradient region, where intrinsic ambipolarity prevents different sizes of electron and ion particle fluxes, the neoclassical ion particle flux in the strong gradient region can be significantly larger than the neoclassical electron particle flux in the presence of sources if the total particle fluxes which include both the turbulent and neoclassical parts obey ambipolarity. Intrinsic ambipolarity (Sugama & Horton 1998; Parra & Catto 2009; Calvo & Parra 2012) does not hold in the strong gradient limit where gradient length scales are of the order of $\rho _p$.

It is also worth pointing out that $\bar {\varGamma }$ and $\bar {\gamma }$ are necessarily of opposite sign if $|\bar {\varGamma }|$ grows as one moves into the transport barrier. An outwards neoclassical ion particle flux requires a negative parallel momentum injection.

6.2. Energy flux

Next, we want to discuss the energy flux equation (6.8). In transport barriers, $\bar {T}$ and $\bar {n}$ decrease. One can use this behaviour to estimate the energy flux in this case. Combining (6.6) and (6.8) to solve for $\partial \bar {T} /\partial \bar {\psi }$ as a function of $\bar {Q}$ and $\bar {\varGamma }$ yields

(6.12)\begin{equation} \frac{\partial{\bar{T} }}{\partial{\bar{\psi} }}=-\underbrace{1.17\frac{|S| ^{3/2}\bar{T} ^{1/2}}{\bar{n} ^2}}_{{>}0} \underbrace{\frac{1}{\Delta\bar{Q}}}_{{\geq} 0} \left\lbrace\bar{Q}-\left[\bar{u} ^2+1.33\bar{T} \frac{H_1(\bar{y} ,\bar{z} )}{G_1(\bar{y} ,\bar{z} )}\right]\bar{\varGamma}\right\rbrace, \end{equation}

where $\Delta \bar {Q}$ was defined in (5.22). Figure 9 shows $\Delta \bar {Q}$ for different values of $\bar {y}$ and $\bar {z}$. It is large for small $\bar {y}$, symmetric in $\bar {y}$ with a maximum at $\bar {y}=0$, and asymmetric in $\bar {z}$ with larger values for $\bar {z}>0$. When $\bar {z}$ increases so does the number of trapped particles. Thus, $\Delta \bar {Q}$ is large when there are many trapped particles.

Figure 9. (a) The quantity $\Delta \bar {Q}$ in (5.22) as a function of $\bar {y}$ for different values of $\bar {z}$. (b) The quantity $\Delta \bar {Q}$ as a function of $\bar {z}$ for different values of $\bar {y}$.

In order to get a negative temperature gradient, the expression in braces in (6.12) must be positive. Thus, we find a lower bound for the energy flux

(6.13)\begin{equation} \bar{Q}\geq\bar{Q}_\text{min}=\left[\bar{u} ^2+1.33\bar{T} \frac{H_1(\bar{y},\bar{z} )}{G_1(\bar{y},\bar{z} )}\right]\bar{\varGamma}. \end{equation}

The factor multiplying $\bar {\varGamma }$ is positive because $\bar {u}^2\geq 0$, $\bar {T}\geq 0$, and $k>0$ and $x\geq |\bar {y}|$ in (5.13) and (5.17). From this, we see that it is not possible to only have neoclassical particle flux and zero neoclassical energy flux. As long as there is neoclassical particle flux, energy will get advected by that particle flux. Thus the energy flux will be in the same direction as the particle flux. The quantity $\bar {Q}_\text {min}$ is shown in figure 10. It is large for large $|\bar {y}|$ and small $\bar {z}$.

Figure 10. (a) The quantity $\bar {Q}_\text {min}$ defined in (6.13) as a function of $\bar {y}$ for different values of $\bar {z}$, where $\bar {u}=0$, $\bar {T}=1$ and $\bar {\varGamma }=1$ (b) The quantity $\bar {Q}_\text {min}$ as a function of $\bar {z}$ for different values of $\bar {y}$, where $\bar {u}=0$, $\bar {T}=1$ and $\bar {\varGamma }=1$.

Surprisingly, a negative density gradient imposes an upper boundary for $\bar {Q}$. From (6.9) it follows that for $\partial \bar {n} /\partial \bar {\psi } <0$

(6.14)\begin{equation} \frac{2\bar{n} }{\bar{T} }\bar{u}>1.17 \frac{|S| ^{3/2}}{\bar{n} \bar{T}^{1/2} } \frac{\bar{Q}-\bar{Q}_\text{min}}{\Delta \bar{Q}}, \end{equation}

and thus we find that in order for $\bar {T}$ and $\bar {n}$ to decay simultaneously, the neoclassical energy flux has to be

(6.15)\begin{equation} \bar{Q}_\text{min}<\bar{Q} <\bar{Q}_\text{min}+1.71 \frac{\bar{n}^2 \bar{u} }{\bar{T} ^{1/2}|S| ^{3/2}}\Delta \bar{Q}. \end{equation}

For zero neoclassical particle flux, the maximum energy flux for decaying density and temperature profiles, is

(6.16)\begin{equation} \bar{Q}_\text{max}=1.71\frac{\bar{n}^2 \bar{u} }{\bar{T} ^{1/2}|S| ^{3/2}}\Delta \bar{Q}. \end{equation}

If the density falls off faster than the temperature in such a way that $\bar {n} ^2/\sqrt {\bar {T} }\rightarrow 0$, which can be expressed as

(6.17)\begin{equation} L_{\bar{n}}<4L_{\bar{T}}, \end{equation}

then the upper bound of the energy flux in (6.15) also decreases unless it is compensated by a stronger growth in $\bar {u}\Delta \bar {Q}/|S|^{3/2}$. In most H-mode pedestals, (6.17) is observed (Viezzer et al. 2016, 2018). It follows, that in order to achieve a growing neoclassical energy flux, it is necessary that $\bar {u}\Delta \bar {Q}/|S|^{3/2}$ increases. Thus, the radial electric field seems to play an important role for the neoclassical energy flux at the top of transport barriers. Note, however, that the result in (6.15) relies strongly on the assumption made in (6.1) between the pressure gradient and the electric field, which is only applicable in the pedestal and not self-consistently derived. A more thorough discussion of this relation will be necessary and we leave it for future work. For now, using (6.2), the estimate (6.15) holds. We already argued in § 6.1 that $\bar {u}$ is positive and growing at the transition from core to pedestal. The quantity $\Delta \bar {Q}$ is large for large $|\bar {z}|$ and small $\bar {y}$ (see figure 9). This is consistent because large $|\bar {z}|$ leads to an increased number of trapped particles. Transport is dominated by trapped particles, so more trapped particles allow for a larger energy flux. Small $\bar {y}$ likewise maximises the number of trapped particles because the trapped region is located close to the maximum of the lowest-order Maxwellian.

In I-mode pedestals, the temperature falls off much faster than the density (Walk et al. 2014). In this case, (6.17) would not necessarily hold and the neoclassical heat flux could grow with a weaker radial electric field than in H-mode.

6.3. Example profiles

To show some example solutions of (6.6)–(6.10), we can take profiles of ion and electron temperature and density loosely based on those measured by Viezzer et al. (2016). With these profiles, we calculate fluxes, velocities and electric potential.

The integration of the mean parallel flow turns out to be very sensitive to the boundary conditions and source terms. Thus, we leave the discussion of the mean parallel flow solutions for future work, and instead only consider cases of known mean parallel flow. The two profiles we discuss for $\bar {V}$ are the ‘high flow’ case and the ‘low flow’ case. Here, ‘high flow’ and ‘low flow’ only refers to the relationship between the mean parallel flow and the gradients of the density, temperature and potential and not to the usual stricter limits that we have discussed at the end of § 5.

For the ‘high flow’ profile, we set

(6.18)\begin{equation} \bar{V}=-\bar{u}. \end{equation}

In this case, there is no friction between trapped and passing particles and the particle flux due to a shift in the Maxwellian is small because $f_M$ is centred around the trapped particle region (see discussion below (5.15)). For the ‘low flow’ profile, we replace condition (6.18) with the usual neoclassical solution (5.26).

The profile of $\bar {u}$ follows directly from assumption (6.9) and consequently $\bar {V}$ is given by (6.18) for the first case or (5.26) for the second case. The quantities $\bar {T}$, $\bar {n}$, $\bar {V}$, and $\bar {u}$ based on realistic profiles or assumptions are presented in figure 11. The input profiles are further discussed in Appendix H.

Figure 11. Input profiles of ion temperature, electron temperature $\bar {T}_e=T_e/T_{0}$ and density based on the profiles reported by Viezzer et al. (2016), as well as the corresponding $\bar {u}$ and $\bar {V}$. The red profile for $\bar {V}$ is the usual neoclassical result for the mean parallel velocity as given by (5.26) and the blue curve is the ‘high flow’ profile as given by (6.18). Vertical dashed lines indicate the position of the top of the pedestal $\bar {\psi }=0.8$ and the point of maximum pressure gradient and minimum radial electric field $\bar {\psi }=0.965$.

The graphs in figures 11 and 12 show the transition between core and pedestal nicely in the sense that at $\bar {\psi } =0$, which corresponds to $\rho _\text {pol}=0.8$ in Viezzer et al. (2016), the profiles of density and temperature are still relatively flat. We see the expected growth of $\bar {u}$ in the strong gradient region starting at $\bar {\psi }=0.8$ (first dashed line in figure 12) which relaxes when the pressure gradient reduces again beyond the dashed line at $\bar {\psi }=0.965$. For $\bar {V}$, we see the difference between ‘high flow’ and standard ‘low flow’ neoclassical theory. The solution for $\bar {V}$ from (6.18) exceeds the standard ‘low flow’ neoclassical result in the pedestal by about a factor of two but becomes as small as the standard ‘low flow’ neoclassical result at the boundary to the core.

Figure 12. Calculated fluxes and poloidally varying potential from the profiles in figure 11. The blue profiles are the solutions with condition (6.18) whereas the red profiles show the solution with the usual neoclassical parallel velocity (5.26). The yellow energy flux is the usual neoclassical result (5.27). Vertical dashed lines highlight the top of the pedestal $\bar {\psi }=0.8$ and the point of maximum pressure gradient and minimum radial electric field $\bar {\psi }=0.965$.

Equation (6.10) gives $\bar {z}$, with which $\bar {\varGamma }$ can be calculated from (6.6). Then, the energy flux can be calculated using (6.8). Lastly, the parallel momentum input that is necessary to sustain the particle flux follows from (6.7). The four graphs for $\bar {\varGamma }$, $\bar {Q}$, $\bar {\gamma }$, and $\bar {\phi }_c$ are presented in figure 12.

The poloidally varying part of the potential is much stronger and changes sign in the pedestal region for $\bar {V}=-\bar {u}$. The neoclassical particle flux, which is close to zero in the core requires parallel momentum input to grow. In the case with condition (6.18), the particle flux and the parallel momentum input are much bigger than for the case with the usual neoclassical mean parallel velocity (5.26). Note that, even for the ‘low flow’ neoclassical mean parallel velocity, the parallel momentum input and the particle flux are non-zero. Interestingly, the neoclassical particle flux and parallel momentum source in the pedestal for (5.26) are of opposite sign to the case with condition (6.18). The energy flux of the ‘high flow’ case matches the standard ‘low flow’ neoclassical result close to the inner boundary but further into the pedestal it grows faster with radius. In the case where we set the parallel velocity to be (5.26), the energy flux is smaller than the usual neoclassical result $Q_\text {neo}$ of (5.27). The prefactor $\bar {n} ^2/\sqrt {\bar {T}}$ in (6.8) decays in the strong gradient region for the example profiles of density and temperature, so (6.17) is satisfied, and the energy flux decays after $\bar {u}$ has reached its maximum. This is consistent with our discussion in § 6.2 and the observation that the energy transport in pedestals reaches significant neoclassical levels only in the middle of a pedestal and not at the top and bottom (Viezzer et al. 2020). If instead we had chosen profiles with a stronger temperature gradient such that $L_{\bar {n}}>4L_{\bar {T}}$, we could have been able to retrieve a growing energy flux throughout the pedestal.

In figure 13 we show the energy fluxes and their respective lower bound (6.13) and upper bound (6.16). In both cases, the energy flux is close to the upper bound in the flat gradient region. The lower bound stays close to zero where the particle flux is small.

Figure 13. Energy flux with upper and lower bounds (6.15) in (a) the ‘high flow’ case, and (b) the ‘low flow’ case.

7. Conclusions

The core is a region of strong turbulent transport. With the transition into a transport barrier such as the pedestal, turbulence gets quenched and we argue that, in order to keep up the total flux, the neoclassical fluxes must increase. This assumption is supported by experiments such as the ones by Viezzer et al. (2018), where it was demonstrated that the heat diffusivity reaches neoclassical levels in the pedestal. This opens the possibility of interaction between turbulent and neoclassical transport which we account for by keeping a source term that represents external particle, momentum and energy injection as well as interaction with turbulence. A random walk estimate was performed to predict the size of this source and to show that trapped particles give the main contribution to particle and energy transport.

We have extended neoclassical theory to transport barriers by choosing gradients to be of the same size as the poloidal gyroradius and expanded in large aspect ratio and low collisionality. A new set of variables, the fixed-$\theta$ variables, were derived from conserved quantities and confirmed that particles are trapped for $v_{\parallel} +u\sim \sqrt {\epsilon }v_t$.

A change of variables to fixed-$\theta$ variables allowed for a convenient reduction of the drift kinetic equation, to which a Maxwellian is the solution to lowest order. We have discussed the trapped–barely passing and freely passing regions in the banana regime. The drift kinetic equation can be solved for the trapped, barely passing and freely passing regions by expanding in $\sqrt {\epsilon }$. The phase space region of trapped and barely passing particles is very narrow for large aspect ratio tokamaks and can be treated as a discontinuity in the freely passing region. The only information needed from the trapped–barely passing region is the jump (4.41) and derivative discontinuity condition (4.48). Additionally, one can find expressions for the poloidal variations of density (4.57) and potential (4.60) which have been observed previously (Theiler et al. 2014; Churchill et al. 2015; Cruz-Zabala et al. 2022). Particles can get trapped on the high-field side because the poloidally varying part of the potential can oppose the magnetic mirror and centrifugal forces. When integrating over velocity space, it is necessary to keep track of particles trapped on either side.

One can take moments of the freely passing particle equation (4.52) using the jump and derivative discontinuity condition to find the particle, parallel momentum and energy conservation equations (5.1), (5.5) and (5.7). From these equations, one can identify the neoclassical particle flux (5.4) and neoclassical energy flux (5.10). We find that the poloidally varying potential affects neoclassical fluxes and that the transport is dominated by trapped particles, which have a parallel velocity close to $-u$. The fluxes match with the usual neoclassical results in the appropriate limits. They equally match with the results for strong density and electric potential gradients derived by Catto et al. (2011) after we account for the missing orbit squeezing factor in the energy flux calculation, which was previously pointed out by Shaing & Hsu (2012). In the limit of small mean velocity gradient and zero poloidal potential, we identify a previously noted discrepancy with Shaing & Hsu (2012), but are otherwise able to reproduce their results.

The parallel momentum equation proves that a parallel momentum source is required to get a non-zero neoclassical particle flux. When there is no external parallel momentum source or sink in the edge (such as impurities or neutral beam injection), this implies that either turbulence does not decay and carries the particle flux throughout the transport barrier or that there is a mechanism by which turbulence supplies parallel momentum to neoclassical transport and the particle flux is indeed partially neoclassical.

For the energy flux, we provided upper and lower bounds in relation to the particle flux to ensure decaying profiles of temperature and density (see (6.15)). The maximum energy flux can be achieved for $\bar {V}+\bar {u}=0$ and large $\bar {z}$. We also found that in pedestals a radially growing radial electric field is needed to obtain a radially growing neoclassical energy flux that substitutes the decreasing turbulent energy flux.

We compared the high flow case $\bar {V}=-\bar {u}$ with the standard low flow neoclassical mean parallel velocity (5.26) to find fluxes for the realistic profiles of temperature and density presented in figure 11, which are similar to those measured by Viezzer et al. (2016). We showed that for $\bar {V}=-\bar {u}$ the non-zero neoclassical particle flux, the energy flux, the mean parallel flow, and the poloidal variation exceed the usual neoclassical values in the strong gradient region. The neoclassical energy flux and especially the neoclassical particle flux are significantly smaller in the low flow case, but non-zero.