Transport barrier and spinning blob dynamics in the tokamak edge

In this work, we investigate the dynamics of plasma blobs in the edge of magnetic confinement devices using a full-f gyrokinetic particle-in-cell code with X-point geometry. In simulations, the evolution of a seeded blob is followed as it approaches a naturally-forming zonal shear layer near the separatrix, where the blob is stabilized by a large spin induced by the self-consistent adiabatic electron response, and blob bifurcation and trapping are observed during the cross-field propagation of blobs. A new theoretical explanation in both the zonal free and zonal shear layer is constructed, where the dominant E×B spin motion is included. A theoretical condition for a transport barrier induced by the interaction between spinning blobs and the zonal shear layer is obtained, and its scaling is verified with simulations. The new theoretical framework, especially the transport barrier, is applicable to explain and predict various experimental phenomena. In particular, the transport barrier condition calculated with experimental parameters demonstrates that the blob radial transport for H mode is smaller than L mode in experiments.

Coherent structures having poloidally and radially localized density structures and extended along the magnetic field, known as plasma blobs, often dominate the cross-field transport of particles and energy in the edge region between the confined plasma and material surface in magnetic confinement devices.The cross-field propagation of plasma blobs is of intrinsic interest for a better understanding of edge transport in future toroidal fusion devices and has drawn significant attention from the physics community [1][2][3][4][5][6][7][8][9][10].Furthermore, we will see that rapidly spinning blobs, considered here, have vortex interactions with a sheared flow layer that are of fundamental interest, reminiscent of similar phenomena in 2D neutral fluids.
In the conventional paradigm [1] of convective blob transport, the charge separation induced by opposite ion and electron vertical drift leads to a dipole potential and such a potential leads to an outward E × B drift.Although much progress has been made in understanding the nonlinear dynamics of these structures [2,3], realistic experimental conditions are complicated and further investigation is warranted.For example, the mushroom shape of blobs [11] due to charge separation has been demonstrated in a linear plasma device; however, blobs in the tokamak have a more stable nearly circular cross-section [12,13] shape.Blobs with negative crossfield velocity (radial velocity in a tokamak) have been observed [12,14].The blob emission in the SOL (scrapeoff layer) region in H mode is less than L mode, which is observed in different tokamaks, such as NSTX [15] and DIII-D [13].
In this work, we investigate the blob dynamics by carrying out simulations using the edge-physics-aimed 5D particle-in-cell (PIC) gyrokinetic code XGC [16,17].The blobs are seeded inside the separatrix using a realistic magnetic equilibrium from a C-Mod H-mode discharge [18].We consider the interaction of blobs born just inside the separatrix with a naturally forming shear layer near the separatrix.In the simulations, blobs with a large amplitude are stable with large spin and small cross-field motion.It is found that under certain conditions the spinning blob bifurcates, or splits, into two distinct parts at the midplane where the inner part is trapped in the core region, and the outer part moves outward radially into the SOL region.Hence, the E r shear layer near the separatrix acts as a transport barrier.A new theory that includes the dominant E × B spin motion in the parallel Ohm's law, vorticity and continuity equations is required.The potential coming from the adiabatic electron response giving rise to the spin is treated as zero-order, whereas the dipole potential is first-order.The blob can be stabilized by spinning and distorted by the large E r shear near the separatrix.The blob's radial motion responds to the zonal electrostatic potential well with the blob being charged positive, consistent with the simulation results.The new theory may explain the negative radial velocity of blobs seen in experiments.The blob bifurcation in the separatrix region is also observed in experiments, and the transport barrier condition obtained here and calculated with experimental parameters demonstrates that the blob radial transport for H mode is smaller than L mode in experiments.

