# Controlling impurity transport in 5D gyrokinetic simulations using a transport barrier

Guillaume Lo-Cascio, Etienne Gravier, Maxime Lesur, Thierry Réveillé, Yanick Sarazin, Xavier Garbet, Kyungtak Lim, Timothé Rouyer, Alejandro Guillevic

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

To cite this version:

Guillaume Lo-Cascio, Etienne Gravier, Maxime Lesur, Thierry Réveillé, Yanick Sarazin, et al. Controlling impurity transport in 5D gyrokinetic simulations using a transport barrier. Nuclear Fusion, 2025, 65 (5), pp.056021. 10.1088/1741-4326/adc5ea. hal-04845964v3

HAL Id: hal-04845964

https://hal.science/hal-04845964v3

Submitted on 9 May 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Controlling impurity transport in 5D gyrokinetic
simulations using a transport barrier

G. Lo-Cascio*¹,², E. Gravier¹, M. Lesur¹,³, T. Réveillé¹, Y. Sarazin⁴, X. Garbet⁴,⁵, K. Lim⁵, T. Rouyer¹, and A. Guillevic⁶

¹IJL, UMR 7198 CNRS, Université de Lorraine, 54000 Nancy, France

²Max-Planck-Institut für Plasmaphysik, Boltzmannstrasse 2, D-85748
Garching bei München, Germany

³Institut Universitaire de France (IUF), Paris, France

⁴CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France

⁵School of Physical and Mathematical Sciences, Nanyang Technological University, 637371 Singapore

⁶Department of Applied Physics, Chalmers University of Technology, Göteborg, Sweden

May 9, 2025

Abstract

The effects of transport barriers on impurity transport, specifically helium (He) and tungsten (W), are investigated using the global, flux-driven, full-F, 5D gyrokinetic code GYSELA. The transport barrier is induced by triggering *E* × *B* shear via an external poloidal momentum source, thereby stabilizing ITG turbulence and reducing outward heat fluxes. These reductions in particle and heat fluxes due to the transport barrier lead to enhanced confinement, thereby steepening the ion temperature profiles and reducing the heat diffusivity of the main ion species (i.e. deuterium). Impurity transport in both turbulent and neoclassical regimes is investigated under various conditions, with and without the transport barrier, along with a reversed density profile for tungsten. The transport barrier is found to reduce outward impurity transport and enhance neoclassical thermal screening due to the steepened temperature profile. However, it also prevents helium from being flushed out, due to an increased inward Banana-Plateau (BP) flux, caused by the poloidal asymmetry of the source. Overall, a transport barrier induced by an *E* × *B* shear proves to be an effective mechanism not only for reducing heat fluxes but also for controlling impurity transport.

<footnote>*Work financed and conducted at IJL, current adress at IPP.</footnote>
# 1 Introduction

Understanding impurity transport in magnetically confined plasmas is crucial for the successful operation of fusion devices, such as tokamaks and stellarators. Even trace amounts of impurities can degrade plasma performance by increasing radiative losses, diluting the fuel, and causing plasma contamination [1]. In this study, we focus on two impurities; helium [2], a by-product of the fusion reaction, and tungsten [3, 4, 5, 6], widely used as a plasma-facing material due to its high melting point and low sputtering yield. Understanding the transport mechanisms of these impurities, both neoclassical and turbulent [7, 8, 9], is essential for maintaining optimal plasma conditions and achieving sustained fusion reactions.

Edge transport barriers (ETB), observed in high-confinement mode (H-mode) plasmas [10, 11, 12, 13], are regions of significantly reduced turbulent transport, leading to improved confinement of heat and particles. These barriers are typically associated with strong poloidal $E \times B$ shear flows [14, 15, 16] that stabilize turbulent modes, such as ion temperature gradient (ITG) instabilities [17, 18]. A strong correlation between the presence and strength of this transport barrier and impurity transport has been identified in Refs. [19, 20, 21], showing a reduced impurity transport within the ETB due to suppressed plasma turbulence.

Neoclassical thermal screening [7, 22], which is expected to play a substantial role in preventing core accumulation of tungsten, can be triggered by having a strong enough temperature gradient of the main ion species compared to its density gradient. Some experimental evidence has been found in JET [23] suggesting this effect already plays an important role in stopping tungsten from accumulating in the core.

Previous studies using gyrokinetic modelling [24] showed a good agreement between the theory [8] and simulations with a strong inward neoclassical flux for heavy impurities such as tungsten, while light impurities are dominated by turbulent transport. The role of poloidal asymmetries has also been extensively investigated [25] and it was found to play a major role for neoclassical transport of impurities. Additionally, the impact of parallel momentum on tungsten transport has also been investigated [26] and concluded that strong toroidal rotation led to increased in-out poloidal asymmetries, enhancing neoclassical transport of impurities as a result. Studies on lighter impurities, such as helium [2], showed a discrepancy between numerical and experimental peak of density profiles, indicating some elements of helium transport is either not modelled correctly or missing. Some additional investigations on the impact of NBI and fast ions [27], polarization drift and ITG-TEM interplay [28] as well as externally applied torque [29] were also performed to improve and complete helium transport modelling. Integrated modelling of impurity transport also shows promising results with a good agreement between ASDEX-U discharges and the ASTRA code [30], with a lower impurity density peaking with a stronger central wave heating scheme.

In the present work, we extend the previous results of light to heavy impurities transport in case of ITG turbulence, as outlined in Ref. [24], by inducing a transport barrier through poloidal $E \times B$ shear injection [31, 32]. We analyse the effects of transport barrier on impurity transport using the 5D full-F gyrokinetic code GYSELA [33]. Both turbulent and neoclassical transports are considered self-consistently, thereby capturing their synergy on the total impurity transport [34]. We compare various cases with and without a transport barrier, as well as scenarios with positive and negative tungsten density gradients, considering that tungsten is typically introduced from the outboard region via wall
sputtering. We also aim to extract the impact of the transport barrier on the diffusion and convection transport coefficients.

The structure of the paper is as follows. Section 2 provides a brief overview of GY-SELA, the vorticity source used for the transport barrier, as well as essential theoretical elements on neoclassical transport. Section 3 presents the numerical setup and parameters chosen to investigate the effects of transport barrier on impurity transport, focusing on both helium and tungsten. Section 4 shows the transport of helium and tungsten, with a particular emphasis on the neoclassical transport of both species. The thermal screening of tungsten is highlighted as a key component of the transport barrier on particle transport. Additionally, an investigation is performed on the reduction of turbulent transport coefficients, namely diffusion and convection, for tungsten when the transport barrier is present, using two different density gradients for tungsten. Finally, Section 5 presents the conclusion and key findings of this study.

# 2 Model and theory

## 2.1 GYSELA

As a full-F gyrokinetic code, GYSELA [33] enables the study of both neoclassical and turbulent transports with a self-consistent model in 5D phase-space $z_s = (x_{GC,s}, v_{G\parallel,s}, \mu_s)$. The code combines the Vlasov equation with the quasi-neutrality equation, written as:

$$B_{\parallel,s}^{*} \partial_t \bar{F}_s + \nabla \cdot (\dot{\mathbf{x}}_{GC,s} B_{\parallel,s}^{*} \bar{F}_s) + \partial_{v_{G\parallel,s}} (\dot{v}_{G\parallel,s} B_{\parallel,s}^{*} \bar{F}_s) = C(\bar{F}_s) + S(\bar{F}_s), \quad (1)$$

$$\underbrace{e N_{0,e} \left( \frac{\phi - \langle \phi \rangle_{FS}}{T_e} \right)}_{\text{Adiabatic e}^{-} \text{ response}} - \sum_{s} \underbrace{Z_s \nabla_{\perp} \cdot \left( \frac{N_{0,s}}{B_0 \omega_{c,s}} \nabla_{\perp} \phi \right)}_{\text{Polarization density}} = \sum_{s} \underbrace{Z_s \int d^3 v_s \mathcal{J} [\bar{F}_s - \bar{F}_{eq,s}]}_{\text{Particle density fluctuations}}. \quad (2)$$

Here, the subscript $s$ denotes the species considered, $\bar{F}_s$ represents the ion distribution function of the species, $\phi$ is the electrostatic potential. The variables $x_{GC,s}$ and $v_{G\parallel,s}$ denote the gyro-center position and parallel velocity, respectively. The infinitesimal volume element in the velocity space $d^3v_s = J_{v,s}d\mu dv_{G\parallel,s}$ is defined through the velocity space Jacobian $J_{v,s} = 2\pi B_{\parallel}^*/M_s$ and the parallel component of the generalized magnetic field $B_{\parallel,s}^* = \mathbf{b} \cdot \mathbf{B}_s^*$, the latter defined through

$$\mathbf{B}_s^* = \mathbf{B} + \frac{M_s}{q_s} v_{G\parallel,s} \nabla \times \mathbf{b}, \quad (3)$$

with $\mathbf{B} = I(\psi) \nabla\varphi + \nabla\varphi \times \nabla\psi$ where $I$ is the current, $\varphi$ the toroidal angle and $\psi$ the poloidal flux function. The mass and charge of the concerned species are respectively written $M_s$ and $q_s$ while $\omega_{c,s} = q_s B_0/M_s$ is the reference cyclotron frequency. $T_e$ and $N_{0,s/e}$ are respectively the initial electron temperature and initial density of the considered species (i.e. $s$ for ion species and $e$ for electrons). The gyro-average operator is written $\mathcal{J}$ while the collision operator is written $C(\bar{F}_s)$. This operator ensures the conservation of energy, momentum, and density [35], while allowing both intra and inter species collisions, which is relevant for multiple kinetic species as explored in this article. However, it should be noted that electron-ion collisions are neglected in the present work due to the assumption of an adiabatic electron response.
The source terms in equation (1), which are limited to heat and vorticity terms in this study, are embedded in

$$S(\bar{F}_s) = S_{H,s} + S_{\Omega,s}, \qquad (4)$$

with their mathematical construction are detailed in [36, 37]. The choice was made to ensure that the vorticity source only injects vorticity without affecting energy or particles, while generating pressure anisotropy as a consequence. The heat source is essential for GYSELA to perform flux-driven simulations and achieve steady states. To inject both energy and poloidal momentum (i.e. vorticity) in the system, the following kinetic sources are used

$$S_{H,s} = \frac{S_0^H S_r^H (r)}{\sqrt{2\pi^{3/2} T_{Sce,H}^{5/2}}} \left[ \bar{v}_{G\parallel}^2 - \frac{1}{2} - \frac{J_{\parallel,B}}{2 - J_{\parallel,B}^2} (2 - \bar{\mu}) (2\bar{v}_{G\parallel} - J_{\parallel,B}) \right] \exp(-\bar{v}_{G\parallel}^2 - \bar{\mu}), \quad (5)$$

