Impact of main ion pressure anisotropy on stellarator impurity transport

TL;DR

This paper investigates how main ion pressure anisotropy affects impurity transport in stellarators, revealing that previously neglected terms related to pressure anisotropy significantly influence impurity dynamics.

Contribution

It introduces the inclusion of next-order impurity-ion collision terms related to ion pressure anisotropy, providing new analytical expressions for impurity flux in stellarator regimes.

Findings

Pressure anisotropy drives impurity transport differently from parallel friction.

Analytical expressions for impurity flux are derived for multiple collisional regimes.

Numerical evaluations show the importance of pressure anisotropy in realistic stellarator plasmas.

Abstract

Main ions influence impurity dynamics through a variety of mechanisms; in particular, via impurity-ion collisions. To lowest order in an expansion in the main ion mass over the impurity mass, the impurity-ion collision operator only depends on the component of the main ion distribution that is odd in the parallel velocity. These lowest order terms give the parallel friction of the impurities with the main ions, which is typically assumed to be the main cause of collisional impurity transport. Next-order terms in the mass ratio expansion of the impurity-ion collision operator, proportional to the component of the main ion distribution that is even in the parallel velocity, are usually neglected. However, in stellarators, the even component of the main ion distribution can be very large. In this article, such next-order terms in the mass ratio expansion of the impurity-ion collision operator are retained, and analytical expressions for the neoclassical radial flux of trace impurities are calculated in the Pfirsch-Schl\"uter, plateau and $1/\nu$ regimes. The new terms provide a drive for impurity transport that is physically very different from parallel friction: they are associated to anisotropy in the pressure of the main ions, which translates into impurity pressure anisotropy. It is argued that main ion pressure anisotropy must be taken into account for a correct description of impurity transport in certain realistic stellarator plasmas. Examples are given by numerically evaluating the analytical expressions for the impurity flux.

Equations

$$\mathcal{E}=\frac{v^2}{2}+\frac{Z_{z}e}{m_z}\phi ,$$

$$\mu =\frac{v_{\perp }^2}{2B},$$

$$v_{||}=\sigma \sqrt{2\left(\mathcal{E}-\mu B-\frac{Z_{z}e}{m_z}\phi \right)},$$

$$\int (\cdot )\text{mbox}{d}^{3}v\equiv \sum _{\sigma }{\int }_{Z_{z}e\phi /m_z}^{\infty }\text{mbox}d\mathcal{E}{\int }_0^{B^{-1}(\mathcal{E}-Z_{z}e\phi /m_z)}\text{mbox}d\mu {\int }_0^{2\pi }\text{mbox}d\varphi \frac{B}{|v_{||}|}(\cdot ).$$

$$\phi (r,\theta ,\zeta )={\phi }_0(r)+{\phi }_1(r,\theta ,\zeta ),$$

$$⟨f⟩(r)=V^{\prime }(r{)}^{-1}{\int }_0^{2\pi }\text{mbox}d\theta {\int }_0^{2\pi }\text{mbox}d\zeta \sqrt{g_{{}_V}}f(r,\theta ,\zeta ),$$

$$V(r_0)={\int }_0^{r_0}\text{mbox}dr{\int }_0^{2\pi }\text{mbox}d\theta {\int }_0^{2\pi }\text{mbox}d\zeta \sqrt{g_{{}_V}}$$

$$\frac{Z_z^2{n}_z}{Z_i^2{n}_i}\sqrt{\frac{m_z}{m_i}}\ll 1,$$

$${\nu }_{zi}=\frac{Z_z^2{Z}_i^2{e}^4{n}_i{m}_i^{1/2}\ln \Lambda }{6\sqrt{2}{\pi }^{3/2}{ϵ}_0^2{m}_z{T}^{3/2}}$$

$$v_{||}\hat{\mathbf{b}}\cdot \nabla F_z+{\mathbf{v}}_d\cdot \nabla F_z=C_{zi}[F_z,F_i],$$

$${\mathbf{v}}_M=\frac{v_{||}^2}{{\Omega }_z}\hat{\mathbf{b}}\times (\hat{\mathbf{b}}\cdot \nabla \hat{\mathbf{b}})+\frac{\mu }{{\Omega }_z}\hat{\mathbf{b}}\times \nabla B$$

$${\mathbf{v}}_E=\frac{1}{B}\hat{\mathbf{b}}\times \nabla ({\phi }_0+{\phi }_1).$$

$$\frac{|\mathbf{E}|}{B{v}_{tz}}\frac{1}{ϵ}\sim \sqrt{\frac{m_z}{m_i}}\frac{{\rho }_{i\ast }}{{ϵ}^2},$$

$$F_z=F_{z0}+F_{z1}+\dots ,$$

$$v_{||}\hat{\mathbf{b}}\cdot \nabla F_{z0}=C_{zi}[F_{z0},F_{Mi}],$$

$$F_{Mi}=n_i(r){\left(\frac{m_i}{2\pi T(r)}\right)}^{3/2}\exp \left(-\frac{m_i{v}^2}{2T(r)}\right).$$

$$C_{zi}[F_{z0},F_{Mi}]={\nu }_{zi}\frac{T}{m_z}{\nabla }_v\cdot \left({\nabla }_v{F}_{z0}+\frac{m_z\mathbf{v}}{T}{F}_{z0}\right).$$

$$-⟨\exp \left(\frac{Z_{z}e({\phi }_0+{\phi }_1)}{T}\right)\int \exp \left(-\frac{m_z{v}^2}{2T}\right){\left|{\nabla }_v\left(\exp \left(\frac{m_z{v}^2}{2T}\right)F_{z0}\right)\right|}^2\text{mbox}{d}^{3}v⟩=0,$$

$${\nabla }_v\left(\exp \left(\frac{m_z{v}^2}{2T}\right)F_{z0}\right)=0.$$

$$v_{||}\hat{\mathbf{b}}\cdot \nabla F_{Mz}=0,$$

$$F_{Mz}(r,\mathcal{E})=\eta (r){\left(\frac{m_z}{2\pi T(r)}\right)}^{3/2}\exp \left(\frac{Z_{z}e{\phi }_0(r)}{T(r)}\right)\exp \left(-\frac{m_z\mathcal{E}}{T(r)}\right).$$

$$n_z:=\int F_{Mz}\text{mbox}{d}^{3}v,$$

$$n_z(r,\theta ,\zeta )=\eta (r)\exp \left(-\frac{Z_{z}e{\phi }_1(r,\theta ,\zeta )}{T(r)}\right).$$

$$v_{||}\hat{\mathbf{b}}\cdot \nabla F_{z1}+{\mathbf{v}}_d\cdot \nabla r{\rm Y}{F}_{Mz}=C_{zi}^{(\ell )}[F_{z1};h_i].$$

$${\rm Y}=\frac{{\eta }^{\prime }}{\eta }+\frac{Z_{z}e{\phi }_0^{\prime }}{T}+\frac{T^{\prime }}{T}\left(\frac{m_z\mathcal{E}}{T}-\frac{Z_{z}e{\phi }_0}{T}-\frac{3}{2}\right)$$

$${\mathbf{v}}_d\cdot \nabla r=\frac{v_{||}}{{\Omega }_z}\nabla \cdot \left(v_{||}\hat{\mathbf{b}}\times \nabla r\right)$$

$$\displaystyle C^{(\ell)}_{zi}[F_{z1};h_{i}]=\nu_{zi}{\cal K}F_{z1}+\nu_{zi}\frac{m_{z}}{T}\left(Av_{||}-\tr(\buildrel\leftrightarrow\over{\mathbf{M}})+\frac{m_{z}}{T}\buildrel\leftrightarrow\over{\mathbf{M}}\,:\mathbf{v}\mathbf{v}\right)F_{Mz},$$

$$\mathcal{K}{F}_{z1}:=\frac{T}{m_z}{\nabla }_v\cdot \left(F_{Mz}{\nabla }_v\left(\frac{F_{z1}}{F_{Mz}}\right)\right).$$

$$A=\frac{3\sqrt{\pi }{T}^{3/2}}{\sqrt{2}{n}_i{m}_i^{3/2}}\int \frac{v_{||}}{v^3}{h}_i(\mathbf{v})\text{mbox}{d}^{3}v,$$

