# Quasi-two-dimensionality of three-dimensional, magnetically dominated,
decaying turbulence

Shreya Dwivedi1,2,† Chandranathan Anandavijayan1,⋆ Pallavi Bhat1,‡
1International Centre for Theoretical Sciences, Tata Institute of Fundamental
Research, Bangalore 560089, India  2Department of Physics, University of
Wisconsin-Madison, Madison, WI 53706, USA

###### Abstract

Decaying magnetohydrodynamic (MHD) turbulence is important in various
astrophysical contexts, including early universe magnetic fields, star
formation, turbulence in galaxy clusters, magnetospheres and solar corona.
Previously known in the nonhelical case of magnetically dominated decaying
turbulence, we show that magnetic reconnection is important also in the fully
helical case and is likely the agent responsible for the inverse transfer of
energy. Again, in the fully helical case, we find that there is a similarity
in power law decay exponents in both 2.5D and 3D simulations. To understand
this intriguing similarity, we investigate the possible quasi-two-
dimensionalization of the 3D system. We perform Minkowski functional analysis
and find that the characteristic length scales of a typical magnetic structure
in the system are widely different, suggesting the existence of local
anisotropies. Finally, we provide a quasi-two-dimensional hierarchical merger
model which recovers the relevant power law scalings. In the nonhelical case,
we show that a helicity-based invariant cannot constrain the system, and the
best candidate is still anastrophy or vector potential squared, which is
consistent with the quasi-two-dimensionalization of the system.

###### keywords:

plasmas – magnetic fields – magnetohydrodynamics(MHD) – turbulence – magnetic
reconnection

†† ${\dagger}$<EMAIL_ADDRESS>
$\star$<EMAIL_ADDRESS>
${\ddagger}$<EMAIL_ADDRESS>
##  1 Introduction

Astrophysical plasmas are typically in a state of turbulence. In the aftermath
of energetic processes that are responsible for generation of the turbulence,
the plasma undergoes relaxation, leading to the decay of velocity and magnetic
fields. Many astrophysical plasmas such as solar and stellar coronae,
accretion disks, and pulsar magnetospheres host magnetic fields of significant
strength and are magnetically dominated (Pevtsov et al., 2003; Blandford et
al., 2017; Narayan et al., 2003). Such magnetic turbulent systems are thought
to relax via magnetic reconnection (Taylor, 1986). Reconnection is a process
by which magnetic fields undergo topological reconfiguration accompanied by
conversion of magnetic energy into heat, radiation and particle acceleration.
These can be highly energetic events and are thought to explain solar flares,
magnetospheric substorms, gamma-ray bursts etc (Ruan et al., 2020; Hesse &
Cassak, 2020; McKinney & Uzdensky, 2012).

More recently, it was shown that magnetic reconnection plays an important role
in 3D MHD decaying turbulence (Bhat et al., 2021; Zhou et al., 2020) further
extending support to the idea of reconnection playing an important role in
relaxation. Some of the astrophysical scenarios are either theoretically
modelled or already well understood to relax via magnetic reconnections. For
example, in the Sun, magnetic fields in the form of flux ropes extending up
into the corona experience twists due to the motion of the foot points in the
photosphere. The Parker model of coronal heating directly employs magnetic
reconnection to relax such twisted magnetic flux ropes (Parker, 1983).
Inevitably, the nonlinear evolution of reconnecting flux ropes leads to
turbulence (Rappazzo et al., 2008), accompanied by growth in magnetic
structure size (Pontin et al., 2011). Another context is the interaction
between solar wind and earth’s magnetotail, which leads to reconnection. This
involves magnetic flux loading and subsequent release of energy resulting in
magnetospheric substorms (McPherron, 1979). Thereafter, the earth’s magnetic
field relaxes to its ground state, and these processes again feature
turbulence (Daughton et al., 2011). In pulsar wind nebulae, the striped wind
and large magnetic hoops are thought to relax via reconnections in decaying
turbulence, leading to nonthermal emission and flares (Zrake, 2016; Zrake &
Arons, 2017). Additionally, the magnetic dissipation due to post-shock
instabilities seems to explain the lower strength of the magnetic field from
observations (Porth et al., 2014).

Magnetic reconnection, as a model for relaxation, can also naturally explain
the phenomenon of inverse energy transfer. Inverse energy transfer occurs
when, as the turbulent system relaxes, energy transfers to larger scales. In
the 3D hydrodynamic case, relaxation of homogeneous and isotropic turbulence
occurs by selective decay, i.e. energy at the smallest scales decays first,
followed by decay of successively larger scales. This process does not lead to
any inverse energy transfer but depicts what is known as the “permanence of
large eddies” (Batchelor & Taylor, 1949; Landau & Lifshitz, 1975; Lesieur &
Ossia, 2000). On the other hand, it is well known that decaying 2D
hydrodynamic turbulence can depict inverse energy transfer because here a
different version of selective decay manifests. Here, one ideal invariant
(kinetic energy) decays more slowly than another (enstrophy). In fact, the
slower decay invariably is due to its condensation at larger scales, thus
being consistent with the original understanding of selective decay (Matthaeus
& Montgomery, 1980; Ting et al., 1986). While permanence of large eddies just
constitutes the continuation of the flow at larger scales (or smaller
wavenumbers that are a part of the forward slope) at similar amplitudes, in
other cases like 2D hydrodynamics or 3D MHD with magnetic helicity, the slower
decaying ideal invariant can cause an actual inverse transfer in energy
(involving increase of energy at larger scales or smaller wavenumbers)
(Mininni & Pouquet, 2013).

In the MHD case, it is easy to observe that in a 2D decaying system, an
inverse transfer emerges as a consequence of reconnection between magnetic
islands (Zhou et al., 2019). A similar process in 3D could lead to larger and
larger magnetic structures at the cost of loss/decay in magnetic energy. The
study in Bhat et al. (2021) (BZL21) shed light on the possibility of magnetic
reconnection driving such a 3D inverse energy transfer. BZL21 were motivated
by the observation that the decay laws in 3D simulations matched with that
seen in 2D. They identified the dynamical timescale in this system of
nonhelical decaying turbulence as the reconnection timescale, $\tau_{rec}$.
This result challenges the conventional notion that the dynamical time scale
in a magnetically dominated system is typically the Alfvenic time scale. Thus,
the inverse transfer phenomenon was shown to occur on a slower time scale.

The inverse transfer phenomenon is useful in astrophysical scenarios where the
coherence length scale of the magnetic field is important as well. The ISM
magnetic fields that first feed onto a central massive object need to be
sufficiently coherent to trigger magnetorotational instability needed for
angular momentum transport in accretion disks (Bhat et al., 2017). In galaxy
clusters, large Faraday rotation measures need coherent magnetic fields
(Subramanian et al., 2006; Bhat & Subramanian, 2013). In IGM, the non-
detection of GeV photons that should have originated from cascading of TeV
photons emanating from blazars again needs coherent magnetic fields for
deflecting charged particles (Neronov & Vovk, 2010; Subramanian, 2016; Hosking
& Schekochihin, 2023). The magnetic fields in ISM, galaxy clusters and IGM, as
mentioned in the examples above, could have benefited from an inverse transfer
in magnetic energy during a turbulent decaying phase (when the dynamo was
inactive).

Direct numerical simulations of magnetically dominated turbulence have shown
evidence for inverse energy transfer even when an ideally conserved quantity
like magnetic helicity was absent (Brandenburg et al., 2015; Zrake, 2014).
BZL21 postulated that anastrophy conservation (in the absence of magnetic
helicity) is responsible for the inverse transfer phenomenon in 3D. They
highlighted the link between anastrophy conservation in the large Lundquist
number limit and magnetic reconnection in 2D and claimed that this carries
over to the case of 3D as well, thus explaining the exact similarity in the
decay law between 2D and 3D. However, Hosking & Schekochihin (2021) introduced
a new integral called the Saffman helicity invariant based on magnetic
helicity fluctuations, which they claimed is responsible for constraining the
decay law (this does not give $t^{-1}$, but gives $t^{-1.18}$ instead). They
extended the BZL21 result of $\tau_{rec}$ (reconnection timescale) being the
dynamical timescale to the helical case as well and found that it leads to a
$t^{-4/7}$ decay law.

In our work here, we examine helical decaying MHD turbulence and find that,
indeed, reconnection is important in this case as well. It makes sense that
reconnection is a common explanation for the phenomenon of structures growing
larger and larger as the decay progresses in both helical and nonhelical
inverse energy transfer. Intriguingly, the 2D(2.5D)-3D correspondence holds in
the helical case as well; the decay law in 2.5D (magnetic helicity cannot be
defined in the strict 2D case) is the same as that seen in 3D. We also assess
the so-called Saffman helicity integral critically and find that it is most
likely ruled out given the 2D-3D correspondence.

To assess the quasi-two-dimensional nature of the magnetically dominated
decaying turbulence, we analyse the system using Minkowski functionals. These
were originally introduced as a tool to study spatial patterns in galaxy
cluster data (Mecke et al., 1994). In cosmology, it is well-known that there
is clustering of galaxies (Peebles, 1980). It was found that the two-point
correlation function for large-scale galaxy clustering is insufficient due to
non-gaussianity. Minkowski functionals, on the other hand, contain information
related to all the correlation functions of the order and above two-point
correlation (Wiegand et al., 2014). In this paper, we use these Minkowski
functionals to examine the morphology of reconnecting fields and the
associated current density structures. This tool was first used in this manner
of assessing shapes in a scalar field in the field of cosmology and structure
formation (Sahni et al., 1998). It has also been previously used in the field
of astrophysical fluid dynamics for analysing magnetic field structures in
simulations of small-scale dynamo (Wilkin et al., 2007; Seta et al., 2020).
And to assess current density structures in strong MHD turbulence (Zhdankin et
al., 2014). We primarily employ Minkowski functionals to quantify the inherent
anisotropy within the system. This approach aids in refining our comprehension
of the quasi-two-dimensional manifestation of the reconnections.

The following section covers the MHD equations and parameters employed in our
direct numerical simulations. Section 3 delves into the analysis of data and
results, with subsection 3.1 dedicated to the helical case and subsection 3.2
focusing on the nonhelical case. Within these subsections, we present evidence
supporting the 2D-3D similarity in decay laws. Subsequent analyses compare
conserved quantities and critically evaluate various integrals believed to
constrain the evolution of these decaying systems. Section 4 presents results
from the Minkowski functional analysis, followed by the introduction of an
analytical quasi-2D hierarchical merger model in Section 5. The concluding
sections, 6 and 7, consist of discussions and overall conclusions
respectively. In the Appendix Sections A and B, we delineate the algorithms
employed for Minkowski functional analysis.

##  2 Numerical Setup

We have run direct numerical simulations of decaying helical and nonhelical
MHD turbulence in 2.5D and 3D cubes, using Pencil-Code111http://pencil-
code.nordita.org, with periodic boundary conditions. We solve the continuity
equation, momentum equation and the uncurled version of induction equations on
a Cartesian $N^{2}$ or $N^{3}$ grid, where $N$ is the number of grid points in
any direction.

