2.1 Flux-friction relation

The formulation of the radial impurity flux in a stellarator in terms of a flux-friction relation was detailed in Sugama & Nishimura (Reference Sugama and Nishimura2002), Braun & Helander (Reference Braun and Helander2010). The flux is decomposed into a sum of contributions, the first due to friction against the background bulk ions and the second the result of the impurity pressure anisotropy,

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{z}=\left\langle \int f_{z}(\boldsymbol{v}_{dz}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r)\,\text{d}^{3}v\right\rangle =\frac{1}{Ze}\left\langle uBR_{z\Vert }+(p_{z\Vert }-p_{z\bot })\frac{\unicode[STIX]{x1D735}_{\Vert }(uB^{2})}{2B}\right\rangle . & & \displaystyle\end{eqnarray}$$

The equilibrium function $u$ satisfies $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}u=-\boldsymbol{b}\times \unicode[STIX]{x1D735}r\boldsymbol{\cdot }\unicode[STIX]{x1D735}(B^{-2})$ . The effect of friction against electrons is small in the electron–ion mass ratio, so it is neglected throughout.

The second term in (2.4) is necessary for the determination of the radial electric field, as it is not intrinsically ambipolar. However, the pressure anisotropy may be expected to be weak for a collisional species, as we discuss below, and the particle flux dominated by the friction drive. This will give rise to a flux very different to that in a tokamak, due to the different structure of the bulk ion distribution function in a low collisionality stellarator. With the linearised, gyroaveraged, collision operator for species $a$ denoted by $C_{a}=\sum _{b}C_{ab}$ , where the sum is over the ion species present, we can compare the magnitude of the flux driven by the parallel friction $R_{z\Vert }=m_{z}\int v_{\Vert }C_{z}(f_{z})\,\text{d}^{3}v$ to that expected due to the species’ pressure anisotropy $p_{z\Vert }-p_{z\bot }$ , by considering the first terms in an expansion of the drift-kinetic equation governing the impurity behaviour.

The expansion is taken as usual with respect to the magnetisation parameter $\unicode[STIX]{x1D70C}_{\ast z}=\unicode[STIX]{x1D70C}_{z}/L$  (Helander Reference Helander2014), where $\unicode[STIX]{x1D70C}_{a}$ is the gyroradius of species $a$ , and $L$ is a characteristic length scale perpendicular to the background magnetic field. We assume $Z\gg 1$ , but not so large to require that $\unicode[STIX]{x1D70C}_{\ast z}$ is higher order with respect to $\unicode[STIX]{x1D70C}_{\ast i}$ . Taking a characteristic parallel length scale $L_{\Vert }$ , which will satisfy $L_{\Vert }>L$ , the ratio of the contributions to the flux in (2.4) is approximately $(p_{z\Vert }-p_{z\bot })/R_{z\Vert }L_{\Vert }$ . Due to the high collisionality the leading-order piece of the expanded distribution function $f_{z}=f_{z0}+f_{z1}+\cdots \,$ will be a Maxwellian, $f_{Mz}=(n_{z}/\unicode[STIX]{x03C0}^{3/2}v_{Tz}^{3})\exp (-v^{2}/v_{Tz}^{2})$ , where the thermal velocity of a species is $v_{Ta}=\sqrt{2T_{a}/m_{a}}$ . The first-order drift-kinetic equation for the distribution function $f_{z1}$ then takes the form

(2.5) $$\begin{eqnarray}\displaystyle C_{z}(f_{z1})=v_{\Vert }\unicode[STIX]{x1D735}_{\Vert }\,f_{z1}+\boldsymbol{v}_{dz}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{Mz}, & & \displaystyle\end{eqnarray}$$

where the independent velocity space coordinates are taken to be $\unicode[STIX]{x1D716}_{z}$ and $\unicode[STIX]{x1D707}_{z}$ . In a subsidiary expansion of (2.5) with respect to collisionality, the pressure anisotropy will appear in first order, as usual for a collisional species (Braun & Helander Reference Braun and Helander2010). We define the collisionality here as $\unicode[STIX]{x1D708}_{\ast ab}=\unicode[STIX]{x1D708}_{ab}/\unicode[STIX]{x1D714}_{ta}=L_{\Vert }/\unicode[STIX]{x1D706}_{mfp}^{ab}$ , where $\unicode[STIX]{x1D714}_{ta}$ is the characteristic transit frequency of species $a$ along the magnetic field, $\unicode[STIX]{x1D708}_{ab}$ represents the characteristic collision frequency between species $a$ and $b$ , and the mean free path $\unicode[STIX]{x1D706}_{mfp}^{ab}=v_{Ta}/\unicode[STIX]{x1D708}_{ab}$ . Comparing the collision and drift terms in (2.5), remembering that the flows of all species are at the diamagnetic level ${\sim}\unicode[STIX]{x1D70C}_{\ast a}v_{Ta}$ , and that a factor $Ze\unicode[STIX]{x1D6F7}_{0}/T$ is introduced through the gradient of $F_{Mz}$ , we may expect $p_{z\Vert }-p_{z\bot }\sim Zp_{z}v_{dz}/\unicode[STIX]{x1D708}_{zz}L\sim Zp_{z}\unicode[STIX]{x1D70C}_{\ast z}/\unicode[STIX]{x1D708}_{\ast zz}$ . The parallel friction between unlike species drives the flux, and for the case of disparate mass ions considered here we may approximate it as $R_{zi\Vert }\sim m_{i}n_{i}(V_{i\Vert }-V_{z\Vert })\unicode[STIX]{x1D708}_{iz}\sim m_{i}n_{i}\unicode[STIX]{x1D70C}_{\ast i}v_{Ti}\unicode[STIX]{x1D708}_{iz}$ , where $V_{i\Vert }$ and $V_{z\Vert }$ are the bulk ion and impurity parallel flows respectively. (The form of the collision operator is discussed in more detail in §§ 2.3–2.4.)

We therefore find that the pressure anisotropy drive will be small when the collisionalities satisfy

(2.6) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D708}_{\ast iz}}\ll \frac{n_{i}}{n_{z}}\sqrt{\frac{m_{i}}{m_{z}}}\unicode[STIX]{x1D708}_{\ast zz}. & & \displaystyle\end{eqnarray}$$

