Nonlinear simulations are carried out for the microtearing mode using particle-based $\delta f$ gyrokinetic simulations for parameters relevant to spherical tokamaks. The present study finds that the microtearing mode can generate significant electron heat flux, which is predominantly carried out by the electromagnetic component of the heat flux with a negligible contribution from the electrostatic component. The mode sustains without the electrostatic component. We observe that the electron heat flux increases with the electron temperature gradient. The heat flux exhibits a very weak dependence on the collisions. It increases with electron *β* initially; however, at very high *β*, the electron heat flux is reduced.

## I. INTRODUCTION

In recent years, the research on the microtearing mode (MTM)^{1–32} in tokamaks has gained renewed interest as a source of heat loss via the electron channel. The advancement of experimental diagnostics and development of comprehensive simulation models aided by the availability of high-performance computation facility have made possible a thorough analysis of the MTM experimentally and numerically.^{19,23,24,29,31} This mode is driven by the electron temperature gradient in the presence of magnetic fluctuations and can produce strong heat flux in the electron channel.^{33–47} The overlapping of current layers on the mode rational surfaces and the electrons moving along the magnetic fields in the presence of the radial temperature gradient, thus, dumps a lot of energy out of the system.^{36,37,40} The mode structure for the electrostatic potential $\varphi $ exhibits odd parity, and the electromagnetic potential $A\u2225$ exhibits even parity. This feature is one of the many characteristics that is used to identify an MTM compared to other drift modes^{52} such as ion temperature gradient mode (ITG),^{53–55} trapped electron mode (TEM),^{56–58} short wavelength ion temperature gradient mode (SWITG),^{59–62} universal drift mode,^{63,64} electron temperature gradient (ETG) mode,^{65,66} and kinetic ballooning mode (KBM).^{49,78} However, one has to be careful while using parity as a tool to identify the MTM because, in the steep gradient region, in tokamaks ITG, ETG, TEM, etc., modes can also exhibit tearing parity.^{50} The occurrence of the MTM is now predicted in all variants of tokamaks, ranging from conventional tokamaks, such as DIII-D,^{11,47} Alcator-C,^{48} JET,^{23,24,29} ASDEX-U,^{5,7,12,19} and JIPPT-IIU,^{22} to spherical tokamaks, such as MAST^{1,2,4,14,15} and NSTX.^{3,8–10,13,16} The reverse field pinch machine such as RFX-mod^{6,17} can also give rise to the microtearing modes. The significance of this mode lies in the fact that while increasing *β* (ratio of kinetic to magnetic pressure, $2\mu 0nT/B2$) appears to debilitate some of the instabilities, e.g., ITG, it destabilizes the microtearing mode, like the KBM. Similarly, the collisions can reduce the growth of some drift modes like TEM by releasing trapped particles to passing particles. In contrast, microtearing modes are found to be destabilized by the presence of collisions under certain conditions. As the fusion community is seeking to build machines with higher and higher *β* and temperature in order to achieve fusion grade plasmas, the microtearing mode can still pose a challenge in achieving such a goal. The microtearing mode is particularly critical for the spherical tokamaks. Spherical tokamaks are another variant of the tokamak, which are compact in design, offer higher *β* by maintaining higher pressure at a given magnetic field, and are economically cheaper compared to the conventional one.^{67} The electron transport is observed to be comparatively stronger in these machines. Owing to high collisionality and high *β*, the microtearing mode might be a plausible candidate for the observed electron heat loss compared to other modes such as the electron temperature gradient mode or trapped electron mode. Several numerical studies relevant to the experimental parameters show that the MTMs can also survive even in the absence of collisions.^{6,17,20,21} Interestingly, MTMs are also predicted to exist before the onset of the edge localized modes^{19} and can exist even in the zero collisional limit.^{20,21} The electrostatic component of the perturbation can be stabilizing or destabilizing depending upon the value of safety factor *q.*^{7,12} The electrostatic component is shown to have little effect on MTM,^{15} whereas it is found to be strongly destabilizing in some other studies.^{4,51} Similarly, the equilibrium EXB shear appears to strongly suppress the turbulence,^{10} but it does not have a substantial effect in the parameter regime considered in Ref. 7. The theories of microtearing mode predict the velocity dependence of collision frequency for the microtearing mode to be unstable;^{36,38–42} however, gyrokinetic analysis shows that the velocity dependence might not be important for the instability in toroidal geometry.^{4,51} The magnetic drifts, which were missed in earlier theories, instead, are shown to play a critical role in causing the MTMs to be unstable.^{4,20,21}

