© The Authors 2024

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 magnetohydrodynamical (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 D_B = δ_jBⁱ, where δ_j = δ/δ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. Methods for analyzing null points from satellite data have also been developed (see Fu et al. 2015; Olshevsky et al. 2020; Zhang et al. 2023).

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, $${\displaystyle {\int }_{P_1}^{P_2}{E}_{\Vert }}\text{d}s\ne 0,$$(1)

where E_‖ is the electric field strength parallel with magnetic field 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₁ and P₂, can “feel” the reconnection and will move to different field lines later.

Three-dimensional reconnection can take place at locations where strong local current concentrates and various types of 3D reconnection have been discovered such as the torsional and shearing reconnection within the fan and 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 reconnec-tion. 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, 2022; Zhong et al. 2021; Guo et al. 2023). Quasi-separatrix layers 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_‖ and J_‖ can reflect the sites of 3D reconnection according to the definition Eq. (1), which has been used in the visualization of turbulent reconnec-tion 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 B = ∇α × ∇β The three vectors ∇α, ∇β, and b̂ = B/ |B| compose a local frame system (α,β, s), where s denotes the arc coordinate along a field line. An (α,β)-pair corresponds to a field line. The gradients of the Euler potentials, ∇α and ∇β, span the local plane normal to B but, generally speaking, they are not orthogonal or normalized.

If setting the vector potential as A = α∇β satisfying the gauge A B = 0, Ohm’s law E = −v × B + N and Faraday’s law produce the evolution rule of the Euler potentials as $$\dot{\alpha }=-\frac{\partial \psi }{\partial \beta }-N^{\beta },$$(2a) $$\dot{\beta }=\frac{\partial \psi }{\partial \alpha }-N^{\alpha },$$(2b)

where N is the general nonideal term, ψ = ϕ + αδβ/δt, ϕ is the electric scalar potential, the dot operator denotes the total derivative of time d/dt, and δψ/δs = −N^s = −E_‖.

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

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 $\text{Ξ}(\alpha ,\beta )\equiv -{\displaystyle {\int }_{P_1}^{P_2}{E}_{\Vert }}(\alpha ,\beta ,s)\text{d}s$, the global differences between the evolutions of Euler potentials at P₁ and P₂ can be expressed by (Hesse et al. 2005) $${\dot{\alpha }|}_{P_2}-{\dot{\alpha }|}_{P_1}=-\frac{\partial \text{Ξ}}{\partial \beta },$$(3a) $${\dot{\beta }|}_{P_2}-{\dot{\beta }|}_{P_1}=\frac{\partial \text{Ξ}}{\partial \alpha }.$$(3b)

According to this equation, if δΞ/δα = δΞ/δβ = 0 (Ξ has an extremal value in the α-β space), P₁ and P₂ will remain on the same field line and thus cannot feel the global reconnection effect. One can notice that δΞ/δα = δΞ/δβ = 0 is equivalent to ∇_⊥Ξ = 0, where ∇_⊥ denotes the gradient perpendicular to 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 Ξ equals the reconnec-tion rate. The extremal-Ξ 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 Ξ.

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 and 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_‖ ≠ 0. In other words, we can use the distribution of E_‖ to exhibit regions where global reconnection effects can happen. Provide that only Joule dissipation is considered and the resistivity η is uniform, E_‖ can be replaced by J_‖ since E_‖ = N · b̂ = ηJ · b̂ = ηJ_‖.

Within the E_‖ ≠ 0 regions, we can further locate field lines with extremal values of Ξ. A sufficient condition for a field line to have an extremal Ξ is that all the E_‖ ≠ 0 positions threaded by this field line satisfy ∇_⊥E_‖ = 0, which can be directly proven by $$0={\displaystyle {\int }_{P_1}^{P_2}{\nabla }_{\perp }}{E}_{\Vert }\text{d}s={\nabla }_{\perp }{\displaystyle {\int }_{P_1}^{P_2}{E}_{\Vert }}\text{d}s={\nabla }_{\perp }\text{Ξ}.$$(4)

In other words, in nonideal regions, locations with ∇_⊥E_‖ = 0 approximately reveal the distributions of extremal-Ξ lines, though the trajectories of extremal-Ξ lines are not necessary to always pass through such locations. E_‖ and ∇_⊥E_‖ 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 δ/δs on both sides of Eq. (2), namely, $$\frac{\partial \dot{\alpha }}{\partial s}=\frac{\partial E_{\Vert }}{\partial \beta }-\frac{\partial N^{\beta }}{\partial s},$$(5a) $$\frac{\partial \dot{\beta }}{\partial s}=-\frac{\partial E_{\Vert }}{\partial \alpha }-\frac{\partial N^{\alpha }}{\partial s},$$(5b)

which implies that the change of $(\dot{\alpha },\dot{\beta })$ along a magnetic field line is governed by two terms, namely, the perpendicular gradient of E_‖ and the parallel gradient of N_⊥ = (N^α, N^β). It can be proved that the line conservation $(\partial \dot{\alpha }/\partial s=\partial \dot{\beta }/\partial s=0)$ is equivalent to B × (∇ × N) = 0, while the condition of flux conservation is ∇ × N = 0 (Hesse & Schindler 1988).

