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

To study the kink/peeling mode at the tokamak plasma edge, a helical coordinate system [3] is used with coordinates $\left(x=r-{{r}_{\text{s}}},y={{l}_{\theta}}-\varepsilon q_{\text{s}}^{-1}{{l}_{\phi}},z={{l}_{\phi}}+\varepsilon q_{\text{s}}^{-1}{{l}_{\theta}}\right)$, where ${{r}_{s}}$ is the plasma minor radius of a rational surface, ${{l}_{\theta}}=r\left(\theta -{{\theta}_{0}}\right)$ is the poloidal arc length starting from the poloidal angle ${{\theta}_{0}}$, ${{l}_{\phi}}={{R}_{0}}\phi$ is the toroidal arc length, ${{R}_{0}}$ is the plasma major radius on the magnetic axis, $\varepsilon =r/{{R}_{0}}\ll 1$ is the inverse aspect ratio, ${{q}_{\text{s}}}=m/n$ is the safety factor on the rational surface, ${{q}_{\text{s}}}\gg 1$ 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, $y=r\left(\theta -{{\theta}_{0}}-q_{\text{s}}^{-1}\phi \right)$, $z={{R}_{0}}\phi +\varepsilon q_{\text{s}}^{-1}r\left(\theta -{{\theta}_{0}}\right)$, ${{l}_{\theta}}\,=\,\,\left(y\,+\,\varepsilon q_{\text{s}}^{-1}z\right)/\left(1\,+\,{{\varepsilon}^{2}}q_{\text{s}}^{-2}\right)\,\approx \,y\,+\,\varepsilon q_{\text{s}}^{-1}z$, ${{l}_{\phi}}=\left(z-\varepsilon q_{\text{s}}^{-1}y\right)/$$\left(1+{{\varepsilon}^{2}}q_{\text{s}}^{-2}\right)\approx z-\varepsilon q_{\text{s}}^{-1}y$. 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 $R={{R}_{0}}+r\cos \theta$ and $Z=r\sin \theta$. The toroidal magnetic field is ${{B}_{\phi}}={{B}_{\phi 0}}{{R}_{0}}/R$ and the poloidal magnetic field is ${{B}_{\theta}}={{B}_{\theta 0}}(1+\varepsilon \Lambda \cos \theta )={{R}^{-1}}{{\partial}_{\text{r}}}\psi$, where $\Lambda ={{\beta}_{\text{p}}}+{{l}_{\text{i}}}/2-1$, ${{\beta}_{\text{p}}}$ is the poloidal beta, ${{l}_{\text{i}}}$ is the internal inductance and $\psi$ is the poloidal magnetic flux. $B={{\left(B_{\phi}^{2}+B_{\theta}^{2}\right)}^{1/2}}$ is the background equilibrium magnetic field with ${{B}_{0}}={{\left(B_{\phi 0}^{2}+B_{\theta 0}^{2}\right)}^{1/2}}$. The safety factor $q(r)=\varepsilon {{B}_{\phi 0}}/{{B}_{\theta 0}}$ is a function of the plasma minor radius.

For the kink mode, the poloidal variation of mode amplitude is unimportant, thus neglected here. The fluctuating generalized vorticity, $\tilde{{\varpi}}\,\equiv {{m}_{\text{i}}}{{e}^{-1}}{{B}^{-2}}{{\nabla}_{\bot}}\centerdot$$\left(e{{n}_{\text{i}0}}{{\nabla}_{\bot}}\tilde{{\varphi}}\,+{{\nabla}_{\bot}}{{\tilde{{p}}\,}_{\text{i}}}\right)={{m}_{\text{i}}}{{n}_{\text{i0}}}{{B}^{-2}}\nabla _{\bot}^{2}{{\tilde{{\varphi}}\,}_{\text{eff}}}$ [27], is the quantity advected by the equilibrium sheared E  ×  B flows, not the fluctuating electrostatic potential, $\tilde{{\varphi}}\,$. Then, the eigenfunction of the mode in this coordinate system can be written in the form

$$\begin{eqnarray}&&\tilde{{\varpi}}\,(t,r,\theta,\phi )=\hat{\varpi}(r)\exp \left[-\text{i}\hat{\omega}t+\text{i}m\left(\theta -{{\theta}_{0}}\right)-\text{i}n\phi \right]\end{eqnarray} \tag{ 5 }$$

where $\hat{\varpi}(r)\equiv {{m}_{\text{i}}}{{e}^{-1}}{{B}^{-2}}{{\nabla}_{\bot}}\centerdot \left(e{{n}_{\text{i0}}}{{\nabla}_{\bot}}{{\hat{\varphi}}_{t=0}}+{{\nabla}_{\bot}}{{\hat{p}}_{\text{i},t=0}}\right)$ is the initial radial amplitude profile of $\tilde{{\varpi}}\,$, i.e. the initial radial distribution of polarization charge. It is only a function of the plasma minor radius. The radial amplitude profile of $\tilde{{\varpi}}\,$ at time t is

$$\begin{eqnarray}&&\hat{\varpi}(r)\exp \left[-\text{i}{{k}_{y}}\left({{v}_{Ey0}}-{{v}_{Ey0\text{s}}}\right)t\right]={{m}_{\text{i}}}{{e}^{-1}}{{B}^{-2}}{{\nabla}_{\bot}}\centerdot \left(e{{n}_{\text{i0}}}{{\nabla}_{\bot}}\hat{\varphi}+{{\nabla}_{\bot}}{{\hat{p}}_{\text{i}}}\right)\end{eqnarray} \tag{ 6 }$$

