Effect of the curvature and the $\beta$ parameter on the nonlinear dynamics of a drift tearing magnetic island

To characterize how the pressure gradient affects the evolution of a magnetic island, linear and nonlinear self consistent numerical simulations have been performed. A grid number of $n_{x}=128$ is chosen for the radial direction and $n_{y}=128$ for the poloidal direction (equivalent to 48 modes in this direction, including dealiasing). The energy conservation relation derived from Eqs. (1–3) is

where $E_{m}=0.5<\!\!\left|\nabla(\psi-\psi_{0})\right|^{2}\!\!>$, $E_{p}=0.5<\!\!p^{2}\!\!>\!\!/C^{2}$ and $E_{k}=0.5<\!\!\left|\nabla\phi\right|^{2}\!\!>$ are respectively the magnetic energy, the pressure energy and the kinetic energy of the fluctuations. The brackets $<.>$ mean here an average over the simulation domain. $S$ is the source term linked to the curvature and the pressure gradient, proportional to the radial pressure flux , $S=-\alpha_{S}<p\partial_{y}\phi>$ with $\alpha_{S}=\frac{v^{\star}}{C^{2}}(1-\kappa_{2})+\kappa_{1}>0$ because $\kappa_{2}<1$. Note that a local flattening of the pressure by radial exchange of pressure cells, gives a fluctuation $\delta S<0$. Moreover, the interchange source term $S$ is not modified by the generation of zonal flow. Figure (3) shows the time evolution of $E_{m}$, $E_{p}$ and $E_{k}$ for the parameters chosen in this work with $\beta=10^{-3}$ as well as the corresponding $E_{m}^{CT}$ and $E_{k}^{CT}$, for a classical tearing mode (i.e $p=0$ and $\kappa_{i}=0$). In comparison with the evolution of a classical tearing mode, four regimes are observed in the nonlinear simulations of a magnetic island in the presence of the interchange term. First, there is a linear regime where the magnetic island is formed. Second, the system reaches a quasi-plateau phase. Then, a transition occurs and as it will be shown later this is linked to the interchange parameters. Finally, the system reaches a new saturated state.

During the first two phases, the evolution of the energies is not strongly affected by the presence of the curvature terms. The evolution of the magnetic island follows closely the time trajectory of an island driven by a tearing instability. However, at $t_{\star}=13200\;\tau_{A}$, a transition occurs. Figures (4) and (5) show snapshots of the fields $\psi$, $\phi$ and $p$ respectively before and after the transition. The two dimensional profiles of the pressure and the electrostatic potential (represented through isocontours) are strongly affected by this transition. After this phase, the structure of the mode changes and a flattening of the pressure is obtained.

In order to understand the origin of the transition and to characterize the structure of the electrostatic potential after the transition, we first assess the importance of small scales. The spectra before the transition (a) and during the transition (b) are shown in Figure (6). An equipartition between the energy of the magnetic flux and the energy of the pressure is observed at large scales $0.8<k_{y}<5$ whereas there is an equipartition between the energy of the pressure and the kinetic energy at small scales $8>k_{y}>14$. We observe that these properties continue to persist when the transition occurs, the energy level of the modes $k_{y}>1.2$ being roughly unaltered. We also note that as apparent in Figure (6) only the large scales are affected at the transition. A detailed analysis shows that the transition occurs when the kinetic energy of the mode $k_{y}=0$ becomes equal to the one of $k_{y}=0.4$ (mode 1). This suggests that the dynamics of the structure of the electrostatic potential is quasi-linear and that the transition occurs when the mode $k_{y}=0$ becomes energetically dominant. To delineate the quasi-linear nature of the magnetic island dynamics, we have performed a simulation with only four poloidal modes and with the same parameters ($\beta=10^{-3}$). Figure (7) shows the time evolution of the energies for this simulation. The comparison with Figure (3) demontrates that one needs only four modes to describe the time evolution of the energies. Therefore the dynamics of the system is the result of quasi linear effects. Further, the time evolutions of the kinetic energies of modes 0, 1 and 2 are presented in the right hand side panel of Figure (7) and compared with the evolution of modes in a classical tearing run, i.e without pressure effects ($p=0$ and $\kappa_{i}=0$). In the latter case, the growth is driven by the mode $1$, the transition does not occur and the mode $0$ is not generated. Conversely, in the case where the pressure effects are included, the transition occurs and the mode $k_{y}=0$, i.e zonal flow, is strongly generated. It is also the first to be amplified exponentially, at the beginning of the transition. This suggests that the transition is linked to a strong amplification of the zonal flow. Nevertheless, when we perform the same run, suppressing artificially the mode 0, i.e the zonal flow, we find that the transition still occurs, roughly at the same time, but with a weaker amplitude.

