An operator-splitting numerical scheme for relativistic magnetohydrodynamics

We describe a novel operator-splitting approach to numerical relativistic magnetohydrodynamics designed to expand its applicability to the domain of ultra-high magnetisation. In this approach, the electromagnetic field is split into the force-free component, governed by the equations of force-free degenerate electrodynamics (FFDE), and the perturbation component, governed by the perturbation equations derived from the full system of relativistic magnetohydrodynamics (RMHD). The combined system of the FFDE and perturbation equations is integrated simultaneously, for which various numerical techniques developed for hyperbolic conservation laws can be used. At the end of every time-step of numerical integration, the force-free and the perturbation components of the electromagnetic field are recombined and the result is regarded as the initial value of the force-free component for the next time-step, whereas the initial value of the perturbation component is set to zero. To explore the potential of this approach, we build a 3rd-order WENO code, which was used to carry out 1D and 2D test simulations. Their results show that this operator-splitting approach allows us to bypass the stiffness of RMHD in the ultra-high-magnetisation regime where the perturbation component becomes very small. At the same time, the cod


Introduction
The strong gravity of astrophysical black holes and neutron stars creates some of the most extreme physical conditions in the Universe, which cannot be achieved in research laboratories.In particular, they naturally develop magnetospheres with extremely high plasma magnetisation.In contrast to the non-relativistic problems, where the magnetisation is well described by the ratio of thermodynamic and magnetic pressures (β = p/p m ), the magnetisation of relativistic plasma is best described by the parameter σ = B2 /4π(e + p), where B is the magnetic field strength as measured in the rest frame of plasma, p is the thermal pressure and e = ρc 2 +e t , where ρ and e t are the proper rest-mass density and thermal energy of plasma, respectively.In the magnetospheres of neutron start σ can reach the values of the order of 10 3 − 10 6 .Unfortunately, modern conservative schemes for relativistic magnetohydrodynamics (RMHD) begin to fail in multi-dimensional problems already when σ is of the order of few.Namely, the conservative variables cannot be converted to physically meaningful primitive variables.Higher values of σ can be handled for isentropic flows, where the energy equation can be eliminated [1].A fix based on the entropy obtained via parallel integration of the adiabatic entropy transport equation can help to bypass crashes for adiabatic flows as well [2].Obviously, this approach would not work for problems involving shocks, like the termination shock of pulsar winds [3], and, more importantly, magnetic reconnection, which involves kinetic and magnetic energy dissipation.Excessive artificial plasma heating due to numerical resistivity is another issue in simulations of highly magnetised flows.
In the limit of infinite magnetisation, σ → ∞, the system of RMHD becomes degenerate.It reduces to the Force-Free Degenerate Electrodynamics, where only two components of the energy-momentum equation are independent [4].It has been suggested, that this is the main reason for the shortcomings of conservative schemes in high-magnetisation problems [5].If so, it would make sense to apply a perturbation approach, where the FFDE solution is a leading order approximation and the plasma inertia enters only equations governing small perturbations to this solution.
A similar issue have been identified some time ago in Newtonian MHD simulations of the collision between highly magnetised (β ≪ 1) Earth's magnetosphere and the solar wind.In that problem, the Earth's largely dipolar magnetic field increases by many orders of magnitude from the collision site to the troposphere, whereas the perturbation of this field remains of about the same magnitude and hence increasingly small relative to the dipolar field on approach to the troposphere.As a result, the perturbation suffers large computational errors.To overcome this problem, Tanaka [6] proposed to separate the strong stationary dipolar field from its perturbation and numerically integrate only the nonlinear equations governing the perturbation.This approach has proven very productive [7].Here we describe a generalisation of this approach, where the zero-order solution is time-dependent.