which appears to have been modified by the differential advection, $\left({{v}_{Ey0}}-{{v}_{Ey0\text{s}}}\right){{\partial}_{y}}$, induced by the sheared E  ×  B flows. Note that ${{k}_{y}}\left({{v}_{Ey0}}-{{v}_{Ey0\text{s}}}\right)$ is r-dependent. Here, $\nabla _{\bot}^{2}\hat{\varphi}=\left[{{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\right)-k_{y}^{2}\right]\hat{\varphi}$.

For kink modes, the unstable region is inside of the rational surface ($x\leqslant 0$) where $q\leqslant {{q}_{\text{s}}}$ [28]. The initial radial amplitude profile of $\tilde{{\varphi}}\,$ is assumed to be a half-Gaussian distribution function centering at the rational surface with a radial decay length of $\lambda$ (the radial mode width), ${{\hat{\varphi}}_{t=0}}(x)={{\hat{\varphi}}_{A}}H(-x)\exp \left(-{{x}^{2}}/2{{\lambda}^{2}}\right)$ with $H(x)$ the step function. Its first and second radial derivatives can be calculated analytically by ${{\partial}_{\text{r}}}{{\hat{\varphi}}_{t=0}}=-{{\hat{\varphi}}_{A}}H(-x)x{{\lambda}^{-2}}\exp \left(-{{x}^{2}}/2{{\lambda}^{2}}\right)$ and ${{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}{{\hat{\varphi}}_{t=0}}\right)=-{{\hat{\varphi}}_{A}}H(-x){{\lambda}^{-2}}\left(1-{{x}^{2}}/{{\lambda}^{2}}+x/r\right)\exp \left(-{{x}^{2}}/2{{\lambda}^{2}}\right)$. The radial decay length $\lambda$ 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}_{\theta}}$ and ${{l}_{\phi}}$ expressions into the mode eigenfunction, we have

$$\begin{eqnarray}&&\tilde{{\varpi}}\,(t,r,y)=\hat{\varpi}(r)\exp \left(-\text{i}\hat{\omega}t+\text{i}{{k}_{y}}y\right)\end{eqnarray} \tag{ 7 }$$

where ${{k}_{y}}={{k}_{\theta}}+{{k}_{\phi}}\varepsilon q_{\text{s}}^{-1}\approx {{k}_{\theta}}$ is the wavenumber in the y direction, ${{k}_{\theta}}=m/r$ is the poloidal wavenumber and ${{k}_{\phi}}=n/{{R}_{0}}$ is the toroidal wavenumber.

The definition of the mode frequency in the laboratory-rest frame, $\hat{\omega}\equiv \omega +{{\omega}_{E}}$, is the same as that in section 3. In a plasma with equilibrium sheared E  ×  B flow, $\hat{\omega}$ is r-dependent, since ${{\omega}_{E}}={{\omega}_{E\text{s}}}+{{k}_{y}}\left({{v}_{Ey0}}-{{v}_{Ey0\text{s}}}\right)$ is r-dependent. However, the mode frequency in the plasma-moving frame, $\omega$, is invariant with the plasma minor radius for a given eigenmode. The E  ×  B shearing rate is given by ${{S}_{v}}\equiv r{{\partial}_{r}}\left({{v}_{Ey0}}/r\right)$ in the cylindrical geometry and ${{S}_{v}}\equiv (r/q){{\partial}_{r}}\left(q{{v}_{Ey0}}/r\right)={{\partial}_{r}}{{v}_{Ey0}}+\left(\hat{s}-1\right){{v}_{Ey0}}/r$ in the toroidal geometry [29]. The shear frequency is defined by ${{\omega}_{\text{s}}}\equiv \lambda {{k}_{y}}{{S}_{v}}$ [30]. Taylor expanding ${{v}_{Ey0}}$ around the rational surface, we have ${{v}_{Ey0}}={{v}_{Ey0\text{s}}}+{{S}_{v}}x$. The radial gradient of $\hat{\omega}$ is ${{\partial}_{\text{r}}}\hat{\omega}={{\partial}_{\text{r}}}{{\omega}_{E}}\approx {{k}_{y}}{{S}_{v}}$.

