Turbulence modelling in neutron star merger simulations

Observations of neutron star mergers have the potential to unveil detailed physics of matter and gravity in regimes inaccessible by other experiments. Quantitative comparisons to theory and parameter estimation require nonlinear numerical simulations. However, the detailed physics of energy and momentum transfer between different scales, and the formation and interaction of small scale structures, which can be probed by detectors, are not captured by current simulations. This is where turbulence enters neutron star modelling. This review will outline the theory and current status of turbulence modelling for relativistic neutron star merger simulations.


Introduction
Neutron star (NS) mergers, including both binary neutron star (BNS) and blackhole neutron-star (BHNS) mergers, are connected to some of the most pressing open questions in nuclear and high-energy astrophysics, such as the origin of elements heavier than iron, the nature of matter at supernuclear densities, and the mechanism producing short-gamma ray bursts (Baiotti and Rezzolla, 2017).Multi-messenger observations of NS mergers enabled by ground-based gravitational wave (GW) detectors LIGO, Virgo, and KAGRA have revolutionised our understanding of these events (Abbott et al, 2018(Abbott et al, , 2019)).However, quantitative theoretical frameworks are needed to turn these extraordinary observations into insight.Here, numerical relativity (NR) has a key role.NR simulations have enabled the development of accurate GW waveform models that have been used to constrain the tidal response of NSs (De et al, 2018;Abbott et al, 2019Abbott et al, , 2018)).NR simulations have also been instrumental for the understanding of the electromagnetic (EM) emission from merging NSs (Hotokezaka et al, 2013;Bauswein et al, 2013;Paschalidis et al, 2015;Kyutoku et al, 2015;Sekiguchi et al, 2015;Foucart et al, 2017;Lehner et al, 2016;Ruiz et al, 2016;Sekiguchi et al, 2016;Radice et al, 2018c;Fujibayashi et al, 2018;Nedora et al, 2019;Vincent et al, 2020;Miller et al, 2019;Ciolfi, 2020;Foucart et al, 2020;Nedora et al, 2021;Combi and Siegel, 2023;Kiuchi et al, 2023a,b).They have enabled joint EM and GW inference providing tight constraints on the equation of state (EOS) of NSs and on the r-process nucleosynthesis yields of the mergers (Margalit and Metzger, 2017;Bauswein et al, 2017;Radice et al, 2018d;Dietrich et al, 2020;Breschi et al, 2022;Ruiz et al, 2021;Rosswog and Korobkin, 2022).
Despite the many successes and the rapid progress in the last twenty years, NS merger simulations still face huge challenges.Models of tidally interacting NSs with an order of magnitude smaller systematic errors will be required in the coming years (Gamba et al, 2021), as current GW detectors improve in sensitivity and with the advent of next-generation ground-based GW detectors, such as LIGO Voyager (Berti et al, 2022), NEMO (Ackley et al, 2020), the Einstein Telescope (Punturo et al, 2010), and Cosmic Explorer (Reitze et al, 2019).Quantitative models of the postmerger signal from BNS mergers will also be required.Calculations of the EM counterparts and of the nucleosynthesis yields from these systems need to be dramatically improved to match the expected reduction in the uncertainties in the rates and in the physics of neutron rich nuclei, which are expected from future observational campaigns and upcoming laboratory experiments (Zappa et al, 2023;Schatz et al, 2022).A major obstacle is that simulations not only need to include more sophisticated microphysics, but they also need to be able to capture the chaotic, turbulent, magnetohydrodynamics (MHD) flows that develop after the merger (Radice et al, 2020).
In the context of BNS mergers, hydrodynamic effects become dominant during the last orbit of the two stars (Radice et al, 2020).In the case of binaries with mass ratio q = M 2 /M 1 ≲ 0.75 the least massive star is totally or partially tidally disrupted and accretes on the primary (Bernuzzi et al, 2020;Perego et al, 2022).At higher mass ratios (more symmetric binaries), the stars remain gravitationally bound and impact each other with radial velocities of the order of 0.1 c.At the time of contact or tidal disruption, the orbital velocity of the two centers of mass approaches ∼0.2 c (Radice et al, 2020).The resulting relativistic shear flow is Kelvin Helmholtz (KH) unstable and becomes turbulent (Price and Rosswog, 2006).The cores of the two stars bounce against each other, maybe multiple times, before finally merging to form a massive NS remnant (RMNS), or collapsing to form a black hole (BH).In the process, they produce shock waves that accelerate a few 10 −3 M ⊙ of neutron rich matter to escape velocity.They also squeeze material out of the remnant that was initially in their contact interface and that was previously stirred by the KH.The expelled matter settles in an accretion disk with mass of up to a few 10 −2 M ⊙ (Radice et al, 2018c).In the context of BHNS mergers, if the NS is disrupted outside of the innermost circular orbit (ISCO) of the BH, the part of the material that is not accreted and not gravitationally unbound eventually circularises to form a massive accretion disk (Kyutoku et al, 2021).Shocks and strong shear flows are generated due to the tidal disruption and the self-intersection of the tidal streams.
In both BNS and BHNS mergers, the immediate outcome of the coalescence of the binary components is the formation of a compact object, a BH or a massive NS, surrounded by a turbulent accretion disk.Turbulence in the disk is also driven by the magneotorational instability (MRI) (Balbus and Hawley, 1998).The redistribution of angular momentum operated by turbulence determines the subsequent evolution of these systems, their multi-messenger emissions, and their nucleosynthetic yields.As we argue in Sect.2, direct numerical simulations of such systems are impossible, due to the enormous separation of scales in these flows.Moreover, there are unsolved open questions concerning the nature of turbulence in hot, dense matter that need to be addressed.
So far we have loosely used turbulence in two senses: as a description of how a fluid transfers energy and momentum between different scales, and as a description of how the fluid processes form ordered structures and flows from disordered behaviour.The key features that can be demonstrated for simple models (such as the incompressible non-relativistic Navier-Stokes equations) include the coupling of power between wavenumbers until the physical viscosity exponentially damps the power, and the universal nature of the Kolmogorov cascade showing the power decaying with the wavenumber as k −5/3 .This has been backed up by experimental and numerical studies in Newtonian and relativistic (Radice and Rezzolla, 2013) cases.These points will be expanded in Sects.2-3.
The implications for numerical modelling are profound.If the numerics can resolve the large wavenumber regime where viscosity acts then all the physical consequences of the theory can be captured.This is the regime of direct numerical simulation (DNS).However, if the numerics can only capture part of the wavenumber range where the power is transferred between scales, and not capture the high wavenumber range, then the numerical approximations will interfere with the physical consequences at high wavenumbers and, through the couplings inherent in the nonlinear terms, eventually on all scales.Correcting for these numerical limitations is the regime of large-eddy simulation (LES).
There are many physical and numerical reasons why more realistic models of neutron stars have more complex mechanisms for transferring energy and momentum between scales.In all cases, however, the fundamental problem remains the same: it is impractical in a single numerical simulation to resolve both the gravitational wave scale and all the small scales to which the physics transfers energy and contains sizeable dynamics.The job of turbulence modelling is to replace the idealized model (that may be correct on all scales) with a practical model (that may be incorrect on very small scales, but captures the correct physics when simulated on practical scales).
In this review we will first discuss, in Sect.2, the properties of neutron star matter and the expected impact on turbulence in mergers.We will then outline the mathematical background for constructing models for turbulent fluids in Sect.3, before looking at current results in relativistic fluids in Sect. 4. Finally we will outline some potential future directions for the field in Sect. 5.

Turbulence in hot neutron star matter
As previously anticipated, in BNS mergers turbulence first appears in the contact region between the stars, which is KH-unstable.The KH instability creates vortices with typical scale of the order of the width of the shear layer, corresponding to the fastest growing mode of the instability.These vortical structures are subject to further hydrodynamical instabilities and form progressively smaller flow structures.This is the well known turbulent cascade.Later, turbulence develops in the remnant accretion disk as a result of the MRI.Although some authors have suggested that convection might operate in the merger remnant, neutrino-radiation simulations appear to exclude it (Radice and Bernuzzi, 2023).However, other mechanisms, such as the Tayler-Spruit dynamo, might transport angular momentum and drive turbulence in the remnant (Margalit et al, 2022).In all these regions it is expected that the end result of the turbulent cascade is the generation of a complex, chaotic flow, with structures spanning a large range of scales.The impact of this in a neutron star merger is sketched in Fig. 1, and discussed through the rest of this section.

