Three-dimensional plasmoid-mediated reconnection and turbulence in Hall magnetohydrodynamics

Magnetic reconnection is a fundamental process that alters the connectivity of magnetic field lines, releasing stored magnetic energy. The magnetic energy released by reconnection is transformed into kinetic, thermal, and non-thermal energy of plasma in explosive events in nature and laboratories, such as geomagnetic substorms, solar flares, coronal mass ejections (CMEs), gamma-ray bursts, and sawtooth crashes in fusion plasmas.^1–7 

Magnetic reconnection occurs at current sheets. In recent years, there has been growing evidence that large-scale reconnection current sheets will likely become fragmented due to the plasmoid instability.^8,9 The plasmoid instability has been found in a wide range of plasma models, including resistive magnetohydrodynamics (MHD),^8–23 Hall MHD,^16,24,25 and fully kinetic particle-in-cell (PIC) simulations.^26–31 

High-resolution two-dimensional (2D) resistive MHD simulations indicate that the plasmoid instability leads to a fast reconnection rate that is nearly independent of the resistivity.^9,14,18 In three-dimensional (3D) reconnection with a guide field, in addition to modes that are symmetric along the out-of-plane direction, oblique tearing modes can become unstable.³² Interaction of oblique tearing modes leads to self-generated turbulent reconnection without the need for external forcing.^28,33–35

Taking into account two-fluid and kinetic physics, the plasmoid instability can trigger even faster collisionless reconnection if the fragmented current sheets reach a kinetic scale $δi$ that corresponds to the ion skin depth $di$ or the ion sound Larmor radius $ρs$⁠, the latter being for cases with a strong guide field.^16,24,27 This effect leads to the consideration of reconnection “phase diagram” in the literature.^16,21,36,37 The phase diagram organizes various types of reconnection in the parameter space of two dimensionless parameters: the Lundquist number $S≡LVA/η$ and the system size to the kinetic scale ratio $Λ≡L/δi$⁠. Here, L is a characteristic length scale of the reconnection layer along the outflow direction, $VA$ is the Alfvén speed of the reconnecting magnetic field component, and $η$ is the magnetic diffusivity. Within resistive MHD, the reconnection current sheet becomes unstable when the Lundquist number S exceeds a critical value $Sc∼104$⁠.^9,38 Once the current sheet becomes unstable, the fragmentation continues until the secondary current sheets become marginally stable. This condition yields the typical widths of the secondary current sheets to scale as $δc∼LSc1/2/S$⁠.¹⁴ Comparing the secondary current sheet width $δc$ with the kinetic scale $δi$ yields a heuristic criterion $Λ<S/Sc1/2$ for the plasmoid-mediated onset of collisionless reconnection.

Based on simple estimates, we can see that this onset criterion is satisfied for numerous types of reconnection events in the solar atmosphere and astrophysical systems.³⁶ For example, consider post-CME current sheets in the corona^39,40 and ultraviolet (UV) bursts in the transition region.^41–43 Typical lengths of post-CME current sheets are on the order of $∼109m$⁠, Lundquist numbers $S∼1014$⁠, and $di∼1m$⁠; hence, the ratio $Λ=L/di∼109$⁠, which is smaller than $S/Sc1/2∼1012$⁠. For UV bursts, the lengths are estimated to be on the order of $105m$⁠, the Lundquist numbers $S∼1010$⁠, and $di∼0.1m$⁠. The ratio $Λ∼106$ is smaller than $S/Sc1/2∼108$⁠. Therefore, plasmoid-mediated onset of collisionless reconnection is potentially important for both types of reconnection events.

Evidence of plasmoids in both types of reconnection events is abundant. In post-CME current sheets, plasmoid-like structures have been regularly observed as moving blobs.^39,40,44 In UV bursts, although plasmoids are not directly visible, their existence may be inferred from the shapes of emission line profiles, which provide information on reconnection outflow through the Doppler shift. The emission line profiles of UV bursts typically have a triangle shape, which is consistent with a highly fluctuating reconnection outflow associated with plasmoids. In contrast, a single-X-line reconnection geometry will predict a double-peak line profile, which is inconsistent with observation.^41,43

