# Impurity Parallel Velocity Gradient instability

J. Bourgeois  Ecole Polytechnique, Institut Polytechnique de Paris, 91128
Palaiseau, France  Interdisciplinary Graduate School of Engineering Sciences
Kyushu University, Fukuoka, Japan  M. Lesur  Université de Lorraine, CNRS,
Institut Jean Lamour, UMR 7198, F-54000 Nancy, France  maxime.lesur@univ-
lorraine.fr G. Cuerva Lazaro  Université de Lorraine, CNRS, Institut Jean
Lamour, UMR 7198, F-54000 Nancy, France  Y. Kosuga  Research Institute for
Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
Research Center for Plasma Turbulence, Kyushu University, Kasuga, Fukuoka
816-8580, Japan

###### Abstract

In magnetized plasmas, a radial gradient of parallel velocity, where parallel
refers to the direction of magnetic field, can destabilise an electrostatic
mode called as Parallel Velocity Gradient (PVG). The theory of PVG has been
mainly developed assuming a single species of ions. Here, the role of
impurities is investigated based on a linear, local analysis, in a
homogeneous, constant magnetic field. To further simplify the analysis, the
plasma is assumed to contain only two ion species - main ions and one impurity
species - while our methodology can be straightforwardly extended to more
species. In the cold-ion limit, retaining polarization drift for both main
ions and impurity ions, and assuming Boltzmanian electrons, the system is
described by 4 fluid equations closed by quasineutrality. The linearized
equations can be reduced to 2 coupled equations: one for the electric
potential, and one for the effective parallel velocity fluctuations, which is
a linear combination of main ion and impurity parallel velocity fluctuations.
This reduced system can be understood as a generalisation of the Hasegawa-Mima
model. With finite radial gradient of impurity parallel flow, the linear
dispersion relation then describes a new instability: the impurity PVG
(i-PVG). Instability condition is described in terms of either the main ion
flow shear, or equivalently, an effective flow shear, which combines main ion
and impurity flow shears. Impurities can have a stabilising or destabilising
role, depending on the parameters, and in particular the direction of main
flow shear against impurity flow shear. Assuming a reasonable value of
perpendicular wavenumber, the maximum growthrate is estimated, depending on
impurity mass, charge, and concentration.

##  1 Introduction

In magnetized plasmas, the Parallel Velocity Gradient (PVG), or Parallel
Velocity Shear (PVS) is a type of Kelvin-Helmholtz fluid instability driven by
a radial gradient of parallel (to magnetic field) plasma flow. It is sometimes
referred to as D’Angelo instability, since D’Angelo investigated a simplified
version (simplified by a radial WKB approximation) for a low-temperature
plasma in a uniform magnetic field [1]. The D’Angelo instability was then
observed in basic plasma experiments [2, 3].

In toroidal magnetic confinement fusion devices, magnetic shear is stabilising
[4]. However, theory indicates that PVG-driven turbulence may be found in the
vicinity of transport barriers [5, 6], in the SOL [7], in the presence of
strong parallel beam injection, and more readily in spherical tokamaks [8, 9].
Even in cases where PVGs are weak, they can have important impact by coupling
with Ion-Temperature-Gradient (ITG) turbulence. Theory predicted subcritical
ITG-PVG turbulence [10, 11, 12, 13], which is consistent with measurements (by
Beam Emission Spectroscopy) on MAST [14]. Strong experimental hints, from
fluctuation level and isotropy of correlation length, indicate a significant
contribution of PVG-driven turbulence in the edge of the CT-6B tokamak plasma
[15]. The universality of PVGs in tokamaks remains an open issue.

Parallel flows are also expected to play major roles in linear magnetized
plasma devices, such as PANTA (Plasma Assembly for Nonlinear Turbulence
Analysis, formerly LMD), where an uphill, near-axis, axial particle flux [16]
has been measured to be consistent with PVG/drift-waves coupling [17], with a
regime transition in quantitative agreement with the theoretical linear
instability threshold [18].

Magnetic confinement fusion plasmas are often contaminated by ion species
other than hydrogen isotopes, which are then called as impurities. These can
include helium, nitrogen, neon, argon, beryllium, carbon, and tungsten. Linear
magnetized plasma experiments can also include impurities. In particular, a
new linear device, called SPEKTRE (Sheath, Plasma Edge Kinetic Turbulence
Radiofrequency Experiment) [19] and currently under construction, is partly
designed to inject various impurities in a controlled manner and investigate
their impact and their dynamics.

Our objective here is to investigate the instability driven by a radial shear
of the parallel flows of multiple species, which we refer to as impurity-PVG
(i-PVG). We consider experimental conditions that are relevant to linear
plasma experiments with cylindrical geometry such as PANTA and SPEKTRE, and
are insightful approximations of tokamak’s edge plasma, where the PVG
instability is of its most significance. Our local, linearized model can be
understood as a modified Hasegawa-Mima model including equilibrium and
perturbed parallel velocities of both main ions and impurities.

In the case of slab geometry with weak magnetic shear, Guo obtained from
gyrokinetic equations a linear dispersion relation for PVG with impurities,
showing the stabilization of the PVG instability by the presence of impurities
and by the increase of electron density gradient [20]. Here we adopt a fluid
approach, and investigate how the i-PVG threshold, frequency and growth-rate
depend on parameters such as impurity charge, mass, concentration and parallel
flow shear. We further assume a cold-ion limit to provide clearer analysis of
the i-PVG. Our fluid approach is supported by this cold ion assumption, plus
by the fact that our study focuses on the linear growth of the instabilities
and ignores the damped states, so that kinetic effects such as Landau damping
can be neglected. Compare to Guo’s previous work, our approach provides a more
flexible analysis in terms of parameter scan for independent ion and impurity
flow shears, and demonstrates that the presence of impurity can be either
stabilizing or destabilizing for PVG. The relationship between linear
properties and plasma parameters turns out to be somewhat complex, with
several non-monotonous dependencies.