Furthermore, it is helpful to evaluate the importance of the E_‖ and N_⊥ terms in Eq. (5). For instance, at locations with ∇_⊥E_‖ = 0, the violation of line conservation is fully determined by the N_⊥ term. Moreover, suppose that the effects of N_⊥ term are ignorable, the local effects would be dominated by the E_‖ term, which means that analysis of E_‖ 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 (∇α, ∇β), 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 r₀ by three orthogonal unit vectors satisfying ê₃ = b̂ (r₀), ê₂ ⊥ ê₃, and ê₁ = ê₂ × ê₃. In this frame, the magnetic field at dr = r − r₀ can be approximated by $$B(\text{d}r)=B_0{\widehat{e}}_3+\text{d}r\cdot \nabla B\left(r_0\right),$$(6)

where B₀ = |B (r₀)| > 0 is the magnetic strength at r₀ and |dr| should be small enough to guarantee the validation of linear-field assumption. ê₁ and ê₂ span the same normal plane dr₃ = 0 as established by ∇α and ∇β.

The second term in Eq. (6) can be further written as two terms, $$\text{d}r\cdot \nabla B\left(r_0\right)=\left[\left(\begin{array}{c}\text{d}{r}_1\\ \text{d}{r}_2\\ 0\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \text{d}{r}_3\end{array}\right)\right]\cdot \nabla B=\left(\begin{array}{c}{B}_{\perp }^1\\ B_{\perp }^2\\ B_{\Vert }\end{array}\right)+\tilde{B},$$(7)

where $B_{\perp }^1\text{\hspace{0.17em}and\hspace{0.17em}}{B}_{\perp }^2$ are two components of the perpendicular magnetic field defined by $$B_{\perp }=M\cdot R,$$(8) $$M=\left(\begin{array}{ll}\partial B^1/\partial R_1\hfill & \partial B^1/\partial R_2\hfill \\ \partial B^2/\partial R_1\hfill & \partial B^2/\partial R_2\hfill \end{array}\right)\equiv \left(\begin{array}{ll}{M}_{11}\hfill & M_{12}\hfill \\ M_{21}\hfill & M_{22}\hfill \end{array}\right),$$(9) $$R=\left(\begin{array}{l}\text{d}{r}_1\hfill \\ \text{d}{r}_2\hfill \end{array}\right)\equiv \left(\begin{array}{l}{R}_1\hfill \\ R_2\hfill \end{array}\right),$$(10)

and B_‖ = R₁ (δB³/δR₁)+ R₂ (δB³/δR₂) is the parallel component. B_⊥ and B_‖ compose the magnetic field lying in the dr₃ = 0 plane, while $\tilde{B}\equiv \text{d}{r}_3\left(\partial B/\partial \text{d}{r}_3\right)$ represents the increments at locations of dr₃ ≠ 0. Hereafter, we analyze the local in-plane magnetic structures on the dr₃ = 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{B}$ and B_‖ is not investigated in detail. Suppose that the magnitudes of |B_‖| and $\left|\tilde{B}\right|$ are small enough compared with |B_⊥|, the magnetic structure within a finite volume near r₀ can be well approximated by B_⊥.

Because B_⊥ is not a real magnetic field satisfying the divergence-free condition, the trace of M, evaluated by tr (M) = M₁₁ + M₂₂ = −δβ³/δdr₃, can be nonzero. Meanwhile, one should notice that B_⊥ vanishes at the frame origin (R = (0,0)^T), where the superscript “T” denotes the transpose operation. In other words, the definition of the local reference frame makes an arbitrary position r₀ a 2D null point of B_⊥. Therefore, its 2D local structures can be fully determined by M, which simplifies the theoretical analysis. Moreover, for discrete fields, the structures of B_⊥ at an arbitrary grid can be directly obtained via M and the computational costs of manipulations such as interpolation and root-finding for searching null points between grids can be saved.

Inspired by the method by Parnell et al. (1996), we rewrite M as $$M=\left(\begin{array}{cc}{M}_{11}& \left(q-J_3\right)/2\\ \left(q+J_3\right)/2& M_{22}\end{array}\right),$$(11)

where q = M₁₂ + M₂₁ and J₃ = J_‖ = M₂₁ − M₁₂. Defining a threshold current $J_{\text{thres}}=\sqrt{{\left(M_{11}-M_{22}\right)}^2+q^2}$, the discriminant of the eigenvalue quadratic |λl − M| = 0, becomes D (M) = 4Det (M) − tr ${(M)}^2=J_3^2-J_{\text{thres}}^2$, where Det (M) is the determinant of M. The eigenvalues are thus $\lambda =(\mathrm{tr}(M)\pm \sqrt{-D(M)})/2$. If |J₃| < J_thres (i.e. D (M) < 0), the eigenvalues of M are distinct real numbers; if |J₃| > J_thres (i.e. D (M) > 0), the eigenvalues of M are two conjugate complex numbers; if |J₃| = J_thres (i.e. D (M) = 0), two eigenvalues are repeated real roots.

The local structures of B_⊥ can be classified into nine types based on the values of tr (M), D (M), and λ (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_‖ ≠ 0 (see Fig. B.1a). Types 2 and 3 are 3D O-types but possess radial components from the trace of M. For different signs of tr (M), field lines of types 2 and 3 show different directions toward r₀ (see Figs. B.1b,c). Types 4 and 5 are respectively source and sink structures dominated by the trace of M (see Figs. B.1d, e). Type 6 has one zero eigenvalue, implying the existence of a neutral line (see Fig. B.1f). 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).

To investigated the geometric properties of B_⊥ , we further decompose M into two parts as $$M=M^{\prime }+T,$$(12)