To model the transition from collisional to collisionless reconnection mediated by the plasmoid instability, particle-in-cell (PIC) simulation provides a first-principles, fully kinetic description.^27,31 However, PIC simulations, especially in 3D, are usually limited to relatively small system sizes due to the computational cost. To simulate systems of larger sizes, this study employs a Hall MHD model incorporating two-fluid effects through the Hall terms in the generalized Ohm's law. Hall MHD is often regarded as a minimal model that captures important aspects of collisionless or weakly collisional reconnection beyond resistive MHD.⁴⁵ Previous 2D Hall MHD studies have demonstrated the onset of Hall reconnection mediated by plasmoid instability.^16,24

From the perspective of plasmoid formation, 2D Hall MHD reconnection is somewhat anomalous: it tends to settle to a single-X-line geometry after expelling all the plasmoids from the reconnection layer. This feature makes it qualitatively different from fully kinetic PIC simulations, where new plasmoids are constantly generated.²⁷ The single-X-line geometry of Hall reconnection also appears to be inconsistent with observations of post-CME current sheets and UV bursts, which indicate that plasmoids continue to form throughout the period of the events. Although non-single-X-line 2D Hall reconnection has been demonstrated, it appears to require a fairly large system size (⁠ $L/di≳5000$⁠) and may only exist in a narrow parameter space of the phase diagram.¹⁶ 

The observed persistence of plasmoid formation throughout solar reconnection events raises an important question regarding whether Hall MHD can adequately describe these phenomena. However, it is important to note that the preference for single-X-line reconnection in Hall MHD may result from the assumed 2D symmetry. Will 3D systems allow the formation of complex 3D structures and significantly change the picture? Moreover, if 3D Hall MHD reconnection geometry deviates from a single X-line, can it support self-generated turbulent reconnection?

To address these questions, we perform 3D Hall MHD simulations of reconnection and compare the results with the corresponding 2D simulations. We aim for large system sizes, up to $L/di=1000$⁠. To mimic natural reconnection events, we start from an initial current sheet that is significantly thicker than $di$ and let it thin down self-consistently. After the plasmoid instability triggers the onset of Hall reconnection, we continue to follow the evolution until it reaches a saturated phase.

This paper is organized as follows: Section II describes the governing equations and the simulation setup in detail. In Sec. III, we present simulation results of reconnection geometry and reconnection rate for various 2D and 3D runs. Section IV looks more deeply into the characteristics of the fully developed reconnection state by examining the energy power spectra and two-point structure functions of kinetic and magnetic fluctuations. Finally, we discuss the implications of our findings for large-scale astrophysical reconnection and conclude in Sec. V.

The normalization of Eqs. (1)–(5) is based on constant reference values of the density $n0$⁠, and the magnetic field $B0$⁠. Lengths are normalized to the system size L, and time is normalized to the global Alfvén time $tA=L/VA$⁠, where $VA=B0/4πn0mi$ and $mi$ is the ion mass. The normalization of physical variables is given by (normalized variables $→$ expressions in Gaussian units): $ρ→ρ/n0mi$⁠, $B→B/B0$⁠, $E→cE/B0VA$⁠, $v→v/VA$⁠, $p→p/n0miVA2$⁠, $J→J/(B0c/4πL)$⁠, and $di→di/L≡mic2/4πn0e2/L$⁠.

The initial plasma density and pressure are both uniform, with $ρ=1$ and $pi=pe=1$ in normalized units. The heat capacity ratio $γ=5/3$ is assumed. Perfectly conducting and free-slip boundary conditions are imposed along both x and y directions. Specifically, we have $v·n̂=0$⁠, $n̂·∇(n̂×v)=0$⁠, and $B·n̂=0$ on the boundaries. Here, $n̂$ is the unit vector normal to the boundary. The z direction is assumed to be periodic.

This model system is numerically solved with a massively parallel HMHD code. The numerical algorithm⁴⁸ approximates spatial derivatives by finite differences with a five-point stencil in each direction, augmented with a fourth-order numerical dissipation for numerical stability. The time-stepping scheme can be chosen from several options, including a second-order accurate trapezoidal leapfrog method and various strong-stability-preserving Runge–Kutta methods.^49,50 We employ the second-order accurate trapezoidal leapfrog method in this study.

