Reduction of blob-filament radial propagation by parallel variation of flows: Analysis of a gyrokinetic simulation

Blob-filaments^1–3 (also called blobs, filaments, intermittent plasma objects, and coherent structures) have a long history of research in plasma physics, dating from observations in the 1980s⁴ to present magnetic fusion devices.^5–10 These coherent structures are of interest for fusion research in that they, or the turbulence from which they emerge, may influence the near-SOL (scrape-off layer) width,^11–17 as well as mid to far-SOL profiles^18–21 and main chamber wall interactions.

The importance of, and intrinsic interest in, blob-filaments has spawned considerable theoretical efforts to understand blob radial velocities.^2,3,22–25 Originally carried out in model two-dimensional geometries using fluid theory,^26–29 the models were later extended to three dimensions,^30–32 and a number of fluid simulation codes now also incorporate increasingly realistic tokamak geometry.^33–36 Recently, kinetic simulation codes have also addressed blob-related research.^37–40 In the present paper, we will continue this thrust by exploring results from the XCG1 code,^41,42 a full-f, 5D particle-in-cell (PIC) gyrokinetic code.

The goal of this paper is to analyze and interpret XGC1 simulation results impacting blob dynamics. In particular, we will explore the effect of toroidal flows that vary in the poloidal (parallel) direction on the parallel structure of blob-filaments and the consequences of this imposed structure for radial blob propagation. Thus, the sheared flows considered here are different from the radial variation of poloidal or toroidal flows that have usually been considered in the literature. However, in a broader sense, our work is at least conceptually related to a very large body of work on the effect of sheared flows on plasma turbulence.^43–45 More specifically, it has been observed experimentally that blob-filaments are distorted and sometimes even split, by the sheared flows that arise from radio frequency induced convective cells in the SOL.⁴⁶ Recently, the effects of sheared poloidal flows on disconnection of blob-filaments from the divertor have been studied.⁴⁷ 

The plan of our paper is as follows: In Sec. II, an overview is presented of the gyrokinetic simulation on which the subsequent analysis is based. This section discusses the magnetic geometry, bulk plasma parameters, and fluctuation statistics. A field line following coordinate system, essential to the filament analysis, is introduced in Sec. III A, followed, in Sec. III B, by the simulation results for the blob-filament structure and motion.

The theory of charge shorting by parallel electron motion is discussed in Sec. IV. Charge shorting occurs when curvature and grad-B induced polarization charges, responsible for blob-filament radial convection, are mitigated by the parallel electron conduction. A heuristic condition for shorting is given for a collisionless system in Sec. IV A, and this result is compared to and reconciled with a similar condition that applies in a collisional fluid model in Sec. IV B.

In Sec. V, the theory is applied to the blob-filaments observed in the simulation. It is shown that they are in the collisionless regime and that they satisfy the condition required for charge shorting. Furthermore, the density and potential fluctuations in the simulation are shown to be nearly in phase, consistent with the dominant role of electron parallel conductivity. Several diagnostic methods were applied to the simulation to measure the radial velocity of the blob-filaments. The blob-filament velocity was found to be smaller than that expected from an ideal analytical estimate that ignores the charge shorting mechanism.

A discussion of some related points is given in Sec. VI, and a summary of this paper and the main conclusions are given in Sec. VII. Three additional topics are discussed in the Appendixes.  Appendix A provides details of the time-lag correlation method used to measure blob-filament motion,  Appendix B discusses the blob-velocity proxy method, and  Appendix C discusses various regimes of the two-region blob model with respect to the simulation.

In our paper, the terms blob-filament, blob, and filament will be used almost interchangeably, with the choice depending on whether the particular emphasis of the discussion is on the perpendicular (blob) or parallel (to B, filament) structure.

In this paper, we analyze an XGC1 simulation of the lower single null Alcator C-Mod EDA high (H) mode discharge #1160930033 with a plasma current of I_p = 1.4 MA, a toroidal magnetic field of 5.4 T, an input power of 5.4 MW, and a record volume averaged pressure exceeding 0.2 MPa.⁴⁸ The geometry is shown in Fig. 1. In addition to the active lower X-point, there is a “virtual” upper X-point, i.e., an X-point outside the analyzed simulation domain whose effect on the near-separatrix geometry is still important. The ion grad-B drift is directed downwards, toward the (lower) X-point.

The simulation was initialized with experimental profiles of density and temperature and reconstructed (kinetic EFIT) magnetic equilibrium. This electrostatic simulation employs drift kinetic electrons and gyrokinetic ions. It implements a logical sheath boundary condition⁴⁹ at the divertor plate (which sets the plasma potential such that electrons with lower parallel velocity are reflected in sufficient numbers to maintain ambipolarity). Two neutral-plasma interaction models are available within XGC1, a built-in model that only treats atoms, and an interface to the more sophisticated DEGAS2 Monte Carlo neutral transport code. The former model is used in the present work and is described in detail elsewhere.⁵⁰ Only half of the toroidal extent of the device was simulated; the rest was reconstructed using periodicity so that only even toroidal mode numbers are present in the simulation. The XGC1 code employs an unstructured mesh in the (R, Z) plane with spatially varying resolution, and this simulation used N_ζ  = 16 toroidal, i.e., (R, Z), planes in the half-torus simulation.

