2.1 Drift-kinetic equation

When the guiding-centre coordinates $\boldsymbol{X}$, the parallel velocity $v_{\Vert }$ and the magnetic moment $\unicode[STIX]{x1D707}=m_{a}v_{\bot }^{2}/(2B)$ are chosen as the phase-space variables, where $B=|\boldsymbol{B}|$ and $v_{\bot }$ is the perpendicular velocity, the guiding-centre equations of motion are given by

(2.1)$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{X}}=v_{\Vert }\boldsymbol{b}+\frac{1}{Z_{a}eB_{\Vert }^{\ast }}\boldsymbol{b}\times (m_{a}v_{\Vert }^{2}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B+Z_{a}e\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}), & \displaystyle\end{eqnarray}$$

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{\Vert }=-\frac{1}{m_{a}v_{\Vert }}\dot{\boldsymbol{X}}\boldsymbol{\cdot }(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B+Z_{a}e\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}), & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{b}=\boldsymbol{B}/B$, $\boldsymbol{B}^{\ast }=\unicode[STIX]{x1D735}\times \boldsymbol{A}^{\ast }$ is the corrected magnetic field with the guiding-centre vector potential $\boldsymbol{A}^{\ast }=\boldsymbol{A}+m_{a}v_{\Vert }\boldsymbol{b}/(Z_{a}e)$, $B_{\Vert }^{\ast }=\boldsymbol{B}^{\ast }\boldsymbol{\cdot }\boldsymbol{b}$ and the overdot denotes the total time derivative ${\dot{A}}\equiv \,\text{d}A/\text{d}t$ of any function $A$. It is assumed that $\unicode[STIX]{x2202}\boldsymbol{A}/\unicode[STIX]{x2202}t=\unicode[STIX]{x2202}\boldsymbol{A}^{\ast }/\unicode[STIX]{x2202}t=0$ in deriving the equations above.

Splitting the distribution function of species $a$ as $f_{a}=f_{a0}+f_{a1}$, where $f_{a1}\ll f_{a0}$, FORTEC-3D solves the equation for $f_{a1}$. The lowest-order distribution function $f_{a0}$ is given by

(2.4)$$\begin{eqnarray}f_{a0}=f_{aM}\text{e}^{-Z_{a}e\unicode[STIX]{x1D6F7}_{1}/T_{a}},\end{eqnarray}$$

with the local Maxwellian

(2.5)$$\begin{eqnarray}f_{aM}=n_{a0}\left(\frac{m_{a}}{2\unicode[STIX]{x03C0}T_{a}}\right)^{3/2}\exp \left[-\frac{m_{a}v_{\Vert }^{2}}{2T_{a}}-\frac{\unicode[STIX]{x1D707}B}{T_{a}}\right].\end{eqnarray}$$

The first-order drift-kinetic equation for $f_{a1}$ is then given by

(2.6)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}f_{a1}}{\unicode[STIX]{x2202}t}+\dot{\boldsymbol{X}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{a1}+\dot{v}_{\Vert }\frac{\unicode[STIX]{x2202}f_{a1}}{\unicode[STIX]{x2202}v_{\Vert }}-C_{TP}(f_{a1})=C_{FP}(f_{a0})+\frac{Z_{a}e}{T_{a}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}_{1}}{\unicode[STIX]{x2202}t}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\dot{\boldsymbol{X}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r\left[\frac{n_{a0}^{\prime }}{n_{a0}}+\frac{Z_{a}e\unicode[STIX]{x1D6F7}_{0}^{\prime }}{T_{a}}+\left(\frac{m_{a}v_{\Vert }^{2}}{2T_{a}}+\frac{\unicode[STIX]{x1D707}B}{T_{a}}-\frac{3}{2}+\frac{Z_{a}e\unicode[STIX]{x1D6F7}_{1}}{T_{a}}\right)\frac{T_{a}^{\prime }}{T_{a}}\right]f_{a0},\end{eqnarray}$$

where $n_{a0}$ and $T_{a}$ are flux functions, $m_{a}$ denotes the mass of the species $a$ and $C_{TP}(f_{a1})$ and $C_{FP}(f_{a0})$ are the test-particle part and the field-particle part of the linearised collision operator, respectively,

(2.7)$$\begin{eqnarray}\displaystyle & \displaystyle C_{TP}(f_{a1})=\mathop{\sum }_{b}C_{TP}^{ab}(f_{a1},f_{b0}), & \displaystyle\end{eqnarray}$$

(2.8)$$\begin{eqnarray}\displaystyle & \displaystyle C_{FP}(f_{a0})=\mathop{\sum }_{b}C_{FP}^{ab}(f_{a0},f_{b1}). & \displaystyle\end{eqnarray}$$

Strictly speaking, however, $f_{a0}$ is approximated with $f_{aM}$ in the linearised collision operator for the present work. Also, the recent attempts to extend the collision operator to include higher-order terms which can be important for impurity transport (Sugama et al. Reference Sugama, Matsuoka, Satake, Nunami and Watanabe2019; Calvo et al. Reference Calvo, Parra, Velasco and García-Regaña2020) are not considered. The equations (2.1)–(2.6) are reduced to the equations for a case without $\unicode[STIX]{x1D6F7}_{1}$ by setting $\unicode[STIX]{x1D6F7}_{1}=0$. A derivation of the drift-kinetic equation (2.6) is provided in appendix A.

Note that in this model, the term corresponding to the so-called $E\times B$ drift

(2.9)$$\begin{eqnarray}\boldsymbol{v}_{E}\equiv \frac{1}{B_{\Vert }^{\ast }}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}\end{eqnarray}$$

is species-dependent due to the denominator

(2.10)$$\begin{eqnarray}B_{\Vert }^{\ast }=B+\frac{m_{a}v_{\Vert }}{Z_{a}e}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}\end{eqnarray}$$

and FORTEC-3D does solve (2.1)–(2.3) as written here so that the phase-space volume is conserved (Littlejohn Reference Littlejohn1983). However, the second term in (2.10) is negligible as far as the low-beta approximation

(2.11)$$\begin{eqnarray}\unicode[STIX]{x1D735}\times \boldsymbol{b}\simeq \frac{\boldsymbol{B}\times \unicode[STIX]{x1D735}B}{B^{2}}\end{eqnarray}$$