##  2 Model

The model is based on a local approximation in a Cartesian basis
(${\boldsymbol{\hat{e_{x}}}}$, ${\boldsymbol{\hat{e_{y}}}}$,
${\boldsymbol{\hat{e_{z}}}}$). The magnetic field is assumed homogeneous and
constant, ${\boldsymbol{B}}=B{\boldsymbol{\hat{e_{z}}}}$. Hereafter, the
parallel and perpendicular subscripts ($\parallel$ and $\perp$) refer to the
direction of the magnetic field ${\boldsymbol{\hat{e_{z}}}}$ and the
perpendicular plane (${\boldsymbol{\hat{e_{x}}}}$,
${\boldsymbol{\hat{e_{y}}}}$). For each species $s$, both equilibrium density
gradient and equilibrium parallel velocity gradient are assumed constant and
in the $x$ direction:

${\boldsymbol{\nabla}}{n_{s,0}}\;=\;\frac{n_{s,0}}{L_{n,s}}{\boldsymbol{\hat{e_{x}}}}$
(1)
${\boldsymbol{\nabla}}{u_{\parallel,s,0}}\;=\;\frac{u_{\parallel,s,0}}{L_{u,s}}{\boldsymbol{\hat{e_{x}}}}$
(2)

where $L_{n,s}$ and $L_{u,s}$ are the density gradient length and parallel
velocity gradient length.

We assume for simplicity that the plasma contains only electrons ($s=e$) and
two ion species (one main ion species $s=i$, and one impurity species $s=z$),
although the model can be straightforwardly generalized to more species. For
each species, mass is noted $m_{s}$, and charge is noted $q_{s}=Z_{s}e$, where
$e$ is the elementary charge ($Z_{e}=-1$ and for simplicity we assume that the
main ions satisfy $Z_{i}=1$).

We adopt the cold ion limit, setting for the wave number an approximate
maximum limit $\sqrt{T_{e}/T_{s}}/\rho_{s}$, with $s=i$ and $z$, and
$\rho_{c,s}$ being the Larmor radius of species $s$. The perpendicular motion
of ions (both main ions and impurity ions) is dominated by $E\times B$ drift,

${\boldsymbol{u_{E}}}\;=\;\frac{{{\boldsymbol{E}}\times{\boldsymbol{B}}}}{B^{2}}$
(3)

and polarization drift,

${\boldsymbol{u_{p,s}}}\;=\;\frac{1}{\omega_{c,s}\,B}\frac{{\mathrm{d}}{\boldsymbol{E_{\perp}}}}{{\mathrm{d}}t}$
(4)

This is valid up to first order in terms of the ratio between wave frequency
and cyclotron frequency, $\omega/\omega_{c,s}$. Here, as verified a
posteriori, $\omega/\omega_{c,i}$ remains as small as a few percents. However,
since $\omega_{c,z}/\omega_{c,i}=Z_{z}m_{i}/m_{z}$, the model requires heavy
impurities to be sufficiently ionized (unless these impurities are in high
enough concentration to bring $\omega$ down to small enough values).

These assumptions yield a system of four fluid equations. Namely, for each of
$s=i$ and $z$, the density equation writes as

$\frac{\partial n_{s}}{\partial
t}\,+\,{{\boldsymbol{{\boldsymbol{\nabla}}{{}_{\perp}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{[}}}n_{s}({\boldsymbol{u_{E}}}+{\boldsymbol{u_{p,s}}})]\,+\,\nabla_{\parallel}(n_{s}u_{\parallel,s})\;=\;0$
(5)

and the parallel flow equation writes as

$\frac{\partial u_{\parallel,s}}{\partial
t}\,+\,{{\boldsymbol{({\boldsymbol{u_{E}}}+{\boldsymbol{u_{p,s}}})}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}{{}_{\perp}}}}}u_{\parallel,s}\,+\,u_{\parallel,s}\nabla_{\parallel}u_{\parallel,s}\;=\;-\frac{q_{s}}{m_{s}}\nabla_{\parallel}\phi$
(6)

Assuming that electron density $n_{e}$ responds to electrostatic fluctuations
$\phi$ according to a Boltzmann distribution, the system of equations is
closed by quasi-neutrality,

$\sum_{s\neq
e}Z_{s}n_{s}\;=\;n_{e}\;=\;n_{0}\,\exp{\left(\frac{e\phi}{T_{e}}\right)}$ (7)

where $T_{e}$ is the temperature.

##  3 Linear analysis

To obtain the linear dispersion relation, densities and parallel velocities
are split into equilibrium and fluctuation parts,
$n_{s}=n_{s,0}+\tilde{n_{s}}$ and
$u_{\parallel,s}=u_{\parallel,s0}+\tilde{u}_{\parallel,s}$. Then, Eqs. (5)-(7)
are linearized,