Although the simulation encompasses the whole (R, Z) plane of the torus, albeit with reduced resolution in the interior core, our analysis in this paper is restricted to the edge and SOL regions shown in Fig. 1. In the radial direction, this analysis domain is defined by normalized poloidal flux in the range ψ = (0.97, 1.02), where ψ = 1 is the separatrix. This corresponds to a range of distances from the separatrix of ΔR = ±4 mm at the outboard midplane, well short of the distance to the main chamber wall. In Fig. 1, the location of the main chamber first wall is shown for general reference; the study of far SOL blobs and wall interactions is deferred to future work. The first wall at locations with ψ < 1.07 was included in the full simulation. In the region with ψ > 1.07, any turbulent or neoclassical potential perturbations were ignored. Thus, the simulation includes neutral effects in the regions where there are blobs, but the self-consistent blob interaction with the wall is not taken into account because of the radial buffer zone.

In the region of interest for this paper, namely, the near-separatrix region, each of the 16 toroidal planes contained 7.6 × 10⁴ nodes for a total of over 1.2 × 10⁶ nodes. The spatial separation of the nodes is about 0.5 mm, which is comparable to the ion gyro-radius at the separatrix. The 16-plane toroidal resolution per half torus is sufficient because the unstructured mesh in the (R, Z) plane is carefully chosen to place nodes nearly along field lines. Field-line-following (FLF) methods used in the analysis are discussed in Sec. III A. The toroidal planes end up determining the parallel resolution of the filament, which varies slowly, while the perpendicular resolution is set by the node spacing.

The (deuterium) plasma conditions in the simulation are roughly characterized by the following parameters on the separatrix at the outboard midplane: density n_e = 1.3 × 10²⁰ m⁻³, electron temperature T_e = 127 eV, ion temperature T_i = 203 eV, and magnetic field B = 4.4 T. The plasma profiles decay rapidly into the SOL. The electron pressure gradient scale length in the near SOL is λ_p ∼ 0.8 mm.

The RMS fluctuation amplitudes peak just outside the separatrix at ∼10% for the relative density fluctuations δn/$n¯$ and ∼20% for the relative potential fluctuation eδΦ/T_e. Here, we define fluctuations as departures from the toroidal (ζ) average, e.g.,

Individual fluctuations can easily reach 50% in the SOL, and infrequent stronger fluctuations occur, i.e., intermittency is present as expected.

The intermittency properties are illustrated in Fig. 2, which used data in the entire range ψ = (0.97, 1.02), poloidal angle θ in the range (−1.7, 2.0), and a swath of toroidal angle of ∼0.48 rad in width. (The statistics are expected to be toroidally axisymmetric.) There is nothing particularly surprising in Fig. 2; rather, it confirms that this XGC1 simulation recovers standard statistical results for SOL intermittency such as a skewed and asymmetric probability distribution function of amplitudes in Fig. 2(a), a mostly increasing normalized density fluctuation level as one moves from the closed to the SOL flux surfaces in Fig. 2(b), and a density fluctuation skewness, which is negative on the closed flux surfaces, passes through zero and increases through positive values as one moves in the SOL in Fig. 2(c). The point where the skewness S is zero has been interpreted as the blob birth location: from this location, holes (S < 0) move inward, while blobs (S > 0) move outward in radius.^2,3 Fluctuations in δn/$n¯$ peak just outside the separatrix near a normalized poloidal flux coordinate value of ψ = 1.01 where much of the blob-filament analysis in this paper will be carried out. This flux surface corresponds to ΔR = 1.5 mm at the outboard midplane.

The intermittency illustrated in Fig. 2 is present in coherent structures that emerge from unstable linear modes. These structures, in the form of blob-filaments, are best viewed and analyzed in field line following (FLF) coordinates because of the large disparity in their spatial scales across and along curved field lines. The FLF coordinates used in this paper are the radial coordinate ψ, the parallel coordinate θ, and a toroidal coordinate ζ₀. An individual field line is labeled by two parameters ψ and ζ₀, where ζ₀ specifies the toroidal angle of the field line at a reference poloidal location θ₀′. Here, this reference location is taken as the location of the X-line. Thus, the transformation from geometric angle (ψ′, θ′, ζ′) coordinates to FLF (ψ, θ, ζ₀) coordinates is defined by

where the equation of a field line is R dζ′/ds = B_ζ/B_θ′ and ds is the distance moved in the poloidal direction while moving through dζ′ in the toroidal direction. The position along the field line is denoted by θ with θ = 0 corresponding to the outboard midplane and increasing in the counterclockwise direction in the (R, Z) plane, as shown in Fig. 1. The angle θ is a geometrical angle referenced to the magnetic axis.

Special care is taken in performing the data analysis by using FLF interpolation. First fluctuation data are obtained on each computational mesh node in the (R, Z) plane for each of the N_ζ  = 16 toroidal planes in the half-torus simulation. We denote these toroidal planes by the discrete index i_ζ = (1, N_ζ). The data are interpolated to create N_ζ continuous functions of δn/$n¯$⁠, one for each toroidal plane. The transform to FLF coordinates is, then, carried out to obtain δn/$n¯$ at an arbitrary value of (ψ, θ, ζ₀) using the following method.

We first note that a given field line (ψ, ζ₀) will intersect each toroidal plane i_ζ at a particular value of θ denoted by θ_i. We identify the two adjacent toroidal planes (i, i + 1) defined by θ_i < θ < θ_i+1. Finally, we, then, interpolate along the field line to obtain the fluctuation at the desired point (ψ, θ, ζ₀).

