Energy distribution and substructure formation in astrophysical MHD simulations

During substructure formation in magnetized astrophysical plasma, dissipation of magnetic energy facilitated by magnetic reconnection affects the system dynamics by heating and accelerating the ejected plasmoids. Numerical simulations are a crucial tool for investigating such systems. In astrophysical simulations, the energy dissipation, reconnection rate and substructure formation critically depend on the onset of reconnection of numerical or physical origin. In this paper, we hope to assess the reliability of the state-of-the-art numerical codes, PLUTO and KORAL by quantifying and discussing the impact of dimensionality, resolution, and code accuracy on magnetic energy dissipation, reconnection rate, and substructure formation. We quantitatively compare results obtained with relativistic and non-relativistic, resistive and non-resistive, as well as two-and three-dimensional setups performing the Orszag-Tang test problem. We find sufficient resolution in each model, for which numerical error is negligible and the resolution does not significantly affect the magnetic energy dissipation and reconnection rate. The non-relativistic simulations show that at sufficient resolution, magnetic and kinetic energies convert to internal energy and heat the plasma. In the relativistic system, energy components undergo mutual conversion during the simulation time, which leads to a substantial increase in magnetic energy at 20% and 90% of the total simulation time of 10 light-crossing times—the magnetic field is amplified by a factor of five due to relativistic shocks. We also show that the reconnection rate in all our simulations is higher than 0 . 1, indicating plasmoid-mediated regime. It is shown that in KORAL simulations more substructures are captured than in PLUTO simulations.


INTRODUCTION
Dissipation processes in astrophysical plasma, including magnetic reconnection (Biskamp 2000), are of fundamental relevance for our understanding of a variety of observed systems, such as solar flares (Giovanelli 1946;Jiang et al. 2021) or magnetic sub-storms in the Earth's magnetosphere (Akasofu 1968;McPherron 1979).The relative motion in plasmas and gas often leads to the formation of shocks.Non-relativistic magnetized shocks in supernovae remnants are possible sources of acceleration of cosmic rays (Chen & Armstrong 1975;Blandford & Ostriker 1978;Bell 1978;van Marle et al. 2017).Energy dissipation in the relativistic regime leads to spectacular displays, such as jets and relativistic ejections from the accretion systems around compact objects (Giannios et al. 2009;Ripperda et al. 2022), or event horizon scale synchrotron emission (Mehlhaff et al. 2020) and flaring (Dexter et al. 2020;Wielgus et al. 2022) in the hot advection-dominated accretion flows.In the context of accretion onto compact objects, understanding dissipation occurring on small spa-★ E-mail: fatima@camk.edu.pltial scales is crucial to finding realistic sub-grid physics prescriptions for global simulations.
Magnetic reconnection is a process by which the magnetic field lines in a plasma break and reconnect, releasing stored energy in the form of heat, particles/plasmoid acceleration, or radiation.Reconnection often occurs spontaneously and is usually associated with the presence of a current sheet, a region where the magnetic field lines become almost antiparallel and the plasma conductivity is finite.The magnetic field lines can break and reconnect due to the tearing instability, which is driven by the pressure of the plasma and the tension of the magnetic field (Coppi et al. 1966;Komissarov et al. 2007;Del Zanna et al. 2016).Spontaneous reconnection is relatively slow, and the rate is determined by the local plasma conditions (Sweet 1958;Baty 2000).Petschek proposed a shock geometry that allows fast magnetic reconnection to occur (Petschek 1964), this may be realized in magnetohydrodynamic (MHD) simulations for large values of (anomalous) resistivity.In systems with the strongly magnetized plasma, Lazarian & Vishniac (1999) state that the reconnection will always occur at some upper limit of the reconnection rate.Another scenario is the forced magnetic reconnection, which occurs due to external perturbation in a turbulent system (Vekstein & Jain 1998;Potter, M. A. et al. 2019;Srivastava et al. 2019).In this scenario, the reconnection rate can be much faster than spontaneous reconnection, as the external forces can overcome the moderating resistances of the plasma.Such turbulent systems can be found in various environments, such as solar wind, the interstellar medium, or the accretion disks around black holes and neutron stars.
In this work, we study energy dissipation and magnetic reconnection in the MHD framework, using a simple example of a vortical system, the Orszag-Tang (OT) vortex (Orszag & Tang 1979), a popular test problem for numerical MHD codes.In such a system, the magnetic field lines stretch and twist thus facilitating the reconnection process.This test has already been performed with state-of-the-art codes like Athena++ (White et al. 2016), BHAC (Olivares Sánchez et al. 2018), and HARM (Gammie et al. 2003).Here, we implement the OT test in two more state-of-the-art codes used in numerical simulations of accretion.We quantitatively compare results obtained with the two codes of our choice at different resolutions and setups in relativistic/non-relativistic, resistive/non-resistive, and two-dimensional (2D) vs. three-dimensional (3D) configurations, to study how much these different aspects impact the obtained results, characterized by the energy balance and reconnection rate.
The well-established codes we selected for the comparison are the widely used, public PLUTO code (Mignone et al. 2007) and the radiative, general relativistic code KORAL (Sądowski et al. 2013(Sądowski et al. , 2014a)).
PLUTO has extensively been used in simulations of magnetospheric star-disk interaction with alpha-viscous disk in Zanni & Ferreira (2009); Čemeljić (2019), with magneto-rotational instability including alpha-dynamo in Flock et al. (2011), jet launching disks in Tzeferacos et al. (2009), accretion-ejection problem in Stepanovs & Fendt (2014), to mention only some.It was also used in the simulations of star-planet magnetospheric interaction, e.g. in Strugarek et al. (2014) and Varela et al. (2018) and related papers.A radiative module was included in simulations of accretion columns in classical T Tauri stars in Colombo et al. (2019).KORAL code is used to study the accreting compact objects in general relativity involving radiation using M1 closure scheme (Sądowski et al. 2013).The code has been used to study the radiative black hole accretion discs (Sądowski et al. 2014b;Sądowski et al. 2017;Lančová et al. 2019;Chael et al. 2019) as well as super-Eddington accretion onto magnetized neutron stars (Abarca et al. 2021).
The paper is organized as follows: in §2 we review the theoretical framework, including the formalism of the MHD equations.The initial conditions in the OT problem in 2D and 3D setups are given in §3.In §4 we discuss the results in different cases.The reconnection rate is studied in §5.In §6 we present the direct comparison of the results in the two codes we used here and we conclude in §7.

