A resistive extension for ideal MHD

TL;DR

This paper introduces a resistive extension to ideal relativistic MHD equations that captures resistivity effects without fine-tuning, improves computational efficiency, and can be broadly applied to various dynamical systems.

Contribution

The authors develop a first-principles resistive extension to ideal MHD that avoids stiffness issues and enhances simulation speed without parameter tuning.

Findings

Emulates fully resistive MHD behaviour across initial data sets

Reduces computational stiffness, enabling faster explicit evolution

Applicable to a wide range of dynamical systems

Abstract

We present an extension to the special relativistic, ideal magnetohydrodynamics (MHD) equations, designed to capture effects due to resistivity. The extension takes the simple form of an additional source term which, when implemented numerically, is shown to emulate the behaviour produced by a fully resistive MHD description for a range of initial data. The extension is developed from first principle arguments, and thus requires no fine tuning of parameters, meaning it can be applied to a wide range of dynamical systems. Furthermore, our extension does not suffer from the same stiffness issues arising in resistive MHD, and thus can be evolved quickly using explicit methods, with performance benefits of roughly an order of magnitude compared to current methods.

Figures

Equations

$${\partial }_t\left(\begin{array}{c}D\\ S^j\\ \tau \\ B^k\end{array}\right)+{\partial }_i\left(\begin{array}{c}D{v}^i\\ S^j{v}^i+p^{\ast }{\delta }^{ij}-b^j{B}^i/W\\ \tau v^i+p^{\ast }{v}^i-b^0{B}^i/W\\ v^i{B}^k-v^k{B}^i\end{array}\right)=0,$$

$$D$$

$$S^j$$

$$\tau$$

$$W$$

$$p^{\ast }$$

$$h^{\ast }$$

$$b^0$$

$$b^i$$

$$b^2$$

$${\partial }_t\left(\begin{array}{c}D\\ S_j\\ \tau \\ B^k\\ E^k\\ \varrho \end{array}\right)+{\partial }_i\left(\begin{array}{c}D{v}^i\\ S_j^i\\ S^i-D{v}^i\\ {ϵ}^{ijk}{E}_j\\ -{ϵ}^{ijk}{B}_j\\ J^i\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ -J^k\\ 0\end{array}\right).$$

$$\left(\begin{array}{c}D\\ S_i\\ \tau \\ J_i\\ S_{ij}\end{array}\right)=\left(\begin{array}{c}\rho W\\ \rho h{W}^2{v}_i+{ϵ}_{ijk}{E}^j{B}^k\\ \rho h{W}^2-p+\frac{1}{2}(E^2+B^2)-\rho W\\ \varrho v_i+W\sigma [E_i+{ϵ}_{ijk}{v}^j{B}^k-(v_k{E}^k)v_i]\\ \rho h{W}^2{v}_i{v}_j+[p+\frac{1}{2}(E^2+B^2)]{\delta }_{ij}-E_i{E}_j-B_i{B}_j\end{array}\right),$$

$${\partial }_t\mathbf{q}+{\partial }_i{\mathbf{f}}^{\mathbf{i}}(\mathbf{q},\overline{\mathbf{q}})$$

$${\partial }_t\overline{\mathbf{q}}+{\partial }_i{\overline{\mathbf{f}}}^i(\mathbf{q},\overline{\mathbf{q}})$$

$$\overline{\mathbf{s}}(\mathbf{q},{\overline{\mathbf{q}}}_0)=0.$$

$$\overline{\mathbf{q}}={\overline{\mathbf{q}}}_0+ϵ{\overline{\mathbf{q}}}_1+{ϵ}^2{\overline{\mathbf{q}}}_2+\mathcal{O}({ϵ}^3),$$

$$\mathbf{A}({\overline{\mathbf{q}}}_0+ϵ{\overline{\mathbf{q}}}_1+\mathcal{O}({ϵ}^2))=\mathbf{A}({\overline{\mathbf{q}}}_0)+ϵ{\partial }_{\overline{\mathbf{q}}}\mathbf{A}(\overline{\mathbf{q}}){|}_{\overline{\mathbf{q}}={\overline{\mathbf{q}}}_0}\cdot {\overline{\mathbf{q}}}_1+\mathcal{O}({ϵ}^2),$$

$${\partial }_t{\overline{\mathbf{q}}}_0+{\partial }_x{\overline{\mathbf{f}}}_0=\text{partialderivative}{\overline{\mathbf{s}}}_0\overline{\mathbf{q}}{\overline{\mathbf{q}}}_1.$$

$$\text{partialderivative}{\overline{\mathbf{q}}}_{0}t=\text{partialderivative}{\overline{\mathbf{q}}}_0\mathbf{q}\text{partialderivative}\mathbf{q}t,$$

$${\partial }_t\mathbf{q}+{\partial }_x{\mathbf{f}}_0={\mathbf{s}}_0,$$

$${\overline{\mathbf{q}}}_1$$

$$\text{partialderivative}\mathbf{q}t+\text{partialderivative}(\mathbf{f}+\hat{\mathbf{f}})x={\mathbf{s}}_0+\hat{\mathbf{s}}+\text{partialderivative}\hat{\mathbf{D}}x.$$

$$\hat{\mathbf{f}}$$

$$\hat{\mathbf{s}}$$

$$\hat{\mathbf{D}}$$

$$\text{partialderivative}\mathbf{q}t+\text{partialderivative}\mathbf{f}x=\text{partialderivative}\hat{\mathbf{D}}x.$$

$$\displaystyle\partialderivative{\bm{q}}{x}+\partialderivative{\bm{f}}{x}=-\partialderivative{x}\bigg{[}\frac{\partial\bm{f}}{\partial\overline{\bm{q}}}\bigg{(}\frac{\partial\overline{\bm{s}}_{0}}{\partial\overline{\bm{q}}}\bigg{)}^{-1}\bigg{(}\frac{\partial\overline{\bm{f}}_{0}}{\partial\bm{q}}-\frac{\partial\overline{\bm{q}}_{0}}{\partial\bm{q}}\frac{\partial\bm{f}}{\partial\bm{q}}\bigg{)}\partialderivative{\bm{q}}{x}\bigg{]},$$

$${\underline{\underline{M}}}_1$$

$${\underline{\underline{M}}}_2$$

$$\text{partialderivative}\mathbf{f}\mathbf{q}=\text{partialderivative}\mathbf{f}\mathbf{w}\text{partialderivative}\mathbf{w}\mathbf{q},$$

$${\underline{\underline{M}}}_1$$

$$\displaystyle=\partialderivative{\bm{f}}{\bm{w}}\bigg{(}\partialderivative{\overline{\bm{s}}_{0}}{\bm{w}}\bigg{)}^{+}$$

$${\underline{\underline{M}}}_2$$

$${\frac{\partial \mathbf{f}}{\partial \mathbf{w}}}^T$$

$$\displaystyle\bigg{(}\frac{\partial\overline{\bm{s}}}{\partial\bm{w}}\bigg{)}^{+}$$

$${\mathcal{A}}_{ij}$$

$${\mathcal{B}}_{ij}$$

Full text

††pubyear: 2019††pagerange: A resistive extension for ideal MHD–References

1 Motivation

