Robust Recovery of Primitive Variables in Relativistic Ideal Magnetohydrodynamics

Modern simulation codes for general relativistic ideal magnetohydrodynamics are all facing a long standing technical problem given by the need to recover fundamental variables from those variables that are evolved in time. In the relativistic case, this requires the numerical solution of a system of nonlinear equations. Although several approaches are available, none has proven completely reliable. A recent study comparing different methods showed that all can fail, in particular for the important case of strong magnetization and moderate Lorentz factors. Here, we propose a new robust, efficient, and accurate solution scheme, along with a proof for the existence and uniqueness of a solution, and analytic bounds for the accuracy. Further, the scheme allows us to reliably detect evolution errors leading to unphysical states and automatically applies corrections for typical harmless cases. A reference implementation of the method is made publicly available as a software library. The aim of this library is to improve the reliability of binary neutron star merger simulations, in particular in the investigation of jet formation and magnetically driven winds.

Modern simulation codes for general relativistic ideal magnetohydrodynamics are all facing a long standing technical problem given by the need to recover fundamental variables from those variables that are evolved in time. In the relativistic case, this requires the numerical solution of a system of nonlinear equations. Although several approaches are available, none has proven completely reliable. A recent study comparing different methods showed that all can fail, in particular for the important case of strong magnetization and moderate Lorentz factors. Here, we propose a new robust, efficient, and accurate solution scheme, along with a proof for the existence and uniqueness of a solution, and analytic bounds for the accuracy. Further, the scheme allows us to reliably detect evolution errors leading to unphysical states and automatically applies corrections for typical harmless cases. A reference implementation of the method is made publicly available as a software library. The aim of this library is to improve the reliability of binary neutron star merger simulations, in particular in the investigation of jet formation and magnetically driven winds.

