# Turbulent transport regimes in the presence of an X-point magnetic
configuration

Fabio Moretti  Nakia Carlevaro  Francesco Cianfrani  Giovanni Montani  ENEA,
Nuclear Department, C.R. Frascati, Via E. Fermi 45, 00044 Frascati (Roma),
Italy  Physics Department, “Sapienza” University of Rome, P.le Aldo Moro 5,
00185 Roma, Italy

###### Abstract

We analyze the transport properties of the two-dimensional electrostatic
turbulence characterizing the edge of a Tokamak device from the study of test
particles motion (passive fluid tracers) following the
$\mathbf{E}\times\mathbf{B}$ drift. We perform statistical tests on the tracer
population in order to assess both the magnitude and the main features of
transport. The role of other physical properties, such as viscosity and
inverse energy cascade in the spectrum, is also considered. We outline that
large scale eddies are responsible for greater transport coefficients, while
the presence of an X-point magnetic field reduces the mean free path of the
particles, however generating a larger outliers population with respect to a
Gaussian profile.

##  1 Introduction

The electrostatic turbulence, characterizing the plasma in the edge region of
a Tokamak device [1, 2, 3], is a fundamental ingredient to determine the
transport properties in a realistic fusion configuration [4]. The first model
able to capture some basic features of the non-linear drift response is the
so-called Hasegawa-Wakatani scenario [5, 6, 7] (see [8, 9, 10, 11] for
extensions), which identifies the source of the turbulent dynamics in the non-
linear coupling between pressure perturbation and electric field fluctuations
[12, 13]. As discussed in [9, 11, 14, 15], the two-dimensional turbulence,
preserving the axisymmetry of the Tokamak configuration, plays an important
role and, to some extent, is an attractor for the three-dimensional non-linear
dynamics.

Here we investigate the influence that a two-dimensional turbulent plasma has
on passive fluid tracers moving inside it according to the
$\mathbf{E}\times\mathbf{B}$ drift velocity, following the same line of
investigation started in [16]. The peculiar and interesting feature of our
analysis consists in the comparison of different physical settings, i.e. the
slab magnetic configuration versus the presence of an X-point profile.
Furthermore, we also distinguish two different turbulent regimes,
corresponding to the presence in their spectral morphology of a condensation
phenomenon (associated to an inverse energy cascade [17, 18]) or in its
absence (when a direct enstrophy cascade dominates). The slab case is also
investigated in the limit of an ideal plasma, in order to understand the role
of the thermodynamical equilibrium reached by the system in this case, as
discussed in [19]. Such an ideal limit will also offer us the possibility of a
comparison of the diffusion coefficient we are able to recover from the tracer
evolution and the corresponding theoretical value provided by the so-called
$K-\epsilon$ model [20, 21].

Specifically, in all these cases, we develop the tracer dynamics (i.e. charged
particles, which can also thought as plasma constituents) under the influence
of the $\mathbf{E}\times\mathbf{B}$ drift flow, defined by the non-linear
dynamics of the electric field living in the plasma. Then, we plot the mean
square displacement (MSD) of the particles with respect to their initial
position versus time. The comparison of the MSD evolution in the different
physical situations allows us to determine the different transport regimes,
within a single scenario and in comparison to other choices of the field set-
up.

The merit of the present study is to fix the following two relevant claims: i)
in the slab case, the presence of the condensation phenomenon for the
turbulence spectrum implies a larger value of the diffusion coefficient, where
it can be recovered, with respect to the case of a direct enstrophy cascade.
Nonetheless, there is an enhancement of the asymptotic value of the MSD for
this former case, estimated in about $50\%$, so that a marked transport
process is recovered. ii) in the comparison between the presence and absence
of an X-point configuration, typical of a real Tokamak device [22], we notice
a smoother transport process, i.e. a smaller asymptotic value of the MSD.

The first result can be interpreted by noticing that the formation of large
sized eddies, in the case of a condensation of the spectrum, causes a growth
of the mean free path for tracers between interactions with vortices and this
results, on average, in a greater amount of transport. The presence of the
X-point and then the emergence of shear magnetic contributions is instead at
the ground of the second claim. In fact, such shear terms reduce the
turbulence of the fluctuating fields and favor a smaller asymptotic MSD.
However, it has to be stressed that, looking at the statistical distribution
of the tracers, we see that, although the variance of the distribution is
consistently smaller than the one in the slab case, the kurtosis of such a
distribution is correspondingly larger. This feature suggests that the X-point
configuration reduces, on average, the transport process, but enhances the
number of outliers that can efficiently escape from the initial position.

It is also remarkable the agreement that we find between the diffusion
coefficient obtained from the tracers analysis and the one we estimate from
the $K-\epsilon$ model, in the case of a slab configuration. This fact
suggests that our investigation is validated in its predictivity by the
capability to reproduce the right turbulence feature of the two-dimensional
model we adopted here. In this respect, we stress that, although we consider a
reduced model for the electric field dynamics (according to the analysis in
[9]), its non-linear features are just those of a fully self-sustained non-
linear drift response. By other words, our reduced model properly captures all
the important properties of the emerging turbulence and the confirmation of
its viability is just in the correct prediction of the diffusion coefficient,
that such a turbulent regime would generate according to the $K-\epsilon$
model.

The manuscript is organized as follows: in Section 2, we present the model of
electrostatic turbulence on which we base our analysis, together with its
reduced version in the toroidal field (slab) case; in Section 3, we describe
the methods implemented to study the tracers transport; in Section 4, we
display the results obtained. Conclusions are drawn in the final section.

##  2 Model of electrostatic turbulence with an X-point magnetic configuration

Here we provide a presentation of the model describing the turbulent dynamics
due to the non-linear low frequency drift response, in the presence of a
background magnetic configuration able to mimic the features of the X-point
region of a Tokamak device. We start by illustrating the local magnetic
equilibrium on which we will subsequently set the analysis of electrostatic
fluctuations. Neglecting toroidal curvature effects, we adopt a set of
Cartesian coordinates $(x,y,z)$ centered in the X-point, i.e. the origin of
our flat grid corresponds to the null of the poloidal magnetic field. The
coordinates $(x,y)$ span the poloidal plane, while $z$ is taken along the
toroidal direction. The magnetic field $\mathbf{B}$ is expressed as

$\mathbf{B}=B_{p}y\mathbf{\hat{e}_{x}}+B_{p}x\mathbf{\hat{e}_{y}}+B_{t}\mathbf{\hat{e}_{z}},$
(1)

where $B_{p}$ and $B_{t}$ are two constants accounting for the strength of the
poloidal and toroidal magnetic field components, respectively. By introducing
the poloidal length scale $L_{p}$, the request of having the magnetic field
toroidal component as the dominant contribution is translated in the condition
$B_{t}\gg B_{p}L_{p}$, which has to be satisfied in all points in which the
magnetic configuration (1) results valid. Further, it is possible to define
the directional versor $\\!\\!\hat{\,\,\textbf{b}}=\frac{\mathbf{B}}{B}$, in
which $B$ represents the coordinate-dependent magnitude of the magnetic field,
namely

$B(x,y)=B_{t}\sqrt{1+\frac{B_{p}^{2}\left(x^{2}+y^{2}\right)}{B_{t}^{2}}}.$
(2)

We stress that we always speak of parallel and perpendicular directions
referring to the versor $\\!\\!\hat{\,\,\textbf{b}}$. For instance, by setting
$B_{p}=0$ it turns out that the parallel direction corresponds to the toroidal
one (namely the $z$ axis), whereas by taking $B_{p}\neq 0$ the versor
$\\!\\!\hat{\,\,\textbf{b}}$ acquires a non-null component on the poloidal
plane. Moreover, the definition of the directional versor
$\\!\\!\hat{\,\,\textbf{b}}$ allows for an unambiguous splitting of the
gradient operator $\bm{\nabla}$ in a parallel and a perpendicular
contribution, i.e. $\bm{\nabla}=\bm{\nabla}_{\perp}+\bm{\nabla}_{\parallel}$,
where the parallel operator is defined as
$\bm{\nabla}_{\parallel}=\\!\\!\hat{\,\,\textbf{b}}(\\!\\!\hat{\,\,\textbf{b}}\cdot\bm{\nabla})$
while the perpendicular component is simply recovered from
$\bm{\nabla}_{\perp}=\bm{\nabla}-\bm{\nabla}_{\parallel}$.

