Recently, the stationary high confinement operations with improved pedestal conditions have been achieved in DIII-D [K. H. Burrell *et al.*, Phys. Plasmas **23**, 056103 (2016)], accompanying the spontaneous transition from the coherent edge harmonic oscillation (EHO) to the broadband MHD turbulence state by reducing the neutral beam injection torque to zero. It is highly significant for the burning plasma devices such as ITER. Simulations about the effects of **E** × **B** shear flow on the quiescent H-mode (QH-mode) are carried out using the three-field two-fluid model in the field-aligned coordinate under the BOUT++ framework. Using the shifted circular cross-section equilibriums including bootstrap current, the results demonstrate that the **E** × **B** shear flow strongly destabilizes low-n peeling modes, which are mainly driven by the gradient of parallel current in peeling-dominant cases and are sensitive to the *E _{r}* shear. Adopting the much more general shape of

**E**×

**B**shear ($\omega E=Er/RB\theta $) profiles, the linear and nonlinear BOUT++ simulations show qualitative consistence with the experiments. The stronger shear flow shifts the most unstable mode to lower-n and narrows the mode spectrum. At the meantime, the nonlinear simulations of the QH-mode indicate that the shear flow in both co- and counter directions of diamagnetic flow has some similar effects. The nonlinear mode interaction is enhanced during the mode amplitude saturation phase. These results reveal that the fundamental physics mechanism of the QH-mode may be shear flow and are significant for understanding the mechanism of EHO and QH-mode.

The high-performance regime (H-mode) is regarded as one of the most promising operational scenarios to achieve sufficient fusion gain in the present tokamaks and for future burning plasma devices such as ITER, due to the excellent high energy confinement at high density. However, the edge-localized modes (ELMs), the quasi-periodic relaxation of pedestal because of the huge free energy held by steep gradients of pressure and parallel current according to the peeling-ballooning (PB) model,^{2} lead to large transient energy and particle flux to the plasma facing components (PFCs) and thereby can result in unacceptable damage to divertor plates and first walls. Considering the detrimental impacts of ELMs on burning plasmas, significant efforts are made on ELM mitigation and control both in theory and in experiments. Providing enough particle transport is the general approach to ELM mitigation and impurity control. The methods can be classified into three main categories:^{3} (1) the intrinsic ELM-free regime (e.g., QH-modes^{4}), (2) triggering frequent smaller ELM pulse (e.g., pellet pacing^{5}), and (3) external means such as resonant magnetic perturbations (RMPs).^{6}

Quiescent H-mode (QH mode) is an ELM-free high confinement plasma and accompanies with the edge harmonic oscillation (EHO), which is conjectured to be a kink/peeling mode driven by current and rotational shear at conditions near but just below the peeling instability limit boundary of the peeling-ballooning model.^{7} The EHO can provide additional energy and particle transport to make it possible to operate steadily for long duration. The features of high energy confinement without ELMs and accessibility at low pedestal collisionality with strong shaping make it attractive for future burning plasma devices. In the edge region, the bootstrap current is dominant and is roughly proportional to the pressure gradient, but is decreased by collisions, which are proportional to the third power of density at a given pressure.^{7} That's to say, the plasma current plays an important role in the QH-mode operation.

The QH-mode was originally achieved on DIII-D in discharges with neutral beam injection (NBI) in the counter direction of plasma current and subsequently also observed in ASDEX-U,^{8} JET,^{9} and JT-60U.^{10} More recently, the QH-mode has been obtained with a strong edge rotation in co-direction of plasma current. It is a significant confirmation of the peeling-ballooning mode theory, which predicts that the QH-mode should exist with both co- and counter plasma rotations and only depend on the strength of the rotation shear.^{11} Burrell *et al.*^{1} report the discovery of stationary operation of the QH-mode in DIII-D by a net-zero NBI torque with impressively improved pedestal condition, while accompanying the spontaneous transition from EHO to broadband MHD turbulence by reducing the NBI torque to zero without external 3D magnetic fields to provide an electromagnetic torque like previous experiment.^{4} Computational simulations of the QH-mode reveal more and more features and improve our understanding of the QH-mode. The results from the M3D-C1 code demonstrate that EHOs can be destabilized in principle with rotation in either direction.^{12} The simulations with JOREK code show that kink peeling modes are the main unstable modes and saturated in the nonlinear phase in QH-mode plasmas.^{13} Consistent with experimental observations during the QH-mode, the simulation results using the NIMROD code illustrate that the particle transport is larger relative to the thermal transport in the QH-mode.^{14} It has great realistic significance for the future fusion devices.

