2 Drift kinetic equation

We begin by dividing the distribution function of each species into even and odd parts in $v_{\Vert }=\unicode[STIX]{x1D70E}|v_{\Vert }|$ ,

(2.1) $$\begin{eqnarray}\displaystyle & f_{a}(\boldsymbol{r},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})=f_{a}^{+}(\boldsymbol{r},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})+f_{a}^{-}(\boldsymbol{r},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707}), & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & f_{a}^{\pm }(\boldsymbol{r},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})={\textstyle \frac{1}{2}}[f_{a}(\boldsymbol{r},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},+1)\pm f_{a}(\boldsymbol{r},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},-1)], & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D716}=m_{a}v^{2}/2+e_{a}\unicode[STIX]{x1D6F7}(\boldsymbol{r})$ denotes the (kinetic $+$  potential) energy and $\unicode[STIX]{x1D707}=m_{a}v_{\bot }^{2}/(2B)$ the magnetic moment. The drift kinetic equation is thus split into the two equations

(2.3) $$\begin{eqnarray}\displaystyle & v_{\Vert }\unicode[STIX]{x1D735}_{\Vert }\,f^{-}=C^{+}(f)-\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & v_{\Vert }\unicode[STIX]{x1D735}_{\Vert }\,f^{+}=C^{-}(f)-\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{-}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{v}_{d}=\boldsymbol{v}_{E}+\boldsymbol{v}_{M}$ is the sum of the $\boldsymbol{E}\times \boldsymbol{B}$ and magnetic drift velocities and $C^{\pm }(f)$ denotes the even and odd parts of the collision operator $C(f)$ . Thus, for the operator (1.1)

(2.5) $$\begin{eqnarray}\displaystyle & C^{+}(f_{1})=\unicode[STIX]{x1D708}{\mathcal{L}}(f_{1}^{+}), & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle C^{-}(f_{1})=\unicode[STIX]{x1D708}\left({\mathcal{L}}(f_{1}^{-})+\frac{mv_{\Vert }}{T}wf_{0}\right). & \displaystyle\end{eqnarray}$$

Here we have omitted the species index, and shall often do so in the following. It is assumed that the electrostatic potential only depends on $\unicode[STIX]{x1D713}$ in leading order – an assumption that will be discussed below. We can thus set $\unicode[STIX]{x1D6F7}=0$ on the particular flux surface we are considering, but must retain a non-zero derivative $\unicode[STIX]{x1D6F7}^{\prime }(\unicode[STIX]{x1D713})$ . Since $\unicode[STIX]{x1D6F7}=0$ , we can use the velocity-space coordinate $v$ instead of $\unicode[STIX]{x1D716}$ . As usual, particle orbits can be divided into trapped and untrapped (or circulating) ones depending on whether the parameter $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D716}$ does, or does not, exceed $1/B_{max}$ , where $B_{max}(\unicode[STIX]{x1D713})$ denotes the maximum value of the magnetic field strength on the surface. An important difference between tokamaks and stellarators is that in the latter the field strength generally reaches $B_{max}$ at a discrete set of points on each surface, whereas $B=B_{max}$ along a circle on the inboard side of a tokamak.

The orbit average of an arbitrary function $g(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})$ is defined by

(2.7) $$\begin{eqnarray}\overline{g}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})=\frac{1}{\unicode[STIX]{x1D70F}_{b}}\int _{l_{1}}^{l_{2}}g^{+}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})\frac{\text{d}l}{|v_{\Vert }|},\end{eqnarray}$$

in the trapped region of phase space, $\unicode[STIX]{x1D706}B_{max}>1$ . Here $l_{1}$ and $l_{2}$ denote consecutive bounce points, $B(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},l_{1})=B(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},l_{2})=1/\unicode[STIX]{x1D706}$ , and

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{b}=\int _{l_{1}}^{l_{2}}\frac{\text{d}l}{|v_{\Vert }|}\end{eqnarray}$$

the bounce time. For circulating particles, $\unicode[STIX]{x1D706}>1/B_{max}$ , the orbit average is defined as

(2.9) $$\begin{eqnarray}\overline{g}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})=\lim _{L\rightarrow \infty }\int _{0}^{L}g(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})\left.\frac{\text{d}l}{v_{\Vert }}\right/\!\!\int _{0}^{L}\frac{\text{d}l}{v_{\Vert }}.\end{eqnarray}$$

Since the integral extends along a field line many times around the torus, the result is independent of $\unicode[STIX]{x1D6FC}$ . It can also be written in terms of the flux-surface average $\langle \cdots \rangle$ as (Helander Reference Helander2014)

(2.10) $$\begin{eqnarray}\overline{g}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})=\left.\left\langle \frac{Bg}{v_{\Vert }}\right\rangle \right/\!\left\langle \frac{B}{v_{\Vert }}\right\rangle .\end{eqnarray}$$

Note that the ‘orbit average’ thus defined is a time average over the trajectories that would result if the drift velocity were neglected in the guiding-centre equations of motion. The actual orbits are more complicated since particles can drift in and out of local trapping wells.

The reason for introducing the orbit average, as done here, is that this operation annihilates the first term in (2.3), which becomes

(2.11) $$\begin{eqnarray}\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+}}=\overline{C^{+}(f)}.\end{eqnarray}$$

The flow velocity of each species within the flux surface is determined by the odd part of the distribution function, which is given by the remainder of (2.3),

(2.12) $$\begin{eqnarray}f^{-}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})=\int _{l_{0}}^{l}[C^{+}(f)-\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+}]\frac{\text{d}l^{\prime }}{v_{\Vert }}+X(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E}),\end{eqnarray}$$

where $l_{0}$ is arbitrary for the time being and $X$ is an integration constant that needs to be determined from the orbit average of (2.4),

(2.13) $$\begin{eqnarray}\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{-}}=\overline{C^{-}(f)}.\end{eqnarray}$$

Our remaining task is to simplify and solve these equations. Of course, since the beginnings of neoclassical transport theory, this has been accomplished before in various approximations. In particular, the $1/\unicode[STIX]{x1D708}$ -regime can be treated in general stellarator geometry, but to our knowledge the bootstrap current has not been calculated before in the $\sqrt{\unicode[STIX]{x1D708}}$ -regime. In order to do so, we must first carefully go through the conventional $1/\unicode[STIX]{x1D708}$ -regime (until § 5), recovering familiar results but carefully identifying the assumptions made on the way.

2.1 Even part of the distribution function

The ratio of the right-hand side to the left-hand side of (2.11) is of the order of $\unicode[STIX]{x1D708}_{\ast }/\unicode[STIX]{x1D70C}_{\ast }$ , where $\unicode[STIX]{x1D70C}_{\ast }=\unicode[STIX]{x1D70C}/L$ is the normalised gyroradius and $\unicode[STIX]{x1D708}_{\ast }=\unicode[STIX]{x1D708}L/v_{T}$ the collisionality. Here $L$ denotes the macroscopic system length (typically the major radius) and $v_{T}=(2T/m)^{1/2}$ the thermal velocity. If $\unicode[STIX]{x1D70C}_{\ast }\ll \unicode[STIX]{x1D708}_{\ast }$ , collisions dominate and the distribution function is forced to be close to a Maxwellian. However, at lower collisionality, $\unicode[STIX]{x1D708}_{\ast }\ll \unicode[STIX]{x1D70C}_{\ast }$ , the plasma will in general not be in local thermodynamic equilibrium, since some particle orbits drift out of the plasma, leaving behind a ‘loss cone’ in velocity space (Helander Reference Helander2014; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017). We are aware of two limits in which the particles are well confined and confinement is restored. One limit is described by Ho & Kulsrud (Reference Ho and Kulsrud1987), who took the $\boldsymbol{E}\times \boldsymbol{B}$ drift,

(2.14) $$\begin{eqnarray}\boldsymbol{v}_{E}=-\frac{1}{B}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}\times \boldsymbol{b}\simeq \frac{1}{B}\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D713}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$

to be much larger than the $\unicode[STIX]{x1D735}B$ and curvature drifts,

(2.15) $$\begin{eqnarray}\boldsymbol{v}_{M}=\frac{v_{\Vert }^{2}}{\unicode[STIX]{x1D6FA}}\boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b})+\frac{\unicode[STIX]{x1D707}}{m\unicode[STIX]{x1D6FA}}\boldsymbol{b}\times \unicode[STIX]{x1D735}B,\end{eqnarray}$$

with $\unicode[STIX]{x1D6FA}$ the gyrofrequency, i.e.

(2.16) $$\begin{eqnarray}\boldsymbol{v}_{E}\gg \boldsymbol{v}_{M}.\end{eqnarray}$$

This limit is consistent with large aspect ratio, but it is also valid in regions of large plasma gradients of a low aspect ratio stellarator. The other possible limit considered in the literature is that of a stellarator sufficiently optimised that it can be considered close to omnigeneous (Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017). In this limit, the orbit average of the radial drift is small compared with its local value,

(2.17a,b ) $$\begin{eqnarray}\boldsymbol{v}_{E}\sim \boldsymbol{v}_{M},\quad (\boldsymbol{v}_{E}+\boldsymbol{v}_{M})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\gg \overline{(\boldsymbol{v}_{E}+\boldsymbol{v}_{M})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}.\end{eqnarray}$$

