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

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.


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 selfsustained 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 recon-nection, how to efficiently locate the reconnection sites and ana-25 lyze their properties from the massive discrete data has become 26 a new challenge (Vlahos & Isliker 2023).

27
A typical location of 3D reconnection is at a magnetic null 28 point, where the magnetic strength vanishes.In the adjacent re-29 gion of a null point, magnetic field forms a fan-spine structure, 30 the local structures of which can be evaluated by the magnetic 31 Jacobian matrix D B = ∂ j B i , where ∂ j ≡ ∂/∂x j and i, j = 1, 2, 3 32 (see Lau & Finn 1990;Parnell et al. 1996).Haynes & Parnell 33 (2007) developed a trilinear method for locating 3D null points 34 in discrete numerical data, which has been widely applied and 35 can also be directly used in turbulent cases.

36
However, 3D reconnection can happen without the presence 37 of null points, thus being more difficult to be analyzed.Schindler 38 et al. (1988) and Hesse & Schindler (1988) proposed the gen-39 eral magnetic reconnection theory and obtained the condition 40 for the "global" effect of reconnection with finite magnetic field, 41 namely, region, P 1 and P 2 , 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 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. 1996Demoulin et al. , 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).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 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) de-veloped the magnetic flux transport method which was recently 110 applied to a kinetic plasma turbulence simulation by Li et al. 111 (2023).

112
For general 3D cases without a strong guide field, local refer-113 ence frames are necessary to analyze small-scale structures.Par-114 nell et al. (2010a) set the direction of a pre-determined separator 115 as the normal direction of the projection plane, on which the 2D 116 projected magnetic field is classified.Kowal et al. (2020) defined 117 a local frame based on the shear tensor of magnetic field, which 118 is further used to determine the properties of current sheets.Re-119 cently, Lapenta (2021) developed a reconnection-site identifica-120 tion method via the electric drift speed, which performs well for 121 data of kinetic simulations.

122
Methods independent with local frames have also been de-123 veloped for the identification of turbulent reconnection struc-124 tures (see Vlahos & Isliker 2023), for instance, the Phase Co-125 herence Index method (Hada et al. 2003), the Partial Variance of 126 Increments (PVI) method (Greco et al. 2017), and the fractal di-127 mension method by box-counting (Isliker et al. 2019).Although 128 these methods can produce useful statistical results, they cannot 129 analyze the magnetic structures at reconnection sites.

130
In this paper, we introduce an efficient method for locating 131 small-scale reconnection structures in 3D discrete data involving 132 turbulent eddies.Both the theoretical basis and the numerical al-133 gorithm of this method are presented in detail.Through two typ-134 ical benchmarks, it is well validated and presents a promise in 135 qualitative and quantitative analysis of 3D reconnection simula-136 tions containing multi-scale structures.

137
The following parts of this paper are organized as follows.138 In Sect.2, we exhibit the theoretical basis of our method.The 139 numerical strategy is presented in Sect.3. We provide a Matlab 140 implementation of this method and introduce its usage in Sect. 4. 141 This method is benchmarked by two typical simulations which 142 are detailed in Sects.5 and 6.The results of this paper are sum-143 marized and discussed in Sect.7.
where N is the general nonideal term, ψ = φ + α∂β/∂t, φ is 159 the electric scalar potential, the dot operator denotes the total 160 derivative of time d/dt, and ∂ψ/∂s = −N s = −E .
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 1 and P 2 , which directly results in the condition Eq. 1.Because N = 0 at P 1 and P 2 , 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 Ξ (α, β) ≡ − E (α, β, s) ds, the global differences between the evolutions of Euler potentials at P 1 and P 2 can be expressed by (Hesse et al. 2005) α| β According to this equation, if ∂Ξ/∂α = ∂Ξ/∂β = 0 (Ξ has an extremal value in the α-β space), P 1 and P 2 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 reconnection 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 Ξ.

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 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 In other words, in nonideal regions, locations with ∇ ⊥ E = 0 approximately reveal the distributions of extremal-Ξ lines, though which implies that the change of α, β along a magnetic field 230 line is governed by two terms, namely, the perpendicular gradi-231 ent of E and the parallel gradient of Hesse & Schindler 1988).

