Analysis of equilibrium and turbulent fluxes across the separatrix in a gyrokinetic simulation

Studies of the plasma edge are indispensable for a variety of tokamak operation and performance issues such as the density and power scrape-off-layer (SOL) widths,^1,3–10 intrinsic rotation and momentum transport,^11–17 edge flows^18–23 and the L-H transition.^24–26 

Out of all the above-mentioned important and inter-related facets of edge physics, the SOL width is crucial in tokamak operation as it governs the power density delivered at the divertor plates. Indeed, a viable fusion reactor would need to have plasma-facing materials that can withstand large heat loads without the need of regular replacement. However, an ITER exhaust power of 100 MW, concentrated on a speculated ∼1 mm SOL width would produce a whopping $5.5 GWm2$ of parallel power density, placing stringent constraints on the lifetime of the divertor materials.²⁷ It is therefore of utmost importance to be able not only to predict the SOL width for a given discharge based on its parameters but also to gain an understanding of the physical mechanisms that control its size. To this end, Eich²⁸ has constructed scaling relations for the power SOL width, λ_q, based on multiple regressions of relevant parameters from various experimental discharges. He has arrived at a scaling law of the form $λq∼Bp−1.19$⁠, where B_p is the poloidal magnetic field, based on which the SOL width of ITER is inferred to be 1 mm. Goldston² has provided theoretical justification for the above scaling by devising a heuristic model that predicts the SOL width. In summary, this model rests on the assumption that the perpendicular ∇B and curvature drift ion flows are balanced exactly by parallel flows with one half of the flux being directed to the divertors at a speed of $Cs2$ and the other half returning upstream via the Pfirsch-Schlütter current. This mechanism is responsible for setting the density SOL width, λ_n. The next assumption is that this density channel fills with turbulent electrons and empties through Spitzer-Härm conductivity, establishing the power SOL width. The heuristic model correctly captures the inverse relation between λ_q and B_p proposed by Eich and can thus, with reasonable accuracy, reproduce the experimental results.

Until recently, simulations of the edge have relied upon reduced fluid models in simplified geometries omitting essential physics. This changed with the introduction of XGC1^26,29–32 which is a full-f, Particle-In-Cell, 5D gyrokinetic code using realistic X-point geometry, optimized to simulate the tokamak edge. The code has gyrokinetic ions and offers an option for drift-kinetic, fluid, or adiabatic electrons. Employing it in a study of ITER SOL width, Chang et al.³ discovered that the simulated width is six times larger than the one predicted by the Eich fit or the Goldston heuristic model, casting doubt on whether the heretofore assumed scaling relations continue to be valid at the ITER regime. The explanation that they provide regarding this disparity has to do with the importance of edge electron turbulence^33–35 in the new regime. More specifically, they claim that the electrostatic-“blobby” edge turbulence is setting the scale of the electron channel. The size of a blob or a streamer (or any other coherent turbulent structure) is, at least partially, related to the width of the edge shearing layer,^36,37 wherein the blobs are being generated. In present day tokamaks, this shearing layer width is relatively small and the density channel is dominated by ion neoclassical flows inside the SOL. On ITER, however, due to the much larger size of the shearing layer, the radial extent of blobs and streamers is bigger and turbulence dominates the density width.

From the above discussion, it is clear that important figures of merit such as the SOL width are the result of a non-linear interaction between particle orbit effects (e.g., neoclassical effects and ion X-point losses^38,39) with the electrostatic-“blobby” edge turbulence. The interaction is through the radial and poloidal edge electric field, established by the ion losses and the sheared radial flows that it generates which can influence both the turbulence and the ion flows. Such interplay can often yield unexpected results and thus merits further study with state-of-the-art codes which might illuminate certain aspects of this complex process. Here, we analyze the results of an XGC1 simulation of DIII-D from the viewpoint of equilibrium and turbulent E × B particle and heat fluxes across the last closed flux surface (LCFS). Within this framework, we attempt to understand the relationship between the direct ion orbit losses and the electrostatic potential that drives these fluxes as well as the effect of the generated shear on the turbulence. Furthermore, we try to measure the impact of turbulence on the electron transport and thereby, test the second assumption of the heuristic model which states that only turbulent electrons are responsible for filling the density SOL channel.

The organization of the paper is as follows: In Sec. II, we give a brief description of the simulation, its geometry, and some relevant parameters. In Sec. III, we provide the definitions of the measured fluxes and the methodology of extracting them from the XGC1 data. We present the equilibrium and turbulent E × B particle and heat fluxes across the separatrix and interpret the observed patterns using the notions of ion X-losses, ion banana drifts, non-adiabaticity, quasineutrality, and shear flows. In Sec. IV, we perform numerical integrations to assess the size of the different fluxes and extricate the contribution of each to the cross-separatrix transport and hence, to the SOL width. In Sec. V, we summarize and draw the conclusions.