Progress in the QH-mode discharge has led to emerging understanding the role of **E **×** B** shear played in the QH-mode operations. During these operations, the pedestal conditions are associated with the edge density changes and altered **E **×** B** shear. In the ELMy H-mode operation, the **E **×** B** shear flow can stabilize the high-n mode and induce the Kelvin-Helmholtz instabilities as well. The performance of the plasma is the result of the competition between these two effects in these ballooning dominant cases.^{15} Here we study the physics mechanism of EHO's features and the role of **E **×** B** shear during the QH-mode.

In our simulations, we use the reduced 3-field nonlinear two-fluid MHD model^{16} in shifted circular cross-section geometry under the BOUT++ framework. This equation set based on the P-B model can be simplified from the complete set of BOUT two-fluid equations^{17} within the flute reduction $k\u2225/k\u22a5\u226a1$ including additional effects of hyper-resistivity,^{16} K-H term.^{15} In our simulations, the ion diamagnetic drift and gyroviscosity are included to get the correct finite Larmour radius (FLR) effect.^{18} The model under the incompressible assumption evolves the perturbations of vorticity ϖ, pressure *P*, and the parallel vector potential $A\u2225$

in which $d/dt=\u2202/\u2202t+VE\xd7B\xb7\u2207,\u2009\kappa 0=b0\xb7\u2207b0,\u2207\u2225F=B\u2202\u2225(F/B),\u2009\u2202\u2225=(b0+b\u0303)\xb7\u2207$ and $b\u0303=\u2207A\u2225\xd7b0/B$. As the equilibrium electric potential is divided into two parts $\Phi 0=\Phi dia0+\Phi V0$, the **E** ×**B** convection flow is split into three parts accordingly $VE\xd7B=VE,dia0+VE,V0+V1$: (I) $VE,dia0=b0\xd7\u2207\Phi dia0/B$, which balances the ion diamagnetic flow $b0\xd7\u2207Pi0/(Zeni0B)$; (II) $VE,V0=b0\xd7\u2207\Phi V0/B$, the net equilibrium flow; (III) $V1=b0\xd7\u2207\varphi /B$, which denotes the perturbed flow and $\varphi $ is the perturbed electric potential. It is very convenient for the study of the effects of equilibrium flow

on the MHD instabilities as we assume that $\Phi V0=\Phi V0(\psi )$ is a flux function.^{15} As the EHOs are thought to be the saturated low-n kink/peeling mode driven by current and the rotational shear, the simulations adopting the experimental shear profiles in peeling dominant situations are vitally important for understanding the mechanism of QH-modes and studying the EHOs.

In our simulations, two serials of equilibrium profiles in JET-like Tokamak geometry with shifted circular cross-section, as shown in Fig. 1, are used in our simulations. The grid size in the radial and poloidal direction is $nx\xd7ny=516\xd7\u200964$ and in the z-direction is *n _{z}* = 16 (for linear simulations) and 128 (for nonlinear simulations), which is sufficient considering the results of numerical convergence study and the computational resource consuming for these cases. The radial simulation domain of the normalized poloidal flux is $\psi n=[0.1,\u20091.4]$. Basic parameters are as follows: minor radius $a=1.2$ m, major radius $R0=3.4$ m, safety factor at 95% flux $q95=2.8$, and toroidal magnetic field $B0=2$ T. Both dimensionless quantities Lundquist number and hyper-Lundquist number are fixed at $S=\mu 0R0VA/\eta =108$ and $SH=\mu 0R03VA/\eta H=1015$, respectively, in the simulations. The linear results demonstrate that the resistivity has a little impact on the instabilities under these equilibriums and the hyper-resistivity, which is also known as anomalous electron viscosity and was studied previously,