This suggests that althouh zonal flows play a predominant role they are not the sole factor responsible for the destabilizing mechanism. Indeed we can show from the analysis of the snapshots of $\phi$, $p$ and $\psi$ given by Figure (8) that the mode 0 plays an important role in the triggering mechanism of the transition phase. We can observe that during the transition, the pressure cells are crossing the resonant surface at the current sheet in both directions, producing at the end a modification of the mode structure. Clearly the potential structure suggests that the crossing results from an advection by the flow. Let us stress that during this phase, the classical tearing picture of an incoming flow from the sheet into the island is no more valid. Between $t=15000\tau_{A}$ and $t=17000\tau_{A}$, $S(t)$ decreases because the reorganisation of the cell is radially equivalent to an exchange of pressure cells with the gradient of the pressure fluctuations being outward. From an energetic point of view, see eq. (9), $\kappa_{1}$ and $\kappa_{2}$ have a negligible effect because the dominant contribution in the interchange source term is linked to $\omega_{\star}/C^{2}$. In Figure (7), the results of a simulation with four modes and without the curvature terms are presented ($\kappa_{1}=\kappa_{2}=0$). We find that the transition does not occur. At least, it does not occur at $2.5t_{\star}$, showing that a more complex mechanism due to curvature terms might be at play. We next investigate the origin of the zonal flow.

We investigate the origin of the zonal flow that is generated whenever $\beta\neq 0$ by considering separately the energy transfer from the tearing mode to the zonal flow and from also the other modes, in particular, the small scales.

Following [8], the equation for the flow energy feeding the mode $k_{m}=m2\pi/L_{y}$ can be written as,

where

Here $\phi_{m}=\hat{\phi}_{m}e^{iky}+\hat{\phi}_{-m}e^{-iky}$, where $\hat{\phi}_{m}$ is the $m$ Fourier component of $\phi$.

In order to understand the origin of the generation of the zonal flow, it is useful to project equation (10) on the mode $m=0$:

The curvature term does not directly feed the zonal flow. In Figure (9), the time evolutions of $T_{0}^{R}$, $T_{0}^{M}$, $T_{0}^{LB}$ and $T_{0}^{KI}$ for the simulation with $\beta=10^{-3}$ are presented. The contributions of the line bending and the kink terms are very weak. However, at the transition, there is a strong generation of the Reynolds stress and the Maxwell stress contributions. To proceed further, it is useful to separate the contributions from zonal flow ($m=0$), the mode $m=1$ and other modes for each of the transfer functions. Let us introduce $\phi_{>1}=\sum_{m>1}\phi_{m}$. We can define three contributions in each transfer function. For instance, in the case of $T^{R}_{m}$, we have:

where $T^{R}_{0m}=\int dx\phi_{m}([\phi_{0},\omega_{0}])_{m}$, $T^{R}_{1m}=\int dx\phi_{m}([\phi_{1},\omega_{1}])_{m}$, and $T^{R}_{>1m}=\int dx\phi_{m}([\phi_{>1},\omega_{>1}])_{m}$. Clearly, by definition $T^{R}_{11}=0$ and $T^{R}_{0m}=T^{M}_{0m}=0$. Let us focus on the energy transfer to zonal flow, neglecting the weak contributions of the line bending and the kink terms. Equation (11) then becomes:

Using the above prescription we have checked that the main contibution to the Reynolds and Maxwell stresses comes from the mode $m=1$ while the contributions of the small scales are weak.

The left panel of Figure (10) presents the transition time and the time where the first quasi plateau saturation occurs for various values of $\beta$. We note that the time corresponding to the first quasi-plateau phase does not depend on the $\beta$ parameter. This is in agreement with the results of Figure (3) which shows the time evolution of the energies for two cases namely the classical tearing case ($\beta=0$) and the $\beta=10^{-3}$ case. However the saturation time depends strongly on the $\beta$ parameter. For a regime where the pressure effects are strong, i.e for a high value of $\beta$, the transition occurs quickly whereas, for low $\beta$ regimes, the transition occurs later. Further, in the right panel of Figure (10), the effect of the interchange parameter $\kappa_{2}$ on the transition is shown for a simulation with four modes and with $\beta=10^{-3}$. The transtion time depends clearly on the $\kappa_{2}$ parameter. It tends to stabilize the first plateau phase.