Now we enumerate the hypotheses on which our model of electrostatic turbulence
is based. We deal with a Hydrogen-like plasma, which satisfies the conditions
of homogeneity and isotropy. We adopt a two-fluid description, assuming quasi-
neutrality (having observed that the turbulence length scale is always much
greater than the Debye length of the plasma we consider), i.e. we set
$\mathcal{N}_{i}=\mathcal{N}_{e}=\mathcal{N}$, being $\mathcal{N}_{i/e}$ the
ion/electron number density. This hypothesis, together with the further
assumption of equal temperature $T$ for ions and electrons, immediately
implies that we can consider a unique pressure for the two fluids, i.e.
$p_{i}=p_{e}=p=\mathcal{N}K_{B}T$, where we use the standard symbol $K_{B}$ to
denote the Boltzmann constant. Moreover, we take the ion polarization drift as
the dominant contribution to the orthogonal motion and we neglect both
diamagnetic effects and parallel ion velocity. Those just listed are the
ordinary assumptions that are necessary to derive the 3D Hasegawa-Wakatani
model [8] describing electrostatic edge turbulence. Our original contribution
is condensed in the last hypothesis, namely we neglect the gradients of the
background quantities. Hence, in this work, we consider $\mathcal{N}=const.$,
$p=const.$ and $T=const.$, so that all the dynamical variables are
perturbations denoted by a symbol $\delta$.

Let us now turn to the construction of the model. The first dynamical equation
we consider is the unique (due to the quasi-neutrality assumption) continuity
equation for both ions and electrons

$\frac{d\delta\mathcal{N}}{dt}-\mathcal{D}\nabla^{2}_{\perp}\delta\mathcal{N}=\frac{1}{e}\bm{\nabla}_{\parallel}\cdot\delta\textbf{J}_{\parallel}\,,$
(3)

where a diffusive term, controlled by the coefficient $\mathcal{D}$, has been
introduced to take into account different regimes of number density transport.
The symbol $\delta\textbf{J}_{\parallel}$ represents the fluctuating parallel
current density and $e$ is the electron charge (we adopt Gaussian units). In
this work we consider the $\textbf{E}\times\textbf{B}$ drift velocity as the
dominant contribution affecting the advective derivative, i.e.
$d/dt=\partial_{t}+\delta\textbf{v}_{\perp}^{e}\cdot\bm{\nabla}$, with

$\displaystyle\delta\textbf{v}_{\perp}^{e}=\frac{c}{B^{2}}\textbf{E}\times\textbf{B}=\frac{c}{B}\;\\!\\!\hat{\,\,\textbf{b}}\times\bm{\nabla}_{\perp}\delta\phi\,=$
$\displaystyle=\frac{c}{B^{2}}\big{[}\left(B_{p}x\partial_{z}\delta\phi-
B_{t}\partial_{y}\delta\phi\right)\hat{\textbf{e}}_{x}+\left(B_{t}\partial_{x}\delta\phi-
B_{p}y\partial_{z}\delta\phi\right)\hat{\textbf{e}}_{y}+$ (4)
$\displaystyle\qquad\qquad+\left(B_{p}y\partial_{y}\delta\phi-
B_{p}x\partial_{x}\delta\phi\right)\hat{\textbf{e}}_{z}\big{]}\,,$

being $c$ the speed of light and $\delta\phi$ the electric field potential.
Next, we write down the momentum balance equation for ions along the
perpendicular direction

$\frac{d\,\delta\textbf{v}_{\perp}^{i}}{dt}=\frac{e}{m_{i}}\left(-\bm{\nabla}_{\perp}\delta\phi+\delta\textbf{v}_{\perp}^{i}\times\textbf{B}/c\right)+\nu\nabla^{2}_{\perp}\,\delta\textbf{v}_{\perp}^{i}\,,$
(5)

where $m_{i}$ and $\nu$ are the ions mass and specific viscosity (here we
neglect the parallel component of the viscous stress). The latter can be
explicitly calculated as

$\displaystyle\nu=\frac{1}{m_{i}\mathcal{N}}\frac{(3/10)\;\mathcal{N}K_{B}T}{\Omega_{i}^{2}\tau_{ii}}\,,\qquad\qquad\tau_{ii}=\frac{3\sqrt{\mathcal{N}}(K_{B}T)^{3/2}}{4\sqrt{\pi}\;e^{4}\mathcal{N}\mathrm{ln}\Lambda_{ii}}\,,$
(6)

in which $\Omega_{i}=eB_{t}/cm_{i}$ is the ion gyro-frequency due to the
toroidal magnetic field, the Coulomb logarithm has been set to
$\mathrm{ln}\Lambda_{ii}=21$ and we have exploited the quasi-neutrality
condition. Due to the fact that we expect $\Omega_{i}$ to be much greater of
the turbulence time-scale inverse, we can assume
$\delta\textbf{v}_{\perp}^{i}=\delta\textbf{v}_{\perp}^{e}+\tilde{\textbf{v}}_{\perp}^{i}$,
with $\tilde{\textbf{v}}_{\perp}^{i}$ a small correction. The latter can be
calculated from a first-order expansion of (5), resulting in

$\tilde{\textbf{v}}_{\perp}^{i}=\frac{cm_{i}}{Be}\Big{(}\frac{d\,}{dt}-\nu\nabla^{2}_{\perp}\Big{)}(\,\\!\\!\hat{\,\,\textbf{b}}\times\delta\textbf{v}_{\perp}^{e})=-\frac{c^{2}m_{i}}{Be}\Big{(}\frac{d\,}{dt}-\nu\nabla^{2}_{\perp}\Big{)}\frac{1}{B}\bm{\nabla}_{\perp}\delta\phi\,.$
(7)

Hence, by splitting the charge conservation equation, namely
$\bm{\nabla}\cdot\delta\textbf{J}=0$, in its parallel and perpendicular
contributions, we obtain

$\bm{\nabla}_{\perp}\cdot\left(\frac{1}{B}\frac{d\,}{dt}\frac{1}{B}\bm{\nabla}_{\perp}\delta\phi\right)-\nu\bm{\nabla}_{\perp}\cdot\left(\frac{1}{B}\nabla_{\perp}^{2}\frac{1}{B}\bm{\nabla}_{\perp}\delta\phi\right)=\frac{1}{\mathcal{N}c^{2}m_{i}}\bm{\nabla}_{\parallel}\cdot\delta\textbf{J}_{\parallel}\,,$
(8)

having noticed that the orthogonal term can be written as
$\delta\textbf{J}_{\perp}\equiv\mathcal{N}e(\delta\textbf{v}_{\perp}^{i}-\delta\textbf{v}_{\perp}^{e})=\mathcal{N}e\,\tilde{\textbf{v}}_{\perp}^{i}$.
The electron momentum balance along the parallel direction results in the
generalized Ohm law

$\displaystyle\delta\textbf{J}_{\parallel}=\frac{\sigma}{\mathcal{N}e}\bm{\nabla}_{\parallel}\delta
p-\sigma\bm{\nabla}_{\parallel}\delta\phi\,,$ (9)

where we introduced the parallel conductivity $\sigma\equiv
1.96\,\mathcal{N}e^{2}/m_{e}\nu_{ei}$, with $\nu_{ei}$ the electron-ion
collision frequency and $m_{e}$ the electron mass. We stress that, given the
assumptions adopted in this analysis, the fluctuating pressure $\delta p$ can
be expressed in terms of the number density via the equation $\delta
p=K_{B}T\delta\mathcal{N}$.

The set (3), (8) and (9) constitutes the system of dynamical equations
characterizing our model. We introduce adimensional coordinates
$\tau\equiv\Omega_{i}t$, $u\equiv(2\pi/L_{p})x$, $v\equiv(2\pi/L_{p})y$,
$w\equiv(2\pi/L_{t})z$ ($L_{t}$ being the toroidal spatial scales, with
$L_{t}\gg L_{p}$), together with the parameters

$\displaystyle\varepsilon=\frac{B_{p}}{B_{t}}\,\left(\frac{L_{p}}{2\pi}\right),\qquad\gamma(u,v)=\sqrt{1+\varepsilon^{2}(u^{2}+v^{2})}$
(10)

accounting for the small ratio of the poloidal component with respect to the
toroidal one and the X-point geometry, respectively. The dimensionless version
of the parallel gradient vector is obtained from
$\textbf{D}_{\parallel}=(L_{p}/2\pi)\bm{\nabla}_{\parallel}$ and its explicit
expression reads

$\displaystyle\textbf{D}_{\parallel}=\Big{(}\frac{\varepsilon}{\gamma}(v\partial_{u}+u\partial_{v})+\frac{L_{p}/L_{t}}{\gamma}\,\partial_{w}\Big{)}\\!\\!\hat{\,\,\textbf{b}}\,.$
(11)

From the latter it is immediate to obtain the perpendicular gradient vector
$\textbf{D}_{\perp}$ and the Laplacian operators $D_{\parallel}^{2}$ and
$D^{2}_{\perp}$, whose lengthy and cumbersome explicit expressions are not
reported here to avoid unnecessary complexity. Setting $\Phi\equiv
e\,\delta\phi/K_{B}T$, $\bar{\mathcal{N}}=\delta\mathcal{N}/\mathcal{N}$ and
$\textbf{Y}_{\parallel}=(2\pi/L_{p})\delta\textbf{J}_{\parallel}/\mathcal{N}e\Omega_{i}$,
we summarize the dynamical equations of our model in dimensionless form as
follows