When both species are collisional, as in Braun & Helander (Reference Braun and Helander2010), this condition is clearly satisfied, even for non-trace impurity levels, and the pressure anisotropy drive is always small. In the mixed collisionality case here, with $Z>1$ , this condition limits how collisionless the bulk ions can be compared to the impurities – otherwise there would be a negligible frictional driving force. It is most readily satisfied for highly charged impurities as the bulk ions become more collisionless. With density and temperature values around $n_{i}\sim 10^{20}~\text{m}^{-3}$ and $T_{i}$ of a few keV in high performance scenarios (Geiger et al. Reference Geiger, Beidler, Feng, Maaßberg, Marushchenko and Turkin2015; Dinklage et al. Reference Dinklage, Sakamoto, Yokoyama, Ida, Baldzuhn, Beidler, Cats, McCarthy, Geiger and Kobayashi2017) giving bulk (hydrogen) ion collisionalities in the range $\unicode[STIX]{x1D708}_{\ast ii}\sim 0.01$ and lower, the relation $\unicode[STIX]{x1D708}_{\ast zz}=\unicode[STIX]{x1D708}_{\ast ii}n_{z}Z^{4}/n_{i}$ for equal species’ temperature indicates that, for example, trace levels of common impurities such as Fe $^{24+}$ and Ar $^{16+}$ satisfy the mixed collisionality regime considered. Further specific examples of potential parameter values can be found in Helander et al. (Reference Helander, Newton, Mollén and Smith2017a ). We assume here that the ordering in (2.6) is satisfied and we will take the dominant drive of the transport in (2.4) to come from the parallel friction. Momentum conservation in collisions then allows us to write the impurity flux in terms of the bulk ion–impurity parallel friction,

(2.7) $$\begin{eqnarray}\displaystyle R_{zi\Vert }=-R_{iz\Vert }=-m_{i}\int v_{\Vert }C_{iz}(f_{i},f_{z})\,\text{d}^{3}v. & & \displaystyle\end{eqnarray}$$

In the next subsections, we develop the expressions for the bulk ion and impurity distribution functions required to evaluate this friction, using model collision operators to treat the low collisionality regimes analytically.

2.2 Bulk ion distribution function

The bulk ion distribution function can be treated throughout the low collisionality regimes of interest here using a recently developed formulation, which was detailed in Helander et al. (Reference Helander, Newton, Mollén and Smith2017a ). The distribution is split into pieces which are even and odd, $f_{i}^{\pm }$ , with respect to the parallel velocity $v_{\Vert }=\unicode[STIX]{x1D70E}|v_{\Vert }|$ , where $\unicode[STIX]{x1D70E}=\pm 1$ . The full bulk ion drift-kinetic equation then splits into two equations,

(2.8) $$\begin{eqnarray}\displaystyle v_{\Vert }\unicode[STIX]{x1D735}_{\Vert }f_{i}^{\mp }=C_{i}^{\pm }(f_{i})-\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{\pm }, & & \displaystyle\end{eqnarray}$$

where $C_{i}^{\pm }(f_{i})$ denotes the even and odd parts of the collision operator $C_{i}(f_{i})$ and the independent coordinates are taken to be $(r,\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716}_{i},\unicode[STIX]{x1D707}_{i},\unicode[STIX]{x1D70E})$ , where $\unicode[STIX]{x1D6FC}$ labels different field lines on the same flux surface and $l$ gives the arc length along the magnetic field.

The orbit average may be introduced, which annihilates the left-hand side of (2.8) and is essentially a time average over the particle trajectory neglecting the drift motion. The parameter $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}_{i}/\unicode[STIX]{x1D716}_{i}$ divides phase space into regions describing particles trapped in the magnetic field structures, for which $\unicode[STIX]{x1D706}>1/B_{\text{max}}$ where $B_{\text{max}}(r)$ is the maximum value of the magnetic field strength on the flux surface, and those able to circulate freely. For circulating particles, the orbit average of an arbitrary function $g$ is defined as

(2.9) $$\begin{eqnarray}\displaystyle \overline{g}(r,\unicode[STIX]{x1D716}_{i},\unicode[STIX]{x1D707}_{i},\unicode[STIX]{x1D70E})=\lim _{L\rightarrow \infty }\int _{0}^{L}g(r,\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716}_{i},\unicode[STIX]{x1D707}_{i})\frac{\text{d}l}{v_{\Vert }}\bigg/\penalty \exhyphenpenalty \int _{0}^{L}\frac{\text{d}l}{v_{\Vert }}. & & \displaystyle\end{eqnarray}$$

This is independent of $\unicode[STIX]{x1D6FC}$ , as the integral extends along a field line so passes many times around the torus on a flux surface, and can also be written in terms of the flux surface average,

(2.10) $$\begin{eqnarray}\displaystyle \overline{g}(r,\unicode[STIX]{x1D716}_{i},\unicode[STIX]{x1D707}_{i},\unicode[STIX]{x1D70E})=\left\langle \frac{Bg}{v_{\Vert }}\right\rangle \bigg/\penalty \exhyphenpenalty \left\langle \frac{B}{v_{\Vert }}\right\rangle . & & \displaystyle\end{eqnarray}$$

In the trapped region the integral is taken between consecutive bounce points, denoted $l_{1}$ and $l_{2}$ at which $B(r,\unicode[STIX]{x1D6FC},l_{1})=B(r,\unicode[STIX]{x1D6FC},l_{2})=1/\unicode[STIX]{x1D706}$ , so

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

where the bounce time $\unicode[STIX]{x1D70F}_{b}=\int _{l_{1}}^{l_{2}}\,\text{d}l/|v_{\Vert }|$ .

The odd piece of the distribution function is needed to evaluate the parallel friction in (2.7). It was determined in Helander et al. (Reference Helander, Parra and Newton2017b ) for a pure plasma, where it was used to evaluate the parallel ion flow. For convenience we outline the arguments here, as we will finally evaluate different velocity space averages of the distribution and account for an impurity species. Formally, the odd piece of the distribution follows from the line integral of the even (2.8),

(2.12) $$\begin{eqnarray}\displaystyle f_{i}^{-}(r,\unicode[STIX]{x1D6FC},l,\unicode[STIX]{x1D716}_{i},\unicode[STIX]{x1D707}_{i},\unicode[STIX]{x1D70E})=\int _{l_{0}}^{l}\left[C_{i}^{+}(f_{i})-\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{+}\right]\frac{\text{d}l^{\prime }}{v_{\Vert }}+X\left(r,\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D716}_{i},\unicode[STIX]{x1D707}_{i},\unicode[STIX]{x1D70E}\right). & & \displaystyle\end{eqnarray}$$

We will return to the definition of $l_{0}$ momentarily. Ruling out the collisional limit for the bulk ions, the odd (2.8) indicates that $f^{+}$ is a function of the constants of the motion, and the integration constant $X$ is then determined by the orbit average

(2.13) $$\begin{eqnarray}\displaystyle \overline{\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{-}}=\overline{C_{i}^{-}(f_{i})}. & & \displaystyle\end{eqnarray}$$

The orbit average of the even equation (2.8) constrains the even piece of the distribution function appearing above,

(2.14) $$\begin{eqnarray}\displaystyle \overline{\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{+}}=\overline{C_{i}^{+}(f_{i})}. & & \displaystyle\end{eqnarray}$$