$$S_{\Omega,s} = \frac{S_0^{\Omega} S_r^{\Omega} (r) B^2}{\sqrt{2\pi^{3/2} T_{Sce,\Omega}^{5/2}}} [2\bar{v}_{G\parallel}^2 - \bar{\mu}] \exp(-\bar{v}_{G\parallel}^2 - \bar{\mu}), \qquad (6)$$

where $S_0^{H,\Omega}$ and $S_r^{H,\Omega}$ denote the source amplitudes and radial profiles respectively, both of which being customizable GYSELA input parameters. The normalization of the parallel velocity $\bar{v}_{G\parallel} = v_{G\parallel}/\sqrt{2T_{Sce}/M_s}$ and magnetic moment $\bar{\mu} = \mu B/T_{Sce}$ is done using a normalized source temperature $T_{Sce}$, taken as 1 hereafter. While the vorticity source injects neither heat nor particles, it injects a marginal amount of parallel momentum alongside vorticity. Pressure anisotropy is, as a result, generated at the source location. This can be seen in equation (6) in the terms $2\bar{v}_{G\parallel}^2 - \bar{\mu}$ which, when integrated in the velocity space, will respectively inject parallel pressure and pump perpendicular pressure<fnref n="1" /> and therefore enhance the local pressure anisotropy $\Pi_{s,\parallel} = P_{s,\parallel} - P_{s,\perp}$. The integral over the velocity space of equations 5 and 6, which are the radial profiles of the fluid sources of heat and vorticity respectively, are shown in Figure 1a. A conservation equation for vorticity can be obtained by integrating the gyro-averaged equation (1)

$$\partial_t W + \partial_\psi K = S_0^\Omega \nabla_\perp^2 S_r^\Omega (r), \qquad (7)$$

with $W = -Z_s \langle \nabla_\perp \cdot \left( \frac{N_{0,s}}{B\omega_{c,s}} \nabla_\perp \phi \right) \rangle_{FS} = q_s \langle \int d^3v_s \mathcal{J}[\bar{F}] \rangle_{FS}$ the vorticity, the vorticity flux is written as $\mathcal{K} = q_s \langle \int d^3v_s \mathcal{J}[(d_t\mathbf{x}_G \cdot \nabla\psi)\bar{F}] \rangle_{FS}$ and $S_0^\Omega \nabla_\perp^2 S_r^\Omega (r)$ the fluid vorticity source. The flux-surface average operation is written $\langle \dots \rangle_{FS} = \iint \dots J_\psi d\theta d\varphi / \iint J_\psi d\theta d\varphi$ with $J_\psi = (\mathbf{B} \cdot \nabla\theta)^{-1}$. The vorticity, written $W$ here, is then completely analogous to the polarization density, hence the effect observed on the $E \times B$ velocity. One may note also that the main drive for the fluid vorticity flux $\mathcal{K}$ is the turbulence itself, which means that we require the system to be in a turbulent state in the first place before using this vorticity source, otherwise the system may not be stable.

However, the source is not poloidally symmetric, as shown in the radial profile the vorticity source for different poloidal angle in Figure 1b. This induces poloidal asymmetry in pressure anisotropy, which in turn affects the neoclassical transport channel, especially the Banana-Plateau component, as discussed in the following subsection.

<footnote n="1">¹We recall that $P_{s,\parallel} = \int d^3v_s \frac{1}{2} \mathcal{U}^2 \bar{F}_s$ with $\mathcal{U} = (v_{G\parallel} - V_{s,\parallel})$ where $V_{s,\parallel}$ is the average fluid velocity, $P_{s,\perp} = \int d^3v_s \mu B \bar{F}_s$ and $P_{s,tot} = \frac{1}{3} P_{s,\parallel} + \frac{2}{3} P_{s,\perp}$.</footnote>
## 2.2 Neoclassical transport

Neoclassical transport is of utmost importance regarding transport of heavy impurities, such as tungsten, in modern tokamaks. They usually are in a high collisionality regime, i.e. Pfirsch-Schlüter, but are expected to be in an intermediate plateau regime for larger future reactors like ITER. Previous studies already emphasized the importance of pressure anisotropy, as well as its poloidal asymmetry regarding transport of heavy impurities [8, 25]. By assuming low impurity concentration $C_s = N_s/N_i \le 10^{-5}$ and main ions in the banana collisionality regime, one can compute the expected neoclassical flux of impurities in steady-state.

Taking the first moment of equation (1) in the steady state, we obtain $\nabla \cdot \Gamma = 0$ (i.e. the conservative form of the continuity equation) from which the flux can be decomposed into parallel and perpendicular components, $\Gamma = \Gamma_\parallel \mathbf{b} + \Gamma_\perp$. The perpendicular flux of impurity can be written as:

$$ \Gamma_{\perp,s} = N_s \mathbf{v}_{E \times B} + \langle \mathbf{v}_D \rangle_v - \nabla \times \left[ \left\langle \frac{\mu_s}{q_s} \mathbf{b} \right\rangle_v \right], \quad (8) $$

with

$$ \mathbf{v}_{E \times B} = \frac{\mathbf{B} \times \nabla \phi}{B^2}, \quad (9) $$

$$ \mathbf{v}_D = \mathbf{v}_\nabla + \mathbf{v}_c = \frac{\mu}{q_s} \frac{\mathbf{B} \times \nabla B}{B^2} + \frac{M_s v_{G\parallel}^2}{q_s B^2} \frac{\mathbf{R}_c \times \mathbf{B}}{R_c^2}, \quad (10) $$

$$ \Gamma_{mag} = -\nabla \times \left[ \left\langle \frac{\mu_s}{q_s} \mathbf{b} \right\rangle_v \right] = \mathbf{j}_{mag}/q_s = \nabla \times \mathbf{M}/q_s. \quad (11) $$

Equations (9–11) are respectively the $E \times B$, gradient plus curvature drifts as well as the magnetization flux $\Gamma_{mag}$ where the bracket $\langle \dots \rangle_v = \int d^3v_s \dots \bar{F}_s$ represents an average over the velocity space with $\bar{F}_s$ the distribution function of the considered impurities species and $\int d^3v \bar{F}_s = N_s$, the density of impurities. The magnetization is a mandatory correction to have the full particle flux<fnref n="2" />. Developing equation (8) by using the CGL tensor $\Pi_s$ notation [38] and adding the parallel component of the flux, one gets

$$ \Gamma_s = \left[ N_s v_{\parallel,GC,s} - \frac{P_{\perp,s}}{q_s B} \mathbf{b} \cdot (\nabla \times \mathbf{b}) \right] \mathbf{b} + N_s \frac{\mathbf{b} \times \nabla \phi}{B} + \frac{\mathbf{b} \times \nabla \cdot \Pi_s}{q_s B}, \quad (12) $$

in which we introduced the CGL tensor $\Pi_{ij,s} = P_{\parallel,s}b_i b_j + P_{\perp,s}(\delta_{ij} - b_i b_j)$. The first term of equation (12), aligned with the magnetic field direction $\mathbf{b}$, represents the parallel particle flux whereas the other terms are respectively the $E \times B$ drift and diamagnetic flux terms. By considering poloidal asymmetries [8, 24, 25] in equation (12) in the form of

$$ \nabla P_{\perp,s} = \frac{\partial P_{\perp,s}}{\partial \theta} \nabla \theta + \frac{\partial P_{\perp,s}}{\partial \psi} \nabla \psi, \quad (13) $$

$$ \nabla \phi = \frac{\partial \phi}{\partial \theta} \nabla \theta + \frac{\partial \phi}{\partial \psi} \nabla \psi, \quad (14) $$

<footnote n="2">²This term does not impact the continuity equation; it is divergence free by nature since $\nabla \cdot (\nabla \times \mathbf{A}) = 0$.</footnote>
one can write the total neoclassical flux-surface-averaged impurity flux across magnetic surfaces as

$$ \mathbf{\Gamma}_s = K_s \mathbf{B} - N_s \Omega_s R^2 \mathbf{\nabla} \varphi + \Pi_{\parallel,s} \frac{\mathbf{b} \times \mathbf{k}}{q_s B} + \frac{\mathbf{b} \times \mathbf{\nabla} \theta}{B} \left( N_s \frac{\partial \phi}{\partial \theta} + \frac{1}{q_s} \frac{\partial P_{\perp,s}}{\partial \theta} \right), \quad (15) $$

with $K_s = \frac{\Gamma_{\parallel,s}}{B} + \frac{IN_s \Omega_s}{B}$, $\Omega_s = \frac{\partial \phi}{\partial \psi} + \frac{1}{N_s q_s} \frac{\partial P_{\perp,s}}{\partial \psi}$, $\Pi_{\parallel,s} = P_{\parallel,s} - P_{\perp,s}$ the pressure anisotropy. By using the incompressibility of the flow $\mathbf{\nabla} \cdot \mathbf{\Gamma}_s = 0$, one can determine the function $K_s$ with a few approximations. By assuming low Mach number for the impurities, heavy impurities ($M_s/M_i \gg 1$) and no local friction force [25], one can rewrite equation (15) as

$$ \langle \mathbf{\Gamma} \cdot \mathbf{\nabla} \psi \rangle_{\text{FS}}^{\text{neo}} = \langle \mathbf{\Gamma} \cdot \mathbf{\nabla} \psi \rangle_{\text{FS}}^{\text{BP}} + \langle \mathbf{\Gamma} \cdot \mathbf{\nabla} \psi \rangle_{\text{FS}}^{\text{PS}}, \quad (16) $$

with each channel detailed as

$$ \langle \mathbf{\Gamma} \cdot \mathbf{\nabla} \psi \rangle_{\text{FS}}^{\text{BP}} = \Gamma^{P_{\perp,s}} + \Gamma^{\Pi_{\parallel,s}}, \quad (17) $$

$$ \langle \mathbf{\Gamma} \cdot \mathbf{\nabla} \psi \rangle_{\text{FS}}^{\text{PS}} = \Gamma^{\nabla N_i} + \Gamma^{u_s} + \Gamma^{\nabla n_s} + \Gamma^{K_s}. \quad (18) $$

The different terms then read

$$ \Gamma^{P_{\perp,s}} = -\frac{I}{q_s} \left\langle \frac{B^2}{N_s} \right\rangle_{\text{FS}}^{-1} \left\langle \frac{\mathbf{B} \cdot \mathbf{\nabla} \theta}{N_s} \left[ \frac{\partial P_{\perp,s}}{\partial \theta} \right] \right\rangle_{\text{FS}}, \quad (19) $$

$$ \Gamma^{\Pi_{\parallel,s}} = -\frac{I}{q_s} \left\langle \frac{B^2}{N_s} \right\rangle_{\text{FS}}^{-1} \left\langle \frac{\mathbf{B} \cdot \mathbf{\nabla} \theta}{N_s} \left[ B \frac{\partial}{\partial \theta} \left( \frac{\Pi_{\parallel,s}}{B} \right) \right] \right\rangle_{\text{FS}}, \quad (20) $$