$\displaystyle\frac{d}{d\tau}\bar{\mathcal{N}}-\bar{\mathcal{D}}D^{2}_{\perp}\bar{\mathcal{N}}$
$\displaystyle=\bm{D}_{\parallel}\cdot\textbf{Y}_{\parallel}\;,$ (12)
$\displaystyle\alpha_{1}\bm{D}_{\perp}\cdot\left(\frac{1}{\gamma}\frac{d\,}{d\tau}\frac{1}{\gamma}\bm{D}_{\perp}\Phi\right)-\alpha_{1}\alpha_{2}\bm{D}_{\perp}\cdot\left(\frac{1}{\gamma}D_{\perp}^{2}\frac{1}{\gamma}\bm{D}_{\perp}\Phi\right)$
$\displaystyle=\bm{D}_{\parallel}\cdot\textbf{Y}_{\parallel}\;,$ (13)
$\displaystyle\alpha_{3}(D_{\parallel}^{2}\bar{\mathcal{N}}-D_{\parallel}^{2}\Phi)$
$\displaystyle=\bm{D}_{\parallel}\cdot\textbf{Y}_{\parallel}\,.$ (14)

The explicit expressions of the dimensionless constants appearing in our
equations are here reported:

$\displaystyle\alpha_{1}=\rho_{i}^{2}\left(\frac{2\pi}{L_{p}}\right)^{2},\quad\alpha_{2}=\frac{\nu}{\Omega_{i}}\left(\frac{2\pi}{L_{p}}\right)^{2},\quad\alpha_{3}=\frac{v_{A}^{2}\rho_{i}^{2}}{\eta_{B}\Omega_{i}}\left(\frac{2\pi}{L_{p}}\right)^{2},\quad\bar{\mathcal{D}}=\frac{\mathcal{D}}{\Omega_{i}}\left(\frac{2\pi}{L_{p}}\right)^{2},$
(15)

where $\rho_{i}^{2}=K_{B}T/m_{i}\Omega_{i}^{2}$ is the ion Larmor radius,
$v_{A}=B_{t}/\sqrt{4\pi\mathcal{N}m_{i}}$ the Alfvén velocity constructed with
the toroidal magnetic field and $\eta_{B}=c^{2}/4\pi\sigma$ the magnetic
diffusivity.

###  2.1 Reduced model: pure toroidal magnetic field

In this section we provide a description of a reduced version of our model,
reproducing the case of a constant toroidal magnetic field, thought as a
predictive paradigm sufficiently close to the X-point. Indeed, by setting
$B_{p}=0$, which in turns implies $\epsilon=0$ and $\gamma=1$, we obtain the
following simplified expressions for the differential operators entering the
dynamical equations of the model:

$\bm{D}_{\parallel}\to\frac{L_{p}}{L_{t}}\partial_{w}\mathbf{\hat{e}_{w}}\,,\quad\bm{D}_{\perp}\to\partial_{u}\mathbf{\hat{e}_{u}}+\partial_{v}\mathbf{\hat{e}_{v}}\,.$
(16)

For what concerns the temporal evolution of the fields involved in our
analysis, the Lagrangian advective derivative, whose general expression can be
deduced from (4), reads

$\displaystyle\frac{d}{d\tau}=\partial_{\tau}+\alpha_{1}(\partial_{u}\Phi\partial_{v}-\partial_{v}\Phi\partial_{u})\,.$
(17)

In this simplified scheme equations (12) and (13) result identical, provided
that the diffusion coefficient is set to $\mathcal{D}=\nu$. Therefore, the
model is condensed in a single equation for the vorticity, namely

$\displaystyle\partial_{\tau}D^{2}_{\perp}\Phi+\alpha_{1}\big{(}\partial_{u}\Phi\partial_{v}D_{\perp}^{2}\Phi-\partial_{v}\Phi\partial_{u}D_{\perp}^{2}\Phi\big{)}=\alpha_{2}\;D_{\perp}^{4}\Phi+\alpha_{3}\;D_{\parallel}^{2}D_{\perp}^{2}\Phi-\frac{\alpha_{3}}{\alpha_{1}}\;D^{2}_{\parallel}\Phi\,$
(18)

and the number density is recovered from the constitutive relation

$\displaystyle\bar{\mathcal{N}}=\alpha_{1}D_{\perp}^{2}\Phi\;.$ (19)

As we showed in a previous work [10], by enforcing periodic boundary
conditions and writing the electric potential in Fourier series

$\displaystyle\Phi(u,v,w)=\sum_{n,\ell,m}\varphi_{n,\ell,m}\;e^{i(nw+\ell
u+mv)}\,,$ (20)

where $(n,\ell,m)$ are integers numbers, the dynamics described by equation
(18) naturally relaxes to a 2D configuration: every non-constant contribution
along the $w$ direction, i.e. every mode with $n\neq 0$, is rapidly suppressed
due to a combined effect of inverse energy cascade and dissipation. Indeed, by
initializing to zero the $n=0$ mode and to random noise all the others $n\neq
0$ modes, it can be observed that, after a short ($\tau\sim 200$) time of
simulation, part of the energy stored in the $n\neq 0$ modes has been
transferred to the $n=0$ mode, while the rest has been dissipated by the
viscous term. Hence, in this work we study the axisymmetric turbulence in the
poloidal plane, by considering a solution $\Phi=\Phi(u,v)$ of the 2D
restriction of (18), namely

$\partial_{\tau}D^{2}_{\perp}\Phi+\alpha_{1}\big{(}\partial_{u}\Phi\partial_{v}D_{\perp}^{2}\Phi-\partial_{v}\Phi\partial_{u}D_{\perp}^{2}\Phi\big{)}=\alpha_{2}\;D_{\perp}^{4}\Phi.$
(21)

It has to be remarked that this equation is equivalent to a 2D Euler equation
for an incompressible fluid with a viscous term in the vorticity
representation [14, 9]. The equivalence becomes evident once the mapping

$\left(v_{x},v_{y}\right)\to\left(-\partial_{v}\Phi,\partial_{u}\Phi\right)$
(22)

is implemented (here $v_{x}$ and $v_{y}$ are the Cartesian components, as
referred to a coordinates system ($x,y$), of the incompressible flow
velocity). Hence, in this analogy with the fluid theory, the electric
potential $\Phi$ plays the role of a stream function, i.e. the curves on which
$\Phi$ is constant correspond to streamlines. Moreover, having recognized this
formal resemblance with the 2D Euler equation allows us to infer a number of
properties of (21) without further analyses. First, the dynamics associated to
the 2D Euler equation and, by extension, to (21), results with great
generality in a turbulent behavior. In addition to this, it can be
demonstrated that, in the ideal inviscid case, there is an infinite number of
conserved quantities associated to such equation. Indeed, when periodic
boundary conditions are imposed on a 2D region of size $A$, it turns out that
the specific kinetic energy

$K=\frac{1}{2A}\int d^{2}x\,v^{2}(x,y)$ (23)

and every other Casimir invariant constructed with a generic function of the
vorticity norm $\omega=|\mathbf{\nabla}\times\mathbf{v}|$, namely

$G_{g}=\frac{1}{2A}\int d^{2}x\,g(\omega),$ (24)

result preserved along the dynamics. Let us focus on a specific realization of
a Casimir invariant, namely the enstrophy

$\Omega=\frac{1}{2A}\int d^{2}x\,\omega^{2}(x,y).$ (25)

Now, if we imagine to expand both the velocity and the vorticity fields in
Fourier series, it is sufficient to invoke Parseval theorem to affirm that the
conserved quantities can be calculated from

$K=\frac{1}{2}\sum_{k}|v_{k}|^{2},\qquad\Omega=\frac{1}{2}\sum_{k}|\omega_{k}|^{2},$
(26)

in which $v_{k}$ and $\omega_{k}$ are the Fourier coefficients of $v$ and
$\omega$ respectively. In the equation above the sums are intended to be
performed on every allowed value of the wavenumber $k$ but, in practice, on
always handle truncated sums in which $k$ ranges between a minimum and a
maximum, i.e. $k_{min}^{2}\leq k^{2}=k_{x}^{2}+k_{y}^{2}\leq k_{max}^{2}$. For
instance, if we assume that the 2D region on which we are integrating Euler
equation is a square with side $L$, the minimum wavenumber $k_{min}$
corresponds to $\frac{2\pi}{L}$, whereas no uncontroversial upper bound can be
imposed on physical grounds in the fluid case, and the cut-off $k_{max}$ is
very often connected to a mere realistic computational time limit.
Nonetheless, it can be shown that $K$ and $\Omega$ are still conserved even
when a a description through a truncated Fourier series is performed, while
many other Casimir invariants result no longer constants. In our case, we
implement an expansion of the electric potential in a 2D Fourier series

$\Phi(u,v)=\sum_{\ell,m}\varphi_{\ell,m}\;e^{i(\ell u+mv)}$ (27)

