Impact of micro-scale stochastic Zonal Flows on the macro-scale V-RMHD modes

TL;DR

This paper presents a stochastic model of micro-scale Zonal Flows and investigates their stabilizing effects on macro-scale MHD modes, revealing potential turbulence-driven stabilization mechanisms.

Contribution

Introduces a stochastic modeling approach for micro-scale ZFs and examines their impact on stabilizing MHD modes, a novel insight into turbulence effects on plasma stability.

Findings

Stochastic ZFs can significantly stabilize (2,1) and (1,1) MHD modes.

Both passive and semi-stochastic models show stabilization effects.

Micro-scale turbulence influences macro-scale plasma stability.

Abstract

A model is developed to simulate micro-scale turbulence driven ZFs, and their impact on the MHD tearing and kink modes is examined. The model is based on a stochastic representation of the micro-scale ZFs with a given Alfv\'{e}n Mach number, M_S. Two approaches were explored: i) passive stochastic model where the ZFs amplitudes are independent of the MHD mode amplitude, and ii) the semi-stochastic model where the amplitudes of the ZFs have a dependence on the amplitude of the MHD mode itself. The results show that the stochastic ZFs can significantly stabilise the (2,1) and (1,1) MHD modes even at very low kinematic viscosity, Pr, where the mode is linearly unstable. Our results, therefore indicate a possible mechanism for stabilisation of the MHD modes via small-scale perturbations in poloidal flow, simulating the turbulence driven ZFs.

Equations

$$\frac{\partial \tilde{W}}{\partial t}+{\mathbf{v}}_0\cdot \nabla \tilde{W}+v_A{\nabla }_{\parallel }{\rho }_s^2{\nabla }_{\perp }^2\bar{\psi }$$

$$=v_A{\rho }_s\frac{1}{r}\frac{\partial \tilde{\psi }}{\partial \theta }\frac{4\pi {\rho }_s}{c{B}_0}{j}_0^{{}^{\prime }}+\frac{v_{th}{\rho }_s}{r}\{\tilde{\psi },{\rho }_s^2{\nabla }_{\perp }^2\tilde{\psi }\}$$

$$+\frac{v_{th}{\rho }_s}{r}\{\tilde{W},\tilde{\varphi }\}+\frac{{\rho }_s^2}{r}\frac{d{W}_S}{dr}\frac{\partial \tilde{\varphi }}{\partial \theta }+\nu {\nabla }_{\perp }^2\tilde{W},$$

$$\frac{\partial \tilde{\psi }}{\partial t}+{\mathbf{v}}_0\cdot \nabla \tilde{\psi }+v_A{\nabla }_{\parallel }\tilde{\varphi }=\frac{v_{th}{\rho }_s}{r}\{\tilde{\psi },\tilde{\varphi }\}+\frac{c^2\eta }{4\pi }{\nabla }_{\perp }^2\tilde{\psi },$$

$$\tilde{W}={\rho }_s^2{\nabla }_{\perp }^2\tilde{\varphi }$$

$$W_S=\nabla \times (\delta v_s{\mathbf{e}}_{\theta })=\frac{1}{r}\frac{d}{dr}(r\delta v_s)$$

Full text

I Introduction

The development of large scale (i.e. machine size) Magneto-Hydro-Dynamic (MHD) instabilities such as $(m=1/n=1)$ internal kink mode Marshall1973 ; VonGoeler1974 ; Ksdomtsev1975 ; Thyagaraja1991a ; Porcelli1996 ; Monticello1986 ; Mendonca2018 ; Kai-Qi2018 and/or $(m=2/n=1)$, Tearing Modes (TM) of which a special case for fusion plasmas is the Neoclassical-Tearing Mode Wilson1996 ; Haye2002 ; Smolyakov1993 ; Hegna1998 ; ZChang1995 ; Carrera1986 ; McDevitt2006 ; Chandra2015 are important limiting factors for establishing a stable confined plasma for long enough time to reach the fusion conditions.