$\displaystyle\frac{D~{}\mathrm{ln}~{}\rho}{Dt}$ $\displaystyle=$
$\displaystyle-\bm{\nabla}\cdot\bm{u},$ (1) $\displaystyle\frac{D\bm{u}}{Dt}$
$\displaystyle=$ $\displaystyle-
c_{s}^{2}\bm{\nabla}~{}\mathrm{ln}~{}\rho+\frac{\bm{J}\times\bm{B}}{\rho}+\frac{\bm{F}_{visc}}{\rho},$
(2) $\displaystyle\frac{\partial\bm{A}}{\partial t}$ $\displaystyle=$
$\displaystyle\bm{u}\times\bm{B}-\eta\mu_{0}\bm{J},$ (3)

where $D/Dt=\partial/\partial t+\bm{u}\cdot\bm{\nabla}$ is the advective
derivative, $\bm{u}$ is the fluid velocity field,
$\bm{B}=\bm{\nabla}\times\bm{A}$ is the magnetic field, with $\bm{A}$ the
magnetic vector potential, and $\bm{J}=\nabla\times\bm{B}/\mu_{0}$ is the
current density, with $\mu_{0}$ the vacuum permeability. The mass density
$\rho$ and the magnetic diffusivity $\eta$ remain constant, the scalar
potential $\Phi=0$ (Weyl gauge), and the viscous force
$\bm{F}_{visc}=\bm{\nabla}\cdot 2\nu\rho\bm{S}$, where $\nu$ is the kinematic
viscosity, and $\bm{S}$ is traceless rate of strain tensor with components
$\bm{S}_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})-\frac{1}{3}\delta_{ij}\bm{\nabla}\cdot\bm{u}$
(commas denote partial derivatives).

A 2.5 simulation is similar to 2D, except the evolving characteristic field
vectors have all three component along $x$, $y$ and $z$. For instance,
$\bm{B}=\nabla\times\bm{A}(x,y)=[B_{x}(x,y),\,B_{y}(x,y),\,B_{z}(x,y)]$.

Helicity | Resolution | $\eta\times 10^{-4}$ | $u_{rms0}$ | $B_{rms0}$ | $S/10^{3}$
---|---|---|---|---|---
nonhel | $1024^{2}$ | 1.3 | 0. | 0.4 | 0.8
nonhel | $2048^{2}$ | 0.34 | 0. | 0.4 | 3.
nonhel | $4096^{2}$ | 0.168 | 0. | 0.4 | 6.
hel | $1024^{2}$ | 1.3 | 0. | 0.4 | 0.8
hel | $2048^{2}$ | 0.34 | 0. | 0.4 | 3.
hel | $4096^{2}$ | 0.168 | 0. | 0.4 | 6.
nonhel | $512^{3}$ | 5.4 | 0. | 0.4 | 0.2
nonhel | $1024^{3}$ | 2.5 | 0. | 0.4 | 0.4
nonhel | $1024^{3}$ | 1.3 | 0. | 0.4 | 0.8
hel | $512^{3}$ | 5.4 | 0. | 0.4 | 0.2
hel | $1024^{3}$ | 2.5 | 0. | 0.4 | 0.4
hel | $1024^{3}$ | 1.3 | 0. | 0.4 | 0.8
Table 1: A summary of all the runs is listed. Helicity refers to whether the
initial random magnetic field is fully helical (given by ’hel’) or nonhelical
(given by ’nonhel’).

The magnetic field is initialized in the Fourier domain. We follow the same
methodology as in BZL21 (the details can be found in there). The initial field
is of random Gaussian nature with a forward slope of $k^{4}$, and peaked at
small scales or large wavenumbers corresponding to the peak wavenumber
$k_{p}\sim\frac{1}{L_{int}}=25$, where $L_{int}$ is the integral length scale
calculated as,

$L_{int}=\frac{\int(2\pi/k)M(k)dk}{\int M(k)dk},$ (4)

where $M(k)$ is the shell averaged 1D magnetic power spectrum. The velocity
field is initialized to be zero. Also $S=(2\pi k_{p}^{-1})\sqrt{2E_{M}}/\eta$
is the Lundquist number, where $E_{M}$ is total magnetic energy. The initial
parameters and configuration of the numerical simulations presented in this
paper are shown in Table 1.

##  3 Results

###  3.1 Reconnection in helical decaying MHD turbulence

Figure 1: Upper panel: evolution curves for $B^{H}_{rms}$ and $u_{rms}$ vs t
from fully helical 2.5D runs with different values for $S$. Lower panel:
collapse of the $B^{H}_{rms}$ evolution curves when plotted against
$t/\tau_{rec}$.

Figure 2: Upper panel: evolution curves for $B^{H}_{rms}$ and $u_{rms}$ vs t
from fully helical 3D runs with different values for $S$. Lower panel:
collapse of the $B^{H}_{rms}$ evolution curves when plotted against
$t/\tau_{rec}$.

It has been well known that in a decaying turbulent system, if the magnetic
fields are helical, then they tend to relax to larger scales. This results in
an inverse cascade of magnetic helicity accompanied by a simultaneous inverse
transfer of magnetic energy (Pouquet et al., 1976). The physical picture
typically discussed involves twists or writhes in the field, relaxing to
larger scales over a period of time. In this context, motivated by our
previous results in BZL21, we explore if magnetic reconnection has a role to
play in helical decaying MHD turbulence as well.

In Figs. 1 and 2, we show results from both 2.5D and 3D simulations. These
were initialized with zero velocity field and random fully helical magnetic
fields peaked at small scales which then decays as the simulations progress.
In particular, we show the evolution curves of the rms value of the helical
part of the magnetic field, defined as,

$B^{H}_{\rm rms}=\sqrt{2E_{M}^{H}},\quad E_{M}^{H}=\frac{1}{2}\int kH(k)dk,$
(5)

where $H(k)$ is the shell-averaged 1D spectrum of the magnetic helicity. In
both the plots for 2.5D and 3D cases in Figs. 1 and 2, the difference between
the top and bottom panels is that in the latter, the time axis is normalized
by the magnetic reconnection time-scale $\tau_{rec}=S^{1/2}\tau_{A}$, where
$\tau_{A}=(2\pi k_{p}^{-1})/\sqrt{2E_{M}})$ is the Alfven time-scale. We
observe that the $B^{H}_{\rm rms}$ evolution curves from runs with different
Lundquist numbers collapse on top of each other. This reveals that the
relevant dynamical time-scale is indeed the reconnection time-scale. This kind
of collapse was first seen in Figs. 10 and 11 in BZL21, in the nonhelical
decaying MHD turbulence case. Although the collapse in 3D is not as smooth as
that in the 2.5D case, it can be inferred that the curves would asymptotically
collapse at significantly higher Lundquist numbers. MHD simulations in 3D at
Lundquist numbers larger than $800$, i.e. at resolutions higher than
$1024^{3}$, could not be conducted because of HPC resource constraints. The
collapsing curves form the first part of the evidence for reconnection playing
a role in decaying helical MHD turbulence. Next we study the slope of the
decay curve.

One can calculate how the magnetic energy scales with time using simple
considerations of power law decay (this is a standard method in hydrodynamics
to derive decay scalings (Kolmogorov, 1941)),

$\frac{dE_{M}}{dt}\sim-\frac{E_{M}}{\tau_{\rm rec}}$ (6)

and the conservation of magnetic helicity, $\int\bm{A}\cdot\bm{B}~{}dV$.
Typically, the system under consideration is initialized to have a spectrum
that rises as $k^{4}$, peaks at large wavenumbers, and falls quickly. Such a
system, as it decays, retains this type of spectral shape with the peak
shifting to the smaller wavenumbers. This indicates that the system, at any
given time, has statistically more structures pertaining to the integral
length scale or scales close to that corresponding to the peak wavenumber,
$k_{p}$. This allows us to simplify and approximate our mathematical
descriptions to follow the dynamics of the fields concentrated at $k_{p}$.
With this in mind, one can simplify the conservation of magnetic helicity as,

$E_{M}k_{p}^{-1}\sim const.$ (7)

Thus, with Eqs. (6) and (7), we obtain $E_{M}\sim t^{-4/7}$ scaling. This
scaling using $\tau_{rec}$ was recently also derived by Hosking & Schekochihin
(2021). We find that this scaling bears out in both 2.5D and 3D simulations as
is seen in Figs. 1 and 2. This further strengthens the idea that magnetic
reconnection drives the decay and leads to the inverse transfer of energy in
this system.

An important point deserving of attention is that the scaling of magnetic
energy with time is the same in both 2.5D and 3D, a fact that has not been
adequately acknowledged. This cannot be a coincidence, given that such a
similarity was also found in the case of nonhelical decaying MHD turbulence in
BZL21. In Fig. 3, we show the evolution curves from both 2.5D and 3D
simulations together to clearly see that the slopes match. The underpinnings
of this similarity between 2.5D and 3D will be discussed in detail in a later
section 3.2.4.

Figure 3: Comparison of evolution curves for $B^{H}_{rms}$ from helical 2.5D
and 3D runs at the same value of $S=800$.

Figure 4: Evolution of $\langle A_{x}^{2}$ \+ $A_{y}^{2}\rangle$ and $\langle
A_{z}^{2}\rangle$ from 2.5D helical run at ${1024}^{2}$ resolution with
$S=800$.

Figure 5: Evolution of $B_{z}$ on 3D domain from 3D helical ($1024^{3}$,
$S=800$) run at times $t=2,12$ and $46$ from left to right.

The reason we explicitly calculated the helical part of the $B_{\rm rms}$ is
that in the 2.5D case, the nonhelical part has a nontrivial evolution curve
and thus total $B_{\rm rms}$ is not dominantly helical throughout the
simulation. In Fig. 3, we show the evolution curves (dashed lines) for the
nonhelical part of the field, $B_{\rm rms}^{NH}=\sqrt{2E_{M}^{NH}}$ computed
by subtracting out the helical part of the energy from the total,
$E_{M}^{NH}=E_{M}-E_{M}^{H}$. It can be seen that, in the 3D case, after an
initial rise around $t\sim 1$, $B_{\rm rms}^{NH}$ gradually decreases until
the end of the simulation. In contrast, in the 2D case, $B_{\rm rms}^{NH}$
shows growth over time. This growth is because even though vector-potential
squared, $\int\langle A_{z}^{2}\rangle dV$ is a conserved quantity in 2D,
since it is actually a 2.5D simulation (and not strictly 2D), the other two
components of the vector potential corresponding to $B_{z}$ are allowed to
grow as is shown in Fig. 4. However, from anti-dynamo theorems, we expect this
growth to cease and $B_{z}$ should also decay ultimately (Shukurov &
Subramanian, 2021). Nonetheless, it is interesting that in the 2.5D case,
nonhelical energy is produced in a system trying to achieve relaxation via
reconnection. So, in fact, the total field decays at a much slower rate. In
the 3D case, nonhelical energy is not produced as much. Thus, the decay of the
total field reflects mainly the evolution of the helical component, which, as
we have mentioned earlier, has the same scaling in time as seen in the 2.5D
case.

In the 2.5D case, the contour plots of the vector potential $A_{z}$ strung
into a movie show the merging of the magnetic islands with time. Localized
current sheets form between two merging islands, leading to magnetic
reconnection and, thus, larger islands. Similarly, plots of magnetic
structures on a 3D domain at a range of time instances for each magnetic field
component $[B_{x},B_{y},B_{z}]$ from helical 3D simulations are strung into a
movie. This movie shows initial small-scale magnetic structures, moving with
the turbulent flow merging with the localized neighbouring like-signed
structures, thereby forming larger structures with the evolution of time as
seen in Fig. 5. The structures merge likely via magnetic reconnection, leading
to an inverse transfer of magnetic energy.

