Analytical model for collisional impurity transport in tokamaks at arbitrary collisionality

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 well suited for fast integrated modelling applications due to its low computational cost.


Introduction
Impurities are an unavoidable and integral element of fusion plasmas. They play both detrimental and beneficial roles across fusion devices, which can potentially be prohibitive or indispensable for their 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]. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
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 modeling of the many coupled physical processes involved in a complete pulse, which must be simultaneously computationally fast and robust enough. In the integrated modeling 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 [6,9,[11][12][13][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 modeled within a robust description of impurities.
Fluid codes that are typically used in the modeling of impurity transport, such as NCLASS [20] and NEOART [21,22] do not include effects associated with 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 modeling. 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 modeling 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 modeling. Section 6 summarizes this work and provides an outlook for our model.

Collisional impurity fluxes
Three physically distinct collisional flux components arise from the force balance equation of impurity species 'z' [26], resulting in a flux-surface-averaged radial particle flux given by where Π z is the viscosity tensor and F z r 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 B = I(ψ)∇φ + ∇φ × ∇ψ.
The first term on the right of equation (1) corresponds to the BP flux. It is generated by neoclassical viscosity and is dominant at low collisionalities, where collisions are not frequent enough to isotropize the pressure tensor. The second term is the PS 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 that is in turn composed of two contributions, respectively proportional to the main ion density and temperature gradients, so that with the index c = PS, BP, CL denoting each flux component and the total collisional flux being given by Γ coll z = c Γ c z . The transport coefficients D c z , K c z and H c z 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, where q is the safety factor, R is the major radius, v ti is the thermal velocity of the main ion, and ε is the local inverse aspect ratio. The trapped particle fraction is given by the approximate expression For each flux component c, both the diffusion coefficient D c z and the convection coefficient K c z 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 c z = ZD c z ), because the diamagnetic flows that are part of the frictional force scale as the pressure gradient of the species over its charge, such that ∇p i /Z i ≫ ∇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 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 TSC 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 modeling 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.

PS component
We begin our model with the PS component of the collisional impurity flux, which is typically assumed to be 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 where ρ L,z is the Larmor radius of the impurity and ν z is the impurity collision frequency, given by ν z = ν zi (1 + A i /A) −1/2 , where ν zi ∝ g is the impurity-main ion collision frequency. C z 0 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 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 In the poloidally symmetric, non-rotating limit, where n = 1 and N = 1 but b ̸ = 1, we have that C G ̸ = 0 and C U = 0, so the second term in H PS z 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 TSC. However, all modifications introduced to that model in this section are consistent with the poloidally asymmetric and rotating expressions in [10,14]. If n ̸ = 1, all three D PS z , K PS z , H PS z are implicitly modified through C G , and H PS z 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 the NEO at arbitrary collisionality.

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 where V z∥ and V i∥ are the parallel fluid velocities, the main ion parallel heat flow, q i∥ , is proportional to the main ion temperature gradient, and its coefficient, C z 0 , controls the collisionality dependence of the PS flux through the dependence of the coefficient of the main ion temperature gradient H PS z 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 z 0 ≈ 1.5. This dependence on collisionality can be extended at intermediate collisionalities by using a more complete expression for C z 0 . Indeed, an expression for this coefficient was calculated numerically in [26,28] (and corrected in [29]), such that it is given by a function C 2 as where the impurity strength parameter is defined as α = (n z Z 2 )/(n i Z 2 i ). This form of C z 0 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 on 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 in 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 which is likewise plotted in figure 2 on the diamond cyan curve. In our working example, the difference between the 'HSW' and 'new' expressions is small because for Ar +18 , f 1 ≈ 0.68 and A i /A = 0.05, whereas for instance for He +2 , f 1 ≈ 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 z 0 allows us to recover the dependence of the TSC at around 0.01 ⩽ g ⩽ 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.

Ion-electron collisional heat exchange
As the collisionality becomes large enough in the deep PS 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 ∆ = 1/ν * i . This expansion was then used to solve the main ion drift-kinetic equation including the ion-electron collision operator The resulting Fülöp-Helander scaling factor h ie (denoted µ in equation (35) of [31]) is where the main ion parallel heat conductivity 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 z 0 coefficient. Therefore, including the effect of ion-electron collisional heat exchange into our model corresponds once again to a transformation in C z 0 , in particular a scaling by 1/h ie . However, similarly to the 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 h ie is proposed as 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 δ z = δ θνii Z 2 to be O (1). Here, δ θ is the ratio of the poloidal ion gyroradius to a radial length scale, assumed to be very small, and the collisionality parameterν 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.
is then introduced to account for the effect of ion-electron collisional heat exchange, derived at large Z, also for low Z impurities. The 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 The effect of ion-electron collisions becomes significant when g 2 ∼ 1/µ ie , typically at g ⩾ 1. This is shown in the dashdotted 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 TSC.
The addition of a charge-dependent factor for the diffusion coefficient, given explicitly in equation (A.3), to more accurately match the NEO results at low Z, concludes our model for the PS component of the impurity flux. The dependence of the model at low collisionalities is only recovered by an appropriate expression for the BP flux, where it is dominant over its PS counterpart.

Classical flux
We conclude this section on the PS flux by also including a subsection on the classical flux, which has close analogies with the PS component that has just been considered. The transport coefficients of the classical flux are given by [14] with the geometric coefficient Completing the collisionality dependence of the classical flux simply corresponds to the use of the modified C z 0 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 ≈ 1 in equation (21) whereas C G ≈ 2ϵ 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 inout asymmetries [6].

BP 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 surfaceaveraged viscous force, which drives the BP flux, can be related to the poloidal particle and heat flows, V zθ and q zθ respectively, in a similar way to friction in equation (11), through where the viscosity coefficients K α jk , with α = i, z, are calculated with the Hirshman-Sigmar formalism [26]. The BP transport coefficients are then given in terms of the K α jk as Similarly to the PS flux, the coefficient of the main ion temperature gradient of the BP flux, H BP z , is the one with a more complex structure. In contrast to 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 BP 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 modeling 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 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 The y i factors are given explicitly in equations (A.4)-(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 where a BP ≈ 1 for Z ⩽ 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 BP and PS 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 bellshaped dependence on the logarithm of the collisionality, peaking at low or intermediate collisionalities and vanishing at high collisionalities. Increasing the 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 3(a) for the K z 11 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 3(b). We use this particular coefficient in figure 3 because the K α jk are proportional to the species density, so in equation (23) we have that 1/(1/K i 11 + 1/K z 11 ) ≈ K z 11 if n z ≪ n i , but the general behavior is common to all K α jk .

Physical aspects and applications
The complete parameter space that was investigated during this work is (g, Z, A, f t ) ∈ [10 −6 , 10 3 ] × [2, 44] × [4, 184] × [0.14, 0.74]. The TSC 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 → 0 and g → ∞ 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 ∼ 10 −5 and a higher end at g ∼ 10. The trapped particle fraction values correspond to radial positions from the core to the 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 are 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 4(a) and (b) 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 4(a). 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 TSC in the case of boron, shown in figure 4(c), is discussed in section 5.1, where an effect that arises in a particular subset of our parameter space is in general analyzed.
In section 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, section 5.3 presents FACIT, a new routine for collisional transport modeling.

Localized increase in the TSC at intermediate collisionalities
As the charge of the impurities increases, the TSC maintains a similar behavior to the one shown in figures 1 and 2: 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 4(c) 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 TSC 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 its coefficient of the main ion temperature gradient, H BP z . This is why the effect is weak for highly charged impurities, which have a sub-dominant BP component of the flux. The form of H BP z 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 BP z . Each peak has a different magnitude and location, depending on Z and f t . It is the magnitude of the positive term in H BP z and its location in collisionality with respect to the transition to a dominant PS flux, which determines the characteristics of the bump at intermediate collisionalities.
This localized increase in the TSC is predicted by NEO ( figure 5(a)) and the model ( figure 5(b)). The subset of the parameter space corresponding to a significant effect of the bump, which is approximately (g, Z, A, f t ) ∈ [10 −2 , 10 0 ] × [2,10] × [4,20] × [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.

Radial particle fluxes with realistic profiles
We now apply our model to the calculation of collisional impurity fluxes using experimental data from an AUG H-mode plasma, namely AUG discharge #38910 at 2.15 s. 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 6(a). Here, ρ 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 6(b) 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 ρ = 0.6, with the PS component becoming relevant at the edge.
The pedestal region of H-mode plasmas is of particular interest in collisional transport modeling because turbulence is suppressed due to the edge transport barrier that is developed, and therefore collisional transport can be dominant. In figure 6(b), 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 ρ = 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 6(a), but it is analyzed in detail in figure 7. Figures 7(a) and (b) 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 that indicate that the tungsten flux in the pedestal of this AUG H-mode is well modeled with respect to NEO.
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.

FACIT tool for transport modeling
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 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 modeling applications that include impurities, with Python and MATLAB standalone versions.
The FACIT routine takes the kinetic profiles, the main ion and impurity charges and masses, and the magnetic equilibrium as inputs, 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, especially for the third option, where an iterative calculation is applied.

Conclusions
We present 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 modeling 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 AUG 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 z 0 , has been shown to greatly influence the behavior of the PS 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 z 0 coefficient, allowing us to recover the results of NEO at high collisionalities. A localized increase in the TSC 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 is 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 well suited for fast integrated modeling applications, and the FACIT routine has been introduced as a tool for collisional impurity transport modeling.
The magnetic geometry can be considered in FACIT. In the case of the PS component, the full flux surface geometry can be used by applying the corresponding B(r, θ) and R(r, θ), as described in [14]. For the BP component, the use of f t instead of ε allows one to capture the geometric dependence more robustly, as also observed in [36]. Likewise, the collisionality parameter g (also used in [26,29]) is defined in terms of physical quantities without introducing an additional ε term, unlike ν * , allowing us to avoid the intrinsic geometric dependence from the definition of ε. 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 → 0 and f t → 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 nontrace impurities. A possible solution to overcome this limitation could consist of 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 closedflux surfaces up to the separatrix, computing flux-surfaceaveraged 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. The 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 [14,38] 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 highpower, 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.

Data availability statement
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 that support the findings of this study are openly available in a Zenodo repository at https://doi.org/ 10.5281/zenodo.6045346 [39].

Acknowledgments
The authors are grateful to Emiliano Fable for having made available the ITER simulations of [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 in 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.
The resolution of the NEO used in the construction of this database was (N θ = 21, N ξ = 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 ξ = 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 ξ 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 ξ for a collisionality scan and a local inverse aspect ratio scan respectively. N θ = N ξ + 2 is used throughout. 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, 4(c) 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 z 0 friction coefficient of the main ion parallel heat flow, introduced in equations (13) and (15), are given by with i = 1, 2, 3 and the coefficients l i,j given by 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 .
The y 2 factor of the (1, 2) plateau-regime viscosity coefficient of the impurity, defined in equation (28), is given by The y 3 factor of the (1, 2) PS-regime viscosity coefficient of the impurity, defined in equation (29), is given by  (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 BP z . 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 such that when y 1 and y 4 are multiplied we obtain a first-order error in fit(y 1 y 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 2 t + w p,2 (Z)f t + w p, 3 (Z) f t + w p,4 (Z), (A.9) with the charge-dependent coefficients 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) = 1.02 − 1.79 × 10 −3 Z 1 + 6.6 × 10 −13 Z 6.66 . (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 θ 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 θ using the parameter ξ c (n θ ) = δ n (n θ , ρ = 0.5) δ n (n θ = 192, ρ = 0.5) , where δ n is defined as the horizontal asymmetry of the impurity density, i.e. n = n z / ⟨n z ⟩ = 1 + δ n cos θ + ∆ n sin θ + O(δ 2 n , ∆ 2 n , δ n ∆ n ), and ∆ n is the vertical asymmetry. The reference value of δ n uses a number of poloidal grid points of n θ = 192 that is prohibitive for integrated modeling applications. Figure B1(a) shows the dependence of ξ c on n θ , where we see that an n θ ≪ 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 θ is shown in figure B1(b). 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 θ .
A default of n θ = 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 θ = 192, while the execution time is approximately 20 milliseconds, which is an order of magnitude faster than when using n θ = 40.