The study of instabilities in the H mode pedestal and understanding of the trigger mechanism of edge localized mode (ELM) are necessary, as H-mode is considered to be a favorite operating regime in fusion devices. Of late, a lot of research has been carried out on the pedestal,^{23,24,29,31} which suggests that the KBM, though still important for pedestal physics, cannot be the sole candidate for producing edge transport. Nonlinear gyrokinetic calculations of the JET-ILW pedestal reveal that the experimental level of transport can be driven by MTM alone, where KBM is found to be in the second stability regime.^{23,24} Adding neoclassical and electron temperature gradient driven transport, it was shown that MTM can reproduce the experimental power balance across most of the pedestal. The use of a novel concept of “fingerprints” showed that microtearing modes along with ETG are more likely candidates for inducing anomalous transport in the H-mode pedestal.^{31} This fact is further supported by the magnetic fluctuation measurements in experiments^{68} (A greater detail on numerical and experimental investigations of MTM in the tokamak pedestal can be found in Ref. 31). More importantly, MTMs are found to be less susceptible to EXB shear suppression, which is strong in the pedestal region. Therefore, this renders MTM a more likely candidate that can exist in the pedestal regime and drive transport. The MTM instability is localized around a narrow region near the mode rational surface, and the fluctuations are mostly magnetic because of which, the effect of shearing is weaker.^{31,69} Thus, there is a growing consensus that the energy transport in the H-mode pedestal is due to the presence of robust MTMs, which can also be the precursor to edge localized modes. More research studies spanning different flavors of fusion devices and simulation methods are required to enhance the knowledge base. The present work aims to contribute to these continuing investigations of MTM in the H-mode pedestal by using a different tool of gyrokinetic simulations, namely, the PIC method, and for a different type of fusion machine, namely, the spherical tokamak. There have been many linear gyrokinetic studies of MTM on the pedestal of the spherical tokamak, for example, Refs. 15 and 16. However, the nonlinear properties of MTM in the pedestal region of spherical tokamaks are yet to be explored. At the same time, more and more PIC gyrokinetic codes are moving toward developing electromagnetic capabilities, but so far, there are not many nonlinear PIC simulations for the MTMs. This is partly due to the complexity pertaining to the kinetic electron physics and high *β* algorithm issues and partly due to the computational cost involved. Thus, this study is perhaps the first nonlinear PIC gyrokinetic simulations for pedestal parameters relevant to spherical tokamaks. Thus, the motivation of the present study is twofold: First, carry out nonlinear PIC simulation for microtearing modes, and second, contribute to the ongoing investigations^{23,24,29,31} on MTM in the H-mode pedestal by undertaking a case relevant to the spherical tokamak pedestal. Nonlinear properties of the MTM for the spherical tokamak pedestal are hitherto not well known. For this purpose, we use a particle-in-cell gyrokinetic code GEM^{70,71} that has been extensively used to simulate both electrostatic^{72–75} and electromagnetic^{76,77} modes including kinetic ballooning mode, kinetic peeling ballooning mode,^{78} and $n=1,\u2009m=2$ tearing mode.^{79,80} This code is extensively benchmarked with other gyrokinetic codes for different parameter regimes.^{81–84} We carry out a detailed investigation of the nonlinear electron heat flux of the microtearing mode and its dependence upon various factors such as the electron temperature gradient, electron *β*, collisions, and EXB shear for pedestal parameters of the spherical tokamak. These factors are considered to be critical for the microtearing mode both linearly and nonlinearly. In our study, we observe that the nonlinear electron heat flux increases with the electron temperature gradient. It increases with increasing electron *β* when *β* is low. However, at very high *β*, the electron heat flux appears to be reduced. We also observe that the nonlinear electron heat flux has a very weak dependence on the collision frequency. The equilibrium EXB shear has a limited effect on suppressing the electron heat flux. The mode survives and can generate substantial electron heat flux without the electrostatic perturbations. The interesting feature of the nonlinear simulations is that the major fraction of electron heat flux is contributed by the electromagnetic perturbations. This paper is organized as follows. Section II delineates a brief introduction of the simulation model GEM and presents the parameters used in the present study. Section III presents several convergence tests with respect to important numerical parameters. In Sec. IV, the important physics results are presented. Finally, the results are summarized in Sec. V.