###  3.2 The case of nonhelical decaying turbulence

Figure 6: Decay rates from 3D helical (upper panel) and nonhelical runs (lower
panel). The integral $I_{H}$ was evaluated using both methods of box-counting
(BC) and correlation-integral (CI) evaluated in a spherical (sph) region with
radius $R=0.1L$. The decay rate of $\langle\bm{A}^{2}\rangle$ is lower than
that of $I_{H}$ in the nonhelical case, indicating that former is the better
conserved quantity.

In BZL21, the study of 3D nonhelical decaying turbulence to investigate the
possibility of magnetic reconnection driving the system was motivated by the
similarity of the numerical result of $E_{M}\sim t^{-1}$ between 2D and 3D. In
the 2D case, it is clear that the system is driven by reconnection.
Additionally, vector-potential squared is not only ideally conserved but also
better conserved than magnetic energy in the limit of large Lundquist numbers.
One can do a similar calculation as we did in section 3.1 using Eq. (6) and
the constraint obtained from the conservation of vector-potential squared or
anastrophy,

$E_{M}k_{p}^{-2}\sim const.,$ (8)

to obtain the scaling law, $E_{M}\sim t^{-1}$.

In BZL21 it was claimed that conservation of anastrophy in Eq. (8), holds in
3D as well. To this end, BZL21 showed numerically that vector-potential
squared is better conserved than magnetic energy in 3D. However, this result
was overlooked in Hosking & Schekochihin (2021) with the claim, "this will be
true for any decay satisfying $B^{\alpha}L\sim const.$ for any $\alpha>0$",
where $L$ is a proxy for $1/k_{p}$ or $L_{int}$. Interestingly, vector
potential squared does not decay at all in the case of helical decaying MHD
turbulence, which corresponds to the case of $\alpha=2$. In the helical case,
$E_{M}\sim t^{-4/7}$ and thus $E_{M}k_{p}^{-2}\sim t^{4/7}$. Indeed, we find
that, both in 2.5D and 3D helical runs, the vector-potential squared grows as
a power law.

In Hosking & Schekochihin (2021), they instead introduced a new quantity,
termed the Saffman helicity invariant,

$I_{H}=\int\langle h(\bm{x})h(\bm{x}+\mbox{\boldmath$r$}{})\rangle
d\mbox{\boldmath$r$}{},$ (9)

where $\langle\rangle$ represents an ensemble average and claimed that this is
the invariant which governs the the scaling laws in the nonhelical case. Their
idea here is that a realistic system of MHD turbulence will contain helical
fields of both polarities, resulting in zero net helicity. The integral
$I_{H}$ which measures the correlation in the magnetic helicity density, is
expected to be ideally conserved and drive the nonhelical decaying system at
large Lundquist numbers. Using this, they obtain the following power law,
$E_{M}\sim t^{-1.18}$, which is actually quite different from $t^{-1}$
measured in past numerical simulations (Zrake, 2014; Berera & Linkmann, 2014;
Brandenburg et al., 2015; Reppin & Banerjee, 2017; Bhat et al., 2021).

We have multiple reasons to conclude that $I_{H}$ is not the underlying
integral constraining nonhelical turbulence decay. We delineate these below.

####  3.2.1 Decay of the invariants/conserved quantities

Figure 7: Evolution of the Saffman integral $I_{H}$ is evaluated by box-
counting (BC) and correlation-integral (CI) methods in a spherical region with
radius $R=0.1L$ in the 3D helical run with $1024^{3}$ resolution and $S=800$.

In the lower panel of Fig. 6, we show the rate of decay of magnetic energy,
vector potential squared $\langle\bm{A}^{2}\rangle$, $I_{H}$ and the total
magnetic helicity, $H_{M}$. We find that while $H_{M}$ is the slowest to
change as expected, the next slowest is, in fact, $\langle\bm{A}^{2}\rangle$.
Here, we have computed $I_{H}$ using the same algorithms as provided in Zhou
et al. (2022), available in Pencil-Code. In both methods of computing $I_{H}$
at $R=0.1L$, its rate of change is initially smaller than that of magnetic
energy but larger compared to $\langle\bm{A}^{2}\rangle$. Eventually, the rate
of change of $I_{H}$ is the same as that of the magnetic energy, and this is
likely because $k_{p}^{-1}\gtrsim R$. Thus, we find that $I_{H}$ is not much
slower than magnetic energy and thus not as robustly conserved as
$\langle\bm{A}^{2}\rangle$. We would like to further critically analyze the
role of $I_{H}$ in decaying MHD turbulence, which we do in the next
subsection. As expected, in the upper panel of Fig. 6, $H_{M}$ has the lowest
decay rate. Notably, $\langle\bm{A}^{2}\rangle$ here actually grows instead of
decaying (as we previously mentioned).

Another interesting thing to note is that in the helical case itself, $I_{H}$
does not decay, but grows with time as can be seen in Fig. 7. This reveals
that the integral $I_{H}$ encodes the correlation scale of magnetic helicity.
As the inverse transfer occurs, the integral scale shifts to larger values,
indicating growth of helical structures to larger scales, which is also
reflected in $I_{H}$.

####  3.2.2 Strength of helicity

Figure 8: Scalings from 3D ($1024^{3}$, $S=800$) runs. Upper panel: helical
case, lower panel: nonhelical case. The ratio of
$B_{rms}^{4}L_{int}^{5}/I_{H}$ (or equivalently, $E_{M}^{2}k_{p}^{-5}/I_{H}$)
$\gg 1$ indicating low strength in magnetic helicity in individual structures
of opposite polarities.

The integral in Eq. (9) can be calculated in two different ways (Zhou et al.,
2022). The one that is used in study of its invariance properties in Fig. 5 in
Hosking & Schekochihin (2021) is to assume ensemble averages are equivalent to
volume averages, where the averages are performed over volumes of size $R$,

$\displaystyle I_{H}(R)$ $\displaystyle=$
$\displaystyle\frac{1}{V_{R}}\left(\int_{V_{R}}d^{3}r\int_{V_{R}}d^{3}x~{}h(x)h(x+r)\right)$
(10) $\displaystyle\equiv$ $\displaystyle\frac{1}{VV_{R}}\int
d^{3}x\left(\int_{V_{R}}d^{3}r~{}h(x+r)\right)^{2}$

The above is called the box counting method in Zhou et al. (2022). The second
way to calculate the integral is to retain the inner correlation integral,

$I_{H}(R)=\frac{1}{V}\left(\int d^{3}r\int_{V_{R}}d^{3}x~{}h(x)h(x+r)\right).$
(11)

The way $I_{H}(R,t)$ varies with $R$ or $t$ in these two ways of calculation
is similar. This is because $R$ is restricted to the range, $L_{int}\ll R\ll
L$, where $L$ is the system size.

An important difference between magnetic helicity or vector-potential squared
integrals and the $I_{H}$ integral is that in the former case, the magnetic
structures are assumed to be mostly limited to a narrow range in wavenumber
$k\sim k_{p}$ (here $1/k_{p}\sim L_{int}$). Hence, the integral over the
volume is superficial, and thus, the scaling does not include a $L_{int}^{3}$
factor. However, in the latter case, the integral over the volume is necessary
to beget the scaling of $E_{M}^{2}k_{p}^{-5}\sim const.$. This complicates the
scaling relation as a simple dependence on a set of scales close to integral
scale or the peak wavenumber $k_{p}$ is not well justified.

However, our main criticism regarding the $I_{H}$ integral is related to the
strength of the helicity in the smaller structures. The strength of helicity
is related to whether the given magnetic field structure is fully helical. The
realizability condition requires, $k|H(k)|/2\leq M(k)$. When $k|H(k)|/2=M(k)$,
the magnetic fields on that scale or wavenumber are fully helical. In helical
decaying MHD turbulence, the reason we can apply the constraint in Eq. (7) is
because the fields in the system are fully helical to begin with, which is the
case in all the helical simulations in our paper as well as the previous ones.
In fact, during the early part of the evolution, the fields remain fully
helical, as we have seen already in Fig. 3, which shows that even though the
nonhelical part grows, it is smaller compared to $B_{\rm rms}^{H}$ by more
than an order of magnitude. However, the gap between the two curves keeps
decreasing monotonically.

Thus, in the fully helical case, we can indeed assume $\bm{A}\cdot\bm{B}\sim
E_{M}L_{int}\sim const.$ is a good constraint. Similarly, in the nonhelical
case, $I_{H}\sim E_{M}^{2}L_{int}^{5}\sim$ const. could be considered to
robustly constrain the system evolution only if the individual magnetic
structures of opposite polarities are fully helical. Our simulations show in
the lower panel of Fig. 8 that indeed the strength of helicity is expectedly
small in nonhelical turbulence and thus $E_{M}^{2}k_{p}^{-5}/I_{H}\gg 1$ (or
$B_{rms}^{4}L_{int}^{5}/I_{H}\gg 1$). This ratio (here $I_{H}$ is measured at
$R=0.1L$) remains larger than an order of magnitude (if we consider the dashed
curve) compared to a factor of $\sim 2$ in the case of the ratio
$B_{rms}L_{int}/A_{rms}$. In the helical run as shown in the upper panel of
Fig. 8, the ratio $B_{rms}^{2}L_{int}/H_{M}\sim 2\pi$ indicating the validity
of usage of $H_{M}$ to constrain the system. Note that the curves
corresponding to the ratios of $B_{rms}^{2}L_{int}/H_{M}$ (in upper panel) and
$B_{rms}L_{int}/A_{rms}$ (lower panel) have initial amplitude around $6$ due
to the $2\pi$ factor in $L_{int}$. In the case of $B_{rms}^{4}L_{int}^{5}\gg
I_{H}$, the volume cancels out and only $(2\pi)^{2}$ remains (due to squaring
of helicity), which should have led to an initial amplitude of $~{}40$ if the
structures were fully helical, but we get $100$ instead. In fact the ratio
involving magnetic helicity, $B_{rms}^{2}L_{int}/H_{M}$, is as large as $\sim
1000$ initially given that the system is nonhelical.

Figure 9: Histograms of the cosine angle between the vector potential and the
magnetic field from 3D ($1024^{3}$, $S=800$) runs. Upper panel: is the helical
case showing a peak at $+1$ with frequency of $\sim 0.1$ and a sharp drop from
there indicating that most of domain has fully helical fields (at all times).
Lower panel: is the nonhelical case which is flat at $t=0$. Peaks at values of
$\pm 1$ show up subsequently but with an incremental factor of only $1.18$
indicating most of the domain is still largely fully nonhelical.

Our argument gains additional support from Fig. 9. Here we present histograms
of the cosine angle between the vector potential and the magnetic field,
offering insights into their alignment. For fully helical fields, we obtain
$\cos{\theta}=1$ and for fully nonhelical fields, $\cos{\theta}=0$. In the
helical case of decaying MHD turbulence, as shown in the upper panel of Fig.
9, the magnetic fields are initialized to be random and fully helical. As
expected, nearly 90% of the total grid points (or counts) fall within the bins
of $\cos{\theta}\in[0.8,1]$, and a steep drop from the peak is observed.
Consequently, the portion of the domain containing nonhelical fields is
minimal, and this trend persists as the system evolves.