with the Fourier coefficients satisfying the reality condition
$\varphi_{-\ell,-m}=\left(\varphi_{\ell,m}\right)^{*}$, where the $*$
indicates complex conjugation. The mapping (22) clearly shows that, in this
plasma analogue of Euler equation, the (dimensionless) kinetic energy
$\bar{K}$ and enstrophy $\bar{\Omega}$ are calculated via

$\bar{K}=\sum_{\ell,m}\left(\ell^{2}+m^{2}\right)|\varphi_{\ell,m}|^{2},\qquad\bar{\Omega}=\sum_{\ell,m}\left(\ell^{2}+m^{2}\right)^{2}|\varphi_{\ell,m}|^{2}.$
(28)

Another aspect that has to be emphasized is that a physical turbulent system
governed by the 2D Euler equation can undergo a process of condensation,
namely a transfer of energy towards the lowest wavenumber allowed by the
geometry and the boundary conditions, corresponding to an inverse cascade of
energy towards the largest spatial scale. It can be shown that the fundamental
parameter separating the regimes of condensation and non-condensation is a
characteristic wavenumber obtained as
$k^{2}_{c}=\frac{\bar{\Omega}}{\bar{K}}$. The inverse of $k_{c}$ corresponds
to the spatial size of the vortices that will dominate the dynamics once the
transient initial phase is ended and equilibrium is reached. Therefore, a
small value of $k_{c}$ indicates that the system will condensate, whereas a
large one will lead to an equilibrium characterized by vortices of many
different sizes exchanging energy without a preferred direction in the $k$
space (we remark that this occurs solely in the ideal inviscid case, given
that the usual viscous term has a form $\propto k^{4}$ in the Fourier space,
causing high-wavenumber modes to be heavily suppressed in the late stages of
the system evolution). The precise value of the threshold between the two
regimes can be empirically determined for each specific physical system under
consideration, the former being dependent on the geometry taken into account
and the choice of length and time scales used to define dimensionless
coordinates. In our case, for instance, it turns out that condensation is
achieved by choosing the initial non-zero Fourier mode such that it results
$k_{c}^{2}\simeq 10$, whereas the non-condensation scenario is reproduced for
$k_{c}^{2}\simeq 60$.

##  3 Setup of the numerical analysis

In this section we provide a description of the numerical methods implemented
to analyze the transport of particles. First, we briefly report the details of
the numerical integration for the electric potential, performed either on
equation (21), when the reduced model is considered, or on the complete set
(12), (13), (14), in the case of an X-point configuration. In both scenarios
our domain of integration is a square in the poloidal plane of side $2\pi$, on
which we impose periodic boundary conditions. The time integration is
performed via a fourth-order Runge-Kutta method, with time step $h=10^{-3}$
and a total simulation time corresponding to $\tau=10^{4}$. The non-linear
terms characterizing the equations of our model are treated by resorting to a
pseudo-spectral approach. We consider a finite number of Fourier modes
$(\ell,m)$ according to the following prescriptions: the lower bound for the
poloidal wavenumber $k$ is provided by the box size, i.e.
$k_{min}=\frac{2\pi}{L_{p}}$, whereas the fluid description of plasma adopted
in this work implies that the smallest spatial scale accessible must be
greater than the ion Larmor radius, so that it results
$k_{max}=\frac{2\pi}{\rho_{i}}$. Hence, we safely consider $\ell$ and $m$
ranging between $-13$ and $13$ with any possible combination taken into
account, excluding the case in which both indices are null, given that the
corresponding term describes a mere constant contribution. This entails a
total number of independent modes equal to $364$. For what concerns the case
with the X-point configuration, we adopted the same numerical scheme presented
in [11] for constant and radially sheared poloidal magnetic field in addition
to the toroidal one, then extended to the X-point case in [23]. This approach
is based on solving the equations for vorticity and perturbed number density
in 3D by deriving at each iteration step the electrostatic potential using a
pre-computed inverse perpendicular Laplacian. An efficient implementation of
this scheme is realized in Python using the functionalities of NumPy [24] and
SciPy [25] packages. In this work, this method is restricted to 2D and to the
equations of the reduced model by imposing the constitutive relation (19) on
initial data. Typical Tokamak parameters are assumed for both the slab and the
X-point cases: $T=100\,$eV, $B_{t}=3\,$T and $\mathcal{N}=5\times
10^{19}\,$m-3 [26], hence $\Omega_{i}\simeq 1.4\times 10^{8}\,$s-1,
$\rho_{i}\simeq 0.048\,$cm and $\nu\simeq 24.9\,$cm2s-1. The poloidal length
scale is set to $L_{p}=1$cm.

Let us now focus on the main goal of this work, namely the study of turbulent
transport from the motion of passive fluid tracers. We initialize a total
number $N=10^{4}$ of tracers randomly distributed on the square of side length
$2\pi$ previously introduced. The tracers motion is determined by the
$\mathbf{E}\times\mathbf{B}$ drift velocity, i.e. we reconstruct the
trajectory of each tracer from the numerical integration of the first-order
dynamical system

$\displaystyle\dot{u}=-\partial_{v}\Phi$ (29)
$\displaystyle\dot{v}=\partial_{u}\Phi,$ (30)

where the overdot indicates derivation with respect to the dimensionless time
$\tau$. The time integration is carried out following the same prescriptions
implemented in the case of the electric potential.

Due to the turbulent nature of the electrostatic potential $\Phi$, the
dynamical system (29)-(30) can be seen as a couple of Langevin equations.
Therefore, we consider the displacements with respect to the initial position
of each tracer at a given time, namely $\Delta u(\tau)=u(\tau)-u(0)$ and
$\Delta v(\tau)=v(\tau)-v(0)$, as a couple of random variables on which we
perform statistical analyses in order to study the underlying 2D distribution
function. Specifically, by computing the second moments
$\sigma_{uu}(\tau)=\left<\left(\Delta u(\tau)\right)^{2}\right>$ and
$\sigma_{vv}(\tau)=\left<\left(\Delta v(\tau)\right)^{2}\right>$, where
$\left<\cdot\right>$ represents an averaging over the ensemble of tracers, it
is possible to calculate the mean squared displacement (MSD) as a function of
time via the formula $\text{MSD}(\tau)=\sigma_{uu}(\tau)+\sigma_{vv}(\tau)$.
This quantity describes the different regimes of transport to which the
tracers population is subject in any phase of the simulation. Indeed, the time
dependence of the MSD can be parameterized as

$\text{MSD}(\tau)\propto\tau^{\gamma}$ (31)

and a value of the exponent $\gamma$ greater (smaller) than $1$ indicates a
supradiffusive (subdiffusive) behavior, whereas the diffusive regime
corresponds to the $\gamma=1$ case. In other words, it is possible to
calculate the turbulent transport coefficient $\mathcal{D}_{T}(\tau)$ as

$\mathcal{D}_{T}(\tau)=\frac{\text{MSD}(\tau)}{4\tau}.$ (32)

A standard diffusive phase is then connected to a linear-in-time behavior of
MSD, thus generating a constant $\mathcal{D}_{T}$. We remark that the quantity
$\mathcal{D}_{T}$ has not a direct connection with the diffusion coefficient
and the viscosity appearing in the dynamical equations of our model, namely
the parameters $\mathcal{D}$ and $\nu$ respectively. The latters play a
significant role when the the calculation of the electric potential is
performed via the integration of the system (12), (13), (14) (or merely of
equation (21) if the reduced model for pure toroidal field is considered). In
this work we are interested in studying the transport properties of the
electric potential itself due to its turbulent nature. To do so we resort to
the calculation of the MSD, which has been linked [27, 28] to the eddy
diffusivity tensor entering the Reynolds-averaged version of Navier-Stokes
equation [29]. In section 4.1 we will a posteriori verify that the turbulent
diffusion coefficient obtained from the passive tracers trajectories analysis
is indeed the same object (or a very good approximation of it) involved in one
of the closure schemes of the averaged Navier-Stokes equation, as for instance
the $K-\epsilon$ model [20].

Since the tracers tend to exit from the simulation region following the
advective $\mathbf{E}\times\mathbf{B}$ field we adopt the following
prescription: we dub as statistically significant only the portion of tracers
which remain at every time of our simulation within a square, with side length
equal to $20\pi$, centered in the origin of our 2D coordinate system. In this
manner we are able to correctly compute the MSD in unwrapped coordinates [30]
and obtain meaningful physical quantities.

Apart from the study of turbulent transport with the methods described above,
we perform normality tests on the 2D distribution function of the random
variables $\Delta u(\tau)$ and $\Delta v(\tau)$, quantifying its departure
from Gaussianity. The symmetry with respect to the central value is
investigated through the computation of a vector-valued measure of skewness,
as given in [31, 32]. This asymmetry estimator is constructed in the following
manner: given a 2D random vector $\mathbf{x}=\left(x_{1},x_{2}\right)^{T}$
with mean $\bm{\mu}$, non-singular covariance matrix $\Sigma$ with elements

$\Sigma=\begin{pmatrix}\sigma_{1,1}&\sigma_{1,2}\\\
\sigma_{1,2}&\sigma_{2,2}\end{pmatrix}$ (33)