In this section, we are going to provide a description of the simulation domain and parameters. The simulation is of the DIII-D⁴⁰ discharge No. 144981, which is an H-mode heated by neutral beam injection. It is initialized with the experimental profiles taken at time 3175 ms. The simulation inputs include experimental kinetic profiles of electron density and temperature (n_e and T_e), ion temperature (T_i), and magnetic equilibrium, from a kinetic EFIT magnetic reconstruction. In Fig. 1, we see a poloidal cross-section of the simulation domain along with the directions of the unit vectors of the cylindrical coordinate system (R, ζ, Z) used in this paper. The origin of the coordinate system in the poloidal cross-section will be taken to be the point (R, Z) = (R_o, Z_o) = (1.67 m, 0). The geometry is lower single null but with a secondary “virtual” upper X-point (outside of the simulation domain). The magnetic field is in the negative $ζ̂$ direction making the ions ∇B-drift towards the X-point and the electrons towards the top. The black dashed line denotes the separatrix. The yellow lines are closed and open flux surfaces, and the black and red continuous lines are showing the poloidal projection of the banana orbits (black for ascending-red for descending parts of the orbit). These have been plotted using contours of toroidal angular momentum. The closed dashed black line close to the center is a potato orbit. The orientation of the coordinate system and the toroidal magnetic field direction are also shown. In the rest of the paper, when we refer to the poloidal angle θ, of a separatrix point, we will mean the angle formed between the line passing through that point and the line passing through the low-field side midplane with both lines passing through the origin. Positive angles are in the counter-clockwise direction.

The total simulation time is 0.16 ms with a time-step of 2.3 × 10⁻⁴ms. Although this time is not enough for the core turbulence to saturate, it is adequate for the saturation of the edge turbulence. The poloidal magnetic field is B_pol = 0.42 T (measured at the outboard midplane, in the positive $êθ$ direction) and the toroidal current is I_p = 1.5 MA. Equilibrium plasma density at the edge was measured to be n ≈ 3.5× 10¹⁹ m⁻³ with a Greenwald density of n_G ≈ 15.6 × 10¹⁹ m⁻³. Sound speed and poloidal flow were estimated from the data to be $Cs≈15×104 ms$ and $Vθ≈4×104 ms$⁠, respectively. Safety factor close to the edge was found to be q₉₅ = 3.7 giving a connection length of about qR ≈ 8.6 m. This connection length, along with an electron thermal velocity of $ute≈5.1×106 ms$ results in a transit frequency of 593 kHz. The collision times were calculated at τ_ei = 3.5 μs and τ_ii = 0.5 ms. The ions and electrons are weakly collisional; thus, the adiabatic invariant μ is conserved over many bounce times.

In Sec. III, the trapped particle orbits are going to play a fundamental role in our interpretation of the flux patterns; here, we give a few relevant numerical estimates regarding them. The fraction of trapped particles at the edge was estimated from the relation $ϵ(1.46−0.46ϵ)≈0.767$ with ϵ being the aspect ratio. The bounce times of deeply trapped particles were calculated to be τ_b,i ≈ 0.16 ms and τ_b,e = 6 μs. Technically, the bounce time of a deeply trapped particle whose banana orbit turning point located at the high field side (HFS) midplane is infinite. The above bounce times are simple estimates found by dividing the length of a poloidal circuit with the magnetic drift speed of each particle and they are presented here to make the point that almost all of the ions have enough time, within the simulation, to complete a banana orbit. With ϵ^3∕2 ≈ 0.2, we calculate $ν⋆i=νiiq95Ruti≈0.12$ and $ν⋆e=νeiq95Rute≈0.48$⁠, we can deduce that at the edge, the electrons are in the plateau regime (ϵ^3∕2 < ν_⋆e < 1), whereas the ions are in the banana transport regime (ν_⋆i < ϵ^3∕2).

One of the goals of this paper is to explore the relationship between equilibrium and turbulent processes and the associated fluxes across the separatrix. For dynamical quantities, such as the plasma density n or velocity v, we define

where $⟨⋯⟩$ denotes an average in either time (t) or toroidal planes (ζ) or both. The time averages considered here are taken over a time interval late in the simulation where a quasi-steady turbulent state has been achieved.

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

where the cross terms $⟨nδv⟩t,ζ$ and $⟨vδn⟩t,ζ$ vanish due to the vanishing of $⟨δn⟩t,ζ$ and $⟨δv⟩t,ζ$⁠.