holds. Then, $\boldsymbol{v}_{E}$ can be approximated as

(2.12)$$\begin{eqnarray}\boldsymbol{v}_{E}\simeq \frac{1}{B}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}\end{eqnarray}$$

and the value of the $E\times B$ drift becomes species-independent. Also, it does not affect the basic property that the lowest-order distribution function does not generate radial particle flux:

(2.13)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{a0}\equiv \left\langle \int \,\text{d}^{3}v(\boldsymbol{v}_{m}+\boldsymbol{v}_{E})\boldsymbol{\cdot }\unicode[STIX]{x1D735}rf_{a0}\right\rangle =0.\end{eqnarray}$$

The verification of (2.13) is provided in appendix A.

2.2 Expressions in the Boozer coordinate system

The Boozer coordinate system (Boozer Reference Boozer1981) is used in FORTEC-3D. In this coordinate system, the covariant representation of the magnetic field is given by

(2.14)$$\begin{eqnarray}\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D704}(r)\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=I(r)\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+G(r)\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}+\unicode[STIX]{x1D6FD}^{\ast }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$

where $\unicode[STIX]{x1D704}$ is the rotational transform, $\unicode[STIX]{x1D713}$ is the toroidal magnetic flux divided by $2\unicode[STIX]{x03C0}$, $I$ is related to the toroidal current inside the flux surface at $r$ and $G$ is related to the poloidal current outside the surface, respectively. Since $\unicode[STIX]{x1D6FD}^{\ast }$ is related to the Pfirsch–Schlüter current (Boozer Reference Boozer1981) and can be estimated as $O(\unicode[STIX]{x1D6FD}/a)$, where $\unicode[STIX]{x1D6FD}$ and $a$ are the plasma beta and the minor radius, respectively, we assume low-$\unicode[STIX]{x1D6FD}$ plasmas and therefore the term $\unicode[STIX]{x1D6FD}^{\ast }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ in (2.14) is ignored in this work. This approximation allows us to put the guiding-centre equations of motion in canonical form. Although we will not write down the equations in canonical form explicitly here, this property has an advantage of ensuring the properties of Hamiltonian systems such as Liouville’s theorem. See § 3.2 of White (Reference White2014) and references therein for a detailed discussion of the construction of the Hamiltonian formulation of the guiding-centre motion. The contravariant representation is given by

(2.15)$$\begin{eqnarray}\boldsymbol{B}=\frac{1}{\sqrt{g}}\frac{\text{d}\unicode[STIX]{x1D713}}{\text{d}r}(\unicode[STIX]{x1D704}\boldsymbol{e}_{\unicode[STIX]{x1D703}}+\boldsymbol{e}_{\unicode[STIX]{x1D701}}),\end{eqnarray}$$

where $\boldsymbol{e}_{j}$ is the unit vector in the $j$-direction and the Jacobian $\sqrt{g}$ is given by

(2.16)$$\begin{eqnarray}\sqrt{g}=\frac{\text{d}\unicode[STIX]{x1D713}}{\text{d}r}\frac{G+\unicode[STIX]{x1D704}I}{B^{2}}=\frac{(\unicode[STIX]{x1D713}^{\prime }G+\unicode[STIX]{x1D712}^{\prime }I)}{B^{2}},\end{eqnarray}$$

with $\unicode[STIX]{x1D712}$ the poloidal magnetic flux divided by $2\unicode[STIX]{x03C0}$.

In this coordinate system, the guiding-centre equations of motion are expressed as follows (a derivation is provided in § A.2):

(2.17)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}r}{\text{d}t}=\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FE}}\left(I\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}-G\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)+\frac{Z_{a}e}{\unicode[STIX]{x1D6FE}}\left(I\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}-G\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right), & \displaystyle\end{eqnarray}$$

(2.18)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}t}=\frac{G}{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}B^{\prime }+Z_{a}e\unicode[STIX]{x1D6F7}^{\prime })+\frac{v_{\Vert }B}{\unicode[STIX]{x1D6FE}}\left(Z_{a}e\unicode[STIX]{x1D712}^{\prime }-\frac{m_{a}v_{\Vert }}{B}G^{\prime }\right), & \displaystyle\end{eqnarray}$$

(2.19)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D701}}{\text{d}t}=-\frac{I}{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}B^{\prime }+Z_{a}e\unicode[STIX]{x1D6F7}^{\prime })+\frac{v_{\Vert }B}{\unicode[STIX]{x1D6FE}}\left(Z_{a}e\unicode[STIX]{x1D713}^{\prime }+\frac{m_{a}v_{\Vert }}{B}I^{\prime }\right), & \displaystyle\end{eqnarray}$$

(2.20)$$\begin{eqnarray}\displaystyle \frac{\text{d}v_{\Vert }}{\text{d}t} & = & \displaystyle -\frac{v_{\Vert }}{B\unicode[STIX]{x1D6FE}}\left[Z_{a}e\left(I\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}-G\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)\unicode[STIX]{x1D6F7}^{\prime }-\unicode[STIX]{x1D707}\left(I\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}-G\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)B^{\prime }\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Z_{a}eB}{m_{a}\unicode[STIX]{x1D6FE}}\left[\left(\frac{m_{a}v_{\Vert }}{Z_{a}eB}G^{\prime }-\unicode[STIX]{x1D712}^{\prime }\right)\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+Z_{a}e\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\left(\frac{m_{a}v_{\Vert }}{Z_{a}eB}I^{\prime }+\unicode[STIX]{x1D713}^{\prime }\right)\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}+Z_{a}e\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\right)\right],\end{eqnarray}$$

where

(2.21)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}\equiv \unicode[STIX]{x1D707}+\frac{m_{a}v_{\Vert }^{2}}{B}\end{eqnarray}$$

and

(2.22)$$\begin{eqnarray}\unicode[STIX]{x1D6FE}\equiv Z_{a}eBB_{\Vert }^{\ast }\sqrt{g}\end{eqnarray}$$

are defined. $\unicode[STIX]{x1D6FF}$ here is not to be confused with the delta function or some ordering parameter. The flux surface average of an arbitrary function $A$ is expressed as

(2.23)$$\begin{eqnarray}\langle A\rangle =\frac{1}{V^{\prime }}\oint \sqrt{g}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}A,\end{eqnarray}$$

where

