Zonal fields as catalysts and inhibitors of turbulence-driven
magnetic islands

D. Villa

Max Planck Institute for Plasma Physics,
Boltzmannstraße 2, 85748, Garching, Germany and
Excellence Cluster ORIGINS, Boltzmannstraße 2, D-85748, Garching, Germany*

N. Dubuit and O. Agullo

Aix-Marseille Université, CNRS, PIIM UMR 7345, Marseille, France

X. Garbet

CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France and
School of Physical and Mathematical Sciences,
Nanyang Technological University, 637371 Singapore

(Dated: December 12, 2024)
Abstract

A novel coalescence process is shown to take place in plasma fluid simulations, leading to the formation of large-scale magnetic islands that become dynamically important in the system. The parametric dependence of the process on the plasma $\beta$ and the background magnetic shear is studied, and the process is broken down at a fundamental level, allowing to clearly identify its causes and dynamics. The formation of magnetic-island-like structures at the spatial scale of the unstable modes is observed quite early in the non-linear phase of the simulation for most cases studied, as the unstable modes change their structure from interchange-like to tearing-like. This is followed by a slow coalescence process that evolves these magnetic structures towards larger and larger scales, adding to the large-scale tearing-like modes that already form by direct coupling of neighbouring unstable modes, but remain sub-dominant without the contribution from the smaller scales through coalescence. The presence of the cubic non-linearities retained in the model is essential in the dynamics of this process. The zonal fields are key actors of the overall process, acting as mediators between the competitive mechanisms from which Turbulence Driven Magnetic Islands can develop. The zonal current is found to slow down the formation of large-scale magnetic islands, acting as an inhibitor, while the zonal flow is needed to allow the system to transfer energy to the larger scales, acting as a catalyst for the island formation process.

Magnetic islands are a well-studied phenomenon in plasma physics, being the product of magnetic reconnection, a process that leads to the re-organization of the magnetic configuration and lowers the overall internal energy of the system. Magnetic reconnection itself is a phenomenon of interest to the whole domain of magnetized plasma physics, as it leads to localized heating and changes in the magnetic field topology, thus strongly affecting the properties and dynamics of the plasma [1-3]. In case the system allows for

it, the newly formed magnetic structures can form subsets of nested flux surfaces that can act as a localized confinement region, the magnetic island itself.

Despite the large volume of work done on magnetic islands [4-8], there are still open questions about the fundamental properties of these structures and the dynamics of reconnection, namely how exactly turbulence drives reconnection, thus forming Turbulence Driven Magnetic Islands (TDMIs), what are the differences in the dynamics for a plasma that undergoes reconnection at different plasma $\beta = \frac{2\mu_0 p}{B^2}$ regimes (where p is

<footnote>* daniele.villa@ipp.mpg.de</footnote>
the plasma pressure, $B$ the magnetic field strength and $\mu_0$ the vacuum permeability constant) and how do zonal fields (i.e. axisymmetric flows and currents) interact with magnetic islands. Of particular interest for this work are the dynamics of turbulent reconnection [7, 9–11], the generation of magnetic islands by turbulence [12–15] and the interaction of magnetic islands with zonal fields [16–19].

In this letter, the generation of magnetic islands by turbulence is shown in systems with varying plasma $\beta$, magnetic shear and instability drive. The relation of the new process to those previously explored in the literature [12, 13] is addressed, and, as will be detailed, this new process comes in addition to, and supercedes, previous processes leading to TDMI formation, in particular in near-marginal regimes relevant to low-$\beta$ plasmas. The core finding of this paper is illustrated in Fig. 1, where the transition of the system from the linear interchange instability to the non-linear formation of TDMIs through coalescence is shown as a sequence of isocontours of the field $\psi$ at different times.

