The physical processes taking place at the separatrix and scrape-off layer regions are crucial for the operation of tokamaks as they govern the interaction of hot plasma with the vessel walls. Numerical modeling of the edge with state-of-the-art codes attempts to elucidate the complex interactions between neoclassical drifts, turbulence, poloidal, and parallel flows that control the physical set-up of the SOL region. Here, we present the post-processing analysis of simulation results from the gyrokinetic code XGC1, comparing and contrasting edge turbulence characteristics from a simulation of the DIII-D tokamak against a simulation of the Alcator C-Mod tokamak. We find that the equilibrium *E* ×* B* flux across the separatrix has a similar poloidal pattern in both discharges, which can be explained by $\u2207B$-drifts and trapped ion excursions. However, collisionality is noted to play a major role in the way that it prevents local charge accumulations from having more global effects in the C-Mod case. In both cases, turbulent electron heat flux is observed to be higher than the ion one and is possibly related to the need of electrons to maintain quasineutrality through the only channel available to them for exiting the confinement. By Fourier analysis, we identify turbulent frequencies and growth rates of the dominant mode in both simulations. In the case of C-Mod, these numbers point to the presence of a drift wave. In the DIII-D case, further linear simulations with the Gene code reveal a trapped electron mode. Furthermore, using a blob detection and tracking tool, we present the amplitude and size distributions of the blobs from both simulations. The amplitude distributions are in qualitative agreement with experimental observations, while the size distributions are consistent with the fact that most of the blobs are not connecting to the divertor plates and suggest that they are generated by the shearing of the turbulent modes.

## I. INTRODUCTION

A critical topic in contemporary tokamak research is the physics of the separatrix region where plasma particles and heat escape into the scrape-off layer (SOL). From the SOL, plasma impacts the divertor target plates and the main chamber walls of the device. In this region, many competitive processes are thought to be important. Poloidal and parallel flows carry plasma exhaust to the divertor, in competition with the cross field transport, both from neoclassical drifts and from plasma turbulence. Some combination of these and other processes establishes the poloidal and radial structures of the SOL plasma and the interaction of that plasma with various material surfaces. Controlling heat fluxes on those surfaces is critical for avoiding material damage; energy distributions and particle fluxes control many interactions of interest for fusion device performance and sustainability, including erosion rates, impurity sputtering, and neutral recycling. For all of these reasons, it is of great importance to develop a deeper understanding of edge and SOL plasma characteristics from experimental data, theory, and numerical modeling.

Qualitatively, the SOL plasma may be divided into near and far SOL regions; frequently, two scale lengths of profiles are observed.^{1} The near-SOL region is typically characterized by steep gradients, particularly in the heat flux channel, and is of primary importance for understanding heat flow to the divertor and plasma–material interactions at the divertor target. The near-SOL heat flux width has been the subject of many experimental, theoretical, and computational studies. Empirical scaling laws have been developed from a large international multi-machine database for both diverted^{2} and inner-wall limited^{3} discharges. Fluid turbulence simulations^{4} have described many features of the observed SOL gradient in the inner-wall limited case. The dominant dependence of the width with the poloidal magnetic field $Bp$ in diverted attached plasmas^{2} was modeled using fluid^{5} and kinetic^{6} physics and was described heuristically based on drifts.^{7} Large-scale computational efforts were undertaken to simulate the heat flux width and its scaling in both fluid^{8} and kinetic^{9} models. Of particular interest is the projection of the heat flux width to ITER, where the scaling at a high poloidal magnetic field $Bp$ becomes critical. This has motivated new high-field experimental results^{10} and recent modeling efforts. A topic of considerable interest is the interplay and competition of neoclassical drift physics and turbulence.^{8,9,11,12} If, as it is argued,^{7} neoclassical drift physics, with its characteristic poloidal Larmor radius width, dominates the Eich scaling for present devices, then as $Bp$ increases and that width narrows, the question arises as to whether a potentially larger turbulence-controlled width will lead to a broader SOL heat flux width. Although we do not address that question specifically in this paper, the interplay of drift–orbit mechanisms and turbulence is a motivating theme of our paper.

Turning to the far SOL, convective transport by coherent turbulent structures known as blob-filaments (also simply as “blobs” or “filaments”) is thought to be important. Observation of coherent structures in magnetic confinement devices has a long history dating back to observations on the Caltech Research Tokamak.^{13} As reviewed in Refs. 14 and 15, blob-filaments, observed in essentially all modern fusion research devices, carry particles, energy, and momentum into and through the SOL. Blob convection results from the *E* ×* B* radial drift of these structures, driven, for example, by interchange forces, which charge polarize the blob.^{16} Of interest for many applications are the net fluxes of particles and energy that are ultimately transported to the chamber walls and the resulting plasma profiles in the SOL. SOL profiles, fluxes, blobs, and statistical properties of SOL turbulence have been studied computationally, mostly with fluid models.^{17–25}

In order to understand and predict the net fluxes driven by blobs in the SOL, the scaling of the radial blob velocity with blob size, background plasma, and magnetic field parameters has received considerable attention. Theoretical predictions of this scaling have been made^{14,15,26–28} and tested against experimental^{29–32,52} and computational results.^{26,33–35} Less well understood, but under active investigation are the formation mechanisms of blobs,^{36,54} their generation rates,^{37,38} and their size distributions.^{39} Once these quantities are determined and understood, the SOL fluxes and profiles may be modeled using statistical methods.^{40,41} The blob size distribution and its relation to the underlying edge and SOL instabilities are topics that will be addressed in this paper.

Our paper presents a post-processing analysis of two simulations carried out with the electrostatic version of the gyrokinetic XGC1 particle-in-cell (PIC) simulation code. Both simulations were first used in Ref. 9 as part of an SOL heat-flux width scaling study. The first one is of a DIII-D discharge that formed the basis of an earlier post-processing analysis.^{12} In that paper, we discussed the development of the electrostatic potential and particle flows near the separatrix that were set up in large part by the neoclassical drift orbit excursions of ions. This mechanism was shown to control the poloidal profiles of the equilibrium $E\xd7B$ flows and fluxes. We showed that ambipolarity in the presence of the ion orbit loss was maintained by the turbulent loss of electrons across the separatrix and that the sheared flows set up by the ion orbit loss were of sufficient strength to affect the turbulence and impact of the poloidal profile of the turbulent fluxes.