(2.24)$$\begin{eqnarray}V^{\prime }=\oint \sqrt{g}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}.\end{eqnarray}$$

The radial particle flux is thus calculated by

(2.25)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{a}=\frac{1}{V^{\prime }}\oint \sqrt{g}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}\int \,\text{d}^{3}v{\dot{r}}f_{a1}.\end{eqnarray}$$

2.3 Two-weight $\unicode[STIX]{x1D6FF}\,f$ method

The two-weight $\unicode[STIX]{x1D6FF}\,f$ scheme (Brunner, Valeo & Krommes Reference Brunner, Valeo and Krommes1999; Wang et al. Reference Wang, Nakajima, Okamoto and Murakami1999) is employed in FORTEC-3D. The microscopic dynamics is simulated by marker particles which evolve according to the equations of motion (2.1)–(2.3). In the two-weight $\unicode[STIX]{x1D6FF}\,f$ method, the distribution function is split into two parts as $f=f_{0}+f_{1}$, where $f_{0}$ is an analytically known part and $f_{1}$ is the deviation from that. In this study, the adiabatic response to $\unicode[STIX]{x1D6F7}_{1}$ is included in the $f_{0}$ part (2.4). See García-Regaña et al. (Reference García-Regaña, Beidler, Kleiber, Helander, Mollén, Alonso, Landreman, Maassberg, Smith and Turkin2017) for a discussion of different choices of $f_{0}$. The species index is omitted for simplicity in this subsection.

The distribution functions are represented by the products of the marker particle distribution $g(\boldsymbol{Z},t)$ and the weight fields $W(\boldsymbol{Z},t)$ and $P(\boldsymbol{Z},t)$, respectively, as

(2.26)$$\begin{eqnarray}\displaystyle & f_{0}(\boldsymbol{Z},t)=P(\boldsymbol{Z},t)g(\boldsymbol{Z},t), & \displaystyle\end{eqnarray}$$

(2.27)$$\begin{eqnarray}\displaystyle & f_{1}(\boldsymbol{Z},t)=W(\boldsymbol{Z},t)g(\boldsymbol{Z},t). & \displaystyle\end{eqnarray}$$

Each marker particle is assigned two weights $w_{i}$ and $p_{i}$, which are defined as the values of the weight fields $W$ and $P$ at the phase position $\boldsymbol{Z}_{i}$, respectively:

(2.28)$$\begin{eqnarray}\displaystyle & w_{i}(t)=W(\boldsymbol{Z}_{i}(t),t), & \displaystyle\end{eqnarray}$$

(2.29)$$\begin{eqnarray}\displaystyle & p_{i}(t)=P(\boldsymbol{Z}_{i}(t),t), & \displaystyle\end{eqnarray}$$

where $i$ is the particle label. The weights also evolve in time following the equations

(2.30)$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{p}}_{i}=\frac{p_{i}}{f_{0}}\dot{\boldsymbol{Z}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{Z}}f_{0}-\frac{Ze}{T}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}_{1}}{\unicode[STIX]{x2202}t}p_{i}, & \displaystyle\end{eqnarray}$$

(2.31)$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{w}}_{i}=-\frac{p_{i}}{f_{0}}\left(\dot{\boldsymbol{Z}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{Z}}-C_{\text{FP}}\right)f_{0}+\frac{Ze}{T}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}_{1}}{\unicode[STIX]{x2202}t}p_{i}, & \displaystyle\end{eqnarray}$$

which are evaluated at $\boldsymbol{Z}=\boldsymbol{Z}_{i}(t)$. We thus can solve the drift-kinetic equation for $f_{1}$ by integrating the evolution equations (2.1)–(2.3), (2.30) and (2.31).

2.4 Evaluation of electrostatic potential

The electrostatic potential can be decomposed as $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}_{0}+\unicode[STIX]{x1D6F7}_{1}$. The first part is a flux function $\unicode[STIX]{x1D6F7}_{0}=\unicode[STIX]{x1D6F7}_{0}(r)$ and generates a radial electric field $E_{r}=-\unicode[STIX]{x1D6F7}_{0}^{\prime }$. The radial electric field $E_{r}$ is determined so that it satisfies the ambipolar condition

(2.32)$$\begin{eqnarray}\mathop{\sum }_{a}Z_{a}\unicode[STIX]{x1D6E4}_{a}=0.\end{eqnarray}$$

In FORTEC-3D, $E_{r}=-\unicode[STIX]{x1D6F7}_{0}^{\prime }$ can be given as a fixed input parameter or it also can be internally solved by the equation (2a) in Satake et al. (Reference Satake, Kanno and Sugama2008).

The second part $\unicode[STIX]{x1D6F7}_{1}=\unicode[STIX]{x1D6F7}_{1}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ is the non-flux function part. It is assumed that $\langle \unicode[STIX]{x1D6F7}_{1}\rangle =0$ in the present study. See appendix B for the justification of this assumption. The non-uniform part $\unicode[STIX]{x1D6F7}_{1}$ is determined by the following argument. From (2.4), the density of species $a$ takes the form

(2.33)$$\begin{eqnarray}n_{a}=n_{a0}\exp (-Z_{a}e\unicode[STIX]{x1D6F7}_{1}/T_{a})+n_{a1}\end{eqnarray}$$

and imposing quasi-neutrality up to the first order yields the expression

(2.34)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{1} & = & \displaystyle \frac{1}{e}\left(\mathop{\sum }_{I}Z_{I}^{2}\frac{n_{I0}}{T_{I}}+\frac{n_{e0}}{T_{e}}\right)^{-1}\mathop{\sum }_{I}Z_{I}n_{I1}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{e}\left(\mathop{\sum }_{I}Z_{I}^{2}\frac{n_{I0}}{T_{I}}+\frac{n_{e0}}{T_{e}}\right)^{-1}\mathop{\sum }_{I}Z_{I}\int \,\text{d}^{3}vf_{I1},\end{eqnarray}$$