To demonstrate the presence of these mechanisms and assess their importance, a reduced 6-field fluid model is considered, developed starting from the Braginskii fluid equations [20, 21], allowing to describe both small-scale phenomena, like interchange instabilities and the deriving turbulence, and large-scale ones, like magnetic islands and zonal fields, on the time scales needed to study (resistive) magnetic reconnection. The model evolves 6 dynamic fields: the magnetic potential $\psi$, the electrostatic potential $\phi$, the electron and ion pressure $p_{e/i}$, the density $n$ and the ionic parallel velocity $u_\|$. Reduction of the system to single-helicity “2D” slab geometry allows studying the fundamentals of the process in a simplified and generic environment, where the only requirement is a strong background magnetic field that allows defining the geometry of the system in reference to a chosen magnetic helicity (that identifies a “resonant position” [4, 22]). The background magnetic field is used to define a parallel and a perpendicular direction and establish drift-ordering expansions [23] (that are used to derive the model), that allow to study slow processes $\omega/\Omega_i \ll 1$, where $\Omega_i$ is the ion cyclotron frequency and $\omega$ the characteristic timescale of the studied phenomenon. The model being a “reduced” model means that fluctuations of the magnetic field in the parallel direction are considered to be higher order corrections to the dynamics, i.e. $\tilde{B}_\|$/$B_0 \ll \tilde{B}_\perp/B_0 < 1$. The $\sim$ indicates fluctuations, the parallel “$\|$” and perpendicular “$\perp$” directions are defined relative to the total magnetic field, and $\partial_z B_z = 0$, with $z$ being the direction of the dominant component
of the background magnetic field. Thus (after normalization) one can write the fluctuations of the magnetic field as $\tilde{B} = \nabla \times (\psi \hat{z})$. In slab geometry, the parallel derivative of a generic function $f$ can be expressed in a convenient form keeping into account the fluctuations of $\tilde{B} = \tilde{B}_{\perp}$ through the Poisson brackets: $\tilde{B} \cdot \nabla f = \{\psi, f\} = (\nabla \psi \times \nabla f) \cdot e_z$ where $\nabla$ is the gradient and $e_z$ the direction of the background out-of-plane magnetic field. The explicit expression for the Poisson bracket in slab geometry is: $\{\psi, f\} = (\partial_x \psi \partial_y f - \partial_x f \partial_y \psi)$. Here, $x$ is the radial direction, $y$ the periodic poloidal direction and any field $f$ can be decomposed according to $f = \sum_{m=-\infty}^{+\infty} f_m(x) \exp(2i\pi my/L_y + \varphi)$.

The full expression of the flux divergence terms is retained in the model, meaning that, in reduced notation, the parallel derivatives are given by the product of 3 fluctuating fields, e.g. $\nabla_\| (pu) = p \nabla_\| u_\| + u_\| \nabla_\| p = p\{\psi, u_\|\} + u_\|\{\psi, p\}$. Such terms, also used in other models [24, 25], are referred to as “cubic terms” [26], and it has been shown [21] that they can significantly impact the dynamics. Velocities in the system are normalized to the Alfvén velocity at the reference position $v_A = \frac{B_0}{\sqrt{\mu_0 mn_0}}$, lengths are normalized to a characteristic but arbitrary perpendicular length scale $L_\perp$, the magnetic field is normalized to the reference value $B_0$ and both pressure and density are also normalized

to the reference values at the resonant magnetic surface. A more detailed description of the model, including simulation parameters, is provided in the supplementary material and/or in the references [21]. This model is implemented in the fluid code AMON [27]. For the simulations of this work, only the interchange instability is present, driven by the coupling of the equilibrium ion pressure gradient to the curvature of the magnetic field. Having a single-helicity system implies that there is only one resonant position, where the magnetic field determines the reference helicity, which allows for clearer analysis of the mechanism. To limit the number of free parameters, the only background gradients in the system are those of the ion pressure and of the magnetic field. Note that in order to have a net 0 perpendicular equilibrium drift one has to impose $\partial_x \phi_{eq} = -\Omega_i \tau_A \rho_*^2 \frac{\partial_x p_{i\,eq}}{n_{eq}}$. The magnetic field is linearly stable to the tearing mode, which is confirmed by computing, at the resonant position $x=0$ the parameter $\Delta' \le -1.9$ [28, 29]. This value of $\Delta'$ is obtained by using a Harris-sheet background magnetic field with sufficiently weak gradient. The resistivity $\eta = 10^{-5}$ implies a typical time-scale for resistive reconnection $\tau_\eta = L_\perp^2/\eta \approx 10^4 \tau_A$ for a typical $L_\perp = 0.5$, meaning that phenomena occurring on time-scales faster than $\tau_\eta$ are due to turbulent dynamics. The simulations are
<figure><img src="image_1.png" /></figure>