We perform two sets of simulations with the ion skin depth $di=0.002$ and 0.001. The plasma resistivity and viscosity $ν$ are both set to a fixed value $5×10−6$⁠. Using the box size as the length scale L, the system size to the ion skin depth ratio $L/di$ is 500 and 1000, respectively. The Lundquist number $S=VAL/η=2×105$ and the magnetic Prandtl number $Prm≡ν/η=1$⁠. Because Hall MHD tends to develop fine structures at small scales, we add a small hyper-resistivity $ηH=10−13$ to smooth fluctuations at grid scales. The initial velocity is seeded with a random noise of amplitude $10−3$ to trigger the plasmoid instability. We perform 2D and 3D simulations with the same parameters and compare the results.

The 3D simulation mesh size is $Nx×Ny×Nz=2000×1000× 2000$ (⁠ $2000×1000$ for 2D), where the grid sizes are uniform along both the x and z directions, and packed along the y direction around the midplane to resolve the reconnection layer better. The grid size along the y direction is $Δy=10−4$ near the midplane (⁠ $y=0$⁠). Moving away from the midplane, the grid size gradually increases and reaches $Δy=0.005$ near the boundaries at $y=±1/2$⁠.

Figure 2 and the associate animation show the time evolution of the 2D simulation with $di=0.002$⁠, corresponding to $L/di=500$⁠. Here, the colormap denotes the reconnection outflow velocity $vx$⁠. The initial thinning of the current sheet is accompanied by the formation of reconnection outflow jets. At around $t=0.6$⁠, several plasmoids start to form in the reconnection layer. Subsequently, the onset of Hall reconnection expels all the plasmoids, and the reconnection geometry settles to a single-X-line configuration at $t=1.5$⁠. The 2D simulation with $di=0.001$ also exhibits qualitatively similar behavior, as shown in the animation in the supplementary material. These two runs demonstrate the strong tendency of 2D Hall MHD to form a single-X-line reconnection geometry, similar to the results of previous studies.^16,24

Three-dimensional (3D) simulations of the same parameters show different behaviors. In the time evolution of the first case with $di=0.002$⁠, shown in Fig. 3, the thinning current sheet first becomes unstable, forming flux ropes. The interaction between flux ropes leads to complex dynamics and the formation of chaotic field lines. However, all the flux ropes are eventually ejected and the magnetic field self-organizes into a nearly single-X-line configuration. The final state is similar to its 2D counterpart, even though a translational symmetry is not imposed. However, for the second case with $di=0.001$ (⁠ $L/di=1000$⁠), shown in Fig. 4, the magnetic configuration does not settle to a single X-line. Instead, it appears to develop a self-generated turbulent state similar to previous resistive MHD and PIC simulations.

To quantify the 3D reconnection rate, we first average the magnetic field along the z direction to obtain the mean magnetic field $B¯$⁠, which now depends only on x and y. We then use the mean field to calculate the reconnection rate in the same manner as the calculation for the 2D reconnection rate. Specifically, because the mean magnetic field $B¯$ is divergence-free, we can construct a mean-field flux function $ψ¯$⁠, subject to the boundary condition $ψ¯=0$ on the conducting-wall boundaries, such that $B¯=Bz¯ẑ+ẑ×∇ψ¯$⁠. We then identify the dominant X-point of the mean-field flux function and its corresponding value of the flux function $ψ¯X$⁠. The reconnection rate is determined by $dψ¯X/dt$⁠.

Figure 5 shows the time histories of reconnection rates for these 2D and 3D runs. The first three cases, all of which end up in a single-X-line configuration, yield faster reconnection rates. In comparison, reconnection in the fourth case, which does not settle to a single-X-line configuration, is slower. For reference, we also add a 3D resistive MHD (i.e., $di=0$⁠) run, which gives a substantially slower reconnection rate.