In the present paper, we extend these results, comparing this DIII-D simulation with a C-Mod simulation of higher collisionality. Our paper has three main goals. One is to inquire into the generality of the DIII-D results described in the preceding paragraph for a higher collisionality, higher B-field device with respect to the relationship of drift–orbit and turbulence effects, sheared flows, and the poloidal profiles of the equilibrium and turbulent fluxes. A second goal is to extend the previous work by studying the linear stability properties of the plasma in the vicinity of the separatrix. For the DIII-D case, where an analysis with the Gene code^{42,43} is possible due to the fact that the dominant mode is localized within closed field lines, we identify the frequency, wavenumber, and driving gradients for the most unstable modes. The final goal of our paper is to investigate the properties of blobs in these simulations, in particular, their amplitude and size distribution, and to attempt to relate these observations to the instabilities and to theoretical results from the blob theory. This work extends earlier studies of blobs in XGC1.^{44}

The plan of our paper is as follows. In Sec. II, we describe and contrast the two simulation discharges under consideration. In Sec. III, we examine flux patterns across the separatrix for these two cases. Linear properties of these cases are investigated in Sec. IV, and in Sec. V, we consider some properties of the blob-filaments that result from these instabilities. Conclusions are given in Sec. VI.

## II. SIMULATIONS

In this section, we provide a description of the two simulations, domains and parameters. The simulations are of the neutral beam heated DIII-D^{45} discharge # 144981 and the RF heated Alcator C-Mod discharge # 1100223023, with both being H-modes. They are initialized with experimental profiles taken at times $3175\u2009ms$ and $1236\u2009ms$, respectively. The simulation inputs include experimental profiles of electron density and temperature (*n _{e}* and

*T*), ion temperature (

_{e}*T*), and magnetic equilibrium, from kinetic Equilibrium Fitting (EFIT) magnetic reconstructions, which can be seen in Fig. 1. The geometry in both cases is lower single null. The magnetic field is in the negative $\zeta \u0302$ direction, in cylindrical $(R,\zeta ,Z)$ coordinates, making the ions $\u2207B$-drift toward the lower X-point and the electrons toward the top (roughly at $90\xb0$ in the poloidal angle). In the rest of the paper, when we refer to the poloidal angle

_{i}*θ*, of a point on the separatrix, we will mean the angle in the

*R*–

*Z*plane of a point measured from the midplane with the magnetic axis (

*R*,

_{o}*Z*) taken as the center, i.e., $\theta =arctan(Z\u2212ZoR\u2212Ro)$. Positive angles are in the counterclockwise direction.

_{o}In Table I, we give some of the physical parameters of the two discharges. Some of them were given as input in the simulations, and some others were directly calculated from the output data. More specifically, we denote the total time of the simulation with *t*; *I _{p}* is the plasma current;

*B*is the toroidal magnetic field strength at the magnetic axis; $B\theta $ is the poloidal magnetic field measured at the outboard midplane separatrix (OMP),

_{tor}*n*; $TiOMP,\u2009TeOMP$ are the quasineutral density of the two species, ion and electron temperature, respectively, measured at OMP;

_{OMP}*n*is the Greenwald density;

_{G}*ρ*is the ion Larmor radius;

_{i}*C*is the sound speed at OMP; $V\theta $ is the poloidal flow at OMP;

_{s}*u*and

_{ti}*u*are the thermal velocities of ions and electrons at OMP, respectively;

_{te}*q*

_{95}is the safety factor;

*qR*is the connection length measured at the surface $\Psi N=0.95,\u2009L\u2225$ is the parallel connection length from the lower divertor plate to the top of the torus, measured at the open field-line region at $\Psi N=1.02,\u2009\u03f5=aR$ is the aspect ratio; $\omega transit=uteqR$ is the transit frequency of electrons;

*f*is the trapped particle fraction at OMP;

_{trapped}*τ*and

_{ii}*τ*are the ion–ion and electron–ion collision times, respectively; and $\lambda mfpi$ and $\lambda mfpe$ are the ion and electron mean free paths, respectively, all evaluated at OMP.

_{ei}Simulation parameters . | ||
---|---|---|

. | DIII-D . | C-Mod . |

$t\u2009(ms)$ | 0.16 | 0.085 |

$dt\u2009(s)$ | 2.3 $\xd710\u22127$ | 9.454 $\xd710\u22128$ |

$Ip\u2009(MA)$ | 1.5 | 0.9 |

$Btor\u2009(T)$ | 2.1 | 5.4 |

$B\theta \u2009(T)$ | 0.42 | 0.806 |

$nOMP\u2009(m\u22123)$ | 3.5 $\xd71019$ | 2.11 $\xd71020$ |

$TiOMP\u2009(eV)$ | 434 | 165 |

$TeOMP\u2009(eV)$ | 160 | 97 |

$nG\u2009(m\u22123)$ | 15.6 $\xd71019$ | 7.16 $\xd71020$ |

$\rho i\u2009(m)$ | 2 $\xd710\u22123$ | 7 $\xd710\u22124$ |

$Cs\u2009(ms)$ | 1.5 $\xd7105$ | 1.07 $\xd7105$ |

$V\theta \u2009(ms)$ | 4 $\xd7104$ | 1.0 $\xd7104\u2009to\u2009(\u22121.5)\xd7103$ |

$uti\u2009(ms)$ | 1.43 $\xd7105$ | 9.1 $\xd7104$ |

$ute\u2009(ms)$ | 5.1 $\xd7106$ | 4.5 $\xd7106$ |

q_{95} | 3.7 | 4.0 |

