Basic Pattern of Three-dimensional Magnetic Reconnection within Strongly Turbulent Current Sheets

Magnetic reconnection occurring in current sheets (CSs) is a fundamental plasma process to drive eruptive phenomena across the universe, during which magnetic energy is rapidly converted to expel fast bulk flows, heat plasmas, and accelerate charged particles (Masuda et al., 1994; Matsumoto et al., 2015; Shukla & Mannheim, 2020; Yang et al., 2024; Huang et al., 2024; Wang et al., 2022a; Ping et al., 2023). Early two-dimensional (2D) reconnection models, including the Sweet-Parker and Petschek models (Parker, 1957; Sweet, 1958; Petschek, 1964), assume a stationary and laminar CS. However, numerous studies over recent decades have highlighted the importance of three-dimensional (3D) effects, revealing that CSs can become highly fragmented and dynamic due to various instabilities, such as the tearing-mode instability (TMI) and Kelvin-Helmholtz instability (KHI), especially in astrophysical plasmas with large magnetic Reynolds numbers (Loureiro et al., 2007; Bhattacharjee et al., 2009; Daughton et al., 2011; Loureiro et al., 2013; Huang & Bhattacharjee, 2016; Kowal et al., 2020; Wang et al., 2023). On the other hand, magnetic reconnection can be inevitably coupled with turbulence (Lazarian & Vishniac, 1999), which not only increases the CS widths but also enhances the reconnection rate (Kowal et al., 2009, 2017; Yang et al., 2020).

Solar flares, the most energetic eruptive phenomena in the solar system, have recently been observed to exhibit numerous turbulent features. Fragmented and dynamic reconnection structures have been identified within the flare CSs that connect coronal mass ejections (CMEs) and flare loops (Warren et al., 2018; Cheng et al., 2018). These structures are likely the fundamental drivers of turbulent flows at the tops of flare loops (Kontar et al., 2011; McKenzie, 2013) and the fine structures observed in flare ribbons (French et al., 2021; Wyper & Pontin, 2021; Corchado Albelo et al., 2024; Thoen Faber et al., 2025). The specific physics of multi-scale 3D reconnection regions constituting turbulent flare CSs, serving as the “Rosetta Stone” of understanding specific energy release mechanisms and interpreting observed fine structures, unfortunately, remains poorly understood.

With the great advances in high-performance computers, high-resolution simulations have become an indispensable tool for investigating 3D turbulent reconnection (Ji et al., 2022). However, accurately identifying and analyzing reconnection regions within a turbulent CS from vast simulation data presents a new challenge, hindering a comprehensive understanding of 3D turbulent reconnection. In weak magnetic turbulence, where the turbulent magnetic field $\delta B$ is much smaller than the background field $B_{0}$, the problem can be simplified using a 2D approximation on planes perpendicular to the background field (Zhdankin et al., 2013; Li et al., 2021, 2023; Dong et al., 2022). In contrast, in strongly turbulent plasmas where $\delta B\geq B_{0}$, as sophisticated methods are absent (Vlahos & Isliker, 2023; Isliker et al., 2019; Kowal et al., 2020; Lapenta, 2021), determining the nature of 3D reconnection regions, including their geometric shapes, reconnection flows, and reconnection rates, is still in the early stage. To tackle this problem, we propose a novel methodology to quantitatively analyze 3D turbulent reconnection. We apply it to the simulation data of the strongly turbulent 3D reconnection during a solar flare, which is self-consistently formed and is comparable with observations in nature (Wang et al., 2023; Ren et al., 2025). Our systematic analysis uncovers a fundamental yet previously unresolved pattern of 3D turbulent reconnection in the flare CSs.

