# Tessellation-based analysis of impurity clustering in the edge plasma of tokamaks

Zetao Lin, Thibault Maurel-Oujia, Benjamin Kadoch, Saddrudin Benkadda, Kai Schneider

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

## To cite this version:

Zetao Lin, Thibault Maurel-Oujia, Benjamin Kadoch, Saddrudin Benkadda, Kai Schneider. Tessellation-based analysis of impurity clustering in the edge plasma of tokamaks. Journal of Plasma Physics, 2025, 91 (1), pp.E30. 10.1017/S0022377824001259 hal-04941575

HAL Id: hal-04941575

https://hal.science/hal-04941575v1

Submitted on 11 Feb 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
Tessellation-based analysis of impurity clustering
in the edge plasma of tokamaks

Zetao Lin<sup>1</sup>, Thibault Maurel-Oujia<sup>1</sup>, Benjamin Kadoch<sup>2</sup>,
Saddrudin Benkadda<sup>3</sup> and Kai Schneider<sup>1,†</sup>

(Received 29 April 2024; revised 17 September 2024; accepted 20 September 2024)

Confinement quality in fusion plasma is influenced significantly by the presence of heavy impurities, which can lead to radiative heat loss and reduced confinement. This study explores the clustering of heavy impurity, i.e. tungsten in edge plasma, using high-resolution direct numerical simulations of the Hasegawa–Wakatani equations. We use the Stokes number to quantify the inertia of impurity particles. It is found that particle inertia will cause spatial intermittency in particle distribution and the formation of large-scale structures, i.e. the clustering of particles. The degrees of clustering are influenced by the Stokes number. To quantify these observations, we apply a modified Voronoi tessellation, which assigns specific volumes to impurity particles. By determining time changes of these volumes, we can calculate the impurity velocity divergence, which allows the clustering dynamics to be assessed. To quantify the clustering statistically, several approaches are applied, such as probability density function (PDF) of impurity velocity divergence and joint PDF of volume and divergence.

**Keywords:** fusion plasma, plasma simulation, plasma flows

**1. Introduction**

Fusion reactors aim to achieve fusion reactions by confining and heating a plasma
mixture of ions and electrons at extremely high temperatures. Impurities in fusion plasmas
can profoundly affect performance and stability. The impurities primarily originate from
the interaction between the intense heat of the plasma and the tokamak's walls. Their
presence can lead to significant radiation losses and make it challenging to maintain
the requisite temperatures for fusion reactions (Wesson & Campbell 2011). These
impurities can also dilute the primary fusion fuel, thereby reducing the overall fusion
reaction rate. In addition, strong radiation on the plasma's edge can drastically reduce
electrical conductivity, leading to disruptions in the current profile inside the q = 2
surface (Stangeby 2000). Moreover, before reaching density limits, plasmas often develop
multifaceted asymmetric radiation from the edge (MARFE) which produces localised

† Email address for correspondence: kai.schneider@univ-amu.fr

<footnote n="1"><sup>1</sup>Aix-Marseille Université, CNRS, I2M, UMR 7373, 13331 Marseille, France</footnote>

<footnote n="2"><sup>2</sup>Aix-Marseille Université, CNRS, IUSTI, UMR 7343, 13453 Marseille, France</footnote>

<footnote n="3"><sup>3</sup>Aix-Marseille Université, CNRS, PIIM, UMR 7345, 13397 Marseille, France</footnote>
radiation, usually at the inner wall or the X-point, caused by strong impurity radiation
near the edge (Stangeby 2000).

The behaviour of plasma flow within reactors, specifically in the edge region of a
tokamak, is influenced by drift-wave turbulence and zonal flows. Our high-resolution
simulations to understand this flow are based on the Hasegawa–Wakatani (HW) model
(Hasegawa & Wakatani 1983), which provides insights into cross-field transport caused
by electrostatic drift waves. The electric field perpendicular to magnetic field lines is
particularly significant because it strongly drives cross-field fluxes, and influences edge
pressure profiles and overall stability (Zhang, Guo & Diamond 2020). The HW equations
are particularly useful for studying plasma turbulence, which is a major cause of energy
loss in tokamak devices. These equations capture the essential physics of drift-wave
turbulence. Despite their simplicity, the HW equations have proven to be a powerful tool
for understanding the complex dynamics of plasma behaviour in tokamak devices (Horton
2012).

For the particles, the Lagrangian perspective has been of interest in recent years
(Bos et al. 2010; Kadoch, Bos & Schneider 2010; Kadoch et al. 2022; Gheorghiu, Militello
& Rasmussen 2024). This perspective delves into transport properties by analysing
the trajectories of multiple tracer particles. Through numerical simulations, we solve
equations about these particle trajectories within particular plasma flow velocity fields,
such as the *E* × *B* field. This method highlights the significant effect of coherent structures
on transport. The interplay between eddy trapping and zonal shear flows results in
non-diffusive transport (del Castillo-Negrete, Carreras & Lynch 2004; van Milligen,
Sanchez & Carreras 2004; Krasheninnikov, D’Ippolito & Myra 2008; Krasheninnikov
2024; Krasheninnikov & Smirnov 2024). Both experimental and numerical findings
confirm the presence of these coherent structures in edge turbulence (Krasheninnikov et al.
2008).

In fusion plasma research, existing studies have focused on passive flow tracers in edge
plasma, using the HW model, without taking inertial effects into account (Futatani et al.
2008a,b; Futatani, Benkadda & del Castillo-Negrete 2009). In the work by Priego et al.
(2005), the inertial effect is considered and a fluid model for impurities is used, wherein the
impurity particle velocity is considered to be the sum of the *E* × *B* and polarisation drifts.
The polarisation drift, accounting for impurity particle inertia, is a higher-order correction
to the *E* × *B* drift velocity and introduces compressible effects. In contrast to a fluid
model, we track each impurity particle individually. In addition, we assume that heavy
impurity particles may ‘lag behind’ plasma flow velocity due to the significant inertia.
Indeed, the role of inertia in particle advection by turbulent flows is well-documented in
fluid dynamics and is known to result in clustering around vortical structures (Provenzale
1999). Our study aims to analyse possible clustering effects for impurities in magnetised
plasma.

