Isotope effect on filament dynamics in fusion edge plasmas
Ole Hauke Heinz Meyer Alexander Kendl

]Isotope effect on filament dynamics in fusion edge plasmas

^1Institut für Ionenphysik und Angewandte Physik, Universität
  Innsbruck, 6020 Innsbruck, Austria 
^2Department of Physics and Technology, UiT - The Arctic University of Norway,
9037 Tromsø, Norway
ole.meyer@uibk.ac.at

The influence of the ion mass on filament propagation in the scrape-off layer
of toroidal magnetised plasmas is analysed for various fusion relevant majority
species, like hydrogen isotopes and helium, on the basis of a computational
isothermal gyrofluid model for the plasma edge. 
Heavy hydrogen isotope plasmas show slower outward filament propagation and
thus improved confinement properties compared to light isotope plasmas,
regardless of collisionality regimes. Similarly, filaments in fully ionised
helium move more slowly than in deutrium. Different mass effects on the
filament inertia through polarisation, finite Larmor radius, and parallel
dynamics are identified.


Keywords: isotope effect, plasma filament, blob, particle transport
	

[
    O H H Meyer^1,2 and A Kendl^1
    
=================================







§ INTRODUCTION


In various tokamak experiments the confinement properties have been shown to
scale favourably with increasing mass of the main (fusion relevant) ion plasma
species, specifically hydrogen isotopes and helium <cit.>.  
The radial cross-field transport of coherent filamentary structures (commonly
denoted “blobs”) in the scrape-off layer (SOL) of tokamaks accounts for a
significant part of particle and heat losses <cit.> to
the plasma facing components. Experimentally the ion mass effect on SOL filament
dynamics has been studied in a simple magnetised torus <cit.>.

Filamentary transport in tokamaks in general is an active subject of studies in
experiments, analytical theory, and by computations in two and three dimensions.

The basic properties of filamentary transport are reviewed in <cit.>. 
Blob propagation results from magnetic drifts that polarise density
perturbations, thus yielding a dipolar electric potential ϕ whose 
resulting B×∇ϕ drift in the magnetic field B 
drives the filaments down the magnetic field gradient and towards the wall. 
The basic physics is illustrated by accounting for the current  
paths involved upon charging of the blob by the diamagnetic current: 
the closure is via perpendicular polarisation currents in the drift plane and
through parallel divergence of the parallel current <cit.>. 

Two-dimensional (2-d) closure schemes are discussed in Ref. <cit.>. 
Depending on parallel resistivity, the dominant closure path features distinct
dynamics: if closure is mainly through the polarisation current, the 2-d
cross-field properties are dominant, leading to a mushroom-cape shaped radial
propagation. For reduced resistivity the closure parallel to the magnetic
field direction in 3-d is dominant, and Boltzmann spinning leads to a more
coherently propagating structure at significantly reduced radial velocity. 

The isotope mass may have influence on the E× B shearing rate in the
edge region <cit.>, and flow shear in the edge region has been suggested
to be a main agent which controls blob formation <cit.>. 
In addition, a finite ion temperature introduces poloidally asymmetric
propagation of blobs <cit.>. The underlying finite Larmor
radius (FLR) effects have been found to contribute to favourable 
isotopic transport scaling of tokamak edge turbulence <cit.>. 

In this work we study the isotopic mass effect on blob filament propagation by
employing an isothermal gyrofluid model so that relevant FLR contributions to
the blob evolution are effectively included, in addition to the mass
dependencies in polarisation and in parallel ion velocities. 





§ GYROFLUID MODEL AND COMPUTATION
 


The present simulations on the isotopic dependence of 3-d filament and 2-d blob
propagation in the edge and SOL of tokamaks are based on the gyrofluid
electromagnetic model introduced by Scott <cit.>.
In the local delta-f isothermal limit the model consists of
evolution equations for the gyrocenter densities n_s and parallel velocities
u_s of electrons and ions, where the index s denotes the species with
s ∈ (e, i):


    _s n_s/ t =     -  u_s  + 𝒦(
    ϕ_s + τ_s n_s),  
    
    
    β  A_/ t + ϵ_s   _s u_s /
      t =     - ( ϕ_s + τ_s n_s) + 2 ϵ_sτ_s𝒦( u_s ) - C J_,


The gyrofluid moments are coupled by the polarisation equation

    ∑_s a_s[ Γ_1 n_s + Γ_0 - 1/τ_sϕ] = 0,

and Ampere's law

    - ^2 A_ = J_ = ∑_s a_s u_s .


The gyroscreened electrostatic potential acting on the ions is given by

    ϕ_s = Γ_1( ρ_s^2 k_⊥^ 2) ϕ_k,

where ϕ_k are the Fourier coefficients of the electrostatic potential.
The gyroaverage operators Γ_0 (b) and Γ_1 (b) = Γ_0^1/2 (b) correspond to
multiplication of Fourier coefficients by I_0(b) e^-b and 
I_0(b/2) e^- b/2, respectively, where I_0 is the modified Bessel
function of zero'th order and b = ρ_s^2 k_⊥^ 2. 
We here use approximate Padé forms with Γ_0 (b) ≈ (1 +
b)^-1 and Γ_1 (b) ≈ (1 + b/2)^-1 <cit.>. 

The perpendicular E×B advective and the parallel derivative operators
for species s are given by 

    _s/ t = / t + δ^-1{ϕ_s,  },


    = / z - δ^-1β{ A_,  },

where we have introduced the Poisson bracket as

    { f, g } = (  f/ x g/ y -  f/ y g/ x).



In local three-dimensional flux tube co-ordinates {x,y,z}, x is a (radial) flux-surface
label, y is a (perpendicular) field line label and z is the position along
the magnetic field line.
In circular toroidal geometry with major radius R, the curvature operator is given by

    𝒦 = ω_B( sin z  / x + cos z  / y),

where ω_B = 2 L_⊥ / R,
and the perpendicular Laplacian is given by

    ^2 = ( ^2/ x^2 + ^2/ y^2).

Flux surface shaping effects <cit.> in more general tokamak or
stellarator geometry on SOL filaments <cit.> are here neglected for simplicity.

Spatial scales are normalised by the drift scale ρ_0 = √(T_
  m_i0)/ B, where T_ is a reference electron temperature, B is
the reference magnetic field strength and m_i0 is a reference ion mass,
for which we use the mass of deuterium m_i0 = m_D. 
The temporal scale is set to by c_0 / L_⊥, where c_0 = √(T_/m_i0),
and L_⊥ is a perpendicular normalisation length (e.g. a generalized
profile gradient scale length), so that δ = ρ_0 / L_⊥ is the drift scale. 
The temporal scale may be expressed alternatively L_⊥ / c_0 = L_⊥ /
(ρ_0 Ω_0) = (δΩ_0)^-1, with the ion-cyclotron
frequency Ω_0 = ρ_0 / c_0.  In the following we employ δ = 0.01.

The main species dependent parameters are

    a_s = Z_s n_s0/n_ 0 ,   τ_s = T_s/Z_s T_,   μ_s = m_s/Z_s m_i0, 
    ρ_s^2 = μ_sτ_sρ_0^2,   ϵ_s = μ_s( q R/L_⊥)^2,

setting the relative concentrations, temperatures, mass ratios and FLR scales
of the respective species. Z_s is the charge state of the species s with
mass m_s and temperature T_s. 

The plasma beta parameter

    β = 4 π p_/B^2( q R/L_⊥)^2,

controls the shear-Alfvén activity, and

    C = 0.51 ν_ L_⊥/c_0m_/m_i 0( q R/L_⊥)^2,

mediates the collisional parallel electron response for Z=1 charged hydrogen
isotopes. The collisional response for other isotopes or ion species is discussed
further below.



 §.§ Parallel boundary conditions


We distinguish between two settings for parallel boundary conditions in 
3-d simulations. In the case of edge simulations a toroidal
closed-flux-surface (CFS) geometry is considered, and quasi-periodic globally
consistent flux-tube boundary conditions in the parallel direction
<cit.> are applied on both state-variables n_, ϕ and flux
variables v_, u_s.   

For SOL simulations, the state variables assume zero-gradient Neumann (sheath)
boundary conditions at the limiter location, and the flux variables are given as 

    u_s |_±π   = p_|_±π = ±Γ_d n_|_±π,
    
    v_   = u_s |_±π - J_ |_±π  = ±Γ_d
    [(Λ + 1) n_|_±π - ϕ|_±π],

at the parallel boundaries z = ±π respectively <cit.>. 
Note that in order to retain the Debye sheath mode in this isothermal model, the Debye current  
J_ |_±π = ±Γ_d (ϕ - Λ T_) is expressed as
J_ |_±π = ±Γ_d (ϕ - Λ n_) and the electron
pressure p_ = n_ T_ is replaced by p_ = n_ <cit.>. 
This edge/SOL set-up and its effects on drift wave turbulence has been
presented in detail by Ribeiro in Refs. <cit.>.

The sheath coupling constant is Γ_d = √((1 + τ_i) / (μ_i ϵ)). 
The floating potential is given by Λ = Λ_0 + Λ_i, where
Λ_0 = log√(m_i0 / (2 π m_)) and Λ_i = log√(μ_i / (1 + τ_i)). 
Here terms with the index i apply only to the ion species. 

The expressions presented here are obtained by considering the finite ion
temperature acoustic sound speed, c_i = √((Z_i T_i + T_) / m_i), 
instead of c_0 in Ref. <cit.>. This results in the additional
Λ_i, and the normalisation scheme yields the extra 
√((1 + τ_i) / μ_i) in Γ_d. 




 §.§ Numerical implementation


Our code TOEFL <cit.> is based on the delta-f isothermal electromagnetic
gyrofluid model <cit.> and uses globally consistent flux-tube
geometry <cit.> with a shifted metric treatment of the coordinates
<cit.> to avoid artefacts by grid deformation. In the SOL region a
sheath boundary condition model is applied <cit.>. 
The electrostatic potential is obtained from the polarisation equation by an
FFT Poisson solver with zero-Dirichlet boundary conditions in the (radial)
x-direction. Gyrofluid densities are adapted at the x-boundaries to ensure
zero vorticity radial boundary conditions for finite ion temperature.  
An Arakawa-Karniadakis scheme is employed for advancing the moment
equations <cit.>. 






§ SCALING LAWS FROM DIMENSIONAL ANALYSIS
 


Blob velocity scalings are commonly deduced from the fluid vorticity
equation. We follow this approach and construct the gyrofluid vorticity
equation to deduce velocity scaling laws. 
The vorticity equation can be obtained upon expressing the gyrocenter ion
density in terms of the electron density and polarisation contribution,  
inserting in the ion gyrocenter density evolution equation and subtracting
the electron density evolution equation <cit.>. 
Up to 𝒪(b) the ion gyrocenter density is

    n_i = n_ - 1/2μ_i^2 p_i - μ_i^2 ϕ,

where the ion pressure is given in terms of the electron particle density p_i = τ_i n_. The gyroaveraged potential for species s up to 𝒪(b) is

    ϕ_s = ϕ +1/2μ_iτ_i^2 ϕ.

Following Ref. <cit.> we obtain

    μ_i∇·/ tϕ^* =  J_ -
    (1 + τ_i) 𝒦 (n_).

Here we have introduced the modified potential ϕ^* = ϕ + p_i.
The vorticity equation is equivalent to the quasi-neutrality statement of
current continuity, ∇·J = 0.  
We identify the divergence of the polarisation current,

    ∇·J_pol = - μ_i ∇·/ tϕ^* ,

and the divergence of the diamagnetic current,

    ∇·J_dia = -(1 + τ_i) 𝒦 (n_).

Blob propagation has in a linearisation of the present gyrofluid model been
analytically analysed by Manz in Ref. <cit.>. 
Therein the dependence of blob velocity 
on the ion isotope mass is in principle present but not explicitly apparent. 
To clarify, we here restate the calculations of Ref. <cit.>, but use the
vorticity eq. (<ref>) with the explicit occurence of μ_i. 
Neglecting parallel currents, employing the blob correspondence
∂_x,∂_y → 1 / σ and /  t →
i γ = i v_b / σ <cit.>, in terms of the blob width σ,  
blob velocity v_b and linear growth rate of the instability γ and
furthermore identifying ϕ = i v_b σ (the radial component of the
electric drift), we get (in normalised units):

    v_b = 1/√(2)√(√(f^2 + g^2) - f),  where  f =
    ( τ_i A/2 σ)^2,  g = 1 + τ_i/μ_iω_B σ A.

A is the initial blob amplitude. In the limit of large blobs, g ≫ f, so that

    v_b ≈1/√(2)√(1 + τ_i/μ_iω_B σ A),

and for smaller blobs satisfying g ≪ f we get

    v_b ≈1 + τ_i/τ_iω_B σ^2/μ_i.

The correspondence with the result of Ref. <cit.> is made explicit upon
renormalising, i.e. letting v_b → v_b / δ c_0, A → A
/ δ, σ→σ / ρ_0, ω_B → 2 L_⊥ / R. 
The limits then are

    v_b ≈√(1 + τ_i/μ_i) c_0 √(2 σ/R A),
     for σ^3 ≫μ_i τ_i^2 A/4 (1 + τ_i) ω_B,

and

    v_b ≈ 2 c_0/μ_i1 + τ_i/τ_i(
    σ/ρ_0)^2,  for σ^3 ≪μ_i
      τ_i^2 A/4 (1 + τ_i) ω_B.

For 2-d computations of sufficiently large blobs we consequently expect v_b
∼ 1 / √(μ_i), whereas for the 3-d model the expected scaling is not
a priori that clear. In Ref. <cit.> 3-d (linear) scaling laws were presented,
where the parallel dynamics was approximated by the Hasegawa-Wakatani closure,
J_ = 1/C^2 (n_ - ϕ). 

In the following we are going to compare reduced 2-d and full 3-d dynamical
blob simulations for various isotope species with the analytical 1 / √(μ_i)-scaling.






§ TWO-DIMENSIONAL BLOB COMPUTATIONS
 



In this section we numerically analyse the dependence ob filament dynamics on
the normalised ion mass μ_i by 
reduced 2-d blob simulations of the isothermal gyrofluid eqs. (<ref>,<ref>).
For the computations in this section we us as parameters:
curvature ω_B = 0.05, drift scale δ = 0.01, grid size L_y = 128
ρ_s, grid points N_x = N_y = 256, initial blob amplitude A = 1 and 
Gaussian blob width σ = 10 ρ_s.  

Fig. <ref> shows contours plots of the electron particle density at
different times of evolution of a seeded blob for several species of cold ions (τ_i=0). 
The initial Gaussian density perturbation n_e (x,y, t=0) = A exp[ -
  (x^2 + y^2) /  σ^2 ] undergoes the familiar transition towards a
mushroom-shaped structure before the blob eventually breaks up due to 
secondary instabilities. This figure illustrates the main point for
the following discussion: lighter isotopes propagate faster than heavier
isotopes. In terms of the (normalising) deuterium mass we consider
μ_H = 1/2, μ_D = 1, μ_T = 3/2 and
μ_He+ = 2 with μ_i = m_i /(Z_i m_D). The species index 
He+ here denotes singly charged
helium-4 with Z_He+≡ 1. The case of (fully ionised) doubly charged helium-4 isotopes
will be discussed further below in context of 3-d simulations in Sec. <ref>.    
Note that the lighter the ion species are, the further the blob is developed
in its radial propagation and evolution at a given snapshot in time.  





For warm ions (τ_i >0) the blob propagation depends on the relative
initialisation of the electron and ion gyrocenter densities.

Commonly, a zero E× B vorticity blob initialisation is assumed where 
n_i (x,y, t=0) = Γ_1^-1 n_ (x,y, t=0): inserting these into the polarisation
eq. (<ref>) results in vorticity Ω≡∇_⊥^2 ϕ = 0. This
initialisation for most parameters leads to an FLR induced rapid development of a
perpendicular propagation component in addition to the radial propagation of
the blob, and thus a pronounced up-down asymmetry in y direction.

Alternatively, the electron and ion gyrocenter densities can be chosen as
equal with n_i (x,y, t=0) = n_e (x,y, t=0), so that ϕ∼ n_e. 
In this case the initial vorticity mostly cancels the FLR asymmetry effect,
and the blob remains more coherent and steady in its radial propagation <cit.>.

The truth may be somewhere in between: as in the experiment blobs are not
“seeded” (in contrast to common simulations), but appear near the separatrix
from E× B drift wave vortices or are sheared off from
poloidal flows, in general some phase-shifted combination of electric
potential and density perturbations will appear. 
For comparison we perform simulations with both of these seeded blob density initialisations.

We note that the x coordinate is effectively pointing radially outwards (in
negative magnetic field gradient direction) at a low-field midplane location
in a tokamak, and the magnetic field here points into the (x,y) plane (e_z
=  e_y × e_x), so that the effective electron diamagnetic
drift direction of poloidal propagation is in the present plots downwards (in
negative y direction).  






Fig. <ref> shows blob propagation for warm ions with τ_i = 1,
initialised with the zero vorticity condition. 
For comparison, we present in Fig. <ref> the propagation for the same
parameters but initialised with equal electron and ion gyrocenter densities. 
Clearly, the latter cases with initial non-zero vorticity  Ω = ∇_⊥^2
ϕ results in faster and more coherent radial propagation, whereas the
zero vorticity cases exhibit significant poloidal translation through the FLR
induced spin-up.
Regardless of initialisation, blobs of light ion species with small μ_i
travel faster and are further developed at a given time compared to heavier species.  





Relevant quantities which determine the intermittent blob related transport properties
of the tokamak SOL are the maximum blob velocity and acceleration.
In Fig. <ref> we present maximum radial center-of-mass velocities
V_x, max and the average radial acceleration A as a function of
the the ion mass parameter μ_i.
The different symbols/colours represent cases with  cold (τ_i = 0, blue
lower curves) and warm  (τ_i = 1) ions, with both types of initial
conditions used on the latter: the zero vorticity condition is depicted in red
(middle curves) and the n_e=n_i  condition in green (upper curves).

It can be seen that the maximum radial blob velocity is slightly larger for n_e=n_i
initialisation due to the mainly radial propagation (left figure), but the
average acceleration is for both cases nearly equal (right figure).

The radial center-of-mass position is given by X_c = [ ∫ x  y  
  x n_ ] / [ ∫ x  y   n_ ].
Taking the temporal derivative gives the radial center-of-mass
velocity, V_x =  X_c /  t. 
The maximum of V_x (t), V_x, max = max{ V_x (t) }
and the corresponding time for the occurence of the maximum, t_max,
then give a measure of the average radial acceleration, 
A = V_x, max / t_max.







Clearly, an inverse dependence of velocities and acceleration on effective
ion mass μ_i can be inferred for all cases. 
For cold ions, the only mass dependence in the present 2-d isothermal
gyrofluid model, lies in the gyrofluid polarisation equation, carrying over
the mass dependence of the polarisation drift in a fluid model. 
As deduced from the basic linear considerations in Sec. <ref>, the
maximum blob velocity scales inversely with the square root of the ion species
or isotope mass: the plotted fits are close to the expected lines V_x,
  max∼μ^-0.5.

From dimensional analysis it follows that the acceleration should scale
according to A ∼γ^2 σ, where γ is the growth-rate 
of the linear instability. For γ∼ 1 / √(μ_i), we expect 
A ∼ 1 / μ_i for cold ions: this is confirmed in Fig. <ref>
(right) where the fitted exponents are close to -1. 
Warm ion simulations also feature μ_i species mass dependence through the FLR
  operators Γ_0(b) and Γ_1(b), where b = ρ_i^2 k_⊥^2 =
  μ_i τ_i ρ_0 k_⊥^2. 


We find that for the parameters at hand, the maximum radial velocity for warm
ions with zero vorticity initialisation is higher compared to cold ions, with
a slightly increased isotopic dependence (seen in an exponent -0.548
compared to -0.500). 

Initialising with non-zero vorticity yields approximately 50 % increased
velocities compared to cold ions, and slightly weakens the isotopic dependence
(expressed by an exponent -0.488). This can be attributed to the
mass dependence in the FLR operators, which we further discuss below in Sec. <ref>.  




§ THREE-DIMENSIONAL FILAMENT COMPUTATIONS



In three dimensions, when the blob extends into an elongated filament along
the magnetic field lines, additional physics enters into the model. 
The basic picture of interchange driving of filaments by charging through ∇
B and curvature drifts to produce a net outward E×B 
propagation still remains valid. However, the total current continuity balance
now also involves parallel currents: 
- ∇·J_pol =  J_ + ∇·J_dia. 

The detailed balance among the current terms determines the overall motion of
the filament. Furthermore, blob filaments in the edge of toroidal magnetised
plasmas generally tend to exhibit ballooning in the unfavourable curvature
region along the magnetic field.  
The parallel gradients in a ballooned blob structure also lead to a parallel
Boltzmann response, mediated mainly through the resistive coupling of (ϕ -
n_) to C J_|| in eq. (<ref>).  
This tends towards (more or less phase shifted) alignment between the electric
potential and the perturbed density, which strongly depends on the
collisionality parameter C.  

For low collisionality, the electric potential in the blob evolves towards
establishment of a Boltzmann relation in phase with the electron density along B, so
that n_∼exp(- ϕ) ∼ϕ. This leads to reduced
radial particle transport, and the resulting spatial alignment of the
potential with the blob density perturbation produces a rotating vortex along
contours of constant density, the so-called Boltzmann spinning <cit.>. 
Large collisionality leads to a delay in the build-up of the potential within
the blob, so that the radial interchange driving can compete with the parallel
evolution, and the perpendicular propagation is similar to the 2-d scenario.

In the following we investigate how 3-d filament dynamics is depending on the
ion mass. Clearly, we expect an impact in addition to the 2-d effects found in
the previous section, since (i) the parallel ion velocity is inversely
dependent on ion mass (but is for any ion species slow compared to the
electron velicity), (ii) the sheath boundary coupling constants are mass
dependent, and (iii) the basic dependence on the ion mass in the polarisation
current will play a more role complicated role compared to the 2-d model.  

For our present study we chose the free computational parameters
basically identical to the 2-d case above: drift scale δ = 0.01,
curvature ω_B=0.05, blob amplitude A=1 and perpendicular blob width σ=10.
The Gaussian width of initial parallel density perturbation is given by
Δ_z=√(32), which represents a slight ballooning with some initial
sheath connection: 

    n_^3D (t = 0, x, y, z) =  n_e ⊥·exp[-(z - z_0)^2/Δ^2]

where z_0 is the parallel reference coordinate at the outboard mid-plane and
n_e ⊥(x,y) is the perpendicular Gaussian initial perturbation
introduced in Sec. <ref>.  
In this section we first focus on zero vorticity initial conditions, non-zero
conditions will be discussed further below.

The perpendicular domain size is L_x = L_y = 128 ρ_s with a grid
resolution of N_x = N_y = 256. The number of parallel grid points is varied
between N_z=8 and 16. The filament simulations have been tested for convergence
with respect to the number of drift planes up to N_z =32: N_z = 8 yields
qualitatively and quantitatively similar results to N_z = 16. 
The plots showing colour cross sections throughout this article are taken from
simulations with N_z = 8, and the presented quantitative results have been
obtained for N_z = 16.  

We here set β = 0 as electromagnetic effects in the SOL are thought to
be of minor importance for the present discussion <cit.>. 
The collisionality parameter is chosen in C = 0.5 - 100 to cover a
likely range of tokamak SOL values. 
Typical values for the collisionality parameter for the SOL in ASDEX Upgrade
L-mode plasmas have in the literature <cit.> been reported as C ∼ 1
- 100, and a reference characteristic collisionality in Ref. <cit.> 
for MAST has been given as C ∼ 2. 




Fig. <ref> illustrates the dependence of filament evolution
with respect to collisionality dependent Boltzmann spinning for warm deuterium ions. 
The case with C = 10 represents the strong Boltzmann spinning (drift wave)
regime, where density and potential perturbations are closely aligned and the
radial filament motion is strongly impeded. 
Increasing the collisionality to C = 100 reduces parallel electron dynamics
and so effectively increases the lag of potential build-up within the density
blob perturbations, so that the Boltzmann spinning is reduced. 
The substantial perpendicular motion component is in this case partly caused
by FLR effects like in the corresponding 2-d case for τ_i=1.
The τ_i contribution to the ion diamagnetic curvature term results in 
enhanced radial driving of the blob compared to cold ion cases.

A measure for blob compactness can be introduced <cit.> by

    I_c (t) = ∫ x ∫ y n_ (x,y,t) h(x,y,t)/∫ x
      ∫ y n_ (x,y,t=0) h(x,y,t=0),

where the Heavyside function h(x,y,t) is defined as

    h (x,y,t) = 1    if  (x - x_max (t))^2 + (y -
    y_max (t))^2 < σ^2,

and zero elsewhere. That is, the integral takes non-zero values for density
contributions located inside a circle of radius σ around its maximum. 
Fig. <ref> quantifies the above observations. In the strong
Boltzmann spinning regime (C = 10, left) the blob retains much of its
initial shape, so that the compactness is higher compared to the weak
Boltzmann spinning regime (C = 100, right), where filaments feature a more 
bean-shaped structure which reduces the compactness measure. 
At the time of measurement (t = 5), the heavier isotopic blobs show slightly more
compactness, which is an indirect result of decreased velocity: at a given
time, the lighter isotopic blobs are further developed and thus less circular. 
The observed trends are similar for cold (τ_i=0) and warm (τ_i=1) ions.

For increased collisionality, the deviation from circularity is more pronounced, as
the mushroom-cape shape is realized. Blobs in light isotopic plasmas are then
again further  developed, i.e. finer scales have emerged at the time of
recording, resulting in a sharper mass dependence of blob compactness compared
to C = 10, where smaller scales are less prominent. 












Fig. <ref> shows filament propagation for cold ions (τ_i=0) and weak
Boltzmann spinning. This can be compared to Fig. <ref> which shows
propagation for warm ions (τ_i=1) and also weak Boltzmann spinning. 
It is observed that there is poloidal propagation also for the cold ion case,
which is a consequence of the non-vanishing Boltzmann spinning that is also present, although greatly reduced, for these high collisionality (C = 100) cases. 
The resulting maximum radial center-of-mass velocites at the outboard midplane
(z = z_0) are shown in Fig. <ref> for weak (right) and strong (left)
Boltzmann spinning. 








The fits of the exponent in μ^α to the simulation data in
Fig. <ref> carry evidence that the additional mass dependences
introduced by the 3-d model via parallel sheath-boundary conditions and
parallel ion velocity dynamics causes the clear deviation from a 1 /
√(μ_i) scaling.  

For high collisionality (and thus reduced Boltzmann spinning), the parallel
current is impeded and the dynamics is more two-dimensional than for lower
collisionalities. The competing nature of the parallel divergence versus
current continuity via the divergence of the polarisation current with
collisionality  is shown in Fig. <ref>:
for each value of collisionality C we compute the isotopic dependence of the
outboard-midplane maximum center-of-mass radial velocity, 

    V_max∼μ_i^α (μ_i),

contained in the scaling exponent, α (μ_i). 
For large values of C the resulting dynamics strongly features
2-d propagation characteristics, since the diamagnetic current is almost
exclusively closed via the polarisation current, which gives the 
1 / √(μ_i) scaling introduced in sec. <ref>. 
Note that the scaling with respect to C cannot be inferred from linear models.



  §.§.§ Non-zero gyrofluid vorticity initialisation
 


So far we have in this section applied the initial condition n_≡Γ_1 n_i
associated with zero initial vorticity. Now the case for non-zero initial
vorticity by the condition n_ = n_i is considered. 

Fig. <ref> shows results for warm ion (τ_i = 1) computations in the strong  
(C = 10, blue) and weak (C = 100, green) Boltzmann spinning regimes.  
(Recall that for τ_i = 0 this discussion is redundant since Γ_1 (τ_i = 0)= 1.)
The left figure depicts the maximum radial center-of-mass velocity, and the right
figure shows the corresponding average acceleration.

Comparing with Fig. <ref> we notice that the resulting filament velocities
are similar to those obtained from zero initial vorticity. We also
find that the isotopic dependence ∼μ_i^α is not significantly altered. 






Recalling the results from 2-d computation in sec. <ref>, we may
conclude that the initialisation is not that important for 3-d numerical
simulations with respect to the maximum radial filament velocity. 
The slight impact of the initial condition on the resulting scaling exponent
for the 2-d case may then be connected to the more prominent mass
dependence in the polarisation current, which is weakened when
parallel currents are taken into account. 




  §.§.§ Comparison of filaments in deuterium and helium plasmas
 


When comparing blob filament propagation in deuterium and in fully ionised
helium-4 plasmas in the present model, the dynamical evolution is identical in
the cold ion limit: in the model parameter μ_i = m_i /(Z_i m_D) the
doubled mass of the helium nucleus excactly cancels with the doubled positive
charge, Z_He = 2. 

Differences are only appearing in warm ion cases. The normalised 
mass ratio is now identical for both species, μ_D = μ_He =
1. The only model parameter that is different, is the helium temperature
ratio, τ_He = T_He / Z_He T_. The species mass
effects thus appears in the combined b ∼μ_i τ_i.
In the following we consider plasmas at equal temperature, T_D =
T_He = 2 T_ such that τ_D = 2 and τ_He = 1. 

The higher charge state of the helium nucleus is also indirectly evident in
the reduced electron-ion collision frequency contained in the C -parameter. 
For electron-ion collisions where the ions are in charge state Z_i we have C
∼α_e ν_, with <cit.> 

    α_e ≈1 + 1.198 Z_i + 0.222 Z_i^2/1 + 2.966 Z_i + 0.753 Z_i^2.

For Z_i = 1 we have α_e ≈ 0.51 and Z_i = 2 gives α_e
≈ 0.43. 

To account for this dependence, we in the following consider two cases: 
(1) equal non-normalised collision frequencies, i.e. C → 0.51 C
for deuterium and C → 0.43 C for helium; (2) using the same C
for both deuterium and helium computations.  

In case (2), setting first the normalised collisionality parameter C = 10
identical for both D and He computations results in V_D = 1.41
δ c_s0 and V_He = 1.14 δ c_s0. 
Setting C = 100 identical for both D and He gives V_D = 3.2 δ
c_s0 and V_He = 2.6 δ c_s0, respectively. 

For case (1) we set the electron-ion collision frequency equal for both
species, so that different C parameters are used according to equation <ref>:
C_D = 10 corresponds to C_He = 8.43, and 
C_D = 100 to C_He = 84.3. 
In these cases, maximal radial He velocities are V_He (C=10) = 1.09
δ c_s0 and V_He (C=100) = 2.46 δ c_s0. 

We find that regardless of how the charge state depency for the  
relative value of the collision parameter is treated, the filaments in
deuterium plasmas move faster than in helium plasmas at identical temperature. 
This is visualized in Fig. <ref> showing filament propagation at equal
electron-ion collision frequency and electron temperature. 






§ CONCLUSIONS
 


We have investigated filament propagation in SOL conditions characteristic for
tokamak fusion devices. Quasi-2-d dynamics is restored in high resistivity
regimes, where the maximum radial blob velocity scales inversely proportional
with the square root of the ion mass. 
In 2-d the diamagnetic current drive is closed solely via the polarisation current,  
yielding this simple characteristic scaling. 

The larger inertia through polarisation of more massive ion species
effectively slows the evolution of filaments, and the maximum radial velocity
occurs later compared to blobs in plasmas with lighter ions.

For non-zero initial vorticity condition, the 2-d warm ion blobs show compact
radial propagation, where the isotopic effect through the mass dependent
FLR terms is slightly less pronounced.

Boltzmann spinning appears in 3-d situations particularly for low
collisionality regimes, and leads to a reduced dependence on the ion isotope
mass. The exponent in the scaling V ∼μ_i^α has been found to be
typically within the range α∈ [-0.1, -0.3] for C < 10, which is a
regime relevant for the edge of most present tokamaks. 

Considering current continuity, the closure via the parallel current
divergence dynamically competes with current loops being closed through the
polarisation current.
For high collisionalities the parallel current is effectively imepeded and 
the polarisation current characteristics dominate the blob evolution,
producing a more 2-d like velocity dependence with respect to ion mass. 
The initial condition has been found to have little influence on the maximum
radial velocity when in 3-d the parallel closure of the current is taken into account. 

For similar ion temperatures and electron-ion collision frequencies, it has
been found that helium filaments travel more slowly compared to deuterium
filaments in both high and low collisionality regimes.

This work was devoted to the identification of isotopic mass effects on
seeded (low amplitude) blob filaments in the tokamak SOL by means of a
delta-f gyrofluid model. 
Naturally, blobs emerge near the separatrix within coupled edge/SOL
turbulence. The dependence of fully turbulent SOL transport on the ion mass
therefore will have to be further studied within a framework that consistently
couples edge and SOL turbulence, preferably through a full-f 3-d gyrofluid
(or gyrokinetic) computational model that does not make any smallness
assumption on the relative amplitude or perturbations compared to the
background <cit.>. 





§ ACKNOWLEDGMENTS

We acknowledge main support by the Austrian Science Fund (FWF) project Y398.
This work has been carried out within the framework of the EUROfusion
Consortium and has received funding from the Euratom research and training
programme 2014-2018 under grant agreement No 633053. The views and opinions
expressed herein do not necessarily reflect those of the European Commission.  







§ REFERENCES


00

Bessenrodt
Bessenrodt-Weberpals M 1993, 33 1205

Hawryluk
Hawryluk R J 1998, Rev. Mod. Phys. 70 537

Liu
Liu B 2016, 56 056012

Xu
Xu Y 2013, 110 265005

Boedo
Boedo J A 2003, Phys. Plasmas 10 1670

LaBombard
LaBombard B 2001, Phys. Plasmas 8 2107

Theiler
Theiler C 2009, 103 065001

DIp_review
D'Ippolito D A, Myra J R and Zweben S J 2011, Phys. Plasmas 18 060501

Krash1
Krasheninnikov S I 2001, Phys. Lett. A 283 368

Krash
Krasheninnikov S I, D'Ippolito D A and Myra J R 2008, J. Plasma Phys. 74, 679

Hahm
Hahm  T S 2013 53 072002

Manz
Manz P 2015, Phys. Plasmas 22 022308

Madsen
Madsen J 2011, Phys. Plasmas 18 112504

Matthias
Wiesenberger M, Madsen J and Kendl A 2014, Phys. Plasmas 21 092301

paper1
Meyer O H H and Kendl  A 2016, 58 115008

scott05b 
Scott B D 2005 Phys. Plasmas 12 102307

dorland93
Dorland W 1993 B 5.3 812–835 

kendl06
Kendl A and Scott B D 2006 Phys. Plasmas 13 012504

kendl03
Kendl A, Scott B D, Ball R and Dewar R L 2003 Phys. Plasmas 10 3684

riva17
Riva F, Lanti E, Jolliet S and Ricci R
2017 Plasmas Phys. Contr. Fusion 59 035001

scott98
Scott B D 1998, Phys. Plasmas 5 2334

Ribeiro05
Ribeiro TT and Scott BD
2005 Plasmas Phys. Contr. Fusion 47 1657

Ribeiro08
Ribeiro TT and Scott BD
2008 Plasmas Phys. Contr. Fusion 50 055007

kendl14
Kendl A 2014 Int. J. Mass Spectrometry 365/366 106–113 

scott01
Scott B D 2001, Phys. Plasmas 8 447

arakawa66
Arakawa A 1966 J. Comput. Phys. 1 119

karniadakis91
Karniadakis G E, Israeli M and Orszag S A 1991 J. Comput. Phys. 97 414

Naulin
Naulin V and Nielsen A 2003, J. Sci. Comput. 25 104

Scott07
Scott B D 2007, Phys. Plasmas 14 102318

Manz13
Manz P 2013, Phys. Plasmas 20 102307

Krash_bc
Krasheninnikov S I 2001, Phys. Lett. A 283 368

Held16
Held 2016, 56 126005

angus12
Angus J R et al 2012 Contrib. Plasma Phys. 52 348

angus14
Angus J R and Umansky V M 2014 Phys. Plasmas 21 012514

Easy1
Easy L 2014, Phys. Plasmas 21, 122515

Hirshman
Hirshman S P 1977, Phys. Fluids 20 589

Kendl15
Kendl A 2015, 57 045012






