Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments

The description of parallel impurity transport can be developed starting from the zero and first order moments of the drift-kinetic equation, that is, the continuity and momentum balance equations [28, 29],

$$\begin{align} \nabla \cdot (n_Z \mathbf{v}_Z) = 0, \end{align} \tag{ 1 }$$

$$\begin{align} & n_Z m_Z \left(\mathbf{v}_Z \cdot \nabla \right) \mathbf{v}_Z + \nabla p_Z + \nabla \cdot \pi_Z \nonumber\\ & \quad + Z e n_Z \left(\nabla \phi + \mathbf{v}_Z \times \mathbf{B}\right) = \mathbf{R}_Z. \end{align} \tag{ 2 }$$

Here, n_Z and $p_Z = n_Z T_\mathrm i$ are the impurity density and pressure, with $T_\mathrm i$ the ion temperature, $\mathbf{v}_Z$ the fluid velocity, π_Z the viscous tensor, and φ the electrostatic potential. The friction force $\mathbf{R}_Z$ includes the friction of the impurity with all the other species. The magnetic field is given by the usual general expression $\mathbf{B} = I \nabla \varphi + \nabla \varphi \times \nabla \psi$, where ψ is the poloidal flux divided by 2π, ϕ is the toroidal angle, and the function I(ψ) is related to the poloidal current inside a magnetic surface. The perpendicular projection of the momentum balance, keeping only the dominant terms provided by pressure gradient and electric force, in combination with the divergence free condition, provides the usual expression of the first order velocity, which relies upon the assumption that the flow velocity is mostly within the flux surface (or, more precisely, that the divergence of the radial flux is much smaller than that of the flux on the magnetic surface)

$$\begin{align} \mathbf{v}_Z = \omega_Z(\psi) R {\hat{\mathbf{e}}}_\phi + \frac{K_Z(\psi)}{n_Z} \mathbf{B}, \end{align} \tag{ 3 }$$

where

$$\begin{align} \omega_Z(\psi) = -\left( \frac{\partial \Phi}{\partial \psi} + \frac{1}{Z e n_Z} \frac{\partial p_Z}{\partial \psi} \right). \end{align} \tag{ 4 }$$

More generally, ω_Z can also include the contributions from the inertial terms from equation (2).

The parallel projection of the momentum balance gives

$$\begin{align} & \frac{m_Z n_Z}{2} \left[ K_Z^2(\psi)\, \nabla_{\parallel} \left(\frac{B^2} {n_Z^2}\right) - {\omega_Z^2}\nabla_{\parallel} R^2 \right] \nonumber\\ & \quad + \nabla_{\parallel}p_Z + \nabla_{\parallel} \cdot \pi_Z +Z e n_Z \nabla_\parallel \Phi = R_{Z\parallel} \end{align} \tag{ 5 }$$

where from left to right we identify the inertia terms connected with poloidal and toroidal centrifugal forces, the pressure gradient term, where the temperature can be considered homogeneous on magnetic surfaces, the viscous term and finally the electric and the friction parallel forces. This is a complete expression, which describes all the collisionality regimes and includes both poloidal and toroidal inertial (centrifugal) terms as well as friction. Starting from this expression, depending on the collisional regime, the adopted ordering and the physical approach, different terms have been included or neglected in various models. In particular, theoretical approaches have been developed with different assumptions on the collisionality regime of the main ions (banana, banana–plateau or Pfirsch–Schlüter), on the collisionality regime and the charge of the impurities (thereby neglecting the viscosity tensor term for impurities in the high collisionality limit), considering impurities only in the trace limit ($\alpha = n_Z Z^2/n_\mathrm i \ll 1$), neglecting inertial terms and thereby centrifugal forces or neglecting the friction force at the lowest order. A comprehensive and more detailed recent review of the different theoretical approaches to parallel impurity transport can be found in [56], section 2.2 and, correspondingly, comparisons of the theoretical predictions of the poloidal variation of a high-Z impurity density have been recently performed on molybdenum measurements in C-Mod [57]. More recent works which are not included in that review are the following [35, 36, 58]. In [58] the role of the Coriolis force is also considered, and a model is developed in which the combined effects of the centrifugal force, the Coriolis force, the poloidal electric field and the ion-impurity friction force on the in–out impurity density asymmetry are investigated. Maget et al [35, 36] extend the work of [32] with the inclusion of the diamagnetic velocity of the impurities in the calculation of the impurity–main ions parallel friction, removing the assumption $Z \gg 1$ and also allowing for the possibility of poloidal inhomogeneities of the main ion density.

The calculation of the friction force plays a critical role, since, as it will be also presented in the next subsection, it is directly connected with both the classical and neoclassical components of the collisional radial flux. In fact, from the momentum balance equation, equation (2), one can see that the perpendicular (radial) particle flux is produced by both perpendicular and parallel friction [29], where the perpendicular friction causes classical collisional transport and parallel friction the neoclassical collisional transport. One of the first examples of the realization of its importance and its direct calculation in the Pfirsch–Schlüter regime can be found in [59]. The friction force between the impurity with charge Z and the (main) ion species 'i' is defined as the first order velocity moment of the collision operator,

$$\begin{align} \mathbf{R}_{Zi} = \int{d^3v \; m_Z \, \mathbf{v} \, C_{Zi}(f_Z)} = - \int{d^3v \; m_i \,\mathbf{v} \, C_{iZ}(f_i)}, \end{align} \tag{ 6 }$$

where the second equality expresses the property of momentum conservation of the Coulomb collisions. Of course, knowledge of the main ion distribution function, for instance by means of a numerical solver, allows the calculation of the friction force by integration in the velocity space. A very common approach considers an expansion of the perturbed distribution function in generalized Laguerre polynomials [28], where the coefficients of the polynomials are directly connected to macroscopic moments of the distribution function, with the first two providing the parallel particle and heat flows. The corresponding expansion of the friction force allows one to express this force as a linear combination of these moments of the distribution function, where the coefficients are the friction coefficients [28]. Analytical expressions for the parallel friction have been derived in various limits, mainly connected with different collisionality regimes of the main ions, in the banana regime [31, 32], in the Pfirsch–Schlüter regime [33] and more recently the plateau regime [34]. The parallel particle and heat flows can be expressed as linear combinations of radial density and temperature gradients as well as the corresponding poloidal flow functions K_Z (ψ). Assuming isothermal flux surfaces, the parallel momentum balance equation becomes an equation for the poloidal distribution of the impurity density,

$$\begin{align} \frac{\nabla_{\parallel} n_Z} {n_Z} & = \frac{m_Z}{2T_\mathrm i} \left[ {\omega_Z^2}\nabla_{\parallel} R^2 - K_Z^2(\psi)\, \nabla_{\parallel} \left(\frac{B^2} {n_Z^2}\right)\right] \nonumber \\ & \quad - \frac{Z e}{T_\mathrm i} \nabla_\parallel \Phi - \frac{m_Z \nu_{zi}}{T_\mathrm i} \left(v_{Z\parallel} - v_{\mathrm i\parallel} + C_Z \frac{2 q_{\mathrm i\parallel}}{5 p_\mathrm i}\right), \end{align} \tag{ 7 }$$

where on the right hand side we can again identify the toroidal and poloidal centrifugal terms, the parallel electric force, and the parallel friction term. In equation (7) we have neglected the viscous tensor term, which becomes small for highly charged impurities, typically in a high collisionality regime [31]. According to the general approach that we have just described, the parallel friction has been expressed as

$$\begin{align} R_{Z\parallel} = - n_Z m_Z \nu_{zi} \left( v_{Z\parallel} - v_{\mathrm i\parallel} + C_Z \frac{2 q_{\mathrm i\parallel}}{5 p_\mathrm i} \right), \end{align} \tag{ 8 }$$