General considerations
Here, we review the basic phenomenology expected for the turbulent flows that develop in NS mergers.The discussion will only be semi-quantitative and our aim is primarily to identify the relevant length and time scales.We will discuss the more precise mathematical formulation in Sect.3.

Turbulent cascade
The starting point for our discussion are the Navier-Stokes equations for Newtonian incompressible flows, where v is the fluid velocity, ρ is the matter density (assumed to be constant), p the pressure, and ν is the kinematic viscosity, assumed to be constant.These equations Fig. 1 A sketch of the different timescales expected in a neutron star merger.In the inspiral phase both the fluid and magnetic power is contained nearly completely in long scales (small wavenumbers k ∝ 1/ℓ).At merger both are disrupted.Shearing instabilities push power into all fluid lengthscales.This rapidly leads to fully developed Kolmogorov turbulence with the power following a ℓ 5/3 or k −5/3 spectrum.The magnetic field power is rapidly disrupted, with the long scale power following a Kazantsev ℓ −3/2 or k 3/2 spectrum.As the evolution progresses the cascade will develop until both fluid and magnetic field power follows the typical ℓ 5/3 cascade.The cascade cuts off below the turbulent lengthscale, with the power falling off exponentially fast due to viscous damping.Physical mechanisms governing this behaviour and best approximations for the values involved are given in Sect. 2. The merger visualization is adapted from Ciolfi, Binary Neutron Star Mergers After GW170817, Front.Astron.Space Sci., 30 June 2020.Reproduced with permission.
do not apply to the relativistic compressible flows present in NS mergers, but they are sufficient to qualitatively describe the turbulence phenomenology.Moreover, on sufficiently small scales, the matter density can be considered to be constant.Similarly, the velocity fluctuations measured in a locally flat frame comoving with a sufficiently small patch of the flow can be expected to be nonrelativistic and subsonic.In such a local frame, Eq. ( 1) provides a good description of the dynamics as long as magnetic stresses can be neglected, as is the case in the early stages of development of the KH instability.We address the impact of magnetic fields in Sect.2.1.2.From the Navier-Stokes equations we can obtain an evolution equation for the specific kinetic energy.To do so, we multiply Eq. (1a) by v to obtain where v is the magnitude of the velocity v. Integrating the last equation over a volume V , contained in the local fluid patch in which Eqs. (1) can be considered to approximately hold, we obtain where we have used the assumption that ρ ≃ const.After integration by parts, the first integral on the RHS can be rewritten as Because turbulence is expected to be statistically isotropic and homogeneous, the first term on the RHS of Eq. ( 4) is expected to vanish.The second term on the RHS of Eq. ( 4) vanishes due to the incompressibility constraint.In a similar way, we find that In the last term of the RHS we have used the notation A : B to denote the inner product between two tensors: The first term on the RHS vanishes because of isotropy and homogeneity.Since V is arbitrary, we conclude that, for isotropic and homogeneous turbulence, the specific kinetic energy satisfies Let ℓ be a characteristic spatial extent for the local patch of fluid under consideration.For example, ℓ could be taken as the diameter of V .Over times t much smaller than the typical global dynamical timescale, but larger than the local eddy turnover time τ (ℓ) = ℓ/v(ℓ), where v(ℓ) is the typical magnitude of the velocity fluctuations on the scale ℓ, the turbulence can be considered stationary (Pope, 2000).By this, we mean that the velocity fluctuations on scales smaller than ℓ are statistically stationary.That is, In other words, over a time t ≃ τ (ℓ), the energy injection from large scales, which is mediated by the nonlinear term (v•∇)v 2 , achieves a balance with the energy dissipation due to viscosity.This fact is well established by experimental and numerical studies of turbulent flows (Frisch, 1995;Pope, 2000).The magnitude of the energy injection or transport rate is estimated as One of the key assumptions in Kolmogorov's theory of turbulence is that v 3 (ℓ)/ℓ ∼ const, independently of ℓ over a large range of length scales, the so-called inertial range.In other words, v 3 (ℓ)/ℓ is independent of the extent of the region V , if that falls within the inertial range.Often Kolmogorov's law is formulated in terms of the velocity spectrum From Kolmogorov's hypothesis v 2 ∼ ℓ 2/3 .This implies that the spectrum has the famous power law scaling Ê ∼ k −5/3 , see, e.g., Frisch (1995) for a more detailed derivation.According to Kolmogorov's phenomenology, turbulence transports kinetic energy from large scales, at which it is injected, to small scales, where it is converted into internal energy.The cascade proceeds until the scale at which the energy transfer rate is balanced by the viscous energy dissipation rate.The latter can be estimated as so the dissipation scale ℓ d can be estimated by equating Eqs. ( 9) and ( 11), giving where we have taken ℓ 0 to be the largest scale of the inertial range for which v 3 (ℓ)/ℓ ∼ const.This outer scale is typically related to a global scale of the system.For example, in the KH-unstable layer in BNS mergers ℓ 0 ≃ 1 km.In the remnant accretion disk ℓ 0 ≃ Hr, where H ≃ 1/3 is the disk thickness and r is the cylindrical radius.
The quantity in parenthesis in Eq. ( 12) is the inverse of the Reynolds number Re, which measures the relative importance of inertial and viscous terms in the Navier-Stokes equations.Experimental studies of transitional flows, where the balance between inertial and viscous effects changes, show that they typically become turbulent when Re ≳ 1000 (Frisch, 1995), although the value of the critical Reynolds number for transition to turbulence is not universal.This means that, for turbulent flows, ℓ d < 0.005 ℓ 0 and the ratio ℓ 0 /ℓ d grows as Re 3/4 .Such a large range of scales significantly limits simulations that resolve all scales of the system, the so called DNS calculations, to flows with modest Re numbers.Indeed, the cost of DNS calculations scales as Re 3 , when using explicit time integration schemes, such as those employed in numerical relativity.As we will see shortly, the Reynolds numbers encountered in NS mergers are extremely large, so DNS simulations are unfeasible for the foreseeable future.

Magnetohydrodynamical effects
In the presence of magnetic fields, the Navier-Stokes equations are modified due to the appearance of the Lorentz force.In CGS units, the momentum equation becomes where B is the magnetic field.The velocity can still be taken as incompressible, ∇•v = 0.The magnetic field is solenoidal and evolves according to the induction equation: where η is the magnetic resistivity.Eqs. ( 14) and ( 15) are the equations of classical, incompressible MHD.Once again, the classical description is not adequate to describe the dynamics of general-relativistic flows.However, these equations can be used to describe the small-scale flow dynamics in NS merger.
The Lorentz force can be decomposed as The first term on the RHS corresponds to the so-called magnetic tension force, while the second term describes the magnetic pressure.This decomposition motivates the introduction of the plasma β parameter which measures the relative contribution of hydrodynamic and magnetic forces.In the bulk of a BNS merger remnant, at typical densities ∼5 × 10 14 g cm −3 , the pressure arising from repulsive interaction between the nucleons is ∼10 34 dyn cm −2 and increases rapidly with density (Hebeler et al, 2013).This implies that equipartition between the magnetic and fluid pressure in the bulk of BNS merger remnants would require fields strengths significantly in excess of 10 17 G.That is more than two orders of magnitude larger than that in the most magnetized NS known (Olausen and Kaspi, 2014).Even in the accretion disk formed in NS merger remnants, where the typical density is "only" ∼10 12 g cm −3 , the pressure is ∼10 30 dyn cm −2 , implying that magnetic pressure forces can be neglected for fields below 10 15 G. Indeed, the plasma β values found in NS merger remnants in GRMHD simulations are typically larger than ∼100 (Kiuchi et al, 2015b;Palenzuela et al, 2022a), with the exception of magnetically dominated regions in the disk corona and in the jet.
The previous discussion would seem to suggest that magnetic stresses can be neglected when considering the cascade of energy to small scales.However, this ignores two important facts.First, magnetic stresses are not isotropic, so they can transport momentum (and angular momentum) in a shear flow.Such transverse momentum exchange is absent for an ideal fluid, so even a small magnetic field can change the dynamics in a qualitative way.Second, there is an important difference between inertial "forces", which dominate the turbulent cascade, and magnetic stresses.At a given scale ℓ, it is always possible to cancel the contribution of the larger scale flow in the (v • ∇)v term of the momentum equation by using a frame in which the mean velocity vanishes.In other words, inertial forces depend only on the velocity fluctuations at a given scale v(ℓ) and not on the absolute velocity.The Lorentz force, instead, is frame invariant, so it cannot be cancelled with a frame transformation (Beresnyak, 2019).
In the absence of viscosity or resistivity, the hydrodynamic turbulent cascade proceeds until a scale ℓ B at which magnetic tension and inertial forces balance: At smaller scales the turbulent cascade changes character.Magnetic tension inhibits the motion of the flow in the plane perpendicular to the magnetic field direction and turbulence becomes strongly anisotropic.MHD turbulence also exhibit an inverse cascade that can create large-scale magnetic structures.A detailed discussion of MHD turbulence is beyond the scope of this review.We refer the interested reader to Beresnyak (2019); Schekochihin (2021) for a discussion of the current state of the field.See also Biskamp (2003); Kulsrud (2004) for textbook introductions.

