Stellarator bootstrap current and plasma flow velocity at low collisionality

TL;DR

This paper provides a comprehensive analysis of bootstrap current and plasma flow velocity in low-collisionality stellarators, highlighting the limitations of conventional formulas and proposing necessary modifications across different collisionality regimes.

Contribution

It offers a unified calculation method for bootstrap current and flow velocity in all low-collisionality regimes in stellarators, extending beyond the traditional $1/ u$-regime.

Findings

Conventional formulas are accurate only in the $1/ u$-regime.

Modifications are needed in the $\sqrt{ u}$-regime.

Finite collisionality correction scales as $ u^{2/5}$.

Abstract

The bootstrap current and flow velocity of a low-collisionality stellarator plasma are calculated. As far as possible, the analysis is carried out in a uniform way across all low-collisionality regimes in general stellarator geometry, assuming only that the confinement is good enough that the plasma is approximately in local thermodynamic equilibrium. It is found that conventional expressions for the ion flow speed and bootstrap current in the low-collisionality limit are accurate only in the $1/\nu$-collisionality regime and need to be modified in the $\sqrt{\nu}$-regime. The correction due to finite collisionality is also discussed and is found to scale as $\nu^{2/5}$.

Equations

$$C_a(f_{a1})\to {\nu }_a\left(\mathcal{L}(f_{a1})+\frac{m_a{v}_{\Vert }}{T_a}{w}_a{f}_{a0}\right),$$

$$\mathcal{L}=\frac{1}{2}\frac{\partial }{\partial \xi }\left(1-{\xi }^2\right)\frac{\partial }{\partial \xi }$$

$$n_a{e}_a({\mathbf{V}}_a\times \mathbf{B}-\nabla \Phi )=\nabla p_a,$$

$${\mathbf{V}}_a(\psi ,\alpha ,l)={\omega }_a\frac{\mathbf{b}\times \nabla \psi }{B}+V_{a\Vert }\mathbf{b},$$

$${\omega }_a(\psi )=\frac{1}{n_a{e}_a}\frac{d{p}_a}{d\psi }+\frac{d\Phi }{d\psi }.$$

$$V_{\Vert }(\psi ,\alpha ,l)={\omega }_a(\psi )\left[u(\psi ,\alpha ,l)+k_a(\psi )\right]B,$$

$$\mathbf{b}\cdot \nabla u=-(\mathbf{b}\times \nabla \psi )\cdot \nabla \left(\frac{1}{B^2}\right),$$

$$f_a(\mathbf{r},ϵ,\mu ,\sigma )=f_a^{+}(\mathbf{r},ϵ,\mu )+f_a^{-}(\mathbf{r},ϵ,\mu ),$$

$$f_a^{\pm }(\mathbf{r},ϵ,\mu )=\frac{1}{2}\left[f_a(\mathbf{r},ϵ,\mu ,+1)\pm f_a(\mathbf{r},ϵ,\mu ,-1)\right],$$

$$v_{\Vert }{\nabla }_{\Vert }{f}^{-}=C^{+}(f)-{\mathbf{v}}_d\cdot \nabla f^{+},$$

$$v_{\Vert }{\nabla }_{\Vert }{f}^{+}=C^{-}(f)-{\mathbf{v}}_d\cdot \nabla f^{-},$$

$$C^{+}(f_1)=\nu \mathcal{L}(f_1^{+}),$$

$$C^{-}(f_1)=\nu \left(\mathcal{L}(f_1^{-})+\frac{m{v}_{\Vert }}{T}w{f}_0\right).$$

$$\overline{g}(\psi ,\alpha ,ϵ,\mu )=\frac{1}{{\tau }_b}{\int }_{l_1}^{l_2}{g}^{+}(\psi ,\alpha ,l,ϵ,\mu )\frac{dl}{|v_{\Vert }|},$$

$${\tau }_b={\int }_{l_1}^{l_2}\frac{dl}{|v_{\Vert }|}$$

$$\overline{g}(\psi,\epsilon,\mu,\sigma)=\lim_{L\rightarrow\infty}\int_{0}^{L}g(\psi,\alpha,l,\epsilon,\mu)\frac{dl}{v_{\|}}\bigg{/}\int_{0}^{L}\frac{dl}{v_{\|}}.$$

$$\overline{g}(\psi,\epsilon,\mu,\sigma)=\left\langle\frac{Bg}{v_{\|}}\right\rangle\bigg{/}\left\langle\frac{B}{v_{\|}}\right\rangle.$$

$$\overline{{\mathbf{v}}_d\cdot \nabla f^{+}}=\overline{C^{+}(f)}.$$

$$f^{-}(\psi ,\alpha ,l,ϵ,\mu ,\sigma )={\int }_{l_0}^l\left[C^{+}(f)-{\mathbf{v}}_d\cdot \nabla f^{+}\right]\frac{d{l}^{\prime }}{v_{\Vert }}+X(\psi ,\alpha ,ϵ,\mu ,\sigma ),$$