This entails the assumptions on the quality of confinement noted in § 1. The ratio of the right to left-hand side of (2.14) is formally of the order $\unicode[STIX]{x1D708}_{\ast i}/\unicode[STIX]{x1D70C}_{\ast i}$ . In the $1/\unicode[STIX]{x1D708}$ regime, collisions are sufficiently dominant that the distribution function is nearly a Maxwellian and the derivation can proceed quite readily (Helander Reference Helander2014). At lower collisionality, orbit drifts can generate loss regions in velocity space, and the plasma is not generally in a local thermodynamic equilibrium (recent treatments of the behaviour of near-quasisymmetric systems can be found, for example, in Helander (Reference Helander2014), Calvo et al. (Reference Calvo, Parra, Velasco and Alonso2015)). Two limits in which confinement can be adequately restored are known (Ho & Kulsrud Reference Ho and Kulsrud1987; Calvo et al. Reference Calvo, Parra, Alonso and Velasco2014, Reference Calvo, Parra, Velasco and Alonso2015, Reference Calvo, Parra, Velasco and Alonso2017) – they were outlined in Helander et al. (Reference Helander, Parra and Newton2017b ) and can be summarised as follows. One is that in which the drift in the radial electric field, $\boldsymbol{v}_{E}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}\times \boldsymbol{b}/B$ , is sufficiently strong compared to the magnetic drift, $\boldsymbol{v}_{M}$ , that the bounce-averaged orbits stay close to a flux surface – this is often consistent with a large aspect ratio system. The other is when the orbit-averaged magnetic drift is small compared to the local value, which is achieved when a stellarator is optimised to be near omnigeneous. In both cases the distribution function is maintained near to Maxwellian, and is constant on a flux surface, thereby making the electrostatic potential a flux surface function, as assumed earlier.

We therefore assume here that either we are in the $1/\unicode[STIX]{x1D708}$ regime, or one of the above low collisionality conditions is satisfied. The even distribution can then be written in the form $f_{i}^{+}=F_{0}+F_{1}$ , where $F_{0}(\unicode[STIX]{x1D716}_{i},r)$ is a Maxwellian, and $F_{1}\ll F_{0}$ , remembering that it is constant along field lines, so is independent of $l$ in the trapped region of phase space, and independent of $\unicode[STIX]{x1D6FC}$ and $l$ in the circulating region. As the averaged drift $\overline{\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r}(\unicode[STIX]{x2202}_{r}F_{0})=0$ in the circulating region, it was argued in Helander et al. (Reference Helander, Parra and Newton2017b ) that $F_{1}$ is small in the circulating region, compared to its value in the trapped region, and we will neglect it. In the $1/\unicode[STIX]{x1D708}$ regime, $F_{1}=0$ also in the trapped region. In the lower collisionality regimes, the orbit average (2.14) requires $\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}F_{1})+\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r}(\unicode[STIX]{x2202}_{r}F_{0})\approx 0$ . (The resolution of the behaviour of the distribution in the trapped–passing boundary layer is required to evaluate the bulk ion transport (Ho & Kulsrud Reference Ho and Kulsrud1987), but is not needed here.) The explicit drift term in (2.12) can then be conveniently written for the low collisionality regimes of interest here as

(2.15) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i}^{+}=\left(\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r-\unicode[STIX]{x1D700}_{t}\overline{\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r}\right)\frac{\unicode[STIX]{x2202}f_{Mi}}{\unicode[STIX]{x2202}r}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D700}_{t}=0$ in the $1/\unicode[STIX]{x1D708}$ regime, and $\unicode[STIX]{x1D700}_{t}=1$ in the trapped region of phase space and zero otherwise in the $\sqrt{\unicode[STIX]{x1D708}}$ regime.

Now we consider the integration constant $X$ . As the odd piece of the distribution function must vanish at a bounce point, if we choose $l_{0}$ in (2.12) to be such a point, then $X=0$ in the trapped region. Therefore, we set:

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

In the circulating region, $X$ is set by the constraint equation (2.13). Using the conservative form of the particle drift, equation (2.3), along with the condition that circulating particles do not drift from their flux surfaces on average, it was shown in detail in Helander et al. (Reference Helander, Parra and Newton2017b ) that this constraint reduces to the following familiar form, for the low collisionality regimes of interest,

(2.17) $$\begin{eqnarray}\displaystyle \left\langle \frac{B}{v_{\Vert }}C_{i}^{-}(f_{i})\right\rangle =0. & & \displaystyle\end{eqnarray}$$

In the next section we introduce a model collision operator which allows the integration constant to be determined explicitly, using this constraint. We will then have, with (2.15), the expression for $f_{i}^{-}$ needed to evaluate the moments giving the bulk ion flow and the impurity flux.

2.3 Bulk ion collision operator

The differences in the bulk ion flow in a pure plasma which result from using different forms of the collision operator to determine the odd piece of the distribution were discussed in Helander et al. (Reference Helander, Parra and Newton2017b ). Similar considerations apply when evaluating particle fluxes via (2.4) and (2.7). It is well known that a momentum conserving collision operator is at least required to maintain the intrinsic ambipolarity of transport driven by friction (see, for example, Sugama & Nishimura (Reference Sugama and Nishimura2002), Maaßberg, Beidler & Turkin (Reference Maaßberg, Beidler and Turkin2009)). Therefore, we adopt here the following description of the bulk ion collisions.

Due to the disparate ion masses, we use a common approximation to the bulk ion–impurity collision operator $C_{iz}$  (Rosenbluth, Hazeltine & Hinton Reference Rosenbluth, Hazeltine and Hinton1972; Helander & Sigmar Reference Helander and Sigmar2002),

(2.18) $$\begin{eqnarray}\displaystyle C_{iz}(f_{i})=\unicode[STIX]{x1D708}_{D}^{iz}(v)\left({\mathcal{L}}(f_{i})+\frac{m_{i}v_{\Vert }V_{z\Vert }}{T}f_{Mi}\right). & & \displaystyle\end{eqnarray}$$

The pitch angle scattering operator ${\mathcal{L}}=(1/2)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}(1-\unicode[STIX]{x1D709}^{2})\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}$ , where $\unicode[STIX]{x1D709}=\cos \unicode[STIX]{x1D703}=v_{\Vert }/v$ is the cosine of the particle pitch angle. With the normalised velocity $x_{a}=v/v_{Ta}$ , the deflection frequency $\unicode[STIX]{x1D708}_{D}^{iz}(v)=3\unicode[STIX]{x03C0}^{1/2}/4\unicode[STIX]{x1D70F}_{iz}x_{i}^{3}=\hat{\unicode[STIX]{x1D708}}_{D}^{iz}/x_{i}^{3}$ and the collision time $\unicode[STIX]{x1D70F}_{iz}=3(2\unicode[STIX]{x03C0})^{3/2}\sqrt{m_{i}}T^{3/2}\unicode[STIX]{x1D716}_{0}^{2}/n_{z}Z^{2}e^{4}\ln \unicode[STIX]{x1D6EC}$ . The parallel impurity flow, $V_{z\Vert }$ , will be determined in the next section. Bulk ion self-collisions are described by an operator with a similar structure (Rosenbluth et al. Reference Rosenbluth, Hazeltine and Hinton1972; Connor Reference Connor1973), that is, a combination of pitch angle scattering and a momentum restoring term,