The $(1,1)$ kink mode arises inside the $q=1$ rational surface when $q$ at the axis is $<1$, where $q$ is the safety factor. This mode is observed to trigger the sawtooth oscillations in the plasma core, influencing the confinement quality. Most notable impact of the sawtooth oscillation, is a slow rise of the order of the resistive time scale, in the plasma temperature in the centre, followed by an abrupt crash (on a 100 microsecond time scale), relaxing the temperature profiles back to the initial values at the beginning of each cycle VonGoeler1974 ; Larche2017 . The loss of temperature in the central plasma region limits the achievable pressure profile. However, the sawtooth crash benefits the plasma core confinement by removing the unwanted impurities from the core. Control of the sawtooth crashes e.g. their amplitude and frequency, hence would be important for optimising plasma performance in a fusion device Campbell1988 ; McClements1996 ; Larche2017 ; Yang2012 ; Mao2001 .

The classical $(2,1)$ TM mode can result in the destruction of the topology of the nested magnetic flux surfaces, where the reconnection of the field lines produces magnetic islands in the confined region. The continuous growth of the islands then destroys magnetic confinement and can produce runaway electrons, and lead to major disruptions. A disruption is an event that induces a rapid loss of thermal energy followed by a quench of the plasma current. Disruptions release large heat loads and electro-magnetic forces on the surrounding machine structures, and can result in a significant damage to the Plasma Facing Components (PFCs).

When the classical $(2,1)$ TM is linearly stable (when $\Delta^{\prime}<0$), it was shown both theoretically and experimentally Wilson1996 ; Hegna1998 ; DeVries1997 , that Neoclassical bootstrap currents (or their absence in the interior of a magnetic island) can lead to a non-linear destabilisation of the island resulting in NTMs. This is a non-linear instability which requires an initial seed island of sufficient threshold amplitude.

Given the important consequences of the MHD instabilities (both linear and non-linear) which generate global deleterious effects, understanding their dynamics and the development of control actuators to avoid them is crucial for current and future fusion machines DeVries2009 ; DeVries2014 ; Lehnen2013 .

Another limiting factor for achieving fusion is the confinement time. This is limited by the existence of a zoo of small-scale electro-magnetic turbulence (due to the ion/electron drift waves) which are ubiquitous in magnetically confined fusion plasmas, and are considered to be responsible for anomalous transport. Although the small-scale turbulence does not result directly in a disruption, it can modify/enhance the large scale MHD modes through non-linear mode coupling and/or the generation of sheared Zonal Flows (ZFs). In a recent review by Ishizawa and co-authors Ishizawa2019 , the findings of an extensive body of work that has been devoted to the study of the multi-scale interactions and the interplay between turbulence and large scale TMs are examined. We will not try to re-examine these works here and will only refer to the a summary: the inverse energy cascade via non-linear parity mixing of the small-scale ballooning modes (such as ITG, ETG, KBM with twisting parity) Ishizawa2015 ; Coppi1977 ; Kadomtsev1971 ; Dorland2000 ; Pueschel2008 or electron Micro-TMs (with tearing parity) Doerk2011 ; Guttenfolder2012 ; moradi2013 , is believed to result in the destabilisation of otherwise stable TMs (e.g. reduce the stability threshold as compared to purely MHD analysis) or further growth of the unstable TMs (i.e. enhancement of the growth rate of the island). Thus, turbulence can create large-scale magnetic islands by the merging of small-scale islands. Other mechanisms such as anomalous current drive by turbulence i.e. much like a small-scale dynamo effect Li2012 , or turbulence-driven polarisation current as the result of turbulence-generated ZFs following the reconnected field lines of the island Ishizawa2013 , have also been suggested as drivers of magnetic seed islands.

The back reaction of the magnetic islands on the turbulence is expressed through a direct energy cascade from MHD modes to turbulence as a result of the non-linear mode coupling observed in EAST tokamak Sun2018 . Moreover the magnetic perturbation of the field lines, especially around the X-point of the island increases the turbulence outside the island DeVries2014 , while the coherent vortex flows generated by the interactions between the turbulence and MHD mode at the boundary can act as a transport barrier that leads to suppression of the drift-wave instabilities inside the island Li2009 ; Ishizawa2013 ; poli2009 ; harriri2015 ; hornsby2015 ; Bardoczi2017 . The enhancement of the turbulence outside, and its reduction inside the island, along with the parallel transport results in increased flux which flattens (not completely) the temperature profile inside the separatrix of the island. This flattening of the temperature profile, changes the resistivity inside the island and further reduces the perturbed bootstrap current which has a destabilising effect on the island. It increases the growth of the island and eventually the degradation of the confinement outside the island equilibrates the profiles inside and outside leading to its saturation at lower plasma stored energy Bardoczi2016 ; Bardoczi2017 .

