# Self-consistent gyrokinetic modeling of turbulent and neoclassical tungsten transport in toroidally rotating plasmas

K. Lim, X. Garbet, Y. Sarazin, E. Gravier, Maxime Lesur, Guillaume

Lo-Cascio, T. Rouyer

► To cite this version:

K. Lim, X. Garbet, Y. Sarazin, E. Gravier, Maxime Lesur, et al., Self-consistent gyrokinetic modeling of turbulent and neoclassical tungsten transport in toroidally rotating plasmas. Physics of Plasmas, 2023, 30 (8), 10.1063/5.0157428. [hal-04180521](https://hal-univ-lorraine.fr/hal-04180521)

HAL Id: hal-04180521

https://hal.univ-lorraine.fr/hal-04180521

Submitted on 12 Aug 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Self-consistent gyrokinetic modelling of turbulent and neoclassical tungsten transport in toroidally rotating plasmas

K. Lim¹*, X. Garbet²,⁴, Y. Sarazin², E. Gravier³, M. Lesur³, G. Lo-Cascio³ and T. Rouyer³

¹École Polytechnique Fédérale de Lausanne (EPFL),
Swiss Plasma Center (SPC), CH-1015 Lausanne, Switzerland

²CEA, IRFM, F-13108 Saint-Paul-Lèz-Durance, France

³Institut Jean Lamour (IJL), UMR 7198 CNRS
- Université de Lorraine, 54000 Nancy, France

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

## Abstract

The effect of toroidal rotation on both turbulent and neoclassical transport of tungsten (W) in tokamaks is investigated using the flux-driven, global, nonlinear 5D gyrokinetic code GYSELA. Nonlinear simulations are carried out with different levels of momentum injection that drive W into the supersonic regime, while the toroidal velocity of the main ions remains in the subsonic regime. The numerical simulations demonstrate that toroidal rotation induces centrifugal forces that cause W to accumulate in the outboard region, generating an in-out poloidal asymmetry. This asymmetry enhances neoclassical inward convection, which can lead to central accumulation of W in cases of strong plasma rotation. The core accumulation of W is mainly driven by inward neoclassical convection. However, as momentum injection continues, roto-diffusion, proportional to the radial gradient of the toroidal velocity, becomes significant and generates outward turbulent flux in the case of ion temperature gradient (ITG) turbulence. Overall, the numerical results from nonlinear GYSELA simulations are in qualitative agreement with the theoretical predictions for impurity transport.

<footnote>*Electronic address: kyungtak.lim@epfl.ch</footnote>
# I. INTRODUCTION

Tungsten (W), a high-Z element, is considered for use in the first-wall plasma-facing components in future tokamaks due to its high melting point, low erosion rate and good thermal conductivity [1]. However, experiments have shown that the accumulation of W in the core region is frequently observed due to neoclassical inward convection, leading to a reduced plasma performance and energy loss through radiation [2, 3]. In particular, even at low W concentration, radiative loss can cause complete extinction of the burning plasma [4]. Consequently, it is crucial to mitigate the accumulation of W in tokamaks to improve the efficiency and performance of the reactor.

Over the past few decades, based on experimental observations of up-down radiative emission [5, 6], there has been a strong emphasis placed on the non-uniform distribution of impurity density across magnetic flux surfaces [7–12]. While conventional neoclassical theory often assumes a uniform density distribution [13, 14], it is now widely accepted that an accurate prediction of neoclassical impurity transport requires a proper estimation of the poloidal asymmetry of impurity density, as it can enhance or reduce neoclassical impurity transport [10, 15–17].

The toroidal rotation of plasma in tokamaks, whether driven intrinsically [18, 19] or externally by neutral beam injection (NBI) heating, often results in a non-uniform distribution of W density. In particular, the centrifugal force generated by toroidal rotation drives impurities towards the outer wall, yielding a strong 'in-out' poloidal asymmetry that leads to the enhancement of neoclassical transport [10, 20–22]. Although the impact of toroidal rotation on turbulent transport is found to be much smaller than on neoclassical transport [21], the turbulent roto-diffusion generated by a radial gradient in the toroidal rotation [23] and the parallel velocity gradient (PVG) instabilities [24] can still play an important role in the case of a strongly rotating plasma.

The impact of toroidal rotation on both turbulent and neoclassical W transport was investigated in Ref. [21]. The main approach used in the study was to compute turbulent and neoclassical transport separately. Turbulent transport was calculated using the gyrokinetic model in the toroidally co-moving frame, whereas neoclassical transport was evaluated using the drift-kinetic code NEO, but without considering frictional effects at the lowest order [17, 25].
The aim of this paper, as an extension of the previous work on impurity transport from light to heavy impurities [26], is to explore the effect of toroidal rotation on turbulent and neoclassical W transport using the full-F gyrokinetic code GYSELA [27]. By employing the multi-species collisional operator implemented in GYSELA [28], we investigate the self-consistent treatment of both neoclassical and turbulent transport channels. The interplay between these channels can have a significant impact on impurity transport and its accumulation in the core of tokamaks [29, 30]. To this end, we examine the impact of sonic (for the impurities) toroidal rotation on turbulent and neoclassical transport in the presence of poloidal asymmetry enhanced by centrifugal forces.

The remainder of the paper is organized as follows. Firstly, in Section II, we introduce the GYSELA code and the source terms for studying impurity transport in a rotating plasma. Section III discusses the results of the nonlinear GYSELA simulations, with a focus on the poloidal asymmetry generated by the external momentum injection. In Section IV, the effect of sonic rotation on both turbulent and neoclassical W transport is explored, highlighting the inward neoclassical convection due to the poloidal asymmetry and the outward roto-diffusion due to the steepening of the velocity gradient. Finally, the paper concludes with a discussion of the main findings and conclusions in Section V.

## II. NUMERICAL MODEL

The GYSELA code is a first-principles, nonlinear, global and flux-driven 5D gyrokinetic code that evolves the gyro-averaged distribution function, $\bar{F}_s$, for each species $s$, in time and phase space $\mathbf{z} = (\mathbf{x}_G, v_{G\parallel}, \mu_G)$, $\mathbf{x}_G$ being the gyro-center position, $v_{G\parallel}$ the parallel velocity and $\mu = mv_{\perp}^2/2B$ the magnetic moment, based on a semi-Lagrangian method for numerical computations [27, 31]. The nonlinear gyrokinetic Vlasov equation is expressed as follows:

$$ \frac{\partial \bar{F}_s}{\partial t} + \frac{1}{B_{\parallel}^*} \nabla_{\mathbf{z}} \cdot (\dot{\mathbf{z}} B_{\parallel}^* \bar{F}_s) = C(\bar{F}_s) + S_{\text{heat}} + S_{\text{mom}} \quad (1) $$

where we use $B_{\parallel}^* = \mathbf{b} \cdot \mathbf{B}_{\parallel}^*$ with the generalized magnetic field $\mathbf{B}_{\parallel}^* = \mathbf{B} + (m_s v_{G\parallel} / e Z_s) \nabla \times \mathbf{b}$, the charge number $Z_s$, the mass $m_s$, the collision operator $C(\bar{F}_s)$ and $S_{\text{heat}}$, $S_{\text{mom}}$ represent the heat and toroidal momentum sources respectively. The evolution of $\bar{F}_s$ is self-consistently coupled to the quasi-neutrality equation to ensure the balance between the total density $N_{eq}$ and the electric potential $\phi$. In the limit of long wavelengths and assuming adiabatic
electrons, the electrostatic quasi-neutrality equation can be written as:

$$e \left( \frac{\phi - \langle \phi \rangle}{T_e} \right) - \sum_s \frac{1}{N_{eq,s}} \nabla_{\perp} \cdot \left( \frac{m_s N_{eq,s}}{e Z_s B^2} \nabla_{\perp} \phi \right) = \sum_s \frac{1}{N_{eq,s}} \int d\mathbf{v} \mathcal{J}(\bar{F}_s - \bar{F}_{eq,s}), \quad (2)$$

where $\langle \dots \rangle$ denotes the flux surface average, $T_e$ is the electron temperature, $\int d\mathbf{v} = \mathcal{J}_v d\mu dv_{G\parallel}$ with $\mathcal{J}_v = (2\pi B_{\parallel}^*/m_s)$ and $\mathcal{J}$ is the gyro-average operator. In the above expressions, the mass $m$, the magnetic field $B$, the density $N$, the temperature $T$ and the electric potential $\phi$ are normalized with respect to the reference values $m_0, B_0, N_0, T_0$ and $T_0/eZ_0$. The length scales are normalized to the reference Larmor radius $\rho_0 = v_{T_0}/\omega_{c0}$, where $v_{T_0} = \sqrt{T_0/m_0}$ and $\omega_{c0} = eZ_0B_0/m_0$. Time is normalized to $\omega_{c0}^{-1}$. All drift velocities are normalized to $v_{T_0}$, while the toroidal velocity is normalized to $v_{T_0,s} = v_{T_0} \sqrt{m_0/m_s}$. Moreover, in the large aspect ratio limit, we assume that the parallel velocity of the particles in GYSELA is equivalent to the toroidal velocity. Under this assumption, the parallel velocity measured in GYSELA can be interpreted as the Mach number, which we define as $M_s = v_{\varphi,s}/v_{th,s}$.

The presence of the heat source $S_{\text{heat}}$ in Eq. (1) is essential to allow flux-driven simulations to reach a statistical steady state. External torque is injected into the plasma via the momentum source $S_{\text{mom}}$. Each of these sources is designed to have a non-vanishing contribution to a single fluid moment equation, related to heat and parallel momentum in the present study.

In this paper, as derived in Ref. [32], we adopt the following form of heat and momentum sources

$$S_{\text{heat}} = \frac{S_E S_r}{\sqrt{2\pi^{3/2} T_{sce}^{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] e^{-\bar{v}_{G\parallel}^2 - \bar{\mu}}, \quad (3)$$

where we use the normalized parallel velocity $\bar{v}_{G\parallel} = v_{G\parallel}/\sqrt{2T_{sce}/m_s}$, the normalized magnetic moment $\bar{\mu} = \mu B/T_{sce}$, $J_{\parallel B} = \sqrt{2T_{sce}}J_{\parallel}/B^2$ with the normalized source temperature $T_{sce}$, the parallel current $J_{\parallel}$, the radial envelope profile $S_r$ and the amplitude of the source $S_E$. Similarly to the heat source, the analytical expression of the momentum source can be recast as:

$$S_{\text{mom}} = \frac{S_{\text{mom}} S_r}{4\pi^{3/2} T_{sce}^2} [2\bar{v}_{G\parallel}(2 - \bar{\mu}) - J_{\parallel B}(1 + 2\bar{v}_{G\parallel}^2 - \bar{\mu})]e^{-\bar{v}_{G\parallel}^2 - \bar{\mu}}, \quad (4)$$

where $S_{\text{mom}}$ is the amplitude of the momentum source that acts as a control parameter of the simulations to achieve different toroidal rotations in the nonlinear simulations. The
radial profile of $S_r$ is composed of two hyperbolic tangents

$$S_r(\rho) = -\frac{1}{2} \left[ \tanh \left( \frac{\rho - \rho_s - 3L_s}{L_s} \right) + \tanh \left( \frac{\rho_s - 3L_s - \rho}{L_s} \right) \right], \quad (5)$$

where we use $\rho = r/a$, the radial position of the source $\rho_s$ and the width of the source $L_s$. Although off-axis injection of NBI is known to provide control over density and current profiles [33], the present study concentrates on central injection of toroidal momentum to achieve consistent gradient profiles with experimental observations.

In a toroidally rotating frame with the velocity $\mathbf{u}$, the particle velocity $\mathbf{v}$ in the laboratory reference frame is expressed as $\mathbf{v} = \Omega R^2 \nabla \varphi + \mathbf{u}$ with the toroidal rotation frequency $\Omega$ and the Hamiltonian is defined as follows:

$$H = \frac{1}{2}mu_{\parallel}^{2} + \mu B + eZ\phi - \frac{m\Omega^{2}R^{2}}{2}, \quad (6)$$

where the electrostatic potential is decomposed into, $\phi = \phi_0(\psi) + \phi_1(\psi, \theta)$, an equipotential part on the flux surface $\phi_0(\psi)$ and a poloidally varying part $\phi_1(\psi, \theta)$ introduced to ensure the quasi-neutrality in the presence of centrifugal forces.

According to the neoclassical theory, developed in Refs. [25, 34], at the lowest order in $\delta = \rho_i/L$, where $\rho_i$ represents the ion Larmor radius and $L$ denotes the radial length scale, the distribution function is a Maxwellian. In addition, the temperature $T_{eq}(\psi)$ and toroidal rotation frequency $\Omega(\psi)$ are flux functions, which are constant on given magnetic flux surfaces, while the density $N_{eq}(\psi, \theta)$ is not a flux function due to the centrifugal force. Consequently, the Maxwellian distribution function in the rotating frame can be expressed as:

$$F_M = N(\psi, \theta) \left( \frac{m}{2\pi T_{eq}(\psi)} \right)^{3/2} \exp \left( -\frac{mu_{\parallel}^2}{2T_{eq}(\psi)} - \frac{\mu B_{eq}(\psi, \theta)}{T_{eq}(\psi)} \right). \quad (7)$$

Using the Hamiltonian derived in Eq. (6), the above expression can be rewritten as:

$$F_M = N_{eq}(\psi) \left( \frac{m}{2\pi T_{eq}(\psi)} \right)^{3/2} \exp \left( -\frac{H}{T_{eq}(\psi)} \right), \quad (8)$$

with the expression of the density distribution in a toroidally rotating plasma

$$N(\psi, \theta) = N_{eq}(\psi) \exp \left( -\frac{eZ\phi_0(\psi)}{T_{eq}(\psi)} - \frac{eZ\phi_1(\psi, \theta)}{T_{eq}(\psi)} + \frac{1}{2} \frac{m\Omega^2 R^2}{T_{eq}(\psi)} \right). \quad (9)$$
In the presence of strong rotation, impurities tend to accumulate in the outboard region due to the centrifugal force, while an induced potential $\phi_1$, which varies poloidally, acts to balance the resulting density asymmetry. These effects are more pronounced for heavy impurity ions because of their large mass and high-Z.

It should be noted that, unlike gyrokinetic formulation in the co-moving frame that requires the inclusion of both the Coriolis drift and centrifugal effect in the derivation of guiding-center equations [35–37], the equations of motion in GYSELA are solved in the laboratory frame. This implies that the effect of toroidal rotation is inherently present in the distribution function described in Eq. (8), without introducing any additional terms in the dynamics of the guiding-center.

### III. NUMERICAL RESULTS

In this section, we describe a set of nonlinear GYSELA simulations to explore the effect of toroidal rotation on impurity transport. To this end, we consider W as the main impurity due to its large mass $A_W = 184$ and high charge state $Z_W = 74$ when fully ionized. However, we restrict our analysis to the partially ionized state $Z_W = 40$ as W is not fully stripped even at the core temperature of tokamaks. We conduct the simulations in a circular geometry with a numerical resolution of $(N_r \times N_\theta \times N_\phi \times N_{v\|} \times N_\mu) = (256 \times 512 \times 32 \times 128 \times 64)$, using a time step of $\Delta t = 16\omega_{ci}^{-1}$. The velocity space parameters are set to $0 < \mu < 20$ and $-7 < v_{\|s} < 7$, to account for both pitch angle scattering and energy diffusion of the collision operator, as well as supersonic impurities [28]. The nonlinear gyrokinetic Vlasov equation defined in Eq. (1) is solved for both main ions and W, while electrons are simplified to follow a Boltzmann response within a flux surface. The initial radial profiles of density and temperature are established at $R/L_n = 2.2$ and $R/L_T = 6$ for both species, triggering ITG turbulence. Other parameters such as normalized gyro-radius $\rho_*= \rho_i/a = 1/190$, inverse aspect ratio $\epsilon = a/R_0 \simeq 0.2$, safety factor $q_{95} \simeq 3$ are kept constant throughout the present study. The radial profile of collisionality is determined by the normalized collisionality $\nu_s^* = \nu_i q R_0 / (v_{T_s} \epsilon^{3/2})$. For the analysis, we set the main ions in the banana regime with $\nu_D^* = 0.1$ and W in the Pfirsch-Schlüter (PS) regime with $\nu_W^* \simeq 23$ near the edge. The W concentration is restricted to the trace limit ($C_W = 10^{-6}$), ensuring that the background turbulence remains unaffected by its presence [38].
In this study, by increasing the value of $S_{mom}$, higher values of Mach number are achieved. The radial profile of the momentum source remains unchanged with $\rho_{s,mom} = 0.1$ and $L_{s,mom} = 0.1$. In addition, the amplitude and radial profile of the heat flux $S_{heat}$ are kept constant throughout the analysis, with $S_E = 0.017$, $\rho_{s,heat} = 0.012$ and $L_{s,heat} = 0.1$. Heat and momentum are injected only for deuterium, indicating that the transfer of heat and kinetic energy occurs via inter-species collisional effects. Figure 1 shows the radial profile of the different source terms for the case with $S_{mom} = 10^{-2}$. No inner boundary condition is applied to cover the magnetic axis ($\rho = 0$), while Dirichlet boundary conditions are chosen at the sepa-ratrix ($\rho = 1$) to satisfy $\phi_{sep} = 0$. A radial diffusion term, $D(\bar{F}_s) = 1/(rB_{\parallel s}^*)\partial_r(D_rB_{\parallel s}^*\partial_r\bar{F}_s)$, with a radial diffusion coefficient $D_r = 0.03$, is applied in the buffer region to mitigate numerical oscillations in the edge region.

To investigate its effect on density distribution and impurity transport, which is typically proportional to $\propto (1+M^2)$ [8], W ions are accelerated to the supersonic regime. Although the Mach number of main ions in fusion plasma is known to be low, typically $M_D < 0.4$, the Mach number of W can exceed the unity ($M_W > 1$) due to the mass ratio $M_W = \sqrt{A_W/A_D}M_D$.

<figure><img src="image_2.png" /><figcaption>FIG. 1: Radial profiles of heat and momentum sources applied to the main ions in GYSELA. No heat and momentum sources are applied to W. The buffer region is implemented to prevent boundary effects from affecting the bulk plasma.</figcaption></figure>

The numerical computation time required to reach a steady state is proportional to the gyro-Bohm scaling $\propto \rho_*^{-3}$ [39], which makes it infeasible to attain a steady state for multi-
species simulations. To overcome this limitation, we first perform a simulation without W impurities for a sufficiently long time ($\sim 10^5 \omega_{ci}^{-1}$) to develop the background turbulence and to saturate the toroidal velocity of main ions. Then, we carry out the second part of the simulation with W impurities for more than $10^5 \omega_{ci}^{-1}$ to get convergence on the toroidal rotation of W through collisional effects. Once the simulations are close to the quasi-steady state, all quantities are analyzed with a time window of $\sim 3000 \omega_{ci}^{-1}$ and are averaged over flux surfaces.

<figure><img src="image_3.png" /><figcaption>FIG. 2: Mach numbers for the main ions (a) and W (b) with different amplitudes of momentum sources plotted as a function of $\rho = r/a$ at $t = 2 \times 10^5 \omega_{ci}^{-1}$.</figcaption></figure>

The radial profiles of the Mach number for both species are illustrated in Figure 2 with different amplitudes of the momentum source. In the absence of the momentum injection $S_{mom} = 0$, the toroidal velocity, produced by turbulent intrinsic rotation [18, 19], exhibits low and flat profile for both species in the counter-current direction, i.e. $M_D \sim -0.02$ and $M_W \sim -0.6$. Once the momentum source is activated, the toroidal rotation increases and converges to a similar shape as the source profile illustrated in Figure 1. W is propelled into the supersonic regime for amplitudes stronger than $S_{mom} = 4 \times 10^{-3}$. In the core region, the value of $M_W$ reaches up to $M_W \simeq 4$ with $S_{mom} = 10^{-2}$. Despite the injection of momentum, the Mach number and its gradient for deuterium remains at low levels, with values of $M_D < 0.4$ and $|dv_{\phi,D}/dr| \simeq 0.4$, respectively. According to the criterion for the PVG instability to occur, which is $|dv_{\phi,D}/dr| > R/(L_n\sqrt{\tau+1})$ with $\tau = T_e/T_i$ [24], the present nonlinear simulations are not subject to PVG instabilities. Similar Mach numbers have been observed in JET experiments with NBI injection [40], where the Mach number
<figure><img src="image_4.png" /></figure>

for deuterium was found to be as high as $M_D \simeq 0.7$, yielding $M_W \simeq \sqrt{A_W/A_D}M_D \simeq 4.5$.

<figure><img src="image_5.png" /><figcaption>FIG. 3: Radial density and temperature profiles of the main ions and W with different values of momentum source measured at $t \simeq 2 \times 10^5 \omega_{ci}^{-1}$. The black solid line represents the initial profile at $t \simeq 10^5 \omega_{ci}^{-1}$ where W is introduced.</figcaption></figure>

The radial profile of the flux-averaged density and temperature for both the main ions and W are displayed in Figure 3. For deuterium, the density and temperature profiles are not significantly affected by the external momentum injection with the considered amplitudes of $2 \times 10^{-3} \le S_{mom} \le 10^{-2}$. The slight increase in $T_D$ is likely due to the balance between the heat source and the outward heat flux, rather than the effect of toroidal rotation. Unlike deuterium, we observe significant changes in W density with increasing the amplitude of $S_{mom}$, while the temperature profiles remain unaffected and saturate at the same level as the main ions, indicating energy transfer through collisional effects. The overall density profile implies a central accumulation in cases of strong rotation.
<figure><img src="image_6.png" /><figcaption>FIG. 4: Radial profile of the root-mean-square (rms) of the perturbed electric potential $\delta\phi = \phi - \langle\phi\rangle_{FS}$ for different values of momentum injection.</figcaption></figure>

Figure 4 describes the radial profile of the root-mean-square (rms) of the perturbed electric potential for different values of momentum injection. In the absence of momentum injection, the turbulence intensity remains below 4%. However, the injection of momentum tends to increase fluctuation levels in the region where the source is applied. Fluctuation levels in the region where $\rho > 0.5$ remain unchanged across different amplitudes, while an enhancement in fluctuation levels is observed in the inner region where ($\rho < 0.5$). This indicates that the injected momentum does not stabilize the turbulence, instead, it results in a slight increase in turbulent levels in the inner region. The background ion temperature gradient (ITG) turbulence consistently remains throughout this study.

Previous works with GYSELA have demonstrated the stabilizing effect on turbulence via the formation of transport barrier generated by a vorticity source term [41, 42]. Furthermore, recent research with the JET-ILW (ITER-like wall) tokamak reports that steep temperature gradients at the plasma periphery can prevent W accumulation in the core region by enhancing outward neoclassical convection [43]. The impact of transport barrier on impurity transport using GYSELA is currently under investigation.
## Poloidal asymmetry of W density in toroidally rotating plasma

Non-uniform distributions of W obtained from the GYSELA simulations of rotating plasma are shown in Figure 5. Three different cases are compared to elucidate the effect of toroidal rotation in driving an enhanced in-out asymmetry. In the absence of momentum injection $S_{mom} = 0$, an up-down asymmetry of W is observed mainly due to the background turbulence [44]. The degree of asymmetry is approximately 20%. As the amplitude of the source term increases, the geometry of asymmetry changes from up-down to in-out due to the centrifugal effects. In particular, the case with $S_{mom} = 10^{-2}$ yields an excess of approximately 40% of W localized in the outboard region.

<figure><img src="image_7.png" /><figcaption>FIG. 5: 2D poloidal section of perturbed W density obtained from the nonlinear GYSELA simulations. Three different cases are represented, i.e. without momentum source $S_{mom} = 0$ (left), the Mach number close to unity $M_W \sim 1$ with $S_{mom} = 4 \times 10^{-3}$ (center) and the high Mach number $M_W \simeq 4$ with $S_{mom} = 10^{-2}$ (right).</figcaption></figure>

An analytical formula for non-uniform density distribution in a toroidally rotating plasma can be obtained by solving the drift-kinetic equation for impurities in the Pfirsch-Schlüter regime [45]. By assuming an axisymmetric system and small inverse aspect ratio ($\epsilon \ll 1$), the poloidal asymmetry of impurity density at the lowest order of $\epsilon$ can be written as:

$$N_W = 1 + N_c \cos \theta + N_s \sin \theta, \qquad (10)$$

with the in-out geometrical factor

$$N_c = 2\epsilon \frac{M_0^2 - (1 + \gamma)g^2}{1 + (1 + \gamma)^2 g^2} \qquad (11)$$
and the up-down geometrical factor

$$N_s = 2\epsilon g \frac{1 + (1 + \gamma)M_0^2}{1 + (1 + \gamma)^2 g^2}. \qquad (12)$$

In the above expressions, the following dimensionless parameters are introduced, such as $M_0^2 = M^2 R_0^2 / R^2$ with

$$M^2 = \frac{m_z \omega^2 R^2}{2T_z} \left( 1 - \frac{Z m_i}{m_z} \frac{T_e}{T_e + T_i} \right) \qquad (13)$$

and

$$g = - \frac{m_i Z^2 I}{e \tau_{iz} \langle \mathbf{B} \cdot \nabla \theta \rangle} \left( \frac{1}{N_i} \frac{\partial N_i}{\partial \psi} - \frac{1}{2} \frac{1}{T_i} \frac{\partial T_i}{\partial \psi} \right), \qquad (14)$$

where we use

$$\gamma = \frac{eL_{\perp} \langle B^2 \rangle u}{T_i}, \qquad (15)$$

$$L_{\perp}^{-1} = -I \left( \frac{1}{N_i} \frac{\partial N_i}{\partial \psi} - \frac{3}{2} \frac{1}{T_i} \frac{\partial T_i}{\partial \psi} \right), \qquad (16)$$

and the poloidal rotation for the main ions

$$u \sim -0.33 \frac{I f_c}{e \langle B^2 \rangle} \frac{\partial T_i}{\partial \psi}. \qquad (17)$$

with the fraction of circulating particles $f_c \simeq 1 - 1.46\sqrt{\epsilon}$.

The first dimensionless parameter, defined in Eq. (13), corresponds to a modified Mach number. The geometrical factors defined in Eqs. (11, 12) scale quadratically with respect to the value of $M$. The second dimensionless parameter, denoted as $g$, characterizes the steepness of the density and temperature gradient profiles for the main ions. In standard neoclassical theory, this parameter is assumed to be small, indicating a weak frictional force. However, its impact becomes dominant when impurities localize in the outboard region, particularly for high-Z impurities ($g \propto Z^2$) or in the pedestal region.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6: (Left) Poloidal asymmetry of W density fluctuation $\tilde{N}_W$ obtained from the GYSELA simulation with $S_{mom} = 10^{-2}$, (center) reconstructed impurity density $\tilde{N}_W$ from Eq. (9) and (right) $\tilde{N}_W$ reconstructed from Eqs. (10-12).</figcaption></figure>

Figure 6 represents the non-uniform distribution of W in a toroidally rotating plasma, obtained from the GYSELA simulation on the left in Figure 6, and its comparison with analytical expressions from Eq. (9) and Eqs. (10-12). The other parameters in the analytical expressions are evaluated from the same GYSELA simulation with $S_{mom} = 10^{-2}$. Both theoretical approaches, namely the Hamiltonian in a rotating frame and the drift-kinetic equation, are in qualitative agreement with the GYSELA simulations regarding the enhanced in-out asymmetry due to the centrifugal effects.

The expression for the density distribution derived in Eq. (9) shows a level of poloidal asymmetry comparable to that of GYSELA, whereas a factor of 1/2 is observed when using the expression derived from the drift-kinetic equation. The main difference between the two expressions is the inclusion of the poloidally varying electric potential $\phi_1(\psi, \theta)$, which arises from toroidal rotation and turbulence. The discrepancies with neoclassical predictions can be attributed to the fact that the expressions in Eqs. (10-12) rely on neoclassical estimates of the poloidal ions velocity [46–49]. These estimates are not always consistent with experimental observations, suggesting that other mechanisms, such as intrinsic rotation induced by turbulence, may be important in this process.

Experimental observations [50–52] have revealed that a strong poloidal asymmetry, which is enhanced by the centrifugal force, can lead to a significant accumulation of W in the plasma core, primarily due to the inward neoclassical convection. Given these observations,
the following section will focus on exploring both neoclassical and turbulent transport of W in toroidally rotating plasma.

## IV. NEOCLASSICAL AND TURBULENT IMPURITY TRANSPORT IN A TOROIDALLY ROTATING PLASMA

The importance of including poloidal asymmetries in neoclassical impurity transport is now well recognized, and experimentally demonstrated that a strong non-uniform distribution in the outboard region can lead to core accumulation [21, 53]. An analytical expression of heavy impurity transport—Pfirsch-Schlüter transport—which takes into account such inhomogeneous poloidal density has been proposed in Ref. [10] assuming the main ions in the banana regime:

$$R\langle\Gamma_W^{\text{neo}} \cdot \nabla\psi\rangle_{\psi} \propto \left[-\frac{R}{L_{ni}} - H_{\text{neo}}\frac{R}{L_{Ti}} + \frac{1}{Z}\frac{R}{L_{nW}}\right] P_A, \quad (18)$$

$$\propto \left[ \left( -\frac{R}{L_{ni}} - H_{\text{sym}}\frac{R}{L_{Ti}} + \frac{1}{Z}\frac{R}{L_{nW}} \right) P_A - 0.33 P_B f_c \frac{R}{L_{Ti}} \right], \quad (19)$$

where $H_{\text{neo}} = H_{\text{sym}} + 0.33(P_B/P_A)f_c$ is the temperature screening coefficient with $H_{\text{sym}} = -1/2$ in absence of poloidal asymmetries. The geometric factors $P_A$ and $P_B$ related to the poloidal asymmetry are defined as

$$P_A = \frac{1}{2\epsilon^2} \frac{\langle B^2 \rangle}{\langle N_W \rangle} \left[ \left\langle \frac{N_W}{B^2} \right\rangle - \left\langle \frac{B^2}{N_W} \right\rangle^{-1} \right], \quad (20)$$

$$P_B = \frac{1}{2\epsilon^2} \frac{\langle B^2 \rangle}{\langle N_W \rangle} \left[ \left\langle \frac{N_W}{B^2} \right\rangle - \left\langle \frac{B^2}{N_W} \right\rangle^{-1} \right]. \quad (21)$$

The values of $P_A$ and $P_B$ depend on the magnitude of poloidal asymmetry, as well as its geometry, namely in-out or up-down. A simple approach to calculate these geometrical factors was proposed in the large aspect ratio limit by considering $b = B/B_0 = 1 - \epsilon \cos\theta$ and using a similar expression to Eq. (10) for $N = 1 + \delta \cos\theta + \Delta \sin\theta$ with $\delta$ (in-out) and $\Delta$ (up-down) asymmetry parameters [10]. By keeping the second order accuracy, the geometrical factors can be rewritten as

$$\left\langle \frac{N_z}{B^2} \right\rangle - \left\langle \frac{B^2}{N_z} \right\rangle^{-1} = 2\epsilon(\epsilon + \delta) + \frac{\delta^2 + \Delta^2}{2}, \quad (22)$$

$$1 - \left\langle \frac{B^2}{N_z} \right\rangle = \epsilon\delta + \frac{\delta^2 + \Delta^2}{2}, \quad (23)$$
where we use $\langle N_Z \rangle = \langle b^2 \rangle = 1$. The radial profile of the simplified geometrical factors is illustrated in Figure 7 from the nonlinear simulations. In case of a uniform density distribution corresponding to the case with $\delta = \Delta = 0$, the geometrical factors are equal to $P_A = 1$ and $P_B = 0$ recovering the standard neoclassical transport derived in Ref. [14]. In line with the previous results from the 2D poloidal W distribution, the toroidal rotation enhances the effect of poloidal asymmetry, especially in the outboard region ($0.7 < \rho < 0.9$) due to the centrifugal effects. The values of $P_A$ and $P_B$ scale with $1/(2\epsilon^2)$, which are usually very large numbers. In the considered simulations, the geometrical factors in Figure 7 are amplified by $1/(2\epsilon^2) \simeq 40$ at $\rho = 0.5$ and $1/(2\epsilon^2) \simeq 12$ at $\rho = 0.9$ yielding $P_A \simeq 4$ and $P_B \simeq 0.5$ at mid-radius. Indeed, a series of recent experimental works highlight that the poloidal asymmetry caused by external heating system, such as NBI and ion cyclotron resonance heating (ICRH) heating systems, can increase the factors $P_A$ and $P_B$ up to 100 and reduce the temperature screening coefficient $H_{neo}$ leading to a core accumulation near the magnetic axis [54].

<figure><img src="image_10.png" /><figcaption>FIG. 7: Radial profile of the geometrical factors defined in Eqs. (22-23) obtained from the nonlinear GYSELA simulations. The green line represents the case with the uniform density distribution.</figcaption></figure>
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_11.png" /><figcaption>FIG. 8: Neoclassical W flux as a function of the $\delta$ (a) and $\Delta$ (b) asymmetry parameters. The reconstructed density (c) $\tilde{N}_W = \delta \cos \theta + \Delta \sin \theta$ is used to calculate the neoclassical flux in Eq. (19), while other terms, i.e. $N_i$, $T_i$, $T_W$, $\nu_{zi}$, are taken from the GYSELA simulations. Flux is normalized to the absolute value of the symmetric case with $\delta = \Delta = 0$.</figcaption></figure>

In general, neoclassical impurity transport depends on various plasma parameters for the main ions and impurities, i.e. $\Gamma_W^{\text{neo}} \propto f(N_i, N_W, T_i, T_W, \nu_{zi})$. However, unlike the standard neoclassical theory, a non-uniform distribution of impurity density has to be included to accurately predict the neoclassical flux. In Figure 8, the effect of inhomogeneous W density on neoclassical W flux is represented by performing a scan of different values in $\delta$ and $\Delta$ from $\tilde{N}_W = \delta \cos \theta + \Delta \sin \theta$. The reconstructed 2D asymmetry of W distribution is then inserted in Eq. (19) to estimate the neoclassical flux, while other terms, measured from the GYSELA simulations, are kept constant. The poloidal asymmetry scan clearly indicates an increase of inward neoclassical flux as the asymmetry grows. In the case of up-down asymmetry (Figure 8b), the trend is found to be symmetric, while the in-out asymmetry study (Figure 8a) reveals a strong impact of impurity localization in the outboard region on
neoclassical transport, even causing a change in the order of magnitude.

Similarly to the neoclassical particle flux defined in Eq. (19), the turbulent W flux can be decomposed into various terms in the trace limit of impurities ($\alpha = N_W Z_W^2 \ll 1$) [55] as

$$ \frac{R\Gamma_W}{N_W} = -D_W \frac{\partial \ln N_W}{\partial r} + C_T \frac{\partial T_W}{\partial r} + C_u \frac{\partial v_{\phi,W}}{\partial r} + C_p \frac{\partial q}{\partial r}, \quad (24) $$

where each term represents turbulent diffusion, thermo-diffusion, roto-diffusion and curvature pinch respectively. The characteristic of these terms depends on the nature of dominant instabilities or magnetic curvature profiles that are present in the plasma. For example, the diffusion coefficient increases (decreases) with impurity charge Z for the case of Trapped-Electron-Mode (TEM) turbulence (ITG turbulence). In addition, thermo-diffusion and roto-diffusion are directed inward (outward) for TEM turbulence (ITG turbulence), while curvature-driven transport is directed inward (outward) for positive (negative) magnetic shear. The same trend was qualitatively reproduced in the recent studies with the bounced-averaged gyrokinetic code TERESA using the quasi-linear approach for the impurity transport [56, 57]. A more comprehensive review on the physical mechanisms underlying each term can be found in Ref. [58].

In the present study, the impact of toroidal rotation on turbulent flux is primarily mediated through the roto-diffusion term. The roto-diffusion term, proportional to the radial gradient of toroidal velocity, is known to be significant for W as the magnitude of the coefficient scales as $C_u \propto A_W/Z_W$.

For the analysis of the nonlinear GYSELA simulations, we adopt the following definitions for turbulent and neoclassical particle flux of W [29, 59],

$$ \langle \Gamma_W^{\text{turb}} \cdot \nabla \psi \rangle_{\psi} = \left\langle \int d^3 v \bar{F}_W \bar{v}_E^{n \neq 0} \cdot \nabla \psi \right\rangle_{\psi}, \quad (25) $$

$$ \langle \Gamma_W^{\text{neo}} \cdot \nabla \psi \rangle_{\psi} = \left\langle \int d^3 v \bar{F}_W (v_{D,W} + \bar{v}_E^{n=0}) \cdot \nabla \psi \right\rangle_{\psi}, \quad (26) $$

where n is the toroidal mode number, $\bar{v}_E$ is the $E \times B$ drift of the guiding-center and $v_D$ is the diamagnetic drift. Note that the axisymmetric component of the $E \times B$ drift, namely $\bar{v}_E^{n=0}$, also contains turbulence-driven terms, sometimes called convective cell [44]. In that respect, the splitting between turbulent and neoclassical contributions contained in these definitions Eqs. (25-26) should not be understood too strictly.
<figure><img src="image_12.png" /><figcaption>FIG. 9: Time cumulative neoclassical (left) and turbulent (right) W flux averaged over the radial position between $\rho = 0.1$ and $\rho = 0.9$ with different magnitudes of the momentum source. The time is set to $t = 0$ when W is injected into the simulations.</figcaption></figure>

Figure 9 illustrates the cumulative turbulent and neoclassical fluxes obtained from GY-SELA simulations at various levels of momentum injection. The fluxes are integrated both in time and volume within the range of $0.1 < \rho < 0.9$. Consistently with the previous findings in Figure 8, the neoclassical inward convection resulting from an enhanced in-out asymmetry of W density is observed to be proportional to the magnitude of toroidal rotation. A recent analysis of W accumulation in KSTAR experiments using the drift-kinetic code NEO [60] revealed that neoclassical inward convection is maximized at a specific toroidal rotation value, beyond which the outward convection increases [51]. The non-monotonic relationship between the toroidal rotation and neoclassical inward convection can be attributed to the different dependencies of the convection terms on the main ion density and temperature profiles, as well as collisionality. An analytical investigation of neoclassical transport shows that the temperature screening of impurities is enhanced at sufficiently high Mach number and low collisionality, while this effect decreases as collisionality increases [61]. In the present study, however, the impact of toroidal rotation on the main profiles is found to be negligible, as depicted in Figure 3, and only a high collisional case with W is considered. This explains the linear relationship observed between neoclassical inward convection and the magnitude of toroidal rotation in Figure 9.

Although the toroidal rotation has a direct impact on neoclassical transport owing to
the presence of strong density inhomogeneity, its influence on turbulent transport in Figure 9 (right) is observed to be relatively small. The enhancement of outward turbulent flux resulting from the increase in toroidal rotation, as depicted in Figure 9, can be attributed to the increased turbulence intensity observed in Figure 4 and the turbulent pinch terms defined in Eq. (24), i.e., roto-diffusion and thermo-diffusion terms. Although both terms are known to drive outward convection in the case of ITG turbulence, the roto-diffusion predominates over the thermo-diffusion term, as the former is proportional to $A_W/Z_W = 4.6$, while the latter is proportional to $1/Z_W = 0.025$. In fact, for heavy impurities with strong toroidal rotation, the magnitude of roto-diffusion can be as significant as the diffusive part [23].

<figure><img src="image_13.png" /><figcaption>FIG. 10: Time evolution of the radial gradient of the toroidal velocity (left), and the time evolution of the turbulent flux of W (right), averaged over the radial position between $\rho = 0.1$ and $\rho = 0.9$ for different values of the injected momentum. The dashed vertical line indicates the time at which the effect of roto-diffusion becomes significant. The time is set to $t = 0$ when W is injected into the simulations.</figcaption></figure>

The impact of the roto-diffusion term on the W turbulent flux becomes apparent when we compare the time evolution of the radial gradient of toroidal velocity and the turbulent flux for W, as depicted in Figure 10. The radial gradient of toroidal rotation gradually increases over time in proportion to the injected momentum. For W, the values of $|dv_{\varphi,W}/dt|$ are observed to exceed unity for amplitudes above $S_{mom} = 4 \times 10^{-3}$ due to its heavy mass, whereas the gradient remains low at $|dv_{\varphi,D}/dt| \simeq 0.4$ for the main ions. As shown in
<figure><img src="image_14.png" /></figure>

Figure 9, the outward turbulent flux emerges at $t = 10^5 \omega_{ci}^{-1}$ for amplitudes greater than $S_{mom} = 6 \times 10^{-3}$. The results demonstrate that the influence of the roto-diffusion term becomes more significant compared to other terms in Eq. (24) as the absolute value of $|dv_{\phi,w}/dr|$ gradually increases. Additionally, it is noteworthy that centrifugal effects also have a non-negligible impact on other convection terms, particularly in the case of heavy impurities [55]. Consequently, the sign and magnitude of the turbulent flux depend on the interplay between different convection terms, which typically results in an outward flux in the case of ITG turbulence.

<figure><img src="image_15.png" /><figcaption>FIG. 11: Core W-accumulation in time for different amplitudes of the toroidal momentum source. The core density $N_{W, core}$ is integrated between $\rho = 0.1$ and $\rho = 0.3$. The strong toroidal velocity leads to an enhanced core accumulation.</figcaption></figure>

Figure 11 presents the W accumulation obtained from the total flux, $\Gamma_W^{tot} = \Gamma_W^{neo} + \Gamma_W^{urb}$, within the radial range of $\rho = 0.1$ and $\rho = 0.3$ for different values of the injected momentum. No significant changes in the core W density are observed up to the time $t = 0.5 \times 10^5 \omega_{ci}^{-1}$, which corresponds to the momentum transfer time from the main ions to W. Following this transfer time, a more significant accumulation of W is observed, approximately 20% of increase in the $S_{mom} = 10^{-2}$ scenario. The central accumulation of W in GYSELA simulations are commonly observed even in the absence of external toroidal rotation [26]. The numerical results obtained in the present study are consistent with experimental evidence from different tokamaks [51, 53, 54] where central NBI heating yields a central deposition of
<figure><img src="image_16.png" /></figure>

W. To prevent the accumulation of W within the core region, the proper use of ICRH and ECRH systems can be helpful in controlling the W within the core region [9, 33, 52, 62–66]. This is achieved by flattening the density profile of main ions and enhancing the temperature gradient of the bulk plasma, which facilitates the outward turbulent diffusion of W and an increase of the temperature screening effect. In addition, the temperature anisotropy that arises from ICRH is known to reduce poloidal asymmetry, thereby decreasing inward neoclassical convection [9, 67].

V. CONCLUSION

In the present paper, the effects of toroidal rotation on both turbulent and neoclassical W transport are addressed using the nonlinear, global, flux-driven 5D gyrokinetic code GYSELA. The injection of toroidal momentum via the external source term drives W to the supersonic regime, while the toroidal velocity of the main ions remains low. By leveraging the multi-species collisional operator and external source terms implemented in GYSELA, both turbulent and neoclassical W transport in a toroidally rotating plasma are treated self-consistently. Unlike the conventional approach, which treats neoclassical and turbulent transport separately, the self-consistent treatment used in this study captures the interplay between these two channels of W transport.

Through a series of nonlinear simulations involving different levels of momentum injection, it has been observed that centrifugal forces cause W to accumulate in the outboard region, creating an in-out poloidal asymmetry. This asymmetry has been found to significantly enhance the neoclassical inward convection, leading to a central accumulation of W in the case of strong plasma rotation. The core W accumulation accelerated by the centrifugal force and the resulting poloidal asymmetry, cannot be offset by the turbulent flux. However, as the momentum injection continues, the roto-diffusion term becomes as significant as the diffusive part once the radial gradient of the rotation velocity exceeds a certain threshold, driving outward turbulent flux in the case of ITG turbulence.

The numerical results obtained in the present study are in qualitative agreement with theoretical predictions. However, for a more comprehensive modelling of W transport in GYSELA, it could be interesting to include the time evolution of the density and temperature profiles of main ions, taking into account their response to other heating systems. These
<figure><img src="image_17.png" /></figure>

additional heating systems are known to mitigate core W accumulation by reducing poloidal W asymmetry and enhancing the temperature screening effect.

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. The simulations presented herein were carried out on the CINECA Marconi supercomputer under FUA35 GSNTITE project.

[1] M. Kaufmann and R. Neu. Tungsten as first wall material in fusion devices. *Fusion Eng. Des.*, **82** 521, (2007).

[2] M.Z. Tokar, J. Rapp, G. Bertschinger, L. Konen, H.R. Koslowski, A. Kramer-Flacken, V. Philipps, U. Samm, and B. Unterberg. Nature of high-Z impurity accumulation in tokamaks. *Nucl. Fusion*, **37** 1691, (1997).

[3] F Köchl, A Loarte, E de la Luna, V Parail, G Corrigan, D Harting, I Nunes, C Reux, F G Rimini, A Polevoi, M Romanelli, and JET Contributors. W transport and accumulation control in the termination phase of JET H-mode discharges and implications for ITER. *Plasma Phys. Control. Fusion*, **60** 074008, (2018).

[4] T. Pütterich, R. Neu, R. Dux, A.D. Whiteford, M.G. O'Mullane, H.P. Summers, and the ASDEX Upgrade Team. Calculation and experimental test of the cooling factor of tungsten. *Nucl. Fusion*, **50** 025012, (2016).

[5] J. L. Terry, E. S. Marmar, K. I. Chen, and H. W. Moos. Observation of poloidal asymmetry in impurity-ion emission due to $\nabla b$ drifts. *Phys. Rev. Lett.*, **39** 1615, (1977).
[6] L C Ingesson, H Chen, P Helander, and M J Mantsinen. Comparison of basis functions in soft x-ray tomography and observation of poloidal asymmetries in impurity density. *Plasma Phys. Control. Fusion*, **42** 161, (2000).

[7] J.A. Wesson. Poloidal distribution of impurities in a rotating tokamak plasma. *Nucl. Fusion*, **37** 577, (1997).

[8] M. Romanelli and M Ottaviani. Effects of density asymmetries on heavy impurity transport in a rotating tokamak plasma. *Plasma Phys. Control. Fusion*, **40** 1767, (1998).

[9] R. Bilato, O. Maj, and C. Angioni. Modelling the influence of temperature anisotropies on poloidal asymmetries of density in the core of rotating plasmas. *Nucl. Fusion*, **54** 072003, (2014).

[10] C. Angioni and P. Helander. Neoclassical transport of heavy impurities with poloidally asymmetric density distribution in tokamaks. *Plasma Phys. Control. Fusion*, **56** 124001, (2014).

[11] J. Dominski, C. S. Chang, R. Hager, P. Helander, S. Ku, and E. S. Yoon. Study of up-down poloidal density asymmetry of high-Z impurities with the new impurity version of XGCa. *J. Plasma Phys.*, **85** 905850510, (2019).

[12] Patrick Maget, Judith Frank, Timothée Nicolas, Olivier Agullo, Xavier Garbet, and Hinrich Lütjens. Natural poloidal asymmetry and neoclassical transport of impurities in tokamak plasmas. *Plasma Phys. Control. Fusion*, **62** 025001, (2020).

[13] F. L. Hinton and R. D. Hazeltine. Theory of plasma transport in toroidal confinement systems. *Rev. Mod. Phys.*, **48** 239, (1976).

[14] S. P. Hirshman and D. J. Sigmar. Neoclassical transport of impurities in tokamak plasmas. *Nucl. Fusion*, **21** 1079, (1981).

[15] F. J. Casson, A. G. Peeters, C. Angioni, Y. Camenen, W. A. Hornsby, A. P. Snodin, and G. Szepesi. Gyrokinetic simulations including the centrifugal force in a rotating tokamak plasma. *Phys. Plasmas*, **17** 102305, (2010).

[16] F.J. Casson, R.M. McDermott, C. Angioni, Y. Camenen, R. Dux, E. Fable, R. Fischer, B. Geiger, P. Manas, L. Menchero, G. Tardini, and the ASDEX Upgrade Team. Validation of gyrokinetic modelling of light impurity transport including rotation in ASDEX Upgrade. *Nucl. Fusion*, **53** 063026, (2013).

[17] E. A. Belli, J. Candy, and C. Angioni. Pfirsch-schlüter neoclassical heavy impurity transport in a rotating plasma. *Plasma Phys. Control. Fusion*, **56** 124002, (2014).
[18] A. Bortolon, B. P. Duval, A. Pochelon, and A. Scarabosio. Observation of spontaneous toroidal rotation inversion in ohmically heated tokamak plasmas. *Phys. Rev. Lett.*, **97** 235003, (2006).

[19] J.E. Rice, A. Ince-Cushman, J.S. deGrassie, L.-G. Eriksson, Y. Sakamoto, A. Scarabosio, A. Bortolon, K.H. Burrell, B.P. Duval, C. Fenzi-Bonizec, M.J. Greenwald, R.J. Groebner, G.T. Hoang, Y. Koide, E.S. Marmar, A. Pochelon, and Y. Podpaly. Inter-machine comparison of intrinsic toroidal rotation in tokamaks. *Nucl. Fusion*, **47** 1618, (2007).

[20] P. Helander. Neoclassical transport in a rotating impure plasma. *Phys. Plasmas*, **5** 1209, (1998).

[21] F J Casson, C Angioni, E A Belli, R Bilato, P Mantica, T Odstrcil, T Pütterich, M Valisa, L Garzotti, C Giroud, J Hobirk, C F Maggi, J Mlynar, M L Reinke, JET EFDA Contributors, and ASDEX-Upgrade Team. Theoretical description of heavy impurity transport and its application to the modelling of tungsten in JET and ASDEX upgrade. *Plasma Phys. Control. Fusion*, **57** 014031, (2015).

[22] S. Breton, F. J. Casson, C. Bourdelle, C. Angioni, E. Belli, Y. Camenen, J. Citrin, X. Garbet, Y. Sarazin, M. Sertoli, and JET Contributors. High Z neoclassical transport: Application and limitation of analytical formulae for modelling JET experimental parameters. *Phys. Plasmas*, **25** 012303, (2018).

[23] Y. Camenen, A. G. Peeters, C. Angioni, F. J. Casson, W. A. Hornsby, A. P. Snodin, and D. Strintzi. Impact of the background toroidal rotation on particle and heat turbulent transport in tokamak plasmas. *Phys. Plasmas*, **16** 012503, (2009).

[24] P. J. Catto, M. N. Rosenbluth, and C. S. Liu. Parallel velocity shear instabilities in an inhomogeneous plasma with a sheared magnetic field. *Phys. Fluids*, **16** 1719, (1973).

[25] F. L. Hinton and S. K. Wong. Neoclassical ion transport in rotating axisymmetric plasmas. *Phys. Fluids*, **28** 3082, (1985).

[26] 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. *Nucl. Fusion*, **61** 046037, (2021).

[27] V. Grandgirard, J. Abiteboul, J. Bigot, T. Cartier-Michaud, N. Crouseilles, G. Dif-Pradalier, Ch. Ehrlacher, D. Esteve, X. Garbet, Ph. Ghendrih, G. Latu, M. Mehrenberger, C. Norscini, Ch. Passeron, F. Rozar, Y. Sarazin, E. Sonnendrücker, A. Strugarek, and D. Zarzoso. A 5D gyrokinetic full-F global semi-lagrangian code for flux-driven ion turbulence simulations.
Comp. Phys. Comm., 207 35, (2016).

[28] P. Donnel, X. Garbet, Y. Sarazin, V. Grandgirard, Y. Asahi, N. Bouzat, E. Caschera, G. Dif-Pradalier, C. Ehrlacher, P. Ghendrih, C. Gillot, G. Latu, and C. Passeron. A multi-species collisional operator for full-F global gyrokinetics codes: Numerical aspects and verification with the GYSELA code. Comput. Phys. Commun., 234 1, (2019).

[29] D. Estève, Y. Sarazin, X. Garbet, V. Grandgirard, S. Breton, P. Donnel, Y. Asahi, C. Bourdelle, G. Dif-Pradalier, C. Ehrlacher, C. Emeriau, Ph. Ghendrih, C. Gillot, G. Latu, and C. Passeron. Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport. Nucl. Fusion, 58 036013, (2018).

[30] Y. Idomura, K. Obrejan, Y. Asahi, and M. Honda. Dynamics of enhanced neoclassical particle transport of tracer impurity ions in ion temperature gradient driven turbulence. Phys. Plasmas, 28 012501, (2021).

[31] X. Garbet, Y. Idomura, L. Villard, and T.H. Watanabe. Gyrokinetic simulations of turbulent transport. Nucl. Fusion, 50 043002, (2010).

[32] Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, Ph. Ghendrih, G. Latu, A. Strugarek, G. Dif-Pradalier, P.H. Diamond, S. Ku, C.S. Chang, B.F. McMillan, T.M. Tran, L. Villard, S. Jolliet, A. Bottino, and P. Angelino. Predictions on heat transport and plasma rotation from global gyrokinetic simulations. Nucl. Fusion, 51 103023, (2011).

[33] R Dux, R Neu, A G Peeters, G Pereverzev, A Mück, F Ryter, J Stober, and ASDEX Upgrade Team. Influence of the heating profile on impurity transport in ASDEX Upgrade. Plasma Phys. Control. Fusion, 45 1815, (2003).

[34] S. K. Wong. Transport of impure plasma with arbitrary toroidal rotation. Phys. Fluids, 30 818, (1987).

[35] A. G. Peeters, D. Strintzi, Y. Camenen, C. Angioni, F. J. Casson, W. A. Hornsby, and A. P. Snodin. Influence of the centrifugal force and parallel dynamics on the toroidal momentum transport due to small scale turbulence in a tokamak. Phys. Plasmas, 16 042310, (2009).

[36] A. J. Brizard. Nonlinear gyrokinetic Vlasov equation for toroidally rotating axisymmetric tokamaks. Phys. Plasmas, 2 459, (1995).

[37] H. Sugama. Modern gyrokinetic formulation of collisional and turbulent transport in toroidally rotating plasmas. Rev. Mod. Plasma Phys., 19, (2017).

[38] M. Lesur, C. Djerroud, K. Lim, E. Gravier, M. Idouakass, J. Moritz, J. Médina, T. Réveillé,
T. Drouot, T. Cartier-Michaud, and X. Garbet. Validity limits of the passive treatment of impurities in gyrokinetic tokamak simulations. *Nucl. Fusion*, **60** 036016, (2020).

[39] R. E. Waltz, J. C. DeBoo, and M. N. Rosenbluth. Magnetic-field scaling of dimensionally similar tokamak discharges. *Phys. Rev. Lett.*, **65** 2390, (1990).

[40] P.C. de Vries, M.-D. Hua, D.C. McDonald, C. Giroud, M. Janvier, M.F. Johnson, T. Tala, K.-D. Zastrow, and JET EFDA Contributors. Scaling of rotation and momentum confinement in JET plasmas. *Nucl. Fusion*, **48** 065006, (2008).

[41] A. Strugarek, Y. Sarazin, D. Zarzoso, J. Abiteboul, A. S. Brun, T. Cartier-Michaud, G. Dif-Pradalier, X. Garbet, Ph. Ghendrih, V. Grandgirard, G. Latu, C. Passeron, and O. Thomine. Unraveling quasiperiodic relaxations of transport barriers with gyrokinetic simulations of tokamak plasmas. *Phys. Rev. Lett.*, **111** 145001, (2011).

[42] 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. *Nucl. Fusion*, **62** 126026, (2022).

[43] A.R. Field, F.J. Casson, D. Fajardo, C. Angioni, C.D. Challis, J. Hobirk, A. Kappatou, Hyun-Tae Kim, E. Lerche, A. Loarte, J. Mailloux, and JET Contributors. Peripheral temperature gradient screening of high-Z impurities in optimised ‘hybrid’ scenario H-mode plasmas in JET-ILW. *Nucl. Fusion*, **63** 016028, (2023).

[44] P Donnel, X Garbet, Y Sarazin, Y Asahi, F Wilczynski, E Caschera, G Dif-Pradalier, P Ghendrih, and C Gillot. Turbulent generation of poloidal asymmetries of the electric potential in a tokamak. *Plasma Phys. Control. Fusion*, **61** 014003, (2018).

[45] T. Fülöp and P Helander. Nonlinear neoclassical transport in a rotating impure plasma with large gradients. *Phys. Plasmas*, **6** 3066, (1999).

[46] Y. B. Kim, P. H. Diamond, and R. J. Groebner. Neoclassical poloidal and toroidal rotation in tokamaks. *Phys. Fluids B*, **3** 2050, (1991).

[47] S. Newton and P. Helander. Neoclassical momentum transport in an impure rotating tokamak plasma. *Phys. Plasmas*, **13** 012505, (2006).

[48] W. M. Stacey, R. W. Johnson, and J. Mandrekas. A neoclassical calculation of toroidal rotation profiles and comparison with DIII-D measurements. *Phys. Plasmas*, **13** 062508, (2006).

[49] S. K. Wong, V. S. Chan, and W. M. Solomon. Poloidal velocity of impurity ions in neoclassical theory. *Phys. Plasmas*, **15** 082503, (2008).
[50] J. Bucalossi, J. Achard, O. Agullo, T. Alarcon, L. Allegretti, H. Ancher, G. Antar, S. Antusch, V. Anzallo, C. Arnas, et al. Operating a full tungsten actively cooled tokamak: overview of WEST first phase of operation. *Nucl. Fusion*, **62** 042007, (2022).

[51] Hyojong Lee, Hyeonjun Lee, Yoon Seong Han, Jiheon Song, E. A. Belli, Wonho Choe, Jisung Kang, Jekil Lee, J. Candy, and Jungpyo Lee. Neoclassical transport analysis of high rotational trace limit tungsten impurities in KSTAR tokamak. *Phys. Plasmas*, **29** 022504, (2022).

[52] Shengyu Shi, Jiale Chen, Clarisse Bourdelle, Xiang Jian, Tomáš Odstrčil, Andrea M. Garofalo, Yunxin Cheng, Yan Chao, Ling Zhang, Yanmin Duan, Mingfu Wu, Fang Ding, Jinping Qian, and Xiang Gao. Understanding core heavy impurity transport in a hybrid discharge on EAST. *Nucl. Fusion*, **62** 066032, (2022).

[53] C. Angioni, F. J. Casson, P. Mantica, T. Pütterich, M. Valisa, E. A. Belli, R. Bilato, C. Giroud, P. Helander, and JET Contributors. The impact of poloidal asymmetries on tungsten transport in the core of JET H-mode plasmas. *Phys. Plasmas*, **22** 055902, (2015).

[54] F.J. Casson, H. Patten, C. Bourdelle, S. Breton, J. Citrin, F. Koechl, M. Sertoli, C. Angioni, Y. Baranov, R. Bilato, E.A. Belli, C.D. Challis, G. Corrigan, A. Czarnecka, O. Ficker, L. Frassinetti, L. Garzotti, M. Goniche, J.P. Graves, T. Johnson, K. Kirov, P. Knight, E. Lerche, M. Mantsinen, J. Mylnar, M. Valisa, and JET contributors. Predictive multi-channel flux-driven modelling to optimise ICRH tungsten control and fusion performance in JET. *Nucl. Fusion*, **60** 066029, (2020).

[55] C. Angioni, F. J. Casson, C. Veth, and A. G. Peeters. Analytic formulae for centrifugal effects on turbulent transport of trace impurities in tokamak plasmas. *Phys. Plasmas*, **19** 122311, (2012).

[56] E. Gravier, M. Lesur, X. Garbet, Y. Sarazin, J. Médina, K. Lim, and M. Idouakass. Diffusive impurity transport driven by trapped particle turbulence in tokamak plasmas. *Phys. Plasmas*, **26** 082306, (2019).

[57] K Lim, E Gravier, M Lesur, X Garbet, Y Sarazin, and J Médina. Impurity pinch generated by trapped particle driven turbulence. *Plasma Phys. Control. Fusion*, **62** 095018, (2020).

[58] C. Angioni. Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments. *Plasma Phys. Control. Fusion*, **63** 073001, (2021).

[59] P Donnel, X Garbet, Y Sarazin, V Grandgirard, N Bouzat, E Caschera, G Dif-Pradalier, P Ghendrih, C Gillot, G Latu, and C Passeron. Neoclassical impurity flux in presence of
turbulent generated poloidal asymmetries and pressure anisotropy. *Plasma Phys. Control. Fusion*, **61** 044006, (2019).

[60] E. A. Belli and J. Candy. Full linearized Fokker-Planck collisions in neoclassical transport simulations. *Plasma Phys. Control. Fusion*, **54** 015015, (2011).

[61] D Fajardo, C Angioni, F J Casson, A R Field, P Maget, P Manas, the ASDEX Upgrade Team, and JET Contributors. Analytical model for the combined effects of rotation and collisionality on neoclassical impurity transport. *Plasma Phys. Control. Fusion*, **65** 035021, (2023).

[62] M. E. Puiatti, M. Valisa, C. Angioni, L. Garzotti, P. Mantica, M. Mattioli, L. Carraro, I. Coffey, C. Sozzi, and JET-EFDA contributors. Analysis of metallic impurity density profiles in low collisionality Joint European Torus H-mode and L-mode plasmas. *Phys. Plasmas*, **13** 042501, (2011).

[63] M. Valisa, L. Carraro, I. Predebon, M.E. Puiatti, C. Angioni, I. Coffey, C. Giroud, L. Lauro Taroni, B. Alper, M. Baruzzo, P. Belo daSilva, P. Buratti, L. Garzotti, D. Van Eester, E. Lerche, P. Mantica, V. Naulin, T. Tala, M. Tsalas, and JET-EFDA contributors. Metal impurity transport control in JET H-mode plasmas with central ion cyclotron radiofrequency power injection. *Nucl. Fusion*, **51** 033002, (2011).

[64] M. Sertoli, C. Angioni, and T. Odstrcil. Parametric dependencies of the experimental tungsten transport coefficients in ICRH and ECRH assisted ASDEX Upgrade H-modes. *Phys. Plasmas*, **24** 112503, (2017).

[65] M Goniche, R J Dumont, V Bobkov, P Buratti, S Brezinsek, C Challis, L Colas, A Czarnecka, P Drewelow, N Fedorczak, J Garcia, C Giroud, M Graham, J P Graves, J Hobirk, P Jacquet, E Lerche, P Mantica, I Monakhov, P Monier-Garbet, M F F Nave, C Noble, I Nunes, T Pütterich, F Rimini, M Sertoli, M Valisa, D Van Eester, and JET Contributors. Ion cyclotron resonance heating for tungsten control in various JET H-mode scenarios. *Plasma Phys. Control. Fusion*, **59** 055001, (2017).

[66] C. Angioni, M. Sertoli, R. Bilato, V. Bobkov, A. Loarte, R. Ochoukov, T. Odstrcil, T. Pütterich, J. Stober, and The ASDEX Upgrade Team. A comparison of the impact of central ECRH and central ICRH on the tungsten behaviour in ASDEX Upgrade H-mode plasmas. *Nucl. Fusion*, **57** 056015, (2017).

[67] R. Bilato, T. Odstrcil, F.J. Casson, C. Angioni, M. Brambilla, Ye.O. Kazakov, and E. Poli. The impact of the ion-cyclotron-resonance location on the poloidal asymmetries of impurity
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

density in an ICRF-heated rotating plasma. Nucl. Fusion, 57 056020, (2017).
<figure><img src="image_18.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_21.png" /></figure>
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_22.png" /></figure>
<figure><img src="image_23.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_24.png" /></figure>
<figure><img src="image_25.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_26.png" /></figure>
<figure><img src="image_27.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_28.png" /></figure>
<figure><img src="image_29.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

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

<figure><img src="image_31.png" /></figure>
<figure><img src="image_32.png" /></figure>

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

<figure><img src="image_33.png" /></figure>
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

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

<figure><img src="image_35.png" /></figure>
<figure><img src="image_36.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

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

<figure><img src="image_38.png" /></figure>
<figure><img src="image_39.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0157428

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