It should be noted parenthetically that a straightforward three-dimensional interpolation of the raw simulation data (i.e., in the form of nodal mesh values of R, Z on each toroidal plane) yielded unsatisfactory results because of the disparity in the perpendicular and parallel structures of the fluctuations.

An example of a blob-filament structure in FLF coordinates is shown in Fig. 3. At a given time, there are many (∼50) such filaments in the SOL of the half torus (0 < ζ₀ < π) of the simulation. In this FLF plane, the magnetic field lines are vertical lines at fixed ζ₀ parameterized by θ. The blob radius (half-width-half-max of δn/$n¯$⁠) in the direction perpendicular to B is about 1.8 mm or about 4ρ_i. In the parallel direction, the filaments extend for about 1 m along the outboard side of the torus from near the top of the device to a location slightly above the X-line. Toroidal flows and the shear in those flows have distorted the blob-filaments and possibly contributed to the X-point disconnection.⁴⁷ The flow can be visualized from the two illustrated frames in Figs. 3(a) and 3(b), with snapshots taken at a time separation of 3.9 μs. (Data were saved every 0.062 μs; thus, there were many other intermediate frames.) The angular velocity may be roughly estimated from these two frames to be about 23 krad/s. This estimate is confirmed and greatly improved upon by the time lag correlation method, as discussed next. Both the time lag correlation method and the blob tracking method discussed in Sec. V were carried out during an interval of 7.9 μs at the end of the 95 μs total simulation time, when the turbulence was in a quasi-steady state.

The time lag correlation method provides a measure of the velocity of structures by searching for two space-time points with maximum correlation. The assigned velocity is, then, v = Δx/Δt, where Δx is the spatial separation of the maximally correlated points and Δt is the time lag. The method is essentially similar to that described in Ref. 51 to analyze experimental gas puff imaging (GPI) data, although some details of the implementation are different. A summary of the method as implemented here is given in  Appendix A. The results of the toroidal angular velocity (ω_ζ0) analysis are shown in Figs. 4 and 5.

Figure 4 shows the radial dependence of ω_ζ0 at two different poloidal locations: the midplane and a point just above the X-line. The angular velocity is roughly constant in the region just inside the separatrix, 0.98 < ψ < 1.00, but deviates strongly from that value in the SOL where it also varies with the poloidal position. The error bars in the figure indicate the standard deviation of multiple measurements including at various ζ₀ locations and over different time intervals (see  Appendix A). In summary, the ω_ζ0 flow is strongly sheared as a result of both radial and poloidal variations in the SOL. This flow shear accounts for the appearance of the sheared filaments in Fig. 3.

The parallel variation of ω_ζ0 is analyzed in more detail in Fig. 5. Again, it can be seen that the flow is more or less uniform in the central region of the filament but becomes strongly sheared at both extremities near the top and bottom of the torus and can even change sign. Part of the reason for the poloidal variation of the flows may be associated with the development of the structure in the electrostatic potential near X-points, e.g., caused by X-point losses.⁵² 

The blobs also move radially, although the radial motion is about two orders of magnitude smaller than the toroidal motion and is, therefore, rather difficult to detect by the correlation method. We will return to the radial motion in Sec. V.

We have seen that spontaneous poloidally sheared flows that develop in the plasma, particularly in the SOL and near the X-points (real and virtual), distort the blob-filaments, causing them to depart slightly from the magnetic field lines. Although the departure is only on the order of the blob perpendicular dimensions, perhaps a few millimeters over a meter or more of field line length, it can play a significant role in the dynamics of the structure. Exploring this is the topic of the present section.

The basic idea is illustrated in Fig. 6, which shows a field line (black) and a filament that departs from the field line. The sketch is in field line following (FLF) coordinates corresponding qualitatively to Fig. 3. Curvature and ∇B drifts charge-polarize the blob-filament in the standard picture of blob dynamics,^1–3 indicated by the red positively charged side of the blob and the blue negatively charged side. If the electron parallel motion along the (black) field line is rapid enough, it will short out the charge polarization. However, the internal charge polarization of the blob is responsible for an internal electric field that propels the blob in the radial direction through its E × B drift. Thus, the shorting mechanism enabled by electron parallel conductivity and the parallel shearing of the filament will suppress the radial motion of the blob. Essentially, the parallel variation of the filament caused by the sheared flows imposes a k_ǁ value on the filament, making it susceptible to electron conductivity. Some different aspects of a similar mechanism in the presence of magnetic shear were also examined in a fluid model in Ref. 53.

The effect of the parallel variation of density in a filament and the subsequent development of a Boltzmann response in the electrostatic potential on the blob-filament dynamics was also investigated in a fluid model.⁵⁴ The resulting mechanisms discussed therein, namely, secondary drift-wave instability of the blob and Boltzmann spinning, are different from the charge shorting mechanism considered here and will be discussed in Sec. VI.