$\frac{\partial\tilde{n_{s}}}{\partial
t}\,+\,{{\boldsymbol{{\boldsymbol{u_{E}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}n_{s,0}\,+\,n_{s,0}{{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{u_{p,s}}}}}}\,+\,n_{s0}\nabla_{\parallel}\tilde{u}_{\parallel,s}\;=\;0$
(8)
$\partial_{t}\tilde{u}_{\parallel,s}+{{\boldsymbol{{\boldsymbol{u_{E}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}u_{\parallel,s,0}\;=\;-\frac{q_{s}}{m_{s}}\nabla_{\parallel}\phi$
(9) $\sum_{s\neq e}Z_{s}\tilde{n_{s}}\;=\;n_{e,0}\,\frac{e\phi}{T_{e}}$ (10)

Here we used the assumption of homogeneous magnetic field to remove the term
in
${{\boldsymbol{{\boldsymbol{\nabla}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{u_{E}}}}}}$.
In addition, we neglected four terms which remain small for typical plasma
parameters:

  1. 1.

The term
${{\boldsymbol{{\boldsymbol{u_{p,s}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}n_{s,0}$
is of order $\omega/\omega_{c,s}$ compared to
${{\boldsymbol{{\boldsymbol{u_{E}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}n_{s,0}$.
Note that, compared to
$n_{s,0}{{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{u_{p,s}}}}}}$,
it is of order $k_{x}/(k_{\perp}^{2}L_{n})$, which, as verified a posteriori,
remains much smaller than unity as long as $\rho_{c,i}/L_{n}\gg 1$.

  2. 2.

Similarly, the term
${{\boldsymbol{{\boldsymbol{u_{p,s}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}u_{\parallel,s,0}$
is of order $\omega/\omega_{c,s}$ compared to
${{\boldsymbol{{\boldsymbol{u_{E}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}u_{\parallel,s,0}$.

  3. 3.

The term $u_{\parallel,s0}\nabla_{\parallel}\tilde{n_{s}}$, compared to
$\partial\tilde{n_{s}}/\partial t$, is of order
$k_{\parallel}u_{\parallel,s0}/\omega$, which must remain small to avoid
strong Landau damping. Note that, compared to
${{\boldsymbol{{\boldsymbol{u_{E}}}}}\,{\mbox{\boldmath$\cdot$}}\,{\boldsymbol{{\boldsymbol{\nabla}}_{\perp}}}}n_{s,0}$,
it is of order $(T_{i}/T_{e})M_{s}^{2}$, where
$M_{s}=u_{\parallel,s0}/v_{T,i}$ is the Mach number for species $s$. To obtain
the latter ordering, we substituted Eq. (20), and assumed that for each
species, density and parallel gradient lengths are comparable. Extension of
this theory for plasmas where the condition $M_{s}^{2}\ll T_{e}/T_{i}$ is not
satisfied, is left for future work, as accounting for the term
$u_{\parallel,s0}\nabla_{\parallel}\tilde{n_{s}}$ expands a parameter space
that is already very large.

  4. 4.

Similarly, the term
$u_{\parallel,s0}\nabla_{\parallel}\tilde{u}_{\parallel,s}$, compared to
$\partial\tilde{u}_{\parallel,s}/\partial t$, is of order
$k_{\parallel}u_{\parallel,s0}/\omega$.

A linear combination reduces the system Eqs. (8)-(10) to two coupled
equations, on the effective variables $\Phi=e\phi/T_{e}=\sum
Z_{s}\tilde{n_{s}}/n_{e,0}$ and $\tilde{\bar{u}}_{\parallel}=\sum
Z_{s}C_{s}\tilde{u}_{\parallel,s}$, where $C_{s}=n_{s,0}/n_{e,0}$ is the
equilibrium concentration of species $s$. The two coupled equations write

$\partial_{t}\left(1-\bar{\rho}^{2}\nabla^{2}_{\perp}\right)\Phi\,+\,v_{\ast,e}\,\partial_{y}\Phi\,+\,\nabla_{\parallel}\tilde{\bar{u}}_{\parallel}\;=\;0$
(11) $\partial_{t}\tilde{\bar{u}}_{\parallel}\,-\,\overline{\nabla
u_{\parallel}}\,\rho_{c,i}v_{T,i}\,\partial_{y}\Phi\;=\;-c_{s}^{2}\nabla_{\parallel}\Phi$
(12)

where

${\overline{\rho}}^{2}\;=\;\sum_{s\neq e}{C_{s}Z_{s}^{2}\rho_{c,s}^{2}}$ (13)

is a linear combination of ion Larmor radii, $v_{\ast,e}=T_{e}/(eBL_{n,e})$ is
the electron diamagnetic drift velocity, $c_{s}$ is the ion-acoustic velocity
including the contribution of impurities,

$c_{s}^{2}\;=\;\sum_{s\neq e}{C_{s}\,Z_{s}^{2}\,\frac{T_{e}}{m_{s}}}$ (14)

(recalling we adopt the cold ion limit), and, most importantly,

$\overline{\nabla u_{\parallel}}\;=\;\sum_{s\neq e}{C_{s}\,Z_{s}\,\nabla
u_{\parallel,s,0}}$ (15)

is a linear combination of main ion and impurity parallel velocity gradients,
which we refer to as the _effective flow shear_. It is the source of free-
energy of the i-PVG instability. Note that throughout this paper, we note the
radial gradient ($\partial/\partial x$) as $\nabla$ (without a subscript). Let
us define the ratio of parallel velocity gradients, $V=\nabla
u_{\parallel,z,0}/\nabla u_{\parallel,i,0}$. Then the effective flow shear can
also be written as $\overline{\nabla
u_{\parallel}}=(C_{i}+V\,C_{Z}\,Z_{Z})\,\nabla u_{\parallel,i,0}$.

The above model can in fact be seen as a modification of the Hasegawa-Mima
model, including the effective flow shear.

Fourier analysis yields the dispersion relation,

$\omega^{2}\,-\,\omega_{\mathrm{iDW}}\,\omega\,-\,\frac{\left(k_{\parallel}c_{s}\right)^{2}}{1\,+\,\left(k_{\bot}\overline{\rho}\right)^{2}}\left(1-\frac{k_{y}}{k_{\parallel}}\frac{\overline{\nabla
u_{\parallel}}}{\overline{\omega_{c}}}\right)\;=\;0$ (16)

where

$\omega_{\mathrm{iDW}}\;=\;\frac{k_{y}\,v_{\ast,e}}{1\,+\,\left(k_{\bot}\overline{\rho}\right)^{2}}$
(17)

is the drift-wave frequency modified by impurities, and

$\overline{\omega_{c}}\;=\;\sum_{s\neq e}{C_{s}\,Z_{s}\,\omega_{c,s}}$ (18)

is a linear combination of main ion and impurity cyclotron frequencies.

The quadratic Eq. (8) includes an unstable solution if and only if the
function

$\Delta(k_{\parallel})\;=\;\left(\frac{\omega_{iDW}}{2}\right)^{2}+\frac{\left(k_{\parallel}c_{s}\right)^{2}}{1\,+\,\left(k_{\bot}\overline{\rho}\right)^{2}}\left(1-\frac{k_{y}}{k_{\parallel}}\frac{\overline{\nabla
u_{\parallel}}}{\overline{\omega_{c}}}\right)$ (19)

is negative.

When the parallel wave number is $k_{\parallel}=k_{\parallel,\mathrm{max}}$,
with

$k_{\parallel,\mathrm{max}}\;=\;\frac{\overline{\nabla
u_{\parallel}}}{2\,\overline{\omega_{c}}}\,k_{y}$ (20)

the function $\Delta$ reaches its minimum value

$\Delta_{min}\;=\;\left(\frac{\omega_{iDW}}{2}\right)^{2}\,\left(1\,-\,L_{n,e}^{2}\,(1+k_{\perp}^{2}\bar{\rho}^{2})\,\frac{\overline{\nabla
u_{\parallel}}^{2}}{\overline{c_{s}}^{2}}\right)$ (21)

Therefore, the linear instability condition can be written in terms of the
effective flow shear, $\left|\overline{\nabla
u_{\parallel}}\right|>\overline{\nabla u_{\parallel}}_{\mathrm{cr}}$ with

$\overline{\nabla
u_{\parallel}}_{\mathrm{cr}}^{2}\;=\;\frac{{(c_{s}/L_{n,e})}^{2}}{1+\left(k_{\bot}\overline{\rho}\right)^{2}}$
(22)

or, alternatively, in terms of the main ion flow shear, $\left|\nabla
u_{\parallel,i,0}\right|>\nabla u_{i,\mathrm{cr}}$ with

$\nabla u_{i,\mathrm{cr}}\;=\;\frac{1}{C_{i}+V\,C_{Z}\,Z_{Z}}\overline{\nabla
u_{\parallel}}_{\mathrm{cr}}$ (23)

When an i-PVG wave is unstable, its frequency is $\omega_{\mathrm{iDW}}/2$ and
its growthrate is $\gamma=\sqrt{-\Delta}$. These expressions are given in Eqs.
(17) and (19). Note that the factor 2 in the frequency is mainly due to the
finite value of $k_{\parallel}$ – the i-PVG frequency does recover the
impurity-modified drift-wave frequency in the limit of $k_{\parallel}=0$ and
$\overline{\nabla u_{\parallel}}=0$.

For a given set of plasma parameters and a given value of $k_{\perp}$, the
maximum growthrate is reached for $k_{x}=0$ and
$k_{\parallel}=k_{\parallel,\mathrm{max}}(k_{y})$, in which case

$\gamma_{max}\;=\;\frac{k_{\perp}\,v_{\ast,e}/2}{1\,+\,k_{\perp}^{2}\bar{\rho}^{2}}\,\left[L_{n,e}^{2}\,(1+k_{\perp}^{2}\bar{\rho}^{2})\,\frac{\overline{\nabla
u_{\parallel}}^{2}}{c_{s}^{2}}\,-\,1\right]^{1/2}$ (24)

Hereafter we focus on that case $k_{x}=0$, which is also more consistent with
our local model.

##  4 Application to typical plasmas

For concision, let us note $M=m_{z}/m_{i}$ the mass ratio, and $V=\nabla
u_{\parallel,z,0}/\nabla u_{\parallel,i,0}$ the ratio of parallel velocity
gradients. Table 1 shows the parameters of a reference case, designed to be
representative of a fusion plasma with highly ionized tungsten. The values
$\nabla u_{\parallel,i,0}$ and $V=2$ are consistent with measurements [21,
22]. Hereafter, we present analyses of cases which are variations of this
reference case. The varying parameters are summarized in Table 2.

| $L_{n}/\rho_{c,i}$ | $\nabla u_{\parallel,i,0}/\omega_{c,i}$ | $V$ | $M$ | $Z$ | $C_{z}$ | $\rho_{c,i}k_{x}$ | $\rho_{c,i}k_{y}$ | $\rho_{c,i}k_{\parallel}$
---|---|---|---|---|---|---|---|---|---
Reference case | 10 | 0.1 | 2 | 184 | 40 | $10^{-3}$ | 0 | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Table 1:  Input parameters for the reference case. | $\nabla u_{\parallel,i,0}/\omega_{c,i}$ | $V$ | $M$ | $Z$ | $C_{z}$ | $\rho_{c,i}k_{y}$ | $\rho_{c,i}k_{\parallel}$
---|---|---|---|---|---|---|---
Fig. 1 | 0.1 | 2 | 12 and 184 | 6 and 40 | 0 and $10^{-3}$ | 0.7 | -0.02 – 0.08
Fig. 2 | 0.1 | 2 | 184 | 40 | $10^{-3}$ | -1 – 1 | -0.1 – 0.1
Fig. 3 | 0.1 | 2 | 12 and 184 | 6 and 40 | 0 and $10^{-3}$ | 0 – 2 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 4 | - | 2 | 0 – 200 | 0 – 200 | $10^{-3}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 5 | 0.1 | -1 – 2 | 5 – 200 | 5 |  $10^{-4}$ – $10^{-2}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 6 | 0.1 | -1 – 2 | 184 | 0 – 50 |  $10^{-4}$ – $10^{-2}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 7 | 0.1 | -1 – 2 | 12 – 184 | 3 – 40 |  $10^{-5}$ – $10^{-2}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 8 | 0.08 | 2 | 12 and 184 | 6 and 40 | 0 and $10^{-3}$ | 0 – 2 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 9 | 0.08 | 2 | 0 – 200 | 0 – 200 | $10^{-3}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 10 | 0.08 | -1 – 2 | 5 – 200 | 5 |  $10^{-4}$ – $10^{-2}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 11 | 0.08 | -1 – 2 | 184 | 0 – 50 |  $10^{-4}$ – $10^{-2}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Fig. 12 | 0.08 | -1 – 2 | 12 – 184 | 3 – 40 |  $10^{-5}$ – $10^{-2}$ | 0.7 | $\rho_{c,i}k_{\parallel,\mathrm{max}}$
Table 2:  Input parameters for each figure. This table includes only
parameters which vary compared to the reference case. In other words, for all
figures, $L_{n}/\rho_{c,i}=10$ and $k_{x}=0$.

###  4.1 Dispersion relation

Figure 1:  Dispersion relation of the i-PVG: frequency and growthrate against
parallel wavenumber, for the reference case (”High-Z Tgst”), as well as a case
with C6+ impurity (”High-Z Carb”). The dispersion relation of the pure PVG is
included for comparison. The curves for the pure PVG are almost hidden behind
the curves for the C6+ impurity. The cross and plus symbols correspond to
linear frequency and growthrate measured in an initial value simulation of the
fluid model Eqs. (5-7).

Figure 1 shows the dispersion relation for the reference case, as well as for
C6+ ($M=12$ and $Z=6$) impurity. The dispersion relation of the pure PVG is
included for comparison. We observe that the dispersion relation is
qualitatively similar for the pure and impure cases, but quantitatively
different. As demonstrated hereafter, differences are more important closer to
the PVG instability threshold, or for higher impurity concentrations.

Let us now verify both the accuracy of the assumptions (i)–(iv) described in
Sec. 3, and the correctness of our linear analysis. We developed an initial
value numerical simulation code to solve the fluid model Eqs. (5-7), which
does not rely on the assumptions (i)–(iv). Time integration is performed by an
explicit RK4 scheme with time-step width $\Delta t=10^{-2}\omega_{c,i}^{-1}$.
All quantities ($\phi$, $n_{s}$, $u_{\parallel,s}$) are described in Fourier
space ($k_{x}$, $k_{y}$) on a 256x256 grid for a fixed $k_{\parallel}$.
Nonlinear terms such as $n_{s}\,u_{\parallel,s}$ are computed in real space
before being transformed back to Fourier space. Fig. 1 includes the frequency
and growthrate extracted from time-traces of complex Fourier component of
$\phi$, during the linear phase. There is good quantitative agreement for the
values of $k_{\parallel}$ such that the growthrate is close to its maximum –
the relative error for the growthrate remains below $0.2\%$ for the range
$\gamma/\omega_{c,i}>0.014$. The relative error can be as high as a few
percents in the range $0.004<\gamma/\omega_{c,i}<0.007$ or a few tens of
percent in the range $\gamma/\omega_{c,i}<0.004$, but these weak modes are
negligible in the linear phase of the instability. The relative error for the
frequency remains below $0.5\%$ for the whole unstable range.

Figure 2:  Growthrate against parallel and azimuthal wavenumbers, for the
reference case. The dashed black line represents the condition
$k_{\parallel}=k_{\parallel,\mathrm{max}}$. The white curve is the instability
threshold.

Figure 2 shows the linear growthrate in the space of wavevector components
($k_{\parallel}$, $k_{y}$). The asymmetry in this space, instability occurring
only for $k_{y}\,k_{\parallel}>0$ (for $\overline{\nabla u_{\parallel}}>0$),
is an important property of the i-PVG. If the sign of $\overline{\nabla
u_{\parallel}}$ was flipped, the unstable region would be in the region
$k_{y}\,k_{\parallel}<0$. As expected, growthrate is maximum for
$k_{\parallel}=k_{\parallel,\mathrm{max}}$. Hereafter, we always set
$k_{\parallel}=k_{\parallel,\mathrm{max}}$.

Figure 3:  Dispersion relation of the i-PVG: frequency and growthrate against
the azimuthal wavenumber $k_{y}$, for the reference case (”High-Z Tgst”), as
well as a case with C6+ impurity (”High-Z Carb”). The dispersion relation of
the pure PVG is included for comparison. The curves for the pure PVG are
almost hidden behind the curves for the C6+ impurity.

Figure 3 displays the maximum growthrate (obtained for
$k_{\parallel}=k_{\parallel,\mathrm{max}}$) against the azimuthal $k_{y}$
wavenumber, as well as the corresponding frequency. Our model does not include
any mechanism for small scale dissipation. As a consequence,
$\gamma_{\mathrm{max}}$ spuriously keeps on increasing for increasing $k_{y}$.
However, our model looses its validity as $\rho_{c,i}k$ approches unity. With
the knowledge of the typical shape of $\gamma_{\mathrm{max}}(k)$ obtained from
more complete models for many instabilities, let us, as a first, qualitative
approach, approximate the physical maximum growthrate by choosing a cut-off at
$\rho_{c,i}k=0.7$. Then, for our reference case, the maximum growthrate is
$\gamma_{\mathrm{max}}\sim 0.02\,\omega_{ci}$.

In summary, the results in this section confirm that the linear properties of
the i-PVG instability are qualitatively the same as that of the pure PVG
instability. For the chosen set of parameters, even the quantitative
differences are small. However, that is for a case where the pure PVG
instability is already quite strong, since $\gamma\approx\omega$. In this
sense, the above results concern a range far from instability threshold. By
contrast, in the next section we focus on the instability threshold itself.

###  4.2 Critical ionic PVG

As explained in Sec. 3, for a given set of parameters ($L_{n}$, $\nabla
u_{\parallel,i,0}$, $V$, $M$, $Z$, $C_{z}$), and bounded value of $k_{y}$, the
i-PVG instability threshold can be expressed either in terms of the effective
flow shear, which is a combination of main ion and impurity flow shears, or,
alternatively, in terms of the main ion flow shear alone. The corresponding
thresholds are given by Eqs. (22) and (23), respectively. The preferred point-
of-view may depend on the experimental conditions which determine preferred
control parameters. Here we focus on the second point-of-view, that of
critical main ion flow shear for a given ratio $V$.

Figure 4:  Critical main ion flow shear against mass and charge ratio.

Figure 4 displays the critical main ion flow shear against mass and charge
ratio, for $V=2$ and $C_{z}=10^{-3}$. It is clear that i-PVG-unstable
conditions are more readily reached when the impurity is massive and highly
ionized. In realistic machines however, radiation from high-mass and high-Z
impurities is a severe issue, so that the applicability of such ranges of
parameters should be carefully studied. On the contrary, at this low
concentration, the threshold for low-mass, low-charge impurities is not so
different than the threshold for high-mass, low-charge impurities, and is
actually close to the threshold of the pure PVG.

Figure 5:  Critical main ion flow shear against mass ratio, for $Z=5$. Left:
fixed $V=2$ and various concentrations. Right: fixed $C_{z}=10^{-3}$ and
various $V$.

Figure 5 illustrates the role of the mass ratio, at fixed $Z=5$, for various
impurity concentrations $C_{z}$ and various flow shear ratios $V$. The
critical main ion flow shear is always decreasing with increasing mass ratio,
but the relationship is nonlinear. Here the impact of $V$ is small, but that
is because $C_{z}$ and $Z$ are both low, such that the i-PVG is actually very
close to the pure PVG.

Figure 6:  Critical main ion flow shear against charge number $Z$, for
$M=184$. Left: fixed $V=2$ and various concentrations. Right: fixed
$C_{z}=10^{-3}$ and various $V$.

Figure 6 illustrates role of the charge number $Z$ of the impurity, at fixed
mass ratio $M=184$ (tungsten), for various impurity concentrations $C_{z}$ and
various flow shear ratios $V$. Higher $Z$ can be either stabilising or
destabilising depending on the parameters, such as $V$. Here, at concentration
$C_{z}=10^{-3}$, it is stabilising for $V=2$, but destabilising for $V=0.5$
and $V=-1$. This stabilising role of impurities is consistent with simple
intuition in the case $V=-1$, since impurity flow counteracts main ion flow in
this case. In the case $V=0.5$, the stabilising effect is more subtle: in this
case the variation in Eq. (23) is dominated by the denominator
$(C_{i}+VC_{Z}Z)$, which can be re-written as $[1+ZC_{Z}(V-1)]$. However,
which term dominates depends on the value of $V$. For example, for $V=1$, the
variation as $Z$ increases is dominated by the decrease of $c_{s}$, so
increasing $Z$ is stabilising. For an intermediate case, $V=0.8$, the
dependency is non-monotonous: increasing $Z$ is stabilising for $Z<28$ and
destabilising for $Z\geq 28$.

Figure 7:  Critical main ion flow shear against impurity concentration. Left:
fixed $V=2$ and various species. Right: various $V$ and the impurity is W40+.

A more intuitive picture can be obtained by varying impurity concentration
while the impurity species and $Z$ are fixed. Figure 7 illustrates the role of
impurity concentration, for various impurity species (C3+, C6+, W5+, W40+) and
various values of $V$. From this point of view, the role of impurity flow is
consistent with simple intuition: if $V>0$, impurity flow adds up to main ion
flow, which is destabilising ; conversely, if $V<0$, impurity flow counters
main ion flow, which is stabilising. Increasing concentration accentuates the
impact – whether stabilising or destabilising – of the impurity on the
instability threshold. This effect may offer new means of control of the PVG
instability by impurity injection.

In summary, the critical main ion flow shear decreases nonlinearly with
increasing mass ratio, and may increase or decrease in a non-trivial manner
depending on all the other parameters.

For our reference case, the instability threshold is $\nabla
u_{i,\mathrm{cr}}=0.075925\,\omega_{c,i}$. In the next subsection, in order to
study the growthrate near the instability threshold, we arbitrarily fix
$\nabla u_{\parallel,i,0}=0.08\,\omega_{c,i}$.

###  4.3 Growthrate

Figure 8:  Dispersion relation of the i-PVG: frequency and growthrate against
azimuthal wavenumber, for the reference case, as well as a case with C6+
impurity (”High-Z Carb”) – same as Fig. 1 but for lower $\nabla
u_{\parallel,i,0}=0.08\,\omega_{c,i}$. The dispersion relation of the pure PVG
is included for comparison. The curves for the pure PVG are almost hidden
behind the curves for the C6+ impurity.

Figure 8 displays the maximum growthrate (obtained for
$k_{\parallel}=k_{\parallel,\mathrm{max}}$) against the azimuthal $k_{y}$
wavenumber, as well as the corresponding frequency, for $\nabla
u_{\parallel,i,0}=0.08\,\omega_{c,i}$. Here, for
$k_{y}<0.75\,\rho_{c,i}^{-1}$, the main ion parallel velocity shear is below
the threshold for pure PVG. As a consequence, the behavior is qualitatively
different compared to Fig. 1 (which is for $\nabla
u_{\parallel,i,0}=0.1\,\omega_{c,i}$, above the pure PVG threshold even for
$k_{y}\rightarrow 0$). There is now an additional bifurcation.

Again, let us approximate the physical maximum growthrate by fixing
$\rho_{c,i}k=0.7$. Then, for our plasma with $W^{40+}$, the maximum growthrate
is $\gamma_{\mathrm{max}}\sim 0.008\,\omega_{c,i}$. Consistently, the maximum
growth-rate in the following figures (9, 10, 11 and 12), where impurity
parameters are varied, is obtained from our model for the wave vector
${\boldsymbol{k}}_{\mathrm{max}}$ with components $k_{x}=0$,
$k_{y}=0.7/\rho_{c,i}$, $k_{\parallel}=k_{\parallel,\mathrm{max}}$.

Figure 9:  Maximum growthrate against mass and charge ratios, for $\nabla
u_{\parallel,i,0}=0.08\,\omega_{c,i}$.

Figure 9 shows our estimated maximum growthrate against mass and charge ratio,
for $V=2$ and $C_{z}=10^{-3}$. Recalling that the real frequency of the mode
is typically $\omega\approx 0.02\,\omega_{c,i}$, we confirm that our choice of
cut-off ($\rho_{c,i}k=0.7$) ensures that $\gamma_{\mathrm{max}}/\omega$
remains significantly below unity, except for unrealistically highly ionized
heavy impurities. Note that here, the main ion parallel velocity shear is
below the threshold for pure PVG (which is $0.08192\,\omega_{c,i}$ for
$k_{y}=0.7\,\rho_{c,i}^{-1}$).

Figure 10:  Maximum growthrate against mass ratio, for $\nabla
u_{\parallel,i,0}=0.08\,\omega_{c,i}$ and $Z=5$. Left: fixed $V=2$ and various
$C_{z}$. Right: fixed $C_{z}=10^{-3}$ and various $V$.

Figure 11:  Maximum growthrate against impurity charge, for $\nabla
u_{\parallel,i,0}=0.08\,\omega_{c,i}$ and $M=184$. Left: fixed $V=2$ and
various concentrations. Right: fixed $C_{z}=10^{-3}$ and various $V$.

Figure 12:  Maximum growthrate against impurity concentration. Left: fixed
$V=2$ and various species. Right: various $V$ and the impurity is W40+.

Figure 10 illustrates the role of the mass ratio on the maximum growthrate,
for fixed $Z=5$ and various $C_{z}$ and $V$. From this perspective, the
instability threshold can be viewed as a threshold in impurity mass, every
other parameters being fixed. Similarly to the effect of mass on the
instability threshold, increasing mass always increases the maximum
growthrate.

Figure 11 illustrates the role of the impurity charge on the maximum
growthrate, for fixed $M=184$ and various $C_{z}$ and $V$. As a rule of thumb,
the growthrate is less sensitive to $Z$ than to $M$. Increasing $Z$ may be
destabilising or stabilising depending on the parameters, and in particular
depending on the flow shear ratio $V$. For $V=0.5$ or $V=-1$, the i-PVG
instability threshold can be viewed as a threshold in impurity charge,
instability occuring only _below_ a critical number.

Figure 12 illustrates the role of impurity concentration on the maximum
growthrate, for various species and $V$. In this point-of-view, the
instability threshold is a threshold in concentration. With our arbitrary
choice of $\nabla u_{\parallel,i,0}=0.08\,\omega_{c,i}=1.05\nabla
u_{i,\mathrm{cr}}$, the concentration threshold is of the order of $10^{-4}$ –
$10^{-3}$ for tunsgten, and of the order of $10^{-3}$ – $10^{-2}$ for carbon.
However, with $\nabla u_{\parallel,i,0}$ closer to $\nabla u_{i,\mathrm{cr}}$,
the concentration threshold can even be much lower. Here, the growthrate
increases with increasing concentration. However, there are also cases where
concentration has a non-monotonous impact [23].

Let us emphasize that figures 11 and 12 are useful to guide experiments which
could test this theory. For example, impurities may be purposely injected in a
linear device, to either excite the i-PVG, or to suppress a pure PVG
(depending on the direction of flow shears). However, the analysis should be
reproduced after measuring the values of parallel flows, since the thresholds
in $Z$ and $C_{Z}$ are sensitive to the value of the main ion flow shear – and
more specifically its closeness to marginal stability.

##  5 Conclusions and discussions

We investigated the role of a radial gradient of the parallel flow of
impurities in a magnetized plasma. The linear dispersion relation indicates
that the impurity-PVG is destabilised when $L_{n,e}^{2}\overline{\nabla
u_{\parallel}}^{2}>c_{s}^{2}/[1+\left(k_{\bot}\overline{\rho}\right)^{2}]$.
Similarly to the pure-PVG, the i-PVG is asymmetric in the space of wavevector
components ($k_{\parallel}$, $k_{y}$): it is unstable only for
$k_{y}\,k_{\parallel}\overline{\nabla u_{\parallel}}>0$.

If the plasma is unstable to the pure-PVG, then impurities can be stabilising
if their flow shear is opposite to the flow shear of main ions. Note that this
does not require flows themselves to be in opposite directions. The role of
impurities is more prominent for higher mass and/or higher charge number. In
general, the instability threshold and the growthrate are more sensitive to
charge number than to mass.

The analytic expression of the maximum growthrate is given in Eq. (24). To
investigate how the maximum growthrate depends on plasma and impurity
parameters (Figs. 10-12), we assumed a perpendicular wavenumber comparable to
the inverse Larmor radius (of main ions).

In the typical case of a hydrogen plasma with a $10^{-3}$ concentration of
$W^{40+}$, tungsten flow shear being twice the hydrogen flow shear, and
assuming that the linear growthrate is maximum for $k_{\perp}\approx
0.7\,\rho_{c,i}^{-1}$, the i-PVG instability is unstable for hydrogen flow
shear above a threshold $\nabla u_{i,\mathrm{cr}}\approx 0.076\,\omega_{c,i}$,
and most unstable for parallel wavenumbers $k_{\parallel}\approx
0.04\,\rho_{c,i}^{-1}$.

Although the present work should in principle capture the main linear
properties of the i-PVG, the reader should keep in mind that it is subject to
several caveats: it is a local analysis, assuming cold ions and Boltzmann
electrons. Our model does not include collisional effects, viscosity,
resistivity. We limited our analysis to two species, but it can be
straightforwardly generalised to more species.

The present work focused on the linear properties of the i-PVG. Although, the
(quasi-)steady-state spectrum of fluctuations may roughly reflect the location
in ${\boldsymbol{k}}$-space of most unstable modes, nonlinear simulations and
analysis are necessary to improve predictions and understanding of steady-
state fluctuations. In particular, the present work cannot address the role of
zonal flows, which can be generated by PVGs [24]. Impurities may play an
interesting role in both the generation and the damping of zonal flows.

Impurities are transported in the radial direction, in general by both
neoclassical and turbulent processes, although one process may dominate
depending on the species and plasma parameters. The efficiency of fusion
depends on the radial distribution of impurities, which results from
transport. For example, slight tungsten contamination of the core yields
prohibitive energy losses. Therefore, it is essential to improve our
understanding of impurity transport. However, the model adopted here assumes a
Boltzmann response of electrons to electric fluctuations, which precludes any
radial transport of electrons, and forces impurity transport to exactly
balance main ion transport. This assumption is reasonable for the presented
linear analysis, but is too strong for providing any conclusion about impurity
transport. Relaxing this assumption is left for future work.

## Acknowledgments

This work was supported by a scholarship by Kyushu University, the joint
research project in RIAM, the Grants-in-Aide for Scientific Research of JSPS
of Japan (JP21H01066), a grant from the Chaire Energies Durables of Ecole
Polytechnique, and funding from the Agence Nationale de la Recherche for the
project GRANUL (ANR-19-CE30-0005).

## References

  * [1] N. D’Angelo, Phys. Fluids 8, 1748 (1965).
  * [2] N. D’angelo and S. V. Goeler, Phys. Fluids 9, 309 (1966).
  * [3] T. Kaneko, H. Tsunoyama, and R. Hatakeyama, Phys. Rev. Lett. 90, 125001 (2003).
  * [4] P. J. Catto, M. N. Rosenbluth, and C. Liu, Phys. Fluids 16, 1719 (1973).
  * [5] D. McCarthy, Phys. Plasmas 9, 2451 (2002).
  * [6] X. Garbet et al., Phys. Plasmas 9, 3893 (2002).
  * [7] J. Drake et al., Nucl. Fusion 32, 1657 (1992).
  * [8] I. Chapman, S. Brown, R. Kemp, and N. Walkden, Nucl. Fusion 52, 042005 (2012).
  * [9] W. Wang et al., Nucl. Fusion 55, 122001 (2015).
  * [10] S. L. Newton, S. C. Cowley, and N. F. Loureiro, Plasma Phys. Control. Fusion 52, 125001 (2010).
  * [11] M. Barnes et al., Phys. Rev. Lett. 106, 175004 (2011).
  * [12] A. Schekochihin, E. Highcock, and S. Cowley, Plasma Phys. Control. Fusion 54, 055011 (2012).
  * [13] E. Highcock et al., Phys. Rev. Lett. 109, 265001 (2012).
  * [14] A. Field et al., Rev. Sci. Instr. 83, 013508 (2012).
  * [15] G. Wang et al., Plasma Phys. Control. Fusion 40, 429 (1998).
  * [16] T. Kobayashi et al., Phys. Plasmas 23, 102311 (2016).
  * [17] S. Inagaki et al., Sci. Rep. 6, 22189 (2016).
  * [18] Y. Kosuga, S.-I. Itoh, and K. Itoh, Plasma Fusion Research 10, 3401024 (2015).
  * [19] F. Brochard et al., in 49th EPS Conference on Plasma Physics (IAEA Publications, Vienna, 2023), p. THS/P5–03.
  * [20] W. Guo, L. Wang, and G. Zhuang, Nucl. Fusion 59, 076012 (2019).
  * [21] D. Testa et al., Phys. Plasmas 9, 243 (2002).
  * [22] B. Grierson et al., Phys. Plasmas 26, 042304 (2019).
  * [23] M. Lesur, J. Bourgeois, and Y. Kosuga, in Fusion Energy 2023 (Proc. 29th Int. Conf. London) (Vienna: IAEA, International Centre, 2023).
  * [24] Y. Kosuga, S.-I. Itoh, and K. Itoh, Phys. Plasmas 24, 032304 (2017).