Numerical study of impurity effects on ion temperature gradient modes in tokamak edge plasmas based on the Euler matrix eigenvalue method

Low-Z impurity injection is frequently used for divertor detachment operations in current tokamaks; however, the impurity effects on the main plasma are yet to be fully understood. In this paper, the impurity effects on the ion temperature gradient (ITG) modes in tokamak edge plasmas are investigated based on the Euler matrix eigenvalue method. The eigen-equations with multiple ion species are established from the fundamental gyrokinetic theory, in which each ion species is treated equally. A novel and efficient gyro-kinetic code is developed for this numerical study, and the code’s availability to examine quasi-linear ITG modes is demonstrated by its comparison with existing results. At the pedestal top parameters in Experimental Advanced Superconducting Tokamak high-β p H-mode plasmas, the ITG mode behavior is investigated in pure deuterium plasmas and with impurities. Impurities can induce destabilizing or stabilizing effects on ITG modes, which are determined by the impurity density scale length. The inwardly peaked impurity density profile tends to reduce the ITG growth rate. The effect strength also increases with the impurity charge concentration. The effects of impurity species, including boron, carbon, neon and argon, are also evaluated. Numerical results show that the strength of destabilizing or stabilizing effect inverses with impurity ion charge at the same effective charge.


Introduction
Nuclear fusion is one of the most attractive alternatives to carbon-dependent energy sources, and the tokamak, a magnetic confinement fusion device, may offer a promising way to utilize nuclear fusion peacefully because of its long energy confinement time [1]. For future advanced fusion devices, such as the International Thermonuclear Experimental Reactor (ITER) [2] and China Fusion Engineering Test Reactor [3], high confinement mode (H-mode) will be taken as the basic operation scheme with a divertor geometry. However, unmitigated steady-state heat loads on divertor target plates severely surpass the material limits, which will restrict the discharge pulse length, affect the component lifetime, and even threaten the device safety. To mitigate the steady-state heat loads on divertor target plates, radiative divertor operation by actively seeding low-Z impurities for divertor detachment is believed to be a promising scenario for future advanced fusion devices [4,5]. However, these low-Z impurities may cross the last closed flux surface and enter the main plasma through perpendicular transport, which may affect the main plasma performance and the core confinement [6]. How to maintain a compatible operation between divertor detachment and high core confinement is one of the key issues in magnetic confinement research.
Radiative divertor experiments in the H-mode have been conducted worldwide, and plentiful experimental phenomena have been obtained. Excessive injected impurities usually enter the main plasma and largely increase the core radiative loss, degrading the core confinement [6]. However, some recent experiments have shown that low-Z impurity injection can promote core confinement under certain conditions. Peaking plasma profile appears in some radiative experiments, which compensate for the negative effects of impurities and help maintain or even promote the main plasma performances, such as strong nitrogen injection for complete detachment in ASDEX-Upgrade [6] and argon seeding at high density in JT-60U [7]. The pedestal dynamics and instabilities are also affected by low-Z impurity seeding, along with the behavior change of edge localized modes (ELMs). In JET-ILW [8], nitrogen seeding upon high-triangularity geometry tends to increase the plasma density and reduce the ELM frequency, which recovers the electron pedestal pressure and enhances global confinement. In DIII-D high-β p plasmas with nitrogen seeding [4], a self-organized synergy between the internal transport barrier (ITB) and the edge transport barrier (ETB) leads to a net gain in energy confinement with H 98 ∼ 1.5. In HL-2A [9], 30% neon seeding by supersonic molecular beam injection (SMBI) converts large ELMs into higher frequency and smaller amplitude fluctuations, whereas pure neon or argon seeding reduces the ELM frequency and improves the plasma confinement moderately. Further research in HL-2A [10] shows that the ion temperature in the edge and core regions rises about 20% ∼ 40% after pure neon or argon SMBI seeding. In contrast, the electron temperature is almost unchanged, wherein the suppression of ion temperature gradient (ITG) modes and profile stiffness may be the key factors. Similar phenomena appear in the radiative divertor experiments with neon seeding on Experimental Advanced Superconducting Tokamak (EAST) [11], whereas argon seeding evidently degrades the core plasma performance [12,13].
The underlying mechanism for impurity injection promoting global confinement is not yet fully understood [6], and the behavior change of microkinetic instabilities may be the key factor. In toroidal plasmas, microkinetic instabilities, such as ITG or trapped electron mode (TEM), play an important role in anomalous transport, which significantly influences the global confinement [14]. Gyrokinetic theory with finite radius effects indicates that the impurity mode driven by impurity ions with outwardly peaked density profiles interacts with the ITG mode [15]: ITG enhances the impurity mode, whereas impurity ions have destabilizing effects on the ITG mode. In DIII-D experiments with a low confinement mode (L-mode) negative central shear discharge, neon seeding significantly reduces long wavelength turbulence and heat and momentum transport with a remarkable increase in the energy confinement time [16]. Simulation works with a gyro-Landau fluid mode GLF23 indicates that the results in DIII-D are primarily attributed to the synergistic effects of impurity-induced enhancement of the E × B shearing rate and growth rate reduction of toroidal drift wave turbulence [17]. Similar phenomena appear in the large helical device stellarator [18] that the injection of submillimetric boron powder grains damps lowfrequency turbulent fluctuations with the increase in plasma confinement. However, in tokamak H-mode plasmas, there is only indirect evidence for turbulence suppression by impurity injection; for instance, remarkable increases in the ion temperature with an almost unchanged electron temperature [9][10][11][12][13], whereas direct turbulence diagnostic evidence is rarely reported. Morever, impurity injection also induces the excitation of electromagnetic turbulence [19] and the stabilizing effect on peeling-ballooning modes [20], which influences the pedestal dynamics and ELM behaviors.
Simulation work is necessary to analyze the effects of impurity injection on microkinetic instabilities in tokamak Hmode plasmas. Global gyro-kinetic particle simulation has achieved great success in single ion species conditions, such as GTC [21]; however, it is too expensive under multiple ion species to explain the experimental phenomena efficiently because particle simulation requires tracing massive particles to reduce the intrinsic noise. Several other methods, such as integral equation, Euler initial value simulation, and Euler matrix eigenvalue solutions, are also suitable to study the quasilinear kinetic instabilities and have the capacity to achieve highly similar results [22]. Based on the integral equation, the gyrokinetic code HD-7 has been developed. It has attained success in the study of kinetic instabilities under the two ion species [23][24][25]. In future fusion devices, more than two ion species need to be considered for microkinetic instabilities, because helium ash is inevitable in deuteriumtritium plasmas with impurity seeding for divertor detachment.
In this paper, quasilinear ITG turbulence is investigated with impurities in EAST edge H-mode plasmas, based on Euler matrix eigenvalue solution. Based on the fundamental gyro-kinetic theory, the eigen-equations are derived in detail using ballooning coordinates, wherein each ion species is treated equally. The availability of the newly developed gyrokinetic code is demonstrated by comparing with existing numerical results. For typical pedestal top parameters in the EAST high-poloidal-beta (β p ) H-mode plasma scenario, the effects of impurity seeding on ITG turbulence will be studied and discussed, including the effects of impurity density gradient scale length, impurity charge concentration, and impurity species. The remainder is organized as follows: in section 2, physical models and eigen-equations are presented. In section 3, numerical results are presented and discussed. Finally, section 4 is devoted to a summary and discussion.