SPECIAL RELATIVISTIC RESISTIVE MHD EQUATIONS
We investigate the energy distribution in astrophysical simulations in the following setups: We begin with presenting the resistive special relativistic MHD equations in Minkowski spacetime, which we then simplify to relativistic ideal MHD and non-relativistic resistive MHD cases.The simulations are performed in the PLUTO and KORAL codes, with the exception of Res-MHD, which is performed in PLUTO alone (KORAL only treats non-resistive MHD equations).
The dynamics of magnetic fluids can be described using the equations of conservation of mass, momentum, and energy, as well as the Maxwell-Faraday, Ampère-Maxwell, and Ohm equations.For a fluid propagating in the laboratory reference frame with bulk velocity  = , the Lorentz factor is defined as Γ = (1 −  2 ) −1/2 , and the fluid four-velocity is  = (Γ, Γ).We denote fluid rest mass density in the fluid frame by  2 , fluid pressure by , fluid internal energy density in the fluid frame by  int , electric field by , and magnetic field by .The  and  fields were redefined to absorb a factor of 1/

√
4 each, so that factors of 1/(4) do not appear in relations such as Eqs. 2 and 3, 7. Furthermore, we define enthalpy density in the fluid frame, momentum density and the total energy density The conservation equations are then where additionally we denote identity matrix with , and the electromagnetic stress tensor with    , hence The Maxwell-Faraday and Ampère-Maxwell equations are respectively, where  is the current density that comes from Ohm's law, where  is the magnetic diffusivity, which is identical to resistivity.The additional condition ∇ •  = 0 from Gauss's law is enforced during the numerical evolution of the magnetic field.
In order to obtain the system of nonrelativistic resistive MHD equations from Eqs. 4-6, we make a number of approximations based on  ≪ 1 and  +  int ≪  assumptions, leading to a following formulation: where the non-relativistic total energy and enthalpy densities are Additionally, Ohm's law in resistive nonrelativistic MHD becomes neglecting the displacement currents (/ = 0) in Eq. 9 to obtain the second equality.
The diffusive time scale,   =  2 / (in conventional units, if  is in cm 2 /s and  is in cm, then   is in seconds) can be compared with the dynamical time scale   = /, where  is the characteristic length scale of the system and  is the characteristic velocity scale.The ratio of the two time scales is known as the magnetic Reynolds number When the typical velocity scale of the system is the Alfvén velocity  A = / √︁ 4, this ratio is called the Lundquist number Astrophysical systems often satisfy the condition  ≫ 1, which is equivalent to  ≫ /  .In such cases, for either relativistic or non-relativistic cases, we can use the ideal MHD approximation As a consequence,  can be readily evaluated and does not need to be evolved with the Ampère-Maxwell equation (Eq.9), simplifying the Maxwell-Faraday equation (Eq.8) for the  field evolution to

ORSZAG-TANG TEST PROBLEM
With implicit inclusion of the most important features of MHD turbulent flow such as energy dissipation and magnetic reconnection (Orszag & Tang 1979;Dahlburg & Picone 1989) the Orszag-Tang vortex is a comprehensive test problem for MHD codes.This problem mostly tests the code performance in simulations with MHD shocks and shock-shock interactions.
We study the energy distribution in different setups by performing the OT test problem simulations using two astrophysical simulation codes: PLUTO (ver. 4.4;Mignone et al. 2007) and KORAL (Sądowski et al. 2014a).The description of our simulations is mostly presented in code units.These are obtained by scaling the equations with fiducial values of certain physical quantities.All velocities are scaled with  0 = , e.g., the statement that  A = 1 in code units signifies that the Alfvén velocity is equal to the speed of light.The density is scaled with some density  0 , the pressure with  0 , and the electromagnetic fields with  0 .The exact value of  0 is immaterial, as long as  0 =  0  2 0 and  0 =  0 √︁ 4 0 .

Two dimensional setup
The simulation is set up in a 2D computational box 0 ≤ ,  ≤ 2 with periodic boundary conditions and the following initial conditions for velocity and magnetic fields (Ripperda et al. 2020): = ṽ(− sin , sin , 0), ( 21) We adopt ṽ = 0.99 0 / √ 2 and B =  0 .The initial density is uniform.In 2D we perform the OT simulations in the range of uniform resolutions from 64 2 to 4096 2 in different setups (Ideal-MHD, Res-MHD, and Rel-MHD), doubling the number of grid points in each dimension to increase the resolution step by step.In 3D we run the Ideal-MHD and Rel-MHD simulations in three resolutions 128 3 , 256 3 , and 512 3 .Only with PLUTO, we run the Res-MHD simulation (in both 2D and 3D) in the resolution 512 3 .Without resistivity, both PLUTO and KORAL are used for Ideal-MHD and Rel-MHD simulations in 2D and 3D setups1 .
All simulations run to the final time  = 10   , where   is the light-crossing time across the typical length in the system.In code units,   = , and we take   = 1.