In both these limits, the bounce-averaged orbits stay close to flux surfaces, orbit confinement is restored and the distribution function becomes approximately Maxwellian and constant on flux surfaces. Moreover, quasineutrality leads to an electrostatic potential $\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D713})$ that is constant on flux surfaces to lowest order (Ho & Kulsrud Reference Ho and Kulsrud1987; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017).

Thus, assuming that the plasma is either in the $1/\unicode[STIX]{x1D708}$ -regime or follows the orderings (2.16) or (2.17), we shall write the solution to (2.11) as

(2.18) $$\begin{eqnarray}f^{+}=F_{0}+F_{1},\end{eqnarray}$$

where $F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D713})$ is a Maxwellian that is constant on flux surfaces and $F_{1}\ll F_{0}$ is constant along magnetic field lines, that is, independent of $l$ for trapped particles, and independent of $\unicode[STIX]{x1D6FC}$ and $l$ for circulating particles. The piece $F_{1}$ is determined by

(2.19) $$\begin{eqnarray}\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}\,\frac{\unicode[STIX]{x2202}F_{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}-\overline{C^{+}(F_{1})}=-\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}\,\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}.\end{eqnarray}$$

The solution to this equation depends on the collisionality and the radial electric field, which enters in $\boldsymbol{v}_{d}$ . Since

(2.20) $$\begin{eqnarray}\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}\,\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}=0\end{eqnarray}$$

in the untrapped part of phase space, there is no ‘drive’ there, and most simple collision operators, including (1.1), will lead to $F_{1}=0$ for circulating particles. The full operator, however, includes a field-particle term that does not have this property (Helander & Sigmar Reference Helander and Sigmar2002), but we shall ignore this complication. The function $F_{1}$ in the circulating domain is small compared to $F_{1}$ in the trapped-particle region if we use the assumptions in (2.16) and (2.17) (Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017).

Equation (2.19), with the approximation (1.1) for the collision operator, is routinely solved by the GSRAKE code (Beidler & DHaeseleer Reference Beidler and DHaeseleer1995), which uses the result to calculate the cross-field transport. This calculation is very rapid compared with codes that attempt to solve the drift kinetic equation without first taking the bounce average.

We consider two low-collisionality limits of (2.19). In the higher-collisionality limit (corresponding to the $1/\unicode[STIX]{x1D708}$ -regime), $C(F_{1})\gg \overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x2202}F_{1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})$ . Thus, $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x2202}F_{1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})\ll \overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}(\unicode[STIX]{x2202}F_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})$ , which in general implies

(2.21) $$\begin{eqnarray}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{1}\ll \boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{0}.\end{eqnarray}$$

In the lower-collisionality limit (corresponding to the $\sqrt{\unicode[STIX]{x1D708}}$ -regime),

(2.22) $$\begin{eqnarray}C(F_{1})\ll \overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x2202}F_{1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}),\end{eqnarray}$$

leading to

(2.23) $$\begin{eqnarray}C(F_{1})\ll \boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{1}\sim \boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{0}.\end{eqnarray}$$

Sections 3–5 below deal with the $1/\unicode[STIX]{x1D708}$ -regime and § 6 with the $\sqrt{\unicode[STIX]{x1D708}}$ -regime.

2.2 Odd part of the distribution function

The odd part of the distribution function is determined by (2.12)–(2.13) and must vanish at the bounce points. This condition alone determines the integration constant $X$ for trapped orbits; specifically, if we choose $l_{0}$ to be a bounce point then $X=0$ . It is thus natural to let $l_{0}$ depend on $\unicode[STIX]{x1D706}$ in the following way:

(2.24) $$\begin{eqnarray}\displaystyle B(l_{0})=\left\{\begin{array}{@{}ll@{}}1/\unicode[STIX]{x1D706},\quad & \unicode[STIX]{x1D706}>1/B_{max},\\ B_{max},\quad & \unicode[STIX]{x1D706}<1/B_{max}.\end{array}\right. & & \displaystyle\end{eqnarray}$$

In order to calculate $X$ for untrapped orbits, we must solve (2.13), which we write as

(2.25) $$\begin{eqnarray}\left\langle \frac{B}{v_{\Vert }}C^{-}(f)\right\rangle =\left\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{Bf^{-}}{v_{\Vert }}\boldsymbol{v}_{d}\right)\right\rangle =\frac{1}{V^{\prime }(\unicode[STIX]{x1D713})}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left\langle \frac{V^{\prime }(\unicode[STIX]{x1D713})Bf^{-}}{v_{\Vert }}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right\rangle ,\end{eqnarray}$$

where $V(\unicode[STIX]{x1D713})$ denotes the volume contained by the flux surface labelled by $\unicode[STIX]{x1D713}$ , $V^{\prime }(\unicode[STIX]{x1D713})=\text{d}V/\text{d}\unicode[STIX]{x1D713}$ , and we have used $\boldsymbol{v}_{d}=(v_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x1D735}\times (v_{\Vert }\boldsymbol{b})$ . On the right-hand side of this expression appears

(2.26) $$\begin{eqnarray}\left\langle \frac{Bf^{-}}{v_{\Vert }}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right\rangle =\left\langle \frac{B\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{v_{\Vert }}\left[\int _{l_{0}}^{l}(C^{+}(f)-\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+})\frac{\text{d}l^{\prime }}{v_{\Vert }}+X\right]\right\rangle ,\end{eqnarray}$$

where the last term vanishes because the constant $X$ is a flux function for untrapped orbits and (Helander Reference Helander2014)

(2.27) $$\begin{eqnarray}\left\langle \frac{B\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{v_{\Vert }}\right\rangle =0.\end{eqnarray}$$

Using (2.18) and the fact that $f^{+}=F_{0}+F_{1}$ does not depend on $\unicode[STIX]{x1D6FC}$ or $l$ in the circulating-particle region, equation (2.25) can be rewritten as

(2.28) $$\begin{eqnarray}\left\langle \frac{B}{v_{\Vert }}C^{-}(f)\right\rangle =\frac{1}{V^{\prime }(\unicode[STIX]{x1D713})}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }(\unicode[STIX]{x1D713})\left\langle \frac{B\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{v_{\Vert }}\int _{l_{0}}^{l}\left(C(F_{1})-\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}f^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)\frac{\text{d}l^{\prime }}{v_{\Vert }}\right\rangle .\end{eqnarray}$$

If we write

(2.29) $$\begin{eqnarray}h(l,\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})=\int _{l_{0}}^{l}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\text{d}l^{\prime }}{v_{\Vert }},\end{eqnarray}$$

then $\unicode[STIX]{x1D735}_{\Vert }h=(\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713})/v_{\Vert }$ and

(2.30) $$\begin{eqnarray}\left\langle \frac{B\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{v_{\Vert }}\frac{\unicode[STIX]{x2202}f^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\int _{l_{0}}^{l}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\text{d}l^{\prime }}{v_{\Vert }}\right\rangle =\frac{\unicode[STIX]{x2202}f^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left\langle \frac{B}{2}\unicode[STIX]{x1D735}_{\Vert }h^{2}\right\rangle =0,\end{eqnarray}$$

since $\langle B\unicode[STIX]{x1D735}_{\Vert }g\rangle$ vanishes for any function $g$ , see equation (21) of Helander (Reference Helander2014). Equation (2.28) thus reduces to

(2.31) $$\begin{eqnarray}\left\langle \frac{B}{v_{\Vert }}C^{-}(f)\right\rangle =\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }\left\langle \frac{B\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{v_{\Vert }}\int _{l_{0}}^{l}C(F_{1})\frac{\text{d}l^{\prime }}{v_{\Vert }}\right\rangle ,\end{eqnarray}$$

where the right-hand side is small compared with the left-hand side, since

(2.32) $$\begin{eqnarray}\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }\left.\left\langle \frac{B\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{v_{\Vert }}\int _{l_{0}}^{l}C(F_{1})\frac{\text{d}l^{\prime }}{v_{\Vert }}\right\rangle \right/\left\langle \frac{B}{v_{\Vert }}C^{-}(f)\right\rangle \sim \frac{\unicode[STIX]{x1D70C}_{\ast }F_{1}}{f^{-}}\sim \frac{F_{1}}{F_{0}}\ll 1,\end{eqnarray}$$

and we have used $f^{-}\sim \unicode[STIX]{x1D70C}_{\ast }F_{0}$ , which is to be verified a posteriori. Furthermore, since we have seen that $F_{1}=0$ for untrapped particles (for the collision operator we employ), we can thus conclude that the equation determining the integration constant $X$ to lowest order is

(2.33) $$\begin{eqnarray}\left\langle \frac{B}{v_{\Vert }}C^{-}(f)\right\rangle =0\end{eqnarray}$$

in all low-collisionality regimes.

Once the integration constant has been obtained, it is straightforward to calculate the flow velocity along the magnetic field of the particle species in question,

(2.34) $$\begin{eqnarray}\displaystyle V_{a\Vert }(l) & = & \displaystyle \frac{1}{n_{a}}\int v_{\Vert }\,f_{a}^{-}d^{3}v\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{n_{a}}\left(-B(l)\int _{l_{max}}^{l}\frac{\text{d}l^{\prime }}{B(l^{\prime })}\int \boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{a}^{+}d^{3}v+\int v_{\Vert }X_{a}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})d^{3}v\right)\nonumber\\ \displaystyle & = & \displaystyle V_{a1}+V_{a2}.\end{eqnarray}$$