Simulation results.
Here, the electrostatic version [17,19] of XGC is used for simulations, which employs drift kinetic electrons and gyrokinetic ions with a logical sheath boundary condition [20] at the divertor plate.For simplicity, the source term is turned off, and no neutral-plasma interaction is included.The geometry and magnetic equilibrium are based on a lower single null Alcator C-Mod EDA high (H) mode discharge [18] #1160930033 with a plasma current of I p = 1.4MA, a toroidal magnetic field on the axis of 5.8T, an input power of 5.4MW, and a record volume averaged pressure exceeding 0.2MPa.We simplify the problem in two ways.First, we set the density and temperature profiles to be uniform to avoid complications of gradient-driven instabilities super-imposed on the seeded blob physics.Second, in order to reduce the required computational resources, we reduce the magnetic field by a factor of 8, resulting in a magnetic field magnitude of 0.6T at the separatrix at the outer-most midplane.(Experimental conditions will be restored in our final scaling estimates.)The ion grad-B drift is directed downwards.The ions and electrons have a uniform temperature profile T e = T i = 60eV and the density profile is uniform n e = n i = 10 19 m −3 .Note, since the gradients of blob density and the profile satisfy ∇lnn b > ∇lnn 0 , the turbulence generated by profile gradients is suppressed.Thus simulations here are used only to investigate the blob physics.
First, we carry out a simulation without blobs as a controlled case.No turbulence is observed.A small amplitude zonal potential develops at the separatrix and the flux surface at the "virtual" upper X-point due to the charge absorption at the divertor.This controlled simulation verifies turbulence is suppressed with uniform density and temperature profiles.We now carry out a simulation where a blob is seeded in the closed-flux region just inside the separatrix at ψ/ψ x = 0.94, where ψ is the poloidal flux and ψ x is the poloidal flux at the separatrix.Our simulations are not intended to model the main pedestal where Reynolds stress and plasma profile gradients play an important role.The radial size of the seeded blobs is 3.71 times the ion gyro radius.The structure of the seeded blob is Gaussian in the radial coordinate and the other binormal coordinate and a cosine function along the field line.The amplitude of density perturbation is an input parameter.When the density perturbation is small, the blob has an obvious charge separation and a radial outward motion; when the density perturbation is large enough, as shown in Fig. 1, the blob structure is stable due to a fast spin induced by the adiabatic electron response.The blob moves radially outward in the closedflux region at a speed much smaller than the spin velocity.This behavior is expected from previous fluid studies in the SOL in simplified geometry [21,22].There is also a small downward drift until the blob's radial motion stops.Then the blob bifurcates into two: one blob stays inside the separatrix and continues to drift downward; the other smaller blob moves across the separatrix.Hence, the separatrix region acts as a transport barrier for the inner blob.
The blob radial velocity v R b can be measured by taking the time derivative of the center-of-mass [23] location at each time step.The radial component of the E × B drift averaged over the blob electron density, v R E , indicates that the radial motion is caused by the average E × B drift.As shown in Fig. 2, v R b and v R E are in reasonable agreement.The dipole potential structure in the vertical direction must still exist since it causes a pure radial E × B radial motion.However, it is not visible in Fig. 1 because the adiabatic potential dominates.In what follows, the dipole potential structure will be treated as a first-order expansion of the zero-order positive potential due to the adiabatic electron response.
Theoretical explanation of results.We start with the usual fluid vorticity equation and the parallel Ohm's law along with the fluid continuity equation, which is used to eliminate density and get J (φ).Note that the spin induced by E × B drift dominates the blob motion, thus the parallel Ohm's law Eq.( 2) and the fluid equations include the E × B and spin motion in leading order.
Here ω pi is the ion plasma frequency, Ω i is the ion gyrofrequency, v E = E × B/B 2 is the E × B drift, φ is the potential, J is the parallel current, κ = b • ∇b is the magnetic curvature, p is the pressure, ν is the collision frequency, n is the total density and ∇ = b 0 • ∇.
Since the adiabatic electron response dominates the zero-order potential and the dipole potential appears in the next order, we solve the vorticity equation for the symmetric zero order (φ 0 , m = 0, k = 0) and first order dipole (φ 1 , m = 1) component of the potential in the local cylindrical blob coordinates where φ ∝ e −imθ+ik z .Here, θ is the local azimuthal angle (θ = 0 in the e R direction at the outboard midplane) relative to the blob center and z indicates the coordinate along the field line.
In zero order, assuming adiabatic electrons, the righthand side of Eq. ( 2) and hence the parallel current vanishes.In the first order, the κ term in Eq. ( 1) acts as a source term for the m = 1 dipole potential.After some algebra that combines the m = 1 first-order vorticity, density, and parallel current equations, one obtains where L 1 is an inertial term given by Here Ω b = c rB dφ0 dr is the spin frequency, v te is the thermal velocity and α ad = Te en0 dn0 dφ0 indicates the adiabaticity.The α ad parameter is introduced for the theoretical description; the simulations employ the full drift-kinetic electron response.
Based on the above equations, we can solve for the propagation of spinning blobs driven by the curvature mechanism which is well known in Eqs. ( 1) and (2) to charge-polarize the blob.In the present formalism, when σ ,1 is real (for example in the collisional limit) Re[φ 1 exp(−imθ)] ∝ cosθ results in an electric field in the e Z direction (up-down in the RZ plane) and hence a radially outward E × B blob propagation velocity.Poloidal blob propagation results from Eqs. ( 3) and ( 4) when the left-hand side of Eq. ( 3) gives φ 1 an imaginary part.In addition to dipole charge polarization, the blob motion is also affected by the background E × B flow, curvature, and grad-B drifts.
Effect of the zonal field.To include the effect of the zonal field and the resulting sheared flows, φ contains zero order terms describing the background zonal flow φ 0f as well as the blob spin φ 0 .The calculation is performed in the frame of the flow at the blob center so that ∂/∂t can be neglected.The resulting equation for the first order response φ 1 +φ 1f from the curvature and flow terms respectively is then The first term on the right-hand side of Eq. ( 5) is the source of the charge separation process, the second term is the source of the zonal field effect.Here local Cartesian coordinates at the outboard midplane are employed where x is radial and y is binormal, v E0f = v E0f e y = (c/B)b × ∇φ 0f and a small flow ordering ∇v E0f ≪ Ω b has been assumed to neglect quadratic products of flow terms and to simplify the operators on the left-hand side.Near the blob center at x = 0, where δ b is the blob radius, assumed to be circular, and L f is the scale length of the flow shear.
The size estimate of the terms on the right-hand side of Eq. ( 5) expresses the sources in terms of the charge polarizing force per unit mass acting on the spinning blob.As is well known [2] the curvature and grad-B effects give rise to the effective gravitational acceleration g ∼ 2c 2 s /R.The factor δn/n takes into account the blob amplitude over the background density [3].
The flow term in square brackets has two parts.The second part can be expressed as eφ 0 ρ 2 s eE 0f /(T e δ 2 b m i ) ∼ (Ω b /Ω i )(eE 0f /m i ) which shows that the force on the blob is proportional to the blob vorticity ωb ≡ Ω b times the zonal electric field force eE 0f /m i , i.e., the blob reacts to the zonal field like a (positive) charge since charge is proportional to vorticity through the Poisson equation.The actual interaction is complicated by the fact that in the flow frame, v E0f and E 0f change signs across the blob center, so that the blob actually reacts to the velocity shear, as will be made explicit in Eq. ( 6).
The first term in square brackets can be rewritten as which shows the interaction of the blob electric field with the vorticity of the flow, ω = v E0f /L f .The extra factor of δ b /L f makes the two terms symmetric noting that one factor of δ b has been pulled outside the square brackets.
Both of these terms show explicitly the mutual interaction of two different sources of vorticity arising from the blob and the shear layer.Blob and shear layer interactions have been studied previously in slab geometry fluid simulations [21,24] and experiments [25,26]; however, the vortex analysis and interpretation turn out to be particularly transparent in the present large spin ordering.
The radial dependence of the zonal field, and the resulting shearing effect on blobs, compete with spin which can stabilize the blob and preserve its coherency.The stabilization criteria for spinning blobs from this competition is Ω b > ∇v E0f ∼ v E0f /L f .The shearing rate can also affect the blob shape when it is enough to distort the blobs as indicated in the simulations.For instance, at 37.36µs in simulations, Ω b ≈ 1.70 × 10 6 s −1 is less than ∇v E0f ≈ 1.95 × 10 6 s −1 , resulting in the blob bifurcation (Fig. 1).Thus, rapidly spinning blobs usually have a more stable shape.Blob stability is widely observed in experiments [12,13], although direct experimental detection of blob spin has so far not been possible.Blob transport barrier.The radial drift velocity induced by dipole response due to the effect of the curvature drift charge separation is always outward, while the radial drift velocity induced by the zonal field effect depends on E 0f and its shear, which can be negative or positive (Fig. 3).Thus, when these two mechanisms compete to determine the radial velocity, a transport barrier for the spinning blobs can be induced.By comparing the flow terms on the right-hand side of Eq. ( 5) to the curvature term, it is possible to derive the condition for a blob transport barrier.Recall that when the left-hand side is real (dissipative) the curvature term drives the blob radially outward, across the separatrix and into the SOL.The condition for the flow terms to dominate the curvature term is where the estimate v E0f = min(L f , δ b )v ′ E0f has been used to cancel one factor of L f or δ b in the final expression.See the comment following in Eq. ( 5).In deriving Eq. ( 6) the condition eφ 0 /T e ∼ δn/n for an adiabatic blob has been employed.As shown in Fig. 4, for the simulation reported here, the inequality in Eq. ( 6) is satisfied within factors of order unity, and the dynamics observed in the simulation agrees well with our theoretical expectation for net radial repulsion by the shear layer and shear-driven distortion and bifurcation of the blob.Furthermore, as discussed subsequently, the transport barrier condition is also relevant for experiments.

