11institutetext: School of Astronomy and Space Science, Nanjing University,
Nanjing 210023, People’s Republic of China
11email: wyulei<EMAIL_ADDRESS>22institutetext: Key Laboratory
for Modern Astronomy and Astrophysics (Nanjing University), Ministry of
Education, Nanjing 210023, People’s Republic of China 33institutetext: Centre
for Mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven, Belgium

# A Method for Determining the Locations and Configurations of Magnetic
Reconnection within 3D Turbulent Plasmas

Yulei Wang 1122 Xin Cheng 1122 Yang Guo 1122 Jinhan Guo 112233 Mingde Ding
1122

###### Abstract

Context. Three-dimensional (3D) reconnection is an important mechanism for
efficiently releasing energy during astrophysical eruptive events, which is
difficult to be quantitatively analyzed especially within turbulent plasmas.

Aims. In this paper, an efficient method for identifying locations and
configurations of 3D reconnection from MHD data is developed.

Methods. This method analyzes the local nonideal electric field and magnetic
structure at an arbitrary position. As only performing algebraical
manipulations on the discrete field data and avoiding computationally
expensive operations like field-line tracing and root-finding, this method
naturally possesses high efficiency. To validate this method, we apply it to
the 3D data from a high-resolution simulation of a Harris-sheet reconnection
and a data-driven simulation of a coronal flux rope eruption.

Results. It is shown that this method can precisely identify the local
structures of discrete magnetic field. Through the information of nonideal
electric field and the geometric attributes of magnetic field, the local
structures of reconnection sites can be effectively and comprehensively
determined. For fine turbulent processes, both qualitative pictures and
quantitative statistical properties of small-scale reconnection structures can
be obtained. For large-scale solar simulations, macro-scale magnetic
structures such as flux ropes and eruption current sheets can also be
recognized.

Conclusions. We develop a powerful method to analyze multi-scale structures of
3D reconnection. It can be applied not only in MHD simulations but also in
kinetic simulations, plasma experiments, and _in-situ_ observations.

###### Key Words.:

3D magnetic reconnection – turbulence – data analysis

##  1 Introduction

Magnetic reconnection is believed to be the driving mechanism of various
astrophysical eruptive phenomena, during which magnetic energy is rapidly
released to heat plasma and accelerate charged particles. The classic theory
of magnetic reconnection is built on two-dimensional (2D) steady models, but
in realistic physical environments, the three-dimensional (3D) effects cannot
be neglected (Priest & Forbes, 2000). Far more complex than 2D scenarios, 3D
reconnection is hard to be investigated analytically, especially for unsteady
reconnection processes coupled with significantly multi-scale and nonlinear
characteristics. Therefore, numerical simulations play vital roles in studies
of 3D reconnection (Ji et al., 2022). With the recent development of large-
scale high-performance computers, high-resolution simulations of both kinetic
and MHD scales have been performed, which proves the formation of self-
sustained turbulence induced by 3D reconnection (see Huang & Bhattacharjee,
2016; Kowal et al., 2020; Zhang et al., 2021; Comisso & Sironi, 2022; Dong et
al., 2022; Wang et al., 2023). Within the turbulent regions, the large-scale
structures keep cascading toward smaller ones following a power-law spectrum,
resulting in complicated patterns and chaotic magnetic structures. Though
high-resolution simulations provide an unprecedented opportunity for
understanding the fine processes of 3D reconnection, how to efficiently locate
the reconnection sites and analyze their properties from the massive discrete
data has become a new challenge (Vlahos & Isliker, 2023).

A typical location of 3D reconnection is at a magnetic null point, where the
magnetic strength vanishes. In the adjacent region of a null point, magnetic
field forms a fan-spine structure, the local structures of which can be
evaluated by the magnetic Jacobian matrix $\mathbf{D}_{B}=\partial_{j}B^{i}$,
where $\partial_{j}\equiv\partial/\partial x_{j}$ and $i,j=1,2,3$ (see Lau &
Finn, 1990; Parnell et al., 1996). Haynes & Parnell (2007) developed a
trilinear method for locating 3D null points in discrete numerical data, which
has been widely applied and can also be directly used in turbulent cases.

However, 3D reconnection can happen without the presence of null points, thus
being more difficult to be analyzed. Schindler et al. (1988) and Hesse &
Schindler (1988) proposed the general magnetic reconnection theory and
obtained the condition for the “global” effect of reconnection with finite
magnetic field, namely,

$\int_{P_{1}}^{P_{2}}E_{\parallel}\mathrm{d}s\neq 0\,,$ (1)

where $E_{\parallel}$ is the electric field strength parallel with magnetic
field $\bf B$, the integration is taken along a magnetic field line, and the
“global” effect refers to the two plasma elements in the ideal region, $P_{1}$
and $P_{2}$, can “feel” the reconnection and will move to different field
lines later.

3D reconnection can take place at locations where strong local current
concentrates and various types of 3D reconnection have been discovered such as
the torsional/shearing reconnection within the fan/spine near a null point,
the separator reconnection, the hyperbolic-flux-tube (HFT) reconnection, and
the braid reconnection (see Pontin & Priest, 2022). Meanwhile, the local
magnetic structures at the reconnection region determine the properties of 3D
reconnection. Priest & Forbes (1989) proposed the singular line (SL)
reconnection which forms X-type magnetic structures on the plane perpendicular
to the local magnetic field. Wilmot-Smith & Priest (2007) provided an example
of 3D reconnection at the center of a flux tube with footpoints rotating
toward different directions. Parnell et al. (2010a) found the O-type
reconnection magnetic structures even along a separator connecting two null
points in their simulations.

To analyze locations of 3D reconnection from the discrete magnetic field data
obtained in simulations, observations, and experiments, both global and local
methods have been developed. The global methods, relying on the field-line
tracing technique, are closely related to the 3D reconnection theory but
mainly suitable for analyzing large-scale laminar magnetic fields. For
example, the quasi-separatrix layer (QSL) method, first proposed by Priest &
Démoulin (1995) and later improved by Titov et al. (2002) and Titov (2007),
has been widely applied to study coronal magnetic field (Démoulin et al.,
1996; Demoulin et al., 1996, 1997; Aulanier et al., 2006; Titov et al., 2009;
Li et al., 2021a; Zhong et al., 2021; Li et al., 2022; Guo et al., 2023). QSLs
are locations where the connectivity of magnetic field, evaluated by the
squashing factor $Q$, changes significantly. Though $Q$-factor is not a direct
indicator of reconnection (Reid et al., 2020), it is useful for investigating
magnetic topologies and exhibiting locations where 3D reconnection might
happen. The calculation of $Q$-factor for discrete data relies on field-line
tracing, which can consume significant computation resources without
appropriate optimizations. Recently, the efficiency of the QSL method has been
improved via the hardware accelerations of GPUs (Zhang et al., 2022). As
another algorithm relying on field-line tracing, the method developed by
Haynes & Parnell (2010) can detect the global laminar magnetic topology
(magnetic skeleton) in numerical results, including null points, spines,
separatrix surfaces, and separators. By using this method, Parnell et al.
(2010b) statistically proved the importance of separator reconnection during
the interaction between emerging and overlying flux. Komar et al. (2013)
proposed a similar method for tracing the separators within discrete data
based on the pre-located null points.

Different from global methods, the local analysis focuses on the differential
structures depending only on the adjacent region of a position and thus
naturally possesses high efficiency. More importantly, it also provides
statistical laws that are vital for understanding the physics of turbulent
reconnection. As a straightforward method, the distributions of
$E_{\parallel}$ and $J_{\parallel}$ can reflect the sites of 3D reconnection
according to the definition Eq. (1), which has been used in the visualization
of turbulent reconnection simulations (see Huang & Bhattacharjee, 2016;
Isliker et al., 2019; Dong et al., 2022).

For 3D systems imposed with a strong uniform background magnetic field, the
methods of locating reconnection structures can be simplified as 2D methods.
For instance, Zhdankin et al. (2013) proposed an algorithm for identifying the
geometrical and physical properties of current sheets. Li et al. (2021b)
developed the magnetic flux transport method which was recently applied to a
kinetic plasma turbulence simulation by Li et al. (2023).

For general 3D cases without a strong guide field, local reference frames are
necessary to analyze small-scale structures. Parnell et al. (2010a) set the
direction of a pre-determined separator as the normal direction of the
projection plane, on which the 2D projected magnetic field is classified.
Kowal et al. (2020) defined a local frame based on the shear tensor of
magnetic field, which is further used to determine the properties of current
sheets. Recently, Lapenta (2021) developed a reconnection-site identification
method via the electric drift speed, which performs well for data of kinetic
simulations.

Methods independent with local frames have also been developed for the
identification of turbulent reconnection structures (see Vlahos & Isliker,
2023), for instance, the Phase Coherence Index method (Hada et al., 2003), the
Partial Variance of Increments (PVI) method (Greco et al., 2017), and the
fractal dimension method by box-counting (Isliker et al., 2019). Although
these methods can produce useful statistical results, they cannot analyze the
magnetic structures at reconnection sites.

In this paper, we introduce an efficient method for locating small-scale
reconnection structures in 3D discrete data involving turbulent eddies. Both
the theoretical basis and the numerical algorithm of this method are presented
in detail. Through two typical benchmarks, it is well validated and presents a
promise in qualitative and quantitative analysis of 3D reconnection
simulations containing multi-scale structures.

The following parts of this paper are organized as follows. In Sect. 2, we
exhibit the theoretical basis of our method. The numerical strategy is
presented in Sect. 3. We provide a Matlab implementation of this method and
introduce its usage in Sect. 4. This method is benchmarked by two typical
simulations which are detailed in Sects. 5 and 6. The results of this paper
are summarized and discussed in Sect. 7.

##  2 Theory

###  2.1 The General Magnetic Reconnection Theory

We first briefly summarize the main concepts of the general magnetic
reconnection theory developed by Hesse & Schindler (1988).

Nonvanishing magnetic field can be (at least locally) expressed by the Euler
potentials as $\bf{B}=\nabla\alpha\times\nabla\beta$. The three vectors
$\nabla\alpha$, $\nabla\beta$, and
$\hat{\mathbf{b}}=\mathbf{B}/\left|\mathbf{B}\right|$ compose a local frame
system $\left(\alpha,\beta,s\right)$, where $s$ denotes the arc coordinate
along a field line. An $\left(\alpha,\beta\right)$-pair corresponds to a field
line. $\nabla\alpha$ and $\nabla\beta$ span the local plane normal to $\bf B$
but, generally speaking, they are not orthogonal or normalized.

If setting the vector potential as $\bf A=\alpha\nabla\beta$ satisfying the
gauge $\bf A\cdot B=0$, Ohm’s law $\bf E=-v\times B+N$ and Faraday’s law
produce the evolution rule of the Euler potentials as

$\displaystyle\dot{\alpha}$
$\displaystyle=-\frac{\partial\psi}{\partial\beta}-N^{\beta}\,,$ (2a)
$\displaystyle\dot{\beta}$
$\displaystyle=\frac{\partial\psi}{\partial\alpha}-N^{\alpha}\,,$ (2b)

where $\bf N$ is the general nonideal term,
$\psi=\phi+\alpha\partial\beta/\partial t$, $\phi$ is the electric scalar
potential, and the dot operator denotes the total derivative of time
$\mathrm{d}/\mathrm{d}t$. Notice that $\partial\psi/\partial
s=-N^{s}=-E_{\parallel}$.

Hesse & Schindler (1988) generalized the definition of reconnection with
finite magnetic field by the violation of line conservation. Supposing two
plasma elements $P_{1}$ and $P_{2}$ outside the nonideal region are connected
by a field line, they will later move to two different field lines if and only
if $\psi$ in Eq. 2 initially takes different values at $P_{1}$ and $P_{2}$,
which directly results in the condition Eq. 1. Because $\bf N=0$ at $P_{1}$
and $P_{2}$, the global effect is independent with the perpendicular
components of $\bf N$, namely, $N^{\alpha}$ and $N^{\beta}$.

Different from 2D reconnection taking place at an X-point, 3D reconnection
happens within a finite nonideal volume (Priest et al., 2003), in which there
exist “special” field lines. After defining a quasi-potential
$\Xi\left(\alpha,\beta\right)\equiv-\int_{P_{1}}^{P_{2}}E_{\parallel}\left(\alpha,\beta,s\right)\mathrm{d}s$,
the global differences between the evolutions of Euler potentials at $P_{1}$
and $P_{2}$ can be expressed by (Hesse et al., 2005)

$\displaystyle\left.\dot{\alpha}\right|_{P_{2}}-\left.\dot{\alpha}\right|_{P_{1}}$
$\displaystyle=-\frac{\partial\Xi}{\partial\beta}\,,$ (3a)
$\displaystyle\left.\dot{\beta}\right|_{P_{2}}-\left.\dot{\beta}\right|_{P_{1}}$
$\displaystyle=\frac{\partial\Xi}{\partial\alpha}\,.$ (3b)