Note that the term with the collision operator does not contribute. In the $1/\unicode[STIX]{x1D708}$ regime, according to (2.21), $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{a}^{+}\simeq \boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{0}$ . If we use this approximation and the expression (e.g. from Landreman & Catto (Reference Landreman and Catto2012) or equation (63) of Helander (Reference Helander2014))

(2.35) $$\begin{eqnarray}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=\frac{m_{a}v^{2}}{e_{a}}\unicode[STIX]{x1D709}(\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{\unicode[STIX]{x1D709}}{B}\right),\end{eqnarray}$$

where the gradient of $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D70E}(1-\unicode[STIX]{x1D706}B)^{1/2}$ is taken at constant $\unicode[STIX]{x1D706}$ , we obtain for the first term

(2.36) $$\begin{eqnarray}\langle V_{a1}B\rangle =-\!\left\langle \frac{m_{a}B^{2}}{n_{a}e_{a}}\int _{0}^{\infty }\frac{\unicode[STIX]{x2202}F_{a0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}2\unicode[STIX]{x03C0}v^{4}\,\text{d}v\int _{l_{max}}^{l}\,\text{d}l^{\prime }\int _{0}^{1/B}(\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{\unicode[STIX]{x1D709}}{B}\right)\text{d}\unicode[STIX]{x1D706}\right\rangle .\end{eqnarray}$$

Evaluating the integrals, we find (Helander et al. Reference Helander, Geiger and Maassberg2011)

(2.37) $$\begin{eqnarray}\langle V_{a1}B\rangle =\frac{T_{a}A_{1a}}{e_{a}}\langle uB^{2}\rangle ,\end{eqnarray}$$

where we have recalled (1.7) and fixed the integration constant therein by taking $u=0$ at the point on the flux surface where $B=B_{max}$ . (In the notation of Nakajima et al. (Reference Nakajima, Okamoto, Todoroki, Nakamura and Wakatani1989) and Helander et al. (Reference Helander, Geiger and Maassberg2011), $g_{2}=-uB^{2}$ .) In addition, we have defined the thermodynamic forces

(2.38) $$\begin{eqnarray}\displaystyle & \displaystyle A_{1a}=\frac{\text{d}\ln p_{a}}{\text{d}\unicode[STIX]{x1D713}}+\frac{e_{a}}{T_{a}}\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle A_{2a}=\frac{\text{d}\ln T_{a}}{\text{d}\unicode[STIX]{x1D713}}. & \displaystyle\end{eqnarray}$$

The approximation $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{a}^{+}\simeq \boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{a0}$ is, however, not appropriate in all collisionality regimes and will be discussed further below, but first we shall solve (2.33) in various cases, still using this approximation, and calculate the corresponding velocities $V_{a2}$ .

3 Lorentz limit

To solve (2.33), we must specify the collision operator. We first consider pure pitch-angle scattering, i.e. we retain only the first term on the right-hand side of (1.1). This applies to electrons in a high- $Z$ plasma and is used in the widespread DKES code (Hirshman et al. Reference Hirshman, Shaing, van Rij, Beasley and Crume1986; van Rij & Hirshman Reference van Rij and Hirshman1989).

In the coordinates $v$ and $\unicode[STIX]{x1D706}$ , the operator ${\mathcal{L}}$ becomes

(3.1) $$\begin{eqnarray}{\mathcal{L}}(f)=\frac{2\unicode[STIX]{x1D709}}{B}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\left(\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right).\end{eqnarray}$$

Our equation (2.33) then becomes

(3.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\unicode[STIX]{x1D706}\left\langle \unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}f^{-}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right\rangle =0,\end{eqnarray}$$

and is readily integrated once, demanding regularity at $\unicode[STIX]{x1D706}=0$ ,

(3.3) $$\begin{eqnarray}\left\langle \unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}f^{-}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right\rangle =0.\end{eqnarray}$$

We substitute (2.12) and note that the term involving the collision operator vanishes in the untrapped part of phase space since $F_{1}=0$ there. Hence we obtain

(3.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}X}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=\frac{1}{\langle \unicode[STIX]{x1D709}\rangle }\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left\langle \unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\int _{l_{max}}^{l}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\text{d}l^{\prime }}{v_{\Vert }(l^{\prime })}\right\rangle ,\end{eqnarray}$$

where $B(l_{max})=B_{max}$ . We thus find

(3.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}X_{a}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=-\frac{m_{a}v}{2e_{a}}\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\frac{\langle g_{4}\rangle }{\langle \unicode[STIX]{x1D709}\rangle },\end{eqnarray}$$

where we have introduced the notation (Nakajima et al. Reference Nakajima, Okamoto, Todoroki, Nakamura and Wakatani1989)

(3.6) $$\begin{eqnarray}g_{4}(\unicode[STIX]{x1D706},l)=\unicode[STIX]{x1D709}\int _{l_{max}}^{l}(\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}^{-1}\,\text{d}l^{\prime },\quad (\unicode[STIX]{x1D706}<1/B_{max}).\end{eqnarray}$$

Knowing $X_{a}$ , we can now calculate

(3.7) $$\begin{eqnarray}V_{a2}=\frac{1}{n_{a}}\int v_{\Vert }X_{a}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})\,\text{d}^{3}v=-\frac{1}{n_{a}}\int v_{\Vert }\unicode[STIX]{x1D706}\frac{\unicode[STIX]{x2202}X_{a}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\,\text{d}^{3}v,\end{eqnarray}$$

which becomes

(3.8) $$\begin{eqnarray}\langle V_{a2}B\rangle =\frac{f_{s}T_{a}A_{1a}}{e_{a}},\end{eqnarray}$$

where we have written

(3.9) $$\begin{eqnarray}f_{s}=\frac{3\langle B^{2}\rangle }{4}\int _{0}^{1/B_{max}}\frac{\langle g_{4}\rangle \unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D706}}{\left\langle \sqrt{1-\unicode[STIX]{x1D706}B}\right\rangle }.\end{eqnarray}$$

The current along $\boldsymbol{B}$ is thus equal to (Helander et al. Reference Helander, Geiger and Maassberg2011)

(3.10) $$\begin{eqnarray}\langle J_{a\Vert }B\rangle =n_{a}e_{a}\langle (V_{a1}+V_{a2})B\rangle =p_{a}A_{1a}(f_{s}+\langle uB^{2}\rangle )\end{eqnarray}$$

in this approximation. Note that this current depends on the radial electric field in $A_{1a}$ but disappears, as it should, from the total current once a sum is taken over all particle species. However, there are perils of using a collision operator that does not conserve momentum. In a tokamak or quasisymmetric stellarator, the resulting transport depends on the radial electric field, but this dependence disappears if a momentum-conserving operator is used. In addition, there is spurious particle transport from self-collisions, which vanishes if momentum is conserved.

4 Rotation velocity of a pure plasma

A better approximation, particularly for the plasma ions ( $a=i$ ), which are assumed to be singly charged, $e_{i}=e$ , is obtained if the momentum-restoring term is retained in the collision operator (1.1). For the bulk ions in a plasma without impurities, we thus take

(4.1) $$\begin{eqnarray}C_{i}^{-}(f_{i})=\unicode[STIX]{x1D708}_{D}^{ii}\left({\mathcal{L}}(f_{i}^{-})+\frac{m_{i}v_{\Vert }}{T_{i}}w_{i}F_{i0}\right),\end{eqnarray}$$

where

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{D}^{ii}=\frac{3\unicode[STIX]{x03C0}^{1/2}}{4\unicode[STIX]{x1D70F}_{ii}}\frac{\unicode[STIX]{x1D719}(x)-G(x)}{x^{3}},\end{eqnarray}$$

denotes the ion–ion pitch-angle scattering frequency, with

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ii}=6\sqrt{2}\unicode[STIX]{x03C0}^{3/2}\unicode[STIX]{x1D716}_{0}^{2}m_{i}^{1/2}T_{i}^{3/2}/e^{4}n_{i}\ln \unicode[STIX]{x1D6EC},\end{eqnarray}$$

where $\ln \unicode[STIX]{x1D6EC}$ is the Coulomb logarithm, $\unicode[STIX]{x1D716}_{0}$ the vacuum permittivity, $x=v/v_{Ti}$ , $v_{Ti}=\sqrt{2T_{i}/m_{i}}$ and

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}(x)=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{x}\text{e}^{-y^{2}}\,\text{d}y, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle G(x)=\frac{\unicode[STIX]{x1D719}(x)-x\unicode[STIX]{x1D719}^{\prime }(x)}{2x^{2}}. & \displaystyle\end{eqnarray}$$

Instead of (3.3) we now obtain

(4.6) $$\begin{eqnarray}\left\langle \unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}f_{i}^{-}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right\rangle =-\frac{m_{i}v}{2T_{i}}\langle w_{i}B\rangle F_{i0},\end{eqnarray}$$

and instead of (3.5)