A key element in understanding these heavy impurities is their tendency for
self-organisation. This self-organisation is seen in clusters of impurities and void regions
(Monchaux, Bourgoin & Cartellier 2010; Lin et al. 2024). It can be mathematically
quantified by deviation from Poissonian statistics (Oujia, Matsuda & Schneider 2020;
Maurel-Oujia, Matsuda & Schneider 2023). The goal is to investigate the statistical
signature of clustering and void regions. Recent results for hydrodynamics turbulence are
promising (Matsuda, Schneider & Yoshimatsu 2021). Finite-time measures to quantify
divergence and the rotation of the particle velocity by determining the volume change rate
of the Voronoi cells and their rotation, respectively, were proposed by Oujia et al. (2020)
and Maurel-Oujia et al. (2023) and applied in the context of hydrodynamic turbulence.
Here we apply them to characterise the dynamics of the self-organisation of heavy
impurities in plasma and assess their clustering and void formation.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Two-dimensional slab geometry. In the tokamak edge, the two-dimensional slab geometry flow configuration is depicted using the HW system. Here, the radial direction is represented by x, whereas y is the poloidal direction. There is an imposed mean plasma density gradient $\nabla n_0$ in the radial direction. The two-dimensional flow is computed within a square domain measuring 64 Larmor radii ($\rho_s$) on each side, as indicated by the green outline. On the right, a vorticity field is displayed.</figcaption></figure>

The paper is structured as follows. In § 2 we outline the theoretical basis used for simulation, tessellation-based analysis method and numerical set-up. Section 3 details our numerical simulation results and tessellation-based analysis results. Finally, § 4 summarises our findings and suggests potential directions for future investigations.

## 2. Models for simulation and analysis method

### 2.1. HW model

The model for our numerical simulations is based on the HW equations for plasma edge turbulence driven by drift-wave instability (Hasegawa & Wakatani 1983). Our focus in this study is on the two-dimensional slab geometry of the HW model; see, e.g., Kadoch et al. (2022). The two-dimensional HW model is a representative paradigm for understanding the nonlinear dynamics of drift-wave turbulence. Figure 1 depicts a representation of the flow configuration.

The magnetic field lines are assumed to be straight and perpendicular to the slab. Ions are considered cold, ignoring temperature gradient effects. The key variables are normalised as detailed in the work done by Bos et al. (2010) and Futatani et al. (2011),

$$x/\rho_s \to x, \quad \omega_{ci}t \to t, \quad e\phi/T_e \to \phi, \quad n_1/n_0 \to n, \qquad (2.1a-d)$$

where $\rho_s$ is the ion Larmor radius at electron temperature $T_e$. It is defined as $\rho_s = \sqrt{T_e/m_i}/\omega_{ci}$, where $\omega_{ci}$ is the ion cyclotron frequency. Here $n_1$ and $n_0$ represent plasma density fluctuation and equilibrium density, whereas $\phi$ indicates electrostatic potential fluctuation. The HW model consists of two equations that describe the evolution of plasma potential and fluctuating plasma density, respectively:

$$\left(\frac{\partial}{\partial t}-\mu_{v} \nabla^{2}\right) \nabla^{2} \phi=\left[\nabla^{2} \phi, \phi\right]+c\left(\phi-n\right), \qquad (2.2)$$

$$\left(\frac{\partial}{\partial t}-\mu_{D} \nabla^{2}\right) n=[n, \phi]-\kappa \frac{\partial \phi}{\partial y}+c\left(\phi-n\right), \qquad (2.3)$$
where $\mu_D$ is cross-field diffusion coefficient and $\mu_v$ is kinematic viscosity. The term $\kappa$, defined as $\kappa = -\partial_x \ln(n_0)$, is a measure of the plasma density gradient. Here $\kappa$ is assumed constant, implying $n_0 \propto \exp(-\kappa x)$. The Poisson bracket is defined as $[A, B] = (\partial A/\partial x)(\partial B/\partial y) - (\partial A/\partial y)(\partial B/\partial x)$. In these equations, the electrostatic potential $\phi$ is the stream function for the $E \times B$ velocity, represented by $u = \nabla^{\perp}\phi$, where $\nabla^{\perp} = (-\partial/\partial y, \partial/\partial x)$. Thus, $u_x = -\partial\phi/\partial y$ and $u_y = \partial\phi/\partial x$ with vorticity $\omega = \nabla^2\phi$.

In the equations, $c$ is the so-called adiabaticity parameter which measures the parallel electron response. It is defined as

$$c = \frac{T_e k_{\parallel}^2}{e^2 n_0 \eta \omega_{ci}}. \qquad (2.4)$$

With $\eta$ as electron resistivity and $k_{\parallel}$ as the effective parallel wavenumber, parameter $c$ determines the phase difference between electrostatic potential and plasma density fluctuations. The model described above is the classical Hasegawa–Wakatani (cHW) model. For $c \gg 1$ (adiabatic limit), the model reduces to the Hasegawa–Mima equation, where electrons follow a Boltzmann distribution. At $c \ll 1$ (hydrodynamic limit), the system reduces to a form that is analogous to the two-dimensional Navier–Stokes equation, where density fluctuations are passively influenced by the $E \times B$ drift. Notably, for $c \sim 1$ (quasi-adiabatic regime), a phase shift between the potentials and densities is observed, enabling particle transport and reflecting a complex interplay between them. To obtain zonal flow, a revised version of the model, known as the modified Hasegawa–Wakatani (mHW) model, can be considered. This modification was introduced in Pushkarev, Bos & Nazarenko (2013) which consists in setting the coupling term $c(\phi - n)$ to zero for modes $k_y = 0$. In this paper, we focus on the quasi-adiabatic regime, i.e. $c = 0.7$ (cHW model), which is relevant to the edge plasma of tokamaks (Bos et al. 2010). For readers interested in the results of other regimes, we refer to Appendix A.

## 2.2. Impurity particles model

Earlier research (Futatani et al. 2008a,b, 2009) assumed that impurity particles, due to no inertia, act as ‘passive tracers’ aligning closely with fluid flow. Yet, for heavier particles, this assumption might not always hold, particularly when the particles have significant mass, leading to notable inertial effects. To address this, the Stokes number ($St$) can be introduced (see, e.g., Oujia et al. 2020), a dimensionless parameter that quantifies particle inertia in fluid flow. The Stokes number is defined as

$$St = \frac{\tau_p}{\tau_\eta}, \qquad (2.5)$$