where, $$M^{\prime }=\left(\begin{array}{cc}p& \left(q-J_3\right)/2\\ \left(q+J_3\right)/2& -p\end{array}\right),$$(13) $$T=\left(\begin{array}{cc}\mathrm{tr}(M)/2& 0\\ 0& \mathrm{tr}(M)/2\end{array}\right),$$(14) $$p=\left(M_{11}-M_{22}\right)/2.$$(15)

M' is a traceless matrix corresponding to the source-free part of B_⊥ which can be treated as a real magnetic field, while T is the source part originating from the 3D effects. Correspondingly, B_⊥ can be written as $B_{\perp }=B_{\perp }^{\prime }+{\nabla }_{p}G$, where $$B_{\perp }^{\prime }=\mathcal{J}{\nabla }_{p}F,$$(16) $$F=\frac{1}{4}\left[\left(q-J_3\right)R_2^2-\left(q+J_3\right)R_1^2\right]+p{R}_1{R}_2,$$(17) $$G=\frac{1}{4}\mathrm{tr}(M)\left(R_1^2+R_2^2\right),$$(18) $${\nabla }_p={\left(\partial /\partial R_1,\partial /\partial R_2\right)}^T,$$(19)

and $$\mathcal{J}=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$$(20)

is the symplectic matrix. Because the second term ∇_pG is isotropic, the anisotropic properties of B_⊥ are completely determined by the first term $B_{\perp }^{\prime }$. As a 2D magnetic field, the integration curves of $B_{\perp }^{\prime }$ are contour lines of its Hamiltonian F.

Following the method by Parnell et al. (1996), we define a new coordinate X = (X₁, X₂)^T satisfying $$X=\left(\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right)R,$$(21)

where the rotation angle θ satisfies tan 2θ = −2p/q. Then F can be transformed to $$F=\frac{1}{4}\left[X_2^2\left(J_{\text{thres }}-J_3\right)-X_1^2\left(J_{\text{thres }}+J_3\right)\right]\text{. }$$(22)

Therefore, the field lines of B' are hyperbolic lines if |J₃| < J_thres (Fig. 1a), are concentric elliptical curves if |J₃| > J_thres (Fig. 1b), and are parallel straight lines if |J₃| = J_thres.