$qR\u2009(m)$ | 8.6 | 3.56 |

$L\u2225\u2009(m)$ | 17.9 | 9.1 |

ϵ | 0.34 | 0.294 |

$\omega transit\u2009(s\u22121)$ | 5.9 $\xd7105$ | 1.26 $\xd7106$ |

f _{trapped} | 77% | 72% |

$\tau ii\u2009(s)$ | 5 $\xd710\u22124$ | 24 $\xd710\u22126$ |

$\tau ei\u2009(s)$ | 3.5 $\xd710\u22126$ | 1.3 $\xd710\u22127$ |

$\lambda mfpi\u2009(m)$ | 71.5 | 2.2 |

$\lambda mfpe\u2009(m)$ | 17.9 | 0.6 |

Simulation parameters . | ||
---|---|---|

. | DIII-D . | C-Mod . |

$t\u2009(ms)$ | 0.16 | 0.085 |

$dt\u2009(s)$ | 2.3 $\xd710\u22127$ | 9.454 $\xd710\u22128$ |

$Ip\u2009(MA)$ | 1.5 | 0.9 |

$Btor\u2009(T)$ | 2.1 | 5.4 |

$B\theta \u2009(T)$ | 0.42 | 0.806 |

$nOMP\u2009(m\u22123)$ | 3.5 $\xd71019$ | 2.11 $\xd71020$ |

$TiOMP\u2009(eV)$ | 434 | 165 |

$TeOMP\u2009(eV)$ | 160 | 97 |

$nG\u2009(m\u22123)$ | 15.6 $\xd71019$ | 7.16 $\xd71020$ |

$\rho i\u2009(m)$ | 2 $\xd710\u22123$ | 7 $\xd710\u22124$ |

$Cs\u2009(ms)$ | 1.5 $\xd7105$ | 1.07 $\xd7105$ |

$V\theta \u2009(ms)$ | 4 $\xd7104$ | 1.0 $\xd7104\u2009to\u2009(\u22121.5)\xd7103$ |

$uti\u2009(ms)$ | 1.43 $\xd7105$ | 9.1 $\xd7104$ |

$ute\u2009(ms)$ | 5.1 $\xd7106$ | 4.5 $\xd7106$ |

q_{95} | 3.7 | 4.0 |

$qR\u2009(m)$ | 8.6 | 3.56 |

$L\u2225\u2009(m)$ | 17.9 | 9.1 |

ϵ | 0.34 | 0.294 |

$\omega transit\u2009(s\u22121)$ | 5.9 $\xd7105$ | 1.26 $\xd7106$ |

f _{trapped} | 77% | 72% |

$\tau ii\u2009(s)$ | 5 $\xd710\u22124$ | 24 $\xd710\u22126$ |

$\tau ei\u2009(s)$ | 3.5 $\xd710\u22126$ | 1.3 $\xd710\u22127$ |

$\lambda mfpi\u2009(m)$ | 71.5 | 2.2 |

$\lambda mfpe\u2009(m)$ | 17.9 | 0.6 |

Here, we should make an observation regarding the collisionality between the two discharges: in the case of DIII-D, we have $\u03f53/2\u22480.2$, and calculating the dimensionless Coulomb collisionality parameters $\nu \u22c6i=\nu iiq95Ruti\u22480.12$ and $\nu \u22c6e=\nu eiq95Rute\u22480.48$, we find that at the edge (OMP), the electrons are in the plateau regime ($\u03f53/2<\nu \u22c6e<1$), whereas the ions are in the banana transport regime ($\nu \u22c6i<\u03f53/2$). For C-Mod though, $\u03f53/2\u22480.16$ and $\nu \u22c6i\u22481.6,\u2009\nu \u22c6e\u22486.1$, which place both ions and electrons in the very collisional, Pfirsch–Schlütter regime. This is also reflected in the very different mean free paths in the two machines. Recalculating the dimensionless collisionalities at the open field-line region, for DIII-D, we get $\nu \u22c6e\u22481$ and $\nu \u22c6i\u22480.25$, which place electrons at the boundary between the plateau and Pfirsch–Schlütter regimes, whereas the ions are barely above the plateau regime. For C-Mod, the situation does not change since at the open field-line region, we obtain $\nu \u22c6e\u224815.6$ and $\nu \u22c6i\u22484.17$. To summarize, these numbers confirm that in the case of C-Mod, we have an overall highly collisional regime for both species, whereas in the case of DIII-D, we have an environment where ions on banana orbits that originate inside the last closed flux surface (LCFS) and overhang it are experiencing a relatively low collisionality even at their maximum excursion.

A second observation is related to the ion poloidal flow: despite the fact that sheared flows are present in both simulations, in the C-Mod case, the ion poloidal flow increases very rapidly approaching the separatrix. There, the shearing rate peaks and the flow get immediately sheared over a very small region, making it impossible to choose a meaningful single value for $V\theta $ at the edge. In the DIII-D case, the variation is much gentler as the flow approaching the separatrix is already strong and most of the shearing happens outside the LCFS. The variation of the poloidal flow near the edge at the outboard midplane regions of the two machines can be seen in Fig. 2.

## III. THE PARTICLE AND HEAT FLUXES

### A. Flux definitions

In this section, we will compare and contrast the particle and heat flux patterns from the two simulations. Our analysis will follow the one in Ref. 12, and some figures from there will be repeated here so that comparisons can be made between the edge of DIII-D and the more collisional edge of C-Mod. First, we start with some definitions

where *n* and *v* are the dynamical quantities of the plasma density and cross field velocity, respectively; and $\u27e8\cdots \u27e9$ denotes an average in either time (*t*) or toroidal planes $(\zeta )$, or both. The time averages considered here are taken over a time interval late in the simulations where a quasi-steady turbulent state has been achieved.

Employing Eq. (1), the fundamental relationship for the fluxes is