and finite third order moments, it is possible to define the standardized
vector
$\mathbf{z}=\left(z_{1},z_{2}\right)^{T}=\Sigma^{-\frac{1}{2}}\left(\mathbf{x}-\bm{\mu}\right)$,
having by construction the null vector as mean value and the two-dimensional
identity as covariance matrix. The vector-valued skewness estimator here
adopted is a function
$\mathbf{s}\left(\mathbf{x}\right):\mathbb{R}^{2}\to\mathbb{R}^{2}$, such that

$\mathbf{s}\left(\mathbf{x}\right)=\mathbf{s}\begin{pmatrix}x_{1}\\\
x_{2}\end{pmatrix}=\begin{pmatrix}\left<z_{1}^{3}\right>+\left<z_{1}^{2}z_{2}\right>\\\
\left<z_{2}^{3}\right>+\left<z_{2}^{2}z_{1}\right>\end{pmatrix}.$ (34)

Each component of the vector $\mathbf{s}=\left(s_{1},s_{2}\right)^{T}$
accounts for an asymmetry estimator on the relative direction. Specifically, a
positive (negative) sign of a component of $\mathbf{s}$ indicates that the
mean value is greater (smaller) than the mode on the corresponding direction,
whereas a null component signals that the section of the 2D distribution
function along that specific direction results symmetrical.

Another normality test we perform is the assessment of the outliers abundance,
i.e. the weight of the distribution tales with respect to a reference Gaussian
with equal variance. This task is carried out through the computation of a
scalar-valued measure of kurtosis, i.e. a function
$\mathcal{K}:\mathbb{R}^{2}\to\mathbb{R}$, whose explicit expression reads

$\mathcal{K}=\frac{1}{\left(\sigma_{1,1}\sigma_{2,2}-\sigma_{1,2}^{2}\right)^{2}}\bigg{[}\sigma_{1,1}^{2}\left<x_{2}^{4}\right>+\sigma_{2,2}^{2}\left<x_{1}^{4}\right>+2\sigma_{1,1}\sigma_{2,2}\left<x_{1}^{2}x_{2}^{2}\right>-\\\
-4\sigma_{1,2}\left(\sigma_{1,1}\left<x_{1}x_{2}^{3}\right>+\sigma_{2,2}\left<x_{1}^{3}x_{2}\right>-\sigma_{1,2}\left<x_{1}^{2}x_{2}^{2}\right>\right)\bigg{]}.$
(35)

This quantity is equal to $8$ in the case of a 2D normal distribution and
values greater (smaller) than $8$ denote heavier (lighter) tails.
Specifically, we speak either of leptokurtic or platykurtic curves referring
to distributions with excess or deficiency of kurtosis, respectively.

##  4 Results

Here we provide an overview of the results achieved from the study of the
fluid passive tracers motion along the advective velocity given by the
$\mathbf{E}\times\mathbf{B}$ drift. The analysis has been carried out in three
distinct scenarios:

  * A1 -

the ideal inviscid case with a pure toroidal magnetic field, i.e. the electric
potential is obtained from the integration of (21) with $\alpha_{2}=0$;

  * A2 -

the viscous case, still with a toroidal magnetic field, i.e. same as in case
A1 but setting $\alpha_{2}\neq 0$;

  * B1 -

assuming a background X-point magnetic configuration with a non-null viscosity
coefficient, i.e. the electric potential is obtained from the integration of
the system (12), (13), (14).

For each of the cases described above, we present the outcome of two different
scenarios, namely the electric field on which we evolve the tracers either
undergo the process of condensation or the latter is absent from the field
evolution at any time. Our goal is therefore to investigate the role of three
distinct physical features, namely viscosity, magnetic field geometry and
condensation, as well as possible interactions between them. In each section
we describe the outcome of our analyses, providing a description of turbulent
transport via the descriptive tools introduced in the previous section, namely
the MSD and the instantaneous transport coefficient. In addition to this, we
show the result of the normality tests performed on the tracers distribution
function.

###  4.1 Case A1

In this section we start to present the results obtained from our analysis of
the turbulent transport induced by equation (21), hence in the absence of
effects due to the X-point magnetic configuration, on fluid passive tracers.
The first scenario we consider is the ideal inviscid setting, obtained by
imposing $\alpha_{2}=0$. Rather than providing a description of a realistic
physical system, the analysis performed in the ideal case constitutes a
benchmark against which a comparison with the viscous version of (21) can be
carried out, in order to precisely assess the effect of viscosity on the
transport phenomenon.

---

Figure 1: Case A1. In the left panel we report the energy spectrum, indicating
with black and red dots the condensation and non-condensation case,
respectively. In the right panels we show the advective
$\mathbf{E}\times\mathbf{B}$ field in the poloidal plane $(u,v)$ for the
condensation (top) and non-condensation case (bottom). Both the energy
spectrum and the advective fields are taken at the end of the simulation time,
namely for $\tau=9500$.

In Fig. 1 we show, for both the condensation and non-condensation cases, the
energy spectrum and the advective $\mathbf{E}\times\mathbf{B}$ field obtained
at the end of the simulation time. As previously stated, we can notice that,
in the condensation case, the energy of the system is properly stored in the
first few modes, while the remaining part of the spectrum follows a profile
$\propto k^{-2}$, as predicted by the theory [17, 18, 19]. For what concerns
the evolution of the system in the absence of condensation, instead, a quasi-
flat spectrum emerges, signaling that eddies of different sizes roughly share
the same amount of energy.

Let us begin the tracer analysis by presenting the MSD, quantifying the mean
departure from the initial position and, in a broad sense, the magnitude of
transport, in the left panel of Fig. 2. The first aspect calling for a comment
is the greater amount of transport observable in the condensation case.
Indeed, as it can be noticed from the plots, the MSD described by the black
curve is roughly $50\%$ larger than the one in red at the end of the
simulation time.

Figure 2: Case A1. The mean squared displacement and the turbulent transport
coefficient are plotted in the left and right panel, respectively. The black
curves describe the condensation case, non-condensation in red.

Moreover, it can be easily noticed that condensation causes the MSD to undergo
different phases characterized by subdiffusive and supradiffusive behaviors,
whereas in the absence of condensation the MSD grows linearly in time. This
aspect can be better emphasized by directly computing the instantaneous
transport coefficient, resulting in the plots presented in the right panel of
Fig. 2. By looking at the red curve, we can observe that the evolution of the
system obtained by initializing the electric potential integration so that
condensation is prevented is basically a quasi-pure diffusive process: indeed,
after a short transient initial phase of supradiffusive transport,
$\mathcal{D}_{T}$ stabilizes around a value $\simeq 0.030$ roughly for
$\tau=1500$ until the end of the simulation time. The same does not hold for
the condensation case: as it can be noticed from the black curve, the system
undergo a series of different transport regimes. In particular, we can observe
that there are two distinct diffusive phases: a first short one between
$\tau\simeq 2000$ and $\tau\simeq 2700$ characterized by a diffusion
coefficient $\simeq 0.08$ followed by a long subdiffusive relaxation towards a
second asymptotic diffusive regime with $\mathcal{D}_{T}\simeq 0.04$. Having
this in mind, we can conclude that, for what concerns the inviscid dynamics in
the case of a pure toroidal magnetic field, condensation acts as an enhancer
of transport, with a maximum value of the instantaneous transport coefficient
$\simeq 0.14$, in contrast to the correspondent maximum obtained in the
absence of condensation amounting to $\simeq 0.035$. In order to give an
explanation of this fact, we have to recall the microscopic interpretation of
a typical diffusive process, i.e. a random walk caused by stochastic impacts.
In the context of turbulent diffusion a single particle gets scattered every
time it interacts with a vortex and the corresponding mean free path is
connected to the typical distance between vortices. As it can be noticed from
Fig. 3, the trajectory of a tracer in the late stages of the simulation is
radically different when the condensation and non-condensation cases are
compared. The large eddies size characterizing the condensation case is
responsible for a greater transport, given that a single tracer is typically
trapped for a while inside a vortex, then gets expelled from it and start
roaming almost unperturbed until the next encounter with another vortex. The
larger amount of transport resulting on average for the whole tracer
population is therefore due to the greater distance between eddies.

Figure 3: Case A1. Portion of the $(u,v)$ poloidal plane with the trajectory
of a single tracer, for the condensation (left panel) and non-condensation
cases (right panel) starting in $\tau=9000$ (red) and ending in $\tau=10000$
(purple).

On the contrary, when the non-condensation case is addressed, we notice that a
single tracer gets scattered much more often or, in other words, the free path
between encounters with vortices is far smaller. Moreover, in the non-
condensation scenario, the interval of time in which a tracer gets trapped
inside a vortex is basically null and the obtained trajectory resembles much
more a proper random walk, hence explaining the quasi-pure diffusive character
of transport observed.

