Theoretical and Global Simulation Analysis of Collisional Microtearing Modes

This study delves into Microtearing Modes (MTMs) in tokamak plasmas, employing advanced simulations within the BOUT++ framework. The research, centering on collisional MTMs influenced by the time-dependent thermal force, enhances our understanding of plasma dynamics. It achieves this through the simplification and linearization of control equations in detailed linear simulations. The study meticulously evaluates various conductivity models, including those proposed by Larakers, Drake, and Hassam, under diverse plasma conditions and collision regimes. A notable achievement of this research is the derivation of a unified dispersion relation that encompasses both MTM and Drift-Alfven Wave (DAW) instabilities. It interestingly reveals that DAW and MTM exhibit instability at different proximities to the rational surface. Specifically, MTMs become unstable near the rational surface but stabilize farther away, whereas the drift-Alfven instability manifests away from the rational surface. Further, the study re-derives MTM dispersion relations based on Ohm’s law and the vorticity equation, providing a thorough analysis of electro-magnetic and electrostatic interactions in tokamaks. Global simulations demonstrate an inverse correlation between MTM growth rates and collisionality, and a direct correlation with temperature gradients. The nonalignment of the rational surface with the peak ω ∗ e stabilizes the MTMs. Nonlinear simulations highlight electron temperature relaxation as the primary saturation mechanism for MTMs


I. INTRODUCTION
The microtearing mode (MTM) is an electromagnetic micro-instability that significantly affects thermal transport and energy confinement in the tokamak plasma.In standard tokamak H-mode discharges, such as those observed in DIII-D 1 , JET 2 , and others, MTMs are frequently observed.
MTM can dominate electron thermal transport in the JET ITER-like wall pedestal under certain conditions 3 , plays a role in determining the behavior of electron temperature during an ELM cycle 4 , and in double barrier operations, it generates a substantial electron thermal transport at the Inner Transport Barrier (ITB) 5 .In the high-β regime, MTMs may become more prominent as KBMs enter their second stable regime, thereby playing a significant role alongside the behavior of KBMs.Therefore, MTMs can become the dominant low-k θ modes as β increase in spherical tokamak devices like NSTX 6,7 .
MTM was initially recognized as a temperature-gradient-driven tearing mode in the ST Tokamak in the 1970s.Using a variant method, Hazeltine et al. 8 first proposed a linear kinetic conductivity model for this type of mode.This kinetic theory was further developed by Drake et al. 9 , with the model being extended to a semi-collisional regime.More recently, Larakers et al. 10,11 have successfully expanded the theory to include the global radial variation effect, providing a comprehensive explanation for the separate bands observed in the MTM frequency spectrum.
BOUT++ is a C++ framework specifically developed for solving initial value partial differential equations 12 related to magnetized plasmas.It ensures a clear delineation between numerical and physics development coding, thereby enhancing the efficiency of the simulation code development process.Successfully applied in tokamak boundary plasma simulations 13,14 , including edge-localized modes, boundary turbulence, and blobs.Plasma models in BOUT++ have been validated with experimental data and results from other codes.
In this work, we have extended the capabilities of the BOUT++ code to include the simulation of the MTM.This paper presents our global simulation of the collisional MTM in a shifted-circle tokamak geometry.We have implemented Hassam's fluid model for the MTM 15 , incorporating electron inertia and a time-dependent thermal force into Ohm's law equation.Both global linear and nonlinear simulations are performed to investigate the MTM characteristics.This paper is organized as follows.In Section II, a brief survey of existing theoretical models on MTM is provided.Section III describes the simulation model and parameters used in this study in detail.Linear simulation results are summarized in Section IV as we explore the mode BOUT-MTM structures and the relation between frequency and growth rate spectrum and collisionality and temperature gradient.The effect of alignment on the growth rate is investigated as well.Nonlinear MTM simulations are presented in Section V. We find the saturation mechanism of MTMs due to temperature profile relaxation and compare the heat transport resulting from E × B convection and magnetic flutter.Finally, Section VI summarizes our findings.

II. THEORY OF MICROTEARING MODE
In this section, we delve into the theory of MTMs, examining its dispersion relation and exploring its connection to classical tearing modes as well as DAW instabilities.