$$\overline{{\mathbf{v}}_d\cdot \nabla f^{-}}=\overline{C^{-}(f)}.$$

$${\mathbf{v}}_E=-\frac{1}{B}\nabla \Phi \times \mathbf{b}\simeq \frac{1}{B}\frac{d\Phi }{d\psi }\mathbf{b}\times \nabla \psi ,$$

$${\mathbf{v}}_M=\frac{v_{\Vert }^2}{\Omega }\mathbf{b}\times (\mathbf{b}\cdot \nabla \mathbf{b})+\frac{\mu }{m\Omega }\mathbf{b}\times \nabla B,$$

$${\mathbf{v}}_E\gg {\mathbf{v}}_M.$$

$${\mathbf{v}}_E\sim {\mathbf{v}}_M,({\mathbf{v}}_E+{\mathbf{v}}_M)\cdot \nabla \psi \gg \overline{({\mathbf{v}}_E+{\mathbf{v}}_M)\cdot \nabla \psi }.$$

$$f^{+}=F_0+F_1,$$

$$\overline{{\mathbf{v}}_d\cdot \nabla \alpha }\frac{\partial F_1}{\partial \alpha }-\overline{C^{+}(F_1)}=-\overline{{\mathbf{v}}_d\cdot \nabla \psi }\frac{\partial F_0}{\partial \psi }.$$

$$\overline{{\mathbf{v}}_d\cdot \nabla \psi }\frac{\partial F_0}{\partial \psi }=0$$

$${\mathbf{v}}_d\cdot \nabla F_1\ll {\mathbf{v}}_d\cdot \nabla F_0.$$

$$C(F_1)\ll \overline{{\mathbf{v}}_d\cdot \nabla \alpha }(\partial F_1/\partial \alpha ),$$

$$C(F_1)\ll {\mathbf{v}}_d\cdot \nabla F_1\sim {\mathbf{v}}_d\cdot \nabla F_0.$$

Full text

1 Introduction

Calculating the neoclassical transport in stellarators is complicated. Only in the limit of high collisionality, the so-called Pfirsch-Schlüter regime, can an accurate analytical calculation be completed using the full, linearised Landau collision operator [Simakov & Helander (2009); Braun & Helander (2010)]. Most plasmas of interest are, however, in the opposite regime of low collisionality, which is divided into sub-regimes of different character, where the transport is proportional to $1/\nu$, $\sqrt{\nu}$ or $\nu$, in order of decreasing collisionality, where $\nu$ is the collision frequency. (When the radial electric field is sufficiently small, there can be a so-called superbanana-plateau regime at collisionalities between those where the transport scales as $1/\nu$ and $\nu$, but this regime is only important in particular magnetic geometries and will be ignored in the present paper.) In practice, the various collisionality regimes are usually not well separated from each other, and the true transport must be calculated numerically [Beidler et al. (2011)].

An exception occurs, formally at least, in quasisymmetric and omnigeneous stellarators, where all collisionless orbits are perfectly confined.111 Quasisymmetric magnetic fields are those in which the field strength is periodic along field lines, $B(l+L)=B(l)$ where $L$ only depends on the magnetic surface. Omnigenous fields are those where the average radial particle drift vanishes. Other definitions of these concepts and how they relate to each other have been given in many places, including Helander (2014). In such systems, the problem of calculating the neoclassical transport can be reduced to the corresponding problem in an axisymmetric tokamak. In quasisymmetric systems, this is obvious from the fact that the orbits, and hence the drift kinetic equation, are “isomorphic” to those in a tokamak [Boozer (1983); Landreman & Catto (2011)], but even in systems that are merely omnigenous it is, less trivially, possible to split the distribution function into two terms, the first of which can be calculated exactly and the second of which is determined by the same drift kinetic equation as in a tokamak [Helander & Nührenberg (2009); Landreman & Catto (2012)]. It is thus possible to calculate the transport analytically insofar as this is possible in a tokamak, i.e., if the exact collision operator for each species $a$ is replaced by an operator of the form

[TABLE]

where $\nu_{a}$ is the collision frequency,

[TABLE]

the Lorentz operator in terms of $\xi=v_{\|}/v$, $w_{a}$ a function of the speed $v$ determined by integral conditions, such as conservation of momentum, and the other symbols have their usual meanings. This operator is not a bad approximation for low collisionality and large aspect ratio, and the problem of calculating the transport in a perfectly omnigenous stellarator can therefore be considered to be solved. The result is, however, of red limited applicability because the deviation from omnigeneity is, in practice, always large enough that the transport has a clear $1/\nu$-regime at low enough collisionality. The deviations from omnigeneity (or quasisymmetry) have to be extremely small in order for tokamak-like levels of neoclassical transport to prevail at low collisionality [Calvo et al. (2013, 2014, 2015, 2017)]. Even in HSX, the most quasisymmetric of stellarators, there is barely a discernible “banana” regime of tokamak-like transport – see Fig.  14 of Beidler et al. (2011).