(2.19) $$\begin{eqnarray}\displaystyle C_{ii}(f_{i})=\unicode[STIX]{x1D708}_{D}^{ii}(v)\left({\mathcal{L}}(f_{i})+\frac{m_{i}v_{\Vert }{\mathcal{V}}_{i\Vert }}{T}f_{Mi}\right). & & \displaystyle\end{eqnarray}$$

The full energy-dependent deflection frequency $\unicode[STIX]{x1D708}_{D}^{ii}(v)=\hat{\unicode[STIX]{x1D708}}_{D}^{ii}[\unicode[STIX]{x1D719}(x_{i})-G(x_{i})]/x_{i}^{3}$ , $\hat{\unicode[STIX]{x1D708}}_{D}^{ii}$ is defined in analogy to $\hat{\unicode[STIX]{x1D708}}_{D}^{iz}$ , the error function $\unicode[STIX]{x1D719}(x)=(2/\sqrt{\unicode[STIX]{x03C0}})\int _{0}^{x}\text{e}^{-y^{2}}\,\text{d}y$ and the Chandrasekhar function $G(x)=[\unicode[STIX]{x1D719}(x)-x\unicode[STIX]{x1D719}^{\prime }(x)]/2x^{2}$ . The momentum restoring coefficient ${\mathcal{V}}_{i\Vert }$ will be set by requiring momentum conservation in bulk ion self-collisions, $\int v_{\Vert }C_{ii}^{-}(f_{i})\,\text{d}^{3}v=0$ . Altogether our model bulk ion collision operator is $C_{i}=C_{ii}+C_{iz}$ , and we introduce the total collision frequency $\unicode[STIX]{x1D708}_{D}^{i}(v)=\unicode[STIX]{x1D708}_{D}^{ii}+\unicode[STIX]{x1D708}_{D}^{iz}$ . The bulk ion flow was evaluated in Helander et al. (Reference Helander, Parra and Newton2017b ) for the case $C_{i}=C_{ii}$ , with $n_{z}=0$ , and as expected many similar steps appear in the derivation here. We highlight throughout the changes introduced by allowing for an impurity species.

For convenience we can set the electrostatic potential to zero on the surface of interest, and use the velocity space coordinate $\unicode[STIX]{x1D706}=v_{\bot }^{2}/v^{2}B$ , which satisfies $\unicode[STIX]{x1D735}_{\Vert }|_{\unicode[STIX]{x1D716},\unicode[STIX]{x1D707}}\unicode[STIX]{x1D706}=0$ . The pitch angle scattering operator can be written as ${\mathcal{L}}=(2\unicode[STIX]{x1D709}/B)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D706}}(\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D706}})$ and $\unicode[STIX]{x1D709}=\pm \sqrt{1-\unicode[STIX]{x1D706}B}$ . The passing region constraint equation (2.17) is then

(2.20) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\unicode[STIX]{x1D706}\left\langle \sqrt{1-\unicode[STIX]{x1D706}B}\frac{\unicode[STIX]{x2202}f_{i}^{-}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right\rangle +\frac{m_{i}v}{2T}\frac{\left\langle (\unicode[STIX]{x1D708}_{D}^{ii}{\mathcal{V}}_{i\Vert }+\unicode[STIX]{x1D708}_{D}^{iz}V_{z\Vert })B\right\rangle }{\unicode[STIX]{x1D708}_{D}^{i}(v)}f_{Mi}=0. & & \displaystyle\end{eqnarray}$$

Integrating over $\unicode[STIX]{x1D706}$ , with $\unicode[STIX]{x1D706}<1/B_{\text{max}}$ , the integration constant vanishes upon requiring regularity at $\unicode[STIX]{x1D706}=0$ . We can now insert the general form for $f_{i}^{-}$ from (2.12), noting that $F_{1}$ is taken to be negligible in the passing region, so the contribution from the term explicitly involving the collision operator vanishes. We thus obtain a simple extension to equation (4.7) of Helander et al. (Reference Helander, Parra and Newton2017b ) to account for the presence of an impurity species,

(2.21) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}X}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=-\frac{m_{i}v}{2T_{i}\left\langle \unicode[STIX]{x1D709}\right\rangle }\left(\frac{\left\langle (\unicode[STIX]{x1D708}_{D}^{ii}{\mathcal{V}}_{i\Vert }+\unicode[STIX]{x1D708}_{D}^{iz}V_{z\Vert })B\right\rangle }{\unicode[STIX]{x1D708}_{D}^{i}(v)}f_{Mi}+\frac{T}{e_{i}}\left\langle g_{4}\right\rangle \frac{\unicode[STIX]{x2202}f_{Mi}}{\unicode[STIX]{x2202}r}\right), & & \displaystyle\end{eqnarray}$$

where the contribution from the drift term in (2.12) gave rise to the known geometry function (Nakajima et al. Reference Nakajima, Okamoto, Todoroki, Nakamura and Wakatani1989; Helander, Geiger & Maaßberg Reference Helander, Geiger and Maaßberg2011)

(2.22) $$\begin{eqnarray}\displaystyle g_{4}(\unicode[STIX]{x1D706},l)=\unicode[STIX]{x1D709}\int _{l_{\text{max}}}^{l}(\boldsymbol{b}\times \unicode[STIX]{x1D735}r)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}^{-1}\,\text{d}l^{\prime }, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D706}<1/B_{\text{max}}$ and $B(l_{\text{max}})=B_{\text{max}}$ . The full form for the integration constant $X$ in the bulk ion distribution equation (2.12) is thus given by (2.21), for $\unicode[STIX]{x1D706}<1/B_{\text{max}}$ , and $X=0$ , for $1/B_{\text{max}}<\unicode[STIX]{x1D706}<1/B_{\text{min}}$ , where $B_{\text{min}}$ is the minimum field strength on the flux surface.

The momentum restoring coefficient, ${\mathcal{V}}_{i\Vert }$ , is determined by momentum conservation in bulk ion self-collisions,

(2.23) $$\begin{eqnarray}\displaystyle 0=m_{i}\int v_{\Vert }C_{ii}^{-}(f_{i})\,\text{d}^{3}v=m_{i}\int v_{\Vert }\unicode[STIX]{x1D708}_{D}^{ii}\left({\mathcal{L}}(f_{i}^{-})+\frac{m_{i}v_{\Vert }{\mathcal{V}}_{i\Vert }}{T}f_{Mi}\right)\text{d}^{3}v, & & \displaystyle\end{eqnarray}$$

as the self-adjoint property of the Lorentz operator gives