Under the drift approximation [28], the leading order convective time derivative is given by ${{\hat{d}}_{\text{t}}}\equiv {{d}_{\text{t}}}+{{\tilde{{v}}\,}_{E\text{r}}}{{\partial}_{\text{r}}}$ with ${{d}_{\text{t}}}\equiv {{\partial}_{\text{t}}}+{{v}_{Ey0}}{{\partial}_{y}}$, ${{\tilde{{v}}\,}_{E\text{r}}}=-{{B}^{-1}}{{\partial}_{y}}\tilde{{\varphi}}\,$ the radial E  ×  B convection velocity and ${{v}_{Ey0}}=-{{E}_{\text{r0}}}/B={{B}^{-1}}{{\partial}_{\text{r}}}{{\varphi}_{0}}$ 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}}{{\partial}_{||}}$, 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 plasma-moving frame of reference. In contrast, the time derivative ${{\partial}_{\text{t}}}$ is calculated in the laboratory-rest frame of reference, thus the frequency is modified by the E  ×  B Doppler shift. Then, we have ${{d}_{\text{t}}}\tilde{{\varpi}}\,=-\text{i}\hat{\omega}\tilde{{\varpi}}\,+\text{i}{{k}_{y}}{{v}_{Ey0}}\tilde{{\varpi}}\,=-\text{i}\tilde{{\varpi}}\,\left(\hat{\omega}-{{\omega}_{E}}\right)=-\text{i}\omega \tilde{{\varpi}}\,$ and ${{\partial}_{\text{t}}}\tilde{{\varpi}}\,=-\text{i}\hat{\omega}\tilde{{\varpi}}\,$. 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. $v_{0}^{2}/v_{\text{Ti}}^{2}\ll 1$ with ${{v}_{\text{Ti}}}\equiv \sqrt{2{{T}_{\text{i0}}}/{{m}_{\text{i}}}}$ the ion thermal velocity, thus the ion inertia term, ${{m}_{\text{i}}}{{n}_{\text{i0}}}{{\mathbf{v}}_{0}}\centerdot \nabla {{\mathbf{v}}_{0}}$, is small in the equilibrium equation. Then, we have the radial force balance equation, ${{\partial}_{\text{r}}}{{p}_{\text{i0}}}/{{Z}_{\text{i}}}e{{n}_{\text{i0}}}-{{E}_{\text{r0}}}=B{{v}_{\text{i}y0}}={{B}_{\theta}}{{v}_{\text{i}\phi 0}}-{{B}_{\phi}}{{v}_{\text{i}\theta 0}}$ [1], where ${{E}_{\text{r0}}}=-{{\partial}_{\text{r}}}{{\varphi}_{0}}$ is the equilibrium radial electric field with ${{\varphi}_{0}}(\psi )$ the stationary electrostatic potential, which is a flux function, ${{v}_{\text{i}\phi 0}}$ and ${{v}_{\text{i}\theta 0}}$ is the toroidal and poloidal ion rotation velocity, respectively. The ion pressure gradient, ${{\partial}_{\text{r}}}{{p}_{\text{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 ${{\mathbf{v}}_{0}}=K(\psi )\mathbf{B}+{{R}^{2}} \Omega (\psi )\nabla \phi$ [31], where $\Omega (\psi )\,\equiv \,\left(\partial {{\varphi}_{0}}+\partial {{p}_{\text{i0}}}/{{Z}_{\text{i}}}e{{n}_{\text{i0}}}\right)/\partial \psi =\,\left({{\partial}_{\text{r}}}{{\varphi}_{0}}+{{\partial}_{\text{r}}}{{p}_{\text{i0}}}/{{Z}_{\text{i}}}e{{n}_{\text{i0}}}\right)/R{{B}_{\theta}}\,=$$B{{v}_{\text{i}y0}}/R{{B}_{\theta}}$ is the ion toroidal rotation frequency, which is a flux function, ${{\mathbf{e}}_{\phi}}=R\nabla \phi$ is the toroidal unit vector. Note that the incompressibility condition is satisfied, $\nabla \centerdot {{\mathbf{v}}_{0}}=K(\psi )\nabla \centerdot \mathbf{B}+\nabla \centerdot \left[{{R}^{2}} \Omega (\psi )\nabla \phi \right]=0$ and the E  ×  B flow is compressible. The parallel flow helps to keep the incompressibility. The total parallel flow is ${{v}_{\text{i}||0}}=K(\psi )B+ \Omega (\psi )R{{B}_{\phi}}/B$ and the perpendicular flow is ${{v}_{\text{i}y0}}= \Omega (\psi )R{{B}_{\theta}}/B$. The total toroidal rotation velocity is ${{v}_{\text{i}\phi 0}}={{B}_{\phi}}K(\psi )+R \Omega (\psi )$ and the poloidal rotation velocity is ${{v}_{\text{i}\theta 0}}={{B}_{\theta}}K(\psi )$. 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 ${{\nabla}_{||}}=B{{\partial}_{||}}\left({{B}^{-1}}\right)$ with ${{\partial}_{||}}=\,\mathbf{b}\,\,\centerdot \nabla +{{B}^{-1}}{{\tilde{{\mathbf{B}}}\,}_{\bot}}\centerdot \nabla \,\approx \,{{\partial}_{z}}\,+\,\varepsilon \left({{q}^{-1}}-q_{\text{s}}^{-1}\right){{\partial}_{y}}+{{\tilde{{B}}\,}_{\text{r}}}{{B}^{-1}}{{\partial}_{\text{r}}}+{{\tilde{{B}}\,}_{y}}{{B}^{-1}}{{\partial}_{y}},$ where ${{\tilde{{\mathbf{B}}}\,}_{\bot}}=\nabla \times {{\tilde{{\mathbf{A}}}\,}_{||}}=-\mathbf{b}\times \nabla {{\tilde{{A}}\,}_{||}}$ is the perpendicular magnetic perturbation. Its radial and y-direction component is ${{\tilde{{B}}\,}_{\text{r}}}={{\partial}_{y}}{{\tilde{{A}}\,}_{||}}=\text{i}{{k}_{y}}{{\tilde{{A}}\,}_{||}}$ and ${{\tilde{{B}}\,}_{y}}=-{{\partial}_{\text{r}}}{{\tilde{{A}}\,}_{||}}$, respectively. With the Ampère's law, we have the parallel current perturbation, ${{\mu}_{0}}{{\tilde{{\mathbf{j}}}\,}_{||}}=\nabla \times {{\tilde{{\mathbf{B}}}\,}_{\bot}}=-\mathbf{b}\nabla _{\bot}^{2}{{\tilde{{A}}\,}_{||}}$. Note that ${{k}_{z}}={{k}_{\theta}}\varepsilon q_{\text{s}}^{-1}-{{k}_{\phi}}=0$, we have the parallel wavenumber along the background equilibrium magnetic field ${{k}_{||}}=\mathbf{b}\centerdot \mathbf{k}={{k}_{\theta}}{{B}_{\theta}}/B-{{k}_{\phi}}{{B}_{\phi}}/B=\varepsilon \left({{q}^{-1}}-q_{\text{s}}^{-1}\right){{k}_{y}}\approx {{k}_{\text{c}}}(m-nq)$, where ${{k}_{\text{c}}}=2\pi /{{L}_{\text{c}}}=1/q{{R}_{0}}$ and ${{L}_{\text{c}}}=2\pi q{{R}_{0}}$ is the connection length. Taylor expanding $q(r)$ around the rational surface, we obtain ${{k}_{||}}\approx -{{k}_{y}}x/{{L}_{\text{s}}}$, its radial gradient ${{\partial}_{\text{r}}}{{k}_{||}}=-{{k}_{y}}/{{L}_{\text{s}}}$ and radial curvature ${{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}{{k}_{||}}\right)=-{{k}_{y}}{{\partial}_{\text{r}}}L_{\text{s}}^{-1}$.

In the helical coordinate system [3], the curvature drift operator becomes ${{B}^{-1}}\mathbf{b}\times \boldsymbol{\kappa }\centerdot \nabla \approx -{{\left({{B}_{0}}{{R}_{0}}\right)}^{-1}}\left({{C}_{\text{n}}}{{\partial}_{y}}+{{C}_{\text{g}}}{{\partial}_{\text{r}}}\right)$ with ${{C}_{n}}=\cos \left(y/r+z/{{R}_{0}}{{q}_{\text{s}}}\right)-\varepsilon \approx \cos \left(\theta -{{\theta}_{0}}\right)-\varepsilon$ and ${{C}_{\text{g}}}=\sin \left(y/r+z/{{R}_{0}}{{q}_{\text{s}}}\right)\approx \sin \left(\theta -{{\theta}_{0}}\right)$ the normal and geodesic curvature, respectively. Here, the average curvature is a constant of order $-\varepsilon$.

The model of kink/peeling modes includes four coupled equations:

$$\begin{eqnarray}&&\nabla \centerdot \mathbf{j}=\nabla \centerdot {{\mathbf{j}}_{\text{pi}}}+\nabla \centerdot {{\mathbf{j}}_{\pi \text{i}}}+\nabla \centerdot {{\mathbf{j}}_{\text{d}}}+\nabla \centerdot {{\mathbf{j}}_{||}}=0\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&e{{n}_{\text{e0}}}{{\partial}_{\text{t}}}{{A}_{||}}={{\partial}_{||}}{{p}_{\text{e}}}-e{{n}_{\text{e0}}}{{\partial}_{||}}\varphi \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\hat{d}}_{\text{t}}}{{p}_{\text{e}}}=0\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{\hat{d}}_{\text{t}}}{{p}_{\text{i}}}=0.\end{eqnarray} \tag{ 11 }$$