Ideal relativistic magnetohydrodynamics
The covariant system of ideal RMHD consists of the Faraday equation the energy-momentum equation the continuity equation and the perfect conductivity condition In these equations, F µν is the Maxwell tensor and * F µν is the Faraday tensor, T νµ is the total stress-energy-momentum tensor, ρ is the rest mass energy of plasma particles, and u µ is the fourvelocity vector of macroscopic motion.In highly ionised plasma the Faraday tensor is Hodge dual to the Maxwell tensor, and hence * where e αβµν = √ −g ϵ αβµν is the Levi-Civita alternating tensor of space-time and ϵ αβµν is the four-dimensional Levi-Civita symbol.The total stress-energy momentum tensor, is the sum of the stress-energy momentum tensor of electromagnetic field and the stress-energy momentum tensor of matter Here p is the thermodynamic pressure, w(p, ρ) is the relativistic enthalpy per unit volume as measured in the rest frame of fluid, and g µν is the metric tensor of space-time.

Force-free and perturbation equations
Let us split the electromagnetic field into the force-free part f µν (0) and the correction part f µν The corresponding expansion for the stress-energy-momentum tensor of the electromagnetic field is where and The differential equations of the leading order system are the FFDE equations Since the four velocity of plasma does not enter these equations, the perfect conductivity condition is written in the form * which ensures the existence of inertial frames where the electric field vanishes.
The perturbation equations equations are obtained from the the full system of RMHD by removing the terms vanishing due to Eqs.13,14 It is important that this system is fully nonlinear with respect to its hydrodynamic component and hence allows the non-linear steepening of magnetosonic waves.In order to close this system we use the perfect conductivity condition By construction, the perturbation equations do not involve terms quadratic in the leading order electromagnetic field, like f µγ (0) f ν (0)γ , which would be dominant in problems with high σ.This is the main goal behind our approach.The terms linear in f µγ (0) in the perturbation stress-energymomentum tensor (12) describe the interaction with the FFDE field.

Numerical scheme
Since the FFDE equations are independent from the perturbation equations, even if initially the electromagnetic component of RMHD solution is close to the FFDE solution, they will eventually diverge.Once the solutions have diverged, the perturbation of the electromagnetic field is no longer small and the splitting becomes useless.To avoid this issue, the leading order and perturbation components of the electromagnetic field are recombined at the end of every full integration time-step, F µν = f µν (0) + f µν (1) and split at the beginning of the next time-step as f µν (0) = F µν , f µν (1) = 0.In other words, the electromagnetic field is first advanced as forcefree and then corrected for the interaction with plasma, making our approach a variant of the operator-splitting method.
We have explored the potential of this approach using a 3rd order WENO scheme for Special Relativistic MHD, similar to the one described in [8].In this explicit scheme, the solution is advanced forward in time using a 3rd order Runge-Kutta method.The FFDE equations 13-14 and the perturbation equations 16-18 are integrated simultaneously, without recombination of the electromagnetic field at the Runge-Kutta sub-steps.To set up the Riemann problems at the cell interfaces, we developed a novel WENO-interpolation which ensures rapid onset of 3rd convergence even for smooth solutions with local extrema.To compute the interface fluxes, we used an HLL Riemann solver [9], with the speed of light as the fastest wave speeds for the FFDE sub-system, and the fast-magnetosonic speeds for the perturbation sub-system.The conserved variables of the FFDE are converted into its primitive variables as in [4], where the energy equation is not utilised at all.The conserved variables of the perturbation sub-system are converted into its primitive variables via the approach similar to the one described in [8].Full details of the scheme will be presented elsewhere.

Test simulations
Here, we present the results of some test simulations.We use relativistic units where c = 1 and 4π does not appear in Maxwell equations, and the equation of state perfect gas w = ρ+Γp/(Γ−1) with the ratio of specific heats Γ = 4/3.

