Dissipative Magnetohydrodynamics for Non-Resistive Relativistic Plasmas:
An implicit second-order flux-conservative formulation with stiff relaxation

In the following we summarize a formalism for describing out-of-equilibrium dynamics allowing for viscous corrections to the system. These can be encoded by the viscous stress tensor $\Pi_{\mu\nu}$, which is added to the fluid’s energy-momentum tensor. In this work, we use the so-called Eckart hydrodynamic frame Eckart (1940) in which the rest-mass current remains $N^{\mu}=\rho u^{\mu}$, all viscous/dissipative effects appear only in the expression for the energy-momentum tensor, and there are no out-of-equilibrium corrections to the energy density.

Considering the number of degrees of freedom, we can see that $T^{\mu\nu}$ and $N^{\mu}$ have a total of 14 (since $T^{\mu\nu}$ is symmetric), but in equilibrium only five of them are actually independent (e.g. $\{e,\rho,u^{\mu}\}$). Accounting for out-of-equilibrium dynamics of the system then amounts to prescribing evolution equations for the remaining 9 degrees of freedom. These can be parameterized in terms of the bulk scalar pressure $\Pi$, the symmetric anisotropic pressure tensor $\pi^{\mu\nu}$, and the heat flux $q^{\mu}$, which obey a distinct set of constraints on the dissipative momenta and stresses, i.e.

Therefore, once these constraints are fulfilled, the set $\{\pi_{\mu\nu},\Pi,q_{\mu}\}$ indeed only has 9 independent degrees of freedom. As we show in this paper, the way these constraints are imposed heavily influences the numerical methods required to study this system. Overall, we can group these contributions into the viscous stress tensor

which is then added to the hydrodynamics equations

The combined stress-energy tensor, including the dissipative contributions, is covariantly conserved

whereas $N^{\mu}$ remains unchanged with dynamics given by (3). We point out that we choose the Eckart frame in this work solely for numerical convenience in the context of astrophysical applications. Clearly, other definitions of the hydrodynamic fields, i.e. other hydrodynamic frames Israel and Stewart (1979), could have been used. For instance, one could have removed the heat flux and introduced a diffusive correction to the baryon current. Alternatively, one may employ other frames in which both terms are present. For a deeper discussion about hydrodynamic frames, and how a judicious choice of hydrodynamic variables can be useful to formulate causal, stable and strongly-hyperbolic theories of viscous relativistic fluids at first-order in derivatives, see Bemfica et al. (2018); Kovtun (2019); Bemfica et al. (2019b); Hoult and Kovtun (2020); Bemfica et al. (2020). Second-order Israel-Stewart-like theories in general hydrodynamic frames have also been investigated in Tsumura and Kunihiro (2010), and more recently in Noronha et al. (2021); Rocha and Denicol (2021).

Having parameterized the viscous degrees of freedom, we need to prescribe a set of causal evolution equations for them. In this work, we use equations of motion obtained from a truncation of suitably defined moments of the relativistic Boltzmann-Vlasov equation for a gas of charged particles (without dipole moment or spin) dynamically coupled to a magnetic field. This approach to obtain the second-order hydrodynamic equations employs a systematic expansion of the Boltzmann equation in powers of the Knudsen and inverse Reynolds numbers¹¹1The Knudsen numbers measure the ratio between microscopic length scales and the macroscopic scales associated with gradients of the conserved quantities. In the approach of Denicol et al. (2012), one denotes quantities such as $\Pi/(e+P)$ as inverse Reynolds numbers, given that they can be seen as relativistic generalizations of the ratio between viscous and inertial forces., and this method was originally applied in Denicol et al. (2012) in the investigation of neutral kinetic systems. The equations of motion, and the transport coefficients obtained within this formulation, are widely used in current heavy-ion collision simulations of the quark-gluon plasma (usually done in the absence of eletromagnetic fields). This formalism was extended to include non-resistive, comoving magnetic field effects using the Boltzmann-Vlasov equation in Denicol et al. (2018). The magnetohydrodynamical equations of motion were obtained in the so-called 14-moment approximation Israel and Stewart (1979) and they essentially provide a generalization of Israel-Stewart fluid dynamics to the case of non-vanishing magnetic fields. The full resistive case (which includes effects from comoving electric and magnetic fields) was worked out later in this approach in Ref. Denicol et al. (2019).