The simulation solves the resistive magnetohydrodynamics (MHD) equations incorporating necessary coronal effects (see Appendix A). The flare CS is $10^{8}\,\mathrm{m}$ in length and initially conforms to the standard flare model (Shibata et al., 1995; Lin et al., 2015), rooted in a high-density, low-temperature chromospheric layer (Wang et al., 2021, 2022b; Ye et al., 2020; Shen et al., 2022). Fast reconnection is initiated by a localized anomalous resistivity in the corona and is subsequently dominated by a low uniform resistivity $\eta_{b}$, corresponding to a background Lundquist number of $S=2\times 10^{5}$. A uniform mesh with a grid size of $\Delta L=26\,\mathrm{km}$ is used in the CS and flare loop top regions to guarantee the accuracy of physical results. Under the influences of tearing mode instability (TMI), kink instability (KI), and Kelvin-Helmholtz instability (KHI) in sequence, the CS finally evolves into a self-sustained strongly turbulent state, indicated by turbulent magnetic energy comparable to background magnetic energy and a spectrum presenting an inertial region with a power-law index close to $-5/3$ (Wang et al., 2023). Notably, this simulation reproduces various key observational features of solar flares (Wang et al., 2023), such as finger-like structures above the flare loop top (Hanneman & Reeves, 2014) and turbulent loop-top regions (McKenzie, 2013). Hereafter, we focus on the data within the CS region $y\in\left[0.45,1\right]$ at the final turbulent moment $t=8.2$. In simulation, the units of length, time, velocity, electric field, magnetic field, and current density are $L_{0}=5\times 10^{7}\,\mathrm{m}$, $t_{0}=114.61\,\mathrm{s}$, $u_{0}=4.36\times 10^{5}\,\mathrm{m\,s^{-1}}$, $E_{0}=8.73\times 10^{2}\,\mathrm{V\,m^{-1}}$, $B_{0}=0.002\,\mathrm{T}$, and $J_{0}=3.18\times 10^{-5}\,\mathrm{A\,m^{-2}}$, respectively.

According to the general magnetic reconnection theory, the condition for 3D reconnection is given by $\Xi\equiv\int E_{\parallel}\mathrm{d}s\neq 0$ (Schindler et al., 1988; Hesse & Schindler, 1988; Hesse et al., 2005; Pontin & Priest, 2022), where $\Xi$ is the quasi-potential, $E_{\parallel}=\eta_{b}\bf{J\cdot B}/\lvert\bf{B}\rvert$ denotes the parallel electric field, $\bf{J}$ is the current density, $\bf{B}$ is the magnetic field, and the integration is carried out along a magnetic field line. In a strongly turbulent state, reconnection regions in the CS characterized by intense $E_{\parallel}$ exhibit a highly chaotic pattern (see Figs.1a and 7). Although 3D reconnection generally occurs within a finite volume (Priest et al., 2003), among all the field lines threading this volume, there exists a special one associated with an extremal quasi-potential $\Xi_{\mathrm{max}}$, which corresponds to the reconnection rate (Hesse et al., 2005; Wyper & Hesse, 2015). For discrete numerical data, the spatial distribution of $\Xi_{\mathrm{max}}$ field lines can be approximately inferred from $E_{\mathrm{max}}$ grid points (see Appendix B), identified by the criteria $\nabla_{\perp}E_{\parallel}\sim 0$ and $\lvert E_{\parallel}\rvert>E_{\mathrm{thres}}$. Here, $\nabla_{\perp}$ denotes the gradient operator perpendicular to the local magnetic field, and $E_{\mathrm{thres}}=5\times 10^{-4}$ serves as a threshold distinguishing regions undergoing reconnection or not. It should be noted that the statistical results presented later are robust to reasonable variations in $E_{\mathrm{thres}}$, as detailed in Appendix E.