$$\frac{\text{buildrel}\leftrightarrow }{\mathbf{M}=-\frac{3\sqrt{\pi }{T}^{5/2}}{2\sqrt{2}{n}_i{m}_z{m}_i^{3/2}}\int {\nabla }_{v^{\prime }}{\nabla }_{v^{\prime }}{\nabla }_{v^{\prime }}{v}^{\prime }\cdot {\nabla }_{v^{\prime }}\left(\frac{h_i({\mathbf{v}}^{\prime })}{F_{Mi}}\right)F_{Mi}\text{mbox}{d}^3{v}^{\prime }.}$$

Full text

March 17, 2024

1 Introduction

Impurity transport is one of the most active research areas in the stellarator community because, as a general rule, impurities are observed to accumulate in the plasma core of these devices [1]. The accumulation of heavy, highly-charged impurities is extremely harmful for confinement, as it leads to fuel dilution [2] and to intolerable energy losses due to radiation [1]. Hence, understanding and controlling impurity transport is essential for the success of the stellarator concept as an alternative to tokamaks in the design of future fusion power plants.

In recent years, remarkable progress has been made on the theory of stellarator impurity transport, both analytical [3, 4, 5, 6] and numerical [7, 8, 9, 10, 11, 12]. In [3], the radial impurity flux is calculated analytically when all the species in the plasma are collisional. In [4], the calculation is given for the case in which impurities are collisional and the main ions have low collisionality. A considerable part of the recent theoretical effort (see [5, 6, 7, 8, 9, 10, 11, 12]) has been related to the extension of neoclassical theory and, consequently, neoclassical codes, to include the effect of the component of the electric field that is tangent to magnetic surfaces in the calculation of radial impurity fluxes. The determination of the radial and tangential components of the electric field is one of the neoclassical mechanisms by which main ion dynamics affects impurity transport. In this article, we will focus on the other obvious mechanism: impurity-ion collisions.

Specifically, we would like to incorporate and quantify a collisional effect usually neglected in the analytical description of neoclassical impurity transport in stellarators: the influence of the pressure anisotropy of the main ions on the dynamics of the impurities via collisions (numerical evidence of the relevance of this effect was given in [13] in neoclassical simulations including the exact Landau collision operator). Let us denote by $m_{i}$ and $m_{z}$ the mass of the main ions and the impurities, respectively, and by $h_{i}$ the deviation of the main ion distribution from a Maxwellian distribution. In order to describe collisions between heavy impurities and main ions, the impurity-ion collision operator is expanded in $\sqrt{m_{i}/m_{z}}\ll 1$ and typically approximated by the lowest-order term, which only depends on the component of $h_{i}$ that is odd in the parallel velocity, $h_{i-}$. Next order terms in the mass ratio expansion depend on the component of $h_{i}$ that is even in the parallel velocity, $h_{i+}$. In tokamaks $h_{i+}$ is small, especially at low main ion collisionality (in this paper, we always assume that the main ions have low collisionality, since this is the situation of interest for fusion purposes), so next-order terms in the $\sqrt{m_{i}/m_{z}}\ll 1$ expansion can be safely dropped. However, in stellarators, $h_{i+}$ can be very large in low collisionality plasmas [14, 15, 16] due to the large main ion pressure anisotropy, which translates into significant impurity pressure anisotropy that in turn can drive non-negligible radial impurity transport. We present an explicit calculation showing this for impurities in the Pfirsch-Schlüter, plateau and $1/\nu$ regimes. In each of these regimes, we will estimate the region of parameter space in which the effect of main ion pressure anisotropy on the impurity flux is relevant.

The rest of the paper is organized as follows. In Section 2, we explain the notation and orderings assumed along the paper. We also introduce the trace impurity drift-kinetic equation, including an explicit expression for the impurity-ion collision operator that keeps next-to-lowest-order terms in the $\sqrt{m_{i}/m_{z}}\ll 1$ expansion, and hence the effect of $h_{i+}$. In Section 3, we give a detailed calculation of the impurity flux when the impurities are in the Pfirsch-Schlüter regime. In Section 4, the case of impurities in the plateau regime is worked out. In Section 5, the neoclassical impurity flux is estimated assuming that the impurities are in the $1/\nu$ regime. In Section 6, we evaluate numerically the expressions derived for the impurity flux and illustrate the effect of main ion pressure anisotropy on neoclassical impurity transport in realistic stellarator plasmas. Section 7 contains the conclusions.

2 Impurity drift-kinetic equation including main ion pressure anisotropy effects

In this section we define notation and explain the orderings assumed in the subsequent calculations. We also present the general form of the impurity-drift kinetic equation that will be solved in different regimes in the next sections. In particular, we provide an explicit expression for the impurity-ion collision operator expanded in $\sqrt{m_{i}/m_{z}}\ll 1$ to sufficiently high order that it includes the influence of the pressure anisotropy of the main ions.

2.1 Phase-space coordinates

Before introducing the drift-kinetic equation, we discuss the coordinates that we will employ to locate points on phase space. As velocity coordinates, we take the total energy per mass unit

[TABLE]

the magnetic moment

[TABLE]

the sign of the parallel velocity $\sigma=v_{||}/|v_{||}|=\pm 1$, with

[TABLE]

and the gyrophase, $\phi$. Here, $v$ is the magnitude of the velocity $\mathbf{v}$, $Z_{z}e$ is the impurity charge, $e$ is the proton charge, $\varphi$ is the electrostatic potential, $B$ is the magnitude of the magnetic field $\mathbf{B}$ and $v_{\perp}^{2}=v^{2}-v_{||}^{2}$. The explicit expression of velocity integrals in these coordinates reads

[TABLE]

As space coordinates we take $r$, $\theta$ and $\zeta$, where $r$, the local minor radius, is a flux surface label, and $\theta$ and $\zeta$ are poloidal and toroidal angles, respectively. We write the electrostatic potential $\varphi$ as

[TABLE]

where we assume $e\varphi_{0}/T\sim 1$ and $|\varphi_{1}|\ll|\varphi_{0}|$ (below, we will be more precise about the ordering assumed for the size of $\varphi_{1}$). The temperatures of the main ions and of the impurities are assumed to be equal and constant on flux surfaces, and denoted by $T(r)$.

In these coordinates, the flux-surface average of any function $f(r,\theta,\zeta)$ is defined by

[TABLE]

where $\sqrt{g_{{}_{V}}}=[(\nabla{r}\times\nabla\theta)\cdot\nabla\zeta]^{-1}>0$ is the volume element and

[TABLE]

is the volume enclosed by the surface labeled by ${r}_{0}$. In this paper, primes stand for derivatives with respect to $r$.

2.2 Impurity drift-kinetic equation

Throughout this paper we use the trace impurity approximation,

[TABLE]

where $Z_{i}$ is the main ion charge number and $n_{i}(r)$ is the lowest-order density of the main ions, which is a flux function, whereas the lowest order impurity density, $n_{z}$, is not a flux function if $Z_{z}e\varphi_{1}/T$ is large enough (see below). When (8) is satisfied, impurity-impurity collisions can be neglected with respect to impurity-ion collisions (impurity-electron collisions can be neglected against impurity-ion collisions due to the small mass of the electrons relative to the mass of the main ions).

We assume that the characteristic impurity gyroradius is much smaller than the characteristic size of the stellarator, $\rho_{z*}:=v_{tz}/\Omega_{z}R_{0}\ll 1$, that the $\mathbf{E}\times\mathbf{B}$ drift, with $\mathbf{E}=-\nabla\varphi$, is small compared to the impurity thermal speed (see below a more detailed discussion on this approximation) and that the typical collision frequency between impurities and ions and the characteristic parallel streaming time are comparable, $\nu_{zi*}:=R_{0}\nu_{zi}/v_{tz}\sim 1$ (we will perform a subsidiary expansion in small and large $\nu_{zi*}$ in what follows). Here $v_{tz}=\sqrt{2T_{z}/m_{z}}$ is the impurity thermal speed, $\Omega_{z}=Z_{z}eB/m_{z}$ is the impurity gyrofrequency, $R_{0}$ is the stellarator major radius,

[TABLE]

is the impurity-ion collision frequency, $\epsilon_{0}$ is the vacuum permittivity and $\ln\Lambda$ is the Coulomb logarithm.

