Quantifying Reconnection and it's Dynamical Role in 2D Magnetic Rayleigh-Taylor Turbulence

TL;DR

This study uses high-resolution simulations to analyze the critical role of magnetic reconnection in the evolution and energy transfer processes of 2D magnetic Rayleigh-Taylor turbulence, revealing its significance in plasma dynamics.

Contribution

It introduces an automated reconnection detection algorithm and provides a comprehensive statistical analysis of reconnection's impact across different magnetic field strengths.

Findings

Reconnection facilitates plume merger and reduces magnetic tension.

Reconnection accounts for up to 80% of magnetic-to-kinetic energy transfer.

Reconnection contributes minimally to magnetic energy dissipation (~3%).

Abstract

Magnetic Rayleigh-Taylor instability (MRTI) governs material transport and mixing in astrophysical and laboratory plasmas under the influence of gravity and magnetic fields. While magnetic reconnection is known to occur during MRTI evolution, its role in the evolution and energy dynamics remains poorly understood. Here, we present a comprehensive analysis of the role of reconnection in the two-dimensional MRTI dynamics, using high-resolution simulations. We establish that reconnection, through facilitating plume merger, relieving magnetic tension, and enabling continued instability growth, forms an essential component for the long-term instability evolution. To quantify the role of reconnection in energy dynamics, we develop a robust automated reconnection detection algorithm and perform a statistical analysis across a range of magnetic field strengths. We find that reconnection accounts for up to $80\%$ of the magnetic-to-kinetic energy transfer in the weak magnetic field regime, while contributing minimally ($\approx 3\%$) to magnetic energy dissipation. Our results establish magnetic reconnection as a critical mechanism that regulates large-scale MRTI dynamics, with implications for astrophysical plasmas and turbulent mixing in magnetized flows.

Equations

$${\partial }_t\mathbf{u}-\nu {\nabla }^2\mathbf{u}=-(\mathbf{u}\cdot \nabla )\mathbf{u}$$

$${\partial }_t\mathbf{B}-\eta {\nabla }^2\mathbf{B}-c_p^2\nabla (\nabla \cdot \mathbf{B})$$

$${\partial }_t\rho -D{\nabla }^2\rho$$

$$\nabla \cdot \mathbf{u}$$

$$\begin{array}{r}\rho ={\rho }_l-\frac{\Delta \rho }{2}\left[\tanh \left(\frac{z-0.45L_z}{0.05}\right)+1\right]+\frac{\Delta \rho }{2}\left[\tanh \left(\frac{z}{0.05}\right)+1\right]\end{array}$$

$$u_z=\sum _{i=1}^{128}{a}_i\sin \left(\frac{2\pi k_{i}x}{L_x}+{\varphi }_i\right)e^{-(z^2/0.01)}.$$

$$B_c=\sqrt{\frac{({\rho }_h-{\rho }_l)g}{2k}}.$$

$$u_z=\sum _{k_i=8}^{100}{a}_i\sin \left(\frac{2\pi k_{i}x}{L_x}+{\varphi }_i\right)e^{-(z^2/0.01)},$$

$$\mathcal{A}={\int }_0^x{B}_z(X,0)\mathrm{d}X-{\int }_{-L_z/2}^Z{B}_x(x,Z)\mathrm{d}Z+C,$$

$$\text{TME}={\int }_0^{L_x}{\int }_{-L_z/2}^{L_z/2}\left(\frac{1}{2}{B}^2-\frac{1}{2}{B}_0^2\right)\mathrm{d}x\mathrm{d}z,$$

$$\text{TKE}={\int }_0^{L_x}{\int }_{-L_z/2}^{L_z/2}\left(\frac{1}{2}\rho u^2\right)\mathrm{d}x\mathrm{d}z.$$

$$D_{TME}={\int }_0^t{\int }_0^{L_x}{\int }_{-L_z/2}^{L_z/2}\eta ({\partial }_j{b}_i{)}^2\mathrm{d}x\mathrm{d}z\mathrm{d}t.$$

$$\mathcal{A}=\frac{B_0}{\kappa }\log [\cosh (\kappa y)+ϵ\cos (\kappa x)],$$

$$\frac{n_{\times }}{\sqrt{R{e}_m}}=\frac{0.15B_c}{B_0}.$$

$$RR=\sum _{i=1}^{n_{\times }}\eta |j_i|.$$

$$n_{\times }\equiv \sqrt{n_1{n}_2}\propto \sqrt{v_{A}h}\propto \sqrt{R{e}_m}.$$

$$E=\frac{\rho v^2}{\tau }=\rho v^3{k}_0.$$

$$k_{\eta }=E^{\alpha }{\rho }^{\beta }{\eta }^{\gamma }\equiv (\rho v^3{k}_0{)}^{\alpha }{\rho }^{\beta }{\eta }^{\gamma }$$

$$L^{-1}=(\rho L^2{\tau }^{-3}{)}^{\alpha }{\rho }^{\beta }(L^2{\tau }^{-1}{)}^{\gamma }$$

$$\alpha +\beta =0;2\alpha +2\gamma =-1;-3\alpha -\gamma =0.$$

$$k_{\eta }=E^{1/4}{\rho }^{-1/4}{\eta }^{-3/4}.$$

$$k_{\eta }=k_0^{1/4}{v}^{3/4}{\eta }^{-3/4}.$$

$$n\propto \frac{Lh}{k_{\eta }{k}_{\eta }}\propto L{h}^{3/4}.$$

Full text

I Introduction

The configuration of a high-density plasma being supported by a low-density plasma in the presence of gravitational and magnetic fields is ubiquitous in many laboratory [1] and astrophysical [2, 3] systems. Perturbations at the interface can lead to the interpenetration and mixing of the two plasmas, which is known as the magnetic Rayleigh-Taylor instability (MRTI). These structures of low-density plasma penetrating high-density plasma or vice versa are referred to as plumes. The region where the plumes penetrate or the two plasmas mix is called the mixing layer, which is typically a turbulent region [4]. The MRTI often plays a crucial role in the mixing or material transportation (through plumes) in numerous systems [2, 3]. For example, mixing of the corona and prominence materials in the solar prominence [5], transportation of the stellar material into the surrounding medium in supernova [6, 7], accretion of the disc material onto the central object in the accretion discs [8].

At the characteristic scales of the instability, plasma in most systems approximately obeys the ideal flux-freezing condition [9], wherein magnetic field lines are tightly coupled to plasma flow. Consequently, perturbations deform both flow streamlines and magnetic field lines, with plasma motion dragging magnetic fields as the instability evolves. When turbulent motions bring oppositely directed magnetic field lines into close proximity, intense current sheets form due to steep magnetic gradients. In these regions, the current becomes strong enough to trigger non-ideal magnetohydrodynamic (MHD) behavior, breaking the flux-freezing condition and altering magnetic topology through magnetic reconnection [10].

Magnetic reconnection within the MRTI mixing layer has been illustrated by numerical simulations [11, 12, 13] and observations [14]. Considering isolated reconnection events, these studies have shown that reconnection drives plasma downflows [11, 14], and leads to plume mergers [12].

Despite these advances, a key missing link in MRTI research is understanding the significance of magnetic reconnection in the large-scale evolution of the instability. Hydrodynamic Rayleigh-Taylor instability studies have established that plume mergers are essential to the instability’s nonlinear development [15, 16]. The role of magnetic reconnection in enabling plume mergers suggests that reconnection may be critical for the overall MRTI dynamics. This raises the question: Is magnetic reconnection an essential component of the long-term evolution of MRTI?

Beyond morphology, magnetic reconnection is also expected to influence energy dynamics. In reconnection regions, magnetic energy converts into kinetic and thermal energies via the Lorentz force and current-sheet dissipation, making reconnection an important source of local plasma heating and flow acceleration [17]. However, in MRTI, gravitational potential energy (GPE) drives the instability by injecting kinetic energy, displacing fluid vertically, and broadening the mixing layer. As the instability evolves, magnetic field lines deform and often bundle around plumes [12], forming dissipative current sheets. Thus, two energy sources coexist: GPE and magnetic reconnection. This motivates the question: Does magnetic reconnection significantly contribute to the energy dynamics of MRTI?