where the subscript $I$ refers to the ion species. In deriving the above expression, $Z_{a}e\unicode[STIX]{x1D6F7}_{1}/T_{a}\ll 1$ is assumed for all species. As we shall see below, this assumption is reasonable for the present study. However, the assumption does not hold for high-$Z$ impurities. Also, the adiabatic response of electrons is assumed, i.e. $n_{e1}=0$. The spatial structure of $\unicode[STIX]{x1D6F7}_{1}$ is thus determined by the superposition of the ion density variations $n_{I1}$, weighted by the charge number $Z_{I}$. Note here that in some studies, $\unicode[STIX]{x1D6F7}_{1}$ is evaluated by considering only the contribution of the main ion. When the contribution of impurities is taken into account as (2.34), the value of $\unicode[STIX]{x1D6F7}_{1}$ usually becomes smaller since the impurity contributions are included in the factor $\sum _{I}Z_{I}^{2}n_{I0}/T_{I}$ in the denominator. In this study, $\unicode[STIX]{x1D6F7}_{1}$ is evaluated during the self-consistent calculation of the set of equations (2.1)–(2.3), (2.6) and (2.34). It has been pointed out that the contribution of kinetic electrons to $\unicode[STIX]{x1D6F7}_{1}$ can be significant, especially in electron-root plasmas (García-Regaña et al. Reference García-Regaña, Estrada, Calvo, Velasco, Alonso, Carralero, Kleiber, Landreman, Mollén and Sánchez2018). However, it is numerically too costly for the present version of FORTEC-3D to solve kinetic electrons with ion species together.

2.5 Evaluation of the flux driven by $\unicode[STIX]{x1D6F7}_{1}$

In the Boozer coordinate system, the radial particle flux driven by $\unicode[STIX]{x1D6F7}_{1}$ through the $E\times B$ drift, $\boldsymbol{v}_{E1}$, is expressed as

(2.35)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{v_{E1}}(r) & = & \displaystyle \frac{1}{V^{\prime }(r)}\oint \sqrt{g}(\boldsymbol{X})\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}\int \,\text{d}^{3}v\frac{1}{B(\boldsymbol{X})B_{\Vert }^{\ast }(\boldsymbol{X},\boldsymbol{v})\sqrt{g}(\boldsymbol{X})}\nonumber\\ \displaystyle & & \displaystyle \times \,\left(I(r)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}(\boldsymbol{X})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}-G(r)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}(\boldsymbol{X})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)f_{a1}(\boldsymbol{X},\boldsymbol{v}),\end{eqnarray}$$

where $\boldsymbol{X}=(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ and $\boldsymbol{v}=(v_{\Vert },\unicode[STIX]{x1D707})$. This quantity can be evaluated in two different ways. One is to simply perform the integral and the other is to use Fourier series expansion. If the low-beta approximation

(2.36)$$\begin{eqnarray}\frac{1}{B(\boldsymbol{X})B_{\Vert }^{\ast }(\boldsymbol{X},\boldsymbol{v})}\simeq \frac{1}{B^{2}(\boldsymbol{X})}\end{eqnarray}$$

is used, the velocity integral can readily be performed and we obtain

(2.37)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{v_{E1}}(r)=\frac{1}{V^{\prime }(r)}\oint \,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}\bar{v}_{E1}(\boldsymbol{X})n_{a1}(\boldsymbol{X}),\end{eqnarray}$$

where we have defined

