Global gyrokinetic simulation of edge coherent mode in EAST

Linear and nonlinear simulations are carried out for the edge coherent mode (ECM) using the global gyrokinetic code GEM based on the EAST experimental parameters. The linear simulation results show that ECM is an electrostatic mode with dominant toroidal mode number n = 18 and frequency about 48 kHz, and propagates along the direction of electron diamagnetic drift, which are consistent with the experimental results. In addition, the density and electron temperature gradients destabilize the mode, while the collision stabilizes the mode. The nonlinear simulation results show that the saturated particle and heat fluxes induced by ECM are mainly due to the perturbed electrostatic ExB drift, and the fluxes of electrons and ions are almost equal. The ECM drives significant outward particle and heat fluxes, thus greatly promoting the maintenance of the long pulse H-mode. The Fourier decomposition of fluxes and potentials demonstrate that the intermediate-n modes of n = 14, 18 grow fastest in the linear phase, while in the nonlinear saturation phase, the low-n modes such as n = 4, 6 dominate and the fluxes are mainly contributed by the mode of n = 10. It is found that zonal flow is not the dominant saturation mechanism of the turbulence. The inverse spectral cascade of turbulence is inevitably observed in the nonlinear saturation process, indicating that it is a more universal turbulence saturation mechanism. It is also found that radial electric field can greatly reduce the turbulence intensity and transport level. From the analyses of frequency and transport channels, it can be concluded that ECM appears to be the collisionless trapped electron mode.