(4.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}X_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=-\frac{m_{i}v}{2T_{i}\langle \unicode[STIX]{x1D709}\rangle }\left(\langle w_{i}B\rangle F_{i0}+\frac{T_{i}\langle g_{4}\rangle }{e}\frac{\unicode[STIX]{x2202}F_{i0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right).\end{eqnarray}$$

The coefficient $w_{i}$ is determined by momentum conservation,

(4.8) $$\begin{eqnarray}0=\int v_{\Vert }C_{i}^{-}(f_{i})\,\text{d}^{3}v=n_{i}w_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}-\int \unicode[STIX]{x1D708}_{D}^{ii}v_{\Vert }\,f_{i}^{-}\text{d}^{3}v,\end{eqnarray}$$

where we have defined the velocity-space average (Helander & Sigmar Reference Helander and Sigmar2002)

(4.9) $$\begin{eqnarray}\{\unicode[STIX]{x1D708}_{D}^{ii}\}=\int \unicode[STIX]{x1D708}_{D}^{ii}\frac{m_{i}v^{2}}{3n_{i}T_{i}}F_{i0}\,\text{d}^{3}v=\frac{\sqrt{2}-\ln \left(1+\sqrt{2}\right)}{\unicode[STIX]{x1D70F}_{ii}}\simeq \frac{0.533}{\unicode[STIX]{x1D70F}_{ii}}.\end{eqnarray}$$

Upon substitution of (2.12) we find

(4.10) $$\begin{eqnarray}\langle w_{i}B\rangle =\frac{1}{n_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\left\langle B^{2}\int _{0}^{\infty }\unicode[STIX]{x1D708}_{D}^{ii}2\unicode[STIX]{x03C0}v^{3}\,\text{d}v\int _{0}^{1/B}\,\text{d}\unicode[STIX]{x1D706}\left(X_{i}-\int _{l_{0}}^{l}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{+}\frac{\text{d}l^{\prime }}{|v_{\Vert }|}\right)\right\rangle ,\end{eqnarray}$$

where we have noted that the part of the distribution function (2.12) containing $C^{+}(f)$ does not contribute to the integral:

(4.11) $$\begin{eqnarray}\int \unicode[STIX]{x1D708}_{D}^{ii}v_{\Vert }\,\text{d}^{3}v\int _{l_{0}}^{l}C^{+}(f)\frac{\text{d}l^{\prime }}{v_{\Vert }}=B(l)\int _{l_{max}}^{l}\frac{\text{d}l^{\prime }}{B(l^{\prime })}\int \unicode[STIX]{x1D708}_{D}^{ii}C^{+}(f)\,\text{d}^{3}v=0,\end{eqnarray}$$

for any collision operator of the form (1.1).

The first of the two terms on the right of (4.10) can be evaluated using (3.5) and is equal to

(4.12) $$\begin{eqnarray}-\frac{1}{n_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\left\langle B^{2}\int _{0}^{\infty }\unicode[STIX]{x1D708}_{D}^{ii}2\unicode[STIX]{x03C0}v^{3}\,\text{d}v\int _{0}^{1/B_{max}}\unicode[STIX]{x1D706}\frac{\unicode[STIX]{x2202}X_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\,\text{d}\unicode[STIX]{x1D706}\right\rangle =f_{c}\langle w_{i}B\rangle +\frac{f_{s}T_{i}}{e}(A_{1i}-\unicode[STIX]{x1D702}A_{2i}),\end{eqnarray}$$

with

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle f_{c}=\frac{3\langle B^{2}\rangle }{4}\int _{0}^{1/B_{max}}\frac{\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D706}}{\left\langle \sqrt{1-\unicode[STIX]{x1D706}B}\right\rangle }, & \displaystyle\end{eqnarray}$$

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\frac{\{\unicode[STIX]{x1D708}_{D}^{ii}(5/2-x^{2})\}}{\{\unicode[STIX]{x1D708}_{D}^{ii}\}}=\frac{5}{2}-\frac{1}{2-\sqrt{2}\ln (1+\sqrt{2})}\simeq 1.17, & \displaystyle\end{eqnarray}$$

The second term on the right of (4.10) becomes

(4.15) $$\begin{eqnarray}-\frac{1}{n_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\left\langle B^{2}\int _{0}^{\infty }\unicode[STIX]{x1D708}_{D}^{ii}2\unicode[STIX]{x03C0}v^{3}\,\text{d}v\int _{0}^{1/B}\,\text{d}\unicode[STIX]{x1D706}\int _{l_{0}}^{l}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{+}\frac{\text{d}l^{\prime }}{|v_{\Vert }|}\right\rangle =\frac{T_{i}\langle uB^{2}\rangle }{e}(A_{1i}-\unicode[STIX]{x1D702}A_{2i}),\end{eqnarray}$$

where we have again approximated $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+}$ by $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{0}$ . As long as this approximation is allowed, equation (4.10) thus implies

(4.16) $$\begin{eqnarray}\langle w_{i}B\rangle =\frac{T_{i}(f_{s}+\langle uB^{2}\rangle )}{e(1-f_{c})}(A_{1i}-\unicode[STIX]{x1D702}A_{2i}),\end{eqnarray}$$

which together with (4.7) determines the integration constant $X_{i}$ and thus the odd part of the ion distribution function.

We are now in a position to calculate the plasma flow velocity along the magnetic field, $\langle V_{i\Vert }B\rangle =\langle (V_{i1}+V_{i2})B\rangle$ . The first term was calculated in (2.37) and the second term is

(4.17) $$\begin{eqnarray}\langle V_{i2}B\rangle =-\frac{2\unicode[STIX]{x03C0}}{n_{i}}\left\langle B^{2}\int _{0}^{\infty }v^{2}\,\text{d}v\int _{0}^{1/B_{max}}\frac{\unicode[STIX]{x2202}X_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D706}\right\rangle =f_{c}\langle w_{i}B\rangle +\frac{f_{s}T_{i}A_{1i}}{e},\end{eqnarray}$$

where we have used (3.5). We thus obtain

(4.18) $$\begin{eqnarray}\langle V_{i\Vert }B\rangle =\frac{T_{i}}{e}\frac{f_{s}+\langle uB^{2}\rangle }{1-f_{c}}(A_{1i}-\unicode[STIX]{x1D702}f_{c}A_{2i}).\end{eqnarray}$$

Apart from the geometric factor, this result is similar to that in tokamaks and was obtained by Landreman & Catto (Reference Landreman and Catto2012) for the case of a perfectly omnigenous stellarator. From our analysis, it is clear that the same result is also valid in the $1/\unicode[STIX]{x1D708}$ -regime of any non-optimised stellarator. A similar conclusion follows from the moment method of Hirshman and Sigmar extended to stellarator geometry by Shaing & Callen (Reference Shaing and Callen1983), in whose paper this expression can also be found. A refined expression is obtained in the next section using a more accurate collision operator. We note, however, that it is only valid if the $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{+}$ can indeed be approximated by $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{i0}$ in (2.12). As we shall see below, the correction that arises in the $\sqrt{\unicode[STIX]{x1D708}}$ -regime can be significant.

In an axisymmetric magnetic field, $\boldsymbol{B}=F(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+\unicode[STIX]{x1D704}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ ,

(4.19) $$\begin{eqnarray}f_{s}=\frac{F}{\unicode[STIX]{x1D704}}\left(f_{c}-\frac{\langle B^{2}\rangle }{B_{max}^{2}}\right),\end{eqnarray}$$

and

(4.20) $$\begin{eqnarray}\langle uB^{2}\rangle =\frac{F}{\unicode[STIX]{x1D704}}\left(\frac{\langle B^{2}\rangle }{B_{max}^{2}}-1\right),\end{eqnarray}$$

so that $f_{s}+\langle uB^{2}\rangle =-(1-f_{c})F/\unicode[STIX]{x1D704}$ . The ion flow velocity (4.18) thus reduces to the usual result (Helander & Sigmar Reference Helander and Sigmar2002)

(4.21) $$\begin{eqnarray}\langle V_{i\Vert }B\rangle =-\frac{F}{\unicode[STIX]{x1D704}}\frac{T_{i}}{e}(A_{1i}-\unicode[STIX]{x1D702}f_{c}A_{2i}).\end{eqnarray}$$

5 Bootstrap current

To calculate the bootstrap current, we also need to determine the distribution function for the electrons. In order to obtain a result that matches that for a standard large aspect ratio tokamak (Rosenbluth, Hazeltine & Hinton Reference Rosenbluth, Hazeltine and Hinton1972; Helander & Sigmar Reference Helander and Sigmar2002), we use a model linearised electron collision operator $C_{e}(f_{e})$ that is more accurate than the operator in (4.1). Rosenbluth et al. (Reference Rosenbluth, Hazeltine and Hinton1972) showed that the operator (1.1) is sufficient in the limit $f_{t}=1-f_{c}\rightarrow 0$ , and it later emerged that improved accuracy can be achieved at finite $f_{t}$ by letting the momentum-restoring coefficient $w$ depend on energy. We therefore employ such an operator as the most accurate one allowing for an explicit analytical solution of the drift kinetic equation, but remind the reader that its accuracy is not perfect (Parker & Catto Reference Parker and Catto2012).

The electron collision operator must include both electron–electron and electron–ion collisions,

(5.1) $$\begin{eqnarray}C_{e}^{-}(f_{e})=C_{ee}(f_{e}^{-})+C_{ei}(f_{e}^{-}).\end{eqnarray}$$

To lowest order in $\sqrt{m_{e}/m_{i}}\ll 1$ , the electron–ion collision operator is

(5.2) $$\begin{eqnarray}C_{ei}(f_{e}^{-})=\unicode[STIX]{x1D708}_{D}^{ei}{\mathcal{L}}(f_{e}^{-})+\frac{\unicode[STIX]{x1D708}_{D}^{ei}m_{e}v_{\Vert }V_{i\Vert }}{T_{e}}F_{e0},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{D}^{ei}=(3\sqrt{\unicode[STIX]{x03C0}}/4\unicode[STIX]{x1D70F}_{ei}x^{3})$ , $x=v/v_{Te}$ , $v_{Te}=\sqrt{2T_{e}/m_{e}}$ and, for singly charged ions, $\unicode[STIX]{x1D70F}_{ei}=\unicode[STIX]{x1D70F}_{ee}=6\sqrt{2}\unicode[STIX]{x03C0}^{3/2}\unicode[STIX]{x1D716}_{0}^{2}m_{e}^{1/2}T_{e}^{3/2}/e^{4}n_{e}\ln \unicode[STIX]{x1D6EC}$ . For $C_{ee}$ , we use a model collision operator (Taguchi Reference Taguchi1988; Helander & Sigmar Reference Helander and Sigmar2002) that has the exact pitch-angle scattering piece of the Landau collision operator, conserves momentum, is self-adjoint, has $f_{e}^{-}=(m_{e}v_{\Vert }V_{e\Vert }/T_{e})F_{e0}$ as the only solution to $C_{ee}(f_{e}^{-})=0$ and in addition, we request that it satisfies

(5.3) $$\begin{eqnarray}\int L_{p}^{(3/2)}(x^{2})v_{\Vert }C_{ee}(L_{q}^{(3/2)}(x^{2})v_{\Vert }F_{e0})\,\text{d}^{3}v=-\frac{K_{pq}n_{e}T_{e}}{m_{e}\unicode[STIX]{x1D70F}_{ee}}\end{eqnarray}$$

for $p,q=0,1,\ldots ,S$ . Here $L_{p}^{(3/2)}(y)$ are Sonine polynomials, and $K_{pq}$ are coefficients that result from calculating $\int L_{p}^{(3/2)}(x^{2})v_{\Vert }C_{ee}(L_{q}^{(3/2)}(x^{2})v_{\Vert }F_{e0})\,\text{d}^{3}v$ using the linearised Landau electron–electron collision operator. The first few Sonine polynomials are

(5.4a-c ) $$\begin{eqnarray}L_{0}^{(3/2)}(y)=1,\quad L_{1}^{(3/2)}(y)=\frac{5}{2}-y,\quad L_{2}^{(3/2)}(y)=\frac{35}{8}-\frac{7y}{2}+\frac{y^{2}}{2},\end{eqnarray}$$

and the corresponding coefficients $K_{pq}$ are

(5.5a-d ) $$\begin{eqnarray}K_{p0}=K_{0q}=0,\quad K_{11}=\sqrt{2},\quad K_{12}=K_{21}=\frac{3\sqrt{2}}{4},\quad K_{22}=\frac{45\sqrt{2}}{16}.\end{eqnarray}$$

We derive an operator that satisfies these conditions in appendix A. For our calculation, it is sufficient to satisfy condition (5.3) for $p,q=0,1,2$ . Thus, we use the model collision operator

(5.6) $$\begin{eqnarray}C_{ee}(f_{e}^{-})=\unicode[STIX]{x1D708}_{D}^{ee}{\mathcal{L}}(f_{e}^{-})+\frac{\unicode[STIX]{x1D708}_{D}^{ee}m_{e}v_{\Vert }}{T_{e}}w_{e}F_{e0},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{D}^{ee}=(3\sqrt{\unicode[STIX]{x03C0}}/4\unicode[STIX]{x1D70F}_{ee})[\unicode[STIX]{x1D719}(x)-G(x)]/x^{3}$ , and $w_{e}$ is not a constant but a function of $v$ given by

(5.7) $$\begin{eqnarray}w_{e}=a_{0}+a_{1}L_{1}^{(3/2)}(x^{2})+a_{2}L_{2}^{(3/2)}(x^{2}).\end{eqnarray}$$

The coefficients $a_{0}$ , $a_{1}$ and $a_{2}$ are determined by the linear set of equations

(5.8) $$\begin{eqnarray}\mathop{\sum }_{q=0}^{2}\unicode[STIX]{x1D615}_{pq}a_{q}=\frac{\unicode[STIX]{x1D70F}_{ee}}{n_{e}}\int \unicode[STIX]{x1D708}_{D}^{ee}L_{p}^{(3/2)}(x^{2})v_{\Vert }\,f_{e}^{-}\,\text{d}^{3}v,\end{eqnarray}$$

where the symmetric matrix $\unicode[STIX]{x1D615}_{pq}$ is given by

(5.9) $$\begin{eqnarray}\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D615}_{00} & \unicode[STIX]{x1D615}_{01} & \unicode[STIX]{x1D615}_{02}\\ \unicode[STIX]{x1D615}_{10} & \unicode[STIX]{x1D615}_{11} & \unicode[STIX]{x1D615}_{12}\\ \unicode[STIX]{x1D615}_{20} & \unicode[STIX]{x1D615}_{21} & \unicode[STIX]{x1D615}_{22}\end{array}\right)=\left(\begin{array}{@{}ccc@{}}0.533 & 0.625 & 0.652\\ 0.625 & 0.0240 & -0.288\\ 0.652 & -0.288 & -0.947\end{array}\right).\end{eqnarray}$$

In appendix A we show how to calculate this matrix.

Using the collision operators in (5.1), (5.2) and (5.6), we find that

(5.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}X_{e}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=-\frac{m_{e}v}{2T_{e}\langle \unicode[STIX]{x1D709}\rangle }\left[\mathop{\sum }_{p=0}^{2}\frac{\unicode[STIX]{x1D708}_{D}^{ee}L_{p}^{(3/2)}(x^{2})\langle a_{p}B\rangle }{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}F_{e0}+\frac{\unicode[STIX]{x1D708}_{D}^{ei}\langle V_{i\Vert }B\rangle }{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}F_{e0}-\frac{T_{e}\langle g_{4}\rangle }{e}\frac{\unicode[STIX]{x2202}F_{e0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right],\end{eqnarray}$$

where the coefficients $\langle a_{p}B\rangle$ are determined from the flux-surface average of (5.8) multiplied by $B$ :

(5.11) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{q=0}^{2}\left(\unicode[STIX]{x1D615}_{pq}-f_{c}\unicode[STIX]{x1D70F}_{ee}\left\{\frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}L_{p}^{(3/2)}L_{q}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}\right\}\right)\langle a_{q}B\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70F}_{ee}\left\{\frac{\unicode[STIX]{x1D708}_{D}^{ee}\unicode[STIX]{x1D708}_{D}^{ei}L_{p}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}\right\}f_{c}\langle V_{i\Vert }B\rangle -\frac{T_{e}\unicode[STIX]{x1D70F}_{ee}}{e}(f_{s}+\langle uB^{2}\rangle )\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,[\{\unicode[STIX]{x1D708}_{D}^{ee}L_{p}^{(3/2)}\}A_{1e}-\{\unicode[STIX]{x1D708}_{D}^{ee}L_{1}^{(3/2)}L_{p}^{(3/2)}\}A_{2e}].\end{eqnarray}$$

Using

(5.12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\left\{{\displaystyle \frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.198}{\unicode[STIX]{x1D70F}_{ee}}},\quad \left\{{\displaystyle \frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}L_{1}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.197}{\unicode[STIX]{x1D70F}_{ee}}},\quad \left\{{\displaystyle \frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}L_{2}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.179}{\unicode[STIX]{x1D70F}_{ee}}},\\ \left\{{\displaystyle \frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}(L_{1}^{(3/2)})^{2}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.458}{\unicode[STIX]{x1D70F}_{ee}}},\quad \left\{{\displaystyle \frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}L_{1}^{(3/2)}L_{2}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.457}{\unicode[STIX]{x1D70F}_{ee}}},\\ \left\{{\displaystyle \frac{(\unicode[STIX]{x1D708}_{D}^{ee})^{2}(L_{2}^{(3/2)})^{2}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.754}{\unicode[STIX]{x1D70F}_{ee}}},\\ \left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ee}\unicode[STIX]{x1D708}_{D}^{ei}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.335}{\unicode[STIX]{x1D70F}_{ee}}},\quad \left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ee}\unicode[STIX]{x1D708}_{D}^{ei}L_{1}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.428}{\unicode[STIX]{x1D70F}_{ee}}},\quad \left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ee}\unicode[STIX]{x1D708}_{D}^{ei}L_{2}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}={\displaystyle \frac{0.473}{\unicode[STIX]{x1D70F}_{ee}}},\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