Three dimensional setup
In order to study the difference between 2D and 3D MHD flows and reconnection, we extend the Orszag-Tang test problem to three dimensions.We set up the simulation in a cubic computational box 0 ≤ (, , ) ≤ 2 with periodic boundary conditions.
For the Rel-MHD simulations, the initial equations are chosen in such a way as to result in a realistic turbulent system, following the definition of a Taylor-Green vortex (Orszag & Tang 1979): where ṽ and B are the same as in the 2D setup.
We find that such initial conditions do not result in a sufficiently turbulent outcome in non-relativistic simulations in 3D, so for Ideal-MHD and Res-MHD simulations in 3D we use different initial conditions, following Mininni et al. (2006): = B(−2 sin 2 + sin , 2 sin  + sin , sin  + sin ), where ṽ = 2 0 and B = 0.8 0 .The initial density is uniform.

ENERGY COMPONENTS IN THE RESULTS
We study the dissipation of magnetic energy and investigate the conversion of energy by following the time evolution of the energy components: the electromagnetic energy density ), the kinetic energy density  K , and internal energy density  int .We study all components in the laboratory frame, thus the kinetic energy and internal energy densities in the relativistic simulations Rel-MHD are computed as follows: Here,  = 4/3 is the polytropic constant.In the non-relativistic limit (simulations Ideal-MHD and Res-MHD) the internal energy density becomes while the kinetic energy density is  as can be seen from Eqs. 13, 14.Another quantity that is a function of space and time is the magnetization defined as  = 2 /( 2 ).
We discuss and compare the averaged energy densities denoted by a bar and computed in 3D through where  is the volume of the simulation box.In 2D the corresponding formula is The results in PLUTO and KORAL simulations are very similar both qualitatively and quantitatively.Unless stated otherwise, we present the PLUTO results.The KORAL results and details of their difference from the PLUTO results are discussed in detail in Section §6.

Ideal-MHD and Res-MHD simulations
In this section, we estimate the numerical dissipation in the simulations and study the effect of resistivity on the evolution of the system.In the left panel of Fig. 1, we plot the time evolution of the averaged squared magnetic field  2 measured in Ideal-MHD simulations for different resolutions.It is clear that at later times the value of  2 increases with an increase in the resolution.This is because in gridbased codes the flux is computed over the surface of every grid cell.In such a calculation there is some amount of computational dissipation, so-called numerical resistivity.Before we study the effect of physical resistivity in simulations, it is important to estimate the numerical dissipation at each resolution and find a reasonable minimal resolution.
We compare the results in non-resistive Ideal-MHD simulations with the Res-MHD simulations set with different physical resistivities ( in Eq. 16), at each resolution 2 .In the right panel of Fig. 1 the results obtained with the resolution 512 2 are shown.We compare  2 of the simulations with  = 0, 10 −4 , 10 −3 , 5 × 10 −3 .The curves corresponding to the Ideal-MHD and Res-MHD simulations with  = 10 −4 are almost overlapping, so at this resolution we estimate the numerical resistivity to be below 10 −4 and conclude that the resolutions higher than 512 2 are reasonably reliable for our simulations with the PLUTO code.
The magnetic energy initially increases and then decreases, forming the hump at 2  in its plot (Fig. 1).This is caused by the compression of a region around a current sheet and subsequent formation of a reconnection layer (at  ≈ 2  ) which then dissipates the magnetic field energy.
In Fig. 2 we show the mass-density plots at  = 2.5  in the simulations Ideal-MHD (numerical resistivity below 10 −4 ) and Res-MHD (physical resisistivity  = 10 −4 ) for the resolution of 4096 2 .In the left panel (Ideal-MHD) we have identified plasmoids (regions of higher density and lower magnetization relative to their surroundings), these are the substructures located in the central region of the simulation box.In the right panel (Res-MHD) the chain of plasmoids does not appear.Similarly, we see no such chain in the simulations with a resistivity larger than 10 −4 .The resistivity of 10 −4 corresponds to the Lundquist number  =   / ≈ 10 4 , with the typical length scale of the system  ≈ 1 and Alfvén velocity  A ≈ 1.This result matches theoretical studies which confirm that the current sheet is plasmoid unstable 3 at  > 10 4 (Loureiro et al. 2007;Ripperda et al. 2020).We also confirm that with a smaller physical resistivity ( < 10 −5 ,  > 10 5 ) some substructures are resolved in the Res-MHD simulations.
We compare the different terms in energy distribution (magnetic energy   , kinetic energy   , internal energy  int , and electric energy   , respectively) in the Res-MHD simulations with  = 5 × 10 −3 and  = 10 −4 (Fig. 3).The first row of this figure shows magnetic energy where the horizontal dashed line, located at   = 0.5, shows the initial value of magnetic energy.We see that with decreasing physical resistivity (from the left panel to the right panel) the rate of magnetic energy decrease becomes smaller.The dissipated magnetic energy converts to the internal energy and heats up plasma as shown in the third row of this plot.We will discuss the energy components in Rel-MHD and Ideal-MHD simulations in the next section.