Modern simulations of astrophysical systems typically involve a number of physical models. For systems that exhibit strong gravitational fields, the equations of general relativity require evolution to describe how the geometry of spacetime develops. As the evolution of the spacetime is coupled to the matter within it, motion of the matter will change the geometry of spacetime, and the geometry of spacetime alters the motion of the matter.

In totality, this system represents a highly non-linear, general relativistic, magneto-hydrodynamics (GRMHD) model. By-and-large, the equations of GRMHD conform to a type of MHD known as the ideal approximation (Font, 2009). For a wide range of systems, the ideal description of GRMHD performs well, and the simplification of perfect electrical conduction, among others, is generally accurate.

There does exist, however, a number of problems in astrophysics that require a description of MHD with non-zero resistivity (Palenzuela et al., 2009; Dionysopoulou et al., 2013; Andersson, Comer & Hawke, 2016). Simulations of binary neutron star mergers (Dionysopoulou, Alic & Rezzolla, 2015), magneto-rotational instabilities (Qian et al., 2016), magnetic reconnection (Lyubarsky, 2005; Mignone et al., 2012) and jet launching from black-hole accretion (Qian, Fendt & Vourellis, 2018) have all shown a sensitive dependence upon the magnitude of the electrical resistivity. And yet, the number of large-scale numerical astrophysics codes that implement resistive GRMHD remains small.

One issue with more physically complex models that prevents the wide spread adoption of the resistive GRMHD equations is the additional computation required in their evolution. Due to the possibly stiff source term of resistive GRMHD (Palenzuela et al., 2009), an implicit integrator, often the IMEX schemes of Pareschi & Russo (2004), is required to keep execution times reasonable (Palenzuela et al., 2009; Dionysopoulou et al., 2013; Aloy & Cordero-Carrión, 2016; Miranda-Aranguren, Aloy & Rembiasz, 2017; Wright & Hawke, 2019). Whilst these schemes allow for a faster evolution than would be possible using traditional explicit schemes for the majority of the parameter space, they still result in at least a factor $5\times$ slow-down compared to conventional, ideal GRMHD models (Dionysopoulou, Alic & Rezzolla, 2015). Furthermore, rewriting a well-tested code-base to evolve the resistive GRMHD equations is no small feat, requiring an increase in the number of fields to evolve, new flux and source vectors and new primitive recovery procedures, along with implementing and testing a new class of integrator.

Model extensions

Sub-grid source terms are beginning to see a number of applications in modelling astrophysical systems. The general principle behind these extensions is to include an additional source term into the equations of motion of the system, aimed at emulating the effects of unresolved fluid motion.

A common application of sub-grid sources is in the modelling of classical turbulence simulations, also known as large-eddy simulations (LES). In LES, the equations of motion are explicitly redefined in terms of resolved and unresolved quantities, and a closure is then assumed that relates the values of the sub-grid scale fields to those that are resolved in the simulation. Using this technique, it is possible to approximate the behaviour that would result from more finely resolved grids. Radice (2017) applied the classical Smagorinsky closure to the equations of general relativistic hydrodynamics for a merger simulation, showing that by modelling the sub-grid scale turbulence, the collapse of the hyper-massive neutron star remnant is altered. Similar work by Vigano & Palenzuela (2019) uses the LES closure approach to model unresolved motion of the non-relativistic MHD equations, with hopes of an extension to the full GRMHD system.

In other work, Giacomazzo et al. (2015) developed a sub-grid model, similar to Palenzuela et al. (2015), that captured how the behaviour of the magnetic fields is altered by unresolved motion. The motivation for such a source term comes from the assumption that the unresolved turbulent motion resulting from the Kelvin-Helmholtz instability (KHI) amplifies the magnetic fields. While this assumption may be justified (Price & Rosswog, 2006; Obergaulinger, Aloy & Müller, 2010; Kiuchi et al., 2015), the free parameters in this model are constrained via local simulations of KHIs (Zrake & Macfadyen, 2013), rather than a first principles argument.

The benefits of these model extensions lie in the huge reduction in computational cost required to evolve the system. These extensions capture physics that occurs in unresolved regions of the domain, regions which would require a prohibitively high resolution to capture with conventional methods—the smallest scales on which the KHI develops during mergers simulations, for example, is on the order $\sim 10$ cm (Palenzuela et al., 2015). Employing these methods allows one to capture the range of relevant physics in a fraction of the time, and with very little change to currently existing numerical codes.

In this paper, we develop an extension to the special relativistic, ideal MHD equations, that captures the resistive effects present in the full, resistive MHD model. This extension, dubbed a resistive extension generated for ideal magneto-hydrodynamics (REGIME), is derived from first principles arguments, and as such requires no fine tuning of parameters for different astrophysical scenarios.

The rest of the paper is laid out as follows. In Section 2 we briefly present the two models of MHD that we are concerned with: special relativistic, ideal and resistive MHD. Section 3 derives the proposed source, details how one should implement this new extension numerically, and discusses its effect on the stability of simulations. Results of the extension are presented in Section 4, and finally, in Section 5, we summarise the findings of the previous sections, discuss how they fit into current astrophysical simulations, and propose the future direction of the project.

2 Magneto-hydrodynamic Models

In this section, we outline two models of MHD that are used in relativistic astrophysics. In order to simplify the numerics in later sections and to test the validity of the method, we will limit ourselves to special relativity. In moving to a general relativistic description, only the form of the equations should change, and so the analysis we perform here should still apply—this includes implementations of GRMHD which evolve the vector potential instead of the magnetic fields (Etienne, Liu & Shapiro, 2010). We will also adopt the Einstein summation convention over repeated indices, where sums are over the three spatial dimensions, $\{x,y,z\}$, and $\delta_{ij}$ and $\epsilon_{ijk}$ are the Kronecker delta and Levi-Civita tensor, respectively. We use units where $\mu_{0}=\epsilon_{0}=1$, such that $c=1$, throughout.

2.1 Ideal MHD

The first model we present is that of ideal MHD, and this is the model that we will extend by means of the sub-grid source. Ideal MHD is the form most commonly used, for example, in merger simulations, due to the relative simplicity of its mathematical form, and the numerical methods employed in solving it. For a more comprehensive description of ideal MHD see Noble & Gammie (2006); Font (2009); Antón et al. (2010); Siegel et al. (2018). We use the notation $\partial_{t}$ and $\partial_{i}$ to represent the partial derivative with respect to time and the $i^{\text{th}}$-spatial coordinate.

The equations of motion, in conservative form, are given as

[TABLE]

where the conserved quantities, $\{D,S^{j},\tau,B^{k}\}$, correspond to the fluid density, specific momentum in the $j^{\text{th}}$-direction, kinetic energy density and magnetic fields, and are related to the primitive quantities, $\{\rho,v^{j},p,B^{k}\}$, namely the rest-mass density, fluid velocity in the lab frame, hydrodynamic pressure and the magnetic field respectively, via

[TABLE]

In addition, we have the relations for the Lorentz factor, $W$, total pressure and specific enthalpy (i.e. including contributions from the magnetic fields), $p^{*}$ and $h^{*}$, and the relations of the components of the magnetic four-vector, $b^{\mu}$, to the lab-frame magnetic field, $B^{i}$:

[TABLE]