A heuristic condition for the shorting mechanism to be effective may be obtained by comparing some time scales. Three time scales are of interest: the time scales of (i) parallel electron motion across the sheared filament τ_eǁ = d_ǁ/v_te, (ii) curvature drift of charges across the filament τ_κ = δ_b/v_κ, and (iii) perpendicular shearing τ_s = δ_b/v_s. Here, d_ǁ is the parallel distance between the two sides of the filament (see Fig. 1), v_te is the electron thermal velocity, δ_b is the perpendicular blob radius, v_κ is the curvature (or grad-B) drift velocity, and v_s is the shearing velocity of the filament, i.e., τ_s is the time for the filament to distort by δ_b in the frame of its average motion. If τ_eǁ < τ_κ, τ_s, then the charge polarization from curvature and grad-B drifts would be shorted out before it can build up. Conversely, the charge shorting mechanism is ineffective in the limit τ_eǁ > τ_κ. The only other possibility is τ_s < τ_eǁ < τ_κ in which case the filament distorts significantly by sheared flow before the electrons can move. In this case, the shorting mechanism is modified because the field line may or may not connect opposite charges any longer or the filament may break up.

The curvature, or equivalently ∇B, drift velocity in a tokamak is estimated as

independent of species for T_e = T_i = T. Here, c_s = (T_e/m_i)^1/2 is the sound speed and ρ_s = c_s/Ω_i is the ion sound radius, where Ω_i is the ion cyclotron frequency. For a fixed shape of the sheared blob-filament like that illustrated in Fig. 6, we have noted above that shorting will occur when τ_eǁ < τ_κ. The entire blob-filament can be E × B drifting both in the radial and binormal directions. This will have no effect on the shorting mechanism since the electrons drift with the structure. However, if the blob-filament is subject to sheared flows, then the structure illustrated in Fig. 6 will distort with time. The relevant shearing velocity is v_s = R Δω_ζ0, where Δω_ζ0 is the change in toroidal angular velocity over a distance comparable to the parallel scale length of the filament.

Neglecting dynamical shearing for the moment, and assuming that the filament is already distorted as shown in Fig. 3, the relevant condition to fulfill for electron shorting of blob propagation in the collisionless kinetic limit is τ_eǁ < τ_κ or equivalently d_ǁc_sρ_s < δ_bRv_te. When T_i > T_e as it is for the blobs in this simulation, the estimate of the curvature velocity in Eq. (3) should be modified by the replacement c_sρ_s → v_tiρ_i. Thus, a more general condition for charge shorting in the collisionless regime takes the form

The preceding inequality is typically easily satisfied since c_s ≪ v_te and ρ_s < δ_b and in the presence of sheared structures d_ǁ < R. This condition, Eq. (4), and importantly the condition for this collisionless estimate to apply, will be compared with simulation parameters in Sec. V.

Since most of the previous simulation work on blobs has been done in the context of fluid models, in this section, a condition for electron shorting is developed in the collisional fluid limit. It will be sufficient to consider an isothermal reduced fluid model, as described, e.g., in Appendix A of Ref. 55.

In Bohm-normalized variables (i.e., time normalized to 1/Ω_i, space to ρ_s, and density to a reference density n₀), the electrostatic potential and plasma density obey the following coupled equations:

where the convective derivative is written as d/dt = ∂/∂t + v·∇ with v=e_z×∇Φ and A_dw is the adiabatic (drift wave) operator,

Here, for any quantity $Q$, ⟨$Q$⟩ ≡ $Q¯$ denotes the zonal or y average part, and in Eq. (7), {$Q$} ≡ $Q−Q¯$ denotes the fluctuating part. In the reduced model, (x, y) are the radial and binormal variables, and B = Be_z. Thus, A_dw enforces a Boltzmann response on fluctuations when the coefficient α_dw is large, in the spirit of the Wakatani–Hasegawa⁵⁶ adiabaticity parameter. Here, the collision frequency ν_ei0 is based on the reference density.

Other coefficients appearing in the model are as follows: the sheath conductivity parameter α_sh = ρ_s0/L_ǁ, where ρ_s0 is the reference sound gyro-radius and the curvature parameter β = 2ρ_s0/R, where R is the effective radius of curvature. In this reduced model, the parallel sheath current J_ǁsh is ultimately given in terms of n and Φ by closure relations in the sheath-limit and conduction-limited (collisional) regimes.

As is well known, the radial blob velocity v_b ∼ ∂Φ/∂y ∼ Φ/δ_b (in Bohm normalized variables) is obtained from the blob polarization potential Φ by balancing either (i) the sheath current α_sh J_ǁ in the sheath-connected limit or (ii) the ion polarization current (d/dt)∇·n∇Φ ∼ v·∇∇·n∇Φ ∼ nΦ²/δ_b⁴ in the inertial limit, with the curvature term β ∂n/∂y ∼ 2n/(Rδ_b).

For a Gaussian shaped blob profile in density, n(y) is even about the blob centroid and ∂n/∂y is odd. This gives rise to an odd or “dipole” potential profile Φ(y) ∼ ∂n/∂y, which, in turn, provides an even poloidal electric field ∂Φ/∂y to propel the radial E × B blob propagation. When α_dw is large enough to dominate the vorticity equation, it instead enforces Φ ∼ ln n, i.e., an even Φ(y) for even n(y) and no net radial propagation. In the language of turbulence theory, the cross-phase between n and Φ is π/2 for small α_dw (high collisionality), which promotes curvature driven radial transport, but the cross-phase tends to 0 for large α_dw (low collisionality, adiabatic electrons), which mitigates that transport channel.

A heuristic condition for adiabatic shorting of the curvature-driven charge polarization in fluid theory is obtained by estimating the potential response arising from Eq. (5),