Ideal-MHD and Rel-MHD simulations
We compare the results of non-relativistic (Ideal-MHD) and relativistic (Rel-MHD) non-resistive MHD simulations in the PLUTO code.The different terms in energy distribution (magnetic energy   , kinetic energy   , internal energy  int , and electric energy   , respectively) are shown in Fig. 4. Panels in the left column show the results for Rel-MHD and in the right column for Ideal-MHD.
The magnetic energy evolution, shown in the first row of panels in Fig. 4, indicates that in simulation Rel-MHD the magnetic energy increases five-fold from the initial value of 0.5 (shown by the black dashed line in both left and right top panels).In the non-relativistic simulation Ideal-MHD there is only a minor initial increase of the magnetic energy followed by a slow decay.
The kinetic energy evolution is presented in the second row of Fig. 4, where a black dashed line is also drawn for reference at the value of 0.5.The kinetic energies were computed using Eqs.27 and 30 for the Rel-MHD and Ideal-MHD simulations, respectively.In Rel-MHD the effect of the Lorentz factor on the kinetic energy leads to an initial value of approximately 0.62, which is higher than the magnetic energy.In contrast, for the Ideal-MHD simulations, the initial value of the kinetic energy is approximately 0.25, half the value of the magnetic energy.Initially, in Rel-MHD the kinetic energy amplifies the magnetic field, while in the nonrelativistic Ideal-MHD case the low value of  K is not enough to amplify the magnetic energy.Thus, in Rel-MHD the effect of kinetic energy on the magnetic energy evolution in the second half of the simulation is significant, causing a secondary increase of  B .In the Ideal-MHD no such effect is observed.
In the third row of panels in Fig. 4, we show the internal energy as computed from Eqs. 28 and 29.Comparison with the first row of panels shows the conversion between magnetic energy and internal energy.
In Rel-MHD, after  ≃ 5  , the large amount of the internal and kinetic energy amplifies the magnetic field.This is visible as the second increase ("hump") in the   curve.Such an outcome in the Rel-MHD simulation offers an explanation for the energy reservoir in magnetized systems like relativistic jets in active galactic nuclei, accretion discs of black holes, and magnetized neutron stars in highenergy astrophysics.In the nonrelativistic Ideal-MHD case, shown in the right panel, the released magnetic energy converts to internal energy and heats up the plasma.In contrast with the relativistic case, the amount of energy in the system is not enough to re-amplify the magnetic field.
The final row of panels in Fig. 4 displays the electric energy, which exhibits a significant evolution in the Rel-MHD simulation.The electric field is a function of magnetic field and velocity (Eq.19), Consequently when the magnetic field is increased around ∼ 2  , the electric energy   also increases.Furthermore, as the system evolves, there is another subsequent increase in   , coinciding with an increase in kinetic energy after 4  .
The sum of all energy components in each of the simulations is conserved over time, as shown in  simulations.Also, the results indicate that Rel-MHD simulation in PLUTO is less dissipative than in KORAL.
Space averaged magnetization in both simulations, Rel-MHD and Ideal-MHD, with the fixed resolution of 512 2 grid cells, is shown as a function of time in Fig. 6.This shows once again how the relativistic system is strongly magnetized and the magnetization increases by the end of the simulation time, while in a non-relativistic simulation, the magnetization does not evolve significantly.

3D simulations
We perform the OT test problem simulations in three dimensions in PLUTO and KORAL with the initial conditions of Eqs. 23, 24 in Rel-MHD simulations and with the initial conditions of Eqs. 25, 26 in Ideal-MHD and Res-MHD simulations.
The time evolution of  2 in the Ideal-MHD simulations is shown in the left panel of Fig. 7.We expect the current sheet to be resolved at time  ≃ 1.5  , because of the increase in magnetic energy discussed in the previous section.We search for the reconnection layers and plasmoids in different slices of the simulation domain at this simulation time.An example of a resulting rest-mass density plot is shown in Fig. 8, which is a slice at  = /2.The plasmoid (in the left panel) is shown at the center of the simulation box, which is zoomed-in at the bottom panel.
We estimate the numerical resistivity at each resolution in Ideal-MHD simulations in 3D by comparing with Res-MHD simulations  for different values of .The plot of  2 with different physical resistivities  = 0, 10 −4 , 10 −3 , 5 × 10 −3 , in the resolution of 512 3 grid cells is shown in Fig. 9 (the method is discussed in §4.1).It is shown that the curves corresponding to  = 10 −4 resistive simulations and the non-resistive Ideal-MHD cases are convergent, so the numerical resistivity in Ideal-MHD simulations with PLUTO at the given resolution is estimated to be ≲ 10 −4 .We expect that at this resolution the current sheets are well resolved.
The rest-mass density plots in the Ideal-MHD simulations (left panel) and resistive Res-MHD simulations with  = 10 −3 (right panel) with the resolution of 512 3 are shown in Fig. 8.The zoomedin frames in the bottom panels show the substructure at the center of each simulation box.From the configuration of the magnetic field which is not shown in this figure, we found that there is a thick current sheet containing a plasmoid in the Ideal-MHD simulation, which is not resolved in the Res-MHD simulation.
The right panel in Fig. 7 shows the time evolution of  2 in the Rel-MHD simulation.It shows that  2 increases to the time  ≃ 7  .At the low resolutions, the magnetic energy drops after this time, but at the high resolution 512 3 , the peak is flattened.We found that at the smaller resolutions, due to the high numerical dissipation, the current sheets are compressed and plasmoids are not resolved.At the high resolution 512 3 , we can see the plasmoid unstable current sheets at different slices in the simulated cubic computational domain.
We show the slice in the rest-mass density at  = /2 in the Rel-MHD simulation with the resolution 512 3 in Fig. 10, with a few magnetic islands in the simulation box 4 .We check the profile of magnetic field components and magnetization in that region.For instance, we take a closer look at one plasmoid located at (, ) = (4.7,3.68).In the right panel, we show the profile of magnetic field components, magnetization, and mass density along the dashed line at  = 3.68 with  ∈ [4, 5.6].The mass density  reaches a local maximum at the position of the plasmoid, while the parallel magnetic field component   , and magnetization  have a minimum local value.Such a profile confirms that there is a plasmoid at this point (Nathanail et al. 2020;Čemeljić et al. 2022).In the same Rel-MHD simulation we made another slice, shown in Fig. A1, through the same simulation box in the  plane at  = 3.68 (where the black dashed line is in Fig. 10).In the top panel we show the reconnection layer and plasmoids.The zoomed plots show the magnetization of the selected reconnection layer.In the next section, we estimate the reconnection rate at this chosen layer.
Using the same method (just described for the 3D Rel-MHD simulation in the last paragraph), we choose the layer shown in Fig. A2 in the 3D Ideal-MHD run.