Note that properly measuring reconnection rates in complex 3D geometries remains a challenge, for which various methods have been proposed (see, e.g., Ref. 30 and references therein). For background configurations exhibiting translation symmetry, our approach of utilizing the mean magnetic field to calculate the reconnection rate provides a straightforward and effective estimate. The obtained reconnection rates qualitatively align with the trends observed in magnetic energy decay rates, $−dEB/dt$⁠, which quantify the rate at which magnetic energy is converted into other energy forms. Here, the magnetic energy $EB=∫B2/2d3x$⁠, integrated over the entire simulation box. Figure 6 presents the magnetic energy decay rates for the same 2D and 3D simulations as in Fig. 5. The results indicate that single-X-line Hall reconnection is more efficient at releasing magnetic energy compared to 3D Hall MHD turbulent reconnection, while resistive MHD exhibits the slowest energy release.

We further investigate whether the reconnection layer evolves to a self-generated turbulent state in the two 3D Hall MHD simulations by examining the energy spectra and two-point structure functions of the kinetic and magnetic fluctuations within the layer. Special treatments are necessary for calculating energy spectra and structure functions due to the strong inhomogeneity and shear in both magnetic and velocity fields. We adopt the procedure of Huang and Bhattacharjee,³³ which is summarized below.

To calculate the kinetic energy spectrum, we first define a weighted velocity field $w≡ρ1/2v$ such that the kinetic energy density is a quadratic form $w2/2$⁠. We then decompose $w$ into the sum of the mean field $w¯$ and the fluctuation $w̃$⁠, where the mean field is defined as the average along the z direction. The total kinetic energy is the sum of the mean-field kinetic energy $Ek¯≡∫w¯2/2 d3x$ and the fluctuation part of the kinetic energy $Ek̃≡∫w̃2/2 d3x$⁠. Likewise, we decompose the magnetic field $B$ into the mean field $B¯$ and the fluctuation $B̃$⁠. The total magnetic energy is the sum of the mean-field magnetic energy $Em¯≡∫B¯2/2 d3x$ and the fluctuation part of the magnetic energy $Em̃≡∫B̃2/2 d3x$⁠. We only consider the fluctuation parts for the calculation of the kinetic and magnetic energy spectra.

Before we calculate the turbulence energy spectra within the reconnection layer, we first multiply the fluctuating fields $B̃$ and $w̃$ by a $C∞$ Planck-taper window function,⁵¹ which equals unity within the range $−0.2≤x≤0.2$ and tapers off smoothly to zero over the ranges where $0.2≤|x|≤0.4$⁠. This step reduces the influence of irrelevant fluctuations from downstream exhaust regions, allowing us to focus on turbulence within the reconnection layer. Next, we compute the discrete Fourier energy spectra of $Ek̃$ and $Em̃$ using the “windowed” variables on each constant-y slice to obtain 2D energy spectra as functions of the wave numbers $kx$ and $kz$⁠. Then, we integrate the 2D energy spectra over the reconnection inflow direction y, which is also the direction of the strongest inhomogeneity, within the range $−0.05≤y≤0.05$⁠. Finally, we calculate one-dimensional (1D) spectra as the functions of $k≡kx2+kz2$ by integrating over the azimuthal direction on the $kx−kz$ plane.

Figure 7 shows the time evolution of the resulting 1D spectra of $Ek̃$ and $Em̃$ for the case of $di=0.002$⁠. Both spectra exhibit qualitatively similar behaviors and remain close to equipartition between the two. When the plasmoid instability occurs at about $t=0.6$⁠, the energy spectra peak at $k≃100$⁠. Subsequently, the energy cascades in both forward and inverse directions to smaller and larger scales and quickly fills a broad range of scales at around $t=1.0$⁠. However, because the reconnection geometry eventually evolves to a single-X-line configuration, the energy spectra do not settle to a quasi-steady state. Instead, the fluctuation parts of the energies peak at around $t=1.0$⁠, then gradually decay as the reconnection geometry becomes increasingly quasi-2D.