where k_y and k_⊥ are typical inverse spatial scales (blob sizes) and γ is a characteristic inverse blob timescale. Here, the sheath term, which is not relevant to this discussion, has been dropped for simplicity. Comparing the curvature and electron parallel streaming terms in the numerator, we see that the phase between δΦ and δn transitions from hydrodynamic-like (δΦ ∝ ik_y δn) to adiabatic-like (δΦ ∝ δn) when α_dw > k_yβ or restoring dimensional variables when

At even larger α_dw ≫ γk_⊥² ≫ k_yβ, the transition to adiabatic electrons is complete. Equation (10) can be written in the form

Fluid theory, and the above condition, is valid for the electrons in the short mean free path limit k_ǁv_te < ν_ei. At the validity boundary of the fluid regime k_ǁv_te ∼ ν_ei and identifying k_ǁ ∼ 1/d_ǁ, Eq. (11) agrees, as it should, with the kinetic result given in Eq. (4). Thus, the charge-shorting mechanism should be present in a 3D fluid turbulence simulation, provided that the required conditions, including Eq. (11), are satisfied.

To decide whether the charge shorting mechanism is active in the XGC1 simulation, we first evaluate the collisionality ratio k_ǁv_te/ν_ei at the location ψ = 1.01 where the relative fluctuation levels maximize. Employing the parameters n_e = 5.6 × 10¹⁹ m⁻³, electron temperature T_e = 40 eV, ion temperature T_i = 155 eV, and B = 4.4 T and taking k_ǁ = 2π/L_ǁb where L_ǁb is the parallel dimension of the blob ∼0.9 m, we find k_ǁv_te/ν_ei ∼ 2.4, which puts the blob in the collisionless kinetic regime. It may also be reasonable to estimate k_ǁ ∼ 1/d_ǁ, where d_ǁ, defined in Fig. 6, is estimated next. This leads to the same conclusion: the blob is at least marginally in the kinetic regime.

A detailed analysis was carried out on the filament shown in Fig. 7. Comparing it with the filaments shown in Fig. 3 also illustrates both the similarity and variability of these structures. The filament shown in Fig. 7 has its centroid at the location of the dashed vertical line at ζ₀ = 1.487. From a cut along this line, the full width half maximum (FWHM) in the θ variable was measured and converted to the distance along the connecting field line. This results in the estimate d_ǁ = 18 cm. For the same blob-filament, the perpendicular size δ_b was measured by converting its half width half maximum (HWHM) in ζ₀ to an equivalent distance in the binormal direction, yielding δ_b = 0.18 cm.

Substituting the preceding estimates into Eq. (4) yields d_ǁρ_sc_s(1+T_i/T_e)/(Rδ_bv_te) = 0.002 ≪ 1, indicating that the condition for charge sorting is easily satisfied. However, in this simulation, the filament shearing is sufficiently pronounced that a simple static picture like that implied in Fig. 6 may not apply. For the parameters given in the first paragraph of this subsection, one finds τ_eǁ = 0.07 μs and τ_κ ∼ 36 μs. From Fig. 5, we estimate the change in toroidal angular velocity to be on the order of Δω_ζ0 ∼ 20 to 40 krad/s, resulting in Δv_ζ0 = RΔω_ζ0 ∼ 18 to 36 km/s and τ_s ∼ 0.10 to 0.05 μs. Thus, the simulation is in the regime τ_s ∼ τ_eǁ ≪ τ_κ where charge shorting should occur but may be modified or even partly mitigated by dynamical shearing.

These considerations suggest a nearly adiabatic electron density response to potential perturbations. In Fig. 8, the electron density and electrostatic potential are overlaid for a field line in the closed surface region. The electron response is indeed nearly adiabatic. The adiabaticity is even more pronounced in the SOL (not shown) where the curves essentially overlay. This nearly adiabatic response mitigates radial blob propagation. Of course, there may be additional reasons for a nearly adiabatic electron density, including the presence of drift-type turbulence with finite k_ǁ. In the nonlinear state, the turbulence and filament shearing are coupled.

To study the actual radial velocity of the structures, several diagnostic methods were employed. First is the time-lag correlation method, described in Sec. III B and in more detail in  Appendix A. As explained there, the time-lag correlation method does not provide a very sensitive diagnostic for the radial velocity in these simulations. The resulting radial blob velocities shown as black data points in Fig. 9 come with large standard deviations on the order of 0.5 km/s. The means are nearly zero but show a slight tendency to decrease moving from the closed surface region into the SOL.

A second velocity diagnostic may be referred to as the velocity proxy method. In this method, we do not attempt to diagnose the actual motion of structures from one time to the next but rather look at the average E × B velocity v_E at locations where δn > 0 (see  Appendix B for details). The results are shown in Fig. 9 as a solid black curve, which peaks at about 300 m/s in the closed surface region and drops to nearly zero in the SOL. The same technique may also be applied to only the most positive parts of the structure δn > δn_th for some positive threshold δn_th. As seen in Fig. 10, the results are not overly sensitive to δn_th but tend to rise somewhat in the edge region as the threshold is raised.

A third method is to measure the blob velocity by a tracking method described in Ref. 40. Briefly, in this method, a region of plasma in the (R, Z) plane with δn/$n¯$ exceeding a threshold (0.1 here) is first identified as a blob in a given frame, i.e., at a given time. A blob with similar features (size, spatial proximity, etc.) is, then, identified in a subsequent frame, and the velocity of the structure is determined directly from its spatial motion in the given time interval. This analysis applied to 128 frames yielded 629 trackable blobs in the range ψ = 1.00 to 1.01 and 798 trackable blobs in the range ψ = 1.01 to 1.02. Outside of these ψ intervals, there were too few trackable blobs for a statistically significant result. The results for the mean blob velocity from this analysis are shown as thick blue bars in Fig. 9. Similar to the blob proxy method in the SOL, the blob radial velocities from the tracking method are nearly zero.