with C_Z an appropriate coefficient and $q_{\mathrm i\parallel}$ the parallel ion heat flow, which can be expressed as proportional to the ion temperature radial gradient. In general, in the limit of $Z \gg 1$, the diamagnetic term can be neglected in the expression of the parallel impurity flow, and the expression for the parallel friction becomes independent of the radial gradients of the impurity. We notice that this does not rigorously correspond to a trace impurity limit $\alpha \ll 1$, but the trace limit still has the similarity of removing the dependence of the parallel friction on the impurity density. If in contrast the dependence of the parallel friction on the impurity density and on its radial gradient is kept through the inclusion of the diamagnetic velocity

$$\begin{align} v_{Z\parallel} = \frac{K_Z(\psi)}{n_Z}B - \frac{I}{B} \left( \frac{\partial \Phi}{\partial \psi} + \frac{1}{Z e n_Z} \frac{\partial p_Z}{\partial \psi} \right), \end{align} \tag{ 9 }$$

then the equation for the parallel distribution of the impurity density also depends on the radial gradients of the impurity density profile through the friction term [35, 36]. In this more complete treatment the poloidal asymmetry of the impurity depends on the radial profile shape. Concomitantly, the radial transport depends on the poloidal asymmetry of the impurity, with the consequence that parallel and perpendicular dynamics are completely connected, and the solution leads to a 'natural' impurity poloidal asymmetry [35, 36] which is consistent with the radial profiles solution of the radial transport equation. In general, the parallel friction force can be expressed through a linear combination of radial gradients of the main ion density and ion temperature, through its dependence on the parallel particle and heat flows. The coefficients of this linear combination, including their poloidal dependence, vary with the collisionality regime. However, the impact of the friction term can be characterized by considering the following dimensionless parameter, first introduced in [31], which is largely equivalent, and more precisely quantifies, the impact of the ordering parameter Δ,

$$\begin{align} g = \frac{q}{\epsilon} \, \frac{m_\mathrm i n_\mathrm i}{e L_\perp (\tau_{iZ}\,n_Z) \langle B \cdot \nabla \theta \rangle}. \end{align} \tag{ 10 }$$

Here  = r/R is the local inverse aspect ratio, and τ_iZ is the ion–impurity collision time

$$\begin{align} \tau_{iZ} = \frac{3(2 \pi)^{1.5} \epsilon_0^2 m_i^{0.5} T_i^{1.5}}{n_Z Z^2 e^4 \textrm{ln} \Lambda}, \end{align} \tag{ 11 }$$

while the corresponding ion–impurity collision frequency is given by $\nu_{iZ} = 3 \sqrt{\pi} v_\mathrm {thi}^3 /(4 \tau_{iz} v^3)$. The characteristic perpendicular length $L_\perp$ is provided by a linear combination of ion density and temperature logarithmic gradients ($R/L_\perp = R/L_{ni}-0.5\,R/L_{Ti}$ and $R/L_\perp = R/L_{ni}$ with main ions in the banana [32] and the Pfirsch–Schlüter [33] regimes respectively). It is easy to observe that this parameter can be rewritten in the equivalent form

$$\begin{align} g = \frac{\rho_{\theta i}}{L_\perp} \frac{qR}{v_\mathrm {thi}} \frac{Z^2}{\tau_{ii}} \simeq \Delta. \end{align} \tag{ 12 }$$

We also observe that friction becomes important when radial logarithmic gradients are large enough in combination with a sufficiently high ion collisionality, and in particular that it is proportional to the impurity charge squared. Therefore it is much stronger for highly charged impurities, whereas it is in general negligible for light impurities, even in the presence of large gradients in the pedestal region. It can also be expected that, at the low collisionalities of a reactor burning plasma, effects coming from friction are very small and only become important close to the separatrix in the peripheral part of the pedestal. In contrast, they can be significant over large portions of the minor radius in present devices for an impurity like W in high density L–mode plasmas, and certainly in the pedestal of present H–mode plasmas, even for light impurities.

The centrifugal term connected with the poloidal rotation is usually negligible and it is regularly neglected in almost the totality of the theoretical models. It has been recently considered in the analysis of ASDEX Upgrade light impurity (boron and nitrogen) density profiles which were concomitantly measured at the low and high field side [60], and has been found to become non–negligible very close to the separatrix, in the external region of the pedestal, while friction plays the dominant role over the entire pedestal.

The centrifugal term connected with the toroidal rotation is usually the dominant term for heavy impurities in the core plasma in the presence of toroidal rotation. It always has to be considered in consistent combination with the electrostatic potential term. The first complete treatment of this term computing the impact on the impurity density and its transport is in [38], starting from the rigorous approach developed in [37]. It has been followed by several other works [32, 36, 39–41]. The essential feature here is that while the centrifugal force pushes impurities on the outboard side, the electrostatic force connected with the different distribution of electrons and ions produced by the same rotation develops a parallel electric field which pushes the impurities on the high field side and slightly opposes the dominant effect of the centrifugal force. Other effects can contribute to the development of a parallel electric field. The most noticeable among these is connected with the presence of an ion species with a temperature anisotropy, as generated by ICRH [42–46] and also by NBI, although usually with weaker effects due to a lower level of anisotropy, but larger density of fast particles [45, 48].

The inhomogeneity of the poloidal distribution of the impurity density strongly affects the neoclassical transport. This can be easily understood, since the net radial flux averaged over the flux surface results from almost equal contributions of opposite sign from the outboard and the inboard sides. Thereby, any modification of the size of the inward or outward components, as produced by an increase of the particle density on the low or high field side, can directly lead to much larger radial fluxes [31–33, 35, 36, 38, 40, 41].

At present there are models and corresponding codes which allow these effects to be computed consistently, although with different levels of assumptions. No complete model is available which takes into account all of the possible sources of poloidal inhomogeneity and consistently computes the consequent radial transport while solving the drift-kinetic equation with a sufficiently complete expression of the collision operator. The drift–kinetic solver NEO [53–55] has been used in many recent applications for impurity transport calculation requiring the inclusion of centrifugal [49–52] and temperature anisotropy [50, 52] effects. This code is based on the Hinton ordering [37] and thereby includes centrifugal effects at the lowest order, but does not include friction effects at the lowest order, thereby it cannot reproduce effects which are instead included in the Helander approach [31], which allows friction effects to be considered at the lowest order. For the comparisons with experimental measurements, complete validation tests could be performed by simultaneous measurements of all of the components that concomitantly impact the parallel transport and determine the poloidal distribution of the impurities. Simultaneous measurements of the toroidal and poloidal flows and of the poloidal distribution of the impurities would provide a more stringent and complete test for the theoretical predictions [57].

As anticipated, and as well known [28, 29], classical and neoclassical perpendicular transport are connected with the perpendicular and the parallel components of the friction force respectively.

The usual approach is to consider the toroidal projection of the momentum balance equation, equation (2). It is of practical use to report here analytical expressions for the classical and neoclassical transport. From equation (2), and considering a rigid toroidal rotation only, with angular velocity Ω, we have

$$\begin{align} & -\frac{m_Z n_Z \Omega^2}{2} \nabla R^2 + \nabla p_Z +\nabla \cdot{\pi_Z} \nonumber\\ & \quad - n_Z Ze \left(-\nabla \phi+ \mathbf{v}_Z \times \mathbf{B}\right) = \mathbf{R_Z}, \end{align} \tag{ 13 }$$

where we have also kept the contribution of the viscous tensor, responsible for the banana-plateau transport.

By taking the scalar product with R²∇ϕ, the toroidal component of this equation allows us to express the radial fluxes in the form

$$\begin{align} \langle \Gamma_Z \cdot \nabla \psi \rangle = \frac{1}{Ze} \left[ I \left \langle \frac{\mathbf{R}_Z\cdot \mathbf{B}}{B^2} \right \rangle + \left \langle \mathbf{B} \cdot \frac{\nabla{\psi} \times {\mathbf{R}_Z}}{B^2} \right \rangle \right]. \end{align} \tag{ 14 }$$