Turbulence in RMNSs
For the typical thermodynamic conditions of matter in NS mergers, momentum transfer due to neutrino diffusion provides the dominant source of viscosity.The kinematic neutrino viscosity ν ν can be estimated as (van den Horn and van Weert, 1984;Guilet et al, 2017) where ϵ ν is the average neutrino energy and ⟨λ ν ⟩ is an energy averaged scattering mean free path for neutrinos.In the conditions relevant for the RMNS this yields (Thompson and Duncan, 1993) where 05 is the proton fraction, with n p , n n being the proton and neutron numbers.In the case of the KH unstable interface between the NSs, a typical velocity is given by the relative velocity between the two stars, v 0 ≃ 0.2 c.A typical length scale is given by the width of the shear layer, ℓ 0 ≃ 1 km, so the Reynolds number can be estimated as Accordingly the dissipation scale can be estimated to be However, neutrinos act as an effective source of viscosity only on scales ℓ ≫ λ ν , the latter being the neutrino mean free path.The neutrino mean free path can be estimated as (Thompson and Duncan, 1993) Since λ ν ≫ ℓ ν , neutrino viscosity is not effective at damping the turbulent fluctuations, which continue to cascade to smaller scales.
The dominant source of viscosity on scales smaller than λ ν is electron scattering, which induces an effective kinematic viscosity1 (Thompson and Duncan, 1993) The Reynolds number of the flow is thus much larger than that implied by Eq. ( 22).A better estimate is Re = The corresponding dissipation scale ℓ d is of the order of one nanometer.
The scale at which MHD effects become dominant is This implies that the back-reaction of the magnetic field becomes dominant at scales that are much larger than the dissipation length scale.

Turbulence in the accretion disk
A similar situation is found in the remnant accretion disk.The inner portion of the accretion disk has densities of the order of 10 12 g cm −3 and temperatures around 10 MeV.In such conditions, nucleons are not degenerate and the neutrino viscosity can be estimated as (Keil et al, 1996) The typical scale of the turbulent motion in the disk is ℓ 0 = Hr, where H ≃ 1/3 is the disk thickness and r ≃ 30 km is the typical orbital radius of the densest part of the disk.Using the orbital velocity we find Re ≃ 10 5 .However, as was the case in the KH-unstable region in the BNS context discussed in Sect.2.2, the corresponding viscous scale has ℓ ν ≪ λ ν .The latter can be estimated using Eq. ( 23) to be λ ν ≃ 165 m.On scales smaller than λ ν the turbulent cascade continues all the way to the electron viscosity scale ℓ e ≃ 5 nm.The associated Reynolds number is ∼2 × 10 16 .If the disk is magnetized, as it is likely to be, then magnetic stresses become dominant at scales smaller than As is the case for massive neutron star remnants, magnetic stresses are likely to substantially alter the hydrodynamic turbulent cascade well before the dissipation range has been reached.

Open problems
There are two main ways in which turbulence impacts the postmerger evolution of binary mergers involving NSs.First, turbulence drives the redistribution of the angular momentum in the remnant, driving accretion and mass ejection.Second, turbulence might strongly amplify the initial magnetic field, which, in turn, can drive relativistic outflows and power jets from these systems.
In the case of BNS mergers, it is known that differential rotation is necessary for the stability of the RMNS.However, the impact of angular momentum transport on the RMNS is an open question.As we discuss in Sect.4, depending on the timescale over which the inner core of the RMNS reaches solid body rotation, angular momentum transport could accelerate, or significantly delay, and even possibly avert, the collapse of the remnant.Understanding the role of turbulence on the life time of the remnant is necessary to interpret any constraints on the lifetime of RMNSs.
Another open question concerns the topology of the magnetic fields amplified by the KH and the MRI in these systems.Large-scale fields, with characteristic field lines whose radius of curvature is comparable to the size of the system, could power ultrarelativistic outflows.It would also determine the evolution of RMNSs that do not collapse even after having achieved solid-body rotation.There is reason to believe that the topology of the magnetic field lines after BH formation cannot be too simple.The reason is that GRMHD simulations of BH accretion show that, in the presence of large scale poloidal flux, magnetic flux is accumulated onto the BH until magnetic pressure arrests the accretion flow (Narayan et al, 2003).This is the so-called magnetically arrested disk (MAD) state.MAD accretion flows generate extremely powerful jets, with energy output E K ∼ M disk c2 .In the NS merger case, M disk = 0.1 M ⊙ , so the resulting jets would have a total kinetic energy of up to ∼10 53 erg, which is much larger than the typical (∼10 50 erg) total energy released by short gamma-ray burst (SGRB) engines (Fong et al, 2015;Mooley et al, 2018;Ghirlanda and others, 2019).Tangled fields, or fields with intermittent polarity, would more naturally explain the energetics of SGRBs.

Turbulence Modelling
The previous section shows that ideal (magneto) hydrodynamics may not be the correct theory on all scales, but is likely to be a reasonable theory on scales that can be resolved in numerical simulations.However, the physical transfer of energy-momentum from resolved scales to shorter scales will mean the ideal model will become steadily less accurate as turbulent effects become important.The modelling challenge becomes that of finding an effective theory that captures the impact of the short-scale physics using only information available at the scales we resolve.This may sound impossible, as the information needed is lost: however, Newtonian approaches suggest a range of possibilities.

Newtonian Reynolds equations
To motivate the modelling in relativity, we start with the standard Newtonian approach, following Pope (2000).Assume the equations of hydrodynamics hold at all scales of interest.Assume that we can resolve some, but not all, of the scales (for reasons of complexity, or numerical resolution, for example).We then want to derive the effective equations of motion at scales that we can resolve.
There are two essential viewpoints.One is that the behaviour on the finest possible scales is essentially unknown, or random.In this case we treat the "true" fine scale behaviour as stochastic, where the flow variables are truly random variables.We then get a mean-field theory on the resolved scales by averaging over all possible realisations of the flow variables.
The second viewpoint is that we know the behaviour on the finest possible scales.In this case, resolving some of the scales means averaging over space and/or time to get a mean-field theory that holds on the resolved scales.This has a number of subtle issues in relativity and so will be looked at later.
To illustrate how this works in the Newtonian case, consider again the constant density incompressible Navier-Stokes equations 2 We now want to construct an averaging or filtering operation for the key fields, such as the velocity v which under the operation gives ⟨v⟩.