These “fluxes” are not technically transport fluxes but rather local density weighted flows. Although collisions are fully accounted for in the simulation, the rapid electron transit time implies that the radial excursions of electrons due to either magnetic or E × B drifts are small and can cancel without causing net transport as is well known from neoclassical theory. In this paper, we will continue to refer to them as fluxes for brevity. As will be seen, they provide a useful diagnostic for understanding basic properties of both the background plasma and the resulting turbulence. Often in the literature, we find the first piece of the rhs of Eq. (2) to be called the equilibrium flux, whereas the second piece is known as the turbulent flux. In this work, we employ the same definition for the turbulent flux; however, for convenience in the analysis, the definition of the equilibrium flux employed in this paper is slightly modified. Our equilibrium flux is given by $Γeq=⟨⟨n⟩ζ⟨v⟩ζ⟩t$⁠. These two definitions of the equilibrium flux are not exactly equivalent. They are connected by the relation: $⟨⟨n⟩ζ⟨v⟩ζ⟩t=⟨n⟩t,ζ⟨v⟩t,ζ−⟨⟨δn⟩ζ⟨δv⟩ζ⟩t$⁠. However, independent calculation of the difference $⟨⟨δn⟩ζ⟨δv⟩ζ⟩t$ has shown that this “ringing” term, associated with temporal fluctuations of the toroidally averaged fields, is three orders of magnitude smaller than the two equilibrium fluxes indicating that, for the level of numerical accuracy assumed in this work, we can take the two definitions to be equivalent.

The definitions for the heat fluxes are similar, with the replacement of the density by the pressure of each species. Because the simulation is quasineutral, we cannot distinguish between the densities or the net particle fluxes of the two species. In general, both magnetic and E × B drifts contribute to the local particle fluxes crossing the separatrix at a particular poloidal position. In principle, the electron and ion E × B drift fluxes can be different due to the larger ion gyroradius. This orbit averaging effect results in the suppression of the ion E × B fluxes by a factor $Γo(b=k⊥2ρi2)=e−bIo(b)$⁠. Close to the separatrix, ρ_i ≈ 2 mm and from power spectra of equilibrium and turbulent fields, $k⊥eq≈61 m−1, k⊥tur≈123 m−1$ which gives a suppression of about 1%–2% for the equilibrium and about 5%–6% for the turbulent E × B fluxes. In the case of the heat fluxes though, a clear distinction between ions and electrons can be made, based on their different temperatures. Finally, we should note that, given the available diagnostics, we lack the capability of calculating the turbulent polarization current and its flux. Nevertheless, in later sections, we estimate its magnitude relative to the rest of the fluxes.

In this section, we describe our method for calculating the various fluxes, using the XGC1 data, which illustrates the limitations of the particular data set and can be useful to practitioners in the field.

The calculation of radial E × B drift velocities from the density and potential fields, whose values we know on the computational nodes, involves derivatives along the R and Z directions. However, a plot of time and toroidally averaged n_e and Φ on the separatrix as the one in Fig. 3 shows that, despite the fact that a characteristic trend can be immediately identified, the values are noisy, especially at the low field side (LFS). Therefore, the strategy of interpolating the computational node values along horizontal and vertical lines and then taking derivatives along those lines to construct the vector velocity can give highly inaccurate results particularly for a sub-dominant component, such as the radial component of the E × B velocity. For this reason, we resorted to expressing the required derivatives as derivatives with respect to the poloidal angle θ. This can be accomplished by writing the magnetic field in terms of the coordinate system (Ψ, ζ, θ). In this coordinate system, the magnetic field is B = RB_ζ∇ζ + ∇Ψ × ∇ζ. Using the formula $J−1=∇Ψ·∇θ×∇ζ$ for the inverse Jacobian, we can turn the expression for the radial component of the E × B velocity, $VEΨ=êΨ·b̂×∇ΦB$⁠, into

For all equilibrium quantities, $∂Φeq∂ζ=0$⁠, since the system has axisymmetry. For turbulent quantities, we can compute such derivatives using the assumption that turbulent structures are field aligned. Demanding B ⋅∇δΦ ≈ 0, we arrive at the equation

which can be readily used. The above expressions allow us to work with derivatives in the poloidal direction for the calculation of the E × B velocity. This procedure has the advantage that we only need to do a single interpolation and smoothing of n_e and Φ around the LCFS and then take the derivative. The results thus obtained are much less noisy.

First, we will describe the fluxes qualitatively and attempt to give an intuitive explanation for their shapes. Afterwards, we will seek to arrive at relevant quantitative conclusions based on their relative sizes.

We start with the spatial distribution of the equilibrium, E × B particle flux normal to, and hence crossing, the separatrix, which is common between the two species, as discussed above.

In Fig. 2, we observe that the magnitude of the particle flux (arrow length and color indicate flux intensity) is not uniform and exhibits a strong poloidal variation, not all of which can be explained by the difference in field strength between the inboard and outboard sides. Certainly, the flux is suppressed at the HFS compared to the LFS. However, at the LFS, the flux seems to be concentrated at two lobes, one outward, starting near the X-point and extending up to a little below the midplane, and one inward, which starts above the midplane and extends near the top of the machine. At the top, near the top X-point, the flux becomes again slightly outward, before it turns again to slightly inward at the HFS all the way to the bottom X-point.