where $\tau_p$ denotes the impurity particle’s relaxation time, it measures the time the particle takes to adapt to fluid flow alterations. A smaller $\tau_p$ means the particle adapts faster. On the other hand, $\tau_\eta$ is the characteristic timescale for turbulence, indicating how fast the fluid flow changes. A smaller $\tau_\eta$ suggests more rapid turbulence changes. A high Stokes number means that a particle’s movement is primarily driven by inertia, causing it to keep its direction despite fluid changes. In previous studies (Futatani et al. 2008a,b, 2009), the impurity particles are treated as passive tracers of the flow, indicating $\tau_p = 0$. However, here we consider heavy impurity particles with $\tau_p > 0$, for instance tungsten.

The heavy impurity particles are modelled as charged point particles. They are considered as test particles, which means that there is no impurity–impurity interaction and impurities have no effect on the plasma dynamics, whereas the plasma flow velocity affects the motion of impurities (one-way coupling). Impurities in our model are assumed
to be cold, just as with the treatment of ions in bulk plasma. Thus, thermal motion of impurities could be neglected for simplicity. The impurity particles experience not only the Lorentz force due to the electric and magnetic fields, but also drag force resulting from momentum transfer with plasma ions. Positive plasma ions transfer momentum to heavy charged impurity particles by interacting with electrostatic potential in the vicinity of the impurity particle. The transfer of momentum from electrons to impurities is negligible due to their small mass. Coulomb scattering theory provides, macroscopically, an approximation of a drag force $F_{drag}$ dependent on the relative velocity between the macroscopic velocity of the plasma ion and the velocity of the impurity particle (Kilgore et al. 1993): $F_{drag} = K_{mt}n_i m_i v_r$, where $v_r$ is the relative velocity between the macroscopic velocity of the plasma ions and the impurity particles, $K_{mt}$ is the momentum transfer coefficient, $n_i$ is the plasma ion density and $m_i$ is the ion mass of the plasma. This drag force tends to reduce the velocity difference, causing the impurity particles to relax to the plasma flow velocity. The relaxation time of the impurity particle $\tau_p$ is the period during which impurity particles adjust to the plasma flow due to the drag force. In our model, we approximate this drag force as $m_p(u_p - v_p)/\tau_p$, which is proportional to the relative velocity between the plasma flow and impurities. The relaxation time can be expressed as $\tau_p = m_p/K_{mt}n_i m_i$. This relation highlights that $\tau_p$ is inversely proportional to the plasma ion density and momentum transfer coefficient, while directly proportional to the impurity mass. Given this relationship, as the ion density $n_i$ increases, the relaxation time decreases, leading to faster coupling of the impurity particles to the plasma flow. Conversely, larger impurity masses result in longer relaxation times, reflecting the slower response of heavier particles to the plasma flow. According to Newton's second law:

$$m_p \frac{d\mathbf{v}_{p,j}}{dt} = \frac{m_p(\mathbf{u}_{p,j} - \mathbf{v}_{p,j})}{\tau_p} + Ze(E(\mathbf{x}_{p,j}) + \mathbf{v}_{p,j} \times \mathbf{B}), \quad (2.6)$$

$$E(\mathbf{x}_{p,j}) = -\nabla\phi(\mathbf{x}_{p,j}). \qquad (2.7)$$

Equation (2.6) describes the forces exerted on each heavy impurity particle, denoted by $j$ ($j = 1, \dots, N_p$, with $N_p$ being the total number of particles). Here $Z$ is charge state of the impurity particle and $e$ is the elementary charge. The variable $\mathbf{v}_{p,j}$ is the $j$th particle's velocity, whereas $\mathbf{u}_{p,j}$, $E(\mathbf{x}_{p,j})$ and $\nabla\phi(\mathbf{x}_{p,j})$ represent the fluid velocity, electric field and gradient of potential at the particle's location $\mathbf{x}_{p,j}$, respectively. Let us note that when $\tau_p = 0$, we recover the situation studied in previous works (Futatani et al. 2008a,b, 2009). In this case, the impurity particles are considered as passive tracers of the flow, i.e. $\mathbf{v}_{p,j} = \mathbf{u}_{p,j}$. Applying the same normalisation as used for the HW equation, we obtain

$$\frac{d\mathbf{v}_{p,j}}{dt} = \frac{(\mathbf{u}_{p,j} - \mathbf{v}_{p,j})}{\tau_p} + \alpha(-\nabla\phi(\mathbf{x}_{p,j}) + \mathbf{v}_{p,j} \times \mathbf{b}), \quad (2.8)$$

where $\alpha = Zm_i/m_p$, $m_i$ is the ion mass of the plasma and $\mathbf{b}$ represents the unit vector along the direction of the magnetic field. Combining with the equation:

$$\frac{d\mathbf{x}_{p,j}}{dt} = \mathbf{v}_{p,j}, \qquad (2.9)$$

it is thus possible to compute the trajectory of the impurity particles.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. (a) Electric potential $\phi$ (stream function) and $10^4$ superimposed impurity particles for $St = 1$ in the case $c = 0.7$ (cHW model). (b) A magnified view with modified Voronoi tessellation. Large cells correspond to void regions, small cells represent clusters.</figcaption></figure>

## 2.3. Modified Voronoi tessellation

Impurity particles satisfy the conservation equation:

$$D_t n_p = -n_p \nabla \cdot \mathbf{v}_p, \qquad (2.10)$$

where $D_t$ is the Lagrangian derivative, $n_p$ is the impurity particle number density and $\mathbf{v}_p$ the particle velocity. From this equation, the divergence of particle velocity $\nabla \cdot \mathbf{v}_p = -D_t n_p/n$ is crucial to understand. A positive divergence ($\nabla \cdot \mathbf{v}_p > 0$) leads to density decrease, whereas a negative divergence, i.e. convergence ($\nabla \cdot \mathbf{v}_p < 0$), causes an increase in density. To quantify the divergence, we employ the modified version of the classical Voronoi tessellation technique, which uses the centre of gravity, instead of the circumcentre of the Delaunay cell, to define Voronoi cell vertices. The purpose of this method is to assign a volume to each impurity particle, enabling the subsequent calculation of the divergence. More details about the classical and modified Voronoi tessellation and their difference can be found in Oujia et al. (2020) and Maurel-Oujia et al. (2023). The modified Voronoi tessellation divides a plane into regions according to the positions of impurity particles, which allows us to assign a cell, called a 'modified Voronoi cell', and thus a corresponding volume $V$ to each particle. Each cell can be considered as the region of influence of an impurity particle. In our two-dimensional scenario, the volume of the cell is determined by calculating its area. Figure 2 shows the modified Voronoi tessellation. In the figure, clustered particles are represented by small cells and void regions are represented by large cells.