where the cross terms $\u27e8n\delta v\u27e9t,\zeta $ and $\u27e8v\delta n\u27e9t,\zeta $ vanish due to the vanishing of $\u27e8\delta n\u27e9t,\zeta $ and $\u27e8\delta v\u27e9t,\zeta $.

As we have pointed out in Ref. 12, these “fluxes” are not technically transport fluxes but rather local density weighted flows. Because of the very fast electron transit time, the radial excursions of electrons due to classical drifts such as *E* ×* B* and magnetic drifts nearly cancel out and negligibly contribute to net particle or heat transport. Therefore, to obtain a transport flux, one would have to integrate those density weighted flows over a flux surface. Here, as we did before, we will keep referring to them as fluxes for the sake of brevity. We recall that the first piece of the RHS of Eq. (2) is known as equilibrium flux and the second as turbulent flux (see Ref. 12 for additional details).

The definitions for the heat fluxes are similar, with the replacement of the density by the pressure of each species. Because the simulations are quasineutral, we cannot distinguish between the densities or the net particle fluxes of the two species. In the case of the heat fluxes though, a clear distinction between ions and electrons can be made, based on their different temperatures.

### B. Flux patterns across the separatrix

The first observation regarding the separatrix fluxes has to do with the degree to which the electrons satisfy the Maxwell–Boltzmann (MB) relation at the edge of the two discharges. On one hand, we see [Fig. 3(a)] that at the low field side (LFS) midplane of DIII-D, there is a significant deviation from MB, which we have ascribed to the fact that the transit frequency, the main frequency of the turbulence, and the collision rate are all of similar order of magnitude. On the other hand, in the C-Mod simulation [Fig. 3(b)], the MB relation is less modified by any other effects. Indeed, the measured frequency of turbulence in C-Mod, which is close to $12\u2009\u2009kHz$, is very low compared to the transit frequency of $1.26\u2009MHz$. Here, we have to remark that even in C-Mod, close to the X-point, the curves of Fig. 3(b) are very sensitive to the exact location where the density and potential are measured. The deviation from the Maxwell–Boltzmann relationship is stronger there, especially inside the separatrix. On the other poloidal locations though, the Maxwell–Boltzmann relationship is well satisfied.

Next, in Fig. 4, we present the *E* ×* B* particle fluxes around the LCFS. Although the DIII-D case [Fig. 4(a)] has been thoroughly analyzed in Ref. 12, we will recall here the main points: the inward flux observed at the HFS is due to Pfirsch–Schlütter flows caused by charge polarization from the opposite magnetic drifts of electrons and ions. The interesting alternating pattern of the inward and outward fluxes at the upper and lower LFS, respectively, is believed to be due to trapped particle excursions. More specifically, trapped ions that exit from the closed flux surfaces at the bottom create a charge hole. The local non-adiabaticity means that the electrons cannot move instantly along the field lines to neutralize the charge build-up. Therefore, a negative potential is created in order to attract more positive charges. The opposite situation takes place at the upper LFS where ions are re-entering. The smaller magnitude of this inward flux compared to that of the outward flux is ascribed to the fact that some ions that leave from the bottom never make it to the top because they get entrained in the parallel flow to the divertors.^{47}

The situation in C-Mod qualitatively resembles that in DIII-D; however, the circulation around X-point seems to dwarf the equilibrium fluxes elsewhere. Again, we have an inward flux at the HFS and an alternating inward and outward flux at the LFS with the same polarity as that in DIII-D. In the C-Mod case though, these features are not as pronounced compared to the very large particle circulation around the X-point. This *E* ×* B* circulation has been predicted and explained in Ref. 48: the downward magnetic drift of the ions, combined with the fact that the poloidal magnetic field at the X-point logarithmically vanishes, results in the loss of counter-traveling ions ($u\u2225<0$) from that point. This ion loss causes a charge and potential build-up in the region. The *E* ×* B* flow resulting from this potential has the direction we see in Fig. 4(c) and can be intuitively understood as the attempt of the plasma to compensate the weak parallel into-the-divertor flow of counter-traveling ions by forming an into-the-divertor *E* ×* B* flow.

These observations are confirmed from Figs. 4(b) and 4(d), where we draw the equipotential lines as a function of the poloidal angle and radial position. These plots cannot capture very rapid variations in such small regions and are meant to show large scale structure formation. The lower X-point is located near $\theta =\u22121.8$. We observe that in the case of closed field lines near the separatrix, the equipotential lines are straight, indicating that closed flux surfaces share the same potential. The situation, however, is very different when we move outside the separatrix or close to the X-point, even in the closed field line region. In the case of DIII-D [Fig. 4(b)], we find closed potential lines in the LFS. Those closed lines should be interpreted as closed contours of potential hills, which is in line with our understanding of charge build-up due to trapped ion excursions at these locations. Similarly, in the case of C-Mod [Fig. 4(d)], there are closed equipotentials outside the separatrix at the LFS, but the most prominent potential hill is at the X-point and is exactly the positive potential hill created by the accumulation of ions that cannot be properly dissipated due to the very weak poloidal motion at this point.^{49}

Even though the *E* ×* B* circulation around the X-point is expected,^{51} there still remains the issue of why it is so strong in C-Mod but relatively weak in DIII-D. A likely explanation has to do with the collisionality regimes of the two devices: in C-Mod, where the collisionality is very high (cf. Table I), any potential build-up remains strongly localized, whereas in DIII-D, which is practically collisionless, a potential perturbation can propagate around the flux surface, influencing the potential around the separatrix in poloidal locations far away from the X-point. In this fashion, in DIII-D, the potential hill at the X-point can be quickly neutralized by fast electrons traveling along the field lines, whereas in C-Mod, the potential hill remains relatively isolated since the spread of the potential perturbation is impeded by the plasma collisionality.