## II. THE $\delta f$ GYROKINETIC SIMULATION MODEL

The present study is carried out using a $\delta f$ simulation code GEM.^{70,71} The GEM code is based on the Particle in Cell method where the 5D gyrokinetic equation is solved in the Lagrangian frame. In the PIC method, marker particles are used to sample, the phase space and the orbits of which are evolved in 5D phase space. The charge and current density that are required to solve the Poisson and Ampère equation are calculated from the marker particles.^{85} The fast parallel motion of the electrons imposes an upper limit on the time steps so that the Courant condition is satisfied. An adjustable split weight scheme is implemented in GEM, which allows us to use larger time steps. In this scheme, a portion of the adiabatic fraction of the perturbed distribution function is calculated separately. In the electromagnetic gyrokinetic formulation, the equation of motion requires calculating $\u2202A\u2225/\u2202t$, which incurs numerical instability with explicit schemes. In order to avoid calculating this term, the coordinate $v\u2225$ is replaced by introducing the canonical momentum $p\u2225=v\u2225+qmA\u2225$. Although this method allows one to avoid calculating $\u2202A\u2225/\u2202t$, the transformation from $v\u2225$ to $p\u2225$, however, introduces a large current term in Ampère's equation. This current term is not physical and should be canceled by the currents carried by the marker particles and weights in the simulation. Any inexact cancelation may lead to inaccuracy as well as numerical instability in certain cases. The cancelation of the large current term in Ampère's equation in GEM is achieved by suitably discretizing this current term in the same manner as the current from marker particles. This facilitates better cancelation of the large current term in Ampère's equation. GEM takes into account full gyrokinetic ions and drift kinetic electrons. A hybrid version with fluid electrons and gyrokinetic ions is also available and used to study low *n* modes and energetic particle physics. The physics modules include collisions, equilibrium flow, arbitrary shaped tokamak equilibrium, and impurity and energetic particle physics.

It is to be noted that the GEM code implements the field line following coordinates. The microinstabilities are typically aligned along the magnetic field in a tokamak with $k\u2225\u226ak\u22a5$. To take this advantage, the field line following coordinate system is used. The field line following coordinates are defined as follows:

where $(r,\u2009\theta ,\zeta )$ are the toroidal coordinates and *R*_{0} and *r*_{0} are the major radius at the magnetic axis and the minor radius at the center of the simulations box. Also, *q*_{0} is the safety factor at location *r*_{0} and $q\u0302(r,\theta )=B\xb7\u2207\zeta /B\xb7\u2207\theta $. This implies that $(x,\u2009y)$ coordinates label the field line, while the *z* coordinate is along the field line. Finally, the particle motion is described by $x\u0307=v\u2192G\xb7\u2207x,\u2009y\u0307=v\u2192G\xb7\u2207y$, and $z\u0307=v\u2192G\xb7\u2207z$. Here, $v\u2192G=v\u2225sb\u0303+vd+vE$ is the guiding center velocity and $vd$ and $vE$ are the magnetic and EXB drifts, respectively. The expression for total heat flux including both electrostatic and electromagnetic components is given by $Qtot=\u222bd3v\delta f12mv2(vE+v\u2225\delta B\u2225B)\xb7\u2207x$, while the heat flux with the electromagnetic component only is calculated as $QEMonly=\u222bd3v\delta f12mv2(v\u2225\delta B\u2225B)\xb7\u2207x$. The GEM code can be run for both local and global simulations. In the present study, we use the local (flux tube) version of GEM. The collisions in the simulations affect electrons only and are incorporated using the Monte Carlo technique. After each time step, a random change in the pitch angle is imparted. The average change in the pitch angle of the colliding particles is determined by two quantities: collision frequency and time step.