The blob tracking code also yields other information about the blob structure. Of particular interest here is the size distribution, which indicates a mean δ_b of 0.23 cm, close to the 0.18 cm of the example shown in Fig. 7.

Of course, the different “blob velocity” diagnostics all measure subtly different things. They are not expected to agree in detail. In particular, the correlation method is the most sensitive to small structures, while the proxy and tracking methods explicitly invoke a threshold on amplitude.

Finally, Fig. 9 shows an analytical estimate of what the blob velocity would be for these parameters in a reduced fluid model where the shorting mechanism is not operative. As shown in  Appendix C, the appropriate estimate is from the inertial or resistive ballooning (RB) regime of the two-region model.²⁴ In this regime, the radial blob velocity v_b is obtained by balancing the ion polarization current (d/dt)∇·n∇Φ (i.e., inertia) with the curvature (charge polarization) term. The inertial regime blob velocity is

where f_b ≡ δn/$n¯$ is a proxy for the perturbation amplitude. For nominal parameters given in the preceding, this evaluates to v_b,in ∼ 1.4 to 3.2 km/s for f_b in the range of 0.1 to 0.5. A more accurate evaluation of Eq. (12) using local plasma parameters and the corresponding rms values of f_b is shown as the red curve in Fig. 9. In principle, for pressure driven perturbations, one should take f_b ≡ δp/$p¯$⁠, which could increase v_b,in even further; however, for the present purpose, the conservative estimate shown in Fig. 9 will suffice.

From these diagnostics and estimates, it can be concluded that the observed blob velocities are indeed much smaller than Eq. (12) would yield where an ideal fluid model is applicable. This is consistent with the predicted shorting mechanism.

It is interesting to note that experimentally observed blob velocities in Alcator C-Mod span a range comparable to the scale of Fig. 9. Using a similar correlation technique to that described here, Zweben⁵¹ found typical radial turbulence velocities for C-Mod in the range 0.1 to < 1 km/s. In different sets of discharges, Grulke⁵⁷ reported 0.5 to 2 km/s using correlation and tracking, with most structures having v_b ∼ 0.5 km/s; Terry⁵⁸ found typical radial turbulence velocities that are less than 1 km/s. These experimental results are roughly consistent with the XGC1 simulation reported here, although the suppression of blob velocities in the SOL is not as strong in the experiments as in the simulation. We can also speculate that the occasional larger velocities reported experimentally in C-Mod (∼1 km/s) may have been filaments that were unaffected by shearing, perhaps because of flow variations throughout the discharge, which is much longer than 95 μs that could be simulated here.

As previously noted, the mechanism for dipole charge shorting, which is the main point of Sec. IV, is different from the mitigation of dipole charges by “Boltzmann spinning.”⁵⁴ The dipole charge shorting condition, given in the collisionless electron limit by Eq. (4) and for collisional fluid electrons in Eq. (11), is applicable when the filament is not aligned along the magnetic field. In this case, the parallel motion of the electrons can respond to the dipole pressure gradient and potential perturbations that are responsible for radial blob convection and thereby slow the convection.

On the other hand, the Boltzmann spinning mechanism is most easily understood when the filament is perfectly aligned along the magnetic field but has a parallel density gradient. Assuming that the parallel electric field and parallel pressure gradient terms dominate Ohm's law, a Maxwell-Boltzmann monopole potential with a parallel variation ΔΦ₁ ∼ ln Δn results (in dimensionless variables assuming isothermal electrons and a parallel filament density variation Δn). It was shown in Ref. 54 using 3D fluid simulations and analytical analysis that Boltzmann spinning develops under these conditions and mitigates the dipole polarization charges when the monopole potential Φ₁ exceeds the dipole potential Φ₂ (within an order unity factor). Boltzmann spinning and polarization charge rotation were also seen in 2D simulations⁵⁵ using a model similar to Eqs. (5) and (6).

Other mechanisms can also give rise to spinning, such as a cross field temperature gradient in sheath-connected blobs,⁵⁹ resulting in the same effect (reduction of radial velocity and a transfer of radial velocity to the poloidal component). While this particular mechanism is not directly relevant here because the blobs are electrically disconnected from the sheaths, any cross field process that links the perpendicular electrostatic potential and temperature or density gradients would also have an analogous effect.

The simultaneous effects of dipole charge shorting for poloidally sheared filaments and spinning have not been studied, but it is clear that blob-filament spinning occurs whenever the electron response allows a monopole potential and that the mitigation condition Φ₁ > Φ₂ should be more easily satisfied if Φ₂ has been reduced by charge shorting. Thus, there may be synergistic effects. Indeed, in the simulation under study, any charges not mitigated by the shorting mechanism would be mitigated by Boltzmann spinning since it can be seen from Fig. 8 that the Boltzmann contribution dominates the total potential.