To summarize, the nonlinear transition results from quasilinear effects. Zonal flow amplification at the transition is due to the energy transfer from the mode $m=1$ to the mode $m=0$ through mainly the Maxwell stress. The curvature term $\kappa_{2}$ linked to the interchange effect does not directly feed the growth of the zonal flow. However, we have shown that this term controls the transition time $t_{\star}$. The transition leads to an effective radial exchange of pressure cells generating an outward mean pressure gradient of fluctuations. The shape of the pressure structure after the transition implies that a diamagnetic velocity $\tilde{\omega}_{\star}$ has been nonlinearly generated, driving a rotation of the island. This driving is of course in competition with the zonal flow. Let us in the next section analyze quantitatively the island poloidal rotation.

In Figure (11), the time evolution of the poloidal position is presented for a nonlinear simulation with $\beta=0.025$. The dynamics of the energies for a simulation with $\beta=0.025$ has the same behaviour as the one obtained in Figure (3) with $\beta=10^{-3}$. However, as shown in the left panel of Figure (10), with such a high value of $\beta$, the transition occurs earlier around $t=5600\tau_{A}$. In Figure (11) a comparison with the island position obtained from the model Eq.(15) is also shown. The derived model is in agreement with the numerical data and the dynamics of the island rotation is recovered. We would like to mention here that the time integration for these results has been performed on a very long time scale compared with the Alfvén time.

Equation (15) shows that the effective island frequency is the result of different contributions. In order to investigate the effect of the interchange terms on the island rotation, Figure (12) presents the evolution of each frequency for a simulation with $\beta=0.025$ setting $\kappa_{1}=\kappa_{2}=0$ (left panel) and for a simulation with $\beta=0.025$ including the curvature/interchange terms (right panel). First, for the two simulations, the contributions to the rotation of the equilibrium magnetic field $L_{\psi_{0}}$ and of the resistivity $L_{\eta}$ are weak. Moreover, Figure (12) shows that the frequency dynamics is not affected by the curvature terms during the linear regime and during the first quasi plateau phase. Actually, during the linear formation of the magnetic island, the rotation is controlled mainly by the linear diamagnetic drift while $\tilde{\omega}^{\star}$ and $\tilde{\omega}_{E\times B}$ are weak. During the first quasi plateau phase, the nonlinear diamagnetic drift and the $E\times B$ poloidal flow are strongly generated, and affect the island rotation. The $E\times B$ poloidal flow is the most important contribution to the frequency during this regime. However, after the first quasi plateau phase, interchange terms affect the dynamics of the frequency. When the interchange terms are switched off, the time evolution of the frequencies is in agreement with previous results found for a drift tearing mode [4]. After the nonlinear generation of the flows, linear and nonlinear diamagnetic drifts cancel each other. As a result, the $E\times B$ poloidal flow controls the effective frequency of the island, $\omega_{t}\sim\tilde{\omega}_{E\times B}$. The right panel of Figure (12) shows that when interchange terms are included, such canceling of the total diamagnetic frequency does not occur anymore after the transition. Hence the total diamagnetic drift then provides the main contribution to the island rotation. However, after the transition, clearly $\partial\omega_{t}/\partial t\sim\partial\omega_{E\times B}/\partial t$. The asymptotic island velocity is enhanced by the curvature terms $\kappa_{1}$ and $\kappa_{2}$.

Nonlinearly, the $E\times B$ flow is generated and affects the rotation of the island. In order to investigate its origin, the flow equation (Eq.1) is projected on the mode $k_{1}$ for the limiting case of $\nabla^{2}_{\perp}\approx\partial^{2}_{x}$. We obtain :

We have defined the $E\times B$ poloidal flow as $\tilde{\omega}_{E\times B}=k_{1}\partial_{x}\phi_{0}$. So multiplying Eq. (16) by $k_{1}$ and averaging over the current sheet $\delta$, we obtain :

where

with $R\left(t\right)$ being the Reynolds stress, $M\left(t\right)$ being the Maxwell stress and $V\left(t\right)$ being the viscosity contribution to the $E\times B$ flow.

