Effects of resistivity and viscosity on dynamic evolution and radial position change of m/n = 3/1 double tearing mode

The effects of the plasma resistivity and viscosity on the dynamic evolution of the m/n = 3/1 double tearing mode (DTM) are studied and analyzed quantitatively using the CLT (Ci-Liu-Ti, which means magnetohydrodynamics in Chinese) code. In this work, we mainly focus on the change in the radial positions and the oscillatory dynamics of the magnetic islands grown on the two rational surfaces. We conduct a systematic investigation on the effect of viscosity on the DTM dynamics, which has rarely been studied before. From the results of the study, it is observed that the time required for entering the explosive phase decreases with decreasing viscosity. In the nonlinear phase, the kinetic energy exhibits an oscillatory behavior due to the magnetic flux injection and magnetic reconnection, and the oscillation amplitude is suppressed for a large viscosity due to dissipation. The effects of the plasma resistivity and viscosity on the change in the radial positions of magnetic islands are systematically explained. The change in the radial positions of magnetic islands occurs in an abrupt growth phase before the kinetic energy reaches its maximum value. Multiple position changes take place with a relatively higher reconnection rate and magnetic flux injection at low viscosity damping. A large range of radial vortices formed as a result of the change in the positions may have a positive effect on the transport.


Introduction
The reversed shear profile of the safety factor (q) is one of the attractive operational scenarios for attaining high performance and steady state operation in a tokamak [1]. A double tearing mode (DTM) is a common phenomenon in the reversed shear profile, which is one of the crucial subjects that should be solved for steady-state operation [1]. The development of the DTM instability can destroy the surface structure of the magnetic field, thus causing the degradation of the plasma confinement and can even lead to disruption [2]. The DTM instability occurs near the rational surfaces where the transport barrier is usually formed [3][4][5]. A steady state with the development of the DTM can form strong radial flows, which can prevent impurities from accumulating in the plasma core and play a role similar to the internal kink mode [6].
One of the important issues in DTM dynamics is the mechanism of nonlinear destabilization, which causes abrupt growth. Wang et al [7] first found that this abrupt or impulsive growth takes place in the early nonlinear stage in the external driving magnetic reconnection with a low resistivity. In the low resistivity regime, the released magnetic flux by magnetic reconnection is less than that by driven-in by the external driving flow in the linear growth phase, which leads to the accumulation of magnetic flux in the vicinity of the reconnection region. A sudden release of the pile-up magnetic flux causes an abrupt growth in the early nonlinear growth. In general, the abrupt growth weakly depends on the resistivity. Ishii et al [8][9][10] proposed a flux structure-driven nonlinear instability due to the triangularity of the magnetic islands around an X point in the cylindrical/toroidal geometry. Similar conclusions were also obtained by Janvier et al [11] in the slab geometry. Wang et al [12,13] found that the mechanism of the abrupt event in the slab reconnection is similar to that found by Wang et al [7].
The abrupt reconnection dynamics are accompanied by a rapid release of magnetic energy. A fast change in the magnetic configuration generally occurs as a result of rapid release of the magnetic energy. Therefore, mode development has been analyzed to understand the physical mechanism of the rapid growth. The phenomenon of the change in the radial positions of the magnetic islands grown on two rational surfaces has been observed from many simulations done so far in different geometries for the m/n = 3/1 DTM [8,9,[12][13][14][15][16][17][18][19] and m/n = 2/1 DTM [20] when other features were studied. For instance, in cylindrical geometry, Ishii et al [8] and Guo et al [18] studied the effect of the distance between the two resonant surfaces on the development of magnetic islands. Their results showed that the system could eventually saturate when all major islands have the same radial size for small distances between the two rational surfaces. The change in the radial positions of two islands occurs when the distance between the two rational surfaces is large in a certain regime of resistivity. Wei and Wang [19] showed the occurrence of the change in the radial positions in the full reconnection process driven by the strongly coupled island when they studied the mechanism causing the oscillatory behavior. Matsumoto et al [20] showed that the change in the radial positions of the islands occurs in collisionless magnetic reconnection of the m/n = 2/1 DTM. However, no systematic explanation for the effects of the resistivity and viscosity on the observed change was given, i.e. a detailed investigation of the conditions responsible for the change in the radial positions needs to be performed.
In this study, a three-dimensional, toroidal, and nonlinear compressible magnetohydrodynamic (MHD) simulation code, namely, CLT (Ci-Liu-Ti, which means magnetohydrodynamics in Chinese), has been used to revisit the dynamic developments of DTMs and an attempt has been made to address the issue of the characteristics and conditions including the change in the radial positions of the magnetic islands. A general rule in predicting the final state obtained from the DTM has been summarized. These results can have an important impact on the range of the magnetic topology transformation and radial flows.

Theoretical method
The simulations of the long-time nonlinear development of the DTM are carried out using the CLT code. The normalized compressible MHD equations used in the CLT code are as follows [6,[21][22][23][24] where ρ, v, p, J, B and E are the plasma density, velocity, plasma pressure, current density, magnetic field, and electric field, respectively. J 0 represents the initial current density and Γ(= 5/3) is the ratio of the specific heat of the plasma. All variables in CLT are normalized as follows: ) → E, and J/(B 00 /µ 0 a) → J where a is the minor radius, v A = B 00 / √ µ 0 ρ 00 is the Alfvén speed, and t A = a/v A is the Alfvén time. B 00 and ρ 00 are the initial magnetic field and plasma density along the magnetic axis, respectively. As in the previous simulations [6,20,21], the resistivity η, the diffusion coefficient D, the perpendicular and parallel thermal conductivity κ ⊥ and κ || , and the viscosity ν, are assumed to be constant in the simulation and are normalized as fol- respectively. In our simulation, the parameters are chosen to be η = 1 × 10 −7 ∼ 1 × 10 −4 , D = 1 × 10 −7 , κ ⊥ = 5 × 10 −6 , κ || = 5 × 10 −2 , and µ = 1 × 10 −7 -1 × 10 −4 , respectively. Because the pressure is not critical for a tearing mode instability, the plasma beta is assumed to β = 0 in our study.
In the reversal sheared configuration, the initial safety factor q-profile for the m/n = 3/1 DTM is as shown in figure 1. The formula used for the q-profile is as follows [25][26][27] where the parameters are set as q c = 2.2864, λ = 1.2, r 0 = 0.6034, δ = 0.4, and A = 0.4459. The distance between the two rational surfaces is 0.28a. In the present simulations, the aspect ratio is chosen to be R/r = 4/1. A uniform mesh with 256 × 32 × 256 (R, φ, Z) is used, and the convergent study is checked. The cylindrical coordinate system (R, φ, Z) is used in CLT, which can avoid the singularity near the r = 0 point that occurs in the toroidal coordinate (ϕ (r), θ, χ). We retain all toroidal and poloidal modes in order to the multi-helicity or mode-coupling MHD physics. The CLT code is accelerated with the graphics processing unit, which dramatically improves its computational ability [28].

Simulation results
The effect of resistivity on the m/n = 3/1 DTM is studied first. Figures 2(a) and (b) show the time evolution of the total kinetic energy E k , for different resistivities when µ = 1 × 10 −5 and µ = 5 × 10 −7 , respectively. In figure 2(a), for η = 1 × 10 −5 , the evolution of E k goes through three nonlinear stages: the linear growth phase, the explosive growth phase and the oscillation phase, as observed from the simulation results for slab geometry [12][13][14][15][29][30][31][32] and the experimental results from the tokamak fusion test reactor [1]. At the low resistivity, the three phases are distinguishable. At large resistivity, the explosive growth phase disappears and steady-state saturation is reached after linear growth, as is the case in figure 2(a) corresponding to η = 5 × 10 −5 . The growth rate of E k in the linear growth phase increases with increasing resistivity in the broader regime, as shown in figure 3(a), which does not follow a fixed scale [15,33]. This trend is similar to different viscosities (1 × 10 −7 -1 × 10 −5 ). At small resistivities, i.e. η ⩽ 1 × 10 −6 , the linear growth rate scales with the resistivity as γ ∼ η 0.6 , which is the same as that from a single resistive tearing instability predicted by theory and demonstrated by numerical calculations [8,15,[34][35][36][37]. This scale law of the linear growth rate is also observed with other slab and cylindrical geometries [11,27,29,[38][39][40][41][42][43]. Figure 3(b) shows the corresponding results as a function of the Prandtl number, Pr = µ/η. After the linear growth phase, the system enters the explosive growth phase, which results from the mutual driving force associated with the island expansion on the two rational surfaces. Since the speed of the island expansion increases with increasing resistivity, the time required for entering the explosive phase decreases, as shown in figure 2, which has also been observed in previous studies with other geometries [8,10,39]. The growth rate of E k in the linear growth phase decreases as the resistivity increases for η > 5 × 10 −5 , as observed in cylindrical plasmas [27]. This is because when the resistivity is large enough, the dissipation can be faster than the magnetic reconnection, resulting in the decrease of the growth of E k . The maximum value of E k reaches almost the same level during the fast reconnection phase, as observed with other geometries [8,10], indicating that the resistivity has little effect on the release of the maximum magnetic energy. Oscillations are observed in the oscillation phase. The oscillation of E k is suppressed to saturation at large resistivity, as shown in figure 2(a).
The effect of the viscosity on the m/n = 3/1 DTM is studied. Figures 4(a) and (b) show the time evolution of the total kinetic energy E k for different viscosities when η = 1 × 10 −5 and η = 1 × 10 −6 , respectively. Similar to the results in figure 2, the time evolution of E k goes through three developmental stages: the linear growth phase, the explosive growth phase, and the oscillation phase. At small viscosities, the two growth phases are distinguishable. At large viscosities, these two phases gradually merge into one, as in the case of µ = 5 × 10 −5 in figure 4(a). The growth rate in the linear growth phase decreases with increasing viscosity due to the dissipation effect: a large viscosity mainly slows down the inflows and outflows, and subsequently the process of magnetic reconnection is reduced [30,33] as seen in figure 5(a), which exhibits an opposite effect between the viscosity and the resistive since the resistivity in general has a destabilization effect on magnetic reconnection as shown in figures 2 and 3. It does not follow a fixed exponential scale law [33] but shows a strong dependence of the linear growth rate on  viscosity in the high-viscosity regime. The linear growth rate scales as γ ∼ µ −0.29 in the regime with high viscosity and low resistivity (µ ⩾ 5 × 10 −6 and η ⩽ 1 × 10 −6 ). This scaling law is very close to γ ∼ µ −1/ 3 and γ ∼ µ −2/ 7 predicted in the theory [39][40][41] and the simulation results with slab and cylindrical geometries [27,[29][30][31][39][40][41]. Figure 5(b) shows the corresponding results as a function of the Prandtl number Pr. The linear growth rate scales as γ ∼ Pr −0.29 for large enough values of Pr and is very close to γ ∼ Pr −1/4 [40] and γ ∼ Pr −1/3 [36] and deviates from γ ∼ Pr −1/6 [44]. For large values of µ, the system takes a long time to enter the explosive nonlinear phase because the large viscous dissipation reduces the growth of the resistive tearing mode [30], as observed with other geometries in previous studies [29,39]. As the viscosity increases, the maximum kinetic energy decreases because the conversion of kinetic energy into thermal energy increases. The kinetic energy exhibits an oscillatory behavior as a function of time and the oscillation amplitude strongly depends on the viscosity. The relative level of oscillations of E k increases with µ until µ = 1 × 10 −5 and then decreases. Such features are relatively difficult to observe from figure 4(b) since few cycles are shown. However, the difference between the first peak and the first valley in the oscillations can be used to show the features. The black solid line segments with the double arrows in figure 4(b) are equivalent to the difference between the first peak and the first valley in the oscillation for the case with µ = 1 × 10 −5 . It is obvious that this difference is larger than the difference corresponding to µ = 5 × 10 −5 and µ = 5 × 10 −6 . The reduction of the oscillation amplitude for large values of µ (µ ⩾ 1 × 10 −5 ) is due to higher dissipation [6,19].
The oscillatory phenomenon in our results is similar to that observed in studies of the single tearing mode [39,45], the m/n = 2/1 DTM [20], and the m/n = 3/1 DTM [19], but the triggering mechanism is different. Coelho [45] observed that a halt in the island growth and a subsequent periodic oscillation reconnection evolution characterized the marginally unstable m/n = 2/1 tearing mode in the case of a cylindrical  quasi-inviscid plasma. It was found that the inertial convection forces significantly affect the island growth rate. Takeda et al [39] suggested that the abrupt growth is due to the island coalescence and the following oscillations are due to the reconnection process influenced by the generation/inhibition dynamics of the external intermittent quadrupolar flow in the case of slab geometry and the highly unstable m = 1 tearing mode. Matsumoto et al [20] found that the exchange of islands takes place in collisionless magnetic reconnection (electron inertia) of the m/n = 2/1 DTM and the secondary reconnection occurs due to the current concentration using the gyrokinetic simulation in the case of cylindrical geometry. They are of a different nature from the observations in this work. Wei and Wang [19] suggested that the first explosive burst is driven by the strongly-coupled island, whereas the following intermittent bursts result from the regeneration of the reversed magnetic shear current profile due to the source term in Ohm's law for cylindrical geometry. Our results suggest that the explosive growth in the low resistivity is resulted from mutual driving force in unstable m/n = 3/1 DTM, and the subsequent smallamplitude oscillations are associated with the restoration of the magnetic flux surfaces owing to the source term in a toroidal geometry, as discussed later, similar to that for the m/n = 2/1 DTM in [6]. In fact, the oscillatory phenomenon, as well as the triggering mechanism in our results, are essentially identical to those presented by Wei and Wang in an earlier study [19], although there are differences in the geometries, source term, etc. It is the classical picture of a system oscillating around some equilibrium state and the system will finally reach a quasi-steady state in the presence of the dissipation (due to viscous and resistive processes).
The effect of resistivity on the mode development was analyzed, following the structure evolution of the magnetic field lines. Figure 6 shows the Poincaré plots of the magnetic field lines with µ = 1 × 10 −5 and η = 5 × 10 −5 at different times. In the linear growth stage, the inner and outer islands grow around each resonance surface at r 1 = 0.360 (at the poloidal angles of π/3, π and 5π/3) and r 2 = 0.632 (at the poloidal  angles of 0, 2π/3 and 4π/3) separately as the conventional tearing mode. As the islands grow bigger, the inner islands are squeezed outward, gradually deformed to a triangular shape, as shown in figure 6(b), and pushed to the X points of the outer islands ( figure 6(c)). With the development of the DTM, the magnetic flux located between the outer and inner islands is gradually reconnected, as seen in figures 6(a)-(c). Until the saturation of the DTM, there are some magnetic field lines un-reconnected (figure 6(d)), which leads to the inner islands being unable to reach the same radial position as the outer islands (r 1 < r 2 ). Thus, the magnetic structure remains in the state seen in figure 6(d) after t = 1500, which is similar to the nonlinear mode saturation for a higher resistivity with slab and cylindrical geometries [8,18].
For decreasing η, such as η = 1 × 10 −5 , the magnetic field structures in the early linear growth phase are very similar to those in the case with η = 5 × 10 −5 as shown in figures 6 and 7(a)-(d). The magnetic pressure generated by piling-up magnetic flux around the X points of the external islands (inner islands) causes the internal islands (external islands) to be pushed outward (inward). At t = 1780.5, the radial positions of the inner and outer islands are r 1 = 0.353 and r 2 = 0.5, respectively. As the DTM develops further, the inner island shrinks in the poloidal direction to push the entire inner island outward, whereas the outer island is pushed inward further. At t = 1869.6, all the magnetic flux between the inner and outer islands is reconnected, and the inner and outer islands reach the same radial position (r 1 = r 2 = 0.412), as shown in figure 7(e). With the further growth of the magnetic islands, the radial positions of the inner and outer islands are, respectively r 1 = 0.529 and r 2 = 0.353 at t = 2136.6, which suggests that the change in the radial positions of the magnetic islands on the two rational surfaces occurs, as seen in figure 7(h). Figure 8 shows the plasma flow pattern on the poloidal cross-section just before and after the change in the radial position of the magnetic islands, corresponding to figures 6(d) and 7(h), respectively. When the position of the magnetic island does not change, a vortex flow is formed in the narrow radial range of the q = 3 rational surfaces, which is just located on the magnetic island chain. The narrow vortex in the radial direction can block the radial transport and thus is conducive to the formation of transport barriers [3,4]. When the position of the magnetic island changes, a large range of vortices are formed in the radial direction, and a strong radial flow is formed locally. Thus, the exchange between core plasma and the plasma outside the q = 3 rational surface can be realized. This radial flow, due to the change in the position of the magnetic islands, degrades the plasma confinement. However, it can prevent impurities from accumulating in the plasma core since the radial flow carries both plasma and impurities from core to the outer region if the impurities are accumulated in the core [6]. It can be seen from figure 8 that the radial flow does not cross the flux surface at r = 4.65, which suggests that the transportation barrier could be formed at the outer boundary r = 4.65.
For a fixed viscosity of µ = 1 × 10 −5 , it is observed that there is no change in the radial positions of the inner and outer islands for η > 1 × 10 −5 , whereas a change is observed for η ⩽ 1 × 10 −5 . The maximum kinetic energy reaches almost the same value, indicating the same magnetic energy released. This can be explained by the results in figure 2, in which an abrupt energy growth is seen in the fast magnetic reconnection phase for η ⩽ 1 × 10 −5 , i.e. the peak becomes sharper by reducing the resistivity. It seems that the interaction of the inner and outer islands is not a sufficient condition for the change in the radial position of the DTM, whether there is a rapid energy explosion in a short time pays a key role in realizing the change in the radial position. In order to clarify the origin of the change in the radial position of the DTM, the evolution of the magnetic fields with a fixed resistivity η = 1 × 10 −5 and various viscosities was also investigated.
The effect of the viscosity on the mode development was analyzed. The time evolution of the mode development with decreasing viscosity was observed to be similar to that with decreasing resistivity, i.e. from no radial position change to radial position change. Even multiple position changes were observed to appear for the low viscosity of µ = 1 × 10 −6 . Figure 9 shows the Poincaré plots of the magnetic field lines with µ = 1 × 10 −6 and η = 1 × 10 −5 at different times. The topologies of the magnetic fields are quite similar to those in figure 7, which appear at different times. The dynamic process in figures 9(a)-(h) is also similar to that in figures 7(a)-(h). The islands grow bigger (figures 9(a)-(d)), and the inner islands are pushed outward as the outer islands are pushed inward (r 1 = 0.294 and r 2 = 0.471 in figure 9(d)). Subsequently, the inner islands shrink in the poloidal direction and are pushed outward, whereas the outer islands are pushed further inward (figures 9(e)-(g)). The inner and outer islands reach the same radial position (r 1 = r 2 = 0.382 in figure 9(e)). Subsequently, the inner and outer magnetic islands change their relative radial positions (r 1 = 0.529 and r 2 = 0.353 in figure 9(h)). The radial positions of the inner and outer islands change once in the early nonlinear stage, as shown in figures 9(d)-(h). In the following oscillation phase, multiple position changes are observed, as shown in figure 10. It eventually evolves toward a quasi-steady 3D equilibrium with magnetic islands after oscillations.
After the explosive reconnection, as shown in figure 4, residual flow is observed since the kinetic energy does not decrease to zero. However, the residual flow is not enough to trigger the following reconnection due to viscosity damping. The new magnetic flux is continuously pumped into the core region due to the source term in Ohm's law for recovering the magnetic flux to trigger another small oscillation. Therefore, the steady-state results from the balance between the newly injected poloidal magnetic flux in the core region and that damped in the reconnection process. The residual flow also indicates that the system never fully recovers during the following small oscillations, which can be seen from the Poincaré plots of the magnetic field lines (figures 9(d) and 10(a)).
For the broader viscosity range (µ = 1 × 10 −7 -1 × 10 −4 ) with a fixed resistivity η = 1 × 10 −5 , there is no change in the radial positions of the inner and outer islands when µ > 1 × 10 −5 , whereas change occurs when µ ⩽ 1 × 10 −5 , and  multiple changes occur when µ ⩽ 1 × 10 −6 , which suggests that the large kinetic energy can help the multiple position changes to occur. This is because the large plasma flow acts as the driving force to develop magnetic reconnection.
The multiple position changes of the inner and outer islands can be clarified by the injection of the magnetic flux, corresponding to the last term in equation (5), and viscosity damping. The magnetic flux injection causes the recovery of the system and viscosity damping causes loss of kinetic energy. Therefore, in order to achieve multiple position changes, a large injection of magnetic flux and low viscosity damping are necessary, which indicates that the multiple position changes take place in the regime with relative high resistivity and low viscosity, which is confirmed by the results of the systematic study shown in figure 11. This figure shows the dynamic behavior of the DTM in the (η, µ) plane. The radial position change is common for η ⩽ 1 × 10 −5 with a broader viscosity regime (µ = 1 × 10 −7 -1 × 10 −4 ). The multiple position changes only occur at relatively large resistivity (η = 1 × 10 −5 ) and low viscosity (η ⩽ 1 × 10 −6 ). Matsumoto et al [20] showed the position exchange of the islands in collisionless magnetic reconnection (electron inertia) using a gyrokinetic simulation considering the non-dissipation limit, which is corresponding to η → 0 and µ → 0 in figure 11.

Conclusion
The effects of the plasma resistivity and viscosity on the dynamic process of the DTM are studied and analyzed quantitatively using the CLT code. We mainly focus on the change in the radial positions and the oscillatory dynamics of the islands originally grown on two rational surfaces.
The evolution of the plasma kinetic energy E k goes through three stages: the linear growth phase, the explosive growth phase, and the oscillation phase. The two growth phases are almost indistinguishable for a large resistivity (η > 1 × 10 −5 ) and/or viscosity (µ > 1 × 10 −5 ). The growth rate in the linear growth phase increases with increasing resistivity and/or decreasing viscosity. Another feature in the linear growth phase is that the time required for entering the explosive phase decreases with increasing resistivity and/or decreasing viscosity. The time required for entering the explosive phase decreases with increasing resistivity because the speed of the island expansion increases in this case.
Since the plasma viscosity generally damps the inflow and outflow of plasma in the magnetic reconnection region, a small viscosity could result in the faster growth of E k and reduce the duration of the linear growth phase. The maximum value of E k is almost the same for different resistivities (with a fixed viscosity) and reaches a high level for low viscosity (with a fixed resistivity) during the growth phase because the main role of the viscosity is to convert the kinetic energy to thermal energy. In the oscillation phase, E k exhibits an oscillatory behavior in time, which is due to the presence of the recovery process resulting from the injection of magnetic flux and the release process associated with magnetic reconnection. The oscillation amplitude exhibits a dependence on η and/or µ. The oscillation of E k can fully be suppressed for a large η, and the oscillation amplitude decreases at a large µ (µ ⩾ 1 × 10 −5 ) because the kinetic energy is dissipated by the viscosity. This study is important because the origin of viscosity is related to the ambient turbulence, which is of significant interest in the fields of space and astrophysics.
A systematic explanation for the effects of plasma resistivity and viscosity on the change in the radial positions of the magnetic islands, originally located on the two rational surfaces, has been given for the first time. From the simulation results of this study, a general rule in predicting the final state of mode development of the DTM has been summarized. The E k characteristics corresponding to changes in the radial positions of the magnetic islands grown on the two rational surfaces is a sharp peak in the explosive growth phase, indicating a rapid energy release in a short time. The E k characteristics suggest that multiple position changes take place with a relatively high reconnection rate and magnetic flux injection at low viscosity damping. The simulation results indicate that the change in the radial positions of the magnetic islands occurs when η ⩽ 1 × 10 −5 and µ ⩽ 1 × 10 −5 , and multiple position change takes place when η ⩽ 1 × 10 −6 and µ = 1 × 10 −5 . A large range of radial vortices formed as a result of the changes in the position can have an important effect on the plasma and impurity transport.
The change in the radial positions of the magnetic islands that have been realized in the explosion stage of m/n = 3/1 DTM instability can lead to a broader spatial distribution of the vortex flow in the radial position. Thus, generating a strong radial flow can downgrade the energy confinement and prevent impurities from accumulating in the plasma core. However, the change in the relative radial positions of the inner and outer islands in the radial direction depends on the resistivity and viscosity. For the reversed shear operation with 2 < q < 3, whether it can provide a method to prevent impurity agglomeration remains to be further studied and discussed.