A. Review and comparison of MTM models
MTM is related to the resistive current layer physics.Different models are characterized by the conductivity, σ , in Ohm's law.
Recently, the classical Hazeltine electric conductivity model 8 has been extended to include radial variation by Larakers et.al 11 .If using a Lorentz Gas collision operator neglecting the electronelectron collision, where ξ = v ∥ /v is the pitch angle, ν ei is the collision rate, ν ei = 16 √ πn e Ze 4 lnΛ/3m 1/2 e v 3 te , where v te is the electron thermal speed, the conductivity has a form of, , where L n = |∇n e |/n e , L T = |∇T e |/T e are characteristic lengths of electron density and temperature gradients, k θ is the poloidal wave number.In the fluid response limit, k ∥ v te ≪ ω, ν ei , the coefficients are, Then Larakers' model in different limits is compared with Drake 16 and Hassam's 17 models.
In the ω ≪ ν ei limit, the coefficients L 11 , L 12 reduce to, After substituting into Eq.(3), the conductivity is, which has the same form as Drake's model and Hassam's model.
In the ω ≫ ν ei limit, following the similar way, the coefficients are After substituting into Eq.(3), the conductivity is, It is the same as Drake's conductivity in the low collisional regime.
In summary, for local assumption around the rational surface, Larakers' conductivity model reduces to Drake's model in both ω ≪ ν ei and ω ≫ ν ei limits, and reduces to Hassam's model in ω ≪ ν ei limit.Larakers' model thus bridges the two limits and also includes k ∥ effects.

B. Dispersion relations of MTM
The dispersion relation of MTM is derived from Ohm's law and vorticity equation (or the quasi-neutral condition) and has already been performed under slab or other geometries by many researchers 8,9,17 .If performed in a slab geometry, after linearization and Fourier transformation, the vorticity equation and Ohm's law take the following forms, where the tilde means perturbed value, x denotes the radial distance from the rational surface, and y is the orthogonal direction on the magnetic surface.
µ 0 is the vacuum permeability.k ∥ = k y x/L s is the parallel wave number, k y is the binomial wave number, approximately k y = m/r, and m is the poloidal mode number, L s is the local magnetic shear length.Using a full collision operator, the conductivity of Hassam's mode is, where, α = 0.8, α′ = 0.54 are constant numbers.These coupled equations can be solved in the inner and outer regions separately, and using a matching condition to obtain the dispersion relation.

BOUT-MTM
In the inner region, the electromagnetic perturbation dominates.The electrostatic potential plays a role in short-circuiting the electric field.When E ∥ ≡ iω A ∥ − ik ∥ φ = 0, a characteristic length can be solved, ∆ E = ω A ∥ L s /k y φ , where L s is the magnetic shear length.Out of the region, |x| > ∆ E , the electric field is short-circuited.If ∆ E ≫ ∆ j , where ∆ j is the width of the current layer determined by σ (x), Eqs.(15) reduces to the equation, The validation of this assumption is 18 .Using a shooting method, the mode structure is obtained and shown in FIG. 1.
However, across the current layer, there is a discontinuity of A ′ ∥ due to the current sheet, and the prime denotes derivative in the radial direction.This discontinuity can be represented by a tearing parameter, calculated under the constant-A ∥ assumption, . The ∆ ′ is solved in the outer region.Approximately A ∥ ∼ e −k y |x| , and ∆ ′ = −2k y 19 .Then the solution of Eq.( 18) can be obtained from an integration within the current layer width, ∆ j 18 .
Furthermore, in general cases, without ∆ E ≫ ∆ j assumption and considering the current drive, J ′ 0∥ , of that classical tearing mode, it is necessary to combine with Eq.( 16) to obtain the dispersion relation in a form, where is the growth rate of classical collisional tearing mode 20 .
We do not repeat the derivation but just show the result here.This shows that MTM could be unstable even when ∆ ′ is negative.
When ω * e ≫ γ c , the influence of γ c on the solution of Eq.( 19) is only a minor correction to its main value.It is applicable under most present tokamak operation parameters.And the main value of the solution can be solved from, 8,11 σ = 0, (20)   This indicates that a dispersion relation σ = 0 can be used to calculate the mode frequency and growth rate.
If performed in local analysis, it becomes evident that the mode transitions from a microtearing mode to a drift Alfvén wave instability with an increase in k ∥ .Let ∇ ⊥ −→ −ik ⊥ , Eqs. ( 15)-( 16) BOUT-MTM are reduced to a dispersion relation, Choosing a set of parameters, k 2 ⊥ η/µ 0 ν ei = 0.05, ω * Te /ν ei = 3ω * ne /ν ei = 0.3, the dispersion relation Eq.( 21) can be solved numerically in slab geometry.As shown in FIG. 2, there are 3 modes.The red curve is recognized as drift Alfvén wave instability.It is unstable with finite k ∥ , in other words, this instability arises at the location where is a little bit shifted away from the rational surface.As shown in FIG.2(b), there is another mode becoming unstable in the small k ∥ regimes and becoming stabilized further away from the rational surface.This instability is recognized as an MTM driven by the time-dependent thermal force.DAW is unstable with finite k ∥ , as a result, DAW and MTM are unstable at different locations to the rational surface.