In our work, the dissipative variables obey the following advection-type relaxation equations²²2We note that the original derivation of the non-resistive equations of motion in Denicol et al. (2018) used the Landau frame Landau and Lifshitz (1987). The Eckart frame version of these equations can be found by taking the non-resistive limit of the equations derived in Denicol et al. (2019).

where we have introduced the symmetric trace-free projector

Above, $\eta$ is the shear viscosity, $\zeta$ is the bulk viscosity, $\kappa$ is the thermal conductivity, $\tau_{\Pi}$ is the bulk relaxation time, $\tau_{q}$ is the relaxation time for heat flux, while $\tau_{\pi}$ is the shear viscosity relaxation time. Besides those coefficients, the equations of motion also contain the additional second-order transport coefficients $\{\delta_{\Pi\Pi},\delta_{qq},\delta_{\pi\pi}\}$ and the new coefficients that determine the coupling between the dissipative quantities and the electromagnetic field tensor, $\delta_{qB}$ and $\delta_{\pi B}$. Expressions for all of these coefficients, including $\delta_{qB}$ and $\delta_{\pi B}$, in terms of the temperature and chemical potential for an ultrarelativistic gas can be found in Denicol et al. (2012, 2019, 2018). In this work, we treat the ratios $\delta_{\Pi\Pi}/\tau_{\Pi}$,$\zeta/\tau_{\Pi}$, $\delta_{qq}/\tau_{q}$, $\kappa/\tau_{q}$, $\delta_{\pi\pi}/\tau_{\pi}$, $\eta/\tau_{\pi}$, $\delta_{\pi B}/\tau_{\pi}$, $\delta_{qB}/\tau_{q}$ as free parameters defined in terms of thermodynamic quantities ($e$ and $\rho$). It is important to remark that in the formulation of Ref. Denicol et al. (2018) the form of the equations of motion remains very close to that found in standard Israel-Stewart theory, with additional terms that couple the fluid variables directly to the magnetic field. Also, we note that in this first work we have not included all the possible second-order terms that appear in Denicol et al. (2012, 2018, 2019). A systematic investigation of their effects will be given in a forthcoming study.

We note that any second-order theory with relaxation equations for the dissipative currents given by (19), (20), and (21) has an asymptotic relativistic Navier-Stokes regime. This can be seen by performing a systematic expansion in gradients, together with the conservation laws Baier et al. (2008). In the absence of magnetic fields, to first order in gradients the system reduces to

which is the standard relativistic generalization of Navier-Stokes theory Eckart (1940). It is important to understand that this limit can never be reached exactly, since it reduces the character of the hydrodynamic evolution equations from hyperbolic to parabolic, rendering the system acausal Pichon (1965). Also, we note that the standard constraints that stem from a linearized analysis of causality and stability, derived in the absence of a magnetic field in Hiscock and Lindblom (1983), automatically hold for the system of equations we use. For a linearized study of stability and causality in the presence of a nonzero magnetic field, see Biswas et al. (2020).

We start out by considering the Israel-Stewart-like equation (19) for the bulk pressure $\Pi$. After a simple algebraic manipulation it can be shown that Eq. (19) is equivalent to

Here we have introduced the shorthand $\theta=\nabla_{\mu}u^{\mu}$. We will provide a detailed treatment of arbitrary viscous relaxation terms in Sec. III.5. We note that depending on the application this equation can be further simplified.