Using these assumptions, the impurity distribution $F_{z}(r,\theta,\zeta,{\cal E},\mu,\sigma)$ can be shown to be independent of the gyrophase $\phi$ to the relevant order [17], and the particle motion can be split into fast thermal motion along magnetic field lines and slow drifts across, giving the trace impurity drift-kinetic equation

[TABLE]

where $C_{zi}[F_{z},F_{i}]$ is the impurity-ion collision term, $F_{i}$ is the main ion distribution and $\mathbf{v}_{d}=\mathbf{v}_{M}+\mathbf{v}_{E}$ is the sum of the magnetic drift

[TABLE]

and the $\mathbf{E}\times\mathbf{B}$ drift

[TABLE]

Of all our assumptions, the small $\mathbf{E}\times\mathbf{B}$ drift approximation is the one that is more likely to fail for highly charged impurities, because for a typical electric field $|\mathbf{E}|\sim|\nabla\varphi_{0}|\sim T/(ea)$, with $a$ the stellarator minor radius, the ratio of the term $B^{-1}(\hat{\mathbf{b}}\times\nabla\varphi_{0})\cdot\nabla F_{z}$ to the parallel streaming term $v_{||}\hat{\mathbf{b}}\cdot\nabla F_{z}$ has a size

[TABLE]

where $\rho_{i*}=v_{ti}/(\Omega_{i}R_{0})$ is the normalized gyroradius of the main ions, $v_{ti}=\sqrt{2T/m_{i}}$ is the main ion thermal speed, $\Omega_{i}=Z_{i}eB/m_{i}$ is the main ion gyrofrequency and $\epsilon=a/R_{0}$ is the stellarator inverse aspect ratio. If $|\mathbf{E}|/(\epsilon Bv_{tz})$ is of order unity, the $\mathbf{E}\times\mathbf{B}$ drift term becomes comparable to the parallel streaming term, and the perpendicular drifts are modified by Coriolis and centrifugal forces [18]. Thus, our expression will only be valid in the cases when the estimate (13) is small, which, for sufficiently heavy impurities, amounts to requiring small radial electric field. We do not consider this a major constraint for the main result of this article because, as we will see, the main ion pressure anisotropy contribution becomes important for small radial electric fields.

Let us expand $F_{z}$ in the small parameter $\rho_{z*}$,

[TABLE]

with $F_{z1}\sim\rho_{z*}F_{z0}$. The lowest order terms in (10) give

[TABLE]

where we are assuming that, to lowest order in the expansion, the main ion distribution is a Maxwellian distribution

[TABLE]

A contains a detailed calculation of the lowest order terms of the linearized impurity-ion collision operator expanded in $\sqrt{m_{i}/m_{z}}$. From the results in A, one can infer (as a particular case) that the collisions of heavy impurities with a Maxwellian background of main ions are given by

[TABLE]

Multiplying (15) (with the collision operator given in (17)) by $\exp(m_{z}{\cal E}/T(r))F_{z0}$, integrating over velocities and flux-surface averaging, we get

[TABLE]

which implies

[TABLE]

That is, $F_{z0}=F_{Mz}$, where $F_{Mz}$ is a Maxwellian distribution with the same temperature as $F_{Mi}$. Going back to (15) and using that $C_{zi}[F_{Mz},F_{Mi}]=0$, we obtain

[TABLE]

which, on an ergodic surface, implies

[TABLE]

The flux function $\eta(r)$ is related to the density of the Maxwellian distribution,

[TABLE]

by the relation

[TABLE]

The next-order terms of the drift-kinetic equation (10) in the $\rho_{z*}$ expansion give an equation for $F_{z1}$,

[TABLE]

Here,

[TABLE]

is a combination of the thermodynamic forces, $h_{i}$ is the non-adiabatic component of the deviation of the main ion distribution from a Maxwellian distribution,

[TABLE]

is the radial drift, $\mathbf{v}_{d}\cdot\nabla{r}=\mathbf{v}_{M}\cdot\nabla{r}+\mathbf{v}_{E}\cdot\nabla{r}$, and $C^{(\ell)}_{zi}[F_{z1};h_{i}]$ is the linearization of $C_{zi}$ around $F_{Mz}$ and $F_{Mi}$,

[TABLE]

where the operator ${\cal K}$ is defined by

[TABLE]

The quantities $A$ and $\buildrel\leftrightarrow\over{\mathbf{M}}$, that depend on the main ions, are given by

[TABLE]

[TABLE]

A detailed derivation of expression (27) can be found in A. It will often be useful to employ that isotropy of the collision operator in velocity space implies

[TABLE]

From (29) and (30) it is clear that only $h_{i-}$ enters $A$ and only $h_{i+}$ enters $\buildrel\leftrightarrow\over{\mathbf{M}}$. Since the trace impurity approximation is made, the calculation of $h_{i}$ and $\varphi_{1}$ does not involve the impurities, and therefore $h_{i}$ and $\varphi_{1}$ enter as known functions in (24). Let us explicitly write equation (24) with the collision operator (27),

[TABLE]

In terms of the solution of the impurity drift-kinetic equation (2.2), the radial impurity flux across the surface $r$ reads

[TABLE]

In this paper we will calculate (33) for several impurity collisionality regimes. Since equation (2.2) is linear, we can treat the effect of its source terms independently. We will denote by $\Gamma_{z,a}$ the radial flux produced by the source terms containing the radial drift and $A$ (i.e., $\Gamma_{z,a}$ is the impurity flux obtained by neglecting main ion pressure anisotropy), and by $\Gamma_{z,b}$ the radial flux produced by the source terms containing $\buildrel\leftrightarrow\over{\mathbf{M}}$ (i.e., $\Gamma_{z,b}$ is the impurity flux due to main ion pressure anisotropy). By comparing the sizes of $\Gamma_{z,a}$ and $\Gamma_{z,b}$ we will learn when main ion pressure anisotropy is non-negligible.

So far, we have not made explicit assumptions about the aspect ratio of the stellarator. In the next subsection, we expand in large aspect ratio and derive the size and scaling with aspect ratio of $A$ and $\buildrel\leftrightarrow\over{\mathbf{M}}$.

2.3 Estimates for the size of $A$ and $\buildrel\leftrightarrow\over{\mathbf{M}}$ in large aspect ratio stellarators

When $\epsilon=a/R_{0}\ll 1$, one can write

[TABLE]

where $B_{0}$ is constant on the flux surface and $B_{1}=O(\epsilon B_{0})$. As for $\Upsilon$, we assume

[TABLE]

The factor $Z_{z}$ in the previous expression is easy to understand looking at (25) and noting that $e\varphi^{\prime}_{0}/T\sim 1/a$. For trace impurities, equilibrium conditions imply $\eta^{\prime}/\eta\sim Z_{z}/a$ as well.

We also need to order $\varphi_{1}$ with respect to $\epsilon$. We take

[TABLE]

which is the size of $\varphi_{1}$ that makes the radial components of the magnetic and $\mathbf{E}\times\mathbf{B}$ drifts comparable. In [5], it was analytically proven that when $Z_{z}e|\varphi_{1}|/T{\scriptstyle{{{}_{\displaystyle>}}\atop{\displaystyle\sim}}}\epsilon$, the effect of $\varphi_{1}$ on the radial neoclassical flux of collisional impurities matters.

It will be important to know the size of $A$ and $\buildrel\leftrightarrow\over{\mathbf{M}}$ when $\epsilon\ll 1$. An analytical expression for $A$ was given in [5] for low collisionality main ions, which implies

[TABLE]

Importantly, from the expression in [5], one can see that

[TABLE]

where $A_{0}$ is constant on the flux surface and $A_{1}(r,\theta,\zeta)=O(\epsilon A_{0})=O(\epsilon^{-1}\rho_{i*}v_{ti})$. An explicit expression for $A_{0}$ will be given in subsection 4.2.1.

We proceed to estimate the size of $\buildrel\leftrightarrow\over{\mathbf{M}}$ for low collisionality main ions. It is useful to employ $v$, $\lambda=v_{\perp}^{2}/(Bv^{2})$, $\sigma$ and $\phi$ as velocity coordinates for the main ions. In these coordinates, the parallel velocity is written as

[TABLE]

