Theoretical and global simulation analysis of collisional microtearing modes Open Access

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,¹ JET,² and others, MTMs are frequently observed. MTM can dominate electron thermal transport in the JET ITER-like wall pedestal under certain conditions,³ plays a role in determining the behavior of electron temperature during an ELM cycle,⁴ and in double barrier operations, it generates a substantial electron thermal transport at the Inner Transport Barrier (ITB).⁵ 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.⁸ first proposed a linear kinetic conductivity model for this type of mode. This kinetic theory was further developed by Drake et al.,⁹ 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¹² 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,¹⁵ 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 Sec. 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 Sec. IV as we explore the mode 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 Sec. 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, Sec. VI summarizes our findings.

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.

In summary, for local assumption around the rational surface, Larakers's conductivity model reduces to Drake's model in both $ω ≪ ν e i$ and $ω ≫ ν e i$ limits, and reduces to Hassam's model in $ω ≪ ν e i$ limit. Larakers's model thus bridges the two limits and also includes $k ∥$ effects.

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, $Δ ′ = [ A ∥ ′ ( 0 + ) − A ∥ ′ ( 0 − ) ] / A ∥ ( 0 )$⁠. The $Δ ′$ is solved in the outer region. Approximately $A ∥ ̃ ∼ e − k y | x |$⁠, and $Δ ′ = − 2 k y$⁠.¹⁹ Then the solution of Eq. (18) can be obtained from an integration within the current layer width, Δ_j.¹⁸ 

In this section, we introduce the control equations, equilibrium profiles, and parameters used in this study.

If the diamagnetic drift effect is considered (i.e., finite $κ$⁠), then the second term in the right-hand side of vorticity Eq. (24) must be included to maintain equilibrium.²⁶ 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, $V E 0 = b 0 × ∇ Φ 0 B 0$⁠, is included to balance equilibrium ion diamagnetic drift for particle conservation.²⁷ Also, E_r shear flow introduces the Doppler shift effect, $∂ ∂ t → ∂ ∂ t + V E 0 · ∇$⁠. Consequently, when referring to the equilibrium ion diamagnetic drift effect, we consider the equilibrium flow induced by equilibrium ion diamagnetic drift, balanced by the E_r shear flow-essentially investigating the impact of E_r shear flow.

In our simulation, we utilize a shifted circular geometry equilibrium generated by CORSICA.²⁸ The equilibrium profile is displayed in Fig. 3. 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 gradients, we maintain a constant density of $10 20 m − 3$ throughout the domain, thereby avoiding the influence of density gradients.

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.²⁹ 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++,³⁰ where $x = Ψ − Ψ 0 , y = θ , z = ζ − ∫ θ 0 θ ν ( Ψ , θ ) d θ$⁠, with ν being the local field-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 , and 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$⁠, and normalization electron temperature $T ¯ e = 840 eV$⁠.

Electrostatic potential $ϕ ̃$ and normalized parallel magnetic vector potential $ψ ̃$ are solved from Eqs. (36) and (37) by the invert-Laplace Laplacian equation solver of BOUT++. This solver is introduced in detail in the  Appendix. In our simulations, we exclude the current driven term by default, $B 0 × ∇ ψ ̃ · ∇ J ∥ 0$⁠, to focus on the temperature gradient-driven MTM.

A linear simulation of n = 8 is performed. Figure 4 displays the poloidal and radial distributions of $ψ ̃ , ϕ ̃$ and $T ̃ e$ 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.

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(b), it is evident that collisionality 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, this finding is a result of our theoretical model, which is only valid in the strongly collisional regime. In general, MTMs are destabilized by collisionality in the presence of a temperature gradient. In the weak collision regime, the growth rate is proportional to collisionality, indicating that collisionality destabilizes MTM,⁵ 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.

The reduction in growth rate is attributed to the non-local shear of E × B flow. As illustrated in Fig. 7, RMS $j ∥$ distributions at outer mid-plane of different collisionality are presented. Following the definition of Yagyu and Numata,³¹ 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.

MTM is unstable only when the rational surface is in good alignment with the peak of the electron diamagnetic frequency, $ω * e$⁠, profile. As pointed out by Larakers and Curie,³² MTM is unstable with $Δ s / x * < 0.3$⁠, where Δ_s is the distance that one rational surface shifts away from the peak of the location of $ω * e$ and $x *$ is the characteristic width of $ω * e$ profile. Once the rational surface is too far away from the peak $ω * e$⁠, MTM cannot survive. The competition between the nonalignment damping effect and gradient-driven effect determines the unstable regime of MTM. The nonalignment is a stabilizing effect. Our conclusion is consistent with Larakers's global calculation,¹¹ Curie's SLiM model,^32,33 and also Hassan's pedestal simulation.³⁴ 

A linear simulation is conducted to confirm this prediction. The equilibrium temperature and safety factor profiles align with those shown in Fig. 3. Toroidal modes are scanned from n = 5 to 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. 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. A higher toroidal mode number leads to a more pronounced $ω * e$ profile since $ω * e ∝ n$⁠, and then increases the drive. The mode structures of the most unstable mode, n = 8, 9, 11, in Fig. 8(b) are presented in Figs. 8(c)–8(e). Our global simulations show that the perturbations peak at the location of peak $ω * e$ but not always at the rational surface. For n = 9, 11 modes, the most unstable m modes offset the rational surface and the symmetry of the mode structure breaks down. As predicted by the local theory, the n = 9 mode is expected to have a larger growth rate than that of n = 8. However, n = 9 has a smaller growth rate despite having a larger $ω * e$⁠. This is due to the non-local nonalignment effect, the linear growth rate of the n = 9 mode decreases.

As the blue dotted line in Fig. 8(b) shows, $J 0 ∥$ destabilizes the nonaligned modes, like n = 9 mode, but has little influence on the aligned n = 8 mode. The curvature effect can stabilize the nonaligned mode as the magenta dotted line shows in Fig. 8(b).

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.

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 $T ̃ e$ 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 equilibrium E_r shear flow, the temperature gradient exhibits a more significant decrease.

To illustrate the MTM thermal transport, we use the case with E_r shear flow as an example. The normalized χ_e of the two components are shown in Fig. 10. The normalized factor $D g B = T e ρ i * / 16 e B 0 = 0.0071 [ m 2 · s − 1 ]$ represents the gyro-Bohm diffusion coefficient, and $ρ i *$ is the ion gyroradius normalized to minor radius. Figures 10(a) and 10(c) present the temporal evolution of the distributions of the two components of χ_e. In comparison with the magnetic flutter component, the electron thermal transport by E × B convection is negligible. This is attributed to the changes in magnetic field topology induced by MTM. It can be illustrated in Fig. 11. In the nonlinear phase, such as $t = 2500 τ A$⁠, the original closed flux surfaces break down due to the strong magnetic 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 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$⁠.

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.