# FINMHD: An Adaptive Finite-element Code for Magnetic Reconnection and Formation of Plasmoid Chains in Magnetohydrodynamics

Hubert Baty

► To cite this version:

Hubert Baty. FINMHD: An Adaptive Finite-element Code for Magnetic Reconnection and Formation of Plasmoid Chains in Magnetohydrodynamics. The Astrophysical Journal Supplement, 2019, 243 (2), pp.23. 10.3847/1538-4365/ab2cd2. hal-02136272

HAL Id: hal-02136272

https://hal.science/hal-02136272v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

# FINMHD: An Adaptive Finite-element Code for Magnetic Reconnection and Formation of Plasmoid Chains in Magnetohydrodynamics

Hubert Baty

Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France
hubert.baty@astro.unistra.fr

Received 2019 April 24; revised 2019 June 21; accepted 2019 June 24; published 2019 July 30

## Abstract

Solving the problem of fast eruptive events in magnetically dominated astrophysical plasmas requires the use of particularly well adapted numerical tools. Indeed, the central mechanism based on magnetic reconnection is determined by a complex behavior with quasi-singular forming current layers enriched by their associated small-scale magnetic islands called plasmoids. A new code is thus presented for the solution of two-dimensional dissipative magnetohydrodynamics (MHD) equations in cartesian geometry specifically developed to this end. A current-vorticity formulation representative of an incompressible model is chosen in order to follow the formation of the current sheets and the ensuing magnetic reconnection process. A finite-element discretization using triangles with quadratic basis functions on an unstructured grid is employed, and implemented via a highly adaptive characteristic-Galerkin scheme. The adaptivity of the code is illustrated on simplified test equations and finally for magnetic reconnection associated with the nonlinear development of the tilt instability between two repelling current channels. Varying the Lundquist number $S$ has allowed us to study the transition between the steady-state Sweet-Parker reconnection regime (for $S \lesssim 10^4$) and the plasmoid-dominated reconnection regime (for $S \gtrsim 10^5$). The implications for the understanding of the mechanism explaining the fast conversion of free magnetic energy in astrophysical environments such as the solar corona are briefly discussed.

**Key words:** magnetic reconnection – magnetohydrodynamics (MHD) – plasmas – stars: coronae – Sun: flares

## 1. Introduction

Magnetic reconnection is a fundamental process in laboratory and space plasmas, which allows the conversion of magnetic energy into bulk flow and heating. The understanding of magnetic reconnection is a major goal of theoretical plasma physics in order to explain explosive events such as solar/stellar flares, coronal mass ejections, and gamma-ray flares in pulsar winds. For collisional plasmas, the magnetohydrodynamic (MHD) model is commonly adopted as an excellent framework (Priest & Forbes 2000). The MHD approximation is a single-fluid description of plasma dynamics, where the electromagnetic properties are taken into account via coupling terms obeying Maxwell's equations. Typically, electrical currents can be generated by via fluid motions producing Lorentz forces on the fluid. The plasma resistivity (inverse of the electrical conductivity) representing the magnetic field dissipation, i.e., $\eta$, is very small in plasmas of interest, but it cannot be neglected because it precisely drives the reconnection process. The classical model of reconnection is based on Sweet-Parker theory in the two-dimensional (2D) resistive MHD framework, in which a steady-state current sheet structure with a small central diffusion layer controls the reconnection between two regions of oppositely directed magnetic fields (Parker 1957; Sweet 1958). However, the Sweet-Parker (hereafter SP) model gives a reconnection rate scaling as $\eta^{1/2}$, which is too small by a few orders of magnitude to explain the fast timescales involved in the previously mentioned eruptive events.

After many years of research on the existence of viable alternative models, it has recently become clear that the problem can be solved by obtaining a new regime of accelerated reconnection when the resistivity is low enough (Loureiro et al. 2007). More explicitly, this is the case in

plasmas having local Lundquist numbers $S$ higher than a critical value $S_c$ that is of order $10^4$, where $S = lV_A/\eta$ is defined by taking $l$ as the characteristic length scale of the current sheet and $V_A$ as the characteristic Alfvén speed based on the local magnetic field just upstream of the layer. As shown by many numerical MHD simulations, this regime is characterized by a current layer broken into a chain of smaller magnetic islands called plasmoids, which are constantly forming, moving, eventually coalescing, and finally being ejected through the outflow boundary. At a given time, the system appears as an aligned layer structure of plasmoids of different sizes, and can be regarded as a statistical steady state with an average reconnection rate that is nearly (or exactly) independent of the resistivity (Bhattacharjee et al. 2009; Samtaney et al. 2009; Huang et al. 2017).

There are, however, some fundamental issues about the growth of the plasmoids and the ensuing reconnection rate that remain unsolved, reflecting our inadequate understanding of the mechanism. For example, a model based on the modal stability analysis of Sweet-Parker forming current sheets, proposed by Comisso et al. (2017), has led to scaling laws for the growth of the plasmoid instability that are not simple power laws of the resistivity (or Lundquist number). Moreover, the fully plasmoid-dominated regime is thus predicted to be obtained when $S$ is higher than another critical Lundquist number $S_T$, called the transitional Lundquist number, which is substantially higher than $S_c$. For example, a value of $S_T \simeq 10^6-10^7$ is estimated for typical solar corona parameters. This transition corresponds to the inequality $\tau_p \lesssim \tau$, where $\tau_p$ and $\tau$ represent the characteristic timescales for the growth of the plasmoids and for the formation of the SP current sheet (that is, for example, an exponentially shrinking process due to some MHD instability) respectively. The latter results are in apparent contradiction with another stability study (Pucci & Velli 2014),
where the linear growth rate of the plasmoids scales as $S^{1/4}$, and where the divergence for an infinite Lundquist number is regularized by considering current sheets as having asymptotic aspect ratios, $l/a \simeq S^{1/3}$ ($a$ being the characteristic width of the current layer), that are smaller than Sweet-Parker ones (i.e., equal to $S^{1/2}$). This second model also predicts constant $\tau_p \simeq l/V_A$ independent of $S$ in the limit of infinite $S$.

Most of the previous numerical studies have used general MHD codes, though not always specifically designed for investigating the magnetic reconnection process. This is the case, for example, for finite-volume-based codes with shock-capturing methods using different Riemann-type solvers in order to handle discontinuities and shocks, such as for example in AMRVAC (see Porth et al. 2014 for code description, and Baty 2017 for an example of application). On the other hand, despite the fact that they are probably the best adapted to such problems, one finds only a few attempts in the literature at using codes based on finite-element methods. The exceptions concern mainly the related topic of the formation of singularities in ideal MHD (Strauss & Longcope 1998; Lankalapalli et al. 2007), but not the reconnection process itself. Conventional codes lack some convergence properties in the low-resistivity regime, which is characterized by a complicated time-dependent bursty dynamics (Keppens et al. 2013). Indeed, even if the global scenario with various evolutionary stages associated with the magnetic reconnection process is generally well captured in the qualitative sense, there is a lack of true convergence properties due to the long-term inherent chaotic dynamics associated with the presence of plasmoids. We thus develop a new code, FINMHD, using a finite-element discretization in order to specifically tackle this problem.

Most of the previous numerical studies on the subject also consider the same initial setup, namely the coalescence instability between two attracting current channels (see Huang et al. 2017 for example). This is a convenient configuration for numerical treatment, although it has the following drawbacks. Indeed, an initial current layer of arbitrary width is generally assumed between the two channels. Second, it is important to use other configurations in order to address the generality of the mechanism, with cases involving the formation of different current sheets and associated timescales. Consequently, we have chosen to consider a different setup, that is the tilt instability between two repelling currents (Richard et al. 1990). Curiously, this latter configuration has been used to test ideal MHD codes (Lankalapalli et al. 2007) or more recently to generate complex magnetic structures to study particle acceleration (Keppens et al. 2014; Ripperda et al. 2017), but not to study magnetic reconnection.

