General-Relativistic Large-Eddy Simulations of Binary Neutron Star Mergers
DDraft version March 8, 2017
Preprint typeset using L A TEX style emulateapj v. 5/2/11
GENERAL-RELATIVISTIC LARGE-EDDY SIMULATIONS OF BINARY NEUTRON STAR MERGERS
David Radice
Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA andDepartment of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA
Draft version March 8, 2017
ABSTRACTThe flow inside remnants of binary neutron star (NS) mergers is expected to be turbulent, becauseof magnetohydrodynamics instability activated at scales too small to be resolved in simulations. Tostudy the large-scale impact of these instabilities, we develop a new formalism, based on the large-eddysimulation technique, for the modeling of subgrid-scale turbulent transport in general relativity. Weapply it, for the first time, to the simulation of the late-inspiral and merger of two NSs. We find thatturbulence can significantly affect the structure and survival time of the merger remnant, as well asits gravitational-wave (GW) and neutrino emissions. The former will be relevant for GW observationof merging neutron stars. The latter will affect the composition of the outflow driven by the mergerand might influence its nucleosynthetic yields. The accretion rate after black-hole formation is alsoaffected. Nevertheless, we find that, for the most likely values of the turbulence mixing efficiency,these effects are relatively small and the GW signal will be affected only weakly by the turbulence.Thus, our simulations provide a first validation of all existing post-merger GW models.
Subject headings:
Gravitational waves — Stars: neutron — Turbulence INTRODUCTIONThe typical outcome of the merger of two neutron stars(NSs) is expected to be the formation of an hypermassiveneutron star (HMNS): a massive NS temporarily sup-ported against gravitational collapse by its fast differen-tial rotation, although prompt black hole (BH) forma-tion might occur for large masses and/or soft equationsof state (EOS) (Baiotti & Rezzolla 2016, and referencestherein). Its survival time and, in general, its properties,are important for the multimessenger signature of NSmergers and for their nucleosynthetic yields. The HMNShas a magnetar-level B-field, and it is a bright sourceof neutrinos (Sekiguchi et al. 2011; Kiuchi et al. 2014).These could drive baryon-rich winds (Dessart et al. 2009;Siegel et al. 2014). The presence of an HMNS could sig-nificantly boost the neutrino annihilation rates at high-latitudes (Richers et al. 2015; Perego et al. 2017) andperhaps contribute to the launching of a relativistic jetand a short γ -ray burst (SGRB) (Nakar 2007). Neutri-nos could also affect the yield and electromagnetic signa-ture of the r-process nucleosynthesis in the binary ejecta(Wanajo et al. 2014; Metzger & Fern´andez 2014). Long-lived massive NSs created in mergers might power the X-ray tails observed in some SGRB (Rowlinson et al. 2013;Lasky et al. 2014; Gao et al. 2016). Finally, gravitationalwaves (GWs) from the HMNS could be used to constrainits EOS (Bauswein & Janka 2012; Takami et al. 2014;Bernuzzi et al. 2015a; Radice et al. 2016a).Despite the rapid recent progress of general-relativisticmagnetohydrodynamics (GRMHD) simulations (Rez-zolla et al. 2011; Kiuchi et al. 2014; Ruiz et al. 2016),the impact of magnetoturbulence on the structure andsurvival time of the HMNS is highly uncertain. Themagnetorotational instability (MRI) (Balbus & Hawley1991) is expected to operate inside the HMNS, drivethe redistribution of angular momentum and affect itslifetime and properties (Duez et al. 2006; Siegel et al.2013). Unfortunately, the fastest growing mode of the MRI in these systems is inaccessible even to the highest-resolution simulations (Kiuchi et al. 2015).A possible way to model the impact of turbulent trans-port of angular momentum in the HMNS would be touse an effective viscosity (Duez et al. 2004). This ap-proach is made difficult by the fact that the Navier-Stokes equations describing relativistic viscous flows areknown to exhibit a number of unphysical pathologies(Hiscock & Lindblom 1985; Kostadt & Liu 2000). Thereare more complex fluid models that do not have theseshortcomings (Andersson & Comer 2007; Rezzolla &Zanotti 2013). However, they are also not entirely with-out problems (Majorana & Motta 1985; Hiscock & Lind-blom 1988), are difficult to implement ( e.g. Takamoto& Inutsuka 2011), and their non-linear properties arepoorly understood. More importantly, they contain alarge number of transport coefficients that have no clas-sical counterpart. These have no clear physical meaningand are not even in principle measurable (Geroch 1995;Lindblom 1996).Here, we propose an alternative approach. Our start-ing point is the observation that turbulence models donot have to be restricted to the class of equations de-scribing fluids with physical viscosity or heat transfer.Instead, we develop an effective model based on a GRextension of the Newtonian large-eddy simulation (LES)framework ( e.g.
Miesch et al. 2015). Our model, whilerecovering the Navier-Stokes equations in the Newtonianlimit, does not correspond to or have the same limita-tions as any relativistic theory of viscous flows.In this
Letter , after a brief description of the GRLESformulation, we present, for the first time, simulationsin full-GR of merging NS with a realistic, tabulated, nu-clear EOS, neutrino cooling, and parametrized turbulenttransport. We show that turbulence could influence thestructure of the HMNS, as well as its GW and neutrinoemissions. On the other hand, for the most realistic val-ues of the turbulent viscosity, these effects appear to be a r X i v : . [ a s t r o - ph . H E ] M a r small and our simulations provide an important confir-mation of a number of previous results where turbulenttransport was not included. FORMULATIONOur starting point is the stress energy tensor of a per-fect fluid T µν = ρhu µ u ν + pg µν , (1)where ρ , h , u µ and g µν are, respectively, density, specificenthalpy, four-velocity, and the spacetime metric.In numerical relativity, spacetime is decomposed inspace-like slices with normal n µ . We decompose T µν withrespect to n µ as T µν = En µ n ν + S µ n ν + S ν n µ + S µν , (2)where E = T µν n µ n ν = ρhW − p , (3) S µ = − γ µα n β T αβ = ρhW v µ , (4) S µν = γ µα γ µβ T αβ = S µ v ν + pγ µν , (5)and γ µν , v µ , p , and W are, respectively, the spatial met-ric, the three-velocity, the pressure, and the Lorentz fac-tor.The equations of energy and momentum conservationare ∂ t (cid:0) √ γS i (cid:1) + ∂ j (cid:104) α √ γ (cid:0) S ij + S i n j (cid:1)(cid:105) = α √ γ (cid:16) S jk ∂ i γ jk + 1 α S k ∂ i β k − E∂ i log α (cid:17) , (6) ∂ t (cid:0) √ γE (cid:1) + ∂ j (cid:104) α √ γ (cid:0) S j + En j (cid:1)(cid:105) = α √ γ (cid:16) K ij S ij − S i ∂ i log α (cid:17) , (7)where α , β i , K ij are, respectively, the lapse function,shift vector, three-metric, and extrinsic curvature. √ γ is the spatial volume element. These equations are thenclosed with an EOS and equations describing the conser-vation of the baryon and lepton numbers.Equations (6) and (7) contain modes at all scales, but,in numerical simulations, only modes resolved with suf-ficiently many grid zones can develop. In essence, anysimulation deals only with a coarse-grained version ofthe hydrodynamics equations. In the LES framework,this observation is made rigorous with the introductionof a linear filtering operation u (cid:55)→ u that removes featuresat scales smaller than a given ∆. Here, we adopt for thefiltering operator the cell-averaging of the finite-volumediscretization of the equations. We leave the investiga-tion of more advanced filters for future work. If we filterEqs. (6) and (7) we obtain ∂ t (cid:0) √ γS i (cid:1) + ∂ j (cid:104) α √ γ (cid:0) S ij + S i n j (cid:1)(cid:105) = α √ γ (cid:16) S jk ∂ i γ jk + 1 α S k ∂ i β k − E∂ i log α (cid:17) , (8) ∂ t (cid:0) √ γE (cid:1) + ∂ j (cid:104) α √ γ (cid:0) S j + En j (cid:1)(cid:105) = α √ γ (cid:16) K ij S ij − S i ∂ i log α (cid:17) . (9) Note that Eqs. (8) and (9) are exact, but are not closed.The reason is that S i v j cannot be expressed only in termsof other coarse-grained quantities. A closure is needed: S i v j = S i v j + τ ij . (10) τ ij is the so-called subgrid-scale turbulent tensor. Similarterms appear in the coarse graining of the baryon and lep-ton number conservation equations, but, for simplicity,we will neglect them here. Simulations usually assume τ ij = 0. Here, instead, we will use τ ij to model small-scale turbulence in merger simulations. To do so, in anal-ogy with the classical Newtonian closure of Smagorinsky(1963), we choose the ansatz τ ij = − ν T ρhW (cid:20) (cid:0) ∇ i v j + ∇ j v i (cid:1) − ∇ k v k γ ij (cid:21) , (11)where ∇ is the covariant derivative compatible with γ ij .The quantity ν T has a dimension of a viscosity. On di-mensional grounds, we are led to assume ν T = (cid:96) mix c s , (12)where c s is the local sound speed, and (cid:96) mix , often calledthe mixing length, is a characteristic length over whichturbulence operates. Note that ν T is not a physical vis-cosity; indeed, its definition depends on the numericalgrid and on the Eulerian observer n µ . This is expected,because the notion of resolved and unresolved scales isobserver dependent in relativity. ν T should be calibratedon the basis of highly-resolved simulations and/or usingself-similarity methods ( e.g. Germano et al. 1991). Weleave this task for future work. For now, we will treat (cid:96) mix as a free parameter. Assuming MRI turbulence, itis then natural to set (cid:96) mix ∼ λ MRI , where (Duez et al.2006) λ MRI ∼ (cid:18) Ω4 rad ms − (cid:19) − (cid:18) B G (cid:19) . (13)Equations (8), (9), (11), and (12), together with theEOS, and the continuity equations are what we refer toas the GRLES equations. We verified, by repeating theanalysis of Hiscock & Lindblom (1985) and numerically,that the GRLES equations are not affected by the samepathologies as the relativistic Navier-Stokes equations. IMPLEMENTATIONWe implement the GRLES equations into the general-relativistic hydrodynamics (GRHD) code
WhiskyTHC (Radice & Rezzolla 2012; Radice et al. 2014a,b). Withour current choice of the filtering operator, this onlyamounts to the inclusion of τ ij in the equations. Wetreat the viscous fluxes in a flux-conservative way andwe self-consistently include the turbulent stress-tensor inthe energy and momentum source terms, as well as inthe calculation of the spacetime geometry.For the simulations presented here, we use the micro-physical EOS of Lattimer & Swesty (1991) with nuclearcompressibility parameter K = 220 MeV. Neutrino cool-ing is treated with the scheme presented in (Radice et al.2016b). For (cid:96) mix , we consider the values 0 (our referencerun), 5, 25, and 50 meters. Over this range, 5 metersis the most likely value for (cid:96) mix given Eq. (13), while 50 t [ ms ] ρ [ g · c m − ] ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m
Figure 1.
Maximum density in the collapse of a differentiallyrotating equilibrium configuration. Turbulent transport of angularmomentum leads to an accelerated collapse. meters might be unphysically large, in the light of thelack of convergence observed in the 17-meter resolutionsimulation of Kiuchi et al. (2015).As a first example, we consider the evolution of anequilibrium configuration constructed with the
RNS code(Stergioulas & Friedman 1995). The initial configura-tion has gravitational mass M (cid:39) . M (cid:12) and angularmomentum J/M (cid:39) . G/c . We use the differentialrotation law of Komatsu et al. (1989), which, in the New-tonian limits reduces toΩ = Ω c (cid:16) (cid:36)R e (cid:17) , (14)where (cid:36) is the cylindrical radius, Ω c is the angular veloc-ity at the center, and R e is the stellar equatorial radius.The resolution for this test is (cid:39)
370 m.We plot the maximum density as a function of timein Fig. 1. As expected on the basis of previous work(Duez et al. 2004), the inclusion of turbulent viscosityresults in the transport of angular momentum leading togravitational collapse. This test shows that
WhiskyTHC is able to capture the effect of turbulent viscosity evenat low resolution. BINARY NEUTRON STAR MERGERSWe consider the last ∼ . M (cid:12) NSs. We already evolved this binary in Bernuzziet al. (2015b), where a description of the properties ofthe initial data is also given. For the evolution, we usethe high-resolution setup of Bernuzzi et al. (2015b), withthe improvements discussed in Radice et al. (2016a). Weperform simulations with resolutions, on the finest re-finement level of ∼
185 m and ∼
246 m. We present re-sults from the high-resolution simulations. In the low-resolution simulations, there are quantitative, but notqualitative differences.We find that turbulent viscosity has a much less ob-vious impact on the evolution of the HMNS than whatcould have been anticipated on the basis of the ideal-ized model in Sec. 3. In the first few milliseconds aftermerger, turbulent transport results in a decrease of thecompactness, as can be seen from the maximum densityevolution (Fig. 2; left panel). Over longer timescales,the behavior is non-linear. The (cid:96) mix -25-m HMNS is themost compact remnant and collapses to a BH ∼
17 ms after merger. The (cid:96) mix -5-m remnant is only slightly lesscompact than that of the reference simulation (cid:96) mix -0-m.BH formation occurs at ∼
20 ms and ∼
22 ms after mergerfor the (cid:96) mix -5-m and (cid:96) mix -0-m binaries, respectively. The (cid:96) mix -50-m HMNS is the least compact and does not col-lapse to a BH within our simulation time. For the modelsthat collapse within our simulation time, we observe theformation of a massive ( ∼ . M (cid:12) ) accretion disk (Fig. 2;right panel). As could have been anticipated, the accre-tion rate is larger for simulations with larger (cid:96) mix .The reason for the different evolutions of the remnantcan be understood from the analysis of its internal struc-ture (Fig. 3). The rotational profile established in theHMNS after the initial, very dynamical, phase is qualita-tively different from that of Eq. (14), as also pointed outby Shibata et al. (2005); Kastaun et al. (2016); Hanauskeet al. (2016); Ciolfi et al. (2017). Consistently with theseprevious studies, we find in the (cid:96) mix -0-m simulation anHMNS composed of a slowly rotating core and a rotation-ally supported massive envelope. As the mixing lengthincreases, the structure of the HMNS is altered due tointerplay between three competing effects. First, angu-lar momentum redistribution spins up the core, reducingits compactness. Second, the loss of angular momentumfrom the massive envelope results in a compression theHMNS. Third, as more kinetic energy is converted intothermal energy by turbulent dissipation, the inner corebecomes hotter and expands because of the increasedpressure. The interplay between these effects is com-plicated by the fact that the angular momentum of theHMNS is not conserved, but is radiated in GWs at a rateproportional to that of the gravitational binding energy(Bernuzzi et al. 2015b). For this reason, as the HMNScontracts, it becomes more bound and at the same timeit looses angular momentum support.The first and third effect are dominant at early times,so the effect of turbulent viscosity is to monotonicallyreduce the compactness in the first few milliseconds af-ter merger. Later, all three effects become important.At this stage, energy and angular momentum lossesto GW play an important role. In the case of the (cid:96) mix -5-m binary, the envelope remains centrifugally sup-ported (Fig. 3; upper-left panel), so the compactness isslightly decreased compared to the reference run withoutturbulence dissipation. For the (cid:96) mix -25-m binary, the ef-fect of turbulent transport is qualitatively similar to the (cid:96) mix -5-m binary at early times. Later, its envelope con-tracts causing the growth of the central density (Fig. 2;left panel) and early BH formation. Finally, in the caseof the (cid:96) mix -50-m run the thermal effect prevails; the hotspots formed in the contact layer at the time of mergersink to the center and enhance the core temperature to ∼
70 MeV. The increased thermal support in the layerswith ∼ · · g · cm − inflates the HMNS. The reducedcompactness, in turn, results in a decrease of the angu-lar momentum loss due to GW and prevents its collapsewithin the simulation time.The changes in the HMNS structure are reflected in itsmultimessenger emissions. The total energy radiated inGWs (Fig. 4; left panel) is closely related to the rate ofincrease of the HMNS compactness. For this reason, atearly times, the amplitude of the signal slightly decreaseswith (cid:96) mix , while, over longer timescales, the behavior isnon-monotonic. The characteristic GW frequency after − t − t mrg [ ms ] ρ [ g · c m − ] ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m t − t BH [ ms ] M d i s k [ − · M (cid:12) ] ‘ mix = ‘ mix = ‘ mix =
25 m
Figure 2.
Maximum density ( left panel ) and (baryonic) disk mass ( right panel ). The disk mass is computed as the total mass outside theapparent horizon. The impact of turbulent mixing on the compactness of the HMNS is non-trivial and non-monotonic. Turbulent angularmomentum transport results in larger accretion rates after BH formation. r [ km ] h Ω i xy [ r a d · m s − ] t − t mrg ’
10 ms ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m r [ km ] h ρ i xy [ g · c m − ] t − t mrg ’
10 ms ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m z [ k m ] ‘ mix = ‘ mix = − − −
10 0 10 20 30 x [ km ] z [ k m ] ‘ mix =
25 m − − −
10 0 10 20 30 x [ km ] ‘ mix =
50 m06121824303642485460 T [ M e V ] ρ [ g · c m − ] t − t mrg ’
10 ms
Figure 3.
Upper panels: angle-averaged angular velocity ( left ) and density ( right ) on the equatorial plane.
Lower panels: temperatureand density in the meridional plane. All data are shown at ∼
10 ms after merger. The white contours in the lower panel are the isodensitycontours for ρ = 10 , , , , , and 5 · g · cm − . Turbulent dissipation leads to angular momentum transport and enhancedthermalization. merger (Fig. 4; right panel) is, instead, only weakly af-fected, with the exception of a slight growth before BHformation, which is a commonly observed feature ( e.g. Radice et al. 2016a).The neutrino emission (Fig. 5) is also influenced by theturbulent dissipation and the consequently higher tem-peratures in the HMNS. The luminosity of neutrinos ofall flavours increases with the mixing length parameterup to (cid:96) mix = 25 m. The luminosity of the (cid:96) mix -50-m sim-ulation is, however, smaller than that of the (cid:96) mix -25-msimulation. This is possibly because, in the (cid:96) mix -50-m HMNS, the maximum of the temperature occurs at thecenter, while for the other models it is off-centered (seealso Kastaun et al. 2016). DISCUSSION AND CONCLUSIONSWe have developed a new framework for the modelingof turbulence in full-GR simulations. Our approach isbased on a relativistic extension of the large-eddy simu-lation technique, which represents the state-of-the-art forturbulence modeling in classical hydrodynamics (Mieschet al. 2015). Our method can naturally exploit turbu- − t − t mrg [ ms ] E G W [ M (cid:12) · c ] ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m − t − t mrg [ ms ] f G W [ k H z ] ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m
Figure 4.
Total energy radiated in GW ( left panel ) and instantaneous GW frequency ( right panel ). The former is smoothed using arunning average with a 0.1-ms window. Turbulent transport can influence the GW luminosity starting from the early post-merger. TheGW instantaneous frequency is, instead, only weakly affected. t − t mrg [ ms ] L X [ · e r g · s − ] X = ν e ‘ mix = ‘ mix = ‘ mix =
25 m ‘ mix =
50 m t − t mrg [ ms ] X = ¯ ν e t − t mrg [ ms ] X = ν x Figure 5.
Electron ( right panel ), anti-electron ( middle panel ) and heavy-lepton ( right panel ) neutrino luminosities. The increasedtemperature of the HMNS due to turbulent dissipation leads to an increase in the neutrino luminosity for all species. This effect seems tobe partially suppressed for the (cid:96) mix -50-m simulation. The sudden drops in the emission for some of the simulations ∼
20 after merger aredue to BH formation. lent closures developed in Newtonian physics, is simpleto implement, robust, and stable.As a first application, we have employed a turbulentviscosity closure to study the effect of angular momen-tum transport and dissipation in NS mergers. We haveperformed, for the first time, general-relativistic large-eddy simulations of merging NS with microphysical nu-clear EOS and neutrino cooling. We have found thatturbulence can modify the structure and collapse timeof the merger remnant. These, in turn, are reflected inthe GW and neutrino emissions from the HMNS. Theaccretion rate after BH formation is also affected.The total energy radiated in GW is the most affectedquantity, since it closely tracks the contraction of theHMNS on its way to the final collapse to BH. In the pres-ence of very efficient turbulent transport, the effectiveviscosity might mask changes in the compactness of theHMNS that would otherwise be attributable to changesin the high-density component of the EOS (Radice et al.2016a). This effect is, however, only modest for moreconservative choices of the turbulent mixing-length pa-rameter. In these cases, turbulence would not signifi-cantly affect the prospect of detecting phase transitionsin the core of the HMNS using GW observations. How-ever, a definitive statement will have to wait until suf- ficiently resolved GRMHD simulations are available toestimate (cid:96) mix .We have also found that the post-merger GW fre-quency is only weakly affected by the effective turbulentviscosity. Thus, our results provide an important val-idation of the several proposed methods relying on itsmeasure to constrain the EOS of dense nuclear matter(Bauswein & Janka 2012; Takami et al. 2014; Bernuzziet al. 2015a). Our results also reaffirm the observation byBernuzzi et al. (2015a) that the post-merger GW peak-frequency is set at the time of merger. Afterwards, thefrequency stays close to constant and is largely insensi-tive to the evolution of the HMNS, with the exception ofthe signature of BH formation.Finally, our results show that the neutrino signal is alsoinfluenced by the turbulent dissipation of kinetic energyinto heat. The increased temperatures and luminosities,especially for the anti-electron neutrinos, will influencethe proton fraction in the outflows and might have aneffect on the resulting nucleosynthetic yields (Wanajoet al. 2014; Metzger & Fern´andez 2014; Foucart et al.2016). Our results strongly suggest that turbulent dis-sipation will have to be included in the next generationof neutrino-radiation-hydrodynamics models of the out-flows from merging NSs.Here, we presented a first application of the GRLESmethod. In the future, on the one hand, we will ex-tend the present study to more binary configurationsand EOS. On the other hand, work is already underwayto develop closures tuned with highly-resolved GRMHDsimulations of HMNSs. Finally, we will extend GRLESto GRMHD and couple it with a subgrid-scale dynamomodel such as the one of Giacomazzo et al. (2015).It is a pleasure to thank S. Bernuzzi for the many stim-ulating discussions on binary neutron star mergers. Ialso thank S. Hild for the ET-D noise curve data andA. Burrows, P. M¨osta, L. Rezzolla, and C. D. Ott fordiscussions. I gratefully acknowledge support from theSchmidt Fellowship and the Sherman Fairchild Founda-tion. The simulations were performed on Stampede NSFXSEDE (TG-PHY160025), and employed computationalresources provided by the TIGRESS high performancecomputer center at Princeton University, which is jointlysupported by the Princeton Institute for ComputationalScience and Engineering (PICSciE) and the PrincetonUniversity Office of Information Technology.
REFERENCESAndersson, N., & Comer, G. L. 2007, Living Rev. Rel., 10, 1Baiotti, L., & Rezzolla, L. 2016, arXiv:1607.03540Balbus, S. A., & Hawley, J. F. 1991, Astrophys. J., 376, 214Bauswein, A., & Janka, H. T. 2012, Phys. Rev. Lett., 108, 011101Bernuzzi, S., Dietrich, T., & Nagar, A. 2015a, Phys. Rev. Lett.,115, 091101Bernuzzi, S., Radice, D., Ott, C. D., et al. 2015b,arXiv:1512.06397Ciolfi, R., Kastaun, W., Giacomazzo, B., et al. 2017,arXiv:1701.08738Dessart, L., Ott, C., Burrows, A., Rosswog, S., & Livne, E. 2009,Astrophys. J., 690, 1681Duez, M. D., Liu, Y. T., Shapiro, S. L., & Shibata, M. 2006,Phys. Rev., D73, 104015Duez, M. D., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2004,Phys. Rev., D69, 104030Foucart, F., O’Connor, E., Roberts, L., et al. 2016, Phys. Rev.,D94, 123016Gao, H., Zhang, B., & L¨u, H.-J. 2016, Phys. Rev., D93, 044065Germano, M., Piomelli, U., Moin, P., & Cabot, W. H. 1991,Physics of Fluids, 3, 1760Geroch, R. P. 1995, J. Math. Phys., 36, 4226Giacomazzo, B., Zrake, J., Duffell, P., MacFadyen, A. I., & Perna,R. 2015, Astrophys. J., 809, 39 Hanauske, M., Takami, K., Bovard, L., et al. 2016,arXiv:1611.07152Hiscock, W. A., & Lindblom, L. 1985, Phys. Rev., D31, 725—. 1988, Submitted to: Phys. Lett. AKastaun, W., Ciolfi, R., & Giacomazzo, B. 2016, Phys. Rev., D94,044060Kiuchi, K., Cerd´a-Dur´an, P., Kyutoku, K., Sekiguchi, Y., &Shibata, M. 2015, Phys. Rev., D92, 124034Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M., & Wada, T.2014, Phys. Rev., D90, 041502Komatsu, H., Eriguchi, Y., & Hachisu, I. 1989,Mon. Not. Roy. Astron. Soc., 237, 355Kostadt, P., & Liu, M. 2000, Phys. Rev., D62, 023003Lasky, P. D., Haskell, B., Ravi, V., Howell, E. J., & Coward,D. M. 2014, Phys. Rev., D89, 047302Lattimer, J. M., & Swesty, F. D. 1991, Nucl. Phys., A535, 331Lindblom, L. 1996, Annals Phys., 247, 1Majorana, A., & Motta, S. 1985, Journal of Non-EquilibriumThermodynamics, 10, doi:10.1515/jnet.1985.10.1.29. https://doi.org/10.1515%2Fjnet.1985.10.1.29https://doi.org/10.1515%2Fjnet.1985.10.1.29