Introduction
High-confinement mode (H-mode) [1] is the basic operation mode for future fusion reactor. The pedestal structure formed in the H-mode plays a key role in the overall confinement of plasma. Due to the strong stiffness of the core plasma, a higher pedestal top helps to improve the core plasma confinement performance. However, the large pressure gradient in the pedestal can drive a number of instabilities [2], including microinstabilities like electrostatic ion temperature gradient modes (ITG) [3,4], electron temperature gradient mode (ETG) [5], trapped electron mode (TEM) [6], and electromagnetic instabilities like kinetic ballooning modes (KBM) [7][8][9], micro tearing modes (MTM) [5,10], and magnetohydrodynamic instabilities like peeling-ballooning mode induced edge localized modes (ELM) [11,12]. These instabilities enhance the heat and particle fluxes transported to the first wall and divertor, severely damage the boundary materials, and affect the evolution of the pedestal structure. Therefore, studying the instabilities in the pedestal is of great significance for realizing high-confinement operation.
So far, various quasi coherent modes have been observed in most experimental devices with frequencies ranging from tens to hundreds of kHz in the plasma pedestal region in H-mode regime. In DIII-D, a set of high frequency coherent modes with the frequency range of 80-250 kHz have been observed in high density quiescent H-mode plasma, which qualitatively matches the characteristics of KBMs [7]. Besides, two distinct bands of fluctuations with different frequencies and propagation directions have been reported, which are lower frequency band  advects in the ion-diamagnetic drift direction and a higher frequency band (200-400 kHz) advects in the electron diamagnetic drift direction [8]. In Alcator C-Mod [9], a quasi coherent fluctuation (QCF) propagated along the direction of ion diamagnetic drift in the plasma frame with frequency of about 300 kHz and toroidal mode number n = 10 has been observed. It is found that the appearance of QCF is related to the evolution of pedestal electron temperature, and it can enhance the edge transport and limit the growth of edge temperature. QCFs have also been observed in DIII-D with a frequency of tens to more than 100 kHz, and its saturation amplitude is consistent with the pedestal temperature and pressure gradient [13]. Modeling and observations indicate that QCFs are consistent with KBMs but can not exclude the MTM as subdominant. In ASDEX upgrade, the high frequency magnetic fluctuations with frequencies above 200 kHz have been found to correlate with the stagnation of the electron temperature pedestal gradient, which are also observed on both the low field side and high field side and might be attributed to a strong peeling part or the presence of non-adiabatic electron response [14]. In addition, a low-frequency magnetic fluctuation with the frequency range of 4 kHz-12 kHz has been observed and may be related to the evolution of pedestal temperature [15]. The density and magnetic fluctuations with frequency of 60-70 kHz have been observed in JET, which indicates that the density fluctuation is dominated by ITG, and the magnetic fluctuation is likely to be MTM [16]. In HL-2A [17], a quasi-coherent mode with the frequency range of 50-100 kHz has been observed, which appears during the ELM free period after L-H transition or inter ELM period at the pedestal, and is driven by the density gradient. In EAST, a lot of coherent modes have also been observed, including the edge electrostatic coherent mode (ECM) with the frequency range of 20-90 kHz [18], broadband fluctuations with the frequency range of 0-200 kHz [19], and magnetic coherent mode with the frequency range of 20-60 kHz [20]. These fluctuations can exist alone or simultaneously during the inter-ELM period [21], and propagate in the electron diamagnetic drift direction in the laboratory frame or plasma frame and correlate with the pedestal density or pressure evolution. Lately, a variety of experimental studies on pedestal localized fluctuations inter-ELM period across several tokamaks have been reviewed and summarized [22]. It is found that the profile recovery phases can be assorted into three categories, so as the pedestal fluctuations.
Meanwhile, a lot of simulation analyses have been done based on these experiments to identify the underlying generation mechanism of these quasi coherent modes and impact on transport as well as pedestal structure formation. Local linear gyrokinetic analysis based on MAST reveals that the dominant microinstability is KBM in the steep gradient region of pedestal, and MTM near the core region [23]. High-n KBM and intermediate-n kinetic version of peeling-ballooning mode near threshold conditions have been found in global electromagnetic gyrokinetic simulations of the edge pedestal of two DIII-D H-mode discharges [24]. Moreover, the linear gyrokinetic simulations of electrostatic [25] and electromagnetic instabilities [26] in DIII-D tokamak pedestal are conducted, which indicates that the dominant instability at the top of the pedestal is ITG, while in the maximum gradient location the most unstable mode is KBM and TEM. Nonlinear gyrokinetic simulations of the JET-ILW pedestal indicate that the bulk of the transport in the steep gradient region is driven by MTM, where KBM is found to be in the second stability regime [4,27]. Adding neoclassical and ETG driven transport, the MTM can reproduce the experimental power balance across most of the pedestal quantitatively. Linear local gyrokinetic simulations based on EAST show that ECM may be the dissipative trapped electron mode (DTEM) [18,28], while the nonlinear fluid simulation indicates that ECM is mainly driven by drift Alfven wave [29]. Recently, a new concept of 'fingerprints' has been proposed to identify the instabilities that cause transport losses in the pedestal from many candidate instabilities [30]. It is mainly based on two points: (a) the mode frequency in the plasma frame; (b) the ratios between the mode-induced diffusivities in various channels. Theoretical analyses and gyrokinetic simulations of two DIII-D discharges reveal that the QCF observed in the experiment is MTM, and the electron energy transport in the pedestal is dominated by MTM and ETG. In addition, an overview of experiments and simulations on the turbulence dynamics during the pedestal evolution between ELMs are given, and the future research challenges are proposed [31]. Lately, a series of experiments and gyrokinetic simulations have been performed using the concept of 'fingerprints' and find that the magnetic fluctuations observed on DIII-D and JET are identified as MTMs [32][33][34][35][36][37]. Although a large number of experimental and theoretical simulation studies on edge quasi coherent modes have been carried out. However, due to the lack of sufficient experimental evidences and more in-depth and accurate theoretical simulation analyses, the physical properties of these quasi-coherent modes have not been fully determined, and their influence on the edge turbulent transport and the evolution of the pedestal structure are not very clear.
ECM is ubiquitous in the H-mode discharge of EAST. Previous simulations mainly focused on ECM in ELM-free regime [18,28], while ECM in inter-ELM regime has not been studied. In the present work, we study the linear characteristics of ECM in small ELMy regime and its influence on turbulent transport through linear and nonlinear simulations using the global gyrokinetic code GEM based on the EAST experimental parameters. The paper is organized as follows. In section 2, the simulation model and setup are presented. The linear simulation results are given in section 3. In section 4, the nonlinear simulation results are performed. Summary is given in section 5.