On the left panel of the Figure (13), the time evolutions of R, M and V are presented for a simulation with $\beta=0.025$. As expected, the stresses are nonlinearly generated at the beginning of the first quasi plateau phase allowing the growth of the $E\times B$ poloidal flow. Except at the end of the linear regime where the Reynolds stress is not yet generated, the viscosity term is very weak and does not play an important role in the generation of the flow. The most important contributions come from the Reynolds and Maxwell stresses which balance each other. There is a strong amplification of the amplitude of the stresses during the transition. On the right panel of Figure (13), a closeup of the temporal dynamics of R, M and V during the transition time evolutions are presented. Note that during this transition, whereas the viscosity term is still weak, the evolutions of the Reynolds and Maxwell stresses are complementary. The two stresses tend to balance each other during this phase where the flow is crossing the separatrices, limiting the level of the generated zonal flow, even if they are both growing in amplitude. Once the transition has occured the amplitude of the mean nonlinear brackets in the vicinity of the separatrix, $R$ and $M$, fall. However, in the new dynamical equilibrium the resulting $E\times B$ poloidal flow persists asymptotically and is driven by the Maxwell stress, as a response of the magnetic structure to the new field distribution.

In order to investigate the origin of the nonlinear diamagnetic drift, we follow the same procedure for the pressure equation (Eq.2). After projection on the mode $k_{1}$, we obtain:

where

$dC$ is the contribution to the nonlinear diamagnetic drift of the divergence of the convective term, $dM$ is the contribution of the divergence of the Maxwell stress and $D$ is the contribution of the diffusivity.

On the left panel of Figure (14), the time evolutions of dC, dM and D are presented. The three contributions to the nonlinear diamagnetic drift are generated at the beginning of the nonlinear regime. During the first quasi plateau phase, dC, dM and D participate actively in the generation and growth of the nonlinear diamagnetic drift $\tilde{\omega}^{\star}$. At the transition there is an amplification of the amplitude of the three contributions. The right panel of Figure (14) presents a closeup of the time evolutions of dC, dM and D during the transition. During this transition, the contribution of the diffusivity does not grow and is relatively weak compared to dC and dM. Note that both the divergence of the convective term and of the Maxwell stress feed the nonlinear diamagnetic drift. We observe that in the first phase of the transition, the nonlinear diamagnetic drift is driven by $dC$, the term linked to the advection of the pressure cells. In the following phase this term is balanced by the Maxwell divergence stress leading to a stabilization of the island dynamics. This shows the importance of the coupling parameter $C$ during the transition. In the saturation phase, the origin of $\tilde{\omega}^{\star}$ comes mainly from the divergence of the convective term, dM and D becoming relatively weak.

In Figure (15), the effect of the $\beta$ parameter on the rotation frequency is presented. The frequencies have been time averaged from $t=15000\tau_{A}$ to $t=30000\tau_{A}$. At low $\beta$, pressure effects are weak and the situation is close the classical tearing situation. The competition between Reynolds and Maxwell stresses produces nonlinearly neither the zonal flow nor a diamagnetic drift. The asymptotic rotation frequency $\omega_{t}$ increases with $\beta$. The nonlinear diamagnetic frequency increases almost linearly with $\beta$, but with a slope lower than the linear one, allowing a global asymptotic drift of the island in the electron diamagnetic direction. Let us note however that the direction of the island rotation depends on the value of the viscosity parameter [4]. For any $\beta$, the zonal flow contribution to the island drift is weaker than the diamagnetic one. Note also that from $\beta\sim 0.015$, $\omega_{E\times B}$ decreases, and this is linked to the transition observed in Figure (16) where the effect of the $\beta$ parameter on the average stresses is shown (from $t=15000\tau_{A}$ to $t=30000\tau_{A}$). The amplification of the stresses leads to the nonlinear generation of the mean flows and hence affects considerably the rotation of the magnetic island. More precisely, it appears that for low $\beta$ values, the amplitudes of the stresses are very weak as in the classical tearing case. This explains why nonlinear diamagnetic drift and nonlinear $E\times B$ poloidal flow do not affect strongly the rotation of the island in those cases. However, for high $\beta$ regimes, the amplitude of the stresses, in particular the Maxwell and convective contribution become more important. At $\beta\sim 0.015$, a transition is observed : first, $dM$ ceases to be neglectable compared to $dC$, second the Reynolds contribution $R$ starts to grow, weakening the global $E\times B$ flow, as observed in Figure (15).