By determining the plasma fuelling, power exhaust, impurities and neutral dynamics, the scrape-off layer (SOL) is critical to the performance of fusion devices. The SOL is the outermost plasma region of a tokamak, characterised by magnetic field lines that intersect the wall, bounded on the inner side by the last closed flux surface (LCFS) and on the outer side by the vessel wall. In the SOL, the turbulent fluctuations are of order unity, meaning there is no separation between equilibrium and fluctuation quantities and no significant separation in their length scales, presenting significant challenges to both simulation and analytical progress. The properties of the turbulence change radially across the SOL in both tokamaks and stellarators, as evidenced for example by an increase in the relative size of fluctuations (Kube et al. Reference Kube, Garcia, Theodorsen, Brunner, Kuang, LaBombard and Terry2018; Niemann et al. Reference Niemann, Jakubowski, Effenberg, Bozhenkov, Cannas, Carralero, Langenberg, Pisano, Rahbarnia and Rudischhauser2020). In the near SOL, fluctuations are not intermittent, whereas in the far SOL their distribution has greater skewness and kurtosis, indicating intermittency (Boedo et al. Reference Boedo, Rudakov, Moyer, McKee, Colchin, Schaffer, Stangeby, West, Allen and Evans2003, Reference Boedo, Myra, Zweben, Maingi, Maqueda, Soukhanovskii, Ahn, Canik, Crocker and D’Ippolito2014; Kuang et al. Reference Kuang, LaBombard, Brunner, Garcia, Kube and Theodorsen2019). This is due to the existence of blobs – high density structures elongated along the magnetic field lines that propagate outwards due to their associated electric field (D’Ippolito, Myra & Zweben Reference D’Ippolito, Myra and Zweben2011). The properties of blobs have been measured extensively in tokamaks (Walkden et al. Reference Walkden, Militello, Harrison, Farley, Silburn and Young2016; Zweben et al. Reference Zweben, Myra, Davis, D’Ippolito, Gray, Kaye, Leblanc, Maqueda, Russell and Stotler2016; Tsui et al. Reference Tsui, Boedo, Myra, Duval, Labit, Theiler, Vianello, Vijvers, Reimerdes and Coda2018), stellarators (Sánchez et al. Reference Sánchez, van Milligen, Newman and Carreras2003), reversed field pinches (Spolaore et al. Reference Spolaore, Antoni, Spada, Bergsåker, Cavazzana, Drake, Martines, Regnoli, Serianni and Vianello2004) and basic plasma experiments (Antar et al. Reference Antar, Krasheninnikov, Devynck, Doerner, Hollmann, Boedo, Luckhardt and Conn2001; Carter Reference Carter2006; Furno et al. Reference Furno, Labit, Fasoli, Poli, Ricci, Theiler, Brunner, Diallo, Graves and Podestà2008a), however, the generation of blobs and prediction of blob-mediated transport remain open questions.
Recently, significant progress has been made in gyrokinetic (Dorf et al. Reference Dorf, Cohen, Compton, Dorr, Rognlien, Angus, Krasheninnikov, Colella, Martin and McCorquodale2012; Chang et al. Reference Chang, Ku, Tynan, Hager, Churchill, Cziegler, Greenwald, Hubbard and Hughes2017; Shi et al. Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017, Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2019; Pan et al. Reference Pan, Told, Shi, Hammett and Jenko2018) and fluid (Dudson et al. Reference Dudson, Umansky, Xu, Snyder and Wilson2009; Tamain et al. Reference Tamain, Bufferand, Ciraolo, Colin, Galassi, Ghendrih, Schwander and Serre2016; Häcker et al. Reference Häcker, Fuchert, Carralero and Manz2018; Paruta et al. Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018; Stegmeir et al. Reference Stegmeir, Coster, Ross, Maj, Lackner and Poli2018; Zhu, Francisquez & Rogers Reference Zhu, Francisquez and Rogers2018; Riva et al. Reference Riva, Militello, Elmore, Omotani, Dudson and Walkden2019) modelling of the SOL. Differences in the near and far SOL have been observed (Halpern & Ricci Reference Halpern and Ricci2017) and blobs have been detected and tracked (Nespoli et al. Reference Nespoli, Furno, Labit, Ricci, Avino, Halpern, Musil and Riva2017a; Paruta et al. Reference Paruta, Beadle, Ricci and Theiler2019). This represents a step forward, building upon earlier single blob studies (Myra, Russell & D’Ippolito Reference Myra, Russell and D’Ippolito2006) and multiblob models (Russell, Myra & D’Ippolito Reference Russell, Myra and D’Ippolito2007; Militello & Omotani Reference Militello and Omotani2016; Walkden et al. Reference Walkden, Wynn, Militello, Lipschultz, Matthews, Guillemaut, Harrison and Moulton2017). Leveraging these achievements, in the present Letter we disentangle the different turbulent mechanisms in the SOL, in particular, the nature of the fluctuations in the near and far SOL, the properties of blobs including their typical size, velocity and generation rate and the parallel transport. Ultimately, this allows us to develop a predictive model for the SOL turbulent transport and density decay lengths. These are key elements towards predicting the heat flux scale length, which is among the most critical issues for the operation of ITER and the design of DEMO (Loarte et al. Reference Loarte, Lipschultz, Kukushkin, Matthews, Stangeby, Asakura, Counsell, Federici, Kallenbach and Krieger2007; Zohm et al. Reference Zohm, Angioni, Fable, Federici, Gantenbein, Hartmann, Lackner, Poli, Porte and Sauter2013; Donné & Morris Reference Donné and Morris2018), and determining the wall recycling, impurity influx and wall erosion. We focus on the tokamak double-null (DN) magnetic configuration. Besides being of interest for DEMO (Wenninger et al. Reference Wenninger, Federici, Albanese, Ambrosino, Bachmann, Barbato, Barrett, Biel, Cavedon and Coster2016), this configuration facilitates the development of simple analytical estimates because the high and low field sides (HFS/LFS) are topologically separated.
Our analytical investigation is based on the results of two-fluid, three-dimensional (3-D) numerical simulations of the SOL dynamics carried out using the GBS code (Ricci et al. Reference Ricci, Halpern, Jolliet, Loizu, Mosetto, Fasoli, Furno and Theiler2012; Halpern et al. Reference Halpern, Ricci, Jolliet, Morales, Mosetto, Musil, Riva, Tran and Wersal2016; Paruta et al. Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018). GBS solves the drift reduced Braginskii equations (Zeiler, Drake & Rogers Reference Zeiler, Drake and Rogers1997) to evolve self-consistently the full profiles of the density, electrical potential, electron temperature and ion and electron parallel velocities with no separation between equilibrium and fluctuations. The simulation domain covers the full toroidal and poloidal angle regime and extends radially from 
${\sim}17\unicode[STIX]{x1D70C}_{s}$ inside the LCFS up to the wall. A source of heat and density within the LCFS mimics heat and plasma outflow from the core. The plasma flows along the field lines while being radially transported due to turbulence until it reaches the wall, which acts as a sink.
In the cold ion, electrostatic limit, the model equations are
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}n}{\text{d}t}=\frac{2}{B}[C(p_{e})-nC(\unicode[STIX]{x1D719})]-\unicode[STIX]{x1D735}_{\Vert }(nv_{\Vert e})+S_{n} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}t}=\frac{2B}{n}C(p_{e})-v_{\Vert i}\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D714}+\frac{B^{2}}{n}\unicode[STIX]{x1D735}_{\Vert }j_{\Vert } & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\text{d}T_{e}}{\text{d}t} & = & \displaystyle \frac{4}{3B}\left[\frac{7}{2}T_{e}C(T_{e})+\frac{T_{e}^{2}}{n}C(n)-T_{e}C(\unicode[STIX]{x1D719})\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{2}{3}\frac{T_{e}}{n}0.71\unicode[STIX]{x1D735}_{\Vert }j_{\Vert }-\frac{2}{3}T_{e}\unicode[STIX]{x1D735}_{\Vert }v_{\Vert e}-v_{\Vert e}\unicode[STIX]{x1D735}_{\Vert }T_{e}+S_{Te}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}v_{\Vert e}}{\text{d}t}=-v_{\Vert e}\unicode[STIX]{x1D735}_{\Vert }v_{\Vert e}+\unicode[STIX]{x1D707}\left(\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D719}-\frac{T_{e}}{n}\unicode[STIX]{x1D735}_{\Vert }n-1.71\unicode[STIX]{x1D735}_{\Vert }T_{e}+\unicode[STIX]{x1D708}j_{\Vert }\right) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}v_{\Vert i}}{\text{d}t}=-v_{\Vert i}\unicode[STIX]{x1D735}_{\Vert }v_{\Vert i}-\frac{1}{n}\unicode[STIX]{x1D735}_{\Vert }p_{e}, & \displaystyle\end{eqnarray}$$ with 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}$ the vorticity and 
$d_{t}\,f=\unicode[STIX]{x2202}_{t}f+\unicode[STIX]{x1D70C}_{\ast }^{-1}[\unicode[STIX]{x1D719},f]$ the convective derivative of field 
$f$. The differential operators are given by 
$\unicode[STIX]{x1D735}_{\Vert }f=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f$, 
$\unicode[STIX]{x1D6FB}_{\Vert }^{2}f=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f)$, 
$[\unicode[STIX]{x1D719},f]=\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\times \unicode[STIX]{x1D735}f)$, 
$C(f)=B[\unicode[STIX]{x1D735}\times (\boldsymbol{b}/B)]\boldsymbol{\cdot }\unicode[STIX]{x1D735}f/2$ and 
$\unicode[STIX]{x1D6FB}_{\bot }^{2}f=\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\boldsymbol{b}\times \unicode[STIX]{x1D735}f)\times \boldsymbol{b}]$ where 
$\boldsymbol{b}=\boldsymbol{B}/B$ is the unit vector of the magnetic field and 
$B$ is its norm. Unless specified otherwise, we normalise 
$n$, 
$T_{e}$, 
$\unicode[STIX]{x1D719}$ and 
$v_{\Vert e,i}$, 
$B$ and 
$t$ to 
$n_{0}$, 
$T_{e0}$, 
$T_{e0}/e$, 
$c_{s0}$, 
$B_{0}$ and 
$R_{0}/c_{s0}$, respectively, where 
$n_{0}$, 
$T_{e0}$ and 
$c_{s0}=\sqrt{T_{e0}/m_{i}}$ are the reference density, electron temperature and sound speed, 
$B_{0}$ and 
$R_{0}$ are the magnetic field strength and major radius at the magnetic axis. Perpendicular lengths are normalised to the ion sonic Larmor radius 
$\unicode[STIX]{x1D70C}_{s0}=c_{s0}/\unicode[STIX]{x1D6FA}_{ci}$, with 
$\unicode[STIX]{x1D6FA}_{ci}=eB_{0}/(cm_{i})$, and parallel lengths normalised to 
$R_{0}$. The normalised resistivity is defined based on the Spitzer resistivity as 
$\unicode[STIX]{x1D708}=m_{e}R_{0}/(1.96m_{i}c_{s0}\unicode[STIX]{x1D70F}_{e})$ where 
$\unicode[STIX]{x1D70F}_{e}$ is the electron collision time, assumed constant across the domain. We define 
$\unicode[STIX]{x1D707}=m_{i}/m_{e}$ and 
$\unicode[STIX]{x1D70C}_{\ast }=\unicode[STIX]{x1D70C}_{s0}/R_{0}$. The electron pressure is denoted 
$p_{e}=nT_{e}$ and dimensionless current 
$j_{\Vert }=n(v_{\Vert i}-v_{\Vert e})$. The coordinates 
$(x,y,z)$ refer to the radial and poloidal directions in units of 
$\unicode[STIX]{x1D70C}_{s0}$ and the toroidal angle in radians respectively. The source terms, 
$S_{n}$ and 
$S_{T_{e}}$, are Gaussian centred at a distance 
$11\unicode[STIX]{x1D70C}_{s0}$ inside the LCFS with half-width half-maximum (HWHM) 
$1.5\unicode[STIX]{x1D70C}_{s0}$ and amplitude 
$1.35n_{0}$. Magnetic presheath boundary conditions (Loizu et al. Reference Loizu, Ricci, Halpern and Jolliet2012) are used at the wall: 
$v_{\Vert i}=\pm \sqrt{T_{e}},v_{\Vert e}=\pm \sqrt{T_{e}}\exp (\unicode[STIX]{x1D706}-\unicode[STIX]{x1D719}/T_{e}),\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}=\mp \sqrt{T_{e}}\unicode[STIX]{x2202}_{x}v_{\Vert i},\unicode[STIX]{x2202}_{x}n=\mp n/\sqrt{T_{e}}\unicode[STIX]{x2202}_{x}v_{\Vert i},\unicode[STIX]{x1D714}=-(\unicode[STIX]{x2202}_{x}v_{\Vert i})^{2}\mp \sqrt{T_{e}},\unicode[STIX]{x2202}_{xx}^{2}v_{\Vert i},\unicode[STIX]{x2202}_{x}T_{e}=0$. The plus and minus signs refer to field lines entering and leaving the wall and 
$\unicode[STIX]{x1D706}\approx 3$. At the inner radial boundary, we use an ad hoc set of boundary conditions: 
$\unicode[STIX]{x2202}_{x}\,f=0$ for all fields 
$f$, except for 
$\unicode[STIX]{x1D714}$ and 
$\unicode[STIX]{x1D719}$, for which we impose 
$\unicode[STIX]{x1D714}=0$ and 
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D706}T_{e}$. The ad hoc inner boundary conditions have no effect in the region outside the source where our analysis is performed.
The axisymmetric magnetic field is based on three infinitely long wires aligned vertically, with the current in the central wire mimicking the plasma current. The upper and lower wires carry a current ten times stronger and are located at a distance 
$2a$ from the central wire with 
$a$ the radius at the wall. We then apply a radial transformation 
$x\rightarrow x-x_{0}$, where 
$x_{0}=0.9a$, to the flux function to obtain a configuration sufficiently circular to fit in the domain. The separatrix is shown in figure 1. Our analysis is independent of the details of the magnetic geometry.
By tuning the poloidal field strength, simulations are run with local safety factor 
$q=(a/\unicode[STIX]{x1D70C}_{\ast })\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}z/\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}y=4.3$, 6.5, 8.6 at the LCFS outer midplane. We also scan the parallel resistivity, 
$\unicode[STIX]{x1D708}=1,0.1$ and 0.01, typical experimental values for the normalised resistivity. Scans over these parameters are chosen since they are the most important controls on the plasma dynamics. We fix 
$\unicode[STIX]{x1D70C}_{\ast }^{-1}=500$ and 
$\unicode[STIX]{x1D707}=200$. The domain size is 
$L_{x}=120$ and 
$L_{y}=2\unicode[STIX]{x03C0}a=800$. Starting from a uniform initial state, the 3-D fields are evolved in time until a steady state is reached where the plasma influx from the source is balanced by losses to the wall and the toroidally averaged fields fluctuate around a constant value. We perform our analysis on this quasi-steady state. Details of the numerical implementation can be found in Paruta et al. (Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018).
A snapshot of 
$n$ and 
$\bar{n}$ on a poloidal plane for the 
$\unicode[STIX]{x1D708}=1$, 
$q=6.5$ simulation is shown in figure 1 (we use an overbar and a tilde to denote the time and toroidal average and the fluctuating component of all quantities, e.g. 
$n=\bar{n}+{\tilde{n}}$). We observe the turbulent LFS and quiescent HFS, experimentally observed in DN configurations (LaBombard et al. Reference LaBombard, Kuang, Brunner, Faust, Mumgaard, Reinke, Terry, Hughes, Walk and Chilenski2016). Removing the interchange instability drive (the curvature term in (0.2)) drastically reduces turbulence on the LFS, indicating that ballooning instabilities are the primary driver of turbulence in this region. Removing the Kelvin–Helmholtz drive (the 
$[\unicode[STIX]{x1D719},\unicode[STIX]{x1D714}]$ term in the convective derivative in (0.2)) suppresses the turbulence on the HFS but has a small impact on the LFS, indicating velocity shear and the Kelvin–Helmholtz instability play a dominant role only in the HFS.
Since most of the heat is exhausted on the LFS, we focus this Letter on the mechanisms determining the density decay length at the outer midplane. The density decay in our simulations cannot be properly described by a single exponential decrease, rather, it can be fitted with two exponentials characterised by a shorter decay length 
$L_{n}$ near the LCFS and longer decay length 
$L_{n}^{\prime }$ in the far SOL (see figure 2). Such a double decay length is a typical observation in a DN configuration, e.g. on C-Mod (LaBombard et al. Reference LaBombard, Kuang, Brunner, Faust, Mumgaard, Reinke, Terry, Hughes, Walk and Chilenski2016) and MAST (Riva et al. Reference Riva, Militello, Elmore, Omotani, Dudson and Walkden2019), as well as in single null (SN) (Carralero et al. Reference Carralero, Siccinio, Komm, Artene, D’Isa, Adamek, Aho-Mantila, Birkenmeier, Brix and Fuchert2017; Kuang et al. Reference Kuang, LaBombard, Brunner, Garcia, Kube and Theodorsen2019) and limited configurations (Horacek et al. Reference Horacek, Pitts, Adamek, Arnoux, Bak, Brezinsek, Dimitrova, Goldston, Gunn and Havlicek2016). It has also been observed in two-fluid simulations (Francisquez, Zhu & Rogers Reference Francisquez, Zhu and Rogers2017). The difference in the scale length is reflected in different turbulent properties in the near and the far SOL. As observed experimentally (Boedo et al. Reference Boedo, Rudakov, Moyer, McKee, Colchin, Schaffer, Stangeby, West, Allen and Evans2003; D’Ippolito et al. Reference D’Ippolito, Myra and Zweben2011; Kuang et al. Reference Kuang, LaBombard, Brunner, Garcia, Kube and Theodorsen2019), the fluctuation distribution is close to Gaussian in the near SOL with increasing skewness and kurtosis, indicative of intermittency, in the far SOL (figure 2). In the following, we identify the two different mechanisms setting 
$L_{n}$ and 
$L_{n}^{\prime }$. We call the width of the inner SOL 
$\unicode[STIX]{x1D6E5}$, this is the distance over which the density decays steeply. We refer to the density, temperature and radial turbulent particle flux at the separatrix by 
$\bar{n}$, 
$\bar{T_{e}}$ and 
$\unicode[STIX]{x1D6E4}$, and at the entrance of the far SOL (a distance 
$\unicode[STIX]{x1D6E5}$ from the separatrix) by 
$\bar{n}^{\prime }$, 
$\bar{T_{e}}^{\prime }$ and 
$\unicode[STIX]{x1D6E4}^{\prime }$, both at the outer midplane.
Figure 1. A snapshot of 
$n$ (a) and 
$\bar{n}$ (b) on a poloidal plane for the 
$\unicode[STIX]{x1D708}=1$, 
$q=6.5$ simulation. The separatrix is shown in black. A comparison of the simulation results and model prediction is shown in (c).
Figure 2. Radial profile of 
$\bar{n}$ at the outer midplane with the fitted exponentials (a), profiles of the skewness and kurtosis (b) and total turbulent transport and transport due to blobs (c). Above are two PDFs of the normalised density fluctuation, 
${\tilde{n}}/\unicode[STIX]{x1D70E}_{n}$, evaluated in the corresponding radius range, where 
$\unicode[STIX]{x1D70E}_{n}$ is the standard deviation evaluated locally.
We start by looking at the near SOL. Since the turbulent radial flux in the near SOL is not intermittent, we estimate the flux based on the development and saturation of a linear instability driven by a background radial gradient in density and temperature. We then match the predicted flux, 
$\unicode[STIX]{x1D6E4}$, to the turbulent flux across the LCFS, 
$\unicode[STIX]{x1D6E4}_{\text{LCFS}}$, to find 
$L_{n}$. It should be noted that although 
${\tilde{n}}(k_{y})$ has a broad spectrum, 
$\unicode[STIX]{x1D6E4}(k_{y})$ has a clear peak (Podestà et al. Reference Podestà, Fasoli, Labit, Furno, Ricci, Poli, Diallo, Müller and Theiler2008).
The turbulent flux can be written 
$\unicode[STIX]{x1D6E4}=\langle \overline{{\tilde{n}}\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D719}}}/B\rangle _{y}$, where the poloidal, 
$y$, average is evaluated over 
$45^{\circ }$ centred around the outer midplane. The density fluctuation can be estimated by noticing that linear instabilities saturate when the gradient of the fluctuations becomes comparable to the background density gradient, hence locally removing the turbulence drive (Ricci, Rogers & Brunner Reference Ricci, Rogers and Brunner2008; Ricci & Rogers Reference Ricci and Rogers2013), that is 
$\unicode[STIX]{x2202}_{x}{\tilde{n}}\sim \unicode[STIX]{x2202}_{x}\bar{n}$ or equivalently 
$k_{x}{\tilde{n}}\sim \bar{n}/L_{n}$, with 
$k_{x}$ the typical radial wavenumber of the perturbation, in agreement with 
${\tilde{n}}/\bar{n}$ in the simulations. We relate 
$\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D719}}$ to 
${\tilde{n}}$ by balancing the leading-order terms in the continuity equation, equation (0.1): 
$\unicode[STIX]{x1D6FE}{\tilde{n}}=\unicode[STIX]{x1D70C}_{\ast }^{-1}\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D719}}\unicode[STIX]{x2202}_{x}\bar{n}/B$ where 
$\unicode[STIX]{x1D6FE}$ is the linear growth rate of the instability driving the transport. The simulation test mentioned above shows that ballooning modes drive turbulence on the LFS. For these modes, as well as for drift waves (Rogers & Dorland Reference Rogers and Dorland2005), non-local linear theory shows that 
$k_{x}=\sqrt{k_{y}/L_{n}}$ (Ricci & Rogers Reference Ricci and Rogers2013). By considering the linear instability that maximises the transport, the turbulent flux 
$\bar{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D70C}_{\ast }\bar{n}(\unicode[STIX]{x1D6FE}/k_{y})_{\text{max}}$ follows.
In order to evaluate 
$(\unicode[STIX]{x1D6FE}/k_{y})_{\max }$, we linearise (0.1)–(0.5) assuming 
$C(f)\sim \unicode[STIX]{x2202}_{y}$, since the flux surfaces are approximately vertical in most of the region we are considering, and neglecting radial variation of the perturbation and poloidal variation of the equilibrium. We take 
$k_{\Vert }=2/q$, where 2 is the minimum parallel mode number, expected from ballooning stability and observed in the simulations. The linearised system that we obtain corresponds to that of the simple magnetic torus geometry (Poli et al. Reference Poli, Ricci, Fasoli and Podesta2008). We take 
$\unicode[STIX]{x1D702}=L_{n}/L_{Te}=0.77$, the theoretically expected value (Ricci et al. Reference Ricci, Rogers and Brunner2008), which is similar to the simulations. Using 
$\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{\text{LCFS}}$ we find numerically the 
$k_{y}$ and 
$L_{n}$ shown in figure 3 for three values of 
$\unicode[STIX]{x1D6E4}_{\text{LCFS}}$. The estimates of 
$L_{n}$ correspond well to the simulations, as shown in figure 1.
To understand the decrease in 
$k_{y}$ and increase in 
$L_{n}$ with increasing 
$\unicode[STIX]{x1D708}$ and 
$q$ (which has been observed experimentally Nespoli et al. Reference Nespoli, Labit, Furno, Horacek, Tsui, Boedo, Maurizio, Reimerdes, Theiler and Ricci2017b), we consider the limit (valid for typical parameters) in which the resistive ballooning mode is dominant: 
$R/L_{n}\gg 1$, 
$\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FE}\ll \unicode[STIX]{x1D708}$ and, to avoid coupling with sound waves and drift waves, 
$k_{\Vert }\ll \unicode[STIX]{x1D6FE}$ and 
$\unicode[STIX]{x1D714}_{\ast }\ll \unicode[STIX]{x1D6FE}$ where 
$\unicode[STIX]{x1D714}_{\ast }=k_{y}R/L_{n}$. Our dispersion relation reduces to 
$\unicode[STIX]{x1D6FE}^{2}-\unicode[STIX]{x1D6FE}_{i}^{2}+\unicode[STIX]{x1D6FE}k_{\Vert }^{2}/(\unicode[STIX]{x1D708}k_{y}^{2})=0$, where 
$\unicode[STIX]{x1D6FE}_{i}=\sqrt{2R(1+\unicode[STIX]{x1D702})/L_{n}}$ is the ideal ballooning growth rate, capturing well the strong transport limit (Halpern et al. Reference Halpern, Ricci, Jolliet, Loizu and Mosetto2014). Expressed in physical units to make explicit the 
$R_{0}$ dependence,
$$\begin{eqnarray}L_{n}=\frac{2}{3}(1+\unicode[STIX]{x1D702})\bar{c_{s}}^{2}\left(\frac{q\bar{n}}{2\unicode[STIX]{x1D6E4}_{\text{LCFS}}}\right)^{4/3}\left(\frac{m_{e}}{m_{i}1.96\unicode[STIX]{x1D70F}_{e}}\right)^{2/3}\unicode[STIX]{x1D6FA}_{ci}^{-4/3}R_{0}^{1/3},\end{eqnarray}$$ which is of the order of 1 mm for typical experimental parameters in TCV and C-Mod (LaBombard et al. Reference LaBombard, Kuang, Brunner, Faust, Mumgaard, Reinke, Terry, Hughes, Walk and Chilenski2016; Vianello et al. Reference Vianello, Carralero, Tsui, Naulin, Agostini, Cziegler, Labit, Theiler, Wolfrum and Aguiam2019). We can use this equation to estimate 
$k_{y}$ by noticing that the peak 
$\unicode[STIX]{x1D6FE}/k_{y}$ occurs approximately where the damping term and ballooning drive are equal. Using 
$\unicode[STIX]{x1D6FE}\sim \unicode[STIX]{x1D6FE}_{i}$, we find 
$k_{y}=(2^{2/3}/3^{1/4})[\bar{n}\unicode[STIX]{x1D6FA}_{ci}^{2}(1.96m_{i}\unicode[STIX]{x1D70F}_{e})/(q^{2}\unicode[STIX]{x1D6E4}_{\text{LCFS}}R_{0}^{2}m_{e})]^{1/3}$ in physical units, the full numerical result is shown in figure 3.
Figure 3. Predicted 
$L_{n}$, 
$k_{y}$, 
$v_{b}$, 
$\unicode[STIX]{x1D6E4}^{\prime }$ and 
$L_{n}^{\prime }$ as a function of 
$\unicode[STIX]{x1D708}$, 
$q$ and 
$\unicode[STIX]{x1D6E4}_{\text{LCFS}}$. All reference quantities are taken at the LCFS.
We now turn to the far SOL, where the fluctuation distribution is heavy tailed, indicating intermittent turbulence and indeed observation of the simulation results reveals the presence of coherent structures of high plasma density (blobs) that propagate outwards due to their self-generated 
$\boldsymbol{E}\times \boldsymbol{B}$ velocity.
We use a pattern recognition algorithm described in Paruta et al. (Reference Paruta, Beadle, Ricci and Theiler2019) to track the blobs (defined here as coherently propagating structures of amplitude greater than 2.5 times the standard deviation of 
$n$) and measure their size, amplitude and velocity. Following Nespoli et al. (Reference Nespoli, Furno, Labit, Ricci, Avino, Halpern, Musil and Riva2017a), we calculate the fraction of the cross-field transport due to blobs by assuming a 2-D Gaussian density distribution of each blob in the poloidal plane with a peak density fluctuation 
$n_{b,i}$ and radial and poloidal HWHM 
$a_{x,i}$ and 
$a_{y,i}$, where 
$i$ is the blob index. The blob flux is calculated by 
$\unicode[STIX]{x1D6E4}_{b}(x,y)=\sum _{i}n_{b,i}v_{b,i}\exp [(x-x_{b,i})^{2}/(2a_{x,i}^{2})+(y-y_{b,i})^{2}/(2a_{y,i}^{2})]$, where the sum is carried out over all blobs and 
$(x_{b,i},y_{b,i})$ are the blobs’ centre of mass. We find that blob transport dominates in the far SOL (figure 2), consistent with the result of Nespoli et al. (Reference Nespoli, Furno, Labit, Ricci, Avino, Halpern, Musil and Riva2017a) and previous experimental works that found blobs to contribute an order unity fraction of the particle flux (Boedo et al. Reference Boedo, Rudakov, Moyer, McKee, Colchin, Schaffer, Stangeby, West, Allen and Evans2003, Reference Boedo, Myra, Zweben, Maingi, Maqueda, Soukhanovskii, Ahn, Canik, Crocker and D’Ippolito2014; D’Ippolito et al. Reference D’Ippolito, Myra and Zweben2011).
We now predict the flux due to blobs using only the near SOL properties. For this purpose, we express the blob flux averaged in the poloidal plane (Russell et al. Reference Russell, Myra and D’Ippolito2007) 
$\unicode[STIX]{x1D6E4}^{\prime }=\langle \unicode[STIX]{x1D6E4}_{b}\rangle _{x,y}=\unicode[STIX]{x1D70E}_{b}f_{b}v_{b}$ in terms of 
$\unicode[STIX]{x1D70E}_{b}$ the average density inside a blob, 
$f_{b}$ the blob packing fraction (ratio of area covered by blobs to total SOL area) and 
$v_{b}$ the average blob velocity. The 
$x$ average is taken from the LCFS to the wall. We address each of these quantities in turn.
We estimate 
$\unicode[STIX]{x1D70E}_{b}=2n_{b}/\text{ln}(2)$, i.e. as the ratio of the average number of particles in a blob, 
$2\unicode[STIX]{x03C0}n_{b}a_{x}a_{y}/\text{ln}(2)$, where we assume blobs have on average a Gaussian shape with HWHM 
$a_{x}$ and 
$a_{y}$ and peak density 
$\bar{n}^{\prime }+n_{b}$, and the average blob area, 
$A_{b}=\unicode[STIX]{x03C0}a_{x}a_{y}$. Since 
$\bar{n}^{\prime }$ decreases radially, 
$n_{b}/\bar{n}^{\prime }$ remains approximately constant over the blob lifetime despite parallel draining, so we combine the definition of a blob and the estimate of the 
${\tilde{n}}/\bar{n}$ in the near SOL to estimate 
$n_{b}\sim 3{\tilde{n}}\sim 3\bar{n}/(L_{n}k_{x})$.
The packing fraction 
$f_{b}=N_{b}A_{b}/(A_{SOL})$, where 
$N_{b}$ is the number of blobs, requires an estimate of the blob size. We observe that the blob size remains approximately constant as the blobs propagate (as observed experimentally in AUG Carralero et al. Reference Carralero, Manz, Aho-Mantila, Birkenmeier, Brix, Groth, Müller, Stroth, Vianello and Wolfrum2015) and that blobs tend to be circular, maximising their Kelvin–Helmholtz stability (Ricci & Rogers Reference Ricci and Rogers2013), so we estimate their size as the geometric mean of the near SOL eddy dimensions 
$a_{x}\sim a_{y}\sim \unicode[STIX]{x03C0}/(2\sqrt{k_{x}k_{y}})$. We infer, therefore, from the results for the near SOL analysis that blob size increases with resistivity, a trend observed both in our simulations (figure 4) and experimentally Vianello et al. (Reference Vianello, Carralero, Tsui, Naulin, Agostini, Cziegler, Labit, Theiler, Wolfrum and Aguiam2019).
Figure 4. Distribution of the vertical size (a) and radial velocity (b) of the blobs at 
$\unicode[STIX]{x1D708}=1$ and 
$\unicode[STIX]{x1D708}=0.01$ with 
$q=6.5$. The mean sizes and velocities are shown with a solid line and the predictions with a dashed line.
We now turn to the estimation of 
$N_{b}$. In steady state, the blob generation and loss rates are equal. Since blobs are generated from instabilities of wavelength 
$2\unicode[STIX]{x03C0}/k_{y}$, we expect the generation rate to be proportional to 
$L_{y}k_{y}/(2\unicode[STIX]{x03C0})$. The generation time scale has previously been proposed as determined by poloidal flow shear (Furno et al. Reference Furno, Labit, Podestà, Fasoli, Müller, Poli, Ricci, Theiler, Brunner and Diallo2008b) or a combination of flow shear and mode phase velocity (D’Ippolito et al. Reference D’Ippolito, Myra and Zweben2011; Fuchert et al. Reference Fuchert, Carralero, Manz, Stroth and Wolfrum2016). In the presence of hot ions, strong 
$\boldsymbol{E}\times \boldsymbol{B}$ flow shear may be present (Zhu, Francisquez & Rogers Reference Zhu, Francisquez and Rogers2017; Paruta et al. Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018). However, in our simulations, we find a flow shear time scale, 
$\unicode[STIX]{x2202}_{x}v_{E\times B}$, almost an order of magnitude less than the observed generation time scale and not to scale with 
$\unicode[STIX]{x1D708}$ and 
$q$, as observed. We reason that blobs are created because the linear instability saturates as the local density gradient is removed and the resulting density perturbation moves outwards without the streamer being refilled from the core, a case which was studied in a basic plasma physics device in Müller et al. (Reference Müller, Theiler, Fasoli, Furno, Labit, Tynan, Xu, Yan and Yu2009). Hence, the generation rate is limited by the time taken for the blob to travel one radial wavelength of the driving instability, 
$4a_{x}/v_{b}$, allowing for the density gradient to be re-established, which is consistent with the simulations. Taking the blob lifetime as the time taken to cross the domain, the loss rate is 
$N_{b}v_{b}/L_{x}$. Hence, 
$N_{b}=4\unicode[STIX]{x03C0}^{2}L_{x}L_{y}/(k_{x}k_{y})$. Using the above relations, we find 
$f_{b}\approx \unicode[STIX]{x03C0}/16$, independent of the SOL parameters. While the universality of 
$f_{b}$ is well supported by the simulation data, the predicted value is an overestimate (likely because we assume all blobs cross the entire radial domain) and a better estimate is 
$f_{b}=0.1$.
The average blob velocity, 
$v_{b}$, is deduced from the average blob size, 
$a_{x}$ and 
$a_{y}$, according to the well-studied size–velocity scaling relations derived using the two-region model (Myra et al. Reference Myra, Russell and D’Ippolito2006; D’Ippolito et al. Reference D’Ippolito, Myra and Zweben2011; Tsui et al. Reference Tsui, Boedo, Myra, Duval, Labit, Theiler, Vianello, Vijvers, Reimerdes and Coda2018; Paruta et al. Reference Paruta, Beadle, Ricci and Theiler2019) (figure 4). The normalised blob velocity, 
$\hat{v}=\text{Im}(\hat{\unicode[STIX]{x1D714}})\hat{a}^{1/2}$, depends on the normalised frequency 
$\hat{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FE}_{b}$ (where 
$\unicode[STIX]{x1D6FE}_{b}^{2}=2\bar{T_{e}}^{\prime }\unicode[STIX]{x1D70C}_{\ast }^{-1}n_{b}/[a_{x}\bar{n}^{\prime }]$ represents the local ballooning drive), which is determined by the dispersion relation 
$1+\hat{\unicode[STIX]{x1D714}}^{2}+\text{i}\hat{\unicode[STIX]{x1D714}}\unicode[STIX]{x1D6E9}/\unicode[STIX]{x1D6EC}=0$, for 
$\unicode[STIX]{x1D6EC}>1$, or 
$1+(1+f^{2})\hat{\unicode[STIX]{x1D714}}^{2}+\text{i}\unicode[STIX]{x1D6E9}\hat{\unicode[STIX]{x1D714}}=0$, for 
$\unicode[STIX]{x1D6EC}<1$, where 
$\unicode[STIX]{x1D6EC}=\unicode[STIX]{x1D708}\bar{n}^{\prime }L_{\Vert 1}^{2}/[L_{\Vert 2}\unicode[STIX]{x1D70C}_{s}]>1$ is the resistivity parameter, 
$L_{\Vert 1}$ is the field line length from the midplane to the region of maximum flux fanning, 
$L_{\Vert 2}$ is the field line length from the region of maximum flux fanning to the wall (Paruta et al. Reference Paruta, Beadle, Ricci and Theiler2019) and 
$f$ is the flux tube fanning (Myra et al. Reference Myra, Russell and D’Ippolito2006). The blob size parameter, 
$\unicode[STIX]{x1D6E9}=\hat{a}^{5/2}$, where 
$\hat{a}=a_{b}/a_{\ast }$ is the normalised blob size, with 
$a_{b}=(2a_{y}/\unicode[STIX]{x03C0})^{4/5}a_{x}^{1/5}$ and 
$a_{\ast }=[2\unicode[STIX]{x1D70C}_{s}^{4}L_{\Vert 2}^{2}n_{b}/(a_{x}\bar{n}^{\prime }\unicode[STIX]{x1D70C}_{\ast })]^{1/5}$. Finally the blob velocity is given by 
$v_{b}=0.5\hat{v}v_{\ast }$, with the reference velocity 
$v_{\ast }=\unicode[STIX]{x1D70C}_{s}[2\unicode[STIX]{x03C0}^{2}a_{x}^{2}\unicode[STIX]{x1D70C}_{s}^{2}\unicode[STIX]{x1D70C}_{\ast }^{2}L_{\Vert 2}n_{b}/(\bar{n}^{\prime }a_{y}^{2})]^{1/5}$ and the factor 0.5, obtained by comparing the scaling with the simulation results, accounting for the fact that our estimate is an upper limit that neglects various mechanisms slowing the blobs (Tsui et al. Reference Tsui, Boedo, Myra, Duval, Labit, Theiler, Vianello, Vijvers, Reimerdes and Coda2018).
Finally, we determine 
$L_{n}^{\prime }$ by balancing the divergence of the blob flux with the divergence of the parallel flow 
$L_{n}^{\prime }=\unicode[STIX]{x1D6E4}^{\prime }L_{\Vert }^{\prime }/(\bar{n}^{\prime }\bar{c_{s}}^{\prime })$ with 
$L_{\Vert }=L_{\Vert 1}+L_{\Vert 2}$, 
$\bar{n}^{\prime }=\bar{n}\exp (-\unicode[STIX]{x1D6E5}/L_{n})$ and 
$\bar{c_{s}}^{\prime }=\bar{c_{s}}\exp (-\unicode[STIX]{x0394}\unicode[STIX]{x1D702}/[2L_{n}])$. We remark that 
$L_{n}^{\prime }$ depends only weakly on 
$\unicode[STIX]{x1D6E5}$ since 
$\unicode[STIX]{x1D6E4}^{\prime }$ also scales approximately with 
$\bar{n}^{\prime }\bar{c_{s}}^{\prime }$. For typical experimental parameters, most blobs are in the 
$\unicode[STIX]{x1D6EC}>\unicode[STIX]{x1D6E9}$ regime, for which
$$\begin{eqnarray}L_{n}^{\prime }=\frac{7.3f_{b}L_{\Vert }\unicode[STIX]{x1D6FA}_{ci}^{0.0167}(R_{0}m_{i}\unicode[STIX]{x1D70F}_{e})^{0.00833}\unicode[STIX]{x1D6E4}_{\text{LCFS}}^{1.04}}{m_{e}^{0.00833}\bar{c_{s}}^{1.05}\bar{n}^{1.04}(1+\unicode[STIX]{x1D702})q^{0.0167}\unicode[STIX]{x1D70C}_{s}^{\prime 0.1}}\end{eqnarray}$$ is of the order of several mm for typical experimental parameters. We note that, since 
$\unicode[STIX]{x1D6EC}\propto R_{0}$ and 
$\unicode[STIX]{x1D6E9}\propto R_{0}^{1.46}$, larger device will likely have blobs in the 
$\unicode[STIX]{x1D6E9}>\unicode[STIX]{x1D6EC}$, 
$\unicode[STIX]{x1D6EC}>1$ regime, for which
$$\begin{eqnarray}L_{n}^{\prime }=\frac{94.8f_{b}L_{\Vert 1}^{2}L_{\Vert }\unicode[STIX]{x1D6E4}_{\text{LCFS}}^{5/3}m_{i}^{17/15}\unicode[STIX]{x1D70F}_{e}^{2/15}\unicode[STIX]{x1D70C}_{s}^{\prime 7/5}\unicode[STIX]{x1D6FA}_{ci}^{19/15}}{R_{0}^{28/15}\bar{n}^{5/3}m_{e}^{17/15}c_{s}^{14/5}\unicode[STIX]{x1D70C}_{s}^{2/5}(1+\unicode[STIX]{x1D702})^{7/5}q^{34/15}}.\end{eqnarray}$$ In figure 3 we show 
$v_{b}$, 
$\unicode[STIX]{x1D6E4}^{\prime }$ and 
$L_{n}^{\prime }$ as a function of 
$\unicode[STIX]{x1D708}$, 
$q$ and 
$\unicode[STIX]{x1D6E4}_{\text{LCFS}}$. We observe that 
$\unicode[STIX]{x1D6E4}^{\prime }$ increases with 
$\unicode[STIX]{x1D708}$ and with 
$q$, primarily due to variation in 
$v_{b}$ and to a lesser extent 
$\unicode[STIX]{x1D70E}_{b}$, as suggested in Russell et al. (Reference Russell, Myra and D’Ippolito2007). The increase in 
$v_{b}$ follows from the 
$\hat{a}-\hat{v}$ scaling. Such an increase has also been observed in gyrofluid simulations (Häcker et al. Reference Häcker, Fuchert, Carralero and Manz2018). The increase in 
$L_{n}^{\prime }$ with resistivity is well documented experimentally (D’Ippolito et al. Reference D’Ippolito, Myra and Zweben2011). The predicted 
$L_{n}^{\prime }$ is compared to the simulation result in figure 1. As for the near SOL, we find good agreement between theory and simulation.
Acknowledgements
The authors thank O. Fevrier, I. Furno, D. Galassi, B. Labit, J. Loizu, C. Theiler, C. Tsui and N. Vianello for useful discussions. The simulations presented herein were carried out in part at the Swiss National Supercomputing Center (CSCS) under the project ID s882 and in part on the CINECA Marconi supercomputer under the GBSedge project. This work was supported in part by the Swiss National Science Foundation, was carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