Previous studies have largely focused on isolated reconnection events [11, 12], but the turbulent nature of the MRTI mixing layer [18, 4, 19, 20] allows the magnetic field lines to contact at multiple locations, leading to simultaneous reconnection events. Such simultaneous multiple reconnection events have been observed in homogeneous turbulence systems [21, 22]. Unlike statistically steady homogeneous turbulence, the MRTI is dynamically evolving: the mixing layer grows approximately quadratically in time, and root-mean-square velocity increases roughly linearly [20]. Hence, the number and nature of reconnection events are expected to vary temporally.

Magnetic field strength is another critical factor influencing reconnection. Increasing the field strength suppresses small-scale perturbations and reduces turbulence intensity within the mixing layer [18, 4, 19, 20], thereby modifying the turbulent dynamics of MRTI. Accordingly, the characteristics of reconnection are expected to vary with magnetic field strength. However, a systematic statistical characterization of this dependence remains lacking. This motivates the question: How do reconnection statistics and dynamics evolve with time and magnetic field strength?

Addressing these questions, this paper aims to provide a comprehensive understanding of magnetic reconnection in MRTI. The qualitative role of reconnection in the large-scale instability evolution and energy dynamics is examined in Sections III and IV, respectively. Section V presents a statistical analysis of reconnection events as functions of time and magnetic field strength. The quantitative contribution of reconnection to the energy dynamics is assessed in Section VI.

II Numerical methodology

In the present study, MRTI is modelled numerically by superimposing a high density plasma $(\rho_{h})$ over a low density plasma $(\rho_{l})$ in the presence of an externally imposed uniform, unidirectional magnetic field $(\mathbf{B_{0}})$ and gravitational field $(\mathbf{g})$ as shown in figure 1(a). The magnetic and gravitational fields are along the horizontal and vertical directions (parallel and perpendicular to the interface), respectively. The numerical modelling was performed using Dedalus [23], an open-source, parallelized computational framework to solve the partial differential equations using the psuedo-spectral method. We solve an approximation of the variable-density non-ideal MHD governing equations [20]

[TABLE]

u, B, $\rho$ represent the velocity, magnetic field, and density, respectively. The bold letter is used to specify a vector quantity. Note that $\nabla\cdot u=0$ is the incompressible limit of the true velocity-divergence of the variable-density approximation [24]. The solenoidal condition for magnetic field ($\mathbf{\nabla}\cdot\mathbf{B}=0$) is ensured through divergence cleaning term ($c_{p}^{2}\mathbf{\nabla}(\mathbf{\nabla}\cdot\mathbf{B})$) [25]. In the present study, the value of $c_{p}$ is set such that the solenoid condition is satisfied to machine precision throughout the simulation. We consider two miscible plasmas with equal and constant density diffusion coefficient ($D$) [26]. The diffusion smooths sharp density gradients and aids mixing of plasma at the grid scale. $\nu,\eta$ are the coefficients of fluid viscosity and magnetic diffusion, respectively. $\nu,\eta,$ and $D$ are set to $10^{-4}$. We use Runge-Kutta scheme of fourth-order for time stepping, integrating the linear terms implicitly and the non-linear terms explicitly (see $\S$IX F of [23]). Dedalus does not use an explicit linear pressure solver. The value of pressure is determined by the divergence-free velocity constraint.

Here, the MRTI is studied in 2D with periodic boundary conditions in both directions. A domain of length $L_{x}{\times}L_{z}{=}4{\times}6$ units $(x:[0,L_{x}],z:[-L_{z}/2,L_{z}/2])$ with resolution of $2048{\times}3072$ is taken. The acceleration due to gravity is taken as 1. All the MRTI cases were run at a density ratio of 3. That is, all simulations have an Atwood number, defined as $\left(\frac{\rho_{h}-\rho_{l}}{\rho_{h}+\rho_{l}}\right)$, of 0.5.

The initial density profile is given by equation 2

[TABLE]

where $\Delta\rho{=}\rho_{h}{-}\rho_{l}$. The profile results in two interfaces, one at $z{=}0$ and the other at $z{=}0.45L_{z}(z=2.7)$ as shown in figure 1(a). The transition between the two densities is made continuous using a hyperbolic tangent profile with a half width of $l{=}0.05$ (see equation 2) spanning 25 grid points. This ensures that the interface is resolved sufficiently throughout the instability evolution despite the thinning of the interface. The above density profile is chosen to ensure periodicity in the $z$-direction.

The MRTI is excited by perturbing the lower interface (at $z=0$) with a vertical velocity of the form,

[TABLE]

The perturbations decay in a Gaussian profile about the interface. The amplitude $(a)$, wave mode $(k)$, and phase $(\phi)$ of the perturbation are chosen between [-0.0125, 0.0125], [1, 128], and [0, $\pi$], respectively. The perturbation has a wide range of wave modes (as suggested by the previous studies [27, 28, 29, 30]) along $x$.

To test the variation of reconnection dynamics with magnetic field strength, we run MRTI simulations over a range of magnetic field strengths between $1\%$ and $25\%B_{c}$. Here, $B_{c}$ refers to the magnetic field strength at which the wave modes parallel to the magnetic field are completely suppressed. $B_{c}$ is given by [31]

[TABLE]

Note that, $B_{c}$ is originally defined for a single-mode undular perturbation. For consistency, we choose the $B_{c}$ based on the smallest mode, $k/L_{x}=1$. We non-dimensionalise the imposed magnetic field strength with $B_{c}$, as a way to characterise the imposed magnetic field strength.

III Reconnection and evolution of MRTI

In $\S$I, we posed the question of whether reconnection is necessary for the long-term evolution of MRTI, analogous to the role of plume mergers in hydrodynamic RTI. To investigate this, we simulate two scenarios: an idealized MRTI where reconnection does not occur ($\S$III.1), and a non-ideal MRTI where reconnection occurs ($\S$III.2). Comparing these two cases, we infer the necessity of reconnection.

III.1 The pseudo-ideal-MRTI case

A truly ideal MRTI, where $D=\nu=\eta=0$, is numerically infeasible in Dedalus due to inherent stability constraints. Alternately, this can be simulated by using low $D,\nu,$ and $\eta$ such that only numerical diffusion drives the non-ideal physics, while ensuring $D,\nu,$ and $\eta$ are large enough to avoid numerical artifacts like ringing. However, this still does not prevent reconnection, the key physics of interest in the current study, since the reconnection occurs due to the numerical diffusion in the codes. In the context of the present study, understanding the role of reconnection, the key aspect is not the magnitude of the diffusion coefficients, but avoiding the reconnection. Therefore, we simulate an analogous ideal-MRTI case where reconnection is prevented by using a carefully designed perturbation (equation 5) and a strong magnetic field. The strong magnetic field ($B_{0}=15\%B_{c}$ in the present study) suppresses small wavelength perturbations [18] ($k>k_{c}\approx 44$ in the present study). As such, no reconnection occurs at the early stages of the instability. The reconnection at the later stages of instability is due to integral length scale plumes, i.e., large wavelength perturbations, which are avoided by limiting the smallest wave mode to 8 in our perturbation. Hence, we choose a perturbation of the form

[TABLE]

where $a_{i}\in[-0.05,0.05]$, $k_{i}\in[8,100]$, and $\phi_{i}\in[0,\pi]$. The values of $\nu,\eta,D,L_{x},L_{z},N_{x},N_{z}$ are same as mentioned in $\S$II.

The evolution of the ideal-analogue-MRTI (along with the magnetic field lines) is shown in figure 2. The magnetic field lines correspond to different values of magnetic potential ($\mathcal{A}$), given by $B(x,z)=\nabla\times\mathcal{A}$, a one-dimensional vector in a 2D problem. We calculate the $\mathcal{A}$ in the post-processing stage as

[TABLE]