RECONNECTION RATE
Magnetic reconnection might occur spontaneously due to the internal MHD instability in a resistive model (Sweet 1958;Petschek 1964) or in the ideal MHD as a kink mode (Baty 2000).In a turbulent system, the external perturbation can cause magnetic reconnection in a so-called forced reconnection, where the plasma is in a state of chaotic and unpredictable motion.The magnetic field lines can become distorted and twisted, leading to reconnection (Vekstein & Jain 1998).
Turbulent systems can be found in various environments, such as in the solar wind, in the interstellar medium, and in the accretion disks of black holes and neutron stars.In these environments, magnetic reconnection can lead to a variety of phenomena, such as the acceleration of particles to high energies, the formation of jets and flares, and the heating of the plasma.An external perturbation in turbulent plasma can accelerate the formation of the X-point, causing a reconnection one order of magnitude faster than spontaneous reconnection.Such a reconnection process is complex and still not well understood, and is an active area of research in astrophysics and plasma physics.There are analytical and numerical studies on forced magnetic reconnection including perturbation in the isolated current sheet (Vekstein & Jain 1998;Potter, M. A. et al. 2019), and a   study searching for the observational signatures of simulated forced reconnection in solar flares (Srivastava et al. 2019).
The OT is a vortex problem, for which turbulence develops during evolution.It is shown in the rest-mass density plots (Figs.A3 and  A4) that the current sheets are not formed in isolation, but are a result of evolution of high-density regions, which are driven together by the evolution of the system.Therefore, fast reconnection is expected in our simulations.
Fig. A5 in Appendix A shows selected reconnection layers in the chosen 2D simulations Ideal-MHD and Rel-MHD.When magnetic field lines reconnect, the magnetic tension acts to shorten the field lines and make a magnetic slingshot, which drives the outflow (plasmoids ejection) from both sides of the layer in the parallel direction (Dahlburg & Norton 1995;Linton et al. 2001).
For a steady-state reconnection, the outflow (from the reconnection area) should be balanced with the inflow (toward the reconnection layer) which is shown with the white arrows in the figure.The ratio of inflow and outflow velocity ( in and  out , respectively), is called the reconnection rate  r =  in / out .
The outflow propagates along the background magnetic field lines with the Alfvén speed  A =  √︁ /( + 1), in conventional units.When  ≃ 10,  A ≃ , the reconnection rate can be approximated with  r =  in /.The magnetization values on both sides of the reconnection layer in all simulations are greater than 8, as demonstrated in Appendix A (Figs. A1, A2 and A5).To compute the reconnection rate we average the inflow velocity of 6 grid cells located on both sides of the layer.The structure of the layer is found by the Harris equilibrium method (Harris 1962;Ripperda et al. 2020).
According to analytical and numerical studies, the reconnection rate in 3D might be both lower or higher than in 2D.The reconnection rate depends on different parameters such as the initial setup, strength of the magnetic field, and turbulence of the system.Čemeljić & Huang (2014) studied magnetic reconnection in 2D and 3D geometries using resistive MHD simulations and found that the reconnection rate in 3D was approximately twice as fast as in 2D.Huang & Bhattacharjee (2016) found that in some cases the 3D reconnection rate can be lower than the 2D reconnection rate due to the complex interplay between the plasmoid instability and the turbulent background.
Our study presents various initial setups in both two and three dimensions (Section 3) that affect the magnetization on both sides of the connection region, which in turn influences the reconnection rate.Our Ideal-MHD simulations result in faster reconnection in 3D than in 2D, while the opposite is observed in the Rel-MHD simulations, where the reconnection rate is slower in 3D.In Fig. 11 we show  r as a function of resolution in the simulations with KORAL simulations.We summarize the results of Fig. 11 as follows.
In 2D setups: 1) Results of the Ideal-MHD simulations show that the resolution does not affect the reconnection rate in the resolutions ≥ 256 2 .We confirm that in the non-relativistic simulations, the current sheet is well resolved in the resolutions ≥ 256 2 (It is also shown in the top panels of Fig. 13 at  ≃ 2.5  that the curves of  2 () at higher resolutions are convergent).In the lower resolutions the reconnection rate changes as a function of resolution  r ≈ 0.04 R−0.7 ( R = /100).
2) Results of the Rel-MHD simulations show that the reconnection rate changes as a function of the resolution as  r ≈ 0.25 R−0.45 in the resolutions ≤ 2048 2 .The current sheets and plasmoids are well resolved in the two highest resolutions.
In both Ideal-MHD and Rel-MHD simulations in the lowest resolutions (64 2 and 128 2 ), the numerical resistivity is much higher than 10 −4 , and the current layer is not resolved.The reconnection rate converges to a constant value at a lower resolution in the Ideal-MHD than in the Rel-MHD simulations.Therefore, in Rel-MHD, it is necessary to increase the resolution with respect to the non-relativistic case to reach a reconnection rate limit that is resolution independent.
In 3D setups, the current sheets are not resolved in the resolution 128 3 .With the higher resolutions 256 3 and 512 3 , we do not see a significant effect of the resolution.In KORAL the lowest value of reconnection rate in 2D simulations at the highest resolution is about 0.1 in the Ideal-MHD and about 0.16 in the Rel-MHD.In 3D simulations, the value of the reconnection rate in the highest resolutions is around 0.3 in both Ideal-MHD and Rel-MHD simulations.
Turning to the resisitive simulations, in Fig. 12 we plot the reconnection rates of Res-MHD runs with  = 10 −4 , 10 −3 , and 5 × 10 −3 in the resolution 512 2 in 2D and 512 3 in 3D.The reconnection rate changes as a function of resistivity, increasing by a factor of about 60% in the 3D case and 30% in the 2D one, as the resistivity changes from 10 −4 to 5 × 10 −3 .This increase is much smaller than the factor 7.07 expected from the Sweet-Parker law (  ∝  1/2 ).The dependence seems to be consistent with 1/log , instead.
Given our fairly low resolution and the small number of points, we cannot make definite claims about the functional form of the reconnection rate.However, the reconnection rate we find is consistent with the dependence on the Lundquist number predicted in Petschek reconnection (  = /log , Petschek 1964).The proportionality constant is  = 0.34 for the 3D simulations,5 and  = 0.10 for the 2D simulations.Here we assumed  A  = 1 and we take logarithms to the base 10 (log ≡ log 10 ).Since our flow is not strongly magnetized nor highly turbulent, the reconnection rate in our resistive simulations is below the rates from Lazarian & Vishniac (1999).