The outline of the paper is as follows. In Section 2, we present our simplified model of 2D incompressible MHD equations written in current-vorticity variables, usually called reduced MHD in the literature. We also define two simplified model equations (advection-diffusion and Burgers equations), which are designed for testing the essential features of our numerical method. This is followed by the description of our characteristic-Galerkin method applied to the simplified equations in Section 3, and the illustration of the efficiency of our adaptive scheme in capturing the formation and evolution of nearly singular structures on three examples in Section 4. In Section 5, we present the FINMHD code and its application to the evolution of the tilt instability. We also investigate different reconnection

regimes by varying the resistivity parameter. Finally, future perspectives using FINMHD, implications for the understanding of the reconnection mechanism in eruptive solar events, and conclusions are reported in Section 6.

## 2. Model Equations

### 2.1. Reduced MHD Model

We first consider the 2D incompressible set of dissipative (viscous and resistive) MHD equations written in flux-vorticity ($\psi-\omega$) scalar variables as follows:

$$ \frac{\partial \omega}{\partial t} + (\mathbf{V} \cdot \nabla)\omega = (\mathbf{B} \cdot \nabla)J + \nu \nabla^2 \omega, \quad (1) $$

$$ \frac{\partial \psi}{\partial t} + (\mathbf{V} \cdot \nabla)\psi = \eta \nabla^2 \psi, \quad (2) $$

$$ \nabla^2 \phi = -\omega, \quad (3) $$

$$ \nabla^2 \psi = -J, \quad (4) $$

where we have introduced the two stream functions, $\phi(x, y)$ and $\psi(x, y)$, defined by $\mathbf{V} = \nabla\phi \wedge e_z$ and $\mathbf{B} = \nabla\psi \wedge e_z$ ($e_z$ being the unit vector perpendicular to the $xOy$ simulation plane). Note that $J$ and $\omega$ are the $z$ components of the current density and vorticity vectors, $\mathbf{J} = \nabla \wedge \mathbf{B}$ and $\omega = \nabla \wedge \mathbf{V}$ respectively (with units using $\mu_0 = 1$). Note that we consider the resistive diffusion via the $\eta \nabla^2 \psi$ term ($\eta$ being assumed uniform for simplicity), and also a viscous term $\nu \nabla^2 \phi$ in a similar way (with $\nu$ being the viscosity parameter). The above definitions results from the choice $\psi \equiv A_z$, where $A_z$ is the $z$ component of the potential vector $\mathbf{A}$ (here $\mathbf{B} = \nabla \wedge \mathbf{A}$). This choice is the one used in Ng et al. (2007) and in Baty & Nishikawa (2016), and different from the one used by Lankalapalli et al. (2007) where the choice $\psi \equiv -A_z$ is made. In the latter case, the two Poisson equations (i.e., Equations (3) and (4)) involve an opposite sign on the right-hand sides. Note that the thermal pressure gradient is naturally absent from our set of equations. The main advantage of the above formulation over a standard one using the velocity and magnetic field vectors ($\mathbf{V}, \mathbf{B}$) as the main variables is the divergence-free property naturally ensured for these two vectors. However, except when spectral/pseudo-spectral methods are employed, the space discretization of the flux-vorticity formulation is known to generate numerical instabilities at cell interfaces due to the evaluation of the third-order spatial derivative term $(\mathbf{B} \cdot \nabla)J$, because $J$ is itself deduced from $\nabla^2 \psi = -J$ (Philip et al. 2007).

Fortunately, another formulation using current-vorticity ($J-\omega$) variables cures the above difficulty. The latter is obtained by taking the Laplacian of the evolution equation for $\psi$, leading to a new equation for $J$ (now the maximum spatial derivative is second order),

$$ \frac{\partial \omega}{\partial t} + (\mathbf{V} \cdot \nabla)\omega = (\mathbf{B} \cdot \nabla)J + \nu \nabla^2 \omega, \quad (5) $$

$$ \frac{\partial J}{\partial t} + (\mathbf{V} \cdot \nabla)J = (\mathbf{B} \cdot \nabla)\omega + \eta \nabla^2 J + g(\phi, \psi), \quad (6) $$

$$ \nabla^2 \phi = -\omega, \quad (7) $$

$$ \nabla^2 \psi = -J, \quad (8) $$
$$ \text{with} \quad g(\phi, \psi) = 2 \left[ \frac{\partial^2 \phi}{\partial x \partial y} \left( \frac{\partial^2 \psi}{\partial x^2} - \frac{\partial^2 \psi}{\partial y^2} \right) - \frac{\partial^2 \psi}{\partial x \partial y} \left( \frac{\partial^2 \phi}{\partial x^2} - \frac{\partial^2 \phi}{\partial y^2} \right) \right]. $$

Note that this second formulation has another advantage over the first one in that the two evolution equations are more symmetric, thus facilitating the numerical matrix calculus (see below) for linear solvers.

## 2.2. Simplified 2D Advection-Diffusion and Viscous Burgers Equations

In order to introduce the basic concepts of our numerical method and test its essential properties, we also consider two simplified 2D time-dependent equations.

First, we consider an advection-diffusion equation for the variable $u(x, y, t)$,

$$ \frac{\partial u}{\partial t} + (V \cdot \nabla)u = \nu \nabla^2 u + f(u), \quad (9) $$

where $V = (V_x, V_y)$ is a given velocity, $\nu$ is the viscosity parameter, and $f(u)$ is a general source term that can be nonlinear.

Second, we consider the following viscous Burgers equation for the variable $u(x, y, t)$:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + V_y \frac{\partial u}{\partial y} = \nu \nabla^2 u + f(u), \quad (10) $$

where $V_y$ is a given velocity along the y direction and $\nu$ is the viscosity parameter. Note that the quasi-linear term involving shocks ($u \frac{\partial u}{\partial x}$) is one-dimensional. This second equation can also be written in a similar way to the previous equation (Equation (9)) with the velocity vector, $V = (u, V_y)$.

# 3. The Adaptive Characteristic-Galerkin Finite-element Method for the Simplified Equations

## 3.1. The Standard Galerkin Method

The standard Galerkin method is obtained by multiplying the relevant equation to be solved by a test function $v_h(x, y)$ defined on a partition $\tau_h$ of the domain $\tau$ (with a mesh of non-overlapping elements of size $h$, which can be triangles for example), and integrating over the whole domain. This leads to the following Galerkin form for the previous test equation:

$$ \left( \frac{\partial u_h}{\partial t} + (V \cdot \nabla)u_h, v_h \right) + \nu(\nabla u_h, \nabla v_h) - (f(u_h), v_h) = 0, \quad (11) $$

where $u_h$ is the weak form solution on $\tau_h$ approximating $u$ on $\tau$, and where we use the definition $(u_h, v_h) = \int_{\tau_h} u_h v_h d\tau$ evaluated as $\sum_{k=1}^N \int_{\tau_k} u_h v_h d\tau$ (N being the number of elements $T_k$ covering the partition $\tau_h$). Note that the Laplacian term is now expressed as a bilinear form, as the result of using the mathematical Green's theorem and assuming homogeneous Dirichlet boundary conditions for simplicity. An additional integral term (in Equation (11)) along the domain boundary must also be added to take into account eventual Neumann conditions, and additional Dirichlet conditions can also be added separately.