The second term, proportional to the perpendicular component of the friction, is the classical transport. The first term, proportional to the parallel component of the friction is the neoclassical transport. In the conventional situation of a poloidally homogeneous impurity density, if the viscous term is neglected, then, since $\langle \mathbf{B} \cdot \nabla \; \rangle = 0$, one has that $\langle \mathbf{R}_Z \cdot \mathbf{B} \rangle = 0$. In the case of poloidally inhomogeneous density, this solubility condition becomes $\langle \mathbf{R}_Z \cdot \mathbf{B} / n_Z \rangle = 0$ [31, 32], with the consequence that the neoclassical transport can be split into the two contributions

$$\begin{align} \langle \Gamma_Z^{\textrm{neo}} \cdot \nabla \psi \rangle & = \frac{I}{Ze} \left[ \left \langle \frac{\mathbf{R}_Z\cdot \mathbf{B}} {n_Z} \right \rangle \; \left \langle \frac{B^2}{n_Z} \right \rangle^{-1} \right.\nonumber\\ & \quad \left. + \left \langle \frac{\mathbf{R}_{Z} \cdot \mathbf{B}}{n_Z} \left( \frac{n_z}{B^2} - \left \langle \frac{B^2}{n} \right \rangle^{-1} \right) \right \rangle \right]. \end{align} \tag{ 15 }$$

Here, the first term on the right hand side is the banana-plateau component, which is substantial only when the viscous term π_Z is important, while the second is the Pfirsch–Schlüter component, which is present independently of the strength of the viscous term. In the limit of highly charged impurities, $Z \gg 1$, which thereby are in a high collisionality regime, the viscous term becomes small compared to the other terms [31], and the dominant contribution is provided by the Pfirsch–Schlüter component.

An analytical expression of the radial flux of Pfirsch–Schlüter impurity transport can then be obtained without making any assumption on the specific poloidal inhomogeneity of the impurity density distribution. Assuming main ions are in the banana limit and with a poloidally homogeneous main ion density, and neglecting the impurity diamagnetic velocity in the friction term, consistent with the limit $Z \gg 1$, this reads [32, 40, 50, 51],

$$\begin{align} \frac{R \langle \Gamma_Z^\mathrm {neo} \rangle}{\langle n_Z \rangle} = \frac{q^2 D_c} {Z} \left[ \left( \frac{1}{Z} \frac{R}{L_{nZ}} \,{-}\, \frac{R}{L_{ni}} {+} \frac{1}{2} \frac{R}{L_{Ti}} \right) P_A \; {-}\,c_0\, P_B \, f_c \, \frac{R}{L_{Ti}} \right]. \end{align} \tag{ 16 }$$

Here the numerical coefficient c₀ is equal to 0.33 in the large aspect ratio limit with an ion–ion collision operator for the main ion distribution function given by pitch-angle scattering combined to a momentum restoring term [32]. The two geometrical coefficients P_A and P_B encapsulate the dependence on the poloidal distribution of the impurity density,

$$\begin{align*} P_A & = \frac{1}{2\epsilon^2} \,\frac{\langle B^2 \rangle}{\langle n_Z \rangle} \left[ \left \langle \frac{n_Z}{B^2} \right \rangle \, \,-\, \left \langle \frac{B^2}{n_Z} \right \rangle ^{-1} \right], \\ P_B & = \frac{1}{2\epsilon^2} \, \frac{\langle B^2 \rangle}{\langle n_Z \rangle} \left[ \frac{\langle n_Z \rangle}{\langle B^2 \rangle} \,-\, \left \langle \frac{B^2}{n_Z} \right \rangle ^{-1} \right]. \end{align*}$$

If the impurity density is poloidally homogeneous, in the large aspect ratio limit, these coefficients are P_A  = 1 and P_B  = 0. In present devices with large NBI heating P_A and P_B can be as large as 100 [61]. The result presented here closely follows the derivation in [32]. Analogous derivations, also in the limit $Z\gg1$, have been obtained with main ions in other collisionality regimes, in particular banana-plateau [34] and Pfirsch–Schlüter [33]. The treatment presented in [32] has been recently generalized to consistently include the diamagnetic velocity of the impurities in the friction term and the poloidal inhomogeneity of the main ions with sources of poloidal asymmetry caused by friction only [35], as well as produced by toroidal rotation or temperature anisotropies of ICRH resonant species [36]. A complete analytical treatment within the drift-kinetic approach valid for general geometry and arbitrary impurity charge has been also obtained in the limit in which both the main ions and the impurities are in the Pfirsch–Schlüter regime [41]. The combined impact of poloidal asymmetries and pressure anisotropies of the impurities on the neoclassical impurity transport has been recently presented in [62] with an analytical derivation, also in relation with the effects that can be produced by turbulence in directly affecting the collisional transport. The aspect of the direct interaction between turbulent and collisional transport will be reviewed in section 5. The analytical model presented in [40] and described by equation (16) has been compared with the predictions of the code NEO [53–55] in [63] for applications to JET plasmas. The NEO code contains a complete treatment of the linearized collision operator, valid in all collisionality regimes, but only assumes in input a poloidal inhomogeneity produced by the centrifugal force from toroidal rotation [38] or a specific form of temperature anisotropy of an ICRH resonant species, following [45], as applied e.g. in [50]. The analytical formula has been found to reproduce reasonably well the NEO numerical results, even when impurities move out of the Pfirsch–Schlüter limit, but requires ions in the banana limit. The strength of the temperature screening terms critically depends on the assumptions made for the calculation of the main ion parallel flow, which is based on an approximated collision operator and can become inappropriate in general geometry. The more general treatment presented in [35, 36], which directly includes friction in the self–consistent determination of the parallel distribution of the main ion and impurity densities, following [32] but without dropping the contribution from the impurity diamagnetic velocity in the friction force, requires main ions in the banana limit. It could be extended to other collisional regimes profiting of the results in [33, 34].

The possibility of developing a sufficiently accurate and fast model for routine integrated modelling applications should be considered a high priority in the direction of obtaining a complete pulse simulator which can be applied for the modelling of an entire plasma discharge with limited computational time. Fluid codes like NCLASS [64] and NEOART [65] based on the Hirshman–Sigmar approach [28] can be used for the calculation of the neoclassical transport of light impurities, but are limited to the most general assumptions of conventional neoclassical ordering and do not include any effect related to the poloidal asymmetry of heavy highly charged impurities. Therefore, in general these codes cannot be applied for realistic computations of the transport of heavy impurities in the presence of any poloidal asymmetry (while of course they are highly reliable for the calculation of heat conductivities of the main species and of the bootstrap current, as well as of the transport of impurities which do not feature any poloidal asymmetry). A set of practical analytical formulae, again in the limit of conventional neoclassical theory ordering, for both the Pfirsch–Schlüter and the banana-plateau contributions can be found here [66], where a useful summary (and the correction of some typos) of the formulae obtained in [28] is provided.

We recall that the entire treatment presented in this section always assumes that $\delta = \rho_{\theta} / L_{\perp} \ll 1$, that is, still follows the local neoclassical ordering. However, following [31] and related works [32–34] as well as more recently [35, 36], it can allow for a non-negligible friction force at the lowest order, that is $g \approx \Delta = (\rho_\theta /L_{\perp}) (\nu_{ii}Z^2 qR)/v_\mathrm {th} \approx 1$. The code NEO [53–55] follows the conventional neoclassical ordering which assumes both $\delta \ll 1$ and $\Delta \ll 1$, however, as already mentioned, it allows for poloidal inhomogeneities of the density distribution as produced by rotation [37, 38] and temperature anisotropies of ICRH resonant species. Considering that friction effects are usually small to negligible for high temperature (low collisionality) plasmas in high confinement regimes, the code NEO is a very practical and reliable tool for accurate applications dedicated to the calculation of the impact of toroidal rotation on heavy and highly charged impurities like W. It however requires considerable computational time for integrated modelling applications in a transport code.