Next, we continue with presenting the poloidal patterns of the turbulent part of the *E* ×* B* flux in the two machines [Figs. 5(a)–5(c)] and relate them to the respective shear rates [Figs. 5(b)–5(d)]. In both machines, the turbulent *E* ×* B* flux is dominant at the LFS (cf. Table II for the linear mode properties of the two simulations there). In the case of DIII-D, this flux is interrupted, at a poloidal location above the outboard midplane, due to the very large local shear,^{50} as we can infer from Fig. 5(b). As we will see in Sec. IV, where we will evaluate the linear properties of the unstable modes, the local shear at this location surpasses the modes' growth rate. In the C-Mod case, on the other hand, we find a relatively uniform magnitude of *E* ×* B* flux, concentrated at the LFS. A visual inspection of Fig. 5(d) shows that the shearing rate is apparently not strong enough to influence the magnitude of the turbulent flux. In Figs. 6(a) and 6(b), we show the radial variation of the shear rate at the outboard midplane plotted together with the unstable mode in both machines. From these plots, we see that eventually, as we approach the separatrix, the shear rate picks up considerably in both machines. One thing to note is that in both cases, the value of the shearing rate strongly varies in a region between $\Psi =0.99$ and $\Psi =1.0$. Although this rapid change in the radial direction makes the plots of velocity shear in Figs. 5(b)–5(d) location sensitive, our main point in this figure is the qualitative structure of the poloidal variation of the shear and flux. An important point, to which we will return later at the section about the blobs in the two machines, is the fact that in C-Mod, the unstable mode peaks outside the LCFS, where $\Omega E\xd7B$ is strongest, whereas the mode of DIII-D, at least at the outboard midplane, is peaking inside the separatrix.

Simulation linear properties . | ||
---|---|---|

. | DIII-D . | C-Mod . |

$\omega lab\u2009(s\u22121)$ | 4 $\xd7106$ | 7.5 $\xd7104$ |

$\omega pl\u2009(s\u22121)$ | 2.5 $\xd7106$ | 3.14 $\xd7106$ |

$\omega \u22c6e\u2009(s\u22121)$ | 1.96 $\xd7106$ | 2.17 $\xd7106$ |

$\gamma GENE\u2009(s\u22121)$ | 2.54 $\xd7105$ | … |

$\gamma LP\u2009(s\u22121)$ | 1.74 $\xd7105$ | 3.5 $\xd7105$ |

$\gamma MHD\u2009(s\u22121)$ | 1.13 $\xd7106$ | 7.8 $\xd7105$ |

$k\u22a5\u2009(m\u22121)$ | 123 | 314 |

$k\u22a5\rho i$ | 0.246 | 0.22 |

$\Omega E\xd7B\u2009(s\u22121)$ | 4.7 $\xd7106$ to 9.4 $\xd7106$ | 2.5 $\xd7105$ to 2 $\xd7107$ |

Simulation linear properties . | ||
---|---|---|

. | DIII-D . | C-Mod . |

$\omega lab\u2009(s\u22121)$ | 4 $\xd7106$ | 7.5 $\xd7104$ |

$\omega pl\u2009(s\u22121)$ | 2.5 $\xd7106$ | 3.14 $\xd7106$ |

$\omega \u22c6e\u2009(s\u22121)$ | 1.96 $\xd7106$ | 2.17 $\xd7106$ |

$\gamma GENE\u2009(s\u22121)$ | 2.54 $\xd7105$ | … |

$\gamma LP\u2009(s\u22121)$ | 1.74 $\xd7105$ | 3.5 $\xd7105$ |

$\gamma MHD\u2009(s\u22121)$ | 1.13 $\xd7106$ | 7.8 $\xd7105$ |

$k\u22a5\u2009(m\u22121)$ | 123 | 314 |

$k\u22a5\rho i$ | 0.246 | 0.22 |

$\Omega E\xd7B\u2009(s\u22121)$ | 4.7 $\xd7106$ to 9.4 $\xd7106$ | 2.5 $\xd7105$ to 2 $\xd7107$ |

The last topic in the flux analysis comparison of the two machines concerns the respective heat fluxes. Because of the different temperatures of the two species, here we can make a distinction between the electron and ion fluxes. The observation made in Ref. 12 was that in DIII-D, the electron turbulent heat flux, shown in Fig. 7(a), was significantly larger than the ion one. This was related to the fact that although ions have many different channels for exiting the confinement, the electrons can only do so through turbulence. Therefore, even if the shearing rate is very strong, electron turbulence survives in order to maintain quasineutrality. When we repeat the same analysis in C-Mod, we see that this holds true as well: in Fig. 7(b), we find the electron turbulent heat flux to be almost everywhere much larger than the ion one, indicating that this might be a common feature of discharges whose species temperatures and scale lengths are of similar size. The only physical location where this stops being true is at the lower X-point with the very strong circulation. This is probably related to the pushing of the ion velocity space loss hole to higher energies due to the emergence of the confining electrostatic field.^{52}

## IV. LINEAR PROPERTIES OF THE SIMULATIONS

In this section, we will describe some properties of the linear modes present in the two simulations with the goal of relating those properties to the blobs that emerge. Of course, in a total-f code simulation, like the ones from XGC1, there can be many modes coexisting and interacting with each other. From here on, we will proceed by approximating this complex interaction of modes as a single unstable mode and find its properties by carrying out a Fourier analysis on XGC1 data. In the case of DIII-D, we also performed linear electrostatic simulations with the gyrokinetic code Gene, assuming that the character of the most unstable linear mode of the Gene simulation would give us some insight into the nature of the unstable mode of the original simulation. As will be seen below, comparison of results from the two methods shows that these assumptions are justified.

In the case of C-Mod, as evidenced by Fig. 2(b), running Gene simulations was not possible because the mode peak (in terms of mode strength $\delta nno$) is outside the last closed flux surface where the Gene code is not applicable. As a result, for this simulation, we have only the results of the Fourier analysis on the XGC1 output. Doing this revealed a turbulent frequency in the lab frame of about $fturb=12\u2009kHz$. The *E* ×* B* flow velocity at the mode location [the blue line in Fig. 2(b)] is of the order of $3\xd7104\u2009ms$, which, combined with a poloidal wave number of $k\theta =314\u2009m\u22121$, results in a plasma frame frequency of the order of the Doppler shift and in the electron diamagnetic direction. Because the poloidal *E* ×* B* varies rapidly and significantly over the region occupied by the mode (see Table I), it is not possible to define a precise value for the plasma frame frequency *ω _{pl}*. The value quoted in Table II is an upper estimate. This estimate is larger than but similar to the electron diamagnetic drift frequency $\omega \u2217e$; an average over the mode width would reduce