The parameters for the present study correspond to the NSTX-like spherical tokamak parameters.^{16} The physical parameters used for this study are given below: $\rho \u2217=0.00705$, aspect ratio 1.36, $r/a=0.93$, elongation $\kappa =1.64$ with $d\kappa /dr=0$, triangularity $\delta =0.40$ with $d\delta /dr=0$, safety factor $q0=7.94$, shear $s\u0302=r/q(dq/dr)=7.63,\u2009\tau =Te/Ti=1.5$, $R/LTi=\u2212R/Ti(dTi/dr)=4.83,\u2009R/Lne=\u2212R/ne(dne/dr)=0.92$, $R/LTe=\u2212R/Te(dTe/dr)=7.74$, electron beta $\beta e=2\mu 0neTe/B2=0.017$, collision frequency $\nu ei(cs/a)=5.58$, and EXB shearing rate $\gamma EXB(cs/a)=0.039$. To reduce the cost of the simulations, we do not consider any impurities.

## III. CONVERGENCE TESTS FOR NONLINEAR RUNS

In this section, we present convergence studies with respect to various numerical parameters. Note that these simulations are extremely expensive in terms of computational resources, and therefore, we restrict ourselves to only a few numbers of nonlinear simulations.

### A. Convergence with respect to the time step

Nonlinear simulation of microtearing mode requires very small time steps compared to the hydrogen ion gyroperiod here, with respect to GEM. This is partly because of the fast electron parallel motion as a result of using kinetic electron physics and partly because of the higher collisionality (compared to collisions in the core of a tokamak) for the parameters used in the simulations. For multimode simulations, for the high $k\u2225$ components, the condition $k\u2225ve\u25b3t\u226a1$ requires to be maintained, which demands very small time steps. Also, the average change in the pitch angle is proportional to the collision frequency times $\u25b3t$, which calls for a smaller $\u25b3t$, when collision frequency is very high, which is the case in the present study. In the following, Fig. 1 depicts the electron heat flux $Qe$ and ion heat flux $Qi$, normalized by $nuTuvu$ with respect to time for different time steps. Here, *n _{u}*,

*T*, and

_{u}*v*are the normalization constants for density, temperature, and velocity, respectively. We consider four different time steps as follows:

_{u}where the time is normalized with proton gyrofrequency Ω_{u}. It is clear from the figures that the simulations are unstable for time steps $\u25b3t\Omega u=0.75$ and greater. However, for time steps $\u25b3t\Omega u=0.50$ and $\u25b3t\Omega u=0.25$, the difference in the mean heat flux in the nonlinear steady state is only 10%. In the subsequent simulations, we consider the time step to be $\u25b3t\Omega u=0.5$, unless mentioned otherwise. It is clear from the figure that ion heat flux is about two orders of magnitude smaller compared to the electron heat flux and, therefore, negligible. Therefore, we shall focus only on the electron heat flux driven by the MTM in the present work.

### B. Convergence with respect to particle numbers

The number of marker particles in a cell plays a vital role in the particle simulations, as it also represents the velocity space resolution. Therefore, it is essential to test if the number of particles per cell is sufficient for robust simulation. In Fig. 2, we show the normalized electron heat flux for different numbers of particles per cell including both ions and electrons. The mean electron heat flux for the case with a total number of 256 marker particles per cell differs by around 12% when compared to the case with a total number of 512 marker particles per cell. Therefore, in all the subsequent simulations, we consider the number of marker particles per cell to be 256.

### C. Convergence with respect to the radial resolution

In Fig. 3, the electron heat fluxes are depicted with respect to time for different radial resolutions. The radial domain is fixed at $lx=0.3a$, where *a* is the minor radius, while the number of radial grid points is increased to $nx=64\u2009128$ and 256. For the latter two cases, the heat flux differs by around 20%. All the simulations presented in Sec. IV correspond to the resolution *n _{x}* = 256.

Based on the above discussed convergence studies, we reiterate here that all the simulations in Sec. IV consider $256\xd764\xd764$ grid points in the $x,\u2009y$, and *z* directions, respectively. The simulations consider a total number of 256 marker particles per cell. The time step considered in the simulations is $dt\Omega u=0.5$. The heat flux is normalized with $nuTuvu$, where *n _{u}*,