It is a common view that in fusion plasmas turbulence is regulated by its self generating sheared $E\times B$ flows (ZFs) Diamond2005 ; Biskamp1983 ; Diamond1984 . The mechanism is explained as generation of the sheared flows via non-linear interactions producing Reynolds-Stresses, which in turn result in shearing of the turbulent eddies, thus reducing their radial correlation length Ham1995 ; Ham2015 . However, it is argued that in high $\beta$ plasma regimes, where $\beta$ is the ratio of plasma kinetic energy to magnetic energy, the ITGs responsible for generating the ZFs are stabilised due to magnetic field bending at high $\beta$, and KBM, TEM or ETG modes arise which do not produce strong ZFs. These micro-scale modes can grow beyond a physically relevant saturation level, unless through mode coupling they transfer their energy to stable modes Ishizawa2015 or form coherent structures such as magnetic islands Ishizawa2019 .

Given the important role of the small-scale turbulence on the stability regime of the large scale MHD modes and their regulation by ZFs, in a more general approach, here we examine the impact of a stochastic ZF following a given probability distribution function (PDF) on the large-scale (1,1) and (2,1) modes. Our aim in this paper is to consider the possibility of dynamic stabilisation of the MHD modes by the stochastic ZFs. We note that it is well-known theoretically and experimentally Thyagaraja1991 ; Thyagaraja1993 ; Yang2012 ; Mao2001 , that in principle, one can dynamically stabilise MHD modes using external applied electro-magnetic perturbations with specific frequencies and amplitudes. We further extend this idea by including a spectrum of random perturbations on small-scales and show that in principle, their associated ZFs can result in stabilisation of the linear and non-linear unstable (1,1) and (2,1) visco-resistive MHD modes.

Our results are based on a series of numerical simulations with the visco-resistive MHD version of the CUTIE code Mendonca2018 ; Chandra2015 where a stochastic ZFs model has been included in the vorticity equation. CUTIE is a non-linear, global, electromagnetic, quasi-neutral, visco-resistive MHD code in a large aspect-ratio cylindrical geometry (neglecting linear toroidal coupling due to curvature effects while including non-linear mode interactions).

Two models were considered: i) passive stochastic ZFs model where the amplitudes of the ZFs are independent of the amplitude of the MHD mode; ii) semi-stochastic coupling model where the amplitudes of the ZFs are proportional to the amplitude of the mode itself. The first model aims to examine the impact of a passive stochastic flow arising from the small-scale turbulence on the MHD stability. The second model aims to represent, in a general manner, the interplay between the MHD mode and the turbulence generated stochastic flows following a direct-cascade phenomenology where the energy is transferred from the large scales to the small-scales which back-react with the mode itself. Our results show that above a certain critical amplitude, the stochastic ZFs have a strong stabilising effect on the (1,1) and (2,1) modes, even at a very low kinematic viscosity (e.g. $Pr=1$) where the modes are linearly unstable.

The paper is arranged as follows: In section II we describe the details of our simple model and the CUTIE code. In section III we present the results of the stochastic passive model or the (2,1) mode. In section IV the results for both (1,1) and (2,1) modes are shown after applying the semi-stochastic ZFs model. Section V contains a summary of our findings and our conclusions.

II Model description

The CUTIE code employes a visco-resistive MHD model of plasma in a periodic cylindrical geometry $(r,\theta,z)$, where $r$ is the radial, $\theta$ is the azimuthal, and $z$ is the axial coordinate Chandra2015 ; Mendonca2018 . The equations for the mean and fluctuating parts of the electro-magnetic fields are split using a Fourier expansion where the mean fields are represented by the $(0,0)$ components, and the remaining terms represent the fluctuating components of the fields. The mean fields are co-evolved with the equations for fluctuating quantities. The fluctuation equations for the vorticity $\tilde{W}$ and poloidal flux function $\tilde{\psi}$ are:

[TABLE]

where vorticity is defined as:

[TABLE]

