# Role of avalanche transport in competing drift wave and interchange turbulence

Philippe Ghendrih, Guilhem Dif-Pradalier, Olivier Panico, Yanick Sarazin, Hugo Bufferand, Guido Ciraolo, Peter Donnel, Nicolas Fedorczak, Xavier Garbet, Virginie Grandgirard, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Philippe Ghendrih, Guilhem Dif-Pradalier, Olivier Panico, Yanick Sarazin, Hugo Bufferand, et al. Role of avalanche transport in competing drift wave and interchange turbulence. Journal of Physics: Conference Series, 2022, 2397, pp.012018. 10.1088/1742-6596/2397/1/012018. [https://doi.org/10.1088/1742-6596/2397/1/012018](https://doi.org/10.1088/1742-6596/2397/1/012018). [hal-04498387](https://HAL.inistre.fr/hal-04498387)

HAL Id: hal-04498387

https://hal.science/hal-04498387

Submitted on 5 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Role of avalanche transport in competing drift wave and interchange turbulence

Philippe GHENDRIH<fnref n="1" />, Guilhem DIF-PRADALIER<fnref n="1" />,
Olivier PANICO<fnref n="1" />,<fnref n="2" />, Yanick SARAZIN<fnref n="1" />, Hugo Bufferand<fnref n="1" />,
Guido CIRAOLO<fnref n="1" />, Peter DONNEL<fnref n="1" />, Nicolas FEDORCZAK<fnref n="1" />,
Xavier GARBET<fnref n="1" />, Virginie GRANDGIRARD<fnref n="1" />,
Pascale HENNEQUIN<fnref n="2" />, Eric SERRE<fnref n="3" />, Patrick TAMAIN<fnref n="1" />

E-mail: philippe.ghendrih@cea.fr

**Abstract.** We complete the 2D 2-fields turbulence model previously used with an interchange-like instability by slightly modifying the parallel loss terms to drive drift wave instabilities. We show that the instability driven by temperature fluctuations of the sheath losses is identical to that of the drift wave turbulence. The linear analysis is performed and used to select control parameters that yield identical maximum growth rates for the interchange alone and drift wave alone instability. Combining the two instabilities doubles the maximum growth rate. The non-linear simulations are used to analyse the SOL width. The simulations allow one to identify a low field side SOL region where interchange and drift wave are unstable and a high field side SOL region where only the drift wave is unstable. The SOL profiles appear exponential in the region close to the source but depart from a simple exponential fall-off in the far SOL. The low field side SOL width is found to be larger in the interchange alone case, slightly smaller when both instabilities are present and finally narrower when only the drift waves. For the high field side SOL, without interchange, the drift wave SOL width is observed to be identical to that on the low field side and larger than that when both instabilities are present. The Sherwood dimensionless parameter, ratio of convective particle flux divided by the diffusive particle flux, is used to compare the efficiency of turbulent transport. The profiles of the Sherwood parameter for time and flux surface averaged transport indicate that turbulent transport is dominant close to the separatrix but is less effective towards the far SOL. The Sherwood parameter evolution, determined with the flux-surface averaged transport, indicates that outward avalanche transport with corrugations governs the case with interchange only. When combining the two instabilities, outward avalanche transport is less pronounced and inward avalanche transport is observed, reducing the overall turbulent transport efficiency. The avalanche transport with drift waves only compared to interchange only is found to be inhibited.

## 1. Introduction

The problem of turbulent transport and its control remains a key issue for the performance of fusion devices. In the core it is detrimental and reduces the confinement [1] while in the SOL it is beneficial by increasing the width of this boundary layer, thereby spreading the heat load on the plasma facing components [2]. Given this importance, a dedicated effort is organised to address this problem. On the one hand sophisticated gyrokinetic codes [3, 4] have been developed to address core transport and a comparable effort is dedicated to the edge and Scrape-Off Layer turbulence [5, 6], mostly using fluid models. Besides such first principle models, reduced models

<figure><img src="image_1.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France.</footnote>

<footnote n="2">² LPP, CNRS, Ecole polytechnique, 91128 Palaiseau, France.</footnote>

<footnote n="3">³ Aix Marseille Univ, CNRS, Centrale Marseille, M2P2, Marseille, France.</footnote>
# Role of avalanche transport in competing drift wave and interchange turbulence

Philippe GHENDRIH<fnref n="1" />, Guilhem DIF-PRADALIER<fnref n="1" />,
Olivier PANICO<fnref n="1" />,<fnref n="2" />, Yanick SARAZIN<fnref n="1" />, Hugo Bufferand<fnref n="1" />,
Guido CIRAOLO<fnref n="1" />, Peter DONNEL<fnref n="1" />, Nicolas FEDORCZAK<fnref n="1" />,
Xavier GARBET<fnref n="1" />, Virginie GRANDGIRARD<fnref n="1" />,
Pascale HENNEQUIN<fnref n="2" />, Eric SERRE<fnref n="3" />, Patrick TAMAIN<fnref n="1" />

E-mail: philippe.ghendrih@cea.fr

**Abstract.** We complete the 2D 2-fields turbulence model previously used with an interchange-like instability by slightly modifying the parallel loss terms to drive drift wave instabilities. We show that the instability driven by temperature fluctuations of the sheath losses is identical to that of the drift wave turbulence. The linear analysis is performed and used to select control parameters that yield identical maximum growth rates for the interchange alone and drift wave alone instability. Combining the two instabilities doubles the maximum growth rate. The non-linear simulations are used to analyse the SOL width. The simulations allow one to identify a low field side SOL region where interchange and drift wave are unstable and a high field side SOL region where only the drift wave is unstable. The SOL profiles appear exponential in the region close to the source but depart from a simple exponential fall-off in the far SOL. The low field side SOL width is found to be larger in the interchange alone case, slightly smaller when both instabilities are present and finally narrower when only the drift waves. For the high field side SOL, without interchange, the drift wave SOL width is observed to be identical to that on the low field side and larger than that when both instabilities are present. The Sherwood dimensionless parameter, ratio of convective particle flux divided by the diffusive particle flux, is used to compare the efficiency of turbulent transport. The profiles of the Sherwood parameter for time and flux surface averaged transport indicate that turbulent transport is dominant close to the separatrix but is less effective towards the far SOL. The Sherwood parameter evolution, determined with the flux-surface averaged transport, indicates that outward avalanche transport with corrugations governs the case with interchange only. When combining the two instabilities, outward avalanche transport is less pronounced and inward avalanche transport is observed, reducing the overall turbulent transport efficiency. The avalanche transport with drift waves only compared to interchange only is found to be inhibited.

## 1. Introduction

The problem of turbulent transport and its control remains a key issue for the performance of fusion devices. In the core it is detrimental and reduces the confinement [1] while in the SOL it is beneficial by increasing the width of this boundary layer, thereby spreading the heat load on the plasma facing components [2]. Given this importance, a dedicated effort is organised to address this problem. On the one hand sophisticated gyrokinetic codes [3, 4] have been developed to address core transport and a comparable effort is dedicated to the edge and Scrape-Off Layer turbulence [5, 6], mostly using fluid models. Besides such first principle models, reduced models

<figure><img src="image_3.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France.</footnote>

<footnote n="2">² LPP, CNRS, Ecole polytechnique, 91128 Palaiseau, France.</footnote>

<footnote n="3">³ Aix Marseille Univ, CNRS, Centrale Marseille, M2P2, Marseille, France.</footnote>
[7, 8, 9] are also contemplated to have fast access to simulations of experiments and forecast experimental scenarios. We propose here an intermediate model providing first principle analysis of turbulent transport but reduced to the simplest 2D, 2 fields model [10, 11] that make such a model very attractive to generate large data bases with scans of control parameters [12]. Such a model is flux driven and characterised by avalanche transport [10, 13], also reported as "blob" transport [14]. Similar transport behaviour has been observed for heat in the gyrokinetic framework [15]. This model has also allowed us to recover staircase confinement [16], first analysed in the gyrokinetic framework [17]. Our claim in such an approach is that physics problems such as turbulent transport must be robust so that a strongly simplified fluid code should capture qualitatively the behaviour of the system, provided the problem at hand is not specifically governed by kinetic properties, such as the problem of the sheath physics [18, 19]. The model that is used is known as the TOKAM2D model [20]. It was developed to investigate interchange like turbulence in the SOL [21, 22]. Recently it has been extended to analyse the interface transport barrier between the edge plasma with closed field lines and the SOL with sheath boundary conditions [23], considering the interchange instability mechanism and zonal flow generation [24]. The model evolves a transported field, typically the density at constant thermal energy, set out of thermodynamic equilibrium either by a prescribed gradient -as used to investigate linear stability or gradient driven turbulence- or by a radially localised source such that the mean flux is imposed and the field gradient an output of turbulent transport. The latter regime has been named flux driven. We show in this paper that the transported field can be the thermal energy, assuming the density to be constant. The charge balance equation for a quasi-neutral plasma yields the vorticity equation, which depends on the electric potential, the second field addressed in the model. The coupling between the two equations is crucial to generate an instability and characterises the mechanism that governs the instability. The aim of the present work is to complete the TOKAM2D framework by incorporating the drift wave instability considering the seminal Hasegawa-Wakatani model [25]. We revisit the derivation of the model to highlight the changes made in the TOKAM2D model to incorporate the parallel loss term that drive the drift wave instability, Section 2. We show in particular that accounting for temperature fluctuations in the parallel sheath loss term, as published in [26], yields an instability drive that is identical to that of the drift wave model, see section 2.2. The linear dispersion relation is then investigated and control parameter values are chosen such that the interchange instability alone or the drift wave instability alone exhibit the same growth rate Section 3. The non-linear simulations with the selected control parameters are presented in Section 4. Finally, Section 5 closes the paper.

## 2. Plasma turbulence model

### 2.1. *Transverse transport model*

Two linear turbulence mechanisms have been used to investigate plasma turbulence with a minimum transport model. The linear drive is either that of drift wave turbulence as investigated in the seminal paper [25], or interchange-like [20]. These models couple two fields and use fluid equations restricting the phase space dimension to the directions transverse to the magnetic field, typically *r* a radial and *θ* an angle coordinate. Parallel transport is taken into account by loss terms that are simplified to only depend on the chosen fields. The transport model then takes the form of a conservation equation for either the density *n* or the scalar pressure *p* = *nT* proportional to the internal energy density, or combining these equations a transport equation for the thermal energy *T*. We address here the transport model in the cold ion limit so that the thermal energy *T* is that of the electrons. The second field of these models is the vorticity *w* and its evolution equation is given by the charge balance equation. The drift expansion is used to determine the transverse transport, at lowest order the electric *v*<sub>*E*</sub> and diamagnetic *v*<sub>*</sub> drift velocities, and at higher order the ion polarisation velocity. In the cold ion limit, the vorticity *W*, is determined by the electric potential *φ*, *W* = ∇<sup>2</sup>*φ*. In this limit the transverse transport model for plasma turbulence is akin to the Rayleigh-B'enard 2D turbulence model [27] for neutral fluids, as underlined in Ref.[11, 28]. Transverse transport is then the competition
[7, 8, 9] are also contemplated to have fast access to simulations of experiments and forecast experimental scenarios. We propose here an intermediate model providing first principle analysis of turbulent transport but reduced to the simplest 2D, 2 fields model [10, 11] that make such a model very attractive to generate large data bases with scans of control parameters [12]. Such a model is flux driven and characterised by avalanche transport [10, 13], also reported as "blob" transport [14]. Similar transport behaviour has been observed for heat in the gyrokinetic framework [15]. This model has also allowed us to recover staircase confinement [16], first analysed in the gyrokinetic framework [17]. Our claim in such an approach is that physics problems such as turbulent transport must be robust so that a strongly simplified fluid code should capture qualitatively the behaviour of the system, provided the problem at hand is not specifically governed by kinetic properties, such as the problem of the sheath physics [18, 19]. The model that is used is known as the TOKAM2D model [20]. It was developed to investigate interchange like turbulence in the SOL [21, 22]. Recently it has been extended to analyse the interface transport barrier between the edge plasma with closed field lines and the SOL with sheath boundary conditions [23], considering the interchange instability mechanism and zonal flow generation [24]. The model evolves a transported field, typically the density at constant thermal energy, set out of thermodynamic equilibrium either by a prescribed gradient -as used to investigate linear stability or gradient driven turbulence- or by a radially localised source such that the mean flux is imposed and the field gradient an output of turbulent transport. The latter regime has been named flux driven. We show in this paper that the transported field can be the thermal energy, assuming the density to be constant. The charge balance equation for a quasi-neutral plasma yields the vorticity equation, which depends on the electric potential, the second field addressed in the model. The coupling between the two equations is crucial to generate an instability and characterises the mechanism that governs the instability. The aim of the present work is to complete the TOKAM2D framework by incorporating the drift wave instability considering the seminal Hasegawa-Wakatani model [25]. We revisit the derivation of the model to highlight the changes made in the TOKAM2D model to incorporate the parallel loss term that drive the drift wave instability, Section 2. We show in particular that accounting for temperature fluctuations in the parallel sheath loss term, as published in [26], yields an instability drive that is identical to that of the drift wave model, see section 2.2. The linear dispersion relation is then investigated and control parameter values are chosen such that the interchange instability alone or the drift wave instability alone exhibit the same growth rate Section 3. The non-linear simulations with the selected control parameters are presented in Section 4. Finally, Section 5 closes the paper.

## 2. Plasma turbulence model

### 2.1. Transverse transport model

Two linear turbulence mechanisms have been used to investigate plasma turbulence with a minimum transport model. The linear drive is either that of drift wave turbulence as investigated in the seminal paper [25], or interchange-like [20]. These models couple two fields and use fluid equations restricting the phase space dimension to the directions transverse to the magnetic field, typically *r* a radial and *θ* an angle coordinate. Parallel transport is taken into account by loss terms that are simplified to only depend on the chosen fields. The transport model then takes the form of a conservation equation for either the density *n* or the scalar pressure *p* = *nT* proportional to the internal energy density, or combining these equations a transport equation for the thermal energy *T*. We address here the transport model in the cold ion limit so that the thermal energy *T* is that of the electrons. The second field of these models is the vorticity *w* and its evolution equation is given by the charge balance equation. The drift expansion is used to determine the transverse transport, at lowest order the electric *v*<sub>*E*</sub> and diamagnetic *v*<sub>*</sub> drift velocities, and at higher order the ion polarisation velocity. In the cold ion limit, the vorticity *W*, is determined by the electric potential *φ*, *W* = ∇<sup>2</sup>*φ*. In this limit the transverse transport model for plasma turbulence is akin to the Rayleigh-B'enard 2D turbulence model [27] for neutral fluids, as underlined in Ref.[11, 28]. Transverse transport is then the competition
between convection and diffusion. Above a threshold in the departure from thermodynamic equilibrium, a bifurcation takes place from purely diffusive transport to convection dominated transport. Let $f$ be the field that is transported, and $Q_\parallel$ and $Q_\perp$ the parallel and transverse flux of $f$. For the transverse flux one then has:

$$Q_\perp = f v_{f,\perp} - D_{f,\perp} \nabla_\perp f \quad (1a)$$

where $v_{f,\perp}$ and $D_{f,\perp}$ are the transverse convective velocity and diffusion coefficient respectively. The evolution of $f$ is then determined by the divergence of the fluxes $Q_\perp$ and $Q_\parallel$ together with local evolution properties such as a localised source term that governs the departure from equilibrium. We shall assume that the diffusion coefficients $D_{f,\perp}$ are constant so that:

$$\nabla_\perp \cdot Q_\perp = \nabla_\perp \cdot f v_{f,\perp} - D_{f,\perp} \nabla_\perp^2 f \quad (1b)$$

The electric and diamagnetic drift velocities are prop ortional to $B \times \nabla g$. Neglecting the terms depending on $\nabla \times B$, one then finds that $\nabla_\perp \cdot f v_{f,\perp} = B \cdot \nabla g \times \nabla(hf) = B g, hf$ therefore proportional to the Poisson bracket of $g$ and $hf$. At lowest order in the drift expansion, one obtains the evolution equation for the electron density $n$:

$$\partial_t n + B \frac{\varphi}{B^2} \frac{\mathbf{i}}{n} + \frac{B}{e} \frac{\mathbf{h}}{B^2} \frac{\mathbf{i}}{n}, p_e - D_\perp \nabla_\perp^2 n = S_n - \nabla_\perp \Gamma_\perp \quad (2a)$$

The charge balance equation, computing the divergence of the transverse current including the ion polarisation current and enforcing the Boussinesq approximation yields the vorticity equation.

$$\partial_t W + \frac{1}{B} \frac{\varphi}{n} \frac{\mathbf{i}}{n} + \frac{B}{n} \frac{\mathbf{h}}{B^2} \frac{\mathbf{i}}{n}, p_e - v_\perp \nabla_\perp^2 W = \frac{1}{n} \nabla_\perp \Gamma_\perp \quad (2b)$$

It is possible to replace the transport equation for the density, taken at constant electron temperature $T_e$, by the transport equation for the electron thermal energy $T_e$ at constant density. The parallel loss term is then typically the parallel divergence of the electron energy flux including parallel heat diffusion a contribution of the parallel convection of thermal energy.

Without changing notation, one the density $n$, magnetic field $B$, length scale $L$, time $t$, electric potential $\varphi$ by $n_0$, $B_0$, $L_0$, $\tau_0$, $T_e/e$ respectively. Given the sound velocity $c_s = T_e/m_i$, and setting $\tau_0 = 1/\Omega_0 = m_i/eB_0$ the inverse ion cyclotron frequency, one defines the normalisation length scale $L_0 = \rho_0 = c_0/\Omega_0$, $\rho_0$, hence the reference ion Larmor radius. One then finds that the diffusion process is normalised by the so-called Bohm diffusion $D_B = T_e/(eB_0) = \rho_0^2\Omega_0$.

$$\partial_t n + \varphi, n - n, \frac{1}{B^2} - D_\perp \nabla_\perp^2 n = S_n - \frac{1}{n_0 \Omega_0} \nabla_\parallel \Gamma_\parallel \quad (3a)$$

$$\partial_t W + \varphi, W + \frac{1}{n} \frac{\mathbf{h}}{B^2}, n - v_\perp \nabla_\perp^2 W = \frac{\rho_0^2}{\Omega_0 \varphi_0} \frac{1}{n} \nabla_\parallel \mathbf{J}_\parallel \quad (3b)$$

This set of equations, coupling the density or the thermal energy to the vorticity, describes 2D plasma turbulence.

## 2.2. Parallel loss term

Two kinds of parallel loss terms can be addressed: that governed by parallel transport or that governed by the sheath boundary conditions as in the Scrape-off Layer of magnetic fusion devices. Let us first address the sheath parallel losses. The boundary conditions come into play when averaging the equations in the parallel direction assuming that all fields in Eq. (3a) and Eq. (3b) are constant along a field line, the so-called flute approximation. The parallel loss term is then
between convection and diffusion. Above a threshold in the departure from thermodynamic equilibrium, a bifurcation takes place from purely diffusive transport to convection dominated transport. Let $f$ be the field that is transported, and $Q_\parallel$ and $Q_\perp$ the parallel and transverse flux of $f$. For the transverse flux one then has:

$$Q_\perp = f v_{f,\perp} - D_{f,\perp} \nabla_\perp f \quad (1a)$$

where $v_{f,\perp}$ and $D_{f,\perp}$ are the transverse convective velocity and diffusion coefficient respectively. The evolution of $f$ is then determined by the divergence of the fluxes $Q_\perp$ and $Q_\parallel$ together with local evolution properties such as a localised source term that governs the departure from equilibrium. We shall assume that the diffusion coefficients $D_{f,\perp}$ are constant so that:

$$\nabla_\perp \cdot Q_\perp = \nabla_\perp \cdot f v_{f,\perp} - D_{f,\perp} \nabla_\perp^2 f \quad (1b)$$

The electric and diamagnetic drift velocities are prop ortional to $B \times \nabla g$. Neglecting the terms depending on $\nabla \times B$, one then finds that $\nabla_\perp \cdot f v_{f,\perp} = B \cdot \nabla g \times \nabla(hf) = B g, hf$ therefore proportional to the Poisson bracket of $g$ and $hf$. At lowest order in the drift expansion, one obtains the evolution equation for the electron density $n$:

$$\partial_t n + B \frac{\phi}{B^2} + \frac{n}{B^2} \frac{\phi}{e} + \frac{B}{e} \frac{\phi}{B^2} p_e - D_\perp \nabla_\perp^2 n = S_n - \nabla_\perp \Gamma_\parallel \quad (2a)$$

The charge balance equation, computing the divergence of the transverse current including the ion polarisation current and enforcing the Boussinesq approximation yields the vorticity equation.

$$\partial_t W + \frac{1}{B} \frac{\phi}{n} \frac{\phi}{B^2} + \frac{B}{n} \frac{\phi}{B^2} p_e - v_\perp \nabla_\perp^2 W = \frac{1}{n} \nabla_\parallel J_\parallel \quad (2b)$$

It is possible to replace the transport equation for the density, taken at constant electron temperature $T_e$, by the transport equation for the electron thermal energy $T_e$ at constant density. The parallel loss term is then typically the parallel divergence of the electron energy flux including parallel heat diffusion a contribution of the parallel convection of thermal energy.

Without changing notation, one the density $n$, magnetic field $B$, length scale $L$, time $t$, electric potential $\phi$ by $n_0$, $B_0$, $L_0$, $\tau_0$, $T_e/e$ respectively. Given the sound velocity $c_s = T_e/m_i$, and setting $\tau_0 = 1/\Omega_0 = m_i/eB_0$ the inverse ion cyclotron frequency, one defines the normalisation length scale $L_0 = \rho_0 = c_0/\Omega_0$, $\rho_0$, hence the reference ion Larmor radius. One then finds that the diffusion process is normalised by the so-called Bohm diffusion $D_B = T_e/(eB_0) = \rho_0^2\Omega_0$.

$$\partial_t n + \frac{\phi}{n} \frac{\phi}{B^2} + \frac{1}{n} \frac{\phi}{B^2} \frac{\phi}{n} - D_\perp \nabla_\perp^2 n = S_n - \frac{1}{n_0 \Omega_0} \nabla_\parallel \Gamma_\parallel \quad (3a)$$

$$\partial_t W + \frac{\phi}{n} W + \frac{1}{n} \frac{\phi}{B^2} \frac{\phi}{n} - v_\perp \nabla_\perp^2 W = \frac{\rho_0^2}{\Omega_0 \phi_0} \frac{1}{n} \nabla_\parallel J_\parallel \quad (3b)$$

This set of equations, coupling the density or the thermal energy to the vorticity, describes 2D plasma turbulence.

## 2.2. Parallel loss term

Two kinds of parallel loss terms can be addressed: that governed by parallel transport or that governed by the sheath boundary conditions as in the Scrape-off Layer of magnetic fusion devices. Let us first address the sheath parallel losses. The boundary conditions come into play when averaging the equations in the parallel direction assuming that all fields in Eq. (3a) and Eq. (3b) are constant along a field line, the so-called flute approximation. The parallel loss term is then
determined by the divergence of an outflux of the form $\nabla_{\parallel}Q_{\parallel}$. The parallel average of the latter yields:

$$ \frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \nabla_{\parallel} Q_{\parallel} = \frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \partial_{s} Q_{\parallel} = \frac{1}{L_{\parallel}} \oint (s = L_{\parallel}) - Q_{\parallel}(s = 0) \quad (4a) $$

We further assume that $s = 0$ corresponds to the stagnation point such that $Q_{\parallel}(s = 0) = 0$ and that the flux leaving the SOL is determined by the sheath conditions, $Q_{\parallel}(s = L_{\parallel}) = Q_{\parallel}^{\text{sheath}}$ so that:

$$ \frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \nabla_{\parallel} Q_{\parallel} = \frac{1}{L_{\parallel}} Q_{\parallel}^{\text{sheath}} \quad (4b) $$

Using the Mach = 1 Bohm boundary condition, grounded wall potential, the electron particle flux to the wall can be determined by the values of the fields at the sheath entrance. This expression depends on the parameter $\Phi_s$ that accounts for the difference of mobility between electrons and ions addressed in a kinetic framework. A strongly simplified description yields $\Phi_s = \frac{1}{2}\log(m_i/m_e)$ therefore a constant depending on the mass ratio, $m_i$ and $m_e$ being the ion and electron mass. For the evolution of the density field Eq. (3a), the loss term is then determined by:

$$ \Gamma_{e|}^{\text{sheath}} = n c_s e^{\Phi_s - e\varphi/T_e} \quad (5a) $$

Linearising this expression for small density $\tilde{n}$ and electric potential $\tilde{\varphi}$ fluctuations, and given $e\varphi = \Phi_s T_e$, one obtains:

$$ \Gamma_{e|}^{\text{sheath}} = \bar{n} \bar{c}_s \left( 1 + \frac{\tilde{n}}{\bar{n}} - \frac{e \tilde{\varphi}}{T_e} \right) \quad (5b) $$

The ion particle outflux at the sheath is $nc_s$ so that the linearised ion flux is $\bar{n}c_s + n/\bar{n}$. The linearised parallel current loss to the wall is then:

$$ j_{e|}^{\text{sheath}} = e \bar{n} \bar{c}_s \frac{e \tilde{\varphi}}{T_e} \quad (5c) $$

One then finds the linear expressions depending on $n$ and $\varphi$ for the parallel outflux that appear in Eq. (3):

$$ \Gamma_{e|}^{\text{sheath}} = n_0 c_0 \frac{n}{n_0} - \frac{e \varphi}{T_0} \quad (6a) $$

$$ j_{e|}^{\text{sheath}} = e n_0 c_0 \frac{e \varphi}{T_0} \quad (6b) $$

When addressing the thermal energy transport at constant density and fluctuating thermal energy $\tilde{T}_e$, the linearised parallel current and relevant energy flux $q_{e|}^{\text{sheath}}$ to the wall depend on $\tilde{T}_e$ and $e\tilde{\varphi}/T_e$:

$$ j_{e|}^{\text{sheath}} = e n \bar{c}_s - \frac{e \tilde{\varphi}}{T_e} - \Phi_s \frac{\tilde{T}_e}{T_e} \quad (7a) $$

$$ q_{e|}^{\text{sheath}} = \gamma_r \bar{n} \bar{c}_s \bar{T}_e + 1 + \frac{\tilde{T}_e}{2 \bar{T}_e} - \frac{e \tilde{\varphi}}{\bar{T}_e} + \Phi_s \frac{\tilde{T}_e}{\bar{T}_e} \quad (7b) $$
determined by the divergence of an outflux of the form $\nabla_{\parallel}Q_{\parallel}$. The parallel average of the latter yields:

$$ \frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \nabla_{\parallel}Q_{\parallel} = \frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \partial_{s}Q_{\parallel} = \frac{1}{L_{\parallel}} \oint (s = L_{\parallel}) - Q_{\parallel}(s = 0) \quad (4a) $$

We further assume that $s = 0$ corresponds to the stagnation point such that $Q_{\parallel}(s = 0) = 0$ and that the flux leaving the SOL is determined by the sheath conditions, $Q_{\parallel}(s = L_{\parallel}) = Q_{\parallel}^{\text{sheath}}$ so that:

$$ \frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \nabla_{\parallel}Q_{\parallel} = \frac{1}{L_{\parallel}} Q_{\parallel}^{\text{sheath}} \quad (4b) $$

Using the Mach = 1 Bohm boundary condition, grounded wall potential, the electron particle flux to the wall can be determined by the values of the fields at the sheath entrance. This expression depends on the parameter $\Phi_s$ that accounts for the difference of mobility between electrons and ions addressed in a kinetic framework. A strongly simplified description yields $\Phi_s = \frac{1}{2}\text{Log}(m_i/m_e)$ therefore a constant depending on the mass ratio, $m_i$ and $m_e$ being the ion and electron mass. For the evolution of the density field Eq. (3a), the loss term is then determined by:

$$ \Gamma_{e|}^{\text{sheath}} = n c_s e^{\Phi_s - e\varphi/T_e} \quad (5a) $$

Linearising this expression for small density $\tilde{n}$ and electric potential $\tilde{\varphi}$ fluctuations, and given $e\varphi = \Phi_s T_e$, one obtains:

$$ \Gamma_{e|}^{\text{sheath}} = \bar{n} \bar{c}_s 1 + \frac{\tilde{n}}{\bar{n}} - \frac{e \tilde{\varphi}}{T_e} \quad (5b) $$

The ion particle outflux at the sheath is $nc_s$ so that the linearised ion flux is $\bar{n}c_s 1 + n/\bar{n}$. The linearised parallel current loss to the wall is then:

$$ j_{e|}^{\text{sheath}} = e \bar{n} \bar{c}_s \frac{e \tilde{\varphi}}{T_e} \quad (5c) $$

One then finds the linear expressions depending on $n$ and $\varphi$ for the parallel outflux that appear in Eq. (3):

$$ \Gamma_{e|}^{\text{sheath}} = n_0 c_0 \frac{n}{n_0} - \frac{e \varphi}{T_0} \quad (6a) $$

$$ j_{e|}^{\text{sheath}} = e n_0 c_0 \frac{e \varphi}{T_0} \quad (6b) $$

When addressing the thermal energy transport at constant density and fluctuating thermal energy $\tilde{T}_e$, the linearised parallel current and relevant energy flux $q_{e|}^{\text{sheath}}$ to the wall depend on $\tilde{T}_e$ and $e\tilde{\varphi}/T_e$:

$$ j_{e|}^{\text{sheath}} = e n \bar{c}_s - \frac{e \tilde{\varphi}}{T_e} - \Phi_s \frac{\tilde{T}_e}{T_e} \quad (7a) $$

$$ q_{e|}^{\text{sheath}} = \gamma_r \bar{n} \bar{c}_s \bar{T}_e + 1 + \frac{\tilde{T}_e}{2 \bar{T}_e} - \frac{e \tilde{\varphi}}{\bar{T}_e} + \Phi_s \frac{\tilde{T}_e}{T_e} \quad (7b) $$
One can determine $\gamma_r$ given the sheath energy transmission factor, typically $\gamma_r = \frac{2}{3}\Phi_s + \frac{1}{3}$. The parallel fluxes for the coupled transport equations are then:

$$j_{\parallel}^{\text{sheath}} = en \oint_0 \frac{e\varphi}{T_0} - \Phi \frac{T_e}{sT_0} \qquad (8a)$$

$$q_{\parallel}^{\text{sheath}} = \gamma \rho \oint_0 \frac{e\varphi}{T_0} + \Phi s + \frac{3}{2} \frac{T_e}{T_0} \qquad (8b)$$

Normalising these loss terms and replacing $n/n_0$ ot $T_e/T_0$ by $f$ to obtain a general expression we find therefore:

$$\frac{1}{L_{\parallel}} \int_{L_{\parallel}}^{L} ds \nabla_{\parallel} Q_{\parallel} = -\sigma_f f + \sigma_f \bar{\varphi} \qquad (9a)$$

$$\frac{1}{L_{\parallel}} \int_{0}^{L_{\parallel}} ds \nabla_{\parallel} j_{\parallel} = \sigma_{\varphi} \varphi - \sigma_{\varphi} \bar{\varphi} \qquad (9b)$$

The $\sigma$ terms are diagonal while the $\bar{\sigma}$ are non-diagonal coupling terms. When addressing the case of density transport $f = n$ with $T_e$ constant, the coefficients are:

$$\sigma_f = \frac{\rho_0}{L_{\parallel}}; \quad \bar{\sigma}_f = \frac{\rho_0}{L_{\parallel}}; \quad \sigma_{\varphi} = \frac{\rho_0}{L_{\parallel}}; \quad \bar{\sigma}_{\varphi} = 0 \qquad (10a)$$

while for the case of thermal energy transport $f = T_e$ with constant density, the coefficients are:

$$\sigma_f = \frac{\rho_0}{L_{\parallel}} \Phi_{s+\frac{3}{2}}; \quad \bar{\sigma}_f = \frac{\rho_0}{L_{\parallel}}; \quad \sigma_{\varphi} = \frac{\rho_0}{L_{\parallel}}; \quad \bar{\sigma}_{\varphi} = \frac{\rho_0}{L_{\parallel}} \Phi_s \qquad (10b)$$

These results hold in the SOL provided the parallel variations are small which allows using the flute approximation and enforces that the parallel losses are governed by the sheath conditions. In the edge plasma or in the high recycling or detached SOL regime this formulation of the parallel loss terms must be modified to account for parallel transport. A crude model can be addressed including collisional friction in the electron momentum balance equation. For adiabatic electrons in the electrostatic regime, the balance equation is then:

$$\nabla_{\parallel} p_e + m_e v n u_{e\parallel} = en \nabla_{\parallel} \varphi \qquad (11a)$$

The friction term is proportional to the collision frequency $v$ and drives the collisional relaxation of the difference between mean electron parallel velocity $u_{e\parallel}$ and mean ion parallel velocity $u_{i\parallel}$ in the limit $u_{i\parallel} \to 0$. We therefore neglect the ion response both to determine both the parallel electron flux and current. In this expression $v$ stands for the electron-ion collision frequency for ions of charge $Z$ and density $n$.

$$v = \frac{4}{3} \frac{\sqrt{2\pi \text{Log}\Lambda}}{Ze^2/(4\pi\epsilon_0)} \left(\frac{n v_{\text{the}}}{T_e^2}\right)^2 \qquad (11b)$$

Neglecting collisions in the limit $v \to 0$ and for constant electron thermal energy $T_e$ then leads to the density dependence used to obtain Eq. (5a). For $v > 0$ the momentum balance equation yields the parallel electron flux driven by the parallel thermodynamic forces.

$$\Gamma_{e\parallel} = n u_{e\parallel} = - \frac{1}{m_e v} \nabla_{\parallel} p_e - en \nabla_{\parallel} \varphi \qquad (11c)$$

Stepping to Fourier space in the parallel direction to determine the divergence of the parallel electron flux, one obtains after linearising the term $n\nabla_{\parallel}\varphi \approx n_0\nabla_{\parallel}\varphi$:

$$-\nabla_{\parallel} \Gamma_{|e} \approx ik_{|e} \hat{\Gamma}_{|e} = - \frac{k^2}{m_e v} \hat{p}_{ek_{|e}} - en_0 \hat{\varphi}_k \qquad (12a)$$
One can determine $\gamma_r$ given the sheath energy transmission factor, typically $\gamma_r = \frac{2}{3}\Phi_s + \frac{1}{3}$. The parallel fluxes for the coupled transport equations are then:

$$j_{\parallel}^{\text{sheath}} = en \oint_0 \frac{e\varphi}{T_0} - \Phi \frac{T_e}{sT_0} \qquad (8a)$$

$$q_{\parallel}^{\text{sheath}} = \gamma \rho \oint_0 \frac{e\varphi}{T_0} + \Phi s + \frac{3}{2} \frac{T_e}{T_0} \qquad (8b)$$

Normalising these loss terms and replacing $n/n_0$ ot $T_e/T_0$ by $f$ to obtain a general expression we find therefore:

$$\frac{1}{L_{\parallel}} \int_{L_{\parallel}}^{L_{\parallel}} ds \nabla_{\parallel} Q_{\parallel} = -\sigma_f f + \sigma_f \bar{\varphi} \qquad (9a)$$

$$\frac{1}{L_{\parallel}} \int_{L_{\parallel}}^{0} ds \nabla_{\parallel} j_{\parallel} = \sigma_{\varphi} \varphi - \sigma_{\varphi} \bar{\varphi} \qquad (9b)$$

The $\sigma$ terms are diagonal while the $\bar{\sigma}$ are non-diagonal coupling terms. When addressing the case of density transport $f = n$ with $T_e$ constant, the coefficients are:

$$\sigma_f = \frac{\rho_0}{L_{\parallel}}; \quad \bar{\sigma}_f = \frac{\rho_0}{L_{\parallel}}; \quad \sigma_{\varphi} = \frac{\rho_0}{L_{\parallel}}; \quad \bar{\sigma}_{\varphi} = 0 \qquad (10a)$$

while for the case of thermal energy transport $f = T_e$ with constant density, the coefficients are:

$$\sigma_f = \frac{\rho_0}{L_{\parallel}} \Phi_{s+\frac{3}{2}}; \quad \bar{\sigma}_f = \frac{\rho_0}{L_{\parallel}}; \quad \sigma_{\varphi} = \frac{\rho_0}{L_{\parallel}}; \quad \bar{\sigma}_{\varphi} = \frac{\rho_0}{L_{\parallel}} \Phi_s \qquad (10b)$$

These results hold in the SOL provided the parallel variations are small which allows using the flute approximation and enforces that the parallel losses are governed by the sheath conditions. In the edge plasma or in the high recycling or detached SOL regime this formulation of the parallel loss terms must be modified to account for parallel transport. A crude model can be addressed including collisional friction in the electron momentum balance equation. For adiabatic electrons in the electrostatic regime, the balance equation is then:

$$\nabla_{\parallel} p_e + m_e v n u_{e\parallel} = en \nabla_{\parallel} \varphi \qquad (11a)$$

The friction term is proportional to the collision frequency $v$ and drives the collisional relaxation of the difference between mean electron parallel velocity $u_{e\parallel}$ and mean ion parallel velocity $u_{i\parallel}$ in the limit $u_{i\parallel} \to 0$. We therefore neglect the ion response both to determine both the parallel electron flux and current. In this expression $v$ stands for the electron-ion collision frequency for ions of charge $Z$ and density $n$.

$$v = \frac{4}{3} \frac{\sqrt{2\pi \text{Log}\Lambda}}{Ze^2/(4\pi\epsilon_0)} \left(\frac{nv_{\text{the}}}{T_e^2}\right)^2 \qquad (11b)$$

Neglecting collisions in the limit $v \to 0$ and for constant electron thermal energy $T_e$ then leads to the density dependence used to obtain Eq. (5a). For $v > 0$ the momentum balance equation yields the parallel electron flux driven by the parallel thermodynamic forces.

$$\Gamma_{e\parallel} = nu_{e\parallel} = - \frac{1}{m_e v} \nabla_{\parallel} p_e - en \nabla_{\parallel} \varphi \qquad (11c)$$

Stepping to Fourier space in the parallel direction to determine the divergence of the parallel electron flux, one obtains after linearising the term $n\nabla_{\parallel}\varphi \approx n_0\nabla_{\parallel}\varphi$:

$$-\nabla_{\parallel} \Gamma_{|e} \approx ik_{|e} \hat{\Gamma}_{|e} = - \frac{k^2}{m_e v} \hat{p}_{ek_{|e}} - en_0 \hat{\varphi}_k \qquad (12a)$$
From this expression one can step directly to the divergence of the electron current:

$$ \nabla \cdot \mathbf{j}_{\parallel} \approx ik_{\parallel} \hat{\mathbf{j}}_{\parallel} = - \frac{k^2 e}{\parallel m_e \hat{\mathbf{v}} \parallel} e_k \cdot e_{n_o} \hat{\phi} \hat{\mathbf{k}}_{\parallel} \quad (12b) $$

When normalised these expressions lead to:

$$ -\nabla \cdot \Gamma_{\parallel} \approx -C \ n - \tilde{\phi} \quad (13a) $$

$$ \nabla \cdot \mathbf{j}_{\parallel} \approx -C \ \tilde{n} - \tilde{\phi} \quad (13b) $$

$$ C = \frac{\frac{m_i}{m_e} k_{\parallel} L_{\parallel}}{L_{\parallel}^2} \frac{\rho_0 L_{\text{coll}}}{L_{\parallel}^2} \quad (13c) $$

For this calculation we have defined $L_{\text{coll}} = c_0/v \propto T_e^2/(nv_{\text{the}})$. Furthermore, we have replaced the Fourier modes by the normalised fluctuations of the electric potential $\tilde{\phi}$ and density $\tilde{n}$. The fluctuations, typically $\tilde{f}$, are such that $k_{\parallel} = 0$, else $C = 0$, and will be approximated by $\tilde{f} = f - \langle f \rangle$. The average performed in the parallel direction is akin to a flux surface average that removes the parallel divergence contribution. Compared to the sheath conductivity parameters $\sigma \approx \rho_0/L_{\parallel}$, either Eq. (10a) or Eq. (10b), one finds that typically $C/\sigma \approx m_i/m_e (k_{\parallel}L_{\parallel})^2 L_{\text{coll}}/L_{\parallel}$. Whenever $C \neq 1$ the transport tends to enforce the standard electron adiabatic transport conditions. In the SOL, the sheath boundary conditions then govern the electron parallel losses. Conversely, when $C \sigma$, the parallel losses imposed by the parallel transport are small and dwarfed by the sheath losses. The regime of interest is such that $C/\sigma \sim 1$. One must then have $m_i/m_e (k_{\parallel}L_{\parallel})^2 \sim (L_{\parallel}/L_{\text{coll}})$ and therefore $k_{\parallel}L_{\parallel} \sim (m_e/m_i)^{1/4} (L_{\parallel}/L_{\text{coll}})^{1/2}$. Only the large wavelength fluctuations can exhibit non adiabatic electron transport and generate parallel transport losses that can compete with the sheath losses. Conversely, the flute approximation and sheath dominated losses in the SOL will hold provided the collisionality is small enough so that $L_{\text{coll}}/L_{\parallel} \sim (m_e/m_i)^{1/2}$.

### 2.3. Model equations for plasma turbulence

Having defined a generic form for the parallel loss terms, we now address the balance equations combining transverse and parallel transport. We first address the density, in Eq. (3), simplifying the notation for the normalised quantities.

$$ \partial_t n + \varphi, n - B, n, \frac{1}{B^2} \cdot \frac{\mathbf{i}}{D} - D \perp \nabla^2 n = S_n - \sigma_n n - \sigma_n \bar{\varphi} \quad (14a) $$

For the particular case of the problem we address in this paper, the contributions that account for the variation of $B$ are small and do not change qualitatively the behaviour of the system. One can then set $B = 1$ in Eq. (14a) to obtain a simplified evolution equation for the density.

$$ \partial_t n + \varphi, n - D \perp \nabla^2 n = S_n - \sigma_n n - \sigma_n \bar{\varphi} \quad (14b) $$

In this expression, the transverse transport of the electron density $n$ is not addressed in the adiabatic framework. Should one enforce that $n$ is a function of $\mathbf{i}$ the electric potential $\varphi$, then the cross-field convection described by the Poisson bracket $\varphi, n$ vanishes altogether, removing both turbulent transport and coupling to the vorticity equation. However, when addressing the parallel transport, the adiabatic electron response has been used, namely the electrons are assumed in mechanical equilibrium on the time scale of the ion transport so that the time derivative of the electron momentum is set to zero. It is to be underlined that although the electron parallel transport is addressed in the adiabatic framework, it does not enforce that the electrons are Boltzmannian and the electron density a function of the electric potential.
From this expression one can step directly to the divergence of the electron current:

$$ \nabla \cdot \mathbf{j}_e \approx ik_\parallel \hat{\mathbf{j}}_e = - \frac{k^2 e}{\| \mathbf{p}_e \|^2} e_k \cdot e_n \hat{\phi}_k \quad (12b) $$

When normalised these expressions lead to:

$$ -\nabla \cdot \Gamma_\parallel \approx -C \ n - \tilde{\phi} \quad (13a) $$

$$ \nabla \cdot \mathbf{j}_e \approx -C \ \tilde{n} - \tilde{\phi} \quad (13b) $$

$$ C = \frac{\frac{m_e}{m_e} k_\parallel L_\parallel}{L_\parallel^2} \frac{\rho_0 L_{coll}}{L_\parallel^2} \quad (13c) $$

For this calculation we have defined $L_{coll} = c_0/v \propto T_e^2/(nv_{the})$. Furthermore, we have replaced the Fourier modes by the normalised fluctuations of the electric potential $\tilde{\phi}$ and density $\tilde{n}$. The fluctuations, typically $\tilde{f}$, are such that $k_\parallel = 0$, else $C = 0$, and will be approximated by $\tilde{f} = f - \langle f \rangle$. The average performed in the parallel direction is akin to a flux surface average that removes the parallel divergence contribution. Compared to the sheath conductivity parameters $\sigma \approx \rho_0/L_\parallel$, either Eq. (10a) or Eq. (10b), one finds that typically $C/\sigma \approx m_e/\tilde{m}_e (k_\parallel L_\parallel)^2 L_{coll}/L_\parallel$. Whenever $C \neq 1$ the transport tends to enforce the standard electron adiabatic transport conditions. In the SOL, the sheath boundary conditions then govern the electron parallel losses. Conversely, when $C \sigma$, the parallel losses imposed by the parallel transport are small and dwarfed by the sheath losses. The regime of interest is such that $C/\sigma \sim 1$. One must then have $m_e/\tilde{m}_e (k_\parallel L_\parallel)^2 \sim (L_\parallel/L_{coll})$ and therefore $k_\parallel L_\parallel \sim (m_e/m_i)^{1/4} (L_\parallel/L_{coll})^{1/2}$. Only the large wavelength fluctuations can exhibit non adiabatic electron transport and generate parallel transport losses that can compete with the sheath losses. Conversely, the flute approximation and sheath dominated losses in the SOL will hold provided the collisionality is small enough so that $L_{coll}/L_\parallel \sim (m_e/m_i)^{1/2}$.

### 2.3. Model equations for plasma turbulence

Having defined a generic form for the parallel loss terms, we now address the balance equations combining transverse and parallel transport. We first address the density, in Eq. (3), simplifying the notation for the normalised quantities.

$$ \partial_t n + \varphi, n - B, n, \frac{1}{B^2} \cdot \frac{\mathbf{i}}{B^2} - D \perp \nabla^2 n = S_n - \sigma_n n - \sigma_n \bar{\varphi} \quad (14a) $$

For the particular case of the problem we address in this paper, the contributions that account for the variation of $B$ are small and do not change qualitatively the behaviour of the system. One can then set $B = 1$ in Eq. (14a) to obtain a simplified evolution equation for the density.

$$ \partial_t n + \varphi, n - D \perp \nabla^2 n = S_n - \sigma_n n - \sigma_n \bar{\varphi} \quad (14b) $$

In this expression, the transverse transport of the electron density $n$ is not addressed in the adiabatic framework. Should one enforce that $n$ is a function of $\mathbf{i}$ the electric potential $\varphi$, then the cross-field convection described by the Poisson bracket $\varphi, n$ vanishes altogether, removing both turbulent transport and coupling to the vorticity equation. However, when addressing the parallel transport, the adiabatic electron response has been used, namely the electrons are assumed in mechanical equilibrium on the time scale of the ion transport so that the time derivative of the electron momentum is set to zero. It is to be underlined that although the electron parallel transport is addressed in the adiabatic framework, it does not enforce that the electrons are Boltzmannian and the electron density a function of the electric potential.
The latter regime is only achieved when the electron parallel current is set to be null. Since the parallel transport is taken into account by the reaction terms on the right hand side of Eq.( 14b), — σ<sub>n</sub>n + σ̄<sub>n</sub>φ, the problem that is addressed is 2D in a plane transverse to the magnetic field, with a radial coordinate and a poloidal angle that labels the field lines. For simplicity, we assume a slab geometry and define the normalised coordinates x = r/ρ<sub>0</sub> and y = aθ/ρ<sub>0</sub>. The evolution equation Eq.( 14b) exhibits a single non-linear term, the divergence of the transverse particle flux written as the Poisson bracket φ, n. To address this term in the linear analysis, we split the density n into a mean n̄ and fluctuating ˜n contributions. The mean is defined as the average on both time and poloidal angle, it therefore only depends on the radial coordinate x. With this separation, the Poisson bracket is the sum of three contributions. The non linear contribution in terms of the fluctuations φ̃, ˜n, and two linearised terms with respect to the fluctuations, φ̃, ˜n = —∂<sub>y</sub> φ̃ ∇<sub>⊥</sub> n̄ and φ̄, ˜n = xφ̄ ∂<sub>y</sub> ˜n. The former linear contribution is of particular importance since it couples the density evolution to the electric potential evolution whatever the right hand side reaction terms. To highlight this term, and replace the source term S in gradient driven regimes, we add the term ∂<sub>y</sub> φ̃ n/L<sub>n</sub> with control parameter 1/L<sub>n</sub> = —∇<sub>⊥</sub> n/n. One then obtains:

$$
\partial_t n + \frac{n}{L_n} \partial_i \varphi + \varphi, n - D \nabla^2 n = S_n - \sigma^n n + \bar{\sigma}^n \varphi \quad (14c)
$$

With the control parameters 1/L_n and S_n one can switch from a flux driven mode with 1/L_n → 0 and a radially localised particle source S_n to a gradient driven mode with S_n = 0 and a finite prescribed value for L_n. We now address the transport of vorticity W = ∇⊥^2 φ.

$$
\partial_t W + \varphi, W - B n, \frac{1}{B^2} - \nu_{\perp} \nabla_{\perp}^2 W = S_W - \bar{\sigma}^\rho n + \sigma^\varphi \varphi \quad (15a)
$$

In the chosen slab geometry, and restricting the poloidal region to the equatorial neighbourhood
on the low field side, the term 1/B = R/RB is proportional to R and therefore to x. The Poisson
bracket n, 1/B² is then proportional to ∂_y n and the order of magnitude of the proportionality
factor is ρ₀/R₀ therefore comparable to ρ₀/L|_‖. The exact calculation depends on the problem
of interest. One either assumes a toroidally localised region, as done when one determines the
parallel losses using the local transport properties, or one performs an average in the parallel
direction so that the parallel losses are governed by the sheath losses. One introduces the control
parameter g so that B n, 1/B² = g∂_y n to encompass these various cases with g ∝ ρ₀/L|_‖. One
furthermore considers physics without a vorticity source so that S_W = 0. The vorticity equation
can then be written in a generic form, analogous to that addressed for the Rayleigh-B'enard
instability.

$$
\partial_t(W) + \varphi, W + g\partial_y n - v_{\perp} \nabla_{\perp}^2(W) = -\sigma_{\varphi} n + \sigma_{\varphi} \varphi \quad (15b)
$$

The system we address to investigate SOL turbulence is therefore:

$$
\partial_t n + \frac{1}{L_n} \partial_y \varphi + \varphi, n - D_{\perp} \nabla^2 n = S_n - \sigma_n n + \sigma_{\varphi} \varphi \quad (16a)
$$

$$
\partial_t W + g \partial_y n + \varphi, W - v_{\perp} \nabla^2 W = -\bar{\sigma}_{\varphi} n + \sigma_{\varphi} \varphi \quad (16b)
$$

The control parameter *g* allows one to address interchange-like plasma turbulence, while the parameter σ<sub>φ</sub> is the key control parameter for drift wave turbulence. For SOL conditions with constant thermal energy, we shall consider the interchange instability assuming σ<sub>n</sub> = 0 and σ̄<sub>φ</sub> = 0. Conversely for drift wave turbulence, we set *g* = 0 and retain σ<sub>n</sub> and σ̄<sub>φ</sub> with σ<sub>n</sub> = σ̄<sub>φ</sub> = σ<sub>n</sub> = σ̄<sub>φ</sub> = C. For SOL conditions the latter case would be similar to the model named "temperature gradient instability induced by conducting end walls" by Berk et al. presented in Ref.[26] provided one replaces *n* by *T*<sub>e</sub> in the present model. One finds therefore
The latter regime is only achieved when the electron parallel current is set to be null. Since the parallel transport is taken into account by the reaction terms on the right hand side of Eq.( 14b), — σ<sub>n</sub>n+σ<sub>n</sub>φ, the problem that is addressed is 2D in a plane transverse to the magnetic field, with a radial coordinate and a poloidal angle that labels the field lines. For simplicity, we assume a slab geometry and define the normalised coordinates x = r/ρ<sub>0</sub> and y = aθ/ρ<sub>0</sub>. The evolution equation Eq.( 14b) exhibits a single non-linear term, the divergence of the transverse particle flux written as the Poisson bracket φ, n. To address this term in the linear analysis, we split the density n into a mean n̄ and fluctuating ˜n contributions. The mean is defined as the average on both time and poloidal angle, it therefore only depends on the radial coordinate x. With this separation, the Poisson bracket is the sum of three contributions. The non linear contribution in terms of the fluctuations φ̃, ˜n, and two linearised terms with respect to the fluctuations, φ̃, ˜n = —∂<sub>y</sub>φ̃ ∇<sub>⊥</sub>n̄ and φ̃, ˜n = xφ̃ ∂<sub>y</sub>˜n. The former linear contribution is of particular importance since it couples the density evolution to the electric potential evolution whatever the right hand side reaction terms. To highlight this term, and replace the source term S in gradient driven regimes, we add the term ∂<sub>y</sub>φ̃ n/L<sub>n</sub> with control parameter 1/L<sub>n</sub> = —∇<sub>⊥</sub>n/n. One then obtains:

$$
\partial_t n + \frac{n}{L_n} \partial_i \varphi + \varphi, n - D \nabla^2 n = S_n - \sigma^n n + \bar{\sigma}^n \varphi \quad (14c)
$$

With the control parameters 1/L_n and S_n one can switch from a flux driven mode with 1/L_n → 0 and a radially localised particle source S_n to a gradient driven mode with S_n = 0 and a finite prescribed value for L_n. We now address the transport of vorticity W = ∇⊥^2 φ.

$$
\partial_t W + \varphi, W - B n, \frac{1}{B^2} - \nu \perp \nabla^2 W = S_W - \bar{\sigma}^\rho n + \sigma^\varphi \varphi \quad (15a)
$$

In the chosen slab geometry, and restricting the poloidal region to the equatorial neighbourhood
on the low field side, the term 1/B = R/RB is proportional to R and therefore to x. The Poisson
bracket n, 1/B² is then proportional to ∂_y n and the order of magnitude of the proportionality
factor is ρ₀/R₀ therefore comparable to ρ₀/L|_‖. The exact calculation depends on the problem
of interest. One either assumes a toroidally localised region, as done when one determines the
parallel losses using the local transport properties, or one performs an average in the parallel
direction so that the parallel losses are governed by the sheath losses. One introduces the control
parameter g so that B n, 1/B² = g∂_y n to encompass these various cases with g ∝ ρ₀/L|_‖. One
furthermore considers physics without a vorticity source so that S_W = 0. The vorticity equation
can then be written in a generic form, analogous to that addressed for the Rayleigh-B'enard
instability.

$$
\partial_t(W) + \varphi, W + g\partial_y n - v_{\perp} \nabla^2(W) = -\sigma_{\bar{\varphi}} n + \sigma_{\varphi} \varphi \quad (15b)
$$

The system we address to investigate SOL turbulence is therefore:

$$
\partial_t n + \frac{1}{L_n} \partial_y \varphi + \varphi, n - D_{\perp} \nabla^2 n = S_n - \sigma_n n + \sigma_{\bar{\varphi}} \varphi \quad (16a)
$$

$$
\partial_t W + g \partial_y n + \varphi, W - v_{\perp} \nabla^2 W = -\sigma_{\bar{\varphi}} n + \sigma_{\varphi} \varphi \quad (16b)
$$

The control parameter *g* allows one to address interchange-like plasma turbulence, while the parameter σ<sub>φ</sub> is the key control parameter for drift wave turbulence. For SOL conditions with constant thermal energy, we shall consider the interchange instability assuming σ<sub>n</sub> = 0 and σ<sub>φ</sub> = 0. Conversely for drift wave turbulence, we set *g* = 0 and retain σ<sub>n</sub> and σ<sub>φ</sub> with σ<sub>n</sub> = σ<sub>φ</sub> = σ<sub>n</sub> = σ<sub>φ</sub> = C. For SOL conditions the latter case would be similar to the model named "temperature gradient instability induced by conducting end walls" by Berk et al. presented in Ref.[26] provided one replaces *n* by *T*<sub>e</sub> in the present model. One finds therefore
that the latter instability is in fact a drift wave instability where the parallel losses are governed by the temperature dependence of the sheath boundary conditions. To investigate the general case with competing interchange and drift waves linear instabilities we shall consider cases with $g/\bar{g} = 0$ and $\bar{\sigma}_\phi/\sigma_\phi = 0$.

### 3. Dispersion relation

To investigate the linear instability we use the system Eq. (16) without enforcing particular values for the $\sigma$ conductivity coefficients. $W$ is the vorticity $W = \nabla \cdot \phi$. We compute the growth rate $\gamma$ with the linearised equations in Fourier space for the gradient driven case, hence for $S_n = 0$ and $1/L_n$ finite.

$$ \gamma \hat{n} + \frac{ik_y}{k_y^n} \hat{\phi}^\wedge + D_\perp k^2 \hat{n} = -\sigma_n \hat{n} + \bar{\sigma}_n \bar{\hat{\phi}} \quad (17a) $$

$$ \gamma \hat{\phi}^\wedge - i \bar{k}^2 g \hat{n} + v_\perp k \hat{\phi}^\wedge = - \bar{k}^2 \bar{\phi}^\wedge + \bar{k}^2 \hat{n} \quad (17b) $$

The dispersion relation determines the condition to achieve a solution different from the trivial solution $\hat{n} = 0$ and $\hat{\phi}^\wedge = 0$, therefore:

$$ \gamma + A_n \quad \gamma + A_\phi \quad -B_n B_\phi = 0 \quad (18a) $$

where:

$$ A_n = Dk^2 + \sigma_n \quad ; \quad A_\phi = vk^2 + \frac{\sigma_\phi}{k^2} \quad (18b) $$

$$ B_n = ik_y \frac{1}{L_n} - \sigma_n \quad ; \quad B_\phi = -ig \frac{k_y}{k^2} \frac{\sigma_\phi}{k^2} \quad (18c) $$

For the specific case $\bar{\sigma}_n = 0$ and $\bar{\sigma}_\phi = 0$, the coefficients $B_n$ and $B_\phi$ are imaginary and therefore $B_n B_\phi$ is real. This particular case corresponds to the linear interchange instability. Conversely, when $g = 0$ and $\bar{\sigma}_\phi = 0$, $B_\phi$ is real and the product $B_n B_\phi$ complex. This property also has an impact on the phase difference of the eigenmodes.

$$ \hat{\phi}^\wedge = \frac{B_\phi}{\gamma + A_\phi} \hat{n} \quad (19) $$

In the interchange case $B_\phi$ is imaginary and $\gamma$ is real so that there is a phase shift of $\pi/2$ between $\hat{n}$ and $\hat{\phi}$. Conversely, when $B_\phi$ is real $\gamma$ is complex and the phase difference departs from $\pi/2$. We now define $A$:

$$ A = \frac{1}{2} A_n + A_\phi \quad (20a) $$

so that the dispersion relation takes the form:

$$ \gamma + A^2 = B_n B_\phi - A_n A_\phi + A^2 = D \quad (20b) $$

The coefficient $D$ can be split into a real part $D_r$ and an imaginary part $D_i$.

$$ D_r = A^2 + \delta \quad ; \quad \delta = \frac{g k_y^2}{L_n k^2} + \frac{\sigma_n \sigma_\phi}{k^2} - A_n A_\phi \quad (20c) $$

$$ D_i = \frac{k_y}{k^2} g \sigma_n - \frac{\sigma_\phi}{L_n} \quad (20d) $$

One can remark that $A \in \mathbb{R}^+$, $A_n A_\phi \in \mathbb{R}^+$, $A - A_n A_\phi \in \mathbb{R}^+$ so that if $B_n B_\phi \in \mathbb{R}^+$ an instability only occurs when $B_n B_\phi \ge A_n A_\phi$. One thus identifies $B_n B_\phi$ as the drive for the instability and
that the latter instability is in fact a drift wave instability where the parallel losses are governed by the temperature dependence of the sheath boundary conditions. To investigate the general case with competing interchange and drift waves linear instabilities we shall consider cases with $g/\bar{g} = 0$ and $\bar{\sigma}_\phi/\sigma_\phi = 0$.

### 3. Dispersion relation

To investigate the linear instability we use the system Eq. (16) without enforcing particular values for the $\sigma$ conductivity coefficients. $W$ is the vorticity $W = \nabla \phi$. We compute the growth rate $\gamma$ with the linearised equations in Fourier space for the gradient driven case, hence for $S_n = 0$ and $1/L_n$ finite.

$$ \gamma \hat{n} + \frac{ik_y}{k_y^n} \hat{\phi}^\wedge + D_\perp k^2 \hat{n} = -\sigma_n \hat{n} + \bar{\sigma}_n \bar{\hat{\phi}} \quad (17a) $$

$$ \gamma \hat{\phi}^\wedge - i \bar{k}^2 g \hat{n} + v_\perp k \hat{\phi}^\wedge = - \bar{k}^2 \bar{\phi}^\wedge + \bar{k}^2 \hat{n} \quad (17b) $$

The dispersion relation determines the condition to achieve a solution different from the trivial solution $\hat{n} = 0$ and $\hat{\phi}^\wedge = 0$, therefore:

$$ \gamma + A_n \quad \gamma + A_\phi \quad -B_n B_\phi = 0 \quad (18a) $$

where:

$$ A_n = Dk^2 + \sigma_n \quad ; \quad A_\phi = vk^2 + \frac{\sigma_\phi}{k^2} \quad (18b) $$

$$ B_n = ik_y \frac{1}{L_n} - \sigma_n \quad ; \quad B_\phi = -ig \frac{k_y}{k^2} \frac{\sigma_\phi}{k^2} \quad (18c) $$

For the specific case $\bar{\sigma}_n = 0$ and $\bar{\sigma}_\phi = 0$, the coefficients $B_n$ and $B_\phi$ are imaginary and therefore $B_n B_\phi$ is real. This particular case corresponds to the linear interchange instability. Conversely, when $g = 0$ and $\bar{\sigma}_\phi = 0$, $B_\phi$ is real and the product $B_n B_\phi$ complex. This property also has an impact on the phase difference of the eigenmodes.

$$ \hat{\phi}^\wedge = \frac{B_\phi}{\gamma + A_\phi} \hat{n} \quad (19) $$

In the interchange case $B_\phi$ is imaginary and $\gamma$ is real so that there is a phase shift of $\pi/2$ between $\hat{n}$ and $\hat{\phi}$. Conversely, when $B_\phi$ is real $\gamma$ is complex and the phase difference departs from $\pi/2$. We now define $A$:

$$ A = \frac{1}{2} A_n + A_\phi \quad (20a) $$

so that the dispersion relation takes the form:

$$ \gamma + A^2 = B_n B_\phi - A_n A_\phi + A^2 = D \quad (20b) $$

The coefficient $D$ can be split into a real part $D_r$ and an imaginary part $D_i$.

$$ D_r = A^2 + \delta \quad ; \quad \delta = \frac{g k_y^2}{L_n k^2} + \frac{\sigma_n \sigma_\phi}{k^2} - A_n A_\phi \quad (20c) $$

$$ D_i = \frac{k_y}{k^2} g \sigma_n - \frac{\sigma_\phi}{L_n} \quad (20d) $$

One can remark that $A \in \mathbb{R}^+$, $A_n A_\phi \in \mathbb{R}^+$, $A - A_n A_\phi \in \mathbb{R}^+$ so that if $B_n B_\phi \in \mathbb{R}^+$ an instability only occurs when $B_n B_\phi \ge A_n A_\phi$. One thus identifies $B_n B_\phi$ as the drive for the instability and
A, A_n and A_φ as the damping processes. In the general case B_nB_φ ∈ C and γ has a real part γ_r that determines the growth rate of the instability and an imaginary part γ_i that determines the mode frequency. We note α + iβ the root of D hence such that α^2 - β^2 = D_r and 2αβ = D_i. One then obtains:

$$ \alpha^2 = \frac{1}{2} D_r + D_r^2 + D_i^2 \quad (21a) $$

so that the largest growth rate is:

$$ \gamma_r = -A + \alpha \quad (21b) $$

with α chosen positive. One recovers here that the coefficient A characterises damping processes since for α = 0, γ_r = -A, therefore negative and leading to an exponentially decaying solution. An instability will occur if α ≥ A, therefore α^2 > A^2. Inspecting Eq. (21a) one finds that D_i^2 is a positive contribution to γ_r and favours the occurrence of an instability. One can also readily notice that if δ > 0, D_r > A^2 so that α^2 > A^2. One finds therefore that the two first terms in the definition of δ Eq. (20c), namely (g/L_n)(k_y^2/k^2) and σ_nσ_φ/k^2 are the driving terms for the instability. One can then proceed the calculation and obtain:

$$ \frac{k_y^2}{k^4} - \frac{g\sigma_n}{L_n} - \frac{\sigma_\phi^2}{L_n} \geq -4A \delta^2 \quad (22a) $$

A sufficient condition for an instability is to verify δ > 0. One can then recover standard cases, that addressed for interchange in the SOL with σ_n = 0 and σ_φ = 0 such that the condition for instability is δ > 0 and therefore for that case:

$$ \frac{g k_y^2}{L_n k^2} > A_n A_\phi \quad (22b) $$

Since A_nA_φ ≥ 0, a threshold for the instability in the departure from thermodynamic equilibrium is required L_n > 0 since g > 0, therefore 1/L_n positive. The linear interchange instability also favours k_y k_x, which maximises the aspect ratio parameter k_y^2/k^2. For a drift wave instability with g = 0 and σ_n = σ_φ = C one obtains:

$$ \frac{C k_y^2}{L_n k^2} \geq 4A^2 k_y^2 A_n A_\phi - \frac{C^2}{k^2} \quad (23a) $$

The marginal instability condition does not depend on the sign of L_n. Should A_n and A_φ not depend on C, then C only accounts for the destabilising effect with a role analogous to g on the left hand side of Eq. (23a) as well as a destabilising role on the right hand side by inhibiting the damping term A_nA_φ. If one now accounts for σ_n = σ_φ = C as done when addressing the standard drift wave instability, one then obtains A_nA_φ = Dvk^4 + C(vk^2 + D) + C^2/k^2 and 2A = (D + v)k^2 + C(1 + 1/k^2). For g = 0, the marginal linear instability condition is then given by:

$$ \frac{1}{|L_n|} > \frac{(D + v)k^4 + C(k^2 + 1)}{C k_y} Dvk^4 + C(vk^2 + D)^{1/2} \quad (23b) $$

One finds therefore that the drift wave instability exhibits a threshold in density gradient and that unlike that for interchange, it does not depend on the sign of the gradient. In that respect interchange can only be unstable in the Low Field side SOL generating a ballooned transport similar to that evidenced experimentally [29, 30, 31, 32]. Conversely, the drift wave turbulence is unstable in both the low and high field side SOL. In Eq. (23b) one finds that at small C the threshold condition is proportional to 1/C while at large C the asymptotic dependence is C^1/2. When scanning C, the threshold condition goes through a minimum, which indicates that for a
A, A_n and A_φ as the damping processes. In the general case B_nB_φ ∈ C and γ has a real part γ_r that determines the growth rate of the instability and an imaginary part γ_i that determines the mode frequency. We note α + iβ the root of D hence such that α^2 - β^2 = D_r and 2αβ = D_i. One then obtains:

$$ \alpha^2 = \frac{1}{2} D_r + D_r^2 + D_i^2 \quad (21a) $$

so that the largest growth rate is:

$$ \gamma_r = -A + \alpha \quad (21b) $$

with α chosen positive. One recovers here that the coefficient A characterises damping processes since for α = 0, γ_r = -A, therefore negative and leading to an exponentially decaying solution. An instability will occur if α ≥ A, therefore α^2 > A^2. Inspecting Eq. (21a) one finds that D_i^2 is a positive contribution to γ_r and favours the occurrence of an instability. One can also readily notice that if δ > 0, D_r > A^2 so that α^2 > A^2. One finds therefore that the two first terms in the definition of δ Eq. (20c), namely (g/L_n)(k_y^2/k^2) and σ_nσ_φ/k^2 are the driving terms for the instability. One can then proceed the calculation and obtain:

$$ \frac{k_y^2}{k^4} - \frac{g\sigma_n}{L_n} - \frac{\sigma_\phi^2}{L_n} \geq -4A \delta^2 \quad (22a) $$

A sufficient condition for an instability is to verify δ > 0. One can then recover standard cases, that addressed for interchange in the SOL with σ_n = 0 and σ_φ = 0 such that the condition for instability is δ > 0 and therefore for that case:

$$ \frac{g k_y^2}{L_n k^2} > A_n A_\phi \quad (22b) $$

Since A_nA_φ ≥ 0, a threshold for the instability in the departure from thermodynamic equilibrium is required L_n > 0 since g > 0, therefore 1/L_n positive. The linear interchange instability also favours k_y k_x, which maximises the aspect ratio parameter k_y^2/k^2. For a drift wave instability with g = 0 and σ_n = σ_φ = C one obtains:

$$ \frac{C k_y^2}{L_n k^2} \geq 4A^2 k_y^2 A_n A_\phi - \frac{C^2}{k^2} \quad (23a) $$

The marginal instability condition does not depend on the sign of L_n. Should A_n and A_φ not depend on C, then C only accounts for the destabilising effect with a role analogous to g on the left hand side of Eq. (23a) as well as a destabilising role on the right hand side by inhibiting the damping term A_nA_φ. If one now accounts for σ_n = σ_φ = C as done when addressing the standard drift wave instability, one then obtains A_nA_φ = Dvk^4 + C(vk^2 + D) + C^2/k^2 and 2A = (D + v)k^2 + C(1 + 1/k^2). For g = 0, the marginal linear instability condition is then given by:

$$ \frac{1}{|L_n|} > \frac{(D + v)k^4 + C(k^2 + 1)}{C k_y} Dvk^4 + C(vk^2 + D)^{1/2} \quad (23b) $$

One finds therefore that the drift wave instability exhibits a threshold in density gradient and that unlike that for interchange, it does not depend on the sign of the gradient. In that respect interchange can only be unstable in the Low Field side SOL generating a ballooned transport similar to that evidenced experimentally [29, 30, 31, 32]. Conversely, the drift wave turbulence is unstable in both the low and high field side SOL. In Eq. (23b) one finds that at small C the threshold condition is proportional to 1/C while at large C the asymptotic dependence is C^1/2. When scanning C, the threshold condition goes through a minimum, which indicates that for a
<figure><img src="image_2.png" /><figcaption>Figure 1: Left hand side, instability growth rate $\gamma_{Dwl} (k_x, k_y)$ with mixed Drift Wave + Interchange conditions, $g = 1.5 \ 10^{-4}$ and $C = 6.0 \ 10^{-5}$ versus the wave numbers $k_x$ and $k_y$. The steps in $k_x$ and $k_y$, $k_{x,min}$ and $k_{y,min}$, are $\pi\rho_*$ with $1/\rho_* = 256$. The black contour is plotted for $\gamma_{Dwl} = 0$. The dashed black line is the line of $k_y$ values at the maximum of the growth rate $\gamma_{Dwl}$ at given $k_x$ and the star symbol locates the value for $k_x = 5k_{x,min}$. Right hand side, maximum growth rate for $k_x \ge 5k_{x,min}$ and all values of $k_y$, blue curve maximum of $\gamma_I(g)$ versus $g$ for Interchange without drift waves, black curve maximum of $\gamma_{DW}(C)$ versus $C$ for Drift Waves without interchange. For $g_c = 1.5 \ 10^{-4}$ and $C_c = 6.0 \ 10^{-5}$, open black square marker and closed blue triangle respectively, the growth rate values are equal $\gamma_I(g_c) = \gamma_{DW}(C_c) = \gamma_c$.</figcaption></figure>

given value of the density gradient $1/|L_n|$ the drift wave turbulence is stable in both the small and large $C$ limit. This property holds for $k/ = 0$. At small values of $k$, the roll-over to the 1/C dependence can only be found for $C$ ($D + \nu k^4$. As for the interchange instability Eq. (22b), the condition Eq. (23b) indicates that the $k_y = 0$ cannot become unstable, the turbulence onset is clearly governed by the loss of poloidal invariance $k_y/ = 0$.

The linear analysis is investigated numerically, varying the wave vectors $k_x$ and $k_y$ as well as the control parameters $g$ and $C$. The control parameters are for all cases $D_\perp = \nu_\perp = 10^{-2}$ and:

* Interchange case, labelled I: $g/ = 0$, $\sigma_n = \sigma_\phi = C_c$ and $\bar{\sigma}_n = \bar{\sigma}_\phi = 0$,

* Drift Wave case, labelled DW: $g = 0$, $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C/ = 0$,

* Drift Wave + Interchange case, labelled Dwl, $g/ = 0$ and $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C/ = 0$.

The wave vectors are stepped respectively by $k_{x,min}$ and $k_{y,min}$ the smallest non-zero chosen value for the wave vectors, here $k_{x,min} = k_{y,min} = \pi\rho_*$ with $1/\rho_* = L_x$ the radial box size normalised by the reference Larmor radius $\rho_0$. For each value of the control parameters, given Eq. (21), one can determine the growth rate $\gamma(k_x, k_y)$, see Figure 1 left hand side. For the chosen example of Dwl linear analysis $g = g_c = 1.5 \ 10^{-4}$ and $C = C_c = 6.0 \ 10^{-5}$. One finds that the growth rate maximum is located at $k_y = 8k_{y,min}$ and $k_x = 0$. For each value of $k_x$ one can determine the maximum of the growth rate when varying $k_y$. The line of these $k_y$ values are plotted on Figure 1 left hand side dash black line. For $k_x = 5k_{x,min}$ one finds that the maximum is at $k_y = 8k_{y,min}$, black star marker. The black contour line on Figure 1 left hand side is that of marginality $\gamma(k_x, k_y) = 0$, on the left hand side, above and below this contour $\gamma(k_x, k_y) < 0$, the modes $k_x, k_y$ are stable. In the following, the reference value chosen to characterise the growth rate $\gamma$ for all values of parameters $g$ and $C$ is the maximum of $\gamma(k_x, k_y)$ for all values of $k_x \ge 5k_{x,min}$ and all values of $k_y$. The characteristic growth rates $\gamma_{DW}(C)$, black curve, and $\gamma_I(g)$, blue curve, are plotted versus $C$ and $g$ on Figure 1 right hand side. One finds that $\gamma_{DW}$
<figure><img src="image_4.png" /><figcaption>Figure 1: Left hand side, instability growth rate $\gamma_{Dwl} (k_x, k_y)$ with mixed Drift Wave + Interchange conditions, $g = 1.5 \ 10^{-4}$ and $C = 6.0 \ 10^{-5}$ versus the wave numbers $k_x$ and $k_y$. The steps in $k_x$ and $k_y$, $k_{x,min}$ and $k_{y,min}$, are $\pi\rho_*$ with $1/\rho_* = 256$. The black contour is plotted for $\gamma_{Dwl} = 0$. The dashed black line is the line of $k_y$ values at the maximum of the growth rate $\gamma_{Dwl}$ at given $k_x$ and the star symbol locates the value for $k_x = 5k_{x,min}$.
Right hand side, maximum growth rate for $k_x \ge 5k_{x,min}$ and all values of $k_y$, blue curve maximum of $\gamma_I(g)$ versus $g$ for Interchange without drift waves, black curve maximum of $\gamma_{DW}(C)$ versus $C$ for Drift Waves without interchange. For $g_c = 1.5 \ 10^{-4}$ and $C_c = 6.0 \ 10^{-5}$, open black square marker and closed blue triangle respectively, the growth rate values are equal $\gamma_I(g_c) = \gamma_{DW}(C_c) = \gamma_c$.</figcaption></figure>

given value of the density gradient $1/|L_n|$ the drift wave turbulence is stable in both the small and large $C$ limit. This property holds for $k/\epsilon = 0$. At small values of $k$, the roll-over to the $1/C$ dependence can only be found for $C/(D+\nu)k^4$. As for the interchange instability Eq. (22b), the condition Eq. (23b) indicates that the $k_y = 0$ cannot become unstable, the turbulence onset is clearly governed by the loss of poloidal invariance $k_y/\epsilon = 0$.

The linear analysis is investigated numerically, varying the wave vectors $k_x$ and $k_y$ as well as the control parameters $g$ and $C$. The control parameters are for all cases $D_\perp = \nu_\perp = 10^{-2}$ and:

* Interchange case, labelled I: $g/\epsilon = 0$, $\sigma_n = \sigma_\phi = C_c$ and $\bar{\sigma}_n = \bar{\sigma}_\phi = 0$,

* Drift Wave case, labelled DW: $g = 0$, $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C/\epsilon = 0$,

* Drift Wave + Interchange case, labelled Dwl, $g/\epsilon = 0$ and $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C/\epsilon = 0$.

The wave vectors are stepped respectively by $k_{x,min}$ and $k_{y,min}$ the smallest non-zero chosen value for the wave vectors, here $k_{x,min} = k_{y,min} = \pi\rho_*$ with $1/\rho_* = L_x$ the radial box size normalised by the reference Larmor radius $\rho_0$. For each value of the control parameters, given Eq. (21), one can determine the growth rate $\gamma(k_x, k_y)$, see Figure 1 left hand side. For the chosen example of Dwl linear analysis $g = g_c = 1.5 \ 10^{-4}$ and $C = C_c = 6.0 \ 10^{-5}$. One finds that the growth rate maximum is located at $k_y = 8k_{y,min}$ and $k_x = 0$. For each value of $k_x$ one can determine the maximum of the growth rate when varying $k_y$. The line of these $k_y$ values are plotted on Figure 1 left hand side dash black line. For $k_x = 5k_{x,min}$ one finds that the maximum is at $k_y = 8k_{y,min}$, black star marker. The black contour line on Figure 1 left hand side is that of marginality $\gamma(k_x, k_y) = 0$, on the left hand side, above and below this contour $\gamma(k_x, k_y) < 0$, the modes $k_x, k_y$ are stable. In the following, the reference value chosen to characterise the growth rate $\gamma$ for all values of parameters $g$ and $C$ is the maximum of $\gamma(k_x, k_y)$ for all values of $k_x \ge 5k_{x,min}$ and all values of $k_y$. The characteristic growth rates $\gamma_{DW}(C)$, black curve, and $\gamma_I(g)$, blue curve, are plotted versus $C$ and $g$ on Figure 1 right hand side. One finds that $\gamma_{DW}$
<figure><img src="image_3.png" /><figcaption>Figure 2: Left hand side, for $k_x = 5k_{x,min}$ growth rate versus $k_y/k_{y,min}$ for Drift Wave (DW), $\gamma_{DW}$ black curve open circles, $g=0$ and $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C_c$, Interchange (I), $\gamma_I$ blue curve open triangles for $g=g_c$, $\sigma_n = \sigma_\phi = C_c$ and $\bar{\sigma}_n = \bar{\sigma}_\phi = 0$ and Drift Wave with Interchange (Dwl), $\gamma_{Dwl}$ red curve closed squares for $g=g_c$, $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C_c$. The control parameters $C_c$ and $g_c$ are chosen so that the maxima of $\gamma_{DW}$ and $\gamma_I$ have the same value. The sum $\gamma_I + \gamma_{DW}$ is observed to peak at a value similar to the peak value of $\gamma_{Dwl}$.</figcaption></figure>

Right hand side, maximum instability growth rate of combined Drift Wave and Interchange instability $\gamma_{Dwl}(C, g)$ in $\log_{10}$ scale, scanning $\log_{10}(C)$ horizontal axis and $\log_{10}(g)$ vertical axis with $g=0$ and $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C=0$. The black-white dashed contour line is for the marginality condition $\gamma_{Dwl}(C, g) = 0$.

becomes negative at large $C$ as expected and decays to the small values of $C$, but will become negative only for $C < 10^{-9}$ given the chosen values of $k_{x,min}$ and $k_{y,min}$. For Interchange, one recovers that $\gamma_I(g)$ increases monotonically with $g$ becoming positive for $g \gtrsim 6 \cdot 10^{-5}$. For the turbulence simulations we select $C = C_c = 6.0 \cdot 10^{-5}$, open black square, and $g = g_c = 1.5 \cdot 10^{-4}$, closed blue triangle, for which the maximum growth rate $\gamma_{DW}(C_c)$ is equal to the maximum growth rate $\gamma_I(g_c)$, see Figure 1 right hand side. The choice that is made $g_c$ and $C_c$ is further analysed of Figure 2 left hand side. The growth rate $\gamma_{DW}$, black curve open circles, and $\gamma_I$, blue curve open triangles, are plotted versus $k_x/k_{x,min}$ for $k_x = 5k_{x,min}$. As chosen, the two growth rates have the same peak value but for different $k_y$, the Drift Wave being shifted towards lower $k_y$ compared to the Interchange case. Combining the two instabilities yields $\gamma_{Dwl}$, red curve closed squares, which close to double growth rate and intermediate peak value of $k_y$. The peak behaviour is found to compare well to $\gamma_{DW}+\gamma_I$ as if the turbulent processes where additive.

The linear stability analysis of Drift waves and Interchange instability mechanisms (Dwl) has been investigated varying the control parameters $C$ and $g$, Figure 2 right hand side where the characteristic growth rate is plotted versus $C$ in log scale lower axis and $g$ in log scale, vertical axis. As observed on Figure 1 right hand side, varying $C$ has little effect. It bounds the unstable region to the right when $C$ becomes large, typically of order $10^{-3}$. This is shown by the vertical dashed black and white line that indicates the marginality condition $\gamma_{Dwl} = 0$. For $g \to 0$, the horizontal dashed black and white line indicates the stabilising effect for small values of $C$. The unstable region governed by the drift wave instability appears to extend for $C$ in the range $10^{-6} \le C \le 10^{-3}$. The main variation of the growth rate is governed by the variation of $g$. Finally the selected conditions $C_c$ and $g_c$, such that the drift wave and interchange instability yield the same value for the characteristic growth rate is indicated by the open black marker. This point lies in a region of rather weak growth rate - $10^{-3}$ agreeing with the peak value indicated on Figure 2 left hand side.
<figure><img src="image_5.png" /><figcaption>Figure 2: Left hand side, for $k_x = 5k_{x,min}$ growth rate versus $k_y/k_{y,min}$ for Drift Wave (DW), $\gamma_{DW}$ black curve open circles, $g=0$ and $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C_c$, Interchange (I), $\gamma_I$ blue curve open triangles for $g=g_c$, $\sigma_n = \sigma_\phi = C_c$ and $\bar{\sigma}_n = \bar{\sigma}_\phi = 0$ and Drift Wave with Interchange (Dwl), $\gamma_{Dwl}$ red curve closed squares for $g=g_c$, $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C_c$. The control parameters $C_c$ and $g_c$ are chosen so that the maxima of $\gamma_{DW}$ and $\gamma_I$ have the same value. The sum $\gamma_I + \gamma_{DW}$ is observed to peak at a value similar to the peak value of $\gamma_{Dwl}$.</figcaption></figure>

Right hand side, maximum instability growth rate of combined Drift Wave and Interchange instability $\gamma_{Dwl}(C, g)$ in $log_{10}$ scale, scanning $log_{10}(C)$ horizontal axis and $log_{10}(g)$ vertical axis with $g=0$ and $\sigma_n = \sigma_\phi = \bar{\sigma}_n = \bar{\sigma}_\phi = C=0$. The black-white dashed contour line is for the marginality condition $\gamma_{Dwl}(C, g) = 0$.

becomes negative at large $C$ as expected and decays to the small values of $C$, but will become negative only for $C < 10^{-9}$ given the chosen values of $k_{x,min}$ and $k_{y,min}$. For Interchange, one recovers that $\gamma_I(g)$ increases monotonically with $g$ becoming positive for $g \gtrsim 6 \cdot 10^{-5}$. For the turbulence simulations we select $C = C_c = 6.0 \cdot 10^{-5}$, open black square, and $g = g_c = 1.5 \cdot 10^{-4}$, closed blue triangle, for which the maximum growth rate $\gamma_{DW}(C_c)$ is equal to the maximum growth rate $\gamma_I(g_c)$, see Figure 1 right hand side. The choice that is made $g_c$ and $C_c$ is further analysed of Figure 2 left hand side. The growth rate $\gamma_{DW}$, black curve open circles, and $\gamma_I$, blue curve open triangles, are plotted versus $k_x/k_{x,min}$ for $k_x = 5k_{x,min}$. As chosen, the two growth rates have the same peak value but for different $k_y$, the Drift Wave being shifted towards lower $k_y$ compared to the Interchange case. Combining the two instabilities yields $\gamma_{Dwl}$, red curve closed squares, which close to double growth rate and intermediate peak value of $k_y$. The peak behaviour is found to compare well to $\gamma_{DW}+\gamma_I$ as if the turbulent processes where additive.

The linear stability analysis of Drift waves and Interchange instability mechanisms (Dwl) has been investigated varying the control parameters $C$ and $g$, Figure 2 right hand side where the characteristic growth rate is plotted versus $C$ in log scale lower axis and $g$ in log scale, vertical axis. As observed on Figure 1 right hand side, varying $C$ has little effect. It bounds the unstable region to the right when $C$ becomes large, typically of order $10^{-3}$. This is shown by the vertical dashed black and white line that indicates the marginality condition $\gamma_{Dwl} = 0$. For $g \to 0$, the horizontal dashed black and white line indicates the stabilising effect for small values of $C$. The unstable region governed by the drift wave instability appears to extend for $C$ in the range $10^{-6} \le C \le 10^{-3}$. The main variation of the growth rate is governed by the variation of $g$. Finally the selected conditions $C_c$ and $g_c$, such that the drift wave and interchange instability yield the same value for the characteristic growth rate is indicated by the open black marker. This point lies in a region of rather weak growth rate - $10^{-3}$ agreeing with the peak value indicated on Figure 2 left hand side.
## 4. Simulation of turbulent SOL widths

The system of two coupled equations Eq. (16) is used to investigate the effect of competing linear drives of plasma turbulence on the SOL width. Simulations are performed with the TOKAM2D code that has been fully verified [33]. The code is pseudo-spectral in both radial and poloidal directions and the time stepping scheme is an order 4 Runge-Kutta. Consequently the simulation domain is periodic in both poloidal y and radial x directions. The poloidal periodicity is consistent with an angle but the radial periodicity requires more attention. A source term is implemented with Gaussian shape in the radial direction, half width $8 \rho_0$ and constant in the poloidal direction. Consequently, the density gradient to the right of the source is negative while the density gradient to the left is positive. At maximum and minimum density, regions with weak density gradients and vanishing radial particle fluxes are expected. The latter minimum region is observed to govern a region of stability that separates the simulation domain in two parts. A negative gradient region, standing for a low field side SOL region with interchange and possibly drift wave instabilities, and a positive density gradient region, standing for a high field side SOL region with only the drift wave instability. The two SOL regions are coupled via the source region, typically the separatrix. The drawback of handling radial periodic conditions is compensated by the numerical efficiency, first the high accuracy of the scheme, but also addressing directly two plasma turbulence cases, that of drift waves with or without interchange. Since the code evolves the Fourier modes, the results of the linear analysis can be directly applied. The mesh stepping is chosen equal to $\rho_0$ in the x and y directions, therefore defining $k_{x,min} = k_{y,min} = \pi\rho_*$ where $1/\rho_*$ is the size of the simulation domain in number of mesh points $N_x$ and $N_y$. For the results presented here we have chosen $N_x = N_y = 256$ ($\rho_* \approx 10^{-3}$). One can show that the amplitude of the source defines the density normalisation and is not a relevant control parameter. For flux driven simulations we set $1/L_n = 0$ so that the available control parameters are $g, \sigma_n = \sigma_\phi$ and $\bar{\sigma}_n = \bar{\sigma}_\phi$ with fixed $D_\perp = v_\perp = 10^{-2}$. As indicated previously, we present in this paper the simulations for the three reference cases listed above, namely Drift Wave (DW), Interchange (I) and allowing for the two instability mechanisms (Dwl).

We first compare the radial profiles of the density averaged over time and poloidal angle. For the Drift Wave case, $n_{DW}$ is symmetric with respect to the mid-box location, Figure 3 left hand side. As expected, the turbulence drives a comparable transport in the low field side and high field side SOL. An exponential decay is observed close to the source with e-folding length $1/\lambda = 0.018$. This value is quite close to that used to perform the linear analysis in the previous Section, $1/L_n = 0.02$. The sharper decay towards the mid-box region is harder to explain since it suggests a reduction of the turbulent transport efficiency despite a larger departure from thermodynamic equilibrium as determined by the density gradient. Let us now analyse the profile of $n_{Dwl}$, red curve open red markers, for the simulation combining both instability mechanisms, Figure 3 right hand side. The profile is no longer symmetric and is clearly broader in the low field side where the interchange mechanism is at play. This imbalance of transport must also modify the source flux that can be expected to be larger for the low field side SOL than for the high field side SOL. Regarding the exponential fall-off, in the low field side SOL close to the source region the e-folding length is typically $1/\lambda \approx 0.011$, and further in the low field side SOL it decreases to $1/\lambda = 0.018$ before dipping to still smaller values. For the high field side SOL, the e-folding length appears to be $1/\lambda \approx 0.029$. More accurate values of the e-folding length are obtained by computing $\nabla_r n/n$ (using the Fourier transform to compute the derivative). This calculation yields directly $1/\lambda$, Figure 4 left hand side. On this figure the data from all three simulations are compared, black curve head up open black triangles for the Drift Wave only instability mechanism, blue curve head down open blue triangles for Interchange as only instability mechanism, and red curve, open red circles, when combining the two instabilities. In some regions one can easily define a nearly constant value of $\nabla_r n/n$, while this becomes less obvious in other regions, in particular for the high field side SOL where $\nabla_r n/n$ is positive. One finds that the width of the low field side SOL is broader with interchange only $1/|\lambda_I| \approx 0.013$, slightly narrower when combining the two instability mechanisms $1/|\lambda_{Dwl}| \approx 0.015$ and still
## 4. Simulation of turbulent SOL widths

The system of two coupled equations Eq. (16) is used to investigate the effect of competing linear drives of plasma turbulence on the SOL width. Simulations are performed with the TOKAM2D code that has been fully verified [33]. The code is pseudo-spectral in both radial and poloidal directions and the time stepping scheme is an order 4 Runge-Kutta. Consequently the simulation domain is periodic in both poloidal y and radial x directions. The poloidal periodicity is consistent with an angle but the radial periodicity requires more attention. A source term is implemented with Gaussian shape in the radial direction, half width $8 \rho_0$ and constant in the poloidal direction. Consequently, the density gradient to the right of the source is negative while the density gradient to the left is positive. At maximum and minimum density, regions with weak density gradients and vanishing radial particle fluxes are expected. The latter minimum region is observed to govern a region of stability that separates the simulation domain in two parts. A negative gradient region, standing for a low field side SOL region with interchange and possibly drift wave instabilities, and a positive density gradient region, standing for a high field side SOL region with only the drift wave instability. The two SOL regions are coupled via the source region, typically the separatrix. The drawback of handling radial periodic conditions is compensated by the numerical efficiency, first the high accuracy of the scheme, but also addressing directly two plasma turbulence cases, that of drift waves with or without interchange. Since the code evolves the Fourier modes, the results of the linear analysis can be directly applied. The mesh stepping is chosen equal to $\rho_0$ in the x and y directions, therefore defining $k_{x,min} = k_{y,min} = \pi\rho_*$ where $1/\rho_*$ is the size of the simulation domain in number of mesh points $N_x$ and $N_y$. For the results presented here we have chosen $N_x = N_y = 256$ ($\rho_* \approx 10^{-3}$). One can show that the amplitude of the source defines the density normalisation and is not a relevant control parameter. For flux driven simulations we set $1/L_n = 0$ so that the available control parameters are $g, \sigma_n = \sigma_\phi$ and $\bar{\sigma}_n = \bar{\sigma}_\phi$ with fixed $D_\perp = v_\perp = 10^{-2}$. As indicated previously, we present in this paper the simulations for the three reference cases listed above, namely Drift Wave (DW), Interchange (I) and allowing for the two instability mechanisms (Dwl).

We first compare the radial profiles of the density averaged over time and poloidal angle. For the Drift Wave case, $n_{DW}$ is symmetric with respect to the mid-box location, Figure 3 left hand side. As expected, the turbulence drives a comparable transport in the low field side and high field side SOL. An exponential decay is observed close to the source with e-folding length $1/\lambda = 0.018$. This value is quite close to that used to perform the linear analysis in the previous Section, $1/L_n = 0.02$. The sharper decay towards the mid-box region is harder to explain since it suggests a reduction of the turbulent transport efficiency despite a larger departure from thermodynamic equilibrium as determined by the density gradient. Let us now analyse the profile of $n_{Dwl}$, red curve open red markers, for the simulation combining both instability mechanisms, Figure 3 right hand side. The profile is no longer symmetric and is clearly broader in the low field side where the interchange mechanism is at play. This imbalance of transport must also modify the source flux that can be expected to be larger for the low field side SOL than for the high field side SOL. Regarding the exponential fall-off, in the low field side SOL close to the source region the e-folding length is typically $1/\lambda \approx 0.011$, and further in the low field side SOL it decreases to $1/\lambda = 0.018$ before dipping to still smaller values. For the high field side SOL, the e-folding length appears to be $1/\lambda \approx 0.029$. More accurate values of the e-folding length are obtained by computing $\nabla_r n/n$ (using the Fourier transform to compute the derivative). This calculation yields directly $1/\lambda$, Figure 4 left hand side. On this figure the data from all three simulations are compared, black curve head up open black triangles for the Drift Wave only instability mechanism, blue curve head down open blue triangles for Interchange as only instability mechanism, and red curve, open red circles, when combining the two instabilities. In some regions one can easily define a nearly constant value of $\nabla_r n/n$, while this becomes less obvious in other regions, in particular for the high field side SOL where $\nabla_r n/n$ is positive. One finds that the width of the low field side SOL is broader with interchange only $1/|\lambda_i| \approx 0.013$, slightly narrower when combining the two instability mechanisms $1/|\lambda_{Dwl}| \approx 0.015$ and still
<figure><img src="image_4.png" /><figcaption>Figure 3: Density profiles in log-scale, the source is located at x = 0 and x = 256 (periodicity condition).</figcaption></figure>

Left hand side, $n_{DW}$ from a simulation with drift wave instability mechanism only, black curve head up open triangles. The profile is symmetric with respect to the mid-box radius at $x = 128$. One can fit an exponential decay in the region closest to the source with e-folding length $\lambda$, $1/\lambda \approx 0.018$.

Right hand side, $n_{Dwl}$ for a simulation combining drift wave and interchange instability mechanisms, red curve open circles. We have identified 3 regions with exponential decay, two for the low field side SOL, $x \le 170$, and one for the high field side SOL, $x \ge 170$. In the low field side SOL, the decay closest to the source exhibits an e-folding length of order $1/\lambda \approx 0.011$, with a sharper decay further radially, $1/\lambda \approx 0.018$. In the high field side the exponential decay away from the source is characterised by an e-folding length $1/\lambda = 0.029$, closer to that expected without turbulent transport $1/\lambda_{diff} = 0.08$.

larger for the drift wave only case $1/\lambda_{Dwl} \approx 0.019$. One finds as previously discussed that the profiles tend to exhibit sharper gradients before the reversal of the slope, the far SOL appear narrower than the SOL close to the separatrix. Regarding the high field *side SOL*, where interchange is stable, one finds that $\nabla r/n$ approaches the value $1/\lambda_{diff} = \sigma_n/D_\perp \approx 0.08$ that would be reached without turbulent transport and only the diffusive transport. The drift wave only case is symmetric, hence $1/|\lambda_{Dwl}| \approx 0.019$. The case combining the two instabilities, appears to generate a narrower SOL with $1/|\lambda_{Dwl}| \approx 0.029$. All these results strongly disagree with the trend one could expect given the analysis of the growth rates. Some effects can be related to expected changes in phase between the density fluctuations and the fluctuations of the radial $E \times B$ drift velocity governed by the ratio $\gamma_i/\gamma_r$. Another point, in line with the importance we give to flux driven conditions, is that the flux that drives the turbulence varies in these simulations, in particular the ratio between the flux to the high field SOL and that to the low field SOL depends on the asymmetry of turbulent transport in these two regions. This could govern a non-linear effect on the SOL width. To illustrate this last aspect, we now present the so-called Sherwood dimensionless parameter for these simulations.

To investigate the transport barriers that appear in global flux driven simulations [34, 35] it is convenient to compute the ratio of the turbulent flux $\Gamma_{turb}$ divided by the total flux $\Gamma_{total} = \Gamma_{diff} + \Gamma_{turb}$ where $\Gamma_{turb} = \langle nv_{r,E \times B} \rangle$ and $\Gamma_{diff} = -D_\perp \langle \nabla_r n \rangle$. In a barrier region $\Gamma_{turb}/\Gamma_{total} \to$ while elsewhere $\Gamma_{turb}/\Gamma_{total} \to$ providing a convenient criterion to analyse the transport barrier dynamics as well as its dynamics [35]. The Sherwood dimensionless parameter $Sh$ for particle transport is similar to the Nusselt number used for heat transport, and is defined as the ratio of the convective (turbulent) flux divided by the diffusive flux, hence $Sh = \Gamma_{turb}/\Gamma_{diff}$.
<figure><img src="image_6.png" /><figcaption>Figure 3: Density profiles in log-scale, the source is located at x = 0 and x = 256 (periodicity condition).</figcaption></figure>

Left hand side, $n_{DW}$ from a simulation with drift wave instability mechanism only, black curve head up open triangles. The profile is symmetric with respect to the mid-box radius at $x = 128$. One can fit an exponential decay in the region closest to the source with e-folding length $\lambda$, $1/\lambda \approx 0.018$.

Right hand side, $n_{Dwl}$ for a simulation combining drift wave and interchange instability mechanisms, red curve open circles. We have identified 3 regions with exponential decay, two for the low field side SOL, $x \le 170$, and one for the high field side SOL, $x \ge 170$. In the low field side SOL, the decay closest to the source exhibits an e-folding length of order $1/\lambda \approx 0.011$, with a sharper decay further radially, $1/\lambda \approx 0.018$. In the high field side the exponential decay away from the source is characterised by an e-folding length $1/\lambda = 0.029$, closer to that expected without turbulent transport $1/\lambda_{diff} = 0.08$.

larger for the drift wave only case $1/\lambda_{DW} \approx 0.019$. One finds as previously discussed that the profiles tend to exhibit sharper gradients before the reversal of the slope, the far SOL appear narrower than the SOL close to the separatrix. Regarding the high field *side SOL*, where interchange is stable, one finds that $\nabla r/n$ approaches the value $1/\lambda_{diff} = \sigma_n/D_\perp \approx 0.08$ that would be reached without turbulent transport and only the diffusive transport. The drift wave only case is symmetric, hence $1/|\lambda_{DW}| \approx 0.019$. The case combining the two instabilities, appears to generate a narrower SOL with $1/\lambda_{DW} \approx 0.029$. All these results strongly disagree with the trend one could expect given the analysis of the growth rates. Some effects can be related to expected changes in phase between the density fluctuations and the fluctuations of the radial $E \times B$ drift velocity governed by the ratio $\gamma_i/\gamma_r$. Another point, in line with the importance we give to flux driven conditions, is that the flux that drives the turbulence varies in these simulations, in particular the ratio between the flux to the high field SOL and that to the low field SOL depends on the asymmetry of turbulent transport in these two regions. This could govern a non-linear effect on the SOL width. To illustrate this last aspect, we now present the so-called Sherwood dimensionless parameter for these simulations.

To investigate the transport barriers that appear in global flux driven simulations [34, 35] it is convenient to compute the ratio of the turbulent flux $\Gamma_{turb}$ divided by the total flux $\Gamma_{total} = \Gamma_{diff} + \Gamma_{turb}$ where $\Gamma_{turb} = \langle nv_{r,E \times B} \rangle$ and $\Gamma_{diff} = -D_\perp \langle \nabla_r n \rangle$. In a barrier region $\Gamma_{turb}/\Gamma_{total} \to$ while elsewhere $\Gamma_{turb}/\Gamma_{total} \to$ providing a convenient criterion to analyse the transport barrier dynamics as well as its dynamics [35]. The Sherwood dimensionless parameter $Sh$ for particle transport is similar to the Nusselt number used for heat transport, and is defined as the ratio of the convective (turbulent) flux divided by the diffusive flux, hence $Sh = \Gamma_{turb}/\Gamma_{diff}$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Left hand side: Profiles of $\nabla n/n$ for the simulation with drift wave instability mechanism only, black curve head-up open triangles, interchange only blue curve head-down open blue triangles and the case combining the two instabilities, red curve open red circles. Regions with near constant values of $\nabla r/n$, indicated by the horizontal dashed lines allow us to define the e-folding lengths discussed in the text. The top horizontal black dashed line is the value of $|\nabla r/n|$ that is expected for diffusive transport without turbulent $E \times B$ convection. Right hand side: Profiles of the Sherwood parameter for the same simulations and same convention, DW black curve with open head-up triangles, I blue curve with blue open head-down triangles, Dwl red curve with red open circles. Here the Sherwood number is computed with the time and y averages of the fluxes. The Sherwood number diverges where $\nabla \times n = 0$, in the source region and towards the mid-box where the density gradient reverses, the latter locations are indicated by the vertical dash-dot lines.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Sherwood number $Sh$ computed with y-averaged fluxes versus radial position $x$, horizontal axis, and time $t$, vertical axis. Data from the simulation "I" with only the interchange instability mechanism.
Left hand side: Plot for the time interval $\Omega_0 \Delta t = 6 \cdot 10^4$, corresponding to 1250 snapshots, Right hand side: Plot for a reduced time window $\Omega_0 \Delta t = 2 \cdot 10^4$ and different colour-scale to highlight the structure with negative values of $Sh$ and the region with $Sh \approx 1$, green colour range. The black plain line is oriented like the typical ballistic propagation of the outward avalanches. The vertical mauve dash-dot lines indicate the radial position where $\langle \Gamma_{\text{diff}} \rangle_{y,t} = 0$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Left hand side: Profiles of $\nabla n/n$ for the simulation with drift wave instability mechanism only, black curve head-up open triangles, interchange only blue curve head-down open blue triangles and the case combining the two instabilities, red curve open red circles. Regions with near constant values of $\nabla r/n$, indicated by the horizontal dashed lines allow us to define the e-folding lengths discussed in the text. The top horizontal black dashed line is the value of $|\nabla r/n|$ that is expected for diffusive transport without turbulent $E \times B$ convection. Right hand side: Profiles of the Sherwood parameter for the same simulations and same convention, DW black curve with open head-up triangles, I blue curve with blue open head-down triangles, Dwl red curve with red open circles. Here the Sherwood number is computed with the time and y averages of the fluxes. The Sherwood number diverges where $\nabla \times n = 0$, in the source region and towards the mid-box where the density gradient reverses, the latter locations are indicated by the vertical dash-dot lines.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Sherwood number $Sh$ computed with y-averaged fluxes versus radial position $x$, horizontal axis, and time $t$, vertical axis. Data from the simulation "I" with only the interchange instability mechanism.
Left hand side: Plot for the time interval $\Omega_0 \Delta t = 6 \cdot 10^4$, corresponding to 1250 snapshots, Right hand side: Plot for a reduced time window $\Omega_0 \Delta t = 2 \cdot 10^4$ and different colour-scale to highlight the structure with negative values of $Sh$ and the region with $Sh \approx 1$, green colour range. The black plain line is oriented like the typical ballistic propagation of the outward avalanches. The vertical mauve dash-dot lines indicate the radial position where $\langle \Gamma_{\text{diff}} \rangle_{y,t} = 0$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Sherwood number $Sh$ computed with y-averaged fluxes versus radial position $x$, horizontal axis, and time $t$, vertical axis. The black plain line is oriented like the typical ballistic propagation of the outward avalanches. The vertical mauve dash-dot lines indicate the radial position where $\langle \Gamma_{diff} \rangle_{y,t} = 0$. Left hand side: Simulation "Dwl" of Drift waves with interchange. Right hand side: simulation "DW" with drift wave only.</figcaption></figure>

It is clearly related to the criterion we have used to analyse the transport barriers and can also be related to the SOL width: $Sh = \lambda^2/\lambda_{diff}^2 - 1$. Depending on the choice made for the fluxes three different values are obtained. Local values of the fluxes depending on radial position $x$, poloidal angle $y$ and time can be used. Very large fluctuations in $x$, $y$ and time are then observed. Regarding the effective radial transport, it is more convenient to consider y-averaged fluxes yielding a time dependent Sherwood profile. This reduces significantly the magnitude of the Sherwood number and retains the dynamics. Finally, one can use the time and y-averaged fluxes yielding a radial profile, with a further reduction in the magnitude of the Sherwood number. The profiles of the Sherwood dimensionless number are plotted on Figure 4 right hand side. This Figure summarizes many points discussed above. First, $Sh$ tends to decrease gradually for $x$ & 64 towards the point with inversion of the density gradient reversal (indicated by the vertical dashed lines). Second, the interchange exhibits a higher value of $Sh$ compared to the case combining the two instabilities. Third, the drift wave case is identical for the high and low field side SOL, when the interchange is the only linear drive it does not contribute to the low field side SOL turbulence and for the high field side SOL the simulation with combined instabilities exhibits less turbulent transport than that with drift wave only.

On Figure 5 are presented the Sherwood number determined with for y-averaged (flux surface average in this 2D model) radial fluxes of particles for the simulation "l" with interchange as only linear drive. On Figure 5 left hand side, the colour-scale is chosen to highlight the avalanches. These exhibit long range propagation from 64 $\le x \le$ 192. For $0 \le x \le$ 64, closer to the source location at $x=0$, corrugations [17, 16], the vertically oriented dark structures, appear to inhibit the long range propagation. The colour scale is changed on Figure 5 right hand side for comparison to the other simulations with drift waves. The time window is reduced and colour-scale symmetric between -250 and +250. The colour-scale also allows one to identify the region with positive Sherwood number and $Sh \approx 1$ that corresponds to the green colour-range. One can then clearly notice that the long range avalanche region extends up to $x \approx 194$ highlighted by the vertical mauve dash dot line that indicates the radial position where the time average density profile reverse, therefore with zero gradient and vanishing diffusive flux $\Gamma_{diff}$. One can remark strongly negative structures of the Sherwood numbers, blue colour-scale, embedded in regions of maximum Sherwood number, with white colour-scale, and located in the region that exhibits the corrugations. The same plot, same time window and same colour-scale, is shown on Figure 6 for the simulation "Dwl" with drift wave and interchange, left hand side, and simulation "DW" with only drift wave. On the left hand side with interchange present,
<figure><img src="image_9.png" /><figcaption>Figure 6: Sherwood number $Sh$ computed with y-averaged fluxes versus radial position $x$, horizontal axis, and time $t$, vertical axis. The black plain line is oriented like the typical ballistic propagation of the outward avalanches. The vertical mauve dash-dot lines indicate the radial position where $\langle \Gamma_{diff} \rangle_{y,t} = 0$. Left hand side: Simulation "Dwl" of Drift waves with interchange. Right hand side: simulation "DW" with drift wave only.</figcaption></figure>

It is clearly related to the criterion we have used to analyse the transport barriers and can also be related to the SOL width: $Sh = \lambda^2/\lambda_{diff}^2 - 1$. Depending on the choice made for the fluxes three different values are obtained. Local values of the fluxes depending on radial position $x$, poloidal angle $y$ and time can be used. Very large fluctuations in $x$, $y$ and time are then observed. Regarding the effective radial transport, it is more convenient to consider y-averaged fluxes yielding a time dependent Sherwood profile. This reduces significantly the magnitude of the Sherwood number and retains the dynamics. Finally, one can use the time and y-averaged fluxes yielding a radial profile, with a further reduction in the magnitude of the Sherwood number. The profiles of the Sherwood dimensionless number are plotted on Figure 4 right hand side. This Figure summarizes many points discussed above. First, $Sh$ tends to decrease gradually for $x$ & 64 towards the point with inversion of the density gradient reversal (indicated by the vertical dashed lines). Second, the interchange exhibits a higher value of $Sh$ compared to the case combining the two instabilities. Third, the drift wave case is identical for the high and low field side SOL, when the interchange is the only linear drive it does not contribute to the low field side SOL turbulence and for the high field side SOL the simulation with combined instabilities exhibits less turbulent transport than that with drift wave only.

On Figure 5 are presented the Sherwood number determined with for y-averaged (flux surface average in this 2D model) radial fluxes of particles for the simulation "I" with interchange as only linear drive. On Figure 5 left hand side, the colour-scale is chosen to highlight the avalanches. These exhibit long range propagation from 64 $x \le 192$. For $0 \le x \le 64$, closer to the source location at $x=0$, corrugations [17, 16], the vertically oriented dark structures, appear to inhibit the long range propagation. The colour scale is changed on Figure 5 right hand side for comparison to the other simulations with drift waves. The time window is reduced and colour-scale symmetric between -250 and +250. The colour-scale also allows one to identify the region with positive Sherwood number and $Sh \approx 1$ that corresponds to the green colour-range. One can then clearly notice that the long range avalanche region extends up to $x \approx 194$ highlighted by the vertical mauve dash dot line that indicates the radial position where the time average density profile reverse, therefore with zero gradient and vanishing diffusive flux $\Gamma_{diff}$. One can remark strongly negative structures of the Sherwood numbers, blue colour-scale, embedded in regions of maximum Sherwood number, with white colour-scale, and located in the region that exhibits the corrugations. The same plot, same time window and same colour-scale, is shown on Figure 6 for the simulation "Dwl" with drift wave and interchange, left hand side, and simulation "DW" with only drift wave. On the left hand side with interchange present,
the pattern of avalanche transport is still apparent but inhibited compared to the case with
only interchange. Corrugation patterns are not visible. Regions with large negative Sherwood
number are more frequent and appear to align according to the avalanche pattern. This suggests
time sequences with large avalanches with outward flux, these being nearly concomitant with
avalanches with large inward fluxes. The reversal of the avalanche direction of transport can
be related to the phase change between the over-dense structures and the location of maximum
radial E × B convection. Finally, in the case with drift wave only turbulence, the avalanche
pattern is clearly inhibited and the magnitude of the Sherwood number is reduced for both the
negative and positive values. The ballistic propagation of the outward avalanches indicated by
the black lines, Figure 5 right hand side and Figure 6 figure: sherwood DW left and right hand
side are found to be quite comparable, 4.6 ρ\*c_s, 4.4 ρ\*c_s and 5.1 ρ\*c_s respectively.

**5. Discussion and conclusion**

In this paper we have revisited the seminal papers on SOL interchange turbulence [10] and
drift wave turbulence [25] to obtain a versatile 2D model to investigate numerically plasma
turbulent transport. We have shown that the instability triggered by the effect of temperature
fluctuations on the sheath conductivity, when linearised, is in fact comparable to drift wave
turbulence. Therefore, provided one extends the model to allow the transported field to be
either the temperature at constant density or the density at constant temperature, one can use
this 2D model with two different instability mechanisms to investigate both the edge plasma
and the SOL. The model then describes the E × B convection of the transported field, here
defined as the density, to the E × B convection of the vorticity, identical in this simplified
version to the Laplacian of the electric potential φ. The physics of interface barriers has been
investigated with such a model [23] modified to handle both closed and open field lines. The
linear analysis indicates that most of the parallel losses participate to damping processes of
the large-scale fluctuations -while dissipative effects are damping processes of the small-scale
fluctuations. Conversely, the buoyancy effect with control parameter *g* and the parallel the
parallel conductivity proportional to the density field with control parameter *C*, govern the
coupling between the two fields and consequently the linear instabilities. The interchange
like instability is shown to depend on the sign of the density gradient making it stable in the
high field side SOL and possibly unstable in the low field side SOL. conversely the drift wave
turbulence does not depend on the sign of the density gradient and can be triggered in both
high and low field SOL. Furthermore, the interchange instability grows with *g* but the drift
wave instability levels-off at modest growth rate and is stabilised at large values of *C*. Using
the linear analysis we have selected a case, such that the growth rate of drift waves and the
interchange growth rate are identical, the growth rate when both instabilities are at play being
roughly the double. The choice defines three cases that have been investigated numerically
for both high and low field SOL. We have used the SOL width to characterise the turbulent
transport. The largest SOL e-folding length is observed with interchange only, and therefore
is restricted to the low field side turbulent transport. The model including both interchange
and drift waves is found to yield a small but comparable width on the low field side, together
with the narrowest SOL width located on the high field side. Finally, the drift wave case
with no buoyancy effect generates the same narrow SOL width in the high and low field SOL.
The Sherwood dimensionless number, ratio of the convective particle transport, therefore the
turbulent particle transport, divided by the diffusive transport, is also used to characterise the
turbulent transport. It indicates, agreeing with the analysis of the density gradient, that for the
chosen parameters the SOL is rather inhomogeneous. The largest e-folding length and turbulent
transport are localised towards the source region (the separatrix), then the amplitude of the
gradients increase significantly while the effectiveness of turbulent transport is reduced. The
ongoing investigation suggests that turbulence with corrugations characterise the broad SOL
region, while long range avalanche transport appear to prevail where the turbulent transport
is reduced and gradient amplitude increase. The connection to the issue of the distance to
marginality, which is most likely equivalent to a weak flux forcing, is presently analysed.
the pattern of avalanche transport is still apparent but inhibited compared to the case with
only interchange. Corrugation patterns are not visible. Regions with large negative Sherwood
number are more frequent and appear to align according to the avalanche pattern. This suggests
time sequences with large avalanches with outward flux, these being nearly concomitant with
avalanches with large inward fluxes. The reversal of the avalanche direction of transport can
be related to the phase change between the over-dense structures and the location of maximum
radial E × B convection. Finally, in the case with drift wave only turbulence, the avalanche
pattern is clearly inhibited and the magnitude of the Sherwood number is reduced for both the
negative and positive values. The ballistic propagation of the outward avalanches indicated by
the black lines, Figure 5 right hand side and Figure 6 figure: sherwood DW left and right hand
side are found to be quite comparable, 4.6 ρ\*c_s, 4.4 ρ\*c_s and 5.1 ρ\*c_s respectively.

**5. Discussion and conclusion**

In this paper we have revisited the seminal papers on SOL interchange turbulence [10] and
drift wave turbulence [25] to obtain a versatile 2D model to investigate numerically plasma
turbulent transport. We have shown that the instability triggered by the effect of temperature
fluctuations on the sheath conductivity, when linearised, is in fact comparable to drift wave
turbulence. Therefore, provided one extends the model to allow the transported field to be
either the temperature at constant density or the density at constant temperature, one can use
this 2D model with two different instability mechanisms to investigate both the edge plasma
and the SOL. The model then describes the E × B convection of the transported field, here
defined as the density, to the E × B convection of the vorticity, identical in this simplified
version to the Laplacian of the electric potential ϕ. The physics of interface barriers has been
investigated with such a model [23] modified to handle both closed and open field lines. The
linear analysis indicates that most of the parallel losses participate to damping processes of
the large-scale fluctuations -while dissipative effects are damping processes of the small-scale
fluctuations. Conversely, the buoyancy effect with control parameter *g* and the parallel the
parallel conductivity proportional to the density field with control parameter *C*, govern the
coupling between the two fields and consequently the linear instabilities. The interchange
like instability is shown to depend on the sign of the density gradient making it stable in the
high field side SOL and possibly unstable in the low field side SOL. conversely the drift wave
turbulence does not depend on the sign of the density gradient and can be triggered in both
high and low field SOL. Furthermore, the interchange instability grows with *g* but the drift
wave instability levels-off at modest growth rate and is stabilised at large values of *C*. Using
the linear analysis we have selected a case, such that the growth rate of drift waves and the
interchange growth rate are identical, the growth rate when both instabilities are at play being
roughly the double. The choice defines three cases that have been investigated numerically
for both high and low field SOL. We have used the SOL width to characterise the turbulent
transport. The largest SOL e-folding length is observed with interchange only, and therefore
is restricted to the low field side turbulent transport. The model including both interchange
and drift waves is found to yield a small but comparable width on the low field side, together
with the narrowest SOL width located on the high field side. Finally, the drift wave case
with no buoyancy effect generates the same narrow SOL width in the high and low field SOL.
The Sherwood dimensionless number, ratio of the convective particle transport, therefore the
turbulent particle transport, divided by the diffusive transport, is also used to characterise the
turbulent transport. It indicates, agreeing with the analysis of the density gradient, that for the
chosen parameters the SOL is rather inhomogeneous. The largest e-folding length and turbulent
transport are localised towards the source region (the separatrix), then the amplitude of the
gradients increase significantly while the effectiveness of turbulent transport is reduced. The
ongoing investigation suggests that turbulence with corrugations characterise the broad SOL
region, while long range avalanche transport appear to prevail where the turbulent transport
is reduced and gradient amplitude increase. The connection to the issue of the distance to
marginality, which is most likely equivalent to a weak flux forcing, is presently analysed.
## Acknowledgements

This work has been carried out thanks to the support of the A*MIDEX project (ANR-11-IDEX-0001 02, TOP project) funded by the ‘Investissements d’Avenir’ French Government program, managed by the French National Research Agency (ANR). This work has been also supported by the French National Research Agency grant SISTEM (ANR-19-CE46-0005-03) and by the French National Research Agency grant AIM4EP (ANR-21-CE30-0018-01). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 633053 and Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

## References

### References

[1] ITER Physics Expert Group 1999 *Nuclear Fusion* **39** 2175-2249 URL https://iopscience.iop.org/article/10.1088/0029-5515/39/12/302

[2] Loarte A, Lipschultz B, Kukushkin A, Matthews G, Stangeby P, Asakura N, Counsell G, Federici G, Kallenbach A, Krieger K, Mahdavi A, Philipps V, Reiter D, Roth J, Strachan J, Whyte D, Doerner R, Eich T, Fundamenski W, Herrmann A, Fenstermacher M, Ghendrih P, Groth M, Kirschner A, Konoshima S, LaBombard B, Lang P, Leonard A, Monier-Garbet P, Neu R, Pacher H, Pegourie B, Pitts R, Takamura S, Terry J, Tsitrone E, the ITPA Scrape-off Layer and Group D 2007 *Nuclear Fusion* **47** S203-S263 URL https://doi.org/10.1088/2F0029-5515/2F47/2F6/2Fs04

[3] Garbet X, Idomura Y, Villard L and Watanabe T 2010 *Nuclear Fusion* **50** 043002 URL http://stacks.iop.org/0029-5515/50/i=4/a=043002

[4] Grandgirard V, Abiteboul J, Bigot J, Cartier-Michaud T, Crouseilles N, Dif-Pradalier G, Ehrlacher C, Esteve D, Garbet X, Ghendrih P, Latu G, Mehrenberger M, Norscini C, Passeron C, Rozar F, Sarazin Y, Sonnendrucker E, Strugarek A and Zarzoso D 2016 *Computer Physics Communications* **207** 35-68 ISSN 0010-4655 URL http://www.sciencedirect.com/science/article/pii/S0010465516301230

[5] Tamain P, Bufferand H, Ciraolo G, Colin C, Galassi D, Ghendrih P, Schwander F and Serre E 2016 *Journal of Computational Physics* **321** 606-623 ISSN 0021-9991 URL http://www.sciencedirect.com/science/article/pii/S0021999116301838

[6] Bufferand H, Bucalossi J, Ciraolo G, Falchetto G, Gallo A, Ghendrih P, Rivals N, Tamain P, Yang H, Giorgiani G, Schwander F, d'Abusco M S, Serre E, Marandet Y and Raghunathan M 2021 *Nuclear Fusion* **61** 116052 URL https://doi.org/10.1088/1741-4326/ac2873

[7] Bourdelle C, Garbet X, Imbeaux F, Casati A, Dubuit N, Guirlet R and Parisot T 2007 *Physics of Plasmas* **14** 112501 (Preprint https://doi.org/10.1063/1.2800869) URL https://doi.org/10.1063/1.2800869

[8] Bufferand H, Bensiali B, Bucalossi J, Ciraolo G, Genesio P, Ghendrih P, Marandet Y, Paredes A, Schwander F, Serre E and Tamain P 2013 *Journal of Nuclear Materials* **438** S445-S448 ISSN 0022-3115 URL http://www.sciencedirect.com/science/article/pii/S0022311513000986

[9] Baschetti S, Bufferand H, Ciraolo G, Ghendrih P, Serre E, Tamain P and the WEST Team 2021 *Nuclear Fusion* **61** 106020 URL https://doi.org/10.1088/1741-4326/ac1e60

[10] Sarazin Y and Ghendrih P 1998 *Physics of Plasmas* **5** 4214-4228 (Preprint https://doi.org/10.1063/1.873157) URL https://doi.org/10.1063/1.873157
## Acknowledgements

This work has been carried out thanks to the support of the A*MIDEX project (ANR-11-IDEX-0001 02, TOP project) funded by the ‘Investissements d’Avenir’ French Government program, managed by the French National Research Agency (ANR). This work has been also supported by the French National Research Agency grant SISTEM (ANR-19-CE46-0005-03) and by the French National Research Agency grant AIM4EP (ANR-21-CE30-0018-01). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 633053 and Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

## References

### References

[1] ITER Physics Expert Group 1999 *Nuclear Fusion* **39** 2175-2249 URL https://iopscience.iop.org/article/10.1088/0029-5515/39/12/302

[2] Loarte A, Lipschultz B, Kukushkin A, Matthews G, Stangeby P, Asakura N, Counsell G, Federici G, Kallenbach A, Krieger K, Mahdavi A, Philipps V, Reiter D, Roth J, Strachan J, Whyte D, Doerner R, Eich T, Fundamenski W, Herrmann A, Fenstermacher M, Ghendrih P, Groth M, Kirschner A, Konoshima S, LaBombard B, Lang P, Leonard A, Monier-Garbet P, Neu R, Pacher H, Pegourie B, Pitts R, Takamura S, Terry J, Tsitrone E, the ITPA Scrape-off Layer and Group D 2007 *Nuclear Fusion* **47** S203-S263 URL https://doi.org/10.1088/2F0029-5515/2F47/2F6/2Fs04

[3] Garbet X, Idomura Y, Villard L and Watanabe T 2010 *Nuclear Fusion* **50** 043002 URL http://stacks.iop.org/0029-5515/50/i=4/a=043002

[4] Grandgirard V, Abiteboul J, Bigot J, Cartier-Michaud T, Crouseilles N, Dif-Pradalier G, Ehrlacher C, Esteve D, Garbet X, Ghendrih P, Latu G, Mehrenberger M, Norscini C, Passeron C, Rozar F, Sarazin Y, Sonnendrucker E, Strugarek A and Zarzoso D 2016 *Computer Physics Communications* **207** 35-68 ISSN 0010-4655 URL http://www.sciencedirect.com/science/article/pii/S0010465516301230

[5] Tamain P, Bufferand H, Ciraolo G, Colin C, Galassi D, Ghendrih P, Schwander F and Serre E 2016 *Journal of Computational Physics* **321** 606-623 ISSN 0021-9991 URL http://www.sciencedirect.com/science/article/pii/S0021999116301838

[6] Bufferand H, Bucalossi J, Ciraolo G, Falchetto G, Gallo A, Ghendrih P, Rivals N, Tamain P, Yang H, Giorgiani G, Schwander F, d'Abusco M S, Serre E, Marandet Y and Raghunathan M 2021 *Nuclear Fusion* **61** 116052 URL https://doi.org/10.1088/1741-4326/ac2873

[7] Bourdelle C, Garbet X, Imbeaux F, Casati A, Dubuit N, Guirlet R and Parisot T 2007 *Physics of Plasmas* **14** 112501 (Preprint https://doi.org/10.1063/1.2800869) URL https://doi.org/10.1063/1.2800869

[8] Bufferand H, Bensiali B, Bucalossi J, Ciraolo G, Genesio P, Ghendrih P, Marandet Y, Paredes A, Schwander F, Serre E and Tamain P 2013 *Journal of Nuclear Materials* **438** S445-S448 ISSN 0022-3115 URL http://www.sciencedirect.com/science/article/pii/S0022311513000986

[9] Baschetti S, Bufferand H, Ciraolo G, Ghendrih P, Serre E, Tamain P and the WEST Team 2021 *Nuclear Fusion* **61** 106020 URL https://doi.org/10.1088/1741-4326/ac1e60

[10] Sarazin Y and Ghendrih P 1998 *Physics of Plasmas* **5** 4214-4228 (Preprint https://doi.org/10.1063/1.873157) URL https://doi.org/10.1063/1.873157
[11] Ghendrih P, Sarazin Y, Attuel G, Benkadda S, Beyer P, Falchetto G, Figarella C, Garbet X, Grandgirard V and Ottaviani M 2003 *Nuclear Fusion* **43** 1013 URL http://stacks.iop.org/0029-5515/43/i=10/a=001

[12] Fedorczak N, Gallo A, Tamain P, Bufferand H, Ciraolo G and Ghendrih P 2018 *Contributions to Plasma Physics* **58** 471–477 (Preprint https://www.onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201700169) URL https://www.onlinelibrary.wiley.com/doi/abs/10.1002/ctpp.201700169

[13] Sarazin Y, Garbet X, Ghendrih P and Benkadda S 2000 *Physics of Plasmas* (1994-present) **7** 1085–1088 URL http://scitation.aip.org/content/aip/journal/pop/7/4/10.1063/1.873947

[14] Krasheninnikov S I 2001 *Physics Letters A* **283** 368–370 ISSN 0375-9601 URL https://www.sciencedirect.com/science/article/pii/S0375960101002523

[15] Sarazin Y, Grandgirard V, Abiteboul J, Allfrey S, Garbet X, Ghendrih P, Latu G, Strugarek A and Dif-Pradalier G 2010 *Nuclear Fusion* **50** 054004 URL http://stacks.iop.org/0029-5515/50/i=5/a=054004

[16] Ghendrih P, Asahi Y, Caschera E, Dif-Pradalier G, Donnel P, Garbet X, Gillot C, Grandgirard V, Latu G, Sarazin Y, Baschetti S, Bufferand H, Cartier-Michaud T, Ciraolo G, Tamain P, Tatali R and Serre E 2018 *Journal of Physics: Conference Series* **1125** 012011 URL http://stacks.iop.org/1742-6596/1125/i=1/a=012011

[17] Dif-Pradalier G, Hornung G, Ghendrih P, Sarazin Y, Clairet F, Vermare L, Diamond P H, Abiteboul J, Cartier-Michaud T, Ehrlacher C, Estève D, Garbet X, Grandgirard V, Gürcan O D, Hennequin P, Kosuga Y, Latu G, Maget P, Morel P, Norscini C, Sabot R and Storelli A 2015 *Phys. Rev. Lett.* **114**(8) 085004 URL https://link.aps.org/doi/10.1103/PhysRevLett.114.085004

[18] Stangeby P 2000 *The Plasma Boundary of Magnetic Fusion Devices* Series in Plasma Physics (IOP Publishing Ltd.)

[19] Bourne E, Munschy Y, Grandgirard V, Mehrenberger M and Ghendrih P 2022 *Non-Uniform Splines for Semi-Lagrangian Kinetic Simulations of the Plasma Sheath* working paper or preprint URL https://hal-cea.archives-ouvertes.fr/cea-03748016

[20] Sarazin Y 1997 *Etude de la turbulence de bord dans les plasmas de tokamaks* Ph.D. thesis Université de Grenoble

[21] Garbet X, Laurent L, Roubin J P and Samain A 1991 *Nuclear Fusion* **31** 967 URL http://stacks.iop.org/0029-5515/31/i=5/a=015

[22] Nedospasov A V 1993 *Physics of Fluids B: Plasma Physics* **5** 3191–3194 (Preprint https://doi.org/10.1063/1.860655) URL https://doi.org/10.1063/1.860655

[23] Norscini C, Cartier-Michaud T, Dif-Pradalier G, Garbet X, Ghendrih P, Grandgirard V and Sarazin Y 2022 *Plasma Physics and Controlled Fusion* **64** 055007 URL https://doi.org/10.1088/1361-6587/ac5a07

[24] Diamond P H et al. 2005 *Plasma Physics and Controlled Fusion* **47** R35–R161 ISSN 1361-6587 URL http://dx.doi.org/10.1088/0741-3335/47/5/R01

[25] Hasegawa A and Wakatani M 1983 *Phys. Rev. Lett.* **50**(9) 682–686 URL http://link.aps.org/doi/10.1103/PhysRevLett.50.682

[26] Berk H L, Ryutov D D and Tsidulko Y A 1991 *Physics of Fluids B: Plasma Physics* **3** 1346–1354 (Preprint https://doi.org/10.1063/1.859699) URL https://doi.org/10.1063/1.859699

[27] Normand C, Pomeau Y and Velarde M G 1977 *Review of Modern Physics* **49**(3) 581–624 URL https://link.aps.org/doi/10.1103/RevModPhys.49.581

[28] Wilczynski F, Hughes D W, Van Loo S, Arter W and Militello F 2019 *Physics of Plasmas* **26** 022510 (Preprint https://doi.org/10.1063/1.5064765) URL https://doi.org/10.1063/1.5064765
[11] Ghendrih P, Sarazin Y, Attuel G, Benkadda S, Beyer P, Falchetto G, Figarella C, Garbet X, Grandgirard V and Ottaviani M 2003 *Nuclear Fusion* **43** 1013 URL http://stacks.iop.org/0029-5515/43/i=10/a=001

[12] Fedorczak N, Gallo A, Tamain P, Bufferand H, Ciraolo G and Ghendrih P 2018 *Contributions to Plasma Physics* **58** 471–477 (Preprint https://www.onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201700169) URL https://www.onlinelibrary.wiley.com/doi/abs/10.1002/ctpp.201700169

[13] Sarazin Y, Garbet X, Ghendrih P and Benkadda S 2000 *Physics of Plasmas* (1994-present) **7** 1085–1088 URL http://scitation.aip.org/content/aip/journal/pop/7/4/10.1063/1.873947

[14] Krasheninnikov S I 2001 *Physics Letters A* **283** 368–370 ISSN 0375-9601 URL https://www.sciencedirect.com/science/article/pii/S0375960101002523

[15] Sarazin Y, Grandgirard V, Abiteboul J, Allfrey S, Garbet X, Ghendrih P, Latu G, Strugarek A and Dif-Pradalier G 2010 *Nuclear Fusion* **50** 054004 URL http://stacks.iop.org/0029-5515/50/i=5/a=054004

[16] Ghendrih P, Asahi Y, Caschera E, Dif-Pradalier G, Donnel P, Garbet X, Gillot C, Grandgirard V, Latu G, Sarazin Y, Baschetti S, Bufferand H, Cartier-Michaud T, Ciraolo G, Tamain P, Tatali R and Serre E 2018 *Journal of Physics: Conference Series* **1125** 012011 URL http://stacks.iop.org/1742-6596/1125/i=1/a=012011

[17] Dif-Pradalier G, Hornung G, Ghendrih P, Sarazin Y, Clairet F, Vermare L, Diamond P H, Abiteboul J, Cartier-Michaud T, Ehrlacher C, Estève D, Garbet X, Grandgirard V, Gürcan O D, Hennequin P, Kosuga Y, Latu G, Maget P, Morel P, Norscini C, Sabot R and Storelli A 2015 *Phys. Rev. Lett.* **114**(8) 085004 URL https://link.aps.org/doi/10.1103/PhysRevLett.114.085004

[18] Stangeby P 2000 *The Plasma Boundary of Magnetic Fusion Devices* Series in Plasma Physics (IOP Publishing Ltd.)

[19] Bourne E, Munschy Y, Grandgirard V, Mehrenberger M and Ghendrih P 2022 *Non-Uniform Splines for Semi-Lagrangian Kinetic Simulations of the Plasma Sheath* working paper or preprint URL https://hal-cea.archives-ouvertes.fr/cea-03748016

[20] Sarazin Y 1997 *Etude de la turbulence de bord dans les plasmas de tokamaks* Ph.D. thesis Université de Grenoble

[21] Garbet X, Laurent L, Roubin J P and Samain A 1991 *Nuclear Fusion* **31** 967 URL http://stacks.iop.org/0029-5515/31/i=5/a=015

[22] Nedospasov A V 1993 *Physics of Fluids B: Plasma Physics* **5** 3191–3194 (Preprint https://doi.org/10.1063/1.860655) URL https://doi.org/10.1063/1.860655

[23] Norscini C, Cartier-Michaud T, Dif-Pradalier G, Garbet X, Ghendrih P, Grandgirard V and Sarazin Y 2022 *Plasma Physics and Controlled Fusion* **64** 055007 URL https://doi.org/10.1088/1361-6587/ac5a07

[24] Diamond P H et al. 2005 *Plasma Physics and Controlled Fusion* **47** R35–R161 ISSN 1361-6587 URL http://dx.doi.org/10.1088/0741-3335/47/5/R01

[25] Hasegawa A and Wakatani M 1983 *Phys. Rev. Lett.* **50**(9) 682–686 URL http://link.aps.org/doi/10.1103/PhysRevLett.50.682

[26] Berk H L, Ryutov D D and Tsidulko Y A 1991 *Physics of Fluids B: Plasma Physics* **3** 1346–1354 (Preprint https://doi.org/10.1063/1.859699) URL https://doi.org/10.1063/1.859699

[27] Normand C, Pomeau Y and Velarde M G 1977 *Review of Modern Physics* **49**(3) 581–624 URL https://link.aps.org/doi/10.1103/RevModPhys.49.581

[28] Wilczynski F, Hughes D W, Van Loo S, Arter W and Militello F 2019 *Physics of Plasmas* **26** 022510 (Preprint https://doi.org/10.1063/1.5064765) URL https://doi.org/10.1063/1.5064765
[29] LaBombard B, Rice J, Hubbard A, Hughes J, Greenwald M, Irby J, Lin Y, Lipschultz B, Marmar E, Pitcher C, Smick N, Wolfe S, Wukitch S and the Alcator Group 2004 *Nuclear Fusion* **44** 1047–1066 URL https://doi.org/10.1088/0029-5515/44/10/001

[30] Asakura N 2007 *Journal of Nuclear Materials* **363-365** 41–51 ISSN 0022-3115 plasma-Surface Interactions-17 URL https://www.sciencedirect.com/science/article/pii/S0022311506006325

[31] Gunn J, Boucher C, Dionne M, Ďuran I, Fuchs V, Loarer T, Nanobashvili I, P'anek R, Pascal J Y, Saint-Laurent F, Stöckel J, Van Rompuy T, Zag'orski R, Ad'amek J, Bucalossi J, Dejarnac R, Devynck P, Hertout P, Hron M, Lebrun G, Moreau P, Rimini F, Sarkissian A and Van Oost G 2007 *Journal of Nuclear Materials* **363-365** 484 – 490 ISSN 0022-3115 plasma-Surface Interactions-17 URL http://www.sciencedirect.com/science/article/pii/S0022311507000943

[32] Dif-Pradalier G, Gunn J, Ciraolo G, Chang C, Chiavassa G, Diamond P, Fedorczak N, Ghendrih P, Isoardi L, Kocan M, Ku S, Serre E and Tamain P 2011 *Journal of Nuclear Materials* **415** S597–S600 ISSN 0022-3115 proceedings of the 19th International Conference on Plasma-Surface Interactions in Controlled Fusion URL https://www.sciencedirect.com/science/article/pii/S0022311510008573

[33] Cartier-Michaud T, Ghendrih P, Sarazin Y, Abiteboul J, Bufferand H, Dif-Pradalier G, Garbet X, Grandgirard V, Latu G, Norscini C, Passeron C and Tamain P 2016 *Physics of Plasmas* **23** 020702 URL http://scitation.aip.org/content/aip/journal/pop/23/2/10.1063/1.4941974

[34] Beyer P, Benkadda S, Fuhr-Chaudier G, Garbet X, Ghendrih P and Sarazin Y 2007 *Plasma Physics and Controlled Fusion* **49** 507 URL http://stacks.iop.org/0741-3335/49/i=4/a=013

[35] Floriani E, Ciraolo G, Ghendrih P, Lima R and Sarazin Y 2013 *Plasma Physics and Controlled Fusion* **55** 095012 URL https://doi.org/10.1088/0741-3335/55/9/095012
[29] LaBombard B, Rice J, Hubbard A, Hughes J, Greenwald M, Irby J, Lin Y, Lipschultz B, Marmar E, Pitcher C, Smick N, Wolfe S, Wukitch S and the Alcator Group 2004 *Nuclear Fusion* **44** 1047–1066 URL https://doi.org/10.1088/0029-5515/44/10/001

[30] Asakura N 2007 *Journal of Nuclear Materials* **363-365** 41–51 ISSN 0022-3115 plasma-Surface Interactions-17 URL https://www.sciencedirect.com/science/article/pii/S0022311506006325

[31] Gunn J, Boucher C, Dionne M, Ďuran I, Fuchs V, Loarer T, Nanobashvili I, P'anek R, Pascal J Y, Saint-Laurent F, Stöckel J, Van Rompuy T, Zag'orski R, Ad'amek J, Bucalossi J, Dejarnac R, Devynck P, Hertout P, Hron M, Lebrun G, Moreau P, Rimini F, Sarkissian A and Van Oost G 2007 *Journal of Nuclear Materials* **363-365** 484 – 490 ISSN 0022-3115 plasma-Surface Interactions-17 URL http://www.sciencedirect.com/science/article/pii/S0022311507000943

[32] Dif-Pradalier G, Gunn J, Ciraolo G, Chang C, Chiavassa G, Diamond P, Fedorczak N, Ghendrih P, Isoardi L, Kocan M, Ku S, Serre E and Tamain P 2011 *Journal of Nuclear Materials* **415** S597–S600 ISSN 0022-3115 proceedings of the 19th International Conference on Plasma-Surface Interactions in Controlled Fusion URL https://www.sciencedirect.com/science/article/pii/S0022311510008573

[33] Cartier-Michaud T, Ghendrih P, Sarazin Y, Abiteboul J, Bufferand H, Dif-Pradalier G, Garbet X, Grandgirard V, Latu G, Norscini C, Passeron C and Tamain P 2016 *Physics of Plasmas* **23** 020702 URL http://scitation.aip.org/content/aip/journal/pop/23/2/10.1063/1.4941974

[34] Beyer P, Benkadda S, Fuhr-Chaudier G, Garbet X, Ghendrih P and Sarazin Y 2007 *Plasma Physics and Controlled Fusion* **49** 507 URL http://stacks.iop.org/0741-3335/49/i=4/a=013

[35] Floriani E, Ciraolo G, Ghendrih P, Lima R and Sarazin Y 2013 *Plasma Physics and Controlled Fusion* **55** 095012 URL https://doi.org/10.1088/0741-3335/55/9/095012