In the simplest case of $\zeta/\tau_{\Pi}=\text{ const}$, or if $\zeta/\tau_{\Pi}$ is advected, i.e. $u^{\mu}\nabla_{\mu}\left(\zeta/\tau_{\Pi}\right)=0$ and $\delta_{\Pi\Pi}=\tau_{\Pi}$, we find

which, remarkably, is free of derivatives on the right-hand side (RHS).

Next we turn to the equation for heat conduction. It can trivially be shown that

where we have used that $u_{\alpha}\nabla_{\mu}u^{\alpha}=0$. The specific form of this equation further allows us to add terms proportional to $u^{\nu}$, since these get projected out by the global $\Delta^{\alpha}_{\nu}$ projection. Along those lines, it will help us to simplify the equations if we add such a terms in the following way

We can now understand the meaning of the projector $\Delta^{\mu}_{\nu}$ in Eq. (29) as follows. Contracting Eq. (29) with $u_{\alpha}$, we see that the sole purpose of the projector is to enforce Eq. (13), such that $q_{\mu}$ is orthogonal to $u^{\mu}$. In practice, the presence of such a projector will lead to a variety of gradient terms, which explicitly impose the constraint (29). Such terms are typically numerically difficult to evaluate and will not only increase the error budget of any numerical solution, but they will also numerically lead to only an approximate enforcement of Eq. (13).
A better approach, which is commonly used in relativistic force-free electrodynamics Alic et al. (2012); Palenzuela (2013), is to instead include the constraint (29) via stiff relaxation Alic et al. (2012). Introducing a new relaxation time $\omega_{q}\ll\tau_{q}$, we can equivalently write

The effect of this term can best be considered in the Navier-Stokes limit (in the absence of magnetic fields), where

Hence, the constraint (13) will indeed be imposed when $\omega_{q}\rightarrow 0$. Although Eq. (31) is more desirable from a numerical point of view than the original Israel-Stewart equation (20), we will need implicit numerical schemes in order to properly compute this term. This will be explained in detail in Sec. IV. Applying the same reordering that leads to a flux-conservative form (26) for the bulk scalar we obtain the final form

For a class of theories, where $\kappa/\tau_{q}=\text{const}$ or advected $\left(u^{\mu}\nabla_{\mu}\left[\kappa/\tau_{q}\right]=0\right)$, and $\delta_{qq}=\tau_{q}$, we find the simple flux-conservative stiff relaxation equation

Nonetheless, our numerical scheme is able to solve systems with arbitrary heat conductivities, see Sec. III.5.

We finally consider the evolution of the shear stress tensor $\pi^{\alpha\beta}$, with the aim of recasting it into flux-conservative form. Starting from Eq. (21), we write

where $2\nabla_{\left(\alpha\right.}u_{\left.\beta\right)}=\nabla_{\alpha}u_{\beta}+\nabla_{\beta}u_{\alpha}$. Similar to the observations made for Eq. (29), we find that the introduction of the trace-free projector $\Delta^{\alpha\beta}_{\lambda\nu}$ was done to ensure the validity of the constraints (14) and (15). In the same spirit of deriving Eq. (31), we can replace the projector by the introduction of a stiff relaxation current. More specifically, we write

where we have introduced the relaxations times $\omega_{\pi},\omega_{\pi\pi}\ll\tau_{\pi}$. Reordering the terms in the same way as for Eqs. (26) and (33), we find

We further note that all terms terms proportional to $u^{\mu}$ are already accounted for in the constraint damping term and will get removed in the $\omega\rightarrow 0$ limit. Remarkably, and differently from the heat conduction and bulk viscosity, this leads to the removal of the $\nabla_{\mu}\left(\eta/\tau_{\pi}\right)$ term from the equations. Hence, no approximation or other special treatment need to be applied to handle non-constant shear viscosity. The final form of the evolution equation then reads

As a final remark we stress the ambiguity in adding or removing constraint terms related to Eq. (14) and (15). While in the continuum limit Eq. (15) will perfectly hold, the discrete version of these equations might behave differently whether or not additional terms proportional to $g^{\mu\nu}$ would be added. One particular example would be the removal of the trace from the shear tensor in Eq. (43). This would lead to a replacement of the principal part of Eq. (43) with