Comparison with experiments.
Comparison between theory, simulations, and experiments is a good tool to verify and validate the theory.We note the following qualitative and quantitative correspondences.(a) The new theoretical framework does not conflict with the dipole potential structure observed in an experimental cross-correction analysis [27], since the dipole potential is also included in the theory.(b) Blob bifurcation in the separatrix region has been clearly observed in C-Mod [28,29] and NSTX [30], which agrees with our simulations and theoretical expectation.(c) Negative radial motion of blobs in TCV [14] and NSTX [12] is potentially due to their interaction with shear layers.(d) The zonal field effect in the edge region is larger in H mode than in L mode.Thus blob emission into the SOL is expected and measured to be lower in H mode than in L mode, for example, in DIII-D [15] and NSTX [13].A 3kHz oscillation in the edge zonal field region of NSTX was found to be anti-correlated [31,32] with blob emission into the SOL.The above two experimental results are consistent with our theoretical prediction.(e) According to the experimental data [28,33] E0f ≈ 2 × 10 6 s −1 for H mode, 5 × 10 5 s −1 for L mode, Ω i ≈ 2.25 × 10 8 rad/s and R ≈ 0.89m in the separatrix region, C TB ≈ 1.09 for H mode and C TB ≈ 0.27 for L mode, which shows a good quantitative agreement.(f) Recently a detailed study of blob-vortex interactions with shear layers in L and H mode has been carried out for the RFX-Mod device.The results [26] are in qualitative agreement with the present theory.Thus, the present theoretical framework for blob structure and motion is consistent with the experimentally observed behavior of blobs.