Since $h_{i}$ and $F_{Mi}$ are independent of the gyrophase, we have

[TABLE]

Employing

[TABLE]

to derive the identity

[TABLE]

valid for any vector $\mathbf{Y}$, we find

[TABLE]

and

[TABLE]

where we have used

[TABLE]

and

[TABLE]

The gyroaverage of (43) and (44) gives

[TABLE]

and

[TABLE]

In this paper, we always assume that the main ions have low collisionality. We can take $h_{i+}$, the component of $h_{i}$ that is even in $v_{||}$, to be zero for passing particles [15, 16, 19]. Trapped trajectories are defined by $\lambda\in[B_{\rm max}^{-1},B_{\rm min}^{-1}]$, where $B_{\rm max}$ and $B_{\rm min}$ are, respectively, the maximum and minimum values of $B$ on the flux surface. Due to (34), the region of phase space corresponding to trapped particles has a size $O(\epsilon B_{0}^{-1})$ in the coordinate $\lambda$. Hence, for trapped particles, $v_{||}\sim\epsilon^{1/2}v_{ti}$, $v\sim v_{\perp}\sim v_{ti}$, and we can approximate

[TABLE]

and

[TABLE]

Finally, for $\epsilon\ll 1$, the matrix $\buildrel\leftrightarrow\over{\mathbf{M}}$ has the form

[TABLE]

where we have used that in coordinates $\{v,\lambda,\phi\}$

[TABLE]

and an integration by parts in $\lambda$ has been performed in order to obtain the form of the first term on the right side of (2.3). Expression (2.3) implies

[TABLE]

If the main ions are in the $1/\nu$ regime, defined by $\nu_{ii*}\ll\epsilon^{3/2}$, then [16]

[TABLE]

where $\nu_{ii*}=R_{0}\nu_{ii}/v_{ti}$ is the main ion collisionality and $\nu_{ii}$ is the ion-ion collision frequency. Then, the size of the matrix $\buildrel\leftrightarrow\over{\mathbf{M}}$ is

[TABLE]

3 Impurities in the Pfirsch-Schlüter regime

The impurities are said to be in the Pfirsch-Schlüter regime when they are collisional; i.e.

[TABLE]

In order to solve (2.2), we write $F_{z1}=H_{a}+H_{b}$, where $H_{a}$ and $H_{b}$ are the solutions of

[TABLE]

and

[TABLE]

An explicit expression for the radial flux due to $H_{a}$ can be found in equation (25) of reference [5]. Assuming the ordering (36) for $\varphi_{1}$, the result of [5] gives

[TABLE]

In this section, we focus on the calculation of $\Gamma_{z,b}:=\Gamma_{z}[H_{b}]$. For that, we must solve (58) in the asymptotic limit $\nu_{zi*}\gg 1$. Subsections 3.1, 3.2, 3.3 and 3.4 are devoted to this. In subsection 3.5, an explicit expression for $\Gamma_{z,b}$ is provided. In subsection 3.6, the size of the ratio $\Gamma_{z,b}/\Gamma_{z,a}$ is estimated.

3.1 Organization of the calculation

We expand $H_{b}$ in powers of $1/\nu_{zi*}\ll 1$,

[TABLE]

with $H_{b}^{(k)}\propto(1/\nu_{zi*})^{k}$. The terms in (58) that scale with $\nu_{zi*}$ yield

[TABLE]

Since

[TABLE]

for any phase-space function $f$, in principle equation (61) might have a solvability condition. However, it does not, because the integral over velocities on the right-hand side of (61) vanishes.

To next order in $1/\nu_{zi*}$, equation (58) gives

[TABLE]

This equation does not have any solvability condition either because $H_{b}^{(0)}$ is even in $v_{||}$ and, therefore, the velocity integral of its right-hand side vanishes.

Finally, from terms in (58) that scale as $1/\nu_{zi*}$, we find

[TABLE]

We will not need to solve (64), but only deal with its solvability condition, which in this case is non-trivial. It is obtained by integrating (64) over velocities, giving

[TABLE]

where $Q(r)$ is a flux function to be computed. It is necessary to carry out the calculation to this order in the $1/\nu_{zi*}\ll 1$ expansion to completely determine $H_{b}^{(0)}$, which gives the dominant contribution of $H_{b}$ to $\Gamma_{z,b}$ in this collisionality regime; that is,

[TABLE]

In subsection 3.2, we solve (61). In subsection 3.3, the solution of (63) is given. In subsection 3.4, we work out (65).

3.2 Solution of (61)

3.2.1 Solution of the homogeneous equation associated to equation (61).

Let us call $H_{b,{\rm hom}}^{(0)}$ the solution of

[TABLE]

If we multiply by $H_{b,{\rm hom}}^{(0)}/F_{Mz}$, integrate over velocities and then integrate by parts, we get

[TABLE]

Every solution of (67) is a solution of (68). The latter is satisfied if and only if

[TABLE]

i.e. if and only if $H_{b,{\rm hom}}^{(0)}$ has the form

[TABLE]

It is easy to check that (70) satisfies (67). The function $N^{(0)}({r},\theta,\zeta)$ is determined by higher order equations.

3.2.2 Particular solution of the inhomogeneous equation (61).

The key for the explicit calculation of $\Gamma_{z,b}$ in the Pfirsch-Schlüter is the diagonalization of the operator ${\cal K}$. In B we find its eigenfunctions and eigenvalues, and explain how to write any function on velocity space as a linear combination of eigenfunctions of ${\cal K}$. Here, like in B, we use spherical coordinates $\{v,\beta,\phi\}$ in velocity space, where $v\in[0,\infty)$, $\beta\in[0,\pi]$ and $\phi\in[0,2\pi)$. Choose a right-handed set of orthonormal vectors $\{\hat{\mathbf{e}}_{1},\hat{\mathbf{e}}_{2},\hat{\mathbf{b}}\}$ and take a point $\mathbf{v}$ of velocity space. Spherical coordinates are defined by

[TABLE]

In addition, some expressions are simpler in terms of

[TABLE]

a variable that we will frequently use in what follows.

Noting that

[TABLE]

denoting by $P_{l}$ the Legendre polynomials and by $L_{p}^{(\alpha)}$ the generalized Laguerre polynomials (see definitions in B), and employing $P_{2}(\cos\beta)=(3\cos^{2}\beta-1)/2$, $L_{0}^{(\alpha)}(x)=1$ and $L_{1}^{(\alpha)}(x)=-x+\alpha+1$, we can rewrite (61) as

[TABLE]

where the right-hand side is now expressed as a linear combination of eigenfunctions of ${\cal K}$ and is, therefore, immediate to solve using (250). Then,

[TABLE]

where $H_{b,{\rm hom}}^{(0)}$ is given in (70) and the particular solution of the inhomogeneous equation, ${\tilde{H}}_{b}^{(0)}$, is

[TABLE]

3.3 Solution of (63)

Again, we need to write the right-hand side of (63) as a linear combination of eigenfunctions of the operator ${\cal K}$. Employing $P_{1}(\cos\beta)=\cos\beta$, $P_{3}(\cos\beta)=(5\cos^{3}\beta-3\cos\beta)/2$, $v_{||}=\sqrt{2T/m_{z}}\,x^{1/2}P_{1}$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

we recast (63) into

[TABLE]

The solution of the homogeneous equation is irrelevant for us and we can set it to zero. The solution of the inhomogeneous equation is

[TABLE]

3.4 Solvability condition (65)

Velocity space integrals in coordinates $x=m_{z}v^{2}/(2T)$ and $\xi:=\cos\beta=v_{||}/v$ read

[TABLE]

Applying this to the left side of (65), and employing the orthogonality relations of the Legrendre and Laguerre polynomials (see (243) and (247)), we get

[TABLE]

Multiplying this equation by $B^{2}/n_{z}$ and flux-surface averaging, we find

[TABLE]

Finally, we rewrite (3.4) as

[TABLE]

This is a magnetic differential equation that can be solved, numerically, determining $N^{(0)}/n_{z}$ up to an irrelevant additive flux function. For $\epsilon\ll 1$, the equation can be solved analytically. This is due to the fact that the right-hand side can be dropped because it is a factor of $\epsilon$ smaller than the second term on the left-hand side. Hence, we can approximate (3.4) by

[TABLE]