The $E_{\mathrm{max}}$ grids represent the smallest numerically resolvable units of reconnection. To investigate the integral properties of reconnecting regions, we employ the region growing algorithm to cluster spatially connected $E_{\mathrm{max}}$ grids, thereby constituting the kernels of reconnection regions (Fig. 1b). Here, “spatially connected” means that, for any given grid, its nearest neighbor resides within the same mesh cell. This approach allows the original complex turbulent reconnection structure shown in Fig. 1a to be decomposed into smaller, more analyzable components. The reconnection kernels exhibit diverse shapes, scales, and $E_{\parallel}$ values (Fig. 1b). Some extend considerable distances in specific directions, resembling the extended X-lines observed in kinetic simulations (Li et al., 2023). It is worth noting that the lower threshold, $E_{\mathrm{thres}}$, ensures all grids within a reconnection kernel have identical signs of $E_{\parallel}$. Among the field lines intersecting the $E_{\mathrm{max}}$ grids, there is one line that can be identified with the maximum value of $\lvert\Xi\rvert$, allowing precise determination of the $\Xi_{\mathrm{max}}$ line for the reconnection region (see the magenta curves in Fig. 1b). Unlike laminar reconnection, the $\Xi_{\mathrm{max}}$ line of a reconnection region may pass multiple reconnection regions (Fig. 1b).

For a given reconnection kernel, there exists an intrinsic reconnection frame (Fig. 1c), consisting of the guide field direction ($\hat{\bf e}_{g}$), the inflow direction ($\hat{\bf e}_{i}$), and the outflow direction ($\hat{\bf e}_{o}$). $\hat{\bf e}_{g}$ approximately aligns with the average magnetic field of $E_{\mathrm{max}}$ grids. $\hat{\bf e}_{i}$ is determined by the local magnetic field structures at all $E_{\mathrm{max}}$ grids and is approximately normal to the surface of the reconnection kernel (see Appendix C). $\hat{\bf e}_{o}$ is set to be perpendicular to both $\hat{\bf e}_{g}$ and $\hat{\bf e}_{i}$. Using this frame, the inflow edges ($\mathcal{E}^{u}_{i}$, $\mathcal{E}^{d}_{i}$) and outflow edges ($\mathcal{E}^{u}_{o}$, $\mathcal{E}^{d}_{o}$) surrounding a reconnection kernel can be identified (see Fig. 1c). Here, the superscripts “u” and “d” denote the upper and down sides relative to $\hat{\bf e}_{i}$ or $\hat{\bf e}_{o}$, respectively. The inflow edges correspond to the nearest isosurfaces with $\lvert E_{\parallel}\rvert=E_{\mathrm{thres}}$ enclosing the reconnection kernel. The outflow edges outline the boundaries of the reconnection region at both ends of the outflow direction, with their widths along the inflow direction determined by the average thickness between inflow edges. These edges enable us to quantitatively investigate reconnection flows and dimensionless reconnection rates of all fragmented reconnection regions. Technical details regarding the determination of the frame and edges are elaborated in Appendix C.

The volumes of reconnection kernels can be represented by their $E_{\mathrm{max}}$ grid numbers, $N_{g}$, which follow a power-law probability distribution function (PDF) decreasing towards larger $N_{g}$ (Fig. 2a). To quantitatively evaluate the shape of a reconnection kernel, we perform principal component analysis (PCA) on its grid coordinates, which outputs the coordinate variances on three principal directions denoted by $\mathcal{L}_{1}$, $\mathcal{L}_{2}$, and $\mathcal{L}_{3}$ in descending order. For reconnection kernels with finite $\mathcal{L}_{1}$, their shapes can be evaluated by $\mathcal{R}_{S}^{\alpha}=\sqrt{\mathcal{L}_{\alpha}/\mathcal{L}_{1}}$, $\alpha=2,\,3$. Figure 2b shows that almost all reconnection kernels satisfy $\mathcal{R}^{3}_{S}<\mathcal{R}^{2}_{S}$, with a large portion having $\mathcal{R}^{3}_{S}$ close to 0, implying that they tend to have 2D patch-like spatial distributions. For a reconnection kernel with $\mathcal{L}_{2}\neq 0$, the characteristic vector $\hat{\bf e}_{c}$ corresponding to the smallest variance $\mathcal{L}_{3}$ is approximately perpendicular to the surface of the reconnection kernel. The mean angle between $\hat{\bf e}_{c}$ and the magnetic field at $E_{\mathrm{max}}$ grids $\mathbf{B}_{gl}$, evaluated by $\theta_{B}=\langle\arccos(\hat{\bf b}_{gl}\cdot\hat{\bf e}_{c})\rangle$, concentrates around $90^{\circ}$ (Fig. 2c), indicating that reconnection kernels tend to form along the local magnetic flux surfaces. Here $\hat{\bf b}_{gl}$ is the unit vector of $\mathbf{B}_{gl}$ and $\langle\cdot\rangle$ in this paper denotes averaging over all $E_{\mathrm{max}}$ grids of a reconnection kernel.