Summary.
In this letter, we investigated blob dynamics by seeding a blob as an initial condition in the comprehensive full-f XGC particle code with gyrokinetic ions, drift-kinetic electrons and realistic Xpoint geometry in the electrostatic limit.We find that spin dominates the blob motion due to the adiabatic electron response.The dipole potential induced by the non-adiabatic part is smaller compared to the non-spinning case but is still responsible for radially outward motion.We take advantage of this ordering to develop a theory where the adiabatic potential is zeroth order and the dipole potential is first order.This theory provides a good explanation of the simulation results and it is able to explain several experimental observations.A transport barrier criteria is obtained from the current theory and the transport barrier condition the simulation and from experimental data is consistent with the behavior of blobs in the simulation and in experiments.The current theoretical framework can be extended and enriched to include kinetic and other effects, enabling more accurate descriptions of the experimental behavior of blobs in magnetically confined devices.
It is noted that this work does not aim to model the buildup of the whole pedestal.This work aims to isolate the spinning blobs and zonal field, and investigate the zonal field effect on the spinning blobs.Although the assumptions such as the blob seeding and flat profiles are used in our simulations, these gyrokinetic simulations with realistic geometry and equilibrium are first principles with the drift kinetic electrons, including both adiabatic and non-adiabatic electron responses.In the future, we will investigate the turbulent interaction of many blobs with realistic magnetic geometry and profiles in different parameter regimes.We will also consider different birth rates in L-mode and H-mode, and their impact on the SOL width, leading to a more complete theoretical transport model.

FIG. 1 .FIG. 2 .
FIG. 1.The blob structure of the electrostatic potential φ in the RZ plane at different times.The left line indicates the flux surface where the blob is initialized and the right line indicates the separatrix surface.For instance, at t = 28.02µs, the blob radial velocity v R b ≈ 76.16m/s, the blob spin velocity vs ≈ 11565.12m/s, the averaged blob perturbed density δn is about 0.51 times of background density, and the adiabaticity eφ/T /(δn/n) 1.

FIG. 3 .
FIG. 3. The time evolution of zonal potential (left) and trajectories of the center-of-mass blob location (right).In the left figure, the solid lines are the zonal potential at different time and the black dashed line indicates the separatrix.In the right figure, the solid line illustrates the time history of the center-of-mass location in the RZ plane; the colored dots indicate the center-of-mass locations at different times; the dashed line indicates the separatrix surface.Note that the zonal potential is due to particle loss in the divertor.