giving

[TABLE]

up to an additive flux function. Recalling (66), (70), (75) and (76), we compute $\Gamma_{z,b}$ in subsection 3.5.

3.5 Explicit expression for $\Gamma_{z,b}$

With the results of previous subsections, we can calculate the right-hand side of (66). To that end, it is convenient to use coordinates $x$ and $\beta$, and make repeated use of expansions in terms of Legendre and generalized Laguerre polynomials (see B). First, note that

[TABLE]

where we have used that, in a magnetohydrodynamic equilibrium $(\nabla\times\mathbf{B})\cdot\nabla r=0$, so that $\hat{\mathbf{b}}\times(\hat{\mathbf{b}}\cdot\nabla\hat{\mathbf{b}})=B^{-1}\hat{\mathbf{b}}\times\nabla B$. Then, we can write

[TABLE]

where

[TABLE]

and

[TABLE]

We start by computing $I_{1}$. It is immediate to find that

[TABLE]

For the second term on the right-hand side of (93), one gets

[TABLE]

which is easy to obtain using (85) and (76). Hence,

[TABLE]

As for $I_{2}$, using (70), (90) and (85), we get

[TABLE]

where the orthogonality relations (243) and (247) have been employed in the last step. Analogous, although a bit lengthier, manipulations allow us to calculate the second term on the right-hand side of (3.5),

[TABLE]

Then, from (98) and (99),

[TABLE]

Finally, for $\epsilon\ll 1$,

[TABLE]

where $\eta$ has been defined in (23) and $B_{0}$ has been defined in (34).

3.6 Comparison of the relative sizes of $\Gamma_{z,a}$ and $\Gamma_{z,b}$

Recalling (53) and (3.5), we find

[TABLE]

And, if the main ions are in the $1/\nu$ regime (see (54)),

[TABLE]

It is convenient to employ the relations

[TABLE]

and

[TABLE]

to write the estimate (103) in terms of $\rho_{z*}$ and $\nu_{zi*}$. Then,

[TABLE]

Let us rewrite the estimate (59) for the size of $\Gamma_{z,a}$ in a slightly different fashion,

[TABLE]

Taking the quotient of (106) and (107) it is easy to find the condition for the main ion pressure anisotropy to matter, assuming that the main ions are in the $1/\nu$ regime and the impurities are in the Pfirsch-Schlüter regime. Since

[TABLE]

the condition is

[TABLE]

where

[TABLE]

only depends on the magnetic geometry (via $\epsilon$), on the impurity charge and on the ratio of main ion to impurity mass. The parameter $\aleph$ will repeatedly appear when we evaluate the relative weight of $\Gamma_{z,a}$ and $\Gamma_{z,b}$ in different collisionality regimes.

The conditions $\nu_{ii*}\ll\epsilon^{3/2}$, $\nu_{zi*}\gg 1$ and $\aleph/\nu_{zi*}^{2}\gg 1$ can be simultaneously satisfied for highly-charged impurities, although the region of parameter space in which this happens is not large. This can be seen by expressing condition $\aleph/\nu_{zi*}^{2}\gg 1$ as

[TABLE]

and noting that $\nu_{zi*}\gg 1$ implies

[TABLE]

4 Impurities in the plateau regime

We deal with the regime

[TABLE]

in which passing particles are collisionless and trapped particles are collisional. This can be seen by comparing the parallel streaming and collision terms in (2.2). For passing particles, $v_{||}\sim v_{tz}$ and the effective collision frequency, $\nu_{zi,{\rm eff}}\sim v_{tz}^{2}\nu_{zi}/v_{||}^{2}\sim\nu_{zi}$, so the collision term is negligible compared to the parallel streaming term in the regime defined by (113). However, for trapped particles, $v_{||}\sim\epsilon^{1/2}v_{tz}$ and $\nu_{zi,{\rm eff}}\sim v_{tz}^{2}\nu_{zi}/v_{||}^{2}\sim\nu_{zi}/\epsilon$, implying that, if (113) is satisfied, the collision term is larger than the parallel streaming term. Although not in the context of impurity transport in stellarators but in the context of plasma transport in tokamaks, a good treatment of the plateau regime can be found in reference [20].

Trapped particles, being collisional, give negligible transport, and we focus on passing particles. It is useful to write the drift-kinetic equation in terms of

[TABLE]

so that (2.2) is recast into

[TABLE]

with

[TABLE]

For passing particles, we can define the transit average of a function $f(r,\theta,\zeta,{\cal E},\mu)$ as

[TABLE]

It is convenient to split (115) into two equations, one whose source is $\overline{S}$ and another one whose source is $S-\overline{S}$. Hence, we take $G_{z}=g_{z}+G_{z,p}$, where $g_{z}$ and $G_{z,p}$ are the solutions of

[TABLE]

and

[TABLE]

The transit average of $S$ is

[TABLE]

The subindex $p$ in $G_{z,p}$ indicates that this piece of the impurity distribution gives the plateau contribution, as we show in subsection 4.2.

4.1 Impurity flux due to $g_{z}$

Radial impurity transport due to $g_{z}$ is small. It is not difficult to show this by expanding $g_{z}=g_{z}^{(0)}+g_{z}^{(1)}+g_{z}^{(2)}+\dots$ in powers of $\nu_{zi*}\ll 1$. To lowest order in this expansion, one finds

[TABLE]

implying that $g_{z}^{(0)}$ is a flux function, which, as a consequence, does not give transport. To next order, we have the equation

[TABLE]

Its transit average,

[TABLE]

determines $g_{z}^{(0)}\sim{\nu_{zi}^{-1}}\overline{S}$. Subtracting (123) to (122), one finds

[TABLE]

From this equation, we deduce that $g_{z}^{(1)}\sim\nu_{zi_{*}}g_{z}^{(0)}$. Note that $g_{z}^{(1)}$ is odd in $v_{||}$ (because $g_{z}^{(0)}$ is even) and then $\Gamma_{z}[g_{z}^{(1)}]=0$. The lowest-order piece of the expansion contributing to $\Gamma_{z}$ is $g_{z}^{(2)}\sim\nu_{zi_{*}}^{2}g_{z}^{(0)}\sim\epsilon^{1/2}\nu_{zi_{*}}^{2}(h_{i+}/F_{Mi})F_{Mz}$. If the main ions are in the $1/\nu$ regime, this implies that

[TABLE]

We will see that this contribution is negligible compared to the one due to $G_{z,p}$.

4.2 Impurity flux due to $G_{z,p}$

4.2.1 Solution of equation (119).

The impurity flux due to $G_{z,p}$ is produced by a small collisional layer around $v_{||}=0$. In order to treat (119), we will employ coordinates $u$ and $\mu$, where now we are denoting the parallel velocity by $u$ (we have changed the notation for the parallel velocity from $v_{||}$ to $u$ to stress that now it is viewed as an independent variable). The total energy per unit mass reads, in terms of these coordinates,

[TABLE]

Denoting by $\hat{G}_{z,p}$ the function $G_{z,p}$ expressed in coordinates $u$ and $\mu$, (119) becomes

[TABLE]

where $\tilde{S}=S-\overline{S}$ is assumed to be expressed in terms of $u$ and $\mu$.

We expect $\hat{G}_{z,p}$ to have large derivatives with respect to $u$ in the small collisional layer of interest. Then, we can take

[TABLE]

By balancing

[TABLE]

we find the size of the layer in the coordinate $u$, $\Delta_{u}$,

[TABLE]

Then, using $\epsilon^{3/2}\ll\nu_{zi*}$ on the parallel streaming term on the left-hand side of (127), we find

[TABLE]

Next, we simplify $\tilde{S}=S-\overline{S}$ (recall (4)). First, employing $\Delta_{u}/v_{tz}\sim\nu_{zi*}^{1/3}\ll 1$, we get

[TABLE]

where $A_{0}$ has been introduced in (38). An explicit expression for $A_{0}$ can be found below, in equation (153).

Second, expanding in $\epsilon$ and noting that in the layer $v_{||}^{2}\ll v_{\perp}^{2}/2$, we have

[TABLE]

where

[TABLE]

Recall that, to lowest order in $\epsilon$,

[TABLE]

and the impurity density is a flux function.

We can rewrite (127) as

[TABLE]

with

[TABLE]