According to this equation, if
$\partial\Xi/\partial\alpha=\partial\Xi/\partial\beta=0$ ($\Xi$ has an
extremal value in the $\alpha$-$\beta$ space), $P_{1}$ and $P_{2}$ will remain
on the same field line and thus cannot feel the global reconnection effect.
Notice that $\partial\Xi/\partial\alpha=\partial\Xi/\partial\beta=0$ is
equivalent to $\nabla_{\perp}\Xi=0$, where $\nabla_{\perp}$ denotes the
gradient perpendicular to $\bf B$. Moreover, Hesse et al. (2005) generalized
the definition of reconnection rate to cases without separators or
separatrices and proved that the extremal value of $\Xi$ equals the
reconnection rate. The extremal-$\Xi$ lines can be treated as generalizations
of magnetic neutral lines or separators. Wyper & Hesse (2015) later
generalized this theory to complex nonideal regions with multiple peaks of
$\Xi$.

###  2.2 Local Analysis of 3D Reconnection

The theory of global reconnection is well-defined but is not convenient to
directly apply for analyzing discrete magnetic data, especially for cases with
strong turbulence. As revealed by various high-resolution simulations, the 3D
reconnection regions are composed of multi-scale current sheets that are far
more complex than theoretical models (see Dong et al., 2022; Wang et al.,
2023; Ye et al., 2023), which significantly enhances the difficulties in
selecting appropriate start/end positions for field-line mappings. Complete
samplings of field lines within the entire system might produce reliable
results, which, however, can hardly be operated as limited by computational
resources and numerical integration errors. Moreover, besides global
information, local reconnection effects and structures are also valuable for
studying the mechanism of 3D reconnection. Therefore, we should seek local
parameters of reconnection.

According to the condition Eq. 1, a _necessary_ condition for a field line to
undergo the global reconnection is that it passes nonideal regions with
$E_{\parallel}\neq 0$. In other words, we can use the distribution of
$E_{\parallel}$ to exhibit regions where global reconnection effects can
happen. Provide that only Joule dissipation is considered and the resistivity
$\eta$ is uniform, $E_{\parallel}$ can be replaced by $J_{\parallel}$ since
$E_{\parallel}=\mathbf{N\cdot\hat{b}}=\eta\mathbf{J\cdot\hat{b}}=\eta
J_{\parallel}$.

Within the $E_{\parallel}\neq 0$ regions, we can further locate field lines
with extremal values of $\Xi$. A _sufficient_ condition for a field line to
have an extremal $\Xi$ is that all the $E_{\parallel}\neq 0$ positions
threaded by this field line satisfy $\nabla_{\perp}E_{\parallel}=0$, which can
be directly proven by

$0=\int_{P_{1}}^{P_{2}}\nabla_{\perp}E_{\parallel}\mathrm{d}s=\nabla_{\perp}\int_{P_{1}}^{P_{2}}E_{\parallel}\mathrm{d}s=\nabla_{\perp}\Xi\,.$
(4)

In other words, in nonideal regions, locations with
$\nabla_{\perp}E_{\parallel}=0$ approximately reveal the distributions of
extremal-$\Xi$ lines, though the trajectories of extremal-$\Xi$ lines are not
necessary to always pass through such locations. $E_{\parallel}$ and
$\nabla_{\perp}E_{\parallel}$ are two useful and easy-to-calculate local
parameters that can provide implications about the locations of global
reconnection effects, though they are not equivalent to the original integral
definitions in Eqs. 1 and 3.

The local reconnection effects, felt by plasma elements inside nonideal
regions, correspond to the local violation of line conservation, which are
blocked out by the global reconnection theory but still worth investigating.
The equation governing the local effects can be obtained by performing
$\partial/\partial s$ on both sides of Eq. 2, namely,

$\displaystyle\frac{\partial\dot{\alpha}}{\partial s}$
$\displaystyle=\frac{\partial E_{\parallel}}{\partial\beta}-\frac{\partial
N^{\beta}}{\partial s}\,,$ (5a)
$\displaystyle\frac{\partial\dot{\beta}}{\partial s}$
$\displaystyle=-\frac{\partial E_{\parallel}}{\partial\alpha}-\frac{\partial
N^{\alpha}}{\partial s}\,,$ (5b)

which implies that the change of $\left(\dot{\alpha},\dot{\beta}\right)$ along
a magnetic field line is governed by two terms, namely, the perpendicular
gradient of $E_{\parallel}$ and the parallel gradient of
$\mathbf{N}_{\perp}=\left(N^{\alpha},N^{\beta}\right)$. It can be proved that
the line conservation ($\partial\dot{\alpha}/\partial
s=\partial\dot{\beta}/\partial s=0$) is equivalent to $\bf
B\times\left(\nabla\times N\right)=0$, while the condition of flux
conservation is $\bf\nabla\times N=0$ (Hesse & Schindler, 1988).

Furthermore, it is helpful to evaluate the importance of the $E_{\parallel}$
and $\bf N_{\perp}$ terms in Eq. 5. For instance, at locations with
$\nabla_{\perp}E_{\parallel}=0$, the violation of line conservation is fully
determined by the $\bf N_{\perp}$ term. Moreover, suppose that the effects of
$\bf N_{\perp}$ term are ignorable, the local effects would be dominated by
the $E_{\parallel}$ term, which means that analysis of $E_{\parallel}$ can
reveal the properties of not only global but also local effects of
reconnection. However, because a quantitative comparison of the two terms
relies on a local flux coordinate frame
$\left(\nabla\alpha,\nabla\beta\right)$, which is difficult to establish for
arbitrary magnetic fields. In Appendix A, we attempt to discuss this problem
and propose a statistical method to compare the effects of the two terms.

###  2.3 Magnetic Structures at Reconnection Sites

As the magnetic null points have been well studied (see Parnell et al., 1996),
here we analyze the local magnetic structures near reconnection locations with
finite magnetic field. To begin with, we define a local frame at the origin
$\mathbf{r}_{0}$ by three orthogonal unit vectors satisfying
$\hat{\mathbf{e}}_{3}=\hat{\mathbf{b}}\left(\mathbf{r}_{0}\right)$,
$\hat{\mathbf{e}}_{2}\perp\hat{\mathbf{e}}_{3}$, and
$\hat{\mathbf{e}}_{1}=\hat{\mathbf{e}}_{2}\times\hat{\mathbf{e}}_{3}$. In this
frame, the magnetic field at $\mathrm{d}\mathbf{r}=\mathbf{r}-\mathbf{r}_{0}$
can be approximated by

$\mathbf{B}\left(\mathrm{d}\mathbf{r}\right)=B_{0}\hat{\mathbf{e}}_{3}+\mathrm{d}\mathbf{r}\cdot\nabla\mathbf{B}\left(\mathbf{r}_{0}\right)\,,$
(6)

where $B_{0}=\left|\mathbf{B}\left(\mathbf{r}_{0}\right)\right|>0$ is the
magnetic strength at $\mathbf{r}_{0}$ and $\left|\mathrm{d}\mathbf{r}\right|$
should be small enough to guarantee the validation of linear-field assumption.
$\hat{\bf{e}}_{1}$ and $\hat{\bf{e}}_{2}$ span the same normal plane
$\mathrm{d}r_{3}=0$ as established by $\nabla\alpha$ and $\nabla\beta$.

The second term in Eq. 6 can be further written as two terms,

$\mathrm{d}\mathbf{r}\cdot\nabla\mathbf{B}\left(\mathbf{r}_{0}\right)=\left[\left(\begin{array}[]{c}\mathrm{d}r_{1}\\\
\mathrm{d}r_{2}\\\ 0\end{array}\right)+\left(\begin{array}[]{c}0\\\ 0\\\
\mathrm{d}r_{3}\end{array}\right)\right]\cdot\nabla\mathbf{B}=\left(\begin{array}[]{c}B_{\perp}^{1}\\\
B_{\perp}^{2}\\\ B_{\parallel}\end{array}\right)+\tilde{\bf B}\,,$ (7)

where $B_{\perp}^{1}$ and $B_{\perp}^{2}$ are two components of the
perpendicular magnetic field defined by

$\displaystyle\mathbf{B}_{\perp}$ $\displaystyle=\mathbf{M}\cdot\mathbf{R}\,,$
(8) $\displaystyle\mathbf{M}$ $\displaystyle=\left(\begin{array}[]{cc}\partial
B^{1}/\partial R_{1}&\partial B^{1}/\partial R_{2}\\\ \partial B^{2}/\partial
R_{1}&\partial B^{2}/\partial
R_{2}\end{array}\right)\equiv\left(\begin{array}[]{cc}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{array}\right)\,,$ (13) $\displaystyle\mathbf{R}$
$\displaystyle=\left(\begin{array}[]{c}\mathrm{d}r_{1}\\\
\mathrm{d}r_{2}\end{array}\right)\equiv\left(\begin{array}[]{c}R_{1}\\\
R_{2}\end{array}\right)\,,$ (18)

and $B_{\parallel}=R_{1}\left(\partial B^{3}/\partial
R_{1}\right)+R_{2}\left(\partial B^{3}/\partial R_{2}\right)$ is the parallel
component. $\mathbf{B}_{\perp}$ and $B_{\parallel}$ compose the magnetic field
lying in the $\mathrm{d}r_{3}=0$ plane, while $\tilde{\bf
B}\equiv\mathrm{d}r_{3}\left(\partial\mathbf{B}/\partial\mathrm{d}r_{3}\right)$
represents the increments at locations of $\mathrm{d}r_{3}\neq 0$. Hereafter,
we analyze the local in-plane magnetic structures on the $\mathrm{d}r_{3}=0$
plane following the ideas of the SL reconnection by Priest & Forbes (1989) and
the O-type separator reconnection by Parnell et al. (2010a). Therefore,
$\tilde{\bf B}$ and $B_{\parallel}$ is not investigated in detail. Suppose
that the magnitudes of $\left|B_{\parallel}\right|$ and $\left|\tilde{\bf
B}\right|$ are small enough compared with $\left|\mathbf{B}_{\perp}\right|$,
the magnetic structure within a finite volume near $\mathbf{r}_{0}$ can be
well approximated by $\mathbf{B}_{\perp}$.

Because $\mathbf{B}_{\perp}$ is not a real magnetic field satisfying the
divergence-free condition, the trace of $\mathbf{M}$, evaluated by
$tr\left(\mathbf{M}\right)=M_{11}+M_{22}=-\partial
B^{3}/\partial\mathrm{d}r_{3}$, can be nonzero. Meanwhile, one should notice
that $\mathbf{B}_{\perp}$ vanishes at the frame origin
($\mathbf{R}=\left(0,0\right)^{T}$), where the superscript “$T$” denotes the
transpose operation. In other words, the definition of the local reference
frame makes an arbitrary position $\mathbf{r}_{0}$ a 2D null point of
$\mathbf{B}_{\perp}$. Therefore, its 2D local structures can be fully
determined by $\bf M$, which simplifies the theoretical analysis. Moreover,
for discrete fields, the structures of $\mathbf{B}_{\perp}$ at an arbitrary
grid can be directly obtained via $\bf M$ and the computational costs of
manipulations like interpolation and root-finding for searching null points
between grids can be saved.

Inspired by the method by Parnell et al. (1996), we rewrite $\mathbf{M}$ as

$\mathbf{M}=\left(\begin{array}[]{cc}M_{11}&\left(q-J_{3}\right)/2\\\
\left(q+J_{3}\right)/2&M_{22}\end{array}\right)\,,$ (19)

where $q=M_{12}+M_{21}$ and $J_{3}=J_{\parallel}=M_{21}-M_{12}$. Defining a
threshold current
$J_{\mathrm{thres}}=\sqrt{\left(M_{11}-M_{22}\right)^{2}+q^{2}}$, the
discriminant of the eigenvalue quadratic
$\left|\lambda\mathbf{I}-\mathbf{M}\right|=0$, becomes
$D\left(\mathbf{M}\right)=4Det\left(\bf M\right)-tr\left(\bf
M\right)^{2}=J_{3}^{2}-J_{\mathrm{thres}}^{2}$, where $Det\left(\bf M\right)$
is the determinant of $\bf M$. The eigenvalues are thus
$\lambda=\left(tr\left(\bf M\right)\pm\sqrt{-D\left(\bf M\right)}\right)/2$.
If $\left|J_{3}\right|<J_{\mathrm{thres}}$ (i.e.
$D\left(\mathbf{M}\right)<0$), the eigenvalues of $\mathbf{M}$ are distinct
real numbers; if $\left|J_{3}\right|>J_{\mathrm{thres}}$ (i.e.
$D\left(\mathbf{M}\right)>0$), the eigenvalues of $\mathbf{M}$ are two
conjugate complex numbers; if $\left|J_{3}\right|=J_{\mathrm{thres}}$ (i.e.
$D\left(\mathbf{M}\right)=0$), two eigenvalues are repeated real roots.