235
Furthermore, it is helpful to evaluate the importance of 236 the E and N ⊥ terms in Eq. 5.For instance, at locations with 237 ∇ ⊥ E = 0, the violation of line conservation is fully determined 238 by the N ⊥ term.Moreover, suppose that the effects of N ⊥ term 239 are ignorable, the local effects would be dominated by the E 240 term, which means that analysis of E can reveal the properties 241 of not only global but also local effects of reconnection.How-242 ever, because a quantitative comparison of the two terms relies 243 on a local flux coordinate frame (∇α, ∇β), which is difficult 244 to establish for arbitrary magnetic fields.In Appendix A, we 245 attempt to discuss this problem and propose a statistical method 246 to compare the effects of the two terms.
where B 0 = |B (r 0 )| > 0 is the magnetic strength at r 0 and |dr| 256 should be small enough to guarantee the validation of linear-field 257 assumption.ê1 and ê2 span the same normal plane dr 3 = 0 as 258 established by ∇α and ∇β.

259
The second term in Eq. 6 can be further written as two terms, 260 where B 1 ⊥ and B 2 ⊥ are two components of the perpendicular magnetic field defined by Inspired by the method by Parnell et al. (1996), we rewrite where q = M 12 + M 21 and ).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 B ⊥ .The star marker " * " denotes complex conjugate.R and C are the sets of real and complex numbers, respectively. Type To investigated the geometric properties of B ⊥ , we further 312 decompose M into two parts as where, M is a traceless matrix corresponding to the source-free part of 315 B ⊥ which can be treated as a real magnetic field, while T is the 316 source part originating from the 3D effects.Correspondingly, B ⊥ 317 can be written as and is the symplectic matrix.Because the second term ∇ p G is 320 isotropic, the anisotropic properties of B ⊥ are completely deter-321 mined by the first term B ⊥ .As a 2D magnetic field, the integra-322 tion curves of B ⊥ are contour lines of its Hamiltonian F.

323
Following the method by Parnell et al. (1996), we define a where the rotation angle θ satisfies tan 2θ = −2p/q.Then F can 326 be transformed to 327 Therefore, the field lines of B are hyperbolic lines if |J 3 | < J thres 328 (Fig. 1a), are concentric elliptical curves if |J 3 | > J thres (Fig. 1b), 329 and are parallel straight lines if eig To quantitatively evaluate the geometric shapes of B ⊥ lines,

331
we define an angle, For |J 3 | < J thres , θ eig is the acute angle spanned by the separatrix 333 lines (see Fig. 1a).For |J 3 | > J thres , θ eig is the angle spanned by 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.

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.

Step 2: Analyze local magnetic structures
First, approximate the 3D magnetic Jacobian matrix D B on each grid via two-order central difference as where l, m, and n are indices of x 1 , x 2 , and x 3 , 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 êlmn 3 = B 1,lmn , B 2,lmn , B and its inverse matrix is T lmn = T lmn T .In the local frame, where ∆L is the spatial step length of the original discrete data.

398
The values of E at these positions can be obtained via linear  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 2 , the second dimension is x1 , 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 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 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 1 , x 2 , and x 3 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 lmn B , êlmn 1 , êlmn 2 , êlmn 3 , 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.

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 × n 2 × n 3 data

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

Benchmark 1: 3D turbulent reconnection within a
Harris sheet

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: 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 /2 and I is the identity matrix.e = p/ (γ − 1)+ρu 2 /2+B 2 /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. where B 0 = 1 and λ = 0.1 is the half-width of the current sheet.

515
The initial mass density and pressure are uniformly set as 1 and 516 0.05, respectively, corresponding to a β = 0.1.We set a uniform 517 background resistivity η = 5×10 −6 to obtain a Lundquist number 518 S = L 0 u 0 /η = 2 × 10 5 .To trigger fast reconnection, a perturba-519 tion on the z-direction of the magnetic vector field is placed at 520 the center of the current sheet (also see Ye et al. 2020), namely, 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 −3 , 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 531 set as 1200, 240, and 96 on x, y, and z directions, respectively.