we find the equations

(5.13) $$\begin{eqnarray}\displaystyle & & \displaystyle \left(\begin{array}{@{}ccc@{}}0.335+0.198f_{t} & 0.428+0.197f_{t} & 0.473+0.179f_{t}\\ 0.428+0.197f_{t} & -0.434+0.458f_{t} & -0.746+0.457f_{t}\\ 0.473+0.179f_{t} & -0.746+0.457f_{t} & -1.701+0.754f_{t}\end{array}\right)\left(\begin{array}{@{}c@{}}\langle a_{0}B\rangle \\ \langle a_{1}B\rangle \\ \langle a_{2}B\rangle \end{array}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{T_{e}}{e}(f_{s}+\langle uB^{2}\rangle )\left[\left(\begin{array}{@{}c@{}}-0.533\\ -0.625\\ -0.652\end{array}\right)A_{1e}+\left(\begin{array}{@{}c@{}}0.625\\ 1.386\\ 1.563\end{array}\right)A_{2e}\right]+\left(\begin{array}{@{}c@{}}0.335\\ 0.428\\ 0.473\end{array}\right)f_{c}\langle V_{i\Vert }B\rangle .\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Solving (5.13), we find

(5.14) $$\begin{eqnarray}\displaystyle \left(\begin{array}{@{}c@{}}\langle a_{0}B\rangle \\ \langle a_{1}B\rangle \\ \langle a_{2}B\rangle \end{array}\right) & = & \displaystyle \left(\begin{array}{@{}c@{}}-1.54+1.34f_{t}-0.22f_{t}^{2}\\ 0.01(1-f_{t}^{2})\\ -0.05+0.07f_{t}-0.02f_{t}^{2}\end{array}\right)\frac{T_{e}}{eD}(f_{s}+\langle uB^{2}\rangle )A_{1e}\nonumber\\ \displaystyle & & \displaystyle +\left(\begin{array}{@{}c@{}}2.72-1.63f_{t}+0.03f_{t}^{2}\\ -0.95-0.20f_{t}+0.21f_{t}^{2}\\ 0.25+0.08f_{t}+0.06f_{t}^{2}\end{array}\right)\frac{T_{e}}{eD}(f_{s}+\langle uB^{2}\rangle )A_{2e}\nonumber\\ \displaystyle & & \displaystyle +\left(\begin{array}{@{}c@{}}1-0.86f_{t}+0.12f_{t}^{2}\\ 0.01f_{t}(1+f_{t})\\ f_{t}(-0.05+0.02f_{t})\end{array}\right)\frac{f_{c}}{D}\langle V_{i\Vert }B\rangle ,\end{eqnarray}$$

where

(5.15) $$\begin{eqnarray}D=1-0.32f_{t}-0.36f_{t}^{2}+0.09f_{t}^{3}.\end{eqnarray}$$

With these results, we can now obtain the electron current. Using (2.34) and (5.10), we find

(5.16) $$\begin{eqnarray}\displaystyle \langle J_{\Vert }B\rangle & = & \displaystyle n_{e}T_{e}(f_{s}+\langle uB^{2}\rangle )A_{1e}-en_{e}\mathop{\sum }_{p=0}^{2}\left\{\frac{\unicode[STIX]{x1D708}_{D}^{ee}L_{p}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}\right\}f_{c}\langle a_{p}B\rangle \nonumber\\ \displaystyle & & \displaystyle +\,en_{e}\left\{\frac{\unicode[STIX]{x1D708}_{D}^{ee}+f_{t}\unicode[STIX]{x1D708}_{D}^{ei}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}\right\}\langle V_{i\Vert }B\rangle ,\end{eqnarray}$$

where

(5.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ee}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}=0.428,\quad \left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ee}L_{1}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}=-0.0563,\\ \left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ee}L_{2}^{(3/2)}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}=-0.0430,\quad \left\{{\displaystyle \frac{\unicode[STIX]{x1D708}_{D}^{ei}}{\unicode[STIX]{x1D708}_{D}^{ee}+\unicode[STIX]{x1D708}_{D}^{ei}}}\right\}=0.572.\end{array}\right\}\end{eqnarray}$$