It is now important to point out that the transport coefficients calculated
above are not related to the diffusion and viscosity coefficients appearing in
the equation of the model, as previously stated in section 3. Indeed, the
object $\mathcal{D}_{T}$ here obtained from the tracers motion analysis
quantifies the rate of transport due to turbulence only, and must be therefore
linked to the turbulent diffusivity tensor emerging from the averaging of
Navier-Stokes equation (to which the basic equation (21) of our reduced model
is equivalent). The expression of the eddies diffusivity tensor in terms of
mean flow quantities is the problem addressed in the closure schemes necessary
to complete the fluid averaged theory. In this work we will consider, as
theoretical paradigm, the $K-\epsilon$ model, from which it is possible to
calculate the turbulent diffusion coefficient as

$\mathcal{D}_{T}=c_{\mu}\mathcal{N}\sqrt{K}L.$ (36)

Here $c_{\mu}$ is a constant characterizing the model and its value is
empirically determined to be $c_{\mu}\simeq 0.09$, while $L$ is the typical
size of eddies in the physical system under consideration. The latter can be
deduced from the data reported above as the square root of the MSD taken at
the beginning of the diffusive phase [16]. For instance, in the condensation
case we have seen that a short diffusive phase begins around $\tau\simeq 2000$
and the corresponding value of the MSD amounts to roughly $610$, hence
$L\simeq 24.7$. The kinetic energy $\bar{K}$ can be calculated from the
initial Fourier mode and the relative amplitudes, from (28), resulting in
$\bar{K}\simeq 1.4\cdot 10^{-3}$. We thus obtain $\mathcal{D}_{T}\simeq
0.082$, reproducing with great precision the value deduced from the data
analysis. Following the same steps we calculate the turbulent diffusion
coefficient in the non-condensation case as $\mathcal{D}_{T}\simeq 0.049$.
Here the discrepancy between theory and simulation slightly grows with respect
to the previous case, but we can nevertheless affirm that the parameter
$\mathcal{D}_{T}$ calculated from the tracers analysis satisfactorily
approximates (at least catching the right order of magnitude) the object
predicted by the closure model.

We now focus on the normality tests, performed in the spirit of highlighting
departures of the 2D probability distribution from a Gaussian profile.

Figure 4: Case A1. Skewness: condensation (left) and non-condensation (right).
The blue curve describes the skewness along the $u$ direction, orange color
for the $v$ direction.

As previously stated, we calculate estimators of the third and fourth moments
of the tracers position distribution, in order to quantify skewness and
kurtosis, whose behaviors are reported in Fig. 4 and Fig. 5, respectively. The
symmetry test conducted through the computation of skewness shows that the
deviation from a Gaussian profile is rather small. The precise outcome of this
test can randomly vary, within a narrow range around $0$, depending on the
specific choice of initialization values for the field integration. In spite
of this, we claim that, with great generality, the skewness test returns
values sufficiently close to $0$ to affirm that the 2D distribution of the
tracers position is characterized by a symmetric shape.

The kurtosis test, instead, shows a tendency which remains stable with respect
to changes in the initialization values of the field integration. Indeed, as
it can be clearly seen from Fig. 5, both in the presence and in the absence of
condensation, a significant departure from a Gaussian profile is highlighted.

Figure 5: Case A1. Kurtosis: the condensation case is described by the black
curve, non-condensation in red. The green line indicates the reference value
for a Gaussian distribution.

Specifically, in both cases we observe that the scalar-valued estimator of
kurtosis adopted in this work takes a value $\simeq 7$, pointing out a smaller
amount of outliers with respect to a Gaussian distribution having the same
variance.

###  4.2 Case A2

We now proceed to present the results obtained from our analysis of transport
induced by equation (21) but, contrarily to the previous section, here we
consider the viscous case by setting $\alpha_{2}\neq 0$. As before, we start
by presenting the spectrum and the advective $\mathbf{E}\times\mathbf{B}$
field, in Fig. 6.

---

Figure 6: Case A2. Left panel: energy spectrum for the condensation (black)
and non-condensation (red) cases. Right panels: the advective
$\mathbf{E}\times\mathbf{B}$ field in the poloidal plane $(u,v)$ for the
condensation (top) and non-condensation case (bottom). Both plots refer to the
end of the simulation time, namely for $\tau=9500$.

Figure 7: Case A2. MSD and $\mathcal{D}_{T}$ for the condensation (black) and
non-condensation (red) cases.

By comparing the plots here reported with their analog in Fig. 1 we notice
that the introduction of viscosity causes the dominance of large scale
structures in the final stages of the system evolution, even in the absence of
condensation. Hence, we expect that, from a certain point of our simulation
time, the condensation and non-condensation cases should roughly overlap, in
contrast to the ideal case presented in the previous section.

We begin to display the results obtained from the tracer analysis by providing
plots of the MSD behavior in time, as given in the left panel of Fig. 7. We
can immediately observe that, either in the presence or absence of
condensation, the total amount of transport has lowered with respect to the
inviscid case. This tendency has a more profound impact in the presence of
condensation: indeed, in this case, the introduction of viscosity is
responsible for a reduction of the MSD at the end of the simulation time of
roughly $25\%$. For what concerns the simulation in the absence of
condensation we notice instead a decrease of the MSD maximum value of about
$10\%$. The introduction of viscosity plays a major role also in modifying the
overall shape of the curves: indeed, while in the inviscid case the MSD
profiles are qualitatively different from each other, we can see that they
acquire a similar behavior when a non-null viscosity is taken into account. It
is the curve describing the non-condensation case that results particularly
modified: from the observation of the MSD profile we can affirm that a number
of phases of anomalous transport are present during the system evolution, in
contrast with the inviscid case in which the MSD grows as a straight line for
almost the entire simulation time.

Figure 8: Case A2. Components of the skewness estimator vector for the
condensation (left) and non-condensation (right) cases. The blue (orange)
curve describes the skewness along the $u$ ($v$) direction.

This fact can be further appreciated from a direct calculation of the
instantaneous transport coefficient, reported in the right panel of Fig. 7.
Let us focus on the red curve describing the behavior in time of the transport
coefficient in the absence of condensation: the qualitative discrepancy of
this curve with the analog obtained in the inviscid case is significant.
Indeed, while in the previous section we noticed a diffusive transport
throughout all stages of the system evolution, here we see that a well-defined
diffusive phase is almost not recognizable. More precisely, the instantaneous
transport coefficient has practically the same behavior as in the inviscid for
very short times, i.e. for $\tau$ between $0$ and $500$: a steep ramp
signaling a highly supradiffusive phase followed by a plateau whose value
results around $0.03$. This similarity can be easily explained given that, on
such short time scale, the effect of viscosity is obviously negligible.
However, for times greater than $\tau=500$ we see that the introduction of the
viscous term in (21) causes a completely different behavior of the
instantaneous transport coefficient. Indeed, we observe that a supradiffusive
regime characterizes the system until $\tau\simeq 4000$, followed by a
subdiffusive phase leading to an asymptotic diffusive regime, matching the
behavior of the condensation case for late times. The curve describing the
latter, in black, overlaps very closely its analog obtained in the absence of
viscosity. We can provide a simple explanation of this fact by observing that
condensation and viscosity qualitatively act in the same manner on the system
morphology: both these physical features are responsible for a dominance of
the large spatial scales on the system energy distribution, leading to a
landscape of large eddies prevailing from a certain stage of the time
evolution until its end (as pointed out by the plots in Fig. 1). This fact has
also an impact on the transport coefficient calculated in the absence of
condensation. In facts we can observe that the maximum value of the transport
coefficient is slightly greater in the viscous case with respect to inviscid
one. This observation confirms that big vortices are, in general, responsible
for an enhancement of transport. However, the effect of viscosity is also to
lower the total energy of the electric field, hence leading to a minor
transport for late stages of the evolution, as signaled by the decrease of the
maximum value of the MSD in both the condensation and the non-condensation
cases.

The introduction of viscosity has a far less significant effect on the
normality of the distribution function describing the tracers positions. As we
can observe from Fig. 8 and Fig. 9, the behavior of skewness and kurtosis
greatly overlaps their analog in the inviscid case. Indeed, as seen in the
previous section, the asymmetry of the distribution remains sufficiently close
to $0$, whereas a slight deviation from a genuine Gaussian behavior concerns
the number of outliers populating the tails of the distribution, the latter
resulting smaller than the standard values implied by normality.

Figure 9: Case A2. The scalar estimator of kurtosis $\mathcal{K}$ as a
function of time. The black (red) curve refers to the condensation (non-
condensation) case. The green line indicates the reference value for a
Gaussian distribution.

We conclude this section by remarking that, contrarily to the ideal case
displayed above, here we do not present the single tracer trajectory for late
times. This is because, as observed from Fig. 6, the introduction of viscosity
makes the condensation and non-condensation cases rather similar in the final
part of the system evolution, as already highlighted.

###  4.3 Case B: Effect of the X-point magnetic configuration

