E × B flow shear drive of the linear low-n modes of EHO in the QH-mode regime

A new model for the edge harmonic oscillations (EHOs) in the quiescent H-mode regime has been developed, which successfully reproduces the recent observations in the DIII-D tokamak. In particular, at high E × B flow shear only a few low- n kink modes remain unstable at the plasma edge, consistent with the EHO behavior, while at low E × B flow shear, the unstable mode spectrum is significantly broadened, consistent with the low- n broadband electromagnetic turbulence behavior. The model is based on a new mechanism for destabilizing low- n kink/peeling modes by the E × B flow shear, which underlies the EHOs, separately from the previously found Kelvin – Helmholtz drive. We find that the differential advection of mode vorticity by sheared E × B flows modifies the 2D pattern of mode electrostatic potential perpendicular to the magnetic field lines, which in turn causes a radial expansion of the mode structure, an increase of field line bending away from the mode rational surface, and a reduction of inertial stabilization. This enhances the kink drive as the parallel wavenumber increases significantly away from the rational surface at the plasma edge where the magnetic shear is also strong. This destabilization is also shown to be independent of the sign of the flow shear, as observed experimentally, and has not been taken into account in previous pedestal linear stability analyses. Verification of the veracity of this EHO mechanism will require analysis of the nonlinear evolution of low- n kink/peeling modes so destabilized in the linear regime.


Introduction
How E × B flow shear influences plasma instability is a long-standing issue in the field of plasma physics. The E × B flow shear was considered responsible for the suppression of plasma edge turbulence [1], leading to a spontaneous transition from a low confinement state (L-mode) to a high confinement state (H-mode) in tokamaks [2]. The instabilities behind the edge plasma turbulence are small in scale (highn, n is the toroidal mode number) [3]. For large-scale (low-n) magnetohydrodynamics (MHD) instabilities, such as kink/ peeling modes [4], whether and how the E × B flow shear can significantly modify these stabilities is still an open question.

E × B flow shear drive of the linear low-n modes of EHO in the QH-mode regime
The kink/peeling mode is believed to be the dominant instability triggering large edge-localized modes (ELMs) at low collisionality in the H-mode plasma edge, a region socalled 'pedestal' where steep density and temperature gradients develop [4]. Large ELMs are a serious concern for next-generation fusion devices because of the intolerable high transient heat loads that they can impose on the plasma-facing wall materials [5]. However, the particle transport driven by ELMs can help prevent impurity accumulation and sustain plasma performance. A robust control of ELMs while retaining good thermal confinement and sufficient edge particle transport for impurity exhaust is therefore essential for fusion energy production.
One potential solution for large ELMs is the quiescent H-mode (QH-mode) regime [6][7][8], in which ELMs are replaced by the edge harmonic oscillations (EHOs). EHOs play an important role in avoiding ELMs and sustaining stationary operations with good thermal confinement and sufficient particle transport for impurity removal. The thermal confinement improvement factor, H 98,y2 , associated with the QH-mode is typically comparable or even higher than that in the type-I ELMy H-mode, and the impurity confinement time in the QH-mode is comparable to that in the type-I ELMy H-mode, suggesting that EHOs may exhaust impurities as effectively as large ELMs. Recent experiments and theories suggested [4,9] that the EHOs might have been saturated lown kink/peeling modes partially destabilized by the E × B flow shear. The EHO saturation was thought due to relaxation of the rotation shear by a braking effect induced by EHO-driven momentum transport [4,10] rather than growing explosively like an ELM. Understanding the physics of EHOs is required in extending the QH-mode operation regime to future fusion reactors.
QH-mode is usually obtained with strong toroidal rotation shear at the plasma edge [6][7][8]. The EHO is thought to be partially driven by the E × B flow shear, as supported by a couple of observations. For example, QH-mode favors counter-I p plasma rotation [8], which reinforces the diamagnetic rotation in determining the E r from radial force balance. Indeed, QH-mode plasmas are found to have a significantly deeper negative E r well than in the ordinary ELMy H-mode plasmas [6][7][8]. A more careful analysis of the E r profile in QH-mode plasmas suggests that it is more likely the E r shear determines the QH-mode access, rather than the impurity carbon ion toroidal rotation [11]. This conclusion was further confirmed in observations where the QH-mode is sustained at nearly zero neutral beam injection (NBI) torque and toroidal rotation, and a strong edge E × B shear is sustained by the counter-I p torque injection from the applied non-axisymmetric non-resonant magnetic fields [11,12].
Very recently, a new QH-mode regime was observed in the DIII-D tokamak at low rotation [13], where a low-n broadband electromagnetic turbulence replaced the coherent EHOs in the pedestal region. Here, a reduced E × B flow shear in the edge pedestal was linked to this change of edge mode behaviors. Low rotation and rotation shear are anticipated in the pedestal region of future fusion plasmas, since the plasma moment of inertia is much larger than that in today's devices, and the applied torque is limited. The QH-mode access under this condition therefore becomes important.
The effects of E × B flow shear on the peeling/ballooning mode stability have been studied numerically. Early pedestal stability analysis with single-fluid ideal-MHD codes ELITE [4] and MINERVA [14] showed that the toroidal rotational shear is stabilizing for high-n modes and destabilizing for low-n modes. The stabilization of high-n modes by toroidal rotational shear has been used to explain ELM mitigation by rotation in experiments. For example, significant reduction of ELM size and increase of ELM frequency was observed in the JT-60U tokamak with increased toroidal rotation in the direction counter to the plasma current (counter-I p ) and more negative radial electric field, E r , at the plasma edge [15]. However, the destabilization of low-n modes was estimated to be relatively small (a few percent or less in growth rate) [4].
Destabilization of the low-n modes was seen in recent MHD simulation using M3D-C1 [16] and BOUT++ [17] codes. A linear eigenmode analysis of the EHO-like low-n modes using M3D-C1 indicated that these modes were destabilized by rotation and/or rotation shear, independently of the rotation direction [16]. This is consistent with the earlier observations of EHOs [8] and supports the view that the EHOs are saturated low-n kink/peeling modes destabilized by edge E × B flow shear [4,9]. However, the underlying mechanisms for destabilization and the associated mode growth rates were not identified.
Simulation study with a two-fluid reduced MHD code BOUT++ suggested that the E × B flow shear destabilizes peeling-ballooning modes mainly through the Kelvin-Helmholtz (K-H) drive [17]. Previously, the K-H instability was shown to be stabilized in the tokamak plasma edge by strong magnetic shear [18]. However, a most recent BOUT++ simulation [19] indicated that the mode remained unstable even after the K-H drive was turned off artificially. These developments therefore beg the question: how could the E × B flow shear destabilize the kink/peeling mode in the absence of the K-H drive?
To address this issue, the effects of E × B flow shear on the linear stability of the low-n kink/peeling mode are studied in detail via a linear eigenvalue analysis guided DIII-D experimental data. We discover that the E × B flow shear can destabilize the low-n kink/peeling mode through a new mechanism, separately from the K-H drive [17]. Our analysis indicates that the differential advection of mode vorticity by sheared E × B flows modifies the 2D pattern of mode electrostatic potential perpendicular to the magnetic field lines, which in turn causes a radial expansion of the mode structures, an increase of field line bending away from the mode rational surface, and a reduction of inertial stabilization. These in turn enhance the kink drive as the parallel wavenumber increases significantly away from the rational surface where the magnetic shear is also strong. This enhanced destabilization by E × B flow shear may well be the underlying mechanism for the generation of EHOs and/or low-n broadband electromagnetic turbulence in the QH-mode pedestal. Our results reproduce the observations that at high E × B flow shear only a few low-n modes remain unstable, consistent with the EHO behavior, while at low E × B flow shear the mode spectrum is significantly broadened, consistent with the low-n broadband electromagnetic turbulence behavior observed recently in the DIII-D tokamak. Furthermore, this destabilization is shown to be independent of the sign of the flow shear, consistent with the experimental observations of EHO in the QH-mode plasmas at both negative and positive rotations [8].
On the other hand, the kink/peeling mode is considered the dominant instability leading to large ELMs at low collisionality, a condition anticipated in the pedestal region of future fusion plasmas [4]. Moreover, low rotation and rotation shear are very likely an essential feature of future fusion plasmas. The E × B flow shear in the pedestal region could therefore be much smaller than that with high-torque injection in today's devices, as indicated in the recent DIII-D experiments [13]. Our results suggest that the kink/peeling instability is strongly dependent on the E × B flow shear at the plasma edge. ELM behaviors in future low-rotation fusion plasmas can therefore be very different from those in today's devices where strong torque injection and strong E × B flow shear prevail at the plasma edge. The required active control and mitigation of the ELMs could therefore be different.
Furthermore, our findings point to a possible mechanism of ELM control via the resonant magnetic perturbations (RMPs). Significant E r profile change at the plasma edge usually accompanies the application of RMPs [20][21][22][23][24][25][26]. Our results suggest that the E × B flow shear may influence ELMs through its influence on the kink/peeling instabilities. Whether the RMP prevents ELMs by inducing E r profile change is therefore of high interest.
In the next section, a schematic 2-dimentional (2D) fluid model is used to illustrate how E × B flow shear induces radial expansion of mode electrostatic potential structure. In the third section, the radial expansion effect is studied analytically. In the fourth section, we propose a new model of kink/ peeling mode in the helical coordinate system [3], derive the linear eigenvalue equation, and calculate dispersion relation using the energy principle. In the fifth section, we conduct a linear eigenvalue analysis of a typical QH-mode discharge measured in DIII-D [13]. Finally, discussions and a summary are presented in the last section.