Using (5.14), we thus obtain

(5.18) $$\begin{eqnarray}\displaystyle \langle J_{\Vert }B\rangle & = & \displaystyle \frac{1.66-1.55f_{t}+0.30f_{t}^{2}}{D}[n_{e}T_{e}(f_{s}+\langle uB^{2}\rangle )A_{1e}+f_{t}en_{e}\langle V_{i\Vert }B\rangle ]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{(1.21-0.69f_{t})f_{c}}{D}n_{e}T_{e}(f_{s}+\langle uB^{2}\rangle )A_{2e},\end{eqnarray}$$

This result has been obtained with a more sophisticated collision operator than the ion flow (4.18). To refine the latter, allowing for an energy-dependent momentum-restoring coefficient,

(5.19) $$\begin{eqnarray}w_{i}=b_{0}+b_{1}L_{1}^{(3/2)}(x^{2})+b_{2}L_{2}^{(3/2)}(x^{2}),\end{eqnarray}$$

we employ an operator similar to (5.6) for ion–ion collisions, which results in algebra very similar to (but simpler than) that for the electron current. Instead of (5.10) one obtains

(5.20) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}X_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=-\frac{m_{i}v}{2T_{i}\langle \unicode[STIX]{x1D709}\rangle }\left[\mathop{\sum }_{p=0}^{2}L_{p}^{(3/2)}(x^{2})\langle b_{p}B\rangle F_{i0}+\frac{T_{i}\langle g_{4}\rangle }{e}\frac{\unicode[STIX]{x2202}F_{i0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right],\end{eqnarray}$$

and in place of the next equation

(5.21) $$\begin{eqnarray}\mathop{\sum }_{q=0}^{2}(\unicode[STIX]{x1D615}_{pq}-f_{c}\unicode[STIX]{x1D70F}_{ii}\{\unicode[STIX]{x1D708}_{D}^{ii}L_{p}^{(3/2)}L_{q}^{(3/2)}\})\langle b_{q}B\rangle =\frac{T_{i}\unicode[STIX]{x1D70F}_{ii}}{e}(f_{s}+\langle uB^{2}\rangle )[\{\unicode[STIX]{x1D708}_{D}^{ii}L_{p}^{(3/2)}\}A_{1i}-\{\unicode[STIX]{x1D708}_{D}^{ii}L_{1}^{(3/2)}L_{p}^{(3/2)}\}A_{2i}].\end{eqnarray}$$

Using the velocity-space averages

(5.22) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\{\unicode[STIX]{x1D708}_{D}^{ii}\}={\displaystyle \frac{0.533}{\unicode[STIX]{x1D70F}_{ii}}},\quad \{\unicode[STIX]{x1D708}_{D}^{ii}L_{1}^{(3/2)}\}={\displaystyle \frac{0.625}{\unicode[STIX]{x1D70F}_{ii}}},\quad \{\unicode[STIX]{x1D708}_{D}^{ii}L_{2}^{(3/2)}\}={\displaystyle \frac{0.652}{\unicode[STIX]{x1D70F}_{ii}}},\\ \{\unicode[STIX]{x1D708}_{D}^{ii}(L_{1}^{(3/2)})^{2}\}={\displaystyle \frac{1.386}{\unicode[STIX]{x1D70F}_{ii}}},\quad \{\unicode[STIX]{x1D708}_{D}^{ii}L_{1}^{(3/2)}L_{2}^{(3/2)}\}={\displaystyle \frac{1.563}{\unicode[STIX]{x1D70F}_{ii}}}\\ \{\unicode[STIX]{x1D708}_{D}^{ii}(L_{2}^{(3/2)})^{2}\}={\displaystyle \frac{2.484}{\unicode[STIX]{x1D70F}_{ii}}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

one finds the following system of equations

(5.23) $$\begin{eqnarray}\displaystyle & & \displaystyle \left(\begin{array}{@{}ccc@{}}0.533f_{t} & 0.625f_{t} & 0.652f_{t}\\ 0.625f_{t} & -1.362+1.386f_{t} & -1.852+1.563f_{t}\\ 0.652f_{t} & -1.852+1.563f_{t} & -3.431+2.484f_{t}\end{array}\right)\left(\begin{array}{@{}c@{}}\langle b_{0}B\rangle \\ \langle b_{1}B\rangle \\ \langle b_{2}B\rangle \end{array}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{T_{i}}{e}(f_{s}+\langle uB^{2}\rangle )\left[\left(\begin{array}{@{}c@{}}0.533\\ 0.625\\ 0.652\end{array}\right)A_{1i}-\left(\begin{array}{@{}c@{}}0.625\\ 1.386\\ 1.563\end{array}\right)A_{2i}\right]\end{eqnarray}$$

and the ion flow

(5.24) $$\begin{eqnarray}\langle BV_{i\Vert }\rangle =\frac{T_{i}}{e}\frac{\langle uB^{2}\rangle +f_{s}}{f_{t}}(A_{1i}-\frac{(1.17-0.89f_{t})f_{c}}{1-1.27f_{t}+0.37f_{t}^{2}}A_{2i}).\end{eqnarray}$$

The final expression for the bootstrap current results from substituting this result in (5.18), giving

(5.25) $$\begin{eqnarray}\displaystyle \langle J_{\Vert }B\rangle & = & \displaystyle n(f_{s}+\langle uB^{2}\rangle )\frac{1.66-1.55f_{t}+0.30f_{t}^{2}}{D}\nonumber\\ \displaystyle & & \displaystyle \times \left[T_{e}A_{1e}-\frac{(0.73-0.42f_{t})f_{c}}{1-0.94f_{t}+0.18f_{t}^{2}}T_{e}A_{2e}+T_{i}A_{1i}-\frac{(1.17-0.89f_{t})f_{c}}{1-1.27f_{t}+0.37f_{t}^{2}}T_{i}A_{2i}\right].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Note that the radial electric field terms present in the thermodynamic forces $A_{1i}$ and $A_{1e}$ cancel, so that the bootstrap current is independent of the electric field in the $1/\unicode[STIX]{x1D708}$ -regime. To lowest order in $f_{t}$ , this expression coincides with that of Shaing et al. (Reference Shaing, Carreras, Dominguez, Lynch and Tolliver1989), who included two instead of three Sonine polynomials in the expansion.Footnote ² It also agrees with the result of Landreman & Catto (Reference Landreman and Catto2012), which it generalises to non-omnigenous stellarators in the $1/\unicode[STIX]{x1D708}$ -regime.

In a large aspect ratio tokamak, the result (5.25) reduces to

(5.26) $$\begin{eqnarray}\langle J_{\Vert }B\rangle =-1.66f_{t}\unicode[STIX]{x1D704}^{-1}F(\unicode[STIX]{x1D713})n(T_{e}A_{1e}+T_{i}A_{1i}-0.73T_{e}A_{2e}-1.17T_{i}A_{2i}),\end{eqnarray}$$

which closely approximates the standard expression (Rosenbluth et al. Reference Rosenbluth, Hazeltine and Hinton1972; Helander & Sigmar Reference Helander and Sigmar2002).

6 Corrections due to tangential drifts

We are now in a position to address our main objective, to calculate the flow velocity of any plasma species that is in the $\sqrt{\unicode[STIX]{x1D708}}$ -regime of collisionality. This is practically always the case for the bulk ions in sufficiently hot stellarator plasmas.

In our analysis above, we have made three approximations. First, the normalised gyroradius $\unicode[STIX]{x1D70C}_{\ast }$ and collisionality $\unicode[STIX]{x1D708}_{\ast }$ have been assumed to be small. Second, it has been assumed that the distribution function is approximately Maxwellian, $f^{+}=F_{0}+F_{1}$ with $F_{1}\ll F_{0}$ . In an omnigenous magnetic field, this would follow from $\unicode[STIX]{x1D70C}_{\ast }\ll 1$ , and in a general stellarator it holds if $\unicode[STIX]{x1D70C}_{\ast }\ll \unicode[STIX]{x1D708}_{\ast }$ , but if $\unicode[STIX]{x1D708}_{\ast }$ is very small it requires that the magnetic field is sufficiently well optimised or that a radial electric field arises to confine the orbits (see the orderings in (2.16) and (2.17)). Our third approximation is that we have, on a couple of occasions, replaced $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+}$ by $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{0}$ in (2.12). This approximation is justified in the $1/\unicode[STIX]{x1D708}$ -regime (see (2.21)) but not in the $\sqrt{\unicode[STIX]{x1D708}}$ -regime. For very low collisionality, the correct ordering is given by (2.23), according to which (2.19) becomes

(6.1) $$\begin{eqnarray}\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}F_{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}+\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}=0.\end{eqnarray}$$

Note that collisions can now be neglected since their frequency is lower than that of the poloidal drift.Footnote ³ Using this equation for $F_{1}$ and employing either ordering (2.16) or ordering (2.17), we find

(6.2) $$\begin{eqnarray}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f^{+}=(\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}})\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}.\end{eqnarray}$$