The primitive quantities along with an equation of state completely define the system. Throughout this paper, we will assume a $\Gamma$-law equation of state, where $p=\rho e(\Gamma-1)$, $e$ is the specific internal energy and $\Gamma$ is the ratio of specific heats, with the specific enthalpy defined as $h=1+e\Gamma$. Whilst the choice of the equation of state will change the form of the source term presented at the end of Section 3.2, the general form of the source term will be unaltered.

Although the system is defined by the primitive quantities, it is the conserved quantities that are evolved explicitly. This means, therefore, that there must be some scheme for transforming from the conserved variables to the primitives such that the fluxes in equation (1) can be computed. The transformations can be found in detail in the aforementioned references.

2.2 Resistive MHD

The second model we present is resistive MHD, and is the model whose behaviour we wish to replicate by means of an extension of ideal MHD. Once again, further details of resistive MHD may be found in Palenzuela et al. (2009); Dionysopoulou et al. (2013). The form of resistive MHD in balance law form reads

[TABLE]

The state vector has picked up four additional fields, three from the electric field and one corresponding to the charge density, $\{E^{i},\varrho\}$ respectively. It is not essential that we evolve the charge density as one can compute it from Gauss’ law, but we choose to do this to avoid difficulties when taking derivatives of near-discontinuous data.

Once again, we relate the conserved vector, the charge density current, $J_{i}$, and the momentum flux, $S_{ij}$, to the primitive quantities with,

[TABLE]

where all quantities are interpreted in the same way as for ideal MHD. In addition to the equation of state, the form of the charge current density is also needed to relate the electric and magnetic fields to the charge density and close the system. We also note the inclusion of the conductivity, $\sigma$, in the definition of the charge current density. This is to be expected, as for resistive MHD there is no assumption that charge flows perfectly (as there is in ideal MHD) and so the conductivity can take any finite, non-negative value. We shall see how this causes difficulties in the next section.

We should also note that the conductivity we present here is scalar, but that in general can be described by a tensor (Dionysopoulou, Alic & Rezzolla, 2015). While including the full tensor description would introduce additional terms into the source of the electric fields, it would not change the analysis that will follow. In this paper, we limit our analysis to a scalar conductivity commonly used in astrophysical settings (Komissarov, 2008; Dionysopoulou, Alic & Rezzolla, 2015; Mohseni et al., 2015; Miranda-Aranguren, Aloy & Rembiasz, 2017).

As is the case for ideal MHD, we also require some algorithm to compute the primitive from the conserved quantities, two examples of which are found within Palenzuela et al. (2009); Dionysopoulou et al. (2013).

Numerical difficulties

The major difference between the two models of MHD that are presented here lie in the form of the source terms. We can see that the source of the electric fields in the resistive model, equation (4), is non-zero and proportional to the conductivity. As mentioned, charge flows without resistance in ideal MHD, and thus the fluid has an infinite conductivity. Clearly then, near the ideal limit when $\sigma\to\infty$, the source term in the resistive MHD equations begins to dominate. This is a common feature of balance laws equations, and can be interpreted as some driving term that acts on a much shorter timescale than the conservation system.

When the timescale that the source acts on is shorter than the timescale of the simulation, $\epsilon\lesssim\Delta t$, the system is said to be stiff, where in this case $\epsilon\propto 1/\sigma$. This results in some difficulties regarding the numerical evolution of the equations. In order to maintain a stable evolution, it is sufficient to either dramatically reduce the size of the timestep used in the simulation, or to employ a set of implicit or semi-implicit time integrators (Pareschi & Russo, 2004). Unfortunately, both of these approaches will increase the execution time of a simulation by around a factor of five or more, depending upon the magnitude of the conductivity.

What would be useful is some source term that captures the behaviour due to finite conductivity, but that avoids the numerical difficulties in the ideal limit of the full, resistive MHD system. The following section will derive such a source term using a Chapman-Enskog-type analysis, but first we introduce the notation that will be used.

We can re-write system (4) in the following, more compact way:

[TABLE]

where we indicate equations which become stiff as $\sigma\to\infty$ with an over-bar. This means that $\overline{\bm{q}}=\{E_{x},E_{y},E_{z}\}$ with the corresponding fluxes, $\overline{\bm{f}^{i}}(\bm{q},\overline{\bm{q}})$, and sources, $\overline{\bm{s}}(\bm{q},\overline{\bm{q}})$, taken from equation (4). The remaining variables are non-stiff in the ideal limit, and denoted $\bm{q}=\{D,S_{x},S_{y},S_{z},\tau,B_{x},B_{y},B_{z}\}$, where we have left out the charge density. We will also denote the vector of primitive variables with $\bm{w}=\{\rho,v_{x},v_{y},v_{z},p,B_{x},B_{y},B_{z},E_{x},E_{y},E_{z},\varrho\}$, where $\varrho$ is the charge density.

As an aside, one might notice that the flux in the charge density evolution may be of the same size as the source of the electric fields, and so the equation can also be stiff. Instead, if we were to multiply this equation by the small timescale, $\epsilon$, we can make the equation non-stiff. Regardless, as the charge density is, in fact, defined by Gauss’s law, the choice of evolving it is purely a numerical one and so it can be excluded from the following expansion.

3 Chapman-Enskog Expansion

In this section we will use the Chapman-Enskog (CE) method of expansion to derive the form of the REGIME source term. The CE expansion was initially used to determine higher order moments for solutions to the Boltzmann equation (Denicol & Noronha, 2018). By expressing the particle distribution function in terms of powers of some small quantity, one can approximate the true solution up to arbitrary order, with each subsequent power encapsulating some new physical phenomenon.

An additional application of this expansion is presented in LeVeque (2002), in which LeVeque demonstrates how a coupled system of potentially stiff balance law equations may be reduced, in a given limit, to a single, modified system. Here, we apply the same logic to ideal and resistive MHD. We consider the solution to ideal MHD to be the equilibrium system of resistive MHD, and derive the form of a perturbation on top of this. As we will see, this perturbation looks like a diffusion term, as is the case in LeVeque’s analysis.

3.1 System Derivation

We begin with equations (6), recalling that the non-stiff and stiff conserved variables are labelled $\bm{q}$ and $\overline{\bm{q}}$ respectively. In order to maintain finite solutions in the ideal limit, we require that $\lim_{\epsilon\to 0}\overline{\bm{s}}(\bm{q},\overline{\bm{q}})=\bm{0}$. We will refer to the solution in the limit $\epsilon\to 0$ as the equilibrium system, as for the models we are interested in here this corresponds to the RHS $=0$, and equation (6a) reduces to ideal MHD. In the case of resistive MHD, this limit corresponds to ideal MHD where $\sigma=\infty$. Therefore, for large but finite conductivities, we can think of the solution of resistive MHD as some perturbation on top of the solution of ideal MHD. This corresponds to small $\epsilon$ when compared to other timescales in the simulation, for example $\Delta t$.

To derive the form of this perturbation, we need to define the equilibrium solution. At this point, we assume that the stiff variables have some values given by $\overline{\bm{q}}_{0}$, and so the stiff source in the equilibrium limit must be

[TABLE]

The Chapman-Enskog expansion involves expressing the stiff variables in terms of a series of increasing powers of the small quantity, $\epsilon$. Expanding the solution for the stiff variables about their equilibrium solution gives