Physical model and eigen-equations
Gyrokinetic theory can describe the low-frequency physics accurately with the assumptions: ω/Ω ci ∼ ρ ti /L ∼ k ∥ /k ⊥ ≪ 1 [26,27], where ω is the instability frequency, Ω ci is the main ion cyclotron frequency, ρ ti is the primary ion cyclotron radius, L is the typical scale length, and k ∥ , k ⊥ are parallel or perpendicular wave vector, respectively. In this paper, we focus on quasilinear ion temperature gradient modes, which exhibit a typical low-frequency electrostatic instability. Moreover, nonlinearity is beyond this study's scope.

Basic dynamics
In inhomogeneous plasmas with multiple ion species, the dynamics of low-frequency electrostatic perturbations can be described by the quasineutrality condition: whereñ e is the electron density perturbation, andñ s is the ion density perturbation of particle species s. Due to the low frequency, the adiabatic response to the electrostatic perturbationφ is suitable for electrons; therefore, the electron density satisfies the condition: where e is the elementary charge; n e , T e are equilibrium values of the electron density and temperature, respectively. For ions of particle species s, the non-adiabatic response should be considered [26], and the perturbed densityñ s in an axisymmetric toroidal geometry, like a tokamak, is described by:ñ where n s , T s is the equilibrium value of the ion density and temperature with the charge z s . J 0 is the Bessel function of zeroth order, coming from the gyrophase average. Ω cs = z s eB/m s c is the cyclotron frequency. h s is the non-adiabatic response. The integral of the velocity space is written aś where ⃗ b is the unit vector of the magnetic field, q is the safety factor, and R is the major radius. The non-adiabatic response h s is determined by solving the gyrokinetic equation [26]: where (7) in which ω, ω Ds , ω * sT and ω * s are the mode frequency, magnetic drift frequency, pressure gradient drift frequency and density gradient drift frequency, respectively.
dr is the magnetic shear and θ k is the ballooning angle parameter. L ns = −(d ln n s /dr) −1 is the density gradient scale length of particle species s, and ϵ ns = L ns /R. L Ts = −(d ln T s /dr) −1 is the temperature gradient scale length, and η s = L ns /L Ts .The components of wave vector sat- . Furthermore, in equation (4), i is the imaginary unit, and F Ms is the Maxwell velocity distribution F Ms = ( ms 2π Ts ) 3/2 exp(− msυ 2 2Ts ).