In the present paper, we calculate the flow of plasma within, rather than across, magnetic flux surfaces. Our principal aim is to calculate the bootstrap current in the $\sqrt{\nu}$-regime, which has not been done analytically before. The collisionality is assumed to be small, but the geometry is kept as general as possible, allowing for a $1/\nu$-regime. We also assume that the electrostatic potential is a flux function to the order of interest, and that the plasma is close to a radially local thermodynamic equilibrium, that is, the distribution function is close to a stationary Maxwellian with density and temperature that are flux functions. In general, for very-low-collisionality regimes (collisionalities corresponding to the $\sqrt{\nu}$ and $\nu$ regimes), trapped particles will drift radially distances of the order of the minor radius of the device, making it impossible to achieve a radially local thermodynamic equilibrium [Helander (2014); Calvo et al. (2017)]. We thus need extra conditions to be satisfied to be able to treat very low collisionalities, and are aware of two different assumptions in the literature that lead to the required condition of locally radial thermodynamic equilibrium: closeness to an omnigeneous configuration [Calvo et al. (2017)], and the condition that the ${\bf E}\times{\bf B}$ drift is much larger than the $\nabla B$ and curvature drifts (valid for large-aspect-ratio stellarators or regions of very large plasma gradients) [Ho & Kulsrud (1987)]. We try to be as general as possible, and we derive equations that are valid for these two possible assumptions.

For a plasma in radially local thermodynamic equilibrium, each species satisfies the leading-order force-balance relation

[TABLE]

where $n_{a}$ denotes the density, $e_{a}$ the charge, ${\bf V}_{a}$ the flow velocity and $p_{a}=n_{a}T_{a}$ the pressure of the species in question. The magnetic field is in Clebsch form, ${\bf B}=\nabla\psi\times\nabla\alpha$, where $\psi$ labels the magnetic surfaces and $\alpha$ the various field lines on each such surface. The density, temperature and electrostatic potential are to lowest order constant on flux surfaces, which forces the flow velocity to be of the form

[TABLE]

where ${\bf b}={\bf B}/B$, $l$ denotes the arc length along the field, and

[TABLE]

Incompressibility, $\nabla\cdot{\bf V}_{a}=0$, then implies

[TABLE]

where the function $u(\psi,\alpha,l)$ satisfies the magnetic differential equation

[TABLE]

and $k_{a}(\psi)$ denotes an arbitrary integration constant. The kinetic theory developed below gives the flux function $k_{a}(\psi)$.

As we shall see, in the $1/\nu$-regime of collisionality the low-collisionality results of Shaing & Callen (1983) and Shaing-1989 (1989Shaing et al.)Shaing, Carreras, Dominguez, Lynch & Tolliver) are recovered, who derived it using the Hirshman-Sigmar moment method [Hirshman & Sigmar (1981); Helander & Sigmar (2002)]. However, we find that this result ceases to be valid in the $\sqrt{\nu}$ or $\nu$-regimes, where the component of the drifts tangential to the flux surface causes the appearance of an additional term that needs to be taken into account. It is well known from numerical results [Beidler et al. (2011); Helander et al. (2011)] that the approach to the low-collisionality limit is slow, and in the penultimate section of the paper we identify a possible reason for this behaviour and derive a scaling for the deviation from the low-collisionality result. Our conclusions are summarised in the final section.

2 Drift kinetic equation

We begin by dividing the distribution function of each species into even and odd parts in $v_{\|}=\sigma|v_{\|}|$,

[TABLE]

[TABLE]

where $\epsilon=m_{a}v^{2}/2+e_{a}\Phi({\bf r})$ denotes the (kinetic + potential) energy and $\mu=m_{a}v_{\perp}^{2}/(2B)$ the magnetic moment. The drift kinetic equation is thus split into the two equations

[TABLE]

[TABLE]

where ${\bf v}_{d}={\bf v}_{E}+{\bf v}_{M}$ is the sum of the ${\bf E}\times{\bf 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)

[TABLE]

[TABLE]

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 $\psi$ in leading order – an assumption that will be discussed below. We can thus set $\Phi=0$ on the particular flux surface we are considering, but must retain a non-zero derivative $\Phi^{\prime}(\psi)$. Since $\Phi=0$, we can use the velocity-space coordinate $v$ instead of $\epsilon$. As usual, particle orbits can be divided into trapped and untrapped (or circulating) ones depending on whether the parameter $\lambda=\mu/\epsilon$ does, or does not, exceed $1/B_{\rm max}$, where $B_{\rm max}(\psi)$ 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_{\rm max}$ at a discrete set of points on each surface, whereas $B=B_{\rm max}$ along a circle on the inboard side of a tokamak.