Illustration of the new physics using a schematic 2D fluid model
In this section, a highly simplified 2D fluid model is used to illustrate the effect of E × B flow shear induced radial expansion of mode electrostatic potential structures. The model is written in the slab geometry with coordinates x y z , , ( ), satisfying the right hand relation, × = e e e x y z , where e is the unit vector, y is the poloidal coordinate, = − x r r s is the radial coordinate and r s is the radial location of a rational magnetic surface. There is a uniform magnetic field B along the z axis and a background equilibrium sheared E × B flow, Equations (1) and (2) are the linearized vorticity equation and ion-pressure equation, respectively. We have separated the fluctuation quantities and equilibrium quantities as usually done in fluid simulations [27], e.g. the generalized vorticity ϖ ϖ ϖ = + ∼ 0 , the ion pressure = + p p p i i 0 i and the electrostatic potential ϕ ϕ ϕ = + 0 . The generalized vorticity is also known as 'polarization charge' [28].
The fluctuating generalized vorticity is defined by [27], where n i0 is the ion density, which is assumed to be a constant in the simulation domain for simplicity, = = × − n n 1 10 m s i0 i0 18 3 , = B b B/ is the unit vector along the background equilibrium magnetic field, is the leading order perpendicular ion flow velocity, ϕ ϕ = +p Z en / is the effective electrostatic potential fluctuation and Z i is the ion charge state. Here, = Z 1 i for singly charged ions. The electrostatic potential fluctuation, ϕ, is solved numerically from the Poisson's equation, The equilibrium generalized vorticity is defined by / / / is the equilibrium perpendicular ion flow velocity, ϕ ϕ = + eff0 0 / p Z en i0 i i0 is the equilibrium effective electrostatic potential and µ = ∇ × j B 0 0 / is the parallel equilibrium current density. Here, positive v iy0 is in the electron diamagnetic direction. The field line curvature is zero in this model, The j v 0 i 0 term of ϖ 0 is usually very small in a tokamak. In this model, the equilibrium parameters are selected to ensure ϖ = ∂ = − m n B v 0 x y 0 i i0 1 i 0 in the simulation domain, so that the K-H drive term, ϖ ∂ v Ex x 0 [17], does not appear in the vorticity equation, equation (1). We avoid involving the K-H drive, since the purpose of this paper is not to study the effects of K-H drive.
The amplitude of generalized vorticity fluctuation, ϖ ∼ , grows exponentially in this model, driven by a linear term on the right-hand side of equation (1) with a constant growth rate γ = 1 MHz 0 . In the meantime, ϖ ∼ is advected by an equilibrium sheared E × B flow in the poloidal direction, which distorts the vorticity structures, as shown in the firstrow plots of figures 3 and 4. We use the velocity difference, Ey v 0 0 s , as the advection velocity, implying that the simulation box is moving with the mode structure and the mode on the rational surface is at rest in the reference frame of the simulation box.
The ion-pressure fluctuation, p i , is driven by the radial , of the background equilibrium ion-pressure gradient and damped by a diffusion term on the right-hand side of equation (2) with a ion thermal conductivity, χ i . In order to ensure the numerical stability, a hyper ion thermal conductivity of × − 1 10 m s 4 2 1 is used, which is much larger than the typical value, ~10 m s 2 / , in the tokamak edge plasma. However, this damping term does not influence the growth of ϖ ∼ , since the equation (1) has been designed to avoid coupling p i into equation (1), so that the time evolution of ϖ ∼ is independent of p i or ϕ. This is similar to the situation of the MHD modes, like the kink mode, where the instability is mainly driven by the equilibrium current-density gradient, rather than the equilibrium ion-pressure gradient directly.
In order to study the effect of E × B flow shear on the mode electrostatic potential structures, we change the E × B shearing rate, S v , and the ion-pressure gradient, ∂ p x i0 , to examine five cases with different equilibrium profiles, as shown in figure 1. Note that some of the panels in figure  . For Cases 1-4, the ion-pressure gradient is the same on the rational surface, the right end of the x domain, as shown in figure 1(d). The ion diamagnetic frequency on the rational surface is ω = ∂ = − * k p eBn 3.14 MHz y x i is i0 s 0s ( )/ . The equilibrium electrostatic-potential profile, as shown in figure 1 on the rational surface and ϕ = 0 0right on the right end, which determines the value of v Ey0s . The equilibrium ion-pressure profile, on the right end, as shown in figure 1(c), is designed to ensure the equilibrium generalized vorticity equal to zero in we set the ion-pressure gradient on the right end equal to zero, , which guarantees no negative ion pressure in the simulation domain. This determines the value of v y i 0 in Case 3. Since ) has been fixed, v y i 0 changes with v Ey0s case by case, as shown in figure 1(e).
The initial conditions of the simulation is is the radial amplitude profile, which is a Gaussian distribution centering at the rational surface with a radial decay length of λ = 1 cm. The initial peak amplitude is set at . Figure 2 shows the time evolution of the peak value of the normalized electrostatic potential fluctuation, ϕ e T e0s / , in 4 µs of linear growth. It grows exponentially with time. By fitting the curves, we obtain the effective linear growth rate, γ, of the electrostatic potential fluctuation. As indicated in figure 2, γ increases with increasing E × B flow shear and is independent of the sign of the shearing rate or the ion-pressure gradient. It exhibits the same growth rate even with flat ionpressure gradient, which clearly demonstrates that the effect is induced by E × B flow shear, not by the ion-pressure gradient. The model has been specifically designed to avoid the influence or disturbance from other factors such as the ion-pressure gradient, so that the unique effect of E × B flow shear on the mode electrostatic potential structures can be clearly isolated and identified. The effective linear growth rate, γ, increases by ~30% from zero shearing rate to a shearing rate of This increase is due to a radial expansion of the mode electrostatic potential structures away from the rational surface, which increases the radial scale length, λ, and thus reduces the effective perpendicular wavenumber, , of the structures, as illustrated by the simulation in figures 3 and 4. Since the linear growth rate of the generalized vorticity fluc- eff , is fixed at γ 0 by the model. The reduction of ⊥ k 2 leads to an increase in the amplitude of the mode electrostatic potential fluctuation, ϕ, and a reduction in the radial curvature of the mode electrostatic potential structures, ϕ ∂ ∂ − r r 1 r r ( ), thus reduces the inertial stabilization  associated with the mode and increases the effective linear growth rate of ϕ. Let us look into the simulation results in more detail. Figures 3 and 4 shows ϖ ∼ (first row), ϕ (second row) and p i (last row) in the five cases at 2 µs and 4 µs during mode linear growth, respectively. We have the following observations. With increasing shearing rate S v from 0 to γ − 0 , the mode electrostatic potential structures expand radially away from the rational surface. This radial expansion is independent of the selection of the diffusion coefficient χ i in the ion-pressure equation, equation (2). The radial expansion increases with time from 2 µs to 4 µs. There is no radial expansion when the flow shear is absent, for Case 1, = S 0 v , even the ionpressure structures are broader than the mode electrostatic potential structures. The radial expansion occurs even without the ion-pressure gradient, ∂ p x i0 , as indicated in Case 5, which indicates that the radial expansion does not need the ion-pressure gradient and is solely induced by the E × B flow shear. The generalized vorticity structures do not expand, since it is designedly driven by a constant linear growth rate, γ 0 , and independent of p i or ϕ, so that we can pin down the origin of the radial expansion.
This radial expansion effect induced by the E × B flow shear will have significant influence on the linear growth of kink mode. As coupling to the shear Alfvén wave associated with field line bending, the kink drive, [28], is proportional to the parallel magnetic vector potential fluctuation, , which is further proportional to the mode electrostatic potential fluctuation, ϕ, and the effective parallel wavenumber, k [28]. Here, ω is the mode frequency in the plasmamoving frame of reference. In ideal-MHD plasmas, the field lines are frozen in the fluid element, thus the radial expansion of the mode electrostatic potential structures will lead to field line bending radially, = ∂ B A y r , significantly away from the rational surface driven by radial E × B convective motion. In a background magnetic field with strong magnetic shear, the angle between the helical mode structures and the helical background field lines increase as the mode structures extend radially away from the rational surface, i.e. the effective parallel wavenumber, ≈ − k k x L y s / , increases with x. Here, ≡ L qR s s 0 / is the magnetic shear length and ≡ ∂ s r q q r ( / ) is the magnetic shear parameter, which are measures of the strength of the magnetic shear, q is the safety factor and R 0 is the plasma major radius on the magnetic axis. Moreover, because of the field line bending, the equilibrium currentdensity gradient has a projection along the magnetic field lines, ∂ B j B r r 0 ( / ), which gives rise to the kink drive at radial locations significantly away from the rational surface, thus enhances the kink drive. Therefore, a radial expansion of the mode electrostatic potential structures increases the effective linear growth rate of the kink mode.
The radial expansion appears to be radially asymmetric about the rational surface. It expands more inwards (towards negative x) with negative shear, < S 0 v , but more outwards (towards positive x) with positive shear, > S 0 v , as shown in figures 3 and 4. The ion-pressure structures also exhibit the same radial asymmetry. The radial asymmetry disappears with flat ion-pressure gradient, as indicated in Case 5, demonstrating that the radial asymmetry originates from the ionpressure gradient. With increasing diffusion coefficient χ i in the ion-pressure equation, equation (2), the radial asymmetry is reduced, because the fluctuation amplitude of p i is reduced and the spatial gradients in the structures are smoothed, which reduces the contribution of ∇ ⊥ p i in the fluctuating generalized vorticity, ϖ ∼ . Another very interesting observation is about the tilting direction of the structures. The generalized vorticity structures are tilted in the same direction as the flow shear. It reverses direction as the sign of the shearing rate, S v , changes. However, the mode electrostatic potential and ionpressure structures appear to be tilted in the opposite direction of the flow shear. It is because that the radial expansion of the mode electrostatic potential structures preferentially arises from the locations where the generalized vor- eff , is minimum, i.e. the radial curvature of effective electrostatic potential is minimum, as shown in figures 3 and 4. The ion-pressure structures are determined by the radial E × B convection, , therefore follow the tilting direction of the mode electrostatic potential structures. Although the initial tilting angle is in the opposite direction with respect to the flow shear. Comparing 2 µs with 4 µs, we find that the tilting angle changes with time and the change is still in the same direction as the flow shear.
Finally, we clarify the fundamental mechanism for the radial expansion of mode electrostatic potential structures. The central reason is that the quantity advected by the equilibrium sheared E × B flows is the polarization charge, ϖ ∼ , not the fluctuating electrostatic potential, ϕ. The sheared flows advect the polarization charge poloidally, introducing a timevarying radial phase shift, which in turn modifies the radial distributions of polarization charge and thereby the radial curvature profiles of the mode electrostatic potential structures, finally leading to the radial expansion.
The detailed physical picture is as follows. For a radially localized mode electrostatic potential structure, the radial curvature, ϕ ∂ x 2 eff , peaks near the rational surface and the low curvature regions are in between the positive and negative potential eddies when there is no flow shear, as shown in figures 3 and 4. The sheared flows poloidally advect the low curvature regions adjacent to the rational surface. When the low curvature regions move to the poloidal locations radially aligned with the high curvature peaks on the rational surface, the mode electrostatic potential will expand radially through those locations simply due to low curvature.
On a longer timescale during mode growth, ω > − t s 1 , the radial expansion would appear periodically with a frequency of ω s , when the mode structure shifts beyond a half-wavelength poloidally, where ω λ ≡ k S y v s is the shear frequency. However, this is prevented in our case, since the linear growth rates of MHD instabilities, such as the kink/peeling modes, are usually much larger than the shear frequency, γ ω s , so that the growth time is too short to allow the structures to shift a half-wavelength poloidally. Thus, the radial expansion continues through the period of mode growth, until nonlinear effects of the kink/peeling modes set in. This is a universal fundamental mechanism, which is applicable to any instabilities with short growth time relative to the shearing time, τ ω < − s 1 , including the kink/peeling modes.