*ω*significantly. Other than that, we estimated the growth rate to be $\gamma =3.5\xd7105\u2009s\u22121$. This estimate comes from fitting an exponential at the initial phase of the simulation and should not be taken to signify the exact linear growth rate of the instability. The start of the simulation is when adjustment of the initial conditions to a non-local neoclassical equilibrium takes place and, unfortunately, there is no way that this process can be separated from the linear instability. Nevertheless, these measurements give us order-of-magnitude estimates about the properties of the linear mode, which we will later relate to the features of the turbulent blobs. For now, we remark that the growth rate is roughly $0.1\xb7\omega pl$ in C-Mod with $k\u22a5\rho i\u22480.2$, both of which are characteristics of drift waves. Similarly, both are very low compared to the transit frequency and collision rate, something that partly explains the fact that the electron response is very close to Maxwell–Boltzman. Moreover, the scale of the growth rate compared to the scale of the shearing rate, just inside the LCFS, explains why there is no suppression there, whereas the rapid increase in $\Omega E\xd7B$ outside the separatrix [cf. Fig. 6(b) and Table II] to a value of about two orders of magnitude larger than

_{pl}*γ*justifies why the mode gets shredded in filaments there.

As we mentioned above, for the DIII-D case, we have results from linear, electrostatic runs of Gene. In Fig. 2(a), we present the shape of the unstable mode in the *R*–*Z* plane, where we have plotted the mode strength $(\delta n)rms$ at the outboard midplane during the quasi-steady turbulent phase of the simulation. There, we see that the mode peaks inside the separatrix; therefore, we can run Gene at this location ($\Psi N=0.97$) and scan a range of parameters. All Gene simulations were local, with the center of the box being the mode peak location. The local parameters of density and electron and ion temperature were used (cf. Table I). Because the profiles vary widely across the mode, the scale lengths used for the simulation were an average of the scale lengths across the mode, weighted by the mode strength, i.e., $1Lx=\u222bdR\u2009(\delta n)2\u2207x\u222bdR\u2009(\delta n)2x$, where $x=n,Ti,Te$, and $(\delta n)2$ is the mode strength. Using the numbers $1Ln=42.45\u2009m\u22121,\u2009\u20091LTi=13.82\u2009m\u22121,1LTe=37.41\u2009m\u22121$, we employed the following scheme for linear runs: We labeled the case with the above parameters as our “base” case and carried out all different combinations of varying the inverse scale lengths between half and double their base value with three values in between. Therefore, we ended up with 243 different runs that demonstrate the response of the linear modes of the system to the variation of the different drivers. The geometric information was specified by the same EFIT file used in the original XGC1 simulation, and all Gene runs were confined to wavenumbers up to $k\u22a5\rho i=2.0$, which is the finest scale that the XGC1 mesh can distinguish.

In Fig. 8, we present the real frequency and linear growth rate of the base case Gene run. Out of the three unstable modes that Gene predicts, the first mode has a growth rate that exactly peaks at the wavenumber that we find from the XGC1 data and has a linear frequency in the same direction (electron diamagnetic) as the one measured in the simulation (Doppler shifted back into the plasma frame). Estimating the growth rate from the XGC1 data by fitting an exponential at the linear phase of the simulation, we find a growth rate of $\gamma =1.74\xd7105\u2009s\u22121$. With the reminder that this number is contaminated by the adjustment of the equilibrium, the Gene result of $\gamma =2.54\xd7105\u2009s\u22121$ is judged to be in satisfactory agreement. Furthermore, this growth rate is lower than the value that $\Omega E\xd7B$ takes above the midplane region, resulting in a suppression of the turbulent *E* ×* B* flux there [cf. Fig. 5(b)].

The response of the linear growth rate and the real frequency of the unstable mode to the variation of the three turbulent drivers, $\u2207n,\u2009\u2207Ti,\u2009\u2207Te$, between half and double their base value can be seen in Fig. 9. Each of the three lines represents five cases where the highlighted driver is varied between its extreme values and the other two scale lengths are kept at their base values. The values of frequency and growth rate recorded are the ones at $k\u22a5\rho i=0.246$. It is evident from the trends that the basic drivers of the mode are the density and electron temperature gradients and that the ion temperature gradient increase seems to have a mild stabilizing effect. We can attribute the latter to the influence of finite Larmor radius (FLR) effects, while the destabilizing effect of $\u2207n$ and $\u2207Te$, along with the direction of propagation ($\omega \u22c6e$) and the ion-scale of the instability, points us to a trapped electron mode. It is worth mentioning that we have also estimated the ideal MHD growth rate, $\gamma MHD=2cs2RLp$, where *L _{p}* is the scale length of the pressure. We found that $\gamma MHD=1.13\xd7106\u2009s\u22121$, which is of the same order of magnitude as the XGC1 measured turbulent frequency and larger than the Gene and XGC1 growth rates. This is an indication that interchange modifications must have a strong influence on the linear mode.

^{51}In Table II, we have gathered all the calculated numbers of the linear properties of the modes to aid the reader.

## V. PROPERTIES OF BLOB FILAMENTS

### A. Blob detection and tracking algorithm

The XGC1 simulations produce a large number of blob filaments, originating close to the edge and propagating into the SOL. In this work, we study the properties of these structures in relation to the unstable modes. We, therefore, developed a blob detection and tracking module in Python for use with the XGC1 data.