FIG. 1. Sequence of isocontours of $\psi$ taken from the simulation with $\beta = 1.28\%$ and $\partial_x B_{eq} = 0.02$ showing the coalescence process taking place from the late linear phase (top) to the non-linear phase (bottom).

thus gradient-driven, though the mode $m = 0$ is self-consistently evolved, allowing for saturation dynamics to take place, should they be present.

Non-linear simulations are run scanning the parameter space as shown in Fig. 2. All simulations are interchange unstable. This identifies a threshold in $\beta$ and $s = \frac{\partial_x B_{eq}|_{x=0}}{B_{eq}|_{x=0}} = \partial_x B_{eq}|_{x=0}$, the magnetic shear, for the formation of TDMIs. The definition of the magnetic shear varies in the literature, here

$B_{eq}$ is the total equilibrium magnetic field, including the constant out-of-plane component, that appears at the denominator when writing $B_y$ as a function of the total magnetic field in a mono-helicity system, hence $|B_{eq}|_{x=0}| = 1$ and $s$ measures the inverse of the length scale of the variation of the “poloidal” magnetic field, but results presented in this paper are independent of the definition of the magnetic shear. Additional simulations have shown that the threshold is not related to the growth rate of the most unstable mode, but rather to the physical mechanisms acting at different shears, since simulations with very similar growth rates of the most unstable mode did or didn’t show the formation of TDMIs depending on magnetic shear.

Notice that the weaker the magnetic shear, the more likely the appearance of TDMIs, reinforcing the point that the structures presented here do not originate from tearing [28]. The non-linear dynamics of the system further weaken the magnetic shear, and non-linear de-stabilization of tearing can also be ruled out.

What is actually meant by “formation of TDMIs” is illustrated in Fig. 3, where the radial structure of the mode $m = 14$ for two simulations with $\beta = 1.28\%$ but different magnetic shears $\partial_x B_{eq} = 0.02$ and $\partial_x B_{eq} = 0.04$ are compared. Since the modes are complex
valued, the plot uses the absolute value of
the function to represent its magnitude, and
the colour scale to represent the phase dif-
ference of the mode at that particular posi-
tion to the phase at the resonant position.
For a mode with interchange parity, in a lin-
ear simulation, the eigenfunction of $\psi$ is odd
(meaning that the real and imaginary parts
are odd), corresponding to an even absolute
value with a change in phase across the res-
onance, while for tearing parity the phase is
constant across the resonance, corresponding
to an even eigenfunction. Since simulations
are non-linear, the words “even” and “odd”
are not to be understood as exact proper-
ties, but refer to the “dominant” component
of the function, which does not necessarily
match the linear one. Focusing on $\psi$, as the
formation of magnetic islands is directly vis-
ible through this field, Fig. 3 refers to $t =$
$1700\tau_A$ (i.e. early in the non-linear phase)
in two simulations where, at that particular
time, the mode $m = 14$ had the highest en-
ergy: the simulation that develops TDMIs ($
\partial_x B_{eq} = 0.02$) has a broad region of weakly
varying phase around the resonance, whereas
the simulation that doesn't form TDMIs has
a clear variation in phase across the reso-
nance. Averaging the absolute values of these
phase differences $\Delta\varphi = (\varphi - \varphi_{res})/(2\pi)$ for
$-1 \le x \le 1$, where the amplitude $A$ of the
modes satisfies $A/A_{max} \ge 0.1$, gives the aver-

age $\langle|\Delta\varphi(x)|\rangle_{x\in[-1,+1]}$ values in Fig. 3, show-
ing a noticeable lower average phase differ-
ence when TDMIs form. Using the phase in-
formation to quantify the parity of the modes
also allows to compensate for the shift of
the resonant position in the non-linear phase,
that would act as the symmetry center to
define the parity of the mode. Performing
this analysis over time and for all modes
shows this behaviour of the average phase dif-
ference to be consistent over time, at least
for all modes with mode-number $m \lesssim m^*$
($m^*$ being the mode-number of the linearly
most unstable mode). Thus, in simulations
where the formation of TDMIs is observed,
modes develop, on average, more tearing-like
structures in their non-linear evolution, even
though they start out with interchange par-
ity.