where $C=\int_{0}^{x}B_{z}(X,0)\mathrm{d}X-\int_{-L_{z}/2}^{0}B_{x}(x,Z)\mathrm{d}Z$ is the gauge constant, chosen from a reference point ($(0,-L_{z}/2)$ in the current study) such that $\mathcal{A}=0$ at the interface ($z=0$) at time $t=0$. Due to the change in value of $B_{z}(0,-L_{z}/2)$, we see an error in the gauge constant by around $2\%$ at the end of the simulation, where the mixing layer reaches the bottom boundary of the domain. The magnetic potential throughout the study is calculated as described above. As the instability grows, the magnetic field lines elongate as shown qualitatively and quantitatively in figures 2, 3, respectively. The magnetic field line initially at $z=0$ deforms early in the instability, providing a representative measure of magnetic tension buildup. We track its temporal length, which correlates with the magnetic tension (figure 3, right axis). As the magnetic field lines elongate and saturate, the magnetic tension also increases and saturates correspondingly. This demonstrates that the length of the magnetic field line is a good tool to estimate the magnetic tension in the system. Note that the correlation between field line length and tension depends on field line location; central lines (e.g., at $z=0$) offer the most representative measure as they evolve from the time $t=0$.

During the instability growth, the turbulent magnetic energy (TME) and the turbulent kinetic energy (TKE) increase equally (see figure 4, $t\in[{\approx}1,{\approx}3.75]$). The TME and TKE are defined as

[TABLE]

[TABLE]

As the plumes evolve to integral length-scales, the energy required to bend the magnetic field lines increases, resulting in an unequal distribution of GPE between TKE and TME. A greater proportion of GPE is converted to magnetic energy. Consequently, the proportion of GPE converted to TKE reduces. This forms a positive feedback loop that leads to a continuous increase in TME (and hence a continued decrease in TKE), see figure 4, $t\gtrapprox 3.75$.

Over time, the magnetic tension becomes strong enough to balance the magnetic buoyancy force and arrest the further growth of instability, as shown in figures 2(c) and (d) where the mixing layer height is approximately the same at $t\approx 7,8$. At this stage, the plumes push against each other and oscillate left and right (figures 2(c), (d)).

To summarize, we showed that the length of the magnetic field line represents the magnetic tension in the system. As the instability evolves, the length of the magnetic field lines and with that the magnetic tension increases. If the tension is not relieved, over time the magnetic tension force becomes strong enough to counter the gravity and arrest further growth of instability. Thus the magnetic fields saturate the instability growth over time.

III.2 The non-ideal-MRTI

To understand how reconnection influences the evolution, we will now consider the case of non-ideal-MRTI. As mentioned before, we change the perturbation to replicate the physics of ideal and non-ideal MRTI physics. Hence, we keep all the parameters the same between $\S$III.1 and $\S$III.2 except for the velocity perturbation. The numerical methodology of the non-ideal-MRTI simulation, initial perturbation, $\nu,\eta,D,L_{x},L_{z},N_{x},N_{z}$ is as explained in $\S$II. For comparison with the ideal-MRTI, we use the same field strength as the ideal-MRTI, $B_{0}{=}15\%B_{c}$. Due to the small $\eta(=10^{-4})$, the field lines evolve approximately with the flow.

Similar to the ideal-MRTI case, the strong magnetic field suppresses most small-scale wave modes, and the mixing layer evolves as laminar plumes. However, the non-linear interaction of large-scale plumes results in the contact of magnetic field lines around them. Due to the non-zero magnetic diffusivity, the magnetic field lines reconnect. Such a reconnection is shown in figures 5(a), (b) using a sample field line ($\mathcal{A}=0.01$). The magnetic reconnection leads to shorter magnetic field lines. The reduction in the length of the field line at $t\approx 5.8$ (when the reconnection occurs) is shown in the inset figure of 6. The shortening of field line length meant that the magnetic tension is relieved (c.f. $\S$III.1). Thus, the magnetic reconnection prevents the domination of TME over TKE. This is illustrated in figure 7 where the system is maintained at approximate energy equipartition throughout the simulation (unlike the ideal-MRTI).

The domination of non-ideal processes at these reconnecting locations results in the transportation of plasma across the reconnecting field lines, allowing the two plumes to coalesce (see figure 5(c)) towards forming a larger plume. The merger process over time is shown in figure 5.

Now, consider the forces acting on the plumes. The plume is under the influence of two forces: $i)$ buoyancy force ($F_{b}=(\rho_{h}-\rho_{l})V_{plume}g$), which tries to displace the plume upward, $ii)$ magnetic tension, which opposes the buoyancy force. The length of the field line $\mathcal{A}=0.01$ decreases from $\approx 20$ to $\approx 17$ due to reconnection at $t\approx 5.6$ (see figure 6). Similarly, the reconnection of field line $\mathcal{A}=0.0$ at $t\approx 6.6$, also shortens its length. On the other hand, the buoyancy force $F_{b}$ is increased due to a larger plume volume. Arresting the growth of the merged plume demands a larger magnetic tension. However, the magnetic tension reduces due to reconnection, making it insufficient to counter the increased buoyancy force due to plume merger. Thus, there is a magnetic tension deficit that leads to further growth of the plume and the instability.

III.3 Section summary

In summary, the ideal MRTI evolves until increasing magnetic tension halts further growth. In contrast, in non-ideal MRTI, reconnection relieves magnetic tension by changing magnetic field topology and reducing the field line length. Further, reconnection allows the merger of the interacting plumes, forming a bigger plume that experiences a greater buoyancy force. The merged plume requires higher magnetic tension to be arrested, but this tension is diminished by reconnection. Thus, magnetic reconnection enables sustained instability growth, demonstrating its critical role in MRTI evolution.

IV Energy dynamics due to reconnection: Qualitative

In $\S$I, we discussed how magnetic reconnection contributes to the energy dynamics of MRTI. To evaluate its significance, we first need to understand how reconnection energy dynamics manifest in the system. In this section, we qualitatively illustrate the conversion of magnetic energy into kinetic and thermal energy during reconnection events. The quantitative evaluation of these processes, which requires a statistical approach, is deferred to $\S$VI.

IV.1 Energy interplay between kinetic and magnetic energy

To understand the interplay between kinetic and magnetic energy, we consider the term $\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}$, which arises in the turbulent kinetic energy (TKE) equation. The term is a product of flow velocity ($\mathbf{u}$) and the Lorentz term ($(\mathbf{B}\cdot\nabla)\mathbf{B}$. Here, $\mathbf{B}$ denotes the total magnetic field, including both the background and perturbations. This term quantifies the work done by magnetic tension on the fluid: it is positive when magnetic energy is transferred to the flow, and negative when kinetic energy is expended against magnetic tension.

Buoyancy forces drive the plumes to rise or fall, deforming the magnetic field and generating magnetic tension that resists this motion. A schematic of the opposing flow and magnetic tension is provided in figure 8 (bottom inset). Due to the opposing nature of flow and magnetic tension, the term $\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}$ would be negative (see the plume heads in figure 8). However, some plume heads exhibit positive values of $\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}$ (one such is marked with a red arrow), indicating that magnetic tension and plume motion are aligned. This alignment likely results from interactions with neighboring plumes. For instance, the larger plume (blue arrow) exerts a downward force on the adjacent smaller plume (red arrow). The figure 5(a) denotes the approximate time instant as the figure 8. From figure 5, we find that the small plume which is initially at $(x,z)\approx(2.1,0.25)$ has a downward motion due to the adjacent large plume. The downward motion of the small plume due to the large plume can be seen from figure 5. Thus, the direction of magnetic tension and the net flow of the smaller plume are the same, making the quantity $\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}$ positive. A schematic of the flow and magnetic tension is shown in figure 8 as an inset at the top.

To demonstrate the influence of reconnection on the energy dynamics, we zoom into a chosen reconnection point (see figure 8(right)). Near the reconnection site, we observe outflows characterized by positive $\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}$, indicating localized conversion of magnetic energy into kinetic energy. The influence of reconnection can be seen to extend far beyond the reconnection point. The injection of kinetic energy into the system improves the fluid mixing locally.

For such distinct dynamics to appear, the reconnecting current sheets should be strong, and the turbulence surrounding the reconnection point should be minimal. But the majority of the reconnection events occur in turbulent regions and due to weak current sheets (evidenced in $\S$V). Hence, such distinct dynamics are uncommon.