[TABLE]

in which each new power of $\epsilon$ represents some higher order perturbation on top of the previous. Using the Taylor expansion of some vector function as

[TABLE]

we can expand the stiff and non-stiff variables, equation (6), about the equilibrium solution. Note, we are considering a 1D system for simplicity—in 3D the following analysis is qualitatively the same, and will be presented in Section 3.4.

Keeping only leading order terms, we get an evolution equation for the equilibrium state of the stiff system up to order $\mathcal{O}(\epsilon^{0})$:

[TABLE]

We will use the notation $\bm{A}_{0}$ to represent the vector $\bm{A}(\bm{q},\overline{\bm{q}}_{0})$, evaluated at the equilibrium point—therefore, $\partial{\bm{s}(\bm{q},\overline{\bm{q}})}/{\partial\overline{\bm{q}}}|_{\overline{\bm{q}}=\overline{\bm{q}}_{0}}\equiv\partial\overline{\bm{s}}_{0}/\partial\overline{\bm{q}}$.

We now wish to determine the form for the first order perturbation of the stiff variables, $\overline{\bm{q}}_{1}$. First, we assume that the equilibrium state of the stiff variables can be fully characterised by the non-stiff variables—that is, $\overline{\bm{q}}_{0}=\overline{\bm{q}}_{0}(\bm{q})$. For the case of resistive and ideal MHD, we know this is true as the electric fields can be expressed as $\bm{E}=-\bm{v}\times\bm{B}$. Therefore, as

[TABLE]

and the evolution of the non-stiff variables, up to zeroth order in $\epsilon$, is

[TABLE]

we can remove the time dependence of the first order perturbation, equation (9), and re-arrange for $\overline{\bm{q}}_{1}$:

[TABLE]

Here, we note that the non-stiff flux evaluated at the equilibrium point is only a function of the non-stiff variables, $\bm{f}(\bm{q},\overline{\bm{q}}_{0}(\bm{q}))\equiv\bm{f}(\bm{q})\equiv\bm{f}_{0}\equiv\bm{f}$.

Upon substituting the form for the perturbation, equation (11), into the expanded evolution equation for the non-stiff variables, and keeping orders up to $\mathcal{O}(\epsilon)$, we can rewrite the perturbed system in an intuitive and compact form:

[TABLE]

where we have absorbed the timescale, $\epsilon$, into the definition of the stiff source, $\bar{\bm{s}}$.

In this form, the modified flux term, modified source term, and the diffusion-like term are

[TABLE]

The hats in equations (13) signify that these terms are perturbations on top of the equilibrium system, coming in at order $\mathcal{O}(\epsilon^{1})$.

Looking at the form for the ideal MHD equations (1), we can see that there is zero source term—$\bm{s}(\bm{q})\equiv 0$. As a result, the modified flux and source terms in the previous equations reduce to zero, and we are left with only the diffusion term being non-zero. This means the equations of motion of ideal MHD are modified only by a single diffusion term111We note that general relativistic, ideal MHD has a non-zero source term, and so this simplification may not be made.,

[TABLE]

We can see that in the ideal limit, which corresponds to $\{\sigma,\epsilon\}\to\{\infty,0\}$, the fact that $\hat{\bm{D}}\propto\epsilon$ means equation (14) reduces to the standard form for ideal MHD, as expected. Then, as the conductivity reduces, larger corrections are made through the diffusive source term to mimic behaviour that should be present in resistive MHD.

Finally, observe how the source term for resistive MHD is proportional to $\sigma$, but that the REGIME source term scales as $\epsilon\propto\sigma^{-1}$. This means that the two models become stiff in opposing limits—near the ideal regime (large $\sigma$) REGIME will be stable as a result of a small source term, and will only become stiff, and potentially unstable, as $\sigma\to 0$. The big benefit of this behaviour is that near the ideal regime we can confidently evolve REGIME with explicit time integrators, knowing that source contributions will remain small. In contrast, in the event of very low conductivities, $\sigma\sim 0$, it will not be sensible to evolve REGIME using implicit schemes, not least because in this regime resistive MHD is likely to be stable with explicit integrators. Because the numerical flux function appears in the diffusion vector (i.e. $\partial_{x}\bm{f}$), and this has a dependence on neighbouring cells, an implicit integration for REGIME would require solving the system for all cells in the domain at once (i.e. an $N_{x}\times N_{y}\times N_{z}\times N_{cons}$ dimensional root-find, where $N_{cons}$ is the size of the conserved vector).

3.2 Numerics

Now we have the form for the source term, our attention turns to implementing it numerically. The interpretation of a term in a fluid’s equations of motion may be intuited by the order of the spatial derivative of the conserved fields. For example, advection is noted by a first order derivative, diffusion by a second order derivative and dispersion via a third order derivative. This is why we call the new piece in equation (14) a diffusion-like term, as it can be re-written to include second order spatial derivatives of the conserved fields:

[TABLE]

recalling that both $\hat{\bm{f}}=\hat{\bm{s}}=\bm{0}$ when considering ideal and resistive MHD in the special relativistic limit.

In order to implement this new piece, it is useful to understand more about what it involves. If we look that the diffusion vector given in equation (13c), the first term is the Jacobian of the non-stiff system with respect to the stiff variables, $\partial_{\overline{\bm{q}}}\bm{f}$, and is therefore a matrix operation—as, indeed, is the second and fourth term. The third and fifth terms correspond to the spatial derivatives of the stiff and non-stiff flux vectors. In fact, the form of the diffusion term in equation (13c) is preferable to that in equation (15), as any numerical flux function that computes $\mathcal{\bm{F}}(\bm{q})=\partial_{x}\bm{f}(\bm{q})$ can be re-used as is for the generation of the source term.

All that is required then is to determine the forms of the two matrices,

[TABLE]

Differentiating with respect to the conserved variables is challenging, as we do not know how to express the flux and source vectors in terms of only the conserved quantities. We can, however, express the conserved, flux and source vectors in terms of the primitive variables, and thus we can differentiate using the primitive quantities. For example,

[TABLE]

with $\bm{w}$ as the vector of primitive variables. In this way, the new form for the $M$ matrices is

[TABLE]

Here, we have used the superscript + to denote the Moore-Penrose pseudoinverse (Barata & Hussein, 2011). As the length of the vectors $\overline{\bm{s}}_{0}$ and $\bm{w}$ are, in general, not the same, the resulting matrix will not be square, nor have a corresponding inverse. As a result, we use the definition of the right-pseudoinverse of a matrix $Q$ as, $Q^{+}=Q^{T}(QQ^{T})^{-1}$. We will use the term inverse to refer to the pseudoinverse henceforth.

Here, we have a choice of how to compute the matrices of interest—that is we can invert them numerically, or try to get the form of the inverted matrix symbolically. Inverting matrices numerically, especially when densely populated, can require a large amount of computation, reducing accuracy as well as slowing down simulations. If the algebraic form of the matrices were at hand, this would lead to a far more efficient simulation, and as we are trying to build a source term to extend ideal MHD with the intention of being faster to evolve than resistive MHD, it is sensible to adopt the performance gains of a purely symbolic source term.