(2.24) $$\begin{eqnarray}\displaystyle {\mathcal{V}}_{i\Vert }=\frac{1}{n_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\int \unicode[STIX]{x1D708}_{D}^{ii}v_{\Vert }f_{i}^{-}\,\text{d}^{3}v. & & \displaystyle\end{eqnarray}$$

Here we have introduced the velocity space average (Hirshman Reference Hirshman1976) for a function of the magnitude of the velocity, $\{F(v)\}=(8/3\sqrt{\unicode[STIX]{x03C0}})\int _{0}^{\infty }F(x)x^{4}\text{e}^{-x^{2}}\,\text{d}x$ , so $\{\unicode[STIX]{x1D708}_{D}^{ii}\}\unicode[STIX]{x1D70F}_{ii}=\sqrt{2}-\ln (1+\sqrt{2})$ . Inserting $f_{i}^{-}$ from (2.12), we see as detailed in Helander et al. (Reference Helander, Parra and Newton2017b ) that the term explicitly containing $C^{+}(f_{i})$ does not contribute when the collision operator is of the form assumed here. The explicit drift term is usefully written in terms of the function $u$ defined in § 2.1, using the projection $\boldsymbol{v}_{di}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r$ of (2.3), and results in the same contribution as in Helander et al. (Reference Helander, Parra and Newton2017b ), with the integration constant in $u$ fixed by taking $u=0$ where $B=B_{\text{max}}$ . The appearance of the impurity flow term in the integration constant here, however, gives an additional contribution compared to equation (4.12) of Helander et al. (Reference Helander, Parra and Newton2017b ),

(2.25) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle -\frac{1}{n_{i}\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\left\langle B^{2}\int _{0}^{\infty }\,\text{d}v\,2\unicode[STIX]{x03C0}v^{2}\unicode[STIX]{x1D708}_{D}^{ii}\int _{0}^{1/B_{\text{max}}}\unicode[STIX]{x1D706}\frac{\unicode[STIX]{x2202}X}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\,\text{d}\unicode[STIX]{x1D706}\right\rangle \nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad =\frac{f_{c}}{\{\unicode[STIX]{x1D708}_{D}^{ii}\}}\left(\left\{\frac{{\unicode[STIX]{x1D708}_{D}^{ii}}^{2}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\left\langle B{\mathcal{V}}_{i\Vert }\right\rangle +\left\{\frac{{\unicode[STIX]{x1D708}_{D}^{iz}}^{2}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\left\langle BV_{z\Vert }\right\rangle \right)+\frac{f_{s}T}{e}(A_{1i}-\unicode[STIX]{x1D702}A_{2i}),\end{eqnarray}$$

where it has been anticipated that we will only need the restoring coefficient in the form $\left\langle B{\mathcal{V}}_{i\Vert }\right\rangle$ . So we find the following modification of equation (6.8) of Helander et al. (Reference Helander, Parra and Newton2017b ) in the presence of an impurity species,

(2.26) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \left\langle B{\mathcal{V}}_{i\Vert }\right\rangle \left[1-\frac{f_{c}}{\left\{\unicode[STIX]{x1D708}_{D}^{ii}\right\}}\left\{\frac{{\unicode[STIX]{x1D708}_{D}^{ii}}^{2}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\right]-\frac{f_{c}}{\left\{\unicode[STIX]{x1D708}_{D}^{ii}\right\}}\left\{\frac{\unicode[STIX]{x1D708}_{D}^{ii}\unicode[STIX]{x1D708}_{D}^{iz}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\left\langle BV_{z\Vert }\right\rangle \nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad =\frac{T}{e_{i}}\left(A_{1i}-\unicode[STIX]{x1D702}A_{2i}\right)\left[f_{s}+\left\langle (u+s)B^{2}\right\rangle \right],\end{eqnarray}$$

which reduces to that expression in the limit of a pure plasma, where $n_{z}\rightarrow 0$ and $\unicode[STIX]{x1D708}_{D}^{i}\rightarrow \unicode[STIX]{x1D708}_{D}^{ii}$ . Here $\unicode[STIX]{x1D702}=\{\unicode[STIX]{x1D708}_{D}^{ii}(5/2-x^{2})\}/\{\unicode[STIX]{x1D708}_{D}^{ii}\}=(5/2)-1/[2-\sqrt{2}\ln (1+\sqrt{2})]$ ,

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

and the term $s$ is zero in the $1/\unicode[STIX]{x1D708}$ regime, and given in the $\sqrt{\unicode[STIX]{x1D708}}$ regime by

(2.28) $$\begin{eqnarray}\displaystyle s(l)=\frac{3}{2}\int _{l_{\text{max}}}^{l}\text{d}l^{\prime }\int _{1/B_{\text{max}}}^{1/B(l^{\prime })}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D709}(l^{\prime })}\overline{\unicode[STIX]{x1D709}\left(\boldsymbol{b}\times \unicode[STIX]{x1D735}r\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{\unicode[STIX]{x1D709}}{B}\right)}. & & \displaystyle\end{eqnarray}$$

Finally, with the assumed quality of confinement described in § 2.2 (that is $F_{0}\approx f_{Mi}$ ) and the model collision operator, equation (2.18), then adopted here, the parallel friction in (2.7) needed to determine the particle flux takes the form

(2.29) $$\begin{eqnarray}\displaystyle R_{zi\Vert }=m_{i}\int \unicode[STIX]{x1D708}_{D}^{iz}(v)v_{\Vert }f_{i}^{-}\,\text{d}^{3}v-\frac{m_{i}n_{i}}{\unicode[STIX]{x1D70F}_{iz}}V_{z\Vert }. & & \displaystyle\end{eqnarray}$$

We therefore now need an expression for the parallel impurity flow and in the next section we consider the impurity distribution function.

2.4 Impurity distribution function

As introduced in § 2.1, the collisional impurity species can be treated by the usual expansion of (2.5) in the small parameter $1/\unicode[STIX]{x1D708}_{\ast zz}$  (Braun & Helander Reference Braun and Helander2010), allowing for $\unicode[STIX]{x1D708}_{zz}\sim \unicode[STIX]{x1D708}_{zi}$ . At order $-1$ , $C_{z}(f_{z1}^{(-1)})=0$ , so the impurity distribution has the form of a perturbed Maxwellian (Helander & Sigmar Reference Helander and Sigmar2002),

(2.30) $$\begin{eqnarray}\displaystyle f_{z1}^{(-1)}=\left[\frac{p_{z}^{(-1)}}{p_{z}}+\frac{m_{z}}{T}v_{\Vert }V_{z\Vert }^{(-1)}+\left(x_{z}^{2}-\frac{5}{2}\right)\frac{T_{z}^{(-1)}}{T}\right]f_{z0}. & & \displaystyle\end{eqnarray}$$

The parallel flow, $V_{z\Vert }^{(-1)}$ , is constrained by momentum conservation in this order

(2.31) $$\begin{eqnarray}\displaystyle \int m_{z}v_{\Vert }C_{zi}(f_{z1}^{(-1)})\,\text{d}^{3}v=0, & & \displaystyle\end{eqnarray}$$

and we take the disparate mass form for the collision operator $C_{zi}$  (Helander & Sigmar Reference Helander and Sigmar2002; Hazeltine & Meiss Reference Hazeltine and Meiss2003),

(2.32) $$\begin{eqnarray}\displaystyle C_{zi}\left(f_{z1}\right)=-\frac{\boldsymbol{R}_{zi}}{m_{z}n_{z}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f_{z0}}{\unicode[STIX]{x2202}\boldsymbol{v}}+\frac{m_{i}n_{i}}{m_{z}n_{z}\unicode[STIX]{x1D70F}_{iz}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\boldsymbol{\cdot }\left[(\boldsymbol{v}-\mathbf{V}_{z})f_{z1}+\frac{T}{m_{z}}\frac{\unicode[STIX]{x2202}f_{z1}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right]. & & \displaystyle\end{eqnarray}$$

Using this in (2.31) gives simply $R_{zi\Vert }^{(-1)}=0$ . Considering the balance in (2.29), this would require an impurity flow $V_{z\Vert }^{(-1)}\sim \unicode[STIX]{x1D70C}_{\ast i}v_{Ti}$ here, that is $V_{z\Vert }^{(-1)}/v_{Tz}\sim \unicode[STIX]{x1D70C}_{\ast z}Z$ . However, by definition of the collisional expansion, $V_{z\Vert }^{(-1)}\sim \unicode[STIX]{x1D70C}_{\ast z}v_{Tz}\unicode[STIX]{x1D708}_{\ast zz}$ , that is $V_{z\Vert }^{(-1)}/v_{Tz}\sim \unicode[STIX]{x1D70C}_{\ast z}Z(\unicode[STIX]{x1D708}_{\ast iz}n_{z}Z/n_{i})$ . For the collisionless ions $\unicode[STIX]{x1D708}_{\ast iz}\ll 1$ , so restricting the impurity density such that $\unicode[STIX]{x1D708}_{\ast iz}n_{z}Z/n_{i}<1$ holds (hence we do not consider a pure ‘impurity’ plasma) the two conditions would give a contradiction. This is resolved by requiring $V_{z\Vert }^{(-1)}=0$ , and so $R_{zi\Vert }^{(0)}$ is found to be the leading-order friction driving the particle flux. (Note this can be compared to the derivation given by Braun & Helander (Reference Braun and Helander2010) for the fully collisional case, where $V_{z\Vert }^{(-1)}=0$ still results.)

The form of the leading-order flow, $V_{z\Vert }\approx V_{z\Vert }^{(0)}$ , may also be found as usual by considering density conservation from the $v_{\Vert }/B$ moment of (2.5) written in conservative form using (2.3), or radial force balance combined with incompressibility of the equilibrium flow in leading order:

(2.33) $$\begin{eqnarray}\displaystyle V_{z\Vert }=\left(\frac{1}{Zen_{z}}\frac{\text{d}p_{z}}{\text{d}r}+\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}r}\right)uB+\frac{K_{z}(r)B}{n_{z}}. & & \displaystyle\end{eqnarray}$$

A constraint on the flux surface function $K_{z}$ is obtained from the Spitzer-type problem for $f_{z1}^{(0)}$ arising at zeroth order in the collisional expansion of (2.5),

(2.34) $$\begin{eqnarray}\displaystyle C_{z}(f_{z1}^{(0)})=v_{\Vert }f_{z0}\left[A_{z1\Vert }^{(-1)}+\left(x_{z}^{2}-\frac{5}{2}\right)A_{z2\Vert }^{(-1)}\right], & & \displaystyle\end{eqnarray}$$

where the parallel driving forces resulting from $f_{z1}^{(-1)}$ are

(2.35a,b ) $$\begin{eqnarray}\displaystyle A_{z1\Vert }^{(-1)}=\frac{\unicode[STIX]{x1D735}_{\Vert }p_{z1}^{(-1)}}{p_{z}}+\frac{Ze}{T}\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D719}_{1}^{(-1)},\quad A_{z2\Vert }^{(-1)}=\frac{\unicode[STIX]{x1D735}_{\Vert }T_{z1}^{(-1)}}{T}. & & \displaystyle\end{eqnarray}$$

Parallel momentum conservation, that is the $m_{z}v_{\Vert }$ moment of (2.34), gives

(2.36) $$\begin{eqnarray}\displaystyle R_{zi\Vert }^{(0)}=n_{z}TA_{z1\Vert }^{(-1)}. & & \displaystyle\end{eqnarray}$$

Upon taking the $B$ -weighted flux surface average, the general property of the divergence of a vector field $\boldsymbol{F}$ ,

(2.37) $$\begin{eqnarray}\displaystyle \left\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{F}\right\rangle =\frac{1}{V^{\prime }(r)}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left\langle V^{\prime }(r)\boldsymbol{F}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r\right\rangle , & & \displaystyle\end{eqnarray}$$

where $V$ is the volume enclosed by a flux surface, annihilates the parallel gradient terms and sets the constraint,

(2.38) $$\begin{eqnarray}\displaystyle \left\langle BR_{zi\Vert }^{(0)}\right\rangle =0. & & \displaystyle\end{eqnarray}$$

This relation was first discussed in the context of transport in the mixed collisionality regime of a tokamak in Hirshman (Reference Hirshman1976). Applying this to (2.29) results in

(2.39) $$\begin{eqnarray}\displaystyle \left\langle BV_{z\Vert }\right\rangle =\frac{T}{Ze}A_{1z}\left\langle uB^{2}\right\rangle +\frac{K_{z}(r)}{n_{z}}\left\langle B^{2}\right\rangle =\frac{\unicode[STIX]{x1D70F}_{iz}}{n_{i}}\left\langle B\int \unicode[STIX]{x1D708}_{D}^{iz}(v)v_{\Vert }f_{i}^{-}\,\text{d}^{3}v\right\rangle . & & \displaystyle\end{eqnarray}$$

We will find in the following section that we do not need to solve explicitly for the function $K_{z}$ to determine the particle flux. Note that, in combination with the choice of the momentum restoring form for the bulk ion collision operator, this leads us to expect that the analytical result obtained here for the impurity flux will be in good agreement with the result of the DKES code when the momentum conserving correction is applied (Hirshman et al. Reference Hirshman, Shaing, van Rij, Beasley and Crume1986; van Rij & Hirshman Reference van Rij and Hirshman1989; Maaßberg et al. Reference Maaßberg, Beidler and Turkin2009). This was seen already in the numerical results presented in Helander et al. (Reference Helander, Newton, Mollén and Smith2017a ).

3.1 Two collisional impurities

There are typically many impurity species present in magnetically confined fusion plasmas. A common situation is one in which there are trace amounts of a particularly heavy impurity, often released from the exhaust region, in a background of an otherwise dominant impurity, which may be released for example from the main walls. The transport of the heavier impurity is of particular importance, as it will be the most difficult to ionise and thus poses the strongest potential source of core radiation losses. The results presented above allow us to make the following interesting observation when both impurity species are taken to be collisional, extending somewhat the analysis presented for the tokamak in Burrell & Wong (Reference Burrell and Wong1981).

We denote the lighter impurity by a subscript $A$ here, with charge $Z_{A}\gg 1$ , and continue to use $z$ for the heavier impurity. Following Braun & Helander (Reference Braun and Helander2010), $V_{z\Vert }^{(-1)}=V_{A\Vert }^{(-1)}=0$ as the species are collisional, and a flow cannot be driven at this order through interaction with the collisionless bulk. Also, as species $A$ is collisional, the radial flux of species $z$ will continue to be dominated by the friction drive, as long as (2.6) is satisfied, so

(3.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{z}=\frac{1}{Ze}\left\langle uBR_{z\Vert }\right\rangle =\frac{1}{Ze}\left\langle uB(R_{zi\Vert }+R_{zA\Vert })\right\rangle . & & \displaystyle\end{eqnarray}$$

Assuming that the bulk ions and species $A$ have disparate masses, $m_{i}\ll m_{A}$ , collisions between them can be modelled by a collision operator analogous to that in (2.18). The contribution to the impurity flux from $R_{zi\Vert }=-R_{iz\Vert }$ can then be determined as a simple extension of the results above – we will again obtain (3.1) and (3.3), but with the flux function $P(r)$ modified such that $\unicode[STIX]{x1D708}_{D}^{i}\rightarrow \unicode[STIX]{x1D708}_{D}^{ii}+\unicode[STIX]{x1D708}_{D}^{iz}+\unicode[STIX]{x1D708}_{D}^{iA}$ and

(3.10) $$\begin{eqnarray}\displaystyle \left\{\frac{{\unicode[STIX]{x1D708}_{D}^{iz}}^{2}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\left\langle BV_{z\Vert }\right\rangle \rightarrow \left\{\frac{{\unicode[STIX]{x1D708}_{D}^{iz}}^{2}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\left\langle BV_{z\Vert }\right\rangle +\left\{\frac{\unicode[STIX]{x1D708}_{D}^{iz}\unicode[STIX]{x1D708}_{D}^{iA}}{\unicode[STIX]{x1D708}_{D}^{i}}\right\}\left\langle BV_{A\Vert }\right\rangle . & & \displaystyle\end{eqnarray}$$

Parallel momentum constraints analogous to (2.38) are obtained similarly for the two collisional species,

(3.11) $$\begin{eqnarray}\displaystyle \left\langle B(R_{zi\Vert }+R_{zA\Vert })\right\rangle =0 & & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle \left\langle B(R_{Ai\Vert }+R_{Az\Vert })\right\rangle =0. & & \displaystyle\end{eqnarray}$$

The first of these allows us to again eliminate $G_{3}$ from (3.4) leaving

(3.13) $$\begin{eqnarray}\displaystyle R_{zi\Vert }=G_{1}(r)\left(uB-\left\langle uB^{2}\right\rangle \frac{B}{\left\langle B^{2}\right\rangle }\right)+G_{2}(r)\left(sB-\left\langle sB^{2}\right\rangle \frac{B}{\left\langle B^{2}\right\rangle }\right)-\frac{B}{\left\langle B^{2}\right\rangle }\left\langle BR_{zA\Vert }\right\rangle . & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Note that the total radial impurity current is

(3.14) $$\begin{eqnarray}\displaystyle J_{\text{imp}}\equiv e_{z}\unicode[STIX]{x1D6E4}_{z}+e_{A}\unicode[STIX]{x1D6E4}_{A}=-\left\langle uB(R_{iz\Vert }+R_{iA\Vert })\right\rangle . & & \displaystyle\end{eqnarray}$$

The disparate mass collision operator adopted will lead to an expression for $R_{iA\Vert }=-R_{Ai\Vert }$ analogous to (3.3). Summing the two constraints in (3.12) gives $\left\langle B(R_{iz\Vert }+R_{iA\Vert })\right\rangle =0$ , which allows all of the unknown flux functions $\left\langle B{\mathcal{V}}_{i\Vert }\right\rangle$ , $K_{A}$ and $K_{z}$ to again be eliminated from the flux, leaving

(3.15) $$\begin{eqnarray}\displaystyle J_{\text{imp}} & = & \displaystyle -\frac{m_{i}p_{i}}{e\unicode[STIX]{x1D70F}_{iz}}\left[\left(\frac{A_{1z}}{Z}+\frac{1}{\unicode[STIX]{x1D701}_{A}}\frac{A_{1A}}{Z_{A}}\right)\left(\left\langle u^{2}B^{2}\right\rangle -\frac{\left\langle uB^{2}\right\rangle ^{2}}{\left\langle B^{2}\right\rangle }\right)\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\left(1+\frac{1}{\unicode[STIX]{x1D701}_{A}}\right)\!\left(A_{1i}-\frac{3}{2}A_{2i}\right)\!\left(\left\langle u(u+s)B^{2}\right\rangle -\left\langle (u+s)B^{2}\right\rangle \frac{\left\langle uB^{2}\right\rangle }{\left\langle B^{2}\right\rangle }\right)\right]\!,\end{eqnarray}$$

where $\unicode[STIX]{x1D701}_{A}=n_{z}Z^{2}/n_{A}Z_{A}^{2}$ . The total impurity current can thus also experience temperature screening, under the conditions described in the previous section.

To form the explicit expression for the flux of the heavier impurity, we still need to determine the combination

(3.16) $$\begin{eqnarray}\displaystyle \left\langle uBR_{zA\Vert }\right\rangle -\frac{\left\langle uB^{2}\right\rangle }{\left\langle B^{2}\right\rangle }\left\langle BR_{zA\Vert }\right\rangle , & & \displaystyle\end{eqnarray}$$

where the friction $R_{zA\Vert }=m_{z}\int v_{\Vert }C_{zA}(f_{z},f_{A})\,\text{d}^{3}v$ contains the linearised collision operator $C_{zA}$ acting on the distribution functions of the two collisional species. These are given to leading order by the solution of (2.34) and the analogous equation for $f_{A1}^{(0)}$ . The solution can be written as an expansion in Sonine polynomials, $L_{\unicode[STIX]{x1D6FC}}^{(3/2)}(x^{2})$ , such that

(3.17) $$\begin{eqnarray}\displaystyle f_{a1}^{(0)}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}=0}^{2}u_{a\unicode[STIX]{x1D6FC}}L_{\unicode[STIX]{x1D6FC}}^{(3/2)}(x_{a}^{2})\frac{m_{a}v_{\Vert }}{T}f_{Ma}. & & \displaystyle\end{eqnarray}$$

With $L_{0}^{(3/2)}(x^{2})=1$ and $L_{1}^{(3/2)}(x^{2})=(5/2)-x^{2}$ , the expansion coefficients may be recognised as $u_{a0}=V_{a\Vert }$ and $u_{a1}=-2q_{a\Vert }/5p_{a}$ , where $q_{a\Vert }$ is the parallel heat flux. Substituting this expansion into $R_{zA\Vert }$ , the integration over the collision operator may be performed directly (Helander & Sigmar Reference Helander and Sigmar2002), and the parallel friction coefficients will depend on the mass ratio of the impurities. We treat the case of disparate impurity masses, $m_{A}\ll m_{z}$ and $Z_{A}\ll Z$ , explicitly here, which may give a good approximation to the experimentally relevant case of a low collisionality bulk H plasma with a main impurity such as C from the walls, and a low density, heavier component, such as Fe. Then

(3.18) $$\begin{eqnarray}\displaystyle R_{zA\Vert }=\frac{m_{i}n_{i}}{\unicode[STIX]{x1D70F}_{iz}}\sqrt{\frac{m_{A}}{m_{i}}}\frac{n_{A}Z_{A}^{2}}{n_{i}}\left(V_{A\Vert }-V_{z\Vert }-\frac{3}{5}\frac{q_{A\Vert }}{p_{A}}+\frac{15}{8}u_{A2}\right). & & \displaystyle\end{eqnarray}$$

The $v_{\Vert }L_{2}$ -moment of (2.34) for species $A$ relates the coefficient $u_{A2}$ to the parallel flows in the disparate mass limit,

(3.19) $$\begin{eqnarray}\displaystyle u_{A2}=\frac{1}{\left(45\sqrt{2}+433\unicode[STIX]{x1D701}_{A}/4\right)}\left[30\unicode[STIX]{x1D701}_{A}\left(V_{z\Vert }-V_{A\Vert }\right)+\frac{2}{5}\left(12\sqrt{2}+69\unicode[STIX]{x1D701}_{A}\right)\frac{q_{A\Vert }}{p_{A}}\right]. & & \displaystyle\end{eqnarray}$$

The parallel species flows in (3.18), $V_{z\Vert }$ and $V_{A\Vert }$ , have the general form of (2.33). In the combination of (3.16), all terms containing the flux functions $K_{z}$ and $K_{A}$ cancel, leaving only contributions to the impurity flux from the radial gradients $A_{1z}$ and $A_{1A}$ . The form of the parallel heat flows in (3.18) can be determined using the $v_{\Vert }\unicode[STIX]{x1D716}_{z}/B$ moment of the conservative form of (2.5) (and the analogous equations for species $A$ and $i$ ), which gives the equation of energy conservation for a species,

(3.20) $$\begin{eqnarray}\displaystyle B\unicode[STIX]{x1D735}_{\Vert }\left(\frac{q_{a\Vert }}{B}-\frac{5}{2}\frac{p_{a}T}{e_{a}}A_{2a}u\right)=\int \frac{m_{a}v^{2}}{2}C_{a}(f_{a1})\,\text{d}^{3}v. & & \displaystyle\end{eqnarray}$$

The energy exchange between species appearing on the right-hand side competes with the parallel heat flux to determine the parallel temperature perturbation on a flux surface. It is typified here for disparate mass species using the second term of (2.32), giving for example

(3.21) $$\begin{eqnarray}\displaystyle \int \frac{m_{z}v^{2}}{2}C_{zi}(f_{z})\,\text{d}^{3}v=\frac{3m_{i}n_{i}}{m_{z}\unicode[STIX]{x1D70F}_{iz}}\left(T_{i}-T_{z}\right), & & \displaystyle\end{eqnarray}$$

between the heaviest impurity and the bulk ions (remember to leading order here the ion temperatures are equal, which leaves only the perturbed temperatures in this expression). The $v_{\Vert }L_{1}$ and $v_{\Vert }L_{2}$ -moments of (2.34) give us, in the disparate mass case, the heavy impurity parallel heat flux

(3.22) $$\begin{eqnarray}\displaystyle q_{z\Vert }=-\frac{125\sqrt{2}}{32}\frac{p_{z}^{2}\unicode[STIX]{x1D70F}_{zz}}{n_{z}m_{z}}A_{z2\Vert }^{(-1)}. & & \displaystyle\end{eqnarray}$$

Thus we see that energy exchange with the low collisionality bulk ions can only be neglected when $1\gg (n_{i}/n_{z})\sqrt{m_{i}/m_{z}}\unicode[STIX]{x1D708}_{\ast zz}\unicode[STIX]{x1D708}_{\ast iz}$ , which cannot be satisfied consistently with the condition (2.6). This arises similarly for energy exchange between the impurity species $A$ and the bulk ions. Energy exchange between the collisional impurity species is dominant when $1\ll (Z^{2}/Z_{A}^{2})\sqrt{m_{A}/m_{z}}\unicode[STIX]{x1D708}_{\ast zz}\unicode[STIX]{x1D708}_{\ast AA}$ , which will always be satisfied. Therefore we take the perturbed temperature of each impurity species to be equal, and set by energy exchange to that of the collisionless bulk ions. The parallel temperature gradients will then be negligible, and so the parallel impurity heat fluxes can be neglected in the expressions above.

The flux of the heaviest impurity can now be constructed from (3.9), (3.13), and (3.18), with (3.19) and the parallel flows just discussed, giving the final form

(3.23) $$\begin{eqnarray}\displaystyle \left\langle \unicode[STIX]{x1D6E4}_{z}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r\right\rangle & = & \displaystyle -\frac{m_{i}p_{i}}{Ze^{2}\unicode[STIX]{x1D70F}_{iz}}\left\{\!\!\left[\left(1+\sqrt{\frac{m_{A}}{m_{i}}}\frac{n_{A}Z_{A}^{2}}{n_{i}}(1-Y)\right)\frac{A_{1z}}{Z}-\sqrt{\frac{m_{A}}{m_{i}}}\frac{n_{A}Z_{A}^{2}}{n_{i}}(1-Y)\frac{A_{1A}}{Z_{A}}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \times \left(\left\langle u^{2}B^{2}\right\rangle -\frac{\left\langle uB^{2}\right\rangle ^{2}}{\left\langle B^{2}\right\rangle }\right)\nonumber\\ \displaystyle & & \displaystyle -\left.\left(A_{1i}-\frac{3}{2}A_{2i}\right)\left(\left\langle u(u+s)B^{2}\right\rangle -\left\langle (u+s)B^{2}\right\rangle \frac{\left\langle uB^{2}\right\rangle }{\left\langle B^{2}\right\rangle }\right)\right\},\end{eqnarray}$$

where $Y=225\unicode[STIX]{x1D701}_{A}/(180\sqrt{2}+433\unicode[STIX]{x1D701}_{A})$ . The net drive from the electric field still vanishes in the $1/\unicode[STIX]{x1D708}$ ( $s=0$ ) regime. The second impurity enhances the flux driven by the impurity density gradient, whilst introducing an oppositely directed component to the flux, when both impurity density gradients have the same sign. The net effect of introducing a second collisional species thus depends on the combination $(Z^{-1}A_{1z}-Z_{A}^{-1}A_{1A})$ , producing an additional outward contribution to the flux when this quantity is negative. Note that the result above does not require that the heaviest impurity $z$ is only present in trace quantities, but does also correctly describe that case.