The method we used roughly follows in the steps of Ref. 52. Before we feed the data into the blob detection algorithm, we submit them to a two-stage preprocessing routine: First, we perform a two-dimensional smoothing over the range of each frame in the *R*–*Z* plane, where we replace the value of each pixel with a weighted average of its nearest neighbors. After that, we also smooth along the field lines since blob filaments are field-aligned structures, replacing the value of each pixel with the weighted average of its toroidal neighboring points that lie on the same field line.

The detection algorithm uses an off-the-shelf Python routine for locating contours above a certain threshold. The distinguishing characteristic that we require from a structure in order to qualify as a blob is that the density perturbation $\delta nn$ exceeds a certain threshold above the background, which, here, we took (arbitrary) to be 20%. After the Python routine tracks all contours above the said threshold, each blob is defined by its peak, that is, given that a set of contours contains the same maximum, we reject all but the largest one. In this way, we avoid double counting as, in the end, each blob contour contains a single, unique peak. Finally, an ellipse is fitted to each blob contour using the algorithm of Ref. 53, and a list of parameters, such as the blob area, peak value and location, length of major and minor axes of the ellipse, tilt angle, etc., are stored in an SQL database for easy retrieval. To give an idea about how the detection module works, we provide Fig. 10, where we see a density perturbation color plot of a time frame from the C-Mod outboard midplane. The blob contours are demarcated by colored lines, which enclose the blob peaks (illustrated by × symbols), and on top of them, we observe the black dashed lines of the fitted ellipse. The few peaks that have no enclosing contours or no fitted ellipses around them are either too small or too close to the edge of the frame so that we cannot find a closed contour. Those blobs are not included in the database and subsequent analysis.

The tracking feature of the module is implemented as follows: For each simulation time step, we take all blobs present at that time and use a comparison algorithm to compare them to each blob present in the following time step. The tracking algorithm computes a score based on the fact that the same blobs would have a radial and poloidal velocity that would follow a roughly normal distribution with means and variances chosen from experimental considerations and that their areas, like their velocities, should also not change dramatically between frames. Using this procedure, with the distribution functions properly calibrated, generally results in an unambiguous $1\u20131$ matching between blobs from adjacent time frames. Occasionally, we are able to track the splitting of one blob into two or more, or the inverse process of blob mergers. In this paper, we will not be presenting results from the tracking feature of the module.

### B. Blob features

Using the blob detection module, filaments are identified in both simulations after the linear phase. In the case of DIII-D, the total number of blobs that fulfill the criteria we set in order to include them in the statistics are 741, detected at the last 290 time steps (135 *μ*s) of the simulation, whereas in the case of C-Mod, it is 1332, detected at the last 500 time steps (94 *μ*s). In both cases, the poloidal range over which the analysis is carried out is roughly $(\u22126.5\xb0,6.5\xb0)$. Although there is some poloidal variation of the blob sizes, the analysis range avoids including in the statistical sample blobs close to the X-point that have been stretched due to flux expansion. Referring back to Figs. 6(a) and 6(b), we can see a possible reason for the larger number of C-Mod blobs in a shorter absolute time frame: The peak of the C-Mod unstable mode is outside the LCFS where the shearing rate is the strongest. This has the effect that the mode is shredded by shear, generating lots of filaments.^{54,55}

We focus on two static blob properties, namely, the amplitude and the size. The blob amplitude is compared to what has been experimentally observed and the blob size set against other characteristic sizes of the problem, and we draw conclusions about the nature of the generated blobs.

In Fig. 11, we give the probability distribution functions for the blob amplitudes in the two simulations. Both plots show the fitted Gaussian kernel density in blue (the y-axis units are in units of probability density). We see that blobs from both simulations have an amplitude that is roughly exponentially distributed, a fact that has been experimentally observed in C-Mod,^{56} MAST,^{41} TCV,^{57} and NSTX.^{52}

The blob size is shown in Fig. 12. To define a size for the blob, we have experimented with various choices for a blob profile (shape). We found that the most reasonable choice for our data was to assume that the blob profile corresponds to the positive part of a sinusoid that has been convected out to the measurement location. (Recall that nonlinearly, curvature-interchange dynamics propel positive fluctuations radially outward and negative fluctuations inwards.) Then, by “blob size,” we define the half width at half maximum (HWHM) of this sinusoid. In Appendix. A, we discuss the procedure by which we reconstruct this size of each blob from the data we have recorded in the database.

In Fig. 12, the histogram of the blob sizes is compared with the Larmor radius *ρ _{s}*, the dimensional characteristic blob scale size $\delta \u2217$, as defined in Ref. 26, $\delta \u2217=\rho s(L\u22252\rho sR)15$, with $L\u2225$ being the connection length and

*R*the major radius, and the HWHM of the linear mode (HWHM

_{LM}), defined by $HWHM=\pi 3k\u22a5$ (assuming it is a simple sine wave), with $k\u22a5$ taken from Table II. We note that even though the shape of these histograms is insensitive to the assumed form of the blob profile, different profile choices can change how the histogram maps to the x-axis. Nevertheless, whichever blob profile is postulated, the peak of the blob size distribution was found to have a scale smaller than that of the linear mode (HWHM

_{LM}), which may be consistent with sheared flows tearing up the linear structures. Also worth noting is the fact that there are very few blobs above the $\delta \u2217$ scale, which is consistent with the fact that none of the observed structures connects to the divertor plate sheaths. Recall here that in the vorticity charge conservation equation, blobs with $\delta <\delta \u2217$ are dominated by inertial (ion polarization drift) currents, while blobs with $\delta >\delta \u2217$ are dominated by the parallel current flow to the sheath.

^{14,26}This is to be expected since XGC1 effectively cuts off currents into the sheath implementing a logical sheath boundary condition, which modifies the sheath potential trying to enforce ambipolar fluxes to the wall.

^{58}