III. FLUID MODEL AND SIMULATION SETUP
In this section, we introduce the control equations, equilibrium profiles, and parameters used in this study.

A. Microtearing mode model for simulation
To minimize the impact of existing fluid-based turbulence model 21 within BOUT++ framework, Hassam's nonlinear fluid model 15 is adapted for the global MTM simulation.In this model, the electron parallel momentum equation includes a time-dependent thermal force ∂ (∇ ∥ T e )/∂t term and thus is slightly different from the original Braginskii model 22 , where η is the parallel resistivity.Replacing A ∥ with B ψ and dividing en e B for both sides of Eq.( 22), this equation becomes, Here we take the semi-electromagnetic assumption 23 and the parallel gradient is written as Because the MTM is driven by temperature gradients, we simplify the model to include the fundamental physics and combine it with the vorticity equation in 6-field model 21 that with finite ion Larmor radius effect (denoted by the gyro-viscosity terms) included.Finally, the control equations are reduced to a three-field model,

BOUT-MTM
Here, b = B/B denotes the direction of magnetic field, V E = b×∇Φ/B is the drift velocity, Φ = φ + Φ 0 is the total electrostatic potential, φ and Φ 0 are the perturbed and equilibrium electrostatic potential respectively, κ = b • ∇b is the magnetic curvature, Ω i is the ion gyro-frequency, µ ∥ , µ ⊥ are the parallel and orthogonal viscosity for numerical stability 24 .Spitzer-Harm resistivity 22 is implemented.
η H is the hyper-resistivity 25 .Other coefficients are set as α = 0.8, α ′ = 0.54, α ′′ = 1.24 17 .C ν ,C T are two free coefficients used for collisionality and gradient dependence study.In the following section, a scan of these coefficients is performed to test the frequency and growth rate spectrum versus collisionality and temperature gradient.Ion parallel velocity is neglected, therefore, the perturbed current is approximately j ∥ = −en e V e∥ .A flux limit closure 21 with the harmonic average of free streaming and conductivity model for heat flux is used, where α H is a weight factor and we call it free streaming parameter, κ FS is the free streaming thermal transport coefficient, and κ SH is the conductive thermal transport coefficient, 23 L ∥ = q L is the parallel correlation length.When α H is close to infinity, the flux limit model will reduce to the conductivity model.In the simulation performed in this article, this factor is chosen to be α H = 0.05 by default.
If the diamagnetic drift effect is considered (i.e., finite κ), the second term in the right-hand side of vorticity Eq.( 24) must be included to maintain equilibrium 26 .Note in our simulations, cold ion model was adapted so that the perturbation of ion density and temperature vanish, making the second term in the right-hand side of Eq.( 27) always zero.At the same time, the E r shear flow, , is included to balance equilibrium ion diamagnetic drift for particle conservation 27 .
Also E r shear flow introduces the Doppler shift effect,

B. Simulation setup
In our simulation, we utilize a shifted circular geometry equilibrium generated by CORSICA 28 .
The equilibrium profile is displayed in FIG. 3 The safety factor profile with 2.6 < q < 3.1 for the simulation domain is intentionally set to be flatter than usual to distinguish rational surfaces clearly for discussions.In contrast to the standard CORSICA-generated equilibrium, our equilibrium features a lower current profile and a stronger magnetic field.It's worth mentioning that Seto et al. 29 have made significant developments in BOUT++ to enable accurate low-n mode simulations, although this version may come with slower code execution speed.So, we do not perform very low-n modes simulation.
A field-aligned coordinate is implemented in BOUT++ 30 , where