Whether or not such modifications will impact the stability of the system will require a more in-depth analysis, which will be provided in forthcoming work.

In the following, we would like to associate a meaning with the magnetic field coupling terms, $\delta_{qB}$ and $\delta_{\pi B}$. To facilitate this, we will draw on an analogy with the Ohm’s law for a single fluid ideal magnetohydrodynamics in the Newtonian regime. To this end, we consider a simple Ohm’s law from Hall magnetohydrodynamics Sturrock and Andrew (1994),

where $\omega_{g}=q/\left(m_{e}cB\right)$ is the electron gyrofrequency, $q$ the electron charge, $m_{e}$ the electron mass and $J^{i}$ the electric current. By fixing the gyration velocity $v_{\perp}$, we can also express this via the inverse Larmor radius $R^{-1}_{L}=\omega_{g}/v_{\perp}$.
Taking the Newtonian limit of Eq. (20), i.e. $\tau_{q}\rightarrow 0$, and going to the rest-frame of the fluid $u^{\mu}=\left(1,0,0,0\right)$, we obtain

Since we only consider non-resistive plamas, we can see that the heat conduction current replaces the electric current compared to Eq. (45). In fact, if we were to consider resistive plasmas Denicol et al. (2019), we would obtain exactly the same expression with $\sigma^{\prime}e^{i}$ on the RHS, where $e^{\mu}$ is the comoving electric field, and $\sigma^{\prime}$ the associated conductivity.

Although Eq. (45) describes an electric and Eq. (46) a heat current, it can be shown that the two are related in a dissipative relativistic fluid description. Indeed, the total heat flux of the theory is given by Denicol et al. (2019)

where $m_{b}$ is the baryon mass and $V_{f}^{\mu}$ is the particle diffusion current, which directly enters the electric current in Maxwell’s equations Denicol et al. (2019). In defining the dissipative system in Eq. (16), we have made a frame choice, to neglect the particle diffusion current, and treat the heat flux vector $q^{\mu}$ directly. We could have additionally chosen to operate in the Landau frame, and instead reexpressed the system in terms of the diffusion current $V_{f}^{\mu}$, being more similar to Eq. (45).

Overall, this comparison allows us to identify the coupling scale with a thermal gyrofrequency

We further find it natural to identify the time scale of collisions $\tau_{e}$ with the relaxation time $\tau_{q}$, as those should be proportional to each other. Physically this implies that the thermal conductivity (in the Navier-Stokes limit) will split into a part parallel $\left(\kappa_{\parallel}\right)$ and perpendicular $\left(\kappa_{\perp}\right.$, $\left.\kappa_{\rm H}\right)$ to the magnetic field, where the latter term is the equivalent of the electric Hall conductivity Denicol et al. (2018). In particular, this would lead to an anisotropic modification of the Navier-Stokes limit as follows

Crucially, in our current prescription these anisotropy effects are controlled by a single parameter $\delta_{qB}$, which can be freely specified. For a more detailed discussion of this effect, see Ref. Denicol et al. (2018).

Although we have so far managed to eliminate most gradient terms from the RHS of the equations of motion, a set of advection and compression terms still remain. We will now, likewise, reformulate them as a set of implicit rate equations in the comoving frame.

We note that because of baryon number conservation, an arbitrary advected scalar $Y$ satisfies

If we want to enforce that $Y$ follows a certain behavior, $Y_{0}$, we may implicitly define an autonomous source term such that

This gives rise to a relaxation current

where $\omega_{Y}$ is a stiff relaxation time scale. In order to enforce the condition (57) we may add this current to the advected part,

In other words, in the local comoving frame of the fluid, this equation is prescribing a rate equation to enforce the damping of $Y$ towards $Y_{0}$ on a (subgrid) time scale $\omega_{Y}$. If we can evaluate the stiff current $\mathcal{I}_{Y}$ numerically, it can be used to replace advective derivatives of the form