The two first terms on the right-hand side of (4.2.1) have a size $O(a^{-1}\rho_{i*}v_{ti}\hat{F}_{Mz})$ and the two last ones have a size $O(\nu_{zi}\epsilon^{1/2}(h_{i+}/F_{Mi})\hat{F}_{Mz})$.

In order to give a completely explicit calculation, we will be more specific about our spatial coordinates. We use Boozer coordinates [21], so that the magnetic field can be simultaneously written as

[TABLE]

and

[TABLE]

where $\iota({r})$ is the rotational transform, $\Psi_{t}({r})$ is the toroidal magnetic flux over $2\pi$, and $I_{t}$ and $I_{p}$ are the toroidal and poloidal currents, respectively. In Boozer coordinates, the volume element is conveniently written in terms of the magnitude of $\mathbf{B}$,

[TABLE]

Note that, to lowest order in $\epsilon$,

[TABLE]

The parallel streaming operator in Boozer coordinates is

[TABLE]

and to lowest order in $\epsilon$,

[TABLE]

In Boozer coordinates, the radial drift reads

[TABLE]

Taking an $\epsilon\ll 1$ expansion and using that in the layer $u^{2}\ll\mu B$, we have

[TABLE]

Equation (4.2.1) becomes

[TABLE]

with

[TABLE]

In (4.2.1), $\hat{F}_{Mz}$, $\hat{\Upsilon}$ and $\mathbf{v}_{d}\cdot\nabla r$ must be understood as given by (134), (137) and (4.2.1).

At this point, it is useful to split (146) into two equations, in analogy to the treatment of Section 3. We take

[TABLE]

where $H_{a}$ and $H_{b}$ are the solutions of

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

It is convenient to write explicitly $A_{0}$, that can be read off from equation (24) of reference [5]. Namely,

[TABLE]

Here, $u$ is the solution of $\hat{\mathbf{b}}\cdot\nabla u=2B^{-3}(\hat{\mathbf{b}}\times\nabla r)\cdot\nabla B$ with vanishing boundary condition for $u$ at the point of the flux surface where $B$ takes the value $B_{\rm max}$, and the flux functions $f_{c}$ and $f_{s}$ are given in C.

Equations (149) and (150) are solved in Fourier space. We write

[TABLE]

with $j=a,b$. Using the Fourier expansions

[TABLE]

and

[TABLE]

and defining the non-dimensional parameter

[TABLE]

equation (154) translates into

[TABLE]

Making the change of variable (we do not change the name of $H_{j,mn}$)

[TABLE]

we find

[TABLE]

The solution to this equation is

[TABLE]

and, employing the variable $u$,

[TABLE]

4.2.2 Explicit expression for the impurity flux in the plateau regime.

We can obtain an explicit expression for $\Gamma_{z}[G_{z,p}]$. Using coordinates $u$ and $\mu$ and to lowest order in $\epsilon$,

[TABLE]

Now, we employ the Fourier expansions (155) and

[TABLE]

Then,

[TABLE]

Here, the approximation (4.2.1) to the radial drift is understood, so that

[TABLE]

with

[TABLE]

[TABLE]

At this point, we note the identity

[TABLE]

which implies that the component of $H_{j,mn}$ that is even in $u$ (the odd component does not give radial transport), $H_{j,mn}^{\rm even}$, is proportional to a delta function for small collisionality,

[TABLE]

Using this and (157) in (165), we get

[TABLE]

where we have assumed $\Psi^{\prime}_{t}>0$. Recall the decomposition (148) and that $\hat{G}_{z,p}$ is $G_{z,p}$ expressed in different coordinates. Clearly, equation (171) is the sum of a contribution from $H_{a}$, that we denote by $\Gamma_{z,a}\equiv\Gamma_{z}[H_{a}]$,

[TABLE]

and a contribution from $H_{b}$, that we denote by $\Gamma_{z,b}\equiv\Gamma_{z}[H_{b}]$,

[TABLE]

We recall that $\tilde{S}_{a}$ and $\tilde{S}_{b}$ are defined in (151) and (4.2.1), respectively.

The integrals over $\mu$ in (172) and (173) can be taken analytically. Doing this and rearranging terms, one can write $\Gamma_{z,a}$ as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and an asterisk denotes complex conjugation.

As for $\Gamma_{z,b}$, we find

[TABLE]

with

[TABLE]

and

[TABLE]

4.2.3 Size of the impurity flux in the plateau regime.

From the expressions for $\Gamma_{z,a}$ and $\Gamma_{z,b}$ found in subsection 4.2.2, we learn that

[TABLE]

and

[TABLE]

In (182), the factor $Z_{z}/\epsilon$ comes from the estimate (35) for $\Upsilon$. Note also that in (151) the term containing $\Upsilon$ has the same typical size as the term proportional to $A_{0}$; hence, in principle, temperature screening is possible.

If the main ions are in the $1/\nu$ regime, then

[TABLE]

Taking the ratio of (184) to (182) we deduce that,

[TABLE]

Hence, when the impurities are in the plateau regime and the main ions are in the $1/\nu$ regime, transport driven by main ion pressure anisotropy is significant if

[TABLE]

where $\aleph$ has been defined in (110).

Finally, note that (recall (125))

[TABLE]

so that transport due to passing particles with large parallel velocity is negligible, as announced above.

5 Impurities in the $1/\nu$ regime

Let us assume that

[TABLE]

In this case, the parallel streaming term is much larger than the collision term in (2.2) for both passing and trapped particles. Therefore, one can average over the motion along magnetic field lines, and it is convenient to choose spatial coordinates in which such averaged is performed in a transparent way. We use $(r,\alpha,l)$, where $\alpha=\theta-\iota\zeta\in[0,2\pi)$ is an angular coordinate labeling magnetic field lines and $l\in[0,l_{\rm max}(r,\alpha))$ is the arc length along the field line. Then, the magnetic field can be written as (cf. equation (138))

[TABLE]

Expanding $F_{z1}=F_{z1}^{(-1)}+F_{z1}^{(0)}+\dots$, with $F_{z1}^{(k)}\sim\nu_{zi*}^{k}\rho_{z*}F_{Mz}$, equation (2.2) gives, to lowest order,

[TABLE]

The distribution $F_{z1}^{(-1)}$ can be taken to be zero for passing particles in large aspect ratio stellarators and/or stellarators with a sufficiently high degree of optimization (see references [15, 16, 19]). For trapped particles, $F_{z1}^{(-1)}(r,\alpha,{\cal E},\mu)$ is found by orbit-averaging the drift-kinetic equation to next order in the $\nu_{zi*}$ expansion,

[TABLE]

where we have used that the term proportional to $A$ in the collision operator (27) is odd in $v_{||}$ and orbit-averages to zero.

When $\varphi_{1}$ is neglected, it is known that terms containing derivatives with respect to the pitch angle coordinate, $\lambda=2\mu/v^{2}$, dominate the differential piece of the collision operator (for us, ${\cal K}$) when $\epsilon\ll 1$ [14, 16], and the collision term simplifies significantly. However, when $\varphi_{1}$ is included, $\lambda$ is not a good coordinate to determine whether a particle is trapped or passing, or to identify the largest piece of the collision operator. In our case, the operator ${\cal K}$ can still be written as a differential operator in a single variable, but this variable is not the pitch angle. We explain this in subsection 5.1.

5.1 Collision operator for large aspect ratio

Let us take velocity coordinates $(\hat{\cal E},\mu,\phi)$, where $\hat{\cal E}$ is defined as follows. We introduce the function $U$,

[TABLE]

and define $U_{\rm M}(r,\mu)$ as the maximum value of $U$ on the magnetic surface $r$ for a given value of $\mu$. Then, $\hat{\cal E}$ is defined by

[TABLE]

The coordinate $\mu$ takes values in $[0,\infty)$. Given a value of $\mu$, the coordinate $\hat{\cal E}$ takes values in the interval $[U_{\rm m}(r,\mu)-U_{\rm M}(r,\mu),\infty)$, where $U_{\rm m}(r,\mu)$ is the minimum value of $U$ on the magnetic surface $r$ for that particular value of $\mu$.

In terms of $\hat{\cal E}$, distinguishing between passing and trapped particles is simple. If $\hat{\cal E}>0$, the particle is passing. The particle is trapped if