*T*, and $vu=Tu/mp$ are the normalizing values for temperature, density, and velocity and

_{u}*m*is the proton mass. Note that $nu=4.6\xd71019/m3,\u2009Tu=1.0\u2009keV$. The bar over

_{p}*Q*in Sec. IV represents the mean value with respect to the time and corresponds to the time average in the nonlinear steady state typically between the time window $t\Omega u=4000$ and the end of the simulations.

_{e}## IV. RESULTS AND DISCUSSION

In this section, we present and discuss the dependence of the microtearing mode on various equilibrium parameters, such as the electron temperature gradient, electron *β*, collisions, and equilibrium EXB shear. We also present a pure electromagnetic microtearing mode. The contribution of the electromagnetic and electrostatic components of the heat flux is also estimated.

### A. Snapshots of potentials in the nonlinear phase

Figure 4 presents the contour plots of the electromagnetic potential of the MTM in the nonlinear phase. The left panel shows the snapshot for $A\u2225$ in the nonlinear phase in the *xy* plane, and the right panel shows the snapshot for $A\u2225$ in the nonlinear phase in the *xz* plane. Note that $A\u2225$ is normalized by $Tu/evu$, where *e* is the electron charge. Similarly, Fig. 5 presents the electrostatic potential $\varphi $ in the *xy* and *xz* planes in the nonlinear phase in the left panel and in the right panel, respectively. Note that $\varphi $ is normalized by $Tu/e$.

### B. Dependence of the MTM on the electron temperature gradient

The free energy source for the microtearing mode is the electron temperature gradient, which gives it a drift mode character. In the presence of electromagnetic perturbation, electrons moving parallel to the magnetic field line also displace in the radial direction. When there exists an electron temperature gradient across the magnetic field, it imparts a thermal force parallel to the magnetic field. This thermal force along the field increases the current perturbation and consequently enhances the electromagnetic perturbation, thus giving rise to MTM instability.^{4–10} This thermal force is proportional to the electron temperature gradient. Therefore, one expects that the instability growth rate and electron heat loss strongly depend upon the electron temperature gradient. However, it is to be noted that these theories are based on simple geometries and without the effect of magnetic drifts. Recent theories^{25,32} indicate that the magnetic drifts and magnetic shear curvature effects can be important for the microtearing mode. In Fig. 6, we plot the mean nonlinear electron heat flux with respect to $R/LTe$, which is a measure of the strength of the electron temperature gradient as $R/LTe=\u2212R\u2009dlnTe/dx$. Therefore, increasing $R/LTe$ indicates an increasing electron temperature gradient or steeper electron temperature profile. It is clear from the figure that with the increasing electron temperature gradient, the electron heat flux increases. This is consistent with previous numerical simulations.^{7,10,12,13,20}

### C. Dependence on electron *β*

The MTM becomes unstable only when *β*, which is the ratio of the kinetic pressure to the magnetic pressure, is above a critical value. Once this threshold value is overcome, the growth rate and consequent electron heat flux of the mode increase with a further increase in *β*. In order to study the effect of *β* on the nonlinear electron heat flux, we carry out several nonlinear simulations with the increasing value of electron *β*. The mean electron heat flux for different values of electron *β* is shown in Fig. 7(a). Note that *β _{e}* is varied consistent with the equilibrium. It is apparent from the figure that the electron heat flux increases with electron

*β*for lower values of

*β*. The increase in the

*β*value increases the $A\u2225$ fluctuation through Ampère's equation, which, in turn, increases the electron heat flux. However, it is observed that the electron heat flux does not increase in a monotonic manner with electron

*β*or electromagnetic perturbation despite MTM being an electromagnetic mode. With an increase in

*β*, the electron heat flux increases initially, but, at higher values of

_{e}*β*, the electron heat flux starts to decrease. An explanation for this behavior can be offered following previous theoretical studies. For higher toroidal mode numbers, the tearing parameter $\u25b3\u2032$ is negative. Therefore, for MTM to become unstable, one requires a finite amount of $\u2207Te$. However, the MTM is the result of the competition between the destabilization produced by the electron temperature gradient and stabilization produced by the magnetic energy, which increases with field line bending.

_{e}^{37,38,40}Therefore, at very high