IV.2 Energy dissipation:

In addition to kinetic energy conversion, reconnection also leads to the dissipation of magnetic energy. Consequently, reconnection regions are also the sites of significant dissipation. In a homogeneous turbulent system, the reconnection sites are major contributors to TME dissipation [21, 22]. However, in MRTI, current sheets also arise from bending of magnetic field line around plume heads. This makes the relative contribution of reconnection to global magnetic energy dissipation less clear.

To understand this, we plot the spatial contours of instantaneous magnetic energy dissipation (figure 9). The magnetic energy dissipation is calculated as

[TABLE]

Upon zooming into the neighborhood of a reconnection point (marked by {\color[rgb]{0,0.88,0}\times}), we observe locally intense dissipation (see figure 9(right)). However, significant dissipation is also observed in non-reconnecting current sheets, particularly near plume heads and flanks, where magnetic fields are highly concentrated due to the opposing nature of flow and magnetic tension. The concentration of these field lines leads to strong current sheets that dissipate magnetic energy significantly.

Since reconnection sites are relatively small in spatial extent compared to other current sheets, their contribution to global magnetic energy dissipation may be limited. A quantitative analysis requires detecting the reconnection regions, which will be discussed in $\S$V.1. Following this, a quantitative analysis of energy dissipation due to reconnection will be made in $\S$VI.1.

IV.3 Section summary

In summary, we qualitatively examined how magnetic reconnection contributes to energy dynamics in MRTI. We illustrated the conversion of magnetic energy into kinetic energy via reconnection-driven outflows, whose influence extends well beyond the reconnection site. However, such clear signatures are rare and typically suppressed by surrounding turbulence. Additionally, we highlighted that while reconnection does lead to magnetic energy dissipation, significant dissipation also occurs in non-reconnecting current sheets due to the natural flow–field interactions in MRTI. A more quantitative treatment follows in Sections V and VI.

V Statistics of MRTI reconnection

To this point, the role of reconnection was elucidated using a sample reconnection point. However, MRTI has simultaneous multiple reconnection events, which change quantitatively with time and magnetic field strength (see $\S$I). Towards a statistical analysis of the reconnection, we propose a reconnection detection algorithm in $\S$V.1. Using the algorithm, we explore the time evolution of the number of reconnection points in $\S$V.2, and characterize their strength in $\S$V.3. In particular, we focus on the fraction of strong reconnection current sheets, since they are expected to have a major influence on the energy dynamics of the system.

V.1 Detection of reconnection points

The first step towards quantifying the reconnection points is to detect them. Detecting the reconnection points could be challenging, particularly in the weak magnetic field case where the mixing layer is turbulent, resulting in numerous reconnection events [21]. Hence, an automated reconnection detection technique is necessary. While a reconnection detection algorithm already exists [21] (see $\S$V.1.1), it is limited to 2D or quasi-2D (3D system with a strong guide field) systems. Hence, we develop a new reconnection detection algorithm that can be extended to 3D, see $\S$V.1.2. However, in this paper, we present the preliminary case of detecting reconnection points for a 2D system.

V.1.1 Previous reconnection detection techniques

In 2D, reconnection typically occurs at magnetic null points, which can exhibit either saddle (X-point) or center (O-point) topologies [32]. Servidio et al. [21] proposed a detection algorithm based on this, originally for homogeneous 2D turbulence. Their method first identifies candidate points where the magnitude of the magnetic field strength is near zero. At these points, the eigenvalues of the Hessian matrix $\left(\partial^{2}\star/\partial x_{i}\partial x_{j}\right)$ of magnetic vector potential $\mathcal{A}$ are calculated to gain insight into the magnetic field topology. The points with positive and negative eigenvalues have saddle topology, i.e., X-points. The above method was also used for the quasi-2D system [33], by detecting the reconnection points at each 2D plane along the guide field direction.

However, this approach is unsuitable for fully 3D configurations, where reconnection does not always occur at null points. Even in 2D, since eigenvalues were used to classify the X- and O-points, we must be as close to the real X or O-point as possible. Hence, this method demands an extremely high resolution or an accurate interpolation technique to map the data to finer scales to locate the X-points.

V.1.2 Present reconnection detection technique

To overcome the limitations of 2D-based methods, we develop a generalized algorithm that detects reconnection points based on local magnetic field configurations. For simplicity, we will demonstrate the algorithm using the Fadeev equilibrium configuration. The Fadeev equilibrium [34, 35] is a continuous arrangement of magnetic islands (plasmoids) with consecutive saddle and center topology leading to the X- and O-points (represented as $\times$ and $\blacksquare$ respectively in figure 10). The Fadeev equilibrium is modelled as

[TABLE]

where $\mathcal{A}$ is the magnetic vector potential, $B_{0}$ is the initial magnetic field strength, and $\kappa$ represents the wave number. $\kappa$ determines the number of plasmoids in the configuration and the number of null points. For $\kappa{=}n\pi,$ we have $2n$ X- and O-points. $\epsilon$ controls the eccentricity of the plasmoid. The limits $\epsilon\rightarrow 0$ and $\epsilon\rightarrow 1$ represent the highly eccentric plasmoid and perfectly circular plasmoid, respectively. The limit $\epsilon=0$ gives the Harris equilibrium, where the magnetic field lines are parallel and have no X- and O-points.

From figure 10(top), we see that the horizontal ($B_{x}$) and vertical ($B_{y}$) components of the magnetic field change sign across the reconnection point along the vertical ($y$) and horizontal ($x$) directions, respectively. To detect such points, we consider a 2D five-point stencil of radius $\Delta x$ around the point $P(i,j)$ (see the solid lines in figure 11), where $\Delta x$ is the grid resolution. We check if $B_{y}$ is changing along $x$ across the point $P$ using condition $B_{y}(i-1,j){\times}B_{y}(i+1,j){<}0$. Similarly, we verify if $B_{x}$ is changing along $y$ using $B_{x}(i,j-1){\times}B_{x}(i,j+1){<}0$. This is schematically shown in figure 11 (a),(b). To avoid the null points due to noise, we discard points of order less than $-10^{-8}$. By filtering the points that satisfy both conditions, we obtain the center and saddle-type null points.

To extract the reconnection points alone, we use the cross gradients of magnetic field components $(\partial_{x}B_{y},\partial_{y}B_{x})$, where $\partial_{\star}\Phi=\partial\Phi/\partial\star$. The product of these cross gradients is positive for X-points and negative for O-type points. The points filtered above that satisfy the condition $\partial_{x}B_{y}{\times}\partial_{y}B_{x}>0$ are taken as reconnection points. The magnetic field configuration at these reconnection points are further verified by repeating the above process with a five-point stencil of radius $2\Delta x$ and $3\Delta x$ (see the dashed lines in figure 11). The points that show saddle nature at all these radii are only considered as X-points.

The above algorithm is suitable when the reconnecting field lines are parallel to the coordinate axes, as shown in figure 10(top). However, due to the turbulent nature of the mixing layer, often the reconnecting field lines are oriented at an angle to the reference frame (see figure 10(bottom)). In such cases, the 2D five-point stencil must first be aligned with the field lines before the above detection algorithm is implemented. Towards this, we first detect the orientation of current at the point ($\mathbf{--}$ in figure 10(bottom)) to the reference coordinate axis (-. line in figure 10(bottom)). Since the current ($J$) is scalar quantity in 2D, we calculate the orientation of the current sheet using $\theta=\tan^{-1}\left(\partial_{x}J/\partial_{z}J\right)$. A $3{\times}3$ box is considered (figure 12(i)) and the box is then oriented along the current sheet $\theta$ as shown in figure 12(ii). The magnetic field components parallel and perpendicular to the current sheet (represented as $B_{\parallel}$ and $B_{\perp}$) are calculated (see figure 12(iii)). We now check if the $B_{\parallel},B_{\perp}$ are changing sign across the point of investigation using $B_{\perp}(i-1,j){\times}B_{\perp}(i+1,j)<0$, $B_{\parallel}(i,j-1){\times}B_{\parallel}(i,j+1)<0$ (figure 11(c)). We then filter the X-point based on the condition $\partial_{\parallel}B_{\perp}{\times}\partial_{\perp}B_{\parallel}>0$. The issues of this method and their solutions were presented in $\S$A. The reconnection detection algorithm was developed as a post-processing code. The current algorithm works only for the 2D case.

