2.1. The numerical code

We developed an algorithm that solves (2.1)–(2.4) by using a combination of finite differences and Runge–Kutta techniques. The 2-D code solves the equation in the ($x, y$) plane assuming periodicity in the $y$ direction. We used classical, second-order, centred finite differences for the spatial differentiation and second-order Runge–Kutta for the temporal evolution of the system. The code adopts a combination of finite differences and spectral methods to solve the Poisson equation, as commonly done in systems with a periodic direction (see below).

The domain consists of a rectangular box of size $L_x \times L_y$ and a cell-centred grid with $N_x$ equidistant points in the $x$-direction and $N_y$ equidistant points in the $y$-direction. In particular, we used a box of $N_x \times N_y = 512 \times 256$ mesh points. The physical region of our simulation domain corresponds to a small rectangle, perpendicular to the mean field, at the boundary of the LCMS. It includes the outer region of the plasma column, in our domain at $x = 0$, extending to the SOL and the near-wall region of a typical fusion device (at $x = L_x$). We normalized these directions to the hybrid gyration radius $\rho _s=c_s/\omega _{ci}$, and in these normalized units, $L_x = 400$ and $L_y = 200$. The radial position corresponding to the LCMS is located at $x_0 = 100$.

We applied the following boundary conditions in the radial direction $x$. At $x = 0$ (the inner boundary) $\varOmega = 0$, and for the gradients ${\partial n}/{\partial x} = {\partial T}/{\partial x} = {\partial \varOmega }/{\partial x} = {\partial \phi }/{\partial x} = 0$. For the outer boundary, at $x=L_x$, we chose $\phi _0 = \varOmega = 0$, $n=T=1$ and again null gradients (Garcia et al. Reference Garcia, Naulin, Nielsen and Rasmussen2005). The quantity $\phi _0$ represents the poloidal average at $L_x$ of the potential.

The initial profile of the density has been chosen as a function with a steep gradient along the radial coordinate corresponding to the LCMS, which analytical form corresponds to $n(x, y) = n_i - \Delta n( 1 + \tanh { (x-x_0)/12 } )$. Here, $n_i$ is the value of the density at the inner boundary of the simulation domain, and $\Delta n$ is a parameter used to modulate the shape of the function, its ‘jump’, so that it always satisfies the outer boundary conditions. In our model, we considered equal initial profiles for both density and temperature. Therefore, we defined the parameter $A$ as $A = n_i = T_i$ and $\Delta A = \Delta n = \Delta T$.

Because of periodicity, $\phi$ and $\varOmega$ can be expressed in terms of the Fourier basis as $\phi (x,y)=\sum _{n=1}^{N_y} \hat {\phi }(x)\exp ({{\rm i} n 2{\rm \pi} /L_y y})=\sum _{n=1}^{N_y} \hat {\phi }(x)\exp ({{\rm i} n k_y y})$ and $\varOmega (x,y)=\sum _{n=1}^{N_y} \hat {\varOmega }(x)\exp ({{\rm i} n k_y y})$ (here, we suppressed the time dependency, for simplicity). Then, by employing a centred finite difference scheme for the second $x$-derivative of $\phi$, from (2.4), one gets the following equation at the $j$th grid point in the $x$-direction:

(2.8)\begin{equation} \frac{\hat{\phi}(x_{j+1})+\hat{\phi}(x_{j-1})-2\hat{\phi}(x_j)}{\Delta x^2}-n^2k_y^2\hat{\phi}(x_j)=\hat{\varOmega}(x_j). \end{equation}

By considering the $x$-boundary conditions discussed above, the previous equation can be written in the form of a linear tridiagonal system for the unknown coefficients $\hat {\phi }(x_j), j=1,\ldots,N_x$. This linear system can be solved through standard linear algebra routines (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1996). Finally, the inverse Fourier transform in $y$ returns the values of $\phi (x_j,y_l)$ at all grid points in the numerical 2-D box.