^{16}is small enough and almost has no effects on the linear growth rate of these cases but is important for the numerical reason in the nonlinear simulation. The toroidal equilibrium module (TEQ) in the CORSICA code

^{19}is utilized to modify the equilibrium and to calculate the bootstrap current according to the Sauter bootstrap current model.

^{20}The equilibrium ‘Serial 1’ in Fig. 1 is with the same pressure profile but various densities (central density $ne0=1,9,20\xd71019m\u22123$) and current profiles, it includes both low-n and high-n modes unstable cases and is much more suitable to investigate the effects of current drive and the

*E*shear. The equilibrium density profiles are given by $ne(\psi n)=ne0*P\u030200.3(\psi n)$, assuming that the ion and electron temperatures are equal ($Te(\psi n)=Ti(\psi n)=P0(\psi n)/2ne(\psi n)$). Due to the QH-mode usually operating near but just below the peeling-ballooning mode stability boundary that sets the ELM limit, the equilibrium ‘Serial 2’ (with green markers in Fig. 1) with much lower pressure and current is adopted to avoid ELMs throughout the nonlinear simulations for the further study of the QH-mode and EHO.

_{r}The P-B model predicts that the edge instabilities are mainly driven by the steep gradient of pressure and parallel current in pedestal. Since the EHO is thought to be the saturated low-n kink/peeling mode and can be destabilized by current and rotational shear, the linear simulations without net flow (K-H term) and current drive term are carried out to check the effects of current drive, using the equilibrium Serial 1 and including the diamagnetic effects. To clearly show the effects of the current drive, we simply chose cases with large growth rate under condition $Er0=2.86Er,dia0$. From Fig. 2 it is evident that after switching off the current drive term (the dashed curves), the growth rate of low-n modes is significantly reduced, while there are much weaker impacts on the high-n modes. After turning off both the current drive and the initial $Er0$ drive, all the low-n (<10) modes are stable. According to the radial force balance, the radial electric field can be expressed as

As the $Er0$ is inversely proportional to the density, the low-n modes (n = 1–10) are much more unstable for low density cases than the high density cases. Figure 2 also demonstrates that the destabilizing effects of current drive on the low density cases are much stronger than the high-density cases in which the ballooning modes are in dominant. From the above metioned linear simulation results, we can easily conjecture the fact that the low-n kink/peeling modes are mainly driven by the currents and the radial electric field shows, the gradient of pressure drives the high-n instabilities, which confirms the theory of peeling-ballooning model to a great extent.

According to the radial force balance Eq. (8), the toroidal rotation will induce the significant change of the radial electric field. In experiments, the neutral beam injection in both co- and counter directions of the plasma current has been used to alter the toroidal velocity $V\varphi i$. The author demonstrated the destabilizing effect of net flow on high-n modes in the ballooning dominant cases.^{15} Here, we qualitatively study the influence of *E _{r}* on the linear instabilities by amplifying

*E*under the peeling dominant regime. As illustrated in Fig. 3, the low density case ($ne0=1\xd71019\u2009m\u22123$, dashed), dominated by low-n peeling modes, is more sensitive to

_{r}*E*and strongly destabilized by

_{r}*E*. But for the high density case ($ne0=20\xd71019\u2009m\u22123$, solid), in which the high-n ballooning modes are in dominant, the

_{r}*E*shows little destabilizing effects on high-n modes. Refer to the results of Ref. 15, the

_{r}*E*shear plays a much more important role than the amplitude of

_{r}*E*on the peeling-ballooning instabilities.

_{r}To date, the QH-mode operations have been achieved by NBI in counter^{21} and co-directions^{11} of plasma current and even created with a net-zero NBI torque^{1} reported recently. The **E **×** B** shear is varied contingently on the combined effects of the power and directions of NBI as indicated from Eq. (8). To be much closer to the actual experimental condition, two general shapes of *E _{r}* profiles with an amplitude close to the DIII-D QH-mode discharge #163518 (∼80 kV/m at 2350 ms, ∼50 kV/m)

^{22}are adopted in the simulations. Larger amplitudes of $E\xd7B$ are chosen to illustrate the effect of