The local structures of $\mathbf{B}_{\perp}$ can be classified into nine types
based on the values of $tr\left(\bf M\right)$, $D\left(\bf M\right)$, and
$\lambda$ (see Table 1). Type 1 has a local X-type structure, being the 3D
generalization of 2D X-points and the location of SL reconnection if
$E_{\parallel}\neq 0$ (see Fig. 12a). Types 2 and 3 are 3D O-types but possess
radial components from the trace of $\bf M$. For different signs of
$tr\left(\bf M\right)$, field lines of types 2 and 3 show different directions
toward $\mathbf{r}_{0}$ (see Fig. 12b, c). Types 4 and 5 are respectively
source and sink structures dominated by the trace of $\bf M$ (see Fig. 12d,
e). Type 6 has one zero eigenvalue, implying the existence of a neutral line
(see Fig. 12f). The other types with zero traces are of 2D configurations,
which are unstable and very rare in 3D evolutions (Priest & Forbes, 2000).
Types 7 and 8 are 2D X and O points, while type 9 are anti-parallel lines
(Parnell et al., 1996).

Table 1: Types of $\mathbf{B}_{\perp}$. The star marker “$*$” denotes complex conjugate. $\mathbb{R}$ and $\mathbb{C}$ are the sets of real and complex numbers, respectively. Type | $tr\left(\bf M\right)$ | $D\left(\bf M\right)$ | $\lambda$ | Description
---|---|---|---|---
1 | $\neq 0$ | $<0$ | $\lambda_{1}\cdot\lambda_{2}<0$ | 3D X
2 | $>0$ | $>0$ | $\lambda_{1}=\lambda_{2}^{*}$ | 3D O (Repelling)
3 | $<0$ | $>0$ | $\lambda_{1}=\lambda_{2}^{*}$ | 3D O (Attracting)
4 | $\neq 0$ | $\leq 0$ | $\lambda_{1},\lambda_{2}>0$ | 3D Repelling
5 | $\neq 0$ | $\leq 0$ | $\lambda_{1},\lambda_{2}<0$ | 3D Attracting
6 | $\neq 0$ | $<0$ | $\lambda_{1}\ or\ \lambda_{2}=0$ | 3D Anti-parallel
7 | $=0$ | $<0$ | $\lambda_{1}=-\lambda_{2}\in\mathbb{R}$ | 2D X
8 | $=0$ | $>0$ | $\lambda_{1}=-\lambda_{2}\in\mathbb{C}$ | 2D O
9 | $=0$ | $=0$ | $\lambda_{1}=\lambda_{2}=0$ | 2D Anti-parallel

To investigated the geometric properties of $\mathbf{B}_{\perp}$, we further
decompose $\mathbf{M}$ into two parts as

$\bf M=M^{\prime}+T\,,$ (20)

where,

$\displaystyle\mathbf{M}^{\prime}$
$\displaystyle=\left(\begin{array}[]{cc}p&\left(q-J_{3}\right)/2\\\
\left(q+J_{3}\right)/2&-p\end{array}\right)\,,$ (23) $\displaystyle\mathbf{T}$
$\displaystyle=\left(\begin{array}[]{cc}tr\left(\bf M\right)/2&0\\\
0&tr\left(\bf M\right)/2\end{array}\right)\,,$ (26) $\displaystyle p$
$\displaystyle=\left(M_{11}-M_{22}\right)/2\,.$ (27)

$\bf M^{\prime}$ is a traceless matrix corresponding to the source-free part
of $\mathbf{B}_{\perp}$ which can be treated as a real magnetic field, while
$\bf T$ is the source part originating from the 3D effects. Correspondingly,
$\mathbf{B}_{\perp}$ can be written as
$\mathbf{B}_{\perp}=\mathbf{B}_{\perp}^{\prime}+\nabla_{p}G$, where

$\displaystyle\mathbf{B}_{\perp}^{\prime}$
$\displaystyle=\mathcal{J}\nabla_{p}F\,,$ (28) $\displaystyle F$
$\displaystyle=\frac{1}{4}\left[\left(q-J_{3}\right)R_{2}^{2}-\left(q+J_{3}\right)R_{1}^{2}\right]+pR_{1}R_{2}\,,$
(29) $\displaystyle G$ $\displaystyle=\frac{1}{4}tr\left(\bf
M\right)\left(R_{1}^{2}+R_{2}^{2}\right)\,,$ (30) $\displaystyle\nabla_{p}$
$\displaystyle=\left(\partial/\partial R_{1},\partial/\partial
R_{2}\right)^{T}\,,$ (31)

and

$\mathcal{J}=\left(\begin{array}[]{cc}0&1\\\ -1&0\end{array}\right)\,\\\ $
(32)

is the symplectic matrix. Because the second term $\nabla_{p}G$ is isotropic,
the anisotropic properties of $\mathbf{B}_{\perp}$ are completely determined
by the first term $\mathbf{B}_{\perp}^{\prime}$. As a 2D magnetic field, the
integration curves of $\mathbf{B}_{\perp}^{\prime}$ are contour lines of its
Hamiltonian $F$.

Following the method by Parnell et al. (1996), we define a new coordinate
$\mathbf{X}=\left(X_{1},X_{2}\right)^{T}$ satisfying

$\mathbf{X}=\left(\begin{array}[]{cc}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\\\ \end{array}\right)\mathbf{R}\,,$ (33)

where the rotation angle $\theta$ satisfies $\tan 2\theta=-2p/q$. Then $F$ can
be transformed to

$F=\frac{1}{4}\left[X_{2}^{2}\left(J_{\mathrm{thres}}-J_{3}\right)-X_{1}^{2}\left(J_{\mathrm{thres}}+J_{3}\right)\right]\,.\\\
$ (34)

Therefore, the field lines of $\mathbf{B}^{\prime}$ are hyperbolic lines if
$\left|J_{3}\right|<J_{\mathrm{thres}}$ (Fig. 1a), are concentric elliptical
curves if $\left|J_{3}\right|>J_{\mathrm{thres}}$ (Fig. 1b), and are parallel
straight lines if $\left|J_{3}\right|=J_{\mathrm{thres}}$.

To quantitatively evaluate the geometric shapes of
$\mathbf{B}_{\perp}^{\prime}$ lines, we define an angle,

$\displaystyle\theta_{\mathrm{eig}}$
$\displaystyle=\begin{cases}\arctan\left(\sqrt{\left(J_{\mathrm{thres}}/J_{3}\right)^{2}-1}\right)\,,&\left|J_{3}\right|<J_{\mathrm{thres}}\,,\\\
\arctan\left(\sqrt{\left(J_{3}/J_{\mathrm{thres}}\right)^{2}-1}\right)\,,&\left|J_{3}\right|>J_{\mathrm{thres}}\,.\end{cases}$
(35)

For $\left|J_{3}\right|<J_{\mathrm{thres}}$, $\theta_{\mathrm{eig}}$ is the
acute angle spanned by the separatrix lines (see Fig. 1a). For
$\left|J_{3}\right|>J_{\mathrm{thres}}$, $\theta_{\mathrm{eig}}$ is the angle
spanned by two minor-axis vertices and one major-axis vertex of an elliptical
curve (see Fig. 1b). According to the definition,
$\theta_{\mathrm{eig}}\in\left[0,90^{\circ}\right]$. If $J_{3}\to
J_{\mathrm{thres}}$, then $\theta_{\mathrm{eig}}\to 0$, which corresponds to
the anti-parallel case. Meanwhile, $\theta_{\mathrm{eig}}\to 90^{\circ}$ means
that the X-type $\mathbf{B}^{\prime}_{\perp}$ tends to be a potential field
satisfying $J_{3}\to 0$ while the O-type $\mathbf{B}^{\prime}_{\perp}$
represents as concentric circle lines governed only by the anti-symmetric part
of $\mathbf{M}^{\prime}$ (the parallel current $J_{3}$);

Figure 1: Cartoons of $\mathbf{B}_{\perp}^{\prime}$ lines of X (a) and O (b)
types. The blue curves depict field lines. The definitions of
$\theta_{\mathrm{eig}}$ are labeled by orange markers.

Because $\theta_{\mathrm{eig}}$ only reflects the property of the source-free
part $\bf M^{\prime}$, we define another parameter to evaluate the importance
of $\bf T$, namely the trace ratio

$R_{\mathrm{tr}}=\|\mathbf{T}\|/\|\mathbf{M}^{\prime}\|\,,$ (36)

where $\|\cdot\|$ is the Euclidean norm. Larger $R_{\mathrm{tr}}$ means larger
effects from the trace part.

##  3 Algorithm

Combining the theory of general reconnection and the analysis of local
magnetic structures, we develop a numerical method to locally analyze 3D
reconnection sites. When introducing this algorithm, we suppose the discrete
data are uniformly sampled in a Cartesian frame, which, however, can be
generalized into arbitrary coordinate systems via transformations.

###  3.1 Step 1: Extract Grids with Large $\left|E_{\parallel}\right|$

As discussed in Sect. 2.2, we can select grids with finite values of
$E_{\parallel}$ to represent the general reconnection sites. But, for
numerical data, $E_{\parallel}$ can hardly be zero. Therefore, one can define
a threshold value of $E_{\mathrm{thres}}>0$ to extract a subset of
$\left|E_{\parallel}\right|>E_{\mathrm{thres}}$ from the entire domain to
reduce the computational costs of the following procedures. However, if the
data scale is acceptable or $E_{\mathrm{thres}}$ is difficult to choose, to
guarantee statistical completeness, one can also skip this step and analyze
all grids.

###  3.2 Step 2: Analyze Local Magnetic Structures

First, approximate the 3D magnetic Jacobian matrix $\mathbf{D}_{B}$ on each
grid via two-order central difference as

$\mathbf{D}_{B}^{lmn}=\left(\frac{\mathbf{B}^{l+1,m,n}-\mathbf{B}^{l-1,m,n}}{2\Delta
x^{1}},\frac{\mathbf{B}^{l,m+1,n}-\mathbf{B}^{l,m-1,n}}{2\Delta
x^{2}},\frac{\mathbf{B}^{l,m,n+1}-\mathbf{B}^{l,m,n-1}}{2\Delta
x^{3}}\right)\,,$ (37)

where $l$, $m$, and $n$ are indices of $x^{1}$, $x^{2}$, and $x^{3}$,
respectively.
$\mathbf{B}^{lmn}=\left(B^{1,lmn},B^{2,lmn},B^{3,lmn}\right)^{T}$ is the
magnetic field at $\mathbf{r}^{lmn}$. Generally speaking, the two-order
spatial difference is accurate enough for most cases, but one can replace it
with higher-order schemes if necessary.

Figure 2: Schematic diagram of the local frame at a grid $\mathbf{r}^{lmn}$.
The purple arrow lines span the local frame. The light blue shade presents the
MPP. The blue dots depict four adjacent grids on the MPP that are one step-
length ($\Delta L$) away from the origin point. The green dots are six
adjacent grids surrounding $\mathbf{r}^{lmn}$ in the frame of original
discrete data.

Second, construct local frame at $\mathbf{r}^{lmn}$ (see Fig. 2). The base
vector on the third direction (the direction of the local magnetic field) is
$\hat{\bf
e}_{3}^{lmn}=\left(B^{1,lmn},B^{2,lmn},B^{3,lmn}\right)/\left|\mathbf{B}^{lmn}\right|$.
$\hat{\bf e}_{2}^{lmn}$ can be any unit vector perpendicular to $\hat{\bf
e}_{3}^{lmn}$. As a convenient choice, we set it as $\hat{\bf
e}_{2}^{lmn}=\left(B^{2,lmn},-B^{1,lmn},0\right)/\sqrt{\left(B^{1,lmn}\right)^{2}+\left(B^{2,lmn}\right)^{2}}$.
For the special case of $B^{1,lmn}=B^{2,lmn}=0$, let $\hat{\bf
e}_{2}^{lmn}=\left(0,1,0\right)$. The last base vector can be directly
obtained by $\hat{\bf e}_{1}^{lmn}=\hat{\bf e}_{2}^{lmn}\times\hat{\bf
e}_{3}^{lmn}$. Hereafter, the local plane spanned by $\hat{\bf e}_{1}^{lmn}$
and $\hat{\bf e}_{2}^{lmn}$ within a cell is called the magnetic projection
plane (MPP).

Third, transform $\mathbf{D}_{B}^{lmn}$ into the local frame. Because both the
original and the local frames are Cartesian frames and the base vector systems
are orthogonal normalized ones, the transformation matrix can be proven to be