Statistical averaging
Assume the probability distribution function of the velocity field is f , which depends on space and time.Write the statistical averaging as We can then write the velocity field in terms of its average and fluctuation as Importantly there is a frame in which the fluctuation averages to zero, ⟨v⟩ = 0.The mean field model will give us the equations of motion in terms of the averaged flow variables, which in this simplified case is the velocity field alone.
A crucial assumption to make is that the underlying fields are continuous.At sufficiently short scales this should be reasonable, as the higher-derivative terms associated with viscosity should prevent wave breaking.This assumption guarantees that statistical averaging commutes with differentiation3 , as, for example, This approach is often associated with Reynolds-Averaged Navier-Stokes (RANS), which is central to the development of turbulence models in the Newtonian theory, but not to relativistic applications.

Filtering
The alternative approach is to assume that we are resolving the behaviour on some lengthscale L and to introduce a filtering kernel G(x ′ ; x) that acts on that lengthscale.A typical example would be the Gaussian kernel G = (4πL The averaging operation is then given by This approach is often associated with LES, and has been the base of most approaches used in relativity.

Reynolds equations
Taking the Navier-Stokes equations, substituting in the average/fluctuation form of the velocity, and averaging, we find the Reynolds equations We see that the Reynolds equations match the Navier-Stokes equations except for the term involving the Reynolds stresses ⟨δv ⊗ δv⟩.These stresses act like an additional source of dissipation, and result from the fine-scale behaviour we are unable to resolve.

The closure problem
When evolving the Reynolds equations, we are assuming we have initial data and boundary conditions for the mean flow ⟨v⟩, but no information on the fluctuations δv.Therefore Eqs. ( 35) are not a closed system until we provide information on the fluctuations.This is the closure problem: we need to provide sufficient information on the fine scale behaviour in order to solve the coarse scale problem.As we are solving only for the coarse scale because working at all scales is impractical, this means giving some summary statistics of the correlations that lead to the Reynolds stresses ⟨δv ⊗ δv⟩.
It is useful to decompose the Reynolds stresses further.A standard choice is to define the turbulent kinetic energy Note that dimensionally k TKE is a specific energy, but in the incompressible case this is often not mentioned.The Reynolds stress can then be written as an isotropic and anisotropic part, We can then rewrite all stress terms in the momentum equation (35a) as Thus the isotropic term 2 3 ρk TKE can be absorbed in a modified pressure.Only the anisotropic term a ij is effective in transporting momentum in novel ways.
This decomposition immediately allows us to write down the first solution to the closure problem.The turbulent-viscosity hypothesis of Boussinesq assumes that the anisotropic stress a ij is directly proportional to the mean rate of strain of the average flow, The turbulent viscosity or eddy viscosity ν T is a scalar coefficient.The Reynolds equations immediately become where the effective viscosity is This hypothesis returns the Reynolds equations to exactly the form of the original incompressible Navier-Stokes equations, only with an effective viscosity coefficient and a modified mean pressure.
Experiments in Newtonian flows have shown that the turbulent-viscosity hypothesis may be extremely inaccurate.Nonetheless, the intuition of introducing modified pressures and viscous coefficients to capture the unresolvable small scale physics remains important in most turbulence models.

Relativity via the equations of motion
The steps in the previous sections can be directly applied to the relativistic Euler equations once we have the equations of motion.We illustrate the procedure following Radice (2020) which uses the Valencia approach.The conserved current (usually thought of as the baryon mass current) J a = mn a , and the stress-energy tensor where e is the total energy density and u a the 4-velocity, are decomposed with respect to the normal to the spatial slice N c .The resulting equations of motion are in balance law form The conserved quantities q = (D, S i , E) are projections of the conserved current and the stress-energy tensor along the normal to the t = const.hypersurfaces N c .The fluxes are which contain nonlinear combinations of the conserved variables, whilst the sources are linear in the conserved variables.
It is worth noting that, by reframing the continuity equation ∇ a J a = 0 in terms of the relativistic material derivative m ṅ = mu a ∇ a n we have and so on scales where the flow is stationary in the sense that ṅ ≃ 0, the flow is incompressible in the sense that ∇ a u a .This gives the link to the assumption of incompressibility in the Newtonian context.
The averaging operation can be directly applied to the equations of motion ( 44), but various assumptions are needed to close the system.In Radice (2020) the first assumption is that metric quantities are unaffected by averaging (as the curvature lengthscale is long compared to the fluid scales).The second assumption is needed to close the nonlinear terms in the fluxes, which can be written The closure terms are τ a b which is the analogue of the Reynolds stress introduced in the Newtonian case, and a new term µ i which describes turbulent mass diffusion.The explicit closure is detailed in Sect.3.6.2.
There are two additional nonlinearities that are less obvious and require closure relations.First, the equation of state relates the pressure p to the other thermodynamic potentials.After averaging, we need a closure relation when the functional form of the equation of state is applied to the averaged quantities, as In most approaches it is assumed that Π = 0 as the turbulent fluctuations are subsonic at sufficiently small scales, so corrections to the pressure are expected to be subdominant.This term is present in the Newtonian compressible case.
The second additional closure is in the recovery of the three-velocity v i from the coarse-grained quantities.This is a purely relativistic effect, which is linked to the "correct" way to average the Lorentz factor, and has been discussed in detail by Carrasco et al (2020); Viganò et al (2020); Duez et al (2020); Radice (2020).The practical approaches suggested there all "work" but explicitly break the underlying 4-covariance of the theory.To retain this underlying symmetry requires a different approach.

Relativity via an observer
The previous sections essentially introduced three separate steps.First, an operation (averaging) was introduced that takes fields that vary on all scales to fields that vary only on long scales.Second, the operation is applied to the "true" equations of motion to create the effective field theory.Finally, additional closure relations are provided to complete the effective field theory.
Covariance can be broken in each of these three steps.In averaging the equations of motion produced by the Valencia approach in Sect.3.3 the normal to the spatial slice N c appears implicitly or explicitly in each step, and leads to particularly problems with (for example) the averaging of the Lorentz factor.An alternative viewpoint is that averaging applied to a single slice inherits the dependencies of the slice.
Discussions of how covariance can be retained are given in Eyink and Drivas (2018) and Celora et al (2021).The approach of Celora et al (2021) will be outlined here.Fundamentally it relies on introducing an observer via their 4-velocity vector U a and constructing the averaging operation in a (Fermi-transported) frame orthogonal to that observer.With respect to that observer it can be shown (rather than assumed) that averaging does not affect metric quantities and also commutes with covariant derivatives.The conserved currents and stress-energy tensor can also be decomposed with respect to the observer.Using notation analogous to Sect.3.3, particularly Eq. ( 46), we will find If we can find a "physical" way of constructing the observer 4-velocity (that is, one independent of the normal to the 3 + 1 slices N c ) then the Valencia approach can be used in the standard way with this new set of equations of motion: a tetrad is introduced with timelike leg along N c , and equations ( 48) are projected with respect to the tetrad.As the total stress-energy tensor expressed with respect to the observer U a is full (not diagonal, as the perfect-fluid stress-energy expressed with respect to the fine-grained 4-velocity is), the resulting equations of motion will look like those of a relativistic non-ideal fluid.
The result should be covariant, but has the equivalent complexities and problems of non-ideal relativistic flows.The closure terms µ a and τ a b are now linked to the particle drift and the non-ideal bulk, shear, and heat-transport terms, all with respect to the observer.

Thermodynamics
When decomposing the stress-energy we will have The terms ẽ, p come purely from projecting the stress-energy tensor with respect to this new 4-velocity.The additional τ a b terms can be interpreted as non-ideal terms, illustrating that stresses in the "true" fluid are, generically, not isotropic with respect to the mean flow.
We again see the subtlety in the standard closure relations -the equation of state.The pressure at the fine scale is given by the EOS in terms of the thermodynamical potentials at the fine scale: for example, p = p(n, e).At the averaged scale we might expect ⟨p⟩ = ⟨p⟩ (⟨n⟩ , ⟨e⟩) or ⟨p⟩ = ⟨p⟩ (⟨n⟩ , ẽ).However, there is no guarantee that the averaging process preserves the functional form of the nonlinear closure that is the equation of state.Any differences could be formally included within a bulk viscous pressure Π.However, this highlights that we need to properly account for the thermodynamics under any averaging procedure, in addition to statistical or spacetime filtering effects.
The averaging procedure identifies, at the coarse-grained level, a number density (say ⟨n⟩) and an energy density ẽ of the mean flow.It also identifies a total pressure p.However, we are interpreting the (full) stress-energy tensor as a non-ideal fluid.That means the total pressure comprises a piece in thermodynamic equilibrium plus the bulk viscous correction.
We are free, at some level, to choose the mean EOS (that gives the equilibrium pressure) as we want.We clearly want to relate this as closely as possible to the "true" (micro-scale) EOS, but the coupling between different fluctuations means there will always be some difference.As an illustration of the differences that can arise when averaging, let us assume that the microscale EOS is barotropic, e = e(n), and consider the Gibbs relation p where the chemical potential is µ = de dn .Applying the split into mean and fluctuating pieces we immediately find The interpretation of this equation is that if we use the same microscale EOS applied to the mean flow variables (i.e., use µ(⟨n⟩)) then the average Gibbs relation is not barotropic.The additional term in square brackets can be interpreted as a second parameter and would usually be thought of as sT , the piece corresponding to entropy and temperature.Alternatively it can be thought of as part of the total pressure.There is no clean distinction here.This argument can be extended to more general EOSs.We can (implicitly) specify the coarse-grained EOS through the thermodynamic potential representing the coarsegrained entropy, s = s(ñ, ẽ).In the same way as above, where we noted that we are free to choose the mean (coarse-grained) EOS as we want, we are also free to choose the coarse-grained entropy as we want.From this we can define the coarse-grained temperature and chemical potential to obey the standard definitions The mean Gibbs relation is then where R p holds the residuals (the differences between the mean flow variables and their averages).At this point we have constructed (by assumption) the entropy of the mean flow, so the natural interpretation of R p is as the bulk viscous pressure Π.It should be clear, however, that there is a direct link between the value of R p and the assumed structure of the entropy of the mean flow.

Numerics
We have seen how (statistical) averaging or (spacetime) filtering leads to equations of motion (in the Newtonian case) that contain dissipative terms.In the relativistic compressible case a similar analysis performed at the level of the stress-energy tensor leads to a full, non-ideal "fluid" stress tensor after averaging or filtering, even when the theory at the fine scale is purely ideal.To consider how practical these averaged or filtered theories are, we should look at how non-ideal relativistic theories are numerically evolved.The starting point is the mean stress tensor ⟨T a b ⟩ and any mean charge currents, ⟨n a ⟩ under consideration.Using the standard identity ∇ a V a = (−g) −1/2 ∂ a (−g) 1/2 V a , the equation of motion for the charge currents will follow as Here Γ n is the net creation rate of the charged species, which may arise purely from the averaging or filtering operation.
For the mean stress tensor, the Valencia approach can be used.A tetrad {e b (j) } is introduced with e b (0) = N b timelike.The contractions V a (j) = ⟨T a b ⟩ e b (j) are therefore vectors, and stress-energy conservation and the above identity imply This is a balance law, evolving the conserved variables (−g) 1/2 T 0 b e b (j) in terms of their associated fluxes (−g) 1/2 T i b e b (j) and the geometric source terms which, involving derivatives only of the chosen tetrad, are algebraic in terms of the matter quantities.
Whilst this will give the equations of motion for a generic filtering operation as those for a non-ideal system, this is not the most practical approach for implementation.As illustrated above, the filtering operation splits the stress-energy into a piece that appears to be ideal with respect to the mean flow, and a non-ideal piece, Again, τ a b is the analogue of the Reynolds stress.This means the equations of motion become The "ideal" part of the equations is implemented in a range of standard codes across the community.The additional terms S are those that need to be provided by the turbulent closure.Therefore, it is standard to extend and re-interpret an existing ideal fluid code, rather than implementing a full turbulence model from scratch.Care should be taken when interpreting the results, particularly on the thermodynamics.There remains a cautionary point to note with the numerical implementation.The numerical error can, via the modified equation approach (Sec.3.6.1 and, (e.g., Warming and Hyett, 1974)), be interpreted as the numerical scheme exactly solving a modified set of equations of motion.This means there is an ambiguity in the "nonideal" behaviour observed in simulation results between effects driven by numerical discretisation error and effects driven by the closure relations attempting to capture short-scale physics.Standard grid convergence tests also need careful interpretation, as the closure terms may sensibly vary depending on the resolved lengthscale as well.

Closure relations
There are three key approaches to the closure problem for relativistic flows that need considering.

Implicit LES
The simplest possible turbulence closure is to neglect the turbulent stresses alltogether, i.e., setting τ ab = 0.In this case, one relies on the intrinsic numerical dissipation of high-resolution shock capturing (HRSC) schemes to model the effects of unresolved turbulence fluctuations.This approach, called implicit large-eddy simulation (ILES), might appear to be unlikely to work at first.Nevertheless, ILES has been tremendously successful in many areas of science and engineering (Grinstein et al, 2007).The theoretical foundation of this method relies on the modified equation approach mentioned in Sec.3.5, which we briefly sketch below.We refer to Grinstein et al (2007) for a more systematic discussion.Consider a second-order accurate HRSC discretization of a conservation law (59) Here, second order accurate means that the numerical solution u h = u + O(h 2 ), where h is some measure of the grid spacing4 .It is possible to show that u h is a third-order accurate approximation to the solution of the modified equation where τ b a (v, h) is a nonlinear function of v which depends on the details of the numerical scheme.In other words, u h = v + O(h 3 ) and τ b a (v, h), sometimes called the "numerical viscosity", can be thought of as representing the leading order contribution to the numerical error.The ILES methodology is based on the observation that the τ a b of widely used numerical schemes, such as the piecewise parabolic method (PPM; Colella and Woodward 1984), are similar to explicit closures τ a b , such as the Smagorinsky closure described in Sec.3.6.2.The main advantage of the ILES method is that it is simple to implement (the only requirement is to choose the components of the numerical solver carefully).The main disadvantage is that ILES requires that a significant fraction of the inertial range to be resolved in order to give convergent results, (e.g., Schmidt et al, 2006;Thornber et al, 2007;Schmidt, 2015;Radice et al, 2015).In practice, even though ILES is used in most published NS merger simulations, none of the ILES simulations have been shown to be in a convergent regime, as discussed in more detail in Sec. 4.

Smagorinsky
Smagorinsky (1963) introduced a closure model for Newtonian atmospheric flows that extends the Boussinesq hypothesis to couple the effective anisotropic stresses to the shear of the mean flow.This model has seen great success in Newtonian flows despite its simplicity.Radice (2017); Radice et al (2020) extended this to relativity, choosing the closure to be Here ⊥ c a = δ c a + N c N a projects into the spatial slice, and vc is the 3-velocity of the mean flow in the spatial slice.
Two key features appear immediately.First, the model explicitly depends on a viscosity scalar ν T similar to the kinematic viscosity introduced in the Newtonian case.This has to be provided based on additional physical arguments.Radice chooses, on dimensional grounds, to link ν T = ℓ mix c s , coupling a lengthscale (the mixing length ℓ mix ) with the local speed of sound of the flow.The mixing length can then be linked to the appropriate turbulent lengthscales for the problem as discussed in Sect. 2.
Second, we see that covariance is explicitly broken as the closure relation depends on the gauge through the normal to the spatial slice N c .In principle this is a problematic feature of the closure model -detailed discussions of covariance are in Eyink and Drivas (2018); Celora et al (2021).At present, however, all practical closure implementations share this feature.

Explicit filtering
In the discussion so far we have first assumed we know the short lengthscale behaviour that we then average/filter, only to discard all that knowledge as soon as we evolve the mean flow equations of motion.This aligns with the assumption that the unresolved behaviour is due to physics which we cannot access -for example, we do not have initial data for the fluctuations.
An alternative approach is to assume that all the closure terms come from the filtering operation, and that operation is explicitly known.For example, we could argue that the numerical lengthscale is much larger than the physical dissipation lengthscales, so the filtering due to numerical discretisation dominates.This approach has been extended to relativity for successively more complex models by Viganò et al (2019); Carrasco et al (2020); Palenzuela et al (2022b).
As in Sect.3.1.2we assume we know the filtering kernel G(x ′ ; x) from which we compute the filtered quantity As the Fourier transform of the convolution is a multiplication, we can write the result and its inverse in frequency space as In principle this allows the original quantity to be computed given the coarse grained quantity and the filter.However, this relies on the inverse of the filter existing and being bounded, which cannot be the case when coarse graining in numerical simulations where information is lost.However, F(G) can be approximated by an expansion with a bounded inverse, and this used to relate filtered quantities (see Stolz and Adams 1999 for links to approximate deconvolution methods).
In particular, using a gradient or Taylor series approximation of F(G) (as in Vlaykov et al 2016; Viganò et al 2019) it can be shown that where C depends on the precise form of the filtering kernel.In particular, for the Gaussian kernel G = (4πL we have that C = L/2, linking the correction terms to the lengthscale of the filter.
With these results, the equations of motion in balance law form can be filtered (assuming, following Viganò et al 2019, that the curvature lengthscale is much greater than the filtering and turbulent lengthscales, meaning the metric terms are pulled out of the filtering) as The analogue of the Reynolds stresses is The product rule in Eq. ( 64) allows these terms to be explicitly computed.
If the terms in the fluxes f (j) were simply related to the conserved variables then the procedure would be complete.However, many of the matter terms, such as the velocities and Lorentz factors, do not have this form.For these other matter variables, generalisations of the product rule ( 64) are needed on a case-by-case basis.
In common with the closure relations in Sect.3.6.2covariance is explicitly broken as the filtering is applied to a spatial slice.In contrast to the closure relations in Sect.3.6.2,there are no tuneable parameters except for the form of the filtering kernel.Typically the kernel will depend on a lengthscale L, but this is usually linked to the numerical grid spacing ∆x, not the physical lengthscales.

Angular momentum transport in RMNSs
The impact of angular momentum transport in RMNSs was first investigated by Duez et al (2004).They studied the evolution of differentially rotating NSs with a 2D axisymmetric code that solved the relativistic Navier-Stokes equations.In other words, these modelled turbulence using the equations for a viscous flow.This approach is conceptually similar to the Smagorinsky closure of Sec.3.6.2.However, an important caveat is that, in relativity, the Navier-Stokes equations are mathematically ill-posed (Hiscock and Lindblom, 1985;Kostädt and Liu, 2000), so their solution can manifest unphysical instabilities, if these are not numerically suppressed.The initial angular velocity profile was chosen according to the so-called "j-constant" differential rotation law (Eriguchi and Mueller, 1985).See the t = 0 line in Fig. 2 for an illustration of the typical angular velocity profiles they considered.They found that these NSs are driven towards solid body rotation within a few viscous timescales.This is expected, because uniformly rotating configurations are the minimum energy equilibrium configuration for stars of fixed rest mass M b and angular momentum J (Hartle and Sharp, 1967).This process is illustrated in Fig. 2, which shows the profile of the angular velocity Ω in the equatorial plane for selected times.Configurations with total mass above the maximum limit for rigidly rotating stars underwent catastrophic collapse to form BHs surrounded by massive accretion disks (see Fig. 3).Configurations with mass below this limit were found to settle to a stationary configuration.
As already mentioned, the mathematical formalism used in this first work was problematic, because the relativistic Navier-Stokes equations do not admit a well posed initial value problem (Hiscock and Lindblom, 1985;Kostädt and Liu, 2000).However, the qualitative findings in Duez et al (2004) were later confirmed with ILES simulations in which viscosity emerged self-consistently from physical MHD stresses (Duez et al, 2006b,a;Siegel et al, 2013), in simulations employing a variant of the Israel-Stewart hyperbolic formulation of dissipative relativistic hydrodynamics (HD) (Shibata et al, 2017), and in general-relativistic large-eddy simulation (GRLES) calculations with a Smagorinsky closure (61) (Radice et al, 2018c;Duez et al, 2020).All these studies confirmed that viscosity drives differentially rotating compact objects towards solid body rotation, although there are quantitative differences in the final rotational profiles between the various approaches (Duez et al, 2020).Several of these studies also confirmed that configurations with mass larger than the rigidly-rotating limit collapse to BHs over a few viscous timescales, while lower mass objects are stable.
That said, while these findings are robust, their applicability to BNS mergers is unclear.The reason is that the rotational profiles of RMNSs are qualitatively different from those predicted by the j-constant law.Rotating stars with j-constant profiles have a maximum in the angular velocity at the centre.Centrifugal support is key to support their cores.RMNSs, instead, have slowly rotating inner cores surrounded by envelopes that are almost entirely supported by rotation (contrast Figs. 4 and 2).This difference was first seen in Shibata et al (2005) and then studied in more detail in Kastaun and Galeazzi (2015); Hanauske et al (2017); Kastaun et al (2017); Ciolfi et al (2017).
Because of this difference, the outcome of the viscous evolution of RMNS is not obvious.This is evident from a comparison of the left and right panels of Fig. 5.While viscosity is found to monotonically decrease the life time of hypermassive differentially rotating NSs with j-constant rotation law, the impact on a RMNS is more complex and large viscosity is found to delay the onset of collapse (Radice, 2017).The counter intuitive behavior for large viscosities is due to two effects.First, since the angular velocity has actually a minimum at the center, viscosity can transport angular momentum into the remnant.Second, if the postmerger viscous timescale is shorter than the typical GW emission timescale (J/ J ∼ 20 ms; Bernuzzi et al 2016), then viscosity can suppress the GW emission (Shibata and Kiuchi, 2017;Radice, 2017).The result is to increase the total amount of angular momentum of the remnant, which would have otherwise been radiated in GWs.On the other hand, as we discuss in Sec.4.2, current GRMHD models suggest that the viscous time might not be sufficiently short to impact the GW emission in a qualitative way.Another reason to exclude the most extreme model shown in Fig. 5 is that large viscosity would produce a massive mildlyrelativistic outflow from the thermalization of the tidal streams of the stars shortly before merger (Radice et al, 2018b).This ejecta is expected to produce bright radio emission when interacting with the interstellar medium on a timescale of months to years.However, such a scenario is in tension with the non-detection of late-time radio emission from SGRBs (Schroeder et al, 2020).
In the more likely scenario in which viscous effects become important only after the remnant has become roughly axisymmetric, viscosity could still either favour or disfavour the collapse.On the one hand, viscosity could drive accretion of the envelope onto the core, thereby triggering its collapse.On the other hand, viscosity could prompt mass ejection and/or transport angular momentum into the core of the RMNS, thereby increasing its centrifugal support.Even in the case of systems with total mass below the uniformly rotating limit, the effect of viscosity is not fully understood.This is because RMNSs have angular momentum in significant excess of the maximum that can be supported by uniformly rotating configurations (Radice et al, 2018a).As such, there is no rigidly-rotating equilibrium configuration towards which systems can evolve while conserving mass and angular momentum.In these cases, the formation of rigidly rotating equilibria is likely preceded by a phase of mass ejection, necessary to remove angular momentum from the system (Radice et al, 2018a).All these considerations also ignore thermal effects, which can also have either a positive or negative impact on the stability of the remnant (Radice et al, 2018a;Hammond et al, 2021).
For concreteness, we show in Fig. 6 an example of a binary highlighting some of the issues discussed so far.The figure shows the evolutionary track of the total angular momentum J and rest mass M b of this particular binary.The binary has a total mass above the maximum limit for a uniformly rotating star, however angular momentum transport operated by viscosity (modelled using the GRLES formalism) and hydrodynamic torques of the remnant on the disk drive significant mass outflow from this binary (Nedora et al, 2019(Nedora et al, , 2021)).A naive extrapolation of the simulation data in time, as well as a simplistic analytic estimate based on angular momentum conservation, suggests that the remnant could settle into a stable configuration.However, the ultimate outcome of the evolution is unknown.The binary could not be evolved until collapse, or until the end of its viscous phase.This was in part because of the large computational costs, but also, and more importantly, because the neutrino scheme used in this simulation was not adequate to capture the diffusion of radiation over a timescale comparable to the neutrino cooling timescale (O(1 s)).Longer simulations, with sophisticated neutrino transport and viscosity, will be required to clarify these issues.

Magnetic field amplification
In a seminal paper, Price and Rosswog (2006) first suggested that in BNS mergers even weak initial magnetic fields could be amplified by the KH instability to magnetar level strengths (∼10 15 G).They employed Newtonian smoothed-particle hydrodynamics (SPH) simulations with no explicit turbulence model and their results could not immediately be confirmed by general relativistic (GR) simulations.This stimulated a vigorous effort in the community to better understand the impact of this instability in mergers.The fact that the KH instability is present was investigated in numerical relativity calculations with ILES (Anderson et al, 2008;Baiotti et al, 2008), but the large field amplification could not be immediately replicated (Giacomazzo et al, 2011).We now know that this was because of the insufficient grid resolution of the simulations.Local ILES simulations by Obergaulinger et al (2010) and Zrake and MacFadyen (2013) found that a substantial fraction of the kinetic energy in the shear layer could be converted into magnetic energy during the merger.We remark that the conversion of 10% of the kinetic energy in the shear layer to magnetic energy would be sufficient to produce fields of up to 10 17 G.The production of magnetarlevel fields is thus plausible.The generation of ultra-strong magnetic fields in mergers was finally confirmed by Kiuchi et al (2015bKiuchi et al ( ,a, 2018)), who performed extremely highresolution, global GRMHD ILES simulations of BNS mergers.The development of the KH instability in one of the simulations of Kiuchi et al (2015a) is reproduced in Fig. 7.
Despite the unprecedented, and so far unmatched, grid resolutions, the simulations of Kiuchi et al. did not achieve convergence.Instead, the saturation strength of the magnetic field after the KHI was found to increase monotonically with resolution.The "low" resolution calculations, which still had resolution higher than most published simulations, only showed modest amplification of the magnetic field, consistent with those of previous studies (e.g., Giacomazzo et al, 2011).As the resolution was increased, both the maximum magnetic field and the growth rate were found to increase, with no sign of saturation (see Fig. 8).They reported fields with strength of up to 10 16 G.The lack of convergence is not surprising.As we have discussed in Sec.2.2, even with a field as large as 10 16 G, the back reaction of the magnetic field on the fluid is expected to interrupt the normal hydrodynamic cascade only on scales of centimeters, well beyond the resolvable scales in the simulations.For an infinitesimally thin shear layer, the fastest growing mode of the KHI is infinitesimal (Biskamp, 2003) and simulations need to capture the magnetic field backreaction scale in order to be converged.In more realistic cases, the fastest growing mode corresponds to the width of the shear layer (Biskamp, 2003).In the context of binary NS mergers the relevant length scale is likely to be the density scale height in the outer layers the NSs ρ/∂ r ρ ∼ 100 m so it might be possible to fully resolve the linear phase of the KHI at extremely high resolutions (∆x ∼ 5 m).
One might wonder whether the stringent resolution requirements for the KH instability can be bypassed by starting simulations with initially large (∼10 16 G) magnetic fields.Unfortunately, the results of Kiuchi et al (2015a) do not support this practice.Indeed, as shown in Fig. 9, they found substantial differences between results of simulations employing initially large fields and re-scaled results from calculations that started with more realistic initial field strengths.Such differences do not manifest themselves in the initial growth of the KH instability, when the magnetic field back reaction is unresolved, but they are apparent in the subsequent evolution.These results suggest that nonlinear MHD effects are important after the merger and that realistic initial field values are required for the simulations to be predictive.This point is also stressed by Aguilera-Miret et al (2023) using results from GRLES simulations with the gradient expansion method discussed in 3.6.3.Kiuchi et al (2018) studied the role of the MRI and of magnetic stresses on the evolution of the RMNSs over a timescale of 30 ms after the merger using ILES simulations.In particular, they measured the effective alpha viscosity induced by magnetic stresses.They found that magnetic stresses are most important in the rotationally supported envelope of the RMNS, peaking at densities ρ < 10 13 g cm −3 , where α ∼ 0.01−0.02.The viscosity at higher densities was found to be significantly smaller, α ∼ 0.0005−0.005.Although they cautioned that their results might not be yet converged, the reported values of α at ρ ∼ 10 13 g cm −3 for different resolutions appear to be converging at first order to ∼0.005.On the basis of these results, the impact of magnetic torques on the stability of the RMNS and on its GW signature could be expected to only manifest over very long timescales.
These results are consistent with the fact that the region with ρ > 10 13 g cm −3 is stable against the MRI, since dΩ dr > 0 (see Fig. 4).Radice (2020) and Radice and Bernuzzi (2023) used the values of the viscosity measured by Kiuchi et al (2018) and Kiuchi et al (2023b) to calibrate their Smagorinsky-type subgrid model.They confirmed that the turbulent viscosity does not have a major impact on the stability of the RMNS.In particular, while differences in the GW amplitude were recorded, the post-merger GW peak-frequency was found to be unaffected.
Convergent results of the magnetic field amplification due to the KH instability could only be achieved with the introduction of subgrid models.Giacomazzo et al (2015) pioneered this approach with a subgrid model inspired by mean-field dynamo theory.Their model reproduced the quick amplification of weak fields to more than 10 16 G.Their simulation showed consistent results across resolutions, in contrast with the direct simulations of Kiuchi et al (2015bKiuchi et al ( ,a, 2018)).However, their subgrid model also predicted large field amplification at the surface of the NSs, driven by the large, unphysical velocity gradients present in the artificial atmosphere of their simulations.As such, while promising, their approach still required fine tuning of the model parameters and it was not fully predictive.A similar mean-field dynamo model has also been recently employed by Most and Quataert (2023); Most (2023) to study the production of flares from long-lived RMNS.
Aguilera-Miret et al ( 2022) and Palenzuela et al (2022a) presented the first simulations with a gradient subgrid model that had no tunable parameters, except for the adopted filter function as outlined in Sect.3.6.3.Their simulations confirmed that the KH instability can amplify weak magnetic fields in the stars prior to merger to magnetar-level fields 10 16 G.The statistical properties of the fields appear to have converged (Palenzuela et al, 2022a).For example, Fig. 10, adapted from Palenzuela et al (2022a), shows the average magnetic field strength in different region of the RMNS.As in Kiuchi et al (2018), the maximum strength and growth rate of the magnetic field were found to initially increase as the grid spacing was reduced.However, and in contrast to the ILES simulations of Kiuchi et al (2018), the GRLES simulations appear to achieve convergence at the highest resolutions.
The maximum strength and ratio of toroidal to poloidal field components shortly after merger was found to be independent of the initial magnetic field in the stars, as long as that was sufficiently weak for the back reaction on the fluid, at the characteristic scale of the eddies generated by the KH instability (∼1 km; Sec. 2), to be negligible (Aguilera-Miret et al, 2022).This is illustrated by Fig. 11, which compares kinetic and magnetic energy spectra obtained with different initial magnetic field configurations.These result are consistent with the expectations for a turbulent dynamo.
An important implication is that the poor knowledge of the interior structure of the magnetic field in NSs prior to merger, (e.g., Braithwaite and Spruit, 2006;Lasky et al, 2011;Ciolfi et al, 2011;Bilous et al, 2019;Sur et al, 2022), might not compromise the predictive power of simulations.
Aguilera-Miret et al (2022) and Aguilera-Miret et al (2023) also found tantalizing evidence for an inverse cascade of the magnetic field.In particular, they reported a growth of the typical magnetic field scale from ∼500 m, immediately after merger, to ∼3.5 km, 110 ms after merger.This is illustrated in Fig. 12, reproduced from their second work, which shows the growth of the field scales.It can also be noticed that the magnetic field reaches equipartition with the kinetic energy at its integral length scale, as expected from dynamo theory.However, the field remained predominantly toroidal and no significant angular momentum or viscous effects, which would be mediated by a poloidal field, were reported (Aguilera-Miret et al, 2022).Palenzuela et al (2022a) also discussed the simulations presented by Aguilera-Miret et al (2022).They reported only modest differences in the amplitude of the GW signal between simulations that resolved the growth of the magnetic fields and those with small post-merger field, see Fig. 13.These studies are currently being extended to include neutrino cooling and heating effects (Palenzuela et al, 2022b;Miravet-Tenés et al, 2022).The changes in the gravitational wave peak frequency were even smaller.This suggests that it will be possible to control for uncertainties due to turbulence in future constraints on the dense matter EOS using the postmerger GW spectrum (Bauswein et al, 2016;Bernuzzi, 2020).
More recently, Kiuchi et al (2023b) presented results of ILES GRMHD simulations with approximate neutrino-transport at extremely high-resolution (with ∆x as small as 12.5 m) and found compelling evidence for an αΩ-dynamo (Brandenburg and Subramanian, 2005) operating in the shear layer at the interface between the star and the disk.Accordingly, poloidal flux is generated from the toroidal field by the turbulent flow (the α-effect) and then converted into toroidal field by the rotation (the Ω-effect).However, Kiuchi et al (2023b) had to employ an artificially large initial field (10 15.5 G), in order to resolve the MRI in this region.As such, even though these results suggests that there is an inverse turbulent cascade, similar to those reported by Aguilera-Miret et al (2022, 2023), it is not clear that these studies are describing the same physics.Indeed, the simulations of Aguilera-Miret et al (2022, 2023) show an inverse cascade for the toroidal field, while those of Kiuchi et al (2023b) shows a dynamo producing a poloidal field.Overall, the results of Kiuchi and collaborators suggest that an αΩdynamo is likely to operate in the remnant, but its implications for the dynamics and multi-messenger signals still need to be understood.

Future directions
Turbulence plays an important role in the dynamics of NS mergers.However, due to the enormous separation between the global scale of the system and the dissipation range, direct numerical simulations capturing the full physics of NS mergers are unfeasible.The GRLES methodology replaces the fundamental equations of GRMHD with an effective theory in which the equations are coarse grained.Small scale effects are included in the form of subgrid-scale closures.
GRLES simulations of NS mergers have now been performed by several groups.These simulations confirmed that the KH instability in binary NS mergers generate strong turbulent magnetic fields.A dynamo process involving the turbulent generation of a poloidal field and its conversion into a toroidal field by the differential rotation has been shown to be active in the exterior layers of the RMNS (Kiuchi et al, 2023b).Within this context, turbulent resistivity is needed to rearrange the field lines and create large scale magnetic structures (Aguilera-Miret et al, 2023).The measured effective viscosity due to magnetic stresses is small.This suggests that turbulence will not jeopardize proposed studies of the dense matter EOS with next-generation GW experiments (Kiuchi et al, 2018;Radice, 2020;Palenzuela et al, 2022a;Breschi et al, 2022).However, more studies are needed to fully quantify this statement.On the other hand, it is expected that turbulence and dynamo action will affect the long-term dynamical evolution of the remnant and leave a significant imprint on the associated multi-messenger emission and nucleosynthesis yields.
The GRLES methodology is very well developed for non-relativistic flow, particularly in the unmagnetized limit.Its application to the relativistic, magnetized flows developing in NS mergers is in its infancy.There are conceptual and practical issues that need to be resolved concerning the GRLES method.In particular, we have discussed the issue of covariance of different formulations and their interpretation.Another issue is model validation.Unlike terrestrial cases, experimental testing and validation is completely impractical.Current validation techniques rely on comparison to fine scale simulations.The potential systematic errors arising from whether these fine scale simulations are truly in the DNS regime, or whether the models transfer to a different spacetime (for example), are currently unknown.Moreover, it is not obvious which observables a subgrid GRLES model should be tuned to.It is now standard for models to be validated against the statistics of a (suitably filtered) DNS model.However, it is typical to validate using statistics of the strain-rate or energy transfer.In the astrophysical context these are unobservable; instead we are interested in observables such as the GW signal or the neutrino lightcurves.Quantifying the systematics from the difference between the validated and observed statistics is necessary future work.
More generally, an underlying problem for the future of NS merger simulations is uncertainty quantification.The many uncertainties within the models include the EOS and transport coefficients and aspects of the initial data.There are formal techniques for studying how this uncertainty propagates to impact on observable quantities (reviewed, for example, by Abgrall and Mishra 2017).However, these techniques suffer from the curse of dimensionality: the high dimension of the parameter space makes the quantification uncertain and inaccurate.The turbulent closures introduced by GRLES methods are necessary but introduce many more parameters.This makes a rigorous understanding of the accuracy with which parameters can be estimated from, for example, observed GW signals, much more complex.
We anticipate that in the next few years there will be substantial progress in our understanding of turbulence in NS mergers.Work is currently ongoing on three fronts 1) developing simulations that combine sophisticated microphysics, MHD, and GRLES, (e.g., Palenzuela et al, 2022b;Zappa et al, 2023); 2) improvement in phenomenological subgrid models using local simulations, (e.g., Miravet-Tenés et al, 2022, 2023); and 3) development of model agnostic, data-driven, subgrid models, (e.g., Brunton et al, 2020;Yuan et al, 2020;Karpov et al, 2022).6 Trajectory of a RMNS simulated with the GRLES method in the angular momentum (J) and rest-mass (M b ) plane.The diamond denotes the location of the system at the time at which the GW timescale J/ J becomes much larger than the viscous and cooling timescales.The green line shows the trajectory predicted with a simple analytic model, while the crosses are extrapolations in time past the end of the simulation.The grey shaded area is the set of all uniformly rotating equilibria.Depending on the outcome of the viscous evolution, this system might collapse to a BH, or end up on one of the stable equilibria in the grey region.From Nedora et al., Astrophys. J. 906:98 (2021).Reproduced with permission.

