© The Authors 2023

1. Introduction

Magnetic reconnection plays a crucial role in a wide range of phenomena in the Universe. For instance, it sparks high-energetic bursts in the accretion disc around the black hole in active galactic nuclei (Liu et al. 2002), it is the basis of thermonuclear power devices, such as the tokamak (Furth et al. 1973), and it strongly affects space weather (Paschmann et al. 1979). On the Sun in particular, this physical process has been shown through numerical experiments to cause several remarkable solar events, such as Ellerman bombs (EBs) and ultraviolet (UV) bursts (e.g., Hansteen et al. 2017, 2019; Danilovic 2017; Nóbrega-Siverio et al. 2017; Peter et al. 2019; Ni et al. 2021), surges and coronal jets (e.g., Yokoyama & Shibata 1995, 1996; Nishizuka et al. 2008; Pariat et al. 2009; Moreno-Insertis & Galsgaard 2013; Archontis & Hood 2013; Fang et al. 2014; Toriumi et al. 2015; Nóbrega-Siverio et al. 2016; Wyper et al. 2016, 2017; Karpen et al. 2017; Luna & Moreno-Insertis 2021; Nóbrega-Siverio & Moreno-Insertis 2022), and flares (e.g., Yokoyama & Shibata 2001; Masson et al. 2009; Cheung et al. 2019; Rempel et al. 2023; Chen et al. 2023), to mention some.

Theoretical reconnection models are commonly divided into two types: slow-reconnection and fast-reconnection. The slow-reconnection model developed by Sweet (1958a,b) and Parker (1957) assumes constant diffusivity over the whole reconnection site and predicts exactly one-half of the inflowing magnetic energy to be converted into heat and the other half into kinetic energy. Nonetheless, the Sweet-Parker model is not efficient enough to reproduce the relatively high reconnection rate observed in flares (e.g., Priest 2014, and references therein). The fast-reconnection model developed by Petschek (1964) instead assumes a diffusion layer limited to a small segment of the boundary layer between the opposing magnetic fields with slow-mode shock waves propagating from the diffusion region. Most of the energy conversion in this model takes place at the shocks, and for a specific heat ratio of $\gamma =\frac{5}{3}$, two-fifths of the inflowing magnetic energy is turned into heat and the remaining three-fifths into kinetic energy. This model predicts a reconnection rate that is high enough to reproduce flares. The Sweet-Parker model and the Petschek model are both steady-state models that assume that the current sheets are stable and do not break. However, reconnection theory has shown that current sheets tend to undergo different resistive instabilities, such as the tearing instability (Furth et al. 1963), causing plasmoids (magnetic islands) to appear and move along the current-flow lines. As a consequence, the reconnection rate and energy conversion rate may deviate from the values predicted analytically with the Sweet-Parker and the Petschek model, and careful analysis is therefore required when studying non-stationary reconnection through numerical simulations.

Mimicking magnetic reconnection processes from a numerical perspective is challenging due to the complex behaviour of the electrical resistivity, η, which appears in Ohm’s law as the ratio of the electric field strength and the current density in the rest frame of the fluid. In the solar atmosphere, this coefficient is commonly derived from kinetic theory of particle collisions and given by Spitzer resistivity (Spitzer 1962). However, under some conditions, such as regions of strong magnetic field gradients, plasma instabilities can affect the dynamics of the charged particles and can cause the resistivity to rise beyond the Spitzer value (Roussev et al. 2002). This effect, known as anomalous resistivity, is also a necessary component to support the theory of dissipation of direct currents (Heyvaerts & Priest 1984) as a significant source of coronal heating because the collisional Spitzer resistivity is too small to dissipate such strong currents (Adamson et al. 2013). In addition, we need to take into account the diffusive effects caused by the discrete nature of numerical codes, which are often significantly greater than those caused by the physical resistivity. Especially in numerical models of the solar atmosphere, regions of large magnetic field gradients require a diffusivity that is much larger than the Spitzer resistivity in order to become numerically resolvable. Because of this, it is common to apply ad hoc terms for anomalous resistivity (Sato & Hayashi 1979; Nordlund & Galsgaard 1995; Roussev et al. 2002; Vögler et al. 2005; Felipe et al. 2010; Adamson et al. 2013; Rempel 2014, 2017; Przybylski et al. 2022) that are set to be large around current sheets in order to dissipate them until they become numerically resolvable, but stay small elsewhere in order to keep the Reynolds and Lundquist numbers relatively high.

For a steady Sweet-Parker- or Petschek-like reconnection model, it is sufficient to use a localised anomalous resistivity model, which means that the resistivity is set to a non-zero value (or to a function of spatial coordinates) in a specific location and zero elsewhere (Innes & Tóth 1999). Non-steady reconnection models with a plasmoid instability can be simulated by using a more adaptive anomalous resistivity model, for instance by enhancing the resistivity when the electron drift velocity or the current density surpass a given threshold value (e.g., Sato & Hayashi 1979), or by applying a fourth-order hyper-diffusive operator consisting of a small global diffusive term and a location-specific diffusion term (e.g., Nordlund & Galsgaard 1995; Gudiksen et al. 2011). However, if the numerical resolution is sufficiently high in areas of strong magnetic field gradients, it is even possible to successfully simulate reconnection with a plasmoid instability without adding any anomalous resistivity terms and only using the actual resistivity in the solar atmosphere (e.g., Ni et al. 2021).

In this paper, three different resistivity models are applied on two numerical experiments for the purpose of analysing their effects on magnetic reconnection. The first experiment mimics a 2D simulation by Syntelis et al. (2019). This enables us to compare our results with already published results that were obtained using a different numerical code. The second experiment simulates a 1D Harris current sheet. We can therefore study the diffusive effects that the resistivity models have in a simple setup.

The structure of the paper is as follows. Section 2 describes the numerical code and the model equations (Sect. 2.1) we used for our experiments, the resistivity models (Sect. 2.2), and the setup for the numerical experiments (Sect. 2.3). Section 3 gives a detailed analysis of the results for the 2D experiment (Sect. 3.1) and the 1D experiment (Sect. 3.2). Finally, Sect. 4 contains a brief discussion of the key results of our study and summarises the conclusions.

2. Numerical model

The simulations of this paper were performed with the Bifrost code (Gudiksen et al. 2011). Bifrost is a massively parallel 3D code that solves the equations of magnetohydrodynamics (MHD) on a staggered grid using a sixth-order differential operator to discretise the spatial derivatives, supported by fifth-order interpolation operators. For the time-stepping, we chose a third-order method (Hyman 1979). The code is modular and can take various physical ingredients into account depending on the experiment.

2.1. Model equations

The model equations for our experiments are given by

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}& =-\mathrm{\nabla }·(\rho \mathbf{u}),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial (\rho \mathbf{u})}{\partial t}& =-\mathrm{\nabla }·(\rho \mathbf{u}\otimes \mathbf{u}-\overline{\overline{\tau }})-\mathrm{\nabla }P+\mathbf{J}\times \mathbf{B}+\rho \mathbf{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{B}}{\partial t}& =-\mathrm{\nabla }\times (-\mathbf{u}\times \mathbf{B}+\overline{\overline{\eta }}\mathbf{J}),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial e}{\partial t}& =-\mathrm{\nabla }·(e\mathbf{u})-P\mathrm{\nabla }·\mathbf{u}+Q_J+Q_V+Q_C,\hfill \end{array}$$(4)

where ρ, u, e, and B are the mass density, fluid velocity, internal energy per unit volume, and the magnetic field, respectively. $\overline{\overline{\tau }}$, P, J, g, $\overline{\overline{\eta }}$, Q_J, Q_V, and Q_C are the viscous stress tensor, gas pressure, electric current density, gravitational acceleration, electrical resistivity tensor, Joule heating, viscous heating, and the Spitzer thermal conductivity term, respectively. Other terms such as non-equilibrium ionisation, ambipolar diffusion, Hall effect, radiative cooling, and optically thin losses are neglected in our experiments. The gravitational term ρg, with g = 0.274 km s⁻², and the Spitzer thermal conductivity term Q_C are only included in the first experiment of this paper (Sect. 2.3.1).

For the equation-of-state, we used the same equation as Syntelis et al. (2019), that is, an electrically neutral ideal gas with a specific heat ratio of $\gamma =\frac{5}{3}$ and a mean molecular weight of μ = 1.2, where P and e are related to the mass density, ρ, and temperature, T, as follows:

$$\begin{array}{cc}& P=\frac{\rho k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}},\hfill \end{array}$$(5)

$$\begin{array}{cc}& e=\frac{P}{(\gamma -1)},\hfill \end{array}$$(6)

where k_B and m_H are the Boltzmann constant and mass of hydrogen, respectively.

2.2. Electrical resistivity models

For the purpose of analysing the effects of the electrical resistivity model on the reconnection in the corona, three different approaches were compared: 1) the default way of handling magnetic resistivity in Bifrost, by means of hyper-diffusion (Gudiksen et al. 2011), hereafter referred to as the Gudiksen-11 model (see Sect. 2.2.1), 2) a resistivity that scales linearly with the current density as was used by Syntelis et al. (2019) for their 2D flux cancellation simulation, which is mimicked in this paper (see Sect. 2.3.1), hereafter referred to as the Syntelis-19 model (see Sect. 2.2.2); and 3) a resistivity that scales quadratically with the electron drift velocity employed by Yokoyama & Shibata (1994) for their simulation of an emerging coronal loop, hereafter referred to as the YS-94 model. Inspired by Sato & Hayashi (1979), the latter resistivity model has been used in several other papers (e.g., Shibata et al. 1992, 1993; Yokoyama & Shibata 1996; Matsumoto et al. 2004).

For later reference, we introduce here the definitions of the Reynolds number, Re, and Lundquist number, S_L,

$$\begin{array}{cc}& Re\equiv \frac{|\mathbf{u}|L_B}{\eta },\hfill \end{array}$$(7)

$$\begin{array}{cc}& S_L\equiv \frac{v_A{L}_B}{\eta },\hfill \end{array}$$(8)

where L_B ≡ (|J|/|B|)⁻¹ is the characteristic length of the magnetic field, and $v_A\equiv |\mathbf{B}|/\sqrt{{\mu }_0\rho }$ is the Alfvén speed of the plasma, where μ₀ is the vacuum permeability.

2.2.1. Gudiksen-11 model

Based on the resistivity model developed by Nordlund & Galsgaard (1995), the Gudiksen-11 resistivity consists of two major terms. The first term is an electrical diffusive speed, U_m, with the x_i component defined by

$$\begin{array}{c}\hfill U_{m,i}={\nu }_1{c}_f+{\nu }_2|u_i|+{\eta }_3\mathrm{\Delta }{x}_i|{\mathrm{\nabla }}_{\perp }{u}_i|,\end{array}$$(9)

where ν₁, ν₂, and η₃ are scaling factors for the fast-mode wave velocity, bulk velocity, and gradients in the velocity perpendicular to the magnetic field, respectively; and $c_f\equiv \sqrt{c_s^2+v_A^2}$ is the fast-mode speed, with the sound speed c_s given by $c_s\equiv \sqrt{\gamma P/\rho }$. In our experiments, we set ν₁ = 0.03, ν₂ = 0.2, and η₃ = 0.2, which are typical values used in Bifrost simulations. In Sect. 3.1.5 we discuss how modifying these free parameters affects the results.

The second term is a positive definite quenching operator defined by

$$\begin{array}{c}\hfill {\mathbb{Q}}_i(g)\equiv \frac{|{\mathrm{\Delta }}_i^{2}g|}{|g|+|{\mathrm{\Delta }}_i^{2}g|/q_{max}},\end{array}$$(10)

where ${\mathrm{\Delta }}_i^2$ is the second-order difference operator in the x_i-direction, g is the first-order derivative (with respect to any spatial coordinate) of any MHD variable, and q_max is the maximum quenching factor. For any perturbation of the wavenumber k, this term quickly approaches q_max as k → ∞ and decreases with k² as k → 0, hence ensuring that perturbations with a wavelength of same order as the grid size are heavily damped, while perturbations with wavelengths that are more than one order of magnitude larger than the grid size are only slightly damped. We used q_max = 8 because this has been empirically shown to work well when Bifrost was used to solve standard test problems.

Thus, the hyper-diffusive resistivity of Bifrost can be written as a diagonal tensor, ${\overline{\overline{\eta }}}_{G11}$, given by

$$\begin{array}{cc}& {\eta }_{G11,xx}=\frac{{\eta }_3}{2}[U_{m,y}\mathrm{\Delta }y{\mathbb{Q}}_y\left(\frac{\partial B_z}{\partial y}\right)+U_{m,z}\mathrm{\Delta }z{\mathbb{Q}}_z\left(\frac{\partial B_y}{\partial z}\right)],\hfill \\ \hfill & {\eta }_{G11,yy}=\frac{{\eta }_3}{2}[U_{m,z}\mathrm{\Delta }z{\mathbb{Q}}_z\left(\frac{\partial B_x}{\partial z}\right)+U_{m,x}\mathrm{\Delta }x{\mathbb{Q}}_x\left(\frac{\partial B_z}{\partial x}\right)],\hfill \\ \hfill & {\eta }_{G11,zz}=\frac{{\eta }_3}{2}[U_{m,x}\mathrm{\Delta }x{\mathbb{Q}}_x\left(\frac{\partial B_y}{\partial x}\right)+U_{m,y}\mathrm{\Delta }y{\mathbb{Q}}_y\left(\frac{\partial B_x}{\partial y}\right)],\hfill \\ \hfill & {\eta }_{G11,xy}={\eta }_{G11,yx}={\eta }_{G11,yz}={\eta }_{G11,zy}={\eta }_{G11,xz}={\eta }_{G11,zx}=0.\hfill \end{array}$$(11)

This resistivity model ensures that the resistive terms in the induction and energy equation become significant only in the regions in which the diffusive velocity is high because of the high fast-mode velocity, advective velocity, or strong magnetic shocks along with strong gradients in the magnetic field, which allow the Reynolds number to stay high outside these regions.

2.2.2. Syntelis-19 model

The Syntelis-19 resistivity, η_S19, is a scalar function given by