To understand the poloidal pattern we just described, we turn to the distributions of density and potential on the LCFS, as shown in Fig. 3. There, we plot the two terms of the adiabatic relation so that we get a measure of the electron departure from adiabaticity at various locations along the LCFS. Moreover, the $ϕTe$ term has the same structure as the potential itself, given that the electrons are almost isothermal on the flux surfaces. From this plot, we immediately recognize the points where the global minimum (ϕ ≈ 120°) and global maximum (ϕ ≈ 270°) are located, which are the top and bottom X-points, respectively. Although the exact location of these two extreme points is a result of a complex interplay between Pfirsch-Schlütter flows and collisional effects,^24,42 it follows from the direction of the magnetic drifts: ions drift downwards and electrons upwards; therefore, the potential arranges itself so that it attracts ions from the core at the top and electrons at the bottom. In both cases, the extremization of the potential follows an analogous increase/decrease in the density. In the case of the bottom X-point, an important factor that would lead to a density build-up is the X-point ion loss. The X-point loss^38,39 of hot ions leads to a lower ion temperature at the LCFS. It should be noted that, in our simulation, the E_r well from the X-loss has already been established. Therefore, the loss cone has moved to higher energies, resulting in hot ion loss. Although the scalar ion pressure is not constant along a field line,⁴³ approximate pressure balance is consistent with the lower temperature to be offset by a density increase. These two extremal points of potential establish an electric field in the HFS region. It is this electric field which in turn creates the E × B flow which is responsible for the inward flux pinch shown in Fig. 2 on the HFS as well as the dominant outward flux on the LFS. The reason why this flux is not as high as the ones we see on the LFS, even if the electric field is stronger is the mitigating effect of the strong magnetic field.

In addition to the global minimum and the global maximum of ϕ, we also have two local ones. Their positions correlate with the two bulges of opposite direction as shown in the LFS in Fig. 2. Indeed, the directions of the fluxes are completely in line with the electric field directions established by the alternating pairs of minimums and maximums of ϕ. However, in contrast to the two global extremal points of potential, it is not immediately clear what is responsible for the two local ones in the LFS. Here, we propose that this potential structure is produced by the orbits of trapped ions which, as shown in Sec. II, constitute roughly 76% of the total ion population near the edge. A trapped ion moving on its banana orbit exits the separatrix at some point below the LFS midplane, leaving behind it a negative charge. If the electrons were completely adiabatic, negative charge would be instantly neutralized by electrons moving along the field lines to prevent the potential build-up. But, as we can readily verify from Fig. 3, there is a significant departure from adiabaticity at the LFS. Therefore, a trapped ion leaving the separatrix will create a negative potential in order to attract more ions. Conversely, the region above the LFS midplane is the region where the trapped ions re-enter the LCFS (at least those that did not get entrained in the parallel SOL flow) generating a positive charge. The non-adiabaticity of the electrons forces the potential to increase, in order to attract them at that location. Thus, the neoclassical ion orbits combined with X-point losses are the dominant mechanism responsible for the structure of the poloidal electric field, and hence the equilibrium radial particle flux in this simulation.

The non-adiabaticity of the electrons at the separatrix seems to be crucial for the exact shape of the edge potential. This non-adiabaticity at the edge can be understood if we compare the time rates of all terms in the electron momentum equation. More specifically, the inertial term, i.e., the transit frequency, was given in Sec. II as 593 kHz. Although large, it is comparable to the electron collision frequency (also given in Sec. II) and the turbulent frequency, which is discussed in Sec. III C 3 and it is found to be 640 kHz. Therefore, both collisions and the inertial term are important in the electron momentum equation, leading to the observed departure from adiabaticity.

The subject of neoclassical direct ion orbit losses and their contribution to the edge flows and electric fields has been treated substantially, both computationally and theoretically.^44–53 Here, we offer a new argument to support our claim that the potential structure at the LFS is due to the excursions of trapped ions on their banana orbits. Assuming that the trapped particles, with orbits within a banana width from the edge, will create either a charge excess or a charge hole at the point of inflection (“banana tip”), we give an estimate for the angular distribution of those banana tips at the LFS.

To get the fraction of trapped particles with banana tips in the range (0, θ), where θ is the angle that the magnetic field is B and zero is the angle at B_min (LFS midplane), we perform the following integral:

where we have changed coordinates from velocity space to the E, μ space, with E being the energy, μ being the adiabatic invariant. Here, σ denotes the sign of the particles parallel velocity.

Plugging in the angular formula for B, $B=Bo1+ϵ cos θ$⁠, valid for circular flux surfaces, we turn this into the cumulative angular distribution for banana tips

Interpreting X(θ) as $X(θ)=∫0θdθ g(θ)$⁠, we get $g(θ)=dX(θ)dθ$⁠. Therefore,

