# Analytical model for collisional impurity transport in tokamaks at arbitrary collisionality

D Fajardo, C Angioni, P Maget, P Manas

► To cite this version:

D Fajardo, C Angioni, P Maget, P Manas. Analytical model for collisional impurity transport in tokamaks at arbitrary collisionality. Plasma Physics and Controlled Fusion, 2022, 64, pp.055017. 10.1088/1361-6587/ac5b4d.cea-04131443

HAL Id: cea-04131443

https://cea.hal.science/cea-04131443

Submitted on 16 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analytical model for collisional impurity transport
in tokamaks at arbitrary collisionality

D. Fajardo<fnref n="1" />, C. Angioni<fnref n="1" />, P. Maget<fnref n="2" /> and P. Manas<fnref n="2" />

E-mail: daniel.fajardo@ipp.mpg.de

Abstract.

The physics governing the collisional transport of impurities in tokamak plasmas can change significantly depending on four main parameters, namely the collisionality, the impurity charge and mass, and the trapped particle fraction, which can vary widely from the core to the edge of a fusion device. We present an analytical model for collisional impurity transport with a consistent dependence on broad scans in these four parameters, showing good agreement with the drift-kinetic code NEO. Radial profiles of collisional fluxes are calculated for different impurity species using ASDEX Upgrade experimental profiles as well as ITER simulated profiles, and they are also compared to NEO. This model is suited for fast integrated modelling applications due to its low computational cost.

Keywords: impurity transport, collisional, tokamak, neoclassical

<footnote n="1">¹ Max-Planck-Institut für Plasmaphysik, D-85748 Garching, Germany</footnote>

<footnote n="2">² CEA, IRFM, F-13108 Saint Paul-lez-Durance, France</footnote>
# 1. Introduction

Impurities are an unavoidable and integral element of fusion plasmas. They play both detrimental and beneficial roles across the fusion device, which can potentially be prohibitive or indispensable for its operation. An accumulation of impurities in the core needs to be avoided, due to the deleterious effects of radiative losses and fuel dilution. On the other hand, a controlled injection of different impurity species for the radiative cooling of the divertor or disruption mitigation systems is expected to be essential for the operation of next generation fusion devices like ITER and DEMO [1, 2].

For the operation of these reactors, it will be necessary to simulate complete plasma discharges from the individual pulse schedules, in order to test not only the physical effects of the planned discharge configuration parameters, but also the control systems and safety of the discharge itself [3]. This requires the integrated modelling of the many coupled physical processes involved in a complete pulse, which must be simultaneously computationally fast and robust enough. In the integrated modelling of a full discharge, the role of impurities is extremely important, as demonstrated for instance in a DEMO simulation of the feedback control of power through the separatrix using Xenon gas puffing and the effect of a disturbance caused by a tungsten flake falling into the plasma [4].

In contrast to the fuel ions and electrons, for which turbulent transport is typically dominant, the collisional transport of impurities is comparable or even dominant, in particular as the charge Z and the mass A of the impurity increase. Impurity species covering a broad range of charges and masses have different effects across the radius of a tokamak. Likewise, the other two parameters that strongly influence collisional transport, namely the collisionality and the trapped particle fraction, can change significantly from core to edge. The collisionality can particularly modify the relative magnitude of the different physical phenomena at interplay in collisional transport [5].

The non-uniformity of the impurity density distribution on the flux surfaces can significantly affect neoclassical transport, typically enhancing it to exceed turbulent transport but possibly also reducing it to classical levels, depending on the localization and magnitude of the asymmetry [6, 7]. These asymmetries can be friction-induced [8, 9, 10] or caused by strong toroidal rotation [11, 12, 9, 13, 6, 14] or by temperature anisotropies introduced by ion cyclotron resonance heating (ICRH) [15, 16]. They have been extensively observed across multiple devices [17, 18, 19], and must be modelled within a robust description of impurities.

Fluid codes that are typically used in the modelling of impurity transport, such as NCLASS

[20] and NEOART [21, 22] do not include effects associated to poloidally asymmetric impurity densities and are therefore not well suited to model high-Z impurities, for which these asymmetries are stronger. The widely-used drift-kinetic code NEO [23, 24, 25], which does include such effects, requires amounts of computation that are prohibitive for integrated modelling. Therefore, a fast analytical model that is able to calculate collisional impurity fluxes accurately over the broad collisional parameter space in a tokamak, including the effects due to poloidally asymmetric impurity densities, would be a valuable element in the integrated modelling of fusion plasmas.

In this work, we complete the collisionality dependence of a recent analytical model for the Pfirsch-Schlüter (PS) flux of impurities with poloidally asymmetric densities, derived in [10, 14]. This allows for its application at arbitrary collisionality, impurity charge and mass, thereby relaxing the condition of having a highly collisional heavy impurity but a collisionless main ion species. The PS flux is complemented by an analytical model for the banana-plateau (BP) flux, which depends additionally on the trapped particle fraction.

The rest of this paper is organized as follows: section 2 provides a general overview of collisional impurity transport. The collisionality dependence of the PS flux is completed in section 3 by introducing appropriate friction coefficients and considering the role of ion-electron collisional heat exchange on the impurity flux. A fully analytical model for the BP flux is obtained in section 4 through a set of new expressions for the viscosity coefficients. Section 5 includes a discussion of a localized decrease in magnitude and even reversal in the temperature screening effect that appears in a small subset of the parameter space, as well as the application of our resulting model to the calculation of radial fluxes of tungsten with ASDEX Upgrade (AUG) experimental data and ITER simulated profiles, and finally the introduction of FACIT, a new routine for collisional impurity transport modelling. Section 6 summarizes this work and provides an outlook for our model.

# 2. Collisional impurity fluxes

Three physically distinct collisional flux components arise from the force balance equation of an impurity species 'z' [26], resulting in a flux-surface-averaged radial particle flux given by

$$ \Gamma_z = - \frac{I}{Ze} \frac{\langle \mathbf{B} \cdot \nabla \cdot \mathbf{\Pi}_z \rangle}{\langle B^2 \rangle} - \frac{I}{Ze} \left\langle \frac{\mathbf{F}_z^\mathrm{r} \cdot \mathbf{B}}{B^2} \right\rangle + \left\langle \frac{\mathbf{F}_z^\mathrm{r}}{Ze} \cdot \frac{\mathbf{B} \times \nabla \psi}{B^2} \right\rangle, \quad (1) $$
where $\Pi_z$ is the viscosity tensor and $F_r^z$ is the frictional force acting on the impurity, assumed to be generated mostly by collisions with the main ion 'i'. The equilibrium axisymmetric magnetic field is parametrized as usual like $\mathbf{B} = I(\psi)\nabla\varphi + \nabla\varphi \times \nabla\psi$.

The first term on the right of equation (1) corresponds to the banana-plateau flux. It is generated by neoclassical viscosity and it is dominant at low collisionalities, where collisions are not frequent enough to isotropize the pressure tensor. The second term is the Pfirsch-Schlüter flux, driven by parallel friction (with respect to the direction of the magnetic field) acting on the guiding center orbits, and it is dominant at high collisionalities. The final term corresponds to the classical flux, caused by perpendicular friction. It is typically much smaller than the BP and PS fluxes, which together constitute the neoclassical flux.

The radial fluxes can be expressed in terms of the thermodynamic gradients through constitutive relations between the frictional and viscous forces and the fluid particle and heat flows, which are mediated by friction and viscosity coefficients. In the Hirshman-Sigmar formalism of [26], these coefficients are calculated by expanding the first order perturbation of the distribution function in a generalized Laguerre polynomial basis, truncating at second order, and using this expansion to solve the corresponding kinetic expressions, closing the fluid system of equations.

A general form for the components of the collisional impurity flux arises from this formalism. In the absence of rotation, each component is given by a diffusive term, proportional to the impurity density gradient, and a convective term which is in turn composed of two contributions, respectively proportional to the main ion density and temperature gradients, so that

$$ \frac{\Gamma_z^c}{n_z} = -D_z^c \partial_r \ln n_z + K_z^c \partial_r \ln n_i + H_z^c \partial_r \ln T_i, \quad (2) $$

with the index c = PS, BP, CL denoting each flux component and the total collisional flux being given by $\Gamma_z^{\text{coll}} = \sum_c \Gamma_z^c$. The transport coefficients $D_z^c$, $K_z^c$ and $H_z^c$ encompass the entire dependence of the collisional flux on the $(g, Z, A, f_t)$ parameter space. Here, g is the collisionality parameter defined as the ratio between the main ion transit time and the ion-ion collision time,