$\mathcal{T}^{lmn}=\left(\begin{array}[]{c}\hat{\bf e}_{1}^{lmn}\\\ \hat{\bf
e}_{2}^{lmn}\\\ \hat{\bf e}_{3}^{lmn}\end{array}\right)\,,$ (38)

and its inverse matrix is $\mathcal{T}^{\prime
lmn}=\left(\mathcal{T}^{lmn}\right)^{T}$. In the local frame,
$\mathbf{D}_{B}^{lmn}$ is transformed to
$\mathbf{D^{\prime}}_{B}^{lmn}=\mathcal{T}^{lmn}\mathbf{D}^{lmn}\mathcal{T^{\prime}}^{lmn}$.

Four, analyze the local magnetic structures via calculating the trace,
discriminant, eigenvalues, $\theta_{\mathrm{eig}}$, and $R_{\mathrm{tr}}$ of
the matrix

$\mathbf{M}^{lmn}=\left(\begin{array}[]{cc}D_{B,11}^{\prime
lmn}&D_{B,12}^{\prime lmn}\\\ D_{B,21}^{\prime lmn}&D_{B,22}^{\prime
lmn}\end{array}\right)\,,$ (39)

as discussed in Sect. 2.3.

###  3.3 Step 3: Locate Reconnection Sites of extremal $E_{\parallel}$

As discussed in Sect. 2.2, we further locate the positions with extremal
values of $E_{\parallel}$. Theoretically, they can be identified by
$\nabla_{\perp}E_{\parallel}=0$. However, numerically,
$\nabla_{\perp}E_{\parallel}$ can hardly vanish and it is also difficult to
define a general small-value threshold for arbitrary cases. Thus, to
approximately locate these sites, we directly judge whether $\mathbf{r}^{lmn}$
is a 2D extremal point of $E_{\parallel}$ on the MPP. To implement this, four
grids adjacent with $\mathbf{r}^{lmn}$ (see Fig. 2) on the MPP is defined as,

$\displaystyle\mathbf{R}^{-0}$ $\displaystyle=\mathbf{r}^{lmn}-\Delta
L\hat{\bf e}_{1}^{lmn}\,,$ (40) $\displaystyle\mathbf{R}^{+0}$
$\displaystyle=\mathbf{r}^{lmn}+\Delta L\hat{\bf e}_{1}^{lmn}\,,$ (41)
$\displaystyle\mathbf{R}^{0-}$ $\displaystyle=\mathbf{r}^{lmn}-\Delta
L\hat{\bf e}_{2}^{lmn}\,,$ (42) $\displaystyle\mathbf{R}^{0+}$
$\displaystyle=\mathbf{r}^{lmn}+\Delta L\hat{\bf e}_{2}^{lmn}\,,$ (43)

where $\Delta L$ is the spatial step length of the original discrete data. The
values of $E_{\parallel}$ at these positions can be obtained via linear
interpolation of original data. If the values of $E_{\parallel}$ have the same
signs at $\mathbf{r}$, $\mathbf{R}^{-0}$, $\mathbf{R}^{+0}$,
$\mathbf{R}^{0-}$, and $\mathbf{R}^{0+}$, and $\left|E_{\parallel}\right|$ at
$\mathbf{r}^{lmn}$ is maximum among the five positions, then we identify
$\mathbf{r}^{lmn}$ as a 2D extreme point of $\left|E_{\parallel}\right|$.

Table 2: Key configuration fields of the Parameters structure. Fields | Description | Values
---|---|---
ARD_AnalyzeAllGrids | Skip Step 1 and analyze all grids | 0: false; 1: true
ARD_ScalarThreshold | Threshold value of $\left|E_{\parallel}\right|$ | Nonnegative real numbers
ARD_ShowThresScalarProfile | Draw a histogram of $\left|E_{\parallel}\right|$ without running ARD | 0: false; 1: true
ARD_FixTrace | Fix the trace of $\mathbf{D}^{lmn}_{B}$ from numerical errors | 0: false; 1: true
ARD_AnalyzeLocalEffects | Analyze local reconnection effects (see Appendix A) | 0: false; 1: true
NumRAMBlock | Number of RAM blocks | Positive integer
OutputType | Select format of output file | -1: no output file; 0: .mat; 1: .csv
OutputDir | Output directory | String
OutputLabel | User-defined label on output filename | String
OutputExtraData | Output ExtraData | 0: false; 1: true

##  4 The ARD function in LoRD Toolkit

This method has been implemented as an easy-to-use function ARD, integrated
into our open-source Matlab toolkit project named LoRD (Locate Reconnection
Distribution). LoRD can be freely downloaded by the GitHub link
‘‘https://github.com/RainthunderWYL/LoRD.git’’. To use the functions in LoRD,
one just needs to add the directory into the Matlab environment by the
following command:

⬇

path(’<Download Directory>/LoRD/matlab’,path);

###  4.1 The ARD function

The ARD function can be executed as follows:

⬇

RDInfo = ARD(B1,B2,B3,x1,x2,x3,Parameters,N1,N2,N3);

Here, B1, B2, and B3 are 3D matrices of discrete magnetic field on three
directions. Notice that they follow the “meshgrid” data model of Matlab,
namely, the first dimension is $x^{2}$, the second dimension is $x^{1}$, and
the third dimension is $x^{3}$. x1, x2, and x3 are 1D coordinate arrays on
three directions. At present, ARD can only process uniform Cartesian mesh. But
users can use interpolation to generate uniform Cartesian mesh inputs from
other complex mesh models. N1, N2, and N3 are three optional inputs with the
same size as the magnetic field data, which passes user-defined data of $\bf
N$ to ARD. Without N1–N3, ARD will suppose the nonideal term is induced by a
Joule dissipation with a constant resistivity and use
$\mathbf{J}=\nabla\times\mathbf{B}$ to replace $\bf N$.

Parameters is a structure containing fields of key configuration parameters as
listed in Table 2. If Parameters.ARD_AnalyzeAllGrids is set as 1, ARD skips
Step 1 (see Sect. 3.1) and analyzes all girds, which should be carefully
considered because of the considerable costs of computational resources. The
threshold value of $\left|E_{\parallel}\right|$ is defined by
Parameters.ARD_ScalarThreshold. To help users determine this value, ARD can
present a histogram of $\left|E_{\parallel}\right|$ without proceeding further
analysis if Parameters.ARD_ShowThresScalarProfile is set as 1. Because the
magnetic Jacobian matrix $\mathbf{D}^{lmn}_{B}$ might have small trace
components resulting from numerical errors, ARD provides an option to clean
the trace by
$\mathbf{D}^{lmn}_{B}-1/3tr\left(\mathbf{D}^{lmn}_{B}\right)\mathbf{I}$, which
is controlled by Parameters.ARD_FixTrace. If Parameters.NumRAMBlock is larger
than 1, ARD divides the original field data into NumRAMBlock blocks and
sequentially analyzes each block to avoid RAM overflow, which is useful for
processing massive magnetic field data.

The output RDInfo is a structure containing two fields, namely, RDInfo.Data
and RDInfo.ExtraData. RDInfo.Data is a $N_{\mathrm{sl}}\times 7$ matrix. Each
row saves the data of an analyzed grid and the columns correspond to 7 key
attributes, including x1, x2, x3, RDType, Is2DExtrema, EigAngle, and
RatioMTrace. x1, x2, and x3 are the $x^{1}$, $x^{2}$, and $x^{3}$ coordinates,
respectively. RDType is an integer number taking values from 1 to 9, marking
the types of $\mathbf{B}_{\perp}$ as listed in Table 1. Is2DExtrema labels
whether the grid is a 2D extreme point of $E_{\parallel}$. EigAngle and
RatioMTrace save the values of $\theta_{\mathrm{eig}}$ and $R_{\mathrm{tr}}$,
respectively. RDInfo.ExtraData is an optional output controlled by
Parameters.OutputExtraData, which has 25 columns containing information of
$\left|\mathbf{B}^{lmn}\right|$, $\mathbf{D}^{\prime lmn}_{B}$, $\hat{\bf
e}_{1}^{lmn}$, $\hat{\bf e}_{2}^{lmn}$, $\hat{\bf e}_{3}^{lmn}$, and
parameters for evaluating local reconnection effects (see Table 3). It should
be noticed that Parameters.OutputExtraData is forcibly set as 1 if
Parameters.ARD_AnalyzeLocalEffects is enabled. ARD can also save the output
RDInfo into “.mat” or “.csv” files, which is set by Parameters.OutputType.
Users can customize the output directory by Parameters.OutputDir and also
attach a label behind the filename via Parameters.OutputLabel.

Table 3: Data saved in RDInfo.ExtraData. $\bf\Gamma$ is defined in Appendix A. Column No. | Index Name | Physical Symbol
---|---|---
1 | B0 | $\left|\mathbf{B}\right|$
2 | DB11 | $D^{\prime}_{B,11}$
3 | DB12 | $D^{\prime}_{B,12}$
4 | DB21 | $D^{\prime}_{B,21}$
5 | DB22 | $D^{\prime}_{B,22}$
6 | DB31 | $D^{\prime}_{B,31}$
7 | DB32 | $D^{\prime}_{B,32}$
8 | DB13 | $D^{\prime}_{B,13}$
9 | DB23 | $D^{\prime}_{B,23}$
10 | DB33 | $D^{\prime}_{B,33}$
11 | e11 | $\hat{\mathrm{e}}_{1,1}$
12 | e12 | $\hat{\mathrm{e}}_{1,2}$
13 | e13 | $\hat{\mathrm{e}}_{1,3}$
14 | e21 | $\hat{\mathrm{e}}_{2,1}$
15 | e22 | $\hat{\mathrm{e}}_{2,2}$
16 | e23 | $\hat{\mathrm{e}}_{2,3}$
17 | e31 | $\hat{\mathrm{e}}_{3,1}$
18 | e32 | $\hat{\mathrm{e}}_{3,2}$
19 | e33 | $\hat{\mathrm{e}}_{3,3}$
20 | DNpara1 | $\partial E_{\parallel}/\partial R_{1}$
21 | DNpara2 | $\partial E_{\parallel}/\partial R_{2}$
22 | Gamma1 | $\Gamma_{1}$
23 | Gamma2 | $\Gamma_{2}$
24 | CurlN | $\left|\nabla\times\mathbf{N}\right|$
25 | BxCurlN | $\left|\mathbf{B}\times\left(\nabla\times\mathbf{N}\right)\right|$

###  4.2 Input API

We implement an input API function for loading the binary data files of
discrete fields to simplify the usage of ARD for users unfamiliar with Matlab.
For example, to load an $n_{1}\times n_{2}\times n_{3}$ data $B^{1}$ from a
binary file with “double” precision named “B1.bin”, one can call the following
commands:

⬇

DIM = [n1,n2,n3];

Precision = ’double’;

LIM = [x1s,x1e,x2s,x2e,x3s,x3e];

[B1,x1,x2,x3] = Tool_LoadData_Bin(’B1.bin’,DIM,Precision,LIM);

The binary file should save the 3D data of $B^{1}$ in a 1D array, in which the
$x^{1}$-direction index changes fastest and the $x^{3}$-direction index
changes slowest. Precision is a string setting the class and size of bits to
read which should be consistent with the style of the binary file. Available
parameters of Precision include ‘double’, ‘float32’, etc.. LIM is an optional
parameter defining the coordinate limitations in three directions. Here, x1s
and x1e define the start and end of $x^{1}$-coordinate, respectively, and x2s,
x2e, x3s, and x3e have similar definitions. Without LIM, this function will
simply assign the outputs x1, x2, and x3 with integer grid indices. The output
variables B1, x1, x2, and x3 can be directly passed to the ARD function.

###  4.3 Standards of Code Upgrade

The functions and configuration parameters of the ARD code introduced above
are based on the current version. With the development of our method, the code
might also be upgraded correspondingly. For example, more output parameters
might be defined to give a more comprehensive picture of reconnection sites
and the Matlab scripts will also be optimized for better efficiency. However,
the code upgrade will follow a backward-compatible standard. To be specific,
the I/O APIs will remain unchanged. If necessary, we will add new
configuration parameters into the Parameters structure or output new
reconnection variables in the RDInfo structure. All new features of our code
will be updated in time on the GitHub website in the future.