Whether this change in structure of the un-
stable modes takes place determines whether
or not the large-scale magnetic islands will
form. What allows this change in structure
is a combination of factors, including the ve-
locity shear, whose role, along with that of
the zonal current, will be expanded upon in
the following section. However, suppressing
the fluctuations (and the equilibrium) of the
mode $\phi_0$ does not prevent the change in par-
ity, so the mechanism requires also the pa-
rameter dependencies illustrated in Fig. 2
and the cubic non-linearities retained in the
model. This latter point has been verified
by running simulations without cubic non-
linearities, that don't show change in par-
ity of the unstable modes. Indeed, in sim-
ulations without cubic terms, one recovers
the results from the literature [10, 13]: di-
rect coupling of the mode $m^*$ with the mode
$m^* \pm 1$ drives non-linearly an $m = 1$ mode
with tearing parity, whose growth rate sat-
isfies $\gamma \approx 2\gamma^*$. However, such modes never
become dynamically important. Only if the
unstable modes (with $m \approx m^*$) change par-
ity, will a clearly identifiable magnetic island
form in the system, after the non-linear en-
ergy transfer across scales takes place. The
small-scale islands that form at the unsta-
ble scales merge (coalesce) into larger islands
over time, allowing the sub-dominant low-
$m$ modes, formed by direct coupling, to be-
come dynamically relevant and form large-
scale magnetic islands. Thus without the
mechanism described here, no large-scale is-
land would be visible in the system (see the
supplementary material for figures illustrating
this process).

The role of the cubic terms (and partially
of $\beta$ and $\partial_x B_{eq}$) can be understood by fo-
cusing on the pressure term in Ohm's law:
$$ \frac{\Omega_i \tau_A \rho_*^2}{n} \{\psi, p_e\} $$
In an interchange unstable sys-
tem, the parity of $\psi_{m*}$ in the linear (and early
non-linear) phase is odd, and non-linear dy-
namics form modes $m=0$ of pressure and

density that also have odd parity (to give
profile flattening). The Poisson bracket of
two odd functions is odd, multiplied by $n_{m=0}$
leads to an even term, at the same $m^*$ of the
linear instability. This changes the parity of
the unstable mode.

It is also important to remark that the phase
change described here happens on much
faster time-scales than the resistive time-
scale (in the late linear phase it's possible to
quantify the average behaviour of the phases
described above), and thus is attributed to
non-linear (turbulent) behaviour. The non-
linear energy transfer also occurs on time-
scales faster than the resistive one, so that
overall the process is quicker than the resis-
tive time-scale.

Summarizing, the novelty of the process de-
scribed here is the fact that the change in
structure of the unstable modes is allowed by
a combination of the background parameters
and the presence of the cubic terms. Through
the above non-linear dynamics, small-scale
TDMI form. They allow the otherwise sub-
dominant large-scale TDMIs formed by the
direct coupling described in the literature
[12, 13, 15] to appear in the system.

These simulations only allowed seeing islands
with $m=2$, while one might expect to
observe islands with $m=1$: for simula-
tions with strong instability drive, the islands
reached the boundaries of the domain, re-
<figure><img src="image_2.png" /></figure>

FIG. 2. Summary of results for 2D slab non-linear runs of linearly unstable interchange where the plasma $\beta$ and background magnetic shear were varied. Green circles represent simulations that showed the formation of magnetic islands and at least the initial phases of coalescence, while simulations indicated with red “X”s didn’t. The mode number $m^*$ and the corresponding growth rate $\gamma^*$ on top of the symbols refer to the most unstable mode. Circled simulations are used in Fig. 3

gardless of the size of the latter, before the mode $m = 1$ formed, while for simulations with weak drive, the non-linear energy transfer slows down progressively, rendering the computational cost of observing the $m = 1$ mode prohibitive. The strong instability case might point to a missing saturation mechanism in our simulations, but this will have to be addressed in future studies.