$$\begin{array}{cc}\hfill {\eta }_{S19}& =\{\begin{array}{cc}{\eta }_0,\hfill & |\mathbf{J}|<J_{\mathrm{crit}}\hfill \\ {\eta }_0+{\eta }_1|\mathbf{J}|/J_{\mathrm{crit}},\hfill & |\mathbf{J}|\ge J_{\mathrm{crit}}.\hfill \end{array}\hfill \end{array}$$(12)

Syntelis et al. (2019) used η₀ = 3.78  ×  10⁻² km² s⁻¹, η₁ = 3.78  ×  10⁻¹ km² s⁻¹, and J_crit = 5.00 × 10⁻⁴ G km⁻¹. In our experiments, we instead chose η₁ = 7.56 km² s⁻¹ in order to obtain approximately the same inflow Alfvén Mach number as when applying the Gudiksen-11 model on the 2D flux cancellation experiment (Sect. 3.1), as well as an average current sheet length similar to that of Syntelis et al. (2019). This change was needed because the MHD solver scheme of Bifrost and the Lare3D code employed by Syntelis et al. (2019) are different. The Lare3D code is a Lagrangian-Eulerian Remap code (Arber et al. 2001).

2.2.3. YS-94 model

The YS-94 resistivity, η_YS94, is defined as

$$\begin{array}{cc}\hfill {\eta }_{{}_{\mathrm{YS}94}}& =\{\begin{array}{cc}0,\hfill & v_d\le v_c\hfill \\ min(\alpha {(\frac{v_d}{v_c}-1)}^2,{\eta }_{max}),\hfill & v_d>v_c,\hfill \end{array}\hfill \end{array}$$(13)

where $v_d=\frac{J}{n_{\mathrm{e}}e}$ is the electron drift velocity, and v_c, α, and η_max are free parameters. Yokoyama & Shibata (1994) used v_c ∈ [4.16 × 10⁻⁷, 8.32 × 10⁻⁶] km s⁻¹, α ∈ [0.20, 2000] km² s⁻¹, and η_max = 2000 km² s⁻¹ (normalisation units and formulae are extracted from Nozawa et al. 1992 and Yokoyama & Shibata 1996).

In our simulations, we used v_c = 8.3 × 10⁻⁶ km s⁻¹, α = 4.0 × 10⁻⁸ km² s⁻¹, and η_max = 2000 km² s⁻¹ in order to obtain a similar inflow Alfvén Mach number in the 2D flux cancellation simulation as when using the other resistivity models. With this, we applied a much lower value of the scaling factor α than Yokoyama & Shibata (1994) used in their study of current sheets located in the convection zone. Our case deals with reconnection in current sheets that are located in the corona, where the density is several orders of magnitude lower. This causes the drift velocity in current sheets to become several orders of magnitude higher. It is therefore logical that a weaker scaling factor between resistivity and drift velocity is needed here. In addition to the resistivity given by Eq. (13), we added a background uniform resistivity of η₀ = 4.00 × 10⁻² km² s⁻¹ when using this model, similar to that of Syntelis-19.

2.3. Numerical experiments

2.3.1. 2D flux cancellation experiment

The first experiment mimics the case 1 simulation by Syntelis et al. (2019), in which reconnection is driven by converging opposite polarities at the solar surface, leading to flux cancellation in a 2D atmosphere. The computational domain was given by x ∈ [ − 30, 30] Mm and z ∈ [0, 30] Mm, and it was discretised over 2048 × 1024 grid points. The initial magnetic field was a superposition of two sources of opposite polarity placed below the photosphere, along with a horizontal uniform background magnetic field. In 2D, the magnetic field strength from one source with a flux of F at a given distance r is F/(πr), with the direction given by unit vector $\widehat{\mathbf{r}}=\mathbf{r}/r$. Thus, the initial magnetic field is given by

$$\begin{array}{c}\hfill \mathbf{B}(x,z,t=0)=\frac{F}{\pi }\frac{{\mathbf{r}}_1}{r_1^2}-\frac{F}{\pi }\frac{{\mathbf{r}}_2}{r_2^2}-B_0\widehat{\mathbf{x}},\end{array}$$(14)

where F = 2500 G Mm is the flux of each source, B₀ = 45 G is the magnetic field strength of the horizontal background field, and

$$\begin{array}{cc}& {\mathbf{r}}_1=(x+d_s)\widehat{\mathbf{x}}+(z-z_0)\widehat{\mathbf{z}},\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\mathbf{r}}_2=(x-d_s)\widehat{\mathbf{x}}+(z-z_0)\widehat{\mathbf{z}},\hfill \end{array}$$(16)

where d_s = 1.8 Mm is the initial half-separation distance between the sources, and z₀ = −0.36 Mm is the height at which the sources are located.

The initial temperature profile of Syntelis et al. (2019), set to mimic the C7 model of Avrett & Loeser (2008), is given by

$$\begin{array}{c}\hfill T(x,z,t=0)=T_{\mathrm{pho}}+\frac{T_{\mathrm{cor}}-T_{\mathrm{pho}}}{2}[tanh\left(\frac{z-z_{\mathrm{cor}}}{w_{\mathrm{tr}}}\right)+1],\end{array}$$(17)

with T_pho = 6109 K and T_cor = 0.61 MK. For the location of the bottom of the corona and the width of the transition region, we used z_cor = 2.31 Mm and w_tr = 0.09 Mm in our simulations. The initial mass density was found by requiring hydrostatic equilibrium, ∂P/∂z = −ρg, and a photospheric density of ρ_pho = 1.67 × 10⁻⁷ g cm⁻³. With P given by the ideal gas law and T given by Eq. (17), the following analytical solution was found:

$$\begin{array}{c}\hfill \rho (x,\stackrel{\sim }{z},t=0)={\rho }_{\mathrm{pho}}{e}^{-2{\chi }_0(\stackrel{\sim }{z}+{\stackrel{\sim }{z}}_c)}{\left(\frac{T_{\mathrm{pho}}+T_{\mathrm{cor}}{e}^{2\stackrel{\sim }{z}}}{T_{\mathrm{pho}}+T_{\mathrm{cor}}{e}^{-2{\stackrel{\sim }{z}}_c}}\right)}^{{\chi }_0-{\chi }_1}\frac{T_{\mathrm{pho}}}{T(\stackrel{\sim }{z})},\end{array}$$(18)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{z}\equiv \frac{z-z_{\mathrm{cor}}}{w_{\mathrm{tr}}},{\stackrel{\sim }{z}}_c\equiv \frac{z_{\mathrm{cor}}}{w_{\mathrm{tr}}},{\chi }_0\equiv \frac{\mu m_{\mathrm{H}}g{w}_{\mathrm{tr}}}{2k_{\mathrm{B}}{T}_{\mathrm{pho}}},{\chi }_1\equiv \frac{\mu m_{\mathrm{H}}g{w}_{\mathrm{tr}}}{2k_{\mathrm{B}}{T}_{\mathrm{cor}}}.\end{array}$$(19)

Initial magnetic field, temperature, and mass density computed from the above equations are shown in Fig. 1. The figure shows that the initial conditions of Syntelis et al. (2019) have indeed been successfully mimicked.

Fig. 1. Initial conditions for the 2D flux cancellation experiment mimicking (Syntelis et al. 2019). Top: map of the temperature with the magnetic field topology superimposed. Middle: vertical component of the magnetic field, Bz, at z = 0. Bottom: stratification of the temperature (black) and mass density (red).

For the bottom boundary conditions, we used a driving mechanism where the horizontal velocity u_x is defined as

$$\begin{array}{c}\hfill u_x(x,z=0,t)=\{\begin{array}{cc}{v}_0(t)& x<0\\ 0& x=0\\ -v_0(t)& x>0\end{array},\end{array}$$(20)

where

$$\begin{array}{c}\hfill v_0(t)=\frac{1}{2}{v}_{max}[tanh\left(\frac{t-t_0}{w}\right)+1],\end{array}$$(21)

v_max = 1 km s⁻¹, t₀ = 10.1 min, and w = 1.4 min; and the magnetic field B is given by

$$\begin{array}{c}\hfill \mathbf{B}(x,z=0,t)=\frac{F}{\pi }\frac{{\mathbf{r}}_{\mathbf{1}}(t)}{r_1^2(t)}-\frac{F}{\pi }\frac{{\mathbf{r}}_{\mathbf{2}}(t)}{r_2^2(t)}-B_0\widehat{\mathbf{x}},\end{array}$$(22)

where

$$\begin{array}{cc}& {\mathbf{r}}_{\mathbf{1}}(t)=(x+d(t))\widehat{\mathbf{x}}+(z-z_0)\widehat{\mathbf{z}},\hfill \end{array}$$(23)

$$\begin{array}{cc}& {\mathbf{r}}_{\mathbf{2}}(t)=(x-d(t))\widehat{\mathbf{x}}+(z-z_0)\widehat{\mathbf{z}},\hfill \end{array}$$(24)

and

$$\begin{array}{c}\hfill d(t)=d_s-(v_{max}\frac{w}{2}[ln(cosh\left(\frac{t-t_0}{w}\right))-ln(cosh\left(\frac{t_0}{w}\right))]+\frac{1}{2}{v}_{max}t).\end{array}$$(25)

In addition, an absorbing layer was applied on u_z, ρ, and e to ensure that waves hitting the boundaries were not reflected. With respect to the top boundary, we set u_x = 0, B to be line-tied to the flow, and applied an absorbing layer for u_z, ρ, and e.

Because Bifrost is designed to use periodic side-boundaries, we superimposed additional terms to the initial and bottom boundary conditions for B, Eqs. (14) and (22), which corresponds to magnetic sources located in neighbouring domains identical to our computational domain. This adjustment had a negligible effect on the central parts of the domain, where the reconnection takes place, but it ensured that the field was horizontal and ∇ ⋅ B-free at the periodic side-boundaries. For u_x, ρ, and e, we also applied an absorbing layer, thus keeping a periodic side-boundary.

As an additional note regarding the boundaries, the Syntelis-19 and YS-94 resistivity models in this experiment were applied within x ∈ [ − 28, 28] Mm ∧ z ∈ [2, 28] Mm. The resistivity was set uniformly to η₀ outside these regions to avoid conflicts near the boundary layers.

2.3.2. 1D Harris current sheet

Our second experiment was a 1D Harris current sheet that was set up in a computational domain of z ∈ [ − 2, 2] Mm and was discretised over 4096 grid points. To keep this experiment relatively simple, we neglected the gravitational term, ρg, and the Spitzer thermal conductivity term, Q_C, when solving Eqs. (1)–(4). The initial condition for the magnetic field was

$$\begin{array}{c}\hfill \mathbf{B}(z,t=0)=B_{0}tanh((z-z_0)/w)\widehat{\mathbf{x}}.\end{array}$$(26)

When we assume a uniform total pressure (the sum of gas pressure and magnetic pressure) with a uniform temperature T(z, t = 0) = T₀, the initial density is given by

$$\begin{array}{c}\hfill \rho (z,t=0)={\rho }_0+\frac{\mu m_{\mathrm{H}}}{k_{\mathrm{B}}{T}_0}\frac{B_0^2}{\sqrt{8\pi }}(1-{tanh}^2((z-z_0)/w)),\end{array}$$(27)

where ρ₀ is the density far away from the current sheet. In our simulations, we used T₀ = 0.61 MK, ρ₀ = 10⁻¹⁵g cm⁻³, and B₀ = 1 G (as well as w = 20 km and z₀ = 0) in order to approximately match the temperature, mass density, and magnetic field strength in the inflow region of the current sheet of the 2D flux cancellation experiment (Sect. 2.3.1). This ensured that the Alfvén velocity and current density in the 1D and 2D experiment were of the same order of magnitude in the regions near the current sheets, which facilitated performing the same comparisons between the same resistivity models in the two experiments.

The boundary condition was handled by applying an absorbing layer for all variables near the two boundaries to ensure that no waves hitting the boundaries were reflected back into the physical domain. The Syntelis-19 and YS-94 resistivities on this experiment were applied within z ∈ [ − 0.5, 0.5] Mm and were set uniformly to η₀ elsewhere to avoid conflicts near the boundary layers.

3. Results

3.1. 2D flux cancellation experiment

3.1.1. Overview

The simulation was run for 40 min, and the results show that the large-scale evolution of the main quantities such as the magnetic field, temperature, and density agrees relatively well with the case 1 simulation of Syntelis et al. (2019) for the three resistivity models.

The two sources of opposite magnetic polarity located immediately below the photosphere move towards each other with the driving velocity given by Eq. (21) until they meet at x = 0 at t = 40 min. Figure 2 shows that the above-lying photospheric polarities do indeed follow the driver very well throughout the simulation time, until they start to slow down after t = 35 min, similar to Syntelis et al. (2019).

Fig. 2. Evolution of the magnetic polarities in the 2D flux cancellation experiment. The black lines show the horizontal position of each source, given by Eq. (25). The blue lines show the horizontal position of the photospheric polarities, that is, the location along z = 0 where Bz reaches its maximum value. The results here are from the simulation with the Gudiksen-11 resistivity model, but nearly identical results are obtained with the other two resistivity models.

As a consequence of the motion of the photospheric polarities towards each other, the null-point, initially located 7.6 Mm above the photosphere, is stretched into a vertical current sheet with a length of up to ∼0.6 Mm. The reconnection site moves slowly downwards along x = 0 during the cancellation phase, that is, from t = 10 min to t = 40 min. Thermal energy from the reconnection is transported outwards from the current sheet along the magnetic field lines and heats up a wide nearly horizontal open reconnection loop above it and a narrow closed reconnection loop below it. The top panels of Fig. 3 show maps of the temperature in the atmosphere at t = 40 min for each resistivity model. The magnetic field topology is superimposed. The bottom panels show the corresponding maps of the mass density in the region surrounding the null-point¹. The resistivity models are indeed capable of producing a large-scale atmospheric response that agrees among the models, except for some differences in terms of final null-point height and maximum temperature. The height of the elongated null-point (here defined as the centre of the current sheet) at t = 40 min lies at 4.05 Mm above the photosphere in the Syntelis-19 case, 4.0 Mm in the YS-94 case, and 3.85 Mm in the Gudiksen-11 case. The maximum temperature in the heated region at this time is 1.49 MK in the Syntelis-19 case, 1.38 MK in the YS-94 case, and 1.78 MK in the Gudiksen-11 case.

Fig. 3. Atmospheric response in the 2D flux cancellation experiment for the three resistivity models (columns) at t = 40 min. Top: maps of the temperature with the magnetic field topology superimposed. Bottom: maps of the mass density around the reconnection site. A movie of the full evolution from t = 0 to t = 40 min of these maps is available online.

A movie of Fig. 3 is available https://www.aanda.org/10.1051/0004-6361/202346447/olm. It shows the evolution of the temperature, magnetic field, and density throughout the whole simulation time for the three cases. While all cases eventually have temperature profiles of the same structural shape, despite some differences in terms of maximum temperature and null-point height, the plasma inside the current sheet behaves notably differently in each case. In the Syntelis-19 model, the current sheet moves steadily downwards without any sign of plasmoid generation. In the other two resistivity models, plasmoids are generated rapidly. The current sheet in the YS-94 case is different from the other two cases by its remarkably lower mass density. In the Gudiksen-11 case, the current sheet coincides with a thin stripe of increased mass density. This is also visible in the Syntelis-19 case, but to a lesser extent.

The Lundquist number at the centre of the current sheet is ∼5 in the Gudiksen-11 case, ∼10 in the Syntelis-19 case, and ∼20 − 100 in the YS-94 case, while the Reynolds number inside the current sheet approaches unity in all three cases (but it is slightly higher in the YS-94 case). At a horizontal distance of 0.1 Mm from the current sheet, the Reynolds and Lundquist numbers are ∼10⁴ or higher in all three models. This is as expected because the resistivity models were scaled so that the simulation was able to obtain roughly the same Alfvén velocities in the inflow region. The plasma-β inside the current sheet reaches maximum values (in the top and bottom points of the current sheet) of ∼2 − 5 in the Gudiksen-11 case, ∼1 in the Syntelis-19 case, and ∼0.5 in the YS-94 case. At a distance of 0.1 Mm from the current sheet, β ∼ 0.1 in all three cases.

To demonstrate that the three resistivity models work differently on the current sheet, maps of the resistivity along x = 0 as function of height relative to the vertical midpoint of the current sheet and time for each resistivity model are shown in Fig. 4. The dashed lines in each panel mark the top and bottom of the current sheet. The relatively smooth behaviour of the resistivity of the Syntelis-19 model agrees well with the fact that the current sheet in this case evolves steadily without any sign of plasmoid instability. Based on this, it is plausible to expect the current sheet in this case to follow a Petschek-like reconnection scheme, especially in terms of energy conversion, which is analysed in Sect. 3.1.4. The resistivity of the Gudiksen-11 and YS-94 models, on the other hand, varies more rapidly in its magnitude due to the frequent plasmoid generation, and therefore we expect the energy conversion rates in these cases to deviate more significantly from the Petschek theory. While the Gudiksen-11 and Syntelis-19 resistivities inside the current sheet mostly stay within the range of 100–1000 km² s⁻¹, the YS-94 resistivity has a lower average value that reaches below 100 km² s⁻¹ within the boundaries of the current sheet. Along with the fact that the diffusive layer is shorter than in the other cases, this explains why the atmosphere in this case has the lowest maximum temperature: the Joule heating scales directly with the resistivity. Although the diffusive layer in the Gudiksen-11 case is of similar size as in the Syntelis-19 case, the average resistivity of the current sheet in the Gudiksen-11 case is slightly higher because the resistivity is enhanced in the plasmoids that appear relatively frequently. This explains why the atmosphere receives the highest amount of heating in the Gudiksen-11 case.

Fig. 4. Evolution of the resistivity η along the current sheet for each resistivity model. The resistivity is measured at x = 0 and is shown as function of the height relative to the current sheet midpoint, zM. The dashed lines mark the top and bottom of the current sheet, which are annotated as Sh and Sl, respectively, in Fig. 5.

3.1.2. Comparison method

We performed the same comparison between simulations and theory as Syntelis et al. (2019) by locating the current sheet and measuring some inflow values near it and comparing them with values predicted from analytical formulae. To demonstrate the localisation of the current sheet and the regions in which the inflow values are measured, Fig. 5 shows maps of the temperature and inverse characteristic length of the magnetic field, L_B, in the surroundings of the null-point with the inflow region delimited by a rectangle of points A, B, C, and D. We defined the current sheet as the oblong vertical region along x = 0 where the characteristic length for magnetic field, L_b, is shorter than a chosen threshold value of 100 km, which is roughly three grid cells because the numerical resolution of the experiments is ∼30 km. The extremes of the current sheet are indicated in the plots with S_h (top) and S_l (bottom). The corresponding current sheet length L_m is measured as the vertical distance between these two points. The index m denotes that it is a numerically measured value. This indexation was applied to several numerically measured values in order to distinguish them from their analytical counterparts. Points A, B, C, and D are defined such that the AB and CD segments form vertical lines parallel to the current sheet at 0.2 Mm to the left and right of the current sheet, respectively. The choice of this location of the line segments was made so that the segments lay within the range in which the analytical formulae for the inflow values used by Syntelis et al. (2019) are valid. We found that placing AB and CD at any horizontal distance between 0.1 and 0.2 Mm was suitable. We used 0.2 Mm to also be consistent with the criterion employed by Syntelis et al. (2019). The figure shows that the inflow rectangle ABCDA does indeed follow the current sheet as it moves downwards throughout the cancellation phase.

Fig. 5. Evolution of the current sheet in the 2D flux cancellation experiment for the Gudiksen-11 resistivity model at different times (columns). Top: temperature maps zoomed in on the region around the elongated current sheet. Bottom: corresponding maps of the inverse characteristic length $L_B^{-1}$. The rectangle marks the region around the current sheet, and the inflow parameters are measured at the line segments AB and CD. Sh and Sl are the top and bottom of the current sheet, respectively.

The inflow magnetic field strength B_im and velocity v_im were measured as the mean absolute value of the magnetic field and the velocity, respectively, along the line segments AB and CD. The Poynting influx Φ_Sim was measured by integrating the Poynting vector component perpendicular to these line segments, S_x = [E × B]_x/μ₀ = E_yB_z/μ₀, over AB and CD. The average density along AB and CD, ρ_im, was also measured because it is needed in the calculations of the analytical estimate for the Poynting influx.

Knowing the numerical measures B_im, v_im, L_m, and Φ_Sim, we compared them with analytical estimates for B_i, v_i, L, and Φ_Si, as derived by Syntelis et al. (2019). The analytical expression for the inflow magnetic field strength B_i is

$$\begin{array}{cc}\hfill B_i(d,d_0,L)& =B_0\sqrt{\frac{d_0}{d}-1}\frac{L}{d_0},\hfill \end{array}$$(28)

where d and d₀ are the source separation distance and the critical source separation distance, respectively. Two different analytical estimates were made for B_i: 1) B_i(d(t),d₀, L_m), based on the source positions with d(t) given by Eq. (25) and $d_0=\frac{2F}{\pi B_0}$; and 2) B_i(d_m(t),d_0m, L_m), based on the photospheric polarity positions, where d_m(t) is the half-separation distance between the photospheric polarities, shown as the blue curve in Fig. 2, and $d_{0m}=\frac{2F_m}{\pi B_0}$, where F_m = 2 200 G Mm is the flux of each photospheric polarity.