The poloidal distribution of the flux caused by these particle excursions is the product of the tip distribution with the actual magnetic flux. In Fig. 4, we can see this flux compared to the E × B fluxes. We point out that indeed, the “trapped” flux tracks qualitatively the equilibrium E × B flux. Two things should be noted: the modelling of the tip distribution is rather simple as it assumes circular flux surfaces and completely ignores the effect of the X-point and the ion losses there. Also, a qualitative similarity between the fluxes is all we can hope for. Indeed, if the equilibrium E × B flux pattern at the LFS is the result of the potential structure due to the trapped particle excursions, as we claim, then, the relationship between their relative sizes is far from obvious.

We proceed with the description of the turbulent E × B flux, as shown in Figs. 4 and 5. There, we observe that particle turbulence is entirely confined to the LFS and has a ballooning-like shape which peaks near the midplane, gets interrupted above it up to roughly 50°, picks up again after that, and recedes off at the top of the machine.

A probable explanation for such a shape is provided by Fig. 6. There, we overlay the plots of the turbulent flux and the E × B shearing rate on the LCFS. We observe that as the shearing rate exceeds a certain value (roughly at 640 kHz) the turbulent flux drops sharply. When Ω_E×B drops below this value, the turbulence rises again. Indeed, we have measured the main frequency of the turbulence in this discharge to be very close to 640 kHz. We can imagine that, therefore, in the absence of shear, the turbulence would have created a ballooning structure in the LFS, peaking at the midplane and tapering off at the top, where the curvature drive is minimal. However, the presence of shear at exactly the right level suppresses it locally creating the resulting pattern.^54–57 

Recall from the previous discussion that the electrostatic potential pattern, and hence the shear, appears to be controlled by neoclassical drift-orbit effects. It is notable that the system has arranged itself so that this equilibrium shear is of sufficient magnitude to influence the turbulence. We have verified that turbulence generated Reynolds stresses play a negligible role in establishing these sheared flows, which, while not unexpected for this H-mode phase of the discharge,¹³ has led us to consider alternative turbulence saturation mechanisms. More specifically, we examined wave breaking (which requires the perturbed radial flow to be larger than the wave phase velocity) and the gradient removal mechanism (where the scale length of the perturbation is comparable to the equilibrium one, leading to profile modification). Order of magnitude estimates show that the size of the perturbations in our simulation is not strong enough to make either mechanism a plausible candidate for turbulence saturation. Even though the above mechanisms seem to be too weak by themselves to provide nonlinear saturation, it is possible that gradient modification could act in conjunction with mean flow shear. Since this equilibrium shear is large enough to impact the wave dynamics, it could suppress the linear growth rate to a low level. In a marginally unstable state, relatively weak nonlinear processes could provide saturation. A study of the linear modes, driving gradients, and the impact of mean flows is the subject of ongoing investigations. One other mechanism which remains to be investigated is mode coupling to damped modes.⁵⁸ 

In addition to shear, it is possible that the observed turbulence pattern may also be related to the existence of curvature driven modes which peak at the top of torus.^59,60 Poloidal flux surface expansion, i.e., the distance between adjacent flux surfaces which is increasing towards the top of the machine, may also play a role. First, it increases the outward velocity and hence the Γ of filamentary turbulent structures in order to maintain their field-aligned character.⁶¹ Second, flux surface expansion reduces E × B shear more quickly than the density or temperature gradients, since the former are proportional to a second radial derivative of the potential and therefore scale more strongly than the gradients which are proportional to a single radial derivative of the profiles.

Finally, given the turbulent frequency, we can evaluate the relative size of the ion polarization flux compared to the turbulent E × B flux. Their ratio is approximately $ωturωci≈0.04$⁠, which means that our inability to calculate the polarization current does not threaten the credibility of our results, at the present level of precision.

In this section, we present the qualitative features of the heat fluxes. Contrary to the E × B particle fluxes where we could not distinguish between the two species (apart from the, almost negligible, gyroaveraging effect of the ion orbits), for the heat fluxes we have a clear separation due to the different temperatures, as shown in Fig. 7. Therefore, in what follows, the calculation of the fluxes still assumes quasineutrality but each of the species has its own temperature. In calculating the heat fluxes, we employ Eq. (2) where in the place of velocity we substitute the pressure P_s = nT_s, with s being denoting the species. Therefore, the heat flux is given by the formula

where again, we call the first term of the rhs as the “equilibrium” and the second as the “turbulent” heat flux. The velocity here is the E × B one.

In Figs. 8 and 9, we present the equilibrium and turbulent E × B fluxes of electrons and ions, respectively, along the LCFS. The shapes of both species equilibrium flux are similar except for the scaling effect of the temperature difference. This can be easily understood from Fig. 7 where we plot the temperatures as a function of the poloidal angle. There, we see that the electrons are very close to being isothermal. Therefore, the qualitative features between the particle equilibrium and electron equilibrium E × B heat fluxes are common.