$$ g = \frac{qR}{v_{ti}\tau_{ii}} = \nu_i^* \epsilon^{3/2}, \quad (3) $$

where q is the safety factor, R is the major radius, $v_{ti}$ is the thermal velocity of the main ion and $\epsilon$ is the local inverse aspect ratio. The trapped particle fraction is given by the approximate expression

$$ f_t = 1 - \frac{(1 - \epsilon)^2}{\sqrt{1 - \epsilon^2 (1 + 1.46\sqrt{\epsilon})}}. \quad (4) $$

For each flux component c, both the diffusion coefficient $D_z^c$ and the convection coefficient $K_z^c$ are positive. For typically negative radial gradients in a tokamak, equation (2) implies an outward diffusive flux and an inward K-convective flux. However, the coefficient of the main ion density gradient is larger than the diffusion coefficient by a factor of the impurity charge ($K_z^c = ZD_z^c$), because the diamagnetic flows that are part of the frictional force scale as the pressure gradient of the species over its charge, such that $\nabla p_i/Z_i \gg \nabla p_z/Z$. Therefore, there is a dominant inward K-convection, which becomes particularly stronger as Z increases.

On the other hand, the coefficient of the main ion temperature gradient can be positive or negative, depending on the values of $(g, Z, A, f_t)$. When $H_z$ is negative, the thermal convection leads to a protective outward flux of the impurity, known as temperature screening. This effect is characterized by the temperature screening coefficient (TSC), defined as

$$ TSC = \frac{H_z}{K_z} = \frac{H_z^{BP} + H_z^{PS} + H_z^{CL}}{K_z^{BP} + K_z^{PS} + K_z^{CL}}. \quad (5) $$

The sign of the TSC indicates the direction of the thermal convection (outwards when it is negative, inwards when it is positive), while its magnitude quantifies its relevance with respect to the K-convection caused by the main ion density gradient, providing a measure of what the ratio between the normalized gradients of the main ion should be in order to have a vanishing collisional flux (assuming a much smaller diffusive flux).

The temperature screening coefficient has a strong, non-monotonic dependence on the collisionality parameter g. This is shown in figure 1, where a collisionality scan of the TSC for Ar$^{+18}$ at mid-radius is made using NEO and NCLASS. In the following sections, we reproduce this dependence with an analytical model that is suited for fast modelling applications, taking the result of NEO (a very complete code which includes the full linearized Fokker-Planck collision operator with multi-species collisional coupling) as the reference.

### 3. Pfirsch-Schlüter component

We begin our model with the Pfirsch-Schlüter component of the collisional impurity flux, which is typically assumed dominant in the case of heavy impurities.

The initial model for the PS flux that we adopt is the analytical treatment of the self-consistent interaction between poloidally asymmetric impurity and main ion densities and the radial flux, derived in [10, 14]. This was done by considering the Helander
ordering of [8] (instead of the conventional neoclassical ordering), maintaining the effect of friction in the parallel force balance equation – which sets the density distributions – while also preserving the diamagnetic velocity of the impurity in the expression for the frictional force. The resulting coupling between the parallel and perpendicular transport equations results in a natural poloidal asymmetry, which is driven by collisions and is present even without rotation or sources of temperature anisotropies. Along with the possibility to include toroidal rotation and the effect of ICRH resonance on the poloidal asymmetries analytically (through the formalism introduced in [16]), this is the main feature of the PS model of [10, 14].

The transport coefficients of this initial model for the PS flux are given by

$$D_z^{PS} = q^2 \rho_{L,z}^2 \nu_z \left( \frac{C_G}{2\epsilon^2} \right), \quad (6)$$

$$K_z^{PS} = \frac{Z}{Z_i} D_z^{PS}, \quad (7)$$

$$H_z^{PS} = \left\{ -\left[ 1 + \frac{Z}{Z_i} (C_0^z - 1) \right] + \frac{C_U}{C_G} \frac{Z}{Z_i} (C_0^z + k_i) \right\} D_z^{PS}, \quad (8)$$

where $\rho_{L,z}$ is the Larmor radius of the impurity and $\nu_z$ is the impurity collision frequency, given by $\nu_z = \nu_{zi}/(1+A_i/A)^{1/2}$, where $\nu_{zi} \propto g$ is the impurity-main ion collision frequency. $C_0^z$ is a coefficient related to friction and $k_i$ is the neoclassical main ion flow coefficient, both of which shall be returned to. The two geometric coefficients

$$C_G = \left\langle \frac{n}{b^2} \right\rangle - \frac{1}{\langle b^2/n \rangle}, \quad (9)$$

$$C_U = \left\langle \frac{n}{N} \right\rangle - \frac{\langle b^2/N \rangle}{\langle b^2/n \rangle}, \quad (10)$$

quantify the effect of the poloidal asymmetry of the magnetic field, which is always present due to the toroidal geometry of the tokamak, and the asymmetries in the density distribution of the impurity and the main ion, through the respective terms $b^2 = B^2/\langle B^2 \rangle$, $n = n_z/\langle n_z \rangle$, $N = n_i/\langle n_i \rangle$.

In the poloidally symmetric, non-rotating limit, where $n = 1$ and $N = 1$ but $b \neq 1$, we have that $C_G \neq 0$ and $C_U = 0$, so the second term in $H_z^{PS}$ vanishes. This is the limit in which we complete the collisionality dependence throughout the following subsections. As shown in figure 1, the model of [10, 14] is not able to reproduce the collisionality dependence of the temperature screening coefficient. However, all modifications to that model introduced in this section are consistent with the poloidally asymmetric and

rotating expressions in [10, 14]. If $n \neq 1$, all three $D_z^{PS}$, $K_z^{PS}$, $H_z^{PS}$ are implicitly modified through $C_G$, and $H_z^{PS}$ is also explicitly modified by the non-zero second term in equation (8). Then, an appropriate expression for $k_i$ must be used to keep the model analytical. We propose the use of the expression constructed in [27], which is accurate with respect to NEO at arbitrary collisionality.

<figure><img src="image_2.png" /><figcaption>Figure 1. Collisionality dependence of the temperature screening coefficient for argon ($Z = 18$, $A = 40$) at mid-radius ($f_t = 0.56$).</figcaption></figure>

### 3.1. Friction coefficients

The Hirshman-Sigmar formalism allows us to calculate a fluid expression for the parallel friction between the impurity and the main ion from its kinetic definition as the first velocity moment of the collision operator. This frictional force is then given by a linear combination of parallel particle and heat flows mediated by friction coefficients [26], such that

$$F_{r\parallel}^{z} = -m_z n_z \nu_z \left( V_{z\parallel} - V_{i\parallel} + C_0^z \frac{2q_{i\parallel}}{5p_i} \right), \quad (11)$$

where $V_{z\parallel}$ and $V_{i\parallel}$ are the parallel fluid velocities, the main ion parallel heat flow, $q_{i\parallel}$, is proportional to the main ion temperature gradient, and its coefficient, $C_0^z$, controls the collisionality dependence of the PS flux through the dependence of the coefficient of the main ion temperature gradient $H_z^{PS}$ on it, as shown in equation (8). The PS model in [10, 14] has a limited dependence on collisionality because it assumes the main ions to be in the deep banana regime but the (heavy) impurity to be highly collisional. These approximations lead to a constant $C_0^z \approx 1.5$.

This dependence on collisionality can be extended at intermediate collisionalities by using a more complete expression for $C_0^z$. Indeed, an expression for this coefficient was calculated numerically in [28,
26] (and corrected in [29]), such that it is given by a function $C_2$ as

$$C_0^z(g, \alpha) = C_2(g, \alpha) = 1.5 - \frac{0.29 + 0.68\alpha}{0.59 + \alpha + 1.34/g^2}, \quad (12)$$

where the impurity strength parameter is defined as $\alpha = (n_z Z^2)/(n_i Z_i^2)$. This form of $C_0^z$ is the one also implemented in NEOART, which is used in the impurity code STRAHL [30], whereas it is not included in NCLASS.

The collisionality dependence of the TSC when this expression is used is shown in figure 2, where it is denoted “HSW” (after Hirshman-Sigmar-Wenzel) and plotted in the dashed purple curve. However, equation (12) was derived by assuming a much heavier impurity than the main ion, and it deviates from NEO for light impurities. Therefore, in order to develop a complete model which is applicable to all impurity charges and masses, we introduce a new parametrization on the impurity charge through two fitted coefficients $f_1$ and $f_2$ (given explicitly in equations (A.1) and (A.2)) and an ion-to-impurity mass ratio. This new form of the $C_2$ function is proposed as