The analytical expression for the inflow velocity is

$$\begin{array}{cc}\hfill v_i(v_0,d_0,L)& =f(d,d_0,z_0)v_0\frac{d_0}{L},\hfill \end{array}$$(29)

where

$$\begin{array}{c}\hfill f(d,d_0,z_0)=1-d_0\frac{z_{max}-z_0}{{(z_{max}-z_0)}^2+d^2}\frac{1}{\sqrt{d_0/d-1}}\end{array}$$(30)

is a flux correction factor, as explained in detail in the appendix of Syntelis et al. (2019), with z_max = 30 Mm as the top of the computational domain. The factor was initially f ≈ 0.72 when d = 1.8 Mm, then approached 1 as d → 0. Again, two analytical estimates were made for the inflow velocity: 1) v_i(v₀(t),d₀, L_m), based on the sources, with v₀(t) given by Eq. (21); and 2) v_i(v_0m(t),d_0m, L_m), based on the photospheric polarities, using f(d_m, d_0m, 0), and where $v_{0m}(t)\equiv {\dot{d}}_m(t)$ is the absolute value of the velocity of the photospheric polarities given by the time derivative of the blue curve in Fig. 2.

The analytical current sheet length is

$$\begin{array}{cc}\hfill L(M_A,d,d_0,v_0,v_{A0})& =\sqrt{f(d,d_0,z_0)}{d}_0\sqrt{\frac{M_{A0}}{M_A}\frac{1}{\sqrt{d_0/d-1}}},\hfill \end{array}$$(31)