which are respectively terms driven by the poloidal asymmetries of perpendicular pressure and pressure anisotropy. Both terms are proportional to $q_s^{-1}$, indicating a low impact on strongly charged ions, such as tungsten. In contrast, weakly charged impurities, like helium, are more significantly affected by these asymmetries. Additionally, since both terms are independent of the collision frequency $\nu_{si}$, they will dominate in low-collisionality regimes, emphasizing their importance for low-mass, low-charge impurities. Consequently, both terms are expected to be dominant for helium.

The other terms are written as

$$ \Gamma^{\nabla N_i} = \frac{I}{q_s} M_s \nu_{si} \left\{ \frac{T_i}{e} \frac{I}{L_{\psi,i}} \left( \left\langle \frac{N_s}{B^2} \right\rangle_{\text{FS}} - \left\langle \frac{B^2}{N_s} \right\rangle_{\text{FS}}^{-1} \right) \right\}, \quad (21) $$

$$ \Gamma^{u_s} = -\frac{I}{q_s} M_s \nu_{si} \left\{ u \left( \left\langle N_s \right\rangle_{\text{FS}} - \frac{\left\langle B^2 \right\rangle_{\text{FS}}}{\left\langle B^2 / N_s \right\rangle_{\text{FS}}} \right) \right\}, \quad (22) $$

$$ \Gamma^{\nabla N_s} = +\frac{I}{q_s} M_s \nu_{si} \left\{ \frac{T_i I}{e} \left( \left\langle \frac{N_s}{B^2 L_{\psi,s}} \right\rangle_{\text{FS}} - \left\langle \frac{1}{L_{\psi,s}} \right\rangle_{\text{FS}} \left\langle \frac{B^2}{N_s} \right\rangle_{\text{FS}}^{-1} \right) \right\}, \quad (23) $$

$$ \Gamma^{K_s} = +\frac{I}{q_s} M_s \nu_{si} \left\{ \frac{\left\langle \frac{B}{N_s} \frac{\partial}{\partial \psi} \left( \frac{\Pi_{\parallel,s}}{q_s B} \right) \right\rangle_{\text{FS}}}{\left\langle B^2 / N_s \right\rangle_{\text{FS}}} - \frac{1}{B} \frac{\partial}{\partial \psi} \left( \frac{\Pi_{\parallel,s}}{q_s B} \right) \right\}, \quad (24) $$
with $u \simeq -0.33 \frac{I}{e^{\langle B^2 \rangle_{FS}}} \partial_{\psi} T_i$, $L_{\psi,s}^{-1} = -\partial_{\psi} P_{\perp,s}/T_i Z_s N_s$ and $L_{\psi,i}^{-1} = \partial_{\psi} \ln P_i - H_{TS} \partial_{\psi} \ln T_i$. The thermal screening factor, $H_{TS}$, is defined here as $H_{TS} = 3/2$, a value valid as long as the main ion species remains in the Banana-Plateau regime [39].

The first two terms, $\Gamma^{\nabla N_i}$ and $\Gamma^{u_s}$, represent convection contributions that are proportional to the ratio $M_s \nu_{ei}/q_s$, making them highly significant for heavy impurities that are not fully ionized and are in a high-collisionality regime, such as tungsten. While $\Gamma^{u_s}$ is particularly sensitive to poloidal asymmetries in density, $\Gamma^{\nabla N_i}$ is more critical as it is proportional to both the density and temperature logarithmic gradients (i.e., gradient lengths $\kappa_{\{N,T\}} = -R \partial_r \log\{N, T\}$) of the main ion species. Specifically, this term is linked to the thermal screening factor $H_{TS}$, which is expected to be important in high-confinement regimes, where the temperature gradient length can substantially exceed the density gradient length, potentially enabling impurity screening.

The term $\Gamma^{\nabla N_s}$ represents a standard diffusion process, which is generally minor unless there is significant impurity peaking in the plasma, like at the pedestal during H-mode operations. On the other hand, $\Gamma^{K_s}$ is associated with the frictional force and pressure anisotropy of impurities. This term becomes essential in the presence of pressure anisotropy, which we induce through a vorticity source to create a local $E \times B$ poloidal shear in the subsequent simulations.

The equations above have been shown to accurately describe neoclassical transport observed in non-linear GYSELA simulations, from light to heavy impurities [24]. The term $\partial_{\theta} \Pi_{\|}$ can be significantly influenced by the vorticity source, as discussed previously. Consequently, we expect the neoclassical impurity flux to be predominantly governed by this term in regions near the applied vorticity source.

Although the anisotropy induced by this vorticity source may not be emulating directly a heating source, its effects are still relevant since heating methods, particularly NBI [40] and ICRH [41, 42, 43], are known to enhance parallel and perpendicular pressure respectively [44] and drive an asymmetrical poloidal distribution of impurities, enhancing their transport as a result. Therefore, the mechanisms involved from experiments, especially the heating systems used, are at least partially retrieved through that vorticity source to the exception of the energy deposition itself.
# 3 Parameters and simulations setup

## 3.1 Simulation parameters

<figure><img src="image_3.png" /><figcaption>Figure 1: (a): Normalized radial profiles of the diffusive buffer (solid blue line), fluid energy source (dashed-dotted green line) and fluid vorticity source (dashed orange line). (b): Normalized fluid vorticity source radial profile for $\theta = 0$ (dashed-dotted gray line), $\theta = \pi$ (dashed purple line) and surface-averaged (solid orange line) for $0.65 \le r/a \le 0.75$, the source "influence" zone where its value is non-zero. The red vertical lines indicate the location of the vorticity source.</figcaption></figure>

Simulation parameters are similar to previous transport barrier studies [32] and summarized in Table 1. The set of simulation presented in this paper consists in a base simulation designed to trigger ITG instabilities with $\eta = \kappa_T/\kappa_N = 3$ on the whole simulation domain, which covers a range of normalized radius from $r/a = 0$ to $r/a = 1.15$. A schematic overview of the simulation branches is provided in Figure 3. A diffusive buffer is applied from $r/a > 0.95$ to damp numerical oscillation at the edge (See Figure 1a for the radial profile of the diffusive buffer). Dirichlet boundary conditions on the electrostatic potential are such that $\phi(r/a = 1.15) = 0$, while $E_r = 0$ is applied at the inner radial boundary, close to the magnetic axis. Since we run the simulation in a flux-driven regime, $N$ and $T$ are allowed to evolve on the whole domain, but all non-axisymmetric fluctuations are forced to zero at both radial boundaries. One must note the magnetic axis is encompassed in the simulation domain, even though the point $r/a = 0$ does not exist to avoid divergence. All quantities are periodic along $\theta$ and $\phi$. A parabolic safety factor profile is chosen with the following equation

$$q(r) = 1.5 + 1.3 \exp [2.5 \ln (r/a)] \quad (25)$$

which creates sufficient magnetic shearing $s = \frac{r}{q} \frac{dq}{dr}$ (see Figure 2a for the radial profile of magnetic shearing and safety factor) to stabilize transverse Kelvin-Helmholtz instability modes [45, 46]. No major effects, such as internal transport barrier triggering, is expected since it is monotonic. The radial profile of the fluid isotropic heat source applied on the main ion species is shown in Figure 1a, green curve, and is shared among all simulations
<table><thead><tr><th>Parameters</th><th>Reference</th><th>Vorticity</th></tr></thead><tbody><tr><td><i>D</i><sup>+</sup> collisionnality</td><td><i>ν</i><sup>*</sup><sub>D+</sub> (<i>r</i>/<i>a</i> = 0.575) = 0.1</td><td></td></tr><tr><td>Main species charge / atomic number</td><td><i>Z</i><sub>i</sub> = 1, <i>A</i><sub>i</sub> = 2</td><td></td></tr><tr><td>Time step</td><td>Δ<i>t</i><i>ω</i><sub>c,0</sub> = 20</td><td></td></tr><tr><td><i>N</i><sub>r</sub> × <i>N</i><sub>θ</sub> × <i>N</i><sub>φ</sub> × <i>N</i><sub><i>v</i><sub>∥</sub></sub> × <i>N</i><sub>μ</sub></td><td>511 × 1024 × 64 × 127 × 51</td><td></td></tr><tr><td>Normalized gyroradius</td><td><i>ρ</i><sup>*</sup> = <i>ρ</i><sub>c,0</sub>/<i>a</i> = 1/200</td><td></td></tr><tr><td>Inverse aspect ratio</td><td>1/<i>ε</i> = <i>R</i><sub>0</sub>/<i>a</i> = 4.4</td><td></td></tr><tr><td>Maximum density gradient<br>for the main species <i>D</i><sup>+</sup></td><td><i>κ</i><sub>N</sub> = <i>R</i><sub>0</sub>/<i>L</i><sub>N</sub> = 2.2</td><td></td></tr><tr><td>Maximum temperature gradient<br>for the main species <i>D</i><sup>+</sup></td><td><i>κ</i><sub>T</sub> = <i>R</i><sub>0</sub>/<i>L</i><sub>T</sub> = 6.6</td><td></td></tr><tr><td>Amplitude of the vorticity source</td><td><i>S</i><sub>0</sub><sup>Ω</sup> = 0</td><td><i>S</i><sub>0</sub><sup>Ω</sup> = 0.005</td></tr><tr><td>Amplitude of the heat source</td><td><i>S</i><sub>0</sub><sup>H</sup> = 0.0085</td><td></td></tr></tbody></table>

Table 1: Simulation parameters used for the main simulations studied in this paper. The sources are only applied to the main ion species, namely deuterium $D^+$. The number of points in the radial ($N_r$) and poloidal ($N_\theta$) directions are, until injection of impurities, kept at half resolution of $N_r \times N_\theta = 255 \times 512$ to reduce computational time. The poloidal momentum (vorticity) source is activated from $t\omega_{c,0} = 120000$ for the vorticity case and disabled in the reference case. See Figure 3 for a better view on the simulation branching.

presented in this paper. Its amplitude is constant in time. The value of $S_0^H = 0.0085$ maintains a relatively constant temperature profile and therefore helps reaching a quasi-steady state faster. Deuterium $D^+$ is the main ion species throughout all simulations presented and is in the banana regime since deuterium collisionality is such that $\nu_{D^+}^* < 1$ (See Figure 2b for the radial profile of collisionality for the different species), with the collisionality defined as [7]

$$ \nu_s^* = \frac{qR_0 4\sqrt{\pi} e^4 \log \Lambda Z_D^2}{3 (4\pi\epsilon_0)^2 T_D^2} \left[ N_D Z_D^2 + \sum_{s \neq D} \sqrt{2}N_s Z_s^2 \frac{1 + M_D/M_s}{(1 + v_{T,s}^2/v_{T_D}^2)^{3/2}} \right], \quad (26) $$

where $v_{T,s} = \sqrt{T_s/M_s}$ is the thermal velocity of the considered species.