(2.38)$$\begin{eqnarray}\bar{v}_{E1}\equiv \frac{1}{B^{2}}\left(I\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}-G\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)\simeq \sqrt{g}\boldsymbol{v}_{E1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r.\end{eqnarray}$$

Thus, if we expand $\bar{v}_{E1}$ and $n_{a1}$ in Fourier series, respectively, and put them into (2.37), only the products of the same modes survive and the equation is formally reduced to the following form:

(2.39)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{v_{E1}}(r)=\frac{1}{V^{\prime }(r)}\mathop{\sum }_{m,n}C_{m,n}\left[\tilde{v}_{E1}^{m,n(c)}(r){\tilde{n}}_{a1}^{m,n(c)}(r)+\tilde{v}_{E1}^{m,n(s)}(r){\tilde{n}}_{a1}^{m,n(s)}(r)\right],\end{eqnarray}$$

where $C_{m,n}$ are normalisation coefficients and $\tilde{v}_{E1}^{m,n(c)}$, $\tilde{v}_{E1}^{m,n(s)}$, ${\tilde{n}}_{a1}^{m,n(c)}$ and ${\tilde{n}}_{a1}^{m,n(s)}$ are the Fourier coefficients of cosine- and sine-$(m,n)$ modes of $\bar{v}_{E1}$ and $n_{a1}$, respectively. In the present article, Fourier coefficients of a quantity $A(\boldsymbol{X})$ are defined as follows:

(2.40)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\tilde{A}^{m,n(c)}(r)={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}^{2}}}\displaystyle \iint \,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}A(\boldsymbol{X})\cos (m\unicode[STIX]{x1D703}-nN\unicode[STIX]{x1D701}),\\ \tilde{A}^{m,n(s)}(r)={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}^{2}}}\displaystyle \iint \,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D701}A(\boldsymbol{X})\sin (m\unicode[STIX]{x1D703}-nN\unicode[STIX]{x1D701}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $N$ is the toroidal period number ($N=10$ for LHD). Note that the flux driven by the magnetic drift, $\boldsymbol{v}_{m}$, cannot be put in the form of (2.37) and therefore also not in the form of (2.39) since $\boldsymbol{v}_{m}$ essentially depends on the velocity variables (and the dependence cannot be approximated away).

Although the first method is more simple and straightforward for evaluating $\unicode[STIX]{x1D6E4}_{v_{E1}}$, this second method has an advantage that it enables us to see how each mode contributes to the additional flux. The $\text{C}^{6+}$ flux will be evaluated by these two different methods below.

3.1 Ambipolar radial electric field

The ambipolar condition has been determined by two different (sets of) numerical codes. One is a set of local codes DKES/PENTA and the other is the global code FORTEC-3D. Both FORTEC-3D and DKES/PENTA determined the ambipolar condition by assuming a pure H plasma. Both solutions are represented in figure 5. The red line corresponds to the global solution and the solid blue line to the local solution. As can be seen, those solutions are completely different. The sign of the global solution is positive. In other words, the global calculation indicates that the plasma is in electron-root, not in ion-root. It is not necessarily surprising that the global calculation yields different solutions from local ones. As mentioned in the introduction, it is known that the magnetic drift (even only the tangential component of it) can significantly change the dependence of neoclassical transport on $E_{r}$ (Matsuoka et al. Reference Matsuoka, Satake, Kanno and Sugama2015; Huang et al. Reference Huang, Satake, Kanno, Sugama and Matsuoka2017; Velasco et al. Reference Velasco, Calvo, Parra and García-Regana2020). In Velasco et al. (Reference Velasco, Calvo, Satake, Alonso, Nunami, Yokoyama, Sato, Estrada, Fontdecaba and Liniers2017), the ambipolar conditions for such impurity hole plasmas are investigated based on the collisionality and it is argued that, even if an impurity hole plasma is in ion-root, the plasma profiles are very close to the condition for $E_{r}$ to become positive. In fact, as represented by the dashed line in figure 5, a positive solution is found by the local codes as well. The solution is stable and physically possible. Its existence suggests that the gap between the local and global solutions is not as large as it appears. Yet, we adopt the ion-root (represented by the solid blue line) as our ‘local solution’ since similar solutions are used in previous local studies. We thus investigate the difference of $\unicode[STIX]{x1D6F7}_{1}$ structure and its impact on impurity transport between the cases, each with the local and the global solutions, respectively.

However, whether the sign of the ambipolar-$E_{r}$ is positive or negative is the most crucial issue in the impurity hole study. If the sign of ambipolar-$E_{r}$ turns out to be positive in real impurity hole plasmas, the scenario is completely changed and the impurity hole phenomenon could be explained by standard neoclassical theory arguments. Nunami et al. (Reference Nunami, Satake, Fujita, Yamaguchi, Matsuoka, Sato, Toda and Sugama2019, Reference Nunami, Nakata, Toda and Sugama2020) have intensively investigated the condition for both solutions by considering the contributions of turbulence and external torque. The result of the study also seems to support the possibility of the positive solution.

There is measurement data on the ambipolar-$E_{r}$ profile in an impurity hole plasma and it indicates that the sign of the ambipolar-$E_{r}$ near the magnetic axis is indeed negative (Ido et al. Reference Ido, Shimizu, Nishiura, Nagaoka, Yokoyama, Ida, Yoshinuma, Toi, Itoh and Nakano2010). However, it also shows that the sign of $E_{r}$ transits to positive at some outer radius. Impurity hole has been observed in plasmas with several different $n$–$T$ profiles (Nakamura et al. Reference Nakamura, Tamura, Yoshinuma, Suzuki, Yoshimura, Kobayashi, Yokoyama, Nunami, Nakata and Nagaoka2017), but the measurement data on ambipolar-$E_{r}$ has been reported only for the case mentioned above. The $n$–$T$ profile used in this study is close to but not the same as that in the measured case. Thus, our global solution for the ambipolar condition does not necessarily contradict the existing experimental data.

Since the following concerns remain, however, we cannot simply accept the electron-root as the correct solution yet. First, the positive solution is, as mentioned above, obtained by a single ion-species calculation and the contribution of impurities to the ambipolar condition is not considered. Second, the $\unicode[STIX]{x1D6F7}_{1}$ effect is not included in the ambipolar condition as well. In contrast, negative solutions similar to our local solution have been obtained by self-consistent calculations for multi-ion-species plasmas including $\unicode[STIX]{x1D6F7}_{1}$ (Mollén et al. Reference Mollén, Landreman, Smith, García-Regaña and Nunami2018). How $\unicode[STIX]{x1D6F7}_{1}$ affects the radial particle fluxes in our electron-root case will be examined below.

In any case, the primary purpose of this study is not to examine the ambipolar condition, but to apply the newly extended code to impurity transport in a real helical configuration. Thus, we hold both possibilities here and investigate two cases corresponding to each solution, respectively. Both cases have the same background parameters (such as the magnetic field configuration, densities and temperatures) except for the radial electric field.

4.1 Ion-root case

In figure 6, the radial profiles of particle fluxes are compared between the cases with and without $\unicode[STIX]{x1D6F7}_{1}$. The red lines represent the fluxes for the case without $\unicode[STIX]{x1D6F7}_{1}$ and the green lines for the case with $\unicode[STIX]{x1D6F7}_{1}$. Before discussing the impact of $\unicode[STIX]{x1D6F7}_{1}$, it is worth commenting on the result without $\unicode[STIX]{x1D6F7}_{1}$. The $\text{H}^{1+}$ flux obtained by our global calculation is smaller by one order than that obtained by local calculations for similar plasmas (Mollén et al. Reference Mollén, Landreman, Smith, García-Regaña and Nunami2018). On the other hand, the values of the $\text{He}^{2+}$ and $\text{C}^{6+}$ fluxes in our global result are of the same order as those in the local study. In this ion-root case, the plasma profile belongs in a parameter region where the discrepancy in the $\text{H}^{1+}$ fluxes between global and local (neglecting the entire $\boldsymbol{v}_{m}$) calculations becomes the largest (Matsuoka et al. Reference Matsuoka, Satake, Kanno and Sugama2015; Huang et al. Reference Huang, Satake, Kanno, Sugama and Matsuoka2017). This significant difference in the main ion flux directly affects the ambipolar condition. Thus, it is natural that the solutions of the ambipolar condition turn out to be qualitatively different between the global and the local simulations and an ion-root is not found by the global simulation.

Figure 6.

Radial particle fluxes of $\text{H}^{1+}$ (a), $\text{He}^{2+}$ (b) and $\text{C}^{6+}$ (c) for the ion-root case without (red) and with (green) $\unicode[STIX]{x1D6F7}_{1}$, respectively.

Figure 7.

From left, structures of $\unicode[STIX]{x1D6F7}_{1}$, $n_{H.1}$, $n_{He,1}$ and $n_{C,1}$ on the flux surface at $r/a=0.2$ for the ion-root case. Figures (a–d) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (e–h) the results with $\unicode[STIX]{x1D6F7}_{1}$. Note that the scale of the colour contour is different for each figure. $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ are poloidal and toroidal angle coordinates in the Boozer coordinate system, respectively, and $N=10$.

Now, let us examine the $\unicode[STIX]{x1D6F7}_{1}$ effect on the particle fluxes. In Mollén et al. (Reference Mollén, Landreman, Smith, García-Regaña and Nunami2018), it is shown that the $\text{C}^{6+}$ flux in ion-root plasmas is inwardly enhanced by $\unicode[STIX]{x1D6F7}_{1}$ while the $\text{H}^{1+}$ and the $\text{He}^{2+}$ fluxes are outwardly driven. A similar result has also been obtained for the $\text{C}^{6+}$ flux in Velasco et al. (Reference Velasco, Calvo, García-Regaña, Parra, Satake and Alonso2018). Figure 6 shows that our global result for the $\text{C}^{6+}$ flux is qualitatively analogous to those local results. That is, the $\text{C}^{6+}$ flux is driven further inwardly, not outwardly. This means that $\unicode[STIX]{x1D6F7}_{1}$ contributes to the impurity accumulation, rather than mitigate it. On the other hand, positive enhancement is seen in the $\text{H}^{1+}$ and $\text{He}^{2+}$ fluxes, especially at outer radii $r/a>0.8$.

Let us look into the structures of $\unicode[STIX]{x1D6F7}_{1}$ and density variation of each ion species to see which results in the modification of the particle fluxes in figure 6. In figure 7, the phase structures of the normalised $\unicode[STIX]{x1D6F7}_{1}$ and the density variation of each ion species on the flux surface at $r/a=0.2$ are shown. Figure 7(a–d) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and figure 7(e–h) the result with $\unicode[STIX]{x1D6F7}_{1}$. Here ‘$\unicode[STIX]{x1D6F7}_{1}$ for the case without $\unicode[STIX]{x1D6F7}_{1}$’ means $\unicode[STIX]{x1D6F7}_{1}$ constructed from the solutions of the drift-kinetic equation neglecting $\unicode[STIX]{x1D6F7}_{1}$ through (2.34). When the difference in $\unicode[STIX]{x1D6F7}_{1}$ between the cases with and without $\unicode[STIX]{x1D6F7}_{1}$ in this sense becomes significant, the necessity of the self-consistent calculation is indicated. Figures 8 and 9 make the same comparison at $r/a=0.5$ and $r/a=0.8$, respectively. The amplitude of $\unicode[STIX]{x1D6F7}_{1}$ is smaller than earlier local results overall and the factor $Z_{a}e\unicode[STIX]{x1D6F7}_{1}/T_{a}$ is at most the order of $10^{-1}$ even for $\text{C}^{6+}$ in this case (for example, $e\unicode[STIX]{x1D6F7}_{1}/T_{i}=1.0\times 10^{-2}$ at $r/a=0.5$ roughly amounts to 30 V of $\unicode[STIX]{x1D6F7}_{1}$ in this case). The comparison of the density structures above implies that $\unicode[STIX]{x1D6F7}_{1}$ is too small to strongly modify the spectrum compositions of the ion species over the flux surface.

Figure 8.

From left, structures of $\unicode[STIX]{x1D6F7}_{1}$, $n_{H.1}$, $n_{He,1}$ and $n_{C,1}$ on the flux surface at $r/a=0.5$ for the ion-root case. Figures (a–d) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (e–h) the results with $\unicode[STIX]{x1D6F7}_{1}$. Note that the scale of the colour contour is different for each figure. $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ are poloidal and toroidal angle coordinates in the Boozer coordinate system, respectively, and $N=10$.

Figure 9.

From left, structures of $\unicode[STIX]{x1D6F7}_{1}$, $n_{H.1}$, $n_{He,1}$ and $n_{C,1}$ on the flux surface at $r/a=0.8$ for the ion-root case. Figures (a–d) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (e–h) the results with $\unicode[STIX]{x1D6F7}_{1}$. Note that the scale of the colour contour is different for each figure. $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ are poloidal and toroidal angle coordinates in the Boozer coordinate system, respectively, and $N=10$.

This can also be seen in figure 10, which represents the Fourier spectra of density variations of the ion species for the cases with and without $\unicode[STIX]{x1D6F7}_{1}$. The density variations are expanded in real trigonometric series. Only leading modes are plotted and reddish colours are assigned to cosine modes and bluish colours to sine modes. Although the Fourier spectra are evaluated on 50 flux surfaces in FORTEC-3D, only the values at $r/a=0.0,0.2,\ldots ,1.0$ are plotted for visualisation purpose. The mode compositions are not largely affected by the inclusion of $\unicode[STIX]{x1D6F7}_{1}$. One point worth noting is that the phase structure of $\text{C}^{6+}$ is qualitatively different from those of the other two species. $\text{H}^{1+}$ or $\text{He}^{2+}$ mainly consists of stellarator symmetric components, which are preserved under the transformation $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})\rightarrow (-\unicode[STIX]{x1D703},-\unicode[STIX]{x1D701})$, as expected for the ions in the $\sqrt{\unicode[STIX]{x1D708}}$ regime (Calvo et al. Reference Calvo, Velasco, Parra, Alonso and García-Regaña2018b). On the other hand, the main components in the density spectrum of $\text{C}^{6+}$, which is in the plateau regime, change their signs under the same transformation.