where M_A is the inflow Alfvén Mach number, and M_A0 ≡ v₀/v_A0 is a hybrid Alfvén Mach number based on the hybrid Alfvén speed $v_{A0}\equiv B_0/\sqrt{{\mu }_0{\rho }_i}$, a quantity introduced by Syntelis et al. (2019) which is based on the external magnetic field B₀ but the inflow mass density ρ_i (therefore “hybrid”). We estimated 1) L(M_Am, d(t),d₀, v₀(t),v_A0m) based on sources, with $v_{A0m}=B_0/\sqrt{{\mu }_0{\rho }_{\mathit{im}}}$, and 2) L(M_Am, d_m(t),d_0m, v_0m(t),v_A0m) based on photospheric polarities.

The analytical Poynting influx is

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{S_i}(M_A,d,d_0,v_0,v_{A0})& =2f^2(d,d_0,z_0)\frac{v_0{B}_0^2}{{\mu }_0}{d}_0\sqrt{d_0/d-1}\frac{M_{A0}}{M_A},\hfill \end{array}$$(32)

where we estimated 1) Φ_Si(M_Am, d(t),d₀, v₀(t),v_A0m) based on the sources and 2) Φ_Si(M_Am, d_m(t),d_0m, v_0m(t),v_A0m) based on the photospheric polarities.

We also calculated the fractions of the Poynting influx that were converted into kinetic energy and into heat. According to Gauss’ theorem, we have

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{S_i}\equiv {\displaystyle {\oint }_C\frac{1}{{\mu }_0}\mathbf{E}\times \mathbf{B}·d\mathbf{C}={\int }_A\frac{1}{{\mu }_0}\mathrm{\nabla }·(\mathbf{E}\times \mathbf{B})\mathrm{d}A,}\end{array}$$(33)

where C is the curve over the points ABCDA, and A is its enclosed area. This simply states that the energy increase in the system equals the energy added into it. The above equation can, with the help of vector calculus as well as Faraday’s law, Ohm’s law, and Ampère’s law, be rewritten as

$$\begin{array}{c}\hfill \left|{\mathrm{\Phi }}_{S_i}\right|\approx {\displaystyle \left|{\int }_A\eta {\mathbf{J}}^{2}dA\right|+|{\int }_A\mathbf{J}·(\mathbf{v}\times \mathbf{B})\mathrm{d}A|,}\end{array}$$(34)

which indeed tells us that the input magnetic energy is converted into heat (first right-hand-side term) and kinetic energy (second right-hand-side term) through reconnection. A third right-hand-side term, ${\int }_A\frac{\partial }{\partial t}\left(\frac{B^2}{2{\mu }_0}\right)$, was neglected here as Syntelis et al. (2019) did the same (we measured this term in our simulations, and it is indeed small compared to the other right-hand side terms in the above equation). To compare the simulated energy conversion with Petschek (1964) theory, we measured the J ⋅ (v × B) term and the Joule heating term integrated over the rectangle A and compared it to three-fifths and two-fifths of the Poynting influx, respectively. For this comparison, we used both the numerical measure Φ_Sim and the analytical estimate Φ_Si.

3.1.3. Inflow magnetic field, velocity, and current sheet length

Figure 6 shows the comparison between the numerical results (solid lines) and the analytical estimates based on the dynamics of the sources (dashed curves) and the photospheric polarities (dash-dotted curves) for the inflow magnetic field (top panels), the inflow velocity (middle panels), and the current sheet length (bottom panels). The quantities shown in the figure were averaged over 100 s to obtain smooth lines, which reduced their rapid fluctuations as a consequence of the non-stationary nature of the current sheet.

Fig. 6. Evolution of the relevant quantities in the 2D flux cancellation experiment for each resistivity model (columns). Top: inflow magnetic field, both numerical measures (solid curves) and analytical estimates (dashed and dash-dotted curves). Middle: inflow velocity. Bottom: length of the current sheet. The quantities are averaged over 100 s to reduce their rapid fluctuations.