Figure 3: The distributions of mass density $\rho$ (a), temperature $T$ (b),
parallel current density $\left|J_{\parallel}\right|$ (c), and 3D X/O-type
grids (d) at $t=9$. For $\rho$ and $T$, we only depict their distributions on
the current-sheet middle plane ($y=0$), the front/back planes ($z=\pm 0.2$),
and the left/right planes ($x=\pm 1$). Panel (c) shows the 3D profile of
$\left|J_{\parallel}\right|$ satisfying $\left|J_{\parallel}\right|>80$. In
panel (d), we plot the grids that (1) satisfy $\left|J_{\parallel}\right|>50$,
(2) have local magnetic structures of 3D X and O types, and (3) are 2D
projected extreme points of $J_{\parallel}$. The X-type (type 1) grids are
depicted by blue dots, while the O-type grids (types 2 and 3) are plotted by
orange ones. The magenta and black curves in panels (c) and (d) are examples
of sheared and twisted field lines, which are traced from the initial sampling
positions near 3D X-type and O-type grids, respectively. The green spheres
mark the regions for sampling start positions of field-line tracing. An
animation of this figure showing the entire evolution from $t=0$ to $9$ is
available.  Figure 4: Four snapshots of 3D X/O-type grids exhibiting the
development of turbulent reconnection. The gray dashed lines approximately
trace the origin of turbulent reconnection regions.

Figure 5: The evolutions of grid numbers for different types of structures.
Only grids with $\left|J_{\parallel}\right|>50$ are counted.  Figure 6: The
number densities of $\left|J_{\parallel}\right|$ (a) , $P_{\mathrm{oh}}$ (b),
$\theta_{\mathrm{eig}}$ (c), and $R_{\mathrm{tr}}$ (d) of the 3D X-type grids
at different moments. Different colors correspond to different moments.
Figure 7: The number densities of $\left|J_{\parallel}\right|$ (a),
$P_{\mathrm{oh}}$ (b), $\theta_{\mathrm{eig}}$ (c), and $R_{\mathrm{tr}}$ (d)
of the 3D X-type grids at $t=9$. The blue curves depict profiles of all 3D
X-types grids, while the purple ones plot the subsets with 2D extrema of
$J_{\parallel}$.

##  5 Benchmark 1: 3D Turbulent Reconnection within a Harris Sheet

###  5.1 Numerical Model

To test the validation of our method, we perform a 3D MHD simulation of
magnetic reconnection within a Harris sheet. We solve the resistive MHD
equations:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho{\bf
u}\right)$ $\displaystyle=0\,,$ $\displaystyle\frac{\partial\left(\rho{\bf
u}\right)}{\partial t}+\nabla\cdot\left(\rho{\bf u}{\bf u}-{\bf BB}+P^{*}{\bf
I}\right)$ $\displaystyle=0\,,$ $\displaystyle\frac{\partial e}{\partial
t}+\nabla\cdot\left[\left(e+P^{*}\right){\bf u}-{\bf B}\left(\bf{B}\cdot{\bf
u}\right)\right]$ $\displaystyle=0\,,$ (44) $\displaystyle\frac{\partial{\bf
B}}{\partial t}-\nabla\times\left(\bf{u}\times\bf{B}\right)$
$\displaystyle=-\nabla\times\left(\eta{\bf J}\right)\,,$ $\displaystyle{\bf
J}$ $\displaystyle=\nabla\times{\bf B}\,,$

where, $\rho$, $\bf{u}$, $p$, $\bf{B}$, $T$, and $\bf{J}$ denote plasma
density, velocity, thermal pressure, magnetic field, temperature, and current
density, respectively. $P^{*}$ equals $p+B^{2}/2$ and $\bf{I}$ is the identity
matrix. $e=p/\left(\gamma-1\right)+\rho u^{2}/2+B^{2}/2$ is the total energy,
where $\gamma=5/3$ is the adiabatic index. All physical quantities are
normalized based on the dimensionless units as listed in Table 4.

Table 4: Dimensionless Units Unit of | Symbol | Value
---|---|---
Mass | $\bar{m}$ | $8.36\times 10^{-25}\,\mathrm{g}$
Number density | $n_{0}$ | $2\times 10^{10}\,\mathrm{cm^{-3}}$
Space | $L_{0}$ | $5\times 10^{9}\,\mathrm{cm}$
Magnetic strength | $B_{0}$ | $20\,\mathrm{G}$
Mass density | $\rho_{0}$ | $1.67\times 10^{-14}\,\mathrm{g\,cm^{-3}}$
Time | $t_{0}$ | $114.61\,\mathrm{s}$
Velocity | $u_{0}$ | $4.36\times 10^{7}\,\mathrm{cm\,s^{-1}}$
Temperature | $T_{0}$ | $1.15\times 10^{7}\,\mathrm{K}$
Pressure | $p_{0}$ | $31.8\,\mathrm{dyn\,cm^{-2}}$
Energy density | $e_{0}$ | $31.8\,\mathrm{erg\,cm^{-3}}$
Current density | $J_{0}$ | $9.54\,\mathrm{statC\,s^{-1}\,cm^{-2}}$

The initial magnetic field forms a force-free Harris sheet, namely,

$\displaystyle B_{x}$ $\displaystyle=B_{0}\tanh\left(y/\lambda\right)\,,$
$\displaystyle B_{y}$ $\displaystyle=0\,,$ (45) $\displaystyle B_{z}$
$\displaystyle=B_{0}/\cosh\left(y/\lambda\right)\,,$

where $B_{0}=1$ and $\lambda=0.1$ is the half-width of the current sheet. The
initial mass density and pressure are uniformly set as $1$ and $0.05$,
respectively, corresponding to a $\beta=0.1$. We set a uniform background
resistivity $\eta=5\times 10^{-6}$ to obtain a Lundquist number
$S=L_{0}u_{0}/\eta=2\times 10^{5}$. To trigger fast reconnection, a
perturbation on the $z$-direction of the magnetic vector field is placed at
the center of the current sheet (also see Ye et al., 2020), namely,

$\tilde{A}_{z}=A_{p}\exp\left(-\frac{x^{2}+y^{2}}{2w_{A}^{2}}\right)\,,$ (46)

where, $A_{p}=0.03$ and $w_{A}=0.1$. The initial velocity field is a small-
amplitude Gaussian thermal noise with a zero mean value and a standard
deviation $\sigma_{u}=10^{-3}$, which works as a seed of the self-sustained
turbulence (also see Huang & Bhattacharjee, 2016).

The simulation domain is $x\in\left[-2.5,2.5\right]$,
$y\in\left[-0.5,0.5\right]$, and $z\in\left[-0.2,0.2\right]$. The
$z$-boundaries are periodic, while the rest are open boundaries. The static
mesh refinement technique is applied to implement a uniform mesh in the
reconnection region and save computational costs. The root level-0 grid
numbers are set as $1200$, $240$, and $96$ on $x$, $y$, and $z$ directions,
respectively. The level-1 grid number doubles in three directions in the
region $0.1<\left|y\right|<0.3$. The level-2 grid, doubling again, is set in
the core reconnection region $\left|y\right|\leq 0.1$, implementing an
effective mesh of $4800\times 960\times 384$ there.

We solve the above system using the Athena++ code (Stone et al., 2020). The
HLLD Riemann solver (Miyoshi & Kusano, 2005), the 2-order piecewise linear
method (PLM), and the 2-order van Leer predictor-corrector scheme are selected
for solving the conservative part of the MHD equations. The resistivity term
is calculated by the explicit operator splitting method. The 2-order RKL2
super-time-stepping algorithm is applied to reduce computational costs (Meyer
et al., 2014). The simulation stops at $t=9$.

###  5.2 Evolution Overview

The evolution presents a standard picture of the Harris-sheet reconnection
(see the animation of Fig. 3). In the beginning, triggered by the perturbation
field $\tilde{A}_{z}$, the center of the current sheet first gets thinner,
where the density, temperature, and current also increase simultaneously.
During this stage, the system shows a translate symmetry on $z$-direction and
tends to the 2D results. After $t=3$, 3D tearing-mode instability starts to
dominate the evolution and the reconnection is significantly boosted. Various
small-scale structures carrying strong currents appear and the distributions
of density and temperature get highly nonuniform in three directions. After
$t=6$, the entire current sheet reaches a fully-developed turbulent pattern
wherein the oblique tearing-mode fluctuations can be recognized (also see
Huang & Bhattacharjee, 2016). Finally, at $t=9$, several relatively large-
scale vortex-like structures containing complex density and temperature
profiles grow near the center of the current sheet (see Fig. 3a and b).

###  5.3 Qualitative Global Picture

As discussed in Sect. 2.1, the regions with strong $E_{\parallel}$ can reflect
the locations of general reconnection, which can be equivalently replaced by
$J_{\parallel}$ since we apply a uniform resistivity. After the reconnection
enters the turbulent stage, the regions with strong $J_{\parallel}$ show
patchy-like patterns (see Fig. 3c).

Based on the algorithm in Sect. 3, we obtain the information of magnetic
structures at strong $J_{\parallel}$ regions. In Fig. 3, to give a clear
picture, we only plot 3D X/O-type grids having 2D projected extrema of
$J_{\parallel}$ (see Fig. 3d). The X-type grids are the locations undergoing
the SL reconnection (Priest & Forbes, 1989). Meanwhile, the O-type grids can
approximately present positions of flux ropes consisting of twisted field
lines.

The distributions of X/O-type grids can exhibit the dynamical evolution of
reconnection regions (see Fig. 4). At $t=3$, several X-lines and O-lines start
to emerge near the center region (Fig. 4a). Later, the X-lines are shattered
by the tearing-mode instability and both X and O-type structures with oblique
patterns are generated simultaneously (Fig. 4b). The newborn reconnection
regions keep extending on the $x$-direction and finally develop complex
turbulent states (Fig. 4c, d).

The local magnetic structures on the MPP can also present global features of a
bundle of field lines near a grid. In Fig. 3c and d, we trace several field
lines starting from a local spherical volume near an X-type grid (see the
magenta curves), which present a global sheared pattern. As another example,
the black field lines, initially sampled within a volume containing O-type
grids, show a typical pattern of flux ropes, which corresponds to the vortex
structure also observed in the density and temperature profiles (Fig. 3).
Similar examples of field lines are also exhibited in Fig. 4.

###  5.4 Quantitative Statistical Analysis

To study the generation probabilities of magnetic structures with different
structures, we count the grid numbers of different types (see Fig. 5). At the
early stage, the magnitudes of $\left|J_{\parallel}\right|$ on all grids are
less than the threshold value and thus no grid is recognized. After $t=1$,
type-1 grids (3D X-type) first appear. Later, their number increases rapidly
and finally reaches a platform after $t=6$. The numbers of type-2 and type-3
grids are almost the same, showing a similar evolution trend as type-1 grids.
The growth of O-type grids corresponds to the formation of twisted flux ropes
and thus represents the development of the tearing-mode instability. The
numbers of type-4 and type-5 are also approximately the same but are at least
2-order less than types 1–3. The type-6 grids, with the 3D anti-parallel
lines, are very rare because zero eigenvalues can hardly appear within
numerical data. However, many grids recognized as types 1–5 are close to anti-
parallel lines, because their $\theta_{\mathrm{eig}}$ is close to zero. This
fact explains the reason why we use $\theta_{\mathrm{eig}}$ as a key parameter
of $\bf B_{\perp}$. Finally, although several grids of 2D types are
identified, their number is small enough to ignore.

Now we focus on the statistical properties of type-1 grids, which represent
the evolution rules of the SL reconnection (see Fig. 6). During $t<3$, though
the relative errors of number density profiles are relatively large because of
the small number of recognized grids, some reconnection information can still
be reflected. To be specific, both $\left|J_{\parallel}\right|$ and
$P_{\mathrm{oh}}=\eta J^{2}$ are small corresponding to a moderate
reconnection; $\theta_{\mathrm{eig}}$ mainly distributes below $10^{\circ}$
meaning the reconnection sites tend to have magnetic structures that are close
to anti-parallel; $R_{\mathrm{tr}}$ is close to 0 implying the 3D effect is
weak. With the development of turbulence, the number density curves of
$\left|J_{\parallel}\right|$, $P_{\mathrm{oh}}$, $\theta_{\mathrm{eig}}$, and
$R_{\mathrm{tr}}$ keep extending toward larger values (see Fig. 6). At $t=9$,
the number density of $\left|J_{\parallel}\right|$ first decreases
monotonically as $\left|J_{\parallel}\right|$ increases to 300, then forms a
platform at $\left|J_{\parallel}\right|\in\left[300,400\right]$, and finally
rapidly damps to zero at $\left|J_{\parallel}\right|\sim 440$ (Fig. 6a). The
distribution of $P_{\mathrm{oh}}$ presents a similar trend with
$\left|J_{\parallel}\right|$ (Fig. 6b). $\theta_{\mathrm{eig}}$ mainly
distributes in the region $\theta_{\mathrm{eig}}<40^{\circ}$ and its peak
number density appears near $5^{\circ}$ (Fig. 6c). The maximum value of
$\theta_{\mathrm{eig}}$ reaches $70^{\circ}$ and there are still considerable
grids as $\theta_{\mathrm{eig}}\to 0$. The number density of $R_{\mathrm{tr}}$
decreases monotonically in the entire domain. After 3D effects get
significant, more grids with larger $R_{\mathrm{tr}}$ emerge (Fig. 6d). It
should be noticed that these results might be model-dependent and should be
reconsidered for different simulations. For instance, the values and behaviors
of $\left|J_{\parallel}\right|$ and $P_{\mathrm{oh}}$ might be sensitive to
the magnitude of $\eta$.