and the vorticity due to the stochastic ZFs is defined as

[TABLE]

In the present problem $T_{0i}=T_{0e}=0.275$keV, and $n_{0}=4\times 10^{14}$cm*-3* are considered uniform and constant, with $T_{0i}$ and $T_{0e}$ being the ion and electron temperatures, and $n_{0}$ is the mean density, respectively. Here, $\rho=r/a$ following the Mendonca2018 , is assumed as $a=10$cm and $R_{0}=100$cm; and in this large aspect ratio limit the linear curvature effects are neglected. The resistivity $\eta$ and viscosity $\nu$ are specified quantities and are kept constant during the simulations. Additionally, $\rho_{s}=\frac{v_{th}}{\omega_{ci}}$, where $v_{th}^{2}=(T_{0i}+T_{0e})/m_{i}$, $\omega_{ci}=(eB_{0}/m_{i}c)$, $e$ is the elementary charge and $m_{i}$ is the ion mass. $v_{A}=B_{0}/(4\pi m_{i}n_{0})^{1/2}=2.18\times 10^{8}$cm/s is the Alfvén velocity, and $\mathbf{v}_{0}$ is the sub-Alfvénic equilibrium sheared flow (both axial and poloidal contributions) which is set to zero in the current study. The equilibrium current density is defined as $j_{0}=-\frac{c}{4\pi}\frac{1}{r}\frac{d}{dr}(r\frac{d\psi_{0}}{dr})$, where $\psi_{0}$ is the poloidal equilibrium flux function.

Alfvén time is set as $\tau_{A}=a/v_{A}$, and resistive and viscous time are defined as: $\tau_{\eta}=\frac{4\pi a^{2}}{c^{2}\eta}$, and $\tau_{\nu}=\frac{a^{2}}{\nu}$. Prandtl number is $Pr=\frac{\tau_{\eta}}{\tau_{\nu}}$, and the Lundquist number is $S=\frac{\tau_{\eta}}{\tau_{A}}$, and the corresponding Reynolds number is defined as: $Re=\frac{\tau_{\nu}}{\tau_{A}}$. $B_{0}=2T$, and $S=10^{6}$ are assumed. The profile of the safety-factor, $q$ ($=rB_{z}/RB_{0\theta}$ where $B_{z}=B_{0}$, and $B_{0\theta}=-d\psi_{0}/dr$), is prescribed as function of $\rho$: $q=q_{0}[1+(\rho/\rho_{0})^{2\lambda}]^{\lambda}$ with $[q_{0}=0.9,\lambda=1,\rho_{0}=0.6]$, and $[q_{0}=1.4,\lambda=2,\rho_{0}=0.74]$ for the (1,1) and (2,1) cases, respectively Mendonca2018 ; Chandra2015 . Since we prescribe both $B_{0}$, $S$ and $q$, the profile of $\eta$ is obtained from the condition that $\eta j_{0}$ is constant.

In the simulations presented we have used 101 radial grid points, 17 poloidal and 9 toroidal Fourier modes. The time step used in the simulations was $\Delta t=2\times 10^{-10}s$, thus, $v_{A}\Delta t/\Delta r=0.44$ is the radial Courant number.

In the case of a passive stochastic ZFs model, the ZF velocity $\delta v_{S}$ is defined as follows: $\delta v_{S}=M_{S}X(\rho,t)v_{A}$ where $M_{S}$ is specified to be the stochastic ZFs Alfvén Mach number, and $X(\rho,t)$ is the stochastic variable is assumed as a white gaussian noise, with zero mean and variance of 1. Here, $X(\rho,t)$ is updated at each radial grid point and at time intervals specified by $\delta t=10^{4}\Delta t$. The time interval is selected as 2$$\mus, which represents the time scale of the micro-turbulence. Therefore, the stochastic variable $X$ is uncorrelated for times larger than $\delta t$, and radial distances greater than $\delta\rho$.

With this prescription, the vorticity Eq. II is no longer a deterministic evolution equation and is a Langevin type stochastic differential equation. In the following we present the results of numerical simulations solving the Eqs. II and II.

III Results for passive stochastic ZFs model