CODE COMPARISON
The codes we used in our simulations, PLUTO and KORAL, rely on solving the MHD equations (given in Section §2) employing the finite volume method.The initial equations are typically formulated in terms of the primitive variables, which include the fluid density, pressure, and velocity, as well as the magnetic field (given in Section §3).To solve the equations using the finite volume method, the computational domain is divided into a grid of cells, each of which contains a set of conserved quantities.These conserved quantities are related to the primitive variables through a set of conversion equations, which are typically derived from the conservation laws of mass, momentum, and energy.Although both PLUTO and KORAL employ the same scheme to calculate conserved quantity fluxes at the boundary of each grid cell, the conversion of primitive to conserved quantities differs between the two codes.PLUTO employs the inversion scheme provided by Mignone et al. (2007), while KORAL uses the 1  inversion scheme outlined in Noble et al. (2006).
We perform simulations of the OT test problem with PLUTO and KORAL codes in the simulations Ideal-MHD and Rel-MHD.The same initial conditions are used in both codes.Here we compare the energy components in the results, the ability of the codes to capture substructures, and the reconnection rates.In Fig. 13 we present the time evolution results for the magnetic energy in the Ideal-MHD and Rel-MHD simulations in PLUTO and KORAL.The value of  2 in the simulations Ideal-MHD slightly increases in KORAL with respect to PLUTO.This difference in the value of  2 is more obvious in the lower resolutions and in the later time steps.In addition, in Fig. 5, we showed that at the identical time steps of Rel-MHD simulation, the residual of the total energy in Rel-MHD in KORAL is typically slightly higher than the one in PLUTO.
To investigate the difference between the codes we plot in Fig. A6 of Appendix A relative differences between KORAL and PLUTO of various quantities.In the Ideal-MHD simulations with sufficient resolution for the small numerical resistivity, both PLUTO and KO-RAL show almost the same numerical dissipation.In the Rel-MHD simulations, the difference between the codes is more pronounced.Also, by comparing the results in Ideal-MHD and Rel-MHD simu-Figure 13.Time evolution of  2 in simulations with different resolutions using PLUTO (left panels) and KORAL (right panels) for the simulations Ideal-MHD (top panels) and Rel-MHD (bottom panels).The value of  2 is slightly higher in the simulations with KORAL.Note: the y-axis is common between left and right panels, and the legend is the same for all panels.lations in Fig. 13, we find that the numerical resistivity is negligible at the largest resolution 4096 2 in the Ideal-MHD simulations (the curves of two larger resolutions overlapping) while in the Rel-MHD simulations, one should increase the resolution to obtain a negligible numerical error.
As mentioned in §4.1, we expect the plasmoid unstable current sheets when there is a hump in  2 plot.We show the rest-mass density plot at  = 2.5  in the simulation Ideal-MHD and  = 9  in the simulation Rel-MHD at the highest resolution 4096 2 in Appendix A, Figs.A3 and A4.These density plots confirm that KORAL is more precise than PLUTO in capturing the substructures.
We compare the reconnection rate in the simulation Rel-MHD in PLUTO and KORAL in Fig. 14.In Fig. A6, we show that the residual relative difference between various quantities in the Ideal-MHD simulation is at the level below 1%, so we only compare   in the Rel-MHD simulation.
We observe that in both 2D and 3D setups the reconnection rate in KORAL simulations is higher than in PLUTO simulations.The magnetization on both sides of the reconnection layer directly affects the reconnection rate (which is discussed in §5), and we showed that in KORAL simulations the magnetic energy (and corresponding magnetization) is higher than in PLUTO simulations.This causes a higher reconnection rate in KORAL simulations compared to PLUTO simulations, as shown in Fig. 14.