Analytical analysis of the radial expansion effect induced by E × B flow shear
In this section, the radial expansion effect of the mode electrostatic potential structures induced by E × B flow shear is studied analytically. Since the fluctuating vorticity, i.e. the polarization charge is the quantity advected by the equilibrium sheared E × B flows. The fluctuating vorticity can be written in the form is the initial radial amplitude profile, ω ω ω ≡ + E is the mode frequency in the laboratory-rest frame of reference, ω ω γ ≡ +i r is the mode frequency in the plasma-moving frame of reference, ω r is the real frequency, γ is the linear growth rate, y Ey s 0 s is the E × B Doppler shift frequency on the rational surface. The initial mode electrostatic potential All other definitions are the same as those in section 2.
The time evolution of the mode electrostatic potential fluctuations can be solved analytically or numerically through radial integrations of the Poisson's equation  Figure 5 clearly demonstrates the radial expansion of mode electrostatic potential structures significantly away from the rational surface, the increase of radial expansion with time and the increase of mode electrostatic potential amplitude with time. Note that the linear growth rate, γ, has been set to zero in this analysis. In addition, the tilt of potential eddies appears to be in the opposite direction of the flow shear and the vorticity tilt, consistent with the results of the 2D fluid simulation in section 2. This analytical study further clarifies that the mechanism underlying the radial expansion of the mode electrostatic potential structures is the poloidal advection of the fluctuating vorticity, i.e. the polarization charge.