1D periodic Alfvén wave
The analytic solution to Alfvén waves of finite amplitude is given in [10].In the exact solution used for the test problem, ρ(x, t) = p(x, t) = 1, σ(x, t) = σ 0 , B x = 0.3B 0 and the tangential magnetic field B t = γ a B 0 e iψ(x,t) , with ψ(x, t) = arcsin(0.3sin(π(x + v a t))), where γ a = (1 − v 2 a ) −1/2 and the wavespeed v a = 0.5.B 0 is the parameter determining the magnetic field strength in the Hoffmann-Teller frame, it is equivalent to the magnetisation parameter σ 0 .Figure 1 illustrates the results for σ 0 = 545 and table 5.1 shows the results of the convergence study based on this problem, confirming the 3rd order accuracy of the scheme.For lower σ 0 , the results are almost as good.The relative L 1 error for B y increases only by the factor of about 4.5 for σ 0 = 5.45 × 10 −4 , suggesting that the approach is suitable for both the high-and low-σ regimes.

1D fast shocks
Here we present two tests, one for the high-σ regime and one for the low-σ regime.In both these cases, the angle between the shock front and the magnetic field in the rest frame of the upstream state is 45 o .The fast magnetosonic Mach number of the shock in the same frame is Shock problems are the most difficult for any numerical scheme.The lack of energy conservation in the FFDE sub-system is a particular source of concern specific to our scheme.The presented tests, however, demonstrate that the code can handle at least some shocks quite well.These and other tests, show good performance for shocks characterised by small variation of the tangential magnetic field.However, when this variation becomes strong the computational errors grow and the code may even crash.One example of these problematic shocks is the above high-σ fast shock when set up in the rest frame of the upstream state.However, the issue may not be as severe as it seems, as illustrated by the strong explosion test considered next.

2D cylindrical explosion
This problem is set on a uniform Cartesian grid with 200 cells in each direction.The computational domain is (−6, 6) × (−6, 6).The initial solution describes a uniform ambient gas at rest with p a = 3 × 10 −5 , ρ a = 10 −4 , B = (B 0 , 0, 0), and v = (0, 0, 0).In the central region with the polar radius r < 0.8 the gas pressure and density are raised to p 0 = 1 and ρ 0 = 0.01 respectively.A tanh-profile is used to introduce a gradual transition between the two states in 0.8 < r < 1.0.This is exactly the same setup as in [12], where B 0 = 1 was the highest value the code could handle (Later, with a somewhat modified code, the author of that paper failed to reproduce this result, as the code crashed.).We have had successful runs for B 0 = 0.01, 0.1, 1, 10, 100, and 1000.For B 0 = 0.1 the results look very close to that in [12] and other papers repeating this test (e.g.[8]).Figure 4 illustrates the results for the extreme cases of B 0 = 0.01 and B 0 = 1000; the corresponding magnetisation in the centre of the explosion, σ 0 = B 2 0 /w(p 0 , ρ 0 ), being σ 0 = 2.5 × 10 −5 and σ 0 = 2.5 × 10 5 respectively.For B 0 = 0.01, the magnetic field has a very small effect on the flow up to t = 4, and so it is almost the same as in the case of unmagnetised fluid.The magnetic field is swept aside by this flow.For B 0 = 1000, the magnetic field is so strong that the magnetic field lines are hardly deflected by the explosion and the fluid velocity is perfectly aligned with the x-axis.