SUMMARY AND CONCLUSIONS
We investigate how the resolution and dimensionality of the simulation setup affect the energy dissipation, substructure formation, and reconnection rate, all of which are critically dependent in astrophysical simulations on the onset of reconnection.We study these effects by performing the Orszag-Tang test problem in the numerical simulation codes PLUTO and KORAL.We perform a quantitative comparison between the results obtained from various setups, including relativistic, non-resistive MHD (Rel-MHD), non-relativistic, non-resistive MHD (Ideal-MHD), non-relativistic, resistive MHD (Res-MHD), in 2D as well as 3D simulations.
First of all, we estimated the numerical resistivity of the simulations in each resolution to find a sufficient resolution in which we can resolve the substructures and study the energy conversion in our simulations.We used PLUTO code in resistive and non-resistive modes (Res-MHD and Ideal-MHD, respectively) in non-relativistic simulations.We show that the numerical resistivity in the resolution 512 2 in both 2D and 3D setups is  ≈ 10 −4 , which is also the limit of the formation of a plasmoid unstable current sheet.
After finding the sufficient resolution for overcoming the effects of numerical resistivity, we study energy conversion in Ideal MHD, Rel-MHD, and Res-MHD simulations.We showed that in Ideal-MHD and Res-MHD simulations magnetic energy converts into internal energy and heats up the plasma.In Ideal-MHD simulation a part of magnetic energy converts to kinetic energy which accelerates the plasmoids out of the reconnection layer.We also show that in Res-MHD simulations, as expected, the magnetic energy dissipation increases with increasing physical resistivity.In higher resistivity cases, there is a corresponding increase in internal energy.
In relativistic simulations, Rel-MHD, we find that the relativistic shocks amplify the magnetic field with the magnetic energy  B increasing by a factor of five at  = 20% of total simulation time.It is also shown that magnetic energy converts into internal and kinetic energies which amplify the magnetic field for the second time during our simulation.The second increase in magnetic energy at  = 90% of total simulation time is coincident with the formation of a set of plasmoid unstable current sheets.
We also compare two state-of-the-art codes, PLUTO and KORAL, in both non-relativistic and relativistic simulations.Our findings indicate that in both Ideal-MHD and Rel-MHD simulations, KORAL simulations show higher magnetic energy,  2 , (implying less magnetic dissipation) compared to PLUTO with the difference more prominent at low resolutions.We show that in resolution 1024 2 , in the Ideal-MHD simulations, the relative difference of relevant quantities in PLUTO and KORAL is less than 10 −2 , while in the Rel-MHD simulations, for some quantities the residual reaches 0.1.In the highest resolution run (4096 2 ), we found that KORAL captures more substructures than PLUTO in both Ideal-MHD and Rel-MHD simulations.We show that the reconnection rate in all simulations in KORAL is higher than that in PLUTO-it is caused by higher magnetization in the reconnection layer region in KORAL.
We study the effect of resolution on the reconnection rate   in our simulations.As expected, numerical resistivity influences the reconnection rate.Increasing the resolution leads to a decrease in both numerical dissipation and reconnection rate.In 2D simulations,   is initially a function of scaled resolution ( R = /100) as  r ≈ 0.04 R−0.7 (Ideal-MHD) and  r ≈ 0.25 R−0.45 (Rel-MHD).In each set of simulations, we find a resolution beyond which the reconnection rate is no longer affected by the resolution, and we find the limiting reconnection rate in this limit: in 2D simulations in KORAL, in the Ideal-MHD runs,  r = 0.1 for resolutions ≥ 512 2 ; in the Rel-MHD,  r ≈ 0.18 for resolutions ≥ 2048 2 .In PLUTO simulations, the reconnection rate is lower than that in KORAL simulations.In PLUTO, in Ideal-MHD   ≈ 0.03, in Rel-MHD   ≈ 0.05.
We conclude that the Rel-MHD simulations should be performed at resolutions at least four times larger than in the non-relativistic Ideal-MHD simulations, to reach a negligible effect of the resolution on the reconnection rate.
In 3D simulations in KORAL the Ideal-MHD and Rel-MHD simulations are not directly comparable since we initialized the velocity and magnetic fields differently.Still, in both setups, the results are remarkably similar, with the effect of resolution on  r not significant in higher resolutions.In both Ideal-MHD and Rel-MHD simulations with resolution 512 3 the reconnection rate  r ≃ 0.3 (Fig. 11).
When comparing the reconnection rate in 2D and 3D setups, it is crucial to consider several parameters, such as the initial setup, the strength and topology of the magnetic field, and the turbulence of the system.In setups with the equivalent magnetization and turbulence levels, we show that the reconnection rate in 3D ideal MHD simulations is lower than that observed in 2D simulations.This trend is particularly notable in relativistic simulations when comparing the 2D and 3D setups.However, in the resistive runs (Res-MHD) the trend is the opposite, the reconnection rate is about a factor of 3 smaller in 2D simulations than in 3D ones.We also show that in the resistive simulations, the reconnection rate seems to be well approximated by a   ∝ 1/log  dependence, reminiscent of Petschek's fast reconnection (Petschek 1964).
The results presented here add to the information needed to evaluate the behavior of numerical MHD codes in different setups.The performance of the codes can be evaluated and compared only with a detailed account of the relation between the substructure formation and the amount of energy in each component.By using the standard Orszag-Tang test, we provided detailed quantitative information on energy components, reconnection rates and substructure formation.Our approach can be followed-and the results compared-for other codes.
A caveat in our work here is that, because of the computational expense, we did not follow the convergence of the results in 3D up to the same resolutions as we did in the 2D setups.The new generation of simulations will unavoidably need such an update in benchmarking.The convergence of vorticity will be addressed in future work.In Fig. A6, we plot the residual quantities   = | KORAL −  PLUTO |/ KORAL ( represents the compared quantity) to clarify the difference between PLUTO and KORAL simulations.The black dashed curves correspond to the Ideal-MHD simulation and the blue solid curves correspond to the Rel-MHD simulation.We compute   in the results with the resolution of 1024 2 , at which the numerical dissipation is small.In the Ideal-MHD simulation, the residuals of magnetic energy   and magnetization  are less than 10 −2 while in the Rel-MHD simulation, the residuals reach 0.1.In the Ideal-MHD simulation, the residual of kinetic energy   is less than 10 −2 , while in the Rel-MHD simulation, it is less than 10 −1 .The residuals of internal energy  int and density  ≡ Γ in the Ideal-MHD simulation are of the order of 10 −4 and in the Rel-MHD simulation they are of the order of 10 −2 .We conclude 10% level code consistency for Rel-MHD and 1% level consistency for non-relativistic Ideal-MHD simulations.