$$\tilde{C}_2(g, Z, A) = \frac{1.5}{1 + f_1(Z)(A_i/A)} - \frac{0.29 + 0.68\alpha}{0.59 + \alpha + (1.34 + f_2(Z))/g^2}, \quad (13)$$

which is likewise plotted in figure 2 in the diamond cyan curve. In our working example the difference between the “HSW” and “new” expressions is small because for $Ar^{+18}$, $f_1 \approx 0.68$ and $A_i/A = 0.05$, whereas for instance for $He^{+2}$, $f_1 \approx 2.3$ and $A_i/A = 0.5$ (in a deuterium plasma). The significant improvement by using the factors for low-Z impurities is shown in figure A3, which reproduces figure 2 for $He^{+2}$ instead of $Ar^{+18}$.

The use of these more complete expressions for $C_0^z$ allows us to recover the dependence of the temperature screening coefficient at around $0.01 \le g \le 3$. However, the limit of very high collisionalities provided by NEO is not recovered by the conventional HSW model and requires the inclusion of an additional physical effect. This effect, which is generated by ion-electron collisions, is presented in the following subsection.

## 3.2. Ion-electron collisional heat exchange

As the collisionality becomes large enough in the deep Pfirsch-Schlüter regime, ion-electron collisions and the energy transfer they generate become non-negligible. This can affect the dynamics of the main ion, subsequently affecting the impurity flux through the ambipolarity condition, which is automatically

fulfilled by the collisional particle fluxes, due to momentum conservation and quasi-neutrality [5].

These ion-electron collisions reduce the main ion parallel heat flow by a factor of $h_{ie}$, defined in equation (14). The derivation of this effect is given in detail in [31], where an auxiliary expansion of the first-order main ion distribution function $f_{i1}$ was made in the smallness of $\Delta = 1/\nu_i^*$. This expansion was then used to solve the main ion drift-kinetic equation including the ion-electron collision operator $C_{ie}[f_{i1}]$, instead of neglecting it in favor of ion-ion collisions $C_{ii}[f_{i1}]$, treating $C_{ie}/C_{ii} \sim O(\Delta)$. The resulting Fülöp-Helander scaling factor $h_{ie}$ (denoted $\mu$ in equation 35 of [31]) is

$$\begin{aligned} h_{ie} &= 1 + 3 \frac{m_e q^2 R^2}{m_i \chi_{i\parallel} \tau_{ei}} \\ &= 1 + \frac{96\sqrt{2}}{125} \frac{1}{Z_i^2} \sqrt{\frac{m_e}{m_i}} \left(\frac{T_i}{T_e}\right)^{3/2} g^2, \end{aligned} \quad (14)$$

where the main ion parallel heat conductivity $\chi_{i\parallel} = (125 T_i \tau_{ii})/(32\sqrt{2} m_i)$ and the definition of $g$ were used.

The main ion parallel heat flow affects the impurity flux through the impurity-main ion frictional force, mediated by the $C_0^z$ coefficient. Therefore, including the effect of ion-electron collisional heat exchange into our model corresponds once again to a transformation in $C_0^z$, in particular a scaling by $1/h_{ie}$. However, similarly to the $\tilde{C}_2$ function, we introduce a parametrization in the impurity charge through a fitted factor $f_3$ (given explicitly in equations (A.1) and (A.2)), such that a new $\tilde{h}_{ie}$ is proposed as

$$\tilde{h}_{ie}(g, Z) = 1 + f_3(Z)\mu_{ie} g^2. \quad (15)$$

This factor arises from the fact that the analytical derivation of [31] is performed by assuming a highly charged impurity, in particular by considering the parameter $\delta_z = \delta_\theta \hat{v}_{ii} Z^2$ to be $O(1)$. Here, $\delta_\theta$ is the ratio of the poloidal ion gyroradius to a radial length scale, assumed to be very small, and the collisionality parameter $\hat{v}_{ii}$ is the ratio of the connection length to the ion-ion mean free path, assumed to be large in the regime where ion-electron collisions are relevant. $\delta_z = O(1)$ then implies $Z = O(\delta_\theta^{-1/2})$, which is large. $f_3(Z)$ is then introduced to account for the effect of ion-electron collisional heat exchange, derived at large $Z$, also for low $Z$ impurities.

The $\mu_{ie} = 96\sqrt{2} \sqrt{m_e/m_i} (T_i/T_e)^{3/2} / (125 Z_i^2)$ coefficient is in general small due to the dependence on the electron to ion mass ratio. The final form of the friction coefficient of the main ion parallel heat flow is thereby given by

$$C_0^z = \frac{\tilde{C}_2(g, Z, A)}{\tilde{h}_{ie}(g, Z)}. \quad (16)$$
The effect of ion-electron collisions becomes significant when $g^2 \sim 1/\mu_{ie}$, typically at $g \ge 1$. This is shown in the dash-dotted pink curve in figure 2, where we can also see that the resulting TSC, shown in the starred orange curve, is able to recover the difference between NEO and the HSW expressions from [26, 29] at high collisionalities. In this deep collisional regime, ion-electron collisional heat exchange leads to an inward thermal convection flux, reversing the sign of the temperature screening coefficient.

The addition of a charge-dependent factor for the diffusion coefficient,

$$D_z^{PS} \rightarrow a_{PS}(Z) D_z^{PS}, \qquad (17)$$

given explicitly in equation (A.3), to more accurately match the NEO results at low Z, concludes our model for the Pfirsch-Schlüter component of the impurity flux. The dependence of the model at low collisionalities is only recovered by an appropriate expression for the banana-plateau flux, where it is dominant over its PS counterpart.

### 3.3. Classical flux

We conclude this section on the PS flux also including a subsection on the classical flux, which has close analogies with the PS component that has been just considered. The transport coefficients of the classical flux are given by [14]

$$D_z^{CL} = \rho_{L,z}^2 \nu_z \left( \frac{C_{G,CL}}{2} \right), \qquad (18)$$

$$K_z^{CL} = \frac{Z}{Z_i} D_z^{CL}, \qquad (19)$$

$$H_z^{CL} = - \left[ 1 + \frac{Z}{Z_i} (C_0^z - 1) \right] D_z^{CL}, \qquad (20)$$

with the geometric coefficient

$$C_{G,CL} = \left\langle \frac{n}{b^2} \right\rangle. \qquad (21)$$

Completing the collisionality dependence of the classical flux simply corresponds to the use of the modified $C_0^z$ coefficient from equation (13), because the same considerations for this coefficient apply for both the parallel and perpendicular friction components.

In the poloidally symmetric limit, the classical flux is lower than the PS flux by the usual factor of $2q^2$, because $C_{G,CL} \approx 1$ in equation (21) whereas $C_G \approx 2\epsilon^2$ in equation (9). However, the classical flux is weakly affected by poloidal asymmetries in the impurity density distribution, and therefore its relative relevance with respect to the PS flux can increase when the latter is suppressed, for instance by weak or intermediate in-out asymmetries [6].

## 4. Banana-Plateau component

Particles in a toroidally-confined magnetized plasma can experience a force due to the anisotropy in the pressure components parallel and perpendicular to the magnetic field, which gives rise to neoclassical viscosity [32]. The parallel surface-averaged viscous force, which drives the banana-plateau flux, can be related to the poloidal particle and heat flows, $V_{z\theta}$ and $q_{z\theta}$ respectively, in a similar way as friction in equation (11), through

$$\langle \mathbf{B} \cdot \nabla \cdot \mathbf{\Pi}_z \rangle = 3 \left\langle \left( \hat{\mathbf{b}} \cdot \nabla B \right)^2 \right\rangle \left[ K_{11}^z V_{z\theta} + \left( K_{12}^z - \frac{5}{2} K_{11}^z \right) \frac{2q_{z\theta}}{5p_z} \right], \qquad (22)$$

where the viscosity coefficients $K_{jk}^\alpha$, with $\alpha = i, z$, are calculated with the Hirshman-Sigmar formalism [26]. The BP transport coefficients are then given in terms of the $K_{jk}^\alpha$ as