*E*shear flow on the instability, as can be seen in Fig. 4, which simulate the

_{r}*E*profiles during the different stages of the QH-mode discharge: the coherent EHO state (high rotation, “high”) was replaced by the broadband MHD turbulence (low rotation, “low”) during the low torque times with an NBI torque ramp.

_{r}To start the simulations without ELMs, the equilibrium Serial 2 displayed in Fig. 1 with much lower pressure and current profile was used in the following linear and nonlinear simulations.

Beginning with the linear cases under the different strength of shear flows, it is suitable to use the amplitude of the net shear flow as the control variable while keeping the similar shape to the low rotation profile as shown in Fig. 4. It is manifest from Fig. 5 that the diminishing net flow, which means the smaller Hahm-Burrell **E **×** B** shearing rate^{23} $\omega E\xd7B=(RBp)2B\u2202\u2202\psi (ErRBp)$, leads to the reduction of the most unstable mode, shifts the peak of $\gamma (n)$ to higher toroidal mode number, and widens the mode spectrum. These linear results are qualitatively consistent with the recent DIII-D experiment results^{1} that the coherent EHO with a few isolated dominant low-n modes is replaced by broadband MHD turbulence with a wider mode spectrum by deliberately reducing the NBI torque slightly to zero during the QH-mode operation.

The bi-spectral analysis is a very useful technique to study the three-wave nonlinear interactions.^{24} The total bi-coherence (bisum) measures the nonlinear interaction of one mode with all the other modes. As is illustrated in Fig. 6(a), in a nonlinear simulation without net flow, the temporal evolution of the mode amplitude entered into the saturation phase after about 1000*τ _{A}*. Figure 6(b) shows the nonlinear interaction among modes at different times marked by same colored dashed lines in Fig. 6(a). As expected, in the saturated phase ($>1000\tau A$) the nonlinear interaction is stronger than the linear phase. It also shows a clear downshift of the peak of the total bi-coherence spectrum from n = 4,5 in linear phase to n = 2,3 in nonlinear phase. For different amplitudes of $E\xd7B$ shear flow, the same trend (bi-sum gap between linear and nonlinear stages) exists.

Figure 7 demonstrates the nonlinear mode spectrum with different strength net flow in the co-direction of the diamagnetic flow. The labels “low” and “high” refer to the profiles used in the simulation, as shown in Fig. 4. The “0.5low” means we reduce the net flow profile of “low” case by half. As can be seen from Fig. 7, in the nonlinear phase the net flow shows the following effects: (1) weaker net flow induces the broader mode spectrum; (2) there are some dominant modes ($n\u223c3$) and enhanced for the cases with stronger net flow (blue and red), but non-obvious dominant mode for the case without net flow (black). The PB-model is used to study the ELM and has recently extended to study the QH-mode in which the dominant toroidal mode number of the EHOs is n = 1–5 typically.^{7} As shown in Fig. 8, for the contour plot of perturbed potential fluctuation versus frequency and time, there are dominant frequencies which the toroidal mode number is 1 in the cases with low and high net flow, and broad turbulence spectrum in case without net flow. The cases with net flow in counter direction of diamagnetic flow have similar features shown in Figs. 7 and 8, while the mode numbers of dominant frequencies are n = 1 and 2 in case with high net flow and it already enters into the turbulence state in case with low net flow. That is because the shear flow can suppress the turbulence, the high n modes are reduced by the net flow as can be seen from Fig. 7. The nonlinear interaction makes the dominant mode number slightly different from the linear results. The dominant mode numbers are sensitive to the *E _{r}* profile which is changed when reversing the net flow direction. It confirms that the QH-mode can be achieved by the net flow in both co- and counter directions of the diamagnetic flow as observed in experiments. It is the strength of shearing, which plays an important role in the QH-mode operation but not the direction of the flow. These nonlinear simulations show qualitative agreement with the recent DIII-D experiments results: when decreasing the flow shear, the mode spectrum becomes broader and the amplitude of single-dominant modes reduced. These results demonstrate the transition from coherent EHO to broadband MHD turbulence state in a certain extent.