The simulations are performed by setting the initial condition, i.e. $q$-profile, to allow for the $(2,1)$ mode to be linearly unstable. In the absence of the stochastic ZFs i.e. $\delta v_{S}=0$, by adjusting the $Pr$ as high as $75$, a non-linearly saturated state can be reached, as seen in Fig. 1 where the time evolution of the $|\tilde{\psi}|_{max}$ is shown (see red solid line).

However, when the stochastic ZFs are included, the $(2,1)$ mode can be stabilised even at significantly lower Prandtl numbers. For example taking $M_{S}=5\times 10^{-3}$ even at $Pr=10$ the mode reaches a saturated state at a lower amplitude, see the black dotted line in Fig. 1.

The relationship between system stability through $Pr$ and $M_{S}$ (always keeping Lundquist number $S$ fixed) shows some interesting features. For $Pr<2$, to stabilise the (2,1) mode the value of $M_{S}$ has to increase. For example, at very low $Pr=1.25$, the stochastic ZFs Mach number, $M_{S}$, has to be increased to $10^{-2}$ in order to avoid the blow up of the mode, see dashed-dotted blue line in Fig. 1. The amplitude of the (2,1) mode in this case, is significantly reduced, and it saturates at a finite value.

In a range of $Pr>4$ and for the same value of $M_{S}$, the mode saturates at lower amplitudes for lower values of $Pr$, see for example green dashed dotted line vs black dotted line in Fig. 1, corresponding to $Pr=4$ and $Pr=10$, respectively. We note that the simulations show that as the $Pr$ is lowered, while keeping the $M_{S}$ fixed, the saturation amplitude is lower. Whilst it would be reasonable to expect that the saturation level to be higher at the lower values of viscosity (i.e. lower $Pr$), the opposite is observed.

A possible explanation for this behaviour is that at high values of $Pr$, the turbulence driven small scales will be strongly dampened by the dissipative effect of the viscosity since the damping effect of viscosity scales as $|k|^{2}$. Therefore the impact of small scale driven ZFs are diminished at these high values of $Pr$. At lower values of $2<Pr<4$, this damping effect is reduced and ZFs are more efficient at stabilising the (2,1) mode.

Figures 2 (a-d) show the corresponding profiles of the $dB_{r}/B_{0}$, $\tilde{W}$, equilibrium current density $j_{0}$, and safety factor $q$, for the case with $Pr=1.25$. The values are compared between $t=8.9$ms and $t=11.3$ms, without and with the stochastic ZFs effects, respectively. As can be seen here, stochastic ZFs can result in a significant stabilisation of the $(2,1)$ mode, comparing the red solid line to the blue dash-dotted line. In the pre ZFs state, $j_{0}$ shows local flattening around $q=2$ resonant surface, which is changed by the ZFs almost back to the linearly unstable original state. Thus, in effect, the ZFs are actually stabilising the (2,1) linearly unstable mode. The vorticity profile in this case shows peaks with higher amplitudes than in the absence of the stochastic ZFs model at $\rho\sim 0.6-0.8$, corresponding to the radial location of the $(2,1)$ mode (see blue dashed line in Fig. 2 (b)). This indicates that the ZFs promote a direct cascade in the vorticity distribution of the (2,1) mode. The free energy of the equilibrium $j_{0}$, is transferred through the ZFs to the high $k$ vorticity of the mode, where viscosity is able to dissipate it, thus resulting in the lowering of the amplitude of the magnetic fluctuations.

Figure 3 illustrates the contour plots (top) and the corresponding spectra (bottom) of the current density fluctuations at the selected times. The poloidal structure of the $(2,1)$ mode remains visible at $t=11.3$ms. However, the amplitude is reduced significantly. This can also be seen in the spectra plots where a strong reduction of the current density fluctuation amplitude is obtained.

To determine the stability dependence of the $(2,1)$ on the amplitude of the ZFs represented by the Mach number, $M_{S}$, a sensitivity scan in $(Pr,M_{S})$-space was performed, while holding all other parameters fixed. The results are shown in Fig. 4. The stability boundary indicates that for $2<Pr$ the minimum ZFs Mach number that is required to stabilise the (2,1) mode is $M_{S}\sim 5\times 10^{-3}$, whilst for $1<Pr<2$ this level has to increase, and for example for a $Pr\sim 1$ it has to be as high as $M_{S}\sim 1.5\times 10^{-2}$. In the very low viscosity regime ($Pr<2$) mode saturation requires an inverse relationship between $Pr$ and $M_{S}$.