$$D_z^{BP} = \frac{3T_i}{2Z^2\epsilon^2 n_z I^2} \left( \frac{1}{1/K_{11}^i + 1/K_{11}^z} \right), \qquad (23)$$

$$K_z^{BP} = \frac{Z}{Z_i} D_z^{BP}, \qquad (24)$$

$$H_z^{BP} = \left[ \frac{Z}{Z_i} \left( \frac{K_{12}^i}{K_{11}^i} - \frac{3}{2} \right) - \left( \frac{K_{12}^z}{K_{11}^z} - \frac{3}{2} \right) \right] D_z^{BP}. \qquad (25)$$

Similarly to the PS flux, the coefficient of the main ion temperature gradient of the BP flux, $H_z^{BP}$, is the one with a more complex structure. Differently from the PS flux, which only depends on the collisionality, charge and mass of the impurity, the BP flux has an additional strong dependence on the trapped particle fraction.

The banana-plateau problem is then reduced to finding an appropriate set of viscosity coefficients. Expressions for the viscosity coefficients in arbitrary collisionality have been derived in [26]. However, they include velocity space averages that require the calculation of computationally-expensive kinetic integrals (as summarized in appendix B of [33]). For a fully analytical model intended for fast modelling applications, bypassing the need to compute these integrals would be advantageous. For the BP flux in our model, we implement instead the approach first derived in [26] and summarized in appendix A of [29]. It consists of solving for the viscosity coefficients in the individual banana, plateau and PS collisionality regimes, evaluating integrals analytically using a Maxwellian distribution function (in contrast to NEO, which uses the full first order perturbed distribution function), and then interpolating the
regimes using a so-called rational approximation, given by

$$K_{jk}^{\alpha} = \frac{K_{jk}^{\alpha, B}}{(1 + K_{jk}^{\alpha, B}/K_{jk}^{\alpha, P})(1 + K_{jk}^{\alpha, P}/K_{jk}^{\alpha, PS})} \quad (26)$$

where the coefficients in the individual regimes are calculated following appendix A of [29]. Naturally, the cumulative effect of these multiple approximations leads to an incorrect behavior with respect to NEO in certain limits of the model. We once again introduce a set of fitted factors to account for these differences, now possibly as functions of the trapped particle fraction in addition to the impurity charge, defined as

$$K_{11}^{z,B} \rightarrow y_1(Z, f_t) K_{11}^{z,B}, \quad (27)$$

$$K_{12}^{z,P} \rightarrow y_2(f_t) K_{12}^{z,P}, \quad (28)$$

$$K_{12}^{z,PS} \rightarrow y_3(Z) K_{12}^{z,PS}, \quad (29)$$

$$K_{12}^{i,B} \rightarrow y_4(Z, f_t) K_{12}^{i,B}. \quad (30)$$

The $y_i$ factors are given explicitly in equations (A.4) to (A.11). Notice that the factor of the viscosity coefficient of the impurity in PS does not depend on the trapped particle fraction, as expected in this regime where the collisions are frequent enough not to allow the particles to complete their orbits.

This new set of viscosity coefficients allows us to complete the dependence of the BP flux on the collisional parameter space. A final charge dependent factor $a_{BP}$ is introduced for the BP coefficient of the main ion temperature gradient, to slightly correct the behavior of this coefficient at very high impurity charges, such that

$$H_z^{BP} \rightarrow a_{BP}(Z) H_z^{BP}, \quad (31)$$

where $a_{BP} \approx 1$ for $Z \le 45$ and $a_{BP} < 1$ for $Z > 45$. High-Z impurities have a BP flux comparable to the PS flux only at very low collisionalities, so $a_{BP}$ is only relevant for instance for tungsten in the core of ITER. This factor is given explicitly in equation (A.12).

The addition of the banana-plateau and Pfirsch-Schlüter components (with a small classical contribution) finalizes our analytical model for collisional impurity transport. The resulting dependence of the TSC on collisionality is shown in the red-square curve in figure 2, where the behavior of NEO is recovered throughout the broad range of collisionalities considered.

The BP flux reduces the magnitude of the temperature screening effect at low collisionalities, although the thermal convective flux is still outwards.

The resulting viscosity coefficients have in general a bell-shaped dependence on the logarithm of the collisionality, peaking at low or intermediate collisionalities and vanishing at high collisionalities. Increasing the

<figure><img src="image_3.png" /><figcaption>Figure 2. Collisionality dependence of the temperature screening coefficient for Ar$^{+18}$ ($Z = 18, A = 40$) at a constant trapped particle fraction $f_t = 0.56$.</figcaption></figure>

charge of the impurity leads to an increase in the magnitude of the peak, but also shifts it to lower collisionalities, as shown in figure 3a for the $K_{11}^z$ viscosity coefficient. On the other hand, increasing the trapped particle fraction decreases the magnitude of the peak and shifts it to higher collisionalities, as shown in figure 3b. We use this particular coefficient in figure 3 because the $K_{jk}^\alpha$ are proportional to the species density, so in equation (23) we have that $1/(1/K_{11}^i + 1/K_{11}^z) \approx K_{11}^z$ if $n_z \ll n_i$, but the general behavior is common to all $K_{jk}^\alpha$.

<figure><img src="image_4.png" /><figcaption>Figure 3. Collisionality dependence of the viscosity coefficient $K_{11}^z$ at a) constant $f_t = 0.56$ with increasing impurity charge, and b) constant $Z = 18, A = 40$ with increasing trapped particle fraction.</figcaption></figure>

## 5. Physical aspects and applications

The complete parameter space that was investigated during this work is $(g, Z, A, f_t) \in [10^{-6}, 10^3] \times [2, 44] \times [4, 184] \times [0.14, 0.74]$. The temperature screening
coefficient becomes independent of collisionality at both low and high enough *g*, saturating to constant values (the diffusion coefficient and therefore also both $K_z$ and $H_z$ are all simply proportional to *g* in these limits). This behavior is present in the analytical expressions, so the model can be applied outside the domain in which it has been fitted: the limits as $g \to 0$ and $g \to \infty$ are well behaved. Such saturation is well within the [$10^{-6}, 10^3$] range selected for *g*. In addition, this collisionality range amply covers the values of *g* in current and planned devices, with the lower end at $g \sim 10^{-5}$ and a higher end at $g \sim 10$. The trapped particle fraction values correspond to radial positions from core to edge of a conventional tokamak. The charge and mass ranges cover impurity species from $He^{+2}$ to $W^{+44}$.

The progressive steps taken in the previous sections in order to have a model with accurate PS and BP fluxes over such a broad parameter space were illustrated in figure 2 for the TSC of $Ar^{+18}$ at mid-radius ($f_t = 0.56$). In order to calculate the particle fluxes, the individual $D_z$, $K_z$ and $H_z$ coefficients must be accurate themselves, for the different impurities that might be present in a tokamak plasma. Figures 4a and 4b show the collisionality dependence of $K_z$ and $H_z$, respectively, for $B^{+5}$ at a lower $f_t = 0.32$. The diffusion coefficient $D_z$ is not shown, because it is simply given by $D_z = K_z/Z$, and therefore it has the same shape as $K_z$ but scaled by a factor of the charge ($Z = 5$ in the case of figure 4a). The collisionality dependence for this low-Z impurity is also well reproduced. Note that both $K_z$ and $H_z$ are plotted by removing a factor of $g$, because otherwise their proportionality to the collision frequency would not allow to see the underlying structure in collisionality. The individual BP and PS components of the flux are also plotted for NCLASS and the model (such a splitting is not possible in NEO).

The particular shape of the temperature screening coefficient in the case of boron, shown in figure 4c, is discussed in subsection 5.1, where an effect that arises in a particular subset of our parameter space is in general analyzed.

In subsection 5.2 we also provide an application of the model to the calculation of radial profiles of the collisional flux, as the entire plasma parameters and magnetic configuration vary realistically from core to edge. Finally, subsection 5.3 presents FACIT, a new routine for collisional transport modelling.

## 5.1. Localized increase in the temperature screening coefficient at intermediate collisionalities

As the charge of the impurities increases, the temperature screening coefficient maintains a similar behavior to the one shown in figures 1 and 2:

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Collisionality dependence of a) the coefficient of the main ion density gradient, $K_z$, b) the coefficient of the main ion temperature gradient, $H_z$, and c) the temperature screening coefficient, TSC=$H_z/K_z$, for boron ($Z = 5$, $A = 10$) at a low trapped particle fraction ($f_t = 0.32$).</figcaption></figure>