[TABLE]

Denote by $\hat{F}_{z1}^{(-1)}$ the distribution $F_{z1}^{(-1)}$ expressed as a function of $(r,\alpha,\hat{\cal E},\mu)$. In coordinates $\hat{\cal E}$, $\mu$ and $\phi$, the operator ${\cal K}$ reads

[TABLE]

Here,

[TABLE]

[TABLE]

and

[TABLE]

with $\mathbf{v}_{\perp}=\mathbf{v}-v_{||}\hat{\mathbf{b}}$ and $B_{M}=B(r,\alpha_{M},l_{M})$, where $\alpha_{M}$ and $l_{M}$ correspond to the point on the magnetic surface where $U=U_{M}$. The Jacobian, $J$, is given by

[TABLE]

Using that $\partial_{\phi}(\hat{F}_{z1}^{(-1)}/F_{Mz})=0$, $\nabla_{v}\hat{\cal E}\cdot\nabla_{v}\phi=0$ and $\nabla_{v}\mu\cdot\nabla_{v}\phi=0$, we find

[TABLE]

In this expression and the expressions that follow in this section, $v_{||}$, $F_{Mz}$, etc. are viewed as functions of the independent coordinates $\hat{\cal E}$ and $\mu$. In these coordinates,

[TABLE]

Recalling (34), (36) and (194), we get the estimate

[TABLE]

for trapped particles. Hence, $\hat{F}_{z1}^{(-1)}$ will typically have large derivatives with respect to $\hat{\cal E}$,

[TABLE]

whereas it will have small derivatives with respect to $\hat{\mu}$,

[TABLE]

Again, using (34) and (36), we find

[TABLE]

Hence, for trapped particles,

[TABLE]

The term containing two derivatives with respect to $\hat{\cal E}$ in (5.1) is larger than the remaining ones by a factor $1/\epsilon$. Noting that, at large aspect ratio, $|\partial_{\hat{\cal E}}F_{Mz}/F_{Mz}|\ll|\partial_{\hat{\cal E}}\hat{F}_{z1}^{(-1)}/\hat{F}_{z1}^{(-1)}|$ and $|\nabla_{v}\hat{\cal E}|^{2}\simeq v_{||}^{2}$, we finally arrive at

[TABLE]

5.2 Impurity flux for impurities in the $1/\nu$ regime

Using (207), equation (191) becomes

[TABLE]

In analogy to what we did in previous sections, we write $F_{z1}^{(-1)}={H_{a}}+{H_{b}}$, where ${H_{a}}$ and ${H_{b}}$ are the solutions of

[TABLE]

and

[TABLE]

These equations can be integrated numerically in a straightforward way.

Noting that

[TABLE]

from (209) we deduce

[TABLE]

where the factor $Z_{z}$ comes from the estimate (35) for $\Upsilon$.

Recalling (53), from (5.2) we get

[TABLE]

Employing (212) and (213) in (33), we reach the estimates

[TABLE]

and

[TABLE]

If the main ions are in the $1/\nu$ regime (recall (54), (104) and (105)), one gets

[TABLE]

The ratio of (216) to (214) gives

[TABLE]

Therefore, when both the impurities and the main ions are in the $1/\nu$ regime, main ion pressure anisotropy must be taken into account if

[TABLE]

Note that $\aleph{\scriptstyle{{{}_{\displaystyle>}}\atop{\displaystyle\sim}}}1$ requires to have highly charged impurities. Since highly charged impurities tend to be collisional, we expect our calculation for the $1/\nu$ regime to be less relevant than the calculations for the plateau and Pfirsch-Schlüter regimes.

6 Numerical evaluation of the analytical results

Figure 1 synthesizes the scalings and typical sizes for $\Gamma_{z}$ derived in previous sections for several collisionality regimes. We have found that there are two essentially different situations as far as the relevance of main ion impurity anisotropy is concerned, distinguished by the value of $\aleph=\epsilon^{3/2}Z_{z}^{2}m_{i}/m_{z}$. If $\aleph\ll 1$, main ion pressure anisotropy is negligible. If $\aleph\gg 1$, main ion pressure anisotropy becomes the main drive for impurity transport. As remarked at the end of Section 5, the condition $\aleph\gtrsim 1$ typically requires impurities with high electric charge, and such impurities are more likely to be in the plateau or Pfirsch-Schlüter regimes than in the $1/\nu$ regime.

Next, we give realistic examples to illustrate the relevance of main-ion-pressure-anisotropy-driven neoclassical impurity transport by numerically evaluating the analytical expressions derived for $\Gamma_{z,a}$ and $\Gamma_{z,b}$ in the plateau and Pfirsch-Schlüter regimes. The expression for $\Gamma_{z,a}$ in the Pfirsch-Schlüter regime can be found in equation (25) of [5]. The expressions for $\Gamma_{z,b}$ in the Pfirsch-Schlüter regime and for $\Gamma_{z,a}$ and $\Gamma_{z,b}$ in the plateau regime have been derived in the present paper (see (3.5) and (171)).

As trace impurity, we take the charge state $Z_{z}=40$ of tungsten and we set $\eta^{\prime}=0$. The main ion species is deuterium. We will perform a scan in the density of the main ions, while keeping constant the rest of the simulation parameters, in particular $n_{i}^{\prime}/n_{i}$. For each point in the scan, we solve the drift-kinetic equation of the main ions using the code KNOSOS [22]. In these examples, we set $\varphi_{1}=0$. We emphasize that KNOSOS solves bounce-averaged equations and can calculate $h_{i+}$, but not $h_{i-}$. This code is employed here to evaluate $\Gamma_{z,b}$, which is the component of the impurity flux that depends on $h_{i+}$.

First, we use an inward-shifted configuration of the Large Helical Device (LHD), characterized by having major radius $R_{0}=3.67\,$m, and we focus on the flux-surface $r/a=0.8$, where $a=0.64\,$m. The magnetic field and profiles are those of the plasma simulated in [5]; in particular, $T=1.3\,$keV, $aT^{\prime}/T=-3.4$ and $an_{i}^{\prime}/n_{i}=-2.6$. The main ion density is scanned so that the collisionality goes from $\nu_{zi*}\sim 10^{-2}$ to $\nu_{zi*}\sim 10$ and we take $\varphi^{\prime}_{0}=0$. When solving the main ion drift-kinetic equation with KNOSOS, we have switched off the component of the magnetic drift that is tangent to the flux surface, whereas the tangential component of the $\mathbf{E}\times\mathbf{B}$ drift vanishes because $\varphi^{\prime}_{0}=0$. An important consequence of not having tangential drifts is that the main ions are in the $1/\nu$ regime for all values of the collisionality in the range selected here (we come back to this below). In figure 2(top), the total flux $\Gamma_{z,a}+\Gamma_{z,b}$ is compared to the flux obtained by neglecting the effect of the main ion pressure anisotropy, showing that neglecting $\Gamma_{z,b}$ would lead to an incorrect result for the impurity flux in most of the collisionality range, and that, actually, $\Gamma_{z}$ is dominated by $\Gamma_{z,b}$ below $\nu_{zi*}\sim 1$. Figure 2(bottom) exhibits the collisionality dependence of $\Gamma_{z,a}$ and $\Gamma_{z,b}$ in the plateau and Pfirsch-Schlüter regimes discussed in this paper. We note that, in this example, $\aleph=1.27$.

In figure 3 we give the results of an analogous calculation on the flux surface $r/a=0.8$ of the standard configuration of W7-X, with $R_{0}=5.51\,$m and $a=0.51\,$m. The trace impurity and main plasma parameters and profiles (normalized to the minor radius) are the same as in the LHD case. Although the results are qualitatively similar, we see that the weight of $\Gamma_{z,b}$ is smaller in W7-X due to its larger aspect ratio compared to LHD. In this example, $\aleph=0.49$.