Linear eigenvalue analysis of kink/peeling modes with E × B flow shear
In this section, we construct a simplified model of kink/ peeling mode in a helical coordinate system [3], derive linear eigenvalue equation, and calculate dispersion relation based on the energy principle [28].

Eigenfunction in a helical coordinate system
To study the kink/peeling mode at the tokamak plasma edge, a helical coordinate system [3] is used with coordinates R 0 is the toroidal arc length, R 0 is the plasma major radius on the magnetic axis, ε = r R 1 0 / is the inverse aspect ratio, = q m n s / is the safety factor on the rational surface, q 1 s at the plasma edge, m and n is the poloidal and toroidal mode number, respectively. On the rational surface, the helical coordinate is aligned with the local magnetic field; while away from the rational surface, magnetic shear causes the magnetic field to tilt locally with respect to the helical coordinate. The coordinates are related to the poloidal and toroidal angles as follows, . The electron diamagnetic direction is the positive direction of y.
For simplicity, a low-β large aspect-ratio tokamak equilibrium is used. The coordinates of the flux surfaces are given by  [27], is the quantity advected by the equilibrium sheared E × B flows, not the fluctuating electrostatic potential, ϕ. Then, the eigenfunction of the mode in this coordinate system can be written in the form is the initial radial amplitude profile of ϖ ∼ , i.e. the initial radial distribution of polarization charge. It is only a function of the plasma minor radius. The radial amplitude profile of ϖ ∼ at time t is which appears to have been modified by the differential advection, . For kink modes, the unstable region is inside of the rational surface (x 0 ⩽ ) where q q s ⩽ [28]. The initial radial amplitude profile of ϕ is assumed to be a half-Gaussian distribution function centering at the rational surface with a radial decay length of λ (the radial mode width), with H x ( ) the step function. Its first and second radial derivatives can be calculated analytically by ϕ ϕ 2 . The radial decay length λ is of ~1 centimeter, which is consistent with the fluctuation measurements of electron density and temperature in the pedestal region of DIII-D, indicating that EHOs are localized in the pedestal steep-pressure-gradient region with a radial width of a few centimeters, as shown by figure 10 in [16]. Substituting θ l and φ l expressions into the mode eigenfunction, we have where is the wavenumber in the y direction, = θ k m r / is the poloidal wavenumber and = φ k n R 0 / is the toroidal wavenumber.
The definition of the mode frequency in the laboratory-rest frame, ω ω ω ≡ + E , is the same as that in section 3. In a plasma with equilibrium sheared E × B flow, ω is r-dependent, However, the mode frequency in the plasma-moving frame, ω, is invariant with the plasma minor radius for a given eigenmode.  [29]. The shear frequency is defined by ω λ ≡ k S y v s [30]. Taylor expanding v Ey0 around the rational surface, we have Under the drift approximation [28], the leading order convective time derivative is given by / the E × B advection velocity in the y direction. Here, positive v Ey0 is in the electron diamagnetic direction. Note that the advection by parallel flow, ∂ v 0 , has been neglected here. This time derivative is calculated in the reference frame moving with the fluid, thus the observed frequency is the mode frequency in the plasmamoving frame of reference. In contrast, the time derivative ∂ t is calculated in the laboratory-rest frame of reference, thus the frequency is modified by the E × B Doppler shift. Then, we have ϖ ωϖ ϖ ϖ ω ω ω ϖ Note that in this analysis we only consider the linear problems, thus all nonlinear terms have been neglected.
The flow in a tokamak plasma is usually subsonic, i.e.
r 0 is the equilibrium radial electric field with ϕ ψ 0 ( ) the stationary electrostatic potential, which is a flux function, φ v i 0 and θ v i 0 is the toroidal and poloidal ion rotation velocity, respectively. The ion pressure gradient, ∂ p r i0 , is usually determined by the ion thermal and particle radial transport. The radial electric field can be modified through changing the ion rotations, especially the toroidal rotation, since the toroidal rotation in a tokamak can be easily controlled by external torque sources such as NBIs.
In a tokamak, a net equilibrium flow can be written as is the ion toroidal rotation frequency, which is a flux function, is the toroidal unit vector.
Note that the incompressibility condition is satisfied, and the E × B flow is compressible. The parallel flow helps to keep the incompressibility. The total parallel flow is