The second step in the Galerkin method concerns the choice of the elements (finite-element functions) in order to cover the domain. In this work, the planar domain is triangulated using the Lagrange basis with polynomial functions of order 1 (usually called $P_1$ elements) or of order 2 ($P_2$ elements). The

variable $u_h$ can thus be expanded as $u_h = \sum_{i=1}^M u_i \lambda_i(x, y)$, where $M$ is the dimension of the approximating space and $\lambda_i(x, y)$ are the usual nodal basis functions. Substituting the above expression into the equation, and replacing the test function by each function of the basis, results in an algebraic system for the unknown coefficient $u_i$. More precisely, for the vector $\bar{u}$ containing the coefficients, we have

$$ M \frac{\partial \bar{u}}{\partial t} + C \bar{u} + \nu S \bar{u} - F \bar{u} = 0, \quad (12) $$

where $M, C$, and $S$ are the mass, convection, and stiffness matrices respectively. $F$ can be obtained by linearizing the source function $f$, and is thus a Jacobian matrix. Thus, $M_{i,j} = \int_{\tau_h} \lambda_i \lambda_j d\tau$, $C_{i,j} = \int_{\tau_h} \lambda_i V \nabla \lambda_j d\tau$, and $S_{i,j} = \int_{\tau_h} \nabla \lambda_i \nabla \lambda_j d\tau$. For the Burgers equation, when implementing the system, the $V_x$ component used in the $C$ matrix must be known (the value of $u_h$ at the previous time step can be used for example).

## 3.2. The Characteristic-Galerkin Method

It is well known that standard Galerkin finite-element methods fail for advection and advection-dominated advection-diffusion problems. Among the different stabilization techniques available, we choose the method of characteristics to discretize the Lagrangian derivative $\frac{\partial}{\partial t} + (V \cdot \nabla)$ in the first two equations (Pironneau 1992). The finite-element software Freefem++ allows us to do this in a simple way (Hecht 2012). Indeed, the discretized version of the above Lagrangian derivative operator (for $u$) can be approximated at time $t_{n+1} = t_n + \Delta t$ ($\Delta t$ being the time step) as

$$ \frac{\partial u}{\partial t} + (V \cdot \nabla)u = \frac{u^{n+1}(x, y, t_{n+1}) - u^n \circ (X^n(x, y, t_{n+1}))}{\Delta t}, \quad (13) $$

where $X^n(x, y, t_{n+1})$ is the space location of a particle obtained by integrating backward (in time) along the characteristic curves from a node (at $t_{n+1}$). Thus, $u^n \circ (X^n(x, y, t_{n+1}))$ is the corresponding value obtained for $u$, which can be interpreted as the particular value at $t_n$, $u^{n'}$. In Freefem++, we use a first-order approximation taking $V^n$ to do this integration and thus deduce $u^{n'}$. Moreover, the corresponding time-discretized matrix formulation similar to Equation (6) can be automatically built in FreeFem++ as

$$ \left( \frac{M}{\Delta t} + \nu S - F \right) \bar{u}^{n+1} = \left( \frac{M}{\Delta t} - C^n \right) \bar{u}^n, \quad (14) $$

where the diffusion is discretized fully implicitly, and $C^n$ is the convection operator using the flow velocity $V^n$ (Hecht 2012).

The previous scheme is only to first order in time, and its extension to second order is not trivial (Rui & Tabata 2002). We propose to use the following two-step predictor scheme:

$$ \left( \frac{M}{\Delta t/2} + \nu S - F \right) \bar{u}^* = \left( \frac{M}{\Delta t/2} - C^n \right) \bar{u}^n, \quad (15) $$

$$ \left( \frac{M}{\Delta t} + \frac{\nu}{2} S - \frac{F}{2} \right) \bar{u}^{n+1} = \left( \frac{M}{\Delta t} - \frac{\nu}{2} S + \frac{F}{2} \right) \bar{u}^n - C^n \bar{u}^*, \quad (16) $$
where the predictor step gives a solution $\bar{u}^*$ at $t + \Delta t/2$, and allows us to estimate a velocity value $V$ and a second-order convection term with $C^*$ operator for the corrector step. At the corrector step, the viscous term also needs to be treated with a second-order semi-implicit (Crank–Nicholson-like) discretization in order to get a second-order discretization.

Finally, a linear solver is chosen in order to solve the systems (Equations (14)–(16)). A large choice of direct and iterative solvers are available in Freefem++.

### 3.3. The Adaptive Mesh

In Freefem++, it is possible to adapt the mesh according to a given function by using the Hessian matrix of this function. In this way, the local error is made nearly uniform by choosing a non-uniform mesh. This procedure is based on the a posteriori error estimates, and is shown to be very efficient by choosing the solution at a given time itself for the function (see the tests below).

## 4. Tests on the Simplified Equations

### 4.1. Advection–Diffusion of a Gaussian Bell

As a first example, we consider a Gaussian bell structure defined by

$$u(x, y) = \exp\left(-\frac{(x - x_c)^2 + (y - y_c)^2}{2\sigma^2}\right), \quad (17)$$

of characteristic width $\sigma = 0.1$, and initially centered around $(x_c = 0.5, y_c = 0)$ in a square domain $[-1:1]^2$. The bell is convected (according to Equation (9) with a zero source term $f=0$) with a given velocity field $V = (-y, x)$ (giving a solid rotation in the clockwise direction), and is also diffusing with a viscosity coefficient $\nu$. The exact solution at $t = t_f$ is

$$u(x, t) = \frac{2\sigma^2}{2\sigma^2 + 4\nu t} \exp\left(-\frac{(x - x_f)^2 + (y - y_f)^2}{2\sigma^2 + 4\nu t}\right). \quad (18)$$

We first apply the first-order (in time) characteristic-Galerkin method using a fixed isotropic mesh (see Figure 1) with $P_1$ elements and integrate during one turn (i.e., $t_f = 2\pi$, $x_f = x_c$, $y_f = y_c$), using a rather small viscosity value $\nu = 1.5 \times 10^{-4}$ and a constant time step $\Delta t = 0.023$. The approximate solution on a crude mesh (2116 triangles) calculated at two different times (corresponding to a quarter of a turn and one full turn) is plotted in Figure 2. A much better approximation is obtained when using adapted meshes for a similar total number of triangles. Indeed, this is illustrated in Figure 3, where the same test is done with an isotropic adapted mesh (i.e., with isogeometric triangles corresponding to a total number varying between 3200 and 2200 at the end of the run) and also with a non-isotropic adapted mesh (with the number of triangles varying between 2000 and 1400). The error is very similar for runs with these two adapted meshes, but the CPU time is advantageous when using non-isogeometric elements because the number of triangles is smaller. Note that for this test, the mesh adaptation is done at each time step, using the Hessian matrix of the solution distribution itself. In Figure 4, we have plotted the $L^2$ error norm as a function of the number of triangles $n_t$. The time step $\Delta t$ is also taken to vary like $1/\sqrt{n_t}$ in order to check the linear dependence expected of the $L^2$ error norm on the solution with the mesh size $h$ of the triangles, as we employ the first-order scheme. Moreover, the results obtained by using our

<figure><img src="image_3.png" /><figcaption>Figure 1. Approximate initial bell solution (color contour map) on a uniform crude mesh with 2116 triangles and $P_1$ elements.</figcaption></figure>