On the other hand, the lower panel in Fig. 9 shows the histogram for the
nonhelical case. In this case, the field is initially random but with no
helicity. Hence, the $t=0$ curve is flat. According to the theory put forward
by Hosking & Schekochihin (2021), the fully nonhelical fields should decay on
a fast time scale, leaving the system with only helical fields. This decay
should have manifested with peaks at both $\cos{\theta}=1$ and
$\cos{\theta}=-1$ and a sharp drop. However, we find that there is a
development of such peaks on a long time scale (compared to the Alfven scale),
and this growth in amplitude is only a factor of about 1.18, which is
insignificant. The growth of such a peak mainly happens because it is easier
for diffusion to kill fields which are less helical. Even at $t=16$, the
percentage of points which lie within the set
$\cos{\theta}\in[-1,-0.8]\cup[0.8,1]$ is only about 25%. Thus, we show that in
the nonhelical system, the strength of helicity is low. Therefore, a helicity-
based integral such as $I_{H}$ likely cannot constrain the evolution of the
system.

Another crucial point to underscore is that even though magnetic helicity
$H_{M}$ is well-conserved in also the nonhelical case (as evident in the lower
panel of Fig. 6), we cannot use it to constrain the evolution of nonhelical
system because the strength of helicity is insignificant. In a similar vein,
$I_{H}$ may function as a well-conserved quantity (with a slightly lower decay
rate compared to magnetic energy) (Zhou et al., 2022), but its effectiveness
in constraining the evolution of a nonhelical system remains questionable.

####  3.2.3 Strictly 2D case of nonhelical decay

The third reason for concluding that $I_{H}$ is not important in decaying
turbulence is related to the observation that, in both 2.5D and 3D, the
helical and nonhelical cases show the same, decay in time, power laws. In the
helical case, it is obviously the magnetic helicity conservation which plays
the critical role in constraining the system’s evolution in both 2.5D and 3D.
Similarly, in the nonhelical case, we expect the same conserved quantity (thus
leading to the same results) to constrain evolution in both 2.5D as well as
3D. A further interesting observation is that, in the nonhelical case, the
same decay power law of $E_{M}\sim t^{-1}$ holds in the strictly 2D case as
well ($B_{z}=A_{x}=A_{y}=0$ throughout the simulation). In the 2D case (as
opposed to 2.5D), magnetic helicity is not well-defined. Thus, any integral
based on magnetic helicity is ruled out, including the Saffman helicity
integral, $I_{H}$. Evolution plots from strictly 2D simulations with
$E_{M}\sim t^{-1}$ were already seen in Fig. 1 of BZL21. There, we also
confirmed that the strictly 2D simulations were similar to 2D simulations
using the reduced-MHD model in Zhou et al. (2019). Thus, only vector-potential
squared can explain the decay power law in the 2D case. And its influence
extends to the 3D case as well, similar to magnetic helicity in the helical
case.

####  3.2.4 Local anisotropy

The question, then, is why there is a match between 2D (or 2.5D) and 3D
results. Such a question was already posed in BZL21. The 2D-3D match has
turned out to be more general, extending to the helical case as well. We
continue to think that a quasi-two-dimensionalization is involved, in the
sense that the individual magnetic reconnection events seem to manifest in a
2D fashion within the 3D domain. In the reconnection literature, a quasi-two-
dimensionalization manifests typically in the presence of a guide field, i.e.
the reconnection rates are similar between a 2D case and a 3D case with a
guide field. We resort to the well-known idea that local mean fields in the
system can provide a local guide field for small-scale fields leading to local
anisotropy (Cho & Vishniac, 2000). We study such possibilities using Minkowski
functionals in the next section.

##  4 Quantifying magnetic field anisotropy using Minkowski Functionals

Minkowski functionals (MFs) are tools which can be used to find the morphology
of structures in a scalar field. It can help identify whether the isosurfaces
in a given dataset are planar or filamentary, etc. The system of interest in
this paper is magnetically dominated decaying MHD turbulence, which does not
contain a guide field. Thus, it is not obvious if there can be anisotropy in
the system, at least locally. But from our analysis in previous sections, we
have provided ample evidence which suggests the possibility of quasi-two-
dimensionalization of the 3D system at the relevant scales where magnetic
reconnections occur. This likely includes the integral length scale (or the
scales corresponding to the peak wavenumber $k_{p}$) and the smaller scales.
Since we posit that reconnection in the 3D system manifests in a 2D fashion,
from the previous literature, we gather that the reconnecting structures need
to be influenced by a local guide field or mean field. Depending on the
physics and the presence of local guide/mean fields, the magnetic structures
can be rendered anisotropic in their morphology. We would like to use
Minkowski functionals to detect and measure this anisotropy.

We can test for quasi-two-dimensionalization by quantifying the spatial
correlation of magnetic fields near reconnection sites in the computational
domain. In such a quasi-2D scenario, the magnetic fields are expected to be
correlated preferentially along one direction over the other two. That is, the
fields should be anisotropic locally around the reconnection sites.

Figure 10: Known as the Blaschke diagram, we show the morphology of different
spheroids as values of $P$ and $F$ varies along the abscissa and the ordinate
respectively.

We first identify reconnection sites and then detect the magnetic fields
around these sites. The necessary condition for reconnection to occur is the
presence of high current density (see Priest (2016) for sufficient
conditions). Around a given high current density site, we analyze the
correlation of magnetic fields by assessing the morphology of isonorm surfaces
of $|\bm{B}|$. The information regarding the shape of the isonorm surfaces is
extracted using Minkowski Functionals (MFs). Mecke et al. (1994) first
introduced MFs in cosmology, and these were later employed by many others for
morphological feature extraction (Refer to Sheth et al. (2003) for a complete
review of MFs). The execution of Minkowski functionals we employ is in Python
from Sheth et al. (2003), who use the shape finders introduced by Sahni et al.
(1998). The four Minkowski functionals are defined as,

$\begin{split}V_{0}=\iiint dV&\quad V_{1}=\frac{1}{6}\iint dS\\\
V_{2}=\frac{1}{6\pi}\iint\left(\kappa_{1}+\kappa_{2}\right)dS&\quad
V_{3}=\frac{1}{4\pi}\iint\kappa_{1}\kappa_{2}\,dS\end{split}$ (12)

where $\kappa_{1}$ and $\kappa_{2}$ are the two principal curvatures. The
volume and surface integrals, respectively, in the equations above, are
specified over closed isosurfaces for a given field. We can extract three
length scales from $V_{0},V_{1},V_{2}$ and $V_{3}$. The typical thickness,
width, and length derived from (12) are,

$T=\frac{V_{0}}{2V_{1}},\quad W=\frac{2V_{1}}{\pi V_{2}},\quad
L=\frac{3V_{2}}{4V_{3}}.$ (13)

These three length scales are ordered such that
$l_{1}=\text{max}(T,W,L),l_{2}=\text{mid}(T,W,L),l_{3}=\text{min}(T,W,L)$.
Next, we can obtain planarity and filamentarity for a given isosurface
structure as follows,

$P=\frac{l_{2}-l_{3}}{l_{2}+l_{3}},\quad F=\frac{l_{1}-l_{2}}{l_{1}+l_{2}}.$
(14)

To both validate our implementation of the MFs and provide an understanding of
the morphology of structures in a phase space governed by $P$ and $F$, we show
the Blaschke diagram in Fig. 10. Here, the abscissa varies over planarity $P$,
and the ordinate varies over filamentarity $F$. It is clear from Eq. (14) that
the range of $P$ and $F$ is $[0,1]$, where zero means no planarity or no
filamentarity. In Fig. 10, we show how the shapes change depending on the
values of ($P,F$). We have a sphere in the lower left corner of the plot that
changes to a plate as $P$ increases (with $F$ constant) in the lower row but
changes to a rugby ball-like shape as $F$ increases (with $P$ constant) in the
left column.

Further, we have also performed the MF analysis separately for fully developed
hydrodynamic turbulence and confirmed the expected results. This can be found
in appendix D.

Since magnetic reconnection occurs near points where there is a significant
value of current density, we implement an algorithm to obtain isolated
clusters of points with the current density value $|\bm{J}|$ exceeding a
threshold of $0.7|\bm{J}|_{\text{max}}$. In each such cluster, we designate
the point with the highest value as the representative point. We designate the
current density isosurface with isovalue of $0.5|\bm{J}|_{\text{max,
cluster}}$ as our fiducial current sheet for that cluster.

Figure 11: Zoomed in version of isosurfaces of the current density $|\bm{J}|$
and the magnetic field $|\bm{B}|$ are shown in green and red respectively from
the nonhelical run of $1024^{3}$ with $S=800$. The domain size shown here is
$72^{3}$. Figure 12: Zoomed in version of isosurfaces of the current density
$|\bm{J}|$ and the magnetic field $|\bm{B}|$ are shown in green and red,
respectively, from the fully helical run of $1024^{3}$ with $S=800$. The
domain size shown here is $116^{3}$.

Figure 13: Nonhelical: Histogram for planarity (left panel) and filamentarity
(right panel) of all the reconnecting magnetic field structures in the
nonhelical $1024^{3}$ simulation with $S=800$, at $t=5$. Also shown in the
black squares is the average volume in each bin. NVC stands for no-volume-
cutoff. The darker shades represent statistics for structures considered above
a certain volume threshold denoted by $\overline{V}/n$ and $n=1$, $5$ or $10$.
The tendency of large structures to have smaller planarity and larger
filamentarity is clear from the graph. The same conclusion can also be reached
by progressively considering structures with larger and larger volumes, as we
have shown in this plot.

Figure 14: Nonhelical: Time evolution of $P$ and $F$ in the early stages of
the decay (when the $B_{rms}$ curve behaviour is close to the theoretical
prediction of $t^{-0.5}$), is shown in the left panel. The time evolution of
length scales of reconnecting magnetic field isosurfaces is shown in the right
panel. The different dashed curves are the least square fit, and they indicate
a scaling close to the expected $l\sim t^{1/2}$.

Next, we would like to obtain a set of magnetic field isosurfaces close to
this current sheet. So, we employ a range of isovalues to do this. To
determine this range of isovalues, consider the following scheme. Around a
current sheet (in an isolated cluster), $|\bm{B}|$ is calculated along rays
emanating from the representative point. Once we obtain the $|\bm{B}|$ profile
around the representative point, we find the maximum and minimum value of
$|\bm{B}|$ nearest to the point. Then we use these extrema values to determine
a range of isovalues of $|\bm{B}|$ given by,

$|\bm{B}|_{\text{max}}-\frac{i}{5}\times\left(\frac{|\bm{B}|_{\text{max}}-|\bm{B}|_{\text{min}}}{4}\right)\\\
\forall i\in{1,2,3,4,5}.$ (15)

Thereafter, Minkowski functionals for $|\bm{J}|$ and $|\bm{B}|$ surfaces are
calculated, and their evolution with time is reported. The details of
algorithm for MFs implementation and the scheme we described briefly for
obtaining isosurfaces is available in Appendix A and B.

###  4.1 Results for Minkowski functionals

For an illustration, we show in Figs. 11 and 12 the isosurfaces of an isolated
current density structure corresponding to a possible reconnection site and
also isosurfaces of the magnetic fields nearby that are potentially
participating in the reconnection. We calculate the MFs for such structures of
current density and the magnetic field throughout the domain of the $1024^{3}$
runs of $S=800$ in both helical and nonhelical cases. We extract the
quantities $l_{1}$, $l_{2}$, $l_{3}$ and the corresponding $P$ and $F$ values.
Below, we study their distributions (histograms) and evolution in time.

####  4.1.1 Nonhelical case