Ideally, the effects of charge shorting on blob motion would best be quantified by comparing simulation results with and without sheared flows satisfying Eq. (4) but with other conditions unchanged as much as possible. This would account for mechanisms in addition to charge shorting that could slow the blob motion. For example, one caveat of our analysis is that neutrals were present in the analyzed simulation; they could, in principle, affect the blob velocity through charge-exchange friction. The large computation expense of running the XGC1 code has so far made such direct comparative runs impractical. Shorter seeded-blob simulations are presently under investigation for this purpose.

Although the dynamics of filaments that are sheared by poloidally varying toroidal flows were studied in the present paper using results from a gyrokinetic simulation, the phenomena of poloidal shearing and charge shorting should, in principle, be accessible to fluid codes in the appropriate parameter regime. Specifically, Eq. (11) must be satisfied, which, in practice, is usually easy, even in the strongly collisional plasmas required for the validity of fluid models. However, for marginally collisional SOL plasmas, it is well known that fluid models may be subject to some inaccuracies in the electron parallel response, and the representation of charge shorting would be more qualitative. Also, some mechanisms for generating sheared flows such as direct particle loss may not be present.

In this paper, we have not addressed the causes for the observed poloidally varying toroidal flows. This is certainly an interesting question but one which is well beyond the scope of this paper. We speculate that X-point ion losses, neoclassical physics, turbulence physics, and the transition from closed to open field lines are involved. All these likely play a role in establishing the electrostatic potential in the edge and SOL regions. Some aspects of the interplay among these mechanisms have been discussed in other works, both for XGC1^14,40,52 and fluid turbulence codes.^47,60

The primary results of this paper are found in Eqs. (4) and (11), Figs. 3–6, and Fig. 9. Figure 3 visually illustrates the shearing of blob filaments within a flux surface in the gyrokinetic simulation. It is caused by poloidal (parallel) variation of toroidal flows. The filament shearing is quantified in Figs. 4 and 5 where it is shown that the shearing is the largest near X-points and in the scrape-off layer. Conversely, shearing is the smallest near the outboard midplane and in the closed surface region just inside the separatrix. The causes for these flow patterns are outside the scope of this paper. Rather, the focus of this paper is on the consequences of such flows, when they occur, on filament shearing and radial propagation. Our focus on the consequences of filament shearing on blob propagation is not meant to imply the absence of other, possibly synergistic or competing mechanisms relevant to blob propagation.

The sketch in Fig. 6 illustrates the fundamental mechanism of charge shorting by parallel electron currents when a filament is sheared by the parallel variation of toroidal flows. This effectively generates a k_||, which enhances the role of the electrons and enables the electron response to approach an adiabatic response. Analytical conditions for when the shorting mechanism is applicable are given for the collisionless kinetic electron limit in Eq. (4) and for collisional fluid theory in Eq. (11). Figure 9 compares the blob velocities diagnosed from the simulation with an analytical estimate of what those velocities would be from an ideal fluid model in the absence of the shorting mechanism. The shorting mechanism inhibits charge polarization of the blob-filament and therefore inhibits the mechanism by which blobs are propelled radially, slowing them down as observed in the simulation. The resulting blob velocities were also found to be broadly consistent with experimental observations.

While simple scaling arguments like those proposed in the original two-region model²⁴ and related models²⁵ give definite values for the blob velocity for given plasma parameters, experimentally and in codes, a statistical distribution of blob radial velocities is invariably observed.^5–7 One reason for this is that the velocity depends on blob size, δ_b, and the turbulence generates a statistical distribution of sizes. Given the present study, it is also of interest to speculate that another reason for the distribution of blob velocities is that the blob-filaments are subject to sheared flows, which may vary throughout a discharge in space and time, on longer time scales than could be simulated here. Blob-filaments may be influenced differently by the charge shorting mechanism depending on the background sheared flows they encounter.

It is possible that the shearing mechanism could be deliberately employed to control blob velocities by biasing. Biasing has been tried experimentally with some success,⁶¹ motivated by other theoretical considerations.⁶² Strongly sheared flows in the SOL may already be present due to self-biasing in experiments that apply RF waves in the ion-cyclotron range of frequencies (ICRF).^46,63,64

In closing, we emphasize that it will be important to continue studies of kinetic and fluid effects on blob-filament dynamics in order to better understand and predict the impact of these turbulent structures on the near-SOL width and on main chamber wall interactions in the far SOL.

This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Fusion Energy Sciences, under Award No. DE-AC02–09CH11466 (sub-contract No. SO15882-C) and U.S. DOE Grant No. DE-FC02–99ER54512, using Alcator C-Mod, a DOE Office of Science User Facility. We thank J. Hughes for providing the Alcator C-Mod geometry and profile data. Discussions with members of the HBPS SciDAC project (Scientific Discovery through Advanced Computing Initiative: High Fidelity Boundary Plasma Simulation) are gratefully acknowledged. We acknowledge the use of computing resources on Summit at the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory through the INCITE program. The digital data for this paper can be found in Ref. 65.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The correlation method provides a measure of the velocity of structures by searching for two space-time points with maximum correlation. The assigned velocity is, then, v = Δx/Δt, where Δx is the spatial separation of the points and Δt is the time lag. However, many detailed considerations enter the actual implementation of this rather simple idea.

The correlation function for the density perturbation δn is defined by