Figure 1 .
Figure 1.The time evolution of  2 in PLUTO simulations: the Ideal-MHD case with different resolutions (left panel), and Res-MHD case with different physical resistivities  for the resolution of 512 2 (right panel).The unit of time is   = /.

Figure 2 .
Figure 2. The rest mass-density at  = 2.5  in the resolution of 4096 2 with PLUTO in the Ideal-MHD simulations (left panels) and Res-MHD simulations with physical resistivity  = 10 −4 (right panels).Plasmoids, zoomed-in in the bottom panels, form only in a case with sufficiently low resistivity, corresponding to a Lundquist number larger than  ∼ 10 4 .

Figure 3 .
Figure 3. Energy distribution in Res-MHD simulations with physical resistivities  = 5 × 10 −3 (left panels) and  = 10 −4 (right panels) at the resolution of 512 2 grid cells.The horizontal black dashed lines in the   and   panels indicate the initial value (0.5) of the magnetic energy   .See the detailed discussion in §4.1.The dissipated magnetic energy heats up the plasma.

Figure 4 .
Figure 4. Energy distribution in Rel-MHD and Ideal-MHD PLUTO simulations at the resolution of 512 2 grid cells are shown in the (left) and (right) panels, respectively.The horizontal black dashed lines in the panels with   and   indicate the initial value of the magnetic energy   = 0.5.See the detailed discussion in §4.2.

Figure 6 .
Figure 6.Magnetization ( =  2 / in code units) in Ideal-MHD and Rel-MHD simulations with PLUTO at the resolution 512 2 grid cells.

Figure 7 .
Figure 7.The time evolution of  2 in 3D simulations with PLUTO for the simulations Ideal-MHD (left panel) and Rel-MHD (right panel).

Figure 8 .
Figure 8.The slice in  = /2 in the simulation box of the rest-mass density  for a vortex at  = 1.5  at the resolution of 512 3 in PLUTO.Left panel: Ideal-MHD.Right panel: Res-MHD with  = 10 −3 .The zoomed-in panels show the current layer in the middle of the simulation boxes.Plasmoids form only in the cases with sufficiently low resistivity, corresponding to a Lundquist number larger than 10 4 ( ≲ 10 −4 ).

Figure 9 .
Figure 9.The time evolution of  2 in 3D Res-MHD simulations with PLUTO, at the resolution 512 3 with different physical resistivities.

Figure 10 .
Figure 10.Left panel: a slice of the rest-mass density at  = /2 in the Rel-MHD simulation in 3D at a resolution 512 3 with PLUTO.The streamlines indicate the magnetic field lines and the white circles show plasmoids.Right panel: the magnetic field components, magnetization, and density profile along the black dashed line at  = 3.68, shown in the left panel.

Figure 11 .
Figure 11.The reconnection rate as a function of resolution in the simulations Ideal-MHD and Rel-MHD with KORAL.

Figure 12 .
Figure 12.The reconnection rate as a function of resistivity for resistivities ≥ 10 −4 1958 in 2D (red circles) and 3D (blue stars) Res-MHD simulations with PLUTO.The change is cosistent with 1/log  dependence,  being the Lundquist number.

Figure 14 .
Figure 14.The reconnection rate as a function of resolution in simulations Rel-MHD in 2D and 3D.Red symbols indicate simulations with KORAL and blue symbols with PLUTO.

Figure A1 .
Figure A1.In the top panel is shown a slice in the rest-mass density at  = 3.68 in the Rel-MHD simulation in 3D with PLUTO, at the resolution 512 3 .The white box shows the reconnection layer contains a plasmoid.The streamlines indicate the magnetic field.In the bottom panel we plot the magnetization of the selected region (white box).

Figure A2 .
Figure A2.In the top panel is shown a slice in the rest-mass density at  =  in the Ideal-MHD simulation in 3D with PLUTO, at the resolution 512 3 .The white box shows the reconnection layer contains a plasmoid.The streamlines indicate the magnetic field.In the bottom panel we plot the magnetization of the selected region (white box).

Figure A3 .
Figure A3.The rest-mass density  for a vortex at  = 2.5  in the Ideal-MHD simulation at the resolution of 4096 2 in PLUTO (left panels) and KORAL (right panels).Due to lower numerical dissipation, KORAL is more precise in capturing the substructure in the simulations.

Figure A4 .
Figure A4.Rest-mass density  for a vortex at  = 9  of the Rel-MHD simulation with the resolution 4096 2 in PLUTO (left panels) and KORAL (right panels).Due to lower numerical dissipation, KORAL is more precise in capturing the substructure.

Figure A5 .
Figure A5.The reconnection layer is visible in the magnetization in simulations with KORAL at the resolution 4096 2 grid cells in the center of the simulation box in the Ideal-MHD simulation (top panel) and along the (0, )-( , 2 ) line in the simulation box in the Rel-MHD simulation (bottom panel).The streamlines show the magnetic field.

Figure A6 .
Figure A6.The residual of quantities ,   = | KORAL −  PLUTO |/ KORAL , in KORAL and PLUTO with the resolution 1024 2 .We show (from top to bottom) the residuals for the magnetic energy,    , magnetization,   , kinetic energy,    , internal energy,   int and density,   .