( ) ( ) and the poloidal rotation velocity is
One can see that only the poloidal rotation has no coupling with the radial electric field in a helical magnetic field.
In the helical coordinate system [3], the gradient along the magnetic field line is given by is the perpendicular magnetic perturbation. Its radial and y-direction

Note that
we have the parallel wavenumber along the background equilibrium magnetic field In the helical coordinate system [3], the curvature drift operator becomes

Model of kink/peeling modes
The model of kink/peeling modes includes four coupled equations: Equation (8) is the quasi-neutrality condition. Equation (9) is the parallel generalized Ohm's law, where the electronion collisions, electron inertia and kinetic effects such as Landau damping and trapped particles have been neglected. Equations (10) and (11) is the simplified electron and ion pressure equation, respectively. After the well-known gyro-viscous cancellation [28], the divergence of the polarization current and gyro-viscosity current is given by which describes the inertia of electrostatic motion. Here, the last term, ϖ ∂ v Er r 0 , is the K-H drive term, responsible for driving the K-H instability [17]. The equilibrium generalized ( ) [17] has been defined in section 2. One can see that the K-H drive in a tokamak is mainly proportional to the radial curvature of the perpendicular flow. However, the K-H instability is usually thought stable in the tokamak plasma edge due to a strong stabilization effect by magnetic shear [18]. The K-H drive will be neglected in the following eigenvalue analysis, since the purpose of this paper is not to study the effects of K-H drive.
We express the fluctuating generalized vorticity, in the quantities normalized by the values on the rational surface, / is the ion density gradient length. The divergence of the diamagnetic drift gives the interchange drive, with the first term representing the kink drive and the second term describing the field line bending [28].
After linearization and normalization, the four-coupled equations become   (13) is the inertia term; the second term is the field line bending term due to coupling to the shear Alfvén wave; the third term is the kink drive term; the fourth term is the K-H drive term and the last term is the interchange drive term. For a kink mode, it experiences only the average curvature which is usually favorable in a tokamak, therefore the interchange drive term is a weak stabilization term with Substituting equation (15) into (14), we have Then, substituting equations (15)-(17) into equation (13), noting that ω is invariant with the plasma minor radius, we arrive at the eigenvalue equation  The diamagnetic stabilization effect [27] has been included, since we are using the fluctuating generalized vorticity. According to the energy principle [28], we multiply the coefficients by Φ r and calculate the radial integral, , over a radial range, r r r a b ⩽ ⩽ , covering the whole pedestal region, e.g. = r a 0.85 a and = r a b with a the plasma minor radius at the last closed flux surface. Then, we obtain a quadratic equation in ω Solving this equation, we finally obtain the dispersion rela- nary part is the required linear growth rate.

Numerical linear eigenvalue analysis using DIII-D experimental data
In this section, we conduct a linear eigenvalue analysis using DIII-D experimental data in a typical QH-mode discharge with toroidal rotation velocity ramping down to zero [13].