Figure 10.

Fourier spectra of $n_{H,1}$ (a,d), $n_{He,1}$ (b,e) and $n_{C,1}$ (c,f) for the ion-root case. Figures (a–c) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (d–f) the results with $\unicode[STIX]{x1D6F7}_{1}$.

Figure 11.

Fourier spectra of $\bar{v}_{E1}$ (a) and $\unicode[STIX]{x1D6F7}_{1}$ (b) for the ion-root case.

The spectrum of $\bar{v}_{E1}$ is plotted in figure 11 as well as the spectrum of $\unicode[STIX]{x1D6F7}_{1}$ itself. For $\text{H}^{1+}$ and $\text{He}^{2+}$, following (2.39), most of the cosine modes in their density spectra couple with the same modes in the $\bar{v}_{E1}$ spectrum having the same signs. These effective couplings result in the outward enhancement of the $\text{He}^{2+}$ flux and $\text{H}^{1+}$ flux at $r/a>0.8$. On the other hand, the sine-$(1,0)$ mode dominates in the $\text{C}^{6+}$ spectrum. The value of this mode peaks around $r/a=0.6$ and the mode couples with $\bar{v}_{E1}$ having the opposite sign and results in the inward enhancement of the particle flux.

In figure 12(a), the radial $\text{C}^{6+}$ flux driven by $v_{E1}$ is plotted as well as the total $\text{C}^{6+}$ flux for the case with $\unicode[STIX]{x1D6F7}_{1}$. The red line represents the total $\text{C}^{6+}$ flux calculated by the integral (1.1). The green line represents the $v_{E1}$-driven part of the $\text{C}^{6+}$ flux calculated by the integral (2.35). The blue line is also the $v_{E1}$-driven part of the $\text{C}^{6+}$ flux, but it is obtained by the sum of the Fourier coefficient products in the form of (2.39). The impact of $\unicode[STIX]{x1D6F7}_{1}$ on the $v_{m}$-driven flux is also shown in figure 12(b). The green line represents the $\text{C}^{6+}$ flux driven by the magnetic drift for the case without $\unicode[STIX]{x1D6F7}_{1}$, i.e. the total $\text{C}^{6+}$ for the case without $\unicode[STIX]{x1D6F7}_{1}$. The blue line represents the $\text{C}^{6+}$ flux driven by $v_{m}$ for the case with $\unicode[STIX]{x1D6F7}_{1}$. It can been seen that the $v_{m}$-driven flux is also enhanced by the inclusion of $\unicode[STIX]{x1D6F7}_{1}$. This is because, as shown in figure 10, $\unicode[STIX]{x1D6F7}_{1}$ amplifies the sine-$(1,0)$ mode in the $\text{C}^{6+}$ spectrum and hence strengthens the opposite-sign couplings between the magnetic drift and the $\text{C}^{6+}$ distribution function as well. Here the sine-$(1,0)$ component of $\boldsymbol{v}_{m}$ represents the grad-$B$ drift by the toroidicity of the magnetic field. Therefore, at least in this case, the presence of $\unicode[STIX]{x1D6F7}_{1}$ is unfavourable for impurity accumulation avoidance.