To compute the divergence of impurity particle velocity $\mathcal{D}$ we analyse particle distributions at two time instants $t^k$ and $t^{k+1} = t^k + \Delta t$, where $\Delta t$ is the time step of the simulation and $k$ is the discrete time index. Each particle distribution snapshot is assessed using modified Voronoi tessellation. The inverse of volume, $1/V$, approximates the local particle number density $n_p$ in the discrete setting. It is proved that (Oujia et al. 2020;
<table><caption>TABLE 1. Simulation parameters for the flow: A, domain area; R, grid resolution; Δt, time step; N<sub>p</sub>, number of impurity particles; μ<sub>D</sub>, diffusion coefficient; μ<sub>ν</sub>, kinematic viscosity; κ ≡ −∂<sub>x</sub> ln(n<sub>0</sub>), a measure of the plasma density gradient; c, adiabaticity parameter.</caption><thead><tr><th>A</th><th>R</th><th>Δt</th><th>N<sub>p</sub></th><th>μ<sub>D</sub></th><th>μ<sub>ν</sub></th><th>κ</th><th>c</th></tr></thead><tbody><tr><td>64 × 64</td><td>1024 × 1024</td><td>5 × 10<sup>-4</sup></td><td>10<sup>6</sup></td><td>5 × 10<sup>-3</sup></td><td>5 × 10<sup>-3</sup></td><td>1</td><td>0.7</td></tr></tbody></table>

<table><caption>TABLE 2. Stokes numbers used in the simulation.</caption><tr><td>St</td><td>0</td><td>0.05</td><td>0.25</td><td>0.5</td><td>1</td><td>2</td><td>5</td><td>50</td></tr></table>

Maurel-Oujia et al. 2023)

$$ \mathcal{D}(\mathbf{v}_p) = -\frac{1}{n_p} D_t n_p = \frac{2 V^{k+1} - V^k}{\Delta t V^{k+1} + V^k} + O\left(\frac{1}{N_p}, \Delta t\right). \quad (2.11) $$

Thus, by measuring the volume change of each particle between $t^k$ and $t^{k+1}$, we can determine the discrete divergence of its velocity $\mathcal{D}(\mathbf{v}_p)$. The velocity divergence of impurity particles measures the rate of volume change for a particle group over time, indicating how much the particle velocity field is spreading out or converging at a particular point in space.

## 2.4. Numerical set-up

The simulation was conducted within a domain that spans an area A of 64 × 64 with periodic boundary condition. The domain was discretised to a resolution R of 1024 × 1024 grid points. The time step Δt was 5 × 10<sup>-4</sup>. The number of impurity particles N<sub>p</sub> was 10<sup>6</sup>. The adiabaticity parameter c is 0.7. The values of specific physical parameters are listed in table 1. The characteristic time scale of turbulence, τ<sub>η</sub>, is defined as 1/√2Z<sub>ms</sub>, with Z<sub>ms</sub> denoting one half of the mean-square vorticity. In the simulation, τ<sub>η</sub> = 0.35 in the statistically steady flow. The simulations are based on the work of Kadoch et al. (2022) where the flow is initialised with Gaussian random fields. Here we start the simulations with a flow already being in the statistically steady state (i.e. after the long transition phase of drift-wave instabilities; see Kadoch et al. 2022) and one million uniformly distributed heavy impurity particles with random velocity are injected.

For the impurity particles, we focus on tungsten <sup>184</sup>W<sup>20+</sup>, a heavy metal that can be made from the material that is exposed to the high heat flux in ASDEX-Upgrade (Krieger, Maier & Neu 1999) and EAST (Yao et al. 2016), expecting the same situation for ITER (Pitts et al. 2009). For <sup>184</sup>W<sup>20+</sup>, α = Zm<sub>i</sub>/m<sub>p</sub> = 0.22. As for Stokes number St = τ<sub>p</sub>/τ<sub>η</sub>, τ<sub>η</sub> = 0.35 in the statistically steady state of flow for c = 0.7 (cHW model), but the exact value of τ<sub>p</sub> for <sup>184</sup>W<sup>20+</sup> is unknown. To resolve this, we explore a range of τ<sub>p</sub> values spanning several orders of magnitude, resulting in Stokes numbers St = 0, 0.05, 0.5, 1, 5, 10 and 50. This parametric approach allows us to systematically explore the dynamics of tungsten impurity, from particles that are tightly coupled to the flow (St ≪ 1) to those
that are essentially independent of the flow (St ≫ 1). The simulation parameters are listed in table 2.

The system of HW equations, specifically (2.2) and (2.3) are solved using a classical pseudospectral method (Canuto et al. 2007). The equations governing the motion of the impurity particles, namely (2.8) and (2.9), are solved using a second-order Runge–Kutta (RK2) scheme and linear interpolation is used for computing the fluid velocity at the particle position. The code of this study builds upon that used in the research done by Kadoch et al. (2022).

## 3. Results

### 3.1. Simulation results

The analysed data set consists of the electric potential $\phi$ (stream function), particle position and velocity data generated by direct numerical simulation (DNS). The contour curves of constant $\phi$ represent the streamlines of the flow field and provide a clear visualisation of fluid flow in two dimensions. Closed contour curves indicate the presence of vortices. Figure 3 illustrates the $\phi$ fields and the behaviour of impurity particles in statistically steady state within the quasi-adiabatic regime ($c = 0.7$, cHW model), where impurity particles are considered as passive tracers without considering inertial effect (St = 0). As shown by figure 3 impurity particles are distributed randomly. Figure 4 shows that with increasing Stokes numbers (St = 0.05, 0.5, 1), particle inertia begins to dominate, causing them to deviate from the fluid streamlines due to Coriolis force induced by the vortices. This deviation causes them to concentrate where the vorticity is minimal, typically at vortex peripheries. For even higher Stokes numbers (St = 5, 10, 50), the impurity particles tend to maintain their trajectories, moving more ballistically and exhibiting reduced clustering. From (2.8), when the coefficient of the flow effect term, $1/\tau_p = 1/(\tau_\eta St)$, is much larger than the coefficient term of the Lorentz force $\alpha$, the particle dynamics are primarily governed by the flow, whereas the Lorentz force plays a relatively minor role. For $St = 0.05, 0.5$ and $1, 1/\tau_p = 1/(0.35St) \gg \alpha = 0.22$, the particles are more strongly influenced by the fluid flow than by the electromagnetic forces. However, as the Stokes number increases to values such as $St = 10$ and 50, the particles are less coupled to the fluid flow and the electromagnetic forces start to play an important role. As shown in figure 4, for case $St = 10$ and 50, the positively charged impurity particles are mostly located in negative potential regions (red region), as expected. The effects of different $\alpha$ values are further explored in Appendix B.