We extract the X-type grids with 2D extrema of $J_{\parallel}$ at $t=9$ to
compare its distribution with that of all X-type grids (Fig. 7). The number
densities of $\left|J_{\parallel}\right|$ and $P_{\mathrm{oh}}$ of 2D extreme
grids almost coincide with that of all X-grids for large values (see the
$\left|J_{\parallel}\right|>300$ region in Fig. 7a and the
$P_{\mathrm{oh}}>0.6$ region in Fig. 7b). But for smaller values, they are
lower than the curves of all X-grids (see Fig. 7a and b). For
$\theta_{\mathrm{eig}}$ (Fig. 7c), compared with all X-type grids, the number
density curve of 2D extreme grids has a similar shape but is about one order
lower. The two curves of $R_{\mathrm{tr}}$ have similar trends (Fig. 7d), but
2D extreme grids tend to concentrate within a smaller value region
($R_{\mathrm{tr}}<0.1$), implying that, near the extremal-$\Xi$ lines, the
variations of local field strength on the direction of guide field
($\hat{\mathbf{e}}_{3}$) is weaker.

The statistical results of 3D O-type grids are also provided in Appendix C.
The rules of $\left|J_{\parallel}\right|$, $P_{\mathrm{oh}}$, and
$R_{\mathrm{tr}}$ of 3D O-type grids are similar to 3D X-type grids (see
Figs.13 and 14). The main difference is that more O-type grids with large
$\theta_{\mathrm{eig}}$ near $90^{\circ}$ are identified, corresponding to
flux ropes with poloidal field lines close to circle shapes (see Figs.13c and
14c).

Figure 8: Distributions of 3D X-type (blue dots) and O-type (orange dots)
grids in a data-driven simulation of a large-scale coronal flux rope eruption.
The plotted grids are 2D projected extreme points of $J_{\parallel}$,
satisfying $\left|J_{\parallel}\right|>11.9\,\mathrm{statC\,s^{-1}\,cm^{-2}}$.
The black curves are field lines traced from the initial sampling positions
near the top of the flux rope as marked by the big green sphere, while the
magenta ones are traced from a small region in the eruption current sheet as
marked by the small green sphere. The distribution of $B_{z}$ is overlapped at
the bottom.

##  6 Benchmark 2: Large-scale Magnetic Structures during a Coronal Flux Rope
Eruption

Besides the statistical analysis of local reconnection structures, our method
can also be applied to recognize large-scale magnetic structures within 3D
simulations. Here we use the data from a data-driven simulation of an erupting
coronal flux rope by Guo et al. (2023) to test our method. Guo et al. (2023)
solved the full MHD equations in which a realistic coronal physical
environment is considered. The simulation initially set a flux rope comparable
with observations and proceeded under the constraint of a bottom condition
obtained from observed magnetic and velocity fields. The simulation reproduced
the observational characteristics of the X1.0 flare on 2021 October 28.
Detailed simulation configurations are introduced in Guo et al. (2023).

We apply our method to analyze the magnetic field data at 15:32 UT. At this
moment, a major flux rope is rising and an eruption current sheet forms under
the flux rope (see Guo et al., 2023, Fig. 3d). As shown by orange dots and
black field lines in Fig. 8, the O-type (types 2 and 3) grids outline the
position of the major flux rope, which is consistent with the results of Guo
et al. (2023). Our method also precisely locates the position of the eruption
current sheet (see Fig. 8), wherein sheared field lines form a current sheet
that plays an important role in driving the eruption of the flux rope. If we
trace a bundle of field lines starting from a local region in the current
sheet, these field lines show a locally sheared pattern and then extend to the
outer layer of the flux rope (see the magenta curves in Fig. 8).

Figure 9: The distributions of 3D X and O-type grids within a 2.5D Harris-
sheet reconnection. The 3D magnetic field has $z$-symmetry and is constructed
from a 2.5D simulation with the same configurations as the 3D Harris sheet in
Sect. 5. This figure shows the result at $t=4$. In panel (a), the blue dots
have the same definition as in Figs. 3 and 4. The black “x” and red “o”
markers depict the positions of 2D X and O points on the $x$-$y$ plane as
determined by the 2D method, while the black and red lines are the
corresponding X and O lines in 3D. The principal reconnection site with
strongest $\left|J_{\parallel}\right|$ is marked by a magenta “x” marker. On
the plane $z=0.2$, the profile of $u_{x}$ is depicted while the magnetic field
lines are shown by gray curves. Panels (b) and (c) are similar to panel (a)
but exhibit results with different thresholds of $\left|J_{\parallel}\right|$
and $\theta_{\mathrm{eig}}$. Panel (d) has the same thresholds as (c) but also
plots 3D O-type grids with orange dots.  Figure 10: The number densities of
$R_{\mathbf{B}_{\perp}}$ (blue curves), $R_{B_{\parallel}}$ (purple curves),
$R_{\tilde{\mathbf{B}}}$ (green curves), and $R_{\nabla\mathbf{B}}$ (black
dashed curves) of all grids with $\left|J_{\parallel}\right|>50$ at different
moments in the Harris-sheet simulation.

##  7 Summary and Discussion

To conclude, we develop an efficient method for analyzing 3D magnetic
reconnection from discrete simulation data. Based on the general reconnection
theory by Hesse & Schindler (1988), we discuss the usage of the nonideal
electric field $E_{\parallel}$ in representing the reconnection sites and
propose that the locations of $\nabla_{\perp}E_{\parallel}=0$ can present the
distribution of extremal-$\Xi$ lines. We perform theoretical analysis on the
local magnetic structure and provide a complete classification including nine
types of magnetic structures (Table 1). By defining two parameters,
$\theta_{\mathrm{eig}}$ and $R_{\mathrm{tr}}$, the geometric properties of the
local magnetic field can be further clarified. We construct an efficient
numerical method, which only performs algebraical manipulations on the
discrete magnetic field, avoiding computationally expensive operations like
field-line tracing and root-finding. This method directly outputs the
classification and geometric parameters of local magnetic fields at arbitrary
grids, which has been implemented in Matlab language and can be freely
obtained on GitHub.

Through two typical numerical benchmarks, namely the 3D Harris-sheet turbulent
reconnection and the coronal flux rope eruption, we show that this method can
precisely identify the local structures of discrete magnetic field (Fig. 12).
Through the strength of the nonideal electric field and the geometric
attributes of magnetic field, the local structures of reconnection sites can
be effectively revealed. We exhibit not only qualitative pictures but also
quantitative statistical results of the 3D turbulent reconnection by use of
the outputs of our method. It is also shown that macro-scale magnetic
structures such as flux ropes and eruption current sheets can be recognized by
our method.

As a generalization of the 2D method for locating reconnection sites, our
method can also determine the classic 2D X/O lines. To clarify this, we run a
2.5D simulation of the Harris-sheet reconnection in Sect. 5 and copy the
simulation data along $z$-grids to construct 3D data with translation symmetry
on $z$-direction. From the perspective of 2D methods, one naturally chooses
the $z$-direction as the “guide field” direction and the $x$-$y$ plane as a
“global” projection plane, on which the 2D null points with X or O-type
structures can be located (see Fig. 9 and also see Huang & Bhattacharjee,
2010; Shen et al., 2011; Lynch et al., 2016; Ye et al., 2019; Wang et al.,
2021). By using our 3D method, the 3D X and O-type grids can also be located
(Fig. 9). Because the fixed global projection plane in the 2D method has been
generalized as the spatially variable local MPPs, more X-type points can be
recognized. After imposing the constraint of 2D extreme $J_{\parallel}$ and
further selecting grids with larger $\left|J_{\parallel}\right|$ and
$\theta_{\mathrm{eig}}$, we find that the resultant distribution of 3D X-grids
gets almost consistent with the 2D X-lines (see Fig. 9a–c). Especially, the
X-line of the principal reconnection site can be well resolved (see the
magenta “x” markers in Fig. 9). Moreover, the 3D O-grids recognized by our
method outline a flux-rope volume full of twisted field lines, while the 2D
method can only locate its axis (O-line) defined on the $x$-$y$ plane.

Different from the global analysis techniques requiring field-line tracing,
our method performs local analysis on discrete grids. Despite the promising
performances in determining global magnetic structures, field-line tracing
techniques can hardly be applied to the analysis of massive discrete data of
turbulence mainly for two reasons. First, unlike the laminar parts of a
magnetic field, the turbulent parts have no clear global topological patterns
but are full of multi-scale chaotic structures. Second, turbulent structures
with smaller scales correspond to larger local gradients which are more
sensitive to numerical truncation errors. The numerical integration for
tracing field lines can further amplify the influences of local errors by
global iterations, which makes the results unreliable.

The technique proposed by Lapenta (2021) is also a local method, which defines
a local frame based on the local electric field $\mathbf{E}$ and electric
drift speed. Because the Lorentz boost speed to eliminate the magnetic field
perpendicular to $\bf E$, namely $\mathbf{v}_{L}=c^{2}\bf{E\times B/E^{2}}$,
can only be lower than the speed of light $c$ in the adjacent region of
reconnection sites, Lapenta (2021) suggested that 3D reconnection sites can be
located by finding locations with $v_{L}<c$, which has been verified within
kinetic simulation data (Lapenta, 2023). But for MHD simulations, the grid
sizes, typically several orders larger than the inertial scales of electrons
and ions, cannot resolve the absolute value of electric field within the
dissipation regions. As a result, when we applied this method to MHD data, no
reconnection site with $v_{L}<c$ can be recognized. In contrast, our method
mainly analyzes the magnetic field and only uses the relative strength of
nonideal electric field as an auxiliary parameter, which is thus independent
of the system scales. Therefore, besides MHD simulations, our method can also
be potentially applied to kinetic simulations, plasma experiments, and even
_in-situ_ observations.

Now we discuss the limitations of our method. First, when analyzing local
magnetic structures, our method depends on a finite magnetic field to
establish a local frame, which gets invalid at magnetic null points. For
simulation data, the magnetic strength on grids can hardly be zeros, which
makes it safe to apply our method. To obtain complete information on 3D
reconnection sites, one still needs to run a 3D null-point analyzing program,
which has been implemented as the ANP function in the LoRD toolkit. However,
because 3D reconnection takes place in the vicinity of a null point (not just
at the null point) where current density and $E_{\parallel}$ concentrate (see
Pontin et al., 2007), our method is still useful near null points.

Second, when analyzing the local magnetic structures near a grid
$\mathbf{r}^{lmn}$, we approximate the surface normal to magnetic field by a
plane (MPP) within a finite volume with the length scale of $\Delta L$, which
validates better if $R_{\nabla\mathbf{B}}\equiv\Delta
L\left\|\mathbf{D}^{\prime lmn}_{B}\right\|/\left|\mathbf{B}^{lmn}\right|\ll
1$. Meanwhile, we mainly use the perpendicular magnetic field on the MPP
($\mathbf{B}_{\perp}$) to determine the local magnetic structures which might
also be affected by $B_{\parallel}$ and $\tilde{\mathbf{B}}$. In applications,
we define three ratios to evaluate the importance of the three terms, namely,
$R_{\mathbf{B}_{\perp}}\equiv\Delta
L\left\|\mathbf{M}^{lmn}\right\|/\left|\mathbf{B}^{lmn}\right|$,
$R_{B_{\parallel}}\equiv\Delta L\sqrt{\left(D^{\prime
lmn}_{B,31}\right)^{2}+\left(D^{\prime
lmn}_{B,32}\right)^{2}}/\left|\mathbf{B}^{lmn}\right|$, and
$R_{\tilde{\mathbf{B}}}\equiv\Delta L\sqrt{\left(D^{\prime
lmn}_{B,13}\right)^{2}+\left(D^{\prime lmn}_{B,23}\right)^{2}+\left(D^{\prime
lmn}_{B,33}\right)^{2}}/\left|\mathbf{B}^{lmn}\right|$. In the Harris-sheet
simulation, for most grids, the magnitude of $R_{\nabla\mathbf{B}}$ is $\sim
0.1$ which implies the MPP approximation is acceptable (see the black dashed
curves in Fig. 10). Moreover, the distribution of $R_{\mathbf{B}_{\perp}}$
almost overlaps with that of $R_{\nabla\mathbf{B}}$ (see the blue curves in
Fig. 10), the peak of $R_{B_{\parallel}}$ is about one order smaller than
$R_{\mathbf{B}_{\perp}}$, and the magnitude of $R_{\tilde{\mathbf{B}}}$ is
much smaller even compared with $R_{B_{\parallel}}$. Therefore, in this
simulation, the variation of magnetic field in a cell is dominated by
$\mathbf{B}_{\perp}$ for most reconnection grids. $R_{\nabla\mathbf{B}}$,
$R_{\mathbf{B}_{\perp}}$, $R_{B_{\parallel}}$, and $R_{\tilde{\mathbf{B}}}$
can be directly calculated by the 1–10 columns of RDInfo.ExtraData (see Table
3), which is helpful for users to evaluate the validations of MPP
approximation and the identified local magnetic structures. We have not added
the four ratios to the output of ARD for the sake of saving storage space.