532
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 2order 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.

Evolution overview
The evolution presents a standard picture of the Harris-sheet reconnection (see the animation of Fig. 3).In the beginning, trig- gered by the perturbation field Ãz , the center of the current sheet tion 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.

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 tearingmode 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 2 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 modeldependent 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  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 reconnection 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 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 fieldline 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 2 E × 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 mag- 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 B z is overlapped at the bottom.netic 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 concentrate (see Pontin et al. 2007), our method is still useful near null points. 2 / B lmn .In the Harrissheet 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 B ⊥ almost overlaps with that of R ∇B (see the blue curves in 796 Fig. 10), the peak of R B is about one order smaller than R B ⊥ , 797 and the magnitude of R B is much smaller even compared with 798 R B .Therefore, in this simulation, the variation of magnetic field 799 in a cell is dominated by B ⊥ for most reconnection grids.R ∇B , 800 R B ⊥ , R B , and R B can be directly calculated by the 1-10 columns 801 of RDInfo.ExtraData (see Table 3), which is helpful for users 802 to evaluate the validations of MPP approximation and the identi-803 fied local magnetic structures.We have not added the four ratios 804 to the output of ARD for the sake of saving storage space.

805
Third, our method provides the reconnection properties of 806 the finest objects in discrete data (i.e., the grids), which, how-807 ever, cannot determine the global properties of small-scale re-808 connection regions containing a series of adjacent grids.Based 809 on the information of grids, it is possible to further recognize 810 local finite-volume reconnection regions through clumping al-811 gorithms, which might produce more useful information includ-812 ing volume-integrated parameters and topological properties.We 813 will attempt to add these functions in the future to make the 814 code more powerful in dealing with multi-scale reconnection 815 processes.816 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 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 J 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 Third, we rewrite the ∂N β /∂s and ∂N α /∂s terms in Eq. 5.
The general magnetic reconnection theory 146 We first briefly summarize the main concepts of the general 147 magnetic reconnection theory developed by Hesse & Schindler 148 (1988).Nonvanishing magnetic field can be (at least locally) ex-149 pressed by the Euler potentials as B = ∇α×∇β.The three vectors 150 ∇α, ∇β, and b = B/ |B| compose a local frame system (α, β, s), 151 where s denotes the arc coordinate along a field line.An (α, β)-152 pair corresponds to a field line.The gradients of the Euler poten-153 tials, ∇α and ∇β, span the local plane normal to B but, generally 154 speaking, they are not orthogonal or normalized.155 If setting the vector potential as A = α∇β satisfying the 156 gauge A • B = 0, Ohm's law E = −v × B + N and Faraday's 157 law produce the evolution rule of the Euler potentials as 158 Magnetic structures at reconnection sites 249As the magnetic null points have been well studied (seeParnell 250  et al. 1996), here we analyze the local magnetic structures near 251 reconnection locations with finite magnetic field.To begin with, 252 we define a local frame at the origin r 0 by three orthogonal unit 253 vectors satisfying ê3 = b (r 0 ), ê2 ⊥ ê3 , and ê1 = ê2 × ê3 .In this 254 frame, the magnetic field at dr = r − r 0 can be approximated by 255