Fig. 2
Fig. 2 Angular velocity in the equatorial plane for a star evolving under the action of viscosity.The time is given in terms of Prot = 2π/Ωc.The viscous timescale for this model is τ vis ≃ 5.5 Prot.From Duez et al., Phys.Rev. D 69:104030 (2004).Reproduced with permission.

Fig. 3 Fig. 4 Fig. 5
Fig. 3 Contours of the rest mass density and velocity profile in the meridional plane for a differentially rotating stars evolving under the action of viscosity.The time is given in terms of Prot = 2π/Ωc.The dark line in the lower right panel denotes the location of the apparent horizon of the formed BH.From Duez et al., Phys.Rev. D 69:104030 (2004).Reproduced with permission.

Fig.
Fig.6Trajectory of a RMNS simulated with the GRLES method in the angular momentum (J) and rest-mass (M b ) plane.The diamond denotes the location of the system at the time at which the GW timescale J/ J becomes much larger than the viscous and cooling timescales.The green line shows the trajectory predicted with a simple analytic model, while the crosses are extrapolations in time past the end of the simulation.The grey shaded area is the set of all uniformly rotating equilibria.Depending on the outcome of the viscous evolution, this system might collapse to a BH, or end up on one of the stable equilibria in the grey region.FromNedora et al., Astrophys.J. 906:98 (2021).Reproduced with permission.