Third, our method provides the reconnection properties of the finest objects
in discrete data (i.e., the grids), which, however, cannot determine the
global properties of small-scale reconnection regions containing a series of
adjacent grids. Based on the information of grids, it is possible to further
recognize local finite-volume reconnection regions through clumping
algorithms, which might produce more useful information like volume-integrated
parameters and topological properties. We will attempt to add these functions
in the future to make the code more powerful in dealing with multi-scale
reconnection processes.

###### Acknowledgements.

The authors would like to thank the anonymous referee for constructive
suggestions and acknowledge J. Chen, P. Fan, and C. Xing for their helpful
discussions. The simulations in this paper were performed in the cluster
system of the High Performance Computing Center (HPCC) of Nanjing University.
This research is supported by the Strategic Priority Research Program of the
Chinese Academy of Sciences (Grant No. XDB0560000), by the National Key R&D
Program of China (Grant No. 2021YFA1600504), and by the Natural Science
Foundation of China (Grant No. 12127901).

## References

  * Aulanier et al. (2006) Aulanier, G., Pariat, E., Démoulin, P., & Devore, C. R. 2006, Sol. Phys., 238, 347
  * Comisso & Sironi (2022) Comisso, L. & Sironi, L. 2022, Astrophys. J. Lett., 936, L27
  * Demoulin et al. (1997) Demoulin, P., Bagala, L. G., Mandrini, C. H., Henoux, J. C., & Rovira, M. G. 1997, Astron. Astrophys., 325, 305
  * Demoulin et al. (1996) Demoulin, P., Henoux, J. C., Priest, E. R., & Mandrini, C. H. 1996, Astron. Astrophys., 308, 643
  * Démoulin et al. (1996) Démoulin, P., Priest, E. R., & Lonie, D. P. 1996, J. Geophys. Res. Sp. Phys., 101, 7631
  * Dong et al. (2022) Dong, C., Wang, L., Huang, Y.-M., et al. 2022, Sci. Adv., 8, eabn7627
  * Greco et al. (2017) Greco, A., Matthaeus, W. H., Perri, S., et al. 2017, Space Sci. Rev., 214, 1
  * Guo et al. (2023) Guo, J. H., Ni, Y. W., Zhong, Z., et al. 2023, Astrophys. J. Suppl. Ser., 266, 3
  * Hada et al. (2003) Hada, T., Koga, D., & Yamamoto, E. 2003, Space Sci. Rev., 107, 463
  * Haynes & Parnell (2007) Haynes, A. L. & Parnell, C. E. 2007, Phys. Plasmas, 14, 82107
  * Haynes & Parnell (2010) Haynes, A. L. & Parnell, C. E. 2010, Phys. Plasmas, 17, 92903
  * Hesse et al. (2005) Hesse, M., Forbes, T. G., & Birn, J. 2005, Astrophys. J., 631, 1227
  * Hesse & Schindler (1988) Hesse, M. & Schindler, K. 1988, J. Geophys. Res. Sp. Phys., 93, 5559
  * Huang & Bhattacharjee (2010) Huang, Y. M. & Bhattacharjee, A. 2010, Phys. Plasmas, 17, 062104
  * Huang & Bhattacharjee (2016) Huang, Y.-M. & Bhattacharjee, A. 2016, Astrophys. J., 818, 20
  * Isliker et al. (2019) Isliker, H., Archontis, V., & Vlahos, L. 2019, Astrophys. J., 882, 57
  * Ji et al. (2022) Ji, H., Daughton, W., Jara-Almonte, J., et al. 2022, Nat. Rev. Phys., 4, 263
  * Komar et al. (2013) Komar, C. M., Cassak, P. A., Dorelli, J. C., Glocer, A., & Kuznetsova, M. M. 2013, J. Geophys. Res. Sp. Phys., 118, 4998
  * Kowal et al. (2020) Kowal, G., Falceta-Gonçalves, D. A., Lazarian, A., & Vishniac, E. T. 2020, Astrophys. J., 892, 50
  * Lapenta (2021) Lapenta, G. 2021, Astrophys. J., 911, 147
  * Lapenta (2023) Lapenta, G. 2023, Nat. Phys.
  * Lau & Finn (1990) Lau, Y.-T. & Finn, J. M. 1990, Astrophys. J., 350, 672
  * Li et al. (2022) Li, H. T., Cheng, X., Guo, J. H., et al. 2022, Astron. Astrophys., 663
  * Li et al. (2021a) Li, T., Priest, E., & Guo, R. 2021a, Proc. R. Soc. A Math. Phys. Eng. Sci., 477, 20200949
  * Li et al. (2021b) Li, T. C., Liu, Y.-H., & Qi, Y. 2021b, Astrophys. J. Lett., 909, L28
  * Li et al. (2023) Li, T. C., Liu, Y.-H., Qi, Y., & Zhou, M. 2023, Phys. Rev. Lett., 131, 85201
  * Lynch et al. (2016) Lynch, B. J., Edmondson, J. K., Kazachenko, M. D., & Guidoni, S. E. 2016, Astrophys. J., 826, 43
  * Meyer et al. (2014) Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2014, J. Comput. Phys., 257, 594
  * Miyoshi & Kusano (2005) Miyoshi, T. & Kusano, K. 2005, J. Comput. Phys., 208, 315
  * Parnell et al. (2010a) Parnell, C. E., Haynes, A. L., & Galsgaard, K. 2010a, J. Geophys. Res. Sp. Phys., 115
  * Parnell et al. (2010b) Parnell, C. E., Maclean, R. C., & Haynes, A. L. 2010b, Astrophys. J. Lett., 725, L214
  * Parnell et al. (1996) Parnell, C. E., Smith, J. M., Neukirch, T., & Priest, E. R. 1996, Phys. Plasmas, 3, 759
  * Pontin et al. (2007) Pontin, D. I., Bhattacharjee, A., & Galsgaard, K. 2007, Phys. Plasmas, 14, 52106
  * Pontin & Priest (2022) Pontin, D. I. & Priest, E. R. 2022, Living Rev. Sol. Phys., 19, 1
  * Priest & Forbes (2000) Priest, E. & Forbes, T. 2000, Magnetic Reconnection: MHD Theory and Applications (Cambridge: Cambridge University Press)
  * Priest & Démoulin (1995) Priest, E. R. & Démoulin, P. 1995, J. Geophys. Res. Sp. Phys., 100, 23443
  * Priest & Forbes (1989) Priest, E. R. & Forbes, T. G. 1989, Sol. Phys., 119, 211
  * Priest et al. (2003) Priest, E. R., Hornig, G., & Pontin, D. I. 2003, J. Geophys. Res. Sp. Phys., 108
  * Reid et al. (2020) Reid, J., Parnell, C. E., Hood, A. W., & Browning, P. K. 2020, Astron. Astrophys., 633, A92
  * Schindler et al. (1988) Schindler, K., Hesse, M., & Birn, J. 1988, J. Geophys. Res. Sp. Phys., 93, 5547
  * Shen et al. (2011) Shen, C., Lin, J., & Murphy, N. A. 2011, Astrophys. J., 737, 14
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, Astrophys. J. Suppl. Ser., 249, 4
  * Titov (2007) Titov, V. S. 2007, Astrophys. J., 660, 863
  * Titov et al. (2009) Titov, V. S., Forbes, T. G., Priest, E. R., Mikić, Z., & Linker, J. A. 2009, Astrophys. J., 693, 1029
  * Titov et al. (2002) Titov, V. S., Hornig, G., & Démoulin, P. 2002, J. Geophys. Res. Sp. Phys., 107, SSH 3
  * Vlahos & Isliker (2023) Vlahos, L. & Isliker, H. 2023, Phys. Plasmas, 30, 40502
  * Wang et al. (2023) Wang, Y., Cheng, X., Ding, M., et al. 2023, Astrophys. J. Lett., 954, L36
  * Wang et al. (2021) Wang, Y., Cheng, X., Ding, M., & Lu, Q. 2021, Astrophys. J., 923, 227
  * Wilmot-Smith & Priest (2007) Wilmot-Smith, A. L. & Priest, E. R. 2007, Phys. Plasmas, 14, 102903
  * Wyper & Hesse (2015) Wyper, P. F. & Hesse, M. 2015, Phys. Plasmas, 22, 42117
  * Ye et al. (2020) Ye, J., Cai, Q., Shen, C., et al. 2020, Astrophys. J., 897, 64
  * Ye et al. (2023) Ye, J., Raymond, J. C., Mei, Z., et al. 2023, Astrophys. J., 955, 88
  * Ye et al. (2019) Ye, J., Shen, C., Raymond, J. C., Lin, J., & Ziegler, U. 2019, Mon. Not. R. Astron. Soc., 482, 588
  * Zhang et al. (2022) Zhang, P., Chen, J., Liu, R., & Wang, C. 2022, Astrophys. J., 937, 26
  * Zhang et al. (2021) Zhang, Q., Guo, F., Daughton, W., Li, H., & Li, X. 2021, Phys. Rev. Lett., 127, 185101
  * Zhdankin et al. (2013) Zhdankin, V., Uzdensky, D. A., Perez, J. C., & Boldyrev, S. 2013, Astrophys. J., 771, 124
  * Zhong et al. (2021) Zhong, Z., Guo, Y., & Ding, M. D. 2021, Nat. Commun., 12, 1

##  Appendix A Analysis of Local Reconnection Effects

To quantitatively compare the $E_{\parallel}$ and $\mathbf{N}_{\perp}$ terms
in Eq. 5, here we deduce a different form of this equation. The basic idea is
transforming the RHS of Eq. 5 from the flux coordinate frame to the local
orthogonal and normalized frame spanned by $\hat{\bf e}_{1}$, $\hat{\bf
e}_{2}$, and $\hat{\bf e}_{3}$. In this section, we use
$\mathbf{x}=\left(x^{1},x^{2},x^{3}\right)^{T}$ and
$\mathbf{x}^{\prime}=\left(x^{\prime 1},x^{\prime 2},x^{\prime
3}\right)\equiv\left(\alpha,\beta,s\right)^{T}$ to denote the coordinates in
the $\left(\hat{\bf e}_{1},\hat{\bf e}_{2},\hat{\bf e}_{3}\right)$ and
$\left(\nabla\alpha,\nabla\beta,\hat{\bf b}\right)$ frames, respectively.
Notice that $\hat{\bf e}_{3}=\hat{\bf b}$, $x^{3}=x^{\prime 3}=s$, and
$\partial/\partial s=\partial/\partial x^{3}$.

First, based on the vector relations

$\displaystyle x^{1}\hat{\bf e}_{1}$ $\displaystyle+x^{2}\hat{\bf
e}_{2}=\alpha\nabla\alpha+\beta\nabla\beta\,,$ (47)
$\displaystyle\nabla\alpha$ $\displaystyle=\left(\hat{\bf
e}_{1}\cdot\nabla\alpha\right)\hat{\bf e}_{1}+\left(\hat{\bf
e}_{2}\cdot\nabla\alpha\right)\hat{\bf e}_{2}\,,$ (48)
$\displaystyle\nabla\beta$ $\displaystyle=\left(\hat{\bf
e}_{1}\cdot\nabla\beta\right)\hat{\bf e}_{1}+\left(\hat{\bf
e}_{2}\cdot\nabla\beta\right)\hat{\bf e}_{2}\,,$ (49)

the transformation matrices between the two frames can be obtained as

$\displaystyle\frac{\partial\mathbf{x}}{\partial\mathbf{x}^{\prime}}$
$\displaystyle=\left(\begin{array}[]{ccc}\hat{\bf
e}_{1}\cdot\nabla\alpha,&\hat{\bf e}_{1}\cdot\nabla\beta,&0\\\ \hat{\bf
e}_{2}\cdot\nabla\alpha,&\hat{\bf e}_{2}\cdot\nabla\beta,&0\\\ 0,&0,&1\\\
\end{array}\right)\,,$ (53)
$\displaystyle\frac{\partial\mathbf{x}^{\prime}}{\partial\mathbf{x}}$
$\displaystyle=\frac{1}{B}\left(\begin{array}[]{ccc}\hat{\bf
e}_{2}\cdot\nabla\beta,&-\hat{\bf e}_{1}\cdot\nabla\beta,&0\\\ -\hat{\bf
e}_{2}\cdot\nabla\alpha,&\hat{\bf e}_{1}\cdot\nabla\alpha,&0\\\ 0,&0,&B\\\
\end{array}\right)\,,$ (57)