V.2 Reconnection events in MRTI

The developed algorithm is used to detect the reconnection points in the MRTI mixing layer for different magnetic field strength cases. From figure 13(left), we see that, for all magnetic field strengths, the number of reconnection points ($n_{\times}$) increases with time. The reason for the increase in $n_{\times}$ is attributed to two factors — $i)$ height of the mixing layer and $ii)$ the turbulent motion of magnetic field lines. As the mixing layer widens and the turbulence increases, we see that the number of reconnection points increases. In the current work, we find that the number of reconnection points scales as the square root of the magnetic Reynolds number, as shown in figure 13(right). Potential reasons for this scaling law are discussed in $\S$VII.1. From figure 13(right) we see that for a given $Re_{m}$, $n_{\times}$ decreases with increasing magnetic field strength. For the current set of simulations, we find that

[TABLE]

To explain the above behaviour, we plot the instantaneous density contours of the mixing layers for different magnetic field strengths along with the detected reconnection points and magnetic field lines (see figure 14). The density contours show that the mixing layer is more turbulent at weaker field strengths. This is due to the selective suppression of small-scale perturbations ($k{\geq}k_{c}$) by the magnetic field in the linear regime [20]. The critical wave number $k_{c}$ is the smallest wave mode that can be suppressed by the parallel magnetic field in the linear regime [31]. From the linear growth rate equation, $k_{c}$ decreases with increasing magnetic field strength [31]. Thus, in the weak magnetic field case, a wide range of waves ($k{\in}[1,k_{c}]$) grow, and the mixing layer is turbulent [4, 19, 20]. The turbulent nature of the mixing leads to frequent interaction of field lines, resulting in an increased number of reconnection points in the weak magnetic field case (see figure 14(left)). Conversely, $k_{c}$ is small in the strong magnetic field case, and only a limited small $k$ modes grow. Therefore, the mixing layer has laminar-like plumes (figure 14(right)). Due to less turbulence, the reconnection points are also low.

V.3 Strength of the reconnection events

In the previous subsection, we quantified the reconnection events. Here, we will delve into their characteristics, particularly the strength of the reconnecting current sheets. This provides a quantitative insight into the influence of reconnection on energy dynamics, like energy conversion to kinetic energy and energy dissipation. Further, the variation of these characteristics with time and imposed magnetic field strength is discussed.

The first step is to segregate the current sheets based on their strength. However, due to the increase in turbulence and magnetic field amplification, the nature of the mixing layer and the characteristics of current sheets change with time. Similarly, due to suppression of perturbations, the nature of the mixing layer (and in turn the characteristics of the current sheets) changes with magnetic field strength. Therefore, a direct comparison of current sheets across time and magnetic field strengths is inappropriate. To make the reconnecting current sheets comparable across different time instants and magnetic field strengths, we need to normalize the current sheet by the instantaneous root mean square current $\left(j_{rms}=\int_{V}|j|\mathrm{d}V\right)$, where $V$ is the volume of mixing layer. As the instability evolves, the height of the mixing layer and the strength of current sheets increase. Capturing these, $j_{rms}$ also increases (see the legend of figure 15), making it an appropriate normalization parameter.

The normalized reconnecting current sheets are segregated into seven bins: current sheets that are weaker than $j_{rms}$ and then on in powers of 2 ($2^{n}$ to $2^{n+1}$, $n\in[0,5]$). Figure 15 shows the histogram of reconnecting current sheets in each range at various time instants. The histogram is plotted for the field strength $B_{0}=1\%B_{c}$. The proportion of strong current sheets (hereon, defined as the current sheets with strength greater than $8j_{rms}$) is observed to be less than $\approx 7\%$ of the total reconnecting current sheets at any well-developed non-linear state ($t>4$). At $t=4,6,8$, the percentage of strong reconnecting current sheets is approximately $6\%,$ $5\%,$ and $7\%$, respectively. The choice of the threshold $>8j_{rms}$ for considering as strong current sheet comes from previous studies where they considered $5j_{rms}$ [33], $10j_{rms}$ [36]. However, the choice is arbitrary, and quantitative results vary accordingly, while the broad conclusions remain the same.

The characteristics of the reconnecting current sheets also vary with the magnetic field strength. Figure 16 shows the distribution of current sheets for MRTI of different magnetic field strengths at a well-evolved non-linear state, $t=8$, see figure 14. The proportion of strong reconnecting sheets is increasing with field strength. The percentage of strong reconnecting sheets is approximately $7\%,$ $8\%,$ $11.5\%,$ $16\%,$ $20.5\%,$ $19.5\%,$ and $50\%$ for $B_{0}=1,$ $2,$ $3,$ $5,$ $7,$ $10,$ and $15$ $(\%B_{c})$ respectively. This trend suggests that while reconnection is more frequent in weak-field, turbulent flows, it is more intense in strong-field, laminarized regimes, reflecting the concentration of magnetic energy into fewer, more coherent structures.

V.4 Section summary

To summarize $\S$V, an algorithm that detects the reconnection points based on magnetic field configuration was developed. Using the algorithm, we quantified the number of reconnection points in the mixing layer as the instability evolves for different magnetic field strengths. For all field strengths, the number of reconnection points ($n_{\times}$) increases with time in a self-similar fashion, following $n_{\times}{\propto}t^{3/2}$. The proportionality constant decreases with increasing magnetic field strength. The large number of reconnection points in the weak field case could be attributed to the higher turbulence in the system. Analysis of current sheet strength revealed that most reconnection events are relatively weak ($<8j_{rms}$), but the fraction of strong events increases with magnetic field strength. The study highlights a key contrast between the reconnection in weak and strong magnetic field cases: while the former is spatially distributed and frequent, the latter is more localized but energetically significant.

VI Energy dynamics due to reconnection: Quantitative

In $\S$III.2 we illustrated that magnetic reconnection controls the TME by relieving magnetic tension. In $\S$IV, we illustrated that the reconnection converts magnetic energy partly to kinetic energy, dissipating the rest. However, a quantitative analysis was deferred due to the lack of precise reconnection point data. Having developed a reconnection detection algorithm (see $\S$V.1), we are equipped to quantify the role of reconnection in the system’s energy dynamics, which is the focus of this section.