In figure 4 we show the result of redoing the W7-X calculations after switching on the component of the magnetic drift tangent to the flux surface in the main ion drift-kinetic equation (the effect of including the tangential drift in the impurity drift-kinetic equation is not studied in this paper). In these simulations, the main ions are not in the $1/\nu$ regime for the whole range of collisionality values. This is evident from figure 4(bottom), where, for sufficiently low values of the collisionality, $\Gamma_{z,b}$ in the ‘plateau’ regime is not independent of $\nu_{zi*}$ and, in the Pfirsch-Schlüter regime, $\Gamma_{z,b}$ does not scale with the inverse of $\nu_{zi*}$. The reason is that, for collisionality low enough for the tangential drifts to be relevant, the distribution function of the main ions ceases to scale with the inverse of the collisionality [15, 16] (the effect would be similar for an $\mathbf{E}\times\mathbf{B}$ drift caused by a relatively small radial electric field). Therefore, the relative weight of $\Gamma_{z,b}$ with respect to $\Gamma_{z,a}$ decreases compared to the case without tangential drifts (see figure 4(top) and figure 3(top)). From figure 4(top) it is clear that, for $\nu_{zi*}\lesssim 1$, neglecting $\Gamma_{z,b}$ would still be incorrect.

7 Conclusions

Main ion dynamics affects impurity transport via two basic mechanisms: determining the electric field and through impurity-ion collisions. In this paper, we focus on the second one.

To lowest order in a mass ratio expansion $\sqrt{m_{i}/m_{z}}\ll 1$, where $m_{i}$ and $m_{z}$ are the main ion and the impurity masses, the impurity-ion collision operator only depends on the component of the main ion distribution that is odd in the parallel velocity, $h_{i-}$. In fluid terms, using this lowest-order collision operator implies that the radial impurity flux, $\Gamma_{z}$, is driven only by the parallel friction of the impurities with the main ions. However, this description might not be complete for stellarators. Whereas terms in the impurity-ion collision operator proportional to the component of the main ion distribution that is even in the parallel velocity, $h_{i+}$, are small in $\sqrt{m_{i}/m_{z}}$, they can end up driving a non-negligible amount of impurity transport because $h_{i+}$ can be large in regimes of main ion low collisionality. Hence, in principle, one can distinguish two collisional drives for stellarator impurity transport with different physical origin: parallel friction and main ion pressure anisotropy. In this article, we have derived the impurity-ion collision operator keeping terms of sufficiently high order in $\sqrt{m_{i}/m_{z}}\ll 1$ to retain the effect of main ion pressure anisotropy. With this collision operator, and using the trace impurity approximation, we have solved the impurity drift-kinetic equation in the Pfirsch-Schlüter, plateau and $1/\nu$ regimes, assuming that the main ions have low collisionality and that the $\mathbf{E}\times\mathbf{B}$ drift term in the impurity drift-kinetic equation is small compared to the parallel streaming term. For each regime, we have deduced the conditions (plasma, impurity and geometric parameters) that determine whether or not main ion pressure anisotropy gives a significant contribution to $\Gamma_{z}$. Finally, we have illustrated the potential importance of this effect by numerically evaluating the analytical expressions of $\Gamma_{z}$ for realistic stellarator plasmas.

The expressions for $\Gamma_{z}$ derived in this article, together with that in [5], and their coupling to the code KNOSOS (in the way described in Section 6) provide a numerical tool for a fast evaluation of neoclassical impurity transport in stellarators, with potential applications to the design of operation scenarios and stellarator optimization.

I. C. and J. L. V. thank H. Sugama and S. Satake for interesting and helpful discussions. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This research was supported in part by grant ENE2015-70142-P, Ministerio de Economía y Competitividad, Spain and grant PGC2018-095307-B-I00, Ministerio de Ciencia, Innovación y Universidades, Spain.

Appendix A Mass ratio expansion of the linealized impurity-ion collision operator

Let us derive expression (27). The Landau operator for collisions between species $s$ and $s^{\prime}$ (see, for example, [20]) is

[TABLE]

where $\mathbf{g}=\mathbf{v}-\mathbf{v}^{\prime}$,

[TABLE]

and $\ln\Lambda$ is the Coulomb logarithm.

Take two Maxwellians with the same temperature,

[TABLE]

and write $F_{s}=F_{Ms}+h_{s}$ and $F_{s^{\prime}}=F_{Ms^{\prime}}+h_{s^{\prime}}$. Recalling that $C_{ss^{\prime}}[F_{Ms},F_{Ms^{\prime}}]=0$ and dropping terms quadratic in the functions $h_{s}$ and $h_{s^{\prime}}$, we obtain the linearized collision operator

[TABLE]

From here on, we select $s=z$ and $s^{\prime}=i$. Hence,

[TABLE]

Let us assume $\sqrt{m_{i}/m_{z}}\ll 1$. First, we focus on $C_{zi}[h_{z},F_{Mi}]$. Expanding in mass ratio,

[TABLE]

so that

[TABLE]

In the second term on the right-hand side of this expression, the leading term (i.e. the first term) on the right-hand side of (224) has given a vanishing contribution due to $\nabla_{v^{\prime}}F_{Ms}(\mathbf{v}^{\prime})=-(m_{s}\mathbf{v}^{\prime}/T)F_{Ms}(\mathbf{v}^{\prime})$ and $\nabla_{v^{\prime}}\nabla_{v^{\prime}}v^{\prime}\cdot\mathbf{v}^{\prime}=0$. It is easy to check that both terms on the right-hand side of (A) have the same size in the mass ratio expansion.

The integrals in (A) can be computed analytically. Noting that

[TABLE]

and

[TABLE]

a straightforward calculation gives

[TABLE]

and

[TABLE]

Hence, to lowest order in the $\sqrt{m_{i}/m_{z}}$ expansion,

[TABLE]

We turn to $C_{zi}[F_{Mz},h_{i}]$. Employing (224),

[TABLE]

In the first term on the right-hand side of (A), it has been enough to keep the contribution of the first term on the right-hand side of (224) because the second term on the right-hand side of (224) always gives a negligible contribution. However, in the second term on the right-hand side of (A), we have kept the contributions coming from the two terms on the right-hand side of (224).

From (A), we easily obtain

[TABLE]

The third term on the right-hand side of (A) is proportional to an odd moment of $h_{i}$. The rest of the terms are proportional to even moments of $h_{i}$ and they have the same typical size. Grouping terms and recalling (230), we finally obtain expression (27).

Appendix B Eigenfunctions of the operator $\cal K$

In this appendix we calculate the eigenfunctions of the operator

[TABLE]

We begin by writing ${\cal K}$ using spherical coordinates $\{v,\beta,\phi\}$ in velocity space, defined in (71) and the paragraph before it. We are interested in the action of ${\cal K}$ on gyrophase-independent functions, $\partial_{\phi}g\equiv 0$. Using this, the formulae

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and the expression

[TABLE]

for the divergence of a vector field, we obtain

[TABLE]

It is useful to rewrite ${\cal K}$ in coordinates $x:=m_{z}v^{2}/(2T)$ and $\beta$. That is,

[TABLE]

We do not change the notation for $h$ or $F_{Mz}$, although in the previous expression we assume that they are written using $x$ instead of $v$. For example,

[TABLE]

Take $h=G_{l}(x)P_{l}(\cos\beta)F_{Mz}$, where $\{P_{l}\,|\,l=0,1,2\dots\}$ are the Legendre polynomials. These polynomials satisfy the differential equations

[TABLE]

and the orthogonality relations

[TABLE]

Using (242), one gets

[TABLE]

with

[TABLE]

Let $L_{p}^{(\alpha)}(x)$, $p=0,1,2,\dots$, $\alpha\in\mathbb{R}$, denote the generalized Laguerre polynomials. They satisfy the differential equations

[TABLE]

and the orthogonality relations

[TABLE]

Try $G_{l}(x)=U_{p,l}(x)$, where $U_{p,l}(x):=x^{l/2}L_{p}^{(l+1/2)}(x)$. Then,

[TABLE]

Employing (246), we deduce that

[TABLE]

Therefore,

[TABLE]

Appendix C Definitions of $f_{c}$ and $f_{s}$

We give here the expressions for $f_{c}$ and $f_{s}$, that enter equation (153). These expressions are taken from [19].

The flux function $f_{c}$ is defined by

[TABLE]

As for $f_{s}$,

[TABLE]

where

[TABLE]

is defined for $\lambda<B_{\rm max}^{-1}$. The integral is taken along the magnetic field line, $l$ is the arc length along the line and $l_{\rm max}$ gives the position at which $B$ takes the value $B_{\rm max}$.