It is clear from the figure that the numerical measures for B_im, v_im, and L_m in each model satisfactorily agree with each other and with the analytical estimates, especially those based on the photospheric response (dash-dotted curves), but they are not identical. The current sheet length in the YS-94 case is slightly shorter than in the other cases, which means that it deviates more strongly from the analytical estimate. The current sheet length in the Syntelis-19 case is similar to that of the Gudiksen-11 case in the first 10 min of the cancellation phase, but it then declines faster. The agreement is best in the Gudiksen-11 model for the numerical measure for L_m and the analytical estimate for L based on photospheric polarities. The inflow velocity in the Syntelis-19 case is more or less the same as in the YS-94 case, both numerically and analytically, while the inflow velocity in the Gudiksen-11 case has a lower maximum value, and the numerical measure and the analytical estimate based on photospheric polarities agree better. The inflow magnetic field in the Gudiksen-11 case has a slightly higher maximum field strength than in the other two cases and simulation and theory agree best, while the field strength in the YS-94 case is weakest and simulation and theory deviate most.

The analytical estimates for L in each model agree very well with each other from t = 15 min and throughout the simulations because the Alfvén Mach number, on which the analytical current sheet length is directly dependent, agrees well. We adjusted the input values of the diffusion scaling parameters of each model (η₁ for the Syntelis-19 model, α for the YS-94 model, and η₃ for the Gudiksen-11 model) on purpose in order to obtain this agreement between the analytical estimates. The analytical estimates for B_i and v_i agree less well when comparing the resistivity models because these estimates depend on the numerical measures for L_m, which are slightly different in each case.

3.1.4. Energy release

Figure 7 shows the energy release in the three models. The quantities here are also averaged over 100 s to reduce their rapid fluctuations. The first row shows the numerical measures of the Poynting influx Φ_Sim (solid line), and the analytical estimates for Φ_Si based on the source positions (dashed curve) and based on the photospheric polarity positions (dash-dotted curve). In all three cases, the numerical measures approach the analytical estimate at t = 13 min, which is approximately the time at which the current sheet length reaches its maximum value. After this time, the numerical Poynting influx stays constant in each case for the next 15 min, instead of increasing, as analytically predicted, before it slowly decreases. These numerical measures roughly follow the same evolution in all the three cases, however, but they reach a slightly lower maximum value in the YS-94 case, and are roughly of same order of magnitude as the analytical estimates based on photospheric polarities.

Fig. 7. Evolution of the energy release in the 2D flux cancellation experiment for each resistivity model (columns). Top: poynting influx, both numerical measures (solid curves) and analytical estimates (dashed and dash-dotted curves). Middle: three-fifths of the numerical measure for the released energy (solid lines), compared to the numerical value for the kinetic energy output (dashed curves) and three-fifths of the analytical estimate for the released energy (dash-dotted curves). Bottom: two-fifths of the numerical measure for the released energy (solid lines), compared to the numerical value for the heat output (dashed curves) and two-fifths of the analytical estimate for the released energy (dash-dotted curves). The quantities are averaged over 100 s to reduce their rapid fluctuations.

The second and third rows show the fraction of the energy that is released through reconnection that is transformed into kinetic energy and thermal energy, respectively, compared to three-fifths and two-fifths, respectively, of the numerical measures and analytical estimates for the Poynting influx. The energy conversion with the Syntelis-19 model is more Petschek-like than with the other two models, with almost exactly three-fifths of the energy input converted into kinetic energy, and slightly less than two-fifths converted into heat. In the Gudiksen-11 model, significantly more than two-fifths of the input energy is converted into heat. It gains more heat than the other two models, and therefore, the agreement between the numerically measured and analytically predicted heat output is best. The YS-94 model deviates most from the Petschek theory: less than one-fifth of the energy is converted into heat. This agrees with Fig. 3, in which the Gudiksen-11 case resulted in the warmest atmosphere. The maximum temperature was almost 0.3 MK higher than in the Syntelis-19 case, while the YS-94 model had the coldest atmosphere with a maximum temperature 0.1 MK lower than in the Syntelis-19 case.

The Syntelis-19 case follows a nearly perfect Petschek-like energy conversion. This agrees with the fact that this simulation has nearly no sign of plasmoid generation in the current sheet, as seen in the movie of Fig. 3. This means that this resistivity model allows the current sheet to undergo Petschek reconnection. In the YS-94 and Gudiksen-11 models, the current sheet undergoes plasmoid-mediated reconnection, which explains why the kinetic and thermal energy released through reconnection is not necessarily equal to three-fifths and two-fifths, respectively, of the input magnetic energy. Still, it is noteworthy that these two cases, while they are plasmoid-mediated, follow completely different energy conversion schemes. While in the Gudiksen-11 case, more of the magnetic energy is converted into heat than predicted with Petschek theory and less into kinetic energy, in the YS-94 case, less magnetic energy is converted into heat and more into kinetic energy. As we described above, this is caused by the significantly stronger diffusive layer in the Gudiksen-11 model than in the YS-94 model, as shown in Fig. 4, where the Gudiksen-11 model clearly has the highest mean resistivity along the centre of the current sheet. The frequency of plasmoids in current sheets as a result of different resistivity models and how this affects the heating of the surrounding plasma will be studied more in detail in an upcoming paper.

3.1.5. Dependence on the choice of diffusion parameters

The results of the above section were obtained by setting the free parameters of the resistivity models to specific values to ensure that the inflow Alfvén speed has roughly the same value in all simulation cases. In this way, we ensured that we solved a very similar physical problem even though we used different numerical approaches. In this section, we study the dependence of the results on an adjustment of these parameters.

For the Gudiksen-11 model (Sect. 2.2.1), we originally used ν₁ = 0.03, ν₂ = 0.2, and η₃ = 0.2. The parameter ν₁ affects the electrical resistivity as well as the viscous terms, and it scales up all the diffusive terms in the MHD equations over the entire computational domain. Therefore, this parameter should be kept as low as possible. It has been shown empirically that ν₁ > 0.02 is needed to obtain stable solutions in several standard test problems to which Bifrost has been applied for a numerical solution (Gudiksen et al. 2011). We studied different choices for this parameter for the 2D flux cancellation experiment and found that ν₁ = 0.03 is a suitable choice because decreasing ν₁ below this value leads to numerical instability in the current sheet, and increasing it much beyond this value will make the whole problem over-diffused.

Furthermore, it has been shown empirically that ν₂ = 0.2 is about the minimum for numerically stable solutions in several standard test problems (Gudiksen et al. 2011). In our case, the length of the current sheet is only slightly affected when this parameter was decreased below that value. However, running the experiment with a higher value of ν₂ led to a reduction of the current sheet length, and therefore, to a considerable deviation between the numerical measures and analytical estimates shown in Figs. 6 and 7.

The only free parameter of the Gudiksen-11 model that is interesting to adjust for our purposes is η₃ because it directly scales the electrical resistivity and has no effect on the viscosity. We tested running the experiment with different values of η₃ and obtained that values below 0.2 are numerically unstable, while values much higher than 0.2 increase the deviation between the numerical measures and the analytical estimates for the inflow values.

The simulation was also run using different values of η₁ for the Syntelis-19 resistivity model. We found that this parameter can be decreased by an order of magnitude from the value used for the results in the above sections without losing numerical stability. However, this reduction of this diffusion parameter causes the current sheet length to be too long compared to the results of Syntelis et al. (2019), thus deviating more from the analytically predicted current sheet length. A further decrease in η₁ will lead to numerical instability. When we instead increase this parameter by an order of magnitude, the current sheet length is too small compared to the analytical estimate. Decreasing the threshold value J_crit has almost the same effect as increasing η₁.

The results obtained with the YS-94 resistivity model seem to be weakly dependent on the scaling parameters: The current sheet length and Poynting influx barely increase when α is decreased by a factor ten. In addition, there is no significant change in the plasmoid behaviour. Decreasing this parameter further causes numerical instability. When the threshold value v_crit is modified, it creates roughly the same effect as adjusting α the opposite way.