To study the mechanism of the nonlinear saturation, we compared the contribution of the quasi-linear part ($\u27e8V1\xb7\u2207Pdc\u27e9rms$) and the nonlinear turbulence part ($\u27e8V1\xb7\u2207P\u27e9rms\u2212\u27e8V1\xb7\u2207Pdc\u27e9rms$) of the convection term ($V1\xb7\u2207P$) in the nonlinear simulations. Here the subscript “*dc*” means surface average and the operation $\u27e8\u27e9rms$ denotes the root mean square of the variable. The results show that the net flow makes the nonlinear turbulence part of the convection term dominant in the nonlinear saturated phase in case with net flow in the co-direction of the diamagnetic flow. In the early nonlinear phase ($700\tau A\u22121000\tau A$), the quasi-linear part is in dominant. After entering into the fully nonlinear phase ($>1000\tau A$), the two parts are split by increasing net flow in co-direction. However, these phenomena are not observed in simulations with net flow in the counter direction of the diamagnetic flow. For these three cases with co-direction net flow, the system enters into saturated phase after about 1000*τ _{A}*, as shown in Fig. 6(a). What roles the quasi-linear and turbulent parts play needs further study.

Both linear and nonlinear simulations are carried out to study the EHO and QH-mode using a 3-Field model under the BOUT++ framework. It is obvious from the linear results that the low-n modes are mainly driven by the parallel current, the radial electric field also destabilizes the low-n modes, and the steep pressure gradient leads to the high-n ballooning instabilities. It also can be inferred from the linear results that the low-n modes are very sensitive to the radial electric field and the net flow shows the destabilizing effects on the low-n kink modes as the theory demonstrates that the QH-mode can be disturbed by rotational shear and current.^{7}

In this manuscript, although we do not simulate the DIII-D discharges directly, the results do show qualitative agreement in trend correlating to the DIII-D experiments which achieved the transition from EHO with dominant toroidal mode to the broadband turbulence state when decreasing the NBI torque to zero. In the cases without net flow, the mode spectrum becomes broader without dominant toroidal modes and the time-frequency analysis shows no frequency in dominance, and these results demonstrate the features of the turbulence state.

Most importantly, the simulations with net flow in both co- and counter directions of the diamagnetic flow are constructed and confirm that both signs of net flow have similar effects on the QH-mode. The state of operation mainly depends on the strength of net flow. These are quite consistent with the experimental facts that the QH-mode can be achieved by the NBI in both co- and counter directions of current. We could come to the conclusion that the shear flow may be the fundamental physics mechanism in the achievement of the QH-mode based on the both experiments and simulation results.

In our simulation, we use the shifted circular cross section geometry and just consider the *E _{r}*'s effect on the instabilities now. The shaping plays an important role

^{1}on the creation of the wide pedestal QH-mode which has only been seen in double null discharges and changes the stability boundary.

^{7}The achievement of the low torque QH-mode experimentally needs strict conditions of density, current profile, and shaping

^{25}and should take the effects of wall conditions,

*E*shape in the account. Hence, for further study of the mechanism of EHO and QH-mode, the Tokamak geometry should be taken into account in the simulations. The additional energy and particle transport without ELM are one of the key features provided by EHO in QH-mode discharges. It is significant to take all these factors into account in the simulations. Whether there is a relationship between energy and particle transports needs more accurate models like 6-field model or the Gyro-Landau fluid model as they are developed to evolve density, temperature and pressure of electron and ion, parallel velocity separately, instead of just evolving the total pressure in the 3-field model. What is the role of net flow in co-direction of diamagnetic flow in the QH-mode acquires further study as well.

_{r}The authors would like to acknowledge G. Q. Li, N. Yan, T. Y. Xia, and Z. B. Guo for the constructive comments. The authors would also like to thank K. H. Burrell and X. Chen for the useful discussions particularly about the *E _{r}* profiles. This work was performed under the auspices of the U.S. DoE by LLNL under Contract No. DE-AC52-7NA27344 and was supported by the ITER-China Program (2013GB111000, 2013GB112006), NSFC (11261140326 and 11375053) and LLNL-JRNL-721760.