B
⊥ and B compose the magnetic field lying in the dr 3 = 0 plane, while B ≡ dr 3 (∂B/∂dr 3 ) represents the increments at locations of dr 3 0. Hereafter, we analyze the local in-plane magnetic structures on the dr 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, B and B is not investigated in detail.Suppose that the magnitudes of B and B are small enough compared with |B ⊥ |, the magnetic structure within a finite volume near r 0 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 11 +M 22 = −∂B 3 /∂dr 3 , 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 0 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.
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 0 (see Fig. B.1b, c).Types 4 and 5 are respectively source and sink structures dominated by the trace of M (see Fig. 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 con-

Fig. 1 .
Fig. 1.Cartoons of B ⊥ lines of X (a) and O (b) types.The blue curves depict field lines.The definitions of θ eig are labeled by orange markers.

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

399
interpolation of original data.If the values of E have the same 400 signs at r, R −0 , R +0 , R 0− , and R 0+ , and E at r lmn is maximum 401 among the five positions, then we identify r lmn as a 2D extreme 402 point of E .

403 4 .
The ARD function in LoRD toolkit 404 This method has been implemented as an easy-to-use func-405 tion ARD, integrated into our open-source Matlab toolkit project 406 named LoRD (Locate Reconnection Distribution) 1 .To use the 407 functions in LoRD, one just needs to add the directory into the 408 Matlab environment by the following command: 409 path( '<Download D i r e c t o r y >/LoRD / m a t l a b ' ,path ); 410 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); d o u b l e ' ; 472 LIM = [x1s ,x1e ,x2s ,x2e ,x3s ,x3e ]; 473 [B1 ,x1 ,x2 ,x3] = Tool_LoadData_Bin ( ' B1 .b i n ' ,DIM ,Precision ,LIM ); 474 The binary file should save the 3D data of B 1 in a 1D array, in 475 which the x 1 -direction index changes fastest and the x 3 -direction 476 index changes slowest.Precision is a string setting the class 477 and size of bits to read which should be consistent with the style 478 of the binary file.Available parameters of Precision include 479 'double', 'float32', etc.. LIM is an optional parameter defin-480 ing the coordinate limitations in three directions.Here, x1s and 481 x1e define the start and end of x 1 -coordinate, respectively, and 482 x2s, x2e, x3s, and x3e have similar definitions.Without LIM, 483 this function will simply assign the outputs x1, x2, and x3 with 484 integer grid indices.The output variables B1, x1, x2, and x3 can be directly passed to the ARD function.

Fig. 3 .
Fig. 3. Distributions of mass density ρ (a), temperature T (b), parallel current density J (c), and 3D X and O-type grids (d) at t = 9.For ρ 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 to 9 is available.

Fig. 4 .
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 .
Fig. 5. Evolutions of grid numbers for different types of structures.Only grids with J > 50 are counted.

Fig. 6 .
Fig. 6.Number densities of J (a) , P oh (b), θ eig (c), and R tr (d) of the 3D X-type grids at different moments.Different colors correspond to different moments.

Fig. 7 .
Fig. 7. Number densities of J (a), P oh (b), θ eig (c), and R 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 .

Fig. 8 .
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 B z is overlapped at the bottom.
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 ∇B ≡ ∆L D lmn B / B lmn 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 B. In applications, we define three ratios to evaluate the importance of the three terms, namely, R B ⊥ ≡ ∆L M lmn / B lmn , R B ≡ ∆L D 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 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 J 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 J and θ eig .Panel (d) has the same thresholds as (c) but also plots 3D O-type grids with orange dots.

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

Fig. C. 2 .
Fig. C.2. Number densities of J (a), P oh (b), θ eig (c), and R 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 .
the trajectories of extremal-Ξ lines are not necessary to always 218 pass through such locations.E and ∇ ⊥ E are two useful and 219 easy-to-calculate local parameters that can provide implications 220 about the locations of global reconnection effects, though they 221 are not equivalent to the original integral definitions in Eqs. 1 222 and 3.The local reconnection effects, felt by plasma elements in-224 side nonideal regions, correspond to the local violation of line 225 conservation, which are blocked out by the global reconnection 226 theory but still worth investigating.The equation governing the 227 local effects can be obtained by performing ∂/∂s on both sides 228 of Eq. 2, namely, 223 3,lmn / B lmn .êlmn 2 can be any unit vector perpendicular to êlmn 3 .As a convenient choice, we set it as 374 êlmn 2 = B 2,lmn , −B 1,lmn , 0 / B 1,lmn 2 + B 2,lmn 2 .For the special 375 case of B 1,lmn = B 2,lmn = 0, let êlmn . 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 lmn B might have small trace components resulting from numerical errors, ARD provides an option to clean the trace by D lmn

Table 2 .
Key configuration fields of the Parameters structure.

Table 3 .
Data saved in RDInfo.ExtraData.Γ is defined in Appendix A.