IV. GLOBAL LINEAR SIMULATION OF MTM
For the linear simulation, in order to easily understand the physics, Eqs.( 24)-( 29) are further reduced.After linearization, the control equations for linear simulation are, where b 0 denotes the direction of the unperturbed magnetic field.And Eqs.( 27) and ( 28) can be written in forms, BOUT-MTM Electrostatic potential φ and normalized parallel magnetic vector potential ψ are solved from Eqs. ( 36)-( 37) by the 'invert-laplace' Laplacian equation solver of BOUT++.This solver is introduced in detail in Appendix A. In our simulations, we exclude the current driven term by default, B 0 × ∇ ψ • ∇J ∥0 , to focus on the temperature gradient driven MTM.

A. Mode structure
A linear simulation of n = 8 is performed.FIG. 4 displays the poloidal and radial distributions of ψ, φ and Te at t = 800τ A in the linear simulation.These mode structures exhibit typical tearing parities, with ψ being even parity and φ odd parity to the rational surface q = 2.75.In FIG. 5, a Poincaré plot reveals the formation of magnetic island chains around rational surfaces, providing compelling evidence for the identification of the MTM.

B. Dependency on collisionality and temperature gradient
The dependence of linear growth rate and frequency on collisionality and temperature gradient is compared with theoretical expectations.In Sec.II B, the conductivity is calculated, Eq. ( 17).The lowest-order dispersion relation is employed, setting σ (ω) = 0, to derive the analytical expressions for growth rate and frequency, This solution from the lowest order dispersion relation can capture the dependence trend.
We perform a scan in C ν with C T = 1.0 for the n = 8 mode to investigate the effects of collisionality.As shown in FIG.6(a), the growth rate decreases as collisionality increases.In comparison with the simulation, the local theory tends to overestimate the linear growth rate.This overestimation is attributed to the theory being local and not considering magnetic shear effects.The declining growth rate in the collisional regime with increasing collisionality is a distinct characteristic that sets it apart from the classical tearing mode.In FIG. 6 has a relatively minor impact on mode frequency in the collisional regime.The scan in C ν indicates that collisionality plays a stabilizing effect for MTM.However, it is valid only in the strong collisional regime.MTM is destabilized by appropriate collisionality and temperature gradient.
In the weak collision regime, the growth rate is proportional to collisionality, indicating that collisionality destabilizes MTM 5 which distinguishes it from drift-wave instabilities.
Another scan is performed in C T using the same approach, as presented in FIG.6(c).As C T increases, the growth rate shows a proportional increase.This is because, with the higher C T , effectively both the thermal force and the time-dependent thermal force increase.These two forces drive the unstable MTM.In other words, when C T is higher, there is more effective free energy from the electron temperature gradient for electron motion and leads to greater instability in the MTM.This provides evidence that the MTM is primarily driven by the temperature gradient, as opposed to the current gradient characteristic for the classical tearing mode.The mode frequency is noticeably influenced by C T , as shown in FIG.6(d).An increase in C T leads to a corresponding increase in the MTM frequency.This observation suggests that the MTM frequency is directly proportional to the electron diamagnetic frequency.
The effects of curvature and ion diamagnetic drift (correspondingly with the E r shear flow induced) are investigated.As illustrated in FIG.6, the magenta stars denote the cases with curvature drive (denoted by the third term on the RHS of Eq.( 33)) included and the green squares denote the cases with both curvature and ion diamagnetic drift effect included.The curvature almost has no influence on the linear growth rate and mode frequency indicating that the MTM is not a curvature-driven mode.
However, when the ion diamagnetic drift effect is considered, the mode frequency and growth rate are well influenced. where Then we obtain a conductivity, Following a similar way, we can obtain the mode frequency and growth rate in the lowest order As illustrated in the previous section, the equilibrium E r flow balances the ion diamagnetic drift flow, V E 0 = ∇P i ×b B 0 .Correspondingly, ω φ 0 = −ω * i = ω * e .In the simulation, the density is constant, therefore, ω φ 0 = ω * Te .Approximately, we obtain the relation, This relation indicates that considering the diamagnetic drift effect, along with the inclusion of equilibrium E r shear flow, results in a mode frequency approximately 1.56 times higher than when the diamagnetic effect is not considered.This relation aligns with the simulation results displayed calculation 11 , Curie's SLiM model 32,33 , and also Hassan's pedestal simulation 34 .
A linear simulation is conducted to confirm this prediction.The equilibrium temperature and The red dotted line denotes the cases without curvature effect and peeling drive effect (related to J ∥0 and denoted by the second term on the RHS of Eq.( 33)), the blue dotted line denotes the cases with peeling drive effect included, and the magenta line denotes the cases with both curvature and peeling drive effect included.
The rational surface of n = 8 almost perfectly aligns with the peak of ω * e , while the rational surface of n = 9 is slightly shifted.Under the discussion in the previous section, when the drive increases, the growth rate is correspondingly higher.And a higher toroidal mode number leads to BOUT-MTM