Finally, as we already mentioned and it is well known, local neoclassical theory orders small the normalized poloidal Larmor radius δ. This ordering can become not appropriate to describe the plasma close to the magnetic axis due to the large orbits and in the edge pedestal region, also due to the strong radial gradients. This leads to the need of developing models and corresponding codes which are not based on the local ordering. Codes with radially global properties, like XGCO [67] and XGCa [68], KEATS [69], PERFECT [70], and IMPGYRO [71] can be applied to these purposes, particularly for the plasma periphery, however they do not include rotational effects. Also global gyrokinetic codes which are usually applied for turbulence simulations but which also have sufficiently sophisticated multi–species collision operators [72–74] can be applied for global neoclassical calculations. Examples of these codes that have recently also addressed the problem of neoclassical transport are GYSELA [75] and GT5D [76]. These codes can also address the problem of the concomitant calculation of neoclassical and turbulent transport, and the direct interaction between the two transport mechanisms. These effects will be reviewed in section 4. Recent studies on global effects [77, 78] have demonstrated a collisional convection which is provided by the combination of an inward pinch connected with the main ion density gradient and a temperature screening term, consistent with the local theory. The strength of the temperature screening has been recently subject of investigations, and it is found in agreement with local theory when global simulations also reproduce the parallel ion heat flux [62].

A general expression of the turbulent impurity particle flux can be formally derived from the gyrokinetic equation, where we also consider the presence of the radial gradient of the toroidal rotation, while still considering the limit in which the impurity toroidal rotation is small compared to its thermal velocity (that is, we neglect centrifugal effects). As already mentioned, a more complete treatment which also includes the impact of the centrifugal effects can be found in [106], where the formulation of the gyrokinetic equation in the plasma co-rotating frame [111] has been adopted.

Here below we derive gyrokinetic expressions of the particle and heat fluxes of a species σ, where we single out the dependence on the particle charge Z. We consider the linearized gyrokinetic equation for the perturbed distribution function

$$\begin{align*} g_{Z} = \delta f_{Z} \,-\; Z e \, \frac{F_{M\,Z}} {T_\mathrm i} \, {\tilde \phi} \end{align*}$$

of a species with charge Z and mass A, and a Maxwellian equilibrium distribution function F_{M  Z} with temperature $T_\mathrm i$. Here ${\tilde \phi}$ is the perturbed electrostatic potential. We consider a single toroidal mode number n = k_y r/q in Fourier space, with r the local minor radius, q the local safety factor and k_y the wave number binormal to the magnetic field, in simple s − α geometry. For simplicity, in the analytic model we consider the electrostatic collisionless limit, neglect centrifugal effects, and focus on the elements of the gyrokinetic equation which are essential to identify the main properties of the diffusive and convective terms. The gyrokinetic equation for the amplitude of the Fourier mode $\hat{g}_{Z k}$, with $g_{Z k} = \hat{g}_{Z k} \; \exp(-i\omega t \,+\, i\textbf{k} \cdot \textbf{x})$, reads

$$\begin{align} & (\omega - \omega_{Gk}) \hat{g}_{Z k} = -\frac{Z T_e}{T_\mathrm i} \nonumber \\ & \quad \times \Bigg\{\omega_{DZk} \Bigg[ \frac{R}{L_{nZ}} + \left( \frac{E}{T_\mathrm i}-\frac{3}{2} \right) \frac{R}{L_{Ti}} \nonumber \\ & \qquad \quad+ \frac{B_T}{B} \frac{v_\parallel}{v_{thZ}}u_Z^{\prime} \Bigg] + \omega \Bigg\} F_{MZ} J_0(k_\perp \rho_Z) \hat{\phi}_k, \end{align} \tag{ 17 }$$

where the generalized gyro–center frequency contains both parallel and drift gyro-center motion

$$\begin{align*} \omega_{G Z k} \; = \; v_\parallel k_\parallel + \omega_{d Z k} \end{align*}$$

and where the parallel wave number $k_\parallel$ is considered to be a number (and not an operator) in this formal derivation. Here $E = m_Z v_\parallel^2/2 + \mu B$ is the energy of the gyro-center, with $\mu = m_Z V_\perp^2/2/B$ the magnetic moment, the fluctuating electrostatic potential is normalized as $\hat{\phi}_k = e \tilde{\phi}_k/T_e$ and the normalized toroidal rotation gradient is $u_Z^{\prime} = -R^2 d \Omega / d r / v_{th Z}$, with Ω the plasma toroidal angular velocity. For simplicity, collisions have been neglected in this analytical description since they do not alter the decomposition of the turbulent particle flux that follows. The gyro–center motion is given by

$$\begin{align} \frac{d\mathbf{X}}{dt} & = v_\parallel \textbf{b} + \frac{m_Z v_\parallel ^2 \textbf{b} \times (\textbf{b} \cdot \nabla) \textbf{b} +\mu \textbf{b} \times \nabla B}{Z e B} \nonumber \\ & \quad +2 \frac{m_Z v_\parallel}{Z e B} \boldsymbol{\Omega}_\perp + \frac{\textbf{b} \times \nabla (\langle \phi \rangle)}{B}, \end{align} \tag{ 18 }$$

where in the order we can recognize the parallel motion, the curvature and grad B drifts, the Coriolis drift due to the plasma angular velocity $\boldsymbol{\Omega}_\perp$, and the E × B drift (non–linear term). Second order terms in the background plasma rotation arising from the centrifugal drift have been omitted (again, general derivation in [106]). In this geometry the expression of the drift frequency is

$$\begin{align*} \omega_{d Z k} & = \textbf{v}_{d Z} \cdot \textbf{k} = - \omega_{D Z k} \, \frac{v_\parallel^2 \,+\, v_\perp^2/2}{v^2_{th\,Z}} \; \nonumber\\ & \quad \times [ \cos \theta + (s \theta - \alpha \sin \theta) \sin \theta], \end{align*}$$

where we have introduced a generic (fluid) drift frequency defined as $\omega_{D Z k} = {k_y T_\mathrm i}/{(Z e B R)}$ and where here $\alpha = -q^2 R d\beta / dr$. By formally solving the gyrokinetic equation and adopting usual definitions of fluxes produced by the fluctuating E × B drift, we find formal analytical expressions of the impurity particle flux

$$\begin{align} \Gamma_{Z k} = \frac{k_y c_s^2}{\Omega_{ci}} \left \langle{\int{d^3v F_{M\,Z} }} {{\frac{\hat{\gamma}_k[R/L_{n Z} +(E/T_\mathrm {i} - 3/2)R/L_{T i} \;+\; ({m_Z v_{thZ} v_\parallel}/{T_\mathrm i}) u_Z^{\prime} - \hat{\omega}_{dZ}]}{(\hat{\omega}_{rk} - \hat{\omega}_{G Z k})^2 + \hat{\gamma}_k^2}\;J_0(k_\perp \rho_s)^2 |\hat{\phi}_k|^2 } }\right \rangle. \end{align} \tag{ 19 }$$

As shown in [109], the impurity transport produced by magnetic flutter is usually negligible, due to the small thermal velocity of heavy impurities, and it will not be specifically discussed here. In contrast, electromagnetic effects can affect the E × B transport, in particular reducing the impurity turbulent diffusion in ion temperature gradient (ITG) turbulence [109, 110].

From the above expression we obtain the usual decomposition of the impurity particle flux,