it saturates to a constant negative value at low collisionalities set by the BP component, decreases to a minimum at intermediate collisionalities, after which the PS component takes over and the coefficient changes sign until it saturates again at very high collisionalities. However, depending on the value of the trapped particle fraction, this is not necessarily the case for low-Z impurities. Figure 4c shows how the TSC of boron changes with collisionality at $f_t = 0.32$, breaking with the aforementioned behavior at intermediate collisionalties.

At low $Z$ and low $f_t$, the temperature screening
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Simultaneous dependence of the temperature screening coefficient on collisionality and trapped particle fraction for boron ($Z = 5, A = 10$) with a) NEO and b) the analytical model.</figcaption></figure>

coefficient exhibits a localized increase at intermediate collisionalities, corresponding to the transition between dominant BP and PS fluxes. Due to the shape of this localized increase, we refer to it as a *bump*. This bump in the TSC can reduce or even change the direction of the thermal convection, taking it from a protective outward screening flux to a deleterious inward flux. The effect of the bump becomes significant at low $Z$ and very low $f_t$, as shown in figure 5, but it disappears as the charge of the impurity increases.

It is important to note that this is an effect associated entirely to the BP flux, in particular to its coefficient of the main ion temperature gradient, $H_z^{BP}$. This is the reason why the effect is weak for highly charged impurities, which have a sub-dominant BP component of the flux. The form of $H_z^{BP}$ presented in equation (25) leads to a difference in viscosity coefficients. Due to the structure of the viscosity

coefficients described in section 4 and shown in figure 3, there is a resonance between a positive and a negative peak in $H_z^{BP}$. Each peak has a different magnitude and location, depending on $Z$ and $f_t$. It is the magnitude of the positive term in $H_z^{BP}$ and its location in collisionality with respect to the transition to a dominant PS flux what determines the characteristics of the bump at intermediate collisionalities.

This localized increase in the TSC is predicted by NEO (figure 5a) and the model (figure 5b). The subset of the parameter space corresponding to a significant effect of the bump, which is approximately $(g, Z, A, f_t) \in [10^{-2}, 10^0] \times [2, 10] \times [4, 20] \times [0.14, 0.33]$, is not very large, and it represents the most challenging combination of parameters for the model to reproduce. The relative error of the model with respect to NEO in $H_z$ (and therefore in the TSC, because the simpler expressions for $D_z$ and $K_z$ have low errors) is below 35 % across all $(g, Z, A, f_t)$, but is significantly lower in most of the parameter space as we exit the bump region.

## 5.2. Radial particle fluxes with realistic profiles

We apply now our model to the calculation of collisional impurity fluxes using experimental data from an ASDEX Upgrade H-mode plasma, namely AUG discharge #38910 at 2.15 seconds. For this, measurements of the electron temperature and density and the main ion temperature profiles are used, along with the reconstruction of the magnetic equilibrium for quantities such as the safety factor profile. We assume a trace impurity, with $n_z = f_z n_i$ and a concentration of $f_z = 10^{-7}$, such that the main ion density profile is obtained from the quasi-neutrality condition. Furthermore, the radial profile of the impurity charge is assumed from a simple average coronal approximation, as a function of the electron temperature. The flux of a tungsten impurity under these conditions is shown in figure 6a. Here, $\rho$ is the mid-plane radius normalized to the minor radius of the device.

We note that in this particular example, there is a dominant BP flux in the core of the device, which is hotter and therefore less collisional, whereas the tungsten flux at the colder, more collisional edge is mostly given by the PS flux. This highlights the importance of an accurate analytical description of both flux components.