For each of the three resistivity models used in this experiment, we observed that the current sheet becomes numerically unstable when the anomalous resistivity is scaled down too strongly. This also shows that the experiment cannot be run without an anomalous resistivity for the given resolution because the current sheet would not be numerically resolvable, unless we were to use a uniform resistivity that is many orders of magnitude greater than the Spitzer resistivity, leading to very unphysical results, or if we were to increase the resolution by several orders of magnitude, causing the experiment to become expensive in terms of compute resources.

3.2. 1D Harris current sheet

In the previous section, we showed that we could use three different resistivity models in a 2D flux cancellation experiment and obtain relatively consistent results in terms of current sheet length and energy release by adjusting the diffusion parameters of each resistivity model. In this section, we begin to study the effects of applying the same resistivity models and parameters to the 1D Harris current sheet experiment introduced in Sect. 2.3.2.

The results of the experiment for the magnetic field B_x, resistivity η, Joule heating Q_J, and temperature T are shown in the first two columns of Fig. 8 at two selected times: one time close to the beginning (0.25 min), and another time at the moment we stopped the simulation (15 min). Even though we applied the same diffusion parameters that ensured relatively consistent results for the 2D flux cancellation experiment, the results for this 1D Harris sheet vary significantly depending on the resistivity model. At t = 0.25 min, the Syntelis-19 model has already had a huge impact in terms of diffusing out the current sheet width and heating up the plasma. The YS-94 model has a significant diffusive effect on the current sheet at t = 15 min, but it is still small compared to the Syntelis-19 model. The Gudiksen-11 model has apparently no diffusive effect on the current sheet with the given values for its free parameters. By fitting the B_x profile to a hyperbolic tangent, tanh(z/w), and finding the width w through the least-squares method, we find that the width of the current sheet, which initially is 20 km, has at t = 15 min increased to 217 km with the Syntelis-19 model and to 30 km with the YS-94 model, but it remained at 20 km with the Gudiksen-11 model. The reason is that the resistivity (second row of the figure) in the Syntelis-19 model is highest: it is up to two orders of magnitude higher than in the YS-94 model. At the end of the simulation, in the Syntelis-19 case, its maximum is ∼15 km² s⁻¹, while for the YS-94 model, it is ∼0.20 km² s⁻¹. The resistivity stays < 0.01 km² s⁻¹ in the Gudiksen-11 model. As a result of this, the Joule heating, as seen in the third row, has a maximum value more than one order of magnitude higher in the Syntelis-19 case than in the YS-94 case at the early stages of the simulation, and then this difference decreases over time as the magnetic field is diffused and the currents are smaller. Since the resistivity is really low for the Gudiksen-11 case, the associated Joule heating in this case is negligible. Consequently, the temperature profile in the current sheet, which is initially uniform with a value of 0.61 MK, has risen to a maximum value above 1.1 MK in the Syntelis-19 case at t = 15 min, but only to 0.69 MK in the YS-94 case. It is unchanged in the Gudiksen-11 case. The large asymmetry seen in the temperature profile for the Syntelis-19 case at t = 15 min is due to the tiny asymmetries in the staggered mesh, which are rapidly magnified by the relatively high diffusivity of this resistivity model (with the given values for the diffusion parameters).

Fig. 8. Evolution of the 1D Harris current sheet. From top to bottom: the magnetic field Bx, the resistivity η, the Joule heating QJ, and the temperature T are plotted as obtained by using the Syntelis-19 (green), YS-94 (blue), and Gudiksen-11 (red) resistivity models. The first and second columns show the results, measured at different times, setting the diffusion parameters to the same values as used in the 2D experiment. The third and fourth columns show the results obtained after adjusting these diffusion parameters to obtain the same behaviour on this 1D Harris sheet for the three resistivity models.

For comparison, the third and fourth columns of Fig. 8 show the results after adjusting the scaling parameter of each resistivity model to ensure that they have roughly the same diffusive effect on this 1D Harris sheet. The new values for the adjusted parameters are η₃ = 1.0 for the Gudiksen-11 model, η₁ = 3.78 × 10⁻³ km² s⁻¹ for the Syntelis-19 model, and α = 2.0  ×  10⁻⁸ km² s⁻¹ for the YS-94 model. With the adjusted values, all three resistivity models diffuse the current sheet out to a final width of ∼26 Mm at t = 15 min. The resistivity at the centre of the current sheet lies at slightly above ∼0.10 km² s⁻¹ in all three cases, causing the final Joule heating profiles to be nearly identical and the final maximum temperature to reach about 0.66 MK in all three cases. One noticeable difference is seen in the resistivity in the regions outside the current sheet, where the magnetic field is nearly constant. The Gudiksen-11 model is nearly an order of magnitude higher than the other two models because the resistivity of this model depends, among other factors, on third derivatives of the magnetic field as well as on the gradients in the velocity perpendicular to the field. This makes it relatively sensitive to tiny perturbations in the current density that are enhanced by the velocity perturbations that arise during the diffusion of the current sheet. However, this enhancement of the resistivity outside the current sheet does not affect the temperature profile at all because the current density, and hence the Joule heating, is here several orders of magnitude lower than at the centre of the current sheet. Additionally, the Lundquist number in the Gudiksen-11 case is above 10⁴ at any distance greater than 0.01 Mm away from the current sheet. This agrees well with the other two resistivity models. This shows indeed that the resistivity outside the current sheet has no effect on the evolution of the plasma.

We have shown that the resistivity models resulted in completely different levels of the diffusive effect when they were applied in this 1D Harris current sheet experiment when the same diffusion parameter values were used that in the 2D flux cancellation experiment gave results that agreed well. We also demonstrated that we can easily adjust the diffusion parameters to obtain roughly the same diffusive behaviour in this relatively simple experiment. The free parameters of the YS-94 and Gudiksen-11 models only needed adjustments within roughly the same order of magnitude to obtain these results, as shown in the second two columns of Fig. 8, but the η₁ value of the Syntelis-19 model needed to be decreased by more than three orders of magnitude. This is due to its direct scaling with the current density, which causes the diffusivity of this resistivity model to be strongly dependent on the magnetic field topology.

4. Discussion

This comparative study of resistivity models has demonstrated that we can use different types of resistivity models in the same numerical experiment and still obtain results that agree relatively well with each other. We successfully mimicked a 2D flux cancellation experiment from Syntelis et al. (2019) and found that using Bifrost’s hyper-diffusive resistivity model (Gudiksen et al. 2011, referred to in this paper as Gudiksen-11) results in a current sheet length that more or less follows the same evolution as when using the current density-proportional resistivity model of the original experiment (Syntelis-19), given the right input values for the diffusion parameters. The magnetic field and velocity measured in the inflow region of the current sheet also develop in a similar way when the experiment is performed with each of these two resistivity models. As a result of this, the Poynting influx evolves similarly in both cases. The energy conversion, on the other hand, follows different schemes in each case. While the energy conversion in the Syntelis-19 case agrees with the Petschek theory, the current sheet in the Gudiksen-11 case undergoes plasmoid-mediated reconnection and a significantly higher portion of the magnetic energy is converted into heat. As a result, the maximum temperature is higher in this last case. The drift velocity-dependent resistivity model (YS-94), previously applied by Yokoyama & Shibata (1994), among others, was also applied for the same experiment. The results obtained when using this resistivity model also agree satisfactorily with the results from the other two resistivity models. The current sheet is slightly shorter and the inflow magnetic field is slightly weaker, however, leading to a significantly lower Poynting influx. Despite undergoing plasmoid-mediated reconnection, a lower portion of the input magnetic energy is converted into heat in this case than in the Petschek-conform Syntelis-19 case, in contrast to the Gudiksen-11 case, in which the conversion rate of magnetic energy to heat is higher. Therefore, the heated region has a lower temperature than in the other two cases. Except for the differences in terms of plasmoid generation and energy conversion, the temperature and mass density profiles of all three cases have a similar structural shape.

Furthermore, we observed that when we numerically solved the same model equations for a 1D Harris current sheet, the results in terms of diffusive rates and Joule heating obtained using each of the three resistivity models were significantly different from each other, given the same input values for the diffusion parameters as in the 2D experiment. Running the same experiment with adjusted values for the diffusion parameters showed that two of these resistivity models, namely Gudiksen-11 and YS-94, needed only adjustments within the same order of magnitude for their scaling parameters in order to obtain the same diffusive rate on the Harris sheet. The scaling parameter η₁ in the Syntelis-19 resistivity model, on the other hand, needed to be scaled down by more than three orders of magnitude from its value applied in the 2D experiment in order to obtain the same diffusive rate in this 1D Harris sheet experiment as the other two resistivity models.