The non-linear dynamics of the system are dominated by the interaction among the structures at the largest scales. These struc-

<figure><img src="image_3.png" /></figure>

FIG. 3. Radial structure and phase differences of the mode $m = 14$ of $\psi$ for two non-linear simulations with $\beta = 1.28\%$ but different magnetic shear as indicated in the legend ($\partial_x B_{eq} = 0.02$ forms TDMIs, $\partial_x B_{eq} = 0.04$ doesn’t). The differences in phase with the reference value at the resonance are computed across the radius, along with the average $\langle|\Delta\varphi|/2\pi\rangle$, over the portion of the radius between the two dashed vertical lines at $x = \pm 1$. Notice that for $\partial_x B_{eq} = 0.02$ the mode has more even (tearing-like) parity. A noticeable lower average phase difference shows up consistently through simulations that form TDMIs for all modes with $m < m^*$.

tures are the TDMI (with poloidal mode-number $m = 2$), the zonal current (i.e. mode $\psi_{m=0}$) and the zonal flow (i.e. mode $\phi_{m=0}$). The fluctuations of the mode $\psi_0$ form a radial region around the resonance where the poloidal magnetic field and its shear vanish. In the non-linear phase, the flow becomes strongly sheared at the resonance, which is expected since that is where the instability is driven, and also at the radial position where the magnetic island has its maximum radial
width. As the island grows in width, this
latter region of sheared flow moves with the
growing island away from the resonant sur-
face, which is indicative of a strong interde-
pendence between the magnetic island and
the localization of the velocity shear.

Shown in Fig. 4 is the time evolution of
the energies of the (selected) modes of $\psi$ and
the zonal flow ($\phi_0$) for the simulation with
$\partial_x B_{eq} = 0.01$ and $\beta = 1.28\%$. It's visible
how in the linear phase the linearly unstable
modes (the curve labeled as $\psi_\Sigma$ is obtained
by summing the energies for all modes with
$10 \le m \le 40$) grow and drive, by non-linear
coupling, the large-scale modes ($\psi_2, \psi_3, \psi_4$).

A period of interaction with the zonal fields,
visible in the quick oscillations of the zonal
current $\psi_0$, enhances the energy transfer from
the small to the large scales. In simula-
tions that did not show the formation of
TDMIs this phase is missing, and the zonal
fields have much “smoother” energy dynam-
ics. Furthermore, in simulations that did not
form TDMIs, the large-scale modes remain at
low energy levels, while here a shift in the en-
ergy distribution is visible quite early in the
non-linear phase. The importance of modes
$m=0$ on the dynamics of coalescence has al-
ready been highlighted also in other contexts
such as in [5]. Late in the simulation the
mode $\psi_2$ becomes dominant, and the large-
scale modes acquire significantly larger en-

ergies than the turbulent scales. One could
expect the mode $\psi_1$ to become dominant over
time, but, as already addressed, this was not
observed.

To further address the role of the zonal fields,
variations on the simulations were run where
the fluctuations of the zonal fields were sup-
pressed either through the use of a stronger
dissipation only to the mode $m=0$ of the
electrostatic potential $\phi$ from the beginning
of the simulations, or by completely eliminat-
ing fluctuations of the mode $m=0$ of $\psi$ or of
$\phi$ from a certain time in the non-linear phase
(once the $m=2$ magnetic island was already
dominant) onwards. When the zonal flow is
suppressed for a simulation that has already
reached the non-linear phase, it is observed
that the coalescence is halted at the current
wave-number, and the island's width stops
evolving, indicating that the zonal flow is re-
sponsible for the transfer of energy at larger
scales (i.e. the coalescence), acting as a cat-
alyst for the formation of TDMIs. This is an
unusual behaviour for zonal flows, that usu-
ally move energy towards the smaller scales
by shearing apart larger eddies [17], but it's
also peculiar that in these simulations the
properties of turbulence aren't much affected
by the disappearance of the zonal flow.