Although an analysis of blob transport is beyond the scope of this paper, the static blob properties in both the DIII-D and C-Mod simulations are seen to be qualitatively similar to experiments and theoretical expectations. Blob amplitude distributions decay rapidly, approximately exponentially, and the blob size distributions cluster between the Larmor radius and the size of the unstable mode, as determined from linear analysis. Thus, the blob properties in the two experiments are basically similar when viewed in terms of relevant physical spatial scales. Not studied here (apart from the observation of shear filamenting the C-Mod mode) but worthy of further investigation is the effect of shearing on the blob formation process.^{54,59}

## VI. SUMMARY AND CONCLUSIONS

We have presented and compared results from the analysis of two XGC1 simulations. Regarding the separatrix fluxes, both DIII-D and C-Mod simulations revealed a similar equilibrium *E* ×* B* poloidal flux pattern. This pattern has previously been ascribed, in the DIII-D case, to a combination of $\u2207B$-drifts and trapped ions exiting and re-entering the closed surface region. We corroborated this conclusion by providing the equipotential contour plots in the $\Psi \u2212\theta $ plane, which reveal potential structures localized at the LCFS that resemble charge accumulations. We showed that this pattern, in principle, holds in C-Mod as well, albeit much attenuated compared to the very strong X-point circulation. The reason for the X-point circulation is the well known X-point loss; however, we have attributed the big difference in the strength of it between the two machines to the rather large difference in their collisionalities. The high collisionality of C-Mod forces the X-point loss potential build-up to be localized, whereas, in the practically collisionless DIII-D case, any potential perturbation will travel around the flux surface and influence the potential at remote poloidal locations. This line of reasoning can also explain why in DIII-D the shear reaches (locally and very close to the separatrix) high enough levels in order to suppress the turbulent flux, but in the C-Mod case, it may remain too low at the separatrix to influence the turbulence. As far as the heat fluxes are concerned, we have shown that in both simulations, the electron turbulent heat flux at the edge is larger than the ion one. This is due to the very small Larmor radius of the electrons, which leaves them only turbulence as a means by which to exit the confinement and maintain quasineutrality.

For both simulations, Fourier analysis has revealed turbulent frequencies and an approximation of the growth rates of the unstable modes. In the case of DIII-D, we also performed a linear analysis using the Gene code. The linear analysis was found to be in close agreement with the measured frequencies and growth rates. The parameter scan that we did in order to find the response of the most unstable mode to the turbulence drivers, along with the frequency direction and spatial scale of the instability, indicates the presence of a TEM mode with strong interchange modifications. For the C-Mod case, because the mode peaks outside the LCFS, we could not perform a Gene simulation. However, the measured frequency, growth rate, and spatial scale of the mode suggest a drift wave.

We also presented the basic features of a blob detection and tracking module that we created for the XGC1 data. The module has the ability to identify individual blobs, fit an ellipse around them, and store relevant information about the filament in a database. It can also track blobs from one time frame to another, measuring their velocities. Here, we focused on two static features of the blobs, the amplitude and size distributions. The distribution of blob amplitudes in both the DIII-D and C-Mod simulations was found to be roughly exponential, in qualitative agreement with the previous experimental observations. Also, in both simulations, the distribution of blob sizes reveals that most of them cluster between the Larmor radius and the size of the unstable mode. This is expected from structures that are created by the shearing of the turbulent mode. Moreover, almost all of the blobs are smaller than the characteristic size $\delta \u2217$ that is relevant for blobs that connect to the divertor and are influenced by the current to the sheath. Therefore, we deduce that in both simulations, blob filaments are dominated by ion polarization currents. In future publications, we would like to explore the dynamic and static blob properties from these and other XGC1 simulations. Such simulations, in combination with information from experimental diagnostics, hold the important promise of a better understanding and predictive capability for plasma–material interactions in future devices.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy Office of Science, Office of Fusion Energy Sciences under Grant Nos. DE-FG02–08ER54954 and DE-SC000801 and by subcontract SO15882-C with PPPL under the U.S. Department of Energy HBPS SciDAC Project No. DE-AC02–09CH11466. The DIII-D U.S. DOE Grant number for the particular shot that was used as a basis for the analysis was DE-FC02–04ER54698 and for the C-Mod one, DE-FC02–99ER54512. The original XGC1 runs used computing resources on Titan at OLCF through the 2015 Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program and the 2016 ALCC (ASCR Leadership Computer Challenge) award.

## DATA AVAILABILITY

Digital data for all the figures of this paper (except Figs. 10 and 13) are openly available in Zenodo, at https://doi.org/10.5281/zenodo.3897755, Ref. 46.

### APPENDIX A: BLOB SIZE EXTRACTED FROM THE DATA

As we mentioned above, the blobs are assumed to have the shape of the positive part of a sinusoid. The blob detecting module records the peak, *p*, and level, *l*, values of the blob in the database. By level value, we mean the $\delta nno$ of the base contour of the blob, i.e., the largest contour that can still be considered part of the blob. In Fig. 13, we illustrate the situation, assuming that the blob is the positive part of the sinusoid oscillation $n=np\u2009sin\u2009(kx)=(p+1)sin\u2009(kx)$. Assuming that we know the size *d* at the level contour, then *d* is related to the wavenumber of the sine wave by $d=k(\pi \u22122\varphi )$, with $\varphi $ being the phase of the oscillation where we find the level contour. From the definition of the sine wave, this phase is given by $\varphi =arcsin(l+1p+1)$. Combining the previous two equations, we get a relationship between *d* and *k*. From that, we can arrive at the equation for the HWHM ($h=\lambda 6$)

Instead of the size *d*, at the database, we have stored the major radius *R _{M}* of the fitted ellipse. To find

*d*, we need to project this length into the binormal direction. This direction is taken to be $e\u0302\chi =b\u0302\xd7e\u0302\psi $, which, after a little manipulation, results in the explicit formula $e\u0302\chi =1|B|Bp(\u2212B\zeta BRe\u0302R\u2212B\zeta BZe\u0302Z+(BR2+BZ2)e\u0302\zeta ),$ with $|B|=BR2+BZ2+B\zeta 2$ and $Bp=BR2+BZ2$. The projected-to-the-binormal blob size is then

where