Regarding the shape of the magnetic field, it was important to set the parameters $\epsilon$ and $\zeta$. Following the literature, their values are usually small, and we set $\epsilon = 0.25$. Since we want to study a toroidal system characterized by a radial gradient of the main magnetic field, we have always assumed a non-null curvature parameter $\zeta$. With this setting, any small perturbation suddenly generates instabilities and turbulence. These instabilities, indeed, are suppressed for $\zeta = 0$, because there is no magnetic gradient that causes the vertical polarization, and hence a radial electric drift. Without this gradient there is no triggering of interchange instability mechanism. This scenario was confirmed by performing a test simulation (not shown here) in which we set $\zeta = 0$.

The collisional coefficients $\nu _n = \nu _T = \nu _\varOmega$ are set to a very low number (${\sim }10^{-2}$), in order to resolve well the structures and to dissipate artificial small-scale fluctuations. In the expression of the linear damping, (2.5), two other parameters have to be defined. In particular, we used $\lambda = 1.6 \times 10^{-4}$ and $\delta = 1$. The parameters $\zeta$ and $\nu$ are varied for each run, as described in the next section.

The large-scale plasma profile cannot be described self-consistently in the framework of our simplified (reduced) Braginskii model. For $x< x_0$, the closed field line region, we therefore maintain constant the mean large-scale profiles, assuming that these would vary on time scales much larger than the typical turbulence dynamical times. In our numerical simulations, a steady-state regime of transport is ensured by maintaining constant time density and temperature poloidal profiles. In this regard, we developed a new driving technique based on an averaging method: we forced the initial poloidal average of the density and temperature fields to remain constant within the LCMS region. The poloidal average of any field is defined as the spatial average over the periodic direction, namely $\langle f \rangle (x) = ({1}/{L_y})\int f(x,y') \,{{\rm d} y}'$. In our model, at each time step, we subtracted the poloidal averages of the plasma density and temperature, replacing them with the initial profiles. This procedure was applied only in the first part of the domain, for $x < x_0$. The method is very efficient and allows us to achieve a stationary configuration of emission/absorption. This novel technique, different from the previous modelling of edge turbulence, is based on maintaining a constant poloidal inner profile. We named this novel forcing procedure the average freezing driving (AFD). In figure 2 we report the behaviour of the poloidal average of the density as a function of time, during a turbulent steady state. A similar approach, applied to other geometries, has been found to give very stable, statistically steady-state simulations of SOL turbulence (Servidio et al. Reference Servidio, Primavera, Carbone, Noullez and Rypdal2008).

Figure 2. Poloidally averaged density at three different times. The initial mean profile remains fixed in the inner region of the simulation domain (RUN$_{{\rm II}}$).

3.1. The Eulerian approach

In order to characterize the dynamics of SOL turbulent plasmas, it is important to define some relevant global quantities. The radial profile of any field, denoted hereafter by a ‘zero’ subscript, is defined as its spatial average over the periodic direction $y$ (as done for the driving technique.) We can define the density and the poloidal flux profiles as $n_0(x,t) = ({1}/{L_y}) \int ^{L_y}_0 n(x,y,t)\,{{\rm d} y}$ and $\nu _0(x,t) = ({1}/{L_y} )\int ^{L_y}_0 ({\partial \phi }/{\partial x})\,{{\rm d} y}$ (Garcia Reference Garcia2001). The deviation from the mean profile of any field is identified as the spatial fluctuation and it is indicated by a tilde, i.e. $\tilde {n}(x,y,t) = n(x,y,t)-n_0(x,t)$. Note that, because of the driving at $x< x_0$, both the poloidal averages of particle density and temperature are constant in time in this inner region. This is an advantage of the AFD, which guarantees stable profiles in time.

To characterize the turbulent state of the system in the outer layer, we studied the evolution in time of the kinetic energy, considering its two components

(3.1 a,b)\begin{equation} K(t) = \int \tfrac{1}{2}(\nabla \tilde{\phi} )^2 \,{\rm d}\boldsymbol{x}, \quad U(t) = \int \tfrac{1}{2}\nu_0^2 \,{\rm d}\boldsymbol{x}, \end{equation}

that is, the kinetic energy in the fluctuating motions and the sheared poloidal flows, respectively. Following Garcia et al. (Reference Garcia, Bian, Paulsen, Benkadda and Rypdal2003), the evolution of these integrals is given by some straightforward manipulations of (2.3) and (3.1a,b), leading to

(3.2)\begin{gather} \frac{{\rm d}K}{{\rm d}t} = \zeta \int nT\widetilde{u_x} \,{\rm d}{\boldsymbol{x}} - \int \nu_0 \frac{\partial}{\partial x}( \widetilde{u_x} \widetilde{u_y})_0 \, {\rm d}{\boldsymbol{x}} - \int \tilde{\phi} \varLambda_\varOmega \,{\rm d}{\boldsymbol{x}}, \end{gather}

(3.3)\begin{gather}\frac{{\rm d}U}{{\rm d}t} = \int \nu_0 \frac{\partial}{\partial x}( \widetilde{u_x} \widetilde{u_y})_0 \,{\rm d}{\boldsymbol{x}} - \int \phi_0 \varLambda_\varOmega \,{\rm d}{\boldsymbol{x}}, \end{gather}

where $\tilde {u}_x$ and $\tilde {u}_y$ are the radial and poloidal components of the electric drift $\boldsymbol {u}_E$, respectively. The first term of (3.2) corresponds to the fluctuating motion change due to the radially advective transport of thermal energy. The second (which appears also in (3.3)) represents the conservative transfer of kinetic energy from fluctuating motion to the sheared flow (i.e. kinetic energy may be transferred between the convective cells and sheared poloidal flows by a radial inhomogeneity in the off-diagonal components of the averaged Reynolds stress tensor). The last term, present in both equations, $\varLambda _\varOmega = \nu _\varOmega \nabla ^2 \varOmega - \sigma _{\varOmega }(x) \varOmega$, represents the diffusive term and the linear damping due to collisions and particle transport. It corresponds to the right-hand side of (2.3). To summarize, the energy transfer rates from thermal energy to the fluctuating motion, and from fluctuating to the sheared flows, can be defined, respectively, as

(3.4 a,b)\begin{equation} F_p(t) = \zeta \int n T \widetilde{u_x} \,{\rm d}\boldsymbol{x}, \quad F_v(t) = \int \frac{\partial \nu_0}{\partial x}( \widetilde{u_x} \widetilde{u_y})_0 \,{\rm d}\boldsymbol{x}. \end{equation}

These fluxes are crucial for the characterization of turbulent transport, as we shall see in the following.

The energies in (3.1a,b) and the fluxes in (3.4a,b) are highly representative of the system dynamics. Hereafter, we show the results for the main RUN$_{{\rm II}}$. The relation between all the above global quantities is better represented in figure 3, where we focused on a single ‘burst event’. In correspondence with blob ejection, there is a rapid growth of the energy transfer rate $F_p$ since the available source of free energy, contained in the mean pressure gradient, is converted into $K$. Then, after a short time, there is a rapid decrease of the latter in favour of the sheared poloidal flow. A peak in the quantity $F_v$ is at the same time observed, although there is a small time shift with the burst in $U$. Overall, these peaks and correspondences are only qualitative, since they emerge from high-level turbulent fluctuations. Both the trend of the kinetic energy and that of the energy transfer rate suggest that the SOL is characterized by the alternation of blob ejection moments and quiet periods.

Figure 3. Kinetic energy and energy transfer rate in a single burst event. Panels represent (a) kinetic energy in the fluctuating motion $K$, (b) the sheared poloidal flow $U$, (c) the energy transfer rate from thermal energy to the fluctuating motion $F_p$ and (d) the transfer rate from fluctuating to sheared flows $F_v$.

Generally, the evolution in time of the kinetic energy of the fluctuating motions is characterized by irregular oscillations with pronounced bursts. This is a manifestation of intermittent, bursty turbulence, typically observed at the edge of many tokamak devices. During such burst events, the formation of blobs takes place. Every time that this fluctuation level is sufficiently large, there is a rapid growth of the energy in the sheared poloidal flows $U$, followed by a suppressed level of $K$. Such a dynamics is the consequence of the exchange of kinetic energy from the fluctuating motion to the sheared flows, as can be evinced from (3.2)–(3.3).

The difference between an emission moment (right) and a quiet period (left) is reported in figure 4 for density, temperature, vorticity and potential, from top to bottom, respectively. In the quiet regime, the plasma is quite homogeneous and isothermal and almost no fluctuations can be observed in the SOL region. There are plumes in the particle density and temperature fields. At the time of a burst event (e–h), it is evident that prominent structures develop in the SOL. During such a burst event, an accumulation of particle density and heat is clearly evident. This ejection is localized in the region of the LCMS, in the SOL, at $x>x_{0}$.

Figure 4. Comparison of the spatial structure of particle density, temperature, vorticity and potential during a quiet period, (a–d), and during a burst event, (e–h).

During a burst event, it is possible to observe the formation of a blob by following the evolution over time of particle density, as can be seen from figure 5. These plots show the rapid development of a structure from the inner region to the outer part of the domain, i.e. towards the open wall. During the emission, the main vortex can also merge with and eat small neighbouring fluctuations.

Figure 5. Spatial evolution of particle density over time, in sequence from (a–f), showing a blob emission from the edge of the plasma column and its propagation towards the outer boundary.

Another useful Eulerian analysis that can help to understand the propagation in space and time of turbulent structures is the analysis of macroscopic variables such as the density at selected regions of the numerical domain (to mimic probe measurements in the SOL). We measured the density and its probability distribution function (p.d.f.) at different points (probes) of the simulation domain. These are 7 probes equally spaced, numbered going from left to right of the domain (so the first probe is situated in the plasma edge region, while the last one is near the wall). These equally spaced probes start at $x_1= 64$, following at $x_j = j \delta$, with $j = 2,\ldots, 7$ with separation $\delta =64$.

Panel $(a)$ of figure 6 represents the change in time of the density at different points of the domain. The signal is stronger at the first probes, because of the higher concentration of density in the inner region. Progressively, the amplitude decreases near the wall, at probe $P_7$. We can see that some events are correlated. These effects are due to the turbulent transport and reveal that blobs structures are moving radially toward the wall. The particle density at the probes within the edge region, $P_1$ and $P_2$, shows chaotic oscillations which appear symmetric concerning the mean value, while all other probe signals reveal predominantly positive fluctuations (higher skewness), as highlighted also from the p.d.f. analysis.

Figure 6. (a) Temporal evolution of the particle density measured by the different probes $P_j$. Note that the signal has been shifted in the $y$-axis to compare, on the same plot, all the measurements; (b) The p.d.f. of the particle density signals, compared for three different probes.

In figure 6(b), we report the p.d.f. of the particle density, at three different probes. The presence of intermittent fluctuations causes a departure of the p.d.f. from the normal distribution, manifesting very high exponential tails (Antar et al. Reference Antar, Krasheninnikov, Devynck, Doerner, Hollmann, Boedo, Luckhardt and Conn2001b; Garcia et al. Reference Garcia, Bian, Paulsen, Benkadda and Rypdal2003; Sattin et al. Reference Sattin, Vianello, Valisa, Antoni and Serianni2005; Killer et al. Reference Killer, Narbutt and Grulke2021; Tatali et al. Reference Tatali, Serre, Tamain, Galassi, Ghendrih, Nespoli, Bufferand, Cartier-Michaud and Ciraolo2021). This is typical of intermittent, extraordinary events. As we can see, at the first probe which is located within the plasma edge, the ${\rm p.d.f.}(n)$ has a nearly symmetric, super-Gaussian shape. At larger distances from the LCMS, due to the intermittent nature of turbulence, it gradually develops higher tails, skewed in the outward direction. This skewness and kurtosis of the distribution is commonly observed in SOL of fusion devices and reflects the high probability of extreme events, as largely investigated in the literature (Antar et al. Reference Antar, Counsell, Yu, Labombard and Devynck2003; van Milligen et al. Reference van Milligen, Sánchez, Carreras, Lynch, LaBombard, Pedrosa, Hidalgo, Gonçalves and Balbín2005; Xu et al. Reference Xu, Jachmich and Weynants2005; Kamataki et al. Reference Kamataki, Itoh, Inagaki, Arakawa, Nagashima, Yamada, Yagi, Fujisawa and Itoh2010).

At this point, we changed the simulation settings to explore the parameter space. By changing $\zeta$ and $A$, we analysed different turbulent regimes and compared all the simulations by evaluating their power spectra and the p.d.f.s of the particle density. For the power spectra, we Fourier transformed the density field $n(x,y,t)$ along the periodic direction, computed the energy spectrum $E(x, k_y, t)=|\tilde {n}(x, k_y, t)|^2$ ($\tilde {n}$ being the Fourier coefficient) and finally obtained the spectrum from an ensemble average both in time and space (in the region $x>x_{0}$).

In RUN$_{\rm I}$, RUN$_{{\rm II}}$ and RUN$_{{\rm III}}$ we have increased progressively the value of $\zeta$. Since this parameter governs the variation along the radial coordinate of the magnetic field, we have effectively increased the gradient of $\boldsymbol {B}$, as can be evinced from figure 1. In figure 7(a) we compare the averaged spectra of the plasma density $n$, as a function of the poloidal mode $k_y$, for this campaign of runs. The turbulent spectrum progressively decreases from RUN$_{\rm I}$ to RUN$_{{\rm III}}$. The spectra all manifest a power-law scaling, consistent with a fluid-like turbulent cascade (Kolmogorov Reference Kolmogorov1941a, Reference Kolmogorovb). We report also the Kolmogorov scaling law as a reference. Even if the scaling is consistent from one run to another, the level of turbulence decreases. This trend reveals that an increase in the value of $\zeta$ leads to more stable configurations, while the system characterized by a small value of $\zeta$ is affected by a stronger turbulent transport. The suppression of turbulent transport is also evident in the p.d.f.s in figure 8(a): the low-$\zeta$ cases are more intermittent.

Figure 7. (a) Comparison between the average energy spectrum of RUN$_{\rm I}$, RUN$_{{\rm II}}$, and RUN$_{{\rm III}}$. The energy spectrum decreases as the $\zeta$ value increases, revealing a more stable configuration in the RUN$_{{\rm III}}$, which is characterized by a higher value of $\zeta$; (b) Comparison between the average energy spectrum of RUN$_{{\rm II}}$, RUN$_{{\rm IV}}$, and RUN$_{{\rm V}}$. In this case, the energy spectrum decreases in the presence of a smaller pressure gradient, reflecting the hypothesis that a higher gradient of density injects more turbulence into the system.

Figure 8. (a) Comparison of the p.d.f. of particle density between RUN$_{\rm I}$, RUN$_{{\rm II}}$, and RUN$_{{\rm III}}$. More pronounced tails are present in the RUN$_{\rm I}$ and RUN$_{{\rm II}}$ rather than in the RUN$_{{\rm III}}$, indicating more intermittency in the former cases; (b) Comparison of the p.d.f. of particle density between RUN$_{{\rm II}}$, RUN$_{{\rm IV}}$, and RUN$_{{\rm V}}$. The latter p.d.f.s are quite similar, revealing that density and temperature gradients do not influence too much the nature of the intermittent transport.

RUN$_{{\rm II}}$, RUN$_{{\rm IV}}$ and RUN$_{{\rm V}}$ are characterized by different values of the initial density and temperature jump $A$, as described in the table. Changing the value of these quantities at the inner boundary of the domain corresponds to a variation of the gradient of pressure that, as a consequence of the inhomogeneity of the magnetic field $\boldsymbol {B}$, drives instabilities. From the spectra in figure 7(b) it is clear that higher density and temperature gradients (higher $\Delta A$) cause more turbulent activity. The intermittent behaviour does not seem to be affected so much by the change of $A$, as can be seen from the p.d.f.s in figure 8(b). This picture confirms that the available source of free energy, driving the convective motions, is contained in the mean pressure gradient. An increase in the free energy leads to an increase in the fluctuating motion and, consequently, to a higher (radial) turbulent transport.

3.2. The Lagrangian approach

In the second part of our analysis, in order to better understand the diffusion of fluid elements in the SOL region and characterize the transport dynamics, a new technique is investigated, based on the Lagrangian approach. We followed in time the motion of fluid tracers in the drift approximation, integrating the trajectory of the plasma element in the numerical domain. In practice, we solved particle integration in our code, following the journey of plasma elements from the central part of the LCMS toward the colder, empty regions at the wall. We solved the fluid equations of motion in which we considered only the major contributions of the cross-field fluid drift – the electric and the diamagnetic drifts.

The equation of motion for the fluid tracers has been derived by using the same Braginskii approximation of our model: we neglected the ion dynamics, concentrating only on the motion of the electrons, and neglected the inertial term in the cross-field drift. The drift velocity then is simply

(3.5)\begin{equation} \boldsymbol{u_\perp} = \boldsymbol{u}_E + \boldsymbol{u}_d = \frac{1}{B} \boldsymbol{\hat{z}} \times \boldsymbol{\nabla} \phi - \frac{1}{nB} \hat{\boldsymbol{z}} \times \boldsymbol{\nabla} p. \end{equation}

Plasma fluid elements, therefore, move with the above velocity, at each time, by obeying to $\dot {\boldsymbol {x}} = \boldsymbol {u_\perp }$. We integrated the trajectories of fluid particles that follow passively the flow. In this regard, we upgraded the numerical code with this Lagrangian tracing technique where the fields have been continuously supplied by the time-varying simulation, and by obtaining the velocity (3.5) at each particle position via a bilinear interpolation algorithm (second order in space).

At the beginning (or at the earlier stages) of each simulation, we placed all the Lagrangian tracers at the LCMS ($x_0$), but with different (random) poloidal positions $y$. We considered the following boundary condition for the tracers: if an element reaches the radial boundary, we removed it from the domain (this is because we are interested in particles within the SOL region); in the poloidal direction, instead, we considered periodic boundaries. The parameters of the Lagrangian integration are summarized in table 1, and their initial motion is represented in figure 9. We show the evolution in time of the trajectories for just 10 fluid elements, to avoid overcrowding. In this figure, one can see how, from the radial position corresponding to the LCMS, elements start their erratic trip in the simulation domain under the effect of turbulence. Some of them are absorbed by the left part of the wall, namely inside the device. Some of them are lost from the outer SOL (right side).

Figure 9. Examples of trajectories of the fluid elements, for RUN$_{{\rm II}}$, starting from the initial position at $x_0$ (here highlighted by a black vertical line). Only early times of such trajectories, for a limited sample, are reported.

To perform a statistical analysis of the turbulent diffusion mechanism, we considered a large sample of fluid elements, $N = 10\,000$, and calculated their displacement from the initial position $\Delta x(t)= x(t)-x_0$, along the $x$ axis. Note that, with this set-up, the maximum value of the displacement is $\Delta x_{\max } = 300$. We verified the convergence of the statistics presented below by increasing systematically the number of particles, as reported later in this section. We noticed that plasma elements in the SOL region experience an irregular random walk.

To make contact with classical diffusion theories, we evaluated the mean square displacement (MSD) of the radial positions $\langle \Delta x^2 \rangle$, where $\langle \bullet \rangle$ represents an ensemble over all the particles. In figure 10, we reported the MSD as a function of time. The behaviour of the MSD is quite interesting, given the fact that we are in a small, finite-size system. In the initial stage of their journey, the particles experience a behaviour consistent with $\tau ^2$, typical of free-streaming motion (Taylor Reference Taylor1921; Taylor & McNamara Reference Taylor and McNamara1971; Huang et al. Reference Huang, Chavez, Taute, Lukić, Jeney, Raizen and Florin2011; Pusey Reference Pusey2011). Here, $\tau$ is simply the difference between $t$ and the injection time $T_{0p}$ (see table 1). After the above initial transient, the curve suggests a diffusive behaviour, $\langle \Delta x^2\rangle \sim \tau$. As in typical plasma turbulence (Pecora et al. Reference Pecora, Servidio, Greco, Matthaeus, Burgess, Haynes, Carbone and Veltri2018), this regime starts when the displacement becomes of the order of the correlation length of turbulence $\lambda _c$. The latter corresponds to the typical size of the largest energy-containing eddies, namely the size of the blobs. For the above simulations, we estimated this length to be, qualitatively, $\lambda _c \sim 20$. The regime of the diffusive transient starts therefore at a characteristic correlation time $\tau _c$, highlighted in the figure with a vertical line. At later times, contrary to unbounded problems of diffusion, the behaviour is influenced by the system size, and the large-scale diffusive regime is lost. At these large times, most of the particles abandoned the domain, either because they are absorbed by the inner boundary or because they are lost at the outer wall.

Figure 10. Temporal evolution of the mean square displacement of the fluid elements, in a double logarithmic scale, for RUN$_{{\rm II}}$. The horizontal (dashed) line represents the square of the correlation length $\lambda _c$, namely the average size of the blobs. The characteristic decorrelation time is reported with a vertical dashed line. The thick (black) line represents the maximum allowed MSD, $\Delta x_{\max }^2$.

To further confirm this diffusive behaviour, we computed the diffusion coefficient, but first, we checked the convergence of our results by increasing systematically the statistics. We divided the simulation times into subdomains. In particular, we divided the particle integration into eight, smaller, subsequent shots, each of which lasts $5\times 10^{4}$ computational times. At each particle restart, we set all the fluid elements at the initial position ($x_0 = 100$) and we studied how the diffusion changes over time. In figure 11 we show the temporal convergence of the MSD: ‘$\langle \Delta x^2\rangle _{[1]}$’ refers to the MSD of the first part of the simulation, the second, $\langle \Delta x^2\rangle _{[1-2]}$, refers to the average of the first and second part of RUN-II, and so on. We continued to enlarge the statistics until we comprehended all the eight time subdomains. We can observe that the curves progressively increase, until they collapse, suggesting a clear convergence of the statistics (note that the convergence is achieved already with an average on just six subsets).

Figure 11. Convergence of the statistical analysis, for RUN$_{{\rm II}}$. Other runs have similar behaviours.

At this point, the diffusion coefficient can be computed as $D \sim \langle \Delta x^2\rangle /(2 \tau )$, for each average shown in figure 11. Similar behaviour has been found for different runs (not shown here). We compared the diffusion coefficient, averaged over all the datasets, for the runs where we had a strong enhancement of turbulence, namely by comparing RUN$_{{\rm II}}$, RUN$_{{\rm IV}}$ and RUN$_{{\rm V}}$. In figure 12 one can see that diffusion is effectively higher for the run with stronger turbulence, namely the case where the density and temperature driving is more intense ($A=3.3$). This preliminary analysis of the diffusion process further confirms that turbulence enhances particle diffusion, even in such a simplified Braginskii model of the SOL.

Figure 12. Comparison of the diffusion coefficients (in a double logarithmic scale) between RUN$_{{\rm II}}$, RUN$_{{\rm IV}}$ and RUN$_{\rm V}$.