In summary, by applying a simple passive stochastic ZFs model we were able to stabilise the $(2,1)$ tearing mode at a very low kinematic viscosity represented by the Prandtl number, $Pr$ for which the mode is linearly unstable. To stabilise the (2,1) mode for a wide range of $Pr>2$, we find that the required Mach number is around $M_{S}=5\times 10^{-3}$ which is of the same order as the ZFs that is expected to be generated by the small scale turbulence.

IV Results for semi-stochastic coupling model

In order to model the dynamics involved in the interaction of the visco-resistive MHD and the turbulence, in a simplified way, we have introduced a “semi-stochastic” coupling model, where the amplitude of the stochastic ZFs is defined as $\delta v_{S}=M_{S}v_{A}\langle(\tilde{j})^{2}+(\tilde{W})^{2}\rangle X(\rho,t)$, with $\langle\dots\rangle$ denoting average over the flux surfaces.

Here, we expect that as the amplitude of the MHD mode increases, the amplitudes of the stochastic ZFs also increase. Note that in our model, the stochasticity properties of the stochastic variable $X(\rho,t)$ are not changed. In this formulation, we attempt to model the MHD mode influence (via a direct cascade) on the interaction of the high $k$ turbulence with the mode itself.

For this model, we have examined two types of MHD instabilities: i) $(2,1)$ tearing mode and ii) $(1,1)$ kink mode. This is done by setting the initial $q$ and $j_{0}$-profiles such that the modes are linearly unstable. The time evolution of the $|\tilde{\psi}|_{max}$ for case (i) is shown in Figs. 5. Here, we examine the effect on the MHD stabilisation for different $M_{S}$ at $Pr=1$.

We start from the saturated phase of the $Pr=75$, $M_{S}=0$ (at $t=8.9$ms) simulation, and continue while applying the semi-stochastic coupling model, and reducing the Prandtl number to $Pr=1$ (for $t>8.9$ms). Figure 5 shows the results of these simulations for a wide range in $M_{S}=10^{-4}-10^{-1}$.

Our results show that a complex dynamic interaction exists between the MHD and ZFs. Initially, as $M_{S}$ is decreased (from $10^{-1}$ to $10^{-3}$), the stabilisation effect of the stochastic ZFs becomes less strong, and therefore, the $(2,1)$ mode saturates at higher amplitudes; compare red solid line at $M_{S}=10^{-1}$ to the green dashed line at $M_{S}=10^{-3}$. However, further reduction of the $M_{S}$ stabilises the $(2,1)$ at a faster rate; see black dotted line in Fig. 5 corresponding to $M_{S}=5\times 10^{-4}$.

By decreasing the $M_{S}$ to $10^{-4}$ increases this initial decay rate even further. However, after reaching a threshold level, it transiently rises to a higher amplitude before eventually saturating to the same amplitude as with $M_{S}=10^{-3}$; see pink dashed-dotted line in Fig. 5 corresponding to $M_{S}=10^{-4}$. A possible reason for this observation is that as the amplitude of the initial seed of ZFs becomes weaker i.e. lowering of the $M_{S}$, the (2,1) mode grows rapidly since at the same time we have significantly reduced the stabilisation effect of $Pr$. As the mode’s growth becomes more rapid, the increase in the amplitude of the seeded ZFs grows (due to the $\langle(\tilde{j})^{2}+(\tilde{W})^{2}\rangle$ factor). Therefore, the stabilisation effect of the ZFs becomes stronger and stronger increasing the damping rate of the $|\tilde{\psi}|_{max}$, as observed in Fig. 5 black dotted line compared to green dashed line.