$$\begin{align} \frac{R \Gamma_Z} {n_Z} = D_{nZ} \frac{R}{L_{n Z}} + D_{T Z} \frac{R}{L_{T Z} } + D_{u Z} u_Z^{\prime} + V_{p Z} \end{align} \tag{ 20 }$$

where the coefficients D_nZ , D_TZ , D_uZ and V_pZ represent impurity turbulent diffusion, thermo-diffusion, roto-diffusion and pure convection respectively. We recall that equation (20) expresses a linear relation in the case of trace impurities, also in nonlinear simulations, however it requires a local description to be rigorously satisfied [112]. It is practical to also define normalized coefficients defined as $C_{TZ} = D_{T Z} / D_{n Z}$, $C_{uZ} = D_{u Z} / D_{n Z}$, and $C_{p Z} = R V_{pZ} / D_{n Z}$ for thermo–diffusion, roto-diffusion and pure convection respectively. These can be more directly computed in quasi-linear models in particular. The different contributions to the flux will be reviewed in the following subsections.

In order to investigate the nature of the turbulent impurity diffusion, it is instructive to compare its kinetic expression with the corresponding expression of the heat conductivity [113]

$$\begin{align} D_{Z\,k} \propto \int {\frac{\hat{E}^{1/2}\, e^{-\hat{E}}}{(\hat{\omega}_{rk} - \hat{\omega}_{di}(\hat{E})/Z)^2 +\hat{\gamma_{k}}^2} \,d\hat{E}} \end{align} \tag{ 21 }$$

and

$$\begin{align} & \chi_{\sigma \, k} \propto \int \nonumber\\ & \ \times {\frac{\hat{E}^{3/2} e^{-\hat{E}} \left[L_{T\sigma}/L_{n\sigma} + \left(\hat{E}-3/2\right) - \left(\hat{\gamma}_k\hat{\omega}_{d\sigma})/(R/L_{T\sigma}\right) \right]}{(\hat{\omega}_{rk} - \hat{\omega}_{d\sigma}(\hat{E}))^2 +\hat{\gamma_{k}}^2} \,d\hat{E} }, \end{align} \tag{ 22 }$$

where the species index σ refers to both electrons e and main ions i. In equations (21) and (22) we have neglected parallel ion motion and we have introduced the energy only dependent drift frequency $\omega_{d\sigma}(\hat{E}) = - \omega_{D\sigma k} \hat{E}$, where $\hat{E} = E/T_\sigma$. Equation (22) provides the expression of the heat conductivity for both main species electrons and ions, with $\hat{\omega}_{de} = -\hat{\omega}_{di} T_e/ T_\mathrm i$. Following [113], we observe that the condition which maximizes the heat conduction of the main species (electrons and main ions) does not correspond to the condition that maximizes the impurity diffusion. Indeed not only the denominators are different ($Z_\sigma = -1$ and 1 for the main species), but also the conduction and diffusion correspond to different energy moments, and thereby the different powers of the energy in the numerators weigh in a different way the impact of the energy dependent resonance provided by the denominator. The role of the relevant resonance produced by the type of turbulence, and thereby the value of the real frequency of the unstable mode resulting from the solution of the dispersion relation is of fundamental importance in turbulent transport. This aspect has been also pointed out in connection with the different size of ion and electron particle diffusion and convection in different turbulence regimes, leading in particular to much larger values of the diffusion and the convection of the ions than those of the electrons in ITG turbulence [114, 115]. Analogously, when considering the ion heat conductivity in the presence of ITG turbulence, we observe that this is maximized in the presence of a large and positive (ion drift direction) real frequency of the unstable mode, as it results from the solution of the dispersion relation for parameters delivering a strong ITG mode. In contrast, at the same ion charge (Z = 1), the particle diffusivity is maximized by smaller values of the real frequency, due to the weaker energy dependence of the numerator. If we consider a highly charged impurity (Z $\gg$ 1), then the drift frequency at the numerator is practically removed by the large value of Z, and the resonant condition is simply obtained when the mode frequency is close to zero. An analogous argument can be developed for trapped electron mode turbulence, where a large and positive mode frequency produces a large heat transport, but again a small impurity particle transport. The maximum reachable value of $D_Z/\chi_\textrm{eff}$ increases with increasing charge at low charge values of the impurity ($Z \leqslant 10$) and already reaches the $Z \gg 1$ asymptotic value for $Z \geqslant 15$. In conclusion, it is evident that conditions for which the impurity diffusion coefficient is large correspond to conditions where the mode real frequency is neither strongly positive or strongly negative, but in fact where the mode frequency approaches zero. These conditions are obtained when the ion and electron heat fluxes are of comparable magnitude, what suggests that the ratio of the impurity diffusion to the main plasma effective heat conductivity $\chi_\textrm{eff}$ is a non–monotonic function of the electron to ion heat flux ratio. This result is indeed confirmed by a set of non–linear gyrokinetic simulations [113] which find that the impurity diffusivity to heat conductivity ratio can vary by more than an order of magnitude and is maximum when $Q_\mathrm e/Q_\mathrm i \simeq 2$. This theoretical result can provide at least qualitative explanation of early observations of a strong increase of the impurity diffusivity in response to the addition of central electron cyclotron resonance heating (ECRH) [116]. It has also found at least qualitative confirmation in a very recent experimental study [117], although the experimentally estimated values of the impurity diffusivity featured variations which are even larger than those which are theoretically predicted. In particular, significantly lower values of the minimum of $D_Z/\chi_\textrm{eff}$ than those predicted have been obtained, whereas maximum values were found to be more consistent with the theoretical predictions and as large as $D_Z/\chi_\textrm{eff} \simeq 10$. A clear demonstration of a drop of the $D_Z / \chi_\textrm{eff}$ ratio at values of $Q_\mathrm e/Q_\mathrm i \gg 1$ has not been obtained so far, but would provide an additional important piece in the theory validation.

Given the common nature of the properties of both the denominator and the numerator in the expressions of the impurity diffusion and convection coefficients and the main species heat conductivities, the property of being maximized when the turbulent electron and ion heat fluxes are of comparable size is also shared by the convective terms [118]. We start now by examining the thermo-diffusion term. The directions of the convective components are given considering usual monotonically increasing safety factor profiles. If the magnetic shear is negative enough or the α parameter large enough to change the sign of the curvature and grad B drift, all of the terms connected with the compression of this drift reverse their direction.

The turbulent thermo-diffusion is an interesting convective term which is analogous to the corresponding term produced by the logarithmic temperature gradient for electron particle transport [79]. The dominant component of this term stems from the toroidal resonance [90, 94], which, analogously to the electron thermo-diffusive convective part, is characterized by having its radial direction connected to the direction of propagation of the turbulence. A detailed physical explanation at the basis of this thermo-diffusion mechanism can be found in [19] and will not be repeated here. As kinetic integral it is given by

$$\begin{align} D_{T Z k} & = \frac{k_y c_s^2}{\Omega_{ci}} \nonumber\\ & \quad \times \left \langle{\int{d^3v F_{M\, Z} \frac{\hat{\gamma}_k (E/T_{i} - 3/2)}{(\hat{\omega}_{rk} - \hat{\omega}_{G Z k})^2 + \hat{\gamma}_k^2}\;J_0(k_\perp \rho_s)^2 |\hat{\phi}_k|^2 } }\right \rangle, \end{align} \tag{ 23 }$$