In the same way, we can also treat gradient terms by noting that

where $\delta^{\nu}_{\mu}$ is the Kronecker symbol. Introducing the current $Z^{\nu}$, we can write in non-covariant form:

where $\Gamma^{\lambda}_{\mu\nu}$ is the Christoffel symbol associated with $g_{\mu\nu}$ and

Here $\omega_{Z}$ is a stiff relaxation timescale.

Applied to our source terms this results in the following set of stiff advection equations

The source terms of those equations are fixed according to

where we have introduced new relaxation time scales $\omega_{\zeta}\,,\omega_{\kappa}\,,\omega_{\eta}\,,\text{and}\,\omega_{\theta}$.

In conclusion, our equations of motion for flux-conservative dissipative magnetohydrodynamics read

Before we proceed, a few comments are in order. All second-order gradient terms $\left(\Pi\theta\,,q^{\nu}\theta\,,\pi^{\mu\nu}\theta\right)$ on the RHS of all equations are treated implicitly via the relaxation equations associated with (57). Although such a treatment seems approximate, compared to the flux-divergence operator on the left-hand side (LHS), it is important to point out that we are primarily interested in near-equilibrium behavior where $\Pi$, $q^{\nu}$, and $\pi^{\mu\nu}$ are small. As such, these specific second-order terms constitute only a minor correction and could even be omitted. This is also consistent with neglecting all other second-order transport terms in these equations, which should otherwise be present Denicol et al. (2012). On the other hand, since first-order gradient terms are important for the evolution, they require more sophisticated numerical methods, as outlined in Sec. IV . Second, we point out that in deriving these equations we have made use of our freedom to remove all terms proportional to $u^{\mu}$ on the RHS. Due to the projected nature of the Israel-Stewart limit, these terms would get removed when projecting the equations into the fluid frame (in the continuum limit). In the stiff-relaxation approach, all these terms are hence implicitly accounted for in the relaxation operator scaling with the fluid four-velocity $u^{\mu}$.

It would be very interesting to investigate if our formulation of the equations of motion can be proven to be hyperbolic in the full nonlinear limit. This would be especially important when coupling our system to Einstein’s equations, which is needed to investigate viscous phenomena in neutron star mergers Shibata et al. (2017); Radice (2017); Alford et al. (2018); Most et al. (2021) or black hole accretion in non-dynamical spacetimes Foucart et al. (2016, 2017). To the best of our knowledge, there is currently no formulation of general-relativistic dissipative magnetohydrodynamics that is proven to be strongly hyperbolic in the nonlinear regime. For instance, such a statement is not known to hold for the elegant approach proposed in Ref. Chandra et al. (2015). Drawing from our experience with the zero magnetic field limit of Israel-Stewart-like equations Bemfica et al. (2019a, 2021), such a nonlinear analysis of hyperbolicity and causality would give rise to important constraints on the values of the dissipative currents, the transport coefficients, and in this case also the magnetic field.

The equations presented in this section are written in a way that seems to be very natural for for performing such a nonlinear investigation. Indeed, here we have taken the first step towards this result as we have been able to rewrite the system as a set of first-order PDEs in flux-conservative form, which can be beneficial in such studies. We remark, however, that proving that our nonlinear system of equations is strongly hyperbolic is a very challenging mathematical task, which we hope to address in the near future. On the other hand, a linearized analysis of our equations can be performed following Ref. Biswas et al. (2020). The only difference would be that the inclusion of the stiff advection equations introduces additional non-hydrodynamic modes to the system parametrized by the new relaxation times $\omega_{\zeta}$, $\omega_{\kappa}$, $\omega_{\eta}$ and $\omega_{\theta}$. Therefore, in the light of Biswas et al. (2020), we leave the investigation of the linear regime of our equations to a future dedicated study.