In order to generate human readable terms for equations (18), we have had to make a number of assumptions. Firstly, we will assume a low velocity limit, in which terms of $\mathcal{O}(v^{2})$ can be ignored. For instance, with this approximation in resistive MHD, the $y$-momentum flux in the $x$-direction reduces to $S_{12}=-E_{x}E_{y}-B_{x}B_{y}$. Secondly, to reduce the number of terms in the inverted matrices, we assume that the fluid only couples weakly to the magnetic field, enforcing this by setting the electric and magnetic fields to zero. As the electric and magnetic fields are still present in the un-inverted matrices, their influence still makes it through to the final source term. As a check of these assumptions and their effect on the system we also implemented the REGIME source term numerically with no approximations. After doing this, we found virtually no difference in the simulation output for mildly relativistic flows, $v\sim 0.5c$. Making these simplifications, we can use a symbolic algebra program to compute the form for the source term. For this task, we used Wolfram Mathematica v11.2 (Wolfram Research, 2018). The full notebooks, along with greater detail, are available through the METHOD GitHub page222https://github.com/AlexJamesWright/METHOD, so we will only show the final forms here.

For the systems of equations presented here, namely ideal and resistive MHD in the special relativistic limit and with the above assumptions made, we get the simple result that $\underline{\underline{M}}_{2}=0$. The reason for this is that $\partial\overline{\bm{q}}_{0}/\partial\bm{w}$ is a $3\times 12$ matrix in which only three entries are non-zero—these correspond to the $\partial E_{x}/\partial E_{x}$, $\partial E_{y}/\partial E_{y}$ and $\partial E_{z}/\partial E_{z}$ terms. These terms are then dotted with $\partial\bm{E}/\partial\bm{S}$, which on account of the weak coupling approximation is zero.

With regards to the $\underline{\underline{M}}_{1}$ matrix, we get

[TABLE]

[TABLE]

where we define the elements of the matrices $\mathcal{A},\mathcal{B}$ and $\mathcal{C}$ as

[TABLE]

the vector $\mathcal{D}$ as

[TABLE]

and the pre-factor, $\alpha$, as

[TABLE]

The pre-factor, $\alpha$, in equation (20) now acts in a similar way to the previous timescale, $\epsilon$. That is, in the ideal limit, $\alpha\propto\sigma^{-4}\to 0$, and the source term tends to zero, recovering ideal MHD. For large but finite conductivities however, the source term will modify the solution.

It should also be noted that in equation (20), the charge density, $\varrho$, appears a number of times. As this term is not explicitly evolved in ideal MHD, we compute its value from Gauss’ law, $\bm{\nabla}\cdot\bm{E}=\varrho$. Furthermore, as the electric fields are not evolved in the ideal prescription, we use the assumption of our previous expansion, equation (7). This condition, that $\overline{\bm{s}}(\bm{q},\overline{\bm{q}}_{0})=0$, implies that there is zero current, $\bm{J}=0$. For this to be valid as $\sigma\to\infty$, it would require that the leading order part, proportional to $\sigma$, is zero, and thus that $\bm{E}=-\bm{v}\times\bm{B}$, which is simply the ideal form for the electric fields.

All that remains to compute the diffusion vector is the $\partial_{x}\overline{\bm{f}}_{0}$ term. This can be done using any existing numerical flux function, which in the code used for this paper is a WENO3 reconstruction (Shu, 1997).

When determining the derivative of the diffusion vector, there are a number of numerical methods available, each of which can have an affect on the accuracy of the solution and the overall stability. For example, employing MINMOD slope limiting to reduce the onset of Gibbs oscillations also improves the accuracy of the solution in most scenarios. Despite this, in Kelvin-Helmholtz simulations, Section 4.2, limiting in this way leads to excessive growth of the magnetic fields, and attempts to identify alternative slope limiting methods have all resulted in similar behaviour. As a result, we simply employ second-order central differencing, which we find gives good results across all simulations and is generally stable at shocks.

3.3 Stability Criterion

The source term we have derived is diffusion-like, containing a second order spatial derivative. Whilst the original equilibrium system is in a strongly hyperbolic form, the addition we have made would be parabolic on its own, i.e. if there were no flux. As a result, we must first understand the limits on the spatio-temporal resolution such that the new system is stable.

Recall our system has the following form:

[TABLE]

In order to proceed, we will make the assumption that in a small region of the domain, $\underline{\underline{M}}_{1}$ is approximately constant. We will also rewrite the stiff flux derivative using the Jacobian, and further assume that it, too, is constant in this region. Finally, if we assume that the resolution constraints set by the hyperbolic system are met, and that we are only interested in the second derivative term, then we can simplify and say

[TABLE]

Here, we make an analogy with the scalar case of the diffusion equation, $\partial_{t}q=k\partial_{x}^{2}q$, in which an explicit scheme is generally considered stable if

[TABLE]

in which $\Delta t$ and $\Delta x$ are the size of the timestep and spatial resolution in a simulation.

For the system, we suggest that $k$ be given by the largest magnitude eigenvalue of $\underline{\underline{K}}$. This, therefore, requires an eigenanalysis of the matrix

[TABLE]

With the help of Mathematica (these scripts can also be found on the GitHub page), it can be shown that the largest eigenvalue is given by

[TABLE]

This result suggests that, assuming a CFL constraint given by $\Delta t/\Delta x=\nu$, the spatial resolution of any simulation using this source term must (roughly) satisfy

[TABLE]

where $\delta$ is some softening factor which we determine in section 4.3. Intuitively, this means there is a limit on how well a simulation may be resolved. For high conductivities, this limit is not an issue as the RHS tends to zero, allowing (almost) arbitrarily fine resolutions. For resistive simulations, however, there is a finest grid resolution for which the system is numerically stable (for a given Courant factor), as the source term begins to become large compared to the hyperbolic part. Of course, one may always reduce the CFL condition and take smaller timesteps in order to achieve higher resolutions.

3.4 Higher dimensions

This section has so far assumed a one-dimensional system to simplify the analysis of the proposed model. In this subsection, we relax that assumption. As the majority of the analysis carries over to higher dimensions in the same manner, we only present the bottom line results.

In general, we will have the non-stiff and stiff systems defined via

[TABLE]

The index $i$ spans the three spatial dimensions, such that $\bm{f}^{i}$ is the non-stiff flux vector in the $i^{th}$ coordinate direction, and $\partial/\partial x^{i}$ is the derivative with respect to the $i^{th}$ coordinate direction.

Applying the same expansion about the equilibrium system, the perturbed system can be written, in analogy with equation (12) such that the small scale is absorbed into the definition of $\hat{\bm{f}},\hat{\bm{s}}$ and $\hat{\bm{D}}$, as

[TABLE]

We note here that each spatial dimension has a corresponding modified flux and diffusion vector. The definitions of all modifications are given by the following expressions:

[TABLE]

where we recall that repeated indices indicate a summation. Clearly, as we move to higher dimensions, the amount of computation required increases correspondingly—3D simulations require the calculation of three diffusion terms, compared to only one for a 1D simulation. We discuss the performance of REGIME in section 4.4.

With regards to the stability analysis, equation (27) becomes,

[TABLE]

where there are now nine matrices,

[TABLE]