The orbit average of an arbitrary function $g(\psi,\alpha,l,\epsilon,\mu,\sigma)$ is defined by

[TABLE]

in the trapped region of phase space, $\lambda B_{\rm max}>1$. Here $l_{1}$ and $l_{2}$ denote consecutive bounce points, $B(\psi,\alpha,l_{1})=B(\psi,\alpha,l_{2})=1/\lambda$, and

[TABLE]

the bounce time. For circulating particles, $\lambda>1/B_{\rm max}$, the orbit average is defined as

[TABLE]

Since the integral extends along a field line many times around the torus, the result is independent of $\alpha$. It can also be written in terms of the flux-surface average $\left\langle\cdots\right\rangle$ as [Helander (2014)]

[TABLE]

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 Eq. (3), which becomes

[TABLE]

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 Eq. (3),

[TABLE]

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 Eq. (4),

[TABLE]

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/\nu$-regime can be treated in general stellarator geometry, but to our knowledge the bootstrap current has not been calculated before in the $\sqrt{\nu}$-regime. In order to do so, we must first carefully go through the conventional $1/\nu$-regime (until Section 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 Eq. (5) is of order $\nu_{\ast}/\rho_{\ast}$, where $\rho_{\ast}=\rho/L$ is the normalised gyroradius and $\nu_{\ast}=\nu 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 $\rho_{\ast}\ll\nu_{\ast}$, collisions dominate and the distribution function is forced to be close to a Maxwellian. However, at lower collisionality, $\nu_{\ast}\ll\rho_{\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 (2014); Calvo et al. (2017)]. We are aware of two limits in which the particles are well confined and confinement is restored. One limit is described by Ho & Kulsrud (1987), who took the ${\bf E}\times{\bf B}$ drift,

[TABLE]

to be much larger than the $\nabla B$ and curvature drifts,

[TABLE]

with $\Omega$ the gyrofrequency, i.e.,

[TABLE]

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. (2017)]. In this limit, the orbit-average of the radial drift is small compared with its local value,

[TABLE]

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 $\Phi(\psi)$ that is constant on flux surfaces to lowest order [Ho & Kulsrud (1987); Calvo et al. (2017)].

Thus, assuming that the plasma is either in the $1/\nu$-regime or follows the orderings (8) or (9), we shall write the solution to Eq. (5) as

[TABLE]

where $F_{0}(\epsilon,\psi)$ 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 $\alpha$ and $l$ for circulating particles. The piece $F_{1}$ is determined by

[TABLE]

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

[TABLE]

in the untrapped part of phase space, there is no “drive” there, and most simple collision operators, including Eq. (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 (2002)], 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 (8) and (9) [Calvo et al. (2017)].

Equation (11), with the approximation (1) for the collision operator, is routinely solved by the GSRAKE code [Beidler & D´Haeseleer (1995)], 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 Eq. (11). In the higher-collisionality limit (corresponding to the $1/\nu$-regime), $C(F_{1})\gg\overline{{\bf v}_{d}\cdot\nabla\alpha}\;(\partial F_{1}/\partial\alpha)$. Thus, $\overline{{\bf v}_{d}\cdot\nabla\alpha}\;(\partial F_{1}/\partial\alpha)\ll\overline{{\bf v}_{d}\cdot\nabla\psi}\,(\partial F_{0}/\partial\psi)$, which in general implies

[TABLE]

In the lower-collisionality limit (corresponding to the $\sqrt{\nu}$-regime),

[TABLE]

leading to

[TABLE]

Sections 3-5 below deal with the $1/\nu$-regime and Section 6 with the $\sqrt{\nu}$-regime.

2.2 Odd part of the distribution function

The odd part of the distribution function is determined by Eqs. (6) - (7) 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 $\lambda$ in the following way:

[TABLE]

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

[TABLE]

where $V(\psi)$ denotes the volume contained by the flux surface labelled by $\psi$, $V^{\prime}(\psi)=dV/d\psi$, and we have used ${\bf v}_{d}=(v_{\|}/\Omega)\nabla\times(v_{\|}{\bf b})$. On the right-hand side of this expression appears

[TABLE]

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

[TABLE]

Using Eq. (10) and the fact that $f^{+}=F_{0}+F_{1}$ does not depend on $\alpha$ or $l$ in the circulating-particle region, Eq. (15) can be rewritten as

[TABLE]

If we write

[TABLE]

then $\nabla_{\|}h=({\bf v}_{d}\cdot\nabla\psi)/v_{\|}$ and

[TABLE]

since $\left\langle B\nabla_{\|}g\right\rangle$ vanishes for any function $g$, see Eq. (21) of Helander (2014). Equation (16) thus reduces to