In the ordering (2.17), the orbit average $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ is negligible compared to $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ , but in the ordering (2.16) these two terms are comparable and $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6F7}^{\prime }(\unicode[STIX]{x1D713})$ is equal to its bounce average. To be as general as possible, we keep $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ . We thus expect that the expressions for the plasma rotation velocity and bootstrap current velocity derived above may not be accurate at collisionalities below that of the $1/\unicode[STIX]{x1D708}$ -regime.

Equations (2.12), (2.23) and (6.2) then imply

(6.3) $$\begin{eqnarray}f^{-}=-\int _{l_{0}}^{l}(\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}})\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\frac{\text{d}l^{\prime }}{v_{\Vert }}+X,\end{eqnarray}$$

where $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ vanishes for circulating particles and $X$ for trapped ones. The odd part of the distribution function has thus acquired a new term, which one expects will tend to partly cancel the term containing $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ . This term vanishes for circulating particles and is small for trapped particles in a nearly omnigenous stellarator, but in general it can make a difference. In particular, if the magnetic geometry is such that the radial drift velocity does not vary much over a bounce orbit, $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\simeq \overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ , then $f^{-}$ (and the plasma flow velocity associated with it) will be considerably smaller than before. This is, for instance, the case in the limit considered in Ho & Kulsrud (Reference Ho and Kulsrud1987) of a stellarator with a large number of local magnetic wells.

Since $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ vanishes for circulating particles, this new term does not affect the analysis leading up to (2.33), but it does modify (2.37), which is replaced by

(6.4) $$\begin{eqnarray}\langle V_{a1}B\rangle =-\frac{1}{n_{a}}\left\langle B^{2}\int _{l_{0}}^{l}\frac{\text{d}l^{\prime }}{B(l^{\prime })}\int (\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}})\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\,\text{d}^{3}v\right\rangle =\frac{T_{a}A_{1a}}{e_{a}}(\langle uB^{2}\rangle +\unicode[STIX]{x1D6FD}),\end{eqnarray}$$

with

(6.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}=\left\langle \frac{3B^{2}}{2}\int _{l_{max}}^{l}\,\text{d}l^{\prime }\int _{1/B_{max}}^{1/B(l^{\prime })}\overline{\unicode[STIX]{x1D709}(\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{\unicode[STIX]{x1D709}}{B}\right)}\frac{\text{d}\unicode[STIX]{x1D706}}{\sqrt{1-\unicode[STIX]{x1D706}B(l^{\prime })}}\right\rangle\end{eqnarray}$$

Whether the new term affects the integration constant $X$ and thus the flow velocity $V_{a2}$ in (2.34) depends on the collision operator. In the case of pure Lorentz scattering, it does not and the current thus becomes equal to

(6.6) $$\begin{eqnarray}\langle BJ_{a\Vert }\rangle =(f_{s}+\langle uB^{2}\rangle +\unicode[STIX]{x1D6FD})p_{a}A_{1a}\end{eqnarray}$$

instead of (3.10).

If the momentum-restoring term is kept in the collision operator (1.1), the integration constant $X$ is still given by (4.7) but is nevertheless affected by the radial electric field through $w_{i}$ . In (4.10), we need to retain $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{1}$ , so that

(6.7) $$\begin{eqnarray}\displaystyle \langle w_{i}B\rangle & = & \displaystyle -\!\frac{1}{n_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\left\langle B^{2}\int _{0}^{\infty }\unicode[STIX]{x1D708}_{D}^{ii}2\unicode[STIX]{x03C0}v^{3}\,\text{d}v\int _{0}^{1/B}\,\text{d}\unicode[STIX]{x1D706}\left(\unicode[STIX]{x1D706}\frac{\unicode[STIX]{x2202}X_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right.\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left.\int _{l_{0}}^{l}(\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}})\frac{\unicode[STIX]{x2202}F_{i0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\frac{\text{d}l^{\prime }}{|v_{\Vert }|}\right)\right\rangle .\end{eqnarray}$$

Evaluating the integrals on the right using (4.7) and (6.5) gives

(6.8) $$\begin{eqnarray}\langle w_{i}B\rangle =\frac{T_{i}(f_{s}+\langle uB^{2}\rangle +\unicode[STIX]{x1D6FD})}{e(1-f_{c})}(A_{1i}-\unicode[STIX]{x1D702}A_{2i}).\end{eqnarray}$$

The piece of the flow $\langle V_{i2}B\rangle$ is still given by (4.17), and the ion flow velocity becomes

(6.9) $$\begin{eqnarray}\langle V_{i\Vert }B\rangle =\frac{T_{i}(f_{s}+\langle uB^{2}\rangle +\unicode[STIX]{x1D6FD})}{e(1-f_{c})}(A_{1i}-\unicode[STIX]{x1D702}f_{c}A_{2i}),\end{eqnarray}$$

instead of (4.18). If the stellarator is well optimised in the sense defined by the ordering (2.17), $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ is small and, as a consequence, the correction $\unicode[STIX]{x1D6FD}$ will be unimportant. However, in other situations, $\unicode[STIX]{x1D6FD}$ will be as large as $f_{s}+\langle uB^{2}\rangle$ , even if the fraction of trapped particles is small. Although the new term containing $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}$ only affects the distribution function in the trapped part of phase space, the response of the circulating particles is modified through the integration constant $X_{i}$ . Plasma rotation is thus expected to be different in the $1/\unicode[STIX]{x1D708}$ - and $\sqrt{\unicode[STIX]{x1D708}}$ -regimes.

Given that the bootstrap current in (5.18) is proportional to $\langle V_{i\Vert }B\rangle$ , there will be a difference in the bootstrap current as well. In particular, the cancellation of the radial electric field does not occur in the $\sqrt{\unicode[STIX]{x1D708}}$ regime, leading to a bootstrap current that depends on $\text{d}\unicode[STIX]{x1D6F7}/\text{d}\unicode[STIX]{x1D713}$ . In plasmas where the electrons are so much hotter than the ions that the radial electric field is determined by the ‘electron root’ of the ambipolarity condition equation, not only the ions, but also the electrons are in the $\sqrt{\unicode[STIX]{x1D708}}$ -regime, and there will be an additional effect of the radial electric field on the bootstrap current.

7 Finite collisionality

In the limit of very small collisionality, the bootstrap current is independent of $\unicode[STIX]{x1D708}_{\ast }$ , but numerical calculations reveal that finite- $\unicode[STIX]{x1D708}_{\ast }$ corrections are important in practice (Beidler et al. Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maassberg, Mikkelsen, Murakami and Schmidt2011). Hinton & Rosenbluth (Reference Hinton and Rosenbluth1973) showed that the first correction to the bootstrap current in a tokamak is proportional to $\unicode[STIX]{x1D708}_{\ast }^{1/2}$ and therefore cannot be calculated by simply continuing the expansion of the drift kinetic equation to higher order in $\unicode[STIX]{x1D708}_{\ast }$ . Instead, a boundary layer forms in phase space between the regions of trapped and circulating orbits, and the correction to the bootstrap current scales as the width of this layer. Essentially the same calculation holds in an omnigenous stellarator (Helander et al. Reference Helander, Geiger and Maassberg2011), partly thanks to the fact that all maxima of the field strength on the same magnetic surface are similar (Cary & Shasharina Reference Cary and Shasharina1997; Parra et al. Reference Parra, Calvo, Helander and Landreman2015). However, the problem is much more complicated in a stellarator of arbitrary geometry, where along a magnetic field line there are local maxima that are arbitrarily close but not equal to the global maximum, $B_{max}$ , making the structure of the boundary layer different from that in an axisymmetric field.

We cannot perform a rigorous calculation like that of Hinton and Rosenbluth for a general stellarator, but we can estimate the width of the boundary layer by considering the statistical properties of the local magnetic field strength maxima. In a tokamak, where all such maxima are equal and equidistant, the estimate is obtained by noting that the effective collision frequency for scattering of $\unicode[STIX]{x1D709}=v_{\Vert }/v$ by an amount $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}$ is $\unicode[STIX]{x1D708}_{eff}\sim \unicode[STIX]{x1D708}/\unicode[STIX]{x0394}\unicode[STIX]{x1D709}^{2}$ . The mean free path for such scattering is thus