Normalization and eigen-equations
For convenience, we introduce a temperature ratio τ s and a charge concentration f s to represent the equilibrium value ratio: Equilibrium densities also satisfy the quasineutrality condition; thus, there are s f s = 1 and 1/L ne = s f s /L ns . The effective charge number is Z eff = s z 2 s n s /n e = s f s z s . By defining υ ti = T i /m i , the cyclotron radius ρ ti = υ ti /Ω ci is used to normalize the wave vector k ⊥ , k θ , which is the same with Xie's notation [22], whereas it has a slight difference ( √ 2) with Dong's definition [23]. Here, the subscript i represents the main ion, which is , usually, the hydrogen isotope with z i = 1 in present tokamaks. The normalized parameters are written with the hat ask ⊥ = k ⊥ ρ ti ,k θ = k θ ρ ti . Frequencies are normalized by the transit frequency, such aŝ ω = ω υti/R . The potential perturbationφ is normalized byφ = eφ/T e . Therefore, equation (2) transfers to: And equation (3) is written as follows: Equation (4) is normalized as: whereω To further simplify the equations, we define: Multiplyω to above equation, and we rewrite equation (13) as follows: Applying equation (15) to equation (10) and (11) with the quasineutrality condition (equation (1)), we could get the potential perturbation equation: With equation (16), we further calculate the integral as follows:ωĝ The boundary conditions areĝ s ,φ → 0 at θ → ±∞. Equations (16) and (17) constitute the eigen-equations aŝ ωX = MX, where X = [ĝ s ,φ] T , which could be used to study the quasilinear growth rate of the quasilinear ITG instability.

Numerical issues
It has been proven that the fourth-order central difference scheme for ∂ ∂θ is much better than the second order scheme [22]. Therefore, at the jth mesh grid of the ballooning angle coordinate space, ∂φ ∂θ j is set as follows: Notice that: where J 1 is a Bessel function of first order. In this paper, the velocity space of different species of ions is set at the same for simplicity, which can also be different. Thus, the vector size of the non-adiabatic responseĝ s for each ion species is N θ × N ∥ × N ⊥ , where N θ , N ∥ , N ⊥ are the grid numbers of θ, υ ∥ , υ ⊥ space mesh, respectively. The vector X = [ĝ s ,φ] T has N total = N s × N θ × N ∥ × N ⊥ + N θ dimensions, where N s is the number of ion species. Therefore, the sparse matrix M has N total × N total dimensions. Each potential at the θ grid is calculated by equations (17) and (18) with the integral obtained approximately by: where ∆υ ∥,i is the width around the ithυ ∥ mesh grid, and ∆υ ⊥,j is the width around the jthυ ⊥ mesh grid. The Euler matrix eigenvalue method treats each ion species equally, and can handle multiple ion species more easily than the integral equation [23] even though the sparse matrix is larger. Due to the rapid development in scientific computing, the eigenvalue calculation of massive sparse matrix is no longer difficult, especially when the matrix is symmetric [28,29]. In the latter calculations, we set a rough solution as initial guess and obtain the high resolution solutions based on the iterative approach, for example, the MATLAB function eigs. In this paper, we focus only on the maximum growth rate; however, multiple branches of quasilinear ITG modes can be obtained at the same time, which will be further studied in the near future. In addition, the eigen-equation form is the same for different ion species, as shown in equations (16) and (17). As a result, the correctness of the derivation process can be easily guaranteed just by simply comparing with the one species condition [22]. The code check work can also be significantly reduced.