*β*, for a given $\u2207Te$, the MTM is weakened due to the stabilizing effect from the field line bending, and hence, the electron heat flux is reduced. In order to investigate whether this is a nonlinear effect or not, we have calculated the linear growth rates for the same equilibrium parameters and plot these in Fig. 7(b) as a function of

*β*. We observe the same trend in the linear regime as well, which shows that the growth rates of the MTM initially increase with increasing

_{e}*β*for smaller values and then growth rates decrease with a further increase in

_{e}*β*. However, this trend might not be observed in the tokamak experiments given the fact that such a roll over occurs at a very high

_{e}*β*and tokamaks may not be operating at such high values of

_{e}*β*. It is to be noted that this observation reminds the second stability regime observed in the case of KBM with increasing

*β*. While such a transition to the second stability regime for the KBM can occur at both low and high magnetic shears, the magnetic shear considered in the present case is high. A more detailed study is required in this regard.

Also, it would be interesting to understand the effect of higher *β* and roll over in heat flux using advanced analytical theories.^{25,28,32} Such an analytical study and comparison with simulation results will be carried out in a future work using simple geometry and parameters.

### D. Effect of equilibrium EXB shearing

The MTM is found to be less susceptible to the EXB shearing^{31,69} owing to the fact that they are localized to a narrow region and the turbulence for MTM is mostly from magnetic fluctuations. This makes MTMs have more robust instability in the pedestal region where the EXB effect is strong. In this section, we investigate the impact of the EXB shearing rate on the MTM turbulence. The results are shown in Fig. 8(a). The mean heat flux is plotted with respect to the EXB shearing rate *γ _{EXB}* normalized by $cs/a$, where

*c*is the ion sound speed and

_{s}*a*is the minor radius. It is clear from the figure that the EXB shearing rate does not have a strong effect on the electron heat flux. In order to further analyze and understand this observation, we carry out linear simulations for the same equilibrium parameters. We calculate growth rates of the MTM corresponding to different toroidal mode numbers and plot the same in Fig. 8(b). In the same plot, we also show the value of the baseline EXB shearing rate. It is evident from the figure that the linear growth rates of the modes are much stronger than the shearing rate, that is, $\gamma \u226b\gamma EXB$. This perhaps explains the observation of the weak effect of EXB shearing on the nonlinear heat flux level.

### E. Dependence on collisions

Previous theories^{36,38–42} underlined the importance of collisions and their velocity dependence for the MTM to be unstable. Therefore, it is usually expected that the MTM grows with increasing collisions, leading to higher electron heat fluxes. While this is observed in some of the gyrokinetic studies^{4,7,10,12} showing the nonmonotonic dependence of the MTM on collisions, there are also reports where the MTM is observed to be independent of collisions^{15,51} and exists even when there are no collisions at all.^{6,20,21} More specifically, the linear MTM is observed to be weakly dependent on collision frequency for pedestal parameters relevant to the spherical tokamak.^{15,51} Therefore, it is interesting to explore the effect of collisions on the MTM turbulence and the resultant nonlinear heat flux. In view of this, we carry out nonlinear simulations for the microtearing mode by varying the collision frequency and calculating the mean electron heat flux. This is shown in Fig. 9, where normalized electron heat fluxes are plotted against the collision frequency *ν _{e}* in units of $cs/a$, where

*c*is the ion sound speed and

_{s}*a*is the minor radius. Note that the collisions here are electron-ion collisions only and they are modeled as pitch angle scattering. Interestingly, we observe that the nonlinear electron heat flux is virtually independent of collision frequency. This is perhaps the first nonlinear simulation of microtearing mode, which shows the weak dependence of the electron heat flux on collisions, although there are several linear studies

^{6,15,20,21,51}showing the independence of the MTM growth rates with respect to the collision frequency. The independence of the electron heat flux with respect to the collision frequency may be explained on the basis that in toroidal device, the role of magnetic drift and trapped particles becomes more important

^{6,15,20,32}for MTM than the collision frequency.

### F. Importance of the electrostatic potential $\varphi $