There are now cross terms in the derivative, but to simplify the analysis we once again take the largest eigenvalue of all $K^{ij}$ to be the value of $k$ as in equation (29). The largest eigenvalue is the same as in the previous section, and so we assume the same form for the resolution criterion, equation (31). As this is an order of magnitude estimation of the stability, we assume that any effects due to the additional terms are captured in $\delta$ to be determined in Section 4.3.

4 Results

In this section we present a series of results relating to the convergence of REGIME with the conductivity, the stability of the proposed model and its performance relative to the resistive model. We then end this section with a discussion on the current limitations of this model.

For the convergence plots presented here, we will consider the results of the resistive MHD model to be the exact solution. Furthermore, any performance results have been optimised to be the fastest execution for a given set of parameters—for example, the runtimes represent the largest timestep (optimum CFL factor) possible for a specific resolution, conductivity and integrator, whilst generating stable simulations and indistinguishable results for a given problem.

4.1 METHOD

All results have been generated using the METHOD code. Details of the numerical schemes used in METHOD can be found in more detail in (Wright & Hawke, 2019). Briefly, the numerical flux is computed using flux splitting, with a third order WENO reconstruction scheme (Shu, 1997). Time integration uses one of two schemes: an operator-split, second-order Runge-Kutta (RK) method (Gottlieb & Shu, 1998); and the second-order implicit-explicit (IMEX) RK scheme from Pareschi & Russo (2004). In order to maintain the divergence constraints imposed by Maxwell’s equations we use the method of hyperbolic divergence cleaning (Dedner et al., 2002).

4.2 Convergence with $\sigma$

We now present results to demonstrate how the new model behaves with varying magnitudes of the conductivity. A number of well-known initial data are evolved, and the results of resistive MHD and REGIME are compared and contrasted. For clarity, we have also evolved some simulations using ideal MHD to highlight its inability in modelling various resistive behaviours, and to show how the REGIME extension to it can capture these behaviours.

Brio-Wu shock tube

The Brio-Wu shock tube test (Brio & Wu, 1988) is a standard numerical fluid problem used to assess how a model behaves when there is discontinuous data. The 1D domain is separated into two regions by a partition which is removed at time $T=0$. The initial data for this problem is given by $(\rho,p,B_{y})=(1,1,0.5)$ for $x<L/2$, and $(0.125,0.1,-0.5)$ for $x\geq L/2$, all other variables are set to zero. To make comparisons with (Palenzuela et al., 2009; Dionysopoulou et al., 2013) we have set $\Gamma=2$. The system is run until $T=0.4$, using $N_{x}=128$ grid points.

The first three columns of figure 1 show the final state of the $y$-direction magnetic field and pressure for ideal MHD (blue dash-dotted), resistive MHD (black dashed) and REGIME (red solid), for varying conductivities. Clearly for the larger conductivities, there is little difference in the output of the models, which is to be expected as both resistive MHD and REGIME limit to ideal MHD. For more resistive simulations (smaller $\sigma$), the resistive models differ more greatly from the ideal MHD results, and REGIME and resistive MHD results are in excellent agreement.

In the final column of figure 1, we calculate the error as the $L_{2}$ norm of the difference between REGIME’s output and resistive MHD, and see how this norm varies with changes in $\sigma$. Clearly, there is a significantly reduced error growth when using REGIME over the ideal MHD model—the diffusion term captures the important features coming from the resistive model.

Self-similar current sheet

As the form of the proposed source term looks much like a diffusion term, we now present a problem in which the solution obeys the diffusion equation. The self-similar current sheet test (Komissarov, 2008) begins with a static fluid in equilibrium, and imposes upon it a magnetic field of the form $(B_{x},B_{y},B_{z})=(0,B(x,t),0)$. This is one of the few MHD tests in which an analytic solution exists, namely,

[TABLE]

where erf is the error function and the magnetic field satisfies

[TABLE]

The system is set up with $(p,\bm{v},\rho)=(50.0,0.0,1.0)$ at $T=1$, using $N_{x}=128$ cells and an adiabatic index of $\Gamma=2$. All simulations are run until the end time, $T=8.0$, with a magnetic field strength of $B_{0}=1$.

An ideal simulation of the current sheet problem should not evolve from the initial condition, which given an infinite resolution would be a step function. Simulations with greater resistivities should expect more diffusion, and so the solution should appear more smeared out. Clearly, for the two final states (left and middle plots) from figure 2, both the resistive MHD and REGIME capture this process, with excellent agreement between both down to conductivities of only $\sigma=20$. The error plot, rightmost in figure 2, shows little error growth between resistive MHD and REGIME.

Resistive magnetic reconnection

Next, to demonstrate the effectiveness of REGIME in higher dimensions, we present the 2D magnetic reconnection problem. Similarly to the current sheet problem, magnetic reconnection is a purely resistive phenomenon—the lack of electrical resistivity in ideal MHD should mean that there is no evolution from the initial condition (excluding effects due to numerical diffusion). Furthermore, this kind of phenomenon should occur extensively during simulations of neutron star mergers, due to the complex interaction between the progenitors’ intense magnetic fields, and so makes a useful test of REGIME in its applicability to these simulations.

Reconnection occurs when two magnetic fields whose vectors are not parallel intersect, forming a new magnetic vector and releasing stored magnetic energy into thermal and kinetic energy. There are multiple configurations of magnetic fields that will lead to some form of reconnection—indeed, understanding magnetic reconnection is so vital for solar physics that there is extensive literature and multiple models to choose from. As we are only interested in demonstrating the effect, we will adopt the first and simplest attempt of understanding reconnection, the Sweet-Parker model (Parker, 1957; Sweet, 1958). One should note, however, that the reconnection rate predicted by the Sweet-Parker model is slow compared to other models, and experimental data. In fact, Loureiro & Uzdensky (2016) claim that “unsatisfactory predictions are obtained for reconnection events in almost all plasmas that one cares to examine."

In the Sweet-Parker model, the reconnection layer is defined to have a thickness of $2\lambda$, and a length of $2L$. The initial magnetic field strength is given by $B_{0}$, with the corresponding Alfven velocity as $V_{A}$. With this, it is possible to estimate the rate of magnetic reconnection, $\mathcal{E}\equiv u_{\text{in}}/u_{\text{out}}$. By conservation-of-mass arguments, we have the relation

[TABLE]

where $u_{\text{in/out}}$ is the velocity of the fluid flowing in/out of the boundary layer. Using Ohm’s law, one can show that the rate of flow in to the boundary layer can be expressed as $u_{\text{in}}=\eta/\lambda$, in which the resistivity is $\eta=\sigma^{-1}$. The Lundquist number, $S$, is defined via $S\equiv LV_{A}/\eta$ and so $\eta/V_{A}=L/S$. Using the fact that $u_{\text{out}}=V_{A}$, one can show that

[TABLE]

The domain, where $x\in[-12.8,12.8]$ and $y\in[-6.4,6.4]$, is set up as follows:

[TABLE]

with all other variables set to zero. The magnetic fields are then perturbed by,

[TABLE]