Descriptions of the experiment and the discharge
Very recently, a new kind of QH-mode regime with increased pedestal pressure and width at low rotation has been discovered serendipitously in the DIII-D tokamak [13]. The low rotation was achieved by slowly ramping down the external injected torque during the plasma discharge through balancing the co-I p and counter-I p oriented external torque injection from neutral beams. A low-n broadband incoherent electromagnetic turbulence, so-called 'broadband MHD' [13], appears in the pedestal region, replacing EHOs. The low rotation leads to a reduced E × B flow shear in the pedestal with an accompanying stabilization of the EHOs and an increase in height and width of the pedestal, which is within the limits allowed by peeling-ballooning mode theory. It could be a very important discovery, since the plasmas in future fusion devices are anticipated to operate at low rotation due to large plasma inertia and limited external torque injection. Stationary operation with improved pedestal conditions but without detrimental ELMs is highly desirable for fusion energy production. Figure 6 shows the time histories of the toroidal rotation velocity at the pedestal top and the toroidal-mode-number spectrogram in a typical discharge #163518 of this experiment. The details of this discharge have been described in [13]. The toroidal rotation velocity and radial electric field profiles were measured by a charge-exchange-recombination spectroscopy (CER) diagnostics on DIII-D [32]. The toroidalmode-number spectrogram is the crosspower between two poloidal magnetic field fluctuation signals measured by two toroidally distributed Mirnov magnetic probes on the outer midplane wall of the DIII-D vacuum vessel. The 30° difference in the toroidal angle allows the determination of the toroidal mode number n up to 5, which is indicated by the color plotted. The color scale on the right hand side of the plot shows the relation between color and mode number. The time point for this analysis is 2. As shown in figure 6, the toroidal rotation velocity at the pedestal top ramps down from ~100 km s − in the counter-I p direction to zero by reducing the counter-I p NBI power. The low-frequency coherent modes are observed in the toroidalmode-number spectrogram at high rotation are the EHOs dominated by modes with toroidal mode number n = 1-4. Integrating the Mirnov signal over time, we have the amplitude of the poloidal magnetic field fluctuation. The n = 2 mode has the highest fluctuation level in the integrated signal at 2.35 s, suggesting that n = 2 is the most unstable mode. The dominant toroidal mode number of EHOs is typically n ⩽ 3 [6][7][8]. In the recent linear eigenmode analysis with M3D-C1 code, n = 2 is the most unstable mode at the experimental rotation level [16], consistent with the dominant EHO mode observed in the experiment. The EHO frequencies decrease with time as the toroidal rotation velocity decreases, mainly due to Doppler shift. The EHOs disappear at ~2.9 s, when the toroidal rotation velocity is reduced to below 50 km s −1 , then replaced by a low-n broadband electromagnetic turbulence, as shown in the toroidal-mode-number spectrogram. The propagation direction of EHOs is in the counter-I p direction, following the plasma rotation with dominantly counter-I p NBI, whereas some low-frequency components of the broadband electromagnetic turbulence follow the plasma current direction instead, probably due to Doppler shift with reduced E × B flows at low NBI torque.

Profile data and equilibrium inputs for the analysis
The analysis focuses on the strong coherent EHOs at 2.35 s in discharge #163518. The radial kinetic profiles as shown in figures 7(e) and ( f ), including electron density n e0 , temperature T e0 and ion density n i0 , temperature T i0 , are obtained by averaging the data points over a 100 ms time window centered at 2.35 s using a Python based profile-fitting program [33]. The equilibria used for profile fitting are produced with the EFIT [34] code based on the magnetic measurements. At each measurement time point, the data points are mapped from physical to poloidal magnetic flux coordinate and then all individually mapped data within the overall averaging time window are collected and fitted to a modified tanh (mtanh) function [35] in the pedestal region and a spline model in the core region except the ion temperature T i0 profile which is fitted to a spline model in the whole region. The electron density n e0 and temper ature T e0 are obtained from Thomson scattering (TS) diagnostics. Impurity (carbon) density n c0 , ion temperature T i0 and radial electric field E r0 are obtained from CER diagnostics. The radial electric field profiles including E 0.5 r0 (red), E r0 (black) and E 2 r0 (blue) are shown in figure 7(c) and their corresponding shearing rates S 0.5 v (red), S v (black) and S 2 v (blue) are shown in figure 7(d). The ion density n i0 profile shown in figure 7(e) is calculated by = − n n n 6 i0 e0 c0 , assuming carbon is the only impurity. The electron (ion) pressure p e0 (p i0 ) shown in figure 7(g) is taken as the product of the electron (ion) density and temperature, respectively.
The experimentally measured total pressure profile is used as a constraint for a kinetic equilibrium reconstruction, which is the sum of the measured electron and ion pressures and the fast ion pressure calculated using the Monte Carlo NUBEAM [36,37] module in the ONETWO [38] transport code. The flux surface averaged toroidal current density in the core plasma is determined from motional Stark effect (MSE) [39] measurements, while in the pedestal region, it is derived from the sum of the parallel bootstrap current given by the Sauter expression [40,41], the neutral beam driven current computed using NUBEAM and the Ohmic current density. The Ohmic current density is determined from the neoclassical expression for the resistivity with a loop voltage, which is assumed spatially constant and adjusted to give a net plasma current matching the measured total plasma current. Iterations are taken in the reconstruction process wherein the pressure profile is remapped to poloidal flux coordinate and the pedestal current density is also recalculated.
The final safety-factor q profile and flux surface averaged current-density j 0 profile from the kinetic equilibrium are shown in figures 7(a) and (b), respectively. According to the q profile, the q = 7 rational surface is located in the steep current-density gradient region where the kink/peeling mode is unstable, marked with purple vertical lines in figure 7, and the current-density gradient at q = 6 rational surface is positive where the mode is stable. Figure 7(h) shows ω A (blue), ω * e (black), ω * i (red) and ω E (blackish green) profiles, where the toroidal mode number n is set to 1.