In the present study, we remove the electrostatic $\varphi $ component of the perturbation and carry out a nonlinear simulation retaining only $A\u2225$ perturbation in order to evaluate the role of electrostatic potential. This is shown in Fig. 10(a) plotted with the simulation that retains both $\varphi $ and $A\u2225$ perturbations. It is clear from the figure that the electrostatic component increases the electron heat flux slightly in the nonlinear steady state. For the case with $\varphi $, however, the simulation quickly grows in the linear phase and enters the nonlinear phase faster than the case without $\varphi $. This is in conformity with the linear results shown in Fig. 10(b) where the growth rates corresponding to modes with different toroidal mode numbers are plotted with and without including $\varphi $. It is clear from the linear analysis that the linear growth rates are found to increase in the presence of the electrostatic perturbations and also consistent with the analytical predictions.^{32,38} Thus, it is interesting to note that although $\varphi $ might have a stronger destabilizing effect linearly, in the nonlinear phase, it contributes little to the total heat flux. The slight reduction in the heat flux without $\varphi $ is consistent with Fig. 11 in Sec. IV G, where the red line is with the total heat flux comprising both EM and electrostatic components and the blue curve is without the electrostatic component. Compared to heat flux from the EM component only in Fig. 11, we found that the heat flux with the pure EM case of Fig. 10 is slightly lower, which means that the reduction in the heat flux in the latter case might partly be coming from the lower magnetic potential.

### G. Comparison of relative contribution of electromagnetic and electrostatic perturbation to the total electron heat flux

In the present study, we estimate the contribution of the electromagnetic component of the electron heat flux to the overall electron heat flux. For this purpose, we plot the total electron heat flux (red curve) along with the electromagnetic component (blue curve) of the electron heat flux in Fig. 11. It is clear from the figure that the electromagnetic heat flux is almost equal to the total electron heat flux. This means that the majority of total electron heat flux comes from the electromagnetic component of the heat flux with little contribution from the electrostatic component. This observation is in conformity with other nonlinear studies^{10,13,24} also, where the significant portion of the electron heat flux is shown to be contributed by the electromagnetic component.

## V. SUMMARY

In this paper, we present detailed linear and nonlinear simulation studies of the microtearing mode for parameters relevant to the edge of the spherical tokamak using a PIC based gyrokinetic code GEM. Some interesting nonlinear characteristics of the microtearing mode are observed. The ion heat flux is negligibly small for the nonlinear microtearing mode compared to the electron heat flux, which is a signature of the MTM.^{10,24} The electron heat flux increases with the increase in the electron temperature gradient. The source of free energy for the MTM is the electron temperature gradient. Therefore, the instability becomes stronger with the increasing electron temperature gradient, which, in turn, enhances the nonlinear turbulence and the resultant electron heat flux. We observe that the electron heat flux changes nonmonotonically with respect to the increase in the electron *β*. The electron heat flux increases with increasing *β* initially for lower values of *β*; however, at very high *β*, the heat flux is found to be reduced. Although a clear understanding is lacking at this moment, however, simpler analytical theories suggest that the field line bending effect might be stronger at higher *β*. A more detailed investigation is required to confirm this. Starting with simple geometries and parameters and comparison with simpler theories might shed light on this effect. This will be pursued in a future work. We find that for the parameters considered, the microtearing mode generated electron heat flux is virtually independent of collision frequency. As suggested in earlier numerical simulations and theoretical studies, curved magnetic shear, magnetic drift, trapped particles, etc., which are dominant effects in the edge of the tokamak, might be more important for the mode than the collisions.^{4,15,20,25,32} The electrostatic potential is found to increase the electron heat flux only slightly when compared to MTM simulations without including $\varphi $ perturbations in the simulation. We also observe that the electron heat flux is mostly contributed by the electromagnetic component, implying that the electrostatic component has very little contribution to the total electron heat flux. This is a characteristic of the microtearing mode observed also in earlier numerical simulations.^{10,24} We finally conclude that a more detailed comparison with analytical theories would be required to understand the physics mechanism starting from a simpler geometry and parameters. Such a study will be carried out in a future work.

## ACKNOWLEDGMENTS

This material is based upon the work supported by United States Department of Energy Grant Nos. DE-FG02-08ER54954 and DE-SC000801. It uses the resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science user facility. We sincerely acknowledge fruitful discussions with Dr. J. Canik and Dr. W. Guttenfelder.