Equation (8) is the quasi-neutrality condition. Equation (9) is the parallel generalized Ohm's law, where the electron-ion 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

$$\begin{eqnarray}&&\nabla \centerdot {{\tilde{{\mathbf{j}}}\,}_{\text{pi}}}+\nabla \centerdot {{\tilde{{\mathbf{j}}}\,}_{\pi \text{i}}}\approx \nabla \centerdot \left[{{B}^{-1}}{{m}_{\text{i}}}{{n}_{\text{i0}}}{{d}_{\text{t}}}\left(\mathbf{b}\times {{\tilde{{\mathbf{v}}}\,}_{\text{i}\bot}}\right)\right]=\frac{1}{B{{\omega}_{\text{ci}}}}\nabla \centerdot \left[e{{n}_{\text{i}0}}{{d}_{\text{t}}}\left({{\tilde{{\mathbf{E}}}\,}_{\bot}}-{{\nabla}_{\bot}}{{\tilde{{p}}\,}_{\text{i}}}/{{Z}_{\text{i}}}e{{n}_{\text{i}0}}\right)\right]\approx -{{d}_{\text{t}}}\tilde{{\varpi}}\,-{{\tilde{{v}}\,}_{E\text{r}}}{{\partial}_{\text{r}}}{{\varpi}_{0}}\end{eqnarray} \tag{ 12 }$$