In Fig. 13, we show the histograms of the measured planarity $P$ and
filamentarity $F$ at a given time, $t=5$ from the $1024^{3}$ run of $S=800$.
In these charts, NVC (no-volume-cutoff) denotes the histogram for the complete
set of structures. The darker shades correspond to histograms derived from a
subset where only structures with a volume surpassing a specific threshold are
considered. The darkest shade represents the subset where only structures with
a volume exceeding the average value $\overline{V}$ are considered. The
planarity histogram (left panel of Fig. 13) peaks within the bin of
$0.1$–$0.2$. With subsets of larger volume structures, the histograms
progressively shrink towards lower values of $P$. This indicates the tendency
of more voluminous structures to have smaller planarity. In the case of
filamentarity histogram (right panel of Fig. 13), the peak is within the bin
of $0.2$–$0.3$. Here, as we limit the set to having more voluminous
structures, the average value of $F$ shifts to the higher side. On the whole,
the structures are more filamentary than planar. These values indicate that
the structures feature between the sphere and the middle ellipsoid in the
first column of the Blaschke diagram in Fig. 10. Thus, their geometry is akin
to prolate spheroids, which are more filamentary than they are planar (unlike
the oblate spheroid). This kind of anisotropy inherent in the magnetic
structure is likely reflective of the local anisotropy within the system.

Next, we study the evolution of the morphology of the structures. In the left
panel of Fig. 14, we show the evolution of the average value of $P$ and $F$
calculated for the magnetic fields. Firstly, both $P$ and $F$ tend to
decrease. This can be understood as the increasing effect of diffusion, which
tends to even out the structures. While the decrease in $F$ is very small
($\sim t^{-0.1}$ given by a straight line fit as mentioned in the legend), the
decrease in $P$ is larger. Here, the $P$ and $F$ values were calculated by
taking the average of the values obtained from all the individual
isosurfaces(detailed steps are in Appendix B). These average $P$ and $F$
values do not have to necessarily match with that resulting from the average
values of $l_{1}$, $l_{2}$ and $l_{3}$ (let us call them $P^{*}$ and $F^{*}$),
and indeed they don’t. However, the trends match, though the indices of power
law decay are different, i.e. $P^{*}\sim t^{-0.25\pm 0.03}$ and $F^{*}\sim
t^{-0.2\pm 0.04}$.

Figure 15: Helical: Histogram for planarity (left panel) and filamentarity
(right panel) of all the reconnecting magnetic field structures in the helical
$1024^{3}$ simulation with $S=800$, at $t=5$. Also shown in the black squares
is the average volume in each bin. NVC stands for no-volume-cutoff. The darker
shades represent statistics for structures considered above a certain volume
threshold denoted by $\overline{V}/n$ and $n=1$, $5$ or $10$. The tendency of
large structures to have smaller planarity and larger filamentarity is clear
from the graph. The same conclusion can also be reached by progressively
considering structures with larger and larger volumes, as we have shown in
this plot.

Figure 16: Helical: Time evolution of $P$ and $F$ in the early stages of the
decay (when the $B_{rms}$ curve behaviour is close to the theoretical
prediction of $t^{-2/7}$), is shown in the left panel. The time evolution of
length scales of reconnecting magnetic field isosurfaces is shown in the right
panel. The different dashed curves are the least square fit, and they indicate
a scaling close to the expected $l\sim t^{4/7}$.

The behaviour of $P^{*}$/$P$ and $F^{*}$/$F$ is roughly consistent with the
power law growth of $l_{1}$, $l_{2}$ and $l_{3}$ as shown in the right panel
of Fig. 14. The thin dashed lines indicate a straight line fit to the
evolution curves, and we infer the respective power-law indices. It can be
seen that indeed $l_{2}$ and $l_{3}$ increase faster than $l_{1}$. This can be
understood as the effect of diffusion, which acts more efficiently on smaller
scales. Thus, the indices increase progressively as we go from $l_{1}$ to
$l_{3}$. Otherwise, the values of the power law indices (within the error
bars) are close to the theoretically expected value of 0.5 (more so for
$l_{1}$ and $l_{2}$ than $l_{3}$, which is most affected by diffusion). It’s
important to note that the accuracy of the straight line fit is influenced by
the chosen time range. Nevertheless, we demonstrate that a fit closely
aligning with the theoretical expectation is not unreasonable, considering the
estimated error bars.

####  4.1.2 Helical case

Next, we consider the iso-surfaces in the decaying helical MHD turbulence. In
Fig. 15, we once again show the histograms of the measured planarity $P$ and
filamentarity $F$ at a given time, $t=5$ from the helical $1024^{3}$ run of
$S=800$. The peak of the planarity histogram, similar to that of the
nonhelical case, peaks within the bin of $0.1-0.2$. However, the filamentarity
histogram peaks at a larger value within the bin of $0.3-0.4$ for more
voluminous structures. Here again, as we restrict to a subset of structures
with larger volumes, the averages shift to lower values for $P$ and higher
values for $F$. As in the nonhelical case, the structures are more filamentary
(even more so) than planar.

In the left panel of Fig. 16, we show the evolution of the average value of
$P$ and $F$ to assess how the morphology of the structures changes. Here, $F$
is nearly constant as time evolves, whereas the value of $P$ shows a decrease
(similar to the nonhelical case) with time. In the right panel of Fig. 16, we
show the evolution of $l_{1}$, $l_{2}$ and $l_{3}$ and we find here the decay
of $l_{2}$ and $l_{3}$ are not all that much faster than $l_{1}$ (unlike the
nonhelical case). We can chalk up this behaviour and constancy of $F$ to the
fact that helical structures tend to be more resilient to diffusion (Bhat et
al., 2014). Here again, we find that the power law indices from the straight
line fit, as shown are close to the theoretically expected value of $\sim
0.57$.

MF analysis was also performed for the current density structures. In Fig. 17,
we show the phase space Blaschke plot for $P$ and $F$ values from all the
different cases of helical MHD, nonhelical MHD and hydrodynamics (no magnetic
fields). Here, the values of ($P,F$) for the current density structures (shown
in a diamond shape) are around $\sim(0.53,0.5)$ in the helical case and
$\sim(0.45,0.37)$ in the nonhelical case. Thus, unlike the case of magnetic
field isosurfaces, the current density structures appear between the bottom
two structures in the middle column of the Blaschke diagram in Fig. 10. Thus,
they are largely sheet-like, reaffirming the understanding of this system as a
quasi-2D system with 2D-like reconnections.

##  5 A quasi-2D hierarchical merger model

Figure 17: Blaschke plot of the combined helical, nonhelical, and hydrodynamic
simulations. The $(P,F)$ of the vorticity of the forced hydrodynamic
simulation has been given as a reference point against which the $(P,F)$ of
the MHD simulations could be compared. Note that there is an outlier in the
$(P,F)$ of the magnetic field in the helical simulation at $(0.71,0.22)$ which
occurs at $t=8$. This is because there was only one current sheet above the
threshold, so we could not obtain enough statistics. We ignore this data point
in the evolution plots.

A hierarchical merger model was first proposed in Zhou et al. (2019) to
understand a 2D system is driven by island mergers due to coalescence
instability. It was useful to derive scaling laws associated with the decay of
magnetic energy and resulting inverse transfer. While BZL21 suggested that the
3D decaying MHD turbulence is akin to the 2D system (also with similar scaling
laws), but they didn’t provide a 3D model to describe the system. A 3D model
was proposed in Zhou et al. (2020), but it was for a system with a uniform
mean field or guide field. We find that their hierarchical model is
incompatible with what we propose below.

In our quasi-2D hierarchical merger model (Q2DHM model), interactions between
magnetic structures happen pairwise. An interaction consists of magnetic
reconnection at the interface, leading to a merger. Such a merger, then,
conserves mass and magnetic flux or magnetic helicity (depending on the
system). The resulting structures are then a size larger but with a lower
strength of magnetic field. Thus, the two ingredients which were required to
solve the Kolmogorov-type calculation using decay equation in Eq. (6) is
available in the Q2DHM model as well : (i) the reconnection timescale being
the driving timescale for the system, (ii) a conserved quantity to constrain
the relationship with the magnetic field and the integral scale in the
problem. However, the advantage of the Q2DHM model over the decay equation in
Eq. (6) is that it allows us to make the distinction between a 2D and a 3D
system and explicitly see the manifestation of quasi-two-dimensionality.

In the model, the mergers happen in discrete stages (given by $n$). A key
assumption is that at any given stage, there is uniformity in magnetic field
structures, i.e. all of them have the same dimensions of $l_{1n}$, $l_{2n}$
and $l_{3n}$. The magnetic field in a structure is given by $B_{n}$.

First, the merger of two structures conserves mass. Assuming the condition of
incompressibility, we can write,

$2l_{1_{n}}l_{2_{n}}l_{3_{n}}=l_{1_{n+1}}l_{2_{n+1}}l_{3_{n+1}}$ (16)

where $n+1$ denotes the next stage (or the next generation). Thus, we have
$l_{1_{n+1}}=2^{1/3}l_{1_{n}}$, $l_{2_{n+1}}=2^{1/3}l_{2_{n}}$ and
$l_{3_{n+1}}=2^{1/3}l_{3_{n}}$.

Second, we demand flux conservation. However, for a system with quasi-two-
dimensionality, the structures reconnect on a plane guided by a local mean
field, then the conservation of flux can be given by,

$B_{n}l_{i_{n}}\sim const.$ (17)

where the $i$ subscript can assume $1,2$ or $3$ depending on the plane of
reconnection. Thus, we have $B_{n+1}=2^{-1/3}B_{n}=2^{-n/3}B_{0}$, where
$B_{0}$ is the initial field. We can now calculate the Lundquist number at the
$n+1$ stage,

$S_{n+1}=\frac{B_{n+1}l_{i_{n+1}}}{\eta}=\frac{2^{-1/3}B_{n}2^{1/3}l_{i_{n}}}{\eta}=S_{n}$
(18)

Since each merger is driven by reconnection, the corresponding timescale is
given by the inverse of the reconnection rate. Since we assume the 2D steady-
state model of Sweet-Parker for reconnection, we obtain,

$\tau_{{rec}_{n+1}}=S_{n+1}^{1/2}\tau_{A_{n+1}}=2^{2/3}\tau_{A_{n}}=2^{2n/3}\tau_{0}$
(19)

The time taken to reach the $n$th generation can be approximated to be
$t_{n}\approx 2^{(2/3)n}\tau_{0}$, for large $n$. Thus we have time
$\tilde{t}=t_{n}/\tau_{0}=2^{(2/3)n}$. Now, we can eliminate $n$ to obtain
time dependence of various quantities,

$\displaystyle B=B_{0}\tilde{t}^{-1/2},\quad l_{1}=l_{1_{0}}\tilde{t}^{1/2}$
(20) $\displaystyle l_{2}=l_{2_{0}}\tilde{t}^{1/2},\quad
l_{3}=l_{3_{0}}\tilde{t}^{1/2}$

Thus, we recover the scalings which otherwise were previously derived using
equation Eq. (6). Note that if we do not assume the quasi-two-dimensionality
for Eq. (17) and instead used $B_{n}l_{i_{n}}l_{j_{n}}\sim const.$, we would
have obtained very different scalings.

In the helical case, instead of Eq. (17), the more important topological
constraint is derived from magnetic helicity. Thus, for each structure, we
consider

$B^{2}_{n}l_{i_{n}}\sim const.$ (21)