In comparison, the time evolution of the energy spectra for the case with $di=0.001$⁠, shown in Fig. 8, is qualitatively similar to the previous case at an early stage. However, because this case does not settle to a single-X-line configuration, the reconnection layer appears to realize a self-generated turbulent state, where the energy spectra become nearly time-independent after $t=1.7$⁠.

During the quasi-steady phase, the kinetic and magnetic energy spectra in Fig. 7 approximately follow the $k−5/3$ power law in the MHD range when $kdi<O(1)$⁠. The spectra steepen in the sub- $di$ range when $kdi>O(1)$⁠. As a reference, we plot the $k−11/3$ power law in the sub- $di$ range, which is the prediction of Galtier and Buchlin.⁵² However, because we do not have a sufficient separation between the ion skin depth $di$ and the dissipation scale, the energy spectra do not exhibit a definitive power law in the sub- $di$ range. The same two power laws are also shown in Fig. 7 for reference, even though the energy spectra do not reach a quasi-steady phase. Notably, the energy spectra in this case also steepen in the sub- $di$ range.

Next, we further investigate the alignment of turbulence eddies with the local magnetic field by calculating two-point structure functions of the kinetic and magnetic fluctuations in terms of the parallel displacement $r∥$ and the perpendicular displacement $r⊥$ relative to the local magnetic field; i.e., $F2w(r∥,r⊥)≡⟨|w(x+r)−w(x)|2⟩$ and $F2B(r∥,r⊥)≡⟨|B(x+r)−B(x)|2⟩$⁠. Structure functions allow us to measure the scale dependence of the anisotropy of turbulence eddies. In MHD turbulence, Goldreich and Sridhar (GS) theory^53,54 predicts that turbulence eddies become increasingly more elongated along the magnetic field at smaller scales. More precisely, GS theory predicts a scale-dependent relationship $k∥∼k⊥2/3$ between the wavenumbers parallel (⁠ $k∥$⁠) and perpendicular (⁠ $k⊥$⁠).

The GS theory relies on the critical balance condition, which assumes that the timescale for the nonlinear energy cascade of an eddy in directions perpendicular to the magnetic field balances the Alfvén wave propagation timescale of the eddy along the field. However, in Hall MHD, the dispersive whistler waves at sub- $di$ scales alter the parallel wave propagation time and may cause different anisotropic scaling relations.⁵⁵ At MHD scales when $k∥di≪1$⁠, the whistler waves become the shear Alfvén wave, and the GS theory may still be applicable.

The GS scaling relation $k∥∼k⊥2/3$ has been confirmed by Cho and Vishniac for homogeneous MHD turbulence using two-point structure function analysis.⁵⁶ In self-generated turbulence within reconnection layers, however, previous resistive MHD studies have shown deviations from the GS scaling relation.^33,34 Here, we explore whether similar behavior is observed in MHD-scale eddies for Hall MHD turbulent reconnection. In particular, we investigate whether our simulations of inhomogeneous reconnection layers in Hall MHD also exhibit deviations from the GS scaling relation.

Here we adopt the procedure of Huang and Bhattacharjee,³³ which is a modification of the procedure of Cho and Vishniac⁵⁶ to accommodate the inhomogeneous background of the reconnection layer. Due to the strong inhomogeneity along the y direction in our system, we calculate the structure functions for each x–z plane, restricting displacements to lie within these individual planes, rather than using the full 3D space. For a pair of points, we define the local magnetic field as the average of the magnetic fields from the two points. To ensure consistency with displacements limited to individual x–z planes, we measure the parallel and perpendicular components of the displacement relative to the in-plane component of the local magnetic field. We compute the structure functions by averaging over $109$ random pairs of points; the x coordinates of these points are limited to the range $−0.25≤x≤0.25$⁠.