adaptive procedure with $P_1$ elements on the same test lead to an error that is substantially lower (divided by two approximately) for the same number of triangles (of order 2000). Note that an isotropic mesh is imposed for this test (see Figure 3), and when a non-isotropic mesh is used the same error is obtained with an even smaller number of triangles (1400 versus 2200). The results are even further ameliorated when using $P_2$ elements instead, as one can see in Figure 4, because the error is divided by a factor of 5 compared to a uniform mesh. A constraint on the maximum edge size $h_{max} = 0.5$ is used, which is not fundamental to our conclusions.

In summary, this means that for a given accuracy, a gain factor of order 100 in the number of elements is obtained thanks to the adaptive scheme with $P_2$ elements in comparison with a uniform mesh calculation with $P_1$ elements. The CPU time is correspondingly much less important. Finally, note that the choice of the solver is not fundamental but is slightly better for CPU time for this test when using a direct one (UMFPACK) rather than an iterative one (GMRES).

### 4.2. Unsteady Advection–Diffusion with Boundary and Internal Layers

We present a challenging test with an unsteady strongly anisotropic solution for the advection equation (Equation (9)) with zero source term $f=0$ (Manzinali et al. 2018). Indeed, we set a condition $u_0(x, y) = 0$ in the domain $[0:1]^2$ except on the boundary, where

$$u_0(x, y) = \begin{cases} 1 & \text{if } x = 0, 0 \le y \le 1, \\ 1 & \text{if } 0 \le x \le 1, y = 1, \\ (\delta - x)/\delta & \text{if } x \le \delta, y = 0, \\ 0 & \text{if } x > \delta, y = 0, \\ (y - 1 + \delta)/\delta & \text{if } x = 1, y \ge 1 - \delta, \\ 0 & \text{if } x = 1, y < 1 - \delta, \end{cases} \quad (19)$$

which is also kept constant during time evolution. A constant velocity field $V = (2, 1)$ is taken with a viscosity coefficient
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Approximate bell solution obtained on a crude uniform mesh with 2116 triangles and $P_1$ elements (see previous figure) for a final time corresponding to a quarter of a turn (left panel) and to a full turn (right panel).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Approximate bell solution obtained for a full turn, using an isotropic adaptive scheme with 2200 (in fact varying between 3200 and 2200) $P_1$ elements (left panel) and a non-isotropic adaptive variant (with the number of elements varying between 1900 and 1400) (right panel). In both cases, the maximum triangle edge size $h_{max}$ is fixed to a value close to 0.5.</figcaption></figure>

$v = 10^{-3}$. We also use $\delta = 7.8125 \times 10^{-3}$. Thus, this problem exhibits boundary layers along $x = 0$ and $y = 1$. A time step $\Delta t = 0.001$ is taken to integrate the equation with our first-order (in time) characteristic-Galerkin method using the efficient scheme with triangles ($P_2$ elements), as described for the previous test.

As time advances, the left boundary layer propagates into the domain, creating an internal layer that finally reaches the right boundary and creates a new boundary layer because of the imposed $u$ value. The initial top boundary layer also reduces progressively in the same time (as illustrated in Figures 5–7). The numerical solution obtained at different times is clearly well captured by our anisotropic adapted mesh, and the quality of our results is similar to that of the best schemes ever developed (see Figures 10–12 in Manzinali et al. 2018).

Indeed, standard Galerkin space discretization would lead to spurious oscillations in such convection-dominated problems, and a fixed mesh would require a prohibitive number of triangles. Our characteristic-Galerkin method with adaptive refined meshes capturing boundary and internal layers is a remedy to this problem. The number of triangles remains relatively low with our scheme, as it varies between 2800 (at the end) and 12,000 (at early time) in this case, with no constraint on isotropy but an arbitrary constraint on the maximum edge size $h_{max} = 0.4$.

### 4.3. Viscous Shock Layer Induced with Burgers Equation

A third (and last) test is proposed in order to check the ability of our method to capture the formation of shock structure. Indeed, we consider now the solution of the viscous Burgers
<figure><img src="image_6.png" /><figcaption>**Figure 4.** $L^2$ error norm obtained for the Gaussian bell rotation (full turn) using uniform meshes as a function of the number of triangles $n_t$ (Unif). A linear variation with $n_t^{0.5}$ is plotted. The $L^2$ errors obtained using the adaptive scheme with isotropic $P_1$ elements (AM P1) and isotropic $P_2$ elements (AM P2) are also reported. Note that the number of triangles and the CPU time are slightly lower for $P_1$ vs. $P_2$. Using a direct solver is also slightly faster.</figcaption></figure>

equation (Equation (10)) with a zero source term $f=0$ and a velocity component $V_y = 1$. Following the initial setup (Chen et al. 2013), we set an initial condition $u_0(x, y) = 1.5 - 2x$ in a $[0:1]^2$ domain, except at the bottom boundary (at $y=0$) where a Dirichlet boundary value $u = 1.5 - 2x$ is taken and kept constant in time. An extremely challenging low viscosity value is also used, $\nu = 1 \times 10^{-5}$. The results obtained with our adaptive method (same scheme as used for the above previous test problem using $P_2$ elements and an adaptive procedure applied at every time step) with a constant time step $\Delta t = 0.004$ are presented in Figures 8 and 9. A shock layer is formed, beginning at $(x=0.75, y=0.5)$ as a consequence of the steepening of the slope solution when propagating upward and rightward, which has been initiated at the bottom layer. As time advances, a final steady-state solution having an oblique shock layer between the previous point and the top-right corner is obtained, as one can see in Figure 8 at $t=11.5$. The final-state structure showing the characteristic curves is nicely reproduced in Figure 9. Depending on time, the number of triangles in our adapted mesh scheme varies between 500 and 17,000 in this case, with no constraint on isotropy but a constraint on the maximum size $h_{\text{max}} = 0.3$.

### 5.1. FINMHD Code

We apply the first-order scheme (Equation (14)) to the MHD equations. This is basically a nonlinear problem, because of the coupling terms in the different equations, that needs to be solved using linear solvers. Different options are possible. We choose the following scheme, where the two (time evolution) equations are solved as a single system as

