© The Authors 2025

1 Introduction

Over the last several decades, the scientific community has witnessed a remarkable improvement in numerical techniques, which has been particularly useful in the astrophysical hydrodynamics scenario for the possibility of modelling the Universe with unprecedented accuracy and of numerically exploring previously inaccessible systems. Between the huge variety of techniques exploited in the astrophysical context, the most widely used surely are those relying on smoothed particle hydrodynamics (SPH; Lucy 1977; Gingold & Monaghan 1977; Price 2008, 2012; Cullen & Dehnen 2010; Read & Hayfield 2012) or on grid-based finite-volume methods (e.g. Berger & Colella 1989; Stone & Norman 1992). Both these approaches have their advantages and disadvantages; for example, while SPH codes have continuous adaptive resolution due to their Lagrangian nature, and great conservation properties, and provide a framework easily coupled with N-body gravity, they suffer from extra computational cost in neighbour finding and some resolution-scale noise associated with particle rearrangement that could cause fluid instabilities (Lecoanet et al. 2016).

Grid-based codes are numerically stable and provide good shock capture; however, they suffer from bad mass and angular momentum conservation, grid alignment effects, large artificial diffusion related to advection errors (Gardiner & Stone 2005), and wasted computational time in empty regions, usually addressed with ‘adaptive mesh refinement’ (AMR) implementations (Berger & Colella 1989). Lately, a number of hybrid methods have tried to get the best of both approaches, as moving mesh schemes based on a Voronoi tessellation of the domain (Springel 2010; Pakmor et al. 2011), particle-in-cell (PIC) simulations (Aunai et al. 2024), and mesh-less methods (Lanson & Vila 2008a,b). The latter have been then implemented in astrophysical codes by (Gaburov & Nitadori 2011, Weighted Particle MHD), (Hopkins 2015; Hopkins & Raives 2016; Hopkins 2016, GIZMO), and (Valentini 2024, OpenGADGET). In particular, GIZMO has been extensively used in many astrophysical scenarios, ranging from the study of active galactic nucleus (AGN) accretion discs and circumbinary disc formation and evolution (Franchini et al. 2022, 2024; Garg et al. 2024) to tidal disruption events (Mainetti et al. 2017), star-forming regions (Grudić et al. 2021; Lupi et al. 2021), and cosmological galaxy formation simulations (Lupi et al. 2017; Hopkins et al. 2018, 2023), proving its accuracy and effectiveness.

All these inquiries were treated with Newtonian physics, yet only recently has a general relativistic hydrodynamics (GRHD) extension been introduced in GIZMO (Lupi 2023). However, an implementation of general relativistic magnetohydrodynamics (GRMHD) in GIZMO was missing. Magnetic fields are ubiquitous in astrophysical plasma and many works have studied the importance of highly magnetized regions in the proximity of compact objects. In fact, the magnetic field is a critical ingredient for the development of magnetorotational instabilities, potentially explaining the angular momentum transfer in AGN accretion discs (Balbus & Hawley 1991; Igumenshchev et al. 2003; Narayan et al. 2003). In addition, magnetically dominated funnels close to the event horizon of single or binary black holes (BHs) appear to be sources of relativistic jets of material, representing important features to look for in electromagnetic (EM) surveys (Palenzuela et al. 2010; Giacomazzo et al. 2012; Bogdanović et al. 2022; Cattorini & Giacomazzo 2024; Fedrigo et al. 2024). Tidal disruption events are other important scenarios in which magnetic fields play an important role; in fact, when a star’s equilibrium gets destroyed by a close-by compact object, the initially confined magnetic field eventually permeates the inflowing stream of gas, affecting its dynamics and contributing to potential signature EM emissions. Finally, a series of works presenting a recent GIZMO simulation (Hopkins et al. 2024c,a,d,b) has studied the formation and evolution of magnetically dominated AGN accretion discs originating in cosmological initial conditions and resolving a dynamical range spanning from ~100 Mpc down to ~100 au in a consistent run. The accretion disc’s innermost region remains, however, unresolved and has been studied using the code H-AMR in Kaaz et al. (2025).

These considerations motivate us to implement a GRMHD scheme in the GIZMO code. Although many other codes that evolve the GRMHD set of equations exist (Gammie et al. 2003; Del Zanna et al. 2007; Giacomazzo & Rezzolla 2007; Mösta et al. 2014; Etienne et al. 2015; Fragile et al. 2019; Kidder et al. 2017; Cipolletta et al. 2020; Liska et al. 2022), none of them do it in a Lagrangian framework. To date, our scheme is the first one to solve the GRMHD equations by exploiting the benefits of a mesh-less scheme. With our implementation, GIZMO can now self-consistently simulate the Universe from cosmological down to BH event horizon scales, including a plethora of physical prescriptions such as cosmological integration, dark matter, radiative heating and cooling (Hopkins 2017), stellar and BH formation and feedback (e.g. Grudić et al. 2021; Guszejnov et al. 2021; Hopkins et al. 2023; Bollati et al. 2024), chemistry, dust, and now also GRMHD effects.

However, a solid implementation of magnetic fields represents a serious challenge due to the need to keep the magnetic field numerical divergence as low as possible. In fact, a non-zero ∇·B can lead to wrong results or to code breaking; a relativistic scheme is even more sensitive to this issue due to the necessity of retrieving the primitive physical quantities from conserved ones via a numerical iterative method (the conservative-to-primitive algorithm). It is easy to understand that a non-zero magnetic field divergence would break the convergence of the latter.

Two methods are commonly employed to control the spurious ∇·B build-up: the constrained transport (CT) scheme (Evans & Hawley 1988), which is able to maintain ∇·B constant to machine precision, and the evolution of the vector potential (Etienne et al. 2012). However, these methods are not straightforward to implement in a mesh-less scheme (Tu et al. 2022). Therefore, we opt here for a hyperbolic divergence cleaning approach (Powell et al. 1999; Dedner et al. 2002), taking inspiration from the Newtonian implementation already present in the original version of GIZMO (Hopkins & Raives 2016). We show that our scheme can pass several special and general relativistic MHD stress tests, keeping the magnetic field divergence value under safe levels.

This paper is organized as follows: in Sect. 2 we recall the GRMHD set of equations, written according to the 3+1 Valencia formulation; in Sect. 3 we present our numerical implementation of the scheme, quickly reviewing the GIZMO mesh-less method, and then focusing on the MHD additions to the GRHD implementation (Lupi 2023); in Sect. 3.3 we explain in detail our divergence cleaning approach; then, in Sect. 4 we show the results from several validation tests, both in special and general relativity. Finally, we draw our conclusions in Sect. 5.

2 GRMHD equations

In this section, we describe the GRMHD set of equations implemented in our code. Our scheme follows the model described in the GRHydro (Mösta et al. 2014) code, implementing the Valencia formulation (Banyuls et al. 1997) of the equations. We also refer the reader to Font (2008) and Mizuno & Rezzolla (2025) for a review on numerical MHD and GRMHD and to Komissarov (1999), Koide et al. (1999), Koide et al. (2000), Koide et al. (2002), De Villiers & Hawley (2003), and Gammie et al. (2003) for pioneering special and general relativistic MHD implementations in multi-dimensional simulation codes. In what follows we use Greek indices for 4D quantities and Latin indices for 3D ones. We also employ natural units, where G = c = M⊙ = 1, and we make use of Einstein notation for repeated indices.

We assume a 3+1 spacetime decomposition in the Arnowitt-Deser-Misner (ADM) formalism (Arnowitt et al. 2008) as

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=(-{\alpha }^2+{\beta }_i{\beta }^i)d{t}^2+2{\beta }_{i}d{x}^{i}dt+{\gamma }_{ij}d{x}^{i}d{x}^j,$$(1)

where g_μν is the 4D metric, α is the lapse function, βⁱ is the shift vector, and γ_ij is the 3D spatial metric. We concentrate on ‘ideal GRMHD’, i.e. we assume infinite conductivity of the fluid and neglect viscosity and heat conduction. This corresponds to imposing u_νF^μν = 0, where u_ν is the four-velocity of the fluid and F^μν is the Faraday tensor. Under this assumption, the electric field as seen by a Eulerian observer can be written as a function of the magnetic field, Bⁱ, as $E^i=-\alpha {\epsilon }^{0ijk}{v}_j{B}_k$, while it vanishes in the fluid rest frame, as is expected for a perfect conductor (Banyuls et al. 1997; Antón et al. 2006; Font 2008).

To describe the fluid, we use the following quantities: ρ is the fluid rest mass density, ε the specific internal energy, P the gas pressure, and $v^i=u^i/W+{\beta }^i/\alpha$ is the fluid three-velocity. In the latter we identify $W=(1-v^i{v}_i{)}^{-1/2}$ as the Lorentz factor, where $v_i{v}^i={\gamma }_{ij}{v}^j{v}^i$. With the addition of the following definitions for the time and spatial components of the magnetic field in the fluid’s rest frame,

$$\begin{array}{rl}& b^0\equiv \frac{W{B}^j{v}_j}{\alpha },\\ & b^i\equiv \frac{B^i}{W}+W(B^j{v}_j)(v^i-\frac{{\beta }^i}{\alpha }),\\ & b^2\equiv \frac{B^j{B}_j}{W^2}+(B^j{v}_j{)}^2,\end{array}$$(2)

we can write the stress-energy tensor of ideal GRMHD as

$$\begin{array}{rl}{T}^{\mu \nu }& =\rho (1+ϵ+\frac{(P+b^2)}{\rho })u^{\mu }{u}^{\nu }+(P+\frac{b^2}{2})g^{\mu \nu }-b^{\mu }{b}^{\nu }\\ & \equiv \rho h^{\ast }{u}^{\mu }{u}^{\nu }+P^{\ast }{g}^{\mu \nu }-b^{\mu }{b}^{\nu },\end{array}$$(3)

where in the second line we have defined the magnetically modified enthalpy, $h^{\ast }\equiv 1+ϵ+(P+b^2)/\rho$, and the total pressure, $P^{\ast }\equiv P+b^2/2$.

The equations of GRMHD are then derived from the mass density conservation law

$$\frac{1}{\sqrt{-g}}{\mathrm{\partial }}_{\mu }(\sqrt{-g}\rho u^{\mu })=0,$$(4)

where g =det(g_μν) is the determinant of the four-metric, the conservation of energy-momentum

$${\mathrm{\nabla }}_{\mu }{T}^{\mu \nu }=0,$$(5)