One of the free parameters of the resistivity model used by Syntelis et al. (2019) requires an adjustment of several orders of magnitude when jumping between these two experiments because the resistivity scales linearly with the current density. This causes the ideal value for the scaling parameter to be strongly dependent of the magnetic field topology of the experiment when a satisfactory result is to be obtained. Moreover, its linear proportionality to the current density causes the resistivity to stay relatively high in relatively large areas around the current sheet. The Lundquist number therefore increases relatively slowly with distance from the current sheet compared to the other two resistivity models that were tested in this paper. Finally, because η scales with the current density, the anomalous resistivity in regions near to magnetic sources needed to be turned off. This resistivity model works in a satisfactory way for several numerical experiments when the scaling parameter is adjusted properly, however.

We observed that the electron drift velocity-dependent resistivity model that was previously used by Yokoyama & Shibata (1994) might be used to obtain results in both experiments of this paper that agree satisfactorily with the corresponding results obtained with Bifrost’s hyper-diffusion model without adjusting the scaling parameter drastically. However, both our experiments dealt with coronal plasma with approximately the same temperature and density as well as similar magnetic field strength. The experiment of Yokoyama & Shibata (1994), on the other hand, which used the same resistivity model to handle current sheets in the upper convection zone, required the scaling parameter to be larger by several orders of magnitude. As the typical electron drift velocity and electron thermal velocity (which typically determines the threshold velocity at which this type of anomalous resistivity is to be activated) differs by several orders of magnitude from the upper convection zone to the upper corona, the ideal values for the free parameters of this resistivity model strongly depend on the local plasma conditions. We were therefore also able to activate the anomalous resistivity of this model only in the coronal region of our 2D experiment (as the scaling parameter was set to handle coronal plasmas) and had to apply a relatively low uniform resistivity below. Despite this, we were fully able to use this resistivity model and obtain results in both our experiments that agreed relatively well with the results obtained with the other two resistivity models, after the free parameters were adjusted properly.

The hyper-diffusive resistivity model of Bifrost (Gudiksen et al. 2011), on the other hand, depends not only on the magnitude of magnetic field gradients, but also on the local fast-mode wave velocity, fluid velocity, and velocity gradients along magnetic field lines. This ensures that the resistivity of this model becomes large only when it is really needed to be large in order to make current sheets numerically resolvable and stay relatively low elsewhere. With a default set of input values for the diffusion parameters, this resistivity model can be applied on anything from coronal plasmas to convection zone plasmas with any type of magnetic field topology without adjusting the parameters drastically. Therefore, this resistivity model does not need to be turned off and replaced by uniform resistivity in specific areas of the computational domain, but can rather be applied on the whole domain.

It is important to point out that several simplifications were made in this study, which is only a rough representation of driven reconnection in the solar atmosphere. For a more detailed study of the reconnection in the Sun, partially ionised effects such as ambipolar diffusion (Zweibel 1989) and the Hall effect (Huang et al. 2011) cannot be ignored, especially when studying the energy balance in the chromosphere (Wargnier et al. 2023) and the heating mechanisms for EBs (Liu et al. 2023) and UV bursts (Ni et al. 2022). These effects also play a significant role in the structure of the inflow current density (Snow et al. 2018), plasmoid formation (Singh et al. 2019; Murtas et al. 2021), and reconnection-driven slow-mode shocks (Hillier et al. 2016). A detailed study of the reconnection rate in plasmoid-mediated reconnection may be performed with high-resolution simulations of a 2D current sheet (Bhattacharjee et al. 2009). More realistic studies of the turbulent energy cascade that occurs in flux ropes generated along the current sheets where the reconnection takes place can be made through high-resolution 3D MHD simulations (Dong et al. 2022) or particle-in-cell simulations (Daughton et al. 2011). We acknowledge that the details of the reconnection physics cannot be revealed through MHD models with anomalous resistivity, and this is not what we attempted to achieve with our study. With the simplifications and assumptions that were made, however, we achieved the insight that three relatively different anomalous resistivity models can be applied on a well-known physical problem to obtain results that agree relatively well with each other. The main gain in knowledge with the hyper-diffusive resistivity model of Bifrost from the results of our experiments is that it is not that strongly dependent on local plasma conditions and magnetic field topology and can therefore be applied on the whole solar atmosphere as well as to upper convection zone in numerical models without using different values for the free parameters in different areas of the computational domain.

Movie

Movie 1 associated with Fig. 3 (results_2dflux_atmospheric_response) Access here

Acknowledgments

This research has been supported by the European Research Council through the Synergy Grant number 810218 (“The Whole Sun”, ERC-2018-SyG) and by the Research Council of Norway through its Centres of Excellence scheme, project number 262622. The simulations were performed on resources provided by Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway. The authors are grateful to the referee for his/her constructive comments to improve the paper.

References

All Figures

	Fig. 1. Initial conditions for the 2D flux cancellation experiment mimicking (Syntelis et al. 2019). Top: map of the temperature with the magnetic field topology superimposed. Middle: vertical component of the magnetic field, Bz, at z = 0. Bottom: stratification of the temperature (black) and mass density (red).
In the text

Fig. 2. Evolution of the magnetic polarities in the 2D flux cancellation experiment. The black lines show the horizontal position of each source, given by Eq. (25). The blue lines show the horizontal position of the photospheric polarities, that is, the location along z = 0 where Bz reaches its maximum value. The results here are from the simulation with the Gudiksen-11 resistivity model, but nearly identical results are obtained with the other two resistivity models.

In the text

	Fig. 3. Atmospheric response in the 2D flux cancellation experiment for the three resistivity models (columns) at t = 40 min. Top: maps of the temperature with the magnetic field topology superimposed. Bottom: maps of the mass density around the reconnection site. A movie of the full evolution from t = 0 to t = 40 min of these maps is available online.
In the text

	Fig. 4. Evolution of the resistivity η along the current sheet for each resistivity model. The resistivity is measured at x = 0 and is shown as function of the height relative to the current sheet midpoint, zM. The dashed lines mark the top and bottom of the current sheet, which are annotated as Sh and Sl, respectively, in Fig. 5.
In the text

Fig. 5. Evolution of the current sheet in the 2D flux cancellation experiment for the Gudiksen-11 resistivity model at different times (columns). Top: temperature maps zoomed in on the region around the elongated current sheet. Bottom: corresponding maps of the inverse characteristic length $L_B^{-1}$. The rectangle marks the region around the current sheet, and the inflow parameters are measured at the line segments AB and CD. Sh and Sl are the top and bottom of the current sheet, respectively.

In the text

	Fig. 6. Evolution of the relevant quantities in the 2D flux cancellation experiment for each resistivity model (columns). Top: inflow magnetic field, both numerical measures (solid curves) and analytical estimates (dashed and dash-dotted curves). Middle: inflow velocity. Bottom: length of the current sheet. The quantities are averaged over 100 s to reduce their rapid fluctuations.
In the text

Fig. 7. Evolution of the energy release in the 2D flux cancellation experiment for each resistivity model (columns). Top: poynting influx, both numerical measures (solid curves) and analytical estimates (dashed and dash-dotted curves). Middle: three-fifths of the numerical measure for the released energy (solid lines), compared to the numerical value for the kinetic energy output (dashed curves) and three-fifths of the analytical estimate for the released energy (dash-dotted curves). Bottom: two-fifths of the numerical measure for the released energy (solid lines), compared to the numerical value for the heat output (dashed curves) and two-fifths of the analytical estimate for the released energy (dash-dotted curves). The quantities are averaged over 100 s to reduce their rapid fluctuations.

In the text

Fig. 8. Evolution of the 1D Harris current sheet. From top to bottom: the magnetic field Bx, the resistivity η, the Joule heating QJ, and the temperature T are plotted as obtained by using the Syntelis-19 (green), YS-94 (blue), and Gudiksen-11 (red) resistivity models. The first and second columns show the results, measured at different times, setting the diffusion parameters to the same values as used in the 2D experiment. The third and fourth columns show the results obtained after adjusting these diffusion parameters to obtain the same behaviour on this 1D Harris sheet for the three resistivity models.

In the text