We focus on the case with $di=0.001$ because it appears to realize a self-generated turbulent state in the reconnection layer. In the following discussion, we present the results from the snapshot at $t=2.1$⁠, but other snapshots during the quasi-steady period show similar behaviors. Figure 9 summarizes the outcome of this analysis. Here, the rows correspond to different x–z slices. The first row shows the results at $y=0$⁠, the second row at $y=0.002$⁠, and the third row at $y=0.005$⁠. The first column of the figure shows the magnitude of the magnetic field fluctuation $B̃$⁠, overlaid with streamlines of the in-plane magnetic field. The magnetic field fluctuations form elongated eddies along the direction of the local magnetic field. The second column shows contours of the structure function $F2B(r∥,r⊥)$⁠, which clearly show elongated eddies along the local magnetic field direction. We have also performed the same analysis for the weighted velocity field fluctuation $w̃$ and the corresponding structure function $F2w(r∥,r⊥)$⁠. The results are qualitatively similar to those of the magnetic field fluctuations; therefore, they are not shown. The third column shows the relations between the semi-major axes $r∥∼k∥−1$ and the semi-minor axes $r⊥∼k⊥−1$ of contours of the structure functions $F2B(r∥,r⊥)$ and $F2w(r∥,r⊥)$⁠. These relations quantify the scale dependence of turbulent eddy anisotropy. The two dashed lines represent the scalings of $k∥∼k⊥$ (scale-independent) and $k∥∼k⊥2/3$ (GS theory), for reference. In the first two rows at $y=0$ and $y=0.002$⁠, relatively close to the mid-plane, the relations between $r∥$ and $r⊥$ appear to be approximately scale-independent. This result is similar to previous resistive MHD results.^33,34 Interestingly, as we move further away from the mid-plane to $y=0.005$ (last row), the GS scaling relation $k∥∼k⊥2/3$ is partially recovered and appears to be consistent with eddies much larger than the $di$ scale.

In conclusion, our simulations show that the occurrence of single-X-line Hall reconnection in 3D is less probable compared to 2D. A case that evolves to a single-X-line geometry can become turbulent with multiple reconnection sites in 3D. Notably, the single-X-line geometry consistently yielded the highest reconnection rate among all studied cases.

Whether 3D Hall MHD reconnection realizes a self-generated turbulent state or a single-X-line configuration depends on the interplay of plasma parameters, including the Lundquist number S and the system size to kinetic scale ratio $Λ=L/di$⁠. This study shows that at a fixed Lundquist number $S=2×105$⁠, a system of $Λ=500$ evolves to a single-X-line geometry, whereas a system of a larger system size $Λ=1000$ evolves to a turbulent reconnection state. However, whether the system with $Λ=1000$ can lead to a single-X-line configuration in even higher Lundquist numbers remains an open question. To comprehensively understand the diverse dynamic behaviors of reconnection across parameter space, a systematic investigation is imperative. Such research will be crucial in advancing our knowledge of reconnection processes in complex plasma systems.

The self-generated turbulent reconnection in 3D Hall MHD exhibits significant qualitative differences compared to resistive MHD. The turbulent reconnection layer in Hall MHD is significantly broader, and the reconnection rate is much faster. Regardless of whether the reconnection geometry settles to a single X-line or becomes turbulent, the hallmarks of Hall reconnection remain visible. Notably, the reconnection exhaust regions open up and form a Petschek-like configuration (see Figs. 3 and 4). We emphasize that the resemblance with the Petschek-like configuration is just geometrical and nothing more because the underlying Hall MHD dynamics is qualitatively different from the resistive MHD dynamics. Nonetheless, our simulation results and future research for even larger system sizes may shed light on recent Parker Solar Probe in situ measurements of reconnection exhausts associated with interplanetary coronal mass ejections and heliospheric current sheets, where bifurcated current sheets resembling Petschek's reconnection model with a pair of slow-shocks or rotational discontinuities have been reported.⁵⁷ 

Many questions remain open, particularly for applications to natural reconnection events where the system sizes are much larger than those considered in this study (i.e., $Λ=L/di⋙1$⁠). Will a larger system size make self-generated turbulent reconnection easier to realize? Do global reconnection rate and geometry depend on the microscopic description of reconnection physics (e.g., Hall MHD or fully kinetic PIC) when the scale separation between L and $di$ is large? In the future, the Hall MHD description of self-generated turbulent reconnection should be further compared with fully-kinetic PIC simulations or more elaborated fluid models such as high-moment multi-fluid models,^58,59 in particular in the regime when the system size is much larger than kinetic scales.