where $B=\hat{\bf b}\cdot\left(\nabla\alpha\times\nabla\beta\right)$ is the
magnetic strength at the origin point.

Second, we rewrite the $\partial E_{\parallel}/\partial\beta$ and $\partial
E_{\parallel}/\partial\alpha$ terms in Eq. 5. Considering that
$\mathbf{E}\cdot\hat{\bf b}=\mathbf{E}\cdot\hat{\bf e}_{3}$, $E_{\parallel}$
is invariant in the new frame. According to Eq. 53, we have
$\partial/\partial\alpha=\left(\hat{\bf
e}_{1}\cdot\nabla\alpha\right)\partial/\partial x^{1}+\left(\hat{\bf
e}_{2}\cdot\nabla\alpha\right)\partial/\partial x^{2}$ and
$\partial/\partial\beta=\left(\hat{\bf
e}_{1}\cdot\nabla\beta\right)\partial/\partial x^{1}+\left(\hat{\bf
e}_{2}\cdot\nabla\beta\right)\partial/\partial x^{2}$, and thus

$\displaystyle\frac{\partial E_{\parallel}}{\partial\alpha}$
$\displaystyle=\frac{\partial E_{\parallel}}{\partial x^{1}}\hat{\bf
e}_{1}\cdot\nabla\alpha+\frac{\partial E_{\parallel}}{\partial x^{2}}\hat{\bf
e}_{2}\cdot\nabla\alpha=\nabla E_{\parallel}\cdot\nabla\alpha\,,$ (58)
$\displaystyle\frac{\partial E_{\parallel}}{\partial\beta}$
$\displaystyle=\frac{\partial E_{\parallel}}{\partial x^{1}}\hat{\bf
e}_{1}\cdot\nabla\beta+\frac{\partial E_{\parallel}}{\partial x^{2}}\hat{\bf
e}_{2}\cdot\nabla\beta=\nabla E_{\parallel}\cdot\nabla\beta\,.$ (59)

Third, we rewrite the $\partial N^{\beta}/\partial s$ and $\partial
N^{\alpha}/\partial s$ terms in Eq. 5. According to Eq. 57, we have

$\displaystyle N^{\alpha}$ $\displaystyle=\frac{1}{B}\left(N^{1}\hat{\bf
e}_{2}\cdot\nabla\beta-N^{2}\hat{\bf e}_{1}\cdot\nabla\beta\right)\,,$ (60)
$\displaystyle N^{\beta}$ $\displaystyle=-\frac{1}{B}\left(N^{1}\hat{\bf
e}_{2}\cdot\nabla\alpha-N^{2}\hat{\bf e}_{1}\cdot\nabla\alpha\right)\,.$ (61)

Therefore,

$\displaystyle\frac{\partial N^{\alpha}}{\partial s}=$
$\displaystyle\left[\hat{\bf e}_{2}\frac{\partial}{\partial
x^{3}}\left(\frac{N^{1}}{B}\right)-\hat{\bf e}_{1}\frac{\partial}{\partial
x^{3}}\left(\frac{N^{2}}{B}\right)\right]\cdot\nabla\beta$
$\displaystyle+\left(\frac{N^{1}}{B}\frac{\partial\hat{\bf e}_{2}}{\partial
x^{3}}-\frac{N^{2}}{B}\frac{\partial\hat{\bf e}_{1}}{\partial
x^{3}}\right)\cdot\nabla\beta$ $\displaystyle+\frac{1}{B}\left(N^{1}\hat{\bf
e}_{2}-N^{2}\hat{\bf e}_{1}\right)\cdot\frac{\partial\nabla\beta}{\partial
x^{3}}\,.$ (62)

Using $\nabla\left(\hat{\bf b}\cdot\nabla\beta\right)=0$, it can be proved
that $\partial\nabla\beta/\partial x^{3}=-\nabla\hat{\bf
b}\cdot\nabla\beta=-\nabla\hat{\bf e}_{3}\cdot\nabla\beta$ and similarly
$\partial\nabla\alpha/\partial x^{3}=-\nabla\hat{\bf e}_{3}\cdot\nabla\alpha$.
As a result, we can define a vector

$\displaystyle{\bf\Gamma}=$ $\displaystyle\hat{\bf
e}_{2}\frac{\partial\left(N^{1}/B\right)}{\partial x^{3}}-\hat{\bf
e}_{1}\frac{\partial\left(N^{2}/B\right)}{\partial
x^{3}}+\frac{N^{1}}{B}\frac{\partial\hat{\bf e}_{2}}{\partial x^{3}}$
$\displaystyle-\frac{N^{2}}{B}\frac{\partial\hat{\bf e}_{1}}{\partial
x^{3}}+\left(-\frac{N^{1}}{B}\hat{\bf e}_{2}+\frac{N^{2}}{B}\hat{\bf
e}_{1}\right)\cdot\nabla\hat{\bf e}_{3}\,,$ (63)

and write Eq. 62 as $\partial N^{\alpha}/\partial
s={\bf\Gamma}\cdot\nabla\beta$. Similarly, we have $\partial
N^{\beta}/\partial s=-{\bf\Gamma}\cdot\nabla\alpha$.

Finally, Eq. 5 becomes

$\displaystyle\frac{\partial\dot{\alpha}}{\partial x^{3}}$
$\displaystyle=\nabla
E_{\parallel}\cdot\nabla\beta+{\bf\Gamma}\cdot\nabla\alpha\,$ (64a)
$\displaystyle\frac{\partial\dot{\beta}}{\partial x^{3}}$
$\displaystyle=-\nabla
E_{\parallel}\cdot\nabla\alpha-{\bf\Gamma}\cdot\nabla\beta.$ (64b)

Because both $\nabla\alpha$ and $\nabla\beta$ are perpendicular to $\hat{\bf
e}_{3}$, only the perpendicular components of $\nabla E_{\parallel}$ and
$\bf\Gamma$, namely, $\nabla_{\perp}E_{\parallel}$ and $\bf\Gamma_{\perp}$,
have effects on the local line conservation. Considering that
$\nabla_{\perp}E_{\parallel}$ and $\bf\Gamma_{\perp}$ can be explicitly
calculated in the $\left(\hat{\bf e}_{1},\hat{\bf e}_{2},\hat{\bf
e}_{3}\right)$ frame, we can use them to qualitatively compare $E_{\parallel}$
and $\bf N_{\perp}$ terms even without knowing $\nabla\alpha$ and
$\nabla\beta$. However, because an absolute comparison still relies on the
flux coordinate frame, $\nabla_{\perp}E_{\parallel}$ and $\bf\Gamma_{\perp}$
only provide a statistical perspective of approximate comparison.

For instance, we compare $\nabla_{\perp}E_{\parallel}$ and $\bf\Gamma_{\perp}$
in the Harris-sheet simulation shown in Sect. 5 (see Fig. 11). At $t=9$, the
peak of $\left|\nabla_{\perp}E_{\parallel}\right|$ locates near $0.1$ while
$\left|\bf\Gamma_{\perp}\right|$ mainly distributes near $0.01$ (see the blue
curves in Fig. 11a). Moreover, during the entire evolution, both the mean
values and standard deviations of $\left|\nabla_{\perp}E_{\parallel}\right|$
are significantly larger than that of $\bf\Gamma_{\perp}$ (see Fig. 11c).
These results show that, statistically speaking, the magnitude of
$\nabla_{\perp}E_{\parallel}$ is larger than $\bf\Gamma_{\perp}$. On the other
hand, we depict the number density of the angle $\theta_{le}$ spanned by
$\nabla_{\perp}E_{\parallel}$ and $\bf\Gamma_{\perp}$, which shows that
$\theta_{le}$ is most likely to be $90^{\circ}$, namely, the two vectors tend
to normal to each other for most of grids (see Fig. 11b). During the
evolution, the mean value of $\theta_{le}$ keeps close to $90^{\circ}$ and its
standard deviation is smaller than that of a uniform distribution, which
implies that the distributions of $\theta_{le}$ tend to concentrate near
$90^{\circ}$.

It should also be noticed that the grids with extremal $E_{\parallel}$
identified by the algorithm introduced in Sect. 3 compose a subset about
$12\%$ of all grids and have smaller values of
$\left|\nabla_{\perp}E_{\parallel}\right|$ (see the purple solid curve in Fig.
11a). Therefore, the algorithm of locating 2D extremal $E_{\parallel}$ is a
good numerical approximation of $\nabla_{\perp}E_{\parallel}=0$.

When using the ARD code, users can enable the calculation of
$\nabla_{\perp}E_{\parallel}$ and $\bf\Gamma_{\perp}$ by setting
ARD_AnalyzeLocalEffects as 1, which costs more RAM and computation resources.
All the derivatives in the expressions of $\nabla_{\perp}E_{\parallel}$ and
$\bf\Gamma_{\perp}$ are first calculated in the original Cartesian frame by
the 2-order central difference method and then transformed into the local
frame by the transformation matrix $\mathcal{T}^{lmn}$ in Eq. 38. ARD also
outputs $\left|\nabla\times\mathbf{N}\right|$ and $\left|\bf
B\times\left(\nabla\times N\right)\right|$ to study the local flux and line
conservations (see Table 3).

Figure 11: Statistical comparisons of $\nabla_{\perp}E_{\parallel}$ and
$\bf\Gamma_{\perp}$ in the Harris-sheet simulation. Panel (a) plots the number
densities of $\left|\nabla_{\perp}E_{\parallel}\right|$ (solid curves) and
$\left|\bf\Gamma_{\perp}\right|$ (dashed curves) at $t=9$. Panel (b) shows the
number densities of the angle $\theta_{le}$ spanned by $\bf\Gamma_{\perp}$ and
$\nabla_{\perp}E_{\parallel}$ at $t=9$. In panels (a) and (b), the blue curves
depict all grids satisfying $\left|J_{\parallel}\right|\geq 50$ while the
purple ones plot the subsets with 2D extremal $E_{\parallel}$. Panel (c)
depicts the time evolutions of the mean values of
$\left|\nabla_{\perp}E_{\parallel}\right|$ (the blue solid curve) and
$\left|\bf\Gamma_{\perp}\right|$ (the orange solid curve). The blue and orange
dashed curves plot the standard deviations of
$\left|\nabla_{\perp}E_{\parallel}\right|$ and
$\left|\bf\Gamma_{\perp}\right|$, respectively. Panel (d) exhibits the time
evolutions of the mean values of $\theta_{le}$ (the solid curve) and its
standard deviation (the dashed curve). The black dashed line denotes the
theoretical standard deviation of a uniform random distribution sampled from
$0^{\circ}$ to $180^{\circ}$. The mean values and standard deviations in
panels (c) and (d) are calculated from all grids satisfying
$\left|J_{\parallel}\right|\geq 50$.

##  Appendix B Examples of Magnetic Structures on the MPP

In Fig. 12, we exhibit the $\mathbf{B}_{\perp}$ lines of types 1–6 as found in
the simulation of 3D Harris-sheet reconnection, which proves that our method
can precisely capture the local magnetic structures.

Figure 12: Six typical local field-line structures of $\mathbf{B}_{\perp}$
(types 1 to 6) as found in the reconnection region at $t=9$. Panels (a)–(f)
depict the field lines of $\mathbf{B}_{\perp}$ on the MPP of selected grid
positions as labeled by black triangle markers in panel (g). Different colors
correspond to different types. The distributions of
$B_{\perp}=\left|\mathbf{B}_{\perp}\right|$ are also plotted with a gray color
map. In panel (g), to give a clear picture, we only plot the grids with
$\theta_{\mathrm{eig}}$ larger than a threshold value. For types 1–3, the
threshold is set as $30^{\circ}$; for types 4–5, the threshold is $5^{\circ}$.

##  Appendix C Quantitative Statistical Results of O-Type Grids in the Harris-
sheet Simulation

Figures 13 and 14 plot the same statistical results as Figs. 6 and 7 but for
3D O-type grids.

Figure 13: The number densities of $\left|J_{\parallel}\right|$ (a) ,
$P_{\mathrm{oh}}$ (b), $\theta_{\mathrm{eig}}$ (c), and $R_{\mathrm{tr}}$ (d)
of the 3D O-type grids at different moments. Different colors correspond to
different moments.  Figure 14: The number densities of
$\left|J_{\parallel}\right|$ (a), $P_{\mathrm{oh}}$ (b),
$\theta_{\mathrm{eig}}$ (c), and $R_{\mathrm{tr}}$ (d) of the 3D O-type grids
at $t=9$. The blue curves depict profiles of all 3D O-types grids, while the
purple ones plot the subsets with 2D extrema of $J_{\parallel}$.