[TABLE]

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

[TABLE]

and we have used $f^{-}\sim\rho_{\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

[TABLE]

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,

[TABLE]

[TABLE]

Note that the term with the collision operator does not contribute. In the $1/\nu$ regime, according to Eq. (12), ${\bf v}_{d}\cdot\nabla f_{a}^{+}\simeq{\bf v}_{d}\cdot\nabla F_{0}$. If we use this approximation and the expression [e.g., from Landreman & Catto (2012) or Eq. (63) of Helander (2014)]

[TABLE]

where the gradient of $\xi=\sigma(1-\lambda B)^{1/2}$ is taken at constant $\lambda$, we obtain for the first term

[TABLE]

Evaluating the integrals, we find [Helander et al. (2011)]

[TABLE]

where we have recalled Eq. (2) and fixed the integration constant therein by taking $u=0$ at the point on the flux surface where $B=B_{\rm max}$. (In the notation of Nakajima et al. (1989) and Helander et al. (2011), $g_{2}=-uB^{2}$.) In addition, we have defined the thermodynamic forces

[TABLE]

[TABLE]

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

3 Lorentz limit

To solve Eq. (17), 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 Eq. (1). This applies to electrons in a high-$Z$ plasma and is used in the widespread DKES code [Hirshman et al. (1986); van Rij & Hirshman (1989)].

In the coordinates $v$ and $\lambda$, the operator $\mathcal{L}$ becomes

[TABLE]

Our equation (17) then becomes

[TABLE]

and is readily integrated once, demanding regularity at $\lambda=0$,

[TABLE]

We substitute Eq. (6) 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

[TABLE]

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

[TABLE]

where we have introduced the notation [Nakajima et al. (1989)]

[TABLE]

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

[TABLE]

which becomes

[TABLE]

where we have written

[TABLE]

The current along $\bf B$ is thus equal to [Helander et al. (2011)]

[TABLE]

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). For the bulk ions in a plasma without impurities, we thus take

[TABLE]

where

[TABLE]

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

[TABLE]

$\ln\Lambda$ the Coulomb logarithm, $\epsilon_{0}$ the vacuum permittivity, $x=v/v_{Ti}$, $v_{Ti}=\sqrt{2T_{i}/m_{i}}$ and

[TABLE]

[TABLE]

Instead of Eq. (20) we now obtain

[TABLE]

and instead of Eq. (21)

[TABLE]

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

[TABLE]

where we have defined the velocity-space average [Helander & Sigmar (2002)]

[TABLE]

Upon substitution of Eq. (6) we find

[TABLE]

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

[TABLE]

for any collision operator of the form (1).

The first of the two terms on the right of Eq. (25) can be evaluated using Eq. (21) and is equal to

[TABLE]

with

[TABLE]

[TABLE]

The second term on the right of Eq. (25) becomes

[TABLE]

where we have again approximated ${\bf v}_{d}\cdot\nabla f^{+}$ by ${\bf v}_{d}\cdot\nabla F_{0}$. As long as this approximation is allowed, Eq. (25) thus implies

[TABLE]

which together with Eq. (24) 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, $\left\langle V_{i\|}B\right\rangle=\left\langle(V_{i1}+V_{i2})B\right\rangle$. The first term was calculated in Eq. (19) and the second term is

[TABLE]

where we have used Eq. (21). We thus obtain

[TABLE]

Apart from the geometric factor, this result is similar to that in tokamaks and was obtained by Landreman & Catto (2012) for the case of a perfectly omnigenous stellarator. From our analysis, it is clear that the same result is also valid in the $1/\nu$-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 (1983), 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 ${\bf v}_{d}\cdot\nabla f_{i}^{+}$ can indeed be approximated by ${\bf v}_{d}\cdot\nabla F_{i0}$ in Eq. (6). As we shall see below, the correction that arises in the $\sqrt{\nu}$-regime can be significant.

In an axisymmetric magnetic field, ${\bf B}=F(\psi)\nabla\varphi+\iota\nabla\varphi\times\nabla\psi$,

[TABLE]

and

[TABLE]

so that $f_{s}+\left\langle uB^{2}\right\rangle=-(1-f_{c})F/\iota$. The ion flow velocity (28) thus reduces to the usual result [Helander & Sigmar (2002)]

[TABLE]

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 et al. (1972); Helander & Sigmar (2002)], we use a model linearised electron collision operator $C_{e}(f_{e})$ that is more accurate than the operator in Eq. (23). Rosenbluth et al. (1972) showed that the operator (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 (2012)).

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

[TABLE]

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

[TABLE]

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

[TABLE]

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_{\|}C_{ee}\left(L_{q}^{(3/2)}(x^{2})v_{\|}F_{e0}\right)\,d^{3}v$ using the linearised Landau electron-electron collision operator. The first few Sonine polynomials are

[TABLE]

