## 1. Introduction

The pedestal, and transport barriers in general, play an important role in tokamak performance (Wagner *et al.* Reference Wagner, Fussmann, Grave, Keilhacker, Kornherr, Lackner, McCormick, Müller, Stäbler and Becker1984; Greenfield *et al.* Reference Greenfield, Schissel, Stallard, Lazarus, Navratil, Burrell, Casper, DeBoo, Doyle and Fonck1997) 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.* Reference Viezzer, Pütterich, Conway, Dux, Happel, Fuchs, McDermott, Ryter, Sieglin and Suttrop2013). Moreover, it has been found that the ion energy transport in pedestals is close to the neoclassical level (Viezzer *et al.* Reference Viezzer, Cavedon, Fable, Laggner, McDermott, Galdon-Quiroga, Dunne, Kappatou, Angioni and Cano-Megias2018). Measurements of H-mode pedestals in Alcator C-Mod (Theiler *et al.* Reference Theiler, Churchill, Lipschultz, Landreman, Ernst, Hughes, Catto, Parra, Hutchinson and Reinke2014; Churchill *et al.* Reference Churchill, Theiler, Lipschultz, Hutchinson, Reinke, Whyte, Hughes, Catto, Landreman and Ernst2015) and Asdex-Upgrade (Cruz-Zabala *et al.* Reference Cruz-Zabala, Viezzer, Plank, McDermott, Cavedon, Fable, Dux, Cano-Megias, Pütterich and van Vuuren2022) 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 Reference Hinton and Hazeltine1976) such as the one by Viezzer *et al.* (Reference Viezzer, Cavedon, Fable, Laggner, McDermott, Galdon-Quiroga, Dunne, Kappatou, Angioni and Cano-Megias2018) 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 Reference Shaing, Hsu and Dominguez1994*a*), and causes an exponential decrease proportional to the radial electric field in plasma viscosity (Shaing *et al.* Reference Shaing, Hsu and Dominguez1994*a*) and radial heat flux (Kagan & Catto Reference Kagan and Catto2010; Shaing & Hsu Reference Shaing and Hsu2012). The mean parallel flow is also affected by a strong radial electric field and can change direction (Kagan & Catto Reference Kagan and Catto2010). 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 Reference Shaing and Hazeltine1992; Shaing, Hsu & Hazeltine Reference Shaing, Hsu and Hazeltine1994*b*).

Combining all these effects, Shaing & Hsu (Reference Shaing and Hsu2012) 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 (Reference Kagan and Catto2008) and Catto *et al.* (Reference Catto, Parra, Kagan, Parker, Pusztai and Landreman2013) have likewise developed extensions to neoclassical theory to allow for stronger density gradients to calculate fluxes. In Kagan & Catto (Reference Kagan and Catto2010) and Catto *et al.* (Reference Catto, Kagan, Landreman and Pusztai2011, Reference Catto, Parra, Kagan, Parker, Pusztai and Landreman2013), 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.* (Reference Landreman, Parra, Catto, Ernst and Pusztai2014), 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.* Reference Viezzer, Cavedon, Fable, Laggner, McDermott, Galdon-Quiroga, Dunne, Kappatou, Angioni and Cano-Megias2018). 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.

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

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

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

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

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

and hence we include a source

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 Reference Sugama and Horton1998; Parra & Catto Reference Parra and Catto2009; Calvo & Parra Reference Calvo and Parra2012) 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.* (Reference Landreman, Parra, Catto, Ernst and Pusztai2014). 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$,

*a–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

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

and

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

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

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

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

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

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

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.* (Reference Shaing, Hsu and Dominguez1994*a*). 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 Reference Shaing and Hazeltine1992). 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.* (Reference Shaing, Hsu and Dominguez1994*a*) (see their (22)), we can deduce that particles are trapped for

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

Equation (3.11) implies that $v_{\parallel f}+u_f\sim \sqrt {|S_f|\epsilon }v_t$, which is consistent with Shaing & Hazeltine (Reference Shaing and Hazeltine1992). 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$

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.

## 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

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

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

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

where

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

where

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

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,

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 (Reference Kagan and Catto2010) and Catto *et al.* (Reference Catto, Parra, Kagan, Parker, Pusztai and Landreman2013). 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

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

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 }$,

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.

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

*a,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

*a,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

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

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),

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

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

*a–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.* (Reference Shaing, Hsu and Dominguez1994*a*). 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

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

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

Thus, the linear collision operator to lowest order is

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

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 )$

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 (A16*a*,*b*), and the fact that $v_{\parallel} \simeq -u$ in the trapped–barely passing region

where we have defined

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

and with that we find

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

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

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

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

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

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

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

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

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

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

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

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$.

### 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

where

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

*a,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 (A16*a*,*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

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

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

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

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

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

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

to write (4.50) as

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

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 5*a*), 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 5*b*). 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.

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

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

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

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