Here, we have $(L_{x},L_{y})=(25.6,12.8)$, with the initial thickness of the boundary layer $\lambda=0.5$. The background density is $\rho_{\infty}=0.2$ with $\rho_{0}=1.0$, a magnetic field strength of $B_{0}=1.0$, and a perturbation size of $\psi_{0}=0.1$. The simulation is then run using $1024\times 512$ cells until $T=100$, with the adiabatic index as $\Gamma=2$. We use periodic boundaries along the $x$-axis and outflow boundaries along the $y$-axis.

In figure 3, we show the final state of the hydrodynamic pressure for a range of conductivities for both resistive MHD and REGIME. The hotspots in the pressure coincide with the formation of two magnetic islands. Clearly, for the simulations with greater conductivity, there is significantly less diffusion occurring—this is seen as a more distinct separation of the two magnetic islands. As the resistivity of the simulations is increased, this separation reduces more quickly. With these plots, one can see no differences in the output of resistive MHD and REGIME. The proposed source term of REGIME captures the diffusive behaviour of the full model exceptionally well.

Next, to determine the rate of diffusion for a given resistivity for both models, we follow Mignone et al. (2012) and compute the current density from Ampere’s law, $\bm{J}=\bm{\nabla}\times\bm{B}$, along the $y$-axis (at $x=0$). As the rate of reconnection can be written $\mathcal{E}=\lambda/L$, where $L=L_{x}$ is the length of the domain in the $x$-direction, we compute the width of the boundary layer, $\lambda$, as the width of the Gaussian profile that best approximates $\bm{J}$ at $T=100.0$.

As the fit for $\sigma=100$ is poor, it has been excluded from the reconnection rate analysis. The results of this analysis are shown in figure 4 for a number of resistivities, and shows a scaling slightly greater than what is expected from the Sweet-Parker analysis for both resistive MHD and REGIME. Clearly the subgrid source term captures this effect, and as a result the reconnection rate of the two models agree remarkably well.

Kelvin-Helmholtz instability

A common type of fluid instability that is believed to play an important role in evolution of magnetic fields during mergers (Kiuchi et al., 2015) and other astrophysical events is the Kelvin-Helmholtz instability (KHI). In this scenario, two differentially flowing fluids are perturbed along the boundary layer between them, resulting in the growth of vortices and a corresponding cascade of kinetic energy from large scales down to the smallest dynamos. The motion of the fluid couples to the magnetic fields, where the folding and twisting increases the fields’ energy density and strength. This is one mechanism by which magnetars are thought to develop such intense magnetic fields.

To investigate this process, we use the initial conditions from (Beckwith & Stone, 2011). We simulate a 2D domain where $x\in[-0.5,0.5]$ and $y\in[-1.0,1.0]$, with the initial $x$-velocity as

[TABLE]

the density as

[TABLE]

and then perturb the $y$-velocity as

[TABLE]

where the shear velocity is $v_{shear}=0.5$, a boundary layer thickness of $a=0.01$, the densities are given by $(\rho_{0},\rho_{1})=(0.55,0.45)$, and a perturbation amplitude of $A_{0}=0.1$ over a length $l=0.1$. The initial pressure is uniform, $p=1.0$, and we impose a perpendicular (to the flow) magnetic field of $B_{z}=0.1$. The adiabatic index is set to $\Gamma=4/3$ and the simulation is evolved until after the end of the linear growth stage, at $T=6.0$, using $(N_{x},N_{y})=(512,1024)$. As in resistive reconnection, we use periodic boundaries in the $x$-direction and outflow boundaries in the $y$-direction.

Whilst there is little-to-no difference between the observed state in the density between ideal MHD, resistive MHD or REGIME, figure 5, we can be more quantitative. As the evolution of the magnetic fields can be affected by the dynamics of this instability, in figure 6 we show their large scale behaviour. In the first plot, we show the average magnetic energy density across the whole domain, and the maximum magnetic field strength in the second.

In ideal MHD, the magnetic field lines lock to the fluid motion, a phenomenon known as magnetic freezing, and as such are efficiently folded within features such as vortices. In resistive MHD however, the magnetic fields are free to move, leading to less efficient folding and a reduced amplification. This effect is clearly captured by REGIME on account of the reduction in the average and maximum magnetic field strength.

We can also determine how REGIME behaves across all scales in the simulation. Kolmogorov predicted (Qian, 1994) what the distribution of energy should look like across all scales for fully developed turbulence, known as his $5/3$-law. Through dimensional arguments alone, he demonstrated that the power spectrum of the kinetic energy density should satisfy

[TABLE]

where $k$ is the wave number of the mode.

Using the procedure laid out in Beckwith & Stone (2011)333The scripts used for this analysis can also be found at the METHOD GitHub page., we compute the kinetic energy density power spectrum, and also the spectrum for the magnetic energy density, comparing the results for ideal and resistive MHD and REGIME. The results are given in figure 7. We can see that for the kinetic energy density spectrum, all models exhibit the predicted behaviour.

For the magnetic energy density power spectrum, however, we can see differences in how the models behave on small scales (large $k$). There is less energy transferred to the smaller scales for resistive MHD than for ideal MHD. This is due to no magnetic freezing, as discussed before.

REGIME seems to be able to capture this change in the dynamics down to a wavenumber of $k\approx 40$, but for larger modes (smaller scales) transfers too much energy to the magnetic fields. The fraction of energy that this represents in this simulation is minute, and as a result the dynamics of the simulation remain relatively unchanged.

4.3 Stability

We now turn our attention to the stability of the new model. In Section 3.3, we argue that there should be some relationship between the resolution and the system’s stability, and that the result should be some maximum possible spatial resolution whilst still giving stable evolutions.

Returning to the Brio-Wu shock tube test, section 4.2, we can see how quickly the simulation can move from a stable region to an unstable region. Figure 8 shows the final state of the density for the Brio-Wu shock tube, where $\sigma=1000$, for four different resolutions. The first two (top row) use $N_{x}=1000$ and $1498$ cells, leading to stable simulations. Adding only two additional cells, however, results in an unstable simulation, as can be seen by the spikes present in the data at $x\approx 0.2$. If we continue to use finer resolutions, the effect becomes more pronounced, and additional artefacts appear.

When performing the stability analysis, we assumed that in moving to higher dimensions it was the same largest eigenvalue that is the main factor in determining stability. To test this, we perform a number of simulations in all dimensions. The initial data for the Orszag-Tang vortex is taken from Beckwith & Stone (2011) for 2D, with the addition of a $z$-perturbation to the fluid velocity of the form $0.2\sin(2\pi z)$ for the three dimensional case.

Figure 9 shows the relationship between the finest possible resolution for a range of conductivities. The one-dimensional Brio-Wu and current sheet problems produce an average $\delta\approx 0.92$ and $0.93$, respectively. When we move to higher dimensions, the $\delta$ value for all problems drops to a value of $\delta\approx 0.3-0.4$, corresponding to a less restrictive stability criterion.

Whilst the reasoning behind this drop in $\delta$ is not fully understood, we believe it comes from the numerics. The fifth plot (bottom left) shows the data for a 2D Brio-Wu test, in which the standard 1D data is rotated about the $z$-axis by $\pi/4$. If the stability was set by the 1D nature of the initial data, a rotation about the $z$-axis should not change the resolution that is possible along the propagating shock wave, and we would expect the same softening factor. In reality, the rotated problem allows for higher resolutions to be achieved, possibly due to the effects of the cross derivatives in the REGIME source term.