and the corresponding coefficients $K_{pq}$ are

[TABLE]

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

[TABLE]

where $\nu_{D}^{ee}=(3\sqrt{\pi}/4\tau_{ee})[\phi(x)-G(x)]/x^{3}$, and $w_{e}$ is not a constant, but a function of $v$ given by

[TABLE]

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

[TABLE]

where the symmetric matrix $N_{pq}$ is given by

[TABLE]

In Appendix A we show how to calculate this matrix.

Using the collision operators in Eqs. (29), (30) and (32), we find that

[TABLE]

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

[TABLE]

Using

[TABLE]

we find the equations

[TABLE]

Solving Eqs. (52), we find

[TABLE]

where

[TABLE]

With these results, we can red now obtain the electron current. Using Eqs. (18) and (34), we find

[TABLE]

where

[TABLE]

Using Eq. (67), we thus obtain

[TABLE]

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

[TABLE]

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

[TABLE]

and in place of the next equation

[TABLE]

Using the velocity-space averages

[TABLE]

one finds the following system of equations

[TABLE]

and the ion flow

[TABLE]

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

[TABLE]

[TABLE]

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/\nu$-regime. To lowest order in $f_{t}$, this expression coincides with that of Shaing-1989 (1989Shaing et al.)Shaing, Carreras, Dominguez, Lynch & Tolliver), who included two instead of three Sonine polynomials in the expansion.222Otherwise there would presumably have been exact agreement. The flux-surface average of the geometry function in Eq. (75b) of Shaing & Callen (1983) which is needed to compare the results can be written in our notation as $\langle\tilde{G}_{b}\rangle=G-\langle uB^{2}\rangle$, where $G=-[f_{c}\langle uB^{2}\rangle+f_{s}]/(1-f_{c})$. It also agrees with the result of Landreman & Catto (2012), which it generalises to non-omnigenous stellarators in the $1/\nu$-regime.

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

[TABLE]

which closely approximates the standard expression [Rosenbluth et al. (1972); Helander & Sigmar (2002)].

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{\nu}$-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 $\rho_{\ast}$ and collisionality $\nu_{\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 $\rho_{\ast}\ll 1$, and in a general stellarator it holds if $\rho_{\ast}\ll\nu_{\ast}$, but if $\nu_{\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 Eqs. (8) and (9)]. Our third approximation is that we have, on a couple of occasions, replaced ${\bf v}_{d}\cdot\nabla f^{+}$ by ${\bf v}_{d}\cdot\nabla F_{0}$ in Eq. (6). This approximation is justified in the $1/\nu$-regime [see Eq. (12)] but not in the $\sqrt{\nu}$-regime. For very low collisionality, the correct ordering is given by Eq. (13), according to which Eq. (11) becomes

[TABLE]

Note that collisions can now be neglected since their frequency is lower than that of the poloidal drift.333The collision operator is however important in regions where $\overline{{\bf v}_{d}\cdot\nabla\alpha}=0$ and in thin boundary layers located at the boundaries between trapped and passing particles and at boundaries between particles trapped in different magnetic wells. The width of this boundary layer is of order $\Delta\xi\sim(\nu/|\Phi^{\prime}(\psi)|)^{1/2}$ Ho & Kulsrud (1987); Helander (2014); Calvo et al. (2017). Using this equation for $F_{1}$ and employing either ordering (8) or ordering (9), we find

[TABLE]

In the ordering (9), the orbit average $\overline{{\bf v}_{d}\cdot\nabla\psi}$ is negligible compared to ${\bf v}_{d}\cdot\nabla\psi$, but in the ordering (8) these two terms are comparable and ${\bf v}_{d}\cdot\nabla\alpha=\Phi^{\prime}(\psi)$ is equal to its bounce average. To be as general as possible, we keep $\overline{{\bf v}_{d}\cdot\nabla\psi}$. 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/\nu$-regime.

Equations (6), (13) and (72) then imply

[TABLE]

where $\overline{{\bf v}_{d}\cdot\nabla\psi}$ 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 ${\bf v}_{d}\cdot\nabla F_{0}$. 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 magnetic geometry is such that the radial drift velocity does not vary much over a bounce orbit, ${\bf v}_{d}\cdot\nabla\psi\simeq\overline{{\bf v}_{d}\cdot\nabla\psi}$, 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 (1987)] of a stellarator with a large number of local magnetic wells.

Since $\overline{{\bf v}_{d}\cdot\nabla\psi}$ vanishes for circulating particles, this new term does not affect the analysis leading up to (17), but it does modify Eq. (19), which is replaced by

[TABLE]

with

[TABLE]

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

[TABLE]

instead of Eq. (22).

If the momentum-restoring term is kept in the collision operator (1), the integration constant $X$ is still given by Eq. (24) but is nevertheless affected by the radial electric field through $w_{i}$. In Eq. (25), we need to retain ${\bf v}_{d}\cdot\nabla F_{1}$, so that

