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 Reference Hasegawa and Wakatani1983). Our focus in this study is on the two-dimensional slab geometry of the HW model; see, e.g., Kadoch et al. (Reference Kadoch, del Castillo-Negrete, Bos and Schneider2022). 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.

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 $\boldsymbol {\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.

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. (Reference Bos, Kadoch, Neffaa and Schneider2010) and Futatani et al. (Reference Futatani, Bos, del Castillo-Negrete, Schneider, Benkadda and Farge2011),

(2.1a)\begin{equation} x/\rho_{{s}} \rightarrow x,\quad \omega_{\mathrm{ci}} t \rightarrow t, \quad e \phi / T_{\rm e} \rightarrow \phi,\quad n_{1} / n_{0} \rightarrow n, \end{equation}

where $\rho _s$ is the ion Larmor radius at electron temperature $T_{\rm e}$. It is defined as ${\rho _s = {\sqrt {T_{\rm e}/m_{\rm i}}}/{\omega _{\rm ci}}}$, where $\omega _{\rm 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:

(2.2)\begin{gather} \left(\frac{\partial}{\partial t}- \mu_{\nu} \nabla^{2}\right) \nabla^{2} \phi=[\nabla^{2} \phi, \phi]+c(\phi-n), \end{gather}

(2.3)\begin{gather}\left(\frac{\partial}{\partial t}-\mu_{\rm D} \nabla^{2}\right) n=[n, \phi]-\kappa \frac{\partial \phi}{\partial y}+c(\phi-n), \end{gather}

where $\mu _{\rm D}$ is cross-field diffusion coefficient and $\mu _{\nu }$ is kinematic viscosity. The term $\kappa$, defined as $\kappa \equiv - \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 $\boldsymbol {E} \times \boldsymbol {B}$ velocity, represented by $\boldsymbol {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

(2.4)\begin{equation} c=\frac{T_{\rm e} k_{{\parallel}}^{2}}{\mathrm{e}^{2} n_{0} \eta \omega_{\rm ci}}. \end{equation}

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 $\boldsymbol {E} \times \boldsymbol {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 (Reference Pushkarev, Bos and Nazarenko2013) 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. Reference Bos, Kadoch, Neffaa and Schneider2010). For readers interested in the results of other regimes, we refer to Appendix A.

2.2 Impurity particles model

Earlier research (Futatani et al. Reference Futatani, Benkadda, Nakamura and Kondo2008a,Reference Futatani, Benkadda, Nakamura, Kondo and Hamaguchib, Reference Futatani, Benkadda and del Castillo-Negrete2009) 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. Reference Oujia, Matsuda and Schneider2020), a dimensionless parameter that quantifies particle inertia in fluid flow. The Stokes number is defined as

(2.5)\begin{equation} St = \frac{\tau_{\rm p}}{\tau_{\eta}}, \end{equation}

where $\tau _{\rm p}$ denotes the impurity particle's relaxation time, it measures the time the particle takes to adapt to fluid flow alterations. A smaller $\tau _{\rm 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. Reference Futatani, Benkadda, Nakamura and Kondo2008a,Reference Futatani, Benkadda, Nakamura, Kondo and Hamaguchib, Reference Futatani, Benkadda and del Castillo-Negrete2009), the impurity particles are treated as passive tracers of the flow, indicating $\tau _{\rm p} = 0$. However, here we consider heavy impurity particles with $\tau _{\rm 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_{\rm drag}$ dependent on the relative velocity between the macroscopic velocity of the plasma ion and the velocity of the impurity particle (Kilgore et al. Reference Kilgore, Daugherty, Porteous and Graves1993): $F_{\rm drag} = K_{\mathrm {mt}} n_{\rm i} m_{\rm i} v_{\rm r}$, where $v_{\rm r}$ is the relative velocity between the macroscopic velocity of the plasma ions and the impurity particles, $K_{\mathrm {mt}}$ is the momentum transfer coefficient, $n_{\rm i}$ is the plasma ion density and $m_{\rm 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 _{\rm 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_{\rm p}(u_{\rm p}-v_{\rm p})/\tau _{\rm p}$, which is proportional to the relative velocity between the plasma flow and impurities. The relaxation time can be expressed as $\tau _{\rm p}={m_{\rm p}}/{K_{\mathrm {mt}} n_{\rm i} m_{\rm i}}$. This relation highlights that $\tau _{\rm 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_{\rm 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:

(2.6)\begin{gather} m_{\rm p} \frac{{\rm d} \boldsymbol{v}_{{\rm p},j}}{{\rm d} t} = \frac{m_{\rm p}(\boldsymbol{u}_{{\rm p},j}- \boldsymbol{v}_{{\rm p},j}) }{\tau_{\rm p}} + Ze(\boldsymbol{E}(\boldsymbol{x}_{{\rm p},j}) + \boldsymbol{v}_{{\rm p},j} \times \boldsymbol{B}), \end{gather}

(2.7)\begin{gather}\boldsymbol{E}(\boldsymbol{x}_{{\rm p},j})={-}\boldsymbol{\nabla} \phi(\boldsymbol{x}_{{\rm p},j}). \end{gather}

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

(2.8)\begin{equation} \frac{{\rm d} \boldsymbol{v}_{{\rm p},j}}{{\rm d} t}=\frac{(\boldsymbol{u}_{{\rm p},j}-\boldsymbol{v}_{{\rm p},j})}{\tau_{\rm p}}+ \alpha (-\boldsymbol{\nabla} \phi(\boldsymbol{x}_{{\rm p},j})+\boldsymbol{v}_{{\rm p},j} \times \boldsymbol{b}),\end{equation}

where $\alpha ={Zm_{\rm i}}/{m_{\rm p}}$, $m_{\rm i}$ is the ion mass of the plasma and $\boldsymbol {b}$ represents the unit vector along the direction of the magnetic field. Combining with the equation:

(2.9)\begin{equation} \frac{{\rm d}\kern0.07em\boldsymbol{x}_{{\rm p},j}}{{\rm d} t} = \boldsymbol{v}_{{\rm p},j},\end{equation}

it is thus possible to compute the trajectory of the impurity particles.

2.3 Modified Voronoi tessellation

Impurity particles satisfy the conservation equation:

(2.10)\begin{equation} D_t n_{\rm p}={-}n_{\rm p} \boldsymbol{\nabla} \boldsymbol{{\cdot}} {\boldsymbol{v}_p},\end{equation}

where $D_t$ is the Lagrangian derivative, $n_{\rm p}$ is the impurity particle number density and $\boldsymbol {v}_{\rm p}$ the particle velocity. From this equation, the divergence of particle velocity $\boldsymbol {\nabla } \boldsymbol {{\cdot }} {\boldsymbol {v}_{\rm p}} = -{D_t n}/{n}$ is crucial to understand. A positive divergence ($\boldsymbol {\nabla } \boldsymbol {{\cdot }} {\boldsymbol {v}_{\rm p}} > 0$) leads to density decrease, whereas a negative divergence, i.e. convergence ($\boldsymbol {\nabla } \boldsymbol {{\cdot }} {\boldsymbol {v}_{\rm 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. (Reference Oujia, Matsuda and Schneider2020) and Maurel-Oujia et al. (Reference Maurel-Oujia, Matsuda and Schneider2023). 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.

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.

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_{\rm p}$ in the discrete setting. It is proved that (Oujia et al. Reference Oujia, Matsuda and Schneider2020; Maurel-Oujia et al. Reference Maurel-Oujia, Matsuda and Schneider2023)

(2.11)\begin{equation} {\mathcal{D}} ({\boldsymbol{v}_{\rm p}}) ={-}\frac{1}{n_{\rm p}} D_t n_{\rm p} = \frac{2}{\Delta t} \frac{V^{k+1}-V^k}{V^{k+1}+V^k} + O\left(\frac{1}{N_{\rm p}}, \Delta t\right).\end{equation}

Thus, by measuring the volume change of each particle between $t^k$ and $t^{k+1}$, we can determine the discreet divergence of its velocity ${\mathcal {D}}({\boldsymbol {v}_{\rm 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 \times 64$ with periodic boundary condition. The domain was discretised to a resolution $R$ of $1024 \times 1024$ grid points. The time step $\Delta t$ was $5\times 10^{-4}$. The number of impurity particles $N_{\rm p}$ was $10^6$. The adiabaticity parameter $c$ is $0.7$. The values of specific physical parameters are listed in table 1. The characteristic time scale of turbulence, $\tau _{\eta }$, is defined as $1 / \sqrt {2 Z_{\rm ms}}$, with $Z_{\rm ms}$ denoting one half of the mean-square vorticity. In the simulation, $\tau _{\eta } = 0.35$ in the statistically steady flow. The simulations are based on the work of Kadoch et al. (Reference Kadoch, del Castillo-Negrete, Bos and Schneider2022) 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. Reference Kadoch, del Castillo-Negrete, Bos and Schneider2022) and one million uniformly distributed heavy impurity particles with random velocity are injected.

Table 1.

Simulation parameters for the flow: $A$, domain area; $R$, grid resolution; $\Delta t$, time step; $N_{\rm p}$, number of impurity particles; $\mu _{\rm D}$, diffusion coefficient; $\mu _{\nu }$, kinematic viscosity; $\kappa \equiv - \partial _{x} \ln (n_{0})$, a measure of the plasma density gradient; $c$, adiabaticity parameter.

For the impurity particles, we focus on tungsten $_{74}^{184}\mathrm {W}^{20+}$, a heavy metal that can be made from the material that is exposed to the high heat flux in ASDEX-Upgrade (Krieger, Maier & Neu Reference Krieger, Maier and Neu1999) and EAST (Yao et al. Reference Yao2016), expecting the same situation for ITER (Pitts et al. Reference Pitts2009). For $_{74}^{184}\mathrm {W}^{20+}$, $\alpha = Zm_{\rm i}/m_{\rm p} = 0.22$. As for Stokes number $St = \tau _{\rm p}/\tau _\eta$, $\tau _\eta = 0.35$ in the statistically steady state of flow for $c = 0.7$ (cHW model), but the exact value of $\tau _{\rm p}$ for $_{74}^{184}\mathrm {W}^{20+}$ is unknown. To resolve this, we explore a range of $\tau _{\rm p}$ 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 \ll 1$) to those that are essentially independent of the flow ($St \gg 1$). The simulation parameters are listed in table 2.

Table 2.

Stokes numbers used in the simulation.

The system of HW equations, specifically (2.2) and (2.3) are solved using a classical pseudospectral method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2007). 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. (Reference Kadoch, del Castillo-Negrete, Bos and Schneider2022).