The ion temperature, on the other hand, does not remain constant along the separatrix. Specifically, it drops near the two X-points indicating that, close to these areas, hotter ions have exited the plasma. High energy ions with the right pitch angle will be in the ion loss cone.^38,39 An interesting feature of the ion temperature is that the average temperatures below and above the outboard midplane are not the same: the region below the midplane has an average temperature which is 4.5 eV higher than the average temperature above. We can understand this difference in the context of the trapped ion banana orbits that we talked about in the particle fluxes section. The ions exiting the flux surface below the midplane are on the inward side of their banana orbit. They travel along a region of higher temperature; therefore, the temperature locally increases. Equivalently, when they re-enter above the midplane, they are coming from outside the LCFS where the temperature has dropped. Indicatively, a typical value for a banana orbit width near the edge is 6 mm and the temperature drop across this distance is about 2 eV, accounting in part for the observed average temperature difference. Moreover, the re-entering ions could have originated from further away in the SOL than a typical banana orbit width, making the above estimate just a lower bound for the real temperature difference. Despite the above observations, the qualitative features of the equilibrium E × B ion heat flux are very similar to those of the equilibrium particle and electron heat E × B fluxes.

The real difference between the two species heat flux behavior comes from the comparison between the relative sizes of equilibrium and turbulent heat fluxes. In the case of the electrons, turbulence dominates; neoclassical electron transport is negligible because of the small electron banana width. Thus, the local equilibrium E × B electron flux shown in Fig. 8 will be further reduced by cancellation from bounce-orbit averaging. Turbulence is the sole mechanism for net electron particle loss across the separatrix. It must compensate for the net loss of ions in order to maintain charge balance. Turbulent electron energy loss follows.

In contrast, for the ions, it seems that the equilibrium process, set up by the neoclassical mechanism of trapped particle orbits, is indeed the dominant mode of ion loss, in this simulation, both for heat and particles. Both species LFS turbulent ballooning shapes are being interrupted in the high shear region. Although the ion turbulent heat flux never recovers and drops to zero, the electron turbulent heat flux picks up again after the initial dip. Moreover, in the region above the midplane, the electron turbulent heat flux is almost twice the size of the equivalent ion one, despite the big difference of the absolute temperature of the two species, with the ions being significantly hotter. This is a clear indication that, in the case of electrons, turbulence is absolutely necessary for the plasma to retain charge balance, whereas in the case of ions, the neoclassical processes are more effective.

In Figs. 10 and 11, we present the effective diffusivities which are defined as $D=−Γ∂n/∂x, χs=−Qs∂Ps/∂x$ and we separate them again into equilibrium and turbulent contributions. The equilibrium diffusivities are presumably driven by neoclassical ion physics as discussed previously, but because of ambipolarity, equilibrium and turbulent diffusivities are fundamentally coupled. A first observation is that both equilibrium and turbulent diffusivities are of a similar order of magnitude, consistent with this coupling.

In Ref. 62, the methods developed for application to pedestal turbulence provide information about the modes present, based on the relative sizes of the transport channels. It is interesting to apply these concepts to our situation for qualitative insight. Here, we find comparable turbulent diffusivities near the midplane in particle, electron heat, and ion heat channels. This is suggestive of an underlying mode with an MHD character, possibly resistive ballooning or Kelvin Helmholtz in this electrostatic simulation. Near the top of the torus, the electron diffusivity is dominant, which may suggest a mode driven primarily by the electron temperature gradient or a mode with an electron drift character. Indeed, we find our turbulence has a frequency very close to the electron diamagnetic drift frequency.

A final observation is that, based on the scale of the above diffusivities, the time scale for further profile evolution is much larger than the total time of the simulation. We have indeed verified that the profiles remain constant throughout, which is a clear indication of achievement of a quasi-steady turbulent state.

Here, we calculate the integrated currents and flux surface averages of the above fluxes and reach some conclusions regarding their relative importance in establishing the SOL width.

First, we start with some definitions. The fluxes of Sec. III are normal to the separatrix components of the total flux, i.e., $ΓN=−BRBp⟨n·VZE×B⟩+BZBp⟨n·VRE×B⟩$⁠. For a quantity A, the flux surface average of it is defined as $⟨A⟩=∫dζdθJA∫dζdθJ$⁠.

For the calculation of such integrals, we must know the Jacobian factor the inverse of which is given by the well-known formula

with r being the distance between (R_o, Z_o) and the point on the flux surface. In the second step of the above calculation, we have replaced the usual formula $∇θ=1rêθ$⁠, which holds in a circular geometry, with $∇θ=1r1+1r2d2rdθ2êθ$⁠, which is an approximation for a general shape flux surface, provided the integrand is non-singular and single valued, which minimizes the integration error. Using this formula, the flux surface average (FSA) of a scalar A is given by