V. GLOBAL NONLINEAR SIMULATION OF MTM
In this section, we conduct an initial nonlinear simulation to explore the nonlinear saturation mechanism of MTM.We will compare radial thermal transport resulting from E × B convection and magnetic flutter.The nonlinear simulation retains all terms except the peeling drive in Eqs.( 24)-( 26) and utilizes the same equilibrium as in the previous section.Also, the E r shear flow is included by default in the following simulations.
Due to the limitations of the collisional Hassam model, accurate calculations of high-n modes are not feasible.Therefore, the nonlinear simulation includes only n = 4, 8, 12 modes while retaining the zonal component of temperature perturbation.As the rational surfaces for these three modes align with the peak ω * e , the nonalignment stabilizing effect is excluded.The initial perturbation is magnetic and takes the form of, where δ x = 0.1, δ y = 0.3, x 0 = 0.5, y 0 = 0.5, k z = 2π/4, ϕ n is a random phase.

A. Saturation of collisional MTM due to profile relaxation
The surface-averaged perturbations of the absolute values of electrostatic potential φ and normalized parallel magnetic vector potential ψ are presented in FIG.9(a).After approximately 1000τ A time, the MTM enters a saturation state.Notably, the temperature profile changes the evolution.We depict the initial temperature profile with a dashed black curve, and the profile after 3000 Alfvén times with a solid blue curve, as shown in FIG.9(b).The corresponding normalized gradient profiles are shown as dotted lines with the same colors.
It is evident that the temperature profile changes, with the peak gradient decreasing occurring around the q = 2.75 rational surface.Consequently, the linear-driven force decreases, leading to a gradual cessation of perturbation growth, ultimately resulting in mode saturation.In our nonlinear simulations, Only the zonal component of temperature perturbation Te is kept, the zonal component of φ and ψ are removed, and correspondingly, zonal flow and zonal field are not included.
We are not going to discuss how they affect the saturation process and will discuss it in another paper.In comparison with the scenario without the equilibrium E r shear flow, perturbations exhibit a reduced linear growth rate but maintain a similar amplitude in the nonlinear phase.Without the The corresponding normalized temperature gradient profiles are plotted in dotted lines with the same color.
equilibrium E r shear flow, the temperature gradient exhibits a more significant decrease.