$$ 
\left\{ 
\begin{aligned}
\frac{\omega^{n+1} - \omega^n \circ X^n}{\Delta t} &= (\mathbf{B}^n \cdot \nabla) J^{n+1} + \nu \nabla^2 \omega^{n+1}, \\
\frac{J^{n+1} - J^n \circ X^n}{\Delta t} &= (\mathbf{B}^n \cdot \nabla) \omega^{n+1} + \eta \nabla^2 J^{n+1} + g(\psi^n, \phi^n),
\end{aligned}
\right. 
\qquad (20) $$

where the two diffusive terms are treated fully implicitly and the nonlinear $g$ term fully explicitly. The two last equations are successively integrated in a separate way, once the solutions for $\omega^{n+1}$ and $J^{n+1}$ are obtained, as

$$ \nabla^2 \phi^{n+1} = -\omega^{n+1}, \qquad (21) $$

$$ \nabla^2 \psi^{n+1} = -J^{n+1}. \qquad (22) $$

The above scheme is only to first order in time, but has been shown to be linearly unconditionally stable (see Appendix A).

A second-order time integrator using a predictor-corrector scheme has also been developed following the idea proposed previously (Equations (15) and (16)):

$$ 
\left\{ 
\begin{aligned}
\frac{\omega^* - \omega^n \circ X^n}{\Delta t / 2} &= (\mathbf{B}^n \cdot \nabla) J^* + \nu \nabla^2 \omega^*, \\
\frac{J^* - J^n \circ X^n}{\Delta t / 2} &= (\mathbf{B}^n \cdot \nabla) \omega^n + \eta \nabla^2 J^* + g(\psi^n, \phi^n),
\end{aligned}
\right. 
\qquad (23) $$

and

$$ \nabla^2 \phi^* = -\omega^*, \qquad (24) $$

$$ \nabla^2 \psi^* = -J^*, \qquad (25) $$

for the predictor step. One must note that the term $\mathbf{B}^n \cdot \nabla$ is evaluated for $\omega^n$ and not for $\omega^*$, because it is known to lead to better stability. The corrector step follows:

$$ 
\left\{ 
\begin{aligned}
\frac{\omega^{n+1} - \omega^n \circ X^*}{\Delta t} &= (\mathbf{B}^* \cdot \nabla) J^* + \frac{\nu}{2}(\nabla^2 \omega^{n+1} + \nabla^2 \omega^n), \\
\frac{J^{n+1} - J^n \circ X^*}{\Delta t} &= (\mathbf{B}^* \cdot \nabla) \omega^* + \frac{\eta}{2}(\nabla^2 J^{n+1} + \nabla^2 J^n) + g(\psi^*, \phi^*),
\end{aligned}
\right. 
\qquad (26) $$

In summary, our scheme nicely reproduces the solution involving a shock despite our scheme not being conservative.

and

$$ \nabla^2 \phi^{n+1} = -\omega^{n+1}, \qquad (27) $$

$$ \nabla^2 \psi^{n+1} = -J^{n+1}. \qquad (28) $$

## 5. FINMHD Code for Tilt Instability and Magnetic Reconnection

Our characteristic-Galerkin method is now applied to the 2D set of reduced MHD equations written in current-vorticity variables (Equations (5)–(8)).

Contrary to the first-order scheme, our second-order variant is linearly conditionally stable. Moreover, the Courant-Friedrichs-Lewy (CFL) limitation (see Appendix A) requires an initial value of the time step that is one order of magnitude smaller than
<figure><img src="image_7.png" /><figcaption>Figure 5. Contour map solution (left panel) of the unsteady advection–diffusion problem (see text) obtained at $t = 0.02$, and corresponding adaptive mesh (right panel).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Same as previous figure but for $t = 0.2$.</figcaption></figure>

the value required for a converged solution obtained by the first-order scheme. This has been numerically checked, leading to the choice of using the first-order scheme in the following magnetic reconnection study.

## 5.2. Magnetic Reconnection using FINMHD

### 5.2.1. Equilibrium Configuration and Numerical Setup for Tilt Instability

The initial magnetic field configuration for tilt instability is a dipole current structure similar to the dipole vortex flow pattern in fluid dynamics (Richard et al. 1990). It consists of two oppositely directed currents embedded in a constant magnetic field (see Figure 10). Contrary to the coalescence instability based on attracting current structures, the two antiparallel currents in the configuration tend to repel. The initial equilibrium is thus defined by taking the following magnetic

flux distribution:

$$ \psi_e(x, y) = \begin{cases} \left(\frac{1}{r} - r\right) \frac{y}{r} & \text{if } r > 1, \\ -\frac{2}{kJ_0(k)} J_1(kr) \frac{y}{r} & \text{if } r \le 1. \end{cases} \qquad (29) $$

And the corresponding current density is

$$ J_e(x, y) = \begin{cases} 0 & \text{if } r > 1, \\ -\frac{2k^2}{kJ_0(k)} J_1(kr) \frac{y}{r} & \text{if } r \le 1, \end{cases} \qquad (30) $$

where $r = \sqrt{x^2 + y^2}$, and $J_0$ and $J_1$ are Bessel functions of order 0 and 1 respectively. Note also that $k$ is the first (non-zero) root of $J_1$, i.e., $k = 3.83170597$. This initial setup is similar to the one used in the previously cited references (Richard et al. 1990; Lankalapalli et al. 2007), and rotated by
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as previous figure but for $t = 0.6$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Contour map solution of the viscous shock layer problem (see text) obtained at two times, $t = 1.5$ (left panel) and the final time $t = 11.5$ (approximate steady-state solution, right panel).</figcaption></figure>

an angle of $\pi/2$ compared to the equilibrium chosen in the other studies (Keppens et al. 2014; Ripperda et al. 2017). Note that the asymptotic (at large $r$) magnetic field strength is unity, and thus defines our normalization. Consequently, our unit time in the rest of this paper will be defined as the Alfvén transit time across the unit distance (i.e., the initial characteristic length scale of the dipole structure). In the usual MHD framework using the flow velocity and magnetic variables, force-free equilibria using an additional vertical (perpendicular to the $x-y$ plane) can be considered (Richard et al. 1990), or non-force-free equilibria can be also ensured through a thermal pressure gradient balancing the Lorentz force (Keppens et al. 2014). In our incompressible reduced MHD model, as thermal pressure is naturally absent, we are not concerned by such a choice.

A stability analysis in the reduced MHD approximation using the energy principle has shown that the linear eigenfunction of the tilt mode is a combination of rotation and outward displacement (Richard et al. 1990). Instead of imposing such a function in order to perturb the initial setup, we have chosen to let the instability develop from the numerical noise. Consequently, an initial zero stream function is assumed, $\psi_e(x, y) = 0$, with zero initial vorticity $\omega_e(x, y) = 0$. Following most of the previously cited studies, this setup and the following computations are done in a square domain $[-3:3]^2$. The boundaries are taken to be relatively far from the center, in order to have a weak effect on the central dynamics. The initial current structure and a few magnetic field lines are reported in Figure 10, with the initial mesh that is adapted using the Hessian matrix of the current density distribution. This arbitrary choice of current density is justified by the fact that, first, the initial dynamics of the tilt instability is
<figure><img src="image_11.png" /><figcaption>Figure 9. Final steady-state solution (see previous figure) using 36 isocontour lines (left panel) and the corresponding mesh grid (right panel).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. Initial configuration for the tilt instability showing the structure of the dipole current density overlaid with magnetic field lines (left panel) and overlaid with the initial grid using the current density to adapt the mesh (right panel).</figcaption></figure>

driven by the current distribution (i.e., this is an ideal current-driven MHD instability), and second, the ensuing magnetic process is controlled by the structure of the current layers.

### 5.2.2. Typical Evolution of the Tilt Instability and Sweet-Parker Reconnection

First, we focus on moderately low values of the resistivity and viscosity, and for simplicity we also assume a fixed Prandtl number, $P_r = \nu/\eta = 1$. The early time evolution of the system corresponding to the tilt instability is well documented (see the previously cited studies). It corresponds to the linear stability analysis, where the pair of oppositely directed currents tend to repel one another, giving rise to a rotation (see Figure 11 for a case with $\eta = 0.0025$) during the linear phase. The sense of rotation (clockwise in Figure 11) is not predetermined and depends only on the numerical noise. This rotation causes two

new regions of enhanced current density at the leading edges of the vortices because of an associated outward component of the linear displacement, taking the form of two bananas. The time history of the maximum current density amplitude (taken over the whole domain) is shown to increase exponentially in time over the equilibrium value (which is equal to approximately 10), as illustrated in Figure 12 for three different runs using $\eta = 0.0025$, 0.001, and 0.0005. We have evaluated the corresponding local Lundquist number at saturation, $S = IV_A/\eta$, where $l$ is the full length of each current layer and $V_A$ is the Alfvén velocity based on the magnetic field measured at the front of the layer. For example, we have estimated that $l = 1.72$ and $V_A = 1.55$ for the case $\eta = 2.5 \times 10^{-3}$, leading to $S \approx 1080$. The values of $l$ and $V_A$ are checked to slightly increase when the resistivity is decreased for $\eta \gtrsim 0.0005$, and they become approximately constant for smaller resistivity values.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Current density at different times during the development of the tilt instability for $\eta = 0.0025$: $t = 6.95$ (left panel) and $t = 8.35$ (right panel).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Time history of the maximum current density amplitude for three resistivity values (i.e., $\eta = 0.0025$, 0.001, and 0.0005) corresponding to the estimated values of the Lundquist number, $S \simeq 1080$, 3400, and 8000 respectively. The horizontal lines indicate the values for the average saturated current densities.</figcaption></figure>

A steady state is obtained with a structure of nearly constant current density subsequently forming, thus driving the reconnection process. For example, an average value of 43 is evaluated for the saturated current density in a run using $\eta = 2.5 \times 10^{-3}$. The choice of the optimal parameters (the initial time step and the adaptive criterion for its current value) necessary for the convergence of the results is reported in Appendix B. For example, the timescale $\tau = (2.6)^{-1} \simeq 0.38$ from the growth rate of the current density is shown to be well reproduced, independently of the resistivity because the tilt mode is an ideal instability. The oscillations in the current density around the average value are due to the sloshing phenomenon, as described in the problem of coalescence of two magnetic islands (Knoll & Chacón 2006), because of the effect of buildup of magnetic pressure in thin current sheets. The detailed structure of the current sheets, the adapted mesh, and corresponding magnetic field lines are visible in Figure 13. Indeed, one can clearly see in the right panel of Figure 13 the few tens of elements covering the width of the current layers. The beginning of the magnetic reconnection mechanism is also well illustrated in the left panel of Figure 13 for a time

corresponding to the early saturation of the current density. The characteristic reconnection time $\tau_r$ can be determined in different equivalent ways. The simplest one is to take the whole duration of the process, and this corresponds to the time elapsed between the saturation (first peak) and the final time when the current density in Figure 12 returns to a very small value corresponding to a new state (free of closed magnetic field lines). A second way is to measure the level of maximum current amplitude $J_{\text{max}}$ during reconnection (horizontal lines in Figure 12), and to multiply it by the resistivity, because $\eta J_{\text{max}}$ is a measure of the reconnection rate (inverse of $\tau_r$) for a steady-state process. We have checked that the reconnection rate for $S \lesssim 10^4$ closely follows an SP scaling, as $J_{\text{max}} \propto S^{1/2}$ (see next subsection).

In the run employing $\eta = 0.0005$, i.e., corresponding to a local Lundquist number $S \simeq 8000$, after two oscillations of the maximum current density at saturation (i.e., $t \simeq 10.7$), a single plasmoid has been observed to form in a transient way (see Figure 14). This plasmoid does not influence the value of the saturated current density and the resulting reconnection rate. This indicates the proximity of the critical Lundquist number as $S_c \simeq 10^4$. The exact value of the critical Lundquist number depends on the numerical noise (Huang et al. 2017).

### 5.2.3. Plasmoid Reconnection Regime

We explore lower resistivity values in the range $[5 \times 10^{-4}: 1 \times 10^{-4}]$, corresponding typically to $S$ values in the range [8000: 40,000]. The time history of the maximum current density is plotted in Figure 15 for four different $S$ values. First, we note that the formation of plasmoids occurs sooner when $S$ is higher. For example, for $S \simeq 13,300$, plasmoids are seen to form close to the first peak ($t \simeq 9$). However, they can appear much earlier for $S \simeq 40,000$ (i.e., at $t \simeq 8.3$). In the latter case, the two current sheets are invaded by a number of plasmoids varying between four and six, in rough agreement with previous studies showing that the number of plasmoids is an increasing function of $S$. This exact number fluctuates during the magnetic reconnection process because new plasmoids are constantly forming, moving, eventually coalescing (giving monster plasmoids), and finally being ejected through the outflow boundaries. These fluctuations also reflect the time oscillations of the maximum current density seen in Figure 15. At a given time, the system thus appears as an aligned layer structure of plasmoids of different sizes, as illustrated in Figure 16.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Zoom-in of current density at saturation ($t = 8.35$) of the tilt instability for $\eta = 0.0025$ (i.e., $S \simeq 1080$), overlaid with magnetic field lines (left panel) and overlaid with the adapted mesh (right panel).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Zoom-in of current density (left current layer) for $S \simeq 8000$ at a relatively late time $t = 10.7$ during magnetic reconnection (see Figure 12) showing the transient formation of one plasmoid (left panel) and the corresponding adapted mesh (right panel).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Time history of the maximum current density for four relatively high Lundquist numbers (cases involving the formation of plasmoids). The Sweet-Parker value expected at saturation for the maximum current density value (i.e., $S^{1/2}$) for the case $S \simeq 40,000$ is plotted for comparison. The letter $P$ indicates the time of disruption of the current sheet due to the plasmoid formation.</figcaption></figure>

<caption>**Figure 16.** Zoom-in of current density at saturation of the tilt instability for $S \simeq 40,000$, overlaid with magnetic field lines. Note that saturated values in the range $[-160:160]$ are used in order to facilitate the visualization of the current structures.</caption>

A characteristic timescale for the growth of plasmoids $\tau_p$ can be deduced from the slopes in Figure 15 for the maximum amplitude of the current density, obtained just before the quasi-steady state involving fluctuations around an average value. This phase also corresponds to a time interval, beginning when the expected SP value is attained and finishing when the plasmoid regime is fully obtained. We have checked that the beginning also corresponds to the early appearance of the

plasmoids (letter $P$ in Figure 15). This timescale clearly decreases with $S$. For example, for the cases with $S \simeq 13,300$ and $S \simeq 20,000$, it is much longer than the characteristic timescale for the growth of the tilt instability (see the slope at a time around $t \simeq 7.5$), which is given by $\tau \simeq 0.38$. This is not the case for $S \simeq 80,000$, because the two timescales now have
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Maximum amplitude at saturation obtained for different Lundquist numbers. Sweet-Parker (SP) and linear scalings are also indicated for comparison.</figcaption></figure>

comparable values $\tau_p \simeq \tau = 0.38$. This scaling law dependence of $\tau_p$ decreasing with $S$ for $S \gtrsim S_c$ agrees with previous numerical studies and also with theoretical models where the growth rate of plasmons for $S \gtrsim S_c$ increases with $S$ (Pucci & Velli 2014; Comisso et al. 2017).

However, another critical Lundquist value ($S_T$) is predicted by Comisso et al. (2017) corresponding to a regime for which $\tau_p \ll \tau$, in disagreement with what is expected from the other model (Pucci & Velli 2014), for which $\tau_p$ remains of the order of the characteristic Alfvén time (i.e., $\tau_p \simeq 1$ in our units) in the limit of high Lundquist number. This important point about the existence of this transitional Lundquist number ($S_T$) will be explored in a future study, because higher $S$ values are required. Indeed, a value of $S_T \simeq 10^6-10^7$ is expected from Comisso et al. (2017).

A significant difference between the two cited models concerns the critical aspect ratio $l/a$ (where $a$ is the width of the current sheet that is thinning due to the evolution of the tilt mode) at which the plasmons begin to disrupt the current layer. Indeed, it is argued that in the limit of infinite $S$ values, a critical value of the aspect ratio $l/a = S^{1/3}$ significantly smaller than the Sweet-Parker one is expected (Pucci & Velli 2014). The investigation of this latter point also requires the use of higher $S$ values.

Finally, we have plotted in Figure 17 the maximum current density obtained at saturation as a function of $S$ for all the runs. The results clearly show a transition between two regimes at a critical Lundquist number $S_c \simeq 10^4$. Indeed, the values for $S \lesssim 10^4$ perfectly follow a Sweet-Parker scaling as $1.3 \times S^{1/2}$, while another linear scaling with $S$ is obtained for $S \gtrsim 10^4$. As concerns the reconnection rate, a constant reconnection rate of $\eta J_{\text{max}} \simeq 0.05$ (independently of $S$) is obtained in the plasmoid regime. This is somewhat surprising because the fully plasmoid-dominated regime is probably not yet obtained for our moderately high $S$ values. The normalized value $\eta J_{\text{max}}/(V_A B_u)$ (where $B_u$ is the local magnetic field amplitude measured upstream of the current sheet), is evaluated to be close to 0.012, in rather good agreement with the value expected from the formula $10^{-2}(1+P_r)^{-1/2}$ (Comisso & Bhattacharjee 2016). This value also agrees with previously reported values of order 0.01 in MHD simulations using a coalescence instability setup (Bhattacharjee et al. 2009; Samtaney et al. 2009).

## 6. Summary, Outlook, and Astrophysical Relevance

In this work, we have presented a new MHD code specifically designed to investigate magnetic reconnection during the formation of quasi-singular current layers. In particular, the focus is on the regime of accelerated magnetic reconnection that is characterized by disruption of the current sheets and the formation of chains of plasmons. A 2D current-vorticity formulation is adopted, which is particularly suitable for our numerical approach. A finite-element method based on the use of quadratic finite-element basis functions is adopted, and is implemented via a characteristic-Galerkin scheme. The complex time-dependent structures are captured by the use of a highly adaptive grid, based on the Hessian matrix of the solution (which is taken to be the current density for magnetic reconnection).

We can summarize our results as follows. First, we have applied our scheme to simplified test problems involving advection-diffusion and Burgers equations with small dissipation coefficients. Indeed, our approach has been demonstrated to be particularly efficient at reproducing the solutions of these numerically challenging problems, where time-dependent quasi-singular structures can form. Second, we have applied our scheme to the 2D set of MHD equations. The code, called FINMHD, has enabled us to follow the development of the tilt instability in a cartesian geometry. Contrary to the widely used coalescence instability, the tilt mode produces in a more self-consistent way two quasi-singular current sheets in the domain. Indeed, the standard setup configuration taken to explore the coalescence problem has a tangential discontinuity, which is usually smoothed in order to avoid numerical difficulties. Thus, an equilibrium having an initial current layer (situated between two attracting current channels) with an arbitrary width is thus chosen (Huang & Bhattacharjee 2010; Huang et al. 2017). In the present problem, such an initial layer is not needed, because the two current sheets are naturally produced as the result of the tilt mode displacement. During the associated shrinking process, an expected steady-state Sweet-Parker reconnection is driven if the Lundquist number $S$ of each layer is below a critical value $S_c \simeq 10^4$. On the other hand, when $S \gtrsim S_c$, the reconnection is accelerated with a reconnection rate that is nearly independent of $S$. In the latter case, the current layers are invaded by plasmons, the number of which is an increasing function of $S$. The maximum amplitude of the current density exhibits some oscillatory behavior with time with an average value $J_{\text{max}}$ that scales as $S$ (contrary to the SP regime where it scales as $S^{1/2}$). The results obtained for $10^4 \lesssim S \lesssim 10^5$ indicate that plasmons can grow on a timescale $\tau_p$ that becomes of the order of the characteristic time for formation of a current sheet $\tau$ ($\tau \simeq 0.38$ for the tilt instability) when $S$ is high enough. The SP aspect ratio of the current sheet ($l/a \simeq S^{1/2}$) is also approximately reached before the disruption time, in agreement with expectations from the theoretical model of Comisso et al. (2017).

As an indication of the performance of our code, a rather modest maximum value of $7 \times 10^4$ triangles was employed to follow the tilt problem for the case with lowest resistivity ($S = 4 \times 10^4$) in our study, to be compared to approximately $10^8$ grid points needed for a traditional finite-volume approach with uniform mesh. Consequently, the use of our code should be particularly promising for investigating magnetic reconnection with high Lundquist number. Future studies employing higher $S$ values are clearly needed in order to investigate the
existence of a second critical Lundquist number, $S_T$, which is expected when the timescale for plasmoid growth $\tau_p$ is much smaller than $\tau$. This is important in order to discriminate between the two previously cited theoretical models. We hope to attain high enough $S$ values of order ~$10^7$ in future studies. The adaptation of the time step using the maximum peak amplitude for the criterion (see Appendix B) should be further improved, because very high peak current amplitudes of order $10^5$ are expected for $S \sim 10^7$, thus leading to very small time steps. Note that the maximum value of current density obtained in this study is only 600 (see Figure 12) and is similar to values obtained previously in conventional MHD codes (Keppens et al. 2014; Ripperda et al. 2017). Moreover, in order to improve the efficiency of the shock-capturing technique, a next step could be the use of a discontinuous Galerkin scheme, following for example the conservative method developed for relativistic magnetic reconnection with high Lundquist number (Zanotti & Dumbser 2011). It is also important to use a wide range of different setups, such as coalescence, tilt, and other instabilities, in order to investigate the dependence on the shrinking process (including the timescale) of the forming current layers.

Such studies are of great relevance in order to understand the mechanism at work so as to explain the fast conversion of free magnetic energy in astrophysical environments such as the solar corona. Given typical values of solar loop parameters in active regions, it is commonly admitted that the amount of available magnetic energy stored is sufficient to feed large flare events. However, observations require two fast timescales, which are the characteristic timescale of the whole process (with typical values of the order of a few minutes) and an even faster timescale representative of an initial explosive-like phase. The reconnection rate of order $0.01V_A B_u$, as obtained in the plasmoid-dominated reconnection regime, is known to give roughly the correct value for the first timescale. On the other hand, the second (explosive) timescale, which is generally called the onset phase, is not clearly identified. The reconnection regime associated with the growth of plasmons, where the reconnection rate increases sharply on the explosive timescale $\tau_p$ (of the order of seconds) would thus account for observations.

The author thanks I. Moufid for fruitful work during his internship that helped improve the choice of the numerical parameters in the FINMHD code. The author also sincerely thanks the anonymous referee for useful suggestions that helped improve the content of the paper.

## Appendix A
Stability Analysis of the Schemes

We consider the stability of a linearized version of Equations (5) and (6), where for the sake of simplicity the linearization is obtained around a uniform equilibrium and the diffusion terms are dropped. In this way, the resulting linear system can be cast in the simple wave-like form

$$ \left\{ \begin{aligned} \frac{\partial \omega_1}{\partial t} &= (\mathbf{B}_0 \cdot \nabla) J_1, \\ \frac{\partial J_1}{\partial t} &= (\mathbf{B}_0 \cdot \nabla) \omega_1, \end{aligned} \right. \qquad (31) $$

where $\mathbf{B}_0$ is the equilibrium magnetic field vector, and $\omega_1$ and $J_1$ are the perturbed vorticity and current density respectively.

Indeed, the previous system is equivalent to

$$ \left\{ \begin{aligned} \frac{\partial^2 J_1}{\partial t^2} &= V_A^2 \nabla^2 J_1, \\ \frac{\partial^2 \omega_1}{\partial t^2} &= V_A^2 \nabla^2 \omega_1, \end{aligned} \right. \qquad (32) $$

where $V_A = B_0$ is the normalized Alfvén speed.

The first-order implicit time-discretized form (following Equation (20)) is

$$ \left\{ \begin{aligned} \frac{\omega_1^{n+1} - \omega_1^n}{\Delta t} &= (\mathbf{B}_0 \cdot \nabla) J_1^{n+1}, \\ \frac{J_1^{n+1} - J_1^n}{\Delta t} &= (\mathbf{B}_0 \cdot \nabla) \omega_1^{n+1}, \end{aligned} \right. \qquad (33) $$

leading to the following amplification matrix,

$$ T_1 = \frac{1}{1 + \alpha^2} \begin{bmatrix} 1 & i\alpha \\ i\alpha & 1 \end{bmatrix}, $$

where $\alpha = kV_A \Delta t$ for a given wavenumber $k$. The associated (complex conjugate) eigenvalues are thus $(1 \pm i\alpha)/(1 + \alpha^2)$ with a modulus equal to $1/(1 + \alpha^2)^{1/2}$ always smaller than one. The scheme is consequently unconditionally stable, and the spectrum of Alfvén waves is more or less damped depending on $\alpha$. For this scheme, the choice of the time step is only determined by the precision required on the linear growth rate of the tilt instability (which is of order unity). Typically, we use a time step $\Delta t = 10^{-2}$, giving a very small second-order damping factor of order $10^{-4}$ on one time step.

The second-order semi-implicit version (following Equations (23) and (26)) is

$$ \left\{ \begin{aligned} \frac{\omega_1^* - \omega_1^n}{\Delta t/2} &= (\mathbf{B}_0 \cdot \nabla) J_1^* \\ \frac{J_1^* - J_1^n}{\Delta t/2} &= (\mathbf{B}_0 \cdot \nabla) \omega_1^n, \end{aligned} \right. \qquad (34) $$

for the predictor step, and

$$ \left\{ \begin{aligned} \frac{\omega_1^{n+1} - \omega_1^n}{\Delta t} &= (\mathbf{B}_0 \cdot \nabla) J_1^* \\ \frac{J_1^{n+1} - J_1^n}{\Delta t} &= (\mathbf{B}_0 \cdot \nabla) \omega_1^*, \end{aligned} \right. \qquad (35) $$

for the corrector step. For this second scheme, the amplification matrix is

$$ T_2 = \begin{bmatrix} 1 - \alpha^2/2 & i\alpha \\ i\alpha(1 - \alpha^2/4) & 1 - \alpha^2/2 \end{bmatrix}. $$

The associated eigenvalues are thus $(1 - \alpha^2/2) \pm i\alpha(4 - \alpha^2)/2$ with a modulus equal to 1. This second scheme is conditionally stable, as $\alpha < 2$ is required. This translates into a necessary CFL condition for stability, $\Delta t < 2/k_M$, where $k_M$ is the maximum numerical wavenumber included in the simulation (typically $2\pi/h_{min}$, with $h_{min}$ being the smallest triangle length scale). For the simulations presented in this work, we have checked that stability is effectively obtained when $\Delta t \lesssim 10^{-3}$.
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Maximum current density amplitude $J_{max}$ (left panel) and maximum vorticity amplitude (right panel) as functions of time for five runs using different values of the initial time step $\Delta t_0$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Maximum current density amplitude measured at the first peak (see previous figure) obtained for different runs as a function of the inverse initial time step $1/\Delta t_0$.</figcaption></figure>

## Appendix B
Convergence of FINMHD with Time Step

The system of reduced MHD equations (Equations (5)–(8)) are integrated using our first-order (in time) characteristic-Galerkin scheme (Equations (20)–(22)) for the tilt instability problem using $\eta = \nu = 0.0025$. Note that $P_2$ elements are necessary in this case, because the maximum spatial derivative in the MHD equations is of second order. In order to save some CPU time, the adaptive mesh procedure is also done only every five time steps. This scheme is checked to be linearly unconditionally stable (see Appendix A). We have employed different time steps in different runs in order to study the convergence of the code. The results are shown in Figure 18 for the maximum amplitude of the current density ($J_{max}$) and vorticity as functions of time, for five different initial time steps $\Delta t_0$. The convergence is shown to be obtained for $\Delta t_0 \simeq 0.01$. We recall that in all runs, the time step is not constant in time, but rather is adapted following $\Delta t = \Delta t_0 J_e / J_{max}$, where $J_e$ and $J_{max}$ are the maximum current density at equilibrium and at current time respectively. This is also confirmed in Figure 19, where is plotted the first peak amplitude of $J_{max}$ as a function of the inverse initial time step $1/\Delta t_0$.

A similar (slightly faster) convergence can be obtained by using the second-order scheme (Equations (23)–(28)), although

that is checked to be only conditionally stable in agreement with the CFL criterion previously derived (see Appendix A). For the simulations presented in this work, we prefer to use the first-order scheme, because we found the stability limitation to be too restrictive (an initial minimum length scale of $10^{-2}$ leading to a maximum initial time step of order $10^{-3}$).

## ORCID iDs

Hubert Baty https://orcid.org/0000-0003-1925-3983

## References

Baty, H. 2017, ApJ, 837, 74

Baty, H., & Nishikawa, H. 2016, MNRAS, 459, 624

Bhattacharjee, A., Huang, Y.-M., Yang, H., & Rogers, B. 2009, PhPl, 16, 112102

Chen, W., Ching-Shan, C., & Chiu-Yen, K. 2013, JCoPh, 234, 452

Comisso, L., & Bhattacharjee, A. 2016, JPIPh, 82, 595820601

Comisso, L., Lingam, M., Huang, Y. M., & Bhattacharjee, A. 2017, ApJ, 850, 142

Hecht, F. 2012, J. Numerical Math., 20, 251

Huang, Y. M., & Bhattacharjee, A. 2010, PhPl, 17, 062104

Huang, Y. M., Comisso, L., & Bhattacharjee, A. 2017, ApJ, 849, 75

Keppens, R., Porth, O., Galsgaard, K., et al. 2013, PhPl, 20, 092109

Keppens, R., Porth, O., & Xia, C. 2014, ApJ, 795, 77

Knoll, D. A., & Chacón, L. 2006, PhPl, 13, 032307

Lankalapalli, S., Flaherty, J. E., Shephard, M. S., & Strauss, H. R. 2007, JChPh, 225, 363

Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007, PhPl, 14, 100703

Manzinali, G., Hachem, E., & Mesri, Y. 2018, CMAME, 340, 864

Ng, C. W., Rosenberg, D., Germaschewski, K., Pouquet, A., & Bhattacharjee, A. 2007, ApJS, 177, 613

Parker, E. N. 1957, JGR, 62, 509

Philip, B., Pernice, M., & Chacon, L. 2007, JCoPh, 227, 8855

Pironneau, O. 1992, CMAME, 100, 117

Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4

Priest, E. R., & Forbes, T. G. 2000, Magnetic Reconnection (Cambridge: Cambridge Univ. Press)

Pucci, F., & Velli, M. 2014, ApJL, 780, L19

Richard, R. L., Sydora, R. D., & Ashour-Abdalla, M. 1990, PhFIB, 2, 488

Ripperda, B., Porth, O., Xia, C., & Keppens, R. 2017, MNRAS, 467, 3279

Rui, H., & Tabata, M. 2002, NuMat, 92, 161

Samtaney, R., Loureiro, N. F., Uzdensky, D. A., Schekochihin, A. A., & Cowley, S. C. 2009, PhRvL, 103, 1050041

Strauss, H. R., & Longcope, D. W. 1998, JChPh, 147, 318

Sweet, P. A. 1958, in IAU Symp. 6, Electromagnetic Phenomena in Cosmical Physics, ed. B. Lehnert (Cambridge: Cambridge Univ. Press), 123

Zanotti, O., & Dumbser, M. 2011, MNRAS, 418, 1004