Electromagnetic drift wave instability in tokamak plasmas with strong pedestal gradient

The linear eigenmode characterizations and the nonlinear turbulence energy spreading of the drift waves in a tokamak plasma with strong pedestal gradient are numerically investigated based on an electromagnetic Landau fluid model. By the linear eigenmode analysis, it is found that the dominant instability in the low β regime is the ion-temperature-gradient (ITGc) mode and the electron drift wave instability (eDWI p ) in the core and edge region with strong density gradient, respectively. Multiple eigenstates of the eDWI p with different peak locations in the poloidal direction can be obtained by the eigenvalue problem solver. The dominant one is the high order eDWI p corresponding to the unconventional ballooning mode structure with multiple peaks in the poloidal position, in contrast to the conventional modes that peak at the outboard mid-plane, and has been verified through initial value simulation. In the high β regime, the dominant eigenmodes in the core and edge region are the conventional and unconventional kinetic ballooning modes respectively. In the nonlinear simulation, an inward turbulence spreading phenomenon during the quasi-saturation phase of the edge turbulence is clearly observed. The inward speed of the turbulence energy front in the high β regime is much faster than that in the low β regime. It is interestingly found that the speed of the turbulence energy front increases with the increase of the plasma β in the low β regime, while it is almost unchanged in the high β regime. It is identified that the turbulence spreading in the low and high β regimes are determined by the nonlinear dynamics and the linear toroidal coupling respectively.


Introduction
Drift wave instabilities, a class of low frequency eigenmodes in the low β plasmas, are important for the understanding of plasma turbulence and the corresponding cross-magnetic field transport of particles, energy and momentum in magnetic fusion devices [1,2]. Here, β = p0 B 2 0 /2µ0 is the ratio of plasma pressure to the magnetic pressure. Therefore, the properties of these drift wave instabilities have attracted much attention over the past several decades.
Based on the previous studies [1,3], the ion-temperaturegradient (ITG) modes, typical drift wave eigenmodes driven by the temperature gradient of the thermal ion species, have been considered as one of the main candidates for the anomalous transport of ion heat in tokamak plasmas. It is well known that both the slab branch and the toroidal branch of ITG modes are kinds of finite threshold instabilities [1,2], i.e. there is a stability threshold for the R/L Ti or η i = L n /L Ti beyond which the ITG or η i modes become unstable, where R is the major radius of tokamak, L Ti and L n are the scale lengths of the temperature gradient and density gradient respectively. It is important to find that the electromagnetic effects in a finite β plasma plays a stabilizing role in the drift wave instabilities [4][5][6], therefore the relatively high β plasma is favorable for the nuclear fusion reaction. However, it is easy to excite the kinetic ballooning modes (KBMs) in toroidal geometry with high β regime [7][8][9][10]. Compared to the understanding of the electrostatic ITG turbulence, the KBM turbulence in the high β plasmas is still rare until now and need to be further studied.
In recent studies, it is found that the properties of the drift wave eigenmodes in the pedestal region of the high confinement plasmas are much different from that in the center region with relatively lower gradient [11][12][13]. For instance, it is found that the most unstable mode in the steep gradient region is not the ground state (the eigenstate that holds the structure of the low degree Hermite polynomial) of the electrostatic drift wave eigenmode [11] or electromagnetic KBM mode [13]. Even though the growth rate of the dominant mode with non-ground eigenstates in the steep gradient regime is larger than that of the dominant mode with ground eigenstate, the unconventional mode structures of the non-ground eigenstate can lead to the significant reduction of the trend of transport level on the background gradient due to the decrease of the effective correlation length. This may provide a possible mechanism for the low to high confinement mode transition without the need of shear flows effect [12]. Moreover, the unconventional ballooning structures for the electron temperature gradient instabilities in the steep gradient region of a JET pedestal are also observed in the gyrokinetic simulations [14].
Another interesting topic of the drift wave turbulence transport in tokamak plasmas is the turbulence spreading between the core and edge region. In the low to high confinement mode transition, the formation of the pedestal with steep gradient in the edge region leads to the suppression of the turbulence transport with the establishment of the edge transport barrier. However, during the back-transition or the high to low confinement mode transition, the fluctuations of turbulence can reach a high level accompanied by the disappearance of the transport barrier [15][16][17]. Thus, with the increased turbulence in the edge region, the turbulence spreading could play an important role in the core-edge coupling.
It is noted that the process of the turbulence spreading is beyond the scope of the local mixing theory. In the paradigm of local mixing, the level of the turbulence transport according to the unstable eigenmodes driven by the local gradient is assessed through the mixing length criteria at a local region related to the mode width. Clearly, for the turbulence spreading in the core-edge coupling system, it is necessary to develop a nonlocal theory including the meso-scale dynamics. To understand the nonlocal phenomena in the tokamak plasmas, Garbet et al [18] first studied the process of turbulence spreading by considering the toroidal coupling and nonlinear coupling effect. It is found that the linear toroidal coupling effect through a convection of the elongated vortex played an important role in the process of turbulence spreading in tokamak plasmas. Chen et al [19] demonstrated that the zonal flow driven by the pump of drift wave can facilitate the turbulence spreading due to the enhanced dispersiveness. Later, Gurcan et al [15] pointed out that the toroidal coupling and the zonalflow induced side-band coupling are not the necessary condition for the ballistic turbulence spreading. Such propagation of the turbulence energy front, proportional to the geometric mean of the growth and diffusion rates of turbulence, is possible via the nonlinear dynamics alone [15,20]. This result is further verified by a gyro-kinetic particle simulations of toroidal ion temperature gradient turbulence spreading from the edge to the core regions of tokamak [16]. Fluid simulation of turbulence spreading in reversed magnetic shear configuration also supported this result [17]. However, until now, there is still few investigation on the topic of electromagnetic effect or finite β effect on the turbulence spreading in the core-edge coupling system.
In this work, the basic features of the electromagnetic drift wave eigenmodes in the edge region with strong pedestal gradient are studied by using a global electromagnetic Landaufluid model [8]. The inward propagations of the electromagnetic turbulence fronts in different β regimes are analyzed in detail. The remainder of this paper is organized as follows. In section 2, the global electromagnetic Landau-fluid model used in this work is briefly introduced. The characterizations of the unstable eigenmode and the inward propagations of the electromagnetic turbulence fronts are analyzed in section 3. Finally, a short summary and discussion are given in section 4.