Having presented the outcome of our simulations in the case of a toroidal
magnetic field, both in the absence and in the presence of viscosity, we now
proceed to display the results obtained when the X-point magnetic geometry is
taken into account. To do so, we resort to the complete set of equations
governing the full model as given in (12), (13) and (14). We remark that here
we will only consider the viscous case in order to faithfully describe a
realistic scenario, having elucidated in the previous sections what are the
main changes occurring when viscosity is turned on starting from an ideal
inviscid setting. In Fig. 10 we display the spectrum and the advective field
obtained in this case. As it can be noticed, the condensation and non-
condensation cases are rather similar in the final stages of the integration,
due to a non-null viscous term in the equations considered. Hence, as in the
previous section, we expect that the tracer dynamics be roughly the same in
the two cases, at least in the asymptotic regime.

---

Figure 10: Case B. Left panel: energy spectrum for the condensation (black)
and non-condensation (red) cases. Right panels: the advective
$\mathbf{E}\times\mathbf{B}$ for the condensation (top) and non-condensation
case (bottom). Both the plots refer to the end of the simulation time, namely
for $\tau=9500$.

Let us begin by analyzing the MSD, depicted in the plots reported in the left
panel of Fig. 11.

Figure 11: Case B. The mean squared displacement and the turbulent transport
coefficient are plotted in the left and right panel, respectively. The black
curves describe the condensation case, non-condensation in red.

Several features of these curves are radically different from the ones
presented in the previous sections. For instance, here the total amount of
transport, deduced from the MSD maximum value, is in this case larger in the
absence of condensation (red curve). Comparing the curves in black and red, it
is easy to see that transport is lowered by an amount of roughly the $20\%$
when condensation is present (black curve). Moreover, while the MSD calculated
in the absence of condensation monotonically increases with time, as every
other MSD curve presented above, the function representing the MSD in the
condensation case presents a maximum for a time $\tau\simeq 5000$, then
slightly decreases in the subsequent phase, relaxing to a constant plateau in
the final stage of the simulation. The value of the MSD at the plateau is
about $10\%$ smaller than the maximum reached at mid-run. In addition to this,
it is important to notice that the introduction of the X-point magnetic
configuration acts as a strong limiter of particle transport. By comparing the
left panels of Fig. 7 and Fig. 11 we observe that the maximum MSD is reduced
by a factor of roughly $1/5$ and $1/3$ in the condensation and non-
condensation cases, respectively. This fact can be easily explained by
noticing that the X-point divides the square constituting our integration
domain into four (quasi) disconnected regions of triangular shape and the
presence of a magnetic shear significantly reduces turbulent transport among
them. Hence, the fraction of tracers initialized in one specific region that
migrates in a different one is very low. It is straightforward to recognize
this as the ultimate cause of the mitigated transport when the X-point
magnetic configuration is taken into account. We now turn to the transport
coefficient, whose plots are given in the right panel of Fig. 11. Let us focus
on the black curve describing the condensation case: by comparing this profile
with the analog in Fig. 7 we see that the introduction of the X-point does not
significantly alter the transport coefficient behavior at a qualitative level.
Indeed, we can still recognize two very short phases, namely a supradiffusive
followed by a subdiffusive regimes, leading to pure diffusive transport
between $\tau\simeq 1750$ and $\tau\simeq 3500$. Then, from around half of the
simulation time until the end, we observe a slow subdiffusive ramp ending in
an asymptotic diffusive regime. At a quantitative level instead, we see that
the presence of the X-point causes a reduction of the transport coefficient
value by a factor between $1/3$ and $1/4$, depending on the particular phase
considered. Let us now focus on the non-condensation case: by comparing the
red curve reported in the right panel of Fig. 11 with its analog plotted in
Fig. 7, we observe profound changes both in the overall morphology of the
curves and in the values reached by the transport coefficient. A first comment
is that it is not possible to recognize a diffusive phase in the initial
stages of the evolution but a quasi-constant plateau appears between
$\tau\simeq 6000$ and $\tau\simeq 8000$. Furthermore, the concavity of the
curve for late times does not signal the arising of a relaxation towards an
asymptotic diffusive regime. For what concerns the values of the transport
coefficient we highlight a reduction even more severe that the one noticed in
the condensation case. Indeed, for the curve obtained in the absence of
condensation we measure a transport coefficient roughly $6$ times smaller than
its analog calculated in the slab case. This being said, we remark that, as
expected from Fig. 10, the late-time features of the tracer dynamics are
roughly similar, regardless the presence or absence of condensation.

Figure 12: Case B. Skewness: condensation (left) and non-condensation (right).
The blue curve describes the skewness along the $u$ direction, orange color
for the $v$ direction. Figure 13: Case B. Kurtosis: the condensation case is
described by the black curve, non-condensation in red. The green line
indicates the reference value for a Gaussian distribution.

Let us now consider the normality tests on the distribution function. We can
observe that, also in this aspect, the X-point magnetic field produced
significant differences with respect to previous cases. First we comment on
the symmetry test conducted via the computation of skewness, whose plots are
given in Fig. 12. It is immediate to notice that the assumption of an X-point
configuration introduces a significant amount of anisotropy. Indeed, as it can
be observed in both the left and right panel, the blue and the orange curves,
describing skewness in the $u$ and $v$ directions, respectively, have two
completely different behaviors. For instance, if we focus on the right panel
we see that the skewness on the $u$ direction is practically negligible. The
same does not hold for the correspondent quantity along the $v$ direction, the
latter presenting a significant deviation from normality. As previously
stated, the precise outcome of this test has a considerable dependence from
the particular choice of initial data for the electric potential simulation,
given that the latter generates a peculiar pattern of vortices affecting the
fluid tracers dynamics and, therefore, a specific realization of anisotropies.

Figure 14: Case B. Snapshot of the tracers population at the end of the
simulation time in the poloidal plane $(u,v)$. The trajectories along which
the outliers are transported much further than the average, here dubbed
“streams”, are highlighted by red boxes.

Last, we display the curves describing the estimator of kurtosis, in Fig. 13.
We outline two main differences with the previous sections. First, the
deviation from the Gaussian value of $8$ is much greater in this case: when we
considered a toroidal magnetic field we calculated a deviation from normality
roughly equal to $1$, whereas here we notice that the departure has a
magnitude up to circa $40$. Moreover, contrarily to previously analyzed
scenarios, the difference between the kurtosis calculated from the tracer data
and the standard Gaussian value of $8$ is, in this case, positive. This
implies that the distribution function of the tracer positions is strongly
leptokurtic, with tails much more heavily populated than a Gaussian with equal
variance. Hence, we observe that the introduction of the X-point configuration
resulted in both a greater confinement, signaled by smaller values of the MSD,
and in the generation of a large population of outliers, whose positions are
much further from the initial ones with respect to the average. These two
features seem, at first glance, somewhat counterintuitive and even
contradictory, but an explanation for both phenomena can be provided. As
previously stated, the decrease of the MSD is due to the change in topology
induced by the X-point, leading to limited and disconnected regions and, on
average, on a mitigated transport. Nevertheless, it is possible to observe
that there are few specific trajectories along which the tracers migrate far
from the initial positions. These particular curves on the poloidal plane are
able to connect different regions thanks to the periodic boundary condition
characterizing our numerical integration. The landscape can be described as a
small number of narrow streams advectively transporting few portions of fluid
with great efficiency. In Fig. 14 we report a snapshot of the tracer
population, taken at the end of the simulation time. From this picture it can
be observed that a small portion of markers is transported much more than the
average along the aforementioned privileged trajectories.

##  5 Concluding remarks

In the analysis above, we performed a detailed investigation about the
influence that the two-dimensional electrostatic turbulence exerts on the test
charged particle motion, according to the $\mathbf{E}\times\mathbf{B}$ drift
velocity, as in [16, 33]. The tracers considered in this work can clearly
represent both plasma constituents as well as other ion species living in the
edge region of a Tokamak device. We discussed the slab magnetic configuration
in the inviscid and in the viscous case and also in correspondence to two
different spectral setups, i.e. with and without condensation phenomenon.
Then, we compared these results with the important case of a magnetic
configuration in which the poloidal field well mimics a real X-point in a
Tokamak equilibrium.

The slab shaped case outlined that an inverse energy cascade (the condensation
process typical of fluid dynamics) induces a more efficient transport
characterized by both diffusive and anomalous phases, in contrast to a direct
enstrophy cascade, i.e. when condensation is absent, the latter being
responsible for a quasi-pure diffusive transport of smaller magnitude. This
result can be interpreted, from a physical point of view, in the capability of
a condensation process to generate large scale eddies, able to inhibit the
stochastic field-particle scattering, while inducing a kind of particle
trapping. Nonetheless, the emergence of larger spatial scales in the
condensation case implies also an enlargement of the mean free path for
particles interacting with vortices, resulting in an overall greater spreading
of the tracer population.