[TABLE]

Evaluating the integrals on the right using (24) and (73) gives

[TABLE]

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

[TABLE]

instead of Eq. (28). If the stellarator is well optimised in the sense defined by the ordering (9), $\overline{{\bf v}_{d}\cdot\nabla\psi}$ is small and, as a consequence, the correction $\beta$ will be unimportant. However, in other situations, $\beta$ will be as large as $f_{s}+\left\langle uB^{2}\right\rangle$, even if the fraction of trapped particles is small. Although the new term containing $\overline{{\bf v}_{d}\cdot\nabla\psi}$ 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/\nu$- and $\sqrt{\nu}$-regimes.

Given that the bootstrap current in Eq. (68) is proportional to $\left\langle V_{i\|}B\right\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{\nu}$ regime, leading to a bootstrap current that depends on $d\Phi/d\psi$. 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{\nu}$-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 $\nu_{\ast}$, but numerical calculations reveal that finite-$\nu_{\ast}$ corrections are important in practice [Beidler et al. (2011)]. Hinton & Rosenbluth (1973) showed that the first correction to the bootstrap current in a tokamak is proportional to $\nu_{\ast}^{1/2}$ and therefore cannot be calculated by simply continuing the expansion of the drift kinetic equation to higher order in $\nu_{\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. (2011)], partly thanks to the fact that all maxima of the field strength on the same magnetic surface are similar [Cary & Shasharina (1997); Parra et al. (2015)]. 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_{\rm 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 $\xi=v_{\|}/v$ by an amount $\Delta\xi$ is $\nu_{\rm eff}\sim\nu/\Delta\xi^{2}$. The mean free path for such scattering is thus

[TABLE]

where

[TABLE]

The distance between two consecutive bounce points close to the trapped-passing boundary is $2\pi qR$ in a tokamak of major radius $R$ and rotational transform $\iota=1/q$, and is therefore equal to the mean free path for scattering by an amount

[TABLE]

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 $(\theta,\varphi)=(0,0)$444Here $\theta$ and $\varphi$ denote field-aligned magnetic coordinates whose origin is chosen to coincide with the magnetic-field maximum., which for simplicity is assumed to be the only point where $B=B_{\rm max}$.555In some stellarators, the maximum field strength is attained in an even number of points, which are stellarator-symmetrically located above and below the midplane. This should not affect the basic scaling derived here. In its vicinity,

[TABLE]

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

[TABLE]

where $N$ is the integer closest to $m\iota$, and if $\iota$ is irrational these values are evenly distributed on the interval $\theta_{m}\in[-\pi,\pi]$, according to a theorem of Hermann Weyl [Weyl (1916); Helander (2014)]. If $\theta_{m}$ is small, the local variation of the field strength along the field line is

[TABLE]

and its local maximum, which is situated where $\partial B/\partial\varphi=0$, is approximately

[TABLE]

with

[TABLE]

The probability that such a local maximum exceeds $B_{\rm max}-\Delta B$ is thus equal to

[TABLE]

On average, one thus expects having to follow the field line a number $N(\Delta B)=\pi\sqrt{\Lambda B_{\rm max}/\Delta B}$ turns around the torus, corresponding to a distance

[TABLE]

until $B$ comes within $\Delta B$ of its maximum value.

If collisional scattering changes $\xi$ by $\Delta\xi$, the value of $B$ at the bounce point will change by $\Delta B/B\sim\Delta\xi$. We obtain an estimate for the width of the collisional boundary layer by equating the corresponding mean free path (75) to the length (76), which gives

[TABLE]

within a factor of order unity. The parameter $\Lambda$ 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 (1997); Helander (2014); Parra et al. (2015)]. For optimized stellarators with small $\Lambda$, the estimate $\Delta\xi\sim\Lambda^{1/5}\nu_{\ast}^{2/5}$ is only valid as long as $L\gg R$, that is, when $\Delta\xi\sim\Delta B/B\ll\Lambda$. For $\Delta\xi\gtrsim\Lambda$, $L\sim R$ and we find a boundary layer width $\Delta\xi\sim\nu_{\ast}^{1/2}$. Then,

[TABLE]

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 Eq. (71) is so slow that it is impossible to say anything about the scaling, see e.g. Figs. 22-24 of Beidler et al. (2011) or Fig. 3 of Kernbichler et al. (2016). At least sometimes, the finite-collisionality correction appears to scale as $\nu_{\ast}^{0.5}$, which is indistinguishable from $\nu_{\ast}^{0.4}$, see Helander et al. (2011).

8 Conclusions