where we introduced the function

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.

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

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

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.* Reference Theiler, Churchill, Lipschultz, Landreman, Ernst, Hughes, Catto, Parra, Hutchinson and Reinke2014; Churchill *et al.* Reference Churchill, Theiler, Lipschultz, Hutchinson, Reinke, Whyte, Hughes, Catto, Landreman and Ernst2015) and ASDEX-Upgrade (Cruz-Zabala *et al.* Reference Cruz-Zabala, Viezzer, Plank, McDermott, Cavedon, Fable, Dux, Cano-Megias, Pütterich and van Vuuren2022) 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

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

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

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

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

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

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

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

where

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

A comparison with (5.7) gives

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)

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

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

and

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}$.

Similarly, $Q$ is

where

and

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.* (Reference Catto, Parra, Kagan, Parker, Pusztai and Landreman2013) and Shaing & Hsu (Reference Shaing and Hsu2012).

In the high flow regime of the usual neoclassical theory (Hinton & Wong Reference Hinton and Wong1985), $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

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

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

where

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

and the Cauchy–Schwarz inequality

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,

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)

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

in agreement with Hinton & Wong (Reference Hinton and Wong1985) and Catto *et al.* (Reference Catto, Bernstein and Tessarotto1987).

We can compare our results with those of Catto *et al.* (Reference Catto, Parra, Kagan, Parker, Pusztai and Landreman2013) 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.* (Reference Catto, Kagan, Landreman and Pusztai2011) and pointed out by Shaing & Hsu (Reference Shaing and Hsu2012). 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 (Reference Shaing and Hazeltine1992) and Shaing & Hsu (Reference Shaing and Hsu2012). We compare our results in the limit $(I/\varOmega )(\partial V_{\parallel} /\partial \psi ) \ll 1$ and $\varGamma =0$ with those of Shaing & Hsu (Reference Shaing and Hsu2012) 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.* Reference Landreman, Parra, Catto, Ernst and Pusztai2014). Other codes such as the axisymmetric versions of XGC (Chang, Ku & Weitzner Reference Chang, Ku and Weitzner2004), Gkeyll (Hakim *et al.* Reference Hakim, Mandell, Bernard, Francisquez, Hammett and Shi2020) and COGENT (Dorf *et al.* Reference Dorf, Cohen, Compton, Dorr, Rognlien, Angus, Krasheninnikov, Colella, Martin and McCorquodale2012) 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 Reference Kagan and Catto2008; McDermott *et al.* Reference McDermott, Lipschultz, Hughes, Catto, Hubbard, Hutchinson, Granetz, Greenwald, LaBombard and Marr2009; Viezzer *et al.* Reference Viezzer, Pütterich, Conway, Dux, Happel, Fuchs, McDermott, Ryter, Sieglin and Suttrop2013). Hence, we impose

which can be written as

We introduce the new, dimensionless quantities

*a–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}

*a–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

*a–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

The parallel momentum equation from (5.6) is

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

The pressure balance equation from (6.2) gives

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

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.* Reference McDermott, Lipschultz, Hughes, Catto, Hubbard, Hutchinson, Granetz, Greenwald, LaBombard and Marr2009). 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

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 8*a*). 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 8*b*). 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 (Reference Landreman and Ernst2012) 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.

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 Reference Sugama and Horton1998; Parra & Catto Reference Parra and Catto2009; Calvo & Parra Reference Calvo and Parra2012) 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

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.

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

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}$.

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$

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

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

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

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.* Reference Viezzer, Fable, Cavedon, Angioni, Dux, Laggner, Bernert, Burckhart, McDermott and Pütterich2016, Reference Viezzer, Cavedon, Fable, Laggner, McDermott, Galdon-Quiroga, Dunne, Kappatou, Angioni and Cano-Megias2018). 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.* Reference Walk, Hughes, Hubbard, Terry, Whyte, White, Baek, Reinke, Theiler and Churchill2014). 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.* (Reference Viezzer, Fable, Cavedon, Angioni, Dux, Laggner, Bernert, Burckhart, McDermott and Pütterich2016). 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

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.

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.* (Reference Viezzer, Fable, Cavedon, Angioni, Dux, Laggner, Bernert, Burckhart, McDermott and Pütterich2016), 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.

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.* Reference Viezzer, Cavedon, Cano-Megias, Fable, Wolfrum, Cruz-Zabala, David, Dux, Fischer and Harrer2020). 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.

## 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.* (Reference Viezzer, Cavedon, Fable, Laggner, McDermott, Galdon-Quiroga, Dunne, Kappatou, Angioni and Cano-Megias2018), 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.* Reference Theiler, Churchill, Lipschultz, Landreman, Ernst, Hughes, Catto, Parra, Hutchinson and Reinke2014; Churchill *et al.* Reference Churchill, Theiler, Lipschultz, Hutchinson, Reinke, Whyte, Hughes, Catto, Landreman and Ernst2015; Cruz-Zabala *et al.* Reference Cruz-Zabala, Viezzer, Plank, McDermott, Cavedon, Fable, Dux, Cano-Megias, Pütterich and van Vuuren2022). 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