Linear eigenvalue analysis using the new model
A linear-stability eigenvalue analysis of the dominant kink/ peeling mode in the edge pedestal is performed for the discharge at t = 2.35 s using the new model developed in section 4. The new mechanism for destabilizing kink/peeling modes by E × B flow shear is the flow-shear-induced radial expansion of the mode structures radially away from the rational surface.
Firstly, we substitute the DIII-D data into the model and solve the eigenvalue equation, equation (20), numerically without radial expansion (no effect of sheared flows), i.e. t = 0 in equation (6). The initial radial amplitude profile of the mode electrostatic potential structures, , is a half-Gaussian distribution function centering at the rational surface with the radial decay length, λ, the only undetermined parameter. We search in the parameter space of λ to maximize the linear growth rate, γ, for each toroidal mode number, n. For a given n, the mode is unstable only when a positive linear growth rate can be found, otherwise the mode is linearly stable and the corresponding parameter λ cannot be determined. This approach is reasonable, since in reality the mode that finally appears in the experiments is the most unstable mode for each n.
The K-H drive term, ω ϖ Φ ∂ k en y ci r 0 i0s / , in expression (19) has been neglected, since it is not the focus of this paper. One can see that there is no S v dependence in expression (19) except through the K-H drive term. Therefore, in the model, without the effect of radial expansion, there is nowhere the sheared E × B flows can influence the solution of the eigenvalue equation. The sheared E × B flows just introduce a Doppler shift, which modifies the mode real frequency, ω ω ω = + E r r , in the laboratory-rest frame of reference.
The stability-analysis result without radial expansion is shown in figure 8 with the magenta curves. The linear growth rate, γ, in figure 8(a) indicates that the modes are unstable in the low to intermediate n range, n = 2-9, consistent with the previous stability-analysis results of kink/peeling mode with the ideal-MHD code ELITE [4,9,16]. The most unstable mode is n = 6, which is considerably higher than the toroidal mode number of the dominant mode observed in the experiment, i.e. n = 2. The peak linear growth rate is γ = 1.34 MHz 0 . The real frequency in the plasma-moving frame of reference, ω π = f 2 r / , in figure 8(b), shows a nearly linear dispersion relation in the ion diamagnetic direction (negative). The radial decay length, λ, of the mode electrostatic potential structures in figure 8(c) is in the range of 3-5 mm, peaking with λ = 5.3 mm at n = 3. The radial amplitude profile of the mode electrostatic potential structures for n = 3 mode is shown in figure 9 with the dark yellow curve, which is a half-Gaussian Nucl. Fusion 57 (2017) 086047 . The radial width of the mode is generally comparable with the pedestal width.
When the radial expansion is taken into account, the radial amplitude profile of the mode electrostatic potential structures, ϕ r ( ), evolves with time according to equation (6), and is no longer a simple half-Gaussian distribution function. It can be solved analytically or numerically through radial integrations of the Poisson's equation with the radial decay length, λ, fixed at the value at t = 0, as described in section 3 by expression (4). Note that here the integrations are conducted in a helical coordinate system, we then have The real and imaginary parts of ϕ r ( ) for n = 3 mode at 2 v (blue), (e) electron density n e0 (black) and ion density n i0 (red), ( f ) electron temperature T e0 (black) and ion temperature T i0 (red), (g) electron pressure p e0 (black) and ion pressure p i0 (red), (h) ω A (blue), ω * e (black), ω * i (red) and ω E (blackish green), where the toroidal mode number n has been set to 1. the shear Alfvén time calculated on the rational surface. Compared with the profile at t = 0, i.e. the half-Gaussian distribution function shown with the dark yellow curve, the real part of the profile at µ = t 0.2 s, shown with the blue curve, appears to be broadened and enhanced in amplitude with a long tail expanding radially significantly away from the rational surface. The imaginary part reflects a tilt of the mode electrostatic potential structures with a radial phase shift, ϕ ϕ r r arctan Im Re [ ( )/ ( )], radially increasing away from the rational surface, as shown in figure 5. This profile change is induced by E × B flow shear, as indicated in expression (21), and it will have significant influence on the linear growth rate of the kink/peeling mode, as shown in figure 8(a).
This theory successfully explains the observations that EHOs always appear with only a few coherent low-n harmonics and that the wavenumber spectrum is significantly broadened at low E × B flow shear. We scan the equilibrium radial electric field, E r0 , by multiplying a factor from half to two times of the measured value, which changes the E × B flow and flow shear simultaneously, as shown in figures 7(c) and (d), but with other parameters fixed. The kink/peeling modes appear to be destabilized at low n, but stabilized at high n and intermediate n by the E × B flow shear. With increasing sheared flows the toroidal mode number of the most unstable mode moves towards low n, i.e. = n 5 with E 0.5 r0 , = n 3 with E r0 and = n 2 with E 2 r0 . The unstable modes with E r0 are n = 2-4, which are generally consistent with the experimental observations in figure 6. The n = 1 mode is stable in this analysis, since it is already stable at t = 0 (no effect of sheared flows). Note that the n = 1 mode is stable only for this particular case and it is unnecessarily stable for other profiles. The peak linear growth rate decreases with increasing sheared flows, which may facilitate the saturation, since a saturation is usually achieved when the nonlinear damping rate balances the linear growth rate. A smaller linear growth rate is easier to be balanced. In addition, the real frequencies (negative) are slightly reduced with increasing sheared E × B flows, as shown in figure 8(b).
In summary, these stability-analysis results indicate that with strong sheared flows, only a few low-n modes remain unstable, and with reduced sheared flows, the toroidal mode number spectrum is broadened, consistent with the observations of coherent EHOs at high rotation and a transition to broadband electromagnetic turbulence at low rotation in the recent DIII-D experiments [13].
To further study the effect of E × B flow shear on the mode stability, the equilibrium radial electric field, E r0 , is shifted radially outwards from negative shear to positive shear on the rational surface step by step, as shown in figure 10, but with other parameters fixed. The results of linear stability analysis are shown in figure 11 with different colors corresponding to different shifts. When E r0 is shifted radially, the E × B  Radial amplitude profile of the mode electrostatic potential structures for n = 3 mode without radial expansion effect, i.e. at t = 0, no effect of sheared flows (dark yellow), which is a half-Gaussian distribution function, The real (blue) and imaginary (red) parts of the radial amplitude profile of the mode electrostatic potential structures for n = 3 mode with radial expansion effect, i.e. at µ τ = ≈ t 0.2 s 0.4 A . The radial electric field profile, E r0 , measured in the DIII-D experiment has been used in this analysis.
shearing rate, S v , on the rational surface is reduced, until the blue curve where the rational surface nearly passes through the flat bottom of the E r0 well. Indeed, the broadest toroidal mode number spectrum is obtained with the smallest S v , as shown with the blue curve in figure 11(a), and as expected, the toroidal mode number spectrum is gradually broadened as S v decreases. With positive shear, as shown with the magenta curve in figures 10 and 11, the toroidal mode number spectrum becomes narrow again, demonstrating that the effect is independent of the sign of the flow shear, consistent with the experimental observations of EHOs at both negative and positive rotations [8].