Modeling equations
A global Landau-fluid model [8,[21][22][23][24][25] is employed in this work, which consists of five dynamic equations for perturbed densityñ e , vorticity ∇ 2 ⊥φ , parallel ion velocityṽ ∥i , parallel magnetic vector potentialÃ ∥ , and ion temperatureT i , Here, the variables with subscript 0 means the time independent equilibrium variables. p i = p i0 a + n 0Ti + T i0ñi and p e = p e0 a + T e0ñe are the ion and electron pressure. τ = T e0 /T i0 is the ratio of electron and ion equilibrium temperatures. According to the Ampère's law, the parallel current can be written asj ∥ = −∇ 2 ⊥Ã ∥ , j ∥0 = ∇ × B 0 /β and j ∥ = j ∥0 a +j ∥ . From the quasi-neutral condition, the perturbed ion and electron densities have the relationñ e ≃ñ i . The circular tokamak geometry (r, θ, ζ) is employed in the numerical simulation, where r is the radius of magnetic surface, θ and ζ are poloidal and toroidal angles, respectively. It follows that the Lagrange Poisson bracket [ f, g] = (∂ r f∂ θ g − ∂ r g∂ θ f )/r, where ϵ = a/R 0 is the inverse aspect ratio, a and R 0 the minor and major radius respectively. The derivatives related to the diamagnetic drift are expressed as The derivatives related to the magnetic drift is expressed as iω Di The variables in the equations are normalized as follows The values at magnetic axis, such as n c , T c , are selected as the references for normalization. v ti = √ T c /m i and ρ i = v ti /ω ci with ω ci the ion cyclotron frequency are the thermal ion velocity and the ion gyro-radius at the axis. Parameter β is defined as a half of the ratio of the thermal ion pressure to magnetic pressure in the center r = 0, i.e. β = n c T c /(B 2 0 /µ 0 ). This model has been employed in the global simulations to study the characteristics of turbulence transport, such as the global structure of zonal flow [8], multi-scale interaction between turbulence and tearing mode [26], the symmetry breaking induced by the electromagnetic fluctuations [27], etc.
Each perturbed variable in equations (1)-(5) can be decomposed as with m and n being the poloidal and toroidal mode numbers. Then we can obtain the linearized equations (1)-(5) by neglecting the nonlinear coupling terms (L k1+k2 ∼ f 1,k1 f * 2,k2 ) related to the Poisson brackets. After writing the time dependence of a perturbed field variable f m,n (r, t) =f m,n (r)exp(Ωt), the eigenvalue problem of the linearized equations can be obtained as where The matrices M m , M d , M p correspond to the operators, respectively, the magnetic drift terms, the diffusivity terms, and the remnant terms in equations (1)- (5). For a fixed toroidal mode number, the eigenvalue Ω and eigenvector [ñ eφṽ∥,iÃ∥Ti ] T of equation (6) can be obtained using an eigenvalue problem solver after giving the boundary conditions. The real and imaginary parts of Ω = ω + iγ are the mode frequency and the growth rate of the eigenmode. The eigenvector can be decomposed into the mode structures of each variable for the unstable mode. This procedure is similar to the previous works in the slab, cylindrical and other geometries [28][29][30][31][32][33][34]. It is noted that due to the contribution of magnetic drift in M m , the eigenmodes with the same toroidal mode and different poloidal modes are coupled together. The matrices in equation (6) are much larger than those in the slab or cylindrical geometries without toroidal coupling effect. Compared to the initial value problem simulation, the strength of the eigenvalue problem simulation is that all the eigenmodes with different growth rates can be obtained, while only the most unstable mode with the fastest growth rate can be obtained in the initial value problem simulation.
In the nonlinear simulation, a two-step predictor-corrector method is used in time advance. The finite difference method in the radial direction and Fourier spectral method in poloidal/toroidal directions are employed for space discretization [8,10]. The radial grid number is set to be 600 to ensure the convergence with strong pedestal gradient. The maximum value of the toroidal mode number is set to be n max = 30 and m is selected to include the resonant perturbations distributing within q(ρ = 0.15) ⩽ m/n ⩽ q(ρ = 1.0), where ρ = r/a. The zero Dirichlet condition is imposed for all f m,n perturbations on the outer boundary. At ρ = 0, ∂ r f m=0,n=0 = 0. For the modes with m ̸ = 0 or n ̸ = 0, f m,n (ρ = 0) = 0. The code has been benchmarked by comparing to other simulation results with similar model [8,10]. To simulate the turbulence in the strong pedestal region, the main radial domain covers the whole radius up to the last closed flux surface, ρ ∈ [0, 1.0], and it is bounded by buffer zones ρ ∈ (1.0, 1.15] with an artificially increased hyper-diffusivity, which represents the region with a pressure sink [35,36]. Other parameters used in the calculations are ϵ = 0.25, ρ i /a = 0.0125, τ = 1. The normalized hyper-diffusivity and resistivity in the main radial domain are set to be D n = D Ω = D v = D T = 5 × 10 −7 m 4 and η = 4 × 10 −5 .

Numerical results
Based on the global electromagnetic Landau-fluid model given in section 2, dependence of the linear growth rates and the mode frequencies of the unstable eigenmodes on the plasma β are calculated by using an eigenvalue solver. Then, the linear eigenmode structures, nonlinear turbulence features and the inward turbulence spreading are analyzed. The profiles used in this work are plotted in figure 1. Here, the strong gradients of density and temperature near the edge region are employed to model the pedestal of the high confinement mode in tokamak plasmas. The density gradient R 0 /L n (r) and temperature gradient R 0 /L T (r) in the edge region are one order of magnitude higher than those in the core region, while the ratio of the scale lengths of density gradient to that of temperature gradient η i (r) = L n (r)/L T (r) = |∂r ln T i0 (r)| |∂r ln n0(r)| is near unity in the core but less than one in the edge region. Here, L n (r) and L T (r) are the profiles of the scale length of density gradient and temperature gradient respectively.

Linear results
The linear growth rates and mode frequencies of each unstable eigenmode for the linearized equations can be obtained by an eigenvalue solver given in section 2. The spectra of the dominant eigenstates with maximum growth rates for each toroidal mode number and different plasma β are displayed in figure 2. Since the electromagnetic effect has a stabilizing impact on the drift wave type instability [8], it is clearly shown in figure 2 that the linear growth rate of the dominant instability decreases with the increasing of the plasma β in the low β region. However, the KBM can be excited in the high β region. Thus, after the plasma β exceeds a critical value, it is found that the linear growth rate increases with the increasing of β.
It is noted that the different eigenstates of eigenmodes with the same toroidal number n in the core and edge regions can co-exist in the global system. Therefore, the curves of the unstable spectra for the dominant instability seem to be discontinuous in figure 2. Roughly speaking, the dominant toroidal mode number for the instabilities in the edge region is relatively smaller than that in the core region, because the eigenmode with relatively large poloidal mode number m = nq(r) in the edge region is more easily damped due to the viscous effect than that in the core region with the same n. In this work, the viscous effect is included in the model as the D U ∇ 4 ⊥ ϕ terms. Since the operator ∇ 4 ⊥ is proportional to the poloidal wave 4 , the viscous effect increases with the increase of the poloidal mode number m. If the perturbation holds large m, the mode is more easily damped by the viscous effect. So the peaks of the curves in the low n and high n regions are for the unstable eigenmode in the edge and core regions respectively.
Moreover, it is found that the mode frequencies of the dominant eigenstates in the low β and edge region is positive (red region in figure 2(c)), i.e. the frequencies are in the electron diamagnetic drift direction. However, the mode frequencies of the dominant eigenstates in the low β and core region are negative or in the ion diamagnetic drift direction. It suggests that the eigenmodes in the low β region propagating in the ion diamagnetic drift direction are the ITG modes, as discussed in the previous [1], which can be stabilized by considering the finite β effect. In the edge region with strong density gradient, the ratio of η i (r) holds a small value, as shown in figure 1(b). Thus, the dominant instability in the edge region is clearly not the ITG modes. In fact, these unstable eigenmodes propagating in the electron diamagnetic drift direction can be classified as the electron drift wave instability (eDWI) in the strong gradient limit. The numerical results that the linear growth rates of this type eigemode are more sensitive to the density gradient in the edge region can support this point.
Furthermore, the eigenvalue solver can obtain all the unstable eigenstates of the system with fixed toroidal number n. To understand the multiple branches of the global eigenmode, figure 3 gives all the eigenvalues of the unstable eigenstates for toroidal mode number n = 5 in the low β = 0.15% and high β = 0.85% regions. The superscript characters c and p in figure 3 represent that the mode structures are located in the core region and edge region, respectively. For the low β case, the dominant instability is eDWI p propagating in the electron diamagnetic drift direction. The dominant, second and third subdominant eigenmode in the low β region are indicated with eDWI p B , eDWI p A and eDWI p C in figure 3(a). The eigen  structures of these modes propagating in the electron diamagnetic drift direction are plotted in figures 4(a)-(c). Clearly, the mode structure of the dominant eigenstate showing with multiple peaks near θ p ≈ ±π/2 is different from the conventional one with single peak located near θ p ≈ 0, where θ p = 0 is at the outboard mid-plane. In fact, the mode structure of the dominant KBM instability in the pedestal region indicated with KBM p with relatively small growth rate also presents unconventional structure with peak near θ p ≈ 3 4 π, as shown in figure 4(d). These unconventional structures are similar to the results in [13]. In the core region, the dominant instability is the conventional ITG mode with finite tilting poloidal angle [37], which is indicated by ITG c in figure 3(a). Both of the growth rates of eDWI p and ITG c decrease with the increasing of plasma β, which can be obtained by comparing the results in figure 3(a) for β = 0.15% and figure 3(b) for β = 0.85%. However, the growth rate of KBM increases with the increasing of plasma β. Mode structures ofφ n=5 for ITG c and KBM c are given in figures 4(e) and (f ). Both of them are similar to the electromagnetic simulation results without edge pedestal region [8].

Nonlinear results
To understand the nonlinear behaviors of the electromagnetic instabilities for the profile with strong pedestal gradient, simulations at the initial value are carried out in the low and high β regimes respectively. Figures 5(a) and (b) give the temporal evolutions of the kinematic energies of the perturbations for low β = 0.15% and high β = 0.85%. Since the perturbations with the same toroidal number n and different poloidal numbers m are coupled together in the toroidal geometry, the kinematic energies (E kin n = ∑ m´| ∇ ⊥ ϕ| 2 ρdρ/2) of the perturbations shown in figure 5 are the summations of the perturbed energies for all the poloidal numbers m with the same toroidal number n. In the linear phase, the growth rates of the dominant unstable modes with different toroidal mode numbers are almost the same with those obtained by the eigenvalue problem solver. Since the growth rate of the dominant unstable mode in the low β = 0.15% regime is smaller than that in the high β = 0.85% regime, the temporal evolution from the linear phase to the quasi-saturation phase in the case of β = 0.15% is slower than that in the case of β = 0.85%. The fluctuations in both cases are located near the edge region with strong gradient. The mode structures are both ballooned, and only the dominant poloidal mode numbers are different for the low and   high β regimes, as shown in figures 6(a) and (e). According to the linear eigenmode analysis in section 3.1, the dominant instabilities are eDWI and KBM propagated in the electron and ion diamagnetic drift directions, respectively, for the low and high β regimes. The unconventional mode structures are also similar to the previous results in the strong gradient regimes [12,13].
Due to the zonal components with m/n = 0/0 excited through the nonlinear coupling effect, the fluctuations are gradually saturated near the edge region. In fact, during the early saturation phase with relatively high level of zonal flow, the inward and outward turbulence propagation can be observed due to the nonlinear coupling effect of the fluctuations in the presence of the zonal flow, as shown in figures 6(b) and (f ), similar to the results in [19]. However, after the turbulence spreading into the center region, the toroidal coupling effect becomes important. The zonal flow modulated ballooning structures can be clearly observed in the center region due to the toroidal ITG and KBM instability in the low and high β regimes, as shown in figures 6(c) and (g). This inward turbulence spreading is somehow consistent with gyrokinetic analyses of the turbulence diffusion in the vicinity of limiters in circular plasmas [38]. The zonal flow structures also have the same inward spreading speed as the turbulence fluctuations, as shown in the shadow regions in figures 5(c) and (d). It is interestingly found that the structure of the zonal flow velocity in the case with β = 0.15% is similar to the E × B staircase reported in [39], as shown in figure 5(c).
The inward spreading speed in the low β regime is found to be slower than that in the high β regime. One possible reason is that the growth rate of KBM c in the high β = 0.85% regime is larger than that of ITG c in the low β = 0.15% regime. Another reason may be that the diffusion rate of turbulence which is proportional to the turbulence intensity in the high β = 0.85% regime is larger than that in the low β = 0.15% case. Both of them seem to be reasonable according to the previous work in the electrostatic turbulence in [15,16] and the electromagnetic turbulence simulation in this work. However, it is found that the well organized steamer structure or ballooning structure can be clearly observed in the 'unfavorable curvature' region during the turbulence inward spreading in the high β regime, as shown in figure 6(g). Thus, the toroidal coupling effect may also play an important role through the elongated structures in the turbulence spreading, especially in the high β regime. Moreover, the turbulence propagations in the poloidal and toroidal directions in the quasi-saturation phase are also different in the low and high β regime. In the low β regime, the turbulence propagates in the electron/ion diamagnetic drift direction with dominant instability of eDWI p and ITG c near the edge and core region. However, in the high β regime, the turbulence propagates in the ion/electron diamagnetic drift direction near the edge and core region, even though both of the dominant instabilities are KBM p and KBM c in the whole region. In fact, in the high β cases, the electromagnetic effect becomes important. The long wavelength drift Alfven type perturbations around the low q resonant surfaces become much higher in the high β case than that in the low β case. The evolution of perturbation of the parallel magnetic vector potentialÃ ∥ (a) is plotted in figures 7(a)-(e). The m/n = 3/1 magnetic perturbation coupling with the m/n = 4/1 magnetic perturbation can be clearly observed in the case with β = 0.85%, as shown in figure 7(a). The reason that the turbulence propagates in the electron diamagnetic drift direction is mainly the excitation of long wavelength Alfven type mode in the high β regime with strong electromagnetic effect [10,26,40], as shown in figure 7. Here, the clockwise and counterclockwise directions  Moreover, the dynamics of the turbulence and zonal flow system can be modulated by the long wavelength electromagnetic perturbations. The change of the magnetic topology near the low q resonant surfaces and the stochastic magnetic field induced by the magnetic perturbation with different helicities can cause the variation of the zonal flow structure induced by the turbulence, as shown in figures 5(c) and (d). Furthermore, the corrugation of the zonal current or zonal field profile can be clearly seen near the q = 3 rational surface due to the generation of the long wavelength electromagnetic perturbation with m/n = 3/1, as shown in the blue line in figure 5(d).
To further investigate the characterizations of the turbulence spreading in the low β and high β regime, the inward spreading speed in different β regimes are compared in detail. As shown in figures 8(a) and (c), in the low β case with β = 0.1% and β = 0.2%, the flux propagates inwards after the quasi-saturation of the edge turbulence. The speed of the inward turbulence spreading increases with the increase of the plasma β in the low β regime. This result is similar to the theoretical model proposed in [15], which indicates that the spreading speed of the turbulence front is proportional to the geometric mean of the growth and diffusion rates of turbulence. While, the result in the low β regime by scanning the β value is a little different from that in the [16] by scanning the temperature gradient. The main reason that the spreading speed of the turbulence front increases with the increasing of the temperature gradient is the increase of the instabilities growth rate. However, in the low β regime, the growth rate usually decreases with the increase of plasma β. The increase of the spreading speed of the turbulence front by increasing plasma β is mainly because of the increase of the diffusion rates of the turbulence proportional to the turbulence amplitude, as shown in figures 8(a) and (c). It is noted that, in this work, the source and sink terms are set to recover the initial profile in a diffusing time scale which is much longer than the turbulence spreading time. Thus, the results in this work may not be suitable to explain the long term event in the tokamak plasmas with strong power injection system.
Moreover, in the high β regime, the spreading speed of the turbulence front seems to be not sensitive to the plasma β, even though both the growth rate and the amplitude of the turbulence increase with the increase of β, as shown in figures 8(b) and (d). This result suggests that the nonlinear diffusion mechanism may not be the dominant reason for the turbulence spreading phenomenon in the high β regime. Based on the perturbation structures during the turbulence inward spreading shown in figures 6(g) and (h), it can be found that the linear toroidal coupling effect acting through the streamer structure [18] plays an important role in the process of the high-β KBM turbulence spreading.

Conclusion and discussion
In this work, both the linear and nonlinear characterization of the electromagnetic drift wave instabilities in a tokamak plasma with strong pedestal gradient are numerically studied based on an electromagnetic Landau fluid model by scanning the plasma β in detail. The main results of this investigation can be summarized as follows.
(1) Under the profiles with strong density and temperature gradients, based on the linear eigenmode analysis, the eDWI p propagating in the electron diamagnetic drift direction is observed in the low β regime. Many eigenstates of the eDWI p with different growth rates and mode frequencies can coexist. The dominant one is the high order state corresponding to the unconventional ballooning mode structure, which is not localized at region of unfavorable curvature or with multiple peaks in poloidal positions. The dominant instability in the core region is the ITG c mode propagating in the ion diamagnetic drift direction, which corresponds to the global toroidal ITG branch with finite tilting poloidal angle [37] in the electrostatic limit. In the high β regime, the dominant eigenmodes in the core and edge regions are both KBM propagating in the ion diamagnetic drift direction. The most unstable KBM in the edge region with strong pedestal gradient is also the unconventional one, while the dominant KBM in the core region is the conventional one. It is worth noting that the nature of the dominant instability strongly depends on the input profiles. The results may be different for the particular experimental scenario, if the experimental profiles are different from the given profiles in this work. (2) In the nonlinear simulation, the inward turbulence spreading from the edge region to the core region can be clearly observed when the edge turbulence enters into the quasisaturation phase. It is found that the inward spreading speed of the turbulence front in the high β regime is much faster than that in the low β regime. The mechanisms of the turbulence inward spreading processes are different in the low and high β regimes. In the low β regime, it is found that the inward spreading speed of the turbulence energy front is proportional to the plasma β. Even though the growth rate of the turbulence decreases with the increase of plasma β, the instantaneous intensity of turbulence or nonlinear diffusion rate during the quasi-saturation stage increases with the increase of plasma β in the electromagnetic global simulation with strong pedestal gradient. The turbulent energy front speed is roughly proportional to the turbulence intensity, which can be roughly related to the nonlinear dynamic induced turbulence spreading. However, in the high β regime, the inward spreading speed of the turbulence front is almost unchanged in a quite wide β range, even though both the growth rate and the instantaneous intensity of turbulence increase with the increase of plasma β. The fast turbulent energy front speed can be related to the linear toroidal coupling induced turbulence spreading. The large scale streamer or ballooning structure induced by the long wavelength electromagnetic modes can be clearly observed during the fast inward turbulence spreading, which can verify the crucial effect of the linear toroidal coupling on the turbulence spreading process in the high β regime.
It should be noted that the heat source [41] which is the necessary condition for the maintenance of the edge strong pressure gradient has not been considered in this work. Thus, the nonlinear results discussed above may be not suitable to explain the long term evolution event with strong power injection. Moreover, the penetration of the boundary error field or resonant magnetic perturbation [42] may also induce large scale electromagnetic perturbations in the core plasmas. The roles of the temperature source and resonant magnetic perturbation in the electromagnetic core-edge turbulence spreading are the objective of our future work.