We also apply our model to ITER simulated profiles, obtained from [34]. Figure 6b shows the profile of the radial collisional flux for tungsten in an ITER-like scenario. In this case, the total flux is once again almost exclusively given by the BP component in the core up to $\rho = 0.6$, with the PS component becoming relevant at the edge.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** a) Radial collisional flux of tungsten in an ASDEX Upgrade H-mode (shot #38910 at 2.15 s). b) Radial collisional flux of tungsten in an ITER-like scenario, magnified between $\rho = 0$ and $\rho = 0.9$ in the inset.</figcaption></figure>

The pedestal region of H-mode plasmas is of particular interest in collisional transport modelling because turbulence is suppressed due to the edge transport barrier that is developed, and therefore collisional transport can be dominant. In figure 6b, the tungsten flux including the pedestal region is shown. Due to the large values of the flux in the pedestal, caused by the stronger gradients of the steep kinetic profiles in this region, visualizing the entire radial profile of the flux becomes difficult. Therefore, its detailed structure throughout most of the radius is presented in the inset with a magnified view up to $\rho = 0.9$. We note two points: first, the model matches the results of NEO for this ITER-like scenario very well, from core to edge and in the pedestal. Second, for the set of simulated profiles from [34], we find an outward radial tungsten flux in the pedestal of ITER

with our model, which is consistent with the results of [35]. In the case of the AUG H-mode, the pedestal region is omitted in figure 6a, but it is analyzed in detail in figure 7. Figures 7a and 7b show the transport coefficients $K_z$ and $H_z$, respectively, for tungsten in a finer discretization in the pedestal. Again, $D_z$ is simply $K_z/Z$. The good agreement between NEO and the model in $K_z$ and $H_z$ and the fact that both codes use the same gradients indicate that the tungsten flux in the pedestal of this AUG H-mode is well modelled with respect to NEO.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Radial profile of a) the coefficient of the main ion density gradient, $K_z$, and b) the coefficient of the main ion temperature gradient, $H_z$, for tungsten at the pedestal of an ASDEX Upgrade H-mode (shot #38910 at 2.15 s).</figcaption></figure>

Throughout these examples, we see that our analytical model is able to reproduce the results of NEO with high accuracy for realistic plasma profiles and magnetic equilibria across different devices.
### 5.3. FACIT tool for transport modelling

The work in [10, 14] focused on the development of an analytical treatment of the self-consistent dependence between the poloidal distribution of the impurity density and the radial PS and classical particle fluxes, and investigated in detail the effect of poloidal asymmetries at different levels of toroidal rotation and ICRH power. In this paper we focused on completing the collisionality dependence of that model, also assuring a correct behavior at low charges, and developing an analytical model for the BP flux through the derivation of a set of analytical viscosity coefficients with a consistent dependence on the complete ($g, Z, A, f_t$) parameter space, thereby completing the collisional flux components.

The Fast and Accurate Collisional Impurity Transport (FACIT) tool was developed through the combination of [10, 14] and the present paper. It consists of a FORTRAN routine that can be implemented in transport codes for integrated modelling applications that include impurities, with Python and MATLAB standalone versions.

The FACIT routine takes as inputs the kinetic profiles, the main ion and impurity charges and masses, and the magnetic equilibrium, and outputs the collisional transport coefficients and fluxes as well as the variables that characterize the poloidal distribution of the impurity density. Toroidal rotation and ICRH-induced temperature anisotropies can also be given to the model as inputs.

Regarding poloidally asymmetric densities, FACIT has three main options:

* the poloidally symmetric limit.

* Considering poloidal asymmetries in circular geometry with a large aspect ratio. This simplified geometry allows us to calculate the self-consistent dependence of the fluxes and the asymmetries in a completely analytical way, as detailed in [14].

* Including poloidal asymmetries in full flux surface shaped geometry.

The computational performance of FACIT is discussed in appendix B, in particular for the third option, where an iterative calculation is applied.

## 6. Conclusions

We have presented an analytical model for collisional impurity transport in tokamaks that is both accurate enough with respect to more complete neoclassical codes, such as the drift-kinetic solver NEO, and fast enough to be employed in integrated modelling applications. This model was successfully compared to NEO over broad scans in the main four parameters that

control collisional transport, namely the collisionality, the impurity charge and mass, and the trapped particle fraction. Our model also shows satisfactory agreement with NEO in the calculation of radial fluxes when realistic plasma profiles from different devices are used, as shown in particular for ASDEX Upgrade and ITER-like data.

This work has also allowed us to identify and analyze several physical elements that take part in the complex dependence of the collisional transport coefficients on the collisional parameter space. The weight of the main ion parallel heat flow on the friction force, characterized by the friction coefficient $C_0^z$, has been shown to greatly influence the behavior of the Pfirsch-Schlüter flux. A previous analytical expression for this coefficient has been extended to allow for lower impurity charges, relaxing the condition of having a much heavier impurity than the main ion. The effect of ion-electron collisional heat exchange has been included in the model through an additional transformation of the $C_0^z$ coefficient, allowing us to recover the results of NEO at high collisionalities. A localized increase in the temperature screening coefficient was identified to occur in a small subset of the parameter space, in particular for low impurity charges and low trapped particle fractions at intermediate collisionalities, where the transition from dominant BP to dominant PS components of the flux takes place. This so-called bump was explained as the resonance between the multiple viscosity coefficients that constitute the BP coefficient of the main ion temperature gradient. While it is reproduced by the model, it remains the combination of parameters where the relative error with respect to NEO is highest.

Even though the modifications that complete the collisionality dependence of the PS model of [10, 14] were done in its poloidally symmetric, non-rotating limit, they are consistent with the poloidally asymmetric and centrifugal terms of those recent works.

The proposed new model is suited for fast integrated modelling applications, and the FACIT routine has been introduced as a tool for collisional impurity transport modelling.

The magnetic geometry can be considered in FACIT. In the case of the Pfirsch-Schlüter component, the full flux surface geometry can be used by applying the corresponding $B(r, \theta)$ and $R(r, \theta)$, as described in Ref. [14]. For the Banana-Plateau component, the use of $f_t$ instead of $\epsilon$ allows one to capture the geometric dependence more robustly, as also observed in [36]. Likewise, the collisionality parameter $g$ (also used in Refs. [26, 29]) is defined in terms of physical quantities without introducing an additional $\epsilon$ term, unlike $\nu_*$, allowing us to avoid the intrinsic
geometric dependence from the definition of $\epsilon$. The obtained formulae have been also tested on magnetic configurations with different types of plasma shape, such as negative triangularity plasmas, still showing good agreement with NEO. In contrast, they have not been systematically compared with NEO in the limit of low aspect ratio (spherical) tokamaks (that is, at $f_t > 0.75$). The formulae are well behaved in the limiting values $f_t \to 0$ and $f_t \to 1$, so in these cases we still expect reasonable applicability, however the same level of agreement with NEO as for conventional tokamaks cannot be guaranteed.

The present model does not allow one to compute the impurity transport produced by collisions with multiple non-trace impurities. A possible solution to overcome this limitation could consist in considering a selected set of experimentally relevant impurities (e.g. helium, a seeding species and tungsten for a device with tungsten walls) and in running the analytical model multiple times in parallel. However, a more complete and consistent approach could be desirable.

FACIT describes collisional impurity transport for closed flux surfaces up to the separatrix, computing flux-surface-averaged fluxes and transport coefficients, and cannot be applied to open field lines. Scrape-off layer (SOL) transport is not included in the model. Physical processes in the SOL are not described by FACIT and would have to be introduced as a boundary condition at the separatrix. Moreover, the model is derived to be consistent with the NEO results. NEO adopts the conventional local ordering of neoclassical transport theory, which orders small the poloidal Larmor radius [37]. Applications in very steep pedestals and very close to the magnetic axis and the separatrix have to consider this limitation, which is common to all local neoclassical transport models.

An analytical treatment of poloidal asymmetries and centrifugal effects for the BP flux, in analogy to the works of Refs. [38, 14] for the PS flux, is still missing. These effects mainly influence heavy impurities, which are more collisional and thus tend to have a stronger PS component, particularly during the normal operation of current devices. However, during high-power, low-density operation of current devices and particularly during normal operation in future, hotter reactors, even heavy impurities like tungsten might be in the BP collisionality regime, and a more complete description of the BP flux will be required.

## Acknowledgements

The authors are grateful to Emiliano Fable for having made available the ITER simulations of Ref. [34]. This work has been carried out within the framework of the EUROfusion Consortium and has received funding

from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work has been performed within the EUROfusion TSVV11 activities. Fruitful discussions with the members of the TSVV11 group are gratefully acknowledged.

## Data availability

The FACIT routine will be made publicly available through a git platform soon. At present, readers can receive the routine upon request to the authors.

The data supporting the results of this work is openly available in a Zenodo repository at https://doi.org/10.5281/zenodo.6045346.

## Appendix A: fitted factors

The set of factors introduced throughout this work were fitted from a database of 2156 NEO calculations, consisting of 56 collisionality scans of 22 points in $g \in \{10^{-6}, 10^3\}$, for 8 species $(Z, A) \in \{(2, 4); (5, 10); (7, 14); (10, 20); (18, 40); (26, 56); (36, 84); (44, 184)\}$ times 7 trapped particle fraction values $f_t \in \{0.08, 0.19, 0.32, 0.41, 0.56, 0.66, 0.73\}$, and 42 trapped particle fraction scans of 22 points in $f_t \in \{0.14, 0.74\}$, for 7 species (same as before except (36, 84)) at 6 constant values of $v_i^* \in 1.5 \times \{10^{-3}, 10^{-2}, 10^{-1}, 1, 10, 100\}$. The magnetic geometry used to construct the database was Miller geometry with elongation ($\kappa = 1.35$) and triangularity ($\delta = 0.11$) obtained from an AUG standard H-mode magnetic equilibrium at mid-radius, namely discharge #36315 at 4.5 s.

The resolution of NEO used in the construction of this database was ($N_\theta = 21, N_\xi = 19, N_x = 10$) for the number of poloidal points, the number of Legendre polynomials in the pitch angle, and the number of energy polynomials respectively, which is high enough to be unpractical for fast applications. The value of $N_\xi = 19$ has been found to be the largest at which our NEO simulations are stable across the entire parameter space, including the more challenging limiting cases of low $g$, low $Z$ and low $f_t$. For consistency, we use the same $N_\xi$ for the entire database. Similarly important considerations include the limited computational time for multi-species calculations (four species: electrons, main ions, and two identical impurity species necessary to extract the transport coefficients from the NEO flux), and the robustness across the large database constructed. Figures A1 and A2 show the convergence on $N_\xi$ for a collisionality scan and a local inverse aspect ratio scan respectively. $N_\theta = N_\xi + 2$ is used throughout.
<figure><img src="image_9.png" /><figcaption>**Figure A1.** Convergence on $N_\xi$ at different values of collisionality for constant $Z = 34$, $A = 184$, $\epsilon = 0.167$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure A2.** Convergence on $N_\xi$ at different values of the local inverse aspect ratio for constant $Z = 18$, $A = 40$, $v_i^* = 0.15$.</figcaption></figure>

The selected discretization allows us to retain a sufficient accuracy of ~10% throughout most of the parameter space while fulfilling the aforementioned constraints. At extremely low collisionalities, the flux can differ more, but because it is proportional to $g$, its absolute value is also very small. The TSC is more quickly converged, as it becomes independent of $g$ at higher $g$, as shown in figures 2, 4c and A3.

The analytical formulae for this set of fitted factors in terms of the impurity charge and the trapped particle fraction are given in the rest of this appendix.

The $f_1$, $f_2$ and $f_3$ factors in the final expression for the $C_0^z$ friction coefficient of the main ion parallel heat flow, introduced in equations (13) and (15), are given by

$$f_i(Z) = \frac{l_{i,1} + l_{i,2}Z}{1 + l_{i,3}Z^{l_{i,4}}}, \qquad (A.1)$$

with $i = 1, 2, 3$ and the coefficients $l_{i,j}$ given by

$$\begin{pmatrix} -6.8 \times 10^5 & 2.5 \times 10^6 & 6.0 \times 10^5 & 1.6 \\ 88.3 & 10.5 & 0.2 & 2.6 \\ -4.4 \times 10^6 & 2.7 \times 10^6 & 5.3 \times 10^6 & 0.8 \end{pmatrix}. \qquad (A.2)$$

The effect of these $f_i$ factors is most relevant at low impurity charge, as shown in figure A3 in a collisionality scan of the TSC for He$^{+2}$.

<figure><img src="image_11.png" /><figcaption>**Figure A3.** Collisionality dependence of the temperature screening coefficient for He$^{+2}$ ($Z = 2, A = 4$) at a constant trapped particle fraction $f_t = 0.66$.</figcaption></figure>

The charge-dependent factor of the Pfirsch-Schlüter diffusion coefficient, defined in equation (17) is given by

$$a_{PS}(Z) = -0.109Z^{-1} + 0.743Z^{0.06}. \qquad (A.3)$$

The $y_1$ factor of the (1,1) banana-regime viscosity coefficient of the impurity, defined in equation (27), is given by

$$y_1(Z, f_t) = w_{1,1}(Z)f_t^2 + w_{1,2}(Z)f_t + w_{1,3}(Z)\sqrt{f_t} + w_{1,4}(Z), \qquad (A.4)$$

with the four charge-dependent coefficients being given by $\mathbf{w} = \mathbf{Sv_z}$, where $\mathbf{w} = (w_{1,1}, w_{1,2}, w_{1,3}, w_{1,4})$, the matrix $\mathbf{S}$ has the coefficients

$$\begin{pmatrix} 1.2 \times 10^{-5} & -1.0 \times 10^{-3} & 0.018 & 1.3 \\ -5.8 \times 10^{-5} & 4.8 \times 10^{-3} & -0.1 & -5.8 \\ 6.0 \times 10^{-5} & -4.8 \times 10^{-3} & 0.1 & 5.7 \\ -1.3 \times 10^{-5} & 1.0 \times 10^{-3} & -0.02 & -1.1 \end{pmatrix}, \qquad (A.5)$$

and $\mathbf{v_z} = (Z^3, Z^2, Z, 1)$.

The $y_2$ factor of the (1,2) plateau-regime viscosity coefficient of the impurity, defined in equation (28), is given by

$$y_2(f_t) = 21.31f_t^3 - 21.88f_t^2 + 7.32f_t + 0.63. \qquad (A.6)$$

The $y_3$ factor of the (1,2) Pfirsch-Schlüter-regime viscosity coefficient of the impurity, defined in equation (29), is given by

$$y_3(Z) = 8.85Z^{-2.98} - 7.96Z^{-1.82} \\ - 9.27Z^{-1.98} + 8.34Z^{-0.82}. \qquad (A.7)$$

The expression for the $y_4$ factor of the (1,2) banana-regime viscosity coefficient of the main ion,
defined in equation (30), is based on a fit of the product of the $y_1$ and $y_4$ values found for each NEO simulation. This was done because $y_4$ appears only in products with $y_1$, in particular in the calculation of $H_z^{BP}$. Therefore, in order to minimize the propagation of errors that would arise from the product of two fits, the fit of $y_4$ is defined as

$$y_4(Z, f_t) = \frac{y_p(Z, f_t)}{y_1(Z, f_t)}, \qquad (A.8)$$

such that when $y_1$ and $y_4$ are multiplied we obtain a first-order error in fit($y_1y_4$) instead of an amplified, second-order error in fit($y_1$)fit($y_4$). In this sense, $y_p$ is given by

$$y_p(Z, f_t) = w_{p,1}(Z)f_t^2 + w_{p,2}(Z)f_t + w_{p,3}(Z)\sqrt{f_t} + w_{p,4}(Z), \qquad (A.9)$$

with the charge-dependent coefficients

$$w_{p,i}(Z) = \frac{v_{i,1} + v_{i,2}Z^{v_{i,3}}}{1 + v_{i,4}Z^{v_{i,5}}}, \qquad (A.10)$$

where $i = 1, 2, 3, 4$ and

$$v_{i,j} = \begin{pmatrix} 0.12 & 0.47 & 0.94 & 0.12 & 1.2 \\ 0.63 & -2.92 & 1.06 & 0.3 & 1.16 \\ -0.69 & 2.86 & 1.12 & 0.34 & 1.2 \\ 1.8 & -1.01 & 0.96 & 0.52 & 1.11 \end{pmatrix}. \qquad (A.11)$$

Finally, the high-Z factor for the BP coefficient of the main ion temperature gradient, defined in equation 31, is given by

$$a_{BP}(Z) = \frac{1.02 - 1.79 \times 10^{-3}Z}{1 + 6.6 \times 10^{-13}Z^{6.66}}. \qquad (A.12)$$

## Appendix B: Computational performance of the FACIT routine

For the self-consistent calculation of radial collisional fluxes and the poloidal asymmetries in the impurity density in full flux surface shaped geometry, FACIT performs iterative calculations that depend on the poloidal discretization $n_\theta$ that is used. For fast applications, a trade-off between poloidal resolution and computation time must be made. This is analyzed in the present appendix, for the FORTRAN routine of FACIT.

Using a typical radial discretization of transport codes of $n_r = 101$, we analyze the convergence of FACIT on $n_\theta$ using the parameter

$$\xi_c(n_\theta) = \frac{\delta_n(n_\theta, \rho = 0.5)}{\delta_n(n_\theta = 192, \rho = 0.5)}, \qquad (B.1)$$

where $\delta_n$ is defined as the horizontal asymmetry of the impurity density, i.e. $n = n_z/\langle n_z \rangle = 1 + \delta_n \cos\theta + \Delta_n \sin\theta + O(\delta_n^2, \Delta_n^2, \delta_n\Delta_n)$, and $\Delta_n$ is the vertical

asymmetry. The reference value of $\delta_n$ uses a number of poloidal grid points of $n_\theta = 192$ that is prohibitive for integrated modelling applications.

Figure B1a shows the dependence of $\xi_c$ on $n_\theta$, where we see that an $n_\theta \ll 192$ can be used and still maintain an acceptable accuracy. The execution time of FACIT (running on an Intel Xeon E5-2680 v3 2.5GHz CPU) for increasing $n_\theta$ is shown in figure B1b. Clearly, a sensible choice on the poloidal grid discretization is important when the full-geometry, asymmetric option of the routine is used. The execution times of the simplified geometry and poloidally symmetric options are likewise shown, where both are naturally much smaller and independent of $n_\theta$.

<figure><img src="image_12.png" /><figcaption>Figure B1. a) Convergence and b) execution time of the FACIT routine as a function of the poloidal discretization for a typical radial discretization of $n_r = 101$ radial points.</figcaption></figure>