2D ABC magnetic grid
The initial configuration of this problem describes a static uniform plasma distribution of density ρ 0 and pressure p 0 and the periodic force-free magnetic field which describes an equidistant grid of twisted magnetic field ropes aligned with the z-axis.In the xy-plane, these ropes appear as magnetic islands of opposite polarities (see the top-left pane of figure 5).Closest islands of the same polarities are separated via an x-point.This equilibrium configuration is unstable, leading to the eventual merger of islands with the same polarity via a collapse of the x-points into current sheets, followed by magnetic reconnection.Here we aimed to reproduce the FFDE and PIC simulations of this process described in [13].
The magnetic reconnection is accompanied by magnetic dissipation.In conservative schemes for RMHD, the total energy is conserved, and the energy of dissipated magnetic field is converted into heat.This is not true for numerical FFDE, where the energy of dissipated magnetic field is lost.When no measures are taken to combat this energy loss, the total energy conservation in this test problem is strongly violated.The approach we used to overcome this problem involves 1) integrating the energy equation of FFDE sub-system together with the momentum equations, 2) computing the difference between the energy obtained this way and the energy computed from the E and B fields of the sub-system, and 3) if it is positive, adding it to the energy of the perturbation sub-system at the end of every timestep.
The parameters of the simulations described here are B 0 = 1.0, p 0 = ρ 0 = 10 −4 .The Cartesian computational domain is (−1, 1) × (−1, 1), with periodic boundary conditions and a 200 × 200 uniform grid.Random noise was used to trigger the instability.In the test simulations, the linear phase of the instability continues up to t = 6.This phase corresponds to the initial plateau of the magnetic energy curve in figure 6.By t = 7, the xpoints separating the islands collapse into current sheets (see the bottom panels of figure 5).This opens the phase of rapid merger of neighbouring magnetic islands accompanied by strong magnetic dissipation.Once the number of the islands in the box reduces to two (of opposing polarities) the dissipation rate slows down significantly.By t = 30 almost 50 percent of the initial magnetic is converted into heat.The timescales of the evolution are similar to that seen in the PIC simulations (see figure 9 in [13]).Although we observe a somewhat smaller fraction of the dissipated magnetic energy, this could be related to the four-times larger size of the computational box in the PIC simulations, allowing it to trace the mergers to larger length scales.The total energy in the periodic box only slightly increases during the simulations.

Conclusions
In this contribution, we have described a rather unusual version of the operator-splitting method tailored to overcome the stiffness of RMHD equations in the regime of ultra-high magnetisation.The results of test simulations show that this approach allows to obtain accurate solutions for a wide range of problems involving smooth flows, shock waves, and magnetic reconnection.It is suitable not only for highly magnetised flows but also for flows with low magnetisation, which is important for many astrophysical applications.

Figure 1 .
Figure1.Alfvén wave test.In the exact solution, the wave moves to the left with the speed v a = 0.5.On both panels, the solid line shows the initial solution (t = 0) and markers show the solutions for t = 1, 2, 3, 4 (In the B y -plots, the solutions for t = 1, 3 and t = 0, 2, 4 are indistinguishable.).The resolution is n=80 and the Courant number Cu = 0.4.

Figure 2 .
Figure 2. 1D slow shock test.The magnetic pressure (left panel) and the gas pressure (right panel) are shown for the initial solution (t = 0, solid lines) and the numerical solution at t = 2 (markers).The exact shock speed is v s = 0.5.The disturbance located near x = −0.3 at t = 2 is caused by the evolution of the numerical shock structure at the beginning of the simulations.

Figure 3 .
Figure 3. Low-σ (left panel) and high-σ (right panel) fast shock tests.The magnetic pressure (left panel) and the gas pressure (right panel) are shown for the initial solution (t = 0, solid lines) and for the numerical solution at t = 2 (markers).The exact shock speed is v s = 0.5 in both the cases.

Figure 5 .
Figure 5. Collapse of the two-dimensional ABC grid.The top panels show the magnetic field lines and B z (left) and σ (right) in the initial solution.The bottom panels show the solution at the onset of rapid merger of magnetic ropes with the same polarity (t = 7): magnetic field lines and B z (left) and p (right).

Figure 6 .
Figure 6.Energy balance in the ABS test.The black solid line shows the total energy in the box in the units of the initial total energy.The red line is the fraction of energy in the electromagnetic field and the green line is the fraction of energy in the plasma.

Table 1 .
Convergence study based on the 1D Alfvén wave test.The numerical and analytic solutions are compared at t = 4.