To quantitatively evaluate the geometric shapes of $B_{\perp }^{\prime }$ lines, we define an angle, $${\theta }_{\text{eig }}=\{\begin{array}{ll}\arctan \left(\sqrt{{\left(J_{\text{thres }}/J_3\right)}^2-1}\right),\hfill & \left|J_3\right|<J_{\text{thres }},\hfill \\ \arctan \left(\sqrt{{\left(J_3/J_{\text{thres }}\right)}^2-1}\right),\hfill & \left|J_3\right|>J_{\text{thres }}.\hfill \end{array}$$(23)

For |J₃| < J_thres, θ_eig is the acute angle spanned by the separatrix lines (see Fig. 1a). For |J₃| > J_thres, θ_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, θ_eig ∈ [0, 90°]. If J₃→ J_thres, then θ_eig → 0, which corresponds to the anti-parallel case. Meanwhile, θ_eig → 90° means that the X-type $B_{\perp }^{\prime }$ tends to be a potential field satisfying J₃ → 0 while the O-type $B_{\perp }^{\prime }$ represents as concentric circle lines governed only by the antisymmetric part of M' (the parallel current J₃);

Because θ_eig only reflects the property of the source-free part M', we define another parameter to evaluate the importance of T, namely the trace ratio $$R_{\text{tr}}=\Vert T\Vert /\Vert M^{\prime }\Vert \text{, }$$(24)

where || · || is the Euclidean norm. Larger R_tr means larger effects from the trace part.

Fig. 1 Cartoons of $B_{\perp }^{\prime }$ lines of X (a) and O (b) types. The blue curves depict field lines. The definitions of θeig are labeled by orange markers.

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 |E_‖|

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

Fig. 2 Schematic diagram of the local frame at a grid rlmn. 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 (ΔL) away from the origin point. The green dots are six adjacent grids surrounding rlmn in the frame of original discrete data.

3.2 Step 2: Analyze local magnetic structures

First, approximate the 3D magnetic Jacobian matrix D_B on each grid via two-order central difference as $$D_B^{lmn}=\left(\frac{B^{l+1,m,n}-B^{l-1,m,n}}{2\text{Δ}{x}^1},\frac{B^{l,m+1,n}-B^{l,m-1,n}}{2\text{Δ}{x}^2},\frac{B^{l,m,n+1}-B^{l,m,n-1}}{2\text{Δ}{x}^3}\right),$$(25)

where l, m, and n are indices of x¹, x², and x³, respectively. B^lmn= (B^1,lmn,B^2,lmn,B^3,lmn)^T is the magnetic field at 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.

Second, construct local frame at r^lmn (see Fig. 2). The base vector on the third direction (the direction of the local magnetic field) is ${\widehat{e}}_3^{lmn}=\left(B^{1,lmn},B^{2,lmn},B^{3,lmn}\right)/\left|B^{lmn}\right|\cdot {\widehat{e}}_2^{lmn}$ can be any unit vector perpendicular to ${\widehat{e}}_3^{lmn}$. As a convenient choice, we set it as ${\widehat{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 ${\widehat{e}}_2^{lmn}=(0,1,0)$. The last base vector can be directly obtained by ${\widehat{e}}_1^{lmn}={\widehat{e}}_2^{lmn}\times {\widehat{e}}_3^{lmn}$. Hereafter, the local plane spanned by ${\widehat{e}}_1^{lmn}\text{and}{\widehat{e}}_2^{lmn}$ within a cell is called the magnetic projection plane (MPP).

Third, transform $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}{l}{\widehat{e}}_1^{lmn}\hfill \\ {\widehat{e}}_2^{lmn}\hfill \\ {\widehat{e}}_3^{lmn}\hfill \end{array}\right),$$(26)

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

Four, analyze the local magnetic structures via calculating the trace, discriminant, eigenvalues, θ_eig, and R_tr of the matrix $$M^{lmn}=\left(\begin{array}{ll}{D}_{B,11}^{\prime lmn}\hfill & D_{B,12}^{\prime lmn}\hfill \\ D_{B,21}^{\prime lmn}\hfill & D_{B,22}^{\prime lmn}\hfill \end{array}\right).$$(27)

Their definitions are exhibited in Sect. 2.3.

3.3 Step 3: Locate reconnection sites of extremal E_‖

As discussed in Sect. 2.2, we further locate the positions with extremal values of E_‖. Theoretically, they can be identified by ∇_⊥E_‖ = 0. However, numerically, ∇_⊥E_‖ 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 r^lmn is a 2D extremal point of E_‖ on the MPP. To implement this, four grids adjacent with r^lmn (see Fig. 2) on the MPP is defined as, $$R^{-0}=r^{lmn}-\text{Δ}L{\widehat{e}}_1^{lmn},$$(28)

$$R^{+0}=r^{lmn}+\text{Δ}L{\widehat{e}}_1^{lmn},$$(29)

$$R^{0-}=r^{lmn}-\text{Δ}L{\widehat{e}}_2^{lmn},$$(30)

$$R^{0+}=r^{lmn}+\text{Δ}L{\widehat{e}}_2^{lmn},$$(31)

where ΔL is the spatial step length of the original discrete data. The values of E_‖ at these positions can be obtained via linear interpolation of original data. If the values of E_‖ have the same signs at r, R⁻⁰, R⁺⁰, R⁰⁻, and R⁰⁺, and |E_‖| at r^lmn is maximum among the five positions, then we identify r^lmn as a 2D extreme point of |E_‖|.

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)¹. 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. It should be noticed that they follow the “meshgrid” data model of Matlab, namely, the first dimension is x², the second dimension is x¹, and the third dimension is x³. 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 N to ARD. Without N1–N3, ARD will suppose the nonideal term is induced by a Joule dissipation with a constant resistivity and use J = ∇ × B to replace 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 |E_‖| is defined by Parameters.ARD_ScalarThreshold. To help users determine this value, ARD can present a histogram of ^||E_‖| without proceeding further analysis if Parameters.ARD_ShowThresScalarProfile is set as 1. Because the magnetic Jacobian matrix $D_B^{lmn}$ might have small trace components resulting from numerical errors, ARD provides an option to clean the trace by $D_B^{lmn}-1/3\text{tr}\left(D_B^{lmn}\right)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_RD × 7 matrix, where N_RD is the number of recognized grids. 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¹, x², and x³ coordinates, respectively. RDType is an integer number taking values from 1 to 9, marking the types of B_⊥as listed in Table 1. Is2DExtrema labels whether the grid is a 2D extreme point of E_‖. EigAngle and RatioMTrace save the values of θ_eig and R_tr, respectively. RDInfo.ExtraData is an optional output controlled by Parameters.OutputExtraData, which has 25 columns containing information of |B^lmn|, $D_B^{\prime lmn},{\widehat{e}}_1^{lmn},{\widehat{e}}_2^{lmn},{\widehat{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.

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₁ × n₂ × n₃ data B¹ 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¹ in a 1D array, in which the x¹ -direction index changes fastest and the x³-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¹-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 input and output 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.

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: $$\begin{array}{cc}\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho u)& =0,\\ \frac{\partial (\rho u)}{\partial t}+\nabla \cdot \left(\rho uu-BB+P^{\ast }I\right)& =0,\\ \frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+P^{\ast }\right)u-B(B\cdot u)\right]& =0,\\ \frac{\partial B}{\partial t}-\nabla \times (u\times B)& =-\nabla \times (\eta J),\\ J& =\nabla \times B,\end{array}$$(32)

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

The initial magnetic field forms a force-free Harris sheet, namely, $$\begin{array}{l}{B}_x=B_0\tanh (y/\lambda ),\\ B_y=0,\\ B_z=B_0/\cosh (y/\lambda ),\end{array}$$(33)

where B₀ = 1 and λ = 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 β = 0.1. We set a uniform background resistivity η = 5 × 10⁻⁶ to obtain a Lundquist number S = L₀u₀/η = 2 × 10⁵. 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),$$(34)

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 σ_u = 10⁻³, which works as a seed of the self-sustained turbulence (also see Huang & Bhattacharjee 2016).

The simulation domain is x ∈ [−2.5, 2.5], y ∈ [−0.5, 0.5], and z ∈ [−0.2, 0.2]. 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 < |y| < 0.3. The level-2 grid, doubling again, is set in the core reconnection region |y| ≤ 0.1, implementing an effective mesh of 4800 × 960 × 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 Ã_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 Figs. 3a and b).

5.3 Qualitative global picture

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