Suppression of the zonal current for a sim-
ulation that has already reached the non-
linear phase does the opposite, accelerating
the transition of energy to larger scales, while strongly suppressing turbulence. This can be understood as the re-establishing of the background magnetic shear, that limits the region where turbulence can develop. Thus the zonal current acts as an inhibitor to the formation of TDMIs, slowing down the coalescence by creating the low-shear region, that favours the accumulation of energy at the turbulent scales.

It has long been known from gyrokinetic analytical calculations [30] and simulations [31] that zonal flow saturation is not dominated by fluid dissipation. Given the important role of zonal flows described above, one should check that our results are not reliant on zonal flow dissipation. Starting simulations anew with varied dissipation for the mode $\phi_0$ shows, as expected, the existence of two regimes of zonal flow saturation (ideal and collisional). Specifically, increase of the dissipation of the zonal flow by a factor 50 doesn't lead to any effect, while a further factor of 5 finally gives a factor $\sim 1/5$ on zonal flow saturation levels. The simulations discussed in this paper are therefore in the ideal regime, where both zonal fields saturate at similar levels regardless of the magnitude of the dissipation.

In conclusion, this paper shows how the formation of TDMIs can happen through a process never observed before that involves

<figure><img src="image_4.png" /></figure>

FIG. 4. Evolution of the energies for the simulation with $\partial_x B_{eq} = 0.01$ and $\beta = 1.28\%$. For the magnetic potential $\psi$ multiple modes are shown to highlight the different phases of the simulation, while for the electrostatic potential $\phi$ only the mode $m = 0$ (i.e. the zonal flow) is shown. The curve labeled as $\psi_{\Sigma}$ is obtained by summing the energies for all modes with $10 \le m \le 40$.

small-scale islands forming at the unstable scales and then merging to progressively larger scales. This happens through a coalescence process that allows the otherwise negligible large-scale islands formed by direct non-linear coupling to become dominant. Coalescence occurs at a later time than the direct coupling of the most unstable modes reported in the literature (e.g. [13]), and, since it requires energy transfer across scales, is slower than the direct coupling (direct coupling takes $\approx 10(\gamma^*)^{-1}$, in Fig. 4 this would be $t \approx 800\tau_A$, while coalescence takes $\ge 50(\gamma^*)^{-1}$), but it modifies the dynamics of the system fundamentally, and, in particular, leads to much larger magnetic islands than otherwise observed. Simulations that form
magnetic islands, unlike those that don't, undergo a change in the radial structure of the magnetic potential $\psi$, rendering the phase of the mode more uniform over the radial coordinate (Fig. 3). The equilibrium magnetic shear of the system, and the presence of the cubic terms in the model, are crucial in allowing this change in parity, and thus the coalescence, to take place. This process is of particular importance for fusion plasmas, since it allows the formation and growth of magnetic islands independent of neoclassical physics, and it might feedback into the dynamics of NTMs. The latter are expected to pose a threat to future devices, operating at higher $\beta$ than existing ones (and than explored here), in particular since the source of the seed islands is, in some cases, unexplained [32]. On the other hand, astrophysical plasmas occur more often in conditions where $\beta \approx 1$, such that pressure fluctuations can play a much larger role on the dynamics of the system, allowing the cubic terms to impact the dynamics much more significantly, making the process described here relevant.

The TDMI is capable of driving a strongly sheared poloidal flow at the position where its separatrix has its maximum radial width, and the presence of the zonal flow is needed for the coalescence process to continue throughout the simulation. The zonal flow

is mostly driven by the turbulence and the magnetic islands, rather than being regulated by dissipation. Given the strong role in the establishment of Internal Transport Barriers [33, 34] played by zonal fields, this might be an unexpected positive outcome of the presence of TDMIs. The zonal current, instead, slows down the coalescence process, favouring the development of small-scale turbulence by creating a region of weak magnetic shear.

## SUPPLEMENTARY MATERIAL

The supplementary material to this letter contains the equations of the model, including normalizations and parameters used to run the simulations. There is also a series of figures illustrating the change of phase and the coalescence through the isocontours of the magnetic potential $\psi$, as well as the crucial role played by the cubic terms, and an explanation of the use of phase information to analyze the parity of the modes.

## ACKNOWLEDGMENTS

The authors would like to thank the referees for their constructive and detailed feedback.