The results of the analysis with an X-point magnetic configuration show a
reduction of the diffusivity by a factor of roughly $1/3$, in comparison to
the pure slab case. This is due to the suppression of the turbulence intensity
by means of the magnetic shear that the X-point configuration unavoidably
induces [1]. However, a very intriguing feature emerges when the tracer
statistical distribution function is analyzed. In fact, while the variance of
such a distribution is, as expected from a lower diffusivity (smaller than in
the slab case), the corresponding kurtosis is now much greater than its
expected Gaussian value (well-recovered in the slab case). This feature has
the important consequence that, if the tracers explore the X-point region, a
small population of them is able to perform an efficient transport, well-above
the corresponding diffusivity.

We see how the present analysis offers intriguing hints on what should be an
optimization of the turbulence features to deal with a well-behaving
(confined) plasma: we should, first of all, reduce the impact of the
condensation phenomenon, so that mainly small scale eddies are formed lowering
the magnitude of transport. Then, it is rather convenient if the particles
near enough to the X-point are prevented to reach the external region by the
creation of a transport barrier. In fact, the population of outliers, which
are able to have larger values of the MSD, should be somehow preventively
confined.

In the real picture of a Tokamak, this could be achieved either by the
introduction of a non-fluctuating radial electric field [34, 35], or by an
efficient transfer of particles from the X-point to the limiter, inhibiting
the condensation of the turbulence spectrum. In the limiter region, the
magnetic shear is much smaller than in the X-point region and the diffusivity
of a non-condensed spectrum could ensure a significant impact on the outgoing
particles.

## References

  * [1] Bruce D Scott.  Tokamak edge turbulence: background theory and computation.  Plasma Physics and Controlled Fusion, 49(7):S25, 2007.
  * [2] A. Stegmeir, A. Ross, T. Body, M. Francisquez, W. Zholobenko, D. Coster, O. Maj, P. Manz, F. Jenko, B. N. Rogers, and K. S. Kang.  Global turbulence simulations of the tokamak edge region with GRILLIX.  Physics of Plasmas, 26(5):052517, May 2019.
  * [3] D.S. Oliveira, T. Body, D. Galassi, C. Theiler, E. Laribi, P. Tamain, A. Stegmeir, M. Giacomin, W. Zholobenko, P. Ricci, H. Bufferand, J.A. Boedo, G. Ciraolo, C. Colandrea, D. Coster, H. de Oliveira, G. Fourestey, S. Gorno, F. Imbeaux, F. Jenko, V. Naulin, N. Offeddu, H. Reimerdes, E. Serre, C.K. Tsui, N. Varini, N. Vianello, M. Wiesenberger, C. Wüthrich, and the TCV Team.  Validation of edge turbulence codes against the tcv-x21 diverted l-mode reference case.  Nuclear Fusion, 62(9):096001, 2022.
  * [4] J. Wesson and D.J. Campbell.  Tokamaks.  International Series of Monographs on Physics. OUP Oxford, 2011.
  * [5] A. Hasegawa and M. Wakatani.  Plasma Edge Turbulence.  Phys. Rev. Lett., 50(9):682–686, 1983.
  * [6] A. Hasegawa and M. Wakatani.  Self-Organization of Electrostatic Turbulence in a Cylindrical Plasma.  Phys. Rev. Lett., 57(14):1581–1584, 1987.
  * [7] A. Hasegawa and K. Mima.  Strong turbulence, self-organization and plasma confinement.  Eur. Phys. J. H, 43:499–521, 2018.
  * [8] D. Biskamp and A. Zeiler.  Nonlinear Instability Mechanism in 3D Collisional Drift-Wave Turbulence.  Phys. Rev. Lett., 74(5):706–709, 1995.
  * [9] Giovanni Montani and Nakia Carlevaro.  On the 3D turbulence regime in a Tokamak plasma edge.  Physica D Nonlinear Phenomena, 451:133774, 2023.
  * [10] Nakia Carlevaro, Giovanni Montani, and Fabio Moretti.  On the effects of tokamak plasma edge symmetries on turbulence relaxation.  Symmetry, 15(9), 2023.
  * [11] F. Cianfrani and G. Montani.  On the edge turbulence in a DTT-like tokamak plasma.  arXiv e-prints, page arXiv:2405.09837, May 2024.
  * [12] B.D. Scott.  Self-Sustained Collisional Drift-Wave Turbulence in a Sheared Magnetic Field.  Phys. Rev. Lett., 65(26):3289–3292, 1990.
  * [13] B.D. Scott.  The nonlinear drift wave instability and its role in tokamak edge turbulence.  New J. Phys., 52:352, 2002.
  * [14] G. Montani, N. Carlevaro, and B. Tirozzi.  On the turbulent behavior of a magnetically confined plasma near the x-point.  Fluids, 7:157, 2022.
  * [15] Giovanni Montani and Fabio Moretti.  Symmetries of the Electromagnetic Turbulence in a Tokamak Edge.  arXiv e-prints, page arXiv:2407.01241, July 2024.
  * [16] L. Scarivaglione, F. Valentini, and S. Servidio.  Simulation campaign of the turbulent diffusion at the edge of fusion devices.  Journal of Plasma Physics, 89(1):905890108, February 2023.
  * [17] Robert H. Kraichnan.  Inertial Ranges in Two‐Dimensional Turbulence.  The Physics of Fluids, 10(7):1417–1423, 07 1967.
  * [18] Robert H. Kraichnan.  Statistical dynamics of two-dimensional flow.  Journal of Fluid Mechanics, 67(1):155–175, 1975.
  * [19] Jr. Seyler, C. E., Yehuda Salu, David Montgomery, and Georg Knorr.  Two‐dimensional turbulence in inviscid fluids or guiding center plasmas.  The Physics of Fluids, 18(7):803–813, 07 1975.
  * [20] W.P Jones and B.E Launder.  The prediction of laminarization with a two-equation model of turbulence.  International Journal of Heat and Mass Transfer, 15(2):301–314, 1972.
  * [21] H. Bufferand, G. Ciraolo, Ph. Ghendrih, Y. Marandet, J. Bucalossi, C. Colin, N. Fedorczak, D. Galassi, J. Gunn, R. Leybros, E. Serre, and P. Tamain.  Interchange turbulence model for the edge plasma in soledge2d-eirene.  Contributions to Plasma Physics, 56(6-8):555–562, 2016.
  * [22] D.D. Ryutov and V.A. Soukhanovskii.  The snowflake divertor.  Phys. Plasmas, 22(11):110901, 2015.
  * [23] F. Cianfrani and G. Montani.  Turbulent features near the X point of a DTT-like tokamak plasma.  In preparation, 2024.
  * [24] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [25] Eric Jones, Travis Oliphant, Pearu Peterson, et al.  SciPy: Open source scientific tools for Python, 2001–.
  * [26] R. Albanese, F. Crisanti, P. Martin, A. Pizzuto, G. Mazzitelli, A.A. Tuccillo, R. Ambrosino, G. Appi, A. Di Gironimo, A. Di Zenobio, A. Frattolillo, G. Granucci, P. Innocente, A. Lampasi, R. Martone, G.M. Polli, G. Ramogida, P. Rossi, S. Sandri, M. Valisa, R. Villari, and V. Vitale.  Design review for the italian divertor tokamak test facility.  Fus. Eng. Des., 146:194–197, 2019.
  * [27] G. I. Taylor.  Diffusion by continuous movements.  Proceedings of the London Mathematical Society, s2-20(1):196–212, 1922.
  * [28] G. K. Batchelor.  Diffusion in a field of homogeneous turbulence: Ii. the relative motion of particles.  Mathematical Proceedings of the Cambridge Philosophical Society, 48(2):345–362, 1952.
  * [29] Osborne Reynolds.  Iv. on the dynamical theory of incompressible viscous fluids and the determination of the criterion.  Philosophical Transactions of the Royal Society of London. (A.), 186:123–164, 1895.
  * [30] H. Qian, M.P. Sheetz, and E.L. Elson.  Single particle tracking. analysis of diffusion and flow in two-dimensional systems.  Biophysical Journal, 60(4):910–921, 1991.
  * [31] K. V. Mardia.  Measures of multivariate skewness and kurtosis with applications.  Biometrika, 57(3):519–530, 1970.
  * [32] Sreenivasa Jammalamadaka, Emanuele Taufer, and Gyorgy Terdik.  On multivariate skewness and kurtosis.  Sankhya A, 83, 08 2020.
  * [33] G. Ciraolo, Ph. Ghendrih, Y. Sarazin, C. Chandre, R. Lima, M. Vittot, and M. Pettini.  Control of test particle transport in a turbulent electrostatic model of the Scrape-Off-Layer.  Journal of Nuclear Materials, 363:550–554, June 2007.
  * [34] Kristel Crombé and Guido Van Oost.  Radial electric fields and transport barriers.  Fusion Science and Technology, 61(2T):169–179, 2012.
  * [35] L.A. Osorio-Quiroga, M. Roberto, I.L. Caldas, R.L. Viana, and Y. Elskens.  Shaping the edge radial electric field to create shearless transport barriers in tokamaks.  Fundamental Plasma Physics, 6:100023, 2023.