Before we delve into the energy dynamics due to reconnection, let us understand the energy flow in the MRTI system. Energy enters the mixing layer as fluid extracts gravitational potential energy (GPE), which increases the total kinetic energy (TKE). This energy extraction is estimated by the quantity $u_{z}\delta\rho g$. A portion of this kinetic energy is transferred into the magnetic field (see blue regions in figure 8), increasing the TME. Some of this magnetic energy is then returned to the fluid (see red regions in figure 8). This exchange between magnetic and kinetic energy is described by the term $\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}$, as discussed in $\S$IV.1. To isolate directional energy transfer, we define $[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{-}$ for flow-to-field energy transfer, and $[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}$ for field-to-flow transfer. We analyze these terms by computing volume-averaged values of $u_{z}\delta\rho g$, $[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{-}$, and $[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}$ over time, as shown in figure 17. The results show that most GPE extraction (panel a) is transferred into the magnetic field (panel b), and about half of this energy is subsequently returned to the fluid (panel c). Consequently, the TKE is modulated by $[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}$, while the net TME is governed by the difference $\left([\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{-}-[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}\right)$.

VI.1 Role of reconnection in magnetic energy dissipation

Reconnection contributes to the energy budget not only by transferring magnetic energy to the fluid but also by dissipating magnetic energy. This dual role was illustrated qualitatively in $\S$IV. Now, having identified reconnection sites, we quantitatively assess their contribution. To estimate dissipation, we focus on the reconnecting current sheets. Below is a step by step procedure of how we estimate energy dissipation due to reconnection:

Define a region of size $60\Delta x\times 60\Delta z$ centered about reconnection point 111Larger regions yielded similar results.. 2. 2.

Rotate the region about the reconnection point (as shown in figure 12) such that the orientation is aligned with current sheet. Orientation of current sheet ($\theta$) is computed as described in $\S$V.1. 3. 3.

Compute the root-mean-square current $\hat{j}_{rms}$ in the box. Note that the $\hat{j}_{rms}$ specified here is different from $j_{rms}$, as the former represents the RMS current calculated from the local rotated box of $60\Delta x\times 60\Delta z$, while the latter is calculated from the whole domain. 4. 4.

Starting from the reconnection point, we trace along and across the sheet until the current drops to $2\hat{j}_{rms}$, defining the sheet’s length and width 222The choice of $2\hat{j}_{rms}$ comes from trial and error of various sampled reconnection points. We found that the threshold of $\hat{j}_{rms}$ includes other regions of the patch beyond the reconnecting current sheet, and threshold larger than $2\hat{j}_{rms}$ are increasingly more stringent, neglecting the outer regions of reconnecting current sheets. Because $\hat{j}_{rms}$ is locally defined, this approach is robust for weak and strong reconnection alike.

The above process is repeated for all reconnecting current sheets. The total energy dissipation is the volumetric average of energy dissipation from all reconnection events.

The TME dissipation is computed as $\eta(\partial_{j}b_{i})^{2}$ for each reconnecting current sheet based on the current sheet’s dimensions. Integrating this dissipation in space and time over all reconnection events yields the total magnetic energy loss attributable to reconnection. We find that reconnection accounts for approximately $3\%$ of the total TME dissipation (see figure 18). Thus, while reconnection clearly contributes to dissipation, its overall role is relatively minor. This conclusion holds when current sheet dimensions are alternatively estimated using the Full Width at Half Maximum (FWHM) method, confirming the robustness of the result.

VI.2 Role of reconnection in the energy transfer from magnetic field to fluid

Estimating the energy transferred from the magnetic field to the fluid is less straightforward, as it depends on the dynamics of reconnection outflows, which spread spatially, and the extent of spatial influence depends on the strength of the current sheet and local turbulence. To obtain a practical estimate, we define a box of size $100\Delta x\times 100\Delta z$ centered at each reconnection point, shown by dashed lines in figure 20 (box size effects are discussed in $\S$VI.2.1). Within each box, we compute $[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}$. Repeating this for all reconnection sites and averaging over volume gives the total contribution of reconnection to field-to-flow energy transfer. This is shown in figure 19. In contrast to energy dissipation, reconnection plays a significant role in energy transfer. At low magnetic field strengths, as much as $80\%$ of the energy transfer from the magnetic field to the fluid occurs near reconnection points, indicating that reconnection plays a dominant role in this regime. This correlates with the higher number of reconnection events observed in weak-field cases. As the magnetic field strength increases, the fraction of energy transfer attributable to reconnection decreases.

VI.2.1 Reconnection box size effect on the dynamics

To quantitatively estimate reconnection’s influence on energy dynamics, we considered a box of $100\Delta x\times 100\Delta z$ in $\S$VI.2. This is to ensure that the region of reconnection influence is well covered. Through figure 20 we evidence that the considered box size encompasses the reconnection influenced region. A larger box will include the dynamics beyond the reconnection, and a smaller box does not capture the dynamics from reconnection fully. From figure 20, the outflow jets are hardly captured from $50\Delta x\times 50\Delta z$ box. Note that the same reconnection point chosen in figure 8(right) is considered here.

Evidently, the choice of box size does influence the magnitude of $\left(\frac{1}{V}\int_{V}[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}dV\right)_{rec}$. But, we acknowledge that considering a single box size may not be applicable for all reconnection points, given the wide range of characteristics of reconnecting current sheets (cf. figures 16, 15). But, choosing a different box size for each reconnection depending on current sheet strength is tedious. Further, a appropriate relation between the current sheet strength and box size is unavailable. To quantitatively estimate the influence of box, we calculate $\left(\frac{1}{V}\int_{V}[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}dV\right)_{rec}$ at a particular time instant (t = 7) for different box sizes (X-axis) and magnetic field strengths. Normalizing with total $\frac{1}{V}\int_{V}[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}dV$ at the considered time instant (t = 7) and magnetic field strengths, we determine the contribution of reconnection to total magnetic to kinetic energy conversion. From figure 21, we see that the box size has a linear relationship with $\frac{\left(\frac{1}{V}\int_{V}[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}dV\right)_{rec}}{\frac{1}{V}\int_{V}[\mathbf{u}\cdot(\mathbf{B}\cdot\nabla)\mathbf{B}]_{+}dV}$.

Note that we use different box sizes to capture energy dissipation ($60\Delta x\times 60\Delta z$) and energy outflow ($100\Delta x\times 100\Delta z$), since the energy dissipation is more dominant at the reconnection point, as shown in figure 9(right), while outflow jet dynamics are more dominant away from reconnection point, as shown in figure 8(right).

VI.3 Reconnection rate

Another characteristic of magnetic reconnection is the rate at which the field lines move through the reconnection point, i.e., how fast the reconnection occurs, called the reconnection rate. For a steady-state 2D reconnecting system, the reconnection rate is calculated based on the electric field ($E=u{\times}B+\eta j$) at the reconnection point. Since the reconnection point in 2D is a magnetic null point ($|B|\rightarrow 0$), the electric field reduces to $\approx\eta j$.

In our system, reconnection occurs at many locations simultaneously, each with different current sheet properties. Consequently, reconnection proceeds at a range of rates—i.e., the system exhibits a distribution of reconnection rates rather than a single characteristic value. These distributions resemble those shown in figures 15 and 16, depending on whether we vary time at fixed magnetic field strength or vice versa, respectively.

Nevertheless, to understand how reconnection rate scales with global parameters such as the magnetic Reynolds number, it is useful to define an aggregate quantity. We compute the total reconnection rate $RR$ at each time instant by summing the reconnection contributions from all $n_{\times}$ detected reconnection events:

[TABLE]

Figure 22 shows the time evolution of $RR$ normalized by the total current dissipation in the domain, i.e., $\int_{V}\eta|j|\mathrm{d}V$, for various magnetic field strengths. Interestingly, the normalized $RR$ remains approximately constant in time, even though the magnetic Reynolds number that drives the field lines to reconnect increases with time. This temporal constancy appears across all magnetic field strengths.

We also observe that $RR$ increases systematically with magnetic field strength. In weak-field regimes, reconnection predominantly involves weak current sheets (cf. figure 16), leading to smaller reconnection rate. As the field strengthens, more intense sheets form, and reconnection becomes correspondingly more vigorous.

VI.4 Section summary

To summarise $\S$VI, we found that the energy extracted from gravitational potential energy (GPE) is largely converted into magnetic energy. The magnetic field, in turn, returns nearly half of this energy to the system as kinetic energy. Reconnection plays a central role in this energy transfer. In weak magnetic field regimes, reconnection accounts for approximately $80\%$ of the total magnetic-to-kinetic energy conversion. Although this fraction decreases with increasing field strength, reconnection remains a significant mediator of energy transfer across all regimes. In contrast, reconnection contributes only marginally to total magnetic energy dissipation ($D_{\text{TME}}$). This limited contribution arises in part because most reconnection events involve weak current sheets, where dissipation is minimal (cf. $D_{\text{TME}}=\eta j^{2}$). Moreover, the system preferentially channels magnetic energy into flow, rather than dissipating. Finally, we found that the magnetic field strength also affects the reconnection rate. As the field strength increases, the total reconnection rate grows. However, for any given field strength, the total reconnected flux summed over all reconnection points remains approximately constant in time.

VII Discussion

VII.1 Variation of number of reconnection points with Reynolds number

In $\S$V.2, we found that the number of reconnection points varies as the square root of the magnetic Reynolds number $Re_{m}$. Here we discuss potential explanations for this scaling. As mentioned in $\S$V.2, the number of reconnection points depends on the size of the mixing layer and the level of turbulence in the mixing layer. For a given turbulence intensity, as the mixing layer widens, the number of locations at which the field lines contact increases. Thus, the number of reconnection events increases with the mixing layer height ($h$), i.e., $n_{1}\propto h$. On the other hand, for a fixed mixing layer height, as the turbulence in the mixing layer increases, the magnetic field lines contact frequently, thereby increasing the number of reconnection events. Thus, we expect the number of reconnection points to increase with the turbulence. The turbulence in the system is often quantified by the turbulent Alfven velocity, defined as $v_{A}{=}\sqrt{\int_{V}\frac{(b_{x}^{2}+b_{z}^{2})}{\rho}dV}$. Thus, $n_{2}\propto v_{A}$. Here, $n_{1}$ and $n_{2}$ are the number of reconnection points at a fixed turbulence intensity and height, respectively.

To understand our assumption of linear relation between $n_{1}$ and $h$, consider a domain of size $L\times h$ that is fully turbulent system with turbulence intensity $v_{A}$. Let $n_{1}$ be the number of reconnection points in this system. Now consider another turbulent domain with the same domain size and turbulence intensity with $n_{1}$ reconnection points. Stacking these two systems and ensuring the total turbulence intensity is same as individual domain, we have a total domain of $L\times 2h$, and the total number of reconnection points would be $2n_{1}$. Thus, $n_{1}$ is linearly proportional to $h$. To understand the linear relation between $n_{2}$ and $v_{A}$, consider another case where we have a fully turbulent system in the domain $L\times h$. As the turbulence intensity in the domain increases from $v_{A}$ to $2v_{A}$, the frequency of magnetic field line interaction becomes twice, and hence the number of reconnection points also increases by 2-fold. Thus, the number of reconnection points is linearly proportional to $v_{A}$.

In the MRTI, as the instability evolves, the mixing layer height and the turbulence increase simultaneously, which is reflected in the increase of magnetic Reynolds number ($Re_{m}{=}v_{A}h/\eta$). The temporal increase of $Re_{m}$ can be seen from figure 13(center). The number of reconnection points in the MRTI mixing layer is given by the geometric mean of the $n_{1},n_{2}$, which gives,

[TABLE]

The latter comes from constant $\eta$. Thus, the number of reconnection points ($n_{\times}$) is expected to scale with magnetic Reynolds number as $\sqrt{Re_{m}}$. This was found to be consistent from our numerical simulations, as shown in figure 13(right).

Kalluri and Hillier [20] showed that $\int_{V}\frac{1}{2}b^{2}\mathrm{d}V\equiv L_{x}h(\frac{1}{2}b^{2})\propto h^{2}$, i.e., $v_{A}{\propto}\sqrt{h}$. Substituting this relation in equation 13, we get $n_{\times}{\propto}h^{3/4}$. From the quadratic growth of mixing layer height [20], $n_{\times}\propto t^{3/2}$. These scaling laws were found to be satisfied for all magnetic field strengths in our simulations. The $t^{3/2}$ growth demonstrates that the number of reconnection events grows in a self-similar fashion in the non-linear MRTI. This signifies that the self-similarity of MRTI is not restricted to flow variables like $\rho u,$ $b,$ TKE, TME, [20] but extends to other parameters like $n_{\times}$, that depend on the flow variables.

Alternatively, a more theoretical proof for $h^{3/4}$ scaling can be obtained from the energy dynamics analysis. The rate of energy injection at the integral length scales ($k_{0}$) is given by

[TABLE]

The latter comes from the scaling $\tau=\frac{1}{k_{0}v}$. From Kolmogorov’s theory, the injected energy trickles down to dissipative length scales. The scale at which the energy dissipates depends on the energy injected and the (magnetic) diffusion coefficient, in the case of reconnection current sheets.

The scale at which energy dissipates predominantly can be written as

[TABLE]

From the dimensional analysis,

[TABLE]

[TABLE]

Solving for $\alpha$, $\beta$, and $\gamma$ we get $\alpha=1/4$, $\beta=-1/4$ and $\gamma=-3/4$. Thus

[TABLE]

Substituting the energy injected rate in the above equation, we get

[TABLE]

To write $n$ in terms of $h$, we write different terms of $k_{\eta}$ in terms of $h$. The energy injection occurs at integral length scales; therefore, $k_{0}$ is of order, $1/h$. From the quadratic growth of mixing layer height with time, $v\propto\sqrt{h}$. With these substitutions, we get $k_{\eta}=h^{1/8}$. Our numerical simulations indeed show that the characteristic length scale of the current sheets is of the order $h^{1/8}$. In figure 23, we evidence this for three magnetic field strength cases.

The number of current sheets of scale $k_{\eta}$ in the mixing layer region ($L\times h$) is given by

[TABLE]

VII.2 Dynamics of reconnection jets

The outbursts from reconnection are often interesting from a turbulence perspective, particularly in the strong magnetic field case, where these are more prominent. We will discuss the dynamics of these structures. From the fluid mixing perspective, we find that the reconnection jets lead to fluid mixing by transporting the diffused matter at the reconnection region into one of the fluids, as shown in the figure 24 where diffuse matter is transported to a lower density fluid. However, a quantitative estimate of the mixing solely due to reconnection is difficult. This is because, in the weak magnetic field case, most fluid is mixed by the turbulence mixing of fluids, making the mixing due to reconnection indistinguishable. In the strong magnetic field case, the intermediate density fluid is predominantly formed due to mass diffusion, and the role of reconnection in fluid mixing seems relatively less prominent.

The other aspect is the dynamics of these out-flowing jets themselves. While the outflow jets are known to exhibit shear instabilities [39], we do not see such secondary KHI in our current simulations. This could be due to the suppression of shear instabilities by the magnetic field [40]. The plume in Figure 20 (right) exhibits little turbulence due to the suppression of small scales by the strong magnetic field. Hence occurrence of secondary instabilities in such conditions is quite difficult. The other possible reason for not seeing secondary instabilities could be due to the inadequate grids to resolve the secondary instabilities. For example, the jets in Figure 20 are resolved by $100\times 10$ grid points approximately.

VII.3 Relevance to 3D MRTI

While the current work explores the role of reconnection in a 2D MRTI system (at fixed magnetic diffusivity, $\eta$), where the reconnection happens in-plane. However, in a 3D system, the reconnection is free to occur in any plane, and reconnection occurs under complex configurations. Hence, one should be wary of extending the results to a 3D MRTI case. We present our hypothesis on the validity of the results from the current work for the 3D MRTI case. In the current study, we found that the number of reconnection points scales according to $\frac{n_{\times}}{\sqrt{Re_{m}}}=\frac{0.15B_{c}}{B_{0}}$. We hypothesize that the scaling between the number of reconnection events and the magnetic Reynolds number $(n_{\times}\propto\sqrt{Re_{m}})$ may hold in 3D, since there are still only three dependent parameters ($h,v_{A},\eta$) in 3D. Also from theoretical scaling perspective in 3D, $k_{\eta}$ would still have $h^{1/8}$ scaling and the only difference is instead of $L$ we would have $L_{x}\times L_{y}$.

While the scaling remains same, the number of reconnection events could increase in 3D, owing to the greater degrees of freedom for the reconnection to occur in 3D. Hence, the proportionality constant $\frac{0.15B_{c}}{B_{0}}$ could change in the 3D case. We expect the qualitative role of reconnection discussed in the long-term evolution of instability ($\S$III.2) and energy dynamics, conversion from magnetic to kinetic and thermal energies ($\S$IV) to be relevant for the 3D MRTI. However, the quantitative significance, discussed in $\S$VI might vary in 3D.

VII.4 Relevance to astrophysical systems

The study presented here is founded on several simplifications, like the 2D, incompressibility, fully ionized charge-neutral plasma, uniform unidirectional magnetic field, and relatively small kinetic and magnetic Reynolds numbers. However, several features present in the simulation are seen in observations of astrophysical systems. For example, Figure 14 shows the formation of isolated blobs of low-density fluid in the high-density fluid due to reconnection, see the blobs at $(x,y)\approx(2.75,0.25)$ in Figures 14(b), 14(c). Such isolated blobs were observed in the Rayleigh Taylor instability of quiescent solar prominences (for example, figure 8 from [41]), which could be potentially explained due to the reconnection. The formation of reconnection outflow jets was demonstrated in the RTI in solar prominences [39], which were also observed in the current simulations.

VIII Conclusions

Magnetic Rayleigh–Taylor instability (MRTI) and magnetic reconnection are ubiquitous in numerous astrophysical and laboratory systems, but the influence of one on the other has remained unclear. The current study brings these two fundamental aspects together and provides a comprehensive understanding of the role of magnetic reconnection in the MRTI through detailed isolated and statistical analyses. The approximate flux-frozen nature of plasma implies that the onset of MRTI demands deformation of not only fluid lines but also the magnetic field lines. Therefore, the non-linear interaction of plumes leads to magnetic reconnection, making reconnection a fundamental and inevitable phenomenon of MRTI. Despite this, the influence of reconnection on MRTI has remained poorly understood. Here, we summarize the key questions addressed in this work and their corresponding answers:

1.   Is magnetic reconnection essential for the long-term instability evolution?

Yes, we establish that the magnetic reconnection is essential for the long-term evolution of instability. In an ideal-MRTI, where reconnection does not occur, as the instability evolves, the magnetic field lines elongate and the magnetic tension grows, eventually become strong enough to counter the gravitational force, the driving force of instability, arresting further growth of instability. In contrast, in the non-ideal case, reconnection changes the magnetic field topology, shortens field lines, and reduces magnetic tension. Reconnection facilitates plume merger by allowing plasma to cross field lines, resulting in larger plumes that experience stronger buoyancy forces. By moderating magnetic tension and promoting plume coalescence, reconnection is crucial for sustaining the instability’s nonlinear growth over long timescales.

2.   Does reconnection play a crucial role in the energy dynamics of MRTI?

Reconnection plays a critical role in the energy dynamics. The majority of the magnetic to kinetic energy transfer is due to reconnection, especially in weak-field regimes, where it accounts for up to $\approx 80\%$ of this transfer. Through energy injection, magnetic reconnection is expected to improve fluid mixing locally. However, reconnection contributes only modestly ($\approx 3\%$) to total magnetic energy dissipation. The majority of the energy dissipation occurs from the non-reconnecting current sheets that form around the plume head and flanks. The reconnection is also a key element in regulating the magnetic tension in the system (cf. figure 7) and ensuring the instability evolves continuously.

3.   How do the reconnection events vary statistically in the MRTI mixing layer?

The MRTI mixing layer exhibits increasingly complex turbulent dynamics, which drive more frequent interactions and reconnection of magnetic field lines over time. It was found that, for any given magnetic field strength, the number of reconnection events increases self-similarly with time ($n_{\times}\propto t^{3/2}$) and magnetic Reynolds number ($n_{\times}\propto Re_{m}^{1/2}$). But, as the magnetic field strength is increased, the number of reconnection events for a given $Re_{m}$ decreases. This is due to the suppression of turbulence in the mixing layer by the magnetic field.

In summary, reconnection emerges as a fundamental mechanism that controls large-scale instability evolution by alleviating magnetic tension and enabling plume mergers. While the number of reconnection events grows self-similarly, the proportion of reconnecting weak current sheets is increasing. Consequently, the dynamical significance of reconnection in energy dissipation is not significant. However, the reconnection plays a crucial role in the transfer of magnetic energy into the system. Thus, the impact of reconnection on MRTI evolution and dynamics is profound.

Acknowledgements.

MTK is supported by the EPSRC Grant No. EP/W523859/1. AH is supported by STFC Research Grant No. ST/R000891/1 and ST/V000659/1. BS is supported by STFC Research Grant No. ST/V000659/1. The computational time was obtained from the University of Exeter High-Performance Computing facility.

Data Availability Statement

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

Author contributions

M.T.K.: reconnection detection algorithm, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing —original draft, review and editing; A.H.: conceptualization, funding acquisition, project administration, resources, supervision, reconnection detection algorithm, writing—review and editing; B.S.: reconnection detection algorithm, writing—review and editing

Appendix A Issues in the reconnection detection algorithm and their solutions

The algorithm mentioned in $\S$V.1 has a few issues which are described below, again using Fadeev equilibrium. The code is tested for Fadeev and Harris equilibrium in a $4{\times}4$ domain $(L_{x}{=}L_{y}{=}[-2.0,2.0])$, with a $400{\times}400$ grid. The problem is tested for different eccentricity $(\epsilon)$ and angle of orientation $(\theta)$, but we only show a sample case here. Figure 25 shows the null points detected by the code for different $\epsilon$ values and orientations of plasmoids. The value of $\kappa$ is $2\pi$ for $\theta=5^{0},$ $85^{0}$, and $3\pi$ for $\theta=45^{0}$.

From the last row of figure 25, at very low $\epsilon$, i.e., when plasmoids are highly elongated, the code gives a chain of points. This is due to the poor resolution of the plasmoids. When the number of grid points is increased, the code detects the location of these points more accurately, as shown in figure 26, where the grid points are increased to $1600\times 1600$ for the same parameters. The need to resolve the plasmoids is a major issue for the algorithm. Often, in a highly turbulent system, the plasmoids are elongated due to strong shear. But, due to limited resolution capability, it is impossible to resolve all plasmoids, and we often end up with a chain of X-points, resulting in false reconnection points.

The reconnection region typically spans over multiple grid points, resulting in a cluster of points (see figure 26). However, all the points represent one null point (or reconnection point). To address this issue, we use a clustering method [42] to cluster these points into a single point. We use the density-based clustering techniques (DBSCAN [43, 44], HDBSCAN [45]) due to the collocation nature of the points. The DBSCAN method requires the user to provide an $eps$ value, which represents the radius of the neighbourhood region. Hence, DBSCAN is suitable for scenarios when the clusters are similar throughout the space. However, in turbulent systems, the reconnecting clusters are often of very different sizes and occur in various sizes. Thus, specifying a value for $eps$ is not possible. Therefore, we use HDBSCAN, where the algorithm determines the clusters by decreasing the $eps$ value. As the $eps$ value decreases, the cluster becomes smaller in size. The part that remains a cluster across the $eps$ value is considered the final cluster. A drawback of the above method is the potential grouping of two or more clusters close to each other into one. This issue is unavoidable due to the choice of autonomous clustering and could also happen in the DBSCAN method (depending on the $eps$ value provided).

Once the clusters are determined, the reconnecting point in each cluster is to be chosen. Two potential methods we can use to choose the point are — $i)$ the centroid of the cluster, $ii)$ the point with the least magnetic field magnitude. For the former, the mean coordinates of each cluster are calculated 333Since the mean need not always be at the grid point, we determine the nearest grid point close to the determined centroid and assign it as the reconnection point.. A problem with this method is the inaccuracy of the $X-$point. If HDBSCAN considers two or more neighbouring clusters as one cluster, the X-point location based on the centroid is inaccurate. This can be avoided with the latter choice, and hence, we choose the least magnetic field magnitude method. While this method reduces, it cannot eliminate false reconnection points.

Appendix B Self-similarity in 2D MRTI

In $\S$V.2, we mentioned that MRTI is self-similar. However, this was not evidenced there since the focus of the section was on reconnection. Here, we evidence the self-similarity of the present MRTI simulations. One of the standard methods to validate MRTI self-similarity is by verifying the quadratic growth of mixing layer height (cf. $h\propto\frac{(\rho_{h}-\rho_{l})gt^{2}}{(\rho_{h}+\rho_{l})},t\gg 1$) [4, 20]. In the self-similar regime, we expect $\frac{(\rho_{h}+\rho_{l})}{((\rho_{h}-\rho_{l})g}\frac{h}{t^{2}}$ to be approximately constant. One of the other ways of confirming self-similarity is by plotting the energy partition ($C_{ep}$), which is the ratio of turbulent magnetic energy to turbulent kinetic energy over time [20]. From figure 27, we find that these parameters are approximately constant over time, evidencing that the system achieved a self-similar state. Sampling the data between $t=6$ and $t=10$, the largest standard deviation during the time frame is less than $1\%$.