Whatever the reason, the effect is to allow a greater resolution in the multi-dimensional case, tending towards a factor $3\times$ for large conductivities. If we therefore take a conservative value for the softening factor of $\delta=1$, we can determine whether the types of astrophysical simulation we are interested in are possible with REGIME. The conductivities one would expect in a warm neutron star crust have values of $\sigma\sim 10^{20-22}$ s -1 (Harutyunyan & Sedrakian, 2016). Such a conductivity would allow resolutions on the order of millimetres444We have $\Delta x>2\delta\nu/\sigma$. When redimensionalising we require a factor of $c^{2}$ on the RHS, thus $\Delta x>2\delta\nu c^{2}/\sigma$., whereas the highest resolution simulations to date (Kiuchi et al., 2014) have resolutions of tens of metres in the central region. Furthermore, the conductivities used in Carrasco et al. (2019) to model the exterior of a magnetar quote conductivities of $\sigma\sim 10^{17}$ s -1, resulting in a maximum resolution of metres. As a result, it appears that REGIME will be suitable for realistic neutron star simulations, and that for physical magnitudes of the conductivity it should be stable, even for the finest resolutions to date.

4.4 Performance

Knowing that REGIME gives sensible results for a range of simulations, and agrees well with resistive MHD in various regimes, we now turn our attention to the relative performance of the two models. The motivation for the model is that resistive MHD simulations generally require implicit time-integration schemes to ensure stable evolutions, and as a result suffer from slow evolution especially when near the ideal MHD limit (large $\sigma$). By instead using the REGIME source term, which is small in this limit compared to the resistive MHD source term, it is possible to perform the evolution using explicit time integrators, and ultimately simulations should benefit in terms of performance.

To determine this, we once again return to the Brio-Wu test, and to the magnetic reconnection problem. The Brio-Wu shock-tube is run using $N_{x}=400$ cells and is run until $T=0.4$. The 2D reconnection problem uses $(N_{x},N_{y})=(128,64)$ run until $T=2.0$. Due to time limitations the 3D reconnection problem uses $(N_{x},N_{y},N_{z})=(128,64,64)$, and is run until $T=0.5$. The problems are evolved using both models for a range of conductivities, and the optimum Courant factor giving the fastest execution, whilst generating similar results, is used. For the resistive model, we present the execution time using both the operator-split RK2 method and the SSP2(222) IMEX (Pareschi & Russo, 2004) integrator for comparison. Conductivities that have no corresponding execution time failed to run even with unreasonably small Courant factors.

Results are given in figure 10 that demonstrate the performance benefits of REGIME over resistive MHD. For the majority of the range of $\sigma$ that we present here, REGIME is the most efficient model to evolve for both problems, being only second to the resistive model (using explicit integrators) for the most resistive set-ups.

With these results, we can see that for any given resistivity there is no need to employ implicit integrators to evolve the system—either the system has a small enough $\sigma$ to allow use of the explicit integrators to evolve resistive MHD, or REGIME will provide a faster execution when implicit schemes are required.

5 Discussion

We present a resistive extension to the relativistic, ideal MHD equations often used in astrophysical simulations. The new source term is derived by expanding equations of motion about an equilibrium solution, resulting in a perturbation term that exhibits diffusive behaviour. The diffusion term scales with the resistivity, $\eta\equiv\sigma^{-1}$, and is numerically non-stiff when the resistive MHD equations are stiff. As a result, near the ideal limit, the source term may be evolved explicitly, resulting in speed-ups in execution times of many factors.

The new source term is demonstrated to produce similar results to resistive MHD in a range of initial set-ups, is able to capture resistive effects near discontinuous data without the onset of Gibbs oscillations, and shows little error growth for smooth, resistive solutions over conductivities spanning many orders of magnitude. For more complex simulations of magnetic reconnection, the correct scaling laws are produced for the reconnection rate.

Because of the nature of the source term, the explicit evolution of the new system results in relatively quick executions. Runtimes can vary dependent upon the conductivity but range from $2$ to $10^{2}\times$ faster than the optimum evolution using either explicit or semi-implicit integrators with resistive MHD. As a result, for many types of simulations, there is no need to employ IMEX integrators to ensure stability with resistive MHD. Near-identical results can be produced in a fraction of the time using the REGIME source term, combined with standard, operator split RK methods.

A limitation on the maximum resolution possible for a given simulation is given, and we show how the new source term can produce instabilities when this criterion is not met. We also give the form to predict at which resolutions such instabilities may occur for a given simulation, and conclude that REGIME should be stable for realistic resistivities in neutron stars.

The only significant difference between REGIME and resistive MHD is seen in the magnetic energy density power spectrum of the Kelvin-Helmholtz instability. Whilst the model still follows Kolmogorov’s 5/3 law for the kinetic energy power spectrum, an analysis of the magnetic energy density power spectrum shows that there are discrepancies between the expected result of resistive MHD and REGIME on small scales. These differences are not enough to affect the large scale behaviour of the magnetic fields, however, and the correct evolution of the maximum and average magnetic field is recovered.

One outstanding challenge for simulations of NS-NS and NS-BH mergers is to accurately model both the highly conducting region within the neutron star, and consistently match the electromagnetic fields to the near-vacuum exterior. Within the neutron star, ideal MHD is expected to be a good approximation. In the exterior, there are two possible limits of the resistive MHD equations: the force-free approximation ($\rho\to 0,\sigma\to\infty$) or the electro-vacuum ($\rho,\sigma\to 0$). In both cases we need to consider the conductivity $\sigma$ varying with space and time. In principle the REGIME approach extends directly to the varying-$\sigma$ case. However, by construction, it will become stiff in the electro-vacuum ($\sigma\to 0$) limit, and hence not practical. Current approaches, as proposed in (Lehner et al., 2012; Dionysopoulou et al., 2013), expect that in the regions around the merger there would be sufficient matter such that the force-free approximation would be valid. In this limit, therefore, REGIME should be sufficient in capturing the highly conducting plasma. If the $\sigma\to 0$ limit is more suitable then it would be more practical to match to a fully resistive model in the exterior, as it would not be numerically stiff in that case, whilst evolving the interior using REGIME. This would be analogous to the approach used in (Lehner et al., 2012).

In Section 3.2, to simplify the form of the inverted matrices in the source, we made the assumption that the fluid only couples weakly to the magnetic fields, and thus set $E=B=0$ for those calculations. In the tests we performed in this paper, this simplification gave only minor differences in the output of REGIME when compared to a full, numerical inversion. One would, however, expect these differences to grow if one modelled a magnetically dominated (low $\beta$) plasma.

Although all simulations have been performed in the special relativistic limit, the techniques we have used are not limited to this alone. A general relativistic extension to REGIME is under way, and will soon allow for more efficient, resistive simulations of neutron star mergers, accretion on to compact objects, and magneto-rotational instabilities.

The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work. We also gratefully acknowledge financial support from the EPSRC Centre for Doctoral Training in Next Generation Computational Modelling grant EP/ L015382/1. Open source software used includes Matplotlib (Hunter, 2007), Seaborn (Waskom, 2014) and CMINPACK (Devernay, 2017).