Simulation model and setup
GEM is a δf electromagnetic gyrokinetic Particle-in-Cell code [38,39]. GEM has been used to study the turbulent transport [40], and lately it extends to the micro-tearing mode [41] and Alfven eigenmode [42]. GEM started as a flux-tube code with adiabatic electron in the electrostatic limit. The splitweight scheme [43] and control-variate method were then used to extend the code to electromagnetic regime with kinetic electron [38]. Later, extension to general magnetic equilibrium was accomplished [39] and the fluid electron model was also improved and implemented in general equilibrium. Recently, the interface module of GEM has been developed and improved, which can be easily used to read the experimental magnetic equilibrium and profile data, and thus simulate and analyze the experimental results. At present, GEM has many features, such as radially global simulation for arbitrary axisymmetric tokamak equilibria, electromagnetic perturbations, multiple ion species, pitch-angle scattering collisional operator and energetic particles with the isotropic slowingdown distribution, etc,. In this paper, we use the global electromagnetic version with the experimental magnetic equilibrium and profiles to study ECM in EAST.
In order to avoid some numerical instabilities, the parallel canonical momentum p ||α = v || + (q α / m α ) A || is used as the coordinate, the so-called p || -formulation. The gyrokinetic equation used in GEM is as follows: where α = i, e, and For low β tokamak plasmas with β << 1, the guide center drift velocity can also be expressed as where only electron collision is considered and a Lorentzian operator is used.
To reduce the particle noise, the δf method is used in GEM.
where ε pα = µ α B + m α p 2 ||α 2, v tα = T α / m α . For ions, the weight w i = δf i / f 0i evolves according to whereε pi = µ i ⇀ v Gi · ∇B + m i p ||iṗ||i . A split-weight scheme is used for electrons, thus the perturbed distribution of electrons can be written as where ε g is a adjustable free parameter. The adiabatic electron response corresponds to ε g = 1. The weight of electrons evolves according to where w e = h/ f 0e , ∂ϕ/ ∂t is obtained from the vorticity equation, which is the time partial derivative of the quasineutrality equation,ṅ where n p = − qi Ti´( ϕ − ⟨ϕ⟩)f 0i d 3 v is the ion polarization density, and perturbed ion and electron densities are At present, only A || component is considered in GEM, so Ampere's equation can be written as where The last two terms on the left of equation (9) are generated by the equilibrium distribution with the parallel canonical momentum, also known as skin terms. Since m i / m e >> 1, only the electron skin term needs considered. This term needs to be exactly canceled by the corresponding part carried by the particles and their weights, and thus can potentially lead to a severe accuracy problem even at moderate β (βm i / m e ⩾ 1), this is the so-called electromagnetic cancellation problem. To solve this problem, methods such as controlled variable method [38], mixed variable method [44] and pull-back transformation method [45] are proposed. In GEM code, this problem is solved by evaluating the skin terms using the same marker particles and the same scattering operation as that used for the perturbed current term in equation (9).
The coordinates used in GEM are field-line following coordinates and (x, y, z) are defined by are the usual toroidal coordinates, R 0 is the major radius at the magnetic axis, r 0 is the minor radius at the center of the simulation domain, q 0 = q (r 0 ) is the safety factor, q = ⃗ B · ∇ζ ⃗ B · ∇θ. In addition, GEM uses the second-order Runge Kutta method to solve the evolution equation of particle weight, and uses the predict-correction method to promote the guiding center motion of particles, while the Fourier series method is used to solve the field equation. Detailed descriptions of the code can be found in [38,39].
The experimental magnetic equilibrium used in GEM is provided in the eqdsk files calculated by EFIT [46] which contains the poloidal flux ψ on the (R, Z) grid of the right-handed cylindrical coordinates (R, Z, ζ), in which R is the radius, Z is the height, and ζ is the toroidal angle. To connect the experimental data to simulations, we need to find each flux surface, label the flux surface by r, and define a mapping between (R, Z) and (r, θ). We have developed and improved the interface module of GEM, which reads the poloidal flux ψ (R, Z) and poloidal current function f (ψ) from eqdsk files, and then obtains the equilibrium magnetic field B (r, θ), R (r, θ), Z (r, θ) and safety factor q (r) required by the simulation through linear interpolation. The plasmas profiles are also transformed from poloidal magnetic flux ψ to minor radius r.
The EAST discharge 54 804 with obvious ECM signal [47] is used for simulation and some typical discharge parameters are shown in figure 1. The L-H transition occurs around 3.3 s. The H-mode is initially obtained by heating with 2.45 GHz lower hybrid current drive and co-current neutral beam injection, and the ELM type at this time is large type I ELMs.
When the counter-current neutral beam injection (Ctr-NBI) was injected into the plasma around 3.72 s, the large ELMs begin to turn into small ELMs. As shown in figure 1(d), an ECM appears after L-H transition, and its spectrum becomes wider after Ctr-NBI was injected. The equilibrium profiles are calculated by EFIT based on this discharge at time 4.5 s, which is a pedestal recovery period in the ELM cycle, as shown in figure 2. The simulation domain is restricted in the edge region of tokamak around 0.899⩽ r/ a ⩽ 0.999. The magnetic field B 0 = 1.67T, R 0 = 1.92m, a = 0.44m, and the realistic ion electron mass ratio is used. The ratio between the plasma pressure and the magnetic pressure at the center of simulation domain is β e = 2µ 0 n 0e T e B 2 0 ≈ 1.7 × 10 −3 , thus electromagnetic effects have to be considered. The number of grid cells in the radial direction is taken to be n x = 128, in toroidal direction n y = 32 and along the field-line coordinate n z = 64. The time step is ∆tω ci = 3, where ω ci is the ion cyclotron frequency. Convergence is achieved with respect to linear and nonlinear simulations. In addition, for simplicity, we do not consider the X-point here, but consider the up-down symmetric magnetic flux commonly used in previous simulations.