If we perform a flux-surface average on the continuity equation, we arrive at the expression

where $V′=∫dζdθJ$ and $∇Ψ=RBpêΨ$⁠. Therefore, the physically relevant quantity is $⟨Γ·∇Ψ⟩$⁠, which is calculated using the explicit formula

We can also define the surface integral of the flux across the LCFS. This should be interpreted as the particle or heat current, i.e., the number of particles per second that cross the separatrix due to the various fluxes. This will be denoted as I_Γ and calculated using

with similar expressions for the heat currents.

In Table I, we present the integrals over the whole separatrix. We have also kept for reference the numbers for flux from the magnetic drifts. In Tables II and III, we give the results for the integrated heat fluxes.

Here, we also define an averaging factor, AF, which is defined as $AF=⟨Γ⟩/max(Γ)$ and shows the degree of “cancellation” of each particular flux. If we assume that the density and potential are constant on a flux surface, then the FSAs of both the equilibrium E × B and the magnetic fluxes can be shown to be identically zero. Close to the edge of the plasma, as shown in Sec. III, the constancy of those quantities on the separatrix is violated; hence, we do not expect the AF to be exactly zero. Indeed, as shown in Table IV, the averaging factor for both the magnetic drift flux and the equilibrium E × B flux are very small but nonzero, indicating that the bulk of these fluxes are returning into the main plasma, whereas, the AF of the turbulent E × B flux is larger. AFs of the heat fluxes are similar.

A few comments about the relative sizes of those integrals are in order: In the particle fluxes and ion heat flux case, the size of the turbulent current is roughly three times larger than the equilibrium one. In the case of the electron heat fluxes though, the turbulent current is almost ten times larger than the equilibrium one indicating the importance of turbulence for the electrons, supporting the main argument of Sec. III.

In this subsection, we explore particle balance in the SOL by comparing the flux of electrons exiting across the separatrix with the exhaust flux in the SOL, which is dominated by the parallel flows. The competition between these cross-field and parallel transport processes sets the density SOL width.

Regarding the separate contribution of each flux to the setting of the SOL width, the ratio of the currents does not tell the full story. Despite the fact that the turbulent flux is fully localized at the outboard midplane, contributing solely to the LFS SOL, the equilibrium fluxes have an important inward component at the HFS. Moreover, as we have mentioned in Sec. III, the particle E × B equilibrium flux cannot take electrons outside of the separatrix. Therefore, we start with the assumption that the current from the turbulent E × B particle flux is the total electron current exiting the separatrix. Now, we take a Gaussian box whose one side is the part of LCFS starting at the point where the E × B turbulent particle flux begins and ending where it ends. The top and bottom of this box are horizontal lines, starting at these two locations and extending outwards, deep inside the SOL. The extent of these lines is left to be determined by the behavior of the normal fluxes on them. The right hand side of the box is unspecified, but we assume that no flux leaks from this side of the box to the wall.

In Figs. 12 and 13, we have plotted the normal component of the particle fluxes at the bottom and top caps of the Gaussian box, respectively. The normal component of the particle flux includes contributions from the E × B flux, the magnetic drift flux and, importantly, the parallel flows. We are not concerned with distinguishing the contributions from equilibrium and turbulent E × B; therefore, we just present the total. An important note here is that, since we do not distinguish between equilibrium and turbulence, in the calculation of the E × B velocity, we take $∂ϕ∂ζ=0$⁠. We believe that the error introduced by such an approximation is minimal, since the E × B contribution to the normal fluxes is only important in a very narrow region close to the separatrix. There, the poloidal magnetic field is almost zero, hence $∂ϕ∂ζ≈B·∇ϕ=0$⁠. Generally, the parallel flux contribution is dominant everywhere and the only place where the other fluxes seem to be relevant is at a very narrow channel close to the LCFS. In both figures, the negative Z direction is downwards, which means that for the top cap, negative fluxes are entering the box, whereas for the bottom cap negative fluxes are leaving.

In the case of the top cap, it seems that there is a well defined point where all fluxes drop to zero, around Ψ_N ≈ 1.16. This seems to be a reasonable place to truncate the integration of the fluxes in order to find the total current crossing the top of the box. The total current thus found is 8.61 × 10²¹ s⁻¹ = 1.38 kA, entering the box.

For the bottom cap, things are not so clear cut. Although the fluxes coming from E × B and magnetic drifts do indeed drop to zero, the parallel flux seems to remain constant. This had to be expected, since there is a significant degree of ionization happening in the SOL. It is very hard to precisely locate a point to truncate the integration, meaning to find where the SOL ends and what is left after that is simply ionization. We believe that a sensible point to stop our integration is at Ψ_N ≈ 1.025 (shown in Fig. 12). The total current found from this choice is 1.17 × 10²² s⁻¹ = 1.87 kA, exiting the box. If we extend the integration further out, this number will increase.