Let us note that in regions of impurity clustering, where the impurity density may locally increase, there is a risk that clustering could lead to deviations from the assumption of quasi-neutrality that is underlying the derivations of the HW model (Naulin et al. 2006). This strict condition, which requires that the impurity density remains smaller than the bulk current divergence on the characteristic drift-wave time scale, cannot always be guaranteed in regions of intense clustering. Nevertheless, as noted in Naulin et al. (2006), considering impurities as test particles is still valuable for understanding tendencies in impurity transport within a given flow field, and the clustering effects observed in our simulations provide insights into impurity dynamics. We further acknowledge that resolving this discrepancy would require a more sophisticated model that accounts for the back-reaction of impurity densities on the plasma fields and their contribution to quasi-neutrality. Such a model would need to self-consistently address the plasma profile problem in the presence of impurities, which is beyond the scope of the present work, but remains a valuable direction for future research.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Electric potential field $\phi$ (stream function) superimposing $10^4$ impurity particles (out of $10^6$) for $St = 0$ in statistically steady state within the quasi-adiabatic regime ($c = 0.7$, cHW model).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Electric potential fields $\phi$ (stream function) superimposing $10^4$ impurity particles (out of $10^6$) for various Stokes numbers in statistically steady state within the quasi-adiabatic regime ($c = 0.7$, cHW model).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. (a) PDF of volume normalised by the mean for different Stokes numbers against the gamma distribution. (b) PDF of the divergence of impurity particle velocity for different Stokes numbers for $^{184}$W$^{20+}$ ($\alpha$ = 0.22) in quasi-adiabatic regime ($c$ = 0.7, cHW model).</figcaption></figure>

## 3.2. Tessellation-based analysis results

We calculate the volume $V$ for each impurity particle, and then plot the probability density function (PDF) of $V/\bar{V}$, where $\bar{V}$ is the average volume, calculated as $\bar{V} = 64^2/N_p$ ($64^2$ is the domain area, $N_p$ is the total number of particles). Figure 5(a) shows the PDF of $V/\bar{V}$.

A random variable $X$ that is gamma-distributed with shape $k$ and scale parameter $\theta$ is denoted $X \sim \Gamma(k, \theta)$. The corresponding PDF is $f(x) = \Gamma(k)^{-1}\theta^{-k}x^{k-1}\exp(-x/\theta)$. For particles that are distributed randomly, the PDF of the volumes follows a gamma distribution (Ferenc & Néda 2007). For $St = 0$, using the maximum-likelihood PDF estimation for fitting, we find that the PDF of the volumes aligns well with the gamma distribution, specifically we have $\Gamma(1.72, 0.58)$. This suggests, as expected, that the impurity particles are distributed randomly when $St = 0$, indicating no clustering of particles. Let us note that the PDF of the volumes for $St = 0$ obtained using classical Voronoi tessellation follows a gamma distribution $\Gamma(7/2, 2/7)$ (Ferenc & Néda 2007). Since we use a different method, the modified Voronoi tessellation, the fitting gives different values, i.e. $\Gamma(1.72, 0.58)$. The PDF of the volume is typically used to classify ‘cluster cells’ and ‘void cells’ (Monchaux et al. 2010). A cell below a certain threshold is considered a cluster cell, whereas a cell above this threshold is categorised as a void cell. In our study, the threshold is $V/\bar{V} \sim 0.5$. As shown by figure 5(a), when $St$ increases, we observe an increase in the number of cluster cells (with $V/\bar{V} < 0.5$). However, once $St$ exceeds $St = 1$, the number of cluster cells begins to decrease. This is also highlighted by figure 4 which shows that the particles are densely packed for $St = 1$.

The PDF of the divergence of impurity particle velocity is displayed in figure 5(b). As shown in figure 5(b), the divergence of the fluid velocity (represented by the divergence of the impurity particle velocity with $St = 0$) is small but not zero. However, as we know, in an incompressible continuous fluid, the divergence of the fluid velocity should be zero. The reason for this discrepancy lies in our computational method, which segments the fluid into separate pieces known as ‘modified Voronoi cells’ as mentioned before. The deformation of a modified Voronoi cell does not correspond to the deformation of a fluid volume in a continuous context. The figure implies that the divergence is attributable to discretisation errors. A negative divergence in a particle velocity field indicates that particles converge towards a point, thus increasing density locally. Conversely, positive divergence shows that particles spread out from a point, reducing the local density. As shown in figure 5(b), for increasing Stokes number, the PDF of divergence $D$ widens, reflecting stronger particles convergence and divergence activities. The symmetry observed in figure 5(b) suggests
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Joint PDF of the volume in log scale and divergence in linear scale for (a) $St = 0.05$, (b) $St = 0.5$, (c) $St = 1$, (d) $St = 5$, (e) $St = 10$ and (f) $St = 50$ in the quasi-adiabatic regime ($c = 0.7$, cHW model).</figcaption></figure>

that positive and negative divergence values are comparable. This balance is due to the conservation of impurity particles.

The joint PDFs of the divergence $D$ and normalised volume $V/V$ show more insight into cluster formation (cf. Figure 6). There is a high probability along $D = 0$, indicating groups of particles move as an incompressible flow. High negative/positive divergence probability in cluster cells ($V/V < 0.5$) may be due to the fact that the probability of finding particles in clusters is higher than in void region as illustrated by figure 5(a). As the Stokes number increases, particles deviate more from the fluid flow, increasing their relative velocity. Therefore, high negative/positive divergence could also be from 'caustics' where the multi-valued particle velocities at a given position cause large velocity differences between neighbouring particles, leading to extreme divergence (Wilkinson & Mehlig 2005). Figure 7 visually supports this by colouring particles based on their divergence $D$. From figure 7, we observe that particles with large negative/positive divergence values are densely packed in clusters, indicating simultaneous particles convergence and divergence.

4. Conclusions