The reconnection inflows and outflows associated with a reconnection kernel can be determined at its inflow and outflow edges (Fig. 1c). Specifically, the inflow velocity on the upper inflow edge $\mathcal{E}^{u}_{i}$ is calculated as $u^{u}_{i}=\hat{\bf e}_{i}\cdot[\mathbf{u}(\mathcal{E}_{i}^{u})-\langle\mathbf{u}\rangle]$, and $u^{d}_{i}$, $u^{u}_{o}$, and $u^{d}_{o}$ at the other three edges can be obtained similarly. Generally speaking, the average velocities of reconnection kernels ($\langle\mathbf{u}\rangle$) have finite values, representing their global motions. We collect the four velocities sampled at the edge grids for all reconnection kernels to obtain their PDFs, which exhibit the flow velocity distributions at the inflow and outflow edges (Fig. 3a). It is found that all four PDFs exhibit significant broadening and approximate symmetry with respect to their peaks (Fig. 3a). Their full widths at half maximums (FWHMs) are much larger than their mean values, indicating a strong turbulent regime. The broadening of PDFs reflects the inherent uncertainties in determining reconnection flows of the fragmented reconnection regions, presenting a distinct nature compared with traditional laminar reconnection with deterministic reconnection flows.

On average, the PDF of the upper inflow velocities drifts towards negative values, while the down inflow velocities tend towards positive ones, reflecting an inward flow along $\hat{\bf e}_{i}$; conversely, the drift directions of the velocities at the upper and down outflow edges are opposite to those of inflows, indicating an outward flow along $\hat{\bf e}_{o}$ (Fig. 3a). The mean inflow and outflow velocities for each reconnection kernel are given by $\langle u_{in}\rangle=\langle u^{u,j}_{i}\rangle-\langle u^{d,j}_{i}\rangle$ and $\langle u_{out}\rangle=\langle u^{u,j}_{o}\rangle-\langle u^{d,j}_{o}\rangle$, respectively. Their joint PDF shows that reconnection kernels are most probably located in the fourth quadrant with $\langle u_{in}\rangle<0$ and $\langle u_{out}\rangle>0$ (Fig. 3b). $\langle u_{in}\rangle<0$ means that the plasmas at inflow edges move towards the reconnection kernel, while $\langle u_{out}\rangle>0$ indicates that the plasmas are expelled out across the outflow edges. This average effect is consistent with the flow patterns surrounding the reconnection region in Sweet-Parker and Petschek models (Priest & Forbes, 2000). However, it is still possible for a reconnection kernel to enter the other three quadrants with a counter-intuitive flow pattern (Fig. 3b), which could be a new feature for 3D turbulent reconnection.

The reconnection rate of each reconnection kernel is given by $\Xi_{\mathrm{max}}$ (Hesse et al., 2005). Its dimensionless value can be obtained by $\mathcal{R}_{p}=\lvert\Xi_{\mathrm{max}}\rvert/(B_{in}V_{Ain}L_{in})$. Here, $B_{in}$ and $V_{Ain}$ are the inflow parameters calculated by the average magnetic field and Alfvénic velocity at the two inflow edges $\mathcal{E}^{u}_{i}$ and $\mathcal{E}^{d}_{i}$, and $L_{in}$, reflecting the length of CS analogous to the 2D model, is set as the standard deviation of the reconnection kernel along $\hat{\bf e}_{o}$. Similar to the velocities, the PDF of $\mathcal{R}_{p}$ also exhibits a strong broadening with the main body in the range of $0.01$–$0.1$ (see the black curve in Fig. 3c).