I. INTRODUCTION
General relativistic magnetohydrodynamic (GRMHD) simulations are an important tool to study many astrophysical scenarios involving neutron stars (NSs) and black holes (BHs). In particular, they represent the leading approach to investigate the dynamics of binary neutron star (BNS) and NS-BH mergers, which are the most important events in the nascent field of multimessenger astrophysics with gravitational wave (GW) sources [1].
Arguably one of the most pressing unsolved problems related to BNS and NS-BH mergers is to find the exact mechanism powering short gamma-ray bursts (SGRBs). The simultaneous detection of the gravitational wave event GW170817 and the SGRB named GRB 170817A [2][3][4], along with the following "afterglow" emission across the entire electromagnetic spectrum (e.g., [1,[5][6][7]), provided compelling evidence that BNS mergers can power SGRBs (e.g., [8][9][10]). In turn, this implies that the remnant object resulting from a BNS merger can launch, at least in some cases, a relativistic jet, which was indeed confirmed for GRB 170817A [9,10]. However, the jet launching mechanism and the nature of the object acting as central engine, either an accreting BH or a massive NS, remain uncertain (e.g., [11]).
Current simulations suggest that a mechanism based on neutrino-antineutrino annihilation would not be powerful enough to explain SGRBs [12,13], reinforcing the alternative idea that the main driver of jet formation should be a strong magnetic field. GRMHD simulations of BNS and NS-BH mergers, while considerably more complex and expensive because of the inclusion magnetic fields, become necessary to properly address the problem. Recent studies in this direction already provided important hints, supporting a scenario where the central engine is an accreting BH [14,15] while disfavoring the massive NS scenario [16]. Nonetheless, a final answer is still missing, and it will be necessary to overcome the tech-nical limitations of present GRMHD codes to ultimately solve the problem.
The merger event GW170817 was also accompanied by the kilonova transient AT 2017gfo, powered by the radioactive decay of heavy r-process elements synthesized within the matter ejected by the merger (e.g., [1,17]). Although this kilonova was observed in unprecedented detail, the interpretation in terms of specific ejecta components and their physical origin is still under debate. Also in this case, numerical relativity simulations represent the ideal approach to fully understand the different mass ejection processes occurring in a BNS (or a NS-BH) merger. Moreover, for some of these ejection processes magnetic fields are likely to play an important role (e.g., [18,19]) and therefore simulations should be performed in GRMHD.
The present work is devoted to a technical but crucial aspect of these simulations that has proven surprisingly difficult, and is motivated by the importance of GRMHD simulations in the context of BNS mergers (see, e.g., [20] for a recent review). Modern evolution codes are based on evolution equations written in form of coupled conservation laws for baryon number density, energy and momentum density including the electromagnetic contributions, and either magnetic field or vector potential. Primitive variables such as matter velocity, density, and pressure, are not directly evolved. Instead, they have to be recovered from the evolved quasi-conserved quantities after each evolution step.
While in Newtonian physics the above recovery is trivial, for the relativistic case one has to numerically solve a system of coupled nonlinear equations. The system involves also the equation of state (EOS), which encodes the behavior of matter up to supranuclear densities by specifying the pressure as function of density and temperature. An additional degree of freedom is the electron fraction, which effectively describes the matter composition and which can only change due to neutrino reactions. Since the EOS is not well constrained theoret-ically or by observation, a crucial requirement is the ability to perform simulations employing arbitrary EOS. This precludes closed-form solutions for the primitive variables, and the system has to be solved numerically. Since the solution is required inside the innermost loop of the evolution, computational efficiency is almost as important as robustness.
Note that most evolution codes make the simplifying assumption of ideal MHD. Although the electrical conductivity in merger remnants is very high, this approximation might not be justified for all aspects of the problem. On the other hand, evolving resistive GRMHD equations introduces even more difficulties (see also [21]). The equations for the primitive variable recovery are also very different for resistive GRMHD. Another complication is that in regions with strong magnetic fields but low mass density, movement of the matter becomes dominated by the field. Treating this "force-free" regime would in principle require different numerical evolution methods (for example, see [22]).
The simpler problem of recovering the primitive variables in relativistic hydrodynamics without magnetic fields is already solved in a robust manner, as described in [23] (also adopted in [24]). For the full problem of ideal GRMHD, several recovery methods with different limitations have been employed in GRMHD evolution codes [25][26][27][28][29][30][31][32]. Older schemes such as [30] are limited to particular analytic prescriptions for the EOS. Newer schemes can in principle work with any EOS, but not all implementations actually allow the use of arbitrary EOS. For a detailed review, we refer to [33].
All of the schemes investigated in [33] were shown to fail in certain regimes. While some of them work well enough in most of the regimes encountered during a merger simulation, even rare primitive recovery failures need to be handled and remain a common hurdle. An additional complication is that not all combinations of values for the evolved variables correspond to physically valid primitive variables. The occurrence of invalid evolved variables due to numerical errors of the evolution needs to be monitored and, if possible, corrected. If the recovery can fail also for valid input, it becomes impossible to reliably assess the overall validity.
In this work, we developed a new recovery scheme with the mathematically proven ability to always find a solution, and which is guaranteed to recognize invalid evolved variables. Furthermore, the scheme provides mathematically derived accuracy bounds. Our scheme is limited to the ideal MHD approximation, but does not introduce problems in the forcefree regime. We provide a reference implementation which is ready to use in any GRMHD evolution code, in the form of a C++ library named RePrimAnd [34]. Our implementation is written completely EOS-agnostic and provides a framework for EOS that can easily be extended. Since our aim is to improve reliability, we perform a comprehensive suite of tests.
The article is organized as follows: In Sec. II, we state the problem, derive the new scheme, prove the existence and uniqueness of a solution, and investigate the expected accuracy. In Sec. III, we discuss possible corrections to invalid evolved variables. In Sec. IV, we present numerical tests of our reference implementation, demonstrating robustness, efficiency, and precision. Here we also compare to other exist-ing schemes. Then, we study error propagation of evolution errors to the primitive variables in Sec. V, identifying potentially problematic regimes. Finally, in Sec. VI we draw our conclusions.

A. Primitive variables
Our scheme is designed for evolution codes which evolve variables defined on a spacelike foliation of spacetime from one timeslice to the next. The hypersurfaces and their normal observers define a frame we will refer to as the Eulerian frame.
We denote the 3-velocity of the fluid with respect to the Eulerian frame as v i , and the corresponding Lorentz factor by W . We will also use a quantity z ≡ W v. The baryon number density in the fluid restframe is denoted as n B . It is common to multiply n B with an arbitrary mass constant m B to define the baryonic mass density ρ = n B m B . The pressure in the fluid restframe is assumed to be isotropic and denoted as P . Denoting the fluid contribution to the total energy density in the fluid restframe as ρ E , we define the specific internal energy We further define a = P/ρ E and the relativistic enthalpy Note that the definitions of and h both depend on the arbitrary choice of m B . The primitive variables we use to describe the electromagnetic field are the electric and magnetic fields as seen by an Eulerian observer. In terms of the Maxwell tensor, where n is the normal to the hypersurfaces of the foliation, and the star denotes the Hodge dual. E µ , B µ are tangential to the hypersurface and thus equivalent to 3-vectors E i , B i . Our scheme neither requires nor provides the fields in the fluid frame, which can be obtained from the above using standard transformations.

B. Equation of State
We assume an equation of state (EOS) of the form The EOS could also depend on further variables, such as the electron fraction, as long as those variables are evolved variables or can be obtained from evolved variables in a trivial way, and can therefore be treated as fixed parameters in the primitive recovery algorithm.
For our purpose, it is also important to specify a validity range for each EOS. The validity range considers both physical and technical constraints. The most important physical constraint is the zero temperature limit for the internal energy. An example of a technical constraint is the range of values available for an EOS given in tabulated form. Currently, our scheme uses an EOS-dependent validity region specified in the following form However, it could easily be adapted to a more general shape in ρ, parameter space. We require that the lower validity bound min (ρ) is the zero-temperature value at the given density. This is the only meaningful choice for any numerical simulation involving cold matter at any time. Our error policy for correcting invalid evolved variables is based on this assumption, as is the proof for guaranteed success of the algorithm.
Our scheme relies on some physical constraints. Causality and thermodynamic stability require where c s is the adiabatic speed of sound, given by We assume positive baryon number density and positive total energy density, which implies We assume that the pressure is positive and further bounded by the total energy density (dominant energy condition), which implies For a given EOS, we also require a positive lower bound h 0 for the relativistic enthalpy h, such that 0 < h 0 ≤ h(ρ, ) over the entire validity region of the EOS. Note we do not assume > 0 or h ≥ 1. The definitions of and h depend on the arbitrary choice of the mass constant m B . Unless the latter is fine-tuned to the average baryon binding energy at low density, nuclear physics EOS often yield slightly negative . In practice, h 0 is of order unity.
By design, our scheme is not confined to any particular equation of state. It only uses the information defined above and does not make any other kind of EOS-specific distinctions or adjustments. For the purpose of testing our scheme, we use two specific EOS as examples: 1. A hybrid EOS given by a cold part and a simple thermal part For the cold part, we employ MS1 EOS from [35] (based on [36]), and we use Γ th = 1.8. This hybridized type of EOS is often used in numerical relativity.
2. The classical ideal gas, given by Here, we use an adiabatic exponent Γ = 2. Pressure and internal energy are zero at zero temperature. We use this unrealistic model, where pressure is only given by thermal effects and degeneracy pressure is ignored, as a corner case for testing our algorithm. Note that in numerical relativity, ideal gas EOS normally refers to a hybrid EOS based on a polytrope with adiabatic exponent Γ = Γ th as the zero-temperature limit.

C. Evolved Variables
Our scheme is intended for numerical evolution codes employing evolution equations for energy, momentum, and baryon number density formulated as a quasi-conservation law. This is the case for finite-volume shock-capturing schemes. The evolved quantities are called conserved variables, although only the baryon number is strictly conserved. The fluid contribution is given bȳ Including also EM contributions, the evolved variables are given by In most formulations, the evolved variables above include the volume element of the spatial metric. Since this factor is available from the spacetime evolution, it is not relevant for the following and was left out of the definitions. In addition to the evolved variables D, τ, S i , we assume that the magnetic field B in the Eulerian frame is either an evolved variable or can be reconstructed from evolved variables without knowledge of the fluid-related primitive variables, such that B is known when our primitive reconstruction scheme is run. This is the case for state of the art methods such as constrained transport schemes or schemes evolving the vector potential.
We only consider evolution codes that assume the ideal MHD limit, which corresponds to the additional condition This precludes the use of our scheme for resistive MHD simulations.
For convenience, we define rescaled variables as We also need to decompose the momentum in parts parallel and normal to the magnetic field The decomposition is undefined for the case of zero magnetic field, but we will exclusively use the product with b 2 , which is always well-defined.

D. Useful Relations
In the following, we collect definitions and analytic relations for later use. First, we define two quantities that will play a central role, Trivially, their ranges are limited to Given the conserved variables and µ, one can directly express the primitive variables analytically. Since the system is over-determined, there are different possible expressions which may disagree if the given µ is inconsistent with the conserved variables. In the latter case, some expressions can diverge. We will use the ambiguity to cast the recovery into a root-finding problem, by expressing the same variables in different ways that only agree for the correct value of µ.
As an intermediate step, we first remove the electromagnetic part of the conserved variables. From Eq. (21) This allows to compute the pure fluid part of the conserved variables. The corresponding quantities relevant for our purpose can be written as We can now express the velocity as v = µr. This expression does however not guarantee that v < 1 for any µ. One way to avoid exceeding the speed of light is to use the quantity z, which yields Although we do not have a closed form expression for z as function of µ, we can use z to obtain a useful upper limit for the velocity, given by After obtaining velocity and Lorentz factor, we can extract rest mass density and specific internal energy using the expressions If ρ and are in the validity range of the EOS, we can now compute pressure P = P (ρ, ) and the enthalpy h(ρ, ). Finally, the following expression for µ itself turned out useful

E. Designing the master function
In the following, we formulate the primitive variable recovery as a root finding problem for a suitable master function. To this end, we employ the following design goals 1. The function should be one-dimensional.
2. It should be continuous.
3. It should always have exactly one root, even for unphysical values of the conserved variables.
4. It should be well behaved in the Newtonian limit.
5. It should be well behaved for zero magnetic field.
6. There should be a known interval which contains the root and on which the function is defined.
7. The root-finding procedure should not require derivatives of the EOS.
For our scheme, we use µ defined in Eq. (25) as the independent variable to solve for, i.e. we will construct a function f (µ) which crosses zero where µ is consistent with the conservative variables. The latter take on the role of fixed parameters. To construct f , we start with Eqs. (29) and (30) and definer Next, we define functions for velocity and Lorentz factor v(µ) = min(µr(µ), v 0 ) ,Ŵ (µ) = 1 where v 0 is the upper velocity limit from Eq. (32). Further, we define rest mass density and specific energy according to Eq. (33) and Eq. (34) provided that the results fall within the validity range of the EOS. Otherwise, the densityρ is adjusted to the closest value within the validity range for ρ, andˆ to the closest value within the validity range for at adjusted densityρ. In Eq. (41), we expressed the term W −1 in a way that prevents large rounding errors in case of small velocities. Using the EOS, we compute the pressure, defininĝ To close the circle, we could now define an enthalpy functionĥ(µ) in terms ofˆ (µ),P (µ) and then express µ itself as functionμ(µ) based onŴ (µ)ĥ(µ). However, we found that this straightforward choice does not yield a function respecting our design goals. One reason is that under extreme conditions, the strong limitations introduced forP ,ˆ ,ρ, andŴ can cause severe kinks in the function. We found a better choice based on Eq. (36). In the latter, we compute the variable ν ≡ h/W in two slightly different ways based on Eqs. (2) and (34), defining The second form computes (1 + ) /W directly from Eq. (34) and uses the range-limitedˆ only when computingâ. Using Eq. (36), we finally obtain our master function f (µ) Examples for the master function resulting for different regimes are shown in Fig. 1. As one can see, the master function is almost linear unless Lorentz factor or magnetization are large, but remains well behaved even then. Note that we only show the root bracketing interval that is constructed in the next section. Beyond this interval, the function can have strong kinks.

F. Existence of solution
In the following, we prove that the master function always has a root, not just for valid evolved variables, but also for invalid ones. To this end, we construct an interval over which where th denotes the difference to the zero temperature case. Velocities are oriented orthogonal to the magnetic field, which is the most difficult case. In the top panel, the independent variable µ is scaled to the initial root bracket µ+, and the function to the maximum value over the interval shown. The lower panel shows the behavior near the rootμ0. the master function changes sign. We start by defining an auxiliary function This is a smooth function which does not require evaluation of the EOS, using only the EOS-specific global lower enthalpy bound h 0 . It is easy to show that f a is strictly increasing and has exactly one root µ + in the interval (0, h −1 0 ]. We will show that µ + constitutes an upper bound for the root of the master function f . For f a (µ + ) = 0, we find where we used r ≥r and the monotonicity of the above expression with respect tor. From Eq. (39), we findv(µ + ) = µ +r (µ + ) ≤ v 0 < 1. Further, f a (µ + ) = 0 implies that Using the definition Eq. (45) ofν, we can write Hence, Using the definition of the master function, Eq. (46), we find that f (µ + ) ≥ 0, and, trivially, f (0) < 0. Since f is continuous, it has at least one root in the interval (0, µ + ]. Conveniently, the interval provides a useful initial bracketing for the root finding algorithm. Although finding µ + requires another numerical root solving, the computation of f a does not require the expensive evaluation of the EOS. Moreover, since f a is smooth, efficient Newton-Raphson methods may be employed. The main reason to expend this effort is to ensure that v ≤ v 0 < 1. Beyond µ + , the cutoff in Eq. (39) can induce a strong kink in the master function, reducing efficiency of the main root finding. With the tight upper limit µ + , the only reason for the cutoff is to make absolutely certain that not even rounding errors in ultra-relativistic cases can lead to not-anumber results. Finally, being able to usev ≤ v 0 simplifies the uniqueness proof in the next section.

G. Uniqueness of solution
Uniqueness of physically valid solutions is obviously important for any evolution scheme based on the conserved variables considered in this work. For the purpose of our recovery scheme, we require in addition that (i) for valid evolved variables, the master function has no additional roots corresponding to invalid solutions, and (ii) for invalid evolved variables, it still has exactly one root. In this section, we will prove that the above conditions are met.
We first compute the derivative of the master function. Differentiation and straightforward algebraic computations yield Sincev is monotonic, and we have shown in Sec. II F that v(µ + ) ≤ v 0 , it follows thatv ≤ v 0 is always satisfied at a root of f . For the case that computed by Eq. (41) is in the valid range of the EOS, we find At a solution, µŴĥ = 1, and we obtain This means that changes of densityρ and specific energyˆ are adiabatic when varying µ near a solution. For the case that computed by Eq. (41) is below the valid range of the EOS,ˆ is set to the lower bound of the validity range, min (ρ), which is the zero temperature limit according to our EOS requirements. Therefore, the variation with µ is also adiabatic. In both of the above cases we can use Eq. (8), obtaining whereĉ s = c s (ρ,ˆ ). At a solution, the derivative of the master function becomes The requirement c s < 1 is therefore sufficient to guarantee uniqueness of the root for all velocities and magnetic fields.
So far, we did not address the corner case where the specific energy is above the valid range of the EOS. In that case,ν is computed from Eq. (44) (this choice is intended to improve the behavior away from the solution). A straightforward computation reveals that uniqueness is always ensured under the condition that The function A is related to the upper validity range max of the EOS. It is defined by the relation and corresponds to the change of specific entropy along the validity boundary, where A = 0 means adiabatic change. The condition given by Eq. (65) does not seem to be very restrictive. In doubt, one can always use a boundary with A = 0 to guarantee uniqueness in all cases. In practice, we encountered no problems using upper validity bounds defined by either constant temperature or constant .

H. Guaranteed accuracy
Since the root of the master function is determined numerically, we require a criterion to stop the iteration once sufficient accuracy is reached. What is sufficient depends on the other errors present in a numerical evolution scheme. We will discuss evolution errors in Sec. V. In this section, we discuss the error propagation of the root finding accuracy to quantify the accuracy of the recovered primitives.
However, we first need to specify how the final result is computed from the outcome of the last root finding iteration. This involves a design decision, since the available variables µ,μ,ν,r,q,v,Ŵ ,ρ,ˆ ,ĥ,P allow to compute the primitives in many different ways, which lead to different error propagation. Here, we useŴ ,ˆ ,ρ,P directly, which turns out to be a good choice in terms of error propagation. To reconstruct the velocity vector, we use the expression which is just Eq. (28) rearranged to avoid degeneracy for the case b = 0. It is easy to verify that the Lorentz factor corresponding to Eq. (67) is exactlyŴ . The conserved densityD computed from the recovered primitives agrees exactly with the original one.
In the following, we only consider the case where the solution is in the validity region of the EOS. For invalid solutions, the accuracy of the solution is less relevant since in this case the cause is the evolution error and the result will either be corrected to the valid range or the simulation aborted. The error introduced by such corrections will be discussed in Sec. V D.
Assuming the root of f (µ) was determined numerically to accuracy of δµ, we now estimate the resulting accuracy of the primitive variables. We already evaluated the relevant derivatives at a root of f in Sec. II G. From those, we obtain where |δv i | denotes the norm given by the 3-metric of the vector δv i . The error in the recovered primitive variables corresponds to an error of conserved variables. We compute this error by insertingρ,ĥ,P ,v i into equations (15)(16)(17)(18)(19)(20) and (21) and then computing the derivative with respect to µ, obtaining the following scaling We find that the accuracy in µ required for a fixed relative error of the primitives increases with increasingly relativistic velocities. A high magnetization, on the other hand, has no impact on the error bounds. It is also worth noting that the error δs of the the specific entropy s is zero (to linear order in δµ) because the variation ofρ,ˆ with respect to µ is adiabatic (see Sec. II G). Finally, we note that the above error bounds do not include numerical rounding errors. Those will be discussed in Sec. IV B.

III. ENFORCING VALIDITY
In typical numerical simulations, the evolved magnetohydrodynamic variables frequently reach an invalid state at some points, mainly due to ordinary numerical error, but also external influences such as gauge pathologies near the center of newly formed black holes. Often, such violations are harmless and can be corrected. Any such correction turns unphysical conditions into regular evolution errors, and obviously different prescriptions will lead to different errors, both in magnitude and in character. Although correcting violations should be regarded as part of the evolution scheme, some basic pointwise corrections can be incorporated into the primitive recovery code, granting it power to change the evolved variables. The following effects cause typical harmless violations: 1. When evolving zero-temperature initial data, arbitrary small evolution errors can lead to evolved variables that correspond to a fluid energy density below the zerotemperature limit.
2. At numerical grid points at the surface of neutron stars moving through vacuum, mass and energy densities during a single timestep can drop by orders of magnitude or even become negative. Although the absolute errors of the conserved variables remain small compared to the global scales of the system, the resulting local error of the specific internal energy and velocity can become huge and lead to an invalid state. The effect is alleviated over time because the errors tend to heat the outermost layer of NS surface, creating a hot atmosphere that reduces the density gradient.
3. During collapse to a black hole, mass density and/or temperature might leave the range covered by the given EOS, arriving at a state that is not unphysical but cannot be evolved further. This typically occurs in regions already inside the horizon or about to be engulfed by a rapidly expanding apparent horizon.
4. The coordinates near a black hole center are strongly stretched for gauges like the puncture gauge, and the surroundings are extremely under-resolved numerically. Under those conditions, all kinds of numerical instabilities can occur for the combined magnetohydrodynamical and spacetime evolution system.

A. Simple corrections
By design, our primitive variable recovery scheme is able to deal also with invalid input. As a side-effect, we obtain a projection onto the valid regime, by simply recomputing the evolved variables from the recovered primitives. As described in Sec. II E, the scheme always yields a pairρ,ˆ such thatˆ is within the validity range of the EOS atρ.
We first consider the important case that the raw value of is below the valid range. In this case, only the recomputed conserved energy τ changes, while S and D stay the same. This can be seen as follows. The only variable through which the adjustment of to the valid range impacts the master function isν. For the case at hand, Eq. (45) impliesν = ν A . Furthermore, the conserved energy τ enters exclusively through Eq. (41). Therefore, if τ is adjusted such that Eq. (41) yields the range-limited value forˆ , we arrive at the same primitive variables without adjustment.
For the case that the energy is above the validity range of the EOS, all recomputed conserved variables can change. One could prevent this by always usingν = ν A , but not without changing the behavior of the master function away from the solution. However, this case is less important, because this correction should only be allowed at low-density fluidvacuum boundaries (NS surfaces) or inside horizons.
In the interior of black holes, it becomes necessary to employ a more lenient error policy than outside. Although physical effects cannot propagate out of the horizon, violations of the constraint equations and gauge effects impact the exterior. Therefore, one cannot allow any runaway instability inside the horizon. For the matter part, this mainly concerns energy and momentum, since the total baryon number is conserved in finite volume schemes (artificial atmosphere corrections aside), and the mass density remains finite. The energy can be limited by allowing the aforementioned correction to the EOS range inside horizons even at high densities.
This leaves the momentum. For pure hydrodynamic simulations, limiting the velocity proved effective to prevent runaway instabilities near the BH center. This was employed for the simulations in [23], by rescaling the velocity to stay within a given limit. For MHD simulations, this approach has a sideeffect. Since the reconstructed electric field depends on the velocity via Eq. (21), it will also change. That might be problematic or not, depending on the evolution scheme. The evolution of the EM field might be problematic in this regime in any case. However, addressing such problems is clearly not inside the scope of the primitive variable recovery, since it operates point-wise and cannot change electric or magnetic fields in any reasonable way.
Another correction often applied is to enforce a minimum mass density, also called artificial atmosphere. There are two motives. One is the wish to use a tabulated EOS that does not include zero density (this might be achieved more consistently by extending the range to zero via analytic expressions). The more fundamental motive is that the hydrodynamic evolution equations break down in vacuum. In purely hydrodynamic simulations, it is common to set the atmosphere velocity to zero with respect the the simulations coordinate system, in order to prevent an unphysical influx of matter. In ideal MHD simulations, the situation is more complicated because the electric field is tied to the velocity via Eq. (21). Therefore, the atmosphere velocity prescription should be the domain of the evolution scheme and not of the primitive recovery.

IV. PERFORMANCE
In the following we assess how well our scheme performs in practice. For this, we devised stand-alone tests that encompass typical parameters encountered in numerical relativity and binary neutron star mergers in particular. Those tests are relevant for the future use of the scheme in actual numerical simulation codes.

A. Code Design
We created a reference implementation for the algorithm described in Sec. II in form of the C++ library RePrimAnd. The library is not tied to any particular evolution framework, allowing the use in arbitrary evolution codes. It also contains a framework providing access to different types of EOS through a generic interface, ensuring that the user code (such as our primitive recovery code) is completely EOS-agnostic. The generic interface also provides the EOS validity ranges and rigidly enforces our EOS requirements (see Sec. II B). The reference implementation is publicly available [34].
In order to find the root of the master function, our implementation uses the TOM748 algorithm [37] provided by the BOOST library. This method keeps the solution bracketed, is guaranteed to converge in a limited amount of steps, and does not make use of function derivatives. The motivation for avoiding derivatives is that, in practice, tabulated EOS tend to have very inaccurate partial derivatives, which is problematic when using a derivative-based root solver such as Newton-Raphson. Our implementation therefore does not make use of the soundspeed or other derivatives. The root is determined to an accuracy specified by prescribing ∆ defined in Eq. (68), where the error δµ is taken as the size of the current tightest bracket of the root.

B. Robustness and Accuracy
To test the algorithm and the implementation, we created a comprehensive testsuite. For the main test, we sample the primitive variable parameter space given by density, temperature/specific energy, magnetization, and velocity. We sample z = W v between 0 and 1000, the magnetization b from 0 to 5, and the specific thermal energy from th = 10 −4 up to 50. For the MS1 EOS, we sample the mass density from 10 6 to 10 15 g/cm 3 . For the ideal gas EOS, the mass density is irrelevant due to the scaling behavior of the EOS. We use two orientations of the velocity, parallel and orthogonal to the magnetic field. The tests are performed both for the ideal gas and for the hybrid EOS, described in Sec. II B.
We verified that the algorithm always, without any exception, succeeds in recovering a solution for valid input. Moreover, we created test cases to assure that input corresponding to energy outside the range possible for a given EOS is correctly classified as such.
To assess the accuracy, we compute the conserved variables from the primitives, apply the primitive recovery algorithm, and compare the result to the original primitives. Further, we compute the conservatives from the recovered primitives and compare to the original conservatives. Our testsuite compares the observed accuracy to the one expected from Eq. (68) to Eq. (72), and Eq. (73). When demanding an accuracy better than ∆ 10 −7 , those bounds are exceeded either for high Lorentz factors or very small and v. We attribute the excess error to various rounding errors.
We identified the most important rounding errors as follows. First, the master function is the difference of two values which can each be expressed only to machine precision. To get the impact on the root, we have to divide by the derivative of the master function, which in this case satisfies f ≤ 1 − v 2 c 2 s . At the same time, we demand an accuracy ∆/W 2 . For the highly relativistic case W = 10 3 and around 16 digits machine precision, this limits ∆ > 10 −10 . If the soundspeed approaches unity, the accuracy is further limited. Second, is computed by subtracting kinetic and magnetic energy density from the total one. If is small compared to these, the cancellation error causes a loss of significant digits. Analyzing Eq. (41), we find additional rounding errors of magnitude z 2 / and b 2 W/ worse than machine precision.
Taking into account both the regular errors predicted by Eq. (68) to Eq. (72) as well as the main rounding errors discussed above, the recovered accuracy is quantitatively within the expected bounds over the whole range of our test cases. Fig. 2 shows the recovered accuracy for the pressure as well as the boundary where the errors caused by rounding start exceeding those caused by root finding.
We do not expect rounding errors to be of practical importance. The rounding errors at low , v are very small and only dominate because the regular errors approach zero. The rounding errors in the ultra-relativistic/highly magnetized regime are still not prohibitive, but will be dominated by the errors of the time evolution, which will be discussed in Sec. V.

C. Efficiency
We measure the computational efficiency of our scheme in terms of calls to the EOS. The motivation is that for a tabulated EOS including thermal and composition degrees of freedom, a single EOS call is likely more expensive than the evaluation of the analytic expressions within our recovery scheme. The worst scenario is when the EOS is tabulated with temperature as one independent variable. Each EOS call then requires its own root solving to convert from to T . Fig. 3 shows how the efficiency varies with specific energy and velocity, either for zero magnetic field, or with a fixed high value for b. Similarly, Fig. 4 shows the efficiency with respect to velocity and magnetization, both for cold and for hot matter. We find that the efficiency does not degrade even for Lorentz factors up 1000 and extremely high magnetizations. The absolute maximum number of calls to the EOS to achieve an accuracy ∆ = 10 −8 is 23, when considering the whole parameter space used in the unit tests (not just the cuts shown in the plots) and both EOS types. The maximum occurs for the ideal gas and only when both > 40 and b > 2, i.e. thermal energies much larger and magnetic energies larger than the rest mass density. Note that the extremes reached in our tests are rather pathological scenarios where the description of matter and the ideal MHD approximation likely break down. For a typical binary neutron star merger simulation, most matter has much lower velocities, temperature and magnetization. In practice, we expect an average number of required calls below 10.

D. Comparison with other schemes
In the following, we compare our scheme to existing ones. We refer to [33] for a comprehensive review and numerical tests of previous schemes. The main characteristics are listed in Table I. One important difference is the number of independent variables. Most of the existing schemes need to solve an equation in two or three unknowns. This is a severe drawback. First, it is difficult to ensure the solution is found. The Newton-Raphson (NR) schemes might not converge. Second, robust but fast schemes that guarantee finding the solution in a limited number of steps only exist for one-dimensional root finding. Third, the recovery schemes based on NR require an initial guess, which is typically taken from the previous timestep during numerical evolution. This makes the methods more unpredictable and more difficult to test, as they do not depend on the conserved variables in a deterministic  way. As two of the existing schemes, our scheme is using one-dimensional root finding. Further, it also makes use of a tight initial bracketing interval proven to contain exactly one solution.
As demonstrated in [33] all of the existing schemes can fail for Lorentz factors 10-1000, depending on the magnetization, and even for small magnetization (compare Fig. 3 therein). For our scheme, the existence and uniqueness of the solution are proven analytically. We also tested our scheme numerically up to extreme magnetization and Lorentz factors 1000, with magnetic pressure not just dominating the fluid pressure, but even the matter total energy density.
Strong magnetization is important for studying the engine of short gamma-ray bursts (SGRBs), where ideally a very low density matter is subject to very strong magnetic fields. This regime is also problematic for the numerical time evolution itself. The ability of our scheme to distinguish reliably between valid and invalid evolved variables is therefore an important advantage.
Comparing efficiency and accuracy to the results in [33] is difficult, since both depend on the tolerance criterion used in the root finding, and because the errors in the different primitive variables are not related in the same way for different recovery methods. For the case of Lorentz factor W = 2 and low magnetization P/P mag = 10 −3 , the recovery accuracy shown in Fig. 2 of [33] seems sufficient for most applications, as is the case for our scheme (compare Fig. 2). However, some existing schemes exhibit isolated points where accuracy degrades intolerably even for this easy scenario. Our scheme has the advantage of well defined theoretical error bounds for each of the primitive variables that hold even for high Lorenz factors and magnetization, and our implementation contains a test suite to verify those over a large parameter space.
Another difference of the schemes is the required form of the EOS. Our scheme is using P = P (ρ, ). It was argued in [33] that it is advantageous if the root function uses an EOS of the form P (ρ, T ), as this does not require another root finding to determine T for each EOS evaluation. However, this is only true if the EOS is a table with T as independent variable. It would be better to simply create a lookup table in terms of the variables provided by the recovery scheme. This allows to chose the most robust recovery procedure without sacrificing speed. For an EOS tabulated in terms of T , however, the scheme in [33] might indeed be faster than ours.
Regarding the efficiency, the different tolerance measures allow only a rough comparison, using the number of root finding iterations shown in Figs. 1 and 3 in [33]. Only the Newman scheme appears to be consistently requiring less than around 10 steps, unless it fails completely. The others need 30 iterations or more in certain regions of parameter space, provided they succeed at all.
Our scheme is guaranteed to converge in a finite number of steps because the root finding algorithm performs bisection steps if needed. The efficiency does not degrade for large Lorentz factors or strong magnetization, in contrast to most other schemes. The worst case scenario for our scheme seems to be extreme temperatures. Even for almost photonic states, it does not require more than 23 EOS calls, however.

V. IMPACT OF NUMERICAL ERROR IN EVOLVED VARIABLES
In the following, we investigate consequences of numerical errors in the evolved variables in conjunction with the corrections of invalid states described in Sec. III. Further, we identify regions in parameter space where the primitive variables are particularly sensitive to errors of the evolved ones. I. Main characteristics of different recovery schemes. We list the independent variables used in the root finding (translated to our notation), the variables for which the EOS needs to be evaluated, whether the scheme can fail, whether the formulation allows a bound on the number of iterations needed for finding the solution, and whether the scheme requires to provide an initial guess for the solution.

A. Newtonian Limit
It is instructive to consider the relation between evolved and primitive variables in the Newtonian limit. Assuming that both kinetic and thermal specific energies are nonrelativistic corresponds to v 1, 1, a 1, h ≈ 1 (for simplicity we chose m B such that h 0 = 1 in this subsection). To leading order in v 2 and , we obtain Taking the Newtonian limit locally does not imply small b. However, if the magnetic field energy is comparable the rest mass density one cannot expect the velocity to stay nonrelativistic during the course of the evolution. It is a plausible assumption that the density of kinetic energy is not much smaller than the magnetic energy density.
, we can also neglect the electric contribution b 2 v 2 ⊥ to . On a side note, it is easy to show that the master function Eq. (46) becomes a linear function in the Newtonian limit (with b 2 1). As x(µ) ≈ 1,r(µ) andq(µ) are independent of µ and equal to the correct values. The same holds in turn forρ,ˆ ,P . Further,ν ≈ĥ ≈ h 0 . Sincer 2 1, the master function becomes f (µ) ≈ µ − 1.
We now turn to the propagation of the evolution error of the variables q, r, D. Even in the Newtonian limit, both v 2 and b 2 contributions can dominate q if is even smaller. Since v is essentially computed from r, computing from the evolved variables suffers from cancellations that amplify the evolution errors. In detail, Once δ / exceeds unity, reconstructing from the evolved variables might lead to larger errors than simply setting it to the zero-temperature value. Assuming some bound for the relative errors of the evolved variables, this corresponds to critical values for b 2 and v 2 .

B. Magnetically Dominated Regime
In the context of magnetohydrodynamic evolution, magnetically dominated refers to the magnetic pressure exceeding the fluid pressure. Increasing the field strength at fixed matter density, the movement of matter becomes constrained along the field lines at some point.
This effect is also reflected in the equations we use for primitive recovery. The relation between total and fluid momentum S ⊥ components perpendicular to the magnetic field isS ⊥ = xS ⊥ , as seen from Eq. (28). The quantity x depends only on µb 2 . In the limit µb 2 1, we find x 1. In that case, the perpendicular part of the evolved momentum is dominated by the electromagnetic part. However, the latter is proportional to v ⊥ in ideal MHD, and also points in the same direction. Therefore, the evolution error of the perpendicular part of momentum is not amplified by cancellations when recovering the fluid velocity orthogonal to the field. Also the parallel part of the evolved momentum is not problematic since it has no electromagnetic contribution.
As discussed for the Newtonian limit, strong magnetic fields are also detrimental to the accuracy of since evolution errors in B are amplified by cancellation. From Eq. (38), we find that the magnetic field contribution to the energy causes strong cancellation error ifq b 2 (1 − v 2 ⊥ ).

C. General case
To quantify the error amplification in the general case, we compute the partial derivatives of the primitive variables with respect to the conserved ones (by means of finite differences). As expected, specific internal energy and pressure exhibit large error amplification in some regimes, while the fluid momentum is well behaved. Fig. 5 shows the behavior of amplification factors of specific energy error in relation to errors in the evolution of q, b, and r, defined as The error of the pressure shows the same qualitative behavior. For a magnetar-strength field B = 10 15 G, we find that a relative evolution error of δb/b = 10 −4 would start to dominate the evolution of at densities of magnitude 10 8 g/cm 3 . The same holds for a relative error δq/q = 10 −4 , which is to be expected since q is dominated by the b 2 contribution. This regime could be relevant for the engine of SGRBs, as a popular model assumes a low-density funnel along the rotation axis of a black hole immersed in a strong magnetic field which is anchored in a surrounding disk. Similarly, the material surrounding a supramassive (i.e. long-lived) neutron star merger remnant could be affected. The consequence of artificial heating could be artificial outflows and increased neutrino luminosity. To assess the potential for spurious winds, we can compare the scales of additional specific energy caused by the evolution error and the specific gravitational binding energy. At a distance 100 km to a M = 2M remnant, we find that the thermal error starts to dominate at densities 10 8 g/cm 3 , again for a fiducial evolution error δb/b = 10 −4 and B = 10 15 G.
On the other hand, the above discussion is overly pessimistic if the increase in thermal energy by physical effects, such as shocks or neutrino absorption, greatly exceeds the one by numerical errors. In other words, the presence of a mild outflow caused by mild heating should be met with more skepticism than stronger outflows caused by prominent heating.

D. Interaction with recovery corrections
Enforcing the evolved variables to stay in the physically valid regime corresponds to a projection onto the validity boundary. Depending on the choice of projection, it is possible that the corrections cause a drift along the boundary, in the worst case with a preferred direction. This is of particular concern for the very frequent correction of limiting the specific energy above the zero temperature value. For our scheme, only the energy density is corrected in that case. Therefore, it does not introduce a drift of the evolved momentum density.
The main effect of limiting above the zero temperature value may be to induce a spurious heating. The reason is that cutting the evolution error distribution from below creates a positive bias until the error distribution has little support below the cut. Of course, the raw error distribution before the correction could already contain a bias. For the idealized case that the evolution errors follow a zero-mean normal distribution around the correct result, we expect the temperature to increase until the thermal energy reaches a level comparable to the width of the error distribution for total energy. Note that excessive artificial heating could reduce the velocity, since the momentum density incorporates a factor h. However, if h is significantly increased non-homogeneously by the errors, the corresponding changes in thermal pressure can be expected to cause gradients and corresponding acceleration.
In the above discussion, we omitted the effect of finite root solving accuracy. Our implementation of the algorithm recomputes all conserved variables from the primitive ones if, and only if, corrections were required. Therefore, the momentum only remains constant to the accuracy of the root solving when applying the correction to the energy. This error is formally bounded by Eq. (73). We cannot predict, however, if the distribution limited by this bound is symmetric or not. Any bias might lead to a cumulative effect over many corrections. For the worst possible scenario where each correction leads to the maximum possible momentum error always pointing along the momentum, a few thousand corrections could add up to intolerable levels.
To assess the actual behavior, we performed a numerical experiment. Instead of using a full numerical evolution, we employ a random walk model representing an evolution error, starting at selected states. After each "evolution" step, we apply the primitive recovery and limit the conserved variables to the allowed regime. The cumulative corrections applied to energy and momentum are monitored.
In this approach, we can prescribe the error distribution. As a worst case example, we use a normal distribution with negative mean for the energy error. Starting at a zero temperature state, this causes frequent corrections to the energy. Note the expected error in the momentum does not depend on the magnitude of the corrections, but this randomized test nevertheless involves different magnitudes. For the root finding accuracy, we use four different values ∆ = 10 −7 , 10 −8 , 10 −9 , 10 −10 .
We find that the average momentum error is orders of magnitude smaller than the limit Eq. (73). Selecting an initial state v = 0.99, b = 2, ρ = 6 × 10 12 g/cm 3 , the momentum errors of individual correction steps approached machine precision levels around ∆ = 10 −9 . We believe that the reason might be that the accuracy increases drastically during the final root finding step, such that the average root error is much smaller than the prescribed maximum. We conclude that cumulative effects of the corrections can likely be neglected. In case of evidence to the contrary, the solution would be to simply not recompute the momentum, sacrificing machine-precision consistency for error reduction.
We also applied the random walk model to states with different combinations of b = {0, 2}, v = {0, 0.99}, th = {0, 10}, perturbing the evolved variables separately with normally distributed relative errors of order 10 −4 . This test confirms that the implementation of the zero-temperature energy correction works as intended. We monitored the behavior of v i , W, , P, µ for the above cases and did not encounter any problematic behavior.

VI. SUMMARY AND DISCUSSION
In this work, we solved the technical problem of primitive variable recovery in relativistic ideal magnetohydrodynamic evolution codes via a new fully reliable scheme. We derived a mathematical proof that the algorithm always finds a valid solution, and that the solution is unique. Moreover, we derived expressions that allow to prescribe the accuracy of the individual primitive variables.
The guaranteed reliability of the new scheme is a big advantage compared to older methods, which are able to handle most of the parameter space encountered in BNS merger scenarios, but may still fail in some cases [33]. Even rare recovery failures are very problematic, since they necessitate manual intervention, and may require repeating parts of the simulation. Recovery failures are practically unpredictable and potentially chaotic (we recall the Newton-Raphson fractal related to convergence properties of a standard root finding procedure). This is aggravated for recovery schemes that rely on an initial guess taken from the previous timestep. The automated approach of using a fixed state (e.g., artificial atmosphere) in case of recovery failure will render simulations unpredictable in practice.
The ability to identify unphysical evolved variables, as well as the nature of the invalidity, is another advantage of our method. All evolution schemes produce unphysical states occasionally, most of which are harmless. However, sometimes invalid input occurs as first symptom of more severe evolution errors. Our method allows to prescribe an error policy and selectively apply corrections based on the nature of the problem. Such corrections (or lack of corrections) should be considered as part of the evolution scheme, but are mentioned in the literature only rarely.
The design of our scheme naturally suggests a particular prescription for correcting slightly unphysical input. We discussed potential cumulative effects of those corrections, and predicted that it will create artificial heating if the matter is close to zero temperature. We also showed that there should be no direct impact on the momentum. We validated this by performing a numerical experiment using random walk perturbations to emulate evolution errors.
Since the implementation of recovery algorithms is a workintensive endeavor, we are making our reference implementation public in form of a well-documented library named RePrimAnd [34], which can be used by any evolution code. In order to be useful in practice, the recovery should not constrain the type of EOS. Therefore, our recovery algorithm is formulated in an EOS-agnostic manner, and the reference implementation contains a generic interface for using arbitrary EOS.
We subjected the code to a comprehensive suite of tests, demonstrating that both the algorithm and the actual implementation are robust up to Lorentz factors and values of magnetization much larger than those relevant for BNS mergers, and beyond the magnetization tested in [33]. We also showed that the scheme is computationally efficient.
While investigating the accuracy of the recovery scheme, we identified regimes where rounding errors are amplified by unavoidable cancellation errors. We quantified the dominant contributions and found that the accuracy measured in our tests is compatible with the predictions. We also found that the rounding errors are irrelevant because the very same cancellation also leads to the amplification of evolution errors. Investigating the error propagation from evolved to primitive variables, we showed that the accuracy of the thermal energy and thermal pressure can severely degrade when evolving strongly magnetized regions of low density.
We believe that our results will be useful in particular for studying the launching mechanism of jets powering SGRBs, as well as the mass ejection processes that are ultimately responsible for kilonova signals. Both astrophysical scenarios involve strongly magnetized matter.