Here, then we have $B_{n+1}=2^{-1/6}B_{n}=2^{-n/6}B_{0}$. This will affect the
the scaling of $S$,

$S_{n+1}=\frac{B_{n+1}l_{i_{n+1}}}{\eta}=\frac{2^{-1/6}B_{n}2^{1/3}l_{i_{n}}}{\eta}=2^{1/6}S_{n}$
(22)

Consequently, instead of Eq. (19), we obtain,

$\tau_{{rec}_{n+1}}=S_{n+1}^{1/2}\tau_{A_{n+1}}=2^{7/12}\tau_{A_{n}}=2^{7n/12}\tau_{0}$
(23)

Note that for calculating $\tau_{rec}$ in both cases, we use the 2D Sweet-
Parker model. We can then, as before, eliminate $n$ to obtain time dependence
of various quantities,

$\displaystyle B=B_{0}\tilde{t}^{-2/7},\quad l_{1}=l_{1_{0}}\tilde{t}^{4/7}$
(24) $\displaystyle l_{2}=l_{2_{0}}\tilde{t}^{4/7},\quad
l_{3}=l_{3_{0}}\tilde{t}^{4/7}$

Thus, we recover the scalings observed in our simulations, which are also
consistent with the Kolmogorov-type calculation.

##  6 Discussions

##### Necessity of high Lundquist number and its connection to quasi-two-
dimensionality

The simulation results of inverse transfer in non-helical MHD decaying
turbulence posed a puzzle upon their emergence, given the absence of magnetic
helicity in this system (Zrake, 2014; Brandenburg et al., 2015). Notably, this
is a high magnetic Reynolds number or high $S$ (available only at a higher
resolution) result, which explains its non-detection in many of the earlier
simulations of nonhelical MHD decaying turbulence (Christensson et al., 2001).
Thus, a comprehensive understanding of this system should duly consider this
important factor. At low $S$ or in a low-resolution simulation, the timescales
of operation for microscopic diffusion can get sufficiently short such that
any effect of local anisotropy can be rendered ineffective. However, a more
significant issue arises when the occurrence of magnetic reconnection is
impeded at sites of large current density due to non-negligible diffusion
(aided by turbulence). One may inquire about the implications for systems that
are not non-helical. Similar considerations should be taken into account in
that scenario as well. Indeed, we observe a steeper decay in lower resolution
simulations, which does not necessarily align with the theoretical calculation
of $E_{M}\sim t^{-4/7}$. Thus, whether nonhelical or helical, in both cases,
the system at large resolution and high Reynolds numbers tends to be more
quasi-two-dimensional. Consequently, it is not surprising when the relevant
reconnection time-scale corresponds to a 2D model, i.e., Sweet-Parker
reconnection. Further, it is desirable to uncover the behaviour of
reconnection with the Prandtl number(Armua et al., 2023).

##### A uniquely MHD result

This quasi-two-dimensionalization leads to exact same decay laws in both 2D
and 3D, which is a uniquely MHD result. Indeed, in hydrodynamics, the 2D and
3D decay laws are different, with different conserved quantities at play. Such
a 2D-3D similarity is not to be expected in hydrodynamics unless there is an
agent for bringing about quasi-two-dimensionalization, which in the MHD case
is naturally manifested by large magnetic field strength (even if the fields
are of turbulent nature and reside at smaller scales). In hydrodynamic
turbulence, such an agent could be strong rotation, and then the relevant
conserved quantities will change back to that relevant in 2D as expected
(Yakhot & Pelz, 1987).

##### Local mean fields

We attribute the quasi-two-dimensionalization to the presence of non-zero
local mean magnetic fields. In the standard theories for steady-state
magnetically dominated MHD turbulence (either weak or strong), a mean-field is
required as the explicit source for anisotropy. The dynamics, then, are
described along the mean field separately from across it. Depending on whether
the turbulence is weak or strong, connections are made between the dynamics
perpendicular and parallel to the mean field (Tobias et al., 2011). In our
system, there is no apriori global mean field (or uniform field), however, we
postulate that as the strong random fields reconnect, a net local mean field
can arise to influence the dynamics locally. This concept aligns with the
conventional understanding in standard MHD turbulence theories that,
ultimately, the local mean field is important (despite the presence of a
uniform mean field) (Cho & Vishniac, 2000). However, in this study, we have
refrained from directly calculating the local mean fields. The reason is that
it is unclear what the scale size for local mean field’s region of influence
needs to be for it to render the reconnection quasi-2D.

##### 3D reconnection

Studies on 3D reconnection are sparse. Most studies either discuss various
kinds of manifestations of reconnection sites or arising current instabilities
(for e.g. kink instability (Pontin, 2011; Landi et al., 2008; Oishi et al.,
2015)). While it is well known that a guide field or mean-field renders the
ensuing reconnection in the system effectively 2D. Another manifestation of
guide field reconnection is when two flux tubes interact at an angle (Linton
et al., 2001). However, a complete and systematic characterization of such
systems is required for understanding within a larger context of turbulent
systems.

##### Anisotropy measurement

Previous studies that have measured anisotropy have used structure function
calculation, taking into account the local mean field in the calculation (Cho
& Vishniac, 2000). However, given the uncertainty in the required coherence
scale for the mean field, we are unable to use similar techniques. Instead, we
use Minkowski functionals to calculate the characteristic length scales of the
magnetic field structures. All the three length scales $l_{1},l_{2}$ and
$l_{3}$ we obtain from this analysis are very different, leading to non-zero
significant planarity and filamentarity. While the difference between $l_{1}$
and $l_{2}$ or $l_{3}$ can be attributed to local mean-field, that between
$l_{2}$ and $l_{3}$ could be due cross helicity considerations and dynamic
alignment (Boldyrev et al., 2009). We have not studied these in this work and
leave it for a future investigation.

##### Issue of gauge

Already, magnetic helicity is a tricky topological quantity to deal with
depending on the nature of the system under study, and there are many issues
to consider in its calculation. Depending on the boundary conditions and other
ingredients like stratification, large-scale shear, etc, there can be magnetic
helicity fluxes which can change the conservation properties dramatically.
Then there is the issue related to gauge, which relates to the meaning and
interpretation of the results in a system with nontrivial (which are less
academic and more realistic) boundaries. But as long as one is in a periodic
system, results can be considered to be reasonably robust.

Similar considerations can be extended to the quantity of vector-potential
squared. This quantity is physically connected to the concept of flux-
freezing. In 2D, magnetic flux conservation is directly related to the
invariance of anastrophy or vector-potential squared. In our decaying
nonhelical MHD turbulence case, anastrophy is particularly useful given that
the system displays quasi-two-dimensionalization. In the absence of mean
magnetic helicity, this quantity displays slower decay time scales compared to
magnetic energy and this is fairly nontrivial as, in the case of non-zero
magnetic helicity, anastrophy tends to increase with time. In other contexts
related to condensates in Yang-Mills theories, the minimization of this
quantity w.r.t the topology of the magnetic field, is associated with a gauge
condition which is invariant. This is used as a probe of phase transition
(Gubarev et al., 2001). Thus, we find that understanding the 3D system with
anastrophy is not unjustified on the grounds of gauge arguments as then those
would equally apply to a more commonly used 3D quantity like magnetic helicity
density.

##### 3D Hierarchical merger models

In our Quasi-2D Hierarchical Merger (Q2DHM) model, we posit that due to quasi-
two-dimensionalization, the magnetic flux conservation constraint is for the
in-plane magnetic field (w.r.t reconnection). This is consistent with the
requirement of vector-potential squared or anastrophy conservation we
discussed above. As mentioned previously, our model is not compatible with
that presented in Zhou et al. (2020) as our first constraint of mass
conservation leads to Eq. (16), which is not the case in theirs. Further, they
draw constraints from strong MHD turbulence of critical balance. The nature of
turbulence in this system is not yet fully understood due to the lack of
probes on the nature of the local mean fields. However, as mentioned in BZL21,
the slope of $k^{-2}$ in the magnetic spectrum is more compatible with that of
weak turbulence.

##  7 Conclusions

In a previous paper, we had shown that magnetic reconnection is important in
decaying nonhelical MHD turbulence. In this work, we showed that the same is
true of decaying fully helical MHD turbulence as well. Total magnetic
helicity, a conserved quantity in the limit of high Lundquist numbers, is used
to constrain the evolution of the helical case. We found that the predicted
theoretical power law decay exponent pans out in both 2.5D and 3D simulations.
We critically analyzed two possibilities of ideal invariants that could
constrain evolution in the nonhelical case, (i) anastrophy or vector potential
squared (which is ideally conserved in 2D or 2.5D), $\langle\bm{A}^{2}\rangle$
(ii) Saffman helicity integral $I_{H}$ (which is undefined in strict 2D). We
showed that $\langle\bm{A}^{2}\rangle$ has a lower decay rate as compared to
$I_{H}$ or the magnetic energy $\langle\bm{B}^{2}\rangle/2$. Next, we argued
that the assumption of the nonhelical system being a composite of fully
helical structures of opposite polarities (where the fully nonhelical part
decays away on fast timescales) is not valid. We provided evidence for this
with (i) a histogram of the cosine angle between vector potential and magnetic
field in Fig. 9, where the concentration at and around $\pm 1$ is low, and
(ii) we showed that the ration of $E_{M}^{2}k_{p}^{-5}/I_{H}\gg 1$. Our final
contention was that the same nonhelical behaviour is also observed in strict
2D case, where $I_{H}$ is basically undefined.

Next, we investigated the possibility of quasi-two-dimensionalization, which
can explain the similarity seen in power law decay exponents between 2D (or
2.5D) and 3D cases. To this end, we performed a Minkowski functional analysis
for the reconnecting magnetic fields. We found that both helical and
nonhelical fields have a planarity $P\sim 0.1$–$0.2$ and filamentarity $F\sim
0.2$–$0.3$, indicating the existence of local anisotropy. Also, we showed that
the fit to power law evolution curves of the characteristic scales, $l_{1}$,
$l_{2}$ and $l_{3}$ displays exponents that are close to the theoretically
expected values of $0.5$ and $0.57$ in nonhelical and helical cases
respectively.

We provided a Quasi-2D Hierarchical Merger (Q2DHM) model based on mass
conservation and either magnetic flux conservation (in the nonhelical case) or
helicity conservation (in the helical case). In this model, we recover the
observed scalings only when we consider (i) a 2D manifestation of flux-
freezing consistent with vector-potential squared conservation in the plane of
reconnection, (ii) and Sweet-Parker model (a 2D model) for reconnection.

## Acknowledgements

SD, CA and PB are thankful to Prof. Kandaswamy Subramanian for his insightful
comments and discussions. SD expresses gratitude to ICTS for hosting her in
the Long Term Visiting Students Program. CA thanks Vinay Kumar and Rajarshi
Chattopadhyay for their valuable help and discussions. We acknowledge project
RTI4001 of the Dept. of Atomic Energy, Govt. of India. SD acknowledges support
from the US Dept. of Energy (grant number DE-SC0018266). The simulations were
performed on the ICTS HPC cluster Contra.