At a given moment $t$, the total reconnection rate contributed by all fragmented reconnected regions is given by $\mathcal{R}_{t}=\sum\lvert\Xi_{\mathrm{max}}\rvert/(B_{0}u_{0}L_{0})$ (Fig. 3d), where the summation is taken over all reconnection kernels. Traditionally, for a CS with approximate translation symmetry along the guide field direction, the reconnection rate $\mathcal{R}_{\bar{\bf B}}$ can be evaluated from the 2D magnetic field averaged along the guide field direction, $\bar{\bf B}=\langle\mathbf{B}\rangle_{z}$ (Wang et al., 2023; Huang & Bhattacharjee, 2016). The mean-field approach relies on identifying the principal X-point of $\bar{\bf B}$, which typically exhibits spatial jumps under turbulent conditions, leading to significant temporal fluctuations (see the gray curve in Fig. 3d). However, by including the contributions from all reconnection fragments, the evolution of $\mathcal{R}_{t}$ shows considerably smaller fluctuations and follows a global trend similar to $\mathcal{R}_{\bar{\bf B}}$ (see the dark curve in Fig. 3d). The curve of $\mathcal{R}_{t}$ exhibits two rising stages corresponding to the development of TMI during $6<t<6.5$ and the growth of turbulence after $t=7$, respectively. Furthermore, $\mathcal{R}_{t}$ is found to be systematically larger than $\mathcal{R}_{\bar{\bf B}}$, which shows that $\mathcal{R}_{\bar{\bf B}}$, as adopted previously, tends to underestimate the total reconnection rate.

To investigate the influences of reconnection kernel scales on statistical results, we categorize the reconnection kernels into three groups: $V_{1}$ ($N_{g}<5$), $V_{2}$ ($5\leq N_{g}<50$), and $V_{3}$ ($N_{g}\geq 50$) (Fig. 2a). The PDF of $\theta_{B}$ remains nearly unchanged across different reconnection kernel sizes (Fig. 2c). For the joint PDFs of $\langle u_{in}\rangle$ and $\langle u_{out}\rangle$, the proportions in the fourth quadrant ($\langle u_{in}\rangle<0$ and $\langle u_{out}\rangle>0$) are similar across all categories, with values of $50.4\%$, $52.2\%$, and $52.5\%$ for $V_{1}$, $V_{2}$, and $V_{3}$, respectively (Fig. 3b). The PDFs of $\mathcal{R}_{p}$ for the three categories exhibit similar trends (Fig. 3c). For $V_{2}$ and $V_{3}$, the PDFs almost overlap; for $V_{1}$, the peak shifts towards smaller values, but 41.7% of the samples still fall within the range of $\left[0.01,0.1\right]$. It should be noted that the reconnection kernels in $V_{1}$ contain fewer than 5 grids, approaching the resolution limit of the simulation, and are likely influenced by numerical errors. Nevertheless, the statistical trends observed in $V_{1}$ remain highly similar to those of the larger reconnection kernels. Therefore, the morphologies, flows, and rates of reconnection kernels can be considered approximately scale-independent.

We also perform a standard 3D box-counting analysis to determine the fractal behavior of the regions covering all reconnection kernels (see Appendix F for technical details). The relationship between the number of covering boxes ($N_{b}$) and their size ($L_{b}$) exhibits a power-law distribution, yielding a fractal dimension of $D_{F}=2.19$. This result suggests that the reconnection kernels exhibit a fractal structure in 3D space with a relatively low filling factor, further supporting the scale-independent quasi-2D patten of reconnection kernels (Fig. 4).

A novel approach is developed to quantitatively analyze 3D reconnection within strongly turbulent flare CS. We recognize multitudes of fragmented reconnection kernels and disclose their basic properties including the morphologies, reconnection in/out-flows, and reconnection rates, as well as the corresponding statistical distributions.