Figure 12.

(a) The total radial $\text{C}^{6+}$ flux, $\unicode[STIX]{x1D6E4}_{C}$, for the case with $\unicode[STIX]{x1D6F7}_{1}$ (red), the $v_{E1}$-driven part of $\unicode[STIX]{x1D6E4}_{C}$ calculated by the integral (green) and the $v_{E1}$-driven part calculated by the products of the Fourier coefficients (blue). (b) The total radial $\text{C}^{6+}$ flux, $\unicode[STIX]{x1D6E4}_{C}$, for the case with $\unicode[STIX]{x1D6F7}_{1}$ (red), the total $\text{C}^{6+}$ flux (driven only by $\boldsymbol{v}_{m}$) for the case without $\unicode[STIX]{x1D6F7}_{1}$ (green) and the $v_{m}$-driven part for the case with $\unicode[STIX]{x1D6F7}_{1}$ (blue).

Figure 13.

Radial particle fluxes of $\text{H}^{1+}$ (a), $\text{He}^{2+}$ (b) and $\text{C}^{6+}$ (c) for the electron-root case without (red) and with (green) $\unicode[STIX]{x1D6F7}_{1}$, respectively.

Figure 14.

From left, structures of $\unicode[STIX]{x1D6F7}_{1}$, $n_{H,1}$, $n_{He,1}$ and $n_{C,1}$ on the flux surface at$r/a=0.2$ for the electron-root case. Figures (a–d) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (e–h) the results with $\unicode[STIX]{x1D6F7}_{1}$. Note that the scale of the colour contour is different for each figure. $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ are poloidal and toroidal angle coordinates in the Boozer coordinate system, respectively, and $N=10$.

Figure 15.

From left, structures of $\unicode[STIX]{x1D6F7}_{1}$, $n_{H,1}$, $n_{He,1}$ and $n_{C,1}$ on the flux surface at $r/a=0.8$ for the electron-root case. Figures (a–d) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (e–h) the results with $\unicode[STIX]{x1D6F7}_{1}$. Note that the scale of the colour contour is different for each figure. $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ are poloidal and toroidal angle coordinates in the Boozer coordinate system, respectively, and $N=10$.

4.2 Electron-root case

In figure 13, the radial particle fluxes with and without $\unicode[STIX]{x1D6F7}_{1}$ for the electron-root case are compared. Regardless of the species, the radial particle flux is outwardly enhanced by $\unicode[STIX]{x1D6F7}_{1}$. The enhancement is largest for $\text{C}^{6+}$ and the existence of an outward convection is indicated. However, the enhancements of the $\text{H}^{1+}$ and $\text{He}^{2+}$ fluxes are not negligible as well. This result implies that the impact of $\unicode[STIX]{x1D6F7}_{1}$ is too large to be ignored in the calculation of the ambipolar condition. Thus, it is expected that an electron-root slightly different from that used in this work would be found by calculation including $\unicode[STIX]{x1D6F7}_{1}$.

The reason for the strong outward enhancement can be seen in the phase structures of density variations. The spatial structures of $\unicode[STIX]{x1D6F7}_{1}$ and $n_{I1}$ at $r/a=0.2$ and $r/a=0.8$ for this case are illustrated in figures 14 and 15, respectively. The amplitude of $n_{I1}$ and hence of $\unicode[STIX]{x1D6F7}_{1}$ are at the same level as those in the ion-root case. Also, the phase structures of $\unicode[STIX]{x1D6F7}_{1}$, $\text{H}^{1+}$ and $\text{He}^{2+}$ seemingly look similar between the ion-root and the electron-root cases. However, the $\text{C}^{6+}$ density structure is distinctively different between the two cases. By comparing figures 9(d,h) and 15(d,h), it can be seen that the $\text{C}^{6+}$ density phase is inverted in the $\unicode[STIX]{x1D703}$-direction.

This can also be seen in terms of Fourier spectra of the density variations and $\bar{v}_{E}$ shown in figures 16 and 17, respectively. As shown in figure 16, the signs of the sine-$(1,0)$ modes in the density spectra are inverted from those in the ion-root case. This is true not only for $\text{C}^{6+}$, but also for $\text{H}^{1+}$ and $\text{He}^{2+}$. Further, this mode is amplified by the inclusion of $\unicode[STIX]{x1D6F7}_{1}$. Figure 17 shows that the sine-$(1,0)$ mode in the $\hat{v}_{E}$ spectrum, which comes from the cosine-$(1,0)$ mode in the $\unicode[STIX]{x1D6F7}_{1}$ spectrum, is negative and is a major component for this case as well. Thus, the inversion of the sine-$(1,0)$ modes in the density spectra leads to strong same-sign couplings and contributes to the outward radial particle fluxes. This sign inversion is not a consequence of the inclusion of $\unicode[STIX]{x1D6F7}_{1}$ but, rather, it seems to be a consequence of the presence of the positive radial electric field.