from which one can directly realize that in the absence on an energy dependent denominator it would be zero, because $\int{d^3v F_{M\,Z} (E/T_{i} - 3/2)} = 0$. Thereby, the existence of this convective term is directly connected with the energy dependence present at the denominator, which is included in the energy dependence of the gyro-center drift. We consider the curvature and grad B drift resonance first. Because the curvature and grad B drift is proportional to 1/Z and thereby reverses direction from electrons to ions, the direction of the turbulent thermodiffusion is also opposite for electrons and ions. Another connected property is that it reverses direction with the type of turbulence, if turbulence is propagating in the ion or the electron drift direction, ITG and TEM respectively for instance. This is produced by the sign reversal of the kinetic dependence $(E/T-3/2)$, which can lead to an overall positive or negative net thermo-diffusive flux depending on whether the resonant real frequency enhances the contribution at low (inward) or high (outward) energies [19]. In the presence of ITG turbulence, the thermo-diffusion of ions and impurities is directed outward, whereas it is directed inward for the electrons. Thereby, in ITG turbulence thermodiffusion leads to a flattening of the impurity density profile with increasing strength of the ion temperature gradient. With increasing impurity charge, the energy dependence of the denominator is removed (or, in other words, the resonant conditions moves at energies which are very far from the $E/T = 3/2$), and thereby the kinetic integral increasingly approaches the unperturbed situation of $\int{d^3v F_{M\,\sigma} (E/T_{\sigma} - 3/2)} = 0$. Indeed, with increasing impurity charge the thermo-diffusion decreases as 1/Z and, compared to the other convective terms, becomes practically negligible for $Z \geqslant 15$, whereas it is independent of the impurity mass. These dependencies can be also obtained with a fluid description [82, 90, 91]. Moving now to the slab resonance provided by the parallel motion term $k_\parallel v_\parallel$ at the denominator, this is charge independent and does not imply a change of sign along the energy axis and thereby provides a contribution which has the same inward direction for both ions and electrons and in both electron and ion propagating turbulence [79]. It is inversely proportional to the impurity mass.

The roto-diffusion term is caused by the presence of a radial gradient of the toroidal rotation. Its existence has been first pointed out within a fluid treatment and confirmed with gyrokinetic simulations [97]. The analysis of its kinetic integral reveals that it is directly connected with the same mechanisms leading to parallel momentum transport

$$\begin{align} D_{u Z k} = \frac{k_y c_s^2}{\Omega_{ci}} \left \langle{\int{d^3v F_{M\,Z} \frac{\hat{\gamma}_k \; ({m_Z v_{thZ} v_\parallel}/{T_i})}{(\hat{\omega}_{rk} {-} \hat{\omega}_{G Z k})^2 {+} \hat{\gamma}_k^2}\;J_0(k_\perp \rho_s)^2 |\hat{\phi}_k|^2 } }\right \rangle. \end{align} \tag{ 24 }$$

It is equal to zero in the absence of a background toroidal rotation (Coriolis drift) or of another symmetry breaking mechanism leading to a finite $k_\parallel$ (which can be also directly provided by a finite rotation gradient $u^{\prime}$). In contrast, when these terms are present, they produce a direct coupling between parallel velocity fluctuations and density fluctuations, leading to radial particle transport. The mechanisms producing roto-diffusion have complete analogy to those which lead to Coriolis momentum pinch [98] and to the intrinsic rotation. In general, the dominant component comes from the presence of a background toroidal rotation in combination with the toroidal (curvature and grad B) resonance, and it is proportional to A/Z. Thereby, it can remain substantial also for heavy and highly charged impurities, differently from the thermo-diffusion term. Its direction also depends on the direction of propagation of the turbulence, and it is directed outward in the case of the ITG turbulence, in the same direction as the thermo-diffusion term. Plasmas characterized by strong NBI heating and developing strong rotation and ion temperature gradients produce thermo-diffusive and roto-diffusive transport contributions which are both directed outward and which concomitantly lead to a flattening of the impurity profile. This predicted dependence has received some qualitative and in some cases even quantitative experimental validation [119–122]. In the limit of the slab resonance only, roto-diffusion is directed outward and is charge and mass independent. It is also modified by neoclassical modifications of the equilibrium distribution function [123], as it will be more specifically discussed in section 4. Finally, also the problem of the impurity transport produced by the parallel velocity gradient instability in slab geometry with weak magnetic shear has been considered with an analytical quasi-linear model [124], where the properties of both diffusion and convection have been separately studied.

We call pure convection a contribution to the particle flux which is not proportional to any gradient of any kinetic profile. It can be identified in the expression for the particle flux equation (19),

$$\begin{align} RV_{p\,Z k} & = \frac{k_y c_s^2}{\Omega_{ci}} \nonumber\\ & \quad \times\left \langle{\int{d^3v F_{M\, Z} \frac{-\hat{\gamma}_k\hat{\omega}_{d Z}}{(\hat{\omega}_{rk} - \hat{\omega}_{G Z k})^2 + \hat{\gamma}_k^2}\;J_0(k_\perp \rho_s)^2 |\hat{\phi}_k|^2 } }\right \rangle. \end{align} \tag{ 25 }$$

The main mechanism of turbulent convection is connected with the inhomogeneity of the magnetic field, leading to the compressibility of the E × B flow. Thereby it is charge and mass independent, and gives a contribution which is directed inward and in simplified geometry in the strong ballooning limit simply equal to C_p  =−2 (where the −2 stems from the sum of the curvature and grad B drift components, each one providing a contribution of −1/R, and we recall that the dimensionless coefficient C_p is given by the ratio of the pure convection to the diffusivity $C_p = RV_p / D_n$, equation (20)). This mechanism is completely analogous to that present for the main plasma [87, 96] and in the past was also interpreted in terms of turbulent equipartition [125]. This E × B compression term is already present in the impurity transport formulation presented in [82], in combination with the thermodiffusive term, and has been also identified in gyrokinetic simulations of He transport [89]. As already anticipated, with reversed magnetic shear this impurity convection reverses sign from inward to outward [126]. An additional contribution to the convection comes from the slab resonance and the parallel dynamics [90, 94]. This component is interestingly reversing direction depending on the direction of propagation of the turbulence, being directed inward in the case of ITG turbulence, and outward in the case of TEM turbulence. It is proportional to the charge to mass ratio Z/A of the impurity and thereby can remain finite also for highly ionized heavy impurities. Analogously to the other terms connected with the slab resonance, it is proportional to the square of the parallel wave number and therefore in general it also increases with decreasing value of the local safety factor. An indication of the impact of this effect on the peaking of the impurity density profile in the central region has been obtained in specific gyrokinetic modelling of ASDEX Upgrade plasmas with central electron heating [127]. Its role has also been confirmed in nonlinear gyrokinetic simulations of He transport, where the He pure convection has been found to move from inward to outward with increasing electron to ion heat flux ratio [112].

Finally, we recall that also in the case of turbulent transport, modifications of the impurity fluxes are produced in the presence of poloidal density asymmetries of the impurities, although these are less important with respect to neoclassical transport. With increasing impurity toroidal rotation, centrifugal effects mainly provide a large positive (outward) addition to the thermodiffusion which is offset by a corresponding large but negative (inward) addition to the roto-diffusion and a small outward component to the pure convection [106]. While the final result can depend on the actual size of these different components, depending on the local parameters, it usually provides a small inward additional net contribution to the convection, leading to a slightly increased logarithmic gradient of the impurity density profile at the low field side. All of these effects can be largely described by means of a simplified model in which the fluxes are computed taking into account these density asymmetries, while the asymmetries are neglected in the actual solution of the dispersion relation [106]. This simplified approach has been more recently included in the QuaLiKiz transport model, where the turbulent transport of impurities has been generalized to the situation of poloidally asymmetric conditions also including the impact of temperature anistropies [128] (specifically in appendix C). The problem of the impact of the presence of a species with temperature anisotropy, such as a resonant minority species with ICRH, on the turbulent impurity transport has been examined in [102–104]. This effect has been also included in the gyrokinetic code GKW [50], where it is applied following the model for impurity poloidal asymmetry proposed in [45].