The authors would like to thank S. Mazzi
and M.J. Pueschel for constructive discussions on the topic.

The project leading to this publication has
received funding from the Excellence Initia-
tive of Aix-Marseille University - A*Midex, a
French "Investissements d'Avenir" program
AMX-19-IET-013. The simulations in this
article were run thanks to the support of EU-
ROfusion and MARCONI-Fusion. This work
has been carried out within the framework
of the EUROfusion Consortium, funded
by the European Union via the Euratom
Research and Training Programme (Grant
Agreement No 101052200 — EUROfusion).
Views and opinions expressed are however

those of the author(s) only and do not
necessarily reflect those of the European
Union or the European Commission. Neither
the European Union nor the European
Commission can be held responsible for
them. This work was partially funded
by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under
Germany's Excellence Strategy - EXC 2094
- 390783311.

The authors have no conflicts to disclose.

The data that support the findings of this
study are available from the corresponding
author upon reasonable request.

[1] DA Falconer, GA Gary, RL Moore, and JG Porter. An assessment of magnetic conditions for strong coronal heating in solar active regions by comparing observed loops with computed potential field lines. *The Astrophysical Journal*, 528(2):1004, 2000.

[2] M Bárta, B Vršnak, and M Karlický. Dynamics of plasmons formed by the current sheet tearing. *Astronomy & Astrophysics*, 477(2):649–655, 2008.

[3] Miho Janvier. Three-dimensional magnetic reconnection and its application to solar flares. *Journal of Plasma Physics*, 83(1), 2017.

[4] FL Waelbroeck. Theory and observations of magnetic islands. *Nuclear Fusion*, 49(10):104025, 2009.

[5] A Poyé, O Agullo, A Smolyakov, S Benkadda, and X Garbet. Dynamics of magnetic islands in large $\delta'$ regimes. *Physics of Plasmas*, 21(2):020705, 2014.

[6] Akihiro Ishizawa, Yasuaki Kishimoto, and Yuji Nakamura. Multi-scale interactions between turbulence and magnetic islands and parity mixture—a review. *Plasma Physics and Controlled Fusion*, 61(5):054006, 2019.

[7] Minjun J Choi. Interaction between a magnetic island and turbulence. *Reviews of*
*Modern Plasma Physics*, 5(1):1–39, 2021.

[8] H. Betar, D. Del Sarto, M. Ottaviani, and A. Ghizzo. Microscopic scales of linear tearing modes: a tutorial on boundary layer theory for magnetic reconnection. *Journal of Plasma Physics*, 88(6):925880601, 2022.

[9] A Lazarian and Ethan T Vishniac. Reconnection in a weakly stochastic field. *The Astrophysical Journal*, 517(2):700, 1999.

[10] A Ishizawa and N Nakajima. Turbulence driven magnetic reconnection causing long-wavelength magnetic islands. *Physics of Plasmas*, 17(7):072308, 2010.

[11] Alex Lazarian, Gregory L Eyink, Amir Jafari, Grzegorz Kowal, Hui Li, Siyao Xu, and Ethan T Vishniac. 3d turbulent reconnection: Theory, tests, and astrophysical implications. *Physics of Plasmas*, 27(1), 2020.

[12] Masatoshi Yagi, Sanae-I Itoh, Kimitaka Itoh, Masafumi Azumi, Patrick H Diamond, Atsushi Fukuyama, and Takayuki Hayashi. Nonlinear drive of tearing mode by microscopic plasma turbulence. *Plasma and Fusion Research*, 2:025–025, 2007.

[13] M Muraglia, O Agullo, S Benkadda, M Yagi, X Garbet, and A Sen. Generation and amplification of magnetic islands by drift interchange turbulence. *Physical review letters*, 107(9):095003, 2011.

[14] O Agullo, M Muraglia, S Benkadda, A Poyé, N Dubuit, and X Garbet. Non-

linear dynamics of turbulence driven magnetic islands. i. theoretical aspects. *Physics of Plasmas*, 24(4):042308, 2017.

[15] N. Dubuit, O. Agullo, M. Muraglia, J. Frank, X. Garbet, and P. Maget. Dynamics of magnetic islands driven by ballooning turbulence. *Physics of Plasmas*, 28(2):022308, 2021.