Benchmark
To demonstrate the availability of the newly developed code for quasilinear ITG modes, necessary comparisons with existing results are conducted carefully. Here, large grid amounts are set to promote the resulting accuracy with N θ = N ∥ = N ⊥ = 61 and |θ| < 2π, |υ ∥ | < 3.5,υ ⊥ < 3.5, though it may lose efficiency.
First, the benchmark work is conducted with the gyrokinetic code HD7 [30] and GENE [31] in one ion species condition with the following parameters: Figure 1 shows that the Matrix method and HD7 agree very well, whereas the growth rates from GENE are slightly larger. A similar discrepancy appears in the comparison with GYRO [32], mainly due to the trapped particle effect, which has not yet been included in HD7 and our code. Second, the results in two ion species conditions are compared with HD7 [30] under the following parameters: s = 0.8, q = 1.5, R/L ne = 5, τ s = T e /T s = 1, R/L Ts = 25(s = i, z). The main ion is hydrogen, and the impurity species is fully ionized carbon. In f z = 0.1 condition, different impurity density gradient scale length L nz is considered with the main ion density gradient scale length L ni calculated from the quasi-neutrality condition: s f s = 1 and 1/L ne = s f s /L ns . Figure 2 shows that the results under different L nz agree well, especially for the maximum growth rates. However, the discrepancy appears obvious when the poloidal wavelength is large (k θ · √ 2ρ ti > 1), in which the growth rates from HD7 are larger. A similar discrepancy also appears in pure hydrogen simulations at large k θ ρ ti , and the reason is not yet clear, which is probably due to the grid convergence in HD7, since the velocity space integral of HD7 using Gauss method is difficult to achieve with high accuracy [22].
The results under the other cases are also compared and all agree well, such as figure 8 in [15], figure 7 in [22] and figure 1 in [32]; however, the absolute values differ slightly with HD7 at large poloidal wavelengths. These comparisons satisfactorily indicate the availability of our newly developed code to study the quasilinear ITG modes in pure plasmas or with impurities under toroidal geometry.

Grid convergence and mode structure
Typical parameters in EAST high-β p H-mode plasma scenario at ρ = 0.85, that is pedestal top [33], is used in later discussions unless otherwise stated:ŝ = 3.5, q = 4.5, R/L ne = 6, R/L Te = 12, τ i = 1. The grid convergence is tested by scanning the grid number in pure deuterium plasmas. The growth rate error (γ − γ 0 )/γ 0 is calculated, in which γ 0 is the growth  rate calculated with N θ = N ⊥ = 61. As shown in figure 3, the error rapidly decreases below 0.5% when grid number N θ or N ∥ surpasses 20, indicating nice grid convergence. Moreover, angle grid convergence is much better than velocity grid convergence, implying that the main error source is the velocity space integral´d 3 υ. Large grid number usually reduces the  computational efficiency. Therefore, n θ = n ∥ = n ⊥ = 41 may be a suitable value to ensure the accuracy and efficiency with the mesh region |θ| < 2, |υ ∥ | < 3.5,υ ⊥ < 3.5, which is set as the default value unless otherwise stated. The standard time consumption for each dot with a single computer kernel is ∼1 min for the one ion species condition and ∼3 min for the latter for two ion species condition.
This potential ϕ(θ) can be used to display the mode structure in ballooning space. Three typical structures are shown in figure 4 with different k θ ρ ti . The eigen-function width is relatively large when k θ ρ ti is small. Moreover, the structure seems to become narrow as k θ ρ ti increases, which is consistent with Li et al's results [24]. The structure width in the ballooning space is a crucial parameter relevant to Landau resonance and is affected at large k θ ρ ti by the impurity density gradient scale length [24], which is beyond the scope of this work and will be further investigated in the near future.