In this work, we have studied the spatial distribution of impurity particles in the edge plasma of tokamaks. High-resolution simulations of the plasma flow were based on the paradigmatic model in edge plasma of tokamaks, namely the HW model, whereas the impurity particles were tracked using a Lagrangian method. Assuming that impurity particles have a non-zero relaxation time ($\tau_p \neq 0$), which means that the particles do
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Spatial distribution of one million particles coloured with the divergence $\mathcal{D}$ for (a) $St = 0.05$, (b) $St = 0.5$, (c) $St = 1$, (d) $St = 5$, (e) $St = 10$ and (f) $St = 50$ in the quasi-adiabatic regime ($c = 0.7$, cHW model).</figcaption></figure>

not immediately adapt to fluctuations in plasma flow, we observed the phenomenon of impurity clustering within the plasma. To gain deeper insights into the self-organisation of impurity particles, we employed a modified Voronoi tessellation technique to compute the divergence of impurity particle velocity $\mathcal{D}$. The main results are as follows.

### (i) The effect of $St$

* As $St$ increases, impurities cluster in low-vorticity regions due to their inability to adapt to rapid changes in plasma flow. At high $St$, impurities move randomly with less clustering.

* With an increase in $St$, there is a higher probability of observing higher values of negative or positive divergence (for $\alpha = 0.22$). This implies that the impurity particle velocity field is converging or diverging faster. These large values of negative or positive divergence are more likely to be found in cluster cells where ($V/V < 0.5$).

### (ii) The effect of $\alpha$

* The influence of $\alpha$ becomes significant only when $St$ is large, which tends to drive the impurity particles into negative electric potential regions.
For future work, it could be interesting to extend the analysis to multi-resolution techniques that would provide a better understanding of the clustering behaviour of impurity particles at different scales; see, e.g., Matsuda *et al.* (2022). This approach could reveal more details about the spatial distribution and dynamics of clusters that are not apparent at a single scale of observation. Another interesting direction is to explore the transport dynamics of heavy impurity particles for different Stokes number by tracking them for a long time. Finally, the limitation of this work is that we assumed that impurities do not affect plasma flow. Although this simplification allowed for a focused study on impurity particle dynamics, it may not capture the full complexity of the interaction between impurities and plasma flow. A possible extension in future studies is to explore the effect of impurity particles on the plasma flow based on the work of Benkadda *et al.* (1996).

## Acknowledgements

Z.L. acknowledges P. Krah for the insightful discussions. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for providing access to its high performance computing resources. The authors thank the referees for constructive comments.

Editor Paolo Ricci thanks the referees for their advice in evaluating this article.

## Declaration of interest

The authors report no conflict of interest.

## Funding

This work was supported by I2M (Z.L., T.M.O. and K.S.); the French Federation for Magnetic Fusion Studies (FR-FCM) and the Eurofusion Consortium, funded by the Euratom Research and Training Programme (Z.L., T.M.O., K.S. and S.B., grant number 633053); and the Agence Nationale de la Recherche (ANR), project CM2E (Z.L., T.M.O., B.K. and K.S., grant number ANR-20-CE46-0010-01). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

# Appendix A. The behaviour of $^{184}_{74}\text{W}^{20+}$ in different flow regimes

## A.1. *Hydrodynamic and adiabatic regime*

We explored the behaviour of $^{184}_{74}\text{W}^{20+}$ in hydrodynamic ($c = 0.01$, cHW model) and adiabatic regime ($c = 2$, cHW model). Figures 8 and 9 show potential fields and impurity behaviour in hydrodynamic and adiabatic regimes. We have the same observation as in the quasi-adiabatic case ($c = 0.7$, cHW model): at $St = 0$, impurities act as passive tracers, uniformly distributed along fluid streamlines; as $St = 0.05$, particles start to cluster in areas of low vorticity and clustering is clear at $St = 1$ ; at $St = 50$, there is less clustering and particles are mostly presented in negative potential areas. Similar to the quasi-adiabatic case, as shown in figures 10(a) and 11(a), as $St$ increases, the number of cluster cells ($V/V < 0.5$) initially increases until $St = 1$. Figures 10(b) and 11(b) demonstrate that the PDF of divergence $D$ broadens with increasing Stokes number.

## A.2. *Zonal flows*

In the mHW model with $c = 2$, zonal flows are present. Figure 12 shows the zonal flows and superimposing impurity particles. Although it may not be clearly visible that particles with $St = 1$ are primarily presented in regions of low vorticity, as indicated by figure 12, it is evident in the vorticity field (figure 13). For $St = 50$, particles are predominantly
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Electric potential fields $\phi$ (stream function) superimposing $10^4$ impurity particles (out of $10^6$) for various Stokes numbers in statistically steady state within the hydrodynamic regime ($c = 0.01$, cHW model).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 9. Electric potential fields $\phi$ (stream function) superimposing $10^4$ impurity particles (out of $10^6$) for various Stokes numbers in statistically steady state within the adiabatic regime ($c = 2$, cHW model).</figcaption></figure>

found in regions of negative potential. As shown in figure 14(a), there are fewer cluster cells ($V/V < 0.5$) compared with other regimes where zonal flows are absent ($c = 0.01, 0.7$ and 2, cHW model). Figure 14(b) illustrates that, in the presence of zonal flows, the likelihood of large negative or positive divergence is reduced compared with other regimes, suggesting that zonal flows tend to prevent the convergence and divergence of particles.

## Appendix B. Analysis of tungsten impurity with different charge state in quasi-adiabatic regime ($c = 0.7$, cHW model)

The clustering of two other tungsten impurities with different charge state is analysed in this appendix. We consider ${}^{184}_{74}\text{W}^{1+}$ ($\alpha = 0.01$) and ${}^{184}_{74}\text{W}^{44+}$ ($\alpha = 0.48$) (Maget et al. 2020). Applying modified Voronoi tessellation, the PDF of volume normalised by the mean ($V/V$) and the PDF of impurity velocity divergence $D$ for different Stokes numbers are plotted in figures 15 and 16.

Figures 15(a) and 16(a) show that for both ${}^{184}_{74}\text{W}^{1+}$ and ${}^{184}_{74}\text{W}^{44+}$, at $St = 0$, the impurity particles exhibit a random distribution, as evidenced by the alignment of the curve with the gamma distribution. This indicates the absence of inertial effects. As $St$ increases, the number of cluster cells ($V/V < 0.5$) increases and then decreases after exceeding $St = 1$. From figure 15(b), for ${}^{184}_{74}\text{W}^{1+}$, it is observed that the PDFs of divergence broaden with
<figure><img src="image_12.png" /></figure>

