General-relativistic hydrodynamics of non-perfect fluids: 3+1 conservative formulation and application to viscous black-hole accretion
MMNRAS , 1–55 (2021) Preprint February 23, 2021 Compiled using MNRAS L A TEX style file v3.0
General-relativistic hydrodynamics of non-perfect fluids: 3+1conservative formulation and application to viscous black-holeaccretion
Michail Chabanov , Luciano Rezzolla , , and Dirk H. Rischke Institut für Theoretische Physik, Goethe-Universität, Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany Frankfurt Institute for Advanced Studies, Ruth-Moufang-Str. 1, 60438 Frankfurt am Main, Germany School of Mathematics, Trinity College, Dublin 2, Ireland
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We consider the relativistic hydrodynamics of non-perfect fluids with the goal of de-termining a formulation that is suited for numerical integration in special-relativisticand general-relativistic scenarios. To this end, we review the various formulations ofrelativistic second-order dissipative hydrodynamics proposed so far and present indetail a particular formulation that is fully general, causal, and can be cast into a 3+1flux-conservative form as the one employed in modern numerical-relativity codes.As an example, we employ a variant of this formulation restricted to a relaxation-type equation for the bulk viscosity in the general-relativistic magnetohydrodynam-ics code
BHAC . After adopting the formulation for a series of standard and non-standard tests in 1+1-dimensional special-relativistic hydrodynamics, we consider anovel general-relativistic scenario, namely, the stationary, spherically symmetric vis-cous accretion onto a black hole. The newly developed solution – which can exhibiteven considerable deviations from the inviscid counterpart – can be used as a testbedfor numerical codes simulating non-perfect fluids on curved backgrounds.
Key words: hydrodynamics, shock waves, accretion
The detection of the first binary neutron-star (BNS) merger event, GW170817 (Abbott et al. 2017a)has provided a new valuable tool to study matter and gravity under extreme conditions. Espe-cially the detection of electromagnetic counterparts in form of a short gamma-ray burst (Abbott © 2021 The Authors a r X i v : . [ g r- q c ] F e b M. Chabanov, L. Rezzolla and D.H. Rischke et al. 2017b) and of a kilonova (Drout et al. 2017; Cowperthwaite et al. 2017) accompanyingGW170817, has made this event an incredibly rich laboratory for physics, providing a number ofconstraints on the equation of state (EOS) of nuclear matter (see, e.g., Margalit & Metzger 2017;Bauswein et al. 2017; Rezzolla et al. 2018; Ruiz et al. 2018; Annala et al. 2018; Radice et al.2018b; Most et al. 2018; De et al. 2018; Abbott et al. 2018; Montaña et al. 2019; Raithel et al.2018; Tews et al. 2018; Malik et al. 2018; Koeppel et al. 2019; Shibata et al. 2019).BNS mergers are highly dynamical and nonlinear phenomena especially during the first fewmilliseconds after merger, (see, e.g., Baiotti & Rezzolla 2017; Paschalidis 2017; Burns 2020, forsome reviews). As has been shown recently, harmonic density oscillations in this violent phasecould be damped significantly due to bulk-viscosity dissipation coming from modified Urca pro-cesses (Alford et al. 2018, 2020; Alford & Haber 2020). Thus, bulk viscosity might lead to modifi-cations in the post-merger gravitational-wave signal. Furthermore, it was shown in high-resolutiongeneral-relativistic magnetohydrodynamic (GRMHD) calculations that the matter after merger isunstable to MHD instabilities, e.g., the Kelvin-Helmholtz instability (Baiotti et al. 2008; Radice &Rezzolla 2012) and the magnetorotational instability (Siegel et al. 2013; Kiuchi et al. 2018). Dueto these instabilities, turbulence can develop and be maintained, which will ultimately influenceBNS-merger observables such as the gravitational-wave signal and the ejected matter.The emergence of turbulence clearly represents a challenge for numerical-relativity simula-tions, which can only be performed with limited resolutions and are normally carried out withresolutions that are well above the scale at which turbulence is physically quenched off. As a re-sult, a number of studies have employed effective shear-viscous models in order to capture theeffects of magneto-turbulent motion in the remnant of BNS mergers. An approach to handle thisproblem has been suggested with the use of general-relativistic large-eddy simulations (LES) inpure hydrodynamics. This method maps effects from numerical calculations with high resolutionto simulations with lower resolution, which would otherwise disappear in an implicit filteringprocedure (Radice 2017; Radice et al. 2018a,c; Radice 2020). The same method has been system-atically extended to the full set GRMHD equations to study the amplification of the magnetic fieldshortly after merger (Viganò et al. 2020; Aguilera-Miret et al. 2020). In this way, it was shown thatby considering the subgrid-scale model for the induction equation only, and suitably increasing(well above the expected value of order one) the phenomenological coefficients associated withthe LES terms, magnetic-field strengths comparable to an high-resolution direct simulation can beachieved.While the robustness and universality – i.e., their effectiveness to capture multi-scale turbu-
MNRAS , 1–55 (2021)
RHD of non-perfect fluids
HL83 ), to serve as our reference
MNRAS000
MNRAS000 , 1–55 (2021)
M. Chabanov, L. Rezzolla and D.H. Rischke second-order theory. We discuss the properties of this formulation and provide a comprehensivecomparison with other and equivalent formulations of relativistic dissipative hydrodynamics thathave been proposed and employed in the literature. In addition, we discuss the implementation ofthis formulation in the general-relativistic magnetohydrodynamics code
BHAC , where it is subjectto a number of tests in special and general relativity.The structure of the paper is as follows. In Section 2 we first provide a brief review of rela-tivistic dissipative hydrodynamics mostly to recall definitions and conventions, while Section 3 isdedicated to the presentation of the set of the
HL83 equations in a general four-dimensional man-ifold. Extending the work of Peitz & Appl (1997) and Peitz & Appl (1999), in Section 4 we derivea 3+1 decomposed version of
HL83 , including all source terms and gradients of fluid variables,which are separately listed in Section 4.3, such that the coupling to a numerical-relativity codebecomes possible. In Section 5 we use this formulation to implement a simplified version for bulkviscosity in the GRMHD code
BHAC (Porth et al. 2017) and perform tests in special relativity inSection 6. As a first but rigorous test in general relativity, we consider in Section 7 the problemof stationary, spherically symmetric viscous accretion onto a black hole. We conclude this workin Section 8, where we also comment on future developments. A number of Appendices is alsoused to provide additional information on the currently available second-order descriptions of rel-ativistic dissipative effects and the differences among them. In addition, details are provided onthe numerical solution of the viscous accretion problem, which is less trivial than may appear atfirst sight.Unless stated otherwise, we use geometrised units, where the speed of light c = 1 and thegravitational constant G = 1 . Greek letters denote spacetime indices, i.e., µ = 0 , , , , whileRoman letters cover spatial indices only, i.e., i = 1 , , . Also, we make use of Einstein’s summa-tion convention and choose the metric signature to be ( − , + , + , +) . Bold symbols such as g referto tensors of generic rank while the symbol ∇ denotes the covariant derivative with respect to g and we use the following definitions for the components of a symmetric and antisymmetric rank-2tensor: T ( µν ) := ( T µν + T νµ ) , T [ µν ] := ( T µν − T νµ ) . A large portion of the theory of relativistic dissipative hydrodynamics dates back to the 70’s and isboth vast and somewhat confusing, with multiple formulations having been produced, often with-out a clear demarcation or obvious guidance on which formulation may be seen as the most robust.
MNRAS , 1–55 (2021)
RHD of non-perfect fluids “First-order” theories are relativistic extensions of the Navier-Stokes(NS) equations, in which the dissipative quantities are directly related to first-order gradients ofthe primary fluid variables appearing in perfect-fluid hydrodynamics, such as e.g., pressure, rest-mass density and fluid velocity . However, already early on (Hiscock & Lindblom 1985; Hiscock& Lindblom 1987), it was shown that a large class of first-order theories suffer from instabilityas well as acausal behaviour due to their parabolic nature (Hiscock & Lindblom 1985; Hiscock &Lindblom 1987; Denicol et al. 2008a). Because of this drawback, “second-order” theories weredeveloped, which extend first-order theories by including terms of second order in quantities whichdescribe deviations from perfect-fluid hydrodynamics. On the one hand, these are dissipative cur-rents like bulk-viscosity pressure, heat current (also referred to as heat flux) or shear-stress tensor.On the other hand, these are gradients of the primary fluid variables, which are in principle inde-pendent from the dissipative currents.While this can be done in different ways, here we focus on three conceptually rather differ-ent approaches. The first approach was suggested originally by Israel (1976) and it extends thedefinition of the entropy current by all terms of second order in the dissipative currents allowedby symmetry. The second-law of thermodynamics then leads to relaxation-type equations for thedissipative currents (see also Hiscock & Lindblom 1983; Muronga 2004; Jaiswal et al. 2013).The second approach, instead, adds to the definitions of the dissipative currents all terms ofsecond order in gradients of the primary fluid variables that are allowed by symmetry. This ap-proach follows in spirit that of a systematic gradient expansion and was first suggested by Baieret al. (2008). Unfortunately, second-order (spatial) gradients usually destroy the hyperbolic natureof a partial differential equation. To counter this problem and obtain a hyperbolic formulation –that can therefore be solved numerically – an amendment to the gradient expansion was proposedby Baier et al. (2008), where gradients of primary fluid variables were replaced by dissipativecurrents, using the first-order Navier-Stokes relation. This effectively leads to a resummation ofgradients and, ultimately, to relaxation-type equations for the dissipative currents. The large majority of these works concentrate on single-fluid scenarios, but valuable work has been done also when considering dissipativehydrodynamics of multi-fluids [see, e.g., Carter 1989; Andersson & Comer 2015, but also Andersson & Comer 2020 for a recent review].MNRAS000
RHD of non-perfect fluids “First-order” theories are relativistic extensions of the Navier-Stokes(NS) equations, in which the dissipative quantities are directly related to first-order gradients ofthe primary fluid variables appearing in perfect-fluid hydrodynamics, such as e.g., pressure, rest-mass density and fluid velocity . However, already early on (Hiscock & Lindblom 1985; Hiscock& Lindblom 1987), it was shown that a large class of first-order theories suffer from instabilityas well as acausal behaviour due to their parabolic nature (Hiscock & Lindblom 1985; Hiscock &Lindblom 1987; Denicol et al. 2008a). Because of this drawback, “second-order” theories weredeveloped, which extend first-order theories by including terms of second order in quantities whichdescribe deviations from perfect-fluid hydrodynamics. On the one hand, these are dissipative cur-rents like bulk-viscosity pressure, heat current (also referred to as heat flux) or shear-stress tensor.On the other hand, these are gradients of the primary fluid variables, which are in principle inde-pendent from the dissipative currents.While this can be done in different ways, here we focus on three conceptually rather differ-ent approaches. The first approach was suggested originally by Israel (1976) and it extends thedefinition of the entropy current by all terms of second order in the dissipative currents allowedby symmetry. The second-law of thermodynamics then leads to relaxation-type equations for thedissipative currents (see also Hiscock & Lindblom 1983; Muronga 2004; Jaiswal et al. 2013).The second approach, instead, adds to the definitions of the dissipative currents all terms ofsecond order in gradients of the primary fluid variables that are allowed by symmetry. This ap-proach follows in spirit that of a systematic gradient expansion and was first suggested by Baieret al. (2008). Unfortunately, second-order (spatial) gradients usually destroy the hyperbolic natureof a partial differential equation. To counter this problem and obtain a hyperbolic formulation –that can therefore be solved numerically – an amendment to the gradient expansion was proposedby Baier et al. (2008), where gradients of primary fluid variables were replaced by dissipativecurrents, using the first-order Navier-Stokes relation. This effectively leads to a resummation ofgradients and, ultimately, to relaxation-type equations for the dissipative currents. The large majority of these works concentrate on single-fluid scenarios, but valuable work has been done also when considering dissipativehydrodynamics of multi-fluids [see, e.g., Carter 1989; Andersson & Comer 2015, but also Andersson & Comer 2020 for a recent review].MNRAS000 , 1–55 (2021)
M. Chabanov, L. Rezzolla and D.H. Rischke
Finally, the third approach considers hydrodynamics as an effective theory in the long-wavelength,low-frequency limit of an underlying microscopic theory. Taking kinetic theory for the latter, it ispossible to derive hydrodynamics from the relativistic Boltzmann equation and the method ofmoments (Israel & Stewart 1979; Betz et al. 2009; Denicol et al. 2012a,b; Jaiswal 2013). Thisapproach leads to an infinite system of coupled equations of motion for the moments of the devia-tion of the single-particle distribution function from local equilibrium. In particular, Denicol et al.(2012b) have proposed a systematic power-counting scheme in Knudsen and inverse Reynoldsnumbers that allows for a truncation of the infinite set of moment equations at an arbitrary orderin these quantities (Denicol et al. 2014).The stability and causality of second-order theories has been studied in a series of works.In particular, stability and causality were analysed in a linear regime using the rest frame of thefluid (Hiscock & Lindblom 1983), with the extension to a moving frame in the case of either bulkviscosity (Denicol et al. 2008a), or of shear viscosity (Pu et al. 2010). In this way, it was found thatstability implies causality and vice-versa, but only if the relaxation times are larger than a certaintimescale. These studies were recently extended to include heat flow (Brito & Denicol 2020) anda magnetic field (Biswas et al. 2020).While most of these works have considered either a flat Minkowski background or a fixedcurved background, recent works have studied the coupling of Einstein equations to second-orderdissipative hydrodynamics. This has been done by Bemfica et al. (2019b) when considering onlythe effects of bulk viscosity and by Bemfica et al. (2020) when including also shear viscosity. Thesesecond-order formulations – that are shown to be causal and to admit unique solutions under ratherreasonable conditions – effectively pave the way for applications in numerical relativity and hencefor the modelling of BNS mergers.We conclude this overview Section by remarking that it has become clear recently that theunstable and acausal behaviour of first-order theories discussed above is actually related to theparticular choice of rest frame of the fluid, for instance the particle frame proposed by Eckart(1940) or the energy-density frame proposed by Landau & Lifshitz (2004). As a result, new causaland stable formulations of first-order theories have been found by choosing different rest frames(Van & Biro 2012; Disconzi et al. 2017; Bemfica et al. 2019a; Kovtun 2019; Hoult & Kovtun2020; Taghinavaz 2020). The prospects of using these theories, that are in principle easier to solvenumerically, are very good but no concrete application to HICs or BNS mergers has been presentedyet.
MNRAS , 1–55 (2021)
RHD of non-perfect fluids In this Section, we briefly review the phenomenological approach to second-order GRDHD lead-ing to the
HL83 formulation in a generic four-dimensional manifold. The first assumption madeis that the rest-mass current J µ and the energy-momentum tensor T µν continue to provide avalid description for fluids which are out of thermodynamical equilibrium. Fluids in such an off-equilibrium state show dissipative effects, which characterize them as non-perfect fluids. Hereafter,we will distinguish perfect fluids from generic (non-perfect) fluids, by using the lower index “PF”.The conservation of rest mass (or some associated conserved number density), energy, andmomentum, is expressed by the five conservation equations ∇ µ J µ = 0 , (1) ∇ µ T µν = 0 , (2)where J µ = J µ PF and T µν = T µν PF . These equations can be complemented by one equation of state(EOS) of the form e PF = e PF ( ρ PF , p PF ) , (3)so that perfect fluids need only five independent variables to be described. However, under moregeneral conditions, i.e., for non-perfect fluids, J µ and T µν = T νµ contain 14 independent vari-ables. The physical meaning of the nine additional degrees of freedom for non-perfect fluids canbe made transparent in a tensor decomposition with respect to the fluid velocity u µ . Choosing theEckart frame (Eckart 1940; Rezzolla & Zanotti 2013), where u µ is the velocity of the flow of restmass, the energy-momentum tensor reads T µν = eu µ u ν + ( p + Π) h µν + 2 q ( µ u ν ) + π µν , (4)while J µ = ρu µ maintains the form of the rest-mass current of a perfect fluid, where ρ is therest-mass density.In Eq. (4), e and p + Π are the energy density and the total isotropic pressure, Π is the bulk-viscosity pressure and h µν := g µν + u µ u ν is the projector orthogonal to u µ , respectively. Further-more, q µ is the heat current, which is orthogonal to the fluid four-velocity, i.e., q µ u µ = 0 . Finally, π µν is the shear-stress tensor with the following properties: it is symmetric π µν = π νµ , purelyspatial π µν u µ = 0 , and trace-free π µµ = 0 . Hereafter we will always refer to rest mass as this is a quantity normally conserved in simulations of neutron stars. However, the rest mass herecan be replaced by any other conserved charge, e.g., baryon number.MNRAS000
HL83 formulation in a generic four-dimensional manifold. The first assumption madeis that the rest-mass current J µ and the energy-momentum tensor T µν continue to provide avalid description for fluids which are out of thermodynamical equilibrium. Fluids in such an off-equilibrium state show dissipative effects, which characterize them as non-perfect fluids. Hereafter,we will distinguish perfect fluids from generic (non-perfect) fluids, by using the lower index “PF”.The conservation of rest mass (or some associated conserved number density), energy, andmomentum, is expressed by the five conservation equations ∇ µ J µ = 0 , (1) ∇ µ T µν = 0 , (2)where J µ = J µ PF and T µν = T µν PF . These equations can be complemented by one equation of state(EOS) of the form e PF = e PF ( ρ PF , p PF ) , (3)so that perfect fluids need only five independent variables to be described. However, under moregeneral conditions, i.e., for non-perfect fluids, J µ and T µν = T νµ contain 14 independent vari-ables. The physical meaning of the nine additional degrees of freedom for non-perfect fluids canbe made transparent in a tensor decomposition with respect to the fluid velocity u µ . Choosing theEckart frame (Eckart 1940; Rezzolla & Zanotti 2013), where u µ is the velocity of the flow of restmass, the energy-momentum tensor reads T µν = eu µ u ν + ( p + Π) h µν + 2 q ( µ u ν ) + π µν , (4)while J µ = ρu µ maintains the form of the rest-mass current of a perfect fluid, where ρ is therest-mass density.In Eq. (4), e and p + Π are the energy density and the total isotropic pressure, Π is the bulk-viscosity pressure and h µν := g µν + u µ u ν is the projector orthogonal to u µ , respectively. Further-more, q µ is the heat current, which is orthogonal to the fluid four-velocity, i.e., q µ u µ = 0 . Finally, π µν is the shear-stress tensor with the following properties: it is symmetric π µν = π νµ , purelyspatial π µν u µ = 0 , and trace-free π µµ = 0 . Hereafter we will always refer to rest mass as this is a quantity normally conserved in simulations of neutron stars. However, the rest mass herecan be replaced by any other conserved charge, e.g., baryon number.MNRAS000 , 1–55 (2021)
M. Chabanov, L. Rezzolla and D.H. Rischke
As a result of the additional degrees of freedom, the bulk-viscosity pressure adds one, theheat current three and the shear-stress tensor five independent variables to the non-perfect energy-momentum tensor. Note that one can define the temperature T and the baryon chemical potential µ by matching the energy density e and the number density n to the corresponding values of afictitious equilibrium state, i.e., e = e PF , n = n PF , such that p = p PF , and then from well-knownthermodynamical relations also T and µ , can be determined from the EOS (3): T = (cid:18) ∂s∂e (cid:19) n , µ = − T (cid:18) ∂s∂n (cid:19) e , s = e + p − µnT , (5)where s denotes the entropy density. We also remark that in the Landau frame (Landau & Lifshitz2004; Rezzolla & Zanotti 2013), where the fluid four-velocity is the timelike eigenvector of theenergy-momentum tensor, the heat current is absent in the tensor decomposition (4), but its threeindependent degrees of freedom reappear in the form of a diffusion current as part of the non-perfect rest-mass current.Hereafter, we will make use of the Eckart frame because it is more intuitive to relate thefluid four-velocity to the motion of particles in astrophysical scenarios. In contrast to that, oneoften chooses the Landau frame in the context of heavy-ion collisions, because baryon chemicalpotentials are typically small in the center of the collision zone.To close the system (1)–(2), we need therefore nine additional equations which determinethe evolution of the dissipative currents Π , q µ and π µν . Following Israel (1976) and Hiscock &Lindblom (1983), such relations are obtained by ensuring the positivity of entropy production.The entropy current depends quadratically on the dissipative currents and reads: S µ = su µ + q µ T − (cid:0) β Π + β q α q α + β π αβ π αβ (cid:1) u µ T + α Π q µ T + α q α π αµ T , (6)where the physical meaning of the coefficients α , α , β , β , β will become clear below. Fromthe second law of thermodynamics, i.e. ∇ µ S µ ≥ , (7) MNRAS , 1–55 (2021)
RHD of non-perfect fluids τ Π ˙Π = Π NS − Π − ζ Π T ∇ µ (cid:18) τ Π u µ ζT (cid:19) + α ζ ∇ µ q µ + γ ζT q µ ∇ µ (cid:16) α T (cid:17) , (8) τ q ˙ q (cid:104) µ (cid:105) = q µ NS − q µ − κT q µ ∇ ν (cid:18) τ q u ν κT (cid:19) + κT (cid:20) α ∇ (cid:104) µ (cid:105) Π + α ∇ ν π ν (cid:104) µ (cid:105) + (1 − γ )Π T ∇ (cid:104) µ (cid:105) (cid:16) α T (cid:17) + (1 − γ ) T π µν ∇ ν (cid:16) α T (cid:17) (cid:21) , (9) τ π ˙ π (cid:104) µν (cid:105) = π µν NS − π µν − ηT π µν ∇ λ (cid:18) τ π u λ ηT (cid:19) + 2 α η ∇ (cid:104) µ q ν (cid:105) + 2 γ ηT q (cid:104) µ ∇ ν (cid:105) (cid:16) α T (cid:17) , (10)where we have introduced two new coefficients, γ and γ , whose existence is due to the ambiguitywhen factoring out the terms which involve the products Π q µ and q α π αµ . Also, note in the expres-sions above the introduction of an operator we will often employ, i.e., the “comoving derivative” ˙ A := ( u · ∇ ) A = u µ ∇ µ A , where A can be an arbitrary tensor field. The comoving derivative isnaturally accompanied by “relaxation times” for the bulk-viscosity pressure, the heat current, andshear-stress tensor, that are respectively defined as τ Π := β ζ , τ q := β κT , τ π := 2 β η . (11)where ζ is the bulk viscosity, κ the heat conductivity, and η the shear viscosity, respectively. All ofthese coefficients are by definition non-negative and they essentially set the timescales over whichnon-equilibrium effects act to push the equations of GRDHD towards the NS solution. Clearly, aninviscid fluid has zero relaxation times.A physically intuitive interpretation of the bulk-viscosity pressure, of the heat current, andof the shear-stress tensor comes from the naive extensions of the non-relativistic NS equations(Landau & Lifshitz 2004), which leads to the NS expressions of these quantities as Π NS = − ζ Θ , (12) q µ NS = − κT (cid:0) ∇ (cid:104) µ (cid:105) ln T + a µ (cid:1) , (13) π µν NS = − ησ µν , (14)Furthermore, we define b (cid:104) µ (cid:105) := h µν b ν the projection of an arbitrary vector b µ onto the directionorthogonal to u µ , while b (cid:104) µν (cid:105) := ( h α ( µ h ν ) β − h µν h αβ ) b αβ denotes the symmetric and trace-free Note that after using the conservation equations (2) for perfect fluids, which read mnha µ = −∇ (cid:104) µ (cid:105) p – where h is the specific enthalpy and m the rest mass of the particles constituting the fluid – and using the Gibbs-Duhem relation, i.e., dp = sdT + ndµ , the NS value of the heatcurrent (13) can also be approximated by the expression q µ NS = ( κT /mh ) ∇ (cid:104) µ (cid:105) ( µ/T ) . This expression is an identity at first order in gradientsof primary fluid variables and in dissipative currents in the sense that if higher-order terms would be used in the form of the acceleration, then termsof order O K or O R or O RK would appear [see Eq. (A4) in Appendix A for the definitions of the symbols O K , O R and O RK ].MNRAS000
RHD of non-perfect fluids τ Π ˙Π = Π NS − Π − ζ Π T ∇ µ (cid:18) τ Π u µ ζT (cid:19) + α ζ ∇ µ q µ + γ ζT q µ ∇ µ (cid:16) α T (cid:17) , (8) τ q ˙ q (cid:104) µ (cid:105) = q µ NS − q µ − κT q µ ∇ ν (cid:18) τ q u ν κT (cid:19) + κT (cid:20) α ∇ (cid:104) µ (cid:105) Π + α ∇ ν π ν (cid:104) µ (cid:105) + (1 − γ )Π T ∇ (cid:104) µ (cid:105) (cid:16) α T (cid:17) + (1 − γ ) T π µν ∇ ν (cid:16) α T (cid:17) (cid:21) , (9) τ π ˙ π (cid:104) µν (cid:105) = π µν NS − π µν − ηT π µν ∇ λ (cid:18) τ π u λ ηT (cid:19) + 2 α η ∇ (cid:104) µ q ν (cid:105) + 2 γ ηT q (cid:104) µ ∇ ν (cid:105) (cid:16) α T (cid:17) , (10)where we have introduced two new coefficients, γ and γ , whose existence is due to the ambiguitywhen factoring out the terms which involve the products Π q µ and q α π αµ . Also, note in the expres-sions above the introduction of an operator we will often employ, i.e., the “comoving derivative” ˙ A := ( u · ∇ ) A = u µ ∇ µ A , where A can be an arbitrary tensor field. The comoving derivative isnaturally accompanied by “relaxation times” for the bulk-viscosity pressure, the heat current, andshear-stress tensor, that are respectively defined as τ Π := β ζ , τ q := β κT , τ π := 2 β η . (11)where ζ is the bulk viscosity, κ the heat conductivity, and η the shear viscosity, respectively. All ofthese coefficients are by definition non-negative and they essentially set the timescales over whichnon-equilibrium effects act to push the equations of GRDHD towards the NS solution. Clearly, aninviscid fluid has zero relaxation times.A physically intuitive interpretation of the bulk-viscosity pressure, of the heat current, andof the shear-stress tensor comes from the naive extensions of the non-relativistic NS equations(Landau & Lifshitz 2004), which leads to the NS expressions of these quantities as Π NS = − ζ Θ , (12) q µ NS = − κT (cid:0) ∇ (cid:104) µ (cid:105) ln T + a µ (cid:1) , (13) π µν NS = − ησ µν , (14)Furthermore, we define b (cid:104) µ (cid:105) := h µν b ν the projection of an arbitrary vector b µ onto the directionorthogonal to u µ , while b (cid:104) µν (cid:105) := ( h α ( µ h ν ) β − h µν h αβ ) b αβ denotes the symmetric and trace-free Note that after using the conservation equations (2) for perfect fluids, which read mnha µ = −∇ (cid:104) µ (cid:105) p – where h is the specific enthalpy and m the rest mass of the particles constituting the fluid – and using the Gibbs-Duhem relation, i.e., dp = sdT + ndµ , the NS value of the heatcurrent (13) can also be approximated by the expression q µ NS = ( κT /mh ) ∇ (cid:104) µ (cid:105) ( µ/T ) . This expression is an identity at first order in gradientsof primary fluid variables and in dissipative currents in the sense that if higher-order terms would be used in the form of the acceleration, then termsof order O K or O R or O RK would appear [see Eq. (A4) in Appendix A for the definitions of the symbols O K , O R and O RK ].MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke projection of an arbitrary tensor b µν orthogonal u µ . It is then easy to show that when the entropycurrent contains only first-order dissipative currents, i.e., β = β = β = α = α = 0 , then(Rezzolla & Zanotti 2013): Π = Π NS , (15) q µ = q µ NS , (16) π µν = π µν NS , (17)As anticipated in Sec. 2, Eqs. (15)–(17) are acausal and fall into the class of equations investigatedby Hiscock & Lindblom (1985) that were shown to be unstable around hydrostatic equilibriumstates. We next derive a general-relativistic 3+1 flux-conservative (or simply “conservative”) formulationof the dissipative-hydrodynamics Eqs. (8) – (10) so as to provide a comprehensive and completeway of including causal dissipative effects in general-relativistic calculations.
As customary in the so-called 3+1 decomposition of spacetime, we decompose the four-dimensionalmanifold such that notions of time and space reappear. This can be achieved by foliating spacetimein terms of a set of non-intersecting spacelike hypersurfaces. On each hypersurface a unit timelikefour-vector field n can be defined such that n is normal to its corresponding hypersurface. As wehave n µ n µ = − , its trajectory can also be interpreted as a reference observer, called “normal” or“Eulerian” observer. In this way, it is possible represent physical laws in terms of projections ei-ther parallel or orthogonal to n . For a detailed description of this formalism we refer to Alcubierre(2008), Gourgoulhon (2012), and Rezzolla & Zanotti (2013).The generic line-element in a 3+1 decomposition can be written as ds = − ( α − β i β i ) dt + 2 β i dx i dt + γ ij dx i dx j , (18)where α is the so-called lapse function, the purely spatial vector β , i.e., β µ = (0 , β i ) T , is the shiftvector and γ ij denotes the components of the purely spatial metric γ , defined as γ µν := g µν + n µ n ν , (19)which acts as the projection operator onto spatial hypersurfaces. Note that the following identities MNRAS , 1–55 (2021)
RHD of non-perfect fluids β i := γ ij β j , √− g = α √ γ , where g := det ( g µν ) and γ := det ( γ ij ) . The componentsof n and its dual are given by: n µ = 1 α (1 , − β i ) T , n µ = ( − α, , , . (20)Having a timelike unit four-vector n , we can use it to decompose any tensor into a part that isparallel ( (cid:107) ) and perpendicular ( ⊥ ) to it. We start from the four-velocity u µ , which can then bewritten as: u µ = u µ (cid:107) + u µ ⊥ := ( − n ν u ν ) n µ + γ µν u ν = W ( n µ + v µ ) , (21)where W := (1 − v i v i ) − / is the Lorentz factor and v µ := (0 , v i ) T the fluid velocity measuredby the normal observer. Analogously, we can decompose the heat current q µ and the shear-stresstensor π µν into components relative to n : q µ = q µ (cid:107) + q µ ⊥ := ( − n ν q ν ) n µ + γ µν q ν , (22) π µν = π µν (cid:107) + π × µν + π νµ × + π µν ⊥ : = (cid:0) n α n β π αβ (cid:1) n µ n ν + (cid:0) − n α γ ν β π αβ (cid:1) n µ + (cid:0) − n β γ µα π αβ (cid:1) n ν + γ µα γ ν β π αβ . (23)Additionally, since the time components of the parallel projections are effectively scalar functions,we treat them as such after the following definitions ˚ q : = − n µ q µ , (24) ˚ π : = n µ n ν π µν , (25) ˚ π λ : = − n µ γ λν π µν . (26)Note that the parallel and perpendicular components are not independent and indeed the followingidentities can be obtained and will be used hereafter ˚ q = v i q i ⊥ , ˚ π i = v j π ij ⊥ , ˚ π = v i ˚ π i = π i ⊥ i . (27)Note that expressions (27) imply that the knowledge of q i ⊥ and π ij ⊥ suffices to fully reconstruct q µ and π µν if the fluid three-velocity v i and the full metric g µν is known.The projections (27) turn out to be helpful in the 3+1 decomposition of the non-perfect energy-momentum tensor (4) T µν = En µ n ν + S µ n ν + S ν n µ + S µν , (28) MNRAS000
RHD of non-perfect fluids β i := γ ij β j , √− g = α √ γ , where g := det ( g µν ) and γ := det ( γ ij ) . The componentsof n and its dual are given by: n µ = 1 α (1 , − β i ) T , n µ = ( − α, , , . (20)Having a timelike unit four-vector n , we can use it to decompose any tensor into a part that isparallel ( (cid:107) ) and perpendicular ( ⊥ ) to it. We start from the four-velocity u µ , which can then bewritten as: u µ = u µ (cid:107) + u µ ⊥ := ( − n ν u ν ) n µ + γ µν u ν = W ( n µ + v µ ) , (21)where W := (1 − v i v i ) − / is the Lorentz factor and v µ := (0 , v i ) T the fluid velocity measuredby the normal observer. Analogously, we can decompose the heat current q µ and the shear-stresstensor π µν into components relative to n : q µ = q µ (cid:107) + q µ ⊥ := ( − n ν q ν ) n µ + γ µν q ν , (22) π µν = π µν (cid:107) + π × µν + π νµ × + π µν ⊥ : = (cid:0) n α n β π αβ (cid:1) n µ n ν + (cid:0) − n α γ ν β π αβ (cid:1) n µ + (cid:0) − n β γ µα π αβ (cid:1) n ν + γ µα γ ν β π αβ . (23)Additionally, since the time components of the parallel projections are effectively scalar functions,we treat them as such after the following definitions ˚ q : = − n µ q µ , (24) ˚ π : = n µ n ν π µν , (25) ˚ π λ : = − n µ γ λν π µν . (26)Note that the parallel and perpendicular components are not independent and indeed the followingidentities can be obtained and will be used hereafter ˚ q = v i q i ⊥ , ˚ π i = v j π ij ⊥ , ˚ π = v i ˚ π i = π i ⊥ i . (27)Note that expressions (27) imply that the knowledge of q i ⊥ and π ij ⊥ suffices to fully reconstruct q µ and π µν if the fluid three-velocity v i and the full metric g µν is known.The projections (27) turn out to be helpful in the 3+1 decomposition of the non-perfect energy-momentum tensor (4) T µν = En µ n ν + S µ n ν + S ν n µ + S µν , (28) MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke where E := n µ n ν T µν = ( e + p + Π) W + 2˚ qW − ( p + Π − ˚ π ) , (29) S µ := − n ν γ µα T αν = ( e + p + Π) W v µ + W (cid:16) ˚ qv µ + q µ ⊥ (cid:17) + ˚ π µ , (30) S µν := γ µα γ ν β T αβ = ( e + p + Π) W v µ v ν + W (cid:16) q µ ⊥ v ν + q ν ⊥ v µ (cid:17) + ( p + Π) γ µν + π µν ⊥ , (31)are: the energy density, the momentum density, and the purely spatial part of the energy-momentumtensor, respectively. Additionally, we will make use of the acceleration of the normal observer andthe extrinsic curvature given by ˆ a µ := n ν ∇ ν n µ = γ µν ∂ ν ln α , (32) K µν := − L n γ µν = − γ µλ ∇ λ n ν = −∇ µ n ν − n µ ˆ a ν , (33)respectively, where L n denotes the Lie derivative along n . Finally, we recall two important four-dimensional tensor identities that will be useful later on to to obtain a flux-conservative formulationfor the conservation of rest mass (1) and energy-momentum (2), namely, ∇ µ J µ = ∂ µ ( √− gJ µ ) √− g , (34) ∇ µ T µν = g νλ (cid:20) ∂ µ ( √− gT µλ ) √− g − T αβ ∂ λ g αβ (cid:21) . (35) We next rewrite the full set of general-relativistic dissipative hydrodynamics, i.e., Eqs. (1), (2),and (8) – (10), and which essentially represent the equations of the
HL83 formulation, in a 3+1flux-conservative form. We start by recalling that a system of partial differential equations is saidto be flux-conservative (or simply “conservative”) if it can be written as ∂ t U + ∂ i F i ( U ) = S , (36)where U is the “state vector” , F i are the “flux vectors” and S is the “source vector” . Employingnow Eq. (34) it is possible to rewrite Eq. (1) as ∂ t ( √ γD ) + ∂ i (cid:2) √ γD (cid:0) αv i − β i (cid:1)(cid:3) = 0 , (37)where D := ραu t = ρW is the conserved rest mass. Similarly, we use equations (28) and (35) toobtain a flux-conservative form of the equations for the conservation of energy and momentum (2) ∂ t ( √ γS j ) + ∂ i (cid:2) √ γ (cid:0) αS ij − β i S j (cid:1)(cid:3) = √ γ (cid:18) αS ik ∂ j γ ik + S i ∂ j β i − E∂ j α (cid:19) , (38) ∂ t [ √ γ ( E − D )] + ∂ i (cid:8) √ γ (cid:2) α (cid:0) S i − v i D (cid:1) − β i ( E − D ) (cid:3)(cid:9) = √ γ (cid:0) αS ij K ij − S j ∂ j α (cid:1) . (39) MNRAS , 1–55 (2021)
RHD of non-perfect fluids E − D rather than simply E in equation (39) because the conservation of E − D is more accurate than that of E only. Comparing Eqs. (36)–(39) we can now read off thecomponents of U , F i and S corresponding to the conservation equations. We will write themdown explicitly in the next Section, after we have reformulated the constitutive equations (8)–(10)in conservative form. In order to do this, we consider the evolution equation for the heat current(9) as an example and extend the treatment to the other equations afterwards. We start with theterm ˙ q (cid:104) µ (cid:105) = h µν u λ ∇ λ q ν = u λ ∇ λ q µ − a ν q ν u µ , (40)where we used the fact that q · u = 0 and have introduced the kinematic acceleration a – not tobe confused with the acceleration of normal observers ˆ a – with components a µ := u λ ∇ λ u µ . Ourgoal is to obtain evolution equations in a conservative form for state variables that are orthogonalto n . Thus, we project Eq. (40) by multiplying with γ iµ : γ iµ ˙ q (cid:104) µ (cid:105) = u λ ∇ λ q i ⊥ − q µ u ν ∇ ν γ iµ − a ν q ν γ iµ u µ = u λ ∂ λ q i ⊥ − G i q − H i q , (41)where, upon using Eq. (33) and the definition of the covariant derivative, we have introduced thenew quantities G i q := Wα (cid:0) K kj v k − ˆ a j (cid:1) q j ⊥ β i + ˚ qW v j K ij − ˚ qW ˆ a i − Wα Γ i j q j ⊥ − Wα Γ ijk q j ⊥ V k , (42) H i q := (cid:16) a ⊥ j q j ⊥ − ˚ a ˚ q (cid:17) W v i . (43)The right-hand sides of Eqs. (42) and (43) also contain other new quantities, namely, the “co-ordinate velocity” V j := u j /u t = αv j − β j and the 3+1 split of the kinematic acceleration a µ := a µ (cid:107) + a µ ⊥ = ˚ an µ + a µ ⊥ , whose explicit components are given by a i ⊥ = A i + W Λ i − W v j K ij , (44) ˚ a = v i a i ⊥ = v i A i + W v i Λ i − W v i v j K ij , (45)with A i := W v j D j (cid:0) W v i (cid:1) , (46) Λ i := 1 α ( ∂ t − L β ) W v i + W ˆ a i . (47)and D j being the fully spatial part of the covariant derivative ∇ , i.e., D j v i := ∂ j v i + Γ ijk v k ,where Γ ijk represent the Christoffel symbols related to the three-metric γ (see Gourgoulhon2012; Rezzolla & Zanotti 2013, for details) Γ ijk := γ il ( ∂ j γ lk + ∂ j γ lk − ∂ l γ jk ) . (48) MNRAS000
RHD of non-perfect fluids E − D rather than simply E in equation (39) because the conservation of E − D is more accurate than that of E only. Comparing Eqs. (36)–(39) we can now read off thecomponents of U , F i and S corresponding to the conservation equations. We will write themdown explicitly in the next Section, after we have reformulated the constitutive equations (8)–(10)in conservative form. In order to do this, we consider the evolution equation for the heat current(9) as an example and extend the treatment to the other equations afterwards. We start with theterm ˙ q (cid:104) µ (cid:105) = h µν u λ ∇ λ q ν = u λ ∇ λ q µ − a ν q ν u µ , (40)where we used the fact that q · u = 0 and have introduced the kinematic acceleration a – not tobe confused with the acceleration of normal observers ˆ a – with components a µ := u λ ∇ λ u µ . Ourgoal is to obtain evolution equations in a conservative form for state variables that are orthogonalto n . Thus, we project Eq. (40) by multiplying with γ iµ : γ iµ ˙ q (cid:104) µ (cid:105) = u λ ∇ λ q i ⊥ − q µ u ν ∇ ν γ iµ − a ν q ν γ iµ u µ = u λ ∂ λ q i ⊥ − G i q − H i q , (41)where, upon using Eq. (33) and the definition of the covariant derivative, we have introduced thenew quantities G i q := Wα (cid:0) K kj v k − ˆ a j (cid:1) q j ⊥ β i + ˚ qW v j K ij − ˚ qW ˆ a i − Wα Γ i j q j ⊥ − Wα Γ ijk q j ⊥ V k , (42) H i q := (cid:16) a ⊥ j q j ⊥ − ˚ a ˚ q (cid:17) W v i . (43)The right-hand sides of Eqs. (42) and (43) also contain other new quantities, namely, the “co-ordinate velocity” V j := u j /u t = αv j − β j and the 3+1 split of the kinematic acceleration a µ := a µ (cid:107) + a µ ⊥ = ˚ an µ + a µ ⊥ , whose explicit components are given by a i ⊥ = A i + W Λ i − W v j K ij , (44) ˚ a = v i a i ⊥ = v i A i + W v i Λ i − W v i v j K ij , (45)with A i := W v j D j (cid:0) W v i (cid:1) , (46) Λ i := 1 α ( ∂ t − L β ) W v i + W ˆ a i . (47)and D j being the fully spatial part of the covariant derivative ∇ , i.e., D j v i := ∂ j v i + Γ ijk v k ,where Γ ijk represent the Christoffel symbols related to the three-metric γ (see Gourgoulhon2012; Rezzolla & Zanotti 2013, for details) Γ ijk := γ il ( ∂ j γ lk + ∂ j γ lk − ∂ l γ jk ) . (48) MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
Note that the vector G i q captures the influence of the choice of spacetime foliation as well asthe curvature of spacetime itself on the transport of heat on each hypersurface. On the other hand,the tensor H i q is a correction term that arises from the projection operator h µν , which ensures thatheat transport occurs always orthogonal to u . Notice that at this point G i q and H i q are not fully3+1 decomposed due to ˚ a , ˚ q and the corresponding Christoffel symbols not yet being fully split.Their full 3+1 decomposition will be given in Section 4.3 as both terms belong to the source terms.Now we exploit the continuity equation in the form (34) in order to modify the first term onthe right-hand side of Eq. (41): α √ γρu λ ∂ λ q i ⊥ = ∂ λ (cid:0) α √ γρq i ⊥ u λ (cid:1) = ∂ t (cid:0) √ γDq i ⊥ (cid:1) + ∂ j (cid:0) √ γDV j q i ⊥ (cid:1) . (49)Using this equation, as well as Eq. (41), and projecting Eq. (9) with γ iµ , we obtain ∂ t (cid:0) √ γDq i ⊥ (cid:1) + ∂ j (cid:0) √ γDV j q i ⊥ (cid:1) = α √ γDτ q W (cid:26) γ iµ q µ NS − q i ⊥ − κT q i ⊥ ∇ ν (cid:18) τ q u ν κT (cid:19) , + κT (cid:20) α γ iµ ∇ (cid:104) µ (cid:105) Π + α γ iµ ∇ ν π ν (cid:104) µ (cid:105) + (1 − γ )Π T γ iµ ∇ (cid:104) µ (cid:105) (cid:16) α T (cid:17) + (1 − γ ) T γ iµ π µν ∇ ν (cid:16) α T (cid:17) (cid:21) , + τ q (cid:2) G i q + H i q (cid:3)(cid:27) . (50)Proceeding in a similar way, we find that the projected version of equation (10) is given by ∂ t (cid:16) √ γDπ ij ⊥ (cid:17) + ∂ k (cid:16) √ γDV i π ij ⊥ (cid:17) = α √ γDτ π W (cid:26) γ iµ γ jν π µν NS − π ij ⊥ − ηT π ij ⊥ ∇ λ (cid:18) τ π u λ ηT (cid:19) , + 2 α ηγ iµ γ jν ∇ (cid:104) µ q ν (cid:105) + 2 γ ηT γ jν γ iµ q (cid:104) µ ∇ ν (cid:105) (cid:16) α T (cid:17) , + τ π (cid:2) G ijπ + H ijπ (cid:3)(cid:27) , (51)with G ijπ := Wα (cid:0) K lk v l − ˆ a k (cid:1) (cid:16) π ki ⊥ β j + π kj ⊥ β i (cid:17) + W v l (cid:0) K li ˚ π j + K lj ˚ π i (cid:1) − W (cid:0) ˆ a i ˚ π j + ˆ a j ˚ π i (cid:1) , − Wα (cid:16) Γ i k π kj ⊥ + Γ j k π ki ⊥ (cid:17) − Wα (cid:16) Γ ikl π kj ⊥ V l + Γ jkl π ki ⊥ V l (cid:17) , (52) H ijπ := W a ⊥ k (cid:16) π ki ⊥ v j + π kj ⊥ v i (cid:17) − W ˚ a (cid:0) ˚ π i v j + ˚ π j v i (cid:1) . (53)Finally, since the bulk-viscosity pressure Π is a scalar quantity, there are no projections thathave to be performed in order to write a 3+1 split version of Eq. (8): ∂ t ( √ γD Π) + ∂ i (cid:0) √ γDV i Π (cid:1) = α √ γDτ Π W (cid:20) Π NS − Π − ζ Π T ∇ µ (cid:18) τ Π u µ ζT (cid:19) , + α ζ ∇ µ q µ + γ ζT q µ ∇ µ (cid:16) α T (cid:17) (cid:21) . (54)In summary, Eqs. (37)–(39), (50), (51) and (54) can be combined into the flux-conservative MNRAS , 1–55 (2021)
RHD of non-perfect fluids U , F i and SU = √ γ DS j E − DD Π Dq j ⊥ Dπ jk ⊥ = √ γ ρW ( e + p + Π) W v j + W (cid:16) ˚ qv j + q j ⊥ (cid:17) + ˚ π j ( e + p + Π) W + 2˚ qW − ( p + Π − ˚ π ) − ρWρW Π ρW q j ⊥ ρW π jk ⊥ , (55) F i = √ γ DV i αS ij − β i S j α ( S i − v i D ) − β i ( E − D ) DV i Π DV i q j ⊥ DV i π jk ⊥ , (56) S = √ γ αS ik ∂ j γ ik + S i ∂ j β i − E∂ j ααS ij K ij − S j ∂ j α ( αD/τ Π W ) (Π NS − Π + ∆ Π )( αD/τ q W ) (cid:16) γ jµ q µ NS − q j ⊥ + ∆ j q + τ q G j q + τ q H j q (cid:17) ( αD/τ π W ) (cid:16) γ jµ γ kν π µν NS − π jk ⊥ + ∆ jkπ + τ π G jkπ + τ π H jkπ (cid:17) , (57) MNRAS000
RHD of non-perfect fluids U , F i and SU = √ γ DS j E − DD Π Dq j ⊥ Dπ jk ⊥ = √ γ ρW ( e + p + Π) W v j + W (cid:16) ˚ qv j + q j ⊥ (cid:17) + ˚ π j ( e + p + Π) W + 2˚ qW − ( p + Π − ˚ π ) − ρWρW Π ρW q j ⊥ ρW π jk ⊥ , (55) F i = √ γ DV i αS ij − β i S j α ( S i − v i D ) − β i ( E − D ) DV i Π DV i q j ⊥ DV i π jk ⊥ , (56) S = √ γ αS ik ∂ j γ ik + S i ∂ j β i − E∂ j ααS ij K ij − S j ∂ j α ( αD/τ Π W ) (Π NS − Π + ∆ Π )( αD/τ q W ) (cid:16) γ jµ q µ NS − q j ⊥ + ∆ j q + τ q G j q + τ q H j q (cid:17) ( αD/τ π W ) (cid:16) γ jµ γ kν π µν NS − π jk ⊥ + ∆ jkπ + τ π G jkπ + τ π H jkπ (cid:17) , (57) MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke where we have introduced the following new quantities: ∆ Π := − ζ Π T ∇ µ (cid:18) τ Π u µ ζT (cid:19) + α ζ ∇ µ q µ + γ ζT q µ ∇ µ (cid:16) α T (cid:17) , (58) ∆ j q := − κT q j ⊥ ∇ ν (cid:18) τ q u ν κT (cid:19) + κT (cid:20) α γ jµ ∇ (cid:104) µ (cid:105) Π + α γ jµ ∇ ν γ jµ ∇ (cid:104) µ (cid:105) (cid:16) α T (cid:17) + (1 − γ ) T γ jµ π µν ∇ ν (cid:16) α T (cid:17) (cid:21) , (59) ∆ jkπ := − ηT π jk ⊥ ∇ λ (cid:18) τ π u λ ηT (cid:19) + 2 α ηγ jµ γ kν ∇ (cid:104) µ q ν (cid:105) + 2 γ ηT γ kµ γ jν q (cid:104) µ ∇ ν (cid:105) (cid:16) α T (cid:17) . (60)The evolution equations (36) for the fifteen components of the state vector (57), with fluxes givenby the vectors (56), and source terms (57) represent the 3+1 flux-conservative formulation of the HL83 system.A few considerations are worth making at this point. First, although the whole solution ofthe set of partial differential equations becomes computationally more expensive, the conversionfrom the conserved to the primitive variables does not gain complexity, at least not beyond what isalready encountered in GRMHD. We recall that such a conversion, which needs to be performednumerically at each grid cell and on each timelevel of the solution, requires the solution of a setof nonlinear equations to obtain the values of the primitive variables from the newly computedconserved ones. Fortunately, the extension of the set of equations in GRDHD does not require newroot-finding steps to obtain the ten new primitive variables, i.e., Π , q i ⊥ , π ij ⊥ , as these are relatedalgebraically with the corresponding components of the state vector U . Second, the way in whichthe new conversion from the conserved to the primitive variables differs from the perfect-fluidcase is in the appearance of the three-velocities v i in the definitions of the conserved variables,i.e., Eqs. (29) – (30), as well as in the contractions with the projected dissipative currents q i ⊥ and π ij ⊥ . While the latter are related algebraically with the corresponding conserved variables,the three-velocity v i still requires a numerical root-finding and hence a continuous update in theroot-finding process. Ultimately, this leads to a numerical matrix inversion or additional root-finding steps in the conversion when compared to the perfect-fluid case. Second, a different choiceof second-order theory will only change the explicit expressions for ∆ Π , ∆ j q and ∆ jkπ . Stateddifferently, ∆ Π , ∆ j q and ∆ jkπ contain off-equilibrium contributions that are not captured by simplerelaxation-type equations and typically contain contributions that are of second and even higherorder. A few examples of other second-order GRDHD equations are presented in Appendix A,where we compare HL83 to other dissipative-hydrodynamics formulations. Finally, we remarkthat the sources for the dissipative quantities in Eqs. (57) are not yet fully 3+1 decomposed andtheir explicit expressions will be presented in the following Section.
MNRAS , 1–55 (2021)
RHD of non-perfect fluids As mentioned above, in order to provide complete expressions for the 3+1, flux-conservative for-mulation of Eqs. (36) we need to obtain right-hand sides that only contain fully spatial quantities,both for the ordinary variables – i.e., ρ , p , v i , α , β i , γ ij , K ij , as well as their spatial (partial) deriva-tives – and for the dissipative ones, as well as their temporal and spatial (partial and covariant)derivatives – i.e., Π , q i ⊥ and π ij ⊥ .To accomplish this second goal, we need to project the corresponding evolution equations –and corresponding right-hand sides – onto spatial hypersurfaces via the metric γ . However, sincemany terms in the source terms involve covariant derivatives [cf., Eqs. (58)–(60)], the resulting3+1 decomposed expressions are inevitably very complicated and lengthy, so that it is difficult toreconstruct the origin of the various source terms. To avoid this, and hence make the calculationsmore transparent and easy to follow, we collect the relevant source terms in the following classes:1. Intrinsically spatial terms.
These are terms containing only quantities which are originallydefined on a spatial hypersurface, i.e., ρ , p , v i , Π , q i ⊥ and π ij ⊥ , as well as spatial partial or covariantderivatives of these quantities. Also part of this class are terms proportional to γ ij and its spatialpartial derivatives. We have marked these terms in green.2. Terms not containing the extrinsic curvature K ij . These are terms involving temporal deriva-tives, spatial derivatives involving α or β i , as well as terms containing projections parallel to theunit normal n , e.g., ˚ q . We have marked these terms in orange.3. Terms containing the extrinsic curvature K ij . These are terms linear in the extrinsic curvature,its trace, K := γ ij K ij , but not containing spatial derivatives of K ij . We have marked these termsin light blue.However, before proceeding to this classification, we recall a number of useful identities thatwill be exploited in the derivation of the source terms and that are associated to gradients of thefluid four-velocity; in doing so we are following in part the convention introduced by Peitz & Appl(1997, 1999). We first recall that the covariant derivative of the four-velocity can be decomposedin terms of tensors that describe its properties in terms of changes in volumes, shape and vorticity,i.e., as (Rezzolla & Zanotti 2013) ∇ µ u ν = ω µν + σ µν + 13 Θ h µν − u µ a ν , (61)where Θ is the “expansion” and is defined as Θ : = ∇ µ u µ = ϑ + Λ − KW , (62)
MNRAS000
MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke and where we have introduced ϑ := D i (cid:0) W v i (cid:1) = ∂ i (cid:0) √ γW v i (cid:1) √ γ , (63) Λ := 1 α ( ∂ t − L β ) W + W v i ˆ a i . (64)Similarly, the “shear tensor” can be written as: σ µν : = ∇ (cid:104) µ u ν (cid:105) = ˚ σn µ n ν + ˚ σ µ n ν + ˚ σ ν n µ + σ µν ⊥ , (65)where σ ij ⊥ = γ iµ γ jν ∇ (cid:104) µ u ν (cid:105) = Σ ij + Λ ij − W K ij , (66)with Σ ij := 12 (cid:2) D i (cid:0) W v j (cid:1) + D j ( W v i ) (cid:3) + 12 W (cid:0) A i v j + A j v i (cid:1) − (cid:0) γ ij + W v i v j (cid:1) ϑ , (67) Λ ij := 12 W (cid:0) Λ i v j + Λ j v i (cid:1) −
13 Λ (cid:0) γ ij + W v i v j (cid:1) , (68) K ij := K ij + W v k (cid:0) K ik v j + K jk v i (cid:1) − K (cid:0) γ ij + W v i v j (cid:1) , (69) ˚ σ i = v j σ ij ⊥ , (70) ˚ σ = γ ij σ ij ⊥ = σ ⊥ ii . (71)Finally, the “kinematic vorticity” has the form: ω µν : = ∇ [ µ u ν ] + u [ µ a ν ] = ˚ ω µ n ν − ˚ ω ν n µ + ω µν ⊥ , (72)where ω ij ⊥ = Ω ij − W (cid:0) Λ i v i − Λ j v i (cid:1) + W v k (cid:0) K ik v j − K jk v i (cid:1) , (73)with Ω ij := D [ i W v j ] − W A [ i v j ] , and ˚ ω i := − n ν γ iλ ω λν = v j ω ij ⊥ . (74) MNRAS , 1–55 (2021)
RHD of non-perfect fluids Listed below are all the spatial source terms [ -th component of the vector S in (57)] appearing inthe evolution equation for the bulk-viscosity pressure, i.e., ∂ t (cid:0) √ γD Π (cid:1) = . . . Π NS = − ζ ( ϑ + Λ − KW ) , (75) ζ Π T ∇ µ (cid:18) τ Π u µ ζT (cid:19) = 12 τ Π Π ϑ + 12 ζ Π T W v i ∂ i (cid:18) τ Π ζT (cid:19) + 12 τ Π ΠΛ + ζ Π T W α ( ∂ t − L β ) (cid:18) τ Π ζT (cid:19) − τ Π Π W K , (76) ∇ µ q µ = D i q i ⊥ + 1 α ( ∂ t − L β ) ˚ q + ˆ a i q i ⊥ − ˚ qK , (77) q µ ∇ µ (cid:16) α T (cid:17) = q i ⊥ ∂ i (cid:16) α T (cid:17) + ˚ qα ( ∂ t − L β ) (cid:16) α T (cid:17) . (78) Similarly, listed below are all the spatial source terms [components - in Eq. (57)] appearing in theevolution equation for the components of the perpendicular heat current, i.e., ∂ t (cid:0) √ γDq i ⊥ (cid:1) = . . .γ iµ q µ NS = − κT (cid:20)(cid:0) γ ij + W v i v j (cid:1) ∂ j ln ( T ) + A i + W α v i ( ∂ t − L β ) ln ( T ) + W Λ i − W K ij v j (cid:21) , (79) κT q i ⊥ ∇ ν (cid:18) τ q u ν κT (cid:19) = 12 τ q q i ⊥ ϑ + 12 κT q i ⊥ W v j ∂ j (cid:16) τ q κT (cid:17) + 12 τ q q i ⊥ Λ + κT q i ⊥ W α ( ∂ t − L β ) (cid:16) τ q κT (cid:17) − τ q q i ⊥ W K , (80) γ iµ h µν ∇ ν Π = (cid:0) γ ij + W v i v j (cid:1) ∂ j Π + W α v i ( ∂ t − L β ) Π , (81) γ iµ ∇ ν π ν (cid:104) µ (cid:105) = D j π ij ⊥ + W v i v k D j π jk ⊥ − W v i D j ˚ π j + 1 α ( ∂ t − L β ) ˚ π i − ˆ a j π ij ⊥ − ˚ π ˆ a i + W v i (cid:20) v k α ( ∂ t − L β ) ˚ π k − α ( ∂ t − L β ) ˚ π − a j ˚ π j − ˚ πv k ˆ a k (cid:21) − K ij ˚ π j − W v i (cid:16) K jk ˚ π j v k − K jk π jk ⊥ − ˚ πK (cid:17) , (82) (1 − γ )Π T γ iµ ∇ (cid:104) µ (cid:105) (cid:16) α T (cid:17) = (1 − γ )Π T (cid:20) γ ij ∂ j + W v i v j ∂ j + W α ( ∂ t − L β ) (cid:21) (cid:16) α T (cid:17) , (83) (1 − γ ) T γ iµ π µν ∇ ν (cid:16) α T (cid:17) = (1 − γ ) T π ij ⊥ ∂ j (cid:16) α T (cid:17) + (1 − γ ) T ˚ π i α ( ∂ t − L β ) (cid:16) α T (cid:17) , (84) MNRAS000
RHD of non-perfect fluids Listed below are all the spatial source terms [ -th component of the vector S in (57)] appearing inthe evolution equation for the bulk-viscosity pressure, i.e., ∂ t (cid:0) √ γD Π (cid:1) = . . . Π NS = − ζ ( ϑ + Λ − KW ) , (75) ζ Π T ∇ µ (cid:18) τ Π u µ ζT (cid:19) = 12 τ Π Π ϑ + 12 ζ Π T W v i ∂ i (cid:18) τ Π ζT (cid:19) + 12 τ Π ΠΛ + ζ Π T W α ( ∂ t − L β ) (cid:18) τ Π ζT (cid:19) − τ Π Π W K , (76) ∇ µ q µ = D i q i ⊥ + 1 α ( ∂ t − L β ) ˚ q + ˆ a i q i ⊥ − ˚ qK , (77) q µ ∇ µ (cid:16) α T (cid:17) = q i ⊥ ∂ i (cid:16) α T (cid:17) + ˚ qα ( ∂ t − L β ) (cid:16) α T (cid:17) . (78) Similarly, listed below are all the spatial source terms [components - in Eq. (57)] appearing in theevolution equation for the components of the perpendicular heat current, i.e., ∂ t (cid:0) √ γDq i ⊥ (cid:1) = . . .γ iµ q µ NS = − κT (cid:20)(cid:0) γ ij + W v i v j (cid:1) ∂ j ln ( T ) + A i + W α v i ( ∂ t − L β ) ln ( T ) + W Λ i − W K ij v j (cid:21) , (79) κT q i ⊥ ∇ ν (cid:18) τ q u ν κT (cid:19) = 12 τ q q i ⊥ ϑ + 12 κT q i ⊥ W v j ∂ j (cid:16) τ q κT (cid:17) + 12 τ q q i ⊥ Λ + κT q i ⊥ W α ( ∂ t − L β ) (cid:16) τ q κT (cid:17) − τ q q i ⊥ W K , (80) γ iµ h µν ∇ ν Π = (cid:0) γ ij + W v i v j (cid:1) ∂ j Π + W α v i ( ∂ t − L β ) Π , (81) γ iµ ∇ ν π ν (cid:104) µ (cid:105) = D j π ij ⊥ + W v i v k D j π jk ⊥ − W v i D j ˚ π j + 1 α ( ∂ t − L β ) ˚ π i − ˆ a j π ij ⊥ − ˚ π ˆ a i + W v i (cid:20) v k α ( ∂ t − L β ) ˚ π k − α ( ∂ t − L β ) ˚ π − a j ˚ π j − ˚ πv k ˆ a k (cid:21) − K ij ˚ π j − W v i (cid:16) K jk ˚ π j v k − K jk π jk ⊥ − ˚ πK (cid:17) , (82) (1 − γ )Π T γ iµ ∇ (cid:104) µ (cid:105) (cid:16) α T (cid:17) = (1 − γ )Π T (cid:20) γ ij ∂ j + W v i v j ∂ j + W α ( ∂ t − L β ) (cid:21) (cid:16) α T (cid:17) , (83) (1 − γ ) T γ iµ π µν ∇ ν (cid:16) α T (cid:17) = (1 − γ ) T π ij ⊥ ∂ j (cid:16) α T (cid:17) + (1 − γ ) T ˚ π i α ( ∂ t − L β ) (cid:16) α T (cid:17) , (84) MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke H i q = W v i (cid:18) A j q j ⊥ + W Λ j q j ⊥ − ˚ qW v j Λ j − ˚ qv j A j − W K kj q k ⊥ v j + 2˚ qW K kj v k v j (cid:19) , (85) G i q = − W Γ ijk q j ⊥ v k − ˚ qW ˆ a i − Wα q j ⊥ ∂ j β i + ˚ qW v j K ij + W K ij q j ⊥ . (86) Finally, listed below are all the spatial source terms [components - in Eq. (57)] appearing inthe evolution equation for the components of the perpendicular shear-stress tensor, i.e., ∂ t (cid:16) √ γDπ ij ⊥ (cid:17) = . . . γ iµ γ jν π µν NS = − η (cid:0) Σ ij + Λ ij − W K ij (cid:1) , (87) ηT π ij ⊥ ∇ λ (cid:18) τ π u λ ηT (cid:19) = 12 τ π π ij ⊥ ϑ + 12 ηT π ij ⊥ W v i ∂ i (cid:18) τ π ηT (cid:19) + 12 τ π π ij ⊥ Λ + ηT π ij ⊥ W α ( ∂ t − L β ) (cid:18) τ π ηT (cid:19) − τ π π ij ⊥ W K , (88) γ iµ γ jν ∇ (cid:104) µ q ν (cid:105) = D ( i q j ) ⊥ + 12 W v i (cid:16) v k D k q j ⊥ + v k D j q k ⊥ (cid:17) + 12 W v j (cid:0) v k D k q i ⊥ + v k D i q k ⊥ (cid:1) − W v i v j A k q k ⊥ − (cid:0) γ ij + W v i v j (cid:1) (cid:0) D k q k ⊥ − A k q k ⊥ (cid:1) + 12 W v i (cid:20) α ( ∂ t − L β ) q j ⊥ + ˚ q ˆ a j − γ jk ∂ k ˚ q (cid:21) + 12 W v j (cid:20) α ( ∂ t − L β ) q i ⊥ + ˚ q ˆ a i − γ ik ∂ k ˚ q (cid:21) − W v i v j (cid:0) W Λ k q k ⊥ − ˚ qW Λ k v k − ˚ qA k v k (cid:1) − (cid:0) γ ij + W v i v j (cid:1) (cid:20) α ( ∂ t − L β ) ˚ q − W Λ k q k ⊥ + ˆ a k q k ⊥ + ˚ qW Λ k v k + ˚ qA k v k (cid:21) − ˚ qK ij − ˚ qW v k (cid:0) K ik v j + K jk v i (cid:1) + 2 W v i v j K kl (cid:0) q k ⊥ v l − ˚ qv k v l (cid:1) + 13 (cid:0) γ ij + W v i v j (cid:1) K kl (cid:0) qW v k v l − W q k ⊥ v l + ˚ qγ kl (cid:1) , (89) γ jν γ iµ q (cid:104) µ ∇ ν (cid:105) (cid:16) α T (cid:17) = (cid:26) (cid:20) q i ⊥ (cid:18) γ jk ∂ k + W v j v k ∂ k + W v j α ( ∂ t − L β ) (cid:19) + q j ⊥ (cid:18) γ ik ∂ k + W v i v k ∂ k + W v i α ( ∂ t − L β ) (cid:19) (cid:21) − (cid:0) γ ij + W v i v j (cid:1) (cid:18) q k ⊥ ∂ k + ˚ qα ( ∂ t − L β ) (cid:19) (cid:27) (cid:16) α T (cid:17) , (90) MNRAS , 1–55 (2021)
RHD of non-perfect fluids H ijπ = W A k (cid:16) π ik ⊥ v j + π jk ⊥ v i (cid:17) + W Λ k (cid:16) π ik ⊥ v j + π jk ⊥ v i (cid:17) − W (cid:0) v l A l + W v l Λ l (cid:1) (cid:0) ˚ π i v j + ˚ π j v i (cid:1) − W K lk v l (cid:16) π ik ⊥ v j + π jk ⊥ v i (cid:17) + 2 W K lk v l v k (cid:0) ˚ π i v j + ˚ π j v i (cid:1) , (91) G ijπ = − W v l (cid:16) Γ ikl π kj ⊥ + Γ jkl π ki ⊥ (cid:17) − W (cid:0) ˚ π i ˆ a j + ˚ π j ˆ a i (cid:1) − Wα (cid:16) π ik ⊥ ∂ k β j + π jk ⊥ ∂ k β i (cid:17) + W v l (cid:0) K li ˚ π j + K lj ˚ π i (cid:1) + W (cid:16) K il π lj ⊥ + K jl π li ⊥ (cid:17) . (92) As a corollary to the lengthy expressions provided above and as a way to help in the actual numeri-cal implementation of Eqs. (36), we provide below also the explicit expressions for the Christoffelsymbols appearing in the source terms (57) Γ = ∂ t ln α + ˆ a i β i − α K ij β i β j , (93) Γ i = ˆ a i − α K ij β j , (94) Γ ij = − α K ij , (95) Γ i = ∂ t β i − β i ∂ t ln α + β j ∂ j β i + 12 γ ij ∂ j α − β i β j ∂ j ln α + Γ ijk β j β k − αK ij β j + 1 α β i K jk β j β k , (96) Γ i j = ∂ j β i − β i ∂ j ln α + Γ ijk β k + 1 α β i K jk β k − αK ij , (97) Γ ijk = Γ ijk + 1 α β i K jk . (98)In summary, the 3+1 flux-conservative formulation of the general-relativistic Eqs. (36), com-bined with the explicit components Eqs. (55)–(57) and Eqs. (58)–(60), provides a complete andready-to-use set of equations for the numerical evaluation of dissipative effects in special-relativisticsimulations of colliding heavy ions as well as in general-relativistic simulations of compact ob-jects. In Appendix A we provide a detailed comparison of the system presented here with that ofother formulations (see Table A2). MNRAS000
RHD of non-perfect fluids H ijπ = W A k (cid:16) π ik ⊥ v j + π jk ⊥ v i (cid:17) + W Λ k (cid:16) π ik ⊥ v j + π jk ⊥ v i (cid:17) − W (cid:0) v l A l + W v l Λ l (cid:1) (cid:0) ˚ π i v j + ˚ π j v i (cid:1) − W K lk v l (cid:16) π ik ⊥ v j + π jk ⊥ v i (cid:17) + 2 W K lk v l v k (cid:0) ˚ π i v j + ˚ π j v i (cid:1) , (91) G ijπ = − W v l (cid:16) Γ ikl π kj ⊥ + Γ jkl π ki ⊥ (cid:17) − W (cid:0) ˚ π i ˆ a j + ˚ π j ˆ a i (cid:1) − Wα (cid:16) π ik ⊥ ∂ k β j + π jk ⊥ ∂ k β i (cid:17) + W v l (cid:0) K li ˚ π j + K lj ˚ π i (cid:1) + W (cid:16) K il π lj ⊥ + K jl π li ⊥ (cid:17) . (92) As a corollary to the lengthy expressions provided above and as a way to help in the actual numeri-cal implementation of Eqs. (36), we provide below also the explicit expressions for the Christoffelsymbols appearing in the source terms (57) Γ = ∂ t ln α + ˆ a i β i − α K ij β i β j , (93) Γ i = ˆ a i − α K ij β j , (94) Γ ij = − α K ij , (95) Γ i = ∂ t β i − β i ∂ t ln α + β j ∂ j β i + 12 γ ij ∂ j α − β i β j ∂ j ln α + Γ ijk β j β k − αK ij β j + 1 α β i K jk β j β k , (96) Γ i j = ∂ j β i − β i ∂ j ln α + Γ ijk β k + 1 α β i K jk β k − αK ij , (97) Γ ijk = Γ ijk + 1 α β i K jk . (98)In summary, the 3+1 flux-conservative formulation of the general-relativistic Eqs. (36), com-bined with the explicit components Eqs. (55)–(57) and Eqs. (58)–(60), provides a complete andready-to-use set of equations for the numerical evaluation of dissipative effects in special-relativisticsimulations of colliding heavy ions as well as in general-relativistic simulations of compact ob-jects. In Appendix A we provide a detailed comparison of the system presented here with that ofother formulations (see Table A2). MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
After having presented a fully general, causal, 3+1 split, and flux-conservative formulation of theequations of general-relativistic dissipative hydrodynamics, i.e., Eqs. (36) with explicit compo-nents (55)–(57), we now turn to a numerical implementation and the strategy that needs to bedeveloped when these equations have to be cast within an already developed GRHD or GRMHDcode. For simplicity, but also because the issues that will be discussed below would apply also formore complicated (and complete) forms of the equations, hereafter we concentrate on a reducedset after neglecting the heat current and the shear-stress tensor in Eqs. (55)–(57), i.e., after setting q i ⊥ = 0 = π ij ⊥ . Furthermore, we will also assume that the off-equilibrium contributions to thesource terms are very small, i.e., ∆ Π (cid:39) , or, equivalently, that Π relaxes towards its NS-valueonly, ignoring corrections coming from terms of order higher than one in Knudsen number. As aresult, the evolution equation for the bulk-viscosity pressure (54) reduces to ∂ t ( √ γD Π) + ∂ i (cid:0) √ γDV i Π (cid:1) = α √ γDτ Π W (Π NS − Π) = − α √ γDτ Π W [ ζ ( ϑ + Λ − KW ) + Π] , (99)which is a relaxation-type equation, describing the evolution of the bulk-viscosity pressure suchthat causality is not violated and stability is guaranteed.The specific implementation discussed here refers to the one made within the Black Hole Ac-cretion Code BHAC (Porth et al. 2017). The code, which solves the equations of GRMHD by meansof a finite-volume approach and high-resolution shock-capturing (HRSC) methods, assumes a sta-tionary but curved background, and it has been employed in a number of studies of accretion ontosupermassive black holes (Mizuno et al. 2018; Event Horizon Telescope Collaboration et al. 2019;Olivares et al. 2020). Of course, in our present implementation all of the electromagnetic fields areset to zero so that the various terms in Eqs. (36) reduce to U = √ γ DS j E − DD Π = √ γ ρW ( e + p + Π) W v j ( e + p + Π) W − ( p + Π) − ρWρW Π , (100) Note that setting to zero the spatial components of the heat current and of the shear-stress tensor implies that also the time components are zero[cf., Eq. (27)]. MNRAS , 1–55 (2021)
RHD of non-perfect fluids F i = √ γ V i DαS ij − β i S j α ( S i − v i D ) − β i ( E − D ) V i D Π , (101) S = √ γ αS ik ∂ j γ ik + S i ∂ j β i − E∂ j α S ik β j ∂ j γ ik + S ij ∂ j β i − S j ∂ j α − ( αD/τ Π W ) [ ζ ( ϑ + Λ − KW ) + Π] . (102)Clearly, the extra terms that need to be handled in Eqs. (99) either involve divergences orpartial derivatives, which can be evaluated using the standard differential operators available within BHAC . For example, the divergence appearing in the term ϑ is discretised at second order as ϑ = ∂ i ( √ γW v i ) √ γ = (cid:82) V ∂ i ( √ γW v i ) dV (cid:82) V √ γdV = 1∆ V (cid:88) i =1 (cid:104)(cid:0) W v i ∆ S i (cid:1) ( x i +∆ x i / − (cid:0) W v i ∆ S i (cid:1) ( x i − ∆ x i / (cid:105) , (103)where the cell volume and cell surfaces are defined respectively as (Porth et al. 2017) ∆ V := (cid:90) V √ γ dV , ∆ S i ( x i +∆ x i / := (cid:90) ( x i +∆ x i / √ γ dS i , (104)Each integral is performed over one cell, where the coordinate volume and surface are given by dV := dx dx dx and dS i := s i dx j (cid:54) = i dx k (cid:54) = i , respectively. The co-vector s i is the i th componentof the unit normal with respect to the boundary of the cell. In addition, we calculate the boundarydata through averages; for example, given a scalar function φ , we compute φ ( x i +∆ x i / := ( φ ( x i +∆ x i ) + φ ( x i )) / , φ ( x i − ∆ x i / := ( φ ( x i − ∆ x i ) + φ ( x i )) / , and ˆ φ := φ ( x i ) .Similarly, the volume-averaged partial derivative of a scalar function φ (e.g., W ) ∂ i φ = (cid:82) V √ γ∂ i φ dV (cid:82) V √ γ dV = 1∆ V (cid:104)(cid:0) φ ∆ S i (cid:1) ( x i +∆ x i / − (cid:0) φ ∆ S i (cid:1) ( x i − ∆ x i / − ˆ φ (cid:0) ∆ S i (cid:1) ( x i +∆ x i / + ˆ φ (cid:0) ∆ S i (cid:1) ( x i − ∆ x i / (cid:105) . (105)Also, the time derivative, such as the one appearing in the term Λ , ( ∂ t − β i ∂ i ) W , is calculated MNRAS000
RHD of non-perfect fluids F i = √ γ V i DαS ij − β i S j α ( S i − v i D ) − β i ( E − D ) V i D Π , (101) S = √ γ αS ik ∂ j γ ik + S i ∂ j β i − E∂ j α S ik β j ∂ j γ ik + S ij ∂ j β i − S j ∂ j α − ( αD/τ Π W ) [ ζ ( ϑ + Λ − KW ) + Π] . (102)Clearly, the extra terms that need to be handled in Eqs. (99) either involve divergences orpartial derivatives, which can be evaluated using the standard differential operators available within BHAC . For example, the divergence appearing in the term ϑ is discretised at second order as ϑ = ∂ i ( √ γW v i ) √ γ = (cid:82) V ∂ i ( √ γW v i ) dV (cid:82) V √ γdV = 1∆ V (cid:88) i =1 (cid:104)(cid:0) W v i ∆ S i (cid:1) ( x i +∆ x i / − (cid:0) W v i ∆ S i (cid:1) ( x i − ∆ x i / (cid:105) , (103)where the cell volume and cell surfaces are defined respectively as (Porth et al. 2017) ∆ V := (cid:90) V √ γ dV , ∆ S i ( x i +∆ x i / := (cid:90) ( x i +∆ x i / √ γ dS i , (104)Each integral is performed over one cell, where the coordinate volume and surface are given by dV := dx dx dx and dS i := s i dx j (cid:54) = i dx k (cid:54) = i , respectively. The co-vector s i is the i th componentof the unit normal with respect to the boundary of the cell. In addition, we calculate the boundarydata through averages; for example, given a scalar function φ , we compute φ ( x i +∆ x i / := ( φ ( x i +∆ x i ) + φ ( x i )) / , φ ( x i − ∆ x i / := ( φ ( x i − ∆ x i ) + φ ( x i )) / , and ˆ φ := φ ( x i ) .Similarly, the volume-averaged partial derivative of a scalar function φ (e.g., W ) ∂ i φ = (cid:82) V √ γ∂ i φ dV (cid:82) V √ γ dV = 1∆ V (cid:104)(cid:0) φ ∆ S i (cid:1) ( x i +∆ x i / − (cid:0) φ ∆ S i (cid:1) ( x i − ∆ x i / − ˆ φ (cid:0) ∆ S i (cid:1) ( x i +∆ x i / + ˆ φ (cid:0) ∆ S i (cid:1) ( x i − ∆ x i / (cid:105) . (105)Also, the time derivative, such as the one appearing in the term Λ , ( ∂ t − β i ∂ i ) W , is calculated MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke using central second-order finite differences and a two-step predictor-corrector time evolution ( ∂ t W ) c ≈ W ( t + ∆ t/ − W ( t − ∆ t/ t , (106)where W ( t − ∆ t/ is saved from the previous timestep and is also used to calculate the predictorstep via first-order backward differences ( ∂ t W ) p ≈ [ W ( t ) − W ( t − ∆ t/ / (∆ t/ , with W ( t +∆ t/ being the Lorentz factor after the predictor step. Finally, for stationary spacetimes, as theone considered here, the trace of the extrinsic curvature can be computed as K = 1 α ∂ i β i + 12 α γ ij β k ∂ k γ ij . (107)Note that the newly evolved bulk-viscosity pressure appears also in the definitions of E , S j , S ij ,which need to be suitably updated. Fortunately, Π only appears as a correction to the equilibriumpressure p in the total pressure of the fluid p t = p + Π . Therefore, we can rewrite equation (100)–(102) in terms of p t , of the total specific enthalpy h t := h + Π ρ = e + p + Π ρ . (108)and of an effective EOS, p t = p t ( ρ, e, Π) . As a result, when considering a simple ideal-gas EOS,where the pressure is given by p = ( e − ρ )( γ − , (109)and γ is the adiabatic index, the new effective EOS would then read p t = ( e − ρ )( γ −
1) + Π , (110)Note that for a generic EOS p = p ( n, e ) and particle rest mass m , the effective sound speed in thepresence of bulk viscosity is given by (Bemfica et al. 2019b) c s,t := (cid:18) ∂p∂e (cid:19) n + 1 mh t (cid:18) ∂p∂n (cid:19) e + ζτ Π ρh t , (111)where n is the number density and ρ = nm . Clearly, c s,t > c s and so the correct characteristicwave-speeds need to be updated in the presence of bulk viscosity.As a result, in the code we make use only of the total pressure, of the total specific enthalpy,and of the corresponding sound speed. In other words, in the recovery of the primitive variablesfrom the conserved ones, we perform the following mapping ρh = ρ + γγ − p → ρh t = ρ + γγ − p t − γ − , (112) p = γ − γ ( ρh − ρ ) → p t = ρ (cid:18) γ − γ (cid:19) ( h t −
1) + Π γ , (113) c s = ( γ − h − h → c s,t = ( γ − h t − h t + ζτ Π ρh t . (114) MNRAS , 1–55 (2021)
RHD of non-perfect fluids
In this Section we present the application of the relativistic dissipative-hydrodynamics Eqs. (100)–(102) described in Section 5 to two 1+1-dimensional tests in special relativity. The time integrationis carried out explicitly by using a two-step predictor-corrector scheme. Furthermore, we employthe so-called Rusanov (Rusanov 1961) or “TVDLF” flux at the cell boundaries, with wave-speedgiven by the expression for c s,t in Eq. (114), although the differences when using instead theexpressions for c s are at most of the order of . for the cases considered here. Furthermore,we use the “minmod” reconstruction scheme to compute state variables at cell boundaries (see,e.g., Rezzolla & Zanotti 2013); for more details on the numerical schemes we refer the interestedreader to Porth et al. (2017). We first consider the time-honoured Bjorken-flow (Bjorken 1983) when bulk viscosity is present,which still represents a well-known and often-used test in relativistic dissipative hydrodynamics(see, e.g., Del Zanna et al. 2013; Inghirami et al. 2016, 2018) and MHD with transverse magneticfields (Roy et al. 2015; Pu et al. 2016). We recall that the Bjorken flow is the idealised repre-sentation of a one-dimensional, longitudinally boost-invariant motion of a fluid, such as the oneproduced in an ultrarelativistic collision of two ions.As usual for the Bjorken-flow scenario, it is convenient to make use of the so-called Milnecoordinates, which are given by cτ := (cid:112) ( ct ) − z , η := artanh (cid:16) zct (cid:17) = 12 ln (cid:18) z/ ( ct )1 − z/ ( ct ) (cid:19) , (115) t = τ cosh η , z = cτ sinh η , (116)where τ and η are defined as the proper time and space-time rapidity, respectively. In our setup,the evolution starts at τ = 1 fm c − and ends at τ = 15 fm c − . The initial conditions are given by ( ρ, p, v, Π ) = (cid:0) − GeV c − fm − , . − , . , . (cid:1) . (117)As the solution does not depend on the coordinate η , all solutions are taken at η = 0 . . Finally, wechoose constant values for ζ and τ Π , as reported in the legend of Fig. 1. The analytic solution for MNRAS000
In this Section we present the application of the relativistic dissipative-hydrodynamics Eqs. (100)–(102) described in Section 5 to two 1+1-dimensional tests in special relativity. The time integrationis carried out explicitly by using a two-step predictor-corrector scheme. Furthermore, we employthe so-called Rusanov (Rusanov 1961) or “TVDLF” flux at the cell boundaries, with wave-speedgiven by the expression for c s,t in Eq. (114), although the differences when using instead theexpressions for c s are at most of the order of . for the cases considered here. Furthermore,we use the “minmod” reconstruction scheme to compute state variables at cell boundaries (see,e.g., Rezzolla & Zanotti 2013); for more details on the numerical schemes we refer the interestedreader to Porth et al. (2017). We first consider the time-honoured Bjorken-flow (Bjorken 1983) when bulk viscosity is present,which still represents a well-known and often-used test in relativistic dissipative hydrodynamics(see, e.g., Del Zanna et al. 2013; Inghirami et al. 2016, 2018) and MHD with transverse magneticfields (Roy et al. 2015; Pu et al. 2016). We recall that the Bjorken flow is the idealised repre-sentation of a one-dimensional, longitudinally boost-invariant motion of a fluid, such as the oneproduced in an ultrarelativistic collision of two ions.As usual for the Bjorken-flow scenario, it is convenient to make use of the so-called Milnecoordinates, which are given by cτ := (cid:112) ( ct ) − z , η := artanh (cid:16) zct (cid:17) = 12 ln (cid:18) z/ ( ct )1 − z/ ( ct ) (cid:19) , (115) t = τ cosh η , z = cτ sinh η , (116)where τ and η are defined as the proper time and space-time rapidity, respectively. In our setup,the evolution starts at τ = 1 fm c − and ends at τ = 15 fm c − . The initial conditions are given by ( ρ, p, v, Π ) = (cid:0) − GeV c − fm − , . − , . , . (cid:1) . (117)As the solution does not depend on the coordinate η , all solutions are taken at η = 0 . . Finally, wechoose constant values for ζ and τ Π , as reported in the legend of Fig. 1. The analytic solution for MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke τ [fm c − ] − . − . − . − . − . . Π [ G e V f m − ] τ Π = 1 . c − , ζ = 0 .
01 GeV fm − analyticΠ NS τ Π = 1 . c − , ζ = 0 .
05 GeV fm − analyticΠ NS Figure 1.
Evolution of the bulk-viscosity pressure in a Bjorken flow. Solid lines show the analytic solution, while dashed lines the solution in theNS approximation. Filled circles report instead the numerical solution from
BHAC . this problem is given by Π ( τ ) = Π ( τ ) exp [ − ( τ − τ ) /τ Π ] + ζcτ Π exp ( − τ /τ Π ) [Ei ( τ /τ Π ) − Ei ( τ /τ Π )] , (118)where Ei ( x ) is the exponential integral function (see, e.g., Del Zanna et al. 2013, for more details).Figure 1 shows the evolution of the bulk-viscosity pressure Π for the initial data given byEq. (117) (the changes in the other hydrodynamical quantities are very small and not particularlyinteresting). As it can be seen, the numerical solutions agree very well with the analytical ones.Also, note that with time, Π converges towards its NS value Π NS = − ζ Θ = − ζ/ ( cτ ) (dashedlines in Fig. 1). At any time during the evolution, the relative difference between the analyticand numerical solution is − at most, confirming the correct implementation of the relativis-tic dissipative-hydrodynamics equations for smooth flows in a flat spacetime. Furthermore, ourchoices for the pair ( ζ, τ ) lie well within the range of applicability for the equations of GRDHD.This can be seen by calculating the effective speed of sound which is given by c s,t c = 13 + ζcτ Π p + Π , (119)Requiring causality, i.e., c s,t /c < , and expressing p through the perfect-fluid solution p ( τ ) = p ( τ ) ( τ o /τ ) / we obtain that τ < (cid:18) cτ Π p ( τ ) ζ (cid:19) / τ , (120)when p (cid:29) | Π | .For our choices ( ζ, τ Π ) = (0 .
01 GeV fm − , . c − ) and ( ζ, τ Π ) = (0 .
05 GeV fm − , . c − ) MNRAS , 1–55 (2021)
RHD of non-perfect fluids | Π | /p remains below 1% and 5%, respectively, up to the time τ = 400 fm c − .Hence, we apply Eq. (120) and find τ <
371 fm c − as well as τ <
110 fm c − , respectively. Bothvalues agree very well with the ones obtained from the exact solution and lie clearly above the endof the simulation at τ = 15 fm c − . We next explore the solution of a shock-tube problem for an ultra-relativistic gas of gluons. Whilethis is a standard 1+1-dimensional test scenario, we here use the same setup implemented by(Bouras et al. 2009a; Gabbana et al. 2020), i.e., we consider the ideal-gas EOS relative to an ultra-relativistic fluid (i.e., γ = 4 / ). Adopting Cartesian coordinates, the spatial domain ranges from x = − . to x = 3 . and the initial discontinuity in pressure and density is located at x = 0 . , while the velocity and bulk-viscosity pressure are assumed to be zero initially. Inother words, the initial conditions are given by ( T, p, v,
Π ) = (cid:0) . k − , .
43 GeV fm − , . , . (cid:1) x < . , (cid:0) . k − , .
33 GeV fm − , . , . (cid:1) x ≥ . . (121)On the other hand, we parametrize the bulk-viscosity coefficient ζ in terms of the entropy densityof the fluid, namely, s = ρ k B m (cid:20) − ln (cid:18) π ρm d F T c (cid:126) k (cid:19)(cid:21) , (122)where d F denotes the number of degrees of freedom and is set to for gluons, such that ζ = 43 k B c (cid:126) ζ s , (123)The coefficient ζ is a non-negative number, for which we choose the values ζ = { . , . , . } to obtain a direct comparison with the data from Gabbana et al. (2020). Note that the equationsdescribing the shock-tube problem with bulk viscosity in one dimension take the same form asthe corresponding equations with shear viscosity; the latter has been investigated in the work ofBouras et al. (2010), as well as more recently by Gabbana et al. (2020). This leads to the mapping ζ = 4 / η between the bulk viscosity employed in this work and the shear viscosity used in Bouraset al. (2010) and Gabbana et al. (2020). Furthermore, to obtain the correct rest-mass density weassume a single-particle rest mass m = 0 . c − , so that ρ (cid:28) e , as required by wanting toconsider an ultra-relativistic limit. Finally, the relaxation time used by Gabbana et al. (2020) is In this Section, to facilitate the comparison with codes designed for describing HICs (Gabbana et al. 2020) we adopt physical units, where k B and (cid:126) are the Boltzmann and reduced Planck constants, respectively.MNRAS000
33 GeV fm − , . , . (cid:1) x ≥ . . (121)On the other hand, we parametrize the bulk-viscosity coefficient ζ in terms of the entropy densityof the fluid, namely, s = ρ k B m (cid:20) − ln (cid:18) π ρm d F T c (cid:126) k (cid:19)(cid:21) , (122)where d F denotes the number of degrees of freedom and is set to for gluons, such that ζ = 43 k B c (cid:126) ζ s , (123)The coefficient ζ is a non-negative number, for which we choose the values ζ = { . , . , . } to obtain a direct comparison with the data from Gabbana et al. (2020). Note that the equationsdescribing the shock-tube problem with bulk viscosity in one dimension take the same form asthe corresponding equations with shear viscosity; the latter has been investigated in the work ofBouras et al. (2010), as well as more recently by Gabbana et al. (2020). This leads to the mapping ζ = 4 / η between the bulk viscosity employed in this work and the shear viscosity used in Bouraset al. (2010) and Gabbana et al. (2020). Furthermore, to obtain the correct rest-mass density weassume a single-particle rest mass m = 0 . c − , so that ρ (cid:28) e , as required by wanting toconsider an ultra-relativistic limit. Finally, the relaxation time used by Gabbana et al. (2020) is In this Section, to facilitate the comparison with codes designed for describing HICs (Gabbana et al. 2020) we adopt physical units, where k B and (cid:126) are the Boltzmann and reduced Planck constants, respectively.MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke given in terms of the bulk viscosity by τ Π = 1516 ζpc . (124)The numerical solution of the shock-tube problem at time t = 3 . c − is shown in Fig. 2,whose upper panels report the behaviour of the pressure normalized to p := 5 .
43 GeV fm − (thetop right panel is a magnification of the top left panel), while the bottom panels show the solutionof the velocity and bulk-viscosity pressure normalized to the fluid pressure. Different lines refereither to solutions obtained with BHAC for different values of ζ , or to solutions obtained with arelativistic lattice-Boltzmann (RLBM) approach (dashed lines) or to solutions of the relativisticBoltzmann equation via the test-particle (RBMTP) approach (dotted lines). Note that the case ζ = 0 . is essentially indistinguishable from an inviscid solution with the precision shown inthe figure and hence can be taken as the perfect-fluid reference.Figure 2 highlights how the strong spatial gradients present in the initial conditions tend tobe washed out by the presence of bulk viscosity and that this smearing of the discontinuities islarger with increasing bulk viscosity. Note that the wave-pattern of perfect-fluid hydrodynamics– which consists of a rarefaction wave and of a shock wave – can still be clearly identified if thebulk viscosity is not too large, i.e., ζ (cid:46) . (see, e.g., the pressure in the upper left panel of Fig.2). Furthermore, the solution behaviour is in good agreement with the results obtained by usingthe RLBM approach. However, in the case of high bulk-viscosity, i.e., ζ = 0 . , the wave-patternof perfect-fluid hydrodynamics is so strongly smeared out that it is difficult to clearly distinguishthe rarefaction wave from the shock wave. This is not surprising, since such large values of thebulk viscosity effectively correspond to a regime of large Knudsen number, which is where thehydrodynamical approach – and hence the formation of shock waves – is expected to fail.Interestingly, and as pointed out by Denicol et al. (2008b) and Bouras et al. (2010), threeadditional discontinuities are present in the high-viscosity, ζ = 0 . , case, two of which can be seenin the upper right panel of Fig. 2, where one is located at the head of the right-propagating shock,while the other near the contact discontinuity of the corresponding inviscid ( ζ = 0 . ) case . Thethird additional discontinuity is located at the head of the left-propagating rarefaction fan which isnot shown here. However, all discontinuities transition smoothly into the wave-pattern of perfect-fluid hydrodynamics at later times; indeed, the pressure jump located near the contact discontinuitydecays to less than 10% of its initial size of (cid:39) . GeV fm − by a time of t (cid:39) . fm c − (see Obviously the contact discontinuity cannot be seen in the pressure profile, but it is apparent in the rest-mass density profile, which is not shownin Fig. 2. MNRAS , 1–55 (2021)
RHD of non-perfect fluids − x [fm]0 . . . . . p / p BHAC ζ = 0 . ζ = 0 . ζ = 0 . ζ = 0 . ζ = 0 . ζ = 0 . . . . . . x [fm]0 . . . . . . . . . p / p − x [fm]0 . . . . . . v [ c ] − x [fm] − . . . . . Π / p . . . . . . . Figure 2.
Solution of the shock-tube test with initial data given by Eq. (121) and at time t = 3 . c − for. Top left: solution of the pressurenormalized to p = 5 .
43 GeV fm − ; top right: same as in top left but zoomed-in at the shock front; bottom left: solution of the three-velocity, bottom right: bulk-viscosity pressure normalised to the fluid pressure. In all panels, the solid lines show the results obtained with BHAC , while thedashed or dotted lines present the results using the relativistic lattice-Boltzmann (RLBM) or the test-particle (RBMTP) approach, respectively. also Bouras et al. 2010, for a more detailed description and a possible explanation). In addition,we find in agreement with Bouras et al. (2009b) that at a time t (cid:39) . fm c − the fluid velocitydownstream of the shock front differs by less than 0.5% from the corresponding fluid velocity ofthe inviscid case. At this point the previously mentioned pressure jump has declined to less than2% of its original size. Because the new discontinuities are absent in the RLBM (dashed line) andRBMTP (dotted line) solutions, their appearance may well indicate a breakdown of dissipativehydrodynamics for such large viscosities, i.e., for such large Knudsen numbers, which is of coursestill correctly described by the microscopic approaches RLBM and RBMTP.In addition, in the high-viscosity case, the BHAC solution underestimates the bulk-viscositypressure at the head of the right-propagating shock with respect to the values computed with theRLBM or RBMTP approaches (see bottom right panel of Fig. 2). While this may again be due to
MNRAS000
MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke the breakdown of the hydrodynamical description, part of the error may also originate from thetruncation of the relativistic evolution equation for Π [we recall that we have set ∆ Π = 0 in Eq.(57)]. As remarked by Bouras et al. (2010), the inclusion of additional source terms, as well as ofa coupling to a heat current, generally yields a better description of the corresponding dissipativecurrent (see Fig. 10 of Bouras et al. 2010). We expect the same to be true here and hence that asmaller deviation would be obtained with a more sophisticated source term for Eq. (99).Overall, Fig. 2 shows that our numerical implementation of the relativistic dissipative-hydrody-namics equations leads to solutions that are in very good agreement with the reference solutionsobtained by the direct solution of the relativistic Boltzmann equation in regimes that are mildlydissipative, i.e., ζ (cid:46) . , and in regimes that are highly dissipative, i.e., ζ (cid:46) . . The relativedifferences remain below 8% for p/p where we only considered the region x ∈ [ − . , .
6] fm . Inthis way we exclude rightmost shock wave and the head of the left-propagating rarefaction wave.Due to the fact that the solution obtained by BHAC reaches its unperturbed initial state at lower | x | than the solutions obtained by RLBM and RBMTP, the relative difference starts growing outsideof [ − . , .
6] fm and can reach values close to 45%. Note that this relative difference dependssolely on the value of p in the unperturbed initial states, for instance if the pressure of the rightinitial state was close to zero then the relative difference would be close to 100% at the rightmostshock wave. Furthermore, we similarly calculate the relative differences for Π /p and obtain valuesbelow a maximum of ∼
60% which is reached at x ∼ . Note that Π /p passes through zerosuch that the corresponding region is excluded, too. In this Section we present a stationary solution of the spherically symmetric equations of GRDHDin a Schwarzschild spacetime. In the following, we return to using geometrised units. For perfectfluids, the fully general-relativistic solution is known as the so-called “Michel solution” (Michel1972) and, together with the “Bondi-Hoyle” solution (Bondi 1952), serves as the reference solutionfor models of accreting nonrotating black holes in spherical symmetry, (see e.g., Nobili et al.1991, and references therein) as well as serves as a testbed for GRHD and GRMHD codes (see,e.g., Hawley et al. 1984; Porth et al. 2017; Weih et al. 2020).The effects of shear viscosity – such as the one arising from turbulent motion – have first beenconsidered by Turolla & Nobili (1989), who however adopted a description in terms of the general-relativistic NS equations. However, already in this simplified setup, Turolla & Nobili (1989) have
MNRAS , 1–55 (2021)
RHD of non-perfect fluids
BHAC in a curved spacetime geometry.Details on the derivation and solution of the ODEs can be found in Appendix B.We assume the fluid to be a mixture of ionised non-relativistic hydrogen coupled to photonsby employing the following EOS (see, e.g., Rezzolla & Zanotti 2013) p = p M (1 + α ) . (125)Here, p M denotes the pressure of the matter component, which is assumed to be an ideal gas, whilethe contribution from the radiation component is fixed using the parameter α . By rearranging Eq.(125), we find that the EOS of the total mixture takes the same form as the EOS for an ideal gashaving an effective adiabatic index γ e : p = ( γ e − e − ρ ) , (126)where γ e = 1 + 2(1 + α ) / [3(1 + 2 α )] . Note that e is the total energy density of the mixture while ρ denotes the rest-mass density of the hydrogen ions. The effective adiabatic index γ e should not beconfused with the generalized adiabatic exponent Γ which instead is defined through a thermo-dynamic relation: Γ := ( ∂ ln p/∂ ln ρ ) ˜ s = (5 / α + 16 α ) / [(3 / α )(1 + α )] , where ˜ s denotes the specific entropy (see, e.g., Mihalas & Mihalas 1984; Rezzolla & Zanotti 2013). Inaddition, the effective adiabatic index γ e can be expressed through the well-known generalizedadiabatic exponent Γ which is defined by Γ − ∂ ln T /∂ ln ρ ) ˜ s = (1 + 4 α ) / (3 / α ) : γ e = 1 + 9 − − . (127)The temperature can be obtained from the ideal-fluid law for the matter component p M = 2( k B /m p ) ρT which yields: T = 12(1 + α ) m p k B pρ , (128)where m p denotes the proton mass. Note the appearance of a factor in the denominator comingfrom the electrons in our charge-neutral plasma. Furthermore, we use a modification of the formulafor radiative bulk viscosity given by Weinberg (1971) and Sawyer (2006) ζ = 4 ζ σ SB T τ mfp (cid:18) − γ e (cid:19) , (129)where σ SB is the Stefan-Boltzmann constant and τ mfp := m p ρ − σ − the mean free time for a MNRAS000
BHAC in a curved spacetime geometry.Details on the derivation and solution of the ODEs can be found in Appendix B.We assume the fluid to be a mixture of ionised non-relativistic hydrogen coupled to photonsby employing the following EOS (see, e.g., Rezzolla & Zanotti 2013) p = p M (1 + α ) . (125)Here, p M denotes the pressure of the matter component, which is assumed to be an ideal gas, whilethe contribution from the radiation component is fixed using the parameter α . By rearranging Eq.(125), we find that the EOS of the total mixture takes the same form as the EOS for an ideal gashaving an effective adiabatic index γ e : p = ( γ e − e − ρ ) , (126)where γ e = 1 + 2(1 + α ) / [3(1 + 2 α )] . Note that e is the total energy density of the mixture while ρ denotes the rest-mass density of the hydrogen ions. The effective adiabatic index γ e should not beconfused with the generalized adiabatic exponent Γ which instead is defined through a thermo-dynamic relation: Γ := ( ∂ ln p/∂ ln ρ ) ˜ s = (5 / α + 16 α ) / [(3 / α )(1 + α )] , where ˜ s denotes the specific entropy (see, e.g., Mihalas & Mihalas 1984; Rezzolla & Zanotti 2013). Inaddition, the effective adiabatic index γ e can be expressed through the well-known generalizedadiabatic exponent Γ which is defined by Γ − ∂ ln T /∂ ln ρ ) ˜ s = (1 + 4 α ) / (3 / α ) : γ e = 1 + 9 − − . (127)The temperature can be obtained from the ideal-fluid law for the matter component p M = 2( k B /m p ) ρT which yields: T = 12(1 + α ) m p k B pρ , (128)where m p denotes the proton mass. Note the appearance of a factor in the denominator comingfrom the electrons in our charge-neutral plasma. Furthermore, we use a modification of the formulafor radiative bulk viscosity given by Weinberg (1971) and Sawyer (2006) ζ = 4 ζ σ SB T τ mfp (cid:18) − γ e (cid:19) , (129)where σ SB is the Stefan-Boltzmann constant and τ mfp := m p ρ − σ − the mean free time for a MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
Model ζ τ [10 − ] low ζ
160 1 . medium ζ ,
000 1 . high ζ ,
000 1 . low τ Π ,
000 0 . high τ Π ,
000 50
Table 1.
Summary of the various models evolved and their corresponding parameters ζ and τ . photon in the mixture. The constant σ T denotes the Thomson scattering cross-section. The dimen-sionless constant ζ is essentially arbitrary and used here to explore the regimes of low and highbulk viscosities. For the relaxation time τ Π we choose the parametrization: τ Π = τ M | ˙ M | (cid:16) r M (cid:17) , (130)where r denotes the circumference radius in Schwarzschild coordinates, M the mass of the blackhole, ˙ M the accretion rate, and τ is a dimensionless parameter to study short and long relaxationtimes. Note that in Eq. (130) the relaxation time increases cubically with radius. Adopting thischoice has turned out to be necessary to prevent the rapid growth of non-equilibrium effects nearthe sonic point, which can yield unphysical solutions.All of the models considered have the sonic point r s at M , a black-hole mass of M = 3 M (cid:12) ,and α = 1 , which then yields γ e = 1 . . Furthermore, we set the constants of motion ˙ M and B from their inviscid values computed at r s = 200 M , i.e., ˙ M = − . and B = − . in code units. Note that we use a polytropic EOS to solve the inviscid case where the polytropicconstant of the fluid mixture is given by k = 2(1 + α ) T ∞ ρ − γ e ∞ with asymptotic values ρ ∞ =2 × − g cm − and T ∞ = 1 . × K . We recall that B is the viscous analogue of the relativisticBernoulli constant − hu t ; (see Appendix B for a definition). Our code units are defined by setting m p /k B = 1 and k = 1 . In these units, a polytropic EOS with polytropic index γ e passing throughthe pair ( ρ ∞ , T ∞ ) is simply given by p = ρ γ e , where T and p are related through Eq. (128). In thisway, the calculation of the corresponding inviscid accretion problem simplifies because withoutdissipative losses the accretion process of an ideal gas is isentropic and can be described by apolytropic EOS.In the following, we consider five different accretion cases, whose parameters are given inTable 1. We employ a grid of 10,000 cells, which ranges from . M to , M using horizon-penetrating Kerr-Schild coordinates and show profiles at a time , M unless stated otherwise.The inviscid solution for the temperature T , together with the high ζ and low τ models are MNRAS , 1–55 (2021)
RHD of non-perfect fluids r [ M ]10 T [ K ] inviscid high ζ low τ Π t [10 M ] − . − . − . − . − . − . − . − . l og ( || δ t || ) Π ρu Figure 3.
Left:
Temperature T as a function of the circumference radius r in units of M and at time , M . Shown are the inviscid solution(filled circles) and the models with high ζ (solid line) and low τ (dashed line), respectively. Note that deviations from the inviscid solution increasetowards the event horizon and that the sonic point is at r s = 200 M . Right: L -norm of the relative time variation δ t for the bulk-viscosity pressure Π (stars), the rest-mass density ρ (crosses), and the primitive fluid velocity u (filled circles) shown as a function of time for the medium ζ modelset using , grid cells. shown in the left panel of Fig. 3. To verify that our calculation reaches a stationary state, we showin the right panel of Fig. 3 the logarithm of the L -norm of the relative time variation δ t for the rest-mass density ρ , the primitive fluid velocity u and the bulk-viscosity pressure Π as a function of timefor the medium ζ model using 20,000 grid cells. For each quantity φ , the relative time variation isdefined as δ t φ ( t ) := 1 − φ ( t − M ) /φ ( t ) . In essence, the right panel of 3 shows that, although thefluid was not stationary in the beginning of the evolution, it reaches an approximately stationarystate for late times. The small but nonzero value of the relative differences at late times is due tolow-amplitude, small-scale oscillations generated at the outer boundary of the numerical domain.We note that after performing a self-convergence test we were able to recover the correct globalconvergence order of BHAC , i.e., ≈ (see Porth et al. 2017). Local self-convergence tests showthat the convergence order is very close to two for small radii, while it exhibits small oscillationsaround two for large radii, with amplitude that increases towards the outer boundary.The left panel of Fig. 4 displays the radial profiles of the temperature for the various viscousmodels when compared to the values obtained from the inviscid solution. As can be seen fromthe solid lines – which refer to the low, medium and high ζ models – the viscous fluid is hotternear the horizon and colder at larger radii than the corresponding inviscid fluid. In particular, forthe high ζ model the temperature at the horizon can be up to ∼ larger than in the inviscidmodel. A similar behaviour can be seen also for the dashed lines, which refer to the low and high τ Π models. Shown instead in the right panel of Fig. 4 is the corresponding comparison in the MNRAS000
Temperature T as a function of the circumference radius r in units of M and at time , M . Shown are the inviscid solution(filled circles) and the models with high ζ (solid line) and low τ (dashed line), respectively. Note that deviations from the inviscid solution increasetowards the event horizon and that the sonic point is at r s = 200 M . Right: L -norm of the relative time variation δ t for the bulk-viscosity pressure Π (stars), the rest-mass density ρ (crosses), and the primitive fluid velocity u (filled circles) shown as a function of time for the medium ζ modelset using , grid cells. shown in the left panel of Fig. 3. To verify that our calculation reaches a stationary state, we showin the right panel of Fig. 3 the logarithm of the L -norm of the relative time variation δ t for the rest-mass density ρ , the primitive fluid velocity u and the bulk-viscosity pressure Π as a function of timefor the medium ζ model using 20,000 grid cells. For each quantity φ , the relative time variation isdefined as δ t φ ( t ) := 1 − φ ( t − M ) /φ ( t ) . In essence, the right panel of 3 shows that, although thefluid was not stationary in the beginning of the evolution, it reaches an approximately stationarystate for late times. The small but nonzero value of the relative differences at late times is due tolow-amplitude, small-scale oscillations generated at the outer boundary of the numerical domain.We note that after performing a self-convergence test we were able to recover the correct globalconvergence order of BHAC , i.e., ≈ (see Porth et al. 2017). Local self-convergence tests showthat the convergence order is very close to two for small radii, while it exhibits small oscillationsaround two for large radii, with amplitude that increases towards the outer boundary.The left panel of Fig. 4 displays the radial profiles of the temperature for the various viscousmodels when compared to the values obtained from the inviscid solution. As can be seen fromthe solid lines – which refer to the low, medium and high ζ models – the viscous fluid is hotternear the horizon and colder at larger radii than the corresponding inviscid fluid. In particular, forthe high ζ model the temperature at the horizon can be up to ∼ larger than in the inviscidmodel. A similar behaviour can be seen also for the dashed lines, which refer to the low and high τ Π models. Shown instead in the right panel of Fig. 4 is the corresponding comparison in the MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke r [ M ] − − − − − − − − − l og ( | − T / T i n v i s c i d | ) low ζ med . ζ high ζ low τ Π high τ Π r [ M ] − − − − − − l og ( | − c s , t / c s | ) . . . . . c s , t [ c ] inviscid Figure 4.
Left:
Relative difference between the temperature T of the viscous models and the corresponding inviscid model. Note that the differencesclearly increase with the increase of the bulk viscosity and that the relative difference changes sign somewhere outside the event horizon. Right:
Same as on the left but for the sound speed. Reported in the inset is the actual value of the total sound speed near the event horizon, which showsconsiderable deviations from the inviscid solution. case of viscous speed of sound when analysed very close to the event horizon. The excess can beup to a factor of almost three for the low τ Π model and a factor of ∼ . for the high ζ model,reaching viscous sound speeds above . and . , respectively (see inset). In general, and as canbe intuitively expected, solutions with high bulk viscosities and low relaxation times tend to havelarger temperatures and viscous sound speeds, and hence larger deviations from the inviscid case.Figure 5 shows radial profiles of the relativistic “inverse Reynolds number” Π /ρh (left panel,solid and dashed lines) and of the bulk-viscosity pressure normalized to the corresponding NSvalue Π / Π NS (right panel, solid and dashed lines). Note that for all models Π /ρh assumes smallbut finite values at large radii, decreases when moving inwards and then increases again sharplyclose to the horizon. The overall magnitude of Π /ρh is very sensitive to the parameter ζ (cf. low , medium , and high ζ models), while τ seems to play a marginal role. On the other hand, τ controlsthe location of the sharp increase (cf. cases medium ζ , low τ Π , and high τ Π ). Note that most modelsapproach their corresponding NS values near the horizon, while the rate at which this happens isagain controlled by τ . This can be seen in the right panel of Fig. 5, where models with different τ show different behaviour. In particular, for the models low ζ , medium ζ , and high ζ , the bulk-viscosity pressure reaches nearly ∼ of the NS value, while the corresponding value for the high τ Π model is considerably smaller and of the order of ∼ . Note also that the low τ Π modelreaches a maximum of ∼ , but not exactly at the horizon; this is most likely a behaviour due MNRAS , 1–55 (2021)
RHD of non-perfect fluids r [ M ] − − − − − − − l og ( Π / ρ h ) r [ M ] − − l og ( Π / Π N S ) low ζ med . ζ high ζ low τ Π high τ Π − − − . . . . Π / Π N S Figure 5.
Left: radial profiles of the ratio of the bulk-viscosity pressure over the enthalpy density (i.e., inverse Reynolds number). Solid lines ofdifferent colours refer to models with low ζ , medium ζ , and high ζ ; dashed lines of different colours refer to models with low τ Π and high τ Π ,respectively. The inset reports the same quantities but near the event horizon. Right: radial profiles of the ratio of the bulk-viscosity pressure overthe corresponding NS value. Solid and dashed lines follow the same convention as in the left panel, and inset reports the same quantities but nearthe event horizon. to a cancellation error near the horizon for the solution obtained by
BHAC , since it is absent in theinitial solution.As a concluding remark we note that there are analogies between the late-time behaviour re-alised in longitudinally expanding fluids, such as the Bjorken flow, and the near-horizon propertiesof the accretion solution considered here. In both cases, in fact, the solution tends to the corre-sponding NS value (for late times in the case of the Bjorken flow and for r ∼ M in the caseof accretion). This behaviour suggests that while the parameter ζ controls the magnitude of first-order non-equilibrium effects – which in the case of the accretion develop mostly in strong-gravityregions – τ controls the degree to which Π approaches its NS value. Of course, these consid-erations are based on the examination of the simplest form of a bulk viscosity. A more extensiveinvestigation of possible initial conditions, transport coefficients and additional source terms in thebulk-viscosity pressure equation will yield a deeper insight into the accretion process of viscousmatter that, as pointed out by Turolla & Nobili (1989), is far from being trivial. After having reviewed the various approaches developed over the years to model relativistic dis-sipative fluids, we have derived a general-relativistic 3+1 flux-conservative formulation of the
MNRAS000
MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke second-order dissipative-hydrodynamics equations first suggested by Israel (1976) and Hiscock& Lindblom (1983) (i.e., the
HL83 set of GRDHD equations). The new set of equations derivedin this way provides a comprehensive and complete way of including causal dissipative effects ingeneral-relativistic calculations.Although a 3+1 formulation of a reduced version of the
HL83 set of GRDHD equations wasalready proposed by Peitz & Appl (1997) and Peitz & Appl (1999) , the work presented hereextends the results of Peitz and Appl in three important ways. First, our set of equations is completeand does not neglect terms that are considered to play a less significant role (see Appendix Afor a comparison of the system presented here with other formulations). Second, the equationspresented are cast into a flux-conservative form suitable for numerical implementation. Finally,also the coupling terms between the different dissipative currents are rewritten a 3+1 form. As aresult, the full system given in Hiscock & Lindblom (1983) can now be readily implemented inmodern numerical-relativity codes and evolved numerically.As a way to test the new set of equations, we have proceeded with the implementation in theGRMHD code BHAC of a reduced version of the equations describing fluids with zero shear andheat currents, and used them against a number of tests in flat and curved background metrics.In the first case, we have first considered the one-dimensional, longitudinally boost-invariantmotion of a viscous fluid, such as the one produced in an ultrarelativistic collision of two ions. Thisexpansion, first proposed by Bjorken for an inviscid fluid, is a standard testbed and can be solvedanalytically. Overall, we find that the numerical solutions agree well with the analytical solutions,with a relative difference that is ∼ − for all the values of viscosity considered. As an additionalflat-spacetime test, we have explored the solution of a shock-tube problem for an ultra-relativisticgas of gluons for different values of the ratio ζ/s , with ζ and s being bulk-viscosity coefficientand the entropy density, respectively. Also in this case, when comparing the solutions with thoseobtained with methods based on the direct solution of the relativistic Boltzmann equation, wefind a very good agreement up to a ratio of ζ/s (cid:46) . , which is already in a regime where adissipative-hydrodynamics framework breaks down.Finally, as a general-relativistic test we have considered for the first time the problem of sta-tionary, spherically symmetric accretion of bulk viscous fluids onto nonrotating black holes withina second-order dissipative-hydrodynamics framework. Starting from initial conditions obtained In Peitz & Appl (1997) and Peitz & Appl (1999), all terms including products of dissipative currents, i.e., Π , q and π , and first-order gradientsof the primary fluid variables are set to zero in Eqs. (8) – (10). MNRAS , 1–55 (2021) RHD of non-perfect fluids
BHAC can deviate from the corresponding inviscid solution withdifferences (cid:46) for the temperature and ∼ for the sound speed. In addition, the bulk-viscosity pressure is highly sensitive to the bulk-viscosity coefficient, exhibiting deviations of upto three orders of magnitude near the event horizon depending on whether the viscosity is large orsmall. We also show that, although BHAC is able to maintain a quasi-stationarity in the solutionand shows global second-order convergence, it does not converge to the reference solution withincreasing grid resolution because of the influence of a finite-size computational domain. Interest-ingly, we note analogies between the late-time behaviour in the Bjorken and accretion flows. Inboth cases, in fact, the solution tends to the corresponding NS value, with this happening at latetimes in the Bjorken flow, and near the horizon in the case of accretion onto a black hole.As a concluding remark we note that although the 3+1 formulation presented here, and itscorresponding discretisation, offers a viable path to the inclusion of non-equilibrium effects ingeneral-relativistic simulations of compact objects, short relaxation times as well as the requiredtemporal and spatial discretization of the source terms in the full system may lead to stiff equations,whose solution will not be feasible with simple explicit schemes. We leave the examination ofmore sophisticated numerical techniques for future work, where mixed explicit-implicit (IMEX)time integrators (see, e.g., Palenzuela et al. 2009; Dionysopoulou et al. 2013; Weih et al. 2020)will be considered.
ACKNOWLEDGEMENTS
We thank Masoud Shokri, Elias Most, Hector Olivares and Lukas Weih for useful discussions. Sup-port comes in part from HGS-HIRe for FAIR; the LOEWE-Program in HIC for FAIR; “PHAROS”,COST Action CA16214; the ERC Synergy Grant “BlackHoleCam: Imaging the Event Horizonof Black Holes” (Grant No. 610058); the Deutsche Forschungsgemeinschaft (DFG, German Re-search Foundation) through the CRC-TR 211 “Strong-interaction matter under extreme condi-tions” - project number 315477589 - TRR 211.
MNRAS000
MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
References
Abbott B. P., et al., 2017a, Phys. Rev. Lett., 119, 161101Abbott B. P., et al., 2017b, Astrophys. J. Lett., 848, L13Abbott B. P., et al., 2018, Physical Review Letters, 121, 161101Afshordi N., Paczy´nski B., 2003, Astrophys. J., 592, 354Aguilera-Miret R., Viganò D., Carrasco F., Miñano B., Palenzuela C., 2020, Phys. Rev. D, 102, 103006Alcubierre M., 2008, Introduction to
Numerical Relativity. Oxford University Press, Oxford, UK,doi:10.1093/acprof:oso/9780199205677.001.0001Alford M. G., Haber A., 2020, arXiv e-prints, p. arXiv:2009.05181Alford M. G., Bovard L., Hanauske M., Rezzolla L., Schwenzer K., 2018, Phys. Rev. Lett., 120, 041101Alford M., Harutyunyan A., Sedrakian A., 2020, Particles, p. arXiv:2006.07975Andersson N., Comer G. L., 2015, Classical and Quantum Gravity, 32, 075008Andersson N., Comer G. L. C., 2020, arXiv e-prints, p. arXiv:2008.12069Annala E., Gorda T., Kurkela A., Vuorinen A., 2018, Phys. Rev. Lett., 120, 172703Baier R., Romatschke P., Son D. T., Starinets A. O., Stephanov M. A., 2008, Journal of High Energy Physics, 2008, 100Baiotti L., Rezzolla L., 2017, Rept. Prog. Phys., 80, 096901Baiotti L., Giacomazzo B., Rezzolla L., 2008, Phys. Rev. D, 78, 084033Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, Astrophys. J. Lett., 850, L34Bemfica F. S., Disconzi M. M., Noronha J., 2019a, Phys. Rev. D, 100, 104020Bemfica F. S., Disconzi M. M., Noronha J., 2019b, Phys. Rev. Lett., 122, 221602Bemfica F. S., Disconzi M. M., Hoang V., Noronha J., Radosz M., 2020, arXiv e-prints, p. arXiv:2005.11632Bernhard J. E., Moreland J. S., Bass S. A., 2019, Nature Physics, 15, 1113Betz B., Henkel D., Rischke D. H., 2009, Progress in Particle and Nuclear Physics, 62, 556Betz B., Denicol G., Koide T., Molnar E., Niemi H., Rischke D., 2011, EPJ Web Conf., 13, 07005Biswas R., Dash A., Haque N., Pu S., Roy V., 2020, Journal of High Energy Physics, 2020, 171Bjorken J. D., 1983, Phys. Rev. D, 27, 140Bondi H., 1952, Mon. Not. R. Astron. Soc., 112, 195Bouras I., Molnar E., Niemi H., Xu Z., El A., Fochler O., Greiner C., Rischke D. H., 2009a, Phys. Rev. Lett., 103, 032301Bouras I., Molnár E., Niemi H., Xu Z., El A., Fochler O., Greiner C., Rischke D. H., 2009b, Nuclear Phys. A, 830, 741Bouras I., Molnar E., Niemi H., Xu Z., El A., Fochler O., Greiner C., Rischke D. H., 2010, Phys. Rev., C82, 024910Brito C. V., Denicol G. S., 2020, Phys. Rev. D, 102, 116009Burns E., 2020, Living Reviews in Relativity, 23, 4Busza W., Rajagopal K., van der Schee W., 2018, Annual Review of Nuclear and Particle Science, 68, 339Carter B., 1989, in Anile A., Choquet-Bruhat Y., eds, Lecture Notes in Mathematics, Vol. 1385, Relativistic Fluid Dynamics. Springer Berlin /Heidelberg, pp 1–64, http://dx.doi.org/10.1007/BFb0084028
Cowperthwaite P. S., et al., 2017, Astrophys. J. Lett., 848, L17De S., Finstad D., Lattimer J. M., Brown D. A., Berger E., Biwer C. M., 2018, Physical Review Letters, 121, 091102Del Zanna L., et al., 2013, European Physical Journal C, 73, 2524Denicol G. S., Kodama T., Koide T., Mota P., 2008a, Journal of Physics G: Nuclear and Particle Physics, 35, 115102Denicol G. S., Kodama T., Koide T., Mota P., 2008b, Phys. Rev. C, 78, 034901Denicol G. S., Molnár E., Niemi H., Rischke D. H., 2012a, The European Physical Journal A, 48, 170Denicol G. S., Niemi H., Molnár E., Rischke D. H., 2012b, Phys. Rev. D, 85, 114047Denicol G., Niemi H., Bouras I., Molnar E., Xu Z., Rischke D., Greiner C., 2014, Phys. Rev. D, 89, 074005Dionysopoulou K., Alic D., Palenzuela C., Rezzolla L., Giacomazzo B., 2013, Phys. Rev. D, 88, 044020Disconzi M. M., Kephart T. W., Scherrer R. J., 2017, International Journal of Modern Physics D, 26, 1750146 MNRAS , 1–55 (2021)
RHD of non-perfect fluids Drout M. R., et al., 2017, Science, 358, 1570Duez M. D., Liu Y. T., Shapiro S. L., Stephens B. C., 2004, Phys. Rev. D, 69, 104030Duez M. D., et al., 2020, Phys. Rev. D, 102, 104050Eckart C., 1940, Phys. Rev., 58, 919Event Horizon Telescope Collaboration et al., 2019, Astrophys. J. Lett., 875, L1Fambri F., Dumbser M., Köppel S., Rezzolla L., Zanotti O., 2018, Mon. Not. R. Astron. Soc., 477, 4543Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018, The Astrophysical Journal, 860, 64Gabbana A., Plumari S., Galesi G., Greco V., Simeoni D., Succi S., Tripiccione R., 2020, Phys. Rev. C, 101, 064904Gourgoulhon E., 2012, 3+1 Formalism in General Relativity. Lecture Notes in Physics, Berlin Springer Verlag Vol. 846, Springer, doi:10.1007/978-3-642-24525-1Hairer E., Wanner G., 1996, Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems, 2nd edn. SpringerHawley J. F., Smarr L. L., Wilson J. R., 1984, Astrophys. J., 277, 296Hebert F., Kidder L. E., Teukolsky S. A., 2018, Phys. Rev. D, 98, 044041Hiscock W. A., Lindblom L., 1983, Annals of Physics, 151, 466Hiscock W. A., Lindblom L., 1985, Phys. Rev. D, 31, 725Hiscock W. A., Lindblom L., 1987, Phys. Rev. D, 35, 3723Hoult R. E., Kovtun P., 2020, Journal of High Energy Physics, 2020, 67Inghirami G., Del Zanna L., Beraudo A., Moghaddam M. H., Becattini F., Bleicher M., 2016, European Physical Journal C, 76, 659Inghirami G., Del Zanna L., Beraudo A., Haddadi Moghaddam M., Becattini F., Bleicher M., 2018, in Journal of Physics Conference Series. p.012043, doi:10.1088/1742-6596/1024/1/012043Israel W., 1976, Annals of Physics, 100, 310Israel W., Stewart J. M., 1979, Annals of Physics, 118, 341Jaiswal A., 2013, Phys. Rev. C, 87, 051901Jaiswal A., Bhalerao R. S., Pal S., 2013, Phys. Rev. C, 87, 021901Jordan D. W., Smith P., 2007, Nonlinear Ordinary Differential Equations : An Introduction for Scientists and Engineers. Oxford University Press,Oxford, UKKiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2018, Phys. Rev. D, 97, 124039Koeppel S., Bovard L., Rezzolla L., 2019, Astrophys. J. Lett., 872, L16Kovtun P., 2019, Journal of High Energy Physics, 2019, 34Landau L. D., Lifshitz E. M., 2004, Fluid Mechanics, Course of Theoretical Physics, Volume 6. Elsevier Butterworth-Heinemann, OxfordMalik T., Alam N., Fortin M., Providência C., Agrawal B. K., Jha T. K., Kumar B., Patra S. K., 2018, Physical Review C, 98, 035804Mandal I., Ray A. K., Das T. K., 2007, Mon. Not. R. Astron. Soc., 378, 1400Margalit B., Metzger B. D., 2017, Astrophys. J. Lett., 850, L19McNelis M., Bazow D., Heinz U., 2021, arXiv e-prints, p. arXiv:2101.02827Michel F. C., 1972, Ap&SS, 15, 153Mihalas D., Mihalas B., 1984, Foundations of radiation hydrodynamicsMizuno Y., et al., 2018, Nature Astronomy,Montaña G., Tolós L., Hanauske M., Rezzolla L., 2019, Phys. Rev. D, 99, 103009Most E. R., Weih L. R., Rezzolla L., Schaffner-Bielich J., 2018, Phys. Rev. Lett., 120, 261103Most E. R., Papenfort L. J., Rezzolla L., 2019, Mon. Not. R. Astron. Soc., 490, 3588Muronga A., 2004, Phys. Rev. C, 69, 034903Nobili L., Turolla R., Zampieri L., 1991, Astrophys. J., 383, 250Olivares H., et al., 2020, MNRAS, 497, 521Palenzuela C., Lehner L., Reula O., Rezzolla L., 2009, Mon. Not. R. Astron. Soc., 394, 1727Paschalidis V., 2017, Classical and Quantum Gravity, 34, 084002MNRAS000
RHD of non-perfect fluids Drout M. R., et al., 2017, Science, 358, 1570Duez M. D., Liu Y. T., Shapiro S. L., Stephens B. C., 2004, Phys. Rev. D, 69, 104030Duez M. D., et al., 2020, Phys. Rev. D, 102, 104050Eckart C., 1940, Phys. Rev., 58, 919Event Horizon Telescope Collaboration et al., 2019, Astrophys. J. Lett., 875, L1Fambri F., Dumbser M., Köppel S., Rezzolla L., Zanotti O., 2018, Mon. Not. R. Astron. Soc., 477, 4543Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018, The Astrophysical Journal, 860, 64Gabbana A., Plumari S., Galesi G., Greco V., Simeoni D., Succi S., Tripiccione R., 2020, Phys. Rev. C, 101, 064904Gourgoulhon E., 2012, 3+1 Formalism in General Relativity. Lecture Notes in Physics, Berlin Springer Verlag Vol. 846, Springer, doi:10.1007/978-3-642-24525-1Hairer E., Wanner G., 1996, Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems, 2nd edn. SpringerHawley J. F., Smarr L. L., Wilson J. R., 1984, Astrophys. J., 277, 296Hebert F., Kidder L. E., Teukolsky S. A., 2018, Phys. Rev. D, 98, 044041Hiscock W. A., Lindblom L., 1983, Annals of Physics, 151, 466Hiscock W. A., Lindblom L., 1985, Phys. Rev. D, 31, 725Hiscock W. A., Lindblom L., 1987, Phys. Rev. D, 35, 3723Hoult R. E., Kovtun P., 2020, Journal of High Energy Physics, 2020, 67Inghirami G., Del Zanna L., Beraudo A., Moghaddam M. H., Becattini F., Bleicher M., 2016, European Physical Journal C, 76, 659Inghirami G., Del Zanna L., Beraudo A., Haddadi Moghaddam M., Becattini F., Bleicher M., 2018, in Journal of Physics Conference Series. p.012043, doi:10.1088/1742-6596/1024/1/012043Israel W., 1976, Annals of Physics, 100, 310Israel W., Stewart J. M., 1979, Annals of Physics, 118, 341Jaiswal A., 2013, Phys. Rev. C, 87, 051901Jaiswal A., Bhalerao R. S., Pal S., 2013, Phys. Rev. C, 87, 021901Jordan D. W., Smith P., 2007, Nonlinear Ordinary Differential Equations : An Introduction for Scientists and Engineers. Oxford University Press,Oxford, UKKiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2018, Phys. Rev. D, 97, 124039Koeppel S., Bovard L., Rezzolla L., 2019, Astrophys. J. Lett., 872, L16Kovtun P., 2019, Journal of High Energy Physics, 2019, 34Landau L. D., Lifshitz E. M., 2004, Fluid Mechanics, Course of Theoretical Physics, Volume 6. Elsevier Butterworth-Heinemann, OxfordMalik T., Alam N., Fortin M., Providência C., Agrawal B. K., Jha T. K., Kumar B., Patra S. K., 2018, Physical Review C, 98, 035804Mandal I., Ray A. K., Das T. K., 2007, Mon. Not. R. Astron. Soc., 378, 1400Margalit B., Metzger B. D., 2017, Astrophys. J. Lett., 850, L19McNelis M., Bazow D., Heinz U., 2021, arXiv e-prints, p. arXiv:2101.02827Michel F. C., 1972, Ap&SS, 15, 153Mihalas D., Mihalas B., 1984, Foundations of radiation hydrodynamicsMizuno Y., et al., 2018, Nature Astronomy,Montaña G., Tolós L., Hanauske M., Rezzolla L., 2019, Phys. Rev. D, 99, 103009Most E. R., Weih L. R., Rezzolla L., Schaffner-Bielich J., 2018, Phys. Rev. Lett., 120, 261103Most E. R., Papenfort L. J., Rezzolla L., 2019, Mon. Not. R. Astron. Soc., 490, 3588Muronga A., 2004, Phys. Rev. C, 69, 034903Nobili L., Turolla R., Zampieri L., 1991, Astrophys. J., 383, 250Olivares H., et al., 2020, MNRAS, 497, 521Palenzuela C., Lehner L., Reula O., Rezzolla L., 2009, Mon. Not. R. Astron. Soc., 394, 1727Paschalidis V., 2017, Classical and Quantum Gravity, 34, 084002MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
Peitz J., Appl S., 1997, Mon. Not. R. Astron. Soc., 286, 681Peitz J., Appl S., 1999, Class. Quantum Grav., 16, 979Porth O., Olivares H., Mizuno Y., Younsi Z., Rezzolla L., Moscibrodzka M., Falcke H., Kramer M., 2017, Computational Astrophysics andCosmology, 4, 1Pu S., Koide T., Rischke D. H., 2010, Phys. Rev. D, 81, 114039Pu S., Roy V., Rezzolla L., Rischke D. H., 2016, Phys. Rev. D, 93, 074022Radice D., 2017, Astrophys. J. Lett., 838, L2Radice D., 2020, Symmetry, 12, 1249Radice D., Rezzolla L., 2012, Astron. Astrophys., 547, A26Radice D., Perego A., Bernuzzi S., Zhang B., 2018a, Mon. Not. R. Astron. Soc., 481, 3670Radice D., Perego A., Zappa F., Bernuzzi S., 2018b, Astrophys. J. Lett., 852, L29Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., Roberts L. F., 2018c, Astrophys. J., 869, 130Raithel C., Özel F., Psaltis D., 2018, Astrophys. J., 857, L23Ray A. K., Bhattacharjee J. K., 2002, Phys. Rev. E, 66, 066303Rezzolla L., Zanotti O., 2013, Relativistic Hydrodynamics. Oxford University Press, Oxford, UK,doi:10.1093/acprof:oso/9780198528906.001.0001Rezzolla L., Most E. R., Weih L. R., 2018, Astrophys. J. Lett., 852, L25Romatschke P., Romatschke U., 2019, Relativistic Fluid Dynamics In and Out of Equilibrium: And Applications to Relativistic Nuclear Collisions.Cambridge Monographs on Mathematical Physics, Cambridge University Press, doi:10.1017/9781108651998Roy V., Pu S., Rezzolla L., Rischke D., 2015, Physics Letters B, 750, 45Ruiz M., Shapiro S. L., Tsokaros A., 2018, Phys. Rev. D, 97, 021501Rusanov V. V., 1961, J. Comput. Math. Phys. USSR, 1, 267Sawyer R. F., 2006, Phys. Rev. D, 74, 043527Shibata M., Hotokezaka K., 2019, Annual Review of Nuclear and Particle Science, 69, 41Shibata M., Kiuchi K., Sekiguchi Y.-i., 2017, Phys. Rev. D, 95, 083005Shibata M., Zhou E., Kiuchi K., Fujibayashi S., 2019, Phys. Rev. D, 100, 023015Siegel D. M., Ciolfi R., Harte A. I., Rezzolla L., 2013, Phys. Rev. D R, 87, 121302Taghinavaz F., 2020, JHEP, 08, 119Tews I., Margueron J., Reddy S., 2018, Physical Review C, 98, 045804Turolla R., Nobili L., 1989, Astrophys. J., 342, 982Van P., Biro T., 2012, Phys. Lett. B, 709, 106Viganò D., Aguilera-Miret R., Carrasco F., Miñano B., Palenzuela C., 2020, Phys. Rev. D, 101, 123019Weih L. R., Olivares H., Rezzolla L., 2020, Mon. Not. R. Astron. Soc., 495, 2285Weinberg S., 1971, Astrophys. J., 168, 175Wolfram Research I., 2020, Mathematica, Version 12.1,
APPENDIX A: REVIEW OF SECOND-ORDER RELATIVISTIC DISSIPATIVEHYDRODYNAMICS
This Appendix serves the scope of providing a guidance in comparing our reference second-orderformulation of general-relativistic dissipative hydrodynamics originally suggested by Hiscock &Lindblom (1983) with other formulations that have appeared in the literature, i.e., the formulations
MNRAS , 1–55 (2021)
RHD of non-perfect fluids
IS79 ), by Denicol et al. (2012a) (hereafter
DMNR12 ) , andby Baier et al. (2008) (hereafter rBRSSS08 ). Since these formulations often adopt different no-tations that makes it hard to compare them, we introduce a generalized notation for the varioustransport coefficients, whose terminology is given in Table A1. Furthermore, since the equationsof HL83 are chosen as our reference equations, we use upper-case letters for transport coefficientsappearing in
HL83 ; the only exception to this rule will be made for the relaxation times. Further-more, since the comparison needs to distinguish terms of first and second order, we introduce thefollowing dimensionless numbers: the Knudsen number
Kn := l micro L macro , (A1)where l micro and L macro are the microscopic lengthscale given by the mean-free-path of the mi-croscopic constituents of the fluid and the macroscopic lengthscale given by the distance overwhich gradients of primary fluid variables appear, respectively. Furthermore, the inverse Reynoldsnumbers are given by R − := | Π | p + e , R − n := | n µ | n , R − π := | π µν | p + e . (A2)Note that in the Eckart frame the inverse Reynolds number R − n is replaced by R − q := | q µ | p + e ∼ O (cid:0) R − n (cid:1) . (A3)Using these dimensionless numbers, we can identify the order of the terms given in Table A2according to the following classification: O (cid:0) R − i , Kn (cid:1) := O , O (cid:0) R − i Kn (cid:1) := O RK , O (cid:0) R − i R − j (cid:1) := O , O (cid:0) Kn (cid:1) := O . (A4)In this way, following the convention of Denicol et al. (2012b), we refer to terms of order O asbeing of first order, while all other terms are referred to as being of second order.With these definitions made, we next proceed to the actual comparison which will take placeby first briefly reviewing each of the formulations considered (Sec. A1–A4) and then proceedswith the actual comparison (Sec. A5). Note that there is another frequently used second-order GRDHD formulation proposed by Denicol et al. (2012b) and often referred to as the“DNMR” formulation. We will not use this formulation for our comparison here since the formulation by Denicol et al. (2012b) is in the Landauframe and we instead consider the formulations adopting the Eckart frame for our comparison.MNRAS000
Kn := l micro L macro , (A1)where l micro and L macro are the microscopic lengthscale given by the mean-free-path of the mi-croscopic constituents of the fluid and the macroscopic lengthscale given by the distance overwhich gradients of primary fluid variables appear, respectively. Furthermore, the inverse Reynoldsnumbers are given by R − := | Π | p + e , R − n := | n µ | n , R − π := | π µν | p + e . (A2)Note that in the Eckart frame the inverse Reynolds number R − n is replaced by R − q := | q µ | p + e ∼ O (cid:0) R − n (cid:1) . (A3)Using these dimensionless numbers, we can identify the order of the terms given in Table A2according to the following classification: O (cid:0) R − i , Kn (cid:1) := O , O (cid:0) R − i Kn (cid:1) := O RK , O (cid:0) R − i R − j (cid:1) := O , O (cid:0) Kn (cid:1) := O . (A4)In this way, following the convention of Denicol et al. (2012b), we refer to terms of order O asbeing of first order, while all other terms are referred to as being of second order.With these definitions made, we next proceed to the actual comparison which will take placeby first briefly reviewing each of the formulations considered (Sec. A1–A4) and then proceedswith the actual comparison (Sec. A5). Note that there is another frequently used second-order GRDHD formulation proposed by Denicol et al. (2012b) and often referred to as the“DNMR” formulation. We will not use this formulation for our comparison here since the formulation by Denicol et al. (2012b) is in the Landauframe and we instead consider the formulations adopting the Eckart frame for our comparison.MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke transport coefficients description τ Π , τ q , τ π relaxation times for bulk viscosity, heat conduction and shear viscosity, respectively δ ij or ∆ ij transport coefficients for gradients of the fluid velocity u l ij or L ij transport coefficients for gradients of the dissipative currents λ ij or Λ ij transport coefficients for the contribution of gradients of ρ and pϕ ij transport coefficients for the contribution of contracted dissipative currents g ij transport coefficients in front of geometrical quantities, i.e. R µνλρ , R µν or R Table A1.
Generalized notation for the transport coefficients. Shown on the left are the newly introduced symbols, while a description of whichfluid field it appears in combination with is shown on the right. The upper index indicates which dissipative current the transport coefficient belongsto, e.g., δ Π is the second transport coefficient expressing the coupling to a gradient of the fluid velocity in the constitutive equation for the bulk-viscosity pressure. Upper-case letters indicate that the corresponding transport coefficient occurs in the HL83 formulation, while lower case lettersindicate that it is set to zero in
HL83 . A1 Israel and Stewart 1979 (IS79)
The equations for the dissipative currents of
IS79 read as follows τ Π ˙Π = Π LNS − Π + δ Π q µ a µ + L Π ∇ µ q µ , (A5) τ q ˙ q (cid:104) µ (cid:105) = q µ NS − q µ + δ q Π a µ + δ q π µν a ν + δ q ω µν q ν + L q ∇ (cid:104) µ (cid:105) Π + L q ∇ ν π (cid:104) µ (cid:105) ν , (A6) τ π ˙ π (cid:104) µν (cid:105) = π µν LNS − π µν + δ π q (cid:104) µ a ν (cid:105) + δ π π λ (cid:104) µ ω ν (cid:105) λ + L π ∇ (cid:104) µ q ν (cid:105) . (A7)Here, Π LNS and π µν LNS denote the NS values of the dissipative currents evaluated in the Landauframe. Equations (A5), (A6) and (A7) correspond to equations (7.1a), (7.1b) and (7.1c) of Israel& Stewart (1979), respectively. This set of equations was derived using the Boltzmann equationfor the single-particle distribution function f , parametrized in the form y ( x µ , p i ) := ln (cid:20) f ( x µ , p i )A (cid:21) , (A8) y ( x µ , p i ) := µT + (cid:15) + (cid:16) u λ T + (cid:15) λ m (cid:17) p λ + 1 m (cid:15) λρ p λ p ρ , (A9) A := 1 + A f (cid:0) x µ , p i (cid:1) , A ∈ {− , , +1 } , (A10)where p = p ( p i , m ) is the on-shell energy of the particles, A = − refers to fermions(bosons) so that A → corresponds to the limit yielding the Boltzmann distribution, and m isthe mass of the particles.The functions (cid:15) , (cid:15) λ and (cid:15) λρ are off-equilibrium corrections. Note that we recover the standardFermi, Bose and Boltzmann local-equilibrium distributions for (cid:15), (cid:15) λ , (cid:15) λρ = 0 . The single-particledistribution f contains 14 independent variables ( T, µ , the three independent components of u µ andthe nine independent components of the functions (cid:15) , (cid:15) λ and (cid:15) λρ ). These are one-to-one matchedto the components of J µ and T µν . The system of conservation equations, corresponding to theequations of motion for the first and second moment of f , is closed by employing the equation ofmotion for the third moment of f , leading to relaxation-type equations for the dissipative currents. MNRAS , 1–55 (2021)
RHD of non-perfect fluids “14-moment approximation” .We should recall that in the derivation of
IS79 two approximations are made. First, they ne-glected second-order terms proportional to ∇ µ ρ and Θ . Following the derivation from the secondlaw of thermodynamics, these terms arise from gradients of the transport coefficients and are gen-erally of the form (see, e.g., the equation for the shear-stress tensor) ∝ π µν u λ I λ and ∝ q (cid:104) µ I ν (cid:105) ,where I µ is a linear combination of gradients of transport coefficients that can be related to gradi-ents of the chosen pair of thermodynamical variables { T, µ } or { ρ, p } . In general, these terms canbe further decomposed by using the first-order version of the conservation laws (1) and (2): ˙ ρ ∝ Θ , (A11) a µ ∝ ∇ (cid:104) µ (cid:105) p , (A12) ˙ e ∝ Θ , (A13)where, again, ˙ A := ( u · ∇ ) A = u µ ∇ µ A . By choosing the pair { ρ, p } as our thermodynamicalvariables and exploiting Eqs. (A11) – (A13), the following relations are valid to second-order π µν u λ I λ ∝ π µν Θ , (A14) q (cid:104) µ I ν (cid:105) ∝ b q (cid:104) µ a ν (cid:105) + b q (cid:104) µ ∇ ν (cid:105) ρ , (A15)where b and b denote scalar functions of the pair { ρ, p } . Hence, the missing terms in IS79 proportional to ∇ µ ρ and Θ are of the form ∝ π µν Θ and ∝ q (cid:104) µ ∇ ν (cid:105) ρ , while the term b q (cid:104) µ a ν (cid:105) couldbe absorbed in the transport coefficient δ π .Second, there are terms of the form ∝ π λ (cid:104) µ σ ν (cid:105) λ and ∝ Π σ µν , appearing for instance in theequation for the shear-stress tensor that are also apparently missing in the IS79 formulation. Theexistence of these terms was first pointed out by Betz et al. (2011), which were subsequentlyincluded in the
DMNR12 formulation (see below).
A2 Hiscock and Lindblom 1983 (HL83)
When using a generalized notation, the evolution equations (8)–(10) for the dissipative currents inthe
HL83 formulation can be written as τ Π ˙Π = Π NS − Π + ∆ Π ΠΘ + L Π ∇ µ q µ + Λ Π Π u µ I µ + Λ Π q µ I µ , (A16) τ q ˙ q (cid:104) µ (cid:105) = q µ NS − q µ + ∆ q q µ Θ + L q ∇ (cid:104) µ (cid:105) Π + L q ∇ ν π (cid:104) µ (cid:105) ν + Λ q q µ u ν I ν + Λ q Π I (cid:104) µ (cid:105) + Λ q π µν I ν , (A17) τ π ˙ π (cid:104) µν (cid:105) = π µν NS − π µν + ∆ π π µν Θ + L π ∇ (cid:104) µ q ν (cid:105) + Λ π π µν u λ I λ + Λ π q (cid:104) µ I ν (cid:105) . (A18) MNRAS000
HL83 formulation can be written as τ Π ˙Π = Π NS − Π + ∆ Π ΠΘ + L Π ∇ µ q µ + Λ Π Π u µ I µ + Λ Π q µ I µ , (A16) τ q ˙ q (cid:104) µ (cid:105) = q µ NS − q µ + ∆ q q µ Θ + L q ∇ (cid:104) µ (cid:105) Π + L q ∇ ν π (cid:104) µ (cid:105) ν + Λ q q µ u ν I ν + Λ q Π I (cid:104) µ (cid:105) + Λ q π µν I ν , (A17) τ π ˙ π (cid:104) µν (cid:105) = π µν NS − π µν + ∆ π π µν Θ + L π ∇ (cid:104) µ q ν (cid:105) + Λ π π µν u λ I λ + Λ π q (cid:104) µ I ν (cid:105) . (A18) MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke where the currents I iµ are given by I µ := ∇ µ ( τ Π /ζT ) , (A19) I µ := ∇ µ ( α /T ) , (A20) I µ := ∇ µ (cid:0) τ q /κT (cid:1) , (A21) I µ := ∇ µ ( α /T ) , (A22) I µ := ∇ µ ( τ π /ηT ) . (A23)In similarity with Sec. A1, we can use the first-order relations (A11)–(A13) to further decom-pose the terms in (A16)–(A18) that include the currents I iµ and obtain Π u µ I µ ∝ ΠΘ , (A24) q µ I µ ∝ b q µ a µ + b q µ ∇ µ ρ , (A25) q µ u ν I ν ∝ q µ Θ , (A26) Π I (cid:104) µ (cid:105) ∝ b Π a µ + b Π ∇ (cid:104) µ (cid:105) ρ , (A27) π µν I ν ∝ b π µν a ν + b π µν ∇ ν ρ , (A28) π µν u λ I λ ∝ π µν Θ , (A29) q (cid:104) µ I ν (cid:105) ∝ b q (cid:104) µ a ν (cid:105) + b q (cid:104) µ ∇ ν (cid:105) ρ , (A30)where, again, the quantities b - b are scalar functions of ρ and p . A3 Denicol et al. 2012a (DMNR12)
The equations for the dissipative currents of the
DMNR12 formulation are τ Π ˙Π = Π NS − Π + ∆ Π ΠΘ + δ Π q µ a µ + δ Π π µν σ µν + L Π ∂ µ q µ + Λ Π q µ I µ , (A31) τ q ˙ q (cid:104) µ (cid:105) = q µ NS − q µ + ∆ q q µ Θ + δ q Π a µ + δ q π µν a ν + δ q ω µν q ν + δ q σ µν q ν + L q ∂ (cid:104) µ (cid:105) Π + L q ∂ ν π (cid:104) µ (cid:105) ν + Λ q Π I (cid:104) µ (cid:105) + Λ q π µν I ν , (A32) τ π ˙ π (cid:104) µν (cid:105) = π µν NS − π µν + ∆ π π µν Θ + δ π q (cid:104) µ a ν (cid:105) + δ π π λ (cid:104) µ ω ν (cid:105) λ + δ π π λ (cid:104) µ σ ν (cid:105) λ + δ π Π σ µν + L π ∂ (cid:104) µ q ν (cid:105) + Λ π q (cid:104) µ I ν (cid:105) . (A33)Equations (A31), (A32) and (A33) correspond to equations (128), (138) and (153) of Denicolet al. (2012a), respectively, when the Eckart frame is chosen as the frame of reference, i.e., when MNRAS , 1–55 (2021)
RHD of non-perfect fluids V µ = 0 and W µ = q µ . Notice that the terms including the current I µ can be decomposed as q µ I µ ∝ b q µ a µ + b q µ ∇ µ ρ , (A34) Π I (cid:104) µ (cid:105) ∝ b Π a µ + b Π ∇ (cid:104) µ (cid:105) ρ , (A35) π µν I ν ∝ b π µν a ν + b π µν ∇ ν ρ , (A36) q (cid:104) µ I ν (cid:105) ∝ b q (cid:104) µ a ν (cid:105) + b q (cid:104) µ ∇ ν (cid:105) ρ , (A37)where I µ := ∂ µ (cid:0) µT (cid:1) and b , b are scalar functions of ρ and p .Note that Eqs. (A31)–(A33) have been derived for a flat spacetime with signature (+ , − , − , − ) ,so that the comoving derivative is ˙ A = u µ ∂ µ A . Hence, in the comparison we will carry out in Sec.A5 the covariant derivatives in our general-relativistic formulation have to be replaced with a sim-ple partial derivative when comparing to DMNR12 .We recall that in the derivation of the
DMNR12 formulation, the full off-equilibrium distribu-tion function is decomposed as f ( x µ , p i ) = f ( x µ , p i ) + δf ( x µ , p i ) , (A38)where f is the local-equilibrium distribution function and δf denotes the deviation from it. Onethen defines the so-called generalized irreducible moments of order r of δf as ρ µ ··· µ l r := (cid:90) gd p (2 π ) p ( E ) r p (cid:104) µ · · · p µ l (cid:105) δf , (A39)where g is the number of internal degrees of freedom and E is defined by p µ =: Eu µ + p (cid:104) µ (cid:105) . Thefull off-equilibrium distribution function can now be expanded in momentum space in a basis ofirreducible tensors p (cid:104) µ · · · p µ l (cid:105) and orthogonal polynomials in E , where the generalized irreduciblemoments appear as coefficients. Some of these moments are directly connected to the dissipativecurrents, e.g., Π = − m ρ = − m (cid:90) gd p (2 π ) p δf , (A40) q µ = ρ µ = (cid:90) gd p (2 π ) p Ep µ δf, (A41) π µν = ρ µν = (cid:90) gd p (2 π ) p p (cid:104) µ p ν (cid:105) δf . (A42)Inserting Eq. (A38) into the special-relativistic version of the Boltzmann equation leads to anequation of motion for δf , δ ˙ f = − ˙ f − E p (cid:104) µ (cid:105) ∂ µ ( f + δf ) . (A43)This equation is then used to evaluate the comoving derivatives of the generalized irreducible This choice introduces a sign difference in some quantities, e.g., π µν NS , which are instead computed with the signature ( − , + , + , +) .MNRAS000
DMNR12 formulation, the full off-equilibrium distribu-tion function is decomposed as f ( x µ , p i ) = f ( x µ , p i ) + δf ( x µ , p i ) , (A38)where f is the local-equilibrium distribution function and δf denotes the deviation from it. Onethen defines the so-called generalized irreducible moments of order r of δf as ρ µ ··· µ l r := (cid:90) gd p (2 π ) p ( E ) r p (cid:104) µ · · · p µ l (cid:105) δf , (A39)where g is the number of internal degrees of freedom and E is defined by p µ =: Eu µ + p (cid:104) µ (cid:105) . Thefull off-equilibrium distribution function can now be expanded in momentum space in a basis ofirreducible tensors p (cid:104) µ · · · p µ l (cid:105) and orthogonal polynomials in E , where the generalized irreduciblemoments appear as coefficients. Some of these moments are directly connected to the dissipativecurrents, e.g., Π = − m ρ = − m (cid:90) gd p (2 π ) p δf , (A40) q µ = ρ µ = (cid:90) gd p (2 π ) p Ep µ δf, (A41) π µν = ρ µν = (cid:90) gd p (2 π ) p p (cid:104) µ p ν (cid:105) δf . (A42)Inserting Eq. (A38) into the special-relativistic version of the Boltzmann equation leads to anequation of motion for δf , δ ˙ f = − ˙ f − E p (cid:104) µ (cid:105) ∂ µ ( f + δf ) . (A43)This equation is then used to evaluate the comoving derivatives of the generalized irreducible This choice introduces a sign difference in some quantities, e.g., π µν NS , which are instead computed with the signature ( − , + , + , +) .MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke moments (A39) and leads to an infinite set of evolution equations for the latter. When this infinitesystem is truncated at the lowest order – corresponding to the 14-moment approximation – it leadsto the
IS79 -like equations, which have been given above by Eqs. (A31)–(A33).However, this procedure is ambiguous because, once the 14-moment approximation is applied,it is possible to obtain a closed evolution equation of the desired dissipative current from any choiceof r in Eq. (A39), for the irreducible moment of tensor rank corresponding to that of the respectivedissipative current . This can be seen from the r -dependence of the transport coefficients in TablesA3 and A4, which reflects the choice of the moment equation. The transport coefficients of IS79 are obtained if one sets r = 3 for the scalar moment, r = 2 for the vector moment and r =1 for the tensor moment. Microphysical properties of the system are encoded in the transportcoefficients obtained from a moment expansion of the collision integral of the Boltzmann equation.Instead of applying the 14-moment approximation, a power-counting scheme in Knudsen andinverse Reynolds numbers is used by Denicol et al. (2012b) to truncate the system. A4 Baier et al. 2008 (rBRSSS08)
The equations for the dissipative currents of the rBRSSS08 formulation read as follows [herewe use the version presented by Romatschke & Romatschke (2019), where they appear as Eq.(2.122)] τ Π ˙Π = Π NS − Π + δ Π ω µν ω µν + λ Π I (cid:104) µ (cid:105) I (cid:104) µ (cid:105) + ϕ Π π µν π µν + ϕ Π Π + g Π R + g Π u µ u ν R µν , (A44) τ π ˙ π (cid:104) µν (cid:105) = π µν NS − π µν + ∆ π π µν Θ + δ π π λ (cid:104) µ ω ν (cid:105) λ + δ π ω (cid:104) µλ ω ν (cid:105) λ + λ π I (cid:104) µ I ν (cid:105) + ϕ π π λ (cid:104) µ π ν (cid:105) λ + g π R (cid:104) µν (cid:105) + g π u λ u ρ R λ (cid:104) µν (cid:105) ρ . (A45)Note that we have defined I µ := ∇ µ ln e and that the heat current is absent because these ex-pressions refer to the Landau frame as the reference frame, where the heat currents are zero bydefinition.As for the previous formulations, we can use the first-order relations (A11)–(A13) to further Stated differently, it is possible to derive a relationship between the dissipative currents Π , q µ , and π µν , which are essentially the irreduciblemoments ρ , ρ µ , and ρ µν , and all the other ρ r , ρ µr , and ρ µνr , for any index r . As a result, there is ambiguity because it is possible to derive anequation of motion for Π , etc. from the equation of motion for ρ r , etc. for any r and not only for r = 0 . In Baier et al. (2008) the equations for the dissipative currents have been derived for the case of conformal fluids. Here we present their completeextension to the case of non-conformal fluids. MNRAS , 1–55 (2021)
RHD of non-perfect fluids I µ as I (cid:104) µ (cid:105) I (cid:104) µ (cid:105) = I µ I µ + u µ I µ u ν I ν ∝ ( b a µ + b ∇ µ ρ )( b a µ + b ∇ µ ρ ) + b Θ ∝ ( b ) a µ a µ + 2 b b a µ ∇ µ ρ + ( b ) ∇ µ ρ ∇ µ ρ + b Θ , (A46) I (cid:104) µ I ν (cid:105) ∝ ( b ) a (cid:104) µ a ν (cid:105) + 2 b b a (cid:104) µ ∇ ν (cid:105) ρ + ( b ) ∇ (cid:104) µ ρ ∇ ν (cid:105) ρ , (A47)where, again, the quantities b , b and b are scalar functions of ρ and p .The rBRSSS08 formulation is based on a systematic expansion in terms of gradients of thefluid variables in equilibrium and of the metric g µν , hence inheriting purely geometric terms in-volving the Ricci tensor R µν and the Ricci scalar R := R µµ . The dissipative currents are thenwritten as a power series in such gradients up to a given order. The individual terms for the shear-stress tensor respect its properties, namely that it is symmetric, trace-free and orthogonal to u .A resummation procedure is then applied to obtain relaxation-type equations, e.g., for theshear-stress tensor the first-order relation π µν = − ησ µν is used to derive u λ ∇ λ π (cid:104) µν (cid:105) ∝ − ηu λ ∇ λ σ (cid:104) µν (cid:105) − b σ µν Θ , (A48)which then leads to hyperbolic equations of motion if the term u λ ∇ λ σ (cid:104) µν (cid:105) is substituted using(A48) in the corresponding parabolic equations of motion. Note that b is a scalar function of ρ and p , and can be absorbed in the definitions of the transport coefficients. In addition, Eq. (A48)is accurate to second order in the gradients. In this way, the previously acausal behaviour is curedand causality is recovered. A5 Comparing different GRDHD formulations
After having introduced and reviewed four different second-order formulations of the equations ofGRDHD, namely,
IS79 , HL83 , DMNR12 and rBRSSS08 , we next proceed with a comparison ofthe sets of evolution equations for the dissipative currents. In order to facilitate the identificationof the transport coefficients in our notation with the transport coefficients in the original notations,we have distinguished between the different currents I iµ , as well as the kinematic acceleration a µ so far. However, as one can see from the first-order relations (A11)–(A13), a further reduction ispossible by using equations (A24)–(A30), (A34)–(A37) and (A46), (A47), respectively. To thisscope, we choose the pair { ρ, p } as our thermodynamical variables and write every gradient asa linear combination of the pair { a µ , ∂ (cid:104) µ (cid:105) ρ } . We apply this reduction in order to compare IS79 , DMNR12 and rBRSSS08 with
HL83 in Table A2. Thus, it is convenient to introduce the new
MNRAS000
MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
Formulation Order bulk-pressure equation heat-current equation shear-tensor equation − + − + − + IS79 O RK ΠΘ q µ Θ ω µν q ν π µν Θ π λ (cid:104) µ ω ν (cid:105) λ q µ I µ Π I (cid:104) µ (cid:105) q (cid:104) µ I ν (cid:105) π µν I ν DMNR12 O RK π µν σ µν ω µν q ν π λ (cid:104) µ ω ν (cid:105) λ σ µν q ν π λ (cid:104) µ σ ν (cid:105) λ Π σ µν rBRSSS08 O RK ΠΘ N/A N/A ∇ (cid:104) µ q ν (cid:105) π λ (cid:104) µ ω ν (cid:105) λ ∇ µ q µ q (cid:104) µ a ν (cid:105) q µ a µ q (cid:104) µ I ν (cid:105) q µ I µ O Θ ω λ (cid:104) µ ω ν (cid:105) λ ω µν ω µν a (cid:104) µ a ν (cid:105) a µ a µ a (cid:104) µ I ν (cid:105) a µ I µ I (cid:104) µ I ν (cid:105) I µ I µ O π µν π µν π λ (cid:104) µ π ν (cid:105) λ Π N/A
R R (cid:104) µν (cid:105) u µ u ν R µν u λ u ρ R λ (cid:104) µν (cid:105) ρ Table A2.
Comparison between different formalisms. Terms that are included in the
HL83 formulation but are missing in the given formulationare coloured in red and listed under ’ − ’. Similarly, those terms that are absent in the HL83 formulation but are included in the other approachesare coloured in blue and listed under ’ + ’. Terms that differ only by a scalar function of ρ and p are not counted as missing, as potential differencescan be absorbed in the definition of the transport coefficients. For each formulation we also use different rows to reflect the classification within thescheme with respect to the order of Knudsen or inverse Reynolds number (A4). Finally, terms without any symbol are not classified because theydo not occur in the corresponding approach. definition I µ := ∇ µ ρ . (A49)Notice that this reduction procedure leads to differences in some transport coefficients, which arecaptured by the previously introduced scalar functions { b i } i ∈ [1 , .Table A2 presents a quick overview of the various formulations considered, whose denomi-nation appears in the first column. The table is written such that terms that are included in the HL83 formulation but are missing in the set we want to compare with are coloured in red andlisted under ’ − ’. Similarly, those terms that are absent in the HL83 formulation but are includedin the other approaches are coloured in blue and listed under ’ + ’. Furthermore, terms that differonly by a scalar function of ρ and p are not counted as missing, as potential differences can beabsorbed in the definition of the transport coefficients. Furthermore, to facilitate the classificationof the various terms within the scheme (A4), i.e., with respect to the order of Knudsen or inverse MNRAS , 1–55 (2021)
RHD of non-perfect fluids O , O RK , O , and O , they are collected on different rows. Finally, termswithout any symbol are not classified because they do not occur in the corresponding approach.Note that the rBRSSS08 formulation includes second-order terms, but does not distinguishbetween ησ µν Π and ζπ µν Θ . However, these are not the only second-order terms present in the rBRSSS08 formulation. Another difference is of course the presence of terms associated with thecurvature of spacetime, which are absent in DMNR12 . Clearly, these terms cannot be classifiedby a certain order within the scheme (A4). However, if g µν is considered as an equilibrium fluidvariable, then these terms can be seen as of second order in Knudsen number. Finally, we recallthat the rBRSSS08 formulation is derived in the Landau frame, where the heat currents are zero.As a corollary to this comparison, we report in Tables A3 and A4 a list of all the transportcoefficients in our and in the original notation, while the definition of all I iµ is given in Table A5. APPENDIX B: DETAILS ON VISCOUS BLACK-HOLE ACCRETION
Given the very limited use and knowledge of the stationary solution of the spherically symmet-ric equations of GRDHD in a Schwarzschild spacetime, we here review the basic mathematicalexpressions and the strategy employed to obtain a numerical solution in the presence of a criticalpoint.
B1 Equations of GRDHD
We recall that the equations of GRDHD are given by Eqs. (1), (2) where J µ = J µ PF , and T µν isgiven by Eq. (4), together with Eq. (99) in the case in which the heat current and the shear-stresstensor are set to zero. By demanding the equations to be stationary and spherically symmetric,i.e., all state variables are functions of the circumference radius r only and the fluid four-velocityhas the form u µ = ( u t , u, , T , it is possible to obtain the following coupled, nonlinear system MNRAS000
We recall that the equations of GRDHD are given by Eqs. (1), (2) where J µ = J µ PF , and T µν isgiven by Eq. (4), together with Eq. (99) in the case in which the heat current and the shear-stresstensor are set to zero. By demanding the equations to be stationary and spherically symmetric,i.e., all state variables are functions of the circumference radius r only and the fluid four-velocityhas the form u µ = ( u t , u, , T , it is possible to obtain the following coupled, nonlinear system MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke
IS79 HL83 DMNR12 rBRSSS08 local source − / ζ V u µE | µ − ζ Θ − ζ r Θ − ζ ∇ ⊥ µ u µ Π NS - κTηβ ∆ µλ α | µ − κT ( ∇ (cid:104) µ (cid:105) ln T + a µ ) κ rq τ rW ψ Wr τ rV β h ∇ (cid:104) µ (cid:105) α N/A q µ NS - − ζ S ∆ α (cid:104) λ ( u E )∆ βµ (cid:105) ( u E ) u Eα | β − ησ µν η r σ µν − ησ µν π µν NS - / ζ V β ζβ τ r Π τ Π τ Π ˙Π κT β κT β τ rW N/A τ q ˙ q (cid:104) µ (cid:105) ζ S β ηβ τ rπ τ π τ π ˙ π (cid:104) µν (cid:105) − / τ Π − τ r Π δ r ΠΠ ∆ Π ΠΘ1 / ζ V a (cid:48) τ r Π τ r Π W δ Π q µ a µ τ r Π λ r Π π δ Π π µν σ µν ξ δ Π ω µν ω µν − / τ q − τ rW δ rWW N/A ∆ q q µ Θ κT a − τ rW τ rq Π /ψ Wr N/A δ q Π a µ κT a − τ rW τ rqπ /ψ Wr N/A δ q π µν a ν κT β = τ q τ rW N/A δ q ω µν q ν − τ rW λ rWW N/A δ q σ µν q ν − / τ π − τ rπ δ rππ − ( τ π d + τ ∗ π ) / ( d −
1) ∆ π π µν Θ2 ζ s a (cid:48) τ rπW τ rπ δ π q (cid:104) µ a ν (cid:105) ζ S β = 2 τ π τ rπ − λ /η δ π π λ (cid:104) µ ω ν (cid:105) λ − τ rπ λ rππ δ π π λ (cid:104) µ σ ν (cid:105) λ τ rπ λ rπ Π δ π Π σ µν τ rπ λ rπ Π λ δ π ω λ (cid:104) µ ω ν (cid:105) λ / ζ V α ζα − τ r Π l r Π W L Π1 ∇ µ q µ κT α κT α l rq Π τ rW /ψ Wr N/A L q ∇ (cid:104) µ (cid:105) Π κT α κT α − l rqπ τ rW /ψ Wr N/A L q ∇ ν π (cid:104) µ (cid:105) ν ζ S a ηα l rπW τ rπ L π ∇ (cid:104) µ q ν (cid:105) Table A3.
Table of transport coefficients I: we here match the transport coefficients used in this work with the corresponding transport coefficientsin the original papers of the
IS79 , HL83 , DMNR12 , and rBRSSS08 formulations. of ordinary differential equations (ODEs) in Schwarzschild coordinates: dρdr = − ρr M/ ( E r ) − Π r/ [( ρh + Π) τ Π u ] − u / E c s,t − u / E , (B1) dudr = ur M/ ( E r ) − Π r/ [( ρh + Π) τ Π u ] − c s,t c s,t − u / E , (B2) d Π dr = − Π (cid:0) c s,t − u / E (cid:1) / ( uτ Π ) + ζ [ M/ ( E r ) − Π r/ [( ρh + Π) τ Π u ] − u / E ] / ( τ Π r ) c s,t − u / E , (B3) dhdr = − ρh + Π ρr (cid:2) c s,t − ( ζ/τ Π − Π) / ( ρh + Π) (cid:3) [ M/ ( E r ) − Π r/ [( ρh + Π) τ Π u ] − u / ( E )] c s,t − u / E , (B4) MNRAS , 1–55 (2021)
RHD of non-perfect fluids IS79 HL83 DMNR12 rBRSSS08 local source − / ζT Λ Π Π u µ I iµ ζγ T τ r Π λ r Π W Λ Π q µ I iµ ξ λ Π I iµ I iµ − / κT Λ q q µ u ν I iν κT (1 − γ ) λ rq Π τ rW /ψ Wr N/A Λ q Π I i (cid:104) µ (cid:105) κT (1 − γ ) λ rqπ τ rW /ψ Wr N/A Λ q π µν I iν − / ηT Λ π π µν u λ I iλ ηγ T − τ rπ λ rπW Λ π q (cid:104) µ I i ν (cid:105) λ λ π I i (cid:104) µ I i ν (cid:105) ξ /η ϕ Π π µν π µν ξ /ζ ϕ Π Π λ /η ϕ π π λ (cid:104) µ π ν (cid:105) λ ξ g Π R ξ g Π u µ u ν R µν κ g π R (cid:104) µν (cid:105) κ ∗ − κ ) g π u λ u ρ R λ (cid:104) µν (cid:105) ρ Table A4.
Table of transport coefficients II: the same as Table A3 but for those terms missing in the previous Table. Note that the definitions of thecurrents I iµ differ for different formalisms and a summary is given in Table A5. IS79 HL83 DMNR12 rBRSSS08
N/A I iµ = ∇ µ ( c i /T ) , i ∈ { , . . . , } I µ = ∂ µ (cid:0) µT (cid:1) I µ = ∇ µ ln e Table A5.
Summary of the various definitions of the currents I iµ given for the four formulations considered here. Note that the IS79 formulationneglects gradients of transport coefficients and that the coefficients c i appearing in the currents of the HL83 formulation are given by c { ,..., } = { τ Π /ζ, α , τ q /κT, α , τ π /η } . where E := u t = − (cid:112) − M /r + u . Furthermore, the system is characterised by two conservedquantities, namely, the mass-accretion rate ˙ M and the “viscous” Bernoulli constant B : ˙ M := 4 πρur , (B5) B := ( ρh + Π) E /ρ . (B6)Note that in the inviscid limit lim Π → B = B PF := h E = hu t , where B PF denotes the relativistic,inviscid Bernoulli constant. Using ˙ M and B , we can express ρ and Π in terms of u , h and r : ρ = ˙ M / (4 πur ) , (B7) Π = ρ B / E − ˙ M h/ (4 πur ) , (B8)so that the viscous speed of sound becomes c s,t = ( γ e − B − EB + ζτ Π π E ur B ˙ M , (B9)
MNRAS000
MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke and Eqs. (B2) and (B4) simplify to dudr = ur M/ ( E r ) − ( B − E h ) r/ ( B τ Π u ) − c s,t c s,t − u / E , (B10) dhdr = − r [( γ e − B − E ) / E + ( B − h E ) / E ] [ M/ ( E r ) − ( B − E h ) r/ ( B τ Π u ) − u / E ] c s,t − u / E . (B11)Equations (B10)–(B11) do not have an analytic solution and their non-trivial numerical solution isdiscussed in more detail in the next section. B2 Singular point analysis and integration
The first step in the solution of Eqs. (B10), (B11) is to choose the location of the sonic point r s ,which is defined as the radial coordinate where u / E = c s,t . This corresponds to the positionwhere the infalling fluid velocity equals that of the speed of sound as observed at rest from spatialinfinity, and where Eqs. (B10), (B11) become singular. Because of that, we choose to start thenumerical integration procedure directly at r s and integrate inward (i.e., supersonic portion) andoutward (i.e., subsonic portion) from there. In this way, we can make use of exact results for thestates and their first derivatives at the sonic point in order to initialise numerical integration. Inthe following we will elaborate on the exact method. Using the inviscid solution with a polytropicconstant (see discussion in Sec. 7), it is possible to obtain values for ˙ M and B PF . In the inviscidcase, all state variables at the sonic point are automatically determined by choosing r s (see Hawleyet al. 1984, for details).Setting B = B PF and using the value for ˙ M found for the inviscid case, we can compute ourinitial values at r s for the viscous case. As in the inviscid case, we demand that the derivatives atthe sonic point constitute a removable singularity. Hence, from equations (B10) and (B11) we findthe conditions M E r − B − E h B rτ Π u − c s,t , (B12) c s,t − u E . (B13)Numerical root finding via Mathematica (Wolfram Research 2020) yields the values of the radialfour-velocity component and of the specific enthalpy at the sonic radius, u s and h s . Because multi-ple solutions exist, we select the solution satisfying the conditions: u s < and h s > , which wasfound to be unique for the choice of constant ζ and τ Π . Next, accurate values for the derivatives du/dr and dh/dr are needed at the sonic point because we would otherwise rely on the evaluationof the right-hand side of Eqs. (B10) and (B11) at r s for the first step of the integration procedure. MNRAS , 1–55 (2021)
RHD of non-perfect fluids u , r , and h can beexpressed in terms of a parameter ξ , i.e., u = u ( ξ ) , r = r ( ξ ) , h = h ( ξ ) , such that when writingEqs. (B10)–(B11) symbolically as du/dr = A/B and dh/dr = C/B , they can be effectivelyrewritten as du/dξ = A , dr/dξ = B , and dh/dξ = C . More specifically, we express Eqs. (B10)–(B11) symbolically as drdξ = r (cid:18) c s,t − u E (cid:19) , (B14) dudξ = u (cid:18) M E r − B − E h B rτ Π u − c s,t (cid:19) , (B15) dhdξ = − (cid:20) ( γ e − B − EE + B − h EE (cid:21) (cid:20) M E r − B − E h B rτ Π u − u E (cid:21) . (B16)Note that ξ does not have a physical interpretation and should be seen simply as a mathematicalparameter. However, with this parameterisation, each of the solutions r ( ξ ) , u ( ξ ) , and h ( ξ ) can bethought as consisting of two branches on either side of the sonic point: i.e., one branch for r < r s and one for r > r s . In this way, and as will be found later, the solution at the sonic point can beobtained in the limit ξ → ±∞ .Hence, Eqs. (B14)–(B16) effectively represent an independent system whose paths in the ( u, h, r ) -space, i.e., the phase space of Eqs. (B14)–(B16), correspond to solutions of the origi-nal system, i.e., Eqs. (B10) and (B11). Thus, the original system can be viewed as differentialequations for the phase paths of Eqs. (B14)– (B16). The solution ( u s , h s ) at the sonic point r s constitutes a solution ( u s , h s , r s ) of Eqs. (B14)–(B16) at the so-called equilibrium point, where du/dξ = dh/dξ = dr/dξ = 0 . Notice the dependence on the specific choice of ζ = ζ ( u, h, r ) and τ Π = τ Π ( u, h, r ) .Equations (B14)–(B16) can can be linearised at the sonic point, i.e., u ≈ u s + δu , h ≈ h s + δh and r ≈ r s + δr , to obtain a local solution (Jordan & Smith 2007). This linearisation procedureand the ansatz δu = δ u exp( λξ ) , δh = δ h exp( λξ ) and δr = δ r exp( λξ ) leads to a an eigenvalueproblem, which can be solved numerically (Wolfram Research 2020). The outcome is a set of threeeigenvalues { λ , λ , λ } and eigenvectors with the properties λ < , λ > and λ ≈ . The This condition follows from the requirement that the perturbation behaves as ∼ e λξ , where λ is a constant eigenvalue that can be either positiveor negative, and that it vanishes at the sonic point; see also below.MNRAS000
RHD of non-perfect fluids u , r , and h can beexpressed in terms of a parameter ξ , i.e., u = u ( ξ ) , r = r ( ξ ) , h = h ( ξ ) , such that when writingEqs. (B10)–(B11) symbolically as du/dr = A/B and dh/dr = C/B , they can be effectivelyrewritten as du/dξ = A , dr/dξ = B , and dh/dξ = C . More specifically, we express Eqs. (B10)–(B11) symbolically as drdξ = r (cid:18) c s,t − u E (cid:19) , (B14) dudξ = u (cid:18) M E r − B − E h B rτ Π u − c s,t (cid:19) , (B15) dhdξ = − (cid:20) ( γ e − B − EE + B − h EE (cid:21) (cid:20) M E r − B − E h B rτ Π u − u E (cid:21) . (B16)Note that ξ does not have a physical interpretation and should be seen simply as a mathematicalparameter. However, with this parameterisation, each of the solutions r ( ξ ) , u ( ξ ) , and h ( ξ ) can bethought as consisting of two branches on either side of the sonic point: i.e., one branch for r < r s and one for r > r s . In this way, and as will be found later, the solution at the sonic point can beobtained in the limit ξ → ±∞ .Hence, Eqs. (B14)–(B16) effectively represent an independent system whose paths in the ( u, h, r ) -space, i.e., the phase space of Eqs. (B14)–(B16), correspond to solutions of the origi-nal system, i.e., Eqs. (B10) and (B11). Thus, the original system can be viewed as differentialequations for the phase paths of Eqs. (B14)– (B16). The solution ( u s , h s ) at the sonic point r s constitutes a solution ( u s , h s , r s ) of Eqs. (B14)–(B16) at the so-called equilibrium point, where du/dξ = dh/dξ = dr/dξ = 0 . Notice the dependence on the specific choice of ζ = ζ ( u, h, r ) and τ Π = τ Π ( u, h, r ) .Equations (B14)–(B16) can can be linearised at the sonic point, i.e., u ≈ u s + δu , h ≈ h s + δh and r ≈ r s + δr , to obtain a local solution (Jordan & Smith 2007). This linearisation procedureand the ansatz δu = δ u exp( λξ ) , δh = δ h exp( λξ ) and δr = δ r exp( λξ ) leads to a an eigenvalueproblem, which can be solved numerically (Wolfram Research 2020). The outcome is a set of threeeigenvalues { λ , λ , λ } and eigenvectors with the properties λ < , λ > and λ ≈ . The This condition follows from the requirement that the perturbation behaves as ∼ e λξ , where λ is a constant eigenvalue that can be either positiveor negative, and that it vanishes at the sonic point; see also below.MNRAS000 , 1–55 (2021) M. Chabanov, L. Rezzolla and D.H. Rischke corresponding set of local solutions is then δu ( ξ ) = ( δ u ) exp( λ ξ ) ; δh ( ξ ) = ( δ h ) exp( λ ξ ) ; δr ( ξ ) = ( δ r ) exp( λ ξ ) , (B17) δu ( ξ ) = ( δ u ) exp( λ ξ ) ; δh ( ξ ) = ( δ h ) exp( λ ξ ) ; δr ( ξ ) = ( δ r ) exp( λ ξ ) , (B18)where [( δ u ) i , ( δ h ) i , ( δ r ) i ] denotes the eigenvector which belongs to the i -th eigenvalue λ i . Now, werequire the linearised solution to pass through ( u s , h s , r s ) exactly, which means that we look forlinearised solutions that fulfil δu i ( ξ ) = δh i ( ξ ) = δr i ( ξ ) = 0 for an arbitrary ξ . Because λ ≈ , itis not possible to fulfil the requirement δu = δh = δr = 0 for any value of ξ ; hence, we neglectthis eigenvalue. Similarly, because λ and λ are of opposite signs, superpositions of both sets arenot allowed for the same reason. Thus, we recover the desired solution in the limit ξ → ∞ for theset (B17) and ξ → −∞ for the set (B18). For each eigenvalue i = 1 , , the derivatives are givenby ( du/dr ) s = ( δ u ) i / ( δ r ) i and dh/dr = ( δ h ) i / ( δ r ) i , respectively, and we select the eigenvaluewhose eigenvector yields ( du/dr ) s > .Finally, we start our numerical integration by calculating the first steps on either side of thesonic point, i.e., at r l = r s − ∆ r and r r = r s + ∆ r , employing forward and backward finite-differences of first order with radial stepsize ∆ r . The integration then proceeds with a fourth-order L -stable singly diagonally implicit Runge-Kutta (SDIRK) method (see Table 6.5. of Hairer& Wanner 1996 and the corresponding Butcher tableau). It is worth noticing that this methodincludes an embedded third-order formula and a continuous solution, both of which can be foundin Hairer & Wanner (1996). We use a simple fixed-point iteration procedure to solve the implicitequation at each stage, i.e., given the solution at a specific radius r n , where r n +1 = r + ∆ r marks the next grid point, we use the solution at r n as an initial guess and iterate on the k i ’s,which are defined by the following general formulas derived assuming the ODE has the form dy/dr = f ( r, y ) and an s -stage Runge-Kutta method y n +1 = y n + ∆ r s (cid:88) i =1 b i k i , (B19) k i = f (cid:32) r n + c i ∆ r, y n + h s (cid:88) j a ij k j (cid:33) , i = 1 , . . . , s , (B20)and where the coefficients { c i , a ij , b i } are given in the form of Butcher tableaus in Table B1. Thenumerical integration is obtained using 700,000 points on a grid that has a higher resolution inthe most delicate portions of the solution, i.e., at the event horizon and at the sonic point. Aftera cubic-spline interpolation, the numerical solution of the ODEs is then used as initial conditionsfor the BHAC evolution.
MNRAS , 1–55 (2021)
RHD of non-perfect fluids / / / / / /
20 17 / − /
25 1 / / / − / /
544 1 /
41 25 / − /
48 125 / − /
12 1 / / − /
48 125 / − /
12 1 / Table B1.
Butcher tableau for a L -stable SDIRK method of order 4 [cf., Eqs. (B19)].This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000