The self-generated turbulence in the reconnection layer also poses new challenges to turbulence theories and future numerical simulations. Our Hall simulation shows that the turbulent state in the vicinity of the reconnection mid-plane does not satisfy the Goldreich–Sridhar (GS) scaling relation of turbulence eddy anisotropy. This result is similar to previous findings in resistive MHD simulations. However, away from the mid-plane, the GS scaling relation is partially recovered for large eddies. This phenomenon has not been observed in resistive MHD simulations and is not predicted by any theory.

Why the GS scaling relation is not satisfied in resistive MHD self-generated turbulent reconnection but is partially recovered in Hall MHD remains an open question. Previously, Huang and Bhattacharjee argued that the discrepancy may be because the reconnection layer is highly inhomogeneous, with the reconnection outflow and the magnetic field being strongly sheared.³³ In Hall MHD turbulent reconnection, the turbulent region is significantly thicker than the resistive MHD counterpart. Because the magnetic field is not as strongly sheared at planes away from the mid-plane, this observation could explain why the GS scaling relation is partially recovered in places away from the mid-plane.

Even within the MHD regime, whether self-generated turbulent reconnection deviates from the GS theory has been unsettled. In a previous study, Kowal et al. showed that the scaling relation for large-scale eddies in turbulent reconnection simulations satisfies the GS scaling, while small-scale eddies are “contaminated” by reconnection and exhibit different scaling behaviors.⁶⁰ However, the initial and boundary conditions in their study differ from those in the present study: their initial condition contains a tangential discontinuity rather than a smooth magnetic field; the strength of the guide field is considerably lower than that of the reconnecting component of the magnetic field; the boundary conditions are periodic in the outflow direction and open in the inflow direction. It is unclear how much of their findings can be ascribed to these distinctions.

Furthermore, Kowal et al. employed full 3D fields and displacements in their calculation of two-point structure functions, while our analysis is restricted to displacements within each x–z slice, utilizing in-plane magnetic field components for determining parallel and perpendicular displacements. The implications of this significant difference warrant further investigation. Our decision to use 2D slices stems from the pronounced inhomogeneity along the inflow (y) direction, which renders displacements in this direction qualitatively distinct from in-plane displacements. Preliminary attempts to calculate two-point structure functions using full 3D fields and displacements yielded results that were challenging to interpret. Future simulations with higher resolution, allowing for sufficient scale separation between the turbulent region thickness and the smallest turbulence eddies, should revisit this analysis to comprehensively assess the consequences of employing 2D vs 3D displacements.

It is worth mentioning that the present study uses a compressible code. The relatively high plasma $β$ (⁠ $≈4$ relative to the reconnecting component of the magnetic field) and the presence of a guide field keep the plasma close to, but not perfectly, incompressible. Throughout the simulations, the root mean square density fluctuation is approximately $3%$⁠, although the density can differ from the mean value by up to $10%$ at some locations. Even though the compressibility of the plasma is relatively low, the compressible waves within the system may introduce subtle effects on the energy cascade and eddy anisotropy.⁶¹ 

In order to assess the effects of compressibility, turbulent reconnection should be investigated with incompressible MHD and Hall MHD, as well as compressible systems with varying plasma $β$⁠. Furthermore, simulations with varying guide field strengths will clarify magnetic shear effects. In a broader context, it is also of great interest to further investigate how Hall and kinetic effects affect plasmoid-mediated energy cascade in homogeneous turbulence, which is an important topic of current research.^62–67 

See the supplementary material for an animation for the 2D simulation with di=0.001 in the reconnection layer, similar to that shown in Fig. 2.

This research was supported by the U.S. Department of Energy (Grant Number DE-SC0021205), the National Science Foundation (Grant Numbers AGS-2301337 and 2209471), and the National Aeronautics and Space Administration (Grant Number 80NSSC18K1285). Computations were performed on facilities at the National Energy Research Scientific Computing Center and the National Center for Atmospheric Research.

The authors have no conflicts to disclose.

Yi-Min Huang: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Resources (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Amitava Bhattacharjee: Conceptualization (equal); Funding acquisition (supporting); Resources (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.