Linear simulation results
In this section, we study the linear characteristics of ECM, such as the mode frequency, dominant toroidal mode number, propagation direction, effect of collision, and driving mechanism. Figure 3 shows the frequency and linear growth rate versus toroidal mode number. Note that only the electron collisional effect is considered, and the effects of uniform collision rate and non-uniform collision rate are compared. The nonuniform collision rate varies with the change of density and  temperature, and its radial distribution is shown in figure 4. The value of uniform collision rate is the collision rate at the center of the simulation domain, i.e. v e / ω ci ≈ 1.07 × 10 −3 , which is also equivalent to the effective collision rate v * e ≈ 1.15, where v * e = 6.92 × 10 −18 qRn e Z ln Λ e T 2 e ε 3/ 2 [48]. The results show that the influence of the non-uniform collision rate on the frequency and growth rate is almost the same as that of the uniform collision rate, indicating that it is sufficient to consider the uniform collision rate. Nevertheless, for accuracy, unless otherwise stated, we consider the non-uniform collision rate in subsequent simulations. Figure 3(b) shows that the dominant toroidal mode number is about n = 18 (k θ ρ i ≈ 0.27), which is consistent with the experimental results that the toroidal mode number of ECM is n = 16-19 [18,47]. As shown in figure 3(a), the frequency of the most unstable mode is about 48 kHz, and propagates along the direction of electron diamagnetic drift (positive frequency in plasma frame), which are also consistent with the experimental results [18,47]. The simulation results indicate that the mode is an electron mode, which may be caused by the kinetic electron effect. We can preliminarily infer that this mode is TEM due to the fact that this mode is a ballooning type electrostatic (see figure 5) electron mode with ion-scale (k θ ρ i ≈ 0.27) from the frequency channel of the concept of 'fingerprints' [30]. Note that ITG can also propagate along the direction of electron diamagnetic drift in the strong gradient regime [49]. Therefore, the propagation direction cannot be used as a decisive criterion to determine the type of one mode, especially for the edge plasma with strong gradient. In order to more accurately judge whether this mode is a TEM, some simulations considering only trapped or passing electrons are carried out. It is noted that the parallel velocity should satisfy the relation v || > 2µ (B max − B (r, θ))/ m for passing electrons. While for the trapped electrons, the parallel velocity should satisfy the relation v || < 2µ (B max − B (r, θ))/ m. It is found that the mode becomes stable (almost unchanged) when only passing electrons (trapped electrons) are retained, indicating that the mode is indeed the TEM. More evidence will be given later to prove that this mode is TEM.
There are two different TEMs [28,50], one is the collisionless trapped electron mode (CTEM), which is mainly driven by electron toroidal precessional resonance without the need of collision. The other is the DTEM, which is mainly driven by temperature gradient and requires large collision. The frequency and growth rate without collision effect are also shown in figure 3. It is found that the mode is stabilized by the collision effect, which further suggests that the mode may be the CTEM. This seems to be different from the previous simulation results that ECM is interpreted as the DTEM [18,28]. It should be noted that previous ECM is observed in ELM-free regime [18], while the ECM studied here is observed in small ELM regime [47]. These two coherent modes may belong to different modes, this requires a detailed comparison of ECM under different parameter conditions and is left for future study. In addition, the pedestal conditions are complex and the relevant turbulence is diverse. The turbulence dominated by different experiments may be different. Moreover, the turbulence in pedestal is very sensitive to parameters [3]. Different plasma profiles, different physical models and even different simulation settings may get different results. This requires very detailed comparative studies, which is beyond the scope of this paper and will be studied in the future. For the most unstable mode n = 18, the ratio of the electron diamagnetic drift frequency to the mode frequency, ω * e / ωis 17, and the ratio between the electron precession frequency and the mode frequency, ω pre,e /ω is approximately 0.44, where ω * e = k θ ρ e v te / L n and ω pre,e = k θ ρ e v te /R 0 . This indicates that the mode is driven by the electron toroidal precessional resonance, which is consistent with the previous theoretical result [50].
Poloidal mode structures of electrostatic potential and parallel magnetic vector potential with toroidal mode number n = 20 are shown in figure 5. Unlike the typical ballooning structure that peaks on the outboard midplane at θ = 0, it peaks near the top and bottom of the plane (θ ≈ ±π/ 2), where θ is the poloidal angle, which is consistent with the previous simulation results [25,26] and may be induced by the large temperature and density gradients in pedestal. In addition, the amplitude of the electrostatic potential is much larger than that of the parallel magnetic vector potential, indicating that the mode is an electrostatic mode, which is also consistent with the experimental results [18,47]. The top of figure 6 shows the mode structure of electrostatic potential with pressure gradient multiplied by coefficients 0.2, 0.6 and nominal experimental values, respectively. The results clearly show that the peak of mode structure moves from θ = 0 to θ ≈ ±π/ 2 or even θ ≈ ±π as the pressure gradient increases, which is consistent with the previous simulation result [25]. The bottom of figure 6 shows that with the increase of the toroidal mode number, the mode structure gradually narrows and moves inward. This is also consistent with the previous theoretical result [51].
Previous theoretical study demonstrates that in the steep density and temperature gradients regime, the growth rate of DTEM is directly proportional to the collision rate for low collision rate v e / εω << 1, while the collision has a stabilizing effect on the DTEM at high collision rate v e / εω >> 1 [52]. To better study the collision effect, the dependence of frequency and growth rate on collision rate with toroidal mode number n = 20 is performed, as shown in figure 7. It is noted that the uniform collision rate is considered and the experimental value is v e / ω ci ≈ 1.07 × 10 −3 . The maximum collision rate is v e / ω ci ≈ 4 × 10 −3 , which is also equivalent to the effective collision rate v * e ≈ 4.30. The results clearly show that collision significantly decreases the growth rate monotonically, even when v e / εω << 1 (see the results shown in the zoom-in plot in figure 7). This is consistent with the previous simulation result [53] and indicates that the mode is indeed the CTEM. In addition, it is found that the effect of collision on the mode frequency is not monotonic. As the collision rate increases, the frequency first increases, then remains unchanged, and finally decreases gradually. The electron collision rate v e and the ratio of plasma pressure to magnetic pressure β e depend on the electron density. In order to simulate a more realistic situation, we further study the properties of the mode for the dependence on  the collision rate and/or β e by changing the electron density while keeping the gradient scale length unchanged, the results are shown in figure 8. Note that the coefficient 1 represents the experimental value. It is found that both β e and collision can stabilize the mode. However, compared with the strong collision stabilization, the electromagnetic stabilization is relatively small, which indicates that the stabilization mainly comes from the collision rather than a combination of collision and electromagnetic effects, and also suggests that this mode is an electrostatic mode. In addition, the mode frequency is mainly affected by electromagnetic effects, which may be induced by the reduction of parallel electric field [26].
It is known that CTEM can be driven by the electron temperature and/or density gradients. To study the influence of   these two driving sources on the mode, the electron temperature and density gradient are multiplied by the coefficients of 0.6, 0.8, 1.0, 1.2 and 1.4, respectively, while other parameters remain unchanged and the results are shown in figure 9. For completeness, the influence of ITG on the mode is also considered. It is found that both electron temperature and density gradients can destabilize the mode, while the ITG stabilize the mode, which is consistent with the theory of CTEM. In addition, it is shown that electron and ITG decreases the mode frequency, while the density gradient increases the frequency first and then has little influence on the frequency. Other parameters, such as the inverse aspect ratio ε, magnetic shearŝ, electron to ion temperature ratio T e / T i and so on also have significant influence on the frequency and growth rate of the CTEM [53]. Since it is not the focus of this paper, we will not study it in detail here.