Summary of new findings
This paper is motivated by recent observations in the DIII-D tokamak that coherent EHOs in QH-mode plasmas at high rotation (and hence high E × B flow shear at the plasma edge) changes into a broadband electromagnetic turbulence in the pedestal region at low rotation. We analyze and identify a new destabilization mechanism for low-n kink/peeling modes in the presence of strong E × B flow shear, separate from the K-H drive [17]. Our analysis indicates that the differ ential advection of mode vorticity by sheared E × B flows modifies the radial distributions of the mode generalized vorticity, [27], the mode polarization charge. This in turn causes a radial expansion of the mode electrostatic potential structures, an increase of field line bending away from the mode rational surface, and a reduction of inertial stabilization. These in turn enhance the kink drive as the parallel wavenumber increases significantly  away from the rational surface where the magnetic shear is also strong. The pedestal fluctuation measurements of electron density and temperature in DIII-D [16] indicate that the EHOs exhibit rather broad radial structure, covering the whole pedestal region or even extending inward beyond the pedestal. The predicted radial expansion of the mode structure is consistent with this observation.
Field line bending can result from the kink drive, i.e. , which in turn is proportional to k , i.e.
ϕ ω = A k / , as shown in [28]. In ideal-MHD, the field lines are frozen in the fluid, the radial expansion of the mode electrostatic potential, ϕ, will lead to field line bending radially, = ∂ B A y r , driven by radial E × B convection away from the rational surface. In the presence of magnetic shear, the angle between the mode structure and the background field line increases as the mode structure extends radially away from the rational surface. That is, the effective parallel wavenumber, ≈ − k k x L y s / , increases with x. Because of the field line bending, the equilibrium current-density gradient has a projection along the magnetic field lines, ∂ B j B r r 0 ( / ). This increases the kink drive at radial locations significantly away from the rational surface. In addition, the radial expansion reduces the radial curvature of the electrostatic potential structures, ϕ ∂ ∂ − r r 1 r r ( ), reduces inertial stabilization, and thereby increases the effective linear growth rate of the kink mode.
The physics effect of the radial expansion is illustrated with 2D fluid model in section 2 and analyzed in section 3. The equilibrium E × B flow shear causes a radial expansion of the mode structure through the poloidal advection of the mode generalized vorticity, which is the polarization charge, ϖ ∼ . This in turn introduces an increasing radial phase shift, modifies the radial distributions of mode polarization charge, modifies the radial curvature profiles of the mode electrostatic potential, and leads to radial mode expansion, as illustrated in figures 3-5.
On a longer timescale during mode growth, ω > − t s 1 , the radial expansion would appear periodically with a frequency of ω λ = k S y v s , when the mode structure shifts beyond a halfwavelength poloidally. However, this is prevented in our case, since the linear growth rates are usually much larger than the shear frequency, γ ω s . Thus, the radial expansion continues through the period of mode growth, until nonlinear effects of the kink/peeling modes set in.
To identify the mechanism of kink/peeling mode destabilization and the resulting variations in mode behavior from different E × B flow shear, kink/peeling mode is analyzed in a helical coordinate system [3] in section 4. A linear eigenvalue analysis is carried out in section 5 for a DIII-D QH-mode discharge (#163518) [13]. The eigenvalue analysis has successfully reproduced the experimental observations that, at high E × B flow shear, only a few low-n modes remain unstable, while at low E × B flow shear, unstable modes in a broad wavenumber spectrum is obtained. These results are consistent with the observation in recent DIII-D experiments [13] of coherent EHOs at high rotation (and hence high E × B flow shear at the plasma edge) and a transition to broadband electro magnetic turbulence at low rotation. MHD simulations indicate that modes with a broad wavenumber spectrum in the linear simulations usually will evolve into broadband turbulence in the nonlinear simulations due to significantly enhanced modemode couplings compared to a narrow wavenumber spectrum [17,19,27]. Our results are further consistent with the kink/ peeling mode properties revealed in simulations using the ideal-MHD codes ELITE [4] and MINERVA [14], and more recently using MHD codes M3D-C1 [16] and BOUT++ [17,19]. Finally, as demonstrated in figure 11, this destabilization effect is shown to be independent of the sign of flow shear, consistent with the experimental observations of EHOs at both negative and positive rotations [8].
The improved understanding of the driving mechanisms of EHOs and broadband electromagnetic turbulence could facilitate the access to the QH-mode regime in future fusion plasmas to prevent the kink/peeling mode induced collapsesthe large ELMs, in subsequent nonlinear development.

Future work
Only linear stability has been analyzed in this paper. The radial mode structure, as shown in figure 9 with the blue curve, as well as the radial mode width in our linear analysis, are consistent with those of EHOs observed by fluctuation diagnostics in experiments, as shown by figure 10 in [16]. However, verification of the veracity of this EHO mechanism will require analysis of the nonlinear evolution of low-n kink/ peeling modes in the pedestal region. These will be addressed in future publications.
For example, the radial expansion of mode structures identified here may increase the radial extent of mode-driven transport as the mode advances into the nonlinear phase. This may provide an explanation for the observation that localized EHOs in the pedestal region may prevent impurity accumulation in the plasma core region by exhausting impurities in the pedestal region. Recently, a careful impurity-transport analysis on DIII-D [42] indicated that the impurity confinement time in QH-mode is comparable to that in ELMy H-mode, suggesting that EHOs may exhaust impurities as effectively as large ELMs.
Furthermore, E × B flow shear break the mode symmetry across the rational surface [43]. This effect has been studied extensively regarding small-scale turbulence [44][45][46][47]. It was suggested that the E × B flow shear could shift the peak amplitude of the mode eigenfunction, ( ) [ ( ) / ], radially away from the rational surface by a distance ξ in a plasma with significant equilibrium density and temperature gradients [43]. The radial distance of shift generally increases with the shear frequency, ξ λ ω γ λ γ ≈ = k S 1 y v s / / / . This can break k versus −k mode symmetry (i.e. a dynamical preference for one sign of k ) in the presence of E × B flow shear, thus introducing a net momentum source needed to cause the so-called 'intrinsic rotation' [43], independently of externally applied torque. Can this effect be applied to the large-scale MHD modes? Our results suggest that the kink/peeling mode could be further destabilized by the symmetry breaking. The strength of this mechanism relative to that of mode radial expansion will be addressed in future publications.