where the data record starts at t = t₁ and stops at t = t₂. First, we choose a time delay, which is a multiple of δt, the time resolution of the dataset, i.e., Δt = n δt, where n is an integer parameter. Then, for a given x, the maxima of C must be located within a given search box. The n = 1 search box is chosen to contain a typical blob as it moves with one δt step. The choice of initial box requires some judgment and experimentation; it is not unique. Typical values used include data up to ±3 spatial grid-points from the reference location. Having determined the initial (n = 1) box, for n larger than 1, the box is expanded by a factor of n, so that blobs do not escape the search box. Thus, the correlation velocities so determined depend on the numerical parameters n and the search box dimensions. Different choices offer competing advantages in terms of velocity scale and resolution.

The velocity resolution is limited by the spatial resolution of the simulation and the largest Δt that can be accommodated without using an unduly large search box. The smallest velocity that can be detected is

where Δx is the grid (or node) spacing. If Δt is made large to improve resolution, then either some fast velocity will be missed or the search box must be increased to the point where adequate spatial resolution and spatial range are lost, with the latter due to the search box beginning to overlap the boundary of the domain.

In the present application, x may be interpreted as an abstract vector in FLF coordinates, x = (ψ, ζ₀). Since the fluctuations are statistically independent of the toroidal coordinate ζ₀, velocities may be determined for a sequence of toroidal positions. In Figs. 4, 5, and 8, we show the mean and standard deviations over ζ₀ and n = 1, … n_max, where the standard deviations are plotted as error bars. Typically, we found good results for n_max in the range of 3 to 5.

The method yields good results for the toroidal velocity because the toroidal resolution in FLF coordinates is very good, and because of periodicity, the search boxes never encounter a domain boundary. On the other hand, the method does not work so well for the radial velocity. The first issue is that the radial velocity is much smaller than the toroidal velocity, making an optimal choice of Δt for a simultaneous measurement of both difficult. This was mitigated by transforming the data into a frame rotating toroidally with the mean toroidal velocity. The second issue is that the radial resolution of the data limits the radial velocity resolution for the reasons just stated to about 1 km/s, and this is only over the middle of the radial domain, where boundary issues do not interfere with the calculation.

In this method, the blob velocities are not determined by direct detection of their motion, but rather, it is assumed that the plasma moves radially across field lines at the local E × B velocity. Because the E × B velocity can have a spatial variation on the scale of the blob itself (see, e.g., Fig. 7), an appropriate average velocity must be defined. Here, we define the blob velocity proxy as

where H is the Heaviside step function, δn_th is a threshold fluctuation level, B_p is the poloidal magnetic field, and the radial E × B velocity in FLF coordinates is given by

An approximation has been made in the final form of Eq. (B2): the contribution of the parallel variation ∝ ∂Φ/∂θ in FLF coordinates is assumed to be negligible specifically,

[In principle, ∂Φ/∂θ|_FLF makes a contribution to v_Eψ because the poloidal electric field in geometric angle coordinates (ψ, θ, ζ) is expressed in terms of both ∂Φ/∂ζ₀ and ∂Φ/∂θ in FLF coordinates.]

The blob velocity proxy given by Eq. (B1) averages v_Eψ over that part of the blob for which the fluctuating density exceeds a threshold value, δn > δn_th. Thus, for δn_th = 0, the average is taken over all locations with positive fluctuations (“blobs”). As δn_th increases, the computed proxy velocity typically increases, as shown in Fig. 9. At large δn_th/δn_rms, the statistics become poor.

It should be kept in mind that the two-region model is rather rudimentary in its representation of X-point geometry, and while useful for understanding blob velocity scaling laws, it is not expected to give accurate order-unity coefficients in those scaling laws or in the regime boundaries. Nevertheless, it is useful to understand the model's implications.

The fundamental parameters of the two-region model are

and the X-point geometry factor ε_x. Here, L_ǁ is the midplane to the divertor connection length ∼8 m and R ∼ 0.89 m is the major radius, and in the X-point region, on the ψ = 1.01 flux surface, we estimate ε_x from the flux expansion ratio to be ε_x ∼ RB_p(X-pt)/RB_p(midplane) ∼ 0.1. Blob and plasma parameters in the SOL from the simulation give n_e ∼ 5.6 × 10¹⁹ m⁻³, T_e ∼ 40 eV, T_i ∼ 155 eV, B = 4.3 T, and δ_b ∼ 0.18 cm, resulting in Λ ∼ 0.4 and Θ ∼ 0.4. This places the blobs on the resistive ballooning - resistive X-point (RB-RX) boundary, Λ = Θ, where the blob velocity is estimated from the RB scaling as

Taking f_b ∼ 0.1 yields v_b ∼ 0.6 km/s similar to Fig. 9 in the SOL. Equation (12) adds the additional curvature drive weighting from finite T_i. Of course, on the RB-RX boundary, the RX scaling will give the same result as Eq. (C3).

The parallel structure of the mode is broadly consistent with these findings. The rms fluctuation level at ψ = 1.01 decays to half of its maximum value at θ = −1.20 and θ = +1.07 in FLF coordinates. At these locations, the factor ε_x defined in the two-region model is 0.88 and 0.95, respectively, hardly different from its midplane value of 1. Thus, the filaments do not experience significant X-point effects such as enhanced inertia and are disconnected from the divertor. This eliminates the relevance of the connected-interchange (CI) regime and the sheath-connected (CS) regime, as expected since the CI regime requires Λ < ε_xΘ and the CS regime requires ε_xΘ > 1. In the simulations, however, the parallel decay and disconnection of the filaments may be assisted by effects other than resistivity, e.g., shearing by the poloidal variation of the background flows.