## References

  * Armua et al. (2023) Armua A., Berera A., Calderón-Figueroa J., 2023, Physical Review E, 107, 055206
  * Batchelor & Taylor (1949) Batchelor G. K., Taylor G. I., 1949, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 195, 513
  * Berera & Linkmann (2014) Berera A., Linkmann M., 2014, Physical Review E, 90, 041003
  * Bhat & Subramanian (2013) Bhat P., Subramanian K., 2013, Monthly Notices of the Royal Astronomical Society, 429, 2469
  * Bhat et al. (2014) Bhat P., Blackman E. G., Subramanian K., 2014, Monthly Notices of the Royal Astronomical Society, 438, 2954
  * Bhat et al. (2017) Bhat P., Ebrahimi F., Blackman E. G., Subramanian K., 2017, Monthly Notices of the Royal Astronomical Society, 472, 2569
  * Bhat et al. (2021) Bhat P., Zhou M., Loureiro N. F., 2021, Monthly Notices of the Royal Astronomical Society, 501, 3074
  * Blandford et al. (2017) Blandford R., Yuan Y., Hoshino M., Sironi L., 2017, Space Science Reviews, 207, 291
  * Boldyrev et al. (2009) Boldyrev S., Mason J., Cattaneo F., 2009, The Astrophysical Journal, 699, L39
  * Brandenburg et al. (2015) Brandenburg A., Kahniashvili T., Tevzadze A. G., 2015, Physical Review Letters, 114, 075001
  * Cho & Vishniac (2000) Cho J., Vishniac E. T., 2000, The Astrophysical Journal, 539, 273
  * Christensson et al. (2001) Christensson M., Hindmarsh M., Brandenburg A., 2001, Physical Review E, 64, 056405
  * Daughton et al. (2011) Daughton W., Roytershteyn V., Karimabadi H., Yin L., Albright B. J., Bergen B., Bowers K. J., 2011, Nature Physics, 7, 539
  * Gubarev et al. (2001) Gubarev F. V., Stodolsky L., Zakharov V. I., 2001, Physical Review Letters, 86, 2220
  * Hesse & Cassak (2020) Hesse M., Cassak P. A., 2020, Journal of Geophysical Research: Space Physics, 125, e2018JA025935
  * Hosking & Schekochihin (2021) Hosking D. N., Schekochihin A. A., 2021, Physical Review X, 11, 041005
  * Hosking & Schekochihin (2023) Hosking D. N., Schekochihin A. A., 2023, Nature Communications, 14, 7523
  * Kolmogorov (1941) Kolmogorov A. N., 1941, Akademiia Nauk SSSR Doklady, 32, 16
  * Landau & Lifshitz (1975) Landau L. D., Lifshitz E. M., 1975, The classical theory of fields
  * Landi et al. (2008) Landi S., Londrillo P., Velli M., Bettarini L., 2008, Physics of Plasmas, 15, 012302
  * Lesieur & Ossia (2000) Lesieur M., Ossia S., 2000, Journal of Turbulence, 1, 007
  * Lewiner et al. (2003) Lewiner T., Lopes H., Vieira A. W., Tavares G., 2003, Journal of Graphics Tools, 8, 1–15
  * Linton et al. (2001) Linton M. G., Dahlburg R. B., Antiochos S. K., 2001, The Astrophysical Journal, 553, 905
  * Lorensen & Cline (1987) Lorensen W. E., Cline H. E., 1987, ACM SIGGRAPH Computer Graphics, 21, 163–169
  * Matthaeus & Montgomery (1980) Matthaeus W. H., Montgomery D., 1980, Annals of the New York Academy of Sciences, 357, 203
  * McKinney & Uzdensky (2012) McKinney J. C., Uzdensky D. A., 2012, Monthly Notices of the Royal Astronomical Society, 419, 573
  * McPherron (1979) McPherron R. L., 1979, Reviews of Geophysics, 17, 657
  * Mecke et al. (1994) Mecke K. R., Buchert T., Wagner H., 1994, Astronomy and Astrophysics, 288, 697
  * Mininni & Pouquet (2013) Mininni P. D., Pouquet A., 2013, Physical Review E, 87, 033002
  * Moffatt et al. (1994) Moffatt H. K., Kida S., Ohkitani K., 1994, Journal of Fluid Mechanics, 259, 241–264
  * Narayan et al. (2003) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, Publications of the Astronomical Society of Japan, 55, L69
  * Neronov & Vovk (2010) Neronov A., Vovk I., 2010, Science, 328, 73
  * Oishi et al. (2015) Oishi J. S., Low M.-M. M., Collins D. C., Tamura M., 2015, The Astrophysical Journal Letters, 806, L12
  * Parker (1983) Parker E. N., 1983, The Astrophysical Journal, 264, 635
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe
  * Pevtsov et al. (2003) Pevtsov A. A., Fisher G. H., Acton L. W., Longcope D. W., Johns-Krull C. M., Kankelborg C. C., Metcalf T. R., 2003, The Astrophysical Journal, 598, 1387
  * Pontin (2011) Pontin D. I., 2011, Advances in Space Research, 47, 1508
  * Pontin et al. (2011) Pontin D. I., Wilmot-Smith A. L., Hornig G., Galsgaard K., 2011, Astronomy & Astrophysics, 525, A57
  * Porth et al. (2014) Porth O., Komissarov S. S., Keppens R., 2014, Monthly Notices of the Royal Astronomical Society, 438, 278
  * Pouquet et al. (1976) Pouquet A., Frisch U., Léorat J., 1976, Journal of Fluid Mechanics, 77, 321–354
  * Priest (2016) Priest E., 2016, in Gonzalez W., Parker E., eds, Astrophysics and Space Science Library, Magnetic Reconnection: Concepts and Applications. Springer International Publishing, Cham, pp 101–142, doi:10.1007/978-3-319-26432-5_3, https://doi.org/10.1007/978-3-319-26432-5_3
  * Rappazzo et al. (2008) Rappazzo A. F., Velli M., Einaudi G., Dahlburg R. B., 2008, The Astrophysical Journal, 677, 1348
  * Reppin & Banerjee (2017) Reppin J., Banerjee R., 2017, Physical Review E, 96, 053105
  * Ruan et al. (2020) Ruan W., Xia C., Keppens R., 2020, The Astrophysical Journal, 896, 97
  * Sahni et al. (1998) Sahni V., Sathyaprakash B. S., Shandarin S. F., 1998, The Astrophysical Journal, 495, L5
  * Seta et al. (2020) Seta A., Bushby P. J., Shukurov A., Wood T. S., 2020, Physical Review Fluids, 5, 043702
  * Sheth et al. (2003) Sheth J. V., Sahni V., Shandarin S. F., Sathyaprakash B. S., 2003, Monthly Notices of the Royal Astronomical Society, 343, 22–46
  * Shukurov & Subramanian (2021) Shukurov A., Subramanian K., 2021, Astrophysical Magnetic Fields: From Galaxies to the Early Universe. Cambridge Astrophysics, Cambridge University Press, Cambridge, doi:10.1017/9781139046657, https://www.cambridge.org/core/books/astrophysical-magnetic-fields/231CE7B2753A87CEB4BFA55AE04412D1
  * Subramanian (2016) Subramanian K., 2016, Reports on Progress in Physics, 79, 076901
  * Subramanian et al. (2006) Subramanian K., Shukurov A., Haugen N. E. L., 2006, Monthly Notices of the Royal Astronomical Society, 366, 1437
  * Taylor (1986) Taylor J. B., 1986, Reviews of Modern Physics, 58, 741
  * Ting et al. (1986) Ting A. C., Matthaeus W. H., Montgomery D., 1986, The Physics of Fluids, 29, 3261
  * Tobias et al. (2011) Tobias S. M., Cattaneo F., Boldyrev S., 2011, MHD Dynamos and Turbulence, doi:10.48550/arXiv.1103.3138, http://arxiv.org/abs/1103.3138
  * Wiegand et al. (2014) Wiegand A., Buchert T., Ostermann M., 2014, Monthly Notices of the Royal Astronomical Society, 443, 241
  * Wilkin et al. (2007) Wilkin S. L., Barenghi C. F., Shukurov A., 2007, Physical Review Letters, 99, 134501
  * Yakhot & Pelz (1987) Yakhot V., Pelz R., 1987, The Physics of Fluids, 30, 1272
  * Zhdankin et al. (2014) Zhdankin V., Boldyrev S., Perez J. C., Tobias S. M., 2014, The Astrophysical Journal, 795, 127
  * Zhou et al. (2019) Zhou M., Bhat P., Loureiro N. F., Uzdensky D. A., 2019, Physical Review Research, 1, 012004
  * Zhou et al. (2020) Zhou M., Loureiro N. F., Uzdensky D. A., 2020, Journal of Plasma Physics, 86, 535860401
  * Zhou et al. (2022) Zhou H., Sharma R., Brandenburg A., 2022, Journal of Plasma Physics, 88, 905880602
  * Zrake (2014) Zrake J., 2014, The Astrophysical Journal, 794, L26
  * Zrake (2016) Zrake J., 2016, The Astrophysical Journal, 823, 39
  * Zrake & Arons (2017) Zrake J., Arons J., 2017, The Astrophysical Journal, 847, 57

##  Appendix A Algorithm for calculating Minkowski Functionals

In section 3.2.4, we have argued that local anisotropy in a 3D turbulent MHD
system can render it quasi-two-dimensional. To characterize this local
anisotropy, we analyze the shape of magnetic field isosurfaces. The required
algorithm must perform the following tasks,

  1. 1.

Extract isosurfaces from 3D volumetric data

  2. 2.

Calculate geometric quantities of these isosurfaces

###  A.1 Iso-surface extraction using marching cubes algorithm

We use the marching cubes algorithm to obtain isosurfaces from the data cube.
The input data is a scalar function $f$ on a 3D cubic grid. The algorithm
samples through each unit cube, hence the name. The algorithm starts by
comparing the value of $f$ at each grid point to the user-provided isosurface
value. Then, with the help of a pre-defined look-up table, triangles are
placed inside a unit cube, thus constituting a triangulation of the surface.
The initial algorithm for the look-up table by Lorensen & Cline (1987) had
face and internal ambiguities and was later modified by Lewiner et al. (2003)
to create a triangulated surface that is topologically consistent. For our
work, we use a Python implementation of the Lewiner et al. (2003) algorithm as
in skimage.measure.marching_cubes module.

###  A.2 Quantifying geometric surfaces using Minkowski Functionals

The marching cubes algorithm provides us with the position vector for vertices
of the triangles and the gradient direction at each vertex of individual
triangles in the order of the flow of the algorithm. To calculate Minkowski
Functionals of the triangulated isosurfaces, we closely follow Sheth et al.
(2003). To separate different connected isosurfaces, we identify the triangles
that share an edge with a given triangle.

On a triangulated isosurface, two triangles can share a common edge. Hence, on
a connected isosurface, each triangle has three other triangles with which it
shares its edges. All the triangles can be viewed together as an undirected
graph. Each node of the graph is a triangle. The three edges that come out of
a node connect to the the three nodes corresponding to the three triangles
that share a common edge with the former. Firstly, we construct a sparse
matrix indicating triangles that share common edges. Then, we use a graph
algorithm implemented in scipy.cs_graph.connected_components to identify
individual connected structures. To avoid inconsistencies, we impose a lower
threshold on the number of triangles making up the isosurfaces. If the number
of triangles is less than 10, we reject those isosurfaces.

We restrict to a particular connected structure and obtain the geometric
quantities leading to MFs as follows:

##### Surface area

The sum of the areas of individual triangles gives the surface area. Using the
position vectors for vertices, we calculate the vectors along the edges of the
triangles. The area of individual triangles is calculated by taking the
absolute value of the cross product of the vectors along any two edges.

##### Volume