It is proved that the reconnection kernels tend to have a patch-like structure aligned with local magnetic flux surface. For the first time, it is shown that, due to strong turbulence, the fragmented reconnection regions produce fluctuated in/out-flows and reconnection rates, represented by the significant broadenings in their PDFs. Despite a strong fluctuation, on average, they most probably induce flows compressing plasmas towards the reconnection kernels and expel exhausts in the outflow direction. The reconnection rates mainly take values of 0.01–0.1, coinciding with the values derived by previous studies for different reconnection configurations (Cassak et al., 2017). More importantly, the statistical laws are found to be approximately independent of the reconnection kernel scales, indicating the existence of a fundamental reconnection pattern. At last, we propose a new strategy to evaluate the reconnection rate of the highly fragmented 3D turbulent reconnection. Compared with the mean-field method which only works for 3D systems with approximate translation symmetry on the guide field direction (Huang & Bhattacharjee, 2016), the new method can be applied to arbitrary 3D systems and estimate the reconnection rate $\mathcal{R}_{t}$ more accurately, thus assessing its intrinsic evolution.

A scenario of 3D reconnection within strongly turbulent flare CSs is suggested and illustrated by Fig. 5. The reconnection kernel with an irregular shape tends to have a 2D patch-like configuration. Owning to the turbulence, the reconnection inflows and outflows are violently perturbed (see the velocity arrows in Fig. 5). The $\Xi_{\mathrm{max}}$ line, a special field line threading the reconnection kernel, satisfies the line conservation (Hesse et al., 2005; Wyper & Hesse, 2015) and approximately moves along with the reconnection kernel (see the magenta curve in Fig. 5). If a field line passes through the diffusion region enclosing the reconnection kernel, its start point flows with the plasma (see the red dot in Fig. 5) and moves a small distance after a short time $\delta t$ but its endpoint might flip a large distance (see the blue curves in Fig. 5 or Fig. 7 of Priest et al. (2003)). This flipping is similar to the slipping of flare loops as frequently observed along flare ribbons (Aulanier et al., 2006, 2007). Thousands of reconnection fragments with various flow patterns and reconnection rates constitute the complex reconnection process in a turbulent CS, but their statistical distributions are independent of scales. The reconnection kernels collectively exhibit a fractal structure in 3D space with a fractal dimension greater than 2. This effectively generalizes the concept of fractal reconnection by Shibata & Tanuma (2001) to 3D turbulent reconnection and provides a quantitative validation.

The quasi-separatrix layers (QSLs) are commonly used to identify where the reconnection may occur (Pontin & Priest, 2022; Pontin et al., 2024), but are not able to pinpoint the reconnection (Reid et al., 2020). In contrast, the methodology developed here not only precisely locates the reconnection sites but also quantitatively determines their physical properties. Furthermore, while our analysis targets the collisional reconnection regime, the method that is developed from general reconnection theory and only requires magnetic field, non-ideal electric field, and velocities as inputs can be applied to kinetic simulations of collisionless reconnection. The resulting statistics of local reconnection properties can aid in initializing physical parameters for small-scale kinetic or hybrid simulations and provide strategies for particle injection in MHD-particle models (Zhang et al., 2024; Drake et al., 2019; Arnold et al., 2021; Sun & Bai, 2023).

The CS for eruptive flares in our simulation has a special configuration, with magnetic field lines rooted in the solar photosphere and extending outward. The guide field within the CS is aligned along the polarity inversion line (Janvier et al., 2014; Xing et al., 2024) and varies spatially and temporally in magnitude (Aulanier et al., 2012; Dahlin et al., 2022). This distinctive configuration plays a key role in determining the development of the turbulent CS and in forming fine flare features observed. Furthermore, despite the turbulent nature of the reconnecting CS varying with time in our simulation, the statistical properties derived above exhibit a limited variation (Fig. 11), suggesting the existence of a unified pattern for 3D turbulent reconnection in flare CSs. Nevertheless, more investigations are still encouraged to justify such a fundamental pattern for other reconnection regimes; the methodology developed here offers a valuable analysis tool.