Consider the difference between the black dotted line ($M_{S}=5\times 10^{-4}$) and the pink dashed-dotted line ($M_{S}=10^{-4}$), in Fig. 5. The rapid transient growth of the mode amplitude due to the significantly lower $M_{S}$ results in a very large rise of the non-linear factor $\langle(\tilde{j})^{2}+(\tilde{W})^{2}\rangle$. As expected, the mode drives the ZFs strongly leading to even faster damping. This results in the pink dashed-dotted line reaching a minimum at about $t\sim 9.2$ms well bellow the black dotted line, in Fig. 5. On the other hand, if the mode decays too rapidly to very low levels, the ZFs amplitude decays with it and so does its stabilisation impact on the (2,1) mode. As soon as the ZFs stability effect disappears, the non-linear growth of the mode rapidly results in an overshoot (see pink dashed-dotted line in Fig. 5 at $t=10.4$ms where the mode has reached reaching its maximum at $t\sim 10.4$ms well above the saturation level of the black dotted line). At this point, the ZFs become large enough to suppress the non-linear growth and stabilise the system to reach a new saturated state which is essentially similar to the previous cases. For a range of initial ZFs seeds with $M_{S}=10^{-4}-10^{-3}$, we observe a non-linear self-organziation interplay between ZFs, $(2,1)$ and its related low (m,n) spectrum that saturates to a similar final state.

Figures 6 (a-d) show a comparison of the profiles of the $dB_{r}/B_{0}$, $\tilde{W}$, equilibrium current density $j_{0}$, and safety factor $q$ between $t=8.9$ms and $t=11.7$ms, without and with the stochastic ZFs, respectively. The equilibrium current profile $j_{0}$ in the final saturated state, is almost back to its original linearly unstable profile, in the case where $M_{S}=10^{-1}$ (see Fig. 6 (c)). However, the saturated mode amplitude is low $dB_{r}/B_{0}\sim 10^{-4}$, thus, the overall ZFs Mach number is of the order of $10^{-3}$.

The vorticity profile in this case shows two distinct peaks with significantly higher amplitudes than in the absence of the stochastic ZFs at $\rho\sim 0.65$ and $\rho\sim 0.75$ (see blue dashed line in Fig. 6 (b)). The impact of this strong shear in the vorticity is observed on the sheared poloidal structure of the current density fluctuations with two peaks developed at $m=2$ surface, see the contour plots shown in Fig. 7(left) and the corresponding spectra plots in Fig. 7(right). Furthermore, a low amplitude $m=5,\;n=4$ mode is found to be active in the plasma core.

In case ii) for the $(1,1)$ mode we start with the Prandtl number of $Pr=30$ (as was reported in Ref. Mendonca2018 ), and after reaching a non-linear saturated state at time $t=2.64$ms, we apply the stochastic ZFs model as described above whilst reducing the Prandtl number to $Pr=1$. As can be seen in Fig. 8, for very small values of the $M_{S}=5\times 10^{-5}$ (pink dashed-dotted line), initially the $(1,1)$ mode grows strongly due to the lower Prandtl number, however as its amplitude increases, so does the amplitude of the stochastic ZFs, which eventually leads to stabilisation and saturation of the mode. We note that during this mode evolution, the equilibrium state also changes ($dj_{0}/dr$) due to the non-linear effects of the mode on itself, the so-called profile-mode interaction Mendonca2018 . After that the mode fluctuates around its original saturated level in the absence of the stochastic ZFs with $Pr=30$. By increasing $M_{S}$, a new saturated state can be reached but at a significantly lower amplitude (see for example, red solid line in Fig. 8 for $M_{S}=10^{-1}$). This dynamic behaviour is closely related to the dynamics that we described in the (2,1) case where the equilibrium current profile $j_{0}$ in the final saturated state, is almost back to its original linearly unstable profile (see Fig. 9 (c)). However, the saturated mode amplitude is low $dB_{r}/B_{0}\sim 10^{-5}$, thus, the overall ZFs Mach number is agin of the order of $10^{-3}$.

Figures 9 (a-d) show a comparison of the profiles of the $dB_{r}/B_{0}$, $\tilde{W}$, equilibrium current density $j_{0}$, and the safety factor $q$ at $t=2.64$ms and $t=4.64$ms, without and with the stochastic ZFs effects at $M_{S}=10^{-1}$, respectively. As can be seen here, addition of the stochastic ZFs results in a significant stabilisation of the $(1,1)$ mode, comparing the red solid line to the blue dash-dotted line in Fig. 9 (a). The distortion of the equilibrium current density due to the $(1,1)$ mode is strongly reduced as the model is applied. The vorticity profile in this case shows a more localised peak with larger amplitude at $\rho\sim 0.15$ (see blue dashed line in Fig. 9 (b)). The impact of this strong shear in the vorticity is observed on the poloidal structure of the current density fluctuations at $m=1$ surface, see the contour plots shown in Fig. 10(Top) and the corresponding spectra plots in Fig. 10(Bottom).