FIGURE 10. (a) PDF of volume normalised by the mean (V/V) and (b) PDF of impurity velocity divergence $\mathcal{D}$ for different Stokes numbers for $^{184}\text{W}^{20+}$ ($\alpha = 0.22$) in hydrodynamic regime ($c = 0.01$, cHW model).

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

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

FIGURE 11. (a) PDF of volume normalised by the mean (V/V) and (b) PDF of impurity velocity divergence $\mathcal{D}$ for different Stokes numbers for $^{184}\text{W}^{20+}$ ($\alpha = 0.22$) in hydrodynamic regime ($c = 2$, cHW model).

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

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

FIGURE 12. Electric potential fields $\phi$ superimposing $10^4$ impurity particles (out of $10^6$) for various Stokes numbers in statistically steady state within zonal flows ($c = 2$, mHW model).

increasing $St$, then the PDFs begin to narrow from $St = 10$ to $St = 50$. From figure 15(b), for $^{184}\text{W}^{44+}$, it is observed that the PDFs of divergence widen as $St$ increases.

Tables 3 and 4 provide a quantitative comparison of how the divergence characteristics change with varying $\alpha$ values, which are indicative of the charge state of the tungsten
<figure><img src="image_17.png" /><figcaption>FIGURE 13. Vorticity field superimposing $10^4$ impurity particles (out of $10^6$) for Stokes number $St = 1$ in statistically steady state within zonal flows ($c = 2$, mHW model).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 14. (a) PDF of volume normalised by the mean (V/V) and (b) PDF of the impurity velocity divergence $\mathcal{D}$ for different Stokes numbers for $^{184}\text{W}^{20+}$ ($\alpha = 0.22$) in zonal flows ($c = 2$, mHW model).</figcaption></figure>

ions ($\alpha = Zm_i/m_p$). Given that negative and positive divergence values are comparable, the overall average divergence approaches zero. Thus, we calculate the mean of positive divergence $\overline{D}_+$, which, by implication, mirrors the behaviour observed for negative divergence. For low Stokes numbers ($St = 0.05, 0.5, 1$), the mean of positive divergence ($\overline{D}_+$) remains relatively consistent regardless of the ionisation state, suggesting that the charge state has a minimal effect on the divergence behaviour of the impurities in this range. However, as we examine higher Stokes numbers ($St = 5, 10$ and $50$), a clear trend emerges: both the variance of divergence ($\sigma_D^2$) and the mean of positive divergence ($\overline{D}_+$) increase with higher $\alpha$ values. This observation implies that at higher $St$ values, the divergence behaviour of the impurities is influenced by their charge state.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 15. (a) PDF of volume normalised by the mean (V/V) and (b) PDF of impurity velocity divergence D for different Stokes numbers for $^{184}\text{W}^{1+}(\alpha = 0.01)$ in the quasi-adiabatic regime (c = 0.7, cHW model).</figcaption></figure>

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

<figure><img src="image_22.png" /><figcaption>FIGURE 16. (a) PDF of volume normalised by the mean (V/V) and (b) PDF of the impurity velocity divergence D for different Stokes numbers for $^{184}\text{W}^{44+}(\alpha = 0.48)$ in the quasi-adiabatic regime (c = 0.7, cHW model).</figcaption></figure>

<table><caption>TABLE 3. Mean of positive divergence $\overline{D}_+$ for different Stokes numbers and different $\alpha$ values.</caption><thead><tr><th>St</th><th>0.05</th><th>0.5</th><th>1</th><th>5</th><th>10</th><th>50</th></tr></thead><tbody><tr><td><span style="text-decoration: overline;">D</span><sub>+</sub>(α = 0.01)</td><td>0.06</td><td>0.4</td><td>1.5</td><td>9.2</td><td>15.1</td><td>9.0</td></tr><tr><td><span style="text-decoration: overline;">D</span><sub>+</sub>(α = 0.22)</td><td>0.06</td><td>0.4</td><td>1.5</td><td>9.5</td><td>17.4</td><td>39.8</td></tr><tr><td><span style="text-decoration: overline;">D</span><sub>+</sub>(α = 0.48)</td><td>0.06</td><td>0.4</td><td>1.4</td><td>10.9</td><td>24.69</td><td>63.0</td></tr></tbody></table>

<table><caption>TABLE 4. Variance of divergence $\sigma_D^2$ for different Stokes numbers and different $\alpha$ values.</caption><thead><tr><th>St</th><th>0.05</th><th>0.5</th><th>1</th><th>5</th><th>10</th><th>50</th></tr></thead><tbody><tr><td>σ<sup>2</sup><sub>D</sub>(α = 0.01)</td><td>0.009</td><td>1.1</td><td>10.7</td><td>1106</td><td>3153</td><td>1087</td></tr><tr><td>σ<sup>2</sup><sub>D</sub>(α = 0.22)</td><td>0.009</td><td>1.0</td><td>10.6</td><td>1214</td><td>4126</td><td>15 902</td></tr><tr><td>σ<sup>2</sup><sub>D</sub>(α = 0.48)</td><td>0.009</td><td>1.0</td><td>9.3</td><td>1632</td><td>7529</td><td>31 319</td></tr></tbody></table>
REFERENCES

BENKADDA, S., GABBAI, P., TSYTOVICH, V.N. & VERGA, A. 1996 Nonlinearities and instabilities of dissipative drift waves in dusty plasmas. *Phys. Rev. E* **53** (3), 2717.

BOS, W.J.T., KADOCH, B., NEFFAA, S. & SCHNEIDER, K. 2010 Lagrangian dynamics of drift-wave turbulence. *Phys. D* **239**, 1269–1277.

CANUTO, C., HUSSAINI, Y.M., QUARTERONI, A. & ZANG, T.A. 2007 *Spectral Methods: Fundamentals in Single Domains*. Springer Science & Business Media.

DEL CASTILLO-NEGRETE, D., CARRERAS, B.A. & LYNCH, V.E. 2004 Fractional diffusion in plasma turbulence. *Phys. Plasmas* **11**, 3854–3864.

FERENC, J.S. & NÉDA, Z. 2007 On the size distribution of Poisson Voronoi cells. *Phys. A* **385** (2), 518–526.

FUTATANI, S., BENKADDA, S. & DEL CASTILLO-NEGRETE, D. 2009 Spatiotemporal multiscaling analysis of impurity transport in plasma turbulence using proper orthogonal decomposition. *Phys. Plasmas* **16** (4), 042506.