which describes the inertia of electrostatic motion. Here, the last term, ${{\tilde{{v}}\,}_{E\text{r}}}{{\partial}_{\text{r}}}{{\varpi}_{0}}$, is the K–H drive term, responsible for driving the K–H instability [17]. The equilibrium generalized vorticity ${{\varpi}_{0}}={{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{m}_{\text{i}}}{{n}_{\text{i0}}}{{B}^{-1}}{{v}_{\text{iy0}}}\right)$ [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, $\tilde{{\varpi}}\,=e{{n}_{\text{i0s}}}\rho _{\text{s}}^{2}\left({{N}_{\text{i}0}}\nabla _{\bot}^{2}\tilde{{ \Phi }}\,+{{\tau}_{\text{i}}}\nabla _{\bot}^{2}{{\tilde{{P}}\,}_{\text{i}}}-L_{\text{ni}}^{-1}{{\partial}_{\text{r}}}\tilde{{ \Phi }}\,\right)$, in the quantities normalized by the values on the rational surface, ${{N}_{\text{i0}}}\equiv {{n}_{\text{i0}}}/{{n}_{\text{i0s}}}$, $\tilde{{ \Phi }}\,\equiv e\tilde{{\varphi}}\,/{{T}_{\text{e0s}}}$, ${{\tau}_{\text{i}}}\equiv {{T}_{\text{i0s}}}/{{T}_{\text{e0s}}}$, ${{\tilde{{P}}\,}_{\text{i}}}\equiv {{\tilde{{p}}\,}_{\text{i}}}/{{n}_{\text{i0s}}}{{T}_{\text{i0s}}}$, ${{L}_{\text{ni}}}\equiv -{{n}_{\text{i0s}}}/{{\partial}_{\text{r}}}{{n}_{\text{i0}}}$ is the ion density gradient length. The divergence of the diamagnetic drift gives the interchange drive, $\nabla \centerdot {{\tilde{{\mathbf{j}}}\,}_{\text{d}}}=-2{{B}^{-1}}{{R}^{-1}}\left({{C}_{\text{n}}}{{\partial}_{y}}+{{C}_{\text{g}}}{{\partial}_{\text{r}}}\right)\,\left({{\tilde{{p}}\,}_{\text{e}}}+{{\tilde{{p}}\,}_{\text{i}}}\right)\,\approx -2{{n}_{\text{i0s}}}{{T}_{\text{e0s}}}{{\left({{B}_{0}}{{R}_{0}}\right)}^{-1}}$$\left({{C}_{\text{n}}}\text{i}{{k}_{y}}\,\,+\,\,{{C}_{\text{g}}}{{\partial}_{\text{r}}}\right)\,\,\left({{Z}_{\text{e}}}{{\tilde{{P}}\,}_{\text{e}}}\,\,+\,\,{{\tau}_{\text{i}}}{{\tilde{{P}}\,}_{\text{i}}}\right)$ [28], where ${{Z}_{\text{e}}}\equiv {{n}_{\text{e0s}}}/{{n}_{\text{i0s}}}$, ${{\tilde{{P}}\,}_{\text{e}}}\equiv {{\tilde{{p}}\,}_{\text{e}}}/{{n}_{\text{e0s}}}{{T}_{\text{e0s}}}$. The system couples to the shear Alfvén wave through the divergence of the parallel current density, $\nabla \centerdot {{\mathbf{j}}_{||}}=B{{\partial}_{||}}\left({{j}_{||}}/B\right)\approx \left(\text{i}{{k}_{y}}{{B}^{-1}}{{\partial}_{\text{r}}}{{j}_{||0}}-\mu _{0}^{-1}\text{i}{{k}_{||}}\nabla _{\bot}^{2}\right){{\tilde{{A}}\,}_{||}}$, 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

$$\begin{eqnarray}&&\frac{i{{d}_{\text{t}}}\tilde{{\varpi}}\,}{e{{n}_{\text{i0s}}}\rho _{\text{s}}^{2}}-\omega _{\text{A}}^{2}k_{||}^{-1}\nabla _{\bot}^{2}{{\tilde{{A}}\,}_{||}}+{{\omega}_{\text{ci}}}{{k}_{y}}\left[\frac{{{\partial}_{\text{r}\,}}{{j}_{||0}}}{e{{n}_{\text{i0s}}}}{{\tilde{{A}}\,}_{||}}+\frac{{{\partial}_{\text{r}}}{{\varpi}_{0}}}{e{{n}_{\text{i0s}}}}\tilde{{ \Phi }}\,\right. \\ &&\left. -\frac{2}{{{R}_{0}}}\langle {{C}_{\text{n}}}-\text{i}k_{y}^{-1}{{C}_{\text{g}}}{{\partial}_{\text{r}}}\rangle \left({{Z}_{\text{e}}}{{\tilde{{P}}\,}_{\text{e}}}+{{\tau}_{\text{i}}}{{\tilde{{P}}\,}_{\text{i}}}\right)\right]=0\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\left({{N}_{\text{e0}}}\omega -{{\omega}_{\ast \text{e}}}\right){{\tilde{{A}}\,}_{||}}={{k}_{||}}\left({{N}_{\text{e0}}}\tilde{{ \Phi }}\,-{{\tilde{{P}}\,}_{\text{e}}}\right)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\omega {{\tilde{{P}}\,}_{\text{e}}}={{\omega}_{\ast \text{e}}}\tilde{{ \Phi }}\,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\omega {{\tilde{{P}}\,}_{\text{i}}}=-\tau _{\text{i}}^{-1}{{\omega}_{\ast \text{i}}}\tilde{{ \Phi }}\,\end{eqnarray} \tag{ 16 }$$

where ${{\omega}_{\text{A}}}\equiv {{k}_{||}}{{v}_{\text{A}}}$ is the Alfvén frequency, ${{v}_{\text{A}}}\equiv B/\sqrt{{{\mu}_{0}}{{m}_{\text{i}}}{{n}_{\text{i0s}}}}$ is the Alfvén velocity on the rational surface, ${{\tilde{{A}}\,}_{||}}\equiv e{{\tilde{{A}}\,}_{||}}/{{T}_{\text{e0s}}}$, $\langle f\rangle ={\oint}^{}B_{\theta}^{-1}f\text{d}{{l}_{\theta}}/{\oint}^{}B_{\theta}^{-1}\text{d}{{l}_{\theta}}$ is the magnetic surface average operator, ${{N}_{\text{e0}}}\equiv {{n}_{\text{e0}}}/{{n}_{\text{e0s}}}$, ${{\omega}_{\ast \text{e}}}\equiv -{{k}_{y}}{{\partial}_{\text{r}}}{{p}_{\text{e0}}}/eB{{n}_{\text{e0s}}}={{k}_{y}}{{T}_{\text{e0s}}}/eB{{L}_{\text{pe}}}$ is the electron diamagnetic frequency, ${{L}_{\text{pe}}}\equiv -{{p}_{\text{e0s}}}/{{\partial}_{\text{r}}}{{p}_{\text{e0}}}$ is the electron pressure gradient length, ${{\omega}_{\ast \text{i}}}\equiv {{k}_{y}}{{\partial}_{\text{r}}}{{p}_{\text{i0}}}/eB{{n}_{\text{i0s}}}=$$-{{k}_{y}}{{T}_{\text{i0s}}}/eB{{L}_{\text{pi}}}$ is the ion diamagnetic frequency and ${{L}_{\text{pi}}}\equiv -{{p}_{\text{i0s}}}/{{\partial}_{\text{r}}}{{p}_{\text{i0}}}$ is the ion pressure gradient length. The electron curvature drift frequency is defined by ${{\omega}_{\text{de}}}\equiv 2{{k}_{y}}{{T}_{\text{e0s}}}/eBR\approx 2{{c}_{\text{s}}}R_{0}^{-1}{{k}_{y}}{{\rho}_{\text{s}}}=2{{\omega}_{\ast \text{e}}}{{L}_{\text{pe}}}/{{R}_{0}}$.

The first term in equation (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 $\langle {{C}_{\text{n}}}-\text{i}k_{y}^{-1}{{C}_{\text{g}}}{{\partial}_{\text{r}}}\rangle \sim -\varepsilon$.

Substituting equation (15) into (14), we have

$$\begin{eqnarray}&&{{\tilde{{A}}\,}_{||}}=\tilde{{ \Phi }}\,{{k}_{||}}/\omega.\end{eqnarray} \tag{ 17 }$$

Then, substituting equations (15)–(17) into equation (13), noting that $\omega$ is invariant with the plasma minor radius, we arrive at the eigenvalue equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left(\omega {{N}_{\text{i0}}}-{{\omega}_{\ast \text{i}}}\right)\omega \left[{{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\rm{\hat \Phi }\right)-\left(k_{\theta}^{2}+k_{\phi}^{2}\right)\rm{\hat \Phi }\right]-{{\omega}^{2}}L_{\text{ni}}^{-1}{{\partial}_{\text{r}}}\rm{\hat \Phi } \\ \quad -\omega _{\text{A}}^{2}\left[{{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\rm{\hat \Phi }\right)-\left(k_{\theta}^{2}+k_{\phi}^{2}+k_{||}^{-1}{{k}_{y}}{{\partial}_{\text{r}}}L_{\text{s}}^{-1}\right)\rm{\hat \Phi }\right] \\ \quad +\rm{\hat \Phi }{{\omega}_{\text{ci}}}{{k}_{y}}\left[\frac{{{k}_{||}}{{\partial}_{\text{r}}}{{j}_{||0}}}{e{{n}_{\text{i0s}}}}+\frac{\omega {{\partial}_{\text{r}}}{{\varpi}_{0}}}{e{{n}_{\text{i0s}}}}-\frac{2}{{{R}_{0}}}\langle {{C}_{\text{n}}}-\text{i}k_{y}^{-1}{{C}_{\text{g}}}{{\partial}_{\text{r}}}\rangle \left({{Z}_{\text{e}}}{{\omega}_{\ast \text{e}}}-{{\omega}_{\ast \text{i}}}\right)\right]=0. \end{array}\end{eqnarray} \tag{ 18 }$$

Note that the radial gradient and curvature of ${{\omega}_{\ast \text{i}}}$ have been neglected. This is a quadratic equation in $\omega$. The coefficients are

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{c}_{1}}={{N}_{\text{i0}}}\left[{{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\rm{\hat \Phi }\right)-\left(k_{\theta}^{2}+k_{\phi}^{2}\right)\rm{\hat \Phi }\right]-L_{\text{ni}}^{-1}{{\partial}_{\text{r}}}\rm{\hat \Phi } \\ {{c}_{2}}=-{{\omega}_{\ast \text{i}}}\left[{{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\rm{\hat \Phi }\right)-\left(k_{\theta}^{2}+k_{\phi}^{2}\right)\rm{\hat \Phi }\right]+\rm{\hat \Phi }{{\omega}_{\text{ci}}}{{k}_{y}}{{\partial}_{\text{r}}}{{\varpi}_{0}}/e{{n}_{\text{i0s}}} \\ {{c}_{3}}=-\omega _{\text{A}}^{2}\left[{{r}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\rm{\hat \Phi }\right)-\left(k_{\theta}^{2}+k_{\phi}^{2}+k_{||}^{-1}{{k}_{y}}{{\partial}_{\text{r}}}L_{\text{s}}^{-1}\right)\rm{\hat \Phi }\right] \\ +\rm{\hat \Phi }{{\omega}_{\text{ci}}}{{k}_{y}}\left[{{k}_{||}}{{\partial}_{\text{r}}}{{j}_{||0}}/e{{n}_{\text{i0s}}}-2R_{0}^{-1}\langle {{C}_{\text{n}}}-ik_{y}^{-1}{{C}_{\text{g}}}{{\partial}_{\text{r}}}\rangle \left({{Z}_{\text{e}}}{{\omega}_{\ast \text{e}}}-{{\omega}_{\ast \text{i}}}\right)\right]. \end{array}\end{eqnarray} \tag{ 19 }$$

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\rm{\hat \Phi }$ and calculate the radial integral, ${{C}_{j}}={\int}_{{{r}_{a}}}^{{{r}_{b}}}r\rm{\hat \Phi }{{c}_{j}}\text{d}r$, over a radial range, ${{r}_{a}}\leqslant r\leqslant {{r}_{b}}$, covering the whole pedestal region, e.g. ${{r}_{a}}=0.85a$ and ${{r}_{b}}=a$ with $a$ the plasma minor radius at the last closed flux surface. Then, we obtain a quadratic equation in $\omega$

$$\begin{eqnarray}&&{{C}_{1}}{{\omega}^{2}}+{{C}_{2}}\omega +{{C}_{3}}=0.\end{eqnarray} \tag{ 20 }$$

Solving this equation, we finally obtain the dispersion relation, $\omega \equiv {{\omega}_{\text{r}}}+\text{i}\gamma =\left(-{{C}_{2}}\pm \sqrt{C_{2}^{2}-4{{C}_{1}}{{C}_{3}}}\right)/2{{C}_{1}}$. Its imaginary part is the required linear growth rate.

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 toroidal-mode-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.35 s as indicated by the vertical line in the plots. Some parameters at this time point, which have been used in the following analysis, are toroidal magnetic field ${{B}_{\phi 0}}=2.009\,\text{T}$ on the magnetic axis, major radius ${{R}_{0}}\text{=1}\text{.788} ~\text{m}$, minor radius $a=0.582\,\text{m}$, plasma internal inductance ${{l}_{\text{i}}}=0.791$ and poloidal beta ${{\beta}_{\text{p}}}=0.968$.

Zoom In Zoom Out Reset image size

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 toroidal-mode-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–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⁻¹, 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.

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}_{\text{e0}}}$, temperature ${{T}_{\text{e0}}}$ and ion density ${{n}_{\text{i0}}}$, temperature ${{T}_{\text{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}_{\text{i0}}}$ profile which is fitted to a spline model in the whole region. The electron density ${{n}_{\text{e0}}}$ and temperature ${{T}_{\text{e0}}}$ are obtained from Thomson scattering (TS) diagnostics. Impurity (carbon) density ${{n}_{\text{c}0}}$, ion temperature ${{T}_{\text{i0}}}$ and radial electric field ${{E}_{\text{r0}}}$ are obtained from CER diagnostics. The radial electric field profiles including $0.5{{E}_{\text{r0}}}$ (red), ${{E}_{\text{r0}}}$ (black) and $2{{E}_{\text{r0}}}$ (blue) are shown in figure 7(c) and their corresponding shearing rates $0.5{{S}_{v}}$ (red), ${{S}_{v}}$ (black) and $2{{S}_{v}}$ (blue) are shown in figure 7(d). The ion density ${{n}_{\text{i0}}}$ profile shown in figure 7(e) is calculated by ${{n}_{\text{i0}}}={{n}_{\text{e0}}}-6{{n}_{\text{c0}}}$, assuming carbon is the only impurity. The electron (ion) pressure ${{p}_{\text{e0}}}$ (${{p}_{\text{i0}}}$) shown in figure 7(g) is taken as the product of the electron (ion) density and temperature, respectively.

Zoom In Zoom Out Reset image size

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 ${{\omega}_{\text{A}}}$ (blue), ${{\omega}_{\ast \text{e}}}$ (black), ${{\omega}_{\ast \text{i}}}$ (red) and ${{\omega}_{E}}$ (blackish green) profiles, where the toroidal mode number n is set to 1.

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, ${{\hat{\varphi}}_{t=0}}(x)={{\hat{\varphi}}_{\text{A}}}H(-x)\exp \left(-{{x}^{2}}/2{{\lambda}^{2}}\right)$, is a half-Gaussian distribution function centering at the rational surface with the radial decay length, $\lambda$, the only undetermined parameter. We search in the parameter space of $\lambda$ to maximize the linear growth rate, $\gamma$, 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 $\lambda$ 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, $\rm{\hat \Phi }{{\omega}_{\text{ci}}}{{k}_{y}}{{\partial}_{\text{r}}}{{\varpi}_{0}}/e{{n}_{\text{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, ${{\hat{\omega}}_{\text{r}}}={{\omega}_{\text{r}}}+{{\omega}_{E}}$, 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, $\gamma$, 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 ${{\gamma}_{0}}=1.34\,\text{MHz}$. The real frequency in the plasma-moving frame of reference, $f={{\hat{\omega}}_{\text{r}}}/2\pi$, in figure 8(b), shows a nearly linear dispersion relation in the ion diamagnetic direction (negative). The radial decay length, $\lambda$, of the mode electrostatic potential structures in figure 8(c) is in the range of 3–5 mm, peaking with $\lambda =5.3\,\text{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 distribution function, ${{\hat{\varphi}}_{t=0}}(x)={{\hat{\varphi}}_{\text{A}}}H(-x)\exp \left(-{{x}^{2}}/2{{\lambda}^{2}}\right)$. The radial width of the mode is generally comparable with the pedestal width.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the radial expansion is taken into account, the radial amplitude profile of the mode electrostatic potential structures, $\hat{\varphi}(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, $\lambda$, 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

$$\begin{eqnarray}&&\hat{\varphi}(r)={\int}_{{{r}_{\text{s}}}}^{r}\left\{{\int}_{{{r}_{\text{s}}}}^{r}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\hat{\varphi}\right)\exp \left[-\text{i}{{k}_{y}}\left({{v}_{Ey0}}-{{v}_{Ey0\text{s}}}\right)t\right]\text{d}r+C\right\}{{r}^{-1}}\text{d}r+{{\hat{\varphi}}_{\text{s}}}.\end{eqnarray} \tag{ 21 }$$

The real and imaginary parts of $\hat{\varphi}(r)$ for n  =  3 mode at $t=0.2\,\mu \text{s}\approx 0.4{{\tau}_{\text{A}}}$ are shown in figure 9, where ${{\tau}_{\text{A}}}\equiv \omega _{\text{A}}^{-1}$ is 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\,\mu \text{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, $\arctan \left[{{\hat{\varphi}}_{Im}}(r)/{{\hat{\varphi}}_{\operatorname{Re}}}(r)\right]$, 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}_{\text{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 $0.5{{E}_{\text{r0}}}$, $n=3$ with ${{E}_{\text{r0}}}$ and $n=2$ with $2{{E}_{\text{r0}}}$. The unstable modes with ${{E}_{\text{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}_{\text{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}_{\text{r0}}}$ is shifted radially, the E  ×  B 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}_{\text{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].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 differential advection of mode vorticity by sheared E  ×  B flows modifies the radial distributions of the mode generalized vorticity, $\tilde{{\varpi}}\,\equiv {{m}_{\text{i}}}{{e}^{-1}}{{B}^{-2}}{{\nabla}_{\bot}}\centerdot \left(e{{n}_{\text{i0}}}{{\nabla}_{\bot}}\tilde{{\varphi}}\,+{{\nabla}_{\bot}}{{\tilde{{p}}\,}_{\text{i}}}\right)$ [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. $\nabla \centerdot {{\mathbf{j}}_{||0}}={{\tilde{{B}}\,}_{\text{r}}}{{\partial}_{\text{r}}}\left({{j}_{||0}}/B\right)\approx \text{i}{{k}_{y}}{{\tilde{{A}}\,}_{||}}{{B}^{-1}}{{\partial}_{\text{r}}}{{j}_{||0}}$, which in turn is proportional to ${{k}_{||}}$, i.e. ${{\tilde{{A}}\,}_{||}}=\tilde{{\varphi}}\,{{k}_{||}}/\omega$, as shown in [28]. In ideal-MHD, the field lines are frozen in the fluid, the radial expansion of the mode electrostatic potential, $\tilde{{\varphi}}\,$, will lead to field line bending radially, ${{\tilde{{B}}\,}_{\text{r}}}={{\partial}_{y}}{{\tilde{{A}}\,}_{||}}$, 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}_{||}}\approx -{{k}_{y}}x/{{L}_{\text{s}}}$, increases with x. Because of the field line bending, the equilibrium current-density gradient has a projection along the magnetic field lines, ${{\tilde{{B}}\,}_{\text{r}}}{{\partial}_{\text{r}}}\left({{j}_{||0}}/B\right)$. 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}^{-1}}{{\partial}_{\text{r}}}\left(r{{\partial}_{\text{r}}}\hat{\varphi}\right)$, 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, $\tilde{{\varpi}}\,$. 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>\omega _{\text{s}}^{-1}$, the radial expansion would appear periodically with a frequency of ${{\omega}_{\text{s}}}=\lambda {{k}_{y}}{{S}_{v}}$, when the mode structure shifts beyond a half-wavelength poloidally. However, this is prevented in our case, since the linear growth rates are usually much larger than the shear frequency, $\gamma \gg {{\omega}_{\text{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 electromagnetic 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 mode-mode 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 collapses—the large ELMs, in subsequent nonlinear development.

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–47]. It was suggested that the E  ×  B flow shear could shift the peak amplitude of the mode eigenfunction, $\hat{\varphi}(x)={{\hat{\varphi}}_{\text{A}}}\exp \left[-{{(x-\xi )}^{2}}/2{{\lambda}^{2}}\right]$, radially away from the rational surface by a distance $\xi$ in a plasma with significant equilibrium density and temperature gradients [43]. The radial distance of shift generally increases with the shear frequency, $\xi /\lambda \approx {{\omega}_{\text{s}}}/\gamma =\lambda {{k}_{y}}{{S}_{v}}/\gamma \ll 1$. 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.