[16] PN Guzdar, RG Kleva, A Das, and PK Kaw. Zonal flow and zonal magnetic field generation by finite $\beta$ drift waves: A theory for low to high transitions in tokamaks. *Physical review letters*, 87(1):015001, 2001.

[17] Patrick H Diamond, SI Itoh, K Itoh, and TS Hahm. Zonal flows in plasma—a review. *Plasma Physics and Controlled Fusion*, 47(5):R35, 2005.

[18] A Fujisawa, K Itoh, A Shimizu, H Nakano, S Ohshima, H Iguchi, K Matsuoka, S Okamura, T Minami, Y Yoshimura, et al. Experimental evidence of a zonal magnetic field in a toroidal plasma. *Physical review letters*, 98(16):165001, 2007.

[19] Ge Dong, Jian Bao, Amitava Bhattacharjee, and Zhihong Lin. Nonlinear saturation of kinetic ballooning modes by zonal fields in toroidal plasmas. *Physics of Plasmas*, 26(1):010701, 2019.

[20] SI Braginskii. Transport processes in a plasma. *Review of Plasma Physics*, 1:216,
1965.

[21] Daniele Villa, Nicolas Dubuit, Olivier Agullo, Alexandre Poyé, Xavier Garbet, and Andrei Smolyakov. Localized compressional self-heating in magnetic islands. *Journal of Plasma Physics*, 88(6):905880613, 2022.

[22] Bruce D. Scott. The mechanism of self-sustainment in collisional drift wave turbulence. *Physics of Fluids B*, 4(8):2468-2494, August 1992.

[23] Bruce Scott. *Turbulence and Instabilities in Magnetised Plasmas*, volume 2 of 2053-2563. IOP Publishing, 2021.

[24] M Giacomin, P Ricci, A Coroado, G Fourestey, D Galassi, E Lanti, D Mancini, N Richart, LN Stenger, and N Varini. The gbs code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary. *Journal of Computational Physics*, page 111294, 2022.

[25] Christoph Pitzal, Andreas Stegmeir, Wladimir Zholobenko, Kaiyu Zhang, and Frank Jenko. Landau-fluid simulations of edge-sol turbulence with grillix. *Physics of Plasmas*, 30(12), 2023.

[26] D Villa, X Garbet, N Dubuit, and O Agullo. On cubic non-linearities in fluid and kinetic models. In *Journal of Physics: Conference Series*, volume 2397, page 012005. IOP Publishing, 2022.

[27] Alexandre Poye. *Dynamique des ilots magnétiques en présence de feuille de courant et en milieu turbulent*. PhD thesis, Aix-Marseille, 2012.

[28] Harold P Furth, John Killeen, and Marshall N Rosenbluth. Finite-resistivity instabilities of a sheet pinch. *The physics of Fluids*, 6(4):459-484, 1963.

[29] N Arcis, DF Escande, and M Ottaviani. Rigorous approach to the nonlinear saturation of the tearing mode in cylindrical and slab geometry. *Physics of plasmas*, 13(5):052305, 2006.

[30] MN Rosenbluth and FL Hinton. Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks. *Physical review letters*, 80(4):724, 1998.

[31] Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang, and R. B. White. Turbulent transport reduction by zonal flows: Massively parallel simulations. *Science*, 281(5384):1835-1837, 1998.

[32] Akihiko Isayama, Go Matsunaga, and Yoichi Hirano. Onset and evolution of m/n= 2/1 neoclassical tearing modes in high-β p mode discharges in jt-60u. *Plasma and Fusion Research*, 8, 2013.

[33] S Benkadda, P Beyer, N Bian, C Figgarella, O Garcia, X Garbet, Ph Ghendrih, Y Sarazin, and PH Diamond. Bursty transport in tokamak turbulence: Role of zonal
flows and internal transport barriers. Nu-
clear fusion, 41(8):995, 2001.

[34] Katsumi Ida and Takaaki Fujita. Internal transport barrier in tokamak and helical plasmas. *Plasma Physics and Controlled Fusion*, 60(3):033001, 2018.