FUTATANI, S., BENKADDA, S., NAKAMURA, Y. & KONDO, K. 2008a Multiscaling dynamics of impurity transport in drift-wave turbulence. *Phys. Rev. Lett.* **100** (2), 025005.

FUTATANI, S., BENKADDA, S., NAKAMURA, Y., KONDO, K. & HAMAGUCHI, S. 2008b Anomalous scaling of impurity transport in drift wave turbulence. *Contrib. Plasma Phys.* **48** (1-3), 111–115.

FUTATANI, S., BOS, W.J.T., DEL CASTILLO-NEGRETE, D., SCHNEIDER, K., BENKADDA, S. & FARGE, M. 2011 Coherent vorticity extraction in resistive drift-wave turbulence: comparison of orthogonal wavelets versus proper orthogonal decomposition. *C. R. Phys.* **12**, 123–131.

GHEORGHIU, T., MILITELLO, F. & RASMUSSEN, J.J. 2024 On the transport of tracer particles in two-dimensional plasma edge turbulence. *Phys. Plasmas* **31** (1), 013901.

HASEGAWA, A. & WAKATANI, M. 1983 Plasma edge turbulence. *Phys. Rev. Lett.* **50**, 682–686.

HORTON, W. 2012 *Turbulent Transport in Magnetized Plasmas*. World Scientific.

KADOCH, B., BOS, W.J. & SCHNEIDER, K. 2010 Origin of Lagrangian intermittency in drift-wave turbulence. *Phys. Rev. Lett.* **105** (14), 145001.

KADOCH, B., DEL CASTILLO-NEGRETE, D., BOS, W.J.T. & SCHNEIDER, K. 2022 Lagrangian conditional statistics and flow topology in edge plasma turbulence. *Phys. Plasmas* **29**, 102301.

KILGORE, M.D., DAUGHERTY, J.E., PORTEOUS, R.K. & GRAVES, D.B. 1993 Ion drag on an isolated particulate in a low-pressure discharge. *J. Appl. Phys.* **73** (11), 7195–7202.

KRASHENINNIKOV, S.I. 2024 On anomalous transport of multi-species plasma associated with the resistive ballooning and resistive drift waves driven turbulence. *Phys. Plasmas* **31** (5), 054502.

KRASHENINNIKOV, S.I., D'IPPOLITO, D.A. & MYRA, J.R. 2008 Recent theoretical progress in understanding coherent structures in edge and SOL turbulence. *J. Plasma Phys.* **74**, 679–717.

KRASHENINNIKOV, S.I. & SMIRNOV, R.D. 2024 Anomalous transport of multi-species edge plasma with the generalized Hasegawa–Wakatani model and the FLR effects. *Phys. Plasmas* **31** (5), 052302.

KRIEGER, K., MAIER, H. & NEU, R. 1999 Conclusions about the use of tungsten in the divertor of ASDEX Upgrade. *J. Nucl. Mater.* **266-269**, 207–216.

LIN, Z., MAUREL-OUJIA, T., KADOCH, B., KRAH, P., SAURA, N., BENKADDA, S. & SCHNEIDER, K. 2024 Synthesizing impurity clustering in the edge plasma of tokamaks using neural networks. *Phys. Plasmas* **31** (3), 032505.

MAGET, P., et al. 2020 Natural poloidal asymmetry and neoclassical transport of impurities in tokamak plasmas. *Plasma Phys. Control. Fusion* **62**, 025001.

MATSUDA, K., SCHNEIDER, K., OUJIA, T., WEST, J., JAIN, S. & MAEDA, K. 2022 Multiresolution analysis of inertial particle tessellations for clustering dynamics. In *Center for Turbulence Research, Proceedings of the Summer Program* 2022, pp. 143–152. Stanford University.

MATSUDA, K., SCHNEIDER, K. & YOSHIMATSU, K. 2021 Scale-dependent statistics of inertial particle distribution in high Reynolds number turbulence. *Phys. Rev. Fluids* **6**, 064304.

MAUREL-OUJIA, T., MATSUDA, K. & SCHNEIDER, K. 2023 Computing differential operators of the particle velocity in moving particle clouds using tessellations. *J. Comput. Phys.* **498**, 112658.

MONCHAUX, R., BOURGOIN, M. & CARTELLIER, A. 2010 Preferential concentration of heavy particles: a Voronoi analysis. *Phys. Fluids* **22** (10), 103304.
NAULIN, V., GARCIA, O.E., PRIEGO, M. & RASMUSSEN, J.J. 2006 The application of passive tracers for investigating transport in plasma turbulence. *Phys. Scr.* **T122**, 129-134.

OUJIA, T., MATSUDA, K. & SCHNEIDER, K. 2020 Divergence and convergence of inertial particles in high Reynolds number turbulence. *J. Fluid Mech.* **905**, A14.

PITTS, R.A., *et al.* 2009 Status and physics basis of the ITER divertor. *Phys. Scr.* **T138**, 014001.

PRIEGO, M., GARCIA, O.E., NAULIN, V. & RASMUSSEN, J.J. 2005 Anomalous diffusion, clustering, and pinch of impurities in plasma edge turbulence. *Phys. Plasmas* **12** (6), 062312.

PROVENZALE, A. 1999 Transport by coherent barotropic vortices. *Annu. Rev. Fluid Mech.* **31** (1), 55-93.

PUSHKAREV, A.V., BOS, W.J.T. & NAZARENKO, S.V. 2013 Zonal flow generation and its feedback on turbulence production in drift wave turbulence. *Phys. Plasmas* **20**, 042304.

STANGEBY, P.C. 2000 *The Plasma Boundary of Magnetic Fusion Devices*. CRC Press.

VAN MILLIGEN, B.P., SANCHEZ, R. & CARRERAS, B.A. 2004 Probabilistic finite-size transport models for fusion: anomalous transport and scaling laws. *Phys. Plasmas* **11**, 2272-2285.

WESSON, J. & CAMPBELL, D.J. 2011 *Tokamaks*. Oxford University Press.

WILKINSON, M. & MEHLIG, B. 2005 Caustics in turbulent aerosols. *Europhys. Lett.* **71** (2), 186-192.

YAO, D.M., *et al.* 2016 Design, R&D and commissioning of EAST tungsten divertor. *Phys. Scr.* **T167**, 014003.

ZHANG, Y., GUO, Z.B. & DIAMOND, P.H. 2020 Curvature of radial electric field aggravates edge magnetohydrodynamics mode in toroidally confined plasmas. *Phys. Rev. Lett.* **125** (25), 255003.