A default of $n_\theta = 20$ in the case of $n_r = 101$ is suggested for FACIT. This value is denoted with the red dotted line in figure B1. At this number of poloidal points, the value of the convergence parameter is within 97 % of its value at $n_\theta = 192$, while the execution time is approximately 20 milliseconds, an order of magnitude faster than when using $n_\theta = 40$.
References

[1] A. Kallenbach et al. “Impurity seeding for tokamak power exhaust: from present devices via ITER to DEMO”. In: *Plasma Physics and Controlled Fusion* 55.12 (Nov. 2013), p. 124041. DOI: 10.1088/0741-3335/55/12/124041. URL: https://doi.org/10.1088/0741-3335/55/12/124041.

[2] E. M. Hollmann et al. “Status of research toward the ITER disruption mitigation system”. In: 22.2 (Nov. 2014), p. 021802. DOI: 10.1063/1.4901251. URL: https://doi.org/10.1063/1.4901251.

[3] W. Treutterer et al. “Concepts of the new ASDEX Upgrade flight simulator”. In: *Fusion Engineering and Design* 146 (Sept. 2019), pp. 1073–1076. DOI: 10.1016/j.fusengdes.2019.02.008. URL: https://doi.org/10.1016/j.fusengdes.2019.02.008.

[4] F. Janky et al. “Simulation of burn control for DEMO using ASTRA coupled with Simulink”. In: *Fusion Engineering and Design* 123 (Nov. 2017), pp. 555–558. DOI: 10.1016/j.fusengdes.2017.04.043. URL: https://doi.org/10.1016/j.fusengdes.2017.04.043.

[5] P. Helander and D.J. Sigmar. *Collisional Transport in Magnetized Plasmas*. 1st ed. Cambridge University Press, 2002. ISBN: 978-0-521-80798-2.

[6] C. Angioni and P. Helander. “Neoclassical transport of heavy impurities with poloidally asymmetric density distribution in tokamaks”. In: *Plasma Physics and Controlled Fusion* 56.12 (Nov. 2014), p. 124001. DOI: 10.1088/0741-3335/56/12/124001. URL: https://doi.org/10.1088/0741-3335/56/12/124001.

[7] C. Angioni. “Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments”. In: 63.7 (May 2021), p. 073001. DOI: 10.1088/1361-6587/abfc9a. URL: https://doi.org/10.1088/1361-6587/abfc9a.

[8] P. Helander. “Bifurcated neoclassical particle transport”. In: *Physics of Plasmas* 5.11 (Nov. 1998), pp. 3999–4004. DOI: 10.1063/1.873121. URL: https://doi.org/10.1063/1.873121.

[9] T. Fülöp and P. Helander. “Nonlinear neoclassical transport in a rotating impure plasma with large gradients”. In: *Physics of Plasmas* 6.8 (Aug. 1999), pp. 3066–3075. DOI: 10.1063/1.873593. URL: https://doi.org/10.1063/1.873593.