Nonlinear simulation results
In this section, we study the nonlinear characteristics of ECM, such as the particle and heat fluxes induced by ECM, the effect of zonal flow on turbulence and the saturation mechanism of turbulence. In order to perform nonlinear simulations, we double the radial and toroidal grid points, and increase the simulation time by 3-4 times. The simulation domain is a half of a torus and the toroidal modes are chosen to be even numbers n = −30, −28, .., −2, 0, 2, 4, .., 28, 30. It should be pointed out that the low-n mode of n = 2 is filtered out in simulations due to the high-n approximation assumed in the field solvers of GEM code [38,39].
The time evolution of the particle and heat fluxes of electrons and ions are shown in figure 10. Note that the particle and heat fluxes are normalized by the gyro-Bohm particle flux Γ GB = n e c s (ρ s /a) 2 with c s = T e / m i , ρ s = c s / ω ci and the gyro-Bohm heat flux Q GB = n e c s T e (ρ s /a) 2 , respectively. It is clearly shown that the saturated particle and heat fluxes caused by ECM are mainly due to the perturbed electrostatic ExB drift, which further confirms that ECM is an electrostatic mode. It is also found that the particle and heat fluxes of electrons and ions are almost equal. Therefore, only the particle and heat fluxes of electrons are analyzed later. In addition, it is shown that ECM can drive significant outward particle and heat fluxes, which exhausts the accumulated energy, fuel particles and impurities, and thus greatly facilitating long pulse H-mode sustainment [18]. It should be noted that the particle and heat fluxes of electrons and ions caused by magnetic fluctuations are inward, except that the electron heat flux is outward. However, compared with the fluxes caused by ExB fluctuations, these fluxes caused by magnetic fluctuations can be ignored. Therefore, in the following manuscript, we only analyze the fluxes induced by ExB fluctuations. The time averaged saturated particle and heat fluxes are 1.51 × 10 20 s −1 m −2 and 22.7 kWm −2 , respectively. Considering that the whole plasma surface area is S ≈ 40m 2 , the estimated heat flux QS ≈ 0.91 MW accounts for about 35% of the net input heating power Q input ≈ 2.58 MW, which is a considerable heat transport. Since the fluxes induced by ECM are mainly due to the perturbed electrostatic ExB drift, and the fluxes of electrons and ions are almost equal, we can conclude that ECM is TEM according to the transport channel of the concept of 'fingerprints' [30].
In order to study the contribution of modes with different toroidal mode number to the fluxes, the toroidal Fourier decomposition of electron heat flux is performed and the time evolution of electron heat flux for each toroidal mode number is shown in figure 11(a). To make the simulation results look more clearer, the higher n modes with smaller magnitudes and less contribution to flux are not shown. It is found that the mode of n = 14 carries the most flux in the early non-linear saturation phase, while the mode of n = 10 dominates the contribution to the flux after t = 5 × 10 4 . In addition, the modes of n = 6, 16 also contribute a lot to the flux. The toroidal Fourier decomposition of other fluxes has similar results. It is noted that the electron heat flux is induced by the perturbed electrostatic ExB drift. The time evolution of electrostatic potential for each toroidal mode number is shown in figure 11(b). Note that the electrostatic potential and parallel magnetic vector potential in figures 11 and 12 are normalized by T u / e and T u / ev u , respectively. It is found that the intermediate-n modes of n = 14, 18 grow fastest in the linear phase, while in the nonlinear saturation phase, the low-n modes like n = 4, 6 dominate. This is not entirely consistent with the flux results shown in figure 11(a) and indicates that the low-n modes such as n = 4, 6 do not dominate the contribution to the flux. In addition, it is found that zonal flow is generated via three wave coupling interaction with a growth rate of twice the linear growth rate of the mode. However, the saturation amplitude of the zonal flow is smaller than that of the dominant modes, indicating that zonal flow may have little effect on the saturation of the turbulence, which is different from the saturation of the ITG turbulence that zonal flow dominate in the non-linear saturation phase. The time averaged saturated electron particle and heat fluxes as well as the potentials versus the toroidal mode number are shown in figure 12. It should be pointed out that selecting the time window is important and can cause large error bars for single modes, which can be seen from figure 11. However, for the total flux it is sure that the selecting the time Figure 11. The time evolution of electrostatic electron heat flux (a) and potential (b) for each toroidal mode number. Note that the higher n modes with smaller magnitudes and less contribution to flux are not shown. window is not that crucial as the decay of one mode is compensated by the growth of other modes. Note that the time window for averaging in figure 12 and subsequent figures is 3.9 × 10 4 ⩽ t ⩽ 6.9 × 10 4 . The results clearly show that the saturation fluxes are mainly contributed by the mode of n = 10. The modes of n = 6, 14, 16 also contribute a lot to the fluxes and the contributions of other modes are relatively small, especially the high-n modes of n > 20. In addition, the potentials are dominated by the modes of n = 4, 6, and the zonal flow and zonal current are smaller compared with the dominant modes. It should be pointed out that the turbulence saturation level of perturbed electrostatic potential energy is about 1% of thermal energy of equilibrium, which indicates the δf method is still valid in this work.
Zonal flows can be driven by the drift wave turbulence through nonlinear Reynolds stress and in turn regulate turbulent transport via flow shearing, and thus reduce the transport level [54]. It is well-known that zonal flows are important in regulating the ITG turbulence [55]. However, the effect of zonal flows on CTEM turbulence is contradictory and sensitive to parameters, such as the electron temperature gradient, electron to ion temperature ratio, and magnetic shear [53,56,57].
To study the influence of zonal flow on turbulent transport, we compared the electron heat flux with and without zonal flow, the results are shown in figure 13. It is found that zonal flow has little effect on the particle flux, but can reduce the heat flux by about 14%, which indicates that zonal flow is not the dominant saturation mechanism of turbulence. This is consistent with the previous simulation results [56]. It is known that the shearing rate should be much larger than the maximum linear growth rate for zonal flow to be the dominant nonlinear saturation mechanism [58]. The shearing rate of the zonal flow can be readily calculated as ω s = d 2 ⟨φ⟩ dx 2 . The maximum shearing rate is ω s / ω ci ≈ 2 × 10 −4 , which is smaller than the maximum linear growth rate γ/ ω ci ≈ 2.91 × 10 −4 , thus the zonal flow is not the dominant saturation mechanism of turbulence for the parameters studied here. Recently, simulation studies have shown that zonal flow has little effect on turbulent saturation in strong gradient regime due to the small eddy size [59], this seems consistent with our present results.
According to the linear simulation results in section 3, the collision has a strong stabilizing effect on the mode. It is necessary to study the nonlinear stabilization effect of collision.
To study the influence of collision on turbulent transport, we compared the saturated electron particle and heat fluxes with and without collision, the results are also shown in figure 13. It is found that collision has little effect on the particle and heat fluxes. This is not very consistent with the linear simulation results. As shown in figure 14, we further compared the saturated electron heat flux versus the toroidal mode number Figure 13. The time averaged saturated electron particle and heat fluxes for three cases. The nominal case represents that both zonal flow and collision are considered in simulation, the w/o coll case represents that collision is not considered in simulation, while the w/o zf case represents that zonal flow is not considered in simulation. between the case with and without collision and found that collision significantly increases the flux contributed by the intermediate-n modes of n = 10,14,16, but greatly reduces the flux contributed by the low-n modes of n = 4,6,8, resulting in relatively small change in the total flux. The effects of ETG and density gradient on turbulent transport are also studied and shown in figure 15. It is found that density gradient has significant influence on the particle and heat fluxes. However, the ETG has relatively small influence on the particle and heat fluxes. This is consistent with the linear simulation results of figure 9, which indicates that density gradient is more important for ECM and we can control the particle and heat fluxes by adjusting the density profile.
Since zonal flow is not the dominant saturation mechanism of ECM. It is thus desirable to find out other possible saturation mechanisms. The nonlinear spectral cascade is universal in turbulence spectral and can affect the nonlinear saturation level and related turbulent transport. In order to study the nonlinear spectral cascade, the toroidal Fourier decomposition of turbulence intensity and electron heat flux are performed, the Figure 15. The time averaged saturated electron particle and heat fluxes for three cases. The nominal case represents that realistic experimental parameters are considered in simulation. The other two cases represent that only 50% of the density gradient value and 50% of the electron temperature gradient value are considered in the simulation, respectively. results are shown in figure 16. To see more clearly, the turbulence intensity and electron heat flux are normalized by its own maximum value at each moment. A clear spectral cascading in the toroidal direction is observed. It is shown that in the linear phase t = 0 ∼ 2 × 10 4 , the turbulence intensity and electron heat flux are dominated by the modes of n = 14,18, which are almost the most unstable mode in the linear simulations. Then, at the early nonlinear saturation phase t = 2 ∼ 4 × 10 4 , the turbulence intensity and electron heat flux are dominated by the modes of n = 14. Finally, at the end of nonlinear saturation phase t = 4 ∼ 6.9 × 10 4 , the turbulence intensity is dominated by the low-n modes such as n = 4, 6, while the electron heat flux is mainly contributed by the mode of n = 10. This indicates that the turbulence energy transfers from short wavelength (high-n) modes to the long wavelength (low-n) modes. This inverse spectral cascading process is considerably universal in the simulations and independent of parameters, which we believe this is more universal saturation mechanism for ECM.
The edge radial electric field is usually large in the pedestal of H-mode plasmas and may be very important for ECM and relevant turbulent transport. However, there is no experimental data of radial electric field, to study the influence of radial electric field on turbulence, we construct an analytical expression of radial electric field, i.e. E r = −E 0 exp −((r − r c )/ ∆r) 2 , where r c = 0.97a, ∆r = 0.015. Note that the maximum value of the edge radial electric field in the H-mode discharge of EAST is more than ten kV m −1 . Here we take E 0 = 9.28 kV m −1 , the profile of radial electric field is shown in figure 17. The nonlinear simulation with the radial electric field is performed and is compared with the previous result without the radial electric field, as shown in figure 18. Note that the maximum shearing rate of radial electric field is ω ExB / ω ci ≈ 2.8 × 10 −3 , which is much larger than the maximum linear growth rate γ/ ω ci ≈ 2.91 × 10 −4 , indicating that radial electric field may play an important role in suppressing   turbulent transport. From figure 18, it is found that the radial electric field has a strong stabilization effect in the linear phase so that the peak value of the turbulence intensity and heat flux is much lower with radial electric field, which is consistent with the linear simulation results that the radial electric field has a significant effect on the linear growth rate of the ECM. In addition, the saturation level of heat flux is reduced by about 54% with radial electric field, which indicates that the radial electric field can greatly reduce the turbulence intensity and transport level. These results confirm the important role of the radial electric field in the edge turbulent transport. It should be pointed out that the radial electric field does not affect the characteristics of ECM and the main conclusions of this paper, and we will further study it with the experimental radial electric field data in the future.