In the present paper, the neoclassical theory of parallel transport in stellarators at low collisionality has been reconsidered. For the $1/\nu$-regime, we have recovered a well-known expression (28) for the ion flow (plasma rotation) and derived a novel explicit formula (68) for the bootstrap current. Both these results, in particular the former, are however found to be sensitive to the influence of tangential drifts and finite collisionality. Tangential drifts are found to affect the parallel ion flow significantly in the $\sqrt{\nu}$-regime, because they reduce the radial excursion of trapped particle orbits from flux surfaces, and the response of the trapped-particle population to the radial pressure and temperature gradients is therefore weakened. The average flow velocity of the trapped particles is thus affected by the electric field, and this change is transmitted to the passing population by collisions. As a result, the plasma rotation velocity changes significantly.

A finite collision frequency affects the boundary layer between the trapped and passing particles slightly differently from that in a tokamak. The width of this boundary layer is set by the requirement that the effective collisional scattering frequency across the layer should be comparable to the bounce frequency between two local magnetic maxima along the same magnetic field line that are high enough to matter for the boundary layer. Whether they do so depends on the width of this layer, and therefore on the collisionality. In a stellarator, the distance along the magnetic field between two relevant field-strength maxima therefore depends on the collision frequency. As a result, the width of the boundary layer is predicted to scale as $\nu_{\ast}^{2/5}$ rather than $\nu_{\ast}^{1/2}$.

Acknowledgment

The first author gratefully acknowledges the splendid hospitality of Oxford University, in particular Merton College, where this work was initiated. We are grateful to Craig Beidler for many interesting discussions and a number of helpful comments on the manuscript. SLN was supported by the Framework grant for Strategic Energy Research (Dnr. 2014-5392) from Vetenskapsrådet.

Appendix A Model collision operator for electron-electron collisions

To calculate the electron distribution function, we need a model linearised electron-electron collision operator $C_{ee}(f_{e}^{-})$ that contains the exact pitch-angle scattering piece of the Landau electron-electron collision operator, that conserves momentum, that is self-adjoint, that has $f_{e}^{-}=(m_{e}v_{\|}V_{e\|}/T_{e})F_{e0}$ as the only solution to $C_{ee}(f_{e}^{-})=0$, and that satisfies (31).

Following Taguchi (1988) and Helander & Sigmar (2002), we propose the model collision operator

[TABLE]

with the coefficients $a_{r}$ determined by the system of equations

[TABLE]

The matrix $N_{pq}$ is determined using conditions (31). For our model collision operator, the left side of Eq. (31) becomes

[TABLE]

Using Eq. (77), we find

[TABLE]

where $(N^{-1})_{rs}$ is the inverse of the matrix $N_{pq}$. Substituting this result into Eq. (31), we find a linear system of equations that determines $(N^{-1})_{rs}$ and hence $N_{pq}$,

[TABLE]

Note that this equation gives a symmetric matrix $N_{pq}=N_{qp}$. We have used Eq. (78) and

[TABLE]

to determine the coefficients of Eq. (33).

The collision operator has been completely determined using condition (31). It is trivial to check that it contains the pitch-angle scattering piece of the Landau collision operator. We proceed to show that it satisfies the other properties postulated above.

A.1 Self-adjointness

The self-adjointness of the collision operator is derived from the fact that $N_{pq}$ is symmetric. Indeed,

[TABLE]

Using (77), we can rewrite this equation as

[TABLE]

proving self-adjointness.

A.2 Conservation of momentum

The coefficients $N_{pq}$ ensure conservation of momentum. Using the property $K_{p0}=K_{0p}=0$ in Eq. (78) gives

[TABLE]

This is a linear system of equations for $\sum_{s=0}^{S}(N^{-1})_{rs}\tau_{ee}\{\nu_{D}^{ee}L_{s}^{(3/2)}\}$ whose solution is

[TABLE]

where $\delta_{pq}$ is the Kronecker delta. Multiplying this equation by $N_{pr}$ and summing over $r$ finally gives

[TABLE]

With this result, we can prove conservation of momentum. The parallel velocity moment of the collision operator is

[TABLE]

or, using Eq. (77),

[TABLE]

Since the coefficients $N_{p0}$ have the property (79), the collision operator conserves momentum.

A.3 Solution to $C_{ee}(f_{e}^{-})=0$

The solution to $C_{ee}(f_{e}^{-})=0$ is

[TABLE]

(recall that $f_{e}^{-}$ is odd in $v_{\|}$ and that it does not depend on gyrophase). We determine the coefficients $a_{p}$ by using Eq. (77). For this particular form of $f_{e}^{-}$, Eq. (77) becomes

[TABLE]

This system of equations implies that $a_{p}=0$ for $p\neq 0$, and that $a_{0}$ is a free parameter since the coefficients $N_{p0}$ satisfy condition (79). Hence $f_{e}^{-}=(m_{e}v_{\|}a_{0}/T_{e})F_{e0}$ is the only solution to $C_{ee}(f_{e}^{-})=0$.