[10] P. Maget et al. “Natural poloidal asymmetry and neoclassical transport of impurities in tokamak plasmas”. In: *Plasma Physics and Controlled Fusion* 62.2 (Nov. 2020), p. 025001. DOI: 10.1088/1361-6587/ab53ab. URL: https://doi.org/10.1088/1361-6587/ab53ab.

[11] F.L. Hinton and S.K. Wong. “Neoclassical ion transport in rotating axisymmetric plasmas”. In: *Physics of Fluids* 28.10 (1985), pp. 3082–3098. DOI: 10.1063/1.865350. URL: https://doi.org/10.1063/1.865350.

[12] S.K. Wong. “Transport of impure plasma with arbitrary toroidal rotation”. In: *Physics of Fluids* 30.3 (1987), pp. 818–829. DOI: 10.1063/1.866333. URL: https://doi.org/10.1063/1.866333.

[13] P. Helander. “Neoclassical transport in a rotating impure plasma”. In: *Physics of Plasmas* 5.4 (Apr. 1998), pp. 1209–1211. DOI: 10.1063/1.872629. URL: https://doi.org/10.1063/1.872629.

[14] P. Maget et al. “An analytic model for the collisional transport and poloidal asymmetry distribution of impurities in tokamak plasmas”. In: 62.10 (Aug. 2020), p. 105001. DOI: 10.1088/1361-6587/aba7f9. URL: https://doi.org/10.1088/1361-6587/aba7f9.

[15] M.L. Reinke et al. “Poloidal variation of high-Z impurity density due to hydrogen minority ion cyclotron resonance heating on Alcator C-Mod”. In: *Plasma Physics and Controlled Fusion* 54.4 (Mar. 2012), p. 045004. DOI: 10.1088/0741-3335/54/4/045004. URL: https://doi.org/10.1088/0741-3335/54/4/045004.

[16] R. Bilato, O. Maj, and C. Angioni. “Modelling the influence of temperature anisotropies on poloidal asymmetries of density in the core of rotating plasmas”. In: *Nuclear Fusion* 54.7 (Apr. 2014), p. 072003. DOI: 10.1088/0029-5515/54/7/072003. URL: https://doi.org/10.1088/0029-5515/54/7/072003.

[17] M.L. Reinke. “Experimental Tests of Parallel Impurity Transport Theory in Tokamak Plasmas”. PhD thesis. Massachusetts Institute of Technology, 2011.

[18] C. Angioni et al. “The impact of poloidal asymmetries on tungsten transport in the core of JET H-mode plasmas”. In: *Physics of Plasmas* 22.5 (May 2015), p. 055902. DOI: 10.1063/1.4919036. URL: https://doi.org/10.1063/1.4919036.
[19] T. Odstrcil et al. "The physics of W transport illuminated by recent progress in W density diagnostics at ASDEX Upgrade". In: 60.1 (Oct. 2017), p. 014003. DOI: 10.1088/1361-6587/aa8690. URL: https://doi.org/10.1088/1361-6587/aa8690.

[20] W.A. Houlberg et al. "Bootstrap current and neoclassical transport in tokamaks of arbitrary collisionality and aspect ratio". In: *Physics of Plasmas* 4.9 (1997), pp. 3230-3242. DOI: 10.1063/1.872465.

[21] A.G. Peeters. "Reduced charge state equations that describe Pfirsch Schlüter impurity transport in tokamak plasma". In: *Physics of Plasmas* 7.1 (2000), pp. 268-275. DOI: 10.1063/1.873812.

[22] R. Dux and A.G. Peeters. "Neoclassical impurity transport in the core of an ignited tokamak plasma". In: *Nuclear Fusion* 40.10 (2000), pp. 1721-1729. DOI: 10.1088/0029-5515/40/10/304.

[23] E.A. Belli and J. Candy. "Kinetic calculation of neoclassical transport including self-consistent electron and impurity dynamics". In: *Plasma Physics and Controlled Fusion* 50.9 (July 2008), p. 095010. DOI: 10.1088/0741-3335/50/9/095010. URL: https://doi.org/10.1088/0741-3335/50/9/095010.

[24] E.A. Belli and J. Candy. "An Eulerian method for the solution of the multi-species drift-kinetic equation". In: *Plasma Physics and Controlled Fusion* 51.7 (June 2009), p. 075018. DOI: 10.1088/0741-3335/51/7/075018. URL: https://doi.org/10.1088/0741-3335/51/7/075018.

[25] E.A. Belli and J. Candy. "Full linearized Fokker-Planck collisions in neoclassical transport simulations". In: *Plasma Physics and Controlled Fusion* 54.1 (Dec. 2011), p. 015015. DOI: 10.1088/0741-3335/54/1/015015. URL: https://doi.org/10.1088/0741-3335/54/1/015015.

[26] S.P. Hirshman and D.J. Sigmar. "Neoclassical transport of impurities in tokamak plasmas". In: *Nuclear Fusion* 21.9 (1981), pp. 1079-1201. DOI: 10.1088/0029-5515/21/9/003.

[27] A. Redl et al. "A new set of analytical formulae for the computation of the bootstrap current and the neoclassical conductivity in tokamaks". In: *Physics of Plasmas* 28.2 (Feb. 2021), p. 022502. DOI: 10.1063/5.0012664. URL: https://doi.org/10.1063/5.0012664.

[28] A. Samain and F. Werkoff. "Diffusion in tokamaks with impurities in the Pfirsch-Schlüter regime". In: 17.1 (Feb. 1977), pp. 53-64. DOI: 10.1088/0029-5515/17/1/006. URL: https://doi.org/10.1088/0029-5515/17/1/006.

[29] K.W. Wenzel and D.J. Sigmar. "Neoclassical analysis of impurity transport following transition to improved particle confinement". In: *Nuclear Fusion* 30.6 (1990), pp. 1117-1127. DOI: 10.1088/0029-5515/30/6/013.

[30] R. Dux. *STRAHL User Manual*. Max-Planck-Institut für Plasmaphysik. Garching bei München, Germany, 2006.

[31] T. Fülöp and P. Helander. "Nonlinear neoclassical transport in toroidal edge plasmas". In: *Physics of Plasmas* 8.7 (July 2001), pp. 3305-3313. DOI: 10.1063/1.1372179. URL: https://doi.org/10.1063/1.1372179.

[32] G.F. Chew, M.L. Goldberger, and F.E. Low. "The Boltzmann equation and the one-fluid hydromagnetic equations in the absence of particle collisions". In: *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* 236.1204 (July 1956), pp. 112-118. DOI: 10.1098/rspa.1956.0116. URL: https://doi.org/10.1098/rspa.1956.0116.

[33] C.E. Kessel. "Bootstrap current in a tokamak". In: *Nuclear Fusion* 34.9 (1994), pp. 1221-1238. DOI: 10.1088/0029-5515/34/9/i04.

[34] E. Fable et al. "The role of the source versus the collisionality in predicting a reactor density profile as observed on ASDEX Upgrade discharges". In: *Nuclear Fusion* 59.7 (2019), p. 076042. DOI: 10.1088/1741-4326/ab1f28.

[35] R. Dux et al. "Transport of tungsten in the H-mode edge transport barrier of ITER". In: *Plasma Physics and Controlled Fusion* 56.12 (Nov. 2014), p. 124003. DOI: 10.1088/0741-3335/56/12/124003. URL: https://doi.org/10.1088/0741-3335/56/12/124003.

[36] O. Sauter, C. Angioni, and Y. R. Lin-Liu. "Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime". In: *Physics of Plasmas* 6.7 (July 1999), pp. 2834-2839. DOI: 10.1063/1.873240. URL: https://doi.org/10.1063/1.873240.

[37] F. L. Hinton and R. D. Hazeltine. "Theory of plasma transport in toroidal confinement systems". In: *Reviews of Modern Physics* 48.2 (Apr. 1976), pp. 239-308. DOI: 10.1103/revmodphys.48.239. URL: https://doi.org/10.1103/revmodphys.48.239.

[38] E.A. Belli, J. Candy, and C. Angioni. "Pfirsch-Schlüter neoclassical heavy impurity transport in a rotating plasma". In: *Plasma Physics and Controlled Fusion* 56.12 (Nov. 2014), p. 124002. DOI: 10.1088/0741-3335/56/12/124002. URL:
https://doi.org/10.1088/0741-3335/56/12/124002.