Discussions on temperature ratio τs = Te/Ts
Limits to the diagnostic capacity make it hard to obtain impurity density and temperature profiles, especially in the edge region. Therefore, parameter scanning work is unavoidable, such as temperature ratio τ z , impurity density gradient scale length L nz , and charge concentration f z .
First, we study the impact of the ion temperature ratio τ i on the growth rate in pure deuterium plasmas, with the same temperature gradient scale length L Ti = L Te and k θ ρ s = 0.5, where ρ s = C s /Ω ci . Figure 5 shows that the growth rate increases rapidly with τ i , reflecting that the main ion temperature exerts a significant effect on the ITG modes. Then, deuterium plasmas with Ne 10+ impurities are considered with τ i = 1, L Tz = L Ti = L Te , f z = 0.1, R/L nz = −10, k θ ρ s = 0.5, in which the effective charge Z eff = s f s z s is 1.9. Figure 5 shows that the impurity temperature ratio τ z only has a slight influence on the ITG modes. Therefore, it is important to obtain the measurement results of main ion temperature to study the effect of impurities on ITG mode in edge plasmas, whereas impurity temperature seems less important.
At the EAST pedestal top, it may be suitable to set the same temperature for different particle species, that is T z = T i = T e , which will be discussed here in detail. In tokamak plasmas, the typical heat transfer time τ kj from particle species j to particle species k can be written as [34]: where m, z, T, n are particle mass, charge, temperature and density, respectively.
where Z is the impurity charge. τ i = T e /T i is close to unity in the edge region for EAST H-mode plasmas, as measured by Thomson scattering diagnostics and charge-exchange recombination spectroscopy [33,35]. Neon is often used in EAST radiative experiments. At the pedestal top, the plasma temperature is about 0.5 KeV and neon is almost ionized completely [36], thus, τ ie /τ zi ∼ 10 m i /m e ≫ 1 is obtained. This implies that the time for main ions and electrons to achieve thermal equilibrium is much longer than the time for impurities and main ions to achieve thermal equilibrium. Therefore, it is convincing for EAST edge plasmas to set τ z = τ i = 1, which is used in later discussions with L Tz = L Ti = L Te unless otherwise stated.

Effects of impurity density gradient scale length Lnz
Typical spectrums with or without neon impurity in deuterium plasmas are shown in figure 6, and the impurity has a significant influence on the growth rate, while only has slight effects on the real frequency. The growth rate curves always rise up to the peak value and then decline with poloidal wavelength. For pure deuterium plasmas, the maximum normalized growth rate is 0.54 at k θ ρ ti = 0.5. With an outwardly peaked neon density profile R/L nz = −10, the growth rate becomes remarkably higher and reaches the maximum value 0.70 at k θ ρ ti = 0.54. With an inwardly peaked neon density profile R/L nz = +10, the growth rate obviously becomes lower and reaches the maximum value 0.33 at k θ ρ ti = 0.44. The slight shift in the growth rate peak position can also be observed in figure 2, in which the underlying reason is not clear.
For real frequencies, the impurity effects are slight and only become obvious at large k θ ρ ti , as shown in figure 6(b). When k θ ρ ti > 0.5, the real frequency becomes slightly higher in R/L nz = ±10 compared to pure deuterium plasma condition. It is worthy to mention that iterative approach to obtain the eigenvalue sometimes has trouble when the growth rate nearly reaches to zero, which can be seen at the end of blue curve (R/L nz = +10). Fortunately, this defect has no effect on our study because only large positive growth rates deserve attention, which means that the mode is unstable and can be easily excited. In the future, other approaches to obtain the eigenvalues will be tested to avoid this problem.
To further study the effects of impurity density profile, the impurity density gradient scale length L nz is scanned with f z = 0.1, k θ ρ ti = 0.5. Figure 7 shows that impurities with outwardly peaked density profile, that is L nz < 0, have destabilizing effects on ITG modes, especially when |L nz | is small. On the other hand, impurities with inwardly peaked density profile, that is L nz > 0, had stabilizing effects, especially when L nz is small. These conclusions are consistent with previous work on tokamak plasmas [15,23,24] or reversedfield pinch plasmas [37].