Summary
In this paper, we present global linear and nonlinear simulation studies of ECM based on the EAST experimental parameters using the gyrokinetic code GEM. Some new and interesting characteristics of the ECM are observed. The linear simulation results show that ECM is an electrostatic mode and peaks in n = 18 with the frequency about 48 KHz, and propagates along the direction of electron diamagnetic drift, which are consistent with the experimental results. It is also shown that density and electron temperature gradients destabilize the mode, while collision stabilize the mode. The nonlinear simulation results show that the saturated particle and heat fluxes induced by ECM are mainly due to the perturbed electrostatic ExB drift, and the fluxes of electrons and ions are almost equal. In addition, it is shown that ECM can drive significant outward particle and heat fluxes, which exhausts the accumulated energy, fuel particles and impurities, and thus greatly promoting the maintenance of the long pulse H mode. The Fourier decomposition of fluxes and potentials show that the modes of n = 14, 18 grow fastest in the linear phase, while the lown modes like n = 4, 6 dominate in the non-linear saturation phase and the saturation fluxes are mainly contributed by the mode of n = 10. It is found that zonal flow is not the dominant saturation mechanism of turbulence for the parameters studied here. The nonlinear spectral cascade analyses imply that the inverse spectral cascading is more universal saturation mechanism for the turbulence since it is fairly universal in the simulations and independent of parameters. It is also found that the radial electric field can greatly reduce the turbulence intensity and transport level. From the analyses of frequency and transport channels, it can be concluded that ECM appears to be the CTEM. It is worth pointing out that various coherent modes co-exist in the pedestal, the unambiguous identification and saturation mechanisms of these modes, and the role played by these modes in ELM bursts and formation of the pedestal through turbulent transport are still far from being resolved [31]. This requires more elaborate experimental and simulation studies, and the present work is only a small step forward.