and the Maxwell equations in the ideal regime

$${\mathrm{\nabla }}_{\nu }{}^{\ast }{F}^{\mu \nu }=0,$$(6)

where ${}^{\ast }{F}^{\mu \nu }=\frac{1}{2}{\epsilon }^{\mu \nu \lambda \sigma }{F}_{\lambda \sigma }$ is the dual of the Faraday tensor. This set of equations can be cast in hyperbolic form:

$$\frac{\mathrm{\partial }\mathbf{\text{U}}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathbf{\text{F}}}^i}{\mathrm{\partial }{x}^i}=\mathbf{\text{S}};$$(7)

here, U is the conserved quantities vector, defined as

$$\mathbf{\text{U}}=\{\begin{array}{ll}D& \equiv \sqrt{\gamma }\rho W,\\ S_j& \equiv \sqrt{\gamma }(\rho h^{\ast }{W}^2{v}_j-\alpha b^0{b}_j),\\ \tau & \equiv \sqrt{\gamma }(\rho h^{\ast }{W}^2-P^{\ast }-(\alpha b^0{)}^2)-D,\\ {\mathcal{B}}^j& \equiv \sqrt{\gamma }{B}^j,\end{array}$$(8)

where γ =det(γ_ij) is the determinant of the spatial metric, F is the flux vector,

$${\mathbf{\text{F}}}^i=\{\begin{array}{ll}& D{\tilde{v}}^i,\\ & S_j{\tilde{v}}^i+\alpha \sqrt{\gamma }{P}^{\ast }{\delta }_j^i-\alpha b_j{\mathcal{B}}^i/W,\\ & \tau {\tilde{v}}^i+\alpha \sqrt{\gamma }{P}^{\ast }{v}^i-{\alpha }^2{b}^0{\mathcal{B}}^i/W,\\ & {\mathcal{B}}^j{\tilde{v}}^i-{\mathcal{B}}^i{\tilde{v}}^j,\end{array}$$(9)

where ${\tilde{v}}^i=\alpha v^i-{\beta }^i$, and S is the sources vector (Duez et al. 2005; Chang & Etienne 2020),

$$\mathbf{\text{S}}=\alpha \sqrt{\gamma }\{\begin{array}{ll}& 0,\\ & \frac{1}{2}{T}^{\alpha \beta }\mathrm{\nabla }{g}_{\alpha \beta },\\ & (T^{00}{\beta }^i{\beta }^j+2T^{0i}{\beta }^j+T^{ij})K_{ij}-(T^{00}{\beta }^i+T^{0i}){\mathrm{\partial }}_i\alpha ,\\ & 0,\end{array}$$(10)

where K_ij is the extrinsic curvature of the metric. We note that the source terms depend on the derivative of the metric and we do not have any modification in form from the GRHD case (note, however, that the stress-energy also considers the magnetic energy-momentum contribution now).

Finally, the time component of Eqs. (6) yields the divergence-free condition of the magnetic field,

$$\frac{1}{\sqrt{\gamma }}{\mathrm{\partial }}_i{\mathcal{B}}^i=0;$$(11)

numerically, this constraint needs to be ensured during the system evolution. We enforce this condition with a ‘divergence cleaning’ technique, which we describe in Sec. 3.3.

An equation of state (EoS) is required to close the system of equations (7). In our work, we always use a Γ-law EoS of the form P = (Γ - 1)ρε, albeit our scheme also allows for different equations of state.

3 Numerical implementation

In this section, we describe the numerical methods used to solve the GRMHD set of Eqs. (7), mainly focusing on the differences between the GRHD version of the scheme, presented in Lupi (2023), and the current implementation. After a brief recall of the GIZMO mesh-free methods, we describe the new volume evolution prescription, the GRMHD Riemann solver, and, finally, the methods used to ensure the divergence-free condition.

3.1 Review of the GIZMO mesh-less methods

The code GIZMO solves the evolution equations in a standard finite-volume Godunov-type style (Hopkins 2015; Hopkins & Raives 2016). The unique trait of the code resides in the domain discretization: the volume is partitioned among elements (which we call either ‘particles’ or ‘cells’ according to need) according to a continuous weighting function, Ψ_i, so that at each meshgenerating point position, x, an infinitesimal volume, dx, can be described by

$$d\mathbf{\text{x}}=\sum _i{\mathrm{\Psi }}_i(\mathbf{\text{x}},h(\mathbf{\text{x}})),$$(12)

where h(x) is the kernel size. The latter is evaluated from the code for each particle to encompass a user-defined number of neighbours. The weighting function, Ψ_i, is defined as

$${\mathrm{\Psi }}_i(\mathbf{\text{x}},h(\mathbf{\text{x}}))=\frac{W({\mathbf{\text{x}}}_i-\mathbf{\text{x}},h(\mathbf{\text{x}}))}{\sum _{j}W({\mathbf{\text{x}}}_j-\mathbf{\text{x}},h(\mathbf{\text{x}}))},$$(13)

where W is any continuous, symmetric kernel function with compact support (for instance, we usually employed a cubic spline kernel function, unless otherwise indicated). The sum at the denominator normalizes the function, Ψ_i, so that the volumes always sum correctly.

As in the GRHD case (Lupi 2023), we decided to solve Eqs. (7) in the lab reference frame, as

$$\frac{\mathrm{d}\mathbf{\text{U}}}{\mathrm{d}t}+\mathrm{\nabla }\cdot (\mathbf{\text{F}}-\mathbf{\text{w}}\otimes \mathbf{\text{U}})=\mathbf{\text{S}},$$(14)

where w is the frame velocity. As is described in the original methods paper (Hopkins 2015), we can discretize Eq. (14) over the discrete elements defined above, by integrating over the domain volume, obtaining

$$\frac{\mathrm{d}(V_i{\mathbf{\text{U}}}_i)}{\mathrm{d}t}+\sum _j{\mathbf{\text{F̃}}}_{ij}\cdot {\mathbf{\text{A}}}_{ij}=V_i{\mathbf{\text{S}}}_i,$$(15)

where F̃_ij is the solution of the Riemann problem between the ith particle and the j-th interacting neighbour, $V_i=\int {\mathrm{\Psi }}_i(\mathbf{\text{x}})d\mathbf{\text{x}}$ is the particle volume, and A_ij is the normal vector to the ‘effective face area’ between i-th and j-th interacting cells.

In principle, the frame velocity, w, in Eq. (14) can be chosen arbitrarily. GIZMO defines two different methods that depend on the assumption made for the frame velocities, as described in the original paper. One assumes the interacting face to be moving with the average fluid velocity of the interacting cells, $w_{ij}=({\tilde{v}}_i+{\tilde{v}}_j)/2$, defining the mesh-less finite volume (MFV) method, which resembles a standard moving-mesh finite volume scheme. The other assumes instead the frame velocity to be that of the contact discontinuity, i.e. the velocity which preserves the mass on both sides of the moving interface. This second choice defines the mesh-less finite mass (MFM) method, in which the particle masses are always preserved throughout the evolution of the system (more on that in Sect. 3.2.1).

3.2 MHD modifications

In this section, we describe the main method modifications from the GRHD version of the scheme, extensively described in Lupi (2023). The first difference resides in the implementation of the Riemann solver, which we describe in the following subsection.

3.2.1 Riemann solver

As in the GRHD case, we implement the 1D HLL Riemann solver (Harten et al. 1983). We find that our single intermediate state HLL solver provides a computationally cheap and sufficiently accurate, although more diffusive, solution to the problem; in the future, we shall also include a multi-state Riemann solver to the scheme, as in the scheme of Mignone & Bodo (2006) or Miyoshi & Kusano (2005), in order to correctly resolve contact discontinuities.

The fluxes given by the solver can be written as

$$\mathbf{\text{F̃}}=\{\begin{array}{l}{\mathbf{\text{F}}}_L-w^{\hat{n}}{\mathcal{U}}_L; w^{\hat{n}}<{\lambda }_{min}\\ {\mathbf{\text{F}}}_{HLL}-w^{\hat{n}}{\mathcal{U}}_{HLL}; {\lambda }_{min}\le w^{\hat{n}}\le {\lambda }_{max}\\ {\mathbf{\text{F}}}_R-w^{\hat{n}}{\mathcal{U}}_R; w^{\hat{n}}>{\lambda }_{max}\end{array},$$(16)

where U_HLL are the intermediate states evaluated as

$${\mathcal{U}}_{HLL}=\frac{{\lambda }_{max}{\mathcal{U}}_R-{\lambda }_{min}{\mathcal{U}}_L+{\mathbf{\text{F}}}_L-{\mathbf{\text{F}}}_R}{{\lambda }_{max}-{\lambda }_{min}}$$(17)

and FHLL are the corresponding fluxes,

$${\mathbf{\text{F}}}_{HLL}=\frac{{\lambda }_{max}{\mathbf{\text{F}}}_L-{\lambda }_{min}{\mathbf{\text{F}}}_R+{\lambda }_{max}{\lambda }_{min}({\mathcal{U}}_R-{\mathcal{U}}_L)}{{\lambda }_{max}-{\lambda }_{min}}.$$(18)

In these equations, $w^{\hat{n}}$ is the face velocity perpendicular to the interface, whereas λ_min and λ_max are the slowest and the fastest wave speeds of the Riemann problem, i.e. the magnetosonic speeds. We evaluate such velocities by solving the quadratic dispersion relation,

$$\begin{array}{r}{\lambda }^2[W^2(V^2-1)-V^2]-2\lambda [W^2{\tilde{v}}^{\hat{n}}(V^2-1)+V^2{\beta }^{\hat{n}}]\\ +[(W{\tilde{v}}^{\hat{n}}{)}^2(V^2-1)+V^2({\alpha }^2{\gamma }^{\hat{n}\hat{n}}-{\beta }^{\hat{n}}{\beta }^{\hat{n}})]=0,\end{array}$$(19)

where $V^2\equiv v_A^2+c_s^2(1-v_A^2)$, with $v_A=\sqrt{b^2/(\rho h^{\ast })}$ the Alfvén velocity and c_s the fluid sound speed.

In order to guarantee the stability of the scheme, the Courant-Friedrichs-Levy (CFL) condition must be fulfilled during the calculation. Even though the solution of λ should be used in this case, for computational efficiency we employ V as the maximum signal velocity. This represents an approximate evaluation of the signal speed, but gives a good enough estimate for our purpose.

To conclude this subsection, we recall that when GIZMO is running in MFM mode, the face velocity, $w^{\hat{n}}$, is assumed to be moving with the zero-mass flux velocity, which for our HLL solver is defined as

$$w^{\hat{n}}=\frac{{\mathbf{\text{F}}}_{D,HLL}}{{\mathcal{D}}_{HLL}}.$$(20)

3.2.2 Time integration and volumes evolution

After the flux computation, we need to evolve the conserved variables in time according to the generalized leap-frog scheme, as is described in Lupi (2023). The only difference with the GRHD scheme lies in the explicit update of the conserved variables through Eq. (15), both for the ‘kick’ procedure on active particles over half a timestep, Δt/2, and the prediction operation over a full timestep, Δt. In fact, in the GRHD scheme, all the conserved quantities (D, S j, τ)_i are proportional to the particle density, D_i, i.e. proportional to the particle mass, m_i, once volume-integrated. This results in just a normalization factor in the conservative-to-primitive inversion algorithm, irrelevant to its solution (Sec. 3.2.3). However, in the volume-integrated GRMHD conserved variables, a proportionality to V_i² appears in the terms containing B_i² in the momentum and energy density definitions, breaking the inversion scheme. For this reason, we modified our scheme, introducing an explicit evolution of the volume, V_i, as

$${\mathbf{\text{U}}}_i^{t+\mathrm{\Delta }t}=\frac{V_i^t[{\mathbf{\text{U}}}_i^t+\mathbf{\text{S}}({\mathbf{\text{U}}}_i^t)\mathrm{\Delta }t]-\sum _j{\mathbf{\text{A}}}_{ij}^t{\mathbf{\text{F̃}}}_{ij}^t\mathrm{\Delta }t}{V_i^{t+\mathrm{\Delta }t}}.$$(21)

This guarantees that the local conservative quantities are consistently passed to the conservative-to-primitive solver.

This procedure boils down to how to explicitly compute the updated volume, $V_i^{t+\mathrm{\Delta }t}$. We decided to estimate it by exploiting the continuity equation as $V_i^{t+\mathrm{\Delta }t}=V_i^t\exp ({\dot{V}}_i\mathrm{\Delta }t/V_i^t)$, where ${\dot{V}}_i=\sum _j{\mathbf{\text{A}}}_{ij}\cdot ({\mathbf{v}}_{part,i}-{\mathbf{w}}_{ij})$; here v_part,i is the i-th particle velocity (which is equal to the fluid velocity in the MFM scheme), w_ij is the velocity of the Riemann problem interface¹, and A_ij is the effective face area between particles i and j², For consistency, the same prescription is employed in the density prediction step, i.e. $D^{t+\mathrm{\Delta }t}=D^t\exp (-\dot{V}\mathrm{\Delta }t/V^t)$.

3.2.3 Other code modifications

In this section, we discuss other code modifications implemented in our GRMHD extension of GIZMO. Despite the use of slope limiters, the process of reconstructing the primitive velocity at the position of the face might yield a resulting fluid velocity slightly larger than unity. In the GRHD version of the code, this was prevented by means of a user-defined limiter on the maximum Lorentz factor allowed by the code. Here instead, we have found it more convenient to reconstruct the quantity Wvⁱ at the face position, as this quantity has physical values ranging in (−∞, ∞), naturally avoiding causality violation. After reconstruction, the corresponding Eulerian, vⁱ, is then easily recovered (Noble et al. 2006) and passed to the Riemann solver. Note that the implementation of a limiter in the evaluation of the Lorentz factor is still available in the code, but it is now optional, and not employed in any of the tests presented in this work.

Because of the presence of the Lorentz factor in the definition of the conservatives variables, primitive quantities cannot be recovered analytically, but only through iterative numerical techniques for the inversion of Eq. (8). Here, we implemented the 2D inversion scheme by Noble et al. (2006), suitably adapted to our scheme, which recovers the GRMHD primitive variables via a Newton-Raphson method. This scheme considers the presence of magnetic fields, can be easily extended to generic EoS, and provides a stable and computationally efficient solution to this inversion problem. We note that keeping the divergence of the magnetic field low is crucial for the correct behaviour of the conservative-to-primitive solver.

In some pathological cases, when the magnetic field contributes to a significant fraction of the total energy of the fluid, the recovery of the specific internal energy, , from the conserved energy, τ, can become a non-trivial numerical task, yielding significant errors on the evolution of the primitive quantities. In our tests, we encountered this problem only in the magnetized Bondi accretion test, when dealing with a strong gravitational field (Sec. 4.2.2). For this reason, we implemented an optional energy-entropy switch, as was suggested in Del Zanna et al. (2007). When active, we evaluate the entropy function $s\equiv P/{\rho }^{\mathrm{\Gamma }}$, which we evolve through the conservation equation

$${\mathrm{\nabla }}_{\mu }(\rho s{u}^{\mu })=0,$$(22)

which can be rewritten as

$${\mathrm{\partial }}_t(Ds)+{\mathrm{\partial }}_i(Ds{\tilde{v}}^i)=0.$$(23)

We note that the last equation only holds for adiabatic flows, in the absence of shocks or any other energy dissipation mechanisms. During the fluid evolution, for each particle the code checks whether the magnetic field contributes to a significant fraction of the total energy and if any shock is present in its interaction with neighbouring particles³. If the first condition is met and shocks are absent, the internal energy is directly recovered from the entropy, 5. Otherwise, if the first condition is false or a shock is present in any direction, we evaluate the internal energy relying on the standard evolution equations (7).

3.3 Divergence cleaning

The ideal GRMHD equations (7) were derived using Eq. (11), i.e. the no-monopole constrain. It is well established, however, that such a constraint is not always guaranteed in numerical schemes. Indeed, the numerical approximations performed during the calculations can easily result in the divergence of the magnetic field significantly differing from zero, leading to an incorrect and unphysical evolution of the conserved quantities or, worse, to code breaking. Keeping the emerging errors in Eq. (11) under control is therefore crucial for a correct evolution of the systems considered. This is usually achieved with numerical techniques such as the constrained transport (CT; Evans & Hawley 1988) or vector potential evolution (Etienne et al. 2012) methods. While successful implementations of these methods in moving-mesh relativistic MHD schemes exist (e.g. He & Tang 2012; Fragile et al. 2019), we opt here for a divergence cleaning technique, to maintain consistency with the Newtonian MHD version of the code (Hopkins & Raives 2016). To date, this is the first implementation of a divergence cleaning method in a GRMHD set of equations solved on a mesh-less scheme.

We now describe our divergence cleaning implementations, starting from the Powell’s ‘eight-wave’ cleaning prescription. As was described in the original paper (Powell et al. 1999), a set of source terms is needed to stabilize the MHD equations against the development of spurious magnetic divergence when the fluxes are evaluated. A special relativistic version of this scheme is analyzed in depth by Wu & Shu (2020). Here, we write a general relativistic analogue, similar to the one presented in Liebling et al. (2010).

We extended the set of equations (7) with additional source terms in the form

$${\mathbf{\text{S}}}_{Powell}=-\alpha ({\mathrm{\partial }}_i{\mathcal{B}}^i)\{\begin{array}{l}0,\\ \frac{B_j}{W^2}+v_j{B}^k{v}_k,\\ B^k{v}_k,\\ v^j-\frac{{\beta }^j}{\alpha }.\end{array}$$(24)

The critical point here resides in the need to consistently evaluate the ${\mathrm{\partial }}_i{\mathcal{B}}^i$ term. As was also reported in the original MHD GIZMO paper (and suggested by Gaburov & Nitadori 2011), the only valuable choice is to compute the divergence term directly from the Riemann solver outputs as

$$(V{\mathrm{\partial }}_k{\mathcal{B}}^k{)}_i=-\sum _j{\overline{\mathcal{B}}}_{ij}^{\perp }|{\mathbf{\text{A}}}_{ij}|,$$(25)

where

$${\overline{\mathcal{B}}}_{ij}^{\perp }=\frac{{\mathcal{B}}_i^{\perp }+{\mathcal{B}}_j^{\perp }}{2}$$(26)

is the average magnetic field perpendicular to the face between each i- j interacting pair of particles. This prescription alone can complete all of our tests without the code breaking. However, some test problems (e.g. the magnetic rotor, Sect. 4.1.4) yield wrong solutions. In order to address these inconsistencies and reduce the errors in the solutions, we added a hyperbolic divergence cleaning on top of the Powell’s prescription. Presented for the first time in Dedner et al. (2002) in its Newtonian formulation and later extended to special and general relativity (Neilsen et al. 2006; Anderson et al. 2006; Liebling et al. 2010; Penner 2011; Mösta et al. 2014; Fambri et al. 2018), the method consists of adding a scalar field, ψ, used to advect divergence errors away from their source location as fast as possible while damping them.

To consistently derive the equations governing the evolution of the ψ field, we started from a modification of Maxwell’s equation (following Mösta et al. 2014),

$${\mathrm{\nabla }}_{\mu }({}^{\ast }{F}^{\mu \nu }+g^{\mu \nu }\hat{\psi })=\sigma n^{\nu }{c}_h\hat{\psi },$$(27)

where $\hat{\psi }=\psi /c_h,n^{\nu }=(1/\alpha ,-{\beta }^j/\alpha )$ is the normal vector to the hypersurface of the 3+1 formulation of spacetime and where we have introduced the advection velocity, ch , and the damping factor, σ. As is suggested in Tricco et al. (2016), we follow the ψ quantity to take into consideration the variability of the advection velocity, c_h, obtaining a better evolution of the energy of the system. We can now solve Eq. (27) for ν = j (spatial components) and ν = 0 (time component), obtaining

$$\{\begin{array}{l}{\mathrm{\partial }}_t{\mathcal{B}}^j+{\mathrm{\partial }}_i({\tilde{v}}^i{\mathcal{B}}^j-{\tilde{v}}^j{\mathcal{B}}^i)=-\alpha \sqrt{\gamma }{\gamma }^{ij}{\mathrm{\partial }}_i\hat{\psi }+{\beta }^j{\mathrm{\partial }}_i{\mathcal{B}}^i,\\ {\mathrm{\partial }}_t\hat{\psi }-{\beta }^i{\mathrm{\partial }}_i\hat{\psi }=-\alpha {\gamma }^{-1/2}{\mathrm{\partial }}_i{\mathcal{B}}^i-\alpha \sigma c_h\hat{\psi }.\end{array}$$(28)

As was suggested in Dedner et al. (2002) and implemented in the MHD version of GIZMO, we considered an additional advection term, $\alpha v^i{\mathrm{\partial }}_i\hat{\psi }$, in the scalar field, ψ, evolution equation, which gives

$$\{\begin{array}{l}{\mathrm{\partial }}_t{\mathcal{B}}^j+{\mathrm{\partial }}_i({\tilde{v}}^i{\mathcal{B}}^j-{\tilde{v}}^j{\mathcal{B}}^i)=-\sqrt{\gamma }({\beta }^j{v}^i+\alpha {\gamma }^{ij}){\mathrm{\partial }}_i\hat{\psi }+{\beta }^j{\mathrm{\partial }}_i{\mathcal{B}}^i,\\ {\mathrm{\partial }}_t\hat{\psi }+{\tilde{v}}^i{\mathrm{\partial }}_i\hat{\psi }=-\alpha {\gamma }^{-1/2}{\mathrm{\partial }}_i{\mathcal{B}}^i-\alpha \sigma c_h\hat{\psi }.\end{array}$$(29)

If we now multiply the second equation by the conserved density, D, and we rearrange the terms using the mass continuity Eq. (7), we can rewrite the evolution equation for ψ̂ as

$${\mathrm{\partial }}_t(D\hat{\psi })+{\mathrm{\partial }}_i({\tilde{v}}^{i}D\hat{\psi })=-\frac{\alpha D{\mathrm{\partial }}_i{\mathcal{B}}^i}{\sqrt{\gamma }}-\alpha \sigma c_{h}D\hat{\psi }.$$(30)

Once averaged over the volume V_i, Eq. (30) yields the evolution of the mass-based ψ̂ field as implemented in our code, which translates into a much easier evaluation of the flux term.

As was done for the ${\mathrm{\partial }}_i{\mathcal{B}}^i$ term in Eq. (26), we need to evaluate the gradient of ψ̂ in a consistent way. In our implementation, this is done by considering the Rusanov fluxes between each pair of interacting particles,

$${\mathbf{\text{F}}}_{HLL}=\frac{1}{2}({\mathbf{\text{F}}}_L+{\mathbf{\text{F}}}_R)-\frac{\tilde{\lambda }}{2}({\mathcal{U}}_R-{\mathcal{U}}_L),$$(31)

for the decoupled (${\mathcal{B}}^{\perp },\hat{\psi }$) system, where λ̃ is the eigenvalue of the 2D problem. Here, $\tilde{\lambda }=\alpha \sqrt{{\gamma }^{\perp \perp }}$, where ${\gamma }^{\perp \perp }$ is the spatial metric component perpendicular to the interface, and from Eq. (31) we obtain

$$\{\begin{array}{l}{\overline{\mathcal{B}}}_{ij}^{\perp }=\frac{1}{2}({\mathcal{B}}_L^{\perp }+{\mathcal{B}}_R^{\perp })-\frac{1}{2}\sqrt{\gamma {\gamma }^{\perp \perp }}({\hat{\psi }}_R-{\hat{\psi }}_L),\\ {\overline{\psi }}_{ij}=\frac{1}{2}({\hat{\psi }}_L+{\hat{\psi }}_R)-\frac{1}{2}(\gamma {\gamma }^{\perp \perp }{)}^{-1/2}({\mathcal{B}}_R^{\perp }-{\mathcal{B}}_L^{\perp }).\end{array}$$(32)

Finally, ${\overline{\mathcal{B}}}_{ij}^{\perp }$ and ${\overline{\psi }}_{ij}$ are used to evaluate the divergence of ${\mathcal{B}}^j$ with Eq. (25) and the gradient of ψ̂ as

$$(V{\mathrm{\partial }}_k\hat{\psi }{)}_i=-\sum _j{\overline{\psi }}_{ij}{\mathbf{\text{A}}}_{ij}.$$(33)

In principle, we are free to choose the Dedner parameters, c_h and σ, arbitrarily. After extensive tests of the scheme, we have opted for defining σ = K/Δx, with K = 0.1 (except for the blast wave tests, in which we use K = 0.75), where Δx is the effective cell size. In a relativistic scheme, the advection velocity is usually set equal to the speed of light, i.e. c_h = 1, since this is the largest value permitted that does not violate causality. This choice, however, would break the hierarchical time-stepping in GIZMO, requiring the entire system to advance at the smallest possible pace. In our implementation, instead, we decided to employ the local magnetosonic speed, v_ms, similar to what is done in the Newtonian MHD version of the code, suitably augmented by a user-defined factor, f , i.e. $c_h=v_{ms}(1+f)/(1+f{v}_{ms}^2)$⁴. In our tests, we typically assumed f = 1 (except for the magnetized TOV equilibrium test where we use f = 2.5), which represents a compromise between computational efficiency and an effective divergence cleaning⁵. We note that, while c_h in Eq. (32) is evaluated from the local quantities entering the Riemann problem for the considered pair of particles, the maximum among all c_h values from the interacting particles pair is employed in Eqs. (29)-(30). We also note that when hyperbolic divergence cleaning is employed, the additional advection speed must be considered when evaluating the system time-step through the CFL condition.

Finally, in order to stabilize the scheme, we need to include two additional source terms in the evolution of the momentum density, S_j, and energy density, τ, yielding a Dedner source vector of the form

$$S_{Dedner}=\{\begin{array}{l}0,\\ {\mathrm{\Xi }}^k[2B_k{v}_j-B_j{v}_k-{\gamma }_{kj}(B^i{v}_i)],\\ {\mathrm{\Xi }}^k[2B_k(1-\frac{1}{2W})-v_k(B^i{v}_i)],\\ {\mathrm{\Xi }}^j,\end{array}$$(34)

where

$${\mathrm{\Xi }}^k\equiv [-\sqrt{\gamma }({\beta }^k{v}^i+\alpha {\gamma }^{ki}){\mathrm{\partial }}_i\hat{\psi }+{\beta }^k{\mathrm{\partial }}_i{\mathcal{B}}^i].$$(35)

This concludes the description of our divergence cleaning procedure.

4 Code validation tests

We next tested the robustness and accuracy of our implementation against some standard tests, both in special and general relativistic MHD. Unless otherwise indicated, all the tests were run with the complete divergence cleaning scheme (with parameters K = 0.1 and f = 1).

4.1 Special relativistic magnetohydrodynamics

4.1.1 Monopole

To test the behaviour and the effectiveness of our divergence cleaning prescription, we performed a simple test in which we initialized a magnetic monopole in a static fluid. The initial setup was similar to the one presented in Mösta et al. (2014): a distribution of 100³ particles initially placed on a Cartesian mesh filled a periodic cubic box of side length 4. The fluid was uniform in density and pressure and it was initially at rest. The magnetic field was set equal to zero everywhere except in a central region where a monopole was introduced in the form

$$B^x=\{\begin{array}{ll}{e}^{-r^2/R_c}-e^{-1}; & r<R_c\\ 0;& r\ge R_c\end{array}$$(36)

where R_c = 0.2. The adiabatic index, Γ, was set to 5/3 during the evolution, we assumed a Minkowski flat background, and we chose the Dedner parameter to be c_h = 1 and either K = (0.1, 1). We report the results of the test performed in MFM mode but we find no difference in the MFV case.

In Fig. 1 we show a 2D slice of the magnetic field divergence evaluated according to Eq. (25) at four times, t = 0, 0.5, 1, 1.5, for both the K = 0.1 and the K = 1 cases. Our scheme can correctly remove spurious monopoles, advecting and damping them in a relatively short time. As was expected, the K = 1 run displays a stronger reduction of ${\mathrm{\partial }}_i{\mathcal{B}}^i$; however, a slightly slower suppression of ψ̂ is recommended so that the evolution of the magnetic field itself is better affected by the corrective terms in Eq. (29). This analysis and the following tests show that higher K values result in inefficient cleaning, while too low K values lead to system instabilities. We also note that even in this simple test, the corrective terms in Eq. (34) are necessary to keep the other evolved quantities stable. For completeness, in Appendix B we show results of the same test performed with the Powell scheme only.

4.1.2 MHD shock tubes

We report here the results of two 1D special relativistic MHD tests. Presented for the first time in Balsara (2001), these tests consist of two-state (left and right) Riemann problems, which are particularly useful for demonstrating the ability of the code to capture shocks. In particular, we performed the relativistic analogue of the Brio-Wu shock tube problem (Brio & Wu 1988), usually called a ‘Balsara1’ test, and a relativistic MHD collision, called a ‘Balsara4’ test.

The initial conditions of the Balsara1 problem consist of a left state (ρ, v^x, v^y, v^z, ε, B^x, B^y, B^z)_L = (1.0, 0, 0, 0, 1.0, 0.5, 1.0, 0) and a right state (ρ, v^x, v^y, v^z, ε, B^x, B^y, B^z)_R = (0.125, 0, 0, 0, 0.8, 0.5, -1.0, 0), with Γ = 2. The discontinuity exists in the x direction. The domain is a periodic 2D box of length 2 in the x direction and 0.05 in the y direction.

To perform the test in the MFM mode, we initialized a distribution of equal-mass particles placed on a Cartesian grid, so that 2297 particles were sampling the x direction. To determine the two different states’ rest-mass densities, the left state had a particle number density eight times higher than the right one. Because of periodicity, two discontinuities are present: one at the centre of the box and one at the edge. We focus here only on the central one, as the other is simply a mirrored version of it, and hence we limit our ‘active’ region to the central half of the box. This translates into a number of active particles in the x direction of ~ 1148. For the MFV case, we discretized the fluid with a closely packed lattice with 2000 particles in the x direction (1000 ‘active’ particles).

In Fig. 2, we plot our solution of the Balsara1 test at t = 0.4 performed in both MFM and MFV modes, using Powell terms only and with the complete cleaning scheme (Powell+Dedner). In the latter, we employ K = 0.1 and f = 1. The exact solution reported in the panels was computed with the exact RMHD Riemann Solver by Giacomazzo & Rezzolla (2006). Our scheme is able to reproduce the correct solution with good accuracy. We recognize some small oscillations on the left of the rarefaction fan and on the right of the fast magnetosonic wave front in all our runs; these artefacts can be controlled with more diffusive slope limiters when performing the reconstruction of the primitive quantities on the interacting faces. In this test, the specific internal energy, , is slightly overestimated at the contact discontinuity; close to the discontinuity, this effect is more severe in the MFM case, but exhibits an overall higher value in the MFV one. At the shock front (at x ' 1.15) oscillations develop in the velocity and the specific internal energy evolution; this effect is probably related to the high volume compression the particles witness when the shock front is approaching. We can attenuate this artefact with more stringent slope limiters, lowerorder reconstructions of the primitive quantities, by employing bigger particle kernels (e.g. the Wendland C4 kernel function) or going to a higher resolution, but some oscillations are always present when particle motion is allowed. Furthermore, initializing a particle configuration without sharp mass gradients relieves these difficulties.

In the bottom panel of Fig. 2, we display the specific internal energy when the problem is initialized with equal-mass particles and evolved with the MFV scheme. When compared with the same problem initialized with equal-volume particles (orange circles in the second panel of Fig. 2), the quantity exhibits a better agreement with the exact solution in the region between the contact wave and the rightgoing slow shock, resembling the MFM solution. Finally, in the same panel, we show results obtained using more stringent slope limiters; while being more diffusive, the solution has significantly fewer numerical artefacts; in fact, both the oscillations on the fast wave and the right slow wave and the overestimation of in front of the contact discontinuity are attenuated.

In a planar shock test, the B^x evolution is a direct indication of the ability of the code to maintain a low magnetic field divergence, as it should remain constant in space and time. We recognize a good evolution of this quantity in all our runs, with only a moderate spike at the contact discontinuity in the MFM mode; however, with the addition of the Dedner divergence cleaning, the B^x evolution is improved overall. Finally, in the MFM runs, the magnetic field witnesses a slow noise build-up in regions yet unaffected by the shock; this is due to the volume discretization of the domain being performed on a 2D Cartesian grid in our initial conditions, which is not optimal for our ${\mathrm{\partial }}_i{\mathcal{B}}^i$ and flux evaluations. Indeed, a better mapping of the simulated volume via a closed packed lattice, as was done in the MFV case, completely removes this artefact⁶. As a check of our implementation, we also performed this test with a fixed particle distribution (similar to a classical finite-volume scheme on a grid). The results, despite not being shown, are in perfect agreement with the exact solution, without any of the numerical artefacts mentioned above, suggesting that their origin has to be ascribed to the particle motion, which reduces numerical diffusivity at a relatively low resolution, but at the same time introduces some spurious oscillations (Hopkins 2015).

For the Balsara4 problem, we initialized a left state (ρ,ν^x,v^y,v^z, ε, B^x, B^y, B^z)_L = (1.0,0.999,0,0,0.15,10,7.0,7.0) and a right state (ρ, v^x, v^y, v^z, ε, B^x, B^y, B^z)_R = (1.0, -0.999, 0, 0, 0.15, 10, -7.0, -7.0) with Γ = 5/3. Since the initial Lorentz factor is W = 22.366, this test is highly relativistic. We filled a 2D (3, 0.05) periodic box with a closely packed lattice so that there were 500 ‘active’ particles sampling the x direction.

In Fig. 3, we plot our solution at time t = 0.4 performed in MFM mode, using either the Powell terms only or with the complete cleaning scheme, alongside the exact solution of the problem. Quantities v_z and B_z are not plotted but exhibit analogue behaviours to v_y and B_y. Our code is able to correctly compute the evolution of the quantities, with only a small overestimation of the specific internal energy and an underestimation of the rest mass density (yielding a correct evaluation of the pressure) in the central regions. The positions of the shock fronts are slightly off, becoming more accurate with the Dedner cleaning scheme. The outer shocks remain instead slightly offset, because of the poor sampling around the strong density jump, which is reflected in a smoothing of the discontinuity across a few neighbours. To verify this hypothesis, we performed a high-resolution run with ~1166 active particles sampling the x direction. The solution yields sharper discontinuities and a better B^x conservation, but still not sufficiently to perfectly resolve the large density jump, which requires an even higher resolution.

Fig. 1 Evolution of the monopole test. 2D slice at z = 0 of the magnetic field divergence at the initial time (left) and at three later times. The top panels display the solution with damping parameter K = 0.1 and the bottom panels the one with K = 1.

Fig. 2 Slice at y = y0 of primitive quantities at time t = 0.4 of the Balsara1 test. The magnetic field divergence in units of the magnetic field intensity is plotted in the seventh panel. Blue squares (Powell) and red triangles (Dedner) mark the MFM solutions, while cyan stars (Powell) and orange circles (Dedner) indicate the MFV ones. In the bottom panel, we display the specific internal energy when the problem is initialized with equal-mass particles and evolved with the MFV scheme; gold diamonds mark the solution employing standard slope limiters, while olive triangles show the results using more diffusive slope limiters. The correct solution, evaluated through the exact Riemann Solver by Giacomazzo & Rezzolla (2006), is plotted with a dashed black line.

Fig. 3 Slice at time t = 0.4 and constant y of the primitive quantities from the Balsara4 test performed in the MFM mode. The magnetic field divergence in units of its intensity is plotted in the bottom panel. Blue squares mark the solution computed with Powell terms only, while red triangles indicate the one computed with the complete divergence cleaning scheme. The high-resolution run is displayed with orange circles. Quantities vz and Bz are not plotted but exhibit analogous behaviours to vy and By. The correct solution, evaluated through the exact Riemann solver by Giacomazzo & Rezzolla (2006), is plotted with a dashed black line.

4.1.3 Loop advection

We next tested our code against the loop advection problem. This test is particularly challenging for grid-based codes due to the necessity to continuously evaluate strong fluxes, resulting from the fluid advection (Gardiner & Stone 2005). The Lagrangian nature of GIZMO, on the other hand, allows it to naturally follow the uniformly moving fluid, reducing the effective fluxes between particles to the minimum. For this set-up, we initialized 2.5 × 10⁵ particles on a Cartesian grid, filling a periodic 2D box of unit length. The fluid was given a uniform density (ρ = 1) and specific internal energy (ε = 4.5), with Γ = 5/3. The magnetic field was set to zero everywhere in the domain, but for a central circular region with radius R_loop = 0.3, in which we set an azimuthal magnetic field loop of the form

$$(B^x,B^y,B^z)=\{\begin{array}{ll}{A}_{\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{p}}(-y/r,x/r,0); & r<R_{loop}\\ (0,0,0);& r\ge R_{loop}\end{array},$$(37)

where A_loop = 10^-3 and $r=\sqrt{x^2+y^2}$ is the cylindrical radius. Finally, we assigned a bulk velocity, (v^x, v^y, v^z) = (1/2, 1/12, 0), to the fluid. We let the system evolve until t = 24 in the MFM mode.

In Fig. 4, we plot B^x and the magnetic pressure, P_mag = b²/2, at the initial and final times. Despite the presence of some smearing, due to the diffusive nature of our HLL solver the magnetic field intensity and shape are well conserved throughout the advection. As has also been observed in many other codes (Tóth & Odstrčil 1996; Stone et al. 2008; Pakmor et al. 2011; Fragile et al. 2019; Fambri et al. 2018; Cipolletta et al. 2020), a magnetic pressure drop in a small central region is present. We note that, during this test, our scheme keeps the magnetic field divergence well controlled, reaching maximum values of $h({\mathrm{\partial }}_i{\mathcal{B}}^i)$ at ~1% of the magnetic field intensity.

Fig. 4 Magnetic field loop advection performed in MFM mode. The left panels display the magnetic field intensity in the x direction multiplied by 103, whereas the right panels show the magnetic pressure Pmag = b2/2 multiplied by 107. Snapshots were taken at the initial time and at time t=24.

Fig. 5 Relativistic magnetic rotor problem at time t = 0.4, performed in MFM mode. We plot 2D maps for (from left to right) the rest mass density, the gas pressure, the magnetic pressure, and the Lorentz factor with magnetic field lines.

4.1.4 Magnetic rotor

As a second multi-dimensional stress test, we performed the relativistic version of the magnetic rotor (Del Zanna et al. 2003). A 2D disc of radius r_d = 0.1 and density ρ_d = 10 was initially rotating with a uniform angular velocity, Ω_d = 9.95. The rest of the domain was filled with a fluid of uniform density, ρ_ext = 1, initially at rest. The gas pressure was P = 1 everywhere (Γ = 5/3) and a magnetic field aligned with the x direction of magnitude B^x = 1 permeated the entire simulated box. This set-up usually represents a really difficult test for numerical schemes, due to the high Lorentz factor, W ≃ 10, the formation of low-density regions in the centre of the rotor, and the presence of strong discontinuities in the initial condition. To slightly alleviate this last point, we applied a smoothed transition on the rest mass density and on the velocity radial profiles between r_d and r_ext = 0.115 of the form f(r) = (r_ext - r)/(r_ext - r_d) so that

$$\rho (r)=\{\begin{array}{ll}{\rho }_d; & r\le r_d\\ 1+9f(r);& r_d<r<r_{ext}\\ {\rho }_0;& r\ge r_{ext}\end{array}$$(38)

and

$$|\mathbf{\text{v}}(r)|=\{\begin{array}{ll}r{\mathrm{\Omega }}_d; & r\le r_d\\ r_d{\mathrm{\Omega }}_{d}f(r);& r_d<r<r_{ext}\\ 0;& r\ge r_{ext}.\end{array}$$(39)

From these, we can compute the Lorentz factor, W(r), and the conserved density, D(r), radial profiles and integrate the cumulative mass function M(R < r). For the MFM run, we first placed ~1.2 × 10⁵ particles on a closely packed lattice and we then stretched their radial position to match the cumulative mass function M(R < r) up to r = r_ext, assuming equal-mass particles. Finally, outside of radius r_ext we initialized a closely packed lattice distribution of particles with uniform density ρ_ext. In the MFV case, we did not perform the radial stretch but we assigned the particle mass to match D(r) = ρ(r) W(r) instead, yielding an equal-volume particle distribution.

Our code can evolve the system without breaking in both MFM and MFV modes. We show our results at time t = 0.4 for the MFM case in Fig. 5. Our solution exhibits good behaviour overall: the centre of the disc is quickly emptied of gas, while an overdensity propagates outside. The presence of the magnetic field slows down the rotation of the fluid, yielding a maximum Lorentz factor of W ≃ 2.48 at the final time; due to the fluxfreeze effect, magnetic field lines are twisted by the fluid in the central region. When compared to results from grid-based codes (e.g. Del Zanna et al. 2003; Duffell & MacFadyen 2011; Mösta et al. 2014), our evolution scheme yields higher maximum values of the rest mass density and of the Lorentz factor, peaked on the oblate shear front. This behaviour is better quantified in Fig. 6, where we plot the most important quantities along horizontal and vertical slices, at the output time t = 0.4. For a more comprehensive view, we display here the results from the MFV run too, alongside the ones computed with fixed particle positions for reference. While our MFV results show a good agreement with the reference run and with the ones presented in Mösta et al. (2014), the MFM run yields higher rest mass density and Lorentz factor peaks by a factor ~1.5 and ~1.3, respectively; these features produce a slightly faster expansion of the rotor. We note that the resolution in regions with higher densities in our MFM run is -2 times higher than the MFV one (and then the resolution used in the referenced grid-based codes), which could therefore better resolve the compression of the fluid⁷.

Fig. 6 Magnetic rotor test 1D slices along y = 0.5 (left panels) and x = 0.5 (right panels) of (from top to bottom) rest mass density, gas pressure, magnetic pressure, Lorentz factor, and magnetic field divergence normalized to the magnetic field intensity, at time t = 0.4. We show results from the MFM run (red triangles), MFV (orange circles), and from a fixed particle position one (black squares) for reference.

4.1.5 Cylindrical blast wave

One of the most demanding special relativistic tests for a multidimensional code is the blast wave explosion. We tested our scheme against the magnetized cylindrical case presented in Komissarov (1999) and later performed as a test for various codes (Del Zanna et al. 2007; Beckwith & Stone 2011; Mösta et al. 2014; Cipolletta et al. 2020). In this problem an initial over-dense and overpressurized gas of ρ_in = 10^-2 and P_in = 1 was confined inside a radius of rin = 0.8. The external medium was instead uniform with rest mass density of ρ_ext = 10^–4 and pressure of P_ext = 3 x 10^-5. Between the inner r_in and an external radius, rext = 1, we applied a transition on the rest mass density of the form

$$\ln \rho (r)=\ln {\rho }_{in}+\left(\frac{r_{in}-r}{r_{in}-r_{ext}}\right)\ln \left(\frac{{\rho }_{ext}}{{\rho }_{in}}\right)$$(40)

and an equivalent one on the pressure profile. A magnetic field, B = (0.1, 0, 0), permeated the entire domain and the adiabatic index was set to Γ = 4/3.

As was done for the magnetic rotor test (Sec. 4.1.4), for the MFM run we initialized ~1.6 × 10⁵ equal-mass particles in a closely packed lattice and we then stretched their radial position to match the integrated mass cumulative function M(R < r) up to r = r_ext. Outside this radius, we did not perform any stretching, initializing the fluid with a uniform density, ρ_ext. For the MFV run, we placed -2.6 × 10⁶ particles on a closely packed lattice and assigned the particle mass to match the ρ(r) profile⁸. We let the system evolve until t = 3 in a periodic box of side length 6, using a Dedner damping parameter, K = 0.75. Our code is able to evolve this difficult test in both MFV and MFM modes, with and without the Dedner divergence cleaning scheme⁹.

We plot the MFM run results in Fig. 7. Our solution exhibits some oscillations in the direction aligned with initial magnetic field; their intensity is initially more relevant to the velocity evolution but later affects the other quantities. This is the same numerical artefact discussed in Sect. 4.1.2 for the Balsara1 test, caused by the strong volume compression of particles ahead of the shock front. Again, this effect can be attenuated by employing a lower-order reconstruction¹⁰ of the magnetic field or a larger kernel. In Fig. 8, we display the pressure distribution resulting from three additional runs. Our MFV run shows fewer oscillations outside the blast wave, due to the smaller volume gradient on the shock in the initial condition, but with a small background pressure enhancement. The equal-mass particle configuration evolved in the MFV mode and using a larger kernel (Wendland C4, Wendland 1995; Dehnen & Aly 2012) yields slightly smaller oscillations in amplitude but with more diffusion. Finally, we show our result evolved with fixed particle positions (as in a fixed grid set-up); this last plot is comparable with results from other grid-based codes (Mösta et al. 2014; Del Zanna et al. 2003).

In Fig. 9, we plot the results evaluated along horizontal and vertical slices. We can recognize a small drop in the magnetic pressure in the MFV equal-volumes run in the direction aligned with the initial magnetic field. The MFM run slightly underestimates the velocity of the shock in the x direction, yielding a lower Lorentz factor. The configuration evolved using the Wendland C4 kernel seems to produce the best solution overall (more similar to the fixed grid case).

Fig. 7 Cylindrical blast wave problem at time t = 3, performed in MFM mode. We plot 2D maps for (from left to right) the gas pressure, the Lorentz factor with magnetic field lines, the magnetic field x component and its y component.

Fig. 8 Pressure distributions in the blast wave problem for three different run: (left) MFV mode, equal-volume particles, using a cubic spline as kernel function; (centre) MFV mode, equal-mass particles, using a Wendland C4 as kernel function; (right) fixed particle positions run.

4.1.6 Spherical blast wave

As a final special relativistic test, we performed the spherical blast wave, an extension of the 2D cylindrical blast wave discussed in Sect. 4.1.5, another very challenging problem for GRMHD codes. We initialized the system as in Cipolletta et al. (2020), i.e. an inner region (r ≤ r_in = 0.8) was filled with a fluid with density, ρ_in = 10^-2, and pressure, P_in = 1. At r ≥ r_ext = 1, the gas density was ρ_ext = 10^-4 and the pressure P_ext = 3 × 10^-5. A smooth transition linked the two regions, as in the cylindrical case. The magnetic field, aligned with the z direction, was set to B = (0, 0, 0.1), and fills the whole domain. The adiabatic index was Γ = 4/3 everywhere. We initialized ~1.8 x 10⁶ equal-mass particles in a stretched closely packed lattice and we let the system evolve until t = 4, with a Dedner damping parameter, K = 0.75.

We plot the results of this test in Fig. 10. Because of the intrinsically larger number of neighbours employed in 3D, both the ‘cell’ volumes and the magnetic field divergence error (Eq. (25)) estimates improve in accuracy. Nonetheless, some spurious oscillations in the Lorentz factor ahead of the shock are still observed, affecting the density and pressure evolution. Despite these artefacts, which we expect can be alleviated by suitably tuning the kernel shape, extension, and the slope limiters employed (as was already discussed), the code is able to evolve this challenging test without breaking and without significantly ruining the solution within the shock front. Moreover, we do not observe any grid alignment effect or diagonal artefacts (usually found in other GRMHD codes (Del Zanna et al. 2003).

4.2 General relativistic magnetohydrodynamics

In this section, we present two general relativistic MHD tests that allow us to determine the robustness of our scheme in non-Minkowskian background spacetimes.

Fig. 9 Cylindrical blast wave test 1D slices along y = 0 (left panels) and x = 0 (right panels) of (from top to bottom) rest mass density, gas pressure, magnetic pressure, Lorentz factor, and magnetic field divergence in units of the magnetic field intensity, at time t = 3. We show results from the MFM run (red triangles), MFV with initially equal volumes particles (orange circles), MFV with initially equal masses particles evolved using the Wendland C4 kernel (cyan stars), and from a fixed-grid one (black squares) for reference.

Fig. 10 2D slices at t = 4 and y = 0 of the spherical blast wave problem, evolved in MFM mode. We report the gas pressure in the left panel and the Lorentz factor with the magnetic field lines distribution overlaid in the right panel.

4.2.1 Magnetized TOV star

As a first test, we considered the evolution of a magnetized, non-rotating Tolman-Oppenheimer-Volkoff (TOV) star (Oppenheimer & Volkoff 1939; Tolman 1939). With this problem, we can test the correct implementation of the terms proportional to the metric functions, $\sqrt{\gamma }$ and α, defined in the ADM formalism. The initial magnetized set-up is presented in many GRMHD codes (Duez et al. 2006; Liu et al. 2008; Mösta et al. 2014; Cipolletta et al. 2020) and represents a perturbation of the pure GRHD configuration, already tested in the relativistic GIZMO paper (Lupi 2023); hence, we used the same initial configuration with ~10⁶ particles, adding a poloidal magnetic field resulting from a purely toroidal vector potential, A_φ, given by

$$A_{\varphi }=A_b{\omega }^{2}max(P-P_{cut},0{)}^{\xi },$$(41)

where $\omega =\sqrt{x^2+y^2}$ is the cylindrical radius. Here we employed ξ = 2, which ensured that the pressure’s derivative was continuous at all radii, A_b = 15, to have a maximum magnetic-to-gas pressure parameter, β^-1 ~ 8 × 10^-4, and P_cut = 0.04P_max, which gives the radius at which the magnetic field goes to zero, where P_max is the maximum initial gas pressure (at the centre of the star). We then projected A_φ on the Cartesian directions as $A_i=\frac{\mathrm{\partial }{x}^{\varphi }}{\mathrm{\partial }{x}^i}{A}_{\varphi }=(-y{A}_{\varphi }/{\omega }^2,x{A}_{\varphi }/{\omega }^2,0)$ and computed the magnetic field as ${\mathcal{B}}^i=\sqrt{\gamma }{n}_{\mu }{\epsilon }^{\mu ijk}{\mathrm{\partial }}_j{A}_k=\sqrt{\gamma }\alpha {\epsilon }^{ijk}{\mathrm{\partial }}_j{A}_k$, obtaining

$${\mathcal{B}}^i=\sqrt{\gamma }\alpha A_b\{\begin{array}{l}\frac{xy}{r}{\mathrm{\partial }}_r\tilde{P};\\ \frac{yz}{r}{\mathrm{\partial }}_r\tilde{P};\\ -(\frac{{\omega }^2}{r}{\mathrm{\partial }}_r\tilde{P}+2\tilde{P});\end{array}$$(42)

where we defined $\tilde{P}\equiv max(P-P_{cut},0{)}^{\xi }$. We plot the resulting magnetic configuration (magnetic pressure distribution and magnetic field lines) on the x - z plane in Fig. 11. This prescription represents a small perturbation of the pure hydrodynamical equilibrium studied in Lupi (2023).

We evolved the TOV star using both MFM and MFV modes, with and without the Dedner cleaning scheme. We note that our solution is stable even when employing the Powell cleaning terms only, although we do not report the results here. In the fiducial case, which employs the complete divergence cleaning scheme, we employed K = 0.1 and f = 2.5; since we were evolving an equilibrium configuration with low magnetization, the estimated signal velocities were very slow, and hence there was a need to overestimate the cleaning speed by a higher factor, f. We note that evolving the system using f = 1 does not corrupt the GRMHD solution and only produces slightly higher ${\mathrm{\partial }}_i{\mathcal{B}}^i$ values.

We let the runs evolve for t ≃ 28t_dyn, where we defined the dynamical time, $t_{dyn}\equiv 1/\sqrt{{\rho }_c}$ (as in the GRHD test, the central rest mass density of the TOV star is ρ_c ≈ 0.129285). In Fig. 12 we show the radial rest mass density profile at three different snapshots, for both the MFV and MFM runs. Our results agree with the exact solution and we do not recognize big differences from the hydrodynamical case presented in Lupi (2023). At the end of the MFV run the central density has lowered by ~2% from its initial value, while we witness a small increase in the MFM run.

To better understand the evolution of the system we plot the central rest mass density and maximum magnetic field intensity evolution in Fig. 13. As was mentioned above, in the MFV case, the ρ_c evolution is comparable with the hydrodynamical case, exhibiting a slow decay over time; the magnetic field displays a similar behaviour: after an initial drop in the central region, the magnetic field intensity witnesses a slow decay, approaching ~30% at the end of the run. As is also studied in Cipolletta et al. (2020), this trend is related to the drop in the central rest mass density and stabilizes with increasing resolution. The initial part of the MFM evolution can be compared with the corresponding pure hydrodynamics run. However, after ~13 dynamical times, the magnetic field witnesses a topological rearrangement and an increase in its maximum intensity by a factor of ~1.6. When using the MFM scheme to evolve an equilibrium configuration, setting the frame velocity in the Riemann problems to cancel out the mass fluxes introduces small perturbations in the fluid velocity. In a pure GRHD scenario, this effect does not pose significant issues, as the quantities readjust correctly without notable consequences. In the GRMHD case, instead, these perturbations significantly affect the evolution of the magnetic fields, leading to the aforementioned increase; as a consequence, the central rest mass density also increases by ~2%, with a delay. In order to assess the long-term effect of this hitch, we let the system evolve up to t ≈ 50t_dyn. Despite the artificial rearrangement of the magnetic field, the configuration does not appear to diverge; the TOV star remains stable and both the B maximum value and the central density eventually (slowly) decay towards their initial value, albeit with a noisier spatial configuration.

Fig. 11 Initial magnetic field configuration of the TOV test. We plot the magnetic pressure with superimposed magnetic field lines. The red circle marks the star radius (defined as the point where P = 10-8 Pmax), while the blue one indicates the radius at which the magnetic field goes to zero. The magnetic field is zero in white-coloured regions.

4.2.2 Magnetized Bondi accretion

As a final test, we considered the magnetized spherical accretion flow onto a non-spinning BH with mass M_BH, as a GRMHD extension of the test performed in Lupi (2023). The solution represents the relativistic generalization of the Bondi accretion solution (Michel 1972; Hawley et al. 1984), with the addition of a radial magnetic field given by

$${\mathcal{B}}^r=\frac{B_0{M}_{BH}^2}{r^2},$$(43)

where r is the coordinate radius and B₀ is a free parameter that controls the magnetic field strength. This spatial configuration of B does not alter the hydrodynamical solution and satisfies the divergence-free constraint (11).

We recall that the steady-state solution can be written as (Liptai & Price 2019)

$$\begin{array}{rl}& u^r(r)=\frac{C_1}{r^2{\zeta }^n(r)};\\ & \rho (r)=K_0{\zeta }^n(r);\\ & ϵ(r)=n\zeta (r),\end{array}$$(44)

where $\zeta \equiv (\mathrm{\Gamma }-1)ϵ$ and n ≡ 1/(Γ - 1). By assuming a critical radius, r_c, where

$$\begin{array}{rl}& u_c\equiv u^r(r_c)=\sqrt{\frac{M_{BH}}{2r_c}};\\ & v_c\equiv v^r(r_c)=\sqrt{\frac{u_c^2}{1-3u_c^2}};\\ & {\zeta }_c\equiv \zeta (r_c)=\frac{n{v}_c^2}{1+n(1-v_c^2)-n^2{v}_c^2},\end{array}$$(45)

one can determine the coefficients,

$$\begin{array}{rl}& C_1=u_c{r}_c^2{\zeta }_c^n;\\ & C_2=(1+(n+1){\zeta }_c{)}^2(1-\frac{2M_{BH}}{r_c}+u_c^2),\end{array}$$(46)

and retrieve ζ(r) by solving

$$C_2=(1+(n+1)\zeta (r){)}^2\{1-\frac{2M_{BH}}{r}+[u^r(r){]}^2\}.$$(47)

Using Eq. (44) and taking into consideration the curved metric tensor, we then computed the Lorentz factor by solving

$$W=\frac{1}{\sqrt{1-{\gamma }_{ij}{v}^i{v}^j}}=\frac{1}{\sqrt{1-{\gamma }_{ij}(\frac{u^i}{W}+\frac{{\beta }^i}{\alpha })(\frac{u^j}{W}+\frac{{\beta }^j}{\alpha })}}$$(48)

and the coordinate velocity, ${\tilde{v}}^i=\alpha u^i/W$. Finally we evaluated $D(r)=\sqrt{\gamma }\rho W$.

Unlike Lupi (2023), here we work with ‘Cartesian’ Kerr-Schild (KS) coordinates (Kerr 1963; Visser 2007; Blakely et al. 2015), as they are horizon-penetrating, i.e. the metric is continuous across the event horizon at r = 2M. In KS coordinates, the general (a ≠ 0) metric can be expressed as

$$g_{\mu \nu }={\eta }_{\mu \nu }+\frac{2M_{BH}{r}^3}{r^4+a^2{z}^2}{l}_a{l}_b,$$(49)

where η_μν is the Minkowski flat metric, a is the spin parameter,

$$r(x,y,z)=\sqrt{\frac{R^2-a^2+\sqrt{(R^2+a^2{)}^2+4a^2{z}^2}}{2}}$$(50)

is the coordinate radius as a function of the spherical radius, $R=\sqrt{x^2+y^2+z^2}$, and

$$l_a=(1,\frac{rx+ay}{r^2+a^2},\frac{ry-ax}{r^2+a^2},\frac{z}{r}).$$(51)

In our test problem, we set a = 0; hence, r = R. Note that, in this coordinate system, the shift vector, β_i = g_i0, differs from zero even for a non-spinning BH, allowing us to directly test the correct implementation of the shift vector-dependent terms in our scheme.

For our initial condition, we considered a critical radius of r_c = 8M, a K₀ value that gives ρ(r_c) = 1/16, and an adiabatic index of Γ = 4/3. We performed runs with two different magnetic field intensities, B₀ = (2.0, 6.25), yielding a magnetic-to-gas pressure parameter of β^-1 ≈ (0.1, 1.1) at the critical radius and β^-1 ≈ (4.7,46) at the BH horizon. We initially placed ~2 × 10⁷ particles on a closely packed lattice ranging from r = 1 M up to r = 100M and then stretched their radial position to match the exact density profile. During the evolution, we set the particle excision radius at r_exc = 1 M, to prevent them from reaching the central singularity.

In the central regions, where the fluid magnetization is stronger, we notice an overestimation of the gas internal energy when Eq (7) is employed, which corrupts the solution (Del Zanna et al. 2007). Since we do not expect shocks in this test, i.e. the entropy should remain constant, in these central regions we activate the energy-entropy switch from Eq. (23), and recover the internal energy from the entropy function, $s\equiv P/{\rho }^{\mathrm{\Gamma }}$.

To avoid numerical issues inside the event horizon, we also dampen the magnetic field and the velocity components via the function

$$\chi =\{\begin{array}{ll}\exp [-\frac{\sqrt{r_{exc}}(r-r_{EH}{)}^2}{(r_{EH}-0.4)}]; & r\le r_{EH}\\ 1;& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$$(52)

where r_exc and r_EH are the excision and the event horizon radii, respectively.

We report in Fig. 14 the results from runs performed in MFV and MFM modes, for both the low and high magnetization cases. We report radial profiles of the rest-mass density, ρ, the radial coordinate velocity, ${\tilde{v}}^r$, the specific internal energy, ε, and the radial component of the conserved magnetic field, B^r , at t = 100M, in the range r ε [2,30]M. Our solutions are in very good agreement with the exact solution. We recognize a small scatter in the radial velocity profile (though more moderate with respect to the hydrodynamical case) and a small overestimation of the same quantity close to r = 2M in the strong magnetization case. Thanks to the horizon penetrating metric, our solution is smooth through r = 2M and does not show any drop in density or specific internal energy near this radius. For completeness purposes, in Appendix A we discuss the same problem evolved without the energy-entropy switch.

Fig. 12 Radial profile of TOV equilibrium at t = 0, t ≈ 16tdyn and t ≈ 28tdyn for the MFM (purple circles) and the MFV (red triangles) run. The dashed black line marks the exact solution. The MFM solution at time t ≈ 50tdyn is displayed as a dashed lime line.

Fig. 13 Evolution of the central TOV rest mass density (top) and the maximum value of the magnetic field (bottom) for the MFM (purple) and MFV (red) run.

Fig. 14 Rest-mass density, ρ, coordinate velocity, ${\tilde{v}}^r$, specific internal energy, ε, and magnetic field, Br, radial profile of the magnetized Bondi accretion test, performed both in the MFM (purple circles) and MFV (orange triangles) modes, at time t = 100M. We plot results of the run performed with low (top row) and high magnetization (bottom row). In each panel, we display the exact solution with a dashed line.

5 Discussions and conclusions

In this paper, we have presented a novel GRMHD scheme implemented in the code GIZMO, which works with the MFV and MFM modes and builds on the GRHD implementation presented in Lupi (2023). The implementation incorporates a general relativistic single-state HLL Riemann solver that considers magnetic fields, a more stable reconstruction of the velocity field at the faces between interacting particles, and an energyentropy switch, which is employed when the fluid is supersonic, strongly magnetized, and isentropic. In the generalized kickdrift-kick scheme, we introduced an explicit time evolution of the conserved quantities that takes into account the particle volume variation during a timestep; in this way, we consistently pass the local, un-densitized evolved variables to the GRMHD conservative-to-primitive routine. To control the magnetic field divergence errors, we implemented a general relativistic version of the Powell eight-wave (Powell et al. 1999) and of the Dedner hyperbolic divergence cleaning (Dedner et al. 2002). We tested and calibrated these routines and verified that the scheme can evolve the system with reasonably low $\mathrm{\nabla }\cdot B$ values.

We performed several benchmark tests, starting with planar relativistic MHD shock tube problems from Balsara (2001) and then testing the magnetic field loop advection and the magnetic rotor set-up from Del Zanna et al. (2003), in two dimensions. As more stringent tests, we performed the magnetized cylindrical (2D) and spherical (3D) blast waves, both in the MFV and the MFM code modes, also employing different kernel functions. From this set of multi-dimensional special relativistic magnetized tests, we conclude that our code can stably evolve various problems yielding good results. We observe no dependence on the geometry of the problem and the magnetic field divergence is always kept under safe levels. However, we recognize some moderate inaccuracies when dealing with violent shocks characterized by a strong density or pressure gradient and high magnetization. These artefacts are enhanced when a sharp particle volume gradient is initialized.

We then challenged our scheme against GRMHD tests, such as the stability of the magnetized TOV and the magnetized spherical accretion onto a non-rotating BH. In the first one, we confirm that the MFV mode is better at maintaining the magnetic field topology intact, even though the code can preserve the star hydrodynamical equilibrium in both modes. For the second test, we implemented the horizon-penetrating KS coordinate system and, employing the energy-entropy switch, we achieved an accurate steady-state solution.

We finally note that our GRMHD implementation can work with any EoS and metric provided by the user, while the ideal fluid EoS, the flat and Kerr metrics (both in Boyer-Lindquist and KS coordinates), are already present in the code. Our scheme will be made public in due course. In forthcoming papers, we shall include dynamic metrics (Combi et al. 2021; Combi & Ressler 2024) and a more sophisticated Riemann solver (Mignone & Bodo 2006; Miyoshi & Kusano 2005; Mignone et al. 2009).

Acknowledgements

Some simulations were run on the Leonardo cluster at CINECA (allocation INF24_teongrav) thanks to the CINECA-INFN agreement. The authors thank the reviewer, Daniel J. Price, for his valuable comments and suggestions that improved the quality of the article. GF thanks Bruno Giaco-mazzo for his advice about the TOV equilibrium test, and Andrea Mignone and Vittoria Berta for helpful discussions. AL acknowledges support from PRIN MUR “2022935STW”.

Appendix A Magnetized Bondi accretion without energy-entropy switch

We report here the results for the magnetized Bondi accretion evolved without the energy-entropy switch discussed in Sec. 3.2.3. The initial set-up is identical to the one presented in Sec. 4.2.2. Here, we evolve the conserved total energy τ through the usual conservation Eq. (7), instead of using the entropy conservation (23). In this case, we let the system evolve until t = 10M and we plot the results in Fig. A.1. Despite the earlier time, we already observe a significant increase (1.5x) of the specific internal energy at the event horizon, in the MFV run. As mentioned in the main text, this effect is due to an incorrect retrieval of the specific internal energy from the conserved energy, when dealing with supersonic flows, strong gravitational fields and relatively high magnetic field pressures b²/2. Despite this hitch, the rest-mass density, the coordinate velocity, and the magnetic field remain consistent with the expected solution overall, but for a small scatter around r = 10M. A subtle underestimation of the radial velocity is found near the event horizon, which is a direct consequence of the incorrect internal energy density estimation. The MFM run yields good results overall, with only a larger scatter in the internal energy and the magnetic field when compared with Fig. (14).

Fig. A.1 Radial profiles of the rest-mass density ρ, the radial coordinate velocity ${\tilde{v}}^r$, the specific internal energy ε, and the radial component of the magnetic field Br of the magnetized Bondi accretion test, performed without the energy-entropy switch. We show the results obtained with the MFV (orange triangles) and MFM mode (purple circles), at time t = 10M in the low-magnetization case. The exact solution is reported in each panel by a dashed black line.

Appendix B Monopole test performed with the Powell scheme only

To better show how our divergence cleaning scheme handles spurious magnetic field divergences, we performed the monopole test from Sec. 4.1.1 employing the Powell terms in Eq. (24) only. The initial set-up is analogous to the one presented in the main text. We let the system evolve until t = 10 and we plot the resulting ${\mathrm{\partial }}_i{\mathcal{B}}^i$ in Fig. B.1, at four different times. At the beginning of the evolution, the Powell scheme is highly effective and quickly damps part of the monopole. After t = 5 it starts to slow down due to the less pronounced B gradients entering Eq. (24). By the final time, the monopole is damped by an order of magnitude and, as expected, we see no advection of the divergence over time. The small ${\mathrm{\partial }}_i{\mathcal{B}}^i$ diffusion along the y-direction is a consequence of the induction and the moment density equations reacting to the presence of an initial non-zero magnetic field along the x-direction. The overall behaviour of this solution was expected; in fact, while the Powell prescription can in principle correct large magnetic field divergences over time, it is better suited to remove the instantaneous formation of spurious errors emerging from the fluxes evaluation. The cleaning of a strong monopole, such as the one in this test, is best handled by the Dedner hyperbolic cleaning.

Fig. B.1 Evolution of the monopole test performed with the Powell scheme only. We plot a 2D slice at z = 0 of the magnetic field divergence at the initial time and at three later times. Note the different colour-map scale with respect to Fig. 1.

References

All Figures

	Fig. 1 Evolution of the monopole test. 2D slice at z = 0 of the magnetic field divergence at the initial time (left) and at three later times. The top panels display the solution with damping parameter K = 0.1 and the bottom panels the one with K = 1.
In the text

Fig. 2 Slice at y = y0 of primitive quantities at time t = 0.4 of the Balsara1 test. The magnetic field divergence in units of the magnetic field intensity is plotted in the seventh panel. Blue squares (Powell) and red triangles (Dedner) mark the MFM solutions, while cyan stars (Powell) and orange circles (Dedner) indicate the MFV ones. In the bottom panel, we display the specific internal energy when the problem is initialized with equal-mass particles and evolved with the MFV scheme; gold diamonds mark the solution employing standard slope limiters, while olive triangles show the results using more diffusive slope limiters. The correct solution, evaluated through the exact Riemann Solver by Giacomazzo & Rezzolla (2006), is plotted with a dashed black line.

In the text

Fig. 3 Slice at time t = 0.4 and constant y of the primitive quantities from the Balsara4 test performed in the MFM mode. The magnetic field divergence in units of its intensity is plotted in the bottom panel. Blue squares mark the solution computed with Powell terms only, while red triangles indicate the one computed with the complete divergence cleaning scheme. The high-resolution run is displayed with orange circles. Quantities vz and Bz are not plotted but exhibit analogous behaviours to vy and By. The correct solution, evaluated through the exact Riemann solver by Giacomazzo & Rezzolla (2006), is plotted with a dashed black line.

In the text

	Fig. 4 Magnetic field loop advection performed in MFM mode. The left panels display the magnetic field intensity in the x direction multiplied by 103, whereas the right panels show the magnetic pressure Pmag = b2/2 multiplied by 107. Snapshots were taken at the initial time and at time t=24.
In the text

	Fig. 5 Relativistic magnetic rotor problem at time t = 0.4, performed in MFM mode. We plot 2D maps for (from left to right) the rest mass density, the gas pressure, the magnetic pressure, and the Lorentz factor with magnetic field lines.
In the text

Fig. 6 Magnetic rotor test 1D slices along y = 0.5 (left panels) and x = 0.5 (right panels) of (from top to bottom) rest mass density, gas pressure, magnetic pressure, Lorentz factor, and magnetic field divergence normalized to the magnetic field intensity, at time t = 0.4. We show results from the MFM run (red triangles), MFV (orange circles), and from a fixed particle position one (black squares) for reference.

In the text

	Fig. 7 Cylindrical blast wave problem at time t = 3, performed in MFM mode. We plot 2D maps for (from left to right) the gas pressure, the Lorentz factor with magnetic field lines, the magnetic field x component and its y component.
In the text

	Fig. 8 Pressure distributions in the blast wave problem for three different run: (left) MFV mode, equal-volume particles, using a cubic spline as kernel function; (centre) MFV mode, equal-mass particles, using a Wendland C4 as kernel function; (right) fixed particle positions run.
In the text

Fig. 9 Cylindrical blast wave test 1D slices along y = 0 (left panels) and x = 0 (right panels) of (from top to bottom) rest mass density, gas pressure, magnetic pressure, Lorentz factor, and magnetic field divergence in units of the magnetic field intensity, at time t = 3. We show results from the MFM run (red triangles), MFV with initially equal volumes particles (orange circles), MFV with initially equal masses particles evolved using the Wendland C4 kernel (cyan stars), and from a fixed-grid one (black squares) for reference.

In the text

	Fig. 10 2D slices at t = 4 and y = 0 of the spherical blast wave problem, evolved in MFM mode. We report the gas pressure in the left panel and the Lorentz factor with the magnetic field lines distribution overlaid in the right panel.
In the text

	Fig. 11 Initial magnetic field configuration of the TOV test. We plot the magnetic pressure with superimposed magnetic field lines. The red circle marks the star radius (defined as the point where P = 10-8 Pmax), while the blue one indicates the radius at which the magnetic field goes to zero. The magnetic field is zero in white-coloured regions.
In the text

	Fig. 12 Radial profile of TOV equilibrium at t = 0, t ≈ 16tdyn and t ≈ 28tdyn for the MFM (purple circles) and the MFV (red triangles) run. The dashed black line marks the exact solution. The MFM solution at time t ≈ 50tdyn is displayed as a dashed lime line.
In the text

	Fig. 13 Evolution of the central TOV rest mass density (top) and the maximum value of the magnetic field (bottom) for the MFM (purple) and MFV (red) run.
In the text

Fig. 14 Rest-mass density, ρ, coordinate velocity, ${\tilde{v}}^r$, specific internal energy, ε, and magnetic field, Br, radial profile of the magnetized Bondi accretion test, performed both in the MFM (purple circles) and MFV (orange triangles) modes, at time t = 100M. We plot results of the run performed with low (top row) and high magnetization (bottom row). In each panel, we display the exact solution with a dashed line.

In the text

Fig. A.1 Radial profiles of the rest-mass density ρ, the radial coordinate velocity ${\tilde{v}}^r$, the specific internal energy ε, and the radial component of the magnetic field Br of the magnetized Bondi accretion test, performed without the energy-entropy switch. We show the results obtained with the MFV (orange triangles) and MFM mode (purple circles), at time t = 10M in the low-magnetization case. The exact solution is reported in each panel by a dashed black line.

In the text

	Fig. B.1 Evolution of the monopole test performed with the Powell scheme only. We plot a 2D slice at z = 0 of the magnetic field divergence at the initial time and at three later times. Note the different colour-map scale with respect to Fig. 1.
In the text