(7.1) $$\begin{eqnarray}L_{eff}=R\unicode[STIX]{x0394}\unicode[STIX]{x1D709}^{2}/\unicode[STIX]{x1D708}_{\ast },\end{eqnarray}$$

where

(7.2) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{\ast }=\frac{R\unicode[STIX]{x1D708}}{v_{T}}.\end{eqnarray}$$

The distance between two consecutive bounce points close to the trapped–passing boundary is $2\unicode[STIX]{x03C0}qR$ in a tokamak of major radius $R$ and rotational transform $\unicode[STIX]{x1D704}=1/q$ , and is therefore equal to the mean free path for scattering by an amount

(7.3) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\sim \sqrt{q\unicode[STIX]{x1D708}_{\ast }}.\end{eqnarray}$$

This is the width of the collisional boundary layer around the trapped–passing boundary: particles within this region do not ‘know’ whether they are trapped since collisional scattering across the boundary occurs on a time scale shorter than that of the bounce motion.

In order to investigate the stellarator case, we need to recognise that the local maxima of $B$ along a field line are in general all unequal. To calculate their distribution, we first choose magnetic coordinates on the flux surface we consider such that the maximum field strength is attained in the point $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=(0,0)$ Footnote ⁴ , which for simplicity is assumed to be the only point where $B=B_{max}$ .Footnote ⁵ In its vicinity,

(7.4) $$\begin{eqnarray}B(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\simeq B_{max}+{\textstyle \frac{1}{2}}(B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}\unicode[STIX]{x1D703}^{2}+2B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D711}}\unicode[STIX]{x1D703}\unicode[STIX]{x1D711}+B_{\unicode[STIX]{x1D711}\unicode[STIX]{x1D711}}\unicode[STIX]{x1D711}^{2}),\end{eqnarray}$$

where derivatives are denoted by subscripts. Let us now follow $B(l)$ along a field line emanating from the point $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=(0,0)$ . If the rotational transform $\unicode[STIX]{x1D704}=n/m$ is a rational number, $B$ will reach $B_{max}$ after $m$ toroidal turns, and if $\unicode[STIX]{x1D704}$ is irrational it will never do so. In any case, it will come close to $B=B_{max}$ only when the field line returns to the neighbourhood of its starting point, i.e. when $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$ are close to zero (or multiples of $2\unicode[STIX]{x03C0}$ , which we may subtract). After $m$ toroidal turns around the torus, the value of the poloidal angle is

(7.5) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{m}=2\unicode[STIX]{x03C0}(m\unicode[STIX]{x1D704}-N),\end{eqnarray}$$

where $N$ is the integer closest to $m\unicode[STIX]{x1D704}$ , and if $\unicode[STIX]{x1D704}$ is irrational these values are evenly distributed on the interval $\unicode[STIX]{x1D703}_{m}\in [-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$ , according to a theorem of Hermann Weyl (Weyl Reference Weyl1916; Helander Reference Helander2014). If $\unicode[STIX]{x1D703}_{m}$ is small, the local variation of the field strength along the field line is

(7.6) $$\begin{eqnarray}B(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\simeq B_{max}+{\textstyle \frac{1}{2}}[B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}_{m}+\unicode[STIX]{x1D704}\unicode[STIX]{x1D711})^{2}+2B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D711}}(\unicode[STIX]{x1D703}_{m}+\unicode[STIX]{x1D704}\unicode[STIX]{x1D711})\unicode[STIX]{x1D711}+B_{\unicode[STIX]{x1D711}\unicode[STIX]{x1D711}}\unicode[STIX]{x1D711}^{2}],\end{eqnarray}$$

and its local maximum, which is situated where $\unicode[STIX]{x2202}B/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$ , is approximately

(7.7) $$\begin{eqnarray}B_{max}^{loc}\simeq B_{max}(1-\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D703}_{m}^{2}),\end{eqnarray}$$

with

(7.8) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}B_{max}=\frac{B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D711}}^{2}-B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}B_{\unicode[STIX]{x1D711}\unicode[STIX]{x1D711}}}{2(\unicode[STIX]{x1D704}^{2}B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}+2\unicode[STIX]{x1D704}B_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D711}}+B_{\unicode[STIX]{x1D711}\unicode[STIX]{x1D711}})}.\end{eqnarray}$$

The probability that such a local maximum exceeds $B_{max}-\unicode[STIX]{x0394}B$ is thus equal to

(7.9) $$\begin{eqnarray}P(B_{max}^{loc}>B_{max}-\unicode[STIX]{x0394}B)=\frac{1}{\unicode[STIX]{x03C0}}\sqrt{\frac{\unicode[STIX]{x0394}B}{\unicode[STIX]{x1D6EC}B_{max}}}.\end{eqnarray}$$

On average, one thus expects having to follow the field line a number $N(\unicode[STIX]{x0394}B)=\unicode[STIX]{x03C0}\sqrt{\unicode[STIX]{x1D6EC}B_{max}/\unicode[STIX]{x0394}B}$ turns around the torus, corresponding to a distance

(7.10) $$\begin{eqnarray}L(\unicode[STIX]{x0394}B)\simeq 2\unicode[STIX]{x03C0}^{2}R\sqrt{\frac{\unicode[STIX]{x1D6EC}B_{max}}{\unicode[STIX]{x0394}B}},\end{eqnarray}$$

until $B$ comes within $\unicode[STIX]{x0394}B$ of its maximum value.

If collisional scattering changes $\unicode[STIX]{x1D709}$ by $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}$ , the value of $B$ at the bounce point will change by $\unicode[STIX]{x0394}B/B\sim \unicode[STIX]{x0394}\unicode[STIX]{x1D709}$ . We obtain an estimate for the width of the collisional boundary layer by equating the corresponding mean free path (7.1) to the length (7.10), which gives

(7.11) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D6EC}^{1/5}\unicode[STIX]{x1D708}_{\ast }^{2/5},\end{eqnarray}$$

within a factor of order unity. The parameter $\unicode[STIX]{x1D6EC}$ vanishes for omnigeneous stellarators since the maximum field strength on each magnetic surface is not reached in a discrete number of points but along a curve that closes poloidally, toroidally or helically around the torus (Cary & Shasharina Reference Cary and Shasharina1997; Helander Reference Helander2014; Parra et al. Reference Parra, Calvo, Helander and Landreman2015). For optimised stellarators with small $\unicode[STIX]{x1D6EC}$ , the estimate $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D6EC}^{1/5}\unicode[STIX]{x1D708}_{\ast }^{2/5}$ is only valid as long as $L\gg R$ , that is, when $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\sim \unicode[STIX]{x0394}B/B\ll \unicode[STIX]{x1D6EC}$ . For $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\gtrsim \unicode[STIX]{x1D6EC}$ , $L\sim R$ and we find a boundary layer width $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D708}_{\ast }^{1/2}$ . Then,

(7.12) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D709}\sim \left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D708}_{\ast }^{1/2}\quad & \text{for }\unicode[STIX]{x1D708}_{\ast }\gtrsim \unicode[STIX]{x1D6EC}^{2}\\ \unicode[STIX]{x1D6EC}^{1/5}\unicode[STIX]{x1D708}_{\ast }^{2/5}\quad & \text{for }\unicode[STIX]{x1D708}_{\ast }\lesssim \unicode[STIX]{x1D6EC}^{2}.\end{array}\right.\end{eqnarray}$$

We thus expect the collisional boundary layer to scale slightly differently in a stellarator as compared with the tokamak.

It is difficult to say whether this prediction is borne out by numerical simulations. In many cases, the approach to the low-collisionality asymptote given by (5.25) is so slow that it is impossible to say anything about the scaling, see e.g. figures 22–24 of Beidler et al. (Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maassberg, Mikkelsen, Murakami and Schmidt2011) or figure 3 of Kernbichler et al. (Reference Kernbichler, Kasilov, Kapper, Martitsch, Nemov, Albert and Heyn2016). At least sometimes, the finite-collisionality correction appears to scale as $\unicode[STIX]{x1D708}_{\ast }^{0.5}$ , which is indistinguishable from $\unicode[STIX]{x1D708}_{\ast }^{0.4}$ , see Helander et al. (Reference Helander, Geiger and Maassberg2011).