## 3.2 Transport barrier and impurity injection

Impurities are introduced as an additional kinetic species in the GYSELA model, requiring an extra Vlasov equation to be solved. The impurities studied are either helium $He^{2+}$ ($A_{He} = 4$, $Z_{He} = 2$) or tungsten $W^{40+}$ ($A_{W} = 184$, $Z_{W} = 40$). The radial profile of collisionality, presented in Figure 2b, shows that helium is in the banana regime, similar to deuterium, whereas tungsten is in an intermediate regime between the plateau and Pfirsch-Schlüter regimes, depending on radial position. The initial radial density profiles of these impurities are shown in Figure 4 (solid lines), where two different tungsten density profiles are investigated. The initial temperature radial profiles of impurities are identical to that of deuterium and are presented in Figure 6 (dashed black line).

Similarly to Ref. [32], once saturation of turbulent modes is reached after $t = 60000\omega_{c,0}^{-1}$, we let the simulation run for an additional $80000\omega_{c,0}^{-1}$ to allow the temperature profile to converge to a quasi-steady state. Letting the profiles evolve, especially the temperature,
<figure><img src="image_4.png" /><figcaption>Figure 2: (a): Safety factor *q* (dashed blue line) and magnetic shear *s* (dash-dotted red line) radial profiles. (b): Radial profile of the collisionality for all species investigated in this paper; *D*⁺ (main ion, blue line), *He*²⁺ (impurity, orange line), *W*⁴⁰⁺ with a negative density gradient (impurity, green line, density increases with the minor radius). Note that for the positive gradient case, meaning when density decreases with the minor radius, collisionality of tungsten is identical to the negative gradient case since initial conditions are very close, leading to identical collisionality profiles. Plateau and Pfirsch-Schlüter limits are also represented in dotted and dashed black lines. The red vertical line indicates the vorticity source position while the red area corresponds to the source "influence" zone presented in Figure 1b, detailing the radial profile of the vorticity source in the radial range 0.65 < *r*/*a* < 0.75.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

Figure 3: Architecture of the simulation runs with the different important timestamps. The reason for this branching procedure is that numerical computation time required to evolve two species in global non-linear gyrokinetic simulations is challenging. Each of these branches is then ran for $50000\omega_{c,0}^{-1}$ to reach a final time of $t = 256000\omega_{c,0}^{-1}$. Note that an additional restart of $6000\omega_{c,0}^{-1}$ has been used to transit from $N_r \times N_\theta = 255 \times 512$ to $N_r \times N_\theta = 511 \times 1024$ and let profiles adapt to the refined grid before adding impurities.

allows the turbulence to settle and fluxes to remain on average constant in time, similarly to reference [24]. Once $t = 140000\omega_{c,0}^{-1}$ is reached, the simulation is split in two distinct branches. One, hereafter called the *reference branch*, does not receive any changes and just continues unaltered whereas the vorticity source is activated in the second branch. The latter is referenced to as the *barrier branch*. Both branches are run for an additional $60000\omega_{c,0}^{-1}$ to reach $t = 200000\omega_{c,0}^{-1}$. No major changes occur in the reference branch while a transport barrier is successfully triggered through the $E \times B$ shear injected in the barrier case, hence the name. This barrier is similar in every way to the one presented in reference [32]. Additional focus and details on the transport barrier and its characteristics can be found in this reference.

To summarize, the vorticity source successfully produces a shear flow and as a result lowers significantly both the turbulence intensity and deuterium heat flux at the source location. The reduction in turbulence intensity is also observed in the core and a small increase in core temperature follows as shown in the radial temperature profile in Figure 6. Additionally, both sides of the created barrier show a steepening of the deuterium temperature radial profile. To check whether this behaviour is recovered in the presented simulations, we compute the radial profile of turbulent heat diffusivity for deuterium normalized to $\chi_{GB} = \rho^{*} \frac{T_{e,0}}{q_i B}$, the gyro-Bohm heat diffusivity of the main ion species (i.e. deuterium), with $Q_{\text{turb,s}} = \left\langle \int v_{E_{n \neq 0}}^r \left( \frac{1}{2}mv_{\|}^2 + \mu B \right) \bar{F}_s d^3 v_s \right\rangle_{\text{FS}}$ the turbulent heat flux, $N_{e,0}$ and $T_{e,0}$ respectively the initial electron density and temperature radial profiles. The diffusivity is averaged over a radial region of $\delta r = 0.1 \times r/a$ with a radial step of $0.05 \times r/a$ and over a time period $\Delta t = 10000\omega_{c,0}^{-1}$ (i.e. from 244000 to 254000$\omega_{c,0}^{-1}$). This moving average along the radial direction allows to smooth out the rapid variations of the temperature gradient and heat flux as well as taking into account the radial extension of
<figure><img src="image_6.png" /><figcaption>Figure 4: Initial radial density profile of $D^+$ (solid blue line), $He^{2+}$ (solid orange-red line), $W^{40+}$ with $\kappa_N = 2.2$ (solid black line) and $\kappa_N = -2.2$ (solid light green line). For $D^+$, $t_{init} = 0\omega_{c,0}^{-1}$ while for the other species $t_i = 206000\omega_{c,0}^{-1}$. The average density profile for the last $10000\omega_{c,0}^{-1}$ is also shown for the reference (dashed lines) and barrier (dotted lines) cases, with a corresponding colour palette (i.e. $He^{2+}$ in orange, $W^{40+}$ with $\kappa_N = 2.2$ in gray and with $\kappa_N = -2.2$ in green). The red vertical line indicates the vorticity source position while the red area corresponds to the source "influence" zone.</figcaption></figure>

the turbulent structures. The turbulent diffusivity is then computed with

$$ \chi_{\text{turb,s}}(r) = - \frac{\langle Q_{\text{turb,s}} \rangle_{FS,\delta r, \Delta t}}{\langle N_s \nabla T_s \rangle_{FS,\delta r, \Delta t}}. \qquad (27) $$

The radial profiles of deuterium diffusivity shown in Figure 5 indicate that in the reference case, the maximum diffusivity is reached at $r/a = 0.7$ at a few gyro-Bohm units. However, in the barrier case, the diffusivity drops by one order of magnitude from $r/a \ge 0.4$, and is unaffected for $r/a \le 0.4$. This is in line with what is observed in [32] and confirms the turbulence and heat flux quench is retrieved here. Although it looks subtle on the radial profile of deuterium temperature (Figure 6), the temperature gradient length $\kappa_{T_D} = -R\partial_r \log(\langle T_D \rangle_{FS,\Delta t})$ is found to be almost doubled in the barrier case compared to the reference case, i.e. from 5 to 9. This increase in temperature peaking is
expected to increase with time and converge to a steeper temperature profile. However, satisfying this condition would require to simulate up to a confinement time, which is not accessible yet.

<figure><img src="image_7.png" /><figcaption>Figure 5: Coarse grained turbulent diffusivity of deuterium for the reference (R, solid line) and barrier (B, dashed line) cases. The red vertical line indicates the position of the vorticity source.</figcaption></figure>

At $t = 206000\omega_{c,0}^{-1}$, both the barrier and reference cases are branched in 3 additional simulations in which we add impurities in a trace regime; one with $He^{2+}$ and two with $W^{40+}$ but with different density gradients. The temperature radial profile for all impurities is the same as the initial one of deuterium with $\kappa_T = 6.6$ as shown in Figure 6 (dashed black line). Impurities are added to both branches (i.e. reference and barrier) in a trace regime with $C_{He^{2+}} \equiv N_{He}/N_D = 2 \times 10^{-5}$ and $C_{W^{40+}} \equiv N_W/N_D = 10^{-6}$, respectively in Banana and intermediate Plateau/Pfirsch-Schlüter collisionality regime (See Figure 2b for the radial profile of collisionality) with similar initial density radial profiles as $D^+$ (See Figure 4). Two density gradients, for which $\kappa_N (r/a = 0.5) = \pm 2.2$ (See Figure 4, green and gray curves), are investigated for tungsten as it is physically expected to come from the outside to contaminate the core afterwards. For the sake of comparison and to study the behaviour of transport relative to density gradient sign and extract particle transport coefficients, we chose to investigate both values with $\kappa_N = \pm 2.2$. Assuming that the turbulent particle flux of impurities can be expressed as

$$ \Gamma_{\text{turb,s}} = -D\nabla N_s + \Gamma_{\text{convection,s}}, \quad (28) $$

both the diffusion coefficient $D$ and the convection component $\Gamma_{\text{convection,s}}$ can be determined at a given radial point. This requires at least two measurements of the turbulent impurity flux $\Gamma_{\text{turb,s}}$, each corresponding to a different density gradient $\nabla N_s$. This will be detailed in section 4.2.3.
# 4 Simulation results

## 4.1 Density and temperature profiles with a transport barrier

The presence of a transport barrier, triggered by $E \times B$ shear, is found to effectively reduce both particle and heat diffusivities [32], thereby affecting the radial density and temperature profiles of both main ions and impurities. Impurity transport, typically driven by the combined effects of multiple factors, is mainly determined by the radial density and temperature profiles [7, 47]. In this section, we present the radial density (Figure 4) and temperature (Figure 6) profiles of both species in the presence of a transport barrier to better understand how it affects the global behavior of both light and heavy impurities.

<figure><img src="image_8.png" /><figcaption>Figure 6: Flux-surface averaged temperature profiles of deuterium (blue lines), helium (orange lines) and tungsten (green lines) for the reference (solid lines) and barrier cases (dashed lines). The black dashed line represents the initial temperature profiles used for all species at their respective injection time.</figcaption></figure>

In the reference case, the helium radial density profile (Figure 4, middle plot, dashed orange line) is hollowing in the core region ($r/a \le 0.6$) and peaks around the mid-radius ($0.6 \le r/a \le 0.8$), suggesting a net outward flux. Activating the vorticity source, however, results in strong peaking within the core ($0.4 \le r/a \le 0.6$) and strong hollowing at the location of the transport barrier ($r/a = 0.7$). A similar pattern is observed for tungsten, which accumulates in the core for both density gradient cases in the absence of a barrier but shows a pronounced hollowing at the barrier location ($r/a = 0.7$) when the barrier is present.

The helium temperature radial profile (Figure 6) is less peaked than deuterium in both reference and barrier cases within the core region ($r/a \le 0.6$), though a slightly stronger peaking is observed in the barrier case than in the reference case. Consequently, the core temperature for helium is, on average, slightly higher in the barrier case, a trend similar to that of deuterium.

In contrast, tungsten exhibits a higher temperature than deuterium over the entire radial range in both reference and barrier cases. In the barrier case, tungsten exhibits a
strong peaking on the inner side of the transport barrier ($0.5 \le r/a \le 0.75$) but lower temperatures outside the barrier region ($r/a \le 0.75$) compared to the reference case.

To better understand the changes in density and temperature profiles, we compute relative variations from the initial profile:

$$ \delta\langle N_s \rangle_{FS,\Delta t} = \frac{\langle N_s \rangle_{FS,\Delta t} - \langle N_s \rangle_{FS,t_0}}{\langle N_s \rangle_{FS,t_0}} \quad \text{and} \quad \delta\langle T_s \rangle_{FS,\Delta t} = \frac{\langle T_s \rangle_{FS,\Delta t} - \langle T_s \rangle_{FS,t_0}}{\langle T_s \rangle_{FS,t_0}}, \quad (29) $$

where the averaged quantities are computed over a time window of $\Delta t = 10000\omega_{c,0}^{-1}$, and other quantities are measured at the time of impurity injection.

Figure 7 illustrates the radial profile of normalized variations in impurity density and temperature across various cases, including different density gradients and the presence of transport barrier. For helium, in the absence of the transport barrier (orange solid line), the particle flux is directed outward, resulting in a hollow profile within the core region ($r/a < 0.6$), and a peaked profile beyond $r/a = 0.6$, consistent with the radial density profiles shown in Figure 4. However, when the transport barrier is triggered, helium accumulates near the barrier forming a density well at $r/a = 0.7$, as observed in the density variations (Figure 7a) and density radial profiles (Figure 4). Overall, the transport barrier causes accumulation within the barrier region, preventing helium from being flushed out of the core region.

<figure><img src="image_9.png" /><figcaption>Figure 7: Density (a) and temperature (b) variations for $He^{2+}$ (orange), $W^{40+}$ with $\kappa_N = 2.2$ (gray) and $\kappa_N = -2.2$ (green) for the reference (solid lines) and barrier (dashed lines) branches. The temperature variations of tungsten are identical regardless of the sign of $\kappa_N$. Additionally, when the transport barrier is established, the density variations are identical between both gradient cases as well.</figcaption></figure>

A similar pattern, albeit less pronounced, is observed in the radial profiles of temperature variations (Figure 7b), where helium temperature variations in the barrier case are slightly higher within $r/a \le 0.6$ than in the reference case. A slightly lower temperature variation appears outside the barrier in the barrier case, potentially linked to reduced heat flux and diffusivity of the main species. As previously discussed, suppression of turbulence results in a lower heat flux, leading to a steeper temperature gradient near the
barrier (Figure 6). This trend observed for deuterium also extends to helium, given that both species primarily rely on turbulent transport.

For tungsten, accumulation in the core is more pronounced in the reference case when $\kappa_N = -2.2$ compared to $\kappa_N = +2.2$, while in the presence of a barrier, radial profile of density variations become nearly identical between the two gradient cases (Figure 7a, green and gray libes). This behaviour may come from the natural tendency of tungsten to accumulate in the core due to convective transport, while reversing the gradient can amplify inward transport via diffusion, potentially creating a synergy between the two mechanisms. Notably, the sign of $\kappa_N$ becomes irrelevant when the barrier is active, leading to similar variations and behaviour to that of helium. Tungsten accumulates on both sides of the barrier, with a pronounced depletion at $r/a = 0.7$ (Figure 4), deeper than for helium as shown on the radial profile of density.

While the sign of the density gradient does not influence tungsten temperature evolution, distinct differences emerge between the reference and barrier cases. In the reference case, tungsten temperature increases more at the edge than in the core, implying an outward-directed heat flux. In the barrier case, the temperature increases by about 40% on the inner side of the barrier ($0.5 \le r/a \le 0.7$), while a lower temperature increase is observed outside the barrier ($r/a > 0.7$) compared to the reference case (Figure 7b). This coincides with a localized reduction in the radial heat flux (Figures 16c, d, and 17c, d), resulting in a steeper temperature profile as thermal energy accumulates within the barrier as shown on the temperature radial profile (Figure 6).

Heat flux colour maps as a function of time and radius are given in Figures 12, 16 and 17.

## 4.2 Impurity particle flux

In GYSELA, one can separate the total flux-surface averaged particle flux in two separate channels as

$$ \Gamma_{\text{tot}}^{\text{GYS}} = \Gamma_{\text{turb}}^{\text{GYS}} + \Gamma_{\text{neo}}^{\text{GYS}}, \qquad (30) $$

where $\Gamma_{\text{neo}}^{\text{GYS}}$ and $\Gamma_{\text{turb}}^{\text{GYS}}$ are the flux-surface averaged neoclassical and turbulent particle fluxes respectively. These two channels can then be written as

$$ \Gamma_{\text{neo}}^{\text{GYS}} = \left\langle \int (v_D^r + v_{E_{n=0}}^r) \bar{F}_s dv \right\rangle_{\text{FS}}, \qquad (31) $$

$$ \Gamma_{\text{turb}}^{\text{GYS}} = \left\langle \int v_{E_{n \neq 0}}^r \bar{F}_s dv \right\rangle_{\text{FS}}. \qquad (32) $$

The neoclassical channel is the sum of the curvature and gradient drift contributions, $v_D^r = \bar{v}_D \cdot \nabla r$, as well as the toroidally axisymmetric $E \times B$ drift contribution $v_{E_{n=0}}^r = \langle \bar{v}_{E \times B} \cdot \nabla r \rangle_{\phi}$, where the n subscript represents the toroidal mode number, n = 0 referring to the toroidally axisymmetrical mode and $\langle \dots \rangle_{\phi} = \int \dots d\phi/L_{\phi}$ the toroidal average with $L_{\phi} = \int d\phi$. Conversely, the turbulent channel is the non axisymmetric $E \times B$ drift contribution $v_{E_{n \neq 0}}^r = \bar{v}_{E \times B} \cdot \nabla r - v_{E_{n=0}}^r$.

The neoclassical channel, defined in equation (31), is usually dominant for heavy impurities and is strongly affected by poloidal asymmetries as shown by previous calculations in section 2.2. Since the vorticity source is poloidally asymmetric and injects
pressure anisotropy, it is expected to play an important role in the neoclassical particle flux profiles, especially for the Banana-Plateau component in equation (17), as evidenced in equation (20).

The turbulent impurity flux, defined in equation (32), can be significantly affected by the transport barrier, as highlighted in previous work [32]. Since light and low-Z impurities are mainly transported by turbulence, and given that the transport barrier can efficiently reduce turbulence levels by an order of magnitude, helium is expected to be strongly affected by this turbulence quench more than tungsten.

In this subsection, particle fluxes are normalized by $\langle N_s \rangle_{FS} / R$ to make the comparison easier between the different species.

### 4.2.1 Helium

For helium (Figure 8), the radial profile of total particle flux switches from outward in the reference case to inward in the barrier case particularly in the inner side of the barrier. Helium on the inner side of the transport barrier (i.e. $r/a \le 0.7$) accumulates toward the core while a thin layer of outward radial flux is present on the outer side of the vorticity source (i.e. $r/a \ge 0.75$). This shows that helium particles cannot escape the plasma core because of the transport barrier just as hinted by the density variation results.

<figure><img src="image_10.png" /><figcaption>Figure 8: Flux-surface averaged total helium particle radial flux for the reference (solid line) and barrier (dashed line) cases averaged over the last $10000\omega_{c,0}^{-1}$ of the simulation branch. The red vertical line locates the vorticity source while the red area represents its influence zone. All the fluxes are normalized to $\langle N_s \rangle_{FS} / R$.</figcaption></figure>

The radial profile of turbulent flux of helium (Figure 9) is, in the reference case, mainly directed outward (positive), and dominates over neoclassical transport. This trend changes completely in the presence of the transport barrier with a quench of the turbulent flux amplitude from $r/a = 0.3$ to $r/a = 0.7$. A sign reversal is even observed between $r/a = 0.7$ to $r/a = 0.8$, in the region where a sign reversal of the radial gradient of density profile is observed in Figure 4. This is understood as the consequence of standard diffusion but directed inward due to this sign reversal.

The radial profile of neoclassical flux (Figure 9) in the reference case exhibits a small (i.e. relative to the turbulent flux) amount of inward transport due to the low collisionality
of helium. However, significant changes occur in the vicinity of the transport barrier (i.e. $r/a = 0.7$) when the barrier is present with a strong dipolar particle flux profile. On both inner and outer sides of the transport barrier, particle fluxes show similar magnitudes to the turbulent flux in the reference case but the inner (outer) side driving an inward (outward) flux is observed, namely with opposite signs.

<figure><img src="image_11.png" /><figcaption>Figure 9: Flux-surface averaged turbulent (dark red lines) and neoclassical (blue lines) helium particle radial flux for the reference (solid lines) and barrier (dashed lines) cases are averaged over the last $10000\omega_{c,0}^{-1}$ of the simulation branch. The red vertical line represents the vorticity source location while the red area represents its influence zone. All the fluxes are normalized to $\langle N_s \rangle_{FS}/R$.</figcaption></figure>

Estimates of the radial profile of neoclassical flux (Figure 10), computed using GY-SELA profiles ($n, P, T$ etc...) with equations (17) and (18), are found to accurately predict neoclassical transport for both the reference and barrier cases by capturing the main characteristics of the profile as well as its amplitude.

As expected for the reference case, the Banana-Plateau (BP hereafter) radial flux contribution dominates over the Pfirsch-Schlüter (PS hereafter) one (Figure 10b) in the reference case due to the low collisionnality of helium. The same observation can be drawn upon the barrier case, with the dipolar structure characteristics arising from the BP flux as well. A more in-depth analysis of the different components of the BP flux reveals the importance of the term $\Gamma^{\Pi_{\parallel,s}}$ (see equation (20)), which accounts for the effect of the poloidal asymmetry of pressure anisotropy $\partial_{\theta}(\Pi_{\parallel,s}) = \partial_{\theta}[(P_{\parallel,s} - P_{\perp,s})/B]$. This is a direct effect of the vorticity source itself as it is poloidally asymmetrical (Figure 1b, representing the radial profile of the fluid vorticity source for different poloidal angle) and drives pressure anisotropy (Equation (6)). Poloidal maps of the pressure anisotropy of impurities $\Pi_{s,\parallel} = P_{s,\parallel} - P_{s,\perp}$ are given in Figure 11 for the reference and barrier cases. Poloidal asymmetries are strongly enhanced in the transport barrier region when the poloidal momentum source is activated, hence the increase observed in $\Gamma^{\Pi_{\parallel}}$.
<figure><img src="image_12.png" /><figcaption>Figure 10: (a): Flux-surface averaged neoclassical helium particle radial flux for the reference (solid lines) and barrier (dashed lines) cases. Light blue lines represent flux computed in GYSELA with equation (31) while the dark blue lines represent the theoretical flux computed using equations (17) and (18) with the GYSELA profiles of density, pressure, etc. (b): Flux-surface averaged BP (i.e. Banana-Plateau, black lines) and PS (i.e. Pfirsch-Schlüter, green lines) components of the theoretical neoclassical fluxes of helium for the reference (solid lines) and barrier (dashed) cases. All the fluxes are averaged over the last $10000\omega_{c,0}^{-1}$ of the simulation branch. The red vertical line represents the vorticity source location while the red area represents its influence zone. All the fluxes are normalized to $R/N$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Pressure anisotropy $\Pi_{s,\parallel} = P_{s,\parallel} - P_{s,\perp}$ poloidal maps for helium (left column) and tungsten for $\kappa_N = -2.2$ (center column) and $\kappa_N = +2.2$ (right column) without (top row) and with (bottom row) the transport barrier. The dotted red circle represents the vorticity source position at $r = 0.7$. Maps are normalized to the maximum absolute value of the reference poloidal map of the impurity considered at the last simulation time (i.e. top row).</figcaption></figure>

The barrier effect observed on helium particle transport is then due to the way the source acts on the distribution function rather than the barrier itself. The temperature variations observed for helium indicates that the barrier has a positive impact on heat losses with a lower helium radial heat flux, just like for deuterium.
<figure><img src="image_14.png" /><figcaption>Figure 12: Total particle flux (color map) as a function of time (y axis) and radius (x axis) for helium, without (a) and with (b) the transport barrier. Turbulent heat flux (color map) as a function of time (y axis) and radius (x axis) for helium, without (c) and with (d) the transport barrier. The black vertical line indicates the transport barrier position.</figcaption></figure>

The time evolution of total radial impurity flux and turbulent flux for helium are illustrated in Figure 12. One can immediately see the effect of transport barrier on the core transport of both particles and heat, which is significantly quenched when the barrier is present. Additionally, a dipolar structure of particle transport is created in the vicinity of the source, with the outer side of the source being positive (i.e. particles are expelled) and the inner side being negative (i.e. particles are transported in) just like presented in the radial profile of radial particle flux shown in Figure 8.

The heat flux follows the same trend, with a strong reduction in heat flux amplitude both in the core and near the source. It is similar to the deuterium diffusivity radial profile shown in Figure 5, which is not surprising as this transport, mainly turbulent, is driven by the background turbulence due to deuterium ITG modes.

### 4.2.2 Tungsten

The total radial profile of radial flux of tungsten (Figure 13) is, as expected in the reference case, mainly inward regardless of the gradient sign, thus leading to accumulation of tungsten in the core. Interestingly, the case with $\kappa_N = -2.2$ shows higher levels of accumulation than with $\kappa_N = 2.2$ on the whole radial domain. Both values of $\kappa_N$ give however the same transport profiles in the barrier case, with a strong dipolar structure around the source location (i.e. $r/a = 0.7$) similarly to the helium case. Knowing that overall turbulence level is quenched in the barrier cases, this hints at a possible inward turbulent diffusion in the reference case with $\kappa_N = -2.2$ compared to $\kappa_N = 2.2$ due to the sign reversal of the density gradient. The barrier creates as a result a strong inward (respectively outward) transport on the inner (outer) side of the vorticity source, keeping tungsten on the outside of the vorticity source and prevent contamination of the core
provided the source of tungsten remains outside the barrier.

<figure><img src="image_15.png" /><figcaption>Figure 13: Flux-surface averaged total tungsten particle radial flux for the reference (solid lines) and barrier (dashed lines) cases for $\kappa_N = +2.2$ (gray lines) and $\kappa_N = -2.2$ (green lines). All the fluxes are averaged over the last $10000\omega_{c,0}^{-1}$ of the simulation branch. The red vertical line represents the vorticity source location while the red area represents its influence zone. All the fluxes are normalized to $\langle N_s \rangle_{FS} / R$.</figcaption></figure>

This is confirmed when comparing the turbulent radial profile of radial flux profile (Figure 14) for $\kappa_N = -2.2$ to $\kappa_N = 2.2$ in the reference case. For $\kappa_N = 2.2$, the turbulent flux seems to be zero on average whereas for $\kappa_N = -2.2$, the turbulent flux is even higher than the neoclassical flux for $r/a > 0.5$. This can be understood as follows: when $\kappa_N = 2.2$, both the diffusion (outward because $\kappa_N > 0$) and convection (inward for tungsten) terms compensate each other whereas for $\kappa_N = -2.2$, both terms are adding up (diffusion is inward when $\kappa_N = -2.2$). In the barrier case, the turbulent flux is on average zero for both gradient values. For the sake of clarity and to avoid redundancy, the turbulent flux of $\kappa_N = 2.2$ in the barrier case as well as the neoclassical fluxes for $\kappa_N = 2.2$ are not shown. The reason is that they are nearly identical to their negative gradients counterpart, the turbulent flux in the reference case for $\kappa_N = 2.2$ being the only noticeable outlier. Hereafter, we will focus on the case with $\kappa_N = -2.2$, since it does not alter the observation made when the transport barrier is present and is more relevant for experimental scenarios.

The radial profile of radial neoclassical flux (Figure 14) is mainly inward (negative) and enhances the accumulation of tungsten in the core. However, it is deeply altered by the presence of the barrier and exhibits a similar dipolar flux profile as the helium one at $r/a = 0.7$. This profile prevents efficiently tungsten from contaminating the core and acts as a proper transport barrier for particles.

Theoretical radial profiles of radial fluxes (Figure 15a) are, for both cases, matching quantitatively the GYSELA fluxes. Looking at the different components of the radial profiles of theoretical radial fluxes (Figure 15b) gives a better understanding of what drives the transport barrier for tungsten. Interestingly, the shearing observed in the neoclassical flux still partly comes from the BP components and the associated $\Gamma^{\Pi_{||},s}$ term, as a direct result of the vorticity source behaviour. However, the main driver for this transport barrier comes from the PS component and more specifically from the main ion species
<figure><img src="image_16.png" /><figcaption>Figure 14: Flux-surface averaged turbulent (dark red and orange lines) and neoclassical (blue lines) tungsten particle radial flux for the reference (solid lines) and barrier (dashed lines) cases. Only the negative gradient case is shown for clarity with addition of the turbulent flux for the positive gradient reference case (orange solid line). All the fluxes are averaged over the last $10000\omega_{c,0}^{-1}$ of the simulation branch. The red vertical line represents the vorticity source location while the red area represents its influence zone. All the fluxes are normalized to $\langle N_s \rangle_{FS}/R$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: (a): Flux-surface averaged neoclassical helium particle radial flux for the $\kappa_N = -2.2$ reference (solid lines) and barrier (dashed lines) cases. Light blue lines represent the flux computed in GYSELA with equation (31) while the dark blue lines represent the theoretical flux computed using equations (17) and (18) with the GYSELA profiles of density, pressure, etc. (b): Flux-surface averaged BP (i.e. Banana-Plateau, black lines) and PS (i.e. Pfirsch-Schlüter, green lines) components of the theoretical neoclassical fluxes of tungsten for the $\kappa_N = -2.2$ reference (solid lines) and barrier (dashed) cases. All the fluxes are averaged over the last $10000\omega_{c,0}^{-1}$ of the simulation branch. The red vertical line represents the vorticity source location while the red area represents its influence zone. All the fluxes are normalized to $\langle N_s \rangle_{FS}/R$.</figcaption></figure>
pressure gradient term $\Gamma^{\nabla N_i}$ (see equation (21)), with $1/L_{\psi,i} = \partial_{\psi} \ln P_i - \frac{3}{2} \partial_{\psi} \ln T_i$. The temperature profile steepening of the main species is then responsible for preventing tungsten to cross $r/a = 0.7$ as a direct result of the transport barrier effect on energy transport rather than the enhanced pressure anisotropy. However, the term $\Gamma^{K_s}$ is unexpectedly low compared to what one could expect from equation (24) as shown in Figure 15. This term, which is proportional to the radial gradient of pressure anisotropy, is moderated by the fact that this anisotropy is in fact spread out poloidally as shown in the 2D poloidal maps of $\Pi_{\parallel,s}$ in Figure 11, in both the reference and barrier cases.

<figure><img src="image_18.png" /><figcaption>Figure 16: Total particle flux (color map) as a function of time (y axis) and radius (x axis) for tungsten with $\kappa_N = +2.2$, without (a) and with (b) the transport barrier. (c) and (d) plots represent the turbulent heat flux of tungsten with $\kappa_N = +2.2$ as a function of time (y axis) and radius (x axis). The black vertical line indicates the transport barrier position.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 17: Total particle flux (color map) as a function of time (y axis) and radius (x axis) for tungsten with $\kappa_N = -2.2$, without (a) and with (b) the transport barrier. (c) and (d) plots represent the turbulent heat flux of tungsten with $\kappa_N = -2.2$ as a function of time (y axis) and radius (x axis). The black vertical line indicates the transport barrier position.</figcaption></figure>

The time evolution of the radial profiles of radial total and turbulent impurity fluxes for tungsten are illustrated in Figures 16 and 17 for $\kappa_N = +2.2$ and $\kappa_N = -2.2$ respectively. Both gradient values lead to similar results as already stated, but also one can note that both particle and turbulent heat fluxes can be positive in the positive gradient case as opposite to the negative one. The total particle flux is then lowered by the transport barrier on the whole radial domain. The dipolar profiles in the vicinity of the source, which we retrieve from the radial profiles of the radial particle fluxes presented in Figure 13, are established quickly after the impurity injection time and increase in amplitude with increasing time.

The heat flux is very similar in behaviour to the total particle flux, with again a strong reduction in heat flux amplitude both in the core and near the source, similarly to the deuterium diffusivity radial profile shown in Figure 5.

### 4.2.3 Tungsten turbulent transport coefficient

The turbulent transport of tungsten becomes comparable to the neoclassical flux, as shown by the radial profiles of radial particle flux (Figure 14) with $\kappa_N = -2.2$, due to the combined effects of diffusion and convection. The turbulent flux in equation (28) can be decomposed into several components [9] in the trace impurity limit as follows:

$$ \frac{R\Gamma_{\text{turb,s}}}{N_s} = D\kappa_{N_s} + D_T\kappa_{T_s} + D_u u'_s + RV_s, \qquad (33) $$

where the terms correspond to turbulent diffusion, thermo-diffusion, roto-diffusion, and pure convection, respectively, with $u'_s = -R^2\partial_r\Omega_{\text{tor}}/v_{\text{th,s}}$ denoting the normalized toroidal
rotation gradient. While the diffusion coefficient *D* is always positive, the sign of *D*<sub>*T*</sub>, *D*<sub>*u*</sub>, and *V*<sub>*s*</sub> vary depending on the nature of the background turbulence (whether ITG or TEM) and the sign of magnetic shear. For instance, thermo-diffusion and roto-diffusion are directed inward in case of TEM turbulence, while outward for ITG turbulence. In addition, positive magnetic shear leads to inward convection, whereas negative magnetic shear results in outward convection [9]. These impurity pinch terms in gyrokinetic modeling have been investigated using the bounce-averaged gyrokinetic code TERESA, yielding qualitatively consistent results [48]. Consequently, the combined effects of these terms can lead to either peaked or hollow impurity density profiles.

We compare the two different density gradients for tungsten, $\kappa_{N,W} = \pm 2.2$, with and without a transport barrier. This approach enables us to distinguish the contribution of the pure diffusion, as thermo-diffusion and roto-diffusion remain the same in both cases since both simulations for tungsten are identical in background turbulence (i.e. *D* is therefore the same), temperature profiles (i.e. identical thermo-diffusion term) as well as parallel velocity profiles (i.e. identical roto-diffusion term). Only the tungsten density profiles (i.e. $N_s$ and subsequently $\kappa_{N_s}$) are different with similar pure convection terms. By focusing on a small radial region in both simulations (i.e. reference cases, $\kappa_{N,s} = 2.2$ and $-2.2$) similarly as for the radial profile of heat diffusivity (See Figure 5), one can expect to extract both diffusion *D* and full convection<fnref n="3" /> $RV_s^F \equiv D_T\kappa_{T_s} + D_uu_s' + RV_s$ terms by performing local linear regressions on $\frac{\Gamma^{turb}}{N_s}$ ($\kappa_{N_s}$). We then extract the slope, *D*, and the intercept, $RV_s^F$, for the considered region. This approach is relevant locally as long as the profiles do not undergo major changes on the radial and time intervals considered.

<figure><img src="image_20.png" /><figcaption>Figure 18: Coarse grained radial profiles of diffusion normalized to the local average gradient length $\langle|\kappa_{N,s}|\rangle_{\delta r} D_W$ (a) and full convection $RV_W^F$ (b) coefficients of tungsten for the reference (blue lines) and barrier (red lines) cases. The coloured areas around the lines represent the error bars. Values are normalized to the local gyro-Bohm diffusivity $\chi_{GB}$.</figcaption></figure>

Figure 18 illustrates the radial profile of diffusion (Figure 18a) and full convection

<footnote n="3">³Here, the superscript *F* refers to the full convection being the sum of the thermo-diffusion, roto-diffusion and pure convection terms.</footnote>
(Figure 18b), both with and without the transport barrier. The diffusion coefficient of tungsten, primarily determined by ITG turbulence in the trace impurity limit, remains positive across the entire minor radius, while convection is directed inward, leading to tungsten accumulation in the core region. The magnitude of convection exceeds that of diffusion in both cases, indicating that the tungsten turbulent flux is mainly driven by the convective term. Once the barrier is activated, both diffusion (D) and convection (RV<sup>F</sup>) are significantly quenched, particularly around the region at r/a ≤ 0.7, but also where the barrier is locally induced. The effects of the barrier are particularly significant on diffusion, reducing the coefficient D by 5 up to an order of magnitude, while convection decreases by almost a factor of three. These results show that the presence of a transport barrier, and more specifically a strong E × B shear, provides an effective means for preventing tungsten accumulation in the core region.

Interestingly, the convection is similar between the reference and barrier cases in the core up to r/a ≈ 0.4 and shows a small sign reversal at r/a = 0.6, a position where density tends to flatten. This indicates that in this area, the tungsten density profile is becoming hollow.

Convection is, in both cases, dominating over diffusion for tungsten. All the turbulent transport coefficients are then strongly reduced on the whole radial domain with a stronger effect in the vicinity of the transport barrier when the latter is present.

# 5 Conclusion and discussion

The effects of transport barrier on both turbulent and neoclassical impurity transport are investigated using global, flux-driven, 5D gyrokinetic code GYSELA simulations. Transport barrier is generated in nonlinear simulations via external poloidal momentum injection, which triggers strong $E \times B$ shear. The induced barrier effectively suppresses the background ITG turbulence, globally reducing particle and heat diffusivity.

The introduction of the $E \times B$ shear through the vorticity source significantly alters both density and temperature profiles of impurities. For helium, the vorticity source induces strong hollowing at the transport barrier location while increasing density peaking in adjacent regions, effectively confining helium within the barrier region and limiting outward flux. Tungsten similarly exhibits a pronounced hollowing at the barrier location, with accumulation on both sides as well. The vorticity source also increases core temperature for helium and creates a sharp temperature gradient for tungsten within the barrier. This modification reflects the reduced turbulent diffusivity of the main ions, enhancing energy confinement and leading to a localized thermal energy accumulation.

The turbulence suppression by the $E \times B$ shear results in a steep ion temperature gradient, which increases the thermal screening term in the neoclassical particle flux. This enhanced thermal screening prevents tungsten from accumulating in the core region. However, the pressure anisotropy, represented by the $\Gamma^{\Pi}_{||}$ term in equation (20), driven by the external vorticity source, often complicates the interpretation of the results. This phenomenon is more pronounced in the case of helium due to its lighter mass, which is not directly observed in experiments with edge transport barriers. Conversely, for heavy impurities like tungsten, the increasing significance of pressure terms, more precisely $\Gamma^{\nabla N_i}$ in equation (21), is more favourable, as it is expected to play an important role in experiments[23] to enhance confinement and reduce significantly tungsten accumulation as well as other heavy impurities.

Considering that the source of tungsten is located at the edge (i.e. walls, divertor),
this may show that if an intense enough $E \times B$ poloidal shear is maintained throughout
the experiment, heavy impurities will not contaminate the core. However, concerns may
be expressed regarding helium since it needs to be flushed out towards the separatrix and
divertor region. Those results suggest that strong $E \times B$ poloidal shearing may conflict
with this specific requirement.

By considering how the transport barrier affects turbulence in the system, particularly
the turbulent intensity, it is possible to understand the processes occurring in the core re-
gion. Reduced turbulent intensity leads to diminished turbulent impurity and heat fluxes,
as both are directly proportional to turbulence intensity. With the reduction of the latter
in the presence of the transport barrier, as previously evidenced [32], the consequent effect
can be easily evaluated. The reduced turbulent diffusivity of deuterium is also a strong
indicator of reduced turbulent activity.

More precisely, the variations of turbulent transport near the transport barrier when
the latter is present still need to be addressed. The main hypothesis to explain the ob-
served behaviour when the transport barrier is active is linked to the density gradient $\kappa_N$.
As presented in Figure 9, the peaks of turbulent particle flux near the transport barrier,
both negative and positive, coincide with the density accumulation and depletion peaks
observed in the radial profile of density variations (Figure 7a). This indicates that the
observed reversal in turbulent flux may be linked to turbulent diffusion. The fact that
both the density gradient and turbulent particle flux reverse their signs when the trans-
port barrier is present provides evidence for this hypothesis. A similar observation can be
made for the slight increase in turbulent flux on the inner side as it aligns with enhanced
turbulent diffusion via a steeper density gradient. However, the reversal of the gradient
sign in this region should be approached with caution. This is due to the influence of neo-
classical transport which, being heavily driven by $\Gamma^{\Pi}_{\parallel,s}$ for helium, modifies the density
profile near the source. The inward turbulent diffusion is then likely an artefact caused
by the vorticity source rather than an effect of the transport barrier itself.

To extend on turbulent transport of tungsten, the effect of the vorticity source is not
limited to its vicinity, as it also reduces the core turbulent transport. More specifically,
thanks to the fact that we have two different density gradient values for tungsten, we have
extracted both the diffusion $D$ and full convection $RV^F$ coefficients for both the reference
and gradient cases. We observe that convection, which is inward in the reference case,
dominates over diffusion and drives additional inward tungsten transport, leading to a
strong tungsten density peaking. The transport barrier leads to a quench of both diffusion
and convection coefficients, strongly reducing overall transport and core accumulation of
tungsten as a result, flattening the density radial profile.

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium,
funded by the European Union via the Euratom Research and Training Programme (Grant
Agreement No. 101052200-EUROfusion). Views and opinions expressed are, however,
those of the author(s) only and do not necessarily reflect those of the European Union or
the European Commission. Neither the European Union nor the European Commission
can be held responsible for them. This work was granted access to the HPC resources of
the EUROfusion High Performance Computer (Marconi-Fusion) under the project FUA36 - TSVV6. The authors would also like to thank Clemente Angioni for the very helpful discussions.

## References

[1] R. Dux, A. Loarte, C. Angioni, D. Coster, E. Fable, and A. Kallenbach, “The interplay of controlling the power exhaust and the tungsten content in ITER,” *Nuclear Materials and Energy*, vol. 12, pp. 28–35, 2017.

[2] A. Kappatou, R. McDermott, C. Angioni, P. Manas, T. Pütterich, R. Dux, E. Viezzer, R. Jaspers, R. Fischer, M. Dunne, et al., “Understanding helium transport: experimental and theoretical investigations of low-Z impurity transport at ASDEX upgrade,” *Nuclear Fusion*, vol. 59, no. 5, p. 056014, 2019.

[3] R. Neu, E. P. Taskforce, J. E. Contributors, A. U. Team, et al., “Preparing the scientific basis for an all metal ITER,” *Plasma Physics and Controlled Fusion*, vol. 53, no. 12, p. 124040, 2011.

[4] G. Federici, C. Bachmann, L. Barucca, C. Baylard, W. Biel, L. Boccaccini, C. Bustreo, S. Ciattaglia, F. Cismondi, V. Corato, et al., “Overview of the DEMO staged design approach in Europe,” *Nuclear fusion*, vol. 59, no. 6, p. 066013, 2019.

[5] P. Rodriguez-Fernandez, A. Creely, M. Greenwald, D. Brunner, S. Ballinger, C. Chrobak, D. Garnier, R. Granetz, Z. Hartwig, N. Howard, et al., “Overview of the SPARC physics basis towards the exploration of burning-plasma regimes in high-field, compact tokamaks,” *Nuclear Fusion*, vol. 62, no. 4, p. 042003, 2022.

[6] P. Maget, P. Manas, R. Dumont, C. Angioni, J. Artaud, C. Bourdelle, F. Casson, L. Colas, P. Devynck, D. Fajardo, et al., “Tungsten accumulation during ion cyclotron resonance heating operation on WEST,” *Plasma Physics and Controlled Fusion*, vol. 65, no. 12, p. 125009, 2023.

[7] S. Hirshman and D. Sigmar, “Neoclassical transport of impurities in tokamak plasmas,” *Nuclear Fusion*, vol. 21, no. 9, p. 1079, 1981.

[8] C. Angioni and P. Helander, “Neoclassical transport of heavy impurities with poloidally asymmetric density distribution in tokamaks,” *Plasma Physics and Controlled Fusion*, vol. 56, no. 12, p. 124001, 2014.

[9] C. Angioni, “Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments,” *Plasma Physics and Controlled Fusion*, vol. 63, no. 7, p. 073001, 2021.

[10] F. Wagner, G. Becker, K. Behringer, D. Campbell, A. Eberhagen, W. Engelhardt, G. Fussmann, O. Gehre, J. Gernhardt, G. v. Gierke, et al., “Regime of improved confinement and high beta in neutral-beam-heated divertor discharges of the ASDEX tokamak,” *Physical Review Letters*, vol. 49, no. 19, p. 1408, 1982.

[11] T. Kobayashi, “The physics of the mean and oscillating radial electric field in the L-H transition: the driving nature and turbulent transport suppression mechanism,” *Nuclear Fusion*, vol. 60, no. 9, p. 095001, 2020.
[12] C. Bourdelle, “Staged approach towards physics-based L-H transition models,” *Nuclear Fusion*, vol. 60, no. 10, p. 102002, 2020.

[13] U. Plank, R. McDermott, G. Birkenmeier, N. Bonanomi, M. Cavedon, G. Conway, T. Eich, M. Griener, O. Grover, P. Schneider, et al., “Overview of L-to H-mode transition experiments at ASDEX Upgrade,” *Plasma Physics and Controlled Fusion*, vol. 65, no. 1, p. 014001, 2022.

[14] H. Biglari, P. Diamond, and P. Terry, “Influence of sheared poloidal rotation on edge turbulence,” *Physics of Fluids B: Plasma Physics*, vol. 2, no. 1, pp. 1–4, 1990.

[15] K. Itoh and S.-I. Itoh, “The role of the electric field in confinement,” *Plasma physics and controlled fusion*, vol. 38, no. 1, p. 1, 1996.

[16] P. Terry, “Suppression of turbulence and transport by sheared flow,” *Reviews of Modern Physics*, vol. 72, no. 1, p. 109, 2000.

[17] P. Guzdar, L. Chen, W. Tang, and P. Rutherford, “Ion-temperature-gradient instability in toroidal plasmas,” *The Physics of Fluids*, vol. 26, no. 3, pp. 673–677, 1983.

[18] G. Lee and P. Diamond, “Theory of ion-temperature-gradient-driven turbulence in tokamaks,” *The Physics of fluids*, vol. 29, no. 10, pp. 3291–3313, 1986.

[19] H. Chen, N. Hawkes, L. Ingesson, N. Peacock, and M. Haines, “Impurity transport with strong and weak internal thermal barriers in JET optimized shear plasmas,” *Nucl. Fusion*, vol. 41, p. 31, jan 2001.

[20] R. Dux, C. Giroud, K.-D. Zastrow, and J. E. contributors, “Impurity transport in internal transport barrier discharges on JET,” *Nucl. Fusion*, vol. 44, p. 260, jan 2004.

[21] S. Mochinaga, N. Kasuya, A. Fukuyama, and M. Yagi, “Effect of toroidal rotation on impurity transport in tokamak improved confinement,” *Nucl. Fusion*, vol. 64, p. 066002, apr 2024.

[22] D. Fajardo, C. Angioni, P. Maget, and P. Manas, “Analytical model for collisional impurity transport in tokamaks at arbitrary collisionality,” *Plasma Physics and Controlled Fusion*, vol. 64, no. 5, p. 055017, 2022.

[23] A. Field, F. Casson, D. Fajardo, C. Angioni, C. Challis, J. Hobirk, A. Kappatou, H.-T. Kim, E. Lerche, A. Loarte, et al., “Peripheral temperature gradient screening of high-Z impurities in optimised ‘hybrid’ scenario H-mode plasmas in JET-ILW,” *Nuclear Fusion*, vol. 63, no. 1, p. 016028, 2022.

[24] K. Lim, X. Garbet, Y. Sarazin, V. Grandgirard, K. Obrejan, M. Lesur, and E. Gravier, “Gyrokinetic modelling of light to heavy impurity transport in tokamaks,” *Nuclear Fusion*, vol. 61, no. 4, p. 046037, 2021.

[25] P. Donnel, X. Garbet, Y. Sarazin, V. Grandgirard, N. Bouzat, E. Caschera, G. Dif-Pradalier, P. Ghendrih, C. Gillot, G. Latu, et al., “Neoclassical impurity flux in presence of turbulent generated poloidal asymmetries and pressure anisotropy,” *Plasma Physics and Controlled Fusion*, vol. 61, no. 4, p. 044006, 2019.
[26] K. Lim, X. Garbet, Y. Sarazin, E. Gravier, M. Lesur, G. Lo-Cascio, and T. Rouyer, "Self-consistent gyrokinetic modelling of turbulent and neoclassical tungsten transport in toroidally rotating plasmas," *Physics of Plasma*, 2023.

[27] P. Manas, A. Kappatou, C. Angioni, R. McDermott, A. U. Team, et al., "Light impurity transport in tokamaks: on the impact of neutral beam fast ions," *Nuclear Fusion*, vol. 60, no. 5, p. 056005, 2020.

[28] D. Palade, "Peaking and hollowness of low-Z impurity profiles: an interplay between ITG and TEM induced turbulent transport," *Nuclear Fusion*, vol. 63, no. 4, p. 046007, 2023.

[29] E. Fable, C. Angioni, A. U. Team, et al., "Impurity outward particle flux from externally applied torque," *Plasma Physics and Controlled Fusion*, vol. 65, no. 11, p. 115007, 2023.

[30] D. Fajardo, C. Angioni, R. Dux, E. Fable, U. Plank, O. Samoylov, G. Tardini, A. U. Team, et al., "Full-radius integrated modelling of ASDEX Upgrade L-modes including impurity transport and radiation," *Nuclear Fusion*, vol. 64, no. 4, p. 046021, 2024.

[31] A. Strugarek, Y. Sarazin, D. Zarzoso, J. Abiteboul, A. Brun, T. Cartier-Michaud, G. Dif-Pradalier, X. Garbet, P. Ghendrih, V. Grandgirard, et al., "Ion transport barriers triggered by plasma polarization in gyrokinetic simulations," *Plasma Physics and Controlled Fusion*, vol. 55, no. 7, p. 074013, 2013.

[32] G. Lo-Cascio, E. Gravier, T. Réveillé, M. Lesur, Y. Sarazin, X. Garbet, L. Vermare, K. Lim, A. Guillevic, and V. Grandgirard, "Transport barrier in 5D gyrokinetic flux-driven simulations," *Nuclear Fusion*, vol. 62, no. 12, p. 126026, 2022.

[33] V. Grandgirard, J. Abiteboul, J. Bigot, T. Cartier-Michaud, N. Crouseilles, G. Dif-Pradalier, C. Ehrlacher, D. Esteve, X. Garbet, P. Ghendrih, et al., "A 5D gyrokinetic full-f global semi-Lagrangian code for flux-driven ion turbulence simulations," *Computer Physics Communications*, vol. 207, pp. 35-68, 2016.

[34] D. Esteve, Y. Sarazin, X. Garbet, V. Grandgirard, S. Breton, P. Donnel, Y. Asahi, C. Bourdelle, G. Dif-Pradalier, C. Ehrlacher, et al., "Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport," *Nuclear Fusion*, vol. 58, no. 3, p. 036013, 2018.

[35] P. Donnel, X. Garbet, Y. Sarazin, V. Grandgirard, Y. Asahi, N. Bouzat, E. Caschera, G. Dif-Pradalier, C. Ehrlacher, P. Ghendrih, et al., "A multi-species collisional operator for full-f global gyrokinetics codes: Numerical aspects and verification with the GYSELA code," *Computer Physics Communications*, vol. 234, pp. 1-13, 2019.

[36] Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, P. Ghendrih, G. Latu, A. Strugarek, G. Dif-Pradalier, P. Diamond, et al., "Predictions on heat transport and plasma rotation from global gyrokinetic simulations," *Nuclear Fusion*, vol. 51, no. 10, p. 103023, 2011.

[37] J. Abiteboul, X. Garbet, V. Grandgirard, S. Allfrey, P. Ghendrih, G. Latu, Y. Sarazin, and A. Strugarek, "Conservation equations and calculation of mean flows in gyrokinetics," *Physics of Plasmas*, vol. 18, no. 8, p. 082503, 2011.
[38] G. Chew, M. Goldberger, and F. Low, “The Boltzmann equation and the one-fluid hydromagnetic equations in the absence of particle collisions,” *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, vol. 236, no. 1204, pp. 112–118, 1956.

[39] D. Fajardo, C. Angioni, F. J. Casson, A. R. Field, P. Maget, P. Manas, J. Contributors, A. U. Team, et al., “Analytical model for the combined effects of rotation and collisionality on neoclassical impurity transport,” *Plasma Physics and Controlled Fusion*, vol. 65, no. 3, p. 035021, 2023.

[40] T. Odstrcil, T. Pütterich, C. Angioni, R. Bilato, A. Gude, M. Odstrcil, A. U. Team, et al., “The physics of W transport illuminated by recent progress in W density diagnostics at ASDEX Upgrade,” *Plasma Physics and Controlled Fusion*, vol. 60, no. 1, p. 014003, 2017.

[41] M. Reinke, I. Hutchinson, J. Rice, N. Howard, A. Bader, S. Wukitch, Y. Lin, D. Pace, A. Hubbard, J. Hughes, et al., “Poloidal variation of high-Z impurity density due to hydrogen minority ion cyclotron resonance heating on Alcator C-Mod,” *Plasma Physics and Controlled Fusion*, vol. 54, no. 4, p. 045004, 2012.

[42] Y. O. Kazakov, I. Pusztai, T. Fülöp, and T. Johnson, “Poloidal asymmetries due to ion cyclotron resonance heating,” *Plasma Physics and Controlled Fusion*, vol. 54, no. 10, p. 105010, 2012.

[43] R. Bilato, T. Odstrcil, F. Casson, C. Angioni, M. Brambilla, Y. O. Kazakov, and E. Poli, “The impact of the ion-cyclotron-resonance location on the poloidal asymmetries of impurity density in an ICRF-heated rotating plasma,” *Nuclear Fusion*, vol. 57, no. 5, p. 056020, 2017.

[44] Z. Qu, M. Fitzgerald, and M. J. Hole, “Analysing the impact of anisotropy pressure on tokamak equilibria,” *Plasma Physics and Controlled Fusion*, vol. 56, no. 7, p. 075007, 2014.

[45] Y. Kim, P. Diamond, H. Biglari, and P. Terry, “Theory of resistivity-gradient-driven turbulence in a differentially rotating plasma,” *Physics of Fluids B: Plasma Physics*, vol. 2, no. 9, pp. 2143–2150, 1990.

[46] B. D. Scott, P. Terry, and P. Diamond, “Saturation of Kelvin-Helmholtz fluctuations in a sheared magnetic field,” *The Physics of fluids*, vol. 31, no. 6, pp. 1481–1491, 1988.

[47] P. Helander and D. J. Sigmar, *Collisional transport in magnetized plasmas*, vol. 4. Cambridge university press, 2005.

[48] K. Lim, E. Gravier, M. Lesur, X. Garbet, Y. Sarazin, and J. Médina, “Impurity pinch generated by trapped particle driven turbulence,” *Plasma Physics and Controlled Fusion*, vol. 62, no. 9, p. 095018, 2020.