As anticipated in the previous subsection, electromagnetic effects can also impact the turbulent transport of impurities. Magnetic flutter usually provides a small to negligible contribution to the total flux, due to the low thermal velocity of the impurities, while the modifications that finite β effects have on the growth rate and frequency of the unstable modes can modify the E × B transport [109, 110]. The behaviour of the impurities when entering a regime with dominant kinetic ballooning modes (KBMs) has been specifically investigated in [110], and found to result in small values of turbulent convection, mainly leading to flat predicted impurity density profiles. Indications that KBM are the dominant instability close to the magnetic axis in H–mode plasmas have been recently obtained by means of gyrokinetic simulations [129], a result that could motivate renewed consideration on the properties of impurity transport in the presence of this type of turbulence. Comparisons among the results of nonlinear gyrokinetic simulations and the predictions of gyrokinetic and/or fluid quasi-linear results have been studied as a function of magnetic shear [93], as well as in TEM and ITG turbulence driven by temperature gradients [130] and in TEM turbulence driven by the density gradient of the main plasma [131]. Consistency between the quasi–linear results and the results from the nonlinear simulations has been found in general although gyro-kinetic quasi-linear results exhibit a slightly stronger dependence on the impurity charge as compared to the nonlinear results. Consistently, in both nonlinear and quasi-linear results, the increase of impurity charge implies an increase and a decrease of the predicted impurity peaking factor in ITG turbulence and TEM turbulence, respectively. This difference can be understood through the opposite sign (outward in ITG and inward in TEM) of the thermodiffusion component, which decreases with increasing charge as presented in section 3.2. Density gradient driven TEM turbulence exhibits impurity density gradients which are weaker than those of the bulk plasma. Consistently with the general properties of impurity turbulent transport, also in these studies [130, 131], the impurity peaking predicted by turbulent transport only remains significantly smaller than the peaking provided by neoclassical transport alone. Further comparisons in the verification of a quasi–linear fluid turbulent transport model against nonlinear gyrokinetic simulations have been performed in [132] where both the electron and the impurity density gradients were consistently obtained from the condition of zero particle fluxes. The low collisionality, turbulence generated, density peaking of the bulk plasma [95] was regularly found to be stronger than that featured by the impurities, confirming the favorable behaviour for bulk plasma and impurity density profiles of low collisionality reactor–like conditions, which are dominated by turbulence. An analytical model for quasi-linear impurity transport has been developed in [100]. This approach appears promising for extensions towards a more complete description for fast applications in transport modelling. Companion studies have been devoted to impurity transport in the presence of TEM [101]. At the transition between ITG and TEM turbulence, more sophisticated effects can occur, which lead to nonlinear results in quantitative disagreement with quasi–linear results that are based on the dominant linear mode only [118]. This property can be understood as a result of the effect of subdominant linear modes, which can determine the dominant properties of the impurity turbulent convection in the case the dominant linear mode has a turbulent convective component which is very small. In these peculiar but experimentally still relevant conditions, quasi-linear models based on the computation of transport combining the effects of both dominant and subdominant linear modes provide results which are closer to the nonlinear results.

The impact of the inclusion of an accurate description of the plasma shape through a numerical magnetic equilibrium solution in comparison with simplified analytical descriptions has been investigated in [133]. Regarding the electromagnetic effects connected with fluctuations of the magnetic field and produced by an increase of the normalized plasma pressure β, with usual ITG modes at sufficiently large mode frequency, the reduction of the growth rate, which is obtained when the threshold for dominant kinetic ballooning modes (KBMs) is approached, leads to a significant decrease of the impurity diffusion coefficient. This reduction is larger than that of the corresponding main ion heat conductivity [109, 110, 134]. This result is also consistent with [113], where the reduction of impurity diffusivity by finite β has been also confirmed in nonlinear gyrokinetic simulations. The effect is practically absent in the case of TEM instabilities, since these are weakly modified by finite β. The impact on the convection is usually weaker, and determined by a more complex combination of multiple effects in the various convective components. In general, however, in the presence of ITG modes at sufficiently large mode frequency, a reduction of the growth rate leads to a stronger relative reduction of the pure convection term with respect to thermo-diffusion and roto-diffusion, which in contrast have a relatively increased effect. Since in these conditions the pure convection is directed inward, whereas both thermo-diffusion and roto-diffusion are directed outward, the reduction of the growth rate can lead to an overall less inward or more outward total convection of the impurity. This effect on the impurity transport has been specifically analyzed also in combination with the stabilization produced by the inclusion of a population of fast ions, like that produced by neutral beam injection heating [134]. As already mentioned, effects connected with the presence of fast ions with temperature anisotropy have been also studied in relation with the consequent development of poloidal asymmetries of the electrostatic potential and the impurity density [102–104]. An overall outward effect in the convection is predicted by these works, which are more specific to the role of ICRH resonant minority ions. The combination of these effects in connection with the presence of beam ions has not been investigated yet to our knowledge. The impact of the consequent poloidal asymmetry of the electrostatic potential on the poloidal distribution is however an exponential function of the impurity charge [45], and therefore becomes negligible for light impurities, for which turbulent transport is usually dominant.

Comparisons of impurity turbulent transport predictions with experimental observations are better performed on light impurities, which have lower levels of collisional transport. First attempts were mainly limited by relatively incomplete theoretical descriptions, in particular still ignoring the presence of roto-diffusion (e.g. [127, 135, 136]), and/or affected by the neglect of the impact that poloidal asymmetries can have on the neoclassical transport component of heavy impurities (e.g. [137]). In one of the first comparisons between dedicated experiments and theoretical predictions, an increase of the central diffusivity of nickel was observed in Tore Supra [138] with increasing electron temperature logarithmic gradient in the central region of the plasma, well beyond the neoclassical levels and consistent with theoretical quasi–linear calculations of impurity diffusivity, as produced by the turbulence destabilization. Additional important steps towards the validation of turbulent transport predictions of impurity transport with experimental observations have been performed in [119, 139] with complementary approaches. The former study compares quasi–linear gyrokinetic predictions with a set of measured density profiles of intrinsic boron, demonstrating the non-negligible importance of the inclusion of roto-diffusion [97] in recovering experimental trends produced by a variation of the electron heating fraction in the presence of background NBI heating, although not reaching a full quantitative agreement. The latter compares, for the first time, global non–linear gyrokinetic simulations with the gyrokinetic code GYRO [140], matching the heat fluxes, with the separate measurements of the diffusion and convection coefficients of a Ca impurity injected in dedicated laser ablation experiments [141]. The comparison provides an unprecedented demonstration that, within the experimental uncertainties, the gyrokinetic predictions were able to reproduce the experimentally estimated values of the impurity transport coefficients. A careful study of the experimental uncertainties was also performed in [121], where the impact of varying plasma toroidal rotation on the measured intrinsic boron density profile was compared to the results of both linear and non–linear gyrokinetic simulations. Within uncertainties a quantitative consistency was possible, although evidences of disagreement at the margin of the uncertainties were observed in the presence of strong NBI heating and strong plasma rotation, with observed hollow boron density profiles being predicted with lower level of profile hollowness. Similar results were also obtained in analogous validation works at JET, consistently documenting the incapability of the gyrokinetic modelling to reproduce the observed hollowness of the intrinsic carbon density profiles [133, 142]. The consistency between predictions from a fluid model and the gyrokinetic results were also verified in these works, and the predicted dependence of the density peaking of injected impurities around mid-radius has been found to be in qualitative agreement with the experimental observations, with an increase of the peaking factor with increasing impurity charge for Z $\leqslant$ 20 and a saturation to constant values for larger values of Z, much below the neoclassical predictions, consistent with a charge independent scaling of impurity transport at large charge, as explained in the previous section.