As we said above, the total electron current crossing into the box from the separatrix has to be the one coming from the turbulent E × B flux given in Table I as 3.02 × 10²¹ s⁻¹ = 0.48 kA. If we compare this number to the total current found from the top and bottom of the boxes, we see that it accounts for 98%. A direct conclusion that can be drawn from such a comparison is that indeed, as presumed in the Goldston model, the density channel is filled up almost entirely with turbulent electron flux. Of course, we had taken the limit of integration of the bottom cap further out, this percentage would have dropped; however, it is reasonable to assume that this result would be severely contaminated by atomic processes that take place outside the LCFS and have nothing to do with turbulence and flows that cross the separatrix.

Here, we need to note that even though we can safely claim that the electron contribution to the SOL channel is almost entirely turbulent, a similar statement cannot be made for the ions. From the data available, we have no way of isolating the turbulent and the neoclassical contributions, since turbulence, in the ions' case, exists on top of a neoclassical background.

If we repeat the previous calculation for electron heat fluxes, we find both the top and bottom caps of the box to have well defined limits where all fluxes go to zero. However, the comparison of the total electron heat current exiting the Gaussian box from the caps to the electron heat current entering from the separatrix shows that the latter is more than twice the former. A plausible explanation for this result is that heat is dissipated in this volume due to various atomic processes, predominantly ionization. A detailed account of the relative strengths of such processes is beyond the scope of this paper.

In this paper, we have analyzed the results of an XGC1 DIII-D simulation focusing on the fluxes exiting the separatrix. We reconstructed the E × B flux distinguishing between the equilibrium and the turbulent part. The poloidal patterns of the two particle fluxes were presented. For the equilibrium flux, we traced its poloidal shape to the LCFS potential. We interpreted this potential structure as the result of ion drifts. Specifically, we could clearly distinguish the impact of the ∇B-drifts and X-point losses on the potential's global maximum and minimum at the bottom X-point and top of the machine, respectively. For the poloidal potential variation at the LFS, we argued that it is due to the trapped ion orbits close to the edge that exit and re-enter the confinement. We constructed a simplified model for the angular distribution function of the banana orbits' inflection points and showed that the resulting “trapped-particle” flux follows the shape of the equilibrium flux, to support our argument. Moreover, we noted that LFS edge electrons exhibit a significant departure from adiabaticity. We provided as reasons for this departure the fact that the turbulent frequency and the collision rate are comparable to the transit frequency. The turbulent flux pattern was also given. It appears to be strongly ballooning but with a significant region of suppression. We contended that this suppression is coming from the ion-neoclassically generated E × B shear which seems to be at the right level to affect the turbulence exactly at that region.

We presented the E × B heat fluxes for ions and electrons at the separatrix. From the difference in size between the turbulent heat fluxes of the two species, it is clear that turbulence is more important for the electrons, since it is their only way of losing confinement. Ions are exiting the separatrix through neoclassical mechanisms; therefore, turbulent electron heat flux has to increase in order to maintain quasineutrality. The electrons were found to be isothermal whereas the ion temperature exhibits significant drops near the X-points, in accordance with the predictions of X-loss theory. The small temperature difference between the upper and lower halves of the LFS can be, in part, explained by ion banana orbits exiting from a hotter and re-entering from a cooler part of the machine. The equilibrium and turbulent diffusivities of those fluxes were calculated, providing insight regarding the present modes. Numerical integrations were performed to calculate the flux-surface-average fluxes and the outgoing particle and heat currents from the separatrix. Finally, we confirmed that anomalous electron diffusion accounts for almost all of the electrons filling the SOL particle channel.

Given the importance of the predictability of the SOL width for modern tokamaks, more state-of-the-art gyrokinetic simulations need to be analyzed to elucidate the complex physical processes that give rise to it. The interaction between the ion drifts and the turbulence seems to be critical for setting the size of the channel and the study of the fluxes and flows at the edge provides illuminating insights into the nature of the system dynamics. In future publications, we will examine simulations from different tokamaks and look into the characteristics of the turbulent “blobs” from this and other simulations. These efforts will be directed towards obtaining a better understanding of the interaction of neoclassical and turbulent physics and the extent to which the present results are applicable to other plasma and device regimes.

I. Keramidas Charidakos would like to acknowledge the useful help with the data from Jugal Chowdhury at the early phase of this project. We acknowledge computing resources on Titan at OLCF through the 2015 INCITE and the 2016 ALCC awards. This work was supported by the U.S. Department of Energy Office of Science, Office of Fusion Energy Sciences under Award No. DE-FG02-97ER54392 and by subcontract SO15882-C with PPPL under the U.S. Department of Energy HBPS SciDAC project.