In summary, by applying a semi-stochastic ZFs model we were able to stabilise the $(2,1)$ tearing mode and the $(1,1)$ kink mode at a very low kinematic viscosity represented by the Prandtl number, $Pr$ for which the modes are linearly unstable. The main difference here with that of the simple stochastic model discussed in the previous section is that, the dynamic stabilisation is more effective even at lower $M_{S}$ and after an initial oscillation the saturation levels for different values of $M_{S}$ are similar due to the interaction between MHD and turbinate.

V Conclusion

A model is developed to simulate micro-scale turbulence driven ZFs, and their impact on the MHD tearing and kink modes is examined. The model is based on a stochastic representation of the micro-scale ZFs with a given Mach number, $M_{S}$. Two approaches were considered i) passive stochastic model where the ZF’s amplitudes are independent of the MHD mode amplitude, and ii) the semi-stochastic model where the amplitudes of the ZFs have a dependence on the amplitude of the MHD mode itself. We have purposely kept the resistive diffusivity ($c^{2}\eta/4\pi$), and the kinematic viscosity ($\nu$) as given functions of $\rho$, and during the time evolution of the system, in order to explore the advective effects of the stochastic ZFs on the vorticity equation. We recognise that in principle, the small-scale turbulence can also modify both of these transport properties for example as it is done in previous two-fluid simulations with CUTIE (Thyagaraja et. al. Thyagaraja2010 ). In neutral fluid dynamics, it has long been recognised (since the time of Prandtl), that turbulence modifies momentum and energy transport properties of the fluid significantly above their Laminar/molecular values.

The results for the approach i) show that the stochastic ZFs can significantly stabilise the (2,1) mode even at very low kinematic viscosity, $Pr$, where the mode is linearly unstable. For very low levels of $Pr<2$ the stochastic ZFs amplitude, described by its Mach number $M_{S}$ has to increase, in order for stabilisation to be effective. However, for $Pr>2$, and at similar values of $M_{S}$, we find that lower values of $Pr$ lead to lower saturation levels of the MHD mode. This is understood to be the result of stronger dissipation of the ZFs at higher $Pr$, which reduces their stabilisation effect on the MHD mode.

In the approach ii) we observe a more complex interaction between the ZFs and the MHD modes (both (2,1) and (1,1) modes), where by decreasing the amplitude of the initial seeded ZFs, i.e. $M_{S}$, the damping rate of the MHD modes actually increases even at linearly unstable $Pr$ levels. This dynamic behaviour is understood to be a result of a predator-prey type of dynamical stabilisation where the initial seeded ZFs are too small to damp the MHD mode. However since the stabilisation effect of the viscosity ($Pr$) has also diminished in these cases, the mode amplitude shoots up. As the amplitude of ZFs are dependent on the energy of the mode, they will grow very fast and back-react on the mode leading to a faster decay rate. For $M_{S}$ values below a threshold, the fast decay of the MHD also kills the ZFs, and hence, it starts to grow again. As the ZFs grow with the mode, at some point they reach a level that can stabilise the MHD mode once more. In our examples, we have not observed any repetition of this cycle, and after this initial interaction between the ZFs and MHD mode, they saturate to an invariant final state (for $10^{-5}<M_{S}<10^{-3}$). The equilibrium current density profile in this final state is close to the initial linearly unstable profile. We note however, that the possibility of a “Hopf bifurcation” with periodic waxing and wining of the mode and turbulence is not ruled out and may occur for different choices of parameters.

Our results, therefore indicate a possible mechanism for stabilisation of the MHD modes via small-scale perturbations in poloidal flow, simulating the turbulence driven ZFs. Generating small-scale perturbations in the vicinity of the MHD modes, without the need for high spatial accuracy, in principle can be achieved by means of Radio Frequency waves (RF). As it is already demonstrated, in fusion plasmas application of RF heating has a stabilisation impact on the core MHD modes such as Sawteeth Campbell1988 ; Larche2017 .