The inconsistency with experimental observations of hollow intrinsic impurity density profiles was also found in more recent studies at JET C–Wall [143] and JET–ILW [144], where the gyrokinetic simulations were found to predict weakly peaked impurity density profiles in conditions where the experimental observation had flat or hollow profiles. The increasing evidence of a disagreement between predictions and observations in these conditions has concomitantly motivated further experimental investigations not only on intrinsic impurities profiles [145] but also on injected impurities with transient transport experiments, separately determining diffusion and convection [122, 146], with comparisons with both quasi-linear gyrofluid and quasi-linear and nonlinear gyrokinetic simulations respectively. In [146], dedicated to high confinement regimes at DIII–D without ELMs, cases of agreement between theoretical predictions and experimental observations where mainly identified in the convective term, whereas evidences of disagreement were reported on the diffusivities. In [122], an ECRH and NBI heated H-mode plasma at ASDEX Upgrade was found to have boron diffusion and convection coefficients quantitatively consistent with the gyrokinetic predictions within experimental uncertainties.

Systematic comparisons over databases of observations in ASDEX Upgrade [145] and JET C-Wall [143] H-mode plasmas provided important informations in the identification of parameter domains in which the disagreement between theory and experiment was particularly large. In JET [143] the disagreement was found to be particularly large at high collisionality. At ASDEX Ugprade a new framework for the evaluation of impurity densities based on measurements from charge exchange recombination spectroscopy was recently developed [147]. In this device the disagreement was documented to be particularly large at high ion temperature logarithmic gradient and rotation gradients [145], correlating with large NBI heating power. The ASDEX Upgrade cases with higher ion temperature and rotation gradients were close to the JET parameter space investigated in [143] and confirmed the overprediction of impurity density peaking found in that previous work.

These evidences of disagreement motivated in particular the investigation of the direct impact of a neoclassical equilibrium on the turbulent impurity transport [123] and, more recently, the impact of beam ions [134]. While the former study has identified only minor effects, which become negligible for light impurities at core plasma conditions, the latter study has clearly identified an important additional ingredient which was missing in the modelling activity performed so far and which leads to a strong reduction or even reversal of the convection at high values of beam ion density and gradients, due to a combination of effects involving both turbulent and neoclassical transport. Further comparisons between theoretical predictions including the impact of beam ions and experimental observations are currently being performed in order to conclusively assess whether the inclusion of the impact of beam ions is the only critical ingredient that was still missing in order to bring theoretical predictions and experimental measurements to agree within uncertainties. The difficulty in quantitatively predicting the convection of the impurity transport has been confirmed by a recent study in C-Mod [148], where theoretical predictions have been compared with calcium impurity density profiles as well as separate diffusive and convective transport coefficients experimentally inferred by means of a fully Bayesian approach. Clear observations of outward impurity convection and hollow density profiles have been obtained which cannot be reproduced by quasi–linear and non–linear gyrokinetic simulations [149]. These results suggest that while beam ions certainly are an important ingredient, they could not be the only element which is still missing in the theoretical description to reliably reproduce the occurrence of hollow density profiles of light impurities. Very recently, theoretical studies with a linear gyrokinetic model solving for multiple eigenmodes have been also dedicated to the possibility that impurity modes [82, 83, 150] produced by hollow impurity density profiles can sustain a turbulent screening effect [151]. If impurity modes are destabilized, what requires sufficient impurity concentration in combination with a sufficiently large reversed impurity density gradient, these modes could produce the dominant turbulent impurity flux, even if the impurity modes are linearly subdominant with respect to the temperature gradient driven instabilities. In these conditions, impurity modes can sustain hollow impurity density profiles in the presence of sufficiently large electron and ion temperature gradients [151]. Transient destabilization of an impurity mode has also been obtained in the transport modelling of impurity laser ablation experiments with a gyrofluid transport model [152], in this case producing a large impurity flux directed inward.

The impact of a change of main ions, such as a hydrogen isotope, on the impurity transport is a topic which has received limited consideration so far in the comparison between theoretical predictions and experimental observations. From the experimental standpoint it is challenging to isolate pure isotopic effects, that is only deriving from a change of mass of the main ion, on the impurity behaviour from effects which are connected with the consequent changes in plasma parameters. These are also produced by a change of the main ions and could be dominant. Direct isotopic effects can be isolated in theoretical studies. An analytical quasi–linear transport model has been applied to investigate the differences in impurity transport produced by a change of main ions, considering hydrogen, deuterium and tritium, in the presence of ITG modes in a simplified sheared slab geometry [153]. The results show that the change of hydrogen isotope has limited effects on the impurity diffusion and some effects on the impurity convection, in particular with a different behaviour for light, low-Z and heavy, high-Z impurities. An increase of main ion mass implies a reduction of outward impurity convection for light impurities, and an increase of the convection for heavy impurities. More studies on this topic could certainly be of interest, where nonlinear effects can also be important, in particular considering initial operation in H and possibly in He in ITER.

Finally, a specific paragraph has to be dedicated to the important topic of the density profiles and the transport of He as an impurity. Past observations of He density profiles as impurity in deuterium plasmas revealed that He density profiles are very similar to the profiles of the main plasma electron density and with He diffusivities which are comparable to the effective heat conductivity $D_\textrm{He} / \chi_\textrm{eff} \sim 1$ [154]. More recently, applying a specifically developed kinetic forward model for the He plume emission [155], the properties of the He density profiles have been characterized in a dedicated study at ASDEX Upgrade [145], which largely confirmed past results, and documented that the shape of the He density profiles follows that of the electron density, but becomes slightly less peaked than the electron density in conditions of large NBI heating fraction. All of the experimental investigations robustly confirm the absence of He accumulation in the plasma center in all the explored confinement regimes, qualitatively consistent with the general concept that transport of low charge impurities is dominated by turbulence and turbulence is theoretically predicted to not produce central accumulation. Recently, from the modelling standpoint, the transport of He has been specifically considered in [89, 112]. The former work focused on the mechanisms leading to the formation of inward convection, with the support of gyrokinetic simulations. The latter work studies the properties of He transport in combination with those of α particles particularly around the operational point provided by the ITER baseline scenario. The He diffusivity is found to be about $D_\textrm{He}/\chi_\textrm{eff} \approx 2 - 3$ in these conditions, which have an electron to ion heat flux ratio around 1. Diffusion coefficients of He which are equal or larger than the effective heat conductivity have been also found in ITG [142] and TEM turbulence simulations [130], with the TEM turbulence providing the larger values. Anomalous values of diffusion and convection of He have been measured in dedicated experiments in MAST [156] L–mode plasmas, with the observed He peaking factor in good agreement with quasi–linear gyrokinetic predictions in conditions of dominant TEM turbulence. More in general, the ratio $D_\textrm{He}/ \chi_\textrm{eff}$ follows the dependence on the electron to ion heat flux ratio which has been pointed out in [113] and reviewed in section 3.1, which also reflects the impact of varying the dominant turbulence conditions from ITG to TEM [114]. Moving from dominant ITG turbulence to conditions of mixed ITG and TEM turbulence is predicted to increase the He diffusivity. Large values of He turbulent diffusion in collisionless TEM turbulence have been recently reported in a dedicated theoretical study, by means of an analytical model [157], with a particularly favorable role of the electron temperature gradient for the removal of He ash. In strong TEM turbulence however the He diffusivity becomes smaller than that of the electrons [114]. Given the central role that the He diffusivity and its ratio to the heat conductivity play in determining the ratio of the He confinement time to the energy confinement time and thereby also the requirements in the He exhaust system in a burning plasma [158, 159], a more complete validation of the theoretical predictions for this critical parameter against accurate experimental measurements should certainly be considered a present priority. From the side of the convective terms of He turbulent transport, in a database of He density profiles at ASDEX Upgrade a systematic underprediction of the observed He density peaking is regularly found [134, 145], when the impact of NBI fast ions is also included. This underprediction of the peaking is also found in a comparison with ³He density profiles at JET [144]. The reasons behind this disagreement, as obtained from a theoretical model which at the same time provides correct quantitative predictions of other light impurity density profiles, are presently still unclear.