Effects of impurity charge concentration fz and impurity species
The strength of impurity destabilizing or stabilizing effect on ITG modes is also determined by the impurity concentration. We scan the Ne 10+ concentration at k θ ρ ti = 0.5 and R/L nz = ±10, as shown in figure 8. With R/L nz = −10, the growth rate increases rapidly with impurity concentration f z . With R/L nz = +10, the growth rate decreases rapidly with f z , and the ITG modes are almost suppressed completely when f z > 0.2. The impurity charge f z is directly associated with the effective charge Z eff = s f s z s ; for instance, Z eff = 1.9 at f z = 0.1 and Z eff = 2.8 at f z = 0.2. In addition, the impurity concentration only has a slight influence on the real frequency. For instance, the real frequency only has <10% decrease when f z = 0.1 compared to pure deuterium plasma conditions. Impurity species may also have an important impact. For instance, different impurity species, including boron, carbon, neon and argon [38], are adopted in EAST experiments for radiative divertor or ELM control with different phenomena. Here, we scan the concentrations of different impurity species under the same parameters k θ ρ ti = 0.5, R/L nz = ±10, as shown in figure 9. It is worth mentioning that argon may not be fully ionized at the EAST pedestal top because argon usually appears as Ar 16+ around 0.5 keV [36]. The numerical results show that under the same f z , different impurity species have similar destabilizing or stabilizing effects on ITG modes.  Obvious differences only appear at high f z , where argon and neon have almost the same effect, while boron has the weakest effect.
However, impurity species charge differs a lot, and the tolerable superior limit of impurity concentration in tokamak plasmas also varies a lot. Instead of impurity charge concentration f z , effective charge Z eff may be a better choice to study the effect on ITG modes, which are directly associated with radiative loss and are an important parameter in tokamak experiments. Considering the same effective charge, the spectrum of different impurity species are shown in figure 10. At Z eff = 2, the charge concentration is fz = 1/4, 1/5, 1/9, 1/15 for B 5+ , C 6+ , Ne 10+ , Ar 16+ respectively. The strength of destabilizing or stabilizing effects on ITG modes follows the order: boron, carbon, neon and argon, inverse with impurity ion charge. This may explain why all tokamaks prefer the inert elements neon and argon for radiative divertor operation instead of other low-Z impurities because the effects of neon and argon on main plasmas is relatively weaker with the same radiative loss.

Summary and discussion
In this paper, a numerical study of impurity effects on quasilinear ITG modes in tokamak edge plasmas is conducted based on the Euler matrix eigenvalue method. For multiple ion species, the derivation of eigen-equations is presented in detail from the fundamental gyrokinetic theory. Based on the MATLAB platform and iterative approach, an efficient gyrokinetic code has been developed in this study. Because each ion species is treated equally, the code has the capacity to handle more than two ion species condition. Comparisons with existing numerical results demonstrate the code's availability to study the quasi-linear ITG modes in toroidal plasmas. Nice grid convergence is also confirmed.
At typical pedestal top parameters (ρ = 0.85) in EAST high-β p H-mode plasmas, the ITG mode behavior is investigated in pure deuterium plasmas or with impurites. Numerical results show that the main ion temperature has a significant influence on the growth rate, whereas the impurity ion temperature seems less important. The destabilizing or stabilizing effect of impurities on ITG modes is mainly determined by the impurity density gradient scale length. The ITG modes are destabilized with L nz < 0, while becomes suppressed with L nz > 0. The strength of destabilizing or stabilizing effects also increases with impurity charge concentration. In addition, the effect strength on ITG modes at the same Z eff follows the order: boron, carbon, neon and argon, inverse with the impurity ion charge.
The numerical results show that the ITG suppression might be the key factor for the phenomenon that the core ion temperature has a remarkable increase with nearly unchanged electron temperature observed in EAST H-mode plasmas with neon seeding [11][12][13], which is similar to that in HL-2A experiments [10]. Nonetheless, it may still be a challenge to directly utilize the newly developed gyrokinetic code to explain EAST experimental phenomena perfectly, because some parameters significantly influencing ITG modes cannot be measured yet in experiments, such as the impurity density gradient scale length. Besides, the slight discrepancy with HD7's results at large poloidal wavelengths also deserves further investigation, and the conclusion in Li et al's results [24] also needs re-examination of the finding that the eigenfunction width of ITG modes at large poloidal wavelengths is affected by impurities. Although ITG modes play an important role in the main plasma, especially in the edge region, the TEM dominates the anomalous transport in the core region in EAST high-β p H-mode plasmas [33]. The dissipative TEM is an important source of electrostatic turbulence in tokamak pedestal plasmas [39], which may be the mechanism for the edge coherent mode observed in EAST H-mode discharges [40]. In the near future, further study of impurity effects on TEM will be conducted by either analytic theory or gyrokinetic simulation.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).