Based on the algorithm in Sect. 3, we obtain the information of magnetic structures at strong J_‖ regions. In Fig. 3, to give a clear picture, we only plot 3D X and O-type grids having 2D projected extrema of J_‖ (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 and 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 (Figs. 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 |J_‖| 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 θ_eig is close to zero. This fact explains the reason why we use θ_eig as a key parameter of B_⊥. 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 |J_‖| and P_oh = ηJ² are small corresponding to a moderate reconnection; θ_eig mainly distributes below 10° meaning the reconnection sites tend to have magnetic structures that are close to anti-parallel; R_tr is close to 0 implying the 3D effect is weak. With the development of turbulence, the number density curves of |J_‖|, P_oh, θ_eig, and R_tr keep extending toward larger values (see Fig. 6). At t =9 , the number density of |J_‖| first decreases monotonically as |J_‖| increases to 300, then forms a platform at |J_‖| ∈ [300, 400], and finally rapidly damps to zero at |J_‖| ~ 440 (Fig. 6a). The distribution of P_oh presents a similar trend with |J_‖|(Fig. 6b). θeig mainly distributes in the region θ_eig < 40^° and its peak number density appears near 5^° (Fig. 6c). The maximum value of θ_eig reaches 70^° and there are still considerable grids as θ_eig → 0. The number density of R_tr decreases monotonically in the entire domain. After 3D effects get significant, more grids with larger R_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 |J_‖| and P_oh might be sensitive to the magnitude of η.

We extract the X-type grids with 2D extrema of J_‖ at t = 9 to compare its distribution with that of all X-type grids (Fig. 7). The number densities of |J_‖| and P_oh of 2D extreme grids almost coincide with that of all X-grids for large values (see the |J_‖| > 300 region in Fig. 7a and the P_oh > 0.6 region in Fig. 7b). But for smaller values, they are lower than the curves of all X-grids (see Figs. 7a and b). For θ_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_tr have similar trends (Fig. 7d), but 2D extreme grids tend to concentrate within a smaller value region (R_tr < 0.1), implying that, near the extremal-Ξ lines, the variations of local field strength on the direction of guide field (ê₃) is weaker.

The statistical results of 3D O-type grids are also provided in Appendix C. The rules of |J_‖|, P_oh, and R_tr of 3D O-type grids are similar to 3D X-type grids (see Figs. C.1 and C.2). The main difference is that more O-type grids with large θ_eig near 90° are identified, corresponding to flux ropes with poloidal field lines close to circle shapes (see Figs. C.1c and C.2c).

Fig. 3 Distributions of mass density p a, temperature T b, parallel current density |J‖| c, and 3D X and O-type grids d at t = 9. For p and T, we only depict their distributions on the current-sheet middle plane (y = 0), the front and back planes (z = ±0.2), and the left and right planes (x = ±1). Panel c shows the 3D profile of |J‖| satisfying |J‖| > 80. In panel d, we plot the grids that (1) satisfy |J‖| > 50, (2) have local magnetic structures of 3D X and O types, and (3) are 2D projected extreme points of J‖. 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–9 is available.

Fig. 4 Four snapshots of 3D X and O-type grids exhibiting the development of turbulent reconnection. The gray dashed lines approximately trace the origin of turbulent reconnection regions.

Fig. 5 Evolutions of grid numbers for different types of structures. Only grids with |J‖| > 50 are counted.

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).

Fig. 6 Number densities of |J‖| (a), Poh (b), θeíg (c), and Rtr (d) of the 3D X-type grids at different moments. Different colors correspond to different moments.

Fig. 7 Number densities of |J‖| (a), Poh (b), θeig (c), and Rtr (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‖.

Fig. 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‖, satisfying |J‖| > 11.9 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 Bz is overlapped at the bottom.

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_‖ in representing the reconnection sites and propose that the locations of ∇_⊥E_‖ = 0 can present the distribution of extremal-Ξ 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, θ_eig and R_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 such as 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. B.1). 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 recon-nection sites, our method can also determine the classic 2D X and 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_‖ and further selecting grids with larger¹ |J_‖| and θ_eig, we find that the resultant distribution of 3D X-grids gets almost consistent with the 2D X-lines (see Figs. 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 E and electric drift speed. Because the Lorentz boost speed to eliminate the magnetic field perpendicular to E, namely v_L = c²E × B/E², 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 υ_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 υ_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 recon-nection 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_‖ concentrate (see Pontin et al. 2007), our method is still useful near null points.

Second, when analyzing the local magnetic structures near a grid r^lmn, we approximate the surface normal to magnetic field by a plane (MPP) within a finite volume with the length scale of δL, which validates better if $R_{\nabla \text{B}}\equiv \text{Δ}L\Vert D_B^{\prime lmn}\Vert /\left|B^{lmn}\right|\ll 1$. Meanwhile, we mainly use the perpendicular magnetic field on the MPP (B_⊥) to determine the local magnetic structures which might also be affected by B_‖ and $\tilde{\text{B}}$. In applications, we define three ratios to evaluate the importance of the three terms, namely, $R_{B_{\perp }}\equiv \text{Δ}L\Vert M^{lmn}\Vert /\left|B^{lmn}\right|,R_{B_{\Vert }}\equiv \text{Δ}L\sqrt{{\left(D_{B,31}^{\prime lm}\right)}^2+{\left(D_{B,32}^{\prime mn}\right)}^2}/\left|B^{lmn}\right|$, and $R_{\tilde{\text{B}}}\equiv \text{Δ}L\sqrt{{\left(D_{B,13}^{\prime lmn}\right)}^2+{\left(D_{B,23}^{\prime lmn}\right)}^2+{\left(D_{B,33}^{\prime lmn}\right)}^2}/\left|B^{lmn}\right|$. In the Harris-sheet simulation, for most grids, the magnitude of R_∇B is ~ 0. 1 which implies the MPP approximation is acceptable (see the black dashed curves in Fig. 10). Moreover, the distribution of $R_{{\text{B}}_{\perp }}$ almost overlaps with that of R_∇B (see the blue curves in Fig. 10), the peak of $R_{B_{\Vert }}$ is about one order smaller than $R_{{\text{B}}_{\perp }}$, and the magnitude of $R_{\tilde{\text{B}}}$ is much smaller even compared with $R_{B_{\Vert }}$. Therefore, in this simulation, the variation of magnetic field in a cell is dominated by B_⊥ for most reconnection grids. R_∇B, $R_{{\text{B}}_{\perp }},R_{B_{\Vert }},\text{and\hspace{0.17em}}{R}_{\widehat{\text{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 recon-nection 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 including 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.

Fig. 9 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 determine d 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 |J‖| is marked by a magenta “x” marker. On the plane z = 0.2, the profile of ux 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 |J‖| and θeig. Panel d has the same thresholds as c but also plots 3D O-type grids with orange dots.

Fig. 10 Number densities of $R_{{\text{B}}_{\perp }}$ (blue curves), $R_{B_{\Vert }}$ (purple curves), $R_{\tilde{\text{B}}}$ (green curves), and R∇B (black dashed curves) of all grids with |J‖| > 50 at different moments in the Harris-sheet simulation.

Movie

Movie associated to Fig.3 Access here

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).

Appendix A Analysis of local reconnection effects

To quantitatively compare the E_‖ and N_⊥ 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 ê_1; ê₂, and ê₃. In this section, we use x = (x¹, x², x³) and x' = (x'¹, x'², x'³) = (α,β, s)^T to denote the coordinates in the (ê₁, ê₂, ê₃) and (∇α, ∇β, b̂) frames, respectively. One can notice that ê₃ = b̂, x³ = x'³ = s, and δ/δs = δ/δx³.

First, based on the vector relations $$x^1{\widehat{e}}_1+x^2{\widehat{e}}_2=\alpha \nabla \alpha +\beta \nabla \beta ,$$(A.1)

$$\nabla \alpha =\left({\widehat{e}}_1\cdot \nabla \alpha \right){\widehat{e}}_1+\left({\widehat{e}}_2\cdot \nabla \alpha \right){\widehat{e}}_2,$$(A.2)

$$\nabla \beta =\left({\widehat{e}}_1\cdot \nabla \beta \right){\widehat{e}}_1+\left({\widehat{e}}_2\cdot \nabla \beta \right){\widehat{e}}_2,$$(A.3)

the transformation matrices between the two frames can be obtained as $$\frac{\partial x}{\partial x^{\prime }}=\left(\begin{array}{ccc}{\widehat{e}}_1\cdot \nabla \alpha ,& {\widehat{e}}_1\cdot \nabla \beta ,& 0\\ {\widehat{e}}_2\cdot \nabla \alpha ,& {\widehat{e}}_2\cdot \nabla \beta ,& 0\\ 0,& 0,& 1\end{array}\right),$$(A.4)

$$\frac{\partial x^{\prime }}{\partial x}=\frac{1}{B}\left(\begin{array}{ccc}{\widehat{e}}_2\cdot \nabla \beta ,& -{\widehat{e}}_1\cdot \nabla \beta ,& 0\\ -{\widehat{e}}_2\cdot \nabla \alpha ,& {\widehat{e}}_1\cdot \nabla \alpha ,& 0\\ 0,& 0,& B\end{array}\right),$$(A.5)

where B = b̂ · (∇α × ∇β) is the magnetic strength at the origin point.

Second, we rewrite the δE_‖/δβ and δE_‖/δα terms in Eq. 5. Considering that E · b̂ = E · ê₃, E_‖ is invariant in the new frame. According to Eq. A.4, we have δ/δα = (ê₁ · ∇α) δ/δx¹ + (ê₂ · ∇α) δ/δx² and δ/δβ = (ê₁ · ∇β) δ/δx¹ + (ê₂ · ∇β) δ/δx², and thus $$\frac{\partial E_{\Vert }}{\partial \alpha }=\frac{\partial E_{\Vert }}{\partial x^1}{\widehat{e}}_1\cdot \nabla \alpha +\frac{\partial E_{\Vert }}{\partial x^2}{\widehat{e}}_2\cdot \nabla \alpha =\nabla E_{\Vert }\cdot \nabla \alpha ,$$(A.6)

$$\frac{\partial E_{\Vert }}{\partial \beta }=\frac{\partial E_{\Vert }}{\partial x^1}{\widehat{e}}_1\cdot \nabla \beta +\frac{\partial E_{\Vert }}{\partial x^2}{\widehat{e}}_2\cdot \nabla \beta =\nabla E_{\Vert }\cdot \nabla \beta .$$(A.7)

Third, we rewrite the δN^β/δs and δN^α/δs terms in Eq. 5. According to Eq. A.5, we have $$N^{\alpha }=\frac{1}{B}\left(N^1{\widehat{e}}_2\cdot \nabla \beta -N^2{\widehat{e}}_1\cdot \nabla \beta \right),$$(A.8)

$$N^{\beta }=-\frac{1}{B}\left(N^1{\widehat{e}}_2\cdot \nabla \alpha -N^2{\widehat{e}}_1\cdot \nabla \alpha \right).$$(A.9)

Therefore, $$\begin{array}{cc}\frac{\partial N^{\alpha }}{\partial s}=& \left[{\widehat{e}}_2\frac{\partial }{\partial x^3}\left(\frac{N^1}{B}\right)-{\widehat{e}}_1\frac{\partial }{\partial x^3}\left(\frac{N^2}{B}\right)\right]\cdot \nabla \beta \\ & +\left(\frac{N^1}{B}\frac{\partial {\widehat{e}}_2}{\partial x^3}-\frac{N^2}{B}\frac{\partial {\widehat{e}}_1}{\partial x^3}\right)\cdot \nabla \beta \\ & +\frac{1}{B}\left(N^1{\widehat{e}}_2-N^2{\widehat{e}}_1\right)\cdot \frac{\partial \nabla \beta }{\partial x^3}.\end{array}$$(A.10)

Using ∇ (b̂ · ∇β = 0, it can be proved that δ∇β/δx³ = −∇b̂ · ∇β = −∇ê₃ · ∇β and similarly δ∇α/δx³ = −∇ê₃ · ∇α. As a result, we can define a vector $$\Gamma ={\widehat{e}}_2\frac{\partial \left(N^1/B\right)}{\partial x^3}-{\widehat{e}}_1\frac{\partial \left(N^2/B\right)}{\partial x^3}+\frac{N^1}{B}\frac{\partial {\widehat{e}}_2}{\partial x^3}$$

$$-\frac{N^2}{B}\frac{\partial {\widehat{e}}_1}{\partial x^3}+\left(-\frac{N^1}{B}{\widehat{e}}_2+\frac{N^2}{B}{\widehat{e}}_1\right)\cdot \nabla {\widehat{e}}_3,$$(A.11)

and write Eq. A.10 as δN^α/δs = Γ · ∇β. Similarly, we have δN^β/δs = −Γ ·∇α.

Finally, Eq. 5 becomes $$\frac{\partial \dot{\alpha }}{\partial x^3}=\nabla E_{\Vert }\cdot \nabla \beta +\Gamma \cdot \nabla \alpha$$(A.12a)

$$\frac{\partial \dot{\beta }}{\partial x^3}=-\nabla E_{\Vert }\cdot \nabla \alpha -\Gamma \cdot \nabla \beta .$$(A.12b)

Because both ∇α and ∇β are perpendicular to ê₃, only the perpendicular components of ∇E_‖ and Γ, namely, ∇_⊥ E_‖ and Γ_⊥ , have effects on the local line conservation. Considering that ∇_⊥ E_‖ and Γ_⊥ can be explicitly calculated in the (ê₁, ê₂, ê₃) frame, we can use them to qualitatively compare ∇E_‖ and N_⊥ terms even withoutknowing ∇α and ∇β. However, because an absolute com-parison still relies on the flux coordinate frame, ∇_⊥ E_‖ and Γ_⊥ only provide a statistical perspective of approximate comparison.

For instance, we compare ∇_⊥ E_‖ and Γ_⊥ in the Harris-sheet simulation shown in Sect. 5 (see Fig. A.1). At t = 9, the peak of |∇_⊥ E_‖| locates near 0.1 while |Γ_⊥| mainly distributes near 0.01 (see the blue curves in Fig. A.1a). Moreover, during the entire evolution, both the mean values and standard deviations of |∇_⊥ E_‖| are significantly larger than that of Γ_⊥ (see Fig. A.1c). These results show that, statistically speaking, the magnitude of ∇_⊥ E_‖ is larger than Γ_⊥ . On the other hand, we depict the number density of the angle θ_le spanned by ∇_⊥ E_‖ and Γ_⊥, which shows that θ_le is most likely to be 90°, namely, the two vectors tend to normal to each other for most of grids (see Fig. A.1b). During the evolution, the mean value of θ_le keeps close to 90° and its standard deviation is smaller than that of a uniform distribution, which implies that the distributions of θ_le tend to concentrate near 90^°.

It should also be noticed that the grids with extremal E_‖ identified by the algorithm introduced in Sect. 3 compose a subset about 12% of all grids and have smaller values of |∇_⊥ E_‖| (see the purple solid curve in Fig. A.1a). Therefore, the algorithm of locating 2D extremal E is a good numerical approximation of ∇_⊥ E_‖=0.^.

When using the ARD code, users can enable the calculation of ∇_⊥ E_‖ and Γ_⊥ by setting ARD_AnalyzeLocalEffects as 1, which costs more RAM and computation resources. All the derivatives in the expressions of ∇_⊥ E_‖ and Γ_⊥ 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 𝒯^lmn in Eq. 26. ARD also outputs |∇ × N| and |B × (∇ × N)| to study the local flux and line conservations (see Table 3).

Fig. A.1 Statistical comparisons of ∇⊥ E‖ and Γ⊥ in the Harris-sheet simulation. Panel (a) plots the number densities of |∇⊥ E‖| (solid curves) and |Γ⊥| (dashed curves) at t = 9. Panel (b) shows the number densities of the angle θle spanned by Γ⊥ and ∇⊥ E‖ at t = 9. In panels (a) and (b), the blue curves depict all grids satisfying |J‖| ≥ 50 while the purple ones plot the subsets with 2D extremal E‖. Panel (c) depicts the time evolutions of the mean values of |∇⊥ E‖| (the blue solid curve) and |Γ⊥| (the orange solid curve). The blue and orange dashed curves plot the standard deviations of |∇⊥ E‖| and |Γ⊥|, respectively. Panel (d) exhibits the time evolutions of the mean values of θ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° to 180°. The mean values and standard deviations in panels (c) and (d) are calculated from all grids satisfying |J‖| ≥ 50.

Appendix B Examples of magnetic structures on the MPP

In Fig. B.1, we exhibit the B_⊥ 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.

Fig. B.1 Six typical local field-line structures of B⊥ (types 1 to 6) as found in the reconnection region at t = 9. Panels (a)–(f) depict the field lines of B⊥ 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⊥ = |B⊥| are also plotted with a gray color map. In panel (g), to give a clear picture, we only plot the grids with θeig larger than a threshold value. For types 1–3, the threshold is set as 30°; for types 4–5, the threshold is 5°.

Appendix C Quantitative statistical results of O-type grids in the Harris-sheet simulation

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

Fig. C.1 Number densities of |J‖| (a) , Poh (b), θeig (c), and Rtr (d) of the 3D O-type grids at different moments. Different colors correspond to different moments.

Fig. C.2 Number densities of |J‖| (a), Poh (b), θeíg (c), and Rtr (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‖.

References

All Tables

All Figures

	Fig. 1 Cartoons of $B_{\perp }^{\prime }$ lines of X (a) and O (b) types. The blue curves depict field lines. The definitions of θeig are labeled by orange markers.
In the text

	Fig. 2 Schematic diagram of the local frame at a grid rlmn. 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 (ΔL) away from the origin point. The green dots are six adjacent grids surrounding rlmn in the frame of original discrete data.
In the text

Fig. 3 Distributions of mass density p a, temperature T b, parallel current density |J‖| c, and 3D X and O-type grids d at t = 9. For p and T, we only depict their distributions on the current-sheet middle plane (y = 0), the front and back planes (z = ±0.2), and the left and right planes (x = ±1). Panel c shows the 3D profile of |J‖| satisfying |J‖| > 80. In panel d, we plot the grids that (1) satisfy |J‖| > 50, (2) have local magnetic structures of 3D X and O types, and (3) are 2D projected extreme points of J‖. 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–9 is available.

In the text

	Fig. 4 Four snapshots of 3D X and O-type grids exhibiting the development of turbulent reconnection. The gray dashed lines approximately trace the origin of turbulent reconnection regions.
In the text

	Fig. 5 Evolutions of grid numbers for different types of structures. Only grids with |J‖| > 50 are counted.
In the text

	Fig. 6 Number densities of |J‖| (a), Poh (b), θeíg (c), and Rtr (d) of the 3D X-type grids at different moments. Different colors correspond to different moments.
In the text

	Fig. 7 Number densities of |J‖| (a), Poh (b), θeig (c), and Rtr (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‖.
In the text

Fig. 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‖, satisfying |J‖| > 11.9 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 Bz is overlapped at the bottom.

In the text

Fig. 9 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 determine d 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 |J‖| is marked by a magenta “x” marker. On the plane z = 0.2, the profile of ux 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 |J‖| and θeig. Panel d has the same thresholds as c but also plots 3D O-type grids with orange dots.

In the text

	Fig. 10 Number densities of $R_{{\text{B}}_{\perp }}$ (blue curves), $R_{B_{\Vert }}$ (purple curves), $R_{\tilde{\text{B}}}$ (green curves), and R∇B (black dashed curves) of all grids with |J‖| > 50 at different moments in the Harris-sheet simulation.
In the text

Fig. A.1 Statistical comparisons of ∇⊥ E‖ and Γ⊥ in the Harris-sheet simulation. Panel (a) plots the number densities of |∇⊥ E‖| (solid curves) and |Γ⊥| (dashed curves) at t = 9. Panel (b) shows the number densities of the angle θle spanned by Γ⊥ and ∇⊥ E‖ at t = 9. In panels (a) and (b), the blue curves depict all grids satisfying |J‖| ≥ 50 while the purple ones plot the subsets with 2D extremal E‖. Panel (c) depicts the time evolutions of the mean values of |∇⊥ E‖| (the blue solid curve) and |Γ⊥| (the orange solid curve). The blue and orange dashed curves plot the standard deviations of |∇⊥ E‖| and |Γ⊥|, respectively. Panel (d) exhibits the time evolutions of the mean values of θ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° to 180°. The mean values and standard deviations in panels (c) and (d) are calculated from all grids satisfying |J‖| ≥ 50.

In the text

Fig. B.1 Six typical local field-line structures of B⊥ (types 1 to 6) as found in the reconnection region at t = 9. Panels (a)–(f) depict the field lines of B⊥ 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⊥ = |B⊥| are also plotted with a gray color map. In panel (g), to give a clear picture, we only plot the grids with θeig larger than a threshold value. For types 1–3, the threshold is set as 30°; for types 4–5, the threshold is 5°.

In the text

	Fig. C.1 Number densities of |J‖| (a) , Poh (b), θeig (c), and Rtr (d) of the 3D O-type grids at different moments. Different colors correspond to different moments.
In the text

	Fig. C.2 Number densities of |J‖| (a), Poh (b), θeíg (c), and Rtr (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‖.
In the text