To calculate the volume, consider one particular triangle on the surface. It
subtends a cone at the grid’s origin $(0,0,0)$. The volume of the triangular
cone is,

$V_{i}=\frac{1}{3}A_{i}(\bm{P}_{c}\cdot\hat{\mbox{\boldmath$n$}}{}),$ (25)

where $A$ is the area of the base triangle, $\hat{\mbox{\boldmath$n$}}{}$ is
its normal, and $\bm{P}_{c}$ is the position vector of the centroid. The
height of the cone is given by
$\mathbf{P}_{c}\cdot\hat{\mbox{\boldmath$n$}}{}$. Now, consider two triangles
$1$ and $2$ on the ‘opposite’ sides of the surface as in Fig. 18. The triangle
$1$’s height is ‘negative’ since the centroid and normal form an obtuse angle.
Hence, the volume is ‘negative’. However, in the case of triangle $2$, the
volume is ‘positive’. It cancels out the negative contribution from triangle
$1$, thereby giving us only the volume inside the closed surface. Note that
diagonally opposite triangles do not necessarily cancel each other’s
contribution exactly. We expect that the calculated value converges to the
actual value in the limit of a large number of triangles for a given surface.
While we described the calculation of the volume subtended at the grid’s
origin, in practice, we use the centroid of the given structure. However, the
above analysis goes through as is.

Figure 18: Calculation of volume by summing contributions due to individual
tetrahedra formed by the triangles and the origin

##### Integrated Mean Curvature (IMC)

The triangulated surface is curved only at the edges of the triangles and not
on the faces. To calculate the curvature at the edges, we fit a cylinder of
radius $r$ between two adjoining triangles as shown in Fig. 19. The tangents
at the point of contact of the two triangles ($i$ and $j$) with the cylinder
form an angle $\phi_{ij}$. The area of this curved surface is
$\phi_{ij}l_{ij}$, where $l_{ij}$ is the length of the common edge. The
principal curvatures of this cylindrical surface is $\kappa_{1}=\frac{1}{r}$
and $\kappa_{2}=0$. The two triangles can be oriented caved in(out),
contributing negatively(positively) to the IMC. We quantify the orientation
$\epsilon$ by calculating the angle between the normal to one of the triangles
and the line joining the centroids of the triangles, $\theta_{ij}$ (please see
Table 2). Hence, the integrated mean curvature of a structure is,

$C=\frac{1}{6\pi}\iint(\kappa_{1}+\kappa_{2})dS=\frac{\epsilon}{6\pi r}\iint
dS\\\ =\sum\frac{\epsilon}{6\pi r}rl_{ij}\phi_{ij}=\sum\frac{\epsilon
l_{ij}\phi_{ij}}{6\pi}.$ (26)

Orientation | $\theta_{ij}$ | $\epsilon$
---|---|---
Convex | $>\frac{\pi}{2}$ | +1
Concave | $<\frac{\pi}{2}$ | -1
Flat | $\frac{\pi}{2}$ | 0
Table 2:  Figure 19: Contribution to IMC from an edge. A cylinder of radius
$r$ is fit between the triangular surfaces. They subtend an angle $\phi_{ij}$,
and the common edge has length $l_{ij}$. The normal at triangle $1$ and the
line joining the triangles $1,2$ subtend an angle $\theta_{ij}$.

##### Gauss curvature

The Gauss curvature is calculated using the Gauss-Bonnet theorem, which
relates it to the Euler characteristic. The Gauss-Bonnet theorem states that

$\int_{M}K\,dS=2\pi\chi(M),$ (27)

where $K$ = $\kappa_{1}\kappa_{2}$. In particular, the above theorem holds
also for polyhedral surfaces. Hence, given a triangulation of the surface and
thereby the number of vertices $V$, edges $E$, and faces $F$, the Euler
characteristic is,

$\chi=V-E+F=V-\frac{3}{2}F+F=V-\frac{F}{2}.$ (28)

We have used the fact that each triangle has three edges, two triangles share
each edge, and thus, $E=\frac{3}{2}F$.

##  Appendix B Identification of sites of reconnection and nearby magnetic
field iso-surfaces

We compute the MFs for $|\bm{B}|$ and $|\bm{J}|$ fields. While computing MFs
for structures throughout the $1024^{3}$ domain is computationally demanding,
it is not necessary as we mainly want to focus on magnetic fields near
reconnection sites. We identify reconnection sites or grid points which have a
large current density value. We quantify this by taking a cutoff value for the
current density, $0.7*|\bm{J}|_{\text{max}}$. This threshold was chosen to
obtain a sufficient number of structures that ensure reliable statistics.

To avoid repetition of processing grid points pertaining to the same current
sheet, we isolate a set by the following method. We first apply the threshold
(mentioned above) to $|\bm{J}|$ to obtain a set of grid-points. Individual
clusters of points belonging to a current sheet are identified by:
(a) Ensuring furthest points in a given cluster are separated at the maximum
by $0.5\ L_{int}$ to prevent the inclusion of neighbouring current sheets,
(b) Ensuring that each cluster is separated from another by atleast $1.5\
L_{int}$.
In this manner, we also avoid double counting of magnetic field isosurfaces.
We choose the point with the maximum value among each cluster as the
representative point for that cluster. We define the current sheet in this
cluster to be the isosurface with isovalue $0.5*|\bm{J}|_{\text{max,
cluster}}$. Henceforth, we focus our attention on individual current sheets.
For a given current sheet, we restrict our domain size to be $L_{int}$ about
its representative point. We do so in order to identify magnetic field
isosurfaces that are in proximity to the current sheet.

To do the MF analysis, we adopt the following scheme to determine, first, the
isovalues for $|\bm{B}|$. We want to map out the radial profile of $|\bm{B}|$
around the current sheet. Ideally, an infinite number of rays with all
possible orientations are needed, but we use 50 rays with random orientations
to get an optimal balance between computational cost and sufficiency. We
parametrize the rays from $[-0.5,0.5]L_{int}$ with $100$ points. We calculate
the $|\bm{B}|$ value along the rays by linear interpolation. The maximum and
the minimum values of the $|\bm{B}|$ along each ray nearest to the
representative point are obtained. Among these extrema, the largest maximum
and the smallest minimum value of $|\bm{B}|$ are calculated and labelled as
$|\bm{B}|_{\text{max}}$ and $|\bm{B}|_{\text{min}}$ respectively. Finally, we
choose the isovalues for the magnetic field in the following manner,

$|\bm{B}|_{\text{max}}-\frac{i}{5}\times\left(\frac{|\bm{B}|_{\text{max}}-|\bm{B}|_{\text{min}}}{4}\right)\\\
\forall i\in{1,2,3,4,5}.$ (29)

This approach of selecting isovalues is to prevent the isosurfaces from being
either minuscule compared to the grid separation or volume-filling. Once the
magnetic isosurfaces are identified in this way, MFs can be calculated. Since
there can be multiple magnetic field isosurfaces, we want to consider the
isosurface closest to the current sheet. So, we calculate the distance between
the current sheet and the $|\bm{B}|$ field isosurfaces222The distance between
two sets, $d(C,D)=min\ d(x_{C},x_{D})$ and choose the isosurface with the
smallest distance to the current sheet.

A particular simulation can be considered one realization from an ensemble of
turbulent configurations. The reconnection sites that we obtain form a sample
of all possible realizations. Having obtained the MFs of $|\bm{B}|$ and
$|\bm{J}|$ isosurfaces, we calculate the sample average and standard
error333It is the standard deviation divided by the square root of the sample
size. of the isosurfaces across all isovalues considered. We then plot the
evolution of sample averages with standard errors as error bars.

##  Appendix C MF analysis of current sheets

###  C.1 Evolution of $P$ and $F$ for current sheets

Here, we display the evolution of the Minkowski Functionals for the current
sheets in Fig. 20 and in Fig. 21 for the helical and nonhelical cases
respectively. We already noted in Fig. 17 that the current isosurfaces in both
the helical and nonhelical cases are sheet-like. We plot the evolution of
$l_{1}$, $l_{2}$ and $l_{3}$ for the current sheets and observe that they
increase with time. This behaviour of the current sheet is consistent with the
fact that the system relaxes to larger and larger scales. We then plot the
evolution of $P$ and $F$. The $P$ and $F$ of the current sheets are almost
constant with time for both helical $(P\sim t^{0.15\pm 0.03},F\sim t^{0.02\pm
0.04})$ and nonhelical $(P\sim t^{0.01\pm 0.03},F\sim t^{-0.11\pm 0.04})$
cases. We also report the linear fits for the length scales $l$ in the plot.

Figure 20: Helical: The time evolution of length scales of current sheet
isosurfaces is shown in the left panel. The different dashed curves are the
least square fit. Time evolution of $P$ and $F$ in the early stages of the
decay (when the $B_{rms}$ curve behaviour is close to the theoretical
prediction of $t^{-2/7}$) is shown in the right panel.

Figure 21: Non-helical: The time evolution of length scales of current sheet
isosurfaces is shown in the left panel. The different dashed curves are the
least square fit. Time evolution of $P$ and $F$ in the early stages of the
decay (when the $B_{rms}$ curve behaviour is close to the theoretical
prediction of $t^{-1/2}$) is shown in the right panel.

##  Appendix D Results from forced hydrodynamic turbulence

We simulated forced hydrodynamic turbulence with a grid size of $512^{3}$,
average forcing wavenumber $k_{f}\approx 1.5$, and zero initial velocity
field. When the turbulence is fully developed (we confirm the $k^{-5/3}$
spectrum), we perform the morphological analysis using MF. At time $t=80$, we
look at the saturation of the velocity field and confirm that turbulence has
fully developed in the system. In 3D hydrodynamic turbulence, the high
vorticity isosurfaces are present as tubular structures (Moffatt et al.,
1994). Hence, we expect structures with isovalues higher than
$\bm{\omega}_{rms}$ at a particular time to be statistically filamentary
(large $F$). A preliminary visual confirmation is seen in Fig. 22. We
calculate MF of $\bm{\omega}$ at an isovalue of $5\bm{\omega}_{\text{rms}}$.
Initially, we plot the histogram of planarity and filamentarity in Fig. 23. We
plot the mean volume in each bin to understand where the structures of various
sizes (volume) contribute. On sequentially removing surfaces with volumes less
than a certain fraction of the total mean volume, the peak of the histogram
shifts to higher values for filamentarity and lower values for planarity.
Hence, larger structures are predominantly filamentarity with an almost
circular cross-section (with normal along the long axis). Hence, the system is
statistically tube-like.

Figure 22: Isosurfaces of vorticity fields at isovalue
$5\bm{\omega}_{\text{rms}}$ of a forced hydrodynamics $512^{3}$ simulations at
$t=80$. From the figure, it is clear that the structures are akin to elongated
tubes.

Figure 23: Histogram for planarity (left panel) and filamentarity (right
panel) of all the vorticity structures at isovalue $5\bm{\omega}_{\text{rms}}$
in the forced hydrodynamics $1024^{3}$ simulation at $t=80$. Also shown in the
black circles is the average volume in each bin. NVC stands for no-volume-
cutoff. The darker shades represent statistics for structures considered above
a certain volume threshold denoted by $\overline{V}/n$ and $n$ = 1, 5 or 10.
The tendency of large structures to have smaller planarity and larger
filamentarity is clear from the graph. The same conclusion can also be reached
by progressively considering structures with larger and larger volumes, as we
have shown in this plot.