Figure 16.

Fourier spectra of $n_{H,1}$ (a,d), $n_{He,1}$ (b,e) and $n_{C,1}$ (c,f) for the electron-root case. Figures (a–c) represent the results without $\unicode[STIX]{x1D6F7}_{1}$ and (d–f) the results with $\unicode[STIX]{x1D6F7}_{1}$.

Figure 17.

Fourier spectra of $\bar{v}_{E1}$ (a) and $\unicode[STIX]{x1D6F7}_{1}$ (b) for the electron-root case with $\unicode[STIX]{x1D6F7}_{1}$.

In figure 18, the contributions of the $v_{E1}$-driven part (a) and the $v_{m}$-driven part (b) to the total radial $\text{C}^{6+}$ flux are separately represented. The red lines in the figures represent the total $\text{C}^{6+}$ flux for the case with $\unicode[STIX]{x1D6F7}_{1}$. The green and blue lines in figure 18(a) represent the $v_{E1}$-driven part of the $\text{C}^{6+}$ flux, each calculated by the integral and the sum of Fourier coefficients, respectively. In figure 18(b), the green line represents the total $\text{C}^{6+}$ flux for the case without $\unicode[STIX]{x1D6F7}_{1}$, which is driven only by $\boldsymbol{v}_{m}$. The blue line is the $v_{m}$-driven flux for the case with $\unicode[STIX]{x1D6F7}_{1}$. As for the ion-root case, the $v_{m}$-driven flux is reinforced by $\unicode[STIX]{x1D6F7}_{1}$ as well. The primary cause of this enhancement is the amplification of the sine-$(1,0)$ mode in the $\text{C}^{6+}$ spectrum by $\unicode[STIX]{x1D6F7}_{1}$ (see figure 16c,f). As the result, the $v_{m}$-driven part becomes at most an order of magnitude larger than the $v_{E1}$-driven part.

Figure 18.

(a) The total radial $\text{C}^{6+}$ flux, $\unicode[STIX]{x1D6E4}_{C}$, for the case with $\unicode[STIX]{x1D6F7}_{1}$ (red), the $v_{E1}$-driven part of $\unicode[STIX]{x1D6E4}_{C}$ calculated by the integral (green) and the $v_{E1}$-driven part calculated by the products of the Fourier coefficients (blue). (b) The total radial $\text{C}^{6+}$ flux, $\unicode[STIX]{x1D6E4}_{C}$, for the case with $\unicode[STIX]{x1D6F7}_{1}$ (red), the total $\text{C}^{6+}$ flux (driven only by $v_{m}$) for the case without $\unicode[STIX]{x1D6F7}_{1}$ (green) and the $v_{m}$-driven part for the case with $\unicode[STIX]{x1D6F7}_{1}$ (blue).

4.3 Discussion of the result

For the ion-root case, a result analogous to previous local results has been obtained: the $\text{C}^{6+}$ flux is driven further inwardly, not outwardly by $\unicode[STIX]{x1D6F7}_{1}$. This enhancement is mainly caused by the couplings of the sine-$(1,0)$ modes between the $\text{C}^{6+}$ density spectrum and the radial drift velocities. $\unicode[STIX]{x1D6F7}_{1}$ amplifies this mode and the opposite-sign couplings of the sine-$(1,0)$ modes result in the further inward flux. For the electron-root case, on the other hand, the radial particle flux is outwardly enhanced regardless of the species mainly due to the same-sign couplings of the sine-$(1,0)$ modes.

The phase structure of $\unicode[STIX]{x1D6F7}_{1}$ is mainly determined by the spatial distribution of main ions. The main components in the $\unicode[STIX]{x1D6F7}_{1}$ spectrum are then determined by the magnetic field configuration. The spectra of $\text{H}^{1+}$ or light impurities such as $\text{He}^{2+}$ in low-collisionality regimes mainly consist of stellarator symmetric components. Thus, stellarator anti-symmetric modes become dominant in the spectrum of $\bar{v}_{E1}$. On the other hand, stellarator anti-symmetric components were found to be dominant in the $\text{C}^{6+}$ density spectrum, and these components tend to couple well with both the magnetic drift and the $E\times B$ drift. The important question here is ‘with which sign do the leading modes, in particular the sine $(1,0)$-modes, of the impurity ions couple with the radial drift velocities?’. An induction from our results implies that the sign of $E_{r}$ is a key factor to determine the signs of the sine-$(1,0)$ modes and, for ion-root plasmas, opposite-sign couplings are realised. To avoid those effective negative couplings or turn them into positive ones, either the phase structure of $\unicode[STIX]{x1D6F7}_{1}$ or the impurity density has to be strongly modified. However, we found that $\unicode[STIX]{x1D6F7}_{1}$ only increases the absolute value of the sine-$(1,0)$ mode in the $\text{C}^{6+}$ density spectrum and does not seem to affect its sign. Increasing the amplitude of $\unicode[STIX]{x1D6F7}_{1}$ is unlikely to make a qualitative difference either. Thus, our conclusion is that it is difficult for $\unicode[STIX]{x1D6F7}_{1}$ to mitigate the inward neoclassical $\text{C}^{6+}$ flux, let alone drive it in the outward direction in ion-root plasmas unless a strong particle or torque source exists.

The argument above ignores the contribution of kinetic electrons and it is possible that their contribution modifies the $\unicode[STIX]{x1D6F7}_{1}$ structure (García-Regaña et al. Reference García-Regaña, Estrada, Calvo, Velasco, Alonso, Carralero, Kleiber, Landreman, Mollén and Sánchez2018). Although it is worth examining, however, it is doubtful that the impact of kinetic electrons can be strong enough to reverse the conclusion. On the other hand, our result does not simply support the electron-root scenario either. In the present study, we have demonstrated that $\unicode[STIX]{x1D6F7}_{1}$ is indeed a necessary factor to study the role of neoclassical transport in the formation of impurity hole, even in a global simulation model. Although it is an important finding that the global simulation alters the ambipolar-$E_{r}$ profile from the prediction by local codes, more detailed and accurate calculation, in which the ambipolar-$E_{r}$ profile is determined without neglecting the effect of $\unicode[STIX]{x1D6F7}_{1}$ on the radial transport, remains for future research.