Fig. 7
Fig.7Profile of density and velocity field on the orbital plane for a magnetized BNS merger simulation with a maximum resolution of ∆x = 17.5 m.Adapted fromKiuchi et al., Phys.Rev. D  92:124034 (2015).Reproduced with permission.

Fig. 10
Fig. 10 Evolution of the average field strength in GRLES simulations of magnetized BNS mergers.The figure shows the magnetic field in the bulk of the remnant (ρ > 10 13 g cm −3 ; top panel) and in the envelope (10 10 g cm −3 < ρ < 10 13 g cm −3 ; bottom panel).The dashed lines are for simulations with no subgrid models, while the solid lines show the results obtained with a subgrid model.The different resolutions are ∆x = 120 m (LR), ∆x = 60 m (MR), and ∆x = 30 m (HR).From Palenzuela et al., Phys.Rev. D 106:023013 (2022).Reproduced with permission.

Fig. 11
Fig. 11 Kinetic (solid) and magnetic (dotted) energy spectra from GRLES simulations of magnetized BNS mergers.From left to right, t − tmrg = 5 ms, 10 ms, 20 ms, 30 ms.The different colors correspond to different initial magnetic field configurations, see Aguilera-Miret et al (2022) for the details.The dots show the integral scales for the kinetic and magnetic spectra.From Aguilera-Miret et al., Astrophys.J. Lett.926:L31 (2022) 2. Reproduced with permission.