B. Electron radial thermal transport
MTM is one candidate for anomalous electron thermal transport in the H-mode discharge.The electron radial thermal transport can be divided into two parts, one is the E × B convection part and the other is the magnetic flutter part.The total heat flux has a form 35 , where p e ≡ n e T e is the electron pressure, B r = b 0 × ∇ ψ is the radial component of the perturbed magnetic field, the bracket ⟨...⟩ denotes the magnetic surface average via the unperturbed surfaces.
The first term at the right-hand side of Eq.( 44 as, To illustrate the MTM thermal transport, we use the case with E r shear flow as an example.perturbation compared to that in the linear phase at t = 600τ A .There forms a magnetic island chain around the q = 2.75 rational surface and a stochastic magnetic field in the major domain between q = 2.67 and q = 2.833 rational surfaces.For the electrons, the parallel conductivity is large, and the parallel thermal transport has a radial component due to the magnetic perturbation.Therefore, the electron thermal transport is enhanced, and the magnetic flutter component dominates the electron thermal transport.In our simulation, we observe an asymmetric spreading of MTM turbulence.There is an evident inward spreading as indicated in FIG.10(a).Such an inward spreading intensifies the thermal transport near the pedestal top.As shown in FIG.10(b), the thermal transport coefficient at ψ N = 0.82 increases after a temporal delay and eventually reaches the same level as that at ψ N = 0.85.

C. The influence of the free streaming parameter α H on electron radial thermal transport
In this section, we discuss how the free streaming parameter α H in the flux limit heat flux model affects the electron thermal transport in the global nonlinear simulations.
The magnetic flutter component dominates electron thermal transport.This component can be written in an explicit form with perturbation and equilibrium parts separated.
Here, the subscript DC denotes the zonal component, and AC denotes the non-zonal component.
Te | DC change the electron temperature profile, ∇ T e0 + Te | DC ≡ ∇T e is the electron temperature gradient at the moment.The first term in the bracket of Eq.( 46) describes the parallel dissipation BOUT-MTM and decreases the magnetic flutter electron thermal transport.The second term dominates the radial thermal transport, it describes an effective radial thermal diffusion when the magnetic filed is perturbed.The third term has little contribution to the heat flux after the magnetic surface average.A magnetic flutter thermal transport coefficient can be defined, Approximately, χ e = χ mf e .The flux-limit κ f l in Eq.( 31) can be rewritten in a form, where α H L ∥ is an effective correlation length along the magnetic field line, and λ mfp = v te /ν ei is the collision mean free path describing a decorrelation length due to collision.Then, a total effective correlation length is defined, Using the equilibrium parameters, L eff is calculated and presented in FIG.12(a).α H L ∥ is comparable to the device size and much smaller than the collision mean free path and determines the L eff .
Therefore, L eff is roughly proportional to the free streaming parameter α H .
Because the second term in the bracket in Eq.( 46) dominates the Q mf er ∝ χ e , electron thermal transport can be rewritten as with coefficient C mf is less than unity.
The electron parallel heat flux q e∥ plays two roles.It doesn't only contribute to a large thermal transport but also dissipates the perturbations.As α H increases, χ f l ≡ v te L eff (α H ) increases correspondingly.But a larger χ f l makes perturbation, like Br , better damped as shown in FIG.12(b).
An increase in α H from 0.03 to 0.15 makes the perturbation Br /B 0 decrease from 1.84 × 10 −4 A thermal transport coefficient in this form indicates that thermal transport can be conducted by a stochastic magnetic field 36,37 However, from the discussion above, we learn that electron thermal transport is strongly dependent on the model of parallel heat flux.Our flux-limit heat flux model can capture some instinct characteristics of microtearing turbulence like the transport scaling law Eq.(51), but quantitatively the results are determined by the α H we chose.

VI. SUMMARY
In our research, we have developed a fluid simulation code for studying the collisional MTM within the BOUT++ framework.We should note that BOUT++ has been updated to include the capability for simulating this mode and make it possible to investigate how MTM interacts with other ion modes like the kinetic ballooning mode.
In our linear simulations, we observe that the MTM exhibits a tearing parity in its mode structure.We explore the dependencies of growth rate and mode frequency on collisionality and temperature gradient.Within the collisional regime, the linear growth rate of the MTM decreases as collisionality increases, while it remains directly proportional to the electron temperature gradient.The E r shear flow causes a reduction in the linear growth rate especially when the collision is strong.Additionally, we investigate the non-local alignment effect, revealing that the MTM becomes unstable only when the rational surface aligns with the peak of the ω * e profile; and the nonalignment effect serves as a stabilizing factor.
In our nonlinear simulations, we find that the saturation of the MTM is primarily due to electron temperature profile relaxation with decreasing temperature gradient around the rational surface.Furthermore, we observe that electron thermal transport from magnetic flutter significantly exceeds that generated by E × B convection due to the change in the topology of the magnetic field by MTM.Evident inward turbulence spreading is observed to intensify the thermal transport near the pedestal top.In our simulations, the electron thermal transport is quantitatively influenced by the heat flux closure or the free streaming parameter α H in the flux-limit heat flux model.
In NSTX H-mode plasmas, the confinement time is observed to be inversely proportional to the electron collisionality 38  Hassam's model is valid in the strong collision regime, and it limits our simulation to reach BOUT-MTM out to the weak collision regime, where ω/ν ei ≈ 1.Unfortunately, in the edge pedestal of modern Tokamak devices, the electron diamagnetic frequency is close to the collision frequency, which means for the MTM ω/ν ei ≈ 1.In the gyro-kinetic simulations of today's Tokamak devices, the most unstable mode is found to be ω/ν ei ≈ 1 39 .We also recognize additional mechanisms for MTM saturation, such as turbulence regulation by Zonal Flows (ZFs) and energy cascading due to nonlinear toroidal mode couplings 40 .Our current BOUT++ simulations primarily focus on the relaxation of the electron temperature profile, partly because ZF calculations are still in development, resulting in their exclusion.Prioritizing the integration of ZFs in future research is essential to fully understand their regulatory role.We also plan to delve into the effects of energy cascading on MTM saturation, acknowledging its potentially significant impact.Our simulations are currently limited by the use of Hassam's model, which assumes ω * e /ν ei ≪ 1, confining our analysis to lower-n modes.Our objective is to advance a fluid model of MTM that can handle arbitrary frequency responses, enabling the inclusion of higher-n modes to extensively assess energy cascading.
In our simulation, we are focused on investigating the unstable microtearing modes.To simplify the analysis, we have considered a limiting case with infinite η e , thereby omitting the complexities associated with the density gradient and its stabilizing effects 31 .This method enables us to pinpoint the maximum growth rate, indicative of the most unstable limit.Moving forward, we plan to integrate the density gradient effects into our research to comprehensively understand their impact on MTM stability.
From the discussion above, we find that the electron thermal transport of tearing mode turbulence in fluid simulations is strongly influenced by the heat flux model/closure.Simulations with the Landau-fluid closure 21,41 and perhaps more sophisticated parallel heat flux models [42][43][44] may provide quantitatively more accurate calculation of electron thermal transport.Other nonlinear effects like energy cascading due to toroidal toroidal mode coupling are believed to relieve the issue.Coupling both ion modes and MTM in one simulation to investigate the thermal transport for both ions and electrons is planned as a future study.

FIG. 1 . 2 δ 2 ,
FIG.1.Diagram of the normalized distribution of absolute parallel magnetic vector potential, A ∥ , calculated from Eq.(18) using a shooting method.The horizontal axis denotes the distance with respect to the rational surface.The parameters are, β = 1%, ω * e /ν ei = 0.18e
FIG. 3. (a)Electron temperature equilibrium profile, (b)safety factor profile, and (c)ω * Te profile normalized to its toroidal mode number.
. Our simulation domain extends from the flux surface ψ N = 0.65 to ψ N = 1.05.To simplify the physics and focus on the impact of temperature BOUT-MTM gradients, we maintain a constant density of 10 20 [m −3 ] throughout the entire domain, thereby avoiding the influence of density gradients.
with ν being the local filed-line pitch.Therefore, e x is along the cross flux direction, and e y , e z are along the parallel and bi-normal directions.For linear simulation, the grid points are NX = 512, NY = 128, NZ = 16 in x, y, and z directions respectively, and for nonlinear simulation, NZ = 32.It is a twist-shift boundary condition in the y direction and periodic in the Z direction.For T e , A j∥ , the boundary condition is the Neumann condition in the inner boundary and the Dirichlet in the outer boundary; for U, it is Dirichlet in both boundaries.Other parameters are: Alfvén time τ A = 8.035 × 10 −7 [s], normalization length L = 3.1574[m], Alfvén speed V A = 3.9295 × 10 6 [m/s], normalization magnetic field B = 2.55[T], normalization electron temperature Te = 840[eV].
FIG. 4. Poloidal distributions of (a) normalized parallel magnetic vector ψ, (b) normalized electrostatic potential φ , and (c) normalized electron temperature perturbation T e for n = 8 MTM.Radial distributions of (d) ψ, (e) φ , and (f) T e at the inner mid-plane.ψ is normalized to L, φ is normalized to B LV A , and Te is normalized to Te .

FIG. 5 .
FIG. 5. Poicare plot in a poloidal slice presents Magnetic island chains in the linear simulation at t = 800τ A .The horizontal axis denotes the poloidal angle, where θ = π is the out mid-plane.The vertical axis denotes the normalized flux surface.The dashed red line indicates the location of q = m/n = 22/8 = 2.75 rational surface.
FIG. 6.(a) Linear growth rate γ and (b) mode frequency versus C ν for C T = 1.0.The blue solid is calculated by Eq.(38);The scatter plots are the simulation results, where the red solid dots denote the simulations without curvature and diamagnetic drift effects, the magenta stars denote the simulations with curvature effect included, and the green squares denote the simulations with both curvature and diamagnetic drift effects included.(c) Linear growth rate and (d) mode frequency versus C T for C ν = 1.0

in FIG. 6 (
FIG.7, RMS j ∥ distributions at outer mid-plan of different collisonality are presented.Following the definition of Yagyu and Numata 31 , the current layer width δ c is defined by the full width at half maximum height of j ∥ where the height is measured from the top to the first local minimal value.The current layer widths of C ν = 0.8, 1.0, 1.4 are δ c /ρ i = 3.22, 3.5, 3.63.MTM is localized around the bottom of the E r well where the local shear is nearly zeros.In such instances, the E r shear flow has minimal impact on the linear growth rate.However, as collisionality increases, the current layer width broadens.Consequently, the E r shear flow could suppress the MTM leading to a decrease in the linear growth rate. safety factor profiles align with those shown in FIG.(3).Toroidal modes are scanned from n = 5 to n = 11.In FIG.8(a), the locations of rational surfaces for different toroidal modes are labeled with various colored vertical lines.Correspondingly, FIG.8(b) displays the growth rate spectrum.

a
FIG. 8. (a) The diagram of ω * e profile and location of rational surfaces of different toroidal mode numbers; (b) the growth rate of different toroidal mode number cases.(c)-(e) presents the mode structure of different poloidal number modes corresponding to n = 8, 9, 11 cases.The dashed vertical red line denotes the location of ω * e peak and the the dashed vertical grey lines denote the location of rational surfaces.The solid black curve represents the safety factor profile.

FIG. 9 .
FIG. 9. (a) The time evolution of the root-mean-square (RMS) of magnetic perturbations is depicted.The solid blue and orange curves represent the evolution of flux surface-averaged magnetic vector potential and electrostatic potential perturbations at the q = 2.75 rational surface, respectively.The dashed curves indicate the perturbation evolution with equilibrium E r shear flow.(b) The electron temperature profiles at different time are illustrated.The dashed black curve represents the initial profile, the blue curve represents the profile after 3000τ A , and the red curve represents the profile after 3000τ A with equilibrium E r shear flow.
FIG. 10.The surface averaged thermal transport coefficient is normalized by the Bohm diffusion coefficient.(a)and (c) represent the temporal evolution of radial distributions of magnetic flutter component and E × B convection respectively.The horizontal dashed white line denotes the evolution of χ e at q = 2.75, ψ N = 0.85 rational surface.Correspondingly, they are plotted in solid lines in panels (b) and (d).

The normalized χ e of
FIG. 11.Poincare plots in the nonlinear simulation at t = 600τ A (top), t = 2500τ A (middle), and t = 8000τ A (bottom), corresponding to the simulation case with E r shear flow.The horizontal dashed red lines show the location of different rational surfaces.

to 4 . 43 ×
FIG. 13.At the q = 2.75 rational surface, (a) temporal evolution of total magnetic flutter electron thermal transport coefficient (blue), χ mf e , for α H = 0.05; (b) The dots in different colors denote the χ mf e calculated at each time-step of different α H cases within t = 2500 − 9500τ A .Black curves in different styles are calculated from Eq.(50).ρ * e is the electron gyro-radius normalized to minor radius; (c) Time averaged thermal transport coefficient (red dots) versus the free streaming parameter α H .
-driven transport being a potential contributor to this scaling, notably in the weakly collisional regime.Contrary to this, our simulations of collisional MTMs show a different trend where higher collisionality stabilizes MTMs, while lower collisionality induces unstable MTMs, leading to increased thermal turbulence transport and decreased confinement.The different results in comparison with the experiment are due to the collision regime.Either very low or high collision stabilizes MTMs.Our model limit the investigation in the strong collisional regime.However, in the opposite regime, MTMs are unstable as collisonality increases, leading to a worse confinement.This finding encourages us to further investigate MTMs in the weakly collisional regime in our future research.
. It is consistent with the discussion in Sec.V B. Simulations for more different α H are performed and the time averaged thermal transport coefficients χ e are presented in FIG.13(c).The maximum of χ e appears in the regime of α H ∈ [0.04, 0.05].A fitted scaling law as indicated by the blue curve in FIG.13(c) is, χ e = 113.2e−23.6αH − 113.3e −30.63αH D gB