Binary Neutron Star Merger Simulations with a Calibrated Turbulence Model

Magnetohydrodynamic (MHD) turbulence in neutron star (NS) merger remnants can impact their evolution and multimessenger signatures, complicating the interpretation of present and future observations. Due to the high Reynolds numbers and the large computational costs of numerical relativity simulations, resolving all the relevant scales of the turbulence will be impossible for the foreseeable future. Here, we adopt a method to include subgrid-scale turbulence in moderate resolution simulations by extending the large-eddy simulation (LES) method to general relativity (GR). We calibrate our subgrid turbulence model with results from very-high-resolution GRMHD simulations, and we use it to perform NS merger simulations and study the impact of turbulence. We find that turbulence has a quantitative, but not qualitative impact on the evolution of NS merger remnants, on their gravitational wave signatures, and on the outflows generated in binary NS mergers. Our approach provides a viable path to quantify uncertainties due to turbulence in NS mergers.


Introduction
Binary neutron star (BNS) mergers are prime targets for the ground-based laser interferometric gravitational-wave (GW) detectors LIGO [1], Virgo [2], and KAGRA [3]. BNS mergers generate loud GW signals and can also power bright electromagnetic (EM) transients [4][5][6][7][8][9][10], as demonstrated by the extraordinary multimessenger observations of GW170817 [11,12]. Finally, BNS mergers can eject neutron rich material which subsequently produces heavy elements, such as gold and uranium, through the r-process [4,[13][14][15]. At the time of writing, one more BNS GW event after GW170817 has been announced by the LIGO/Virgo collaboration (LVC): GW190425 [16,17]. However, several more candidates have been reported and are currently being analyzed by the LVC [18]. Many more detections are expected in the next years as GW observatories improve their sensitives and as more facilities are added to the global network of detectors [19].
Submitted to Symmetry, pages 1 -20 www.mdpi.com/journal/symmetry Numerical relativity (NR) simulations are the only tool able to study the dynamics of BNS mergers in the strong field regime and its connection to the multimessenger signals they produce. State-of-the-art NR simulations include a microphysical treatment of dense matter, the impact of weak reactions and neutrino radiation, and magnetic effects [63][64][65][66]. Even though modern simulations ostensibly include all of the physics believed to determine the outcome of BNS mergers, the long-term evolution of binaries after merger remains poorly known, e.g., [67]. Leading sources of uncertainty are connected to our limited knowledge of the behavior of matter at extreme densities and temperatures, the crudeness with which neutrino radiation is treated in the simulations, and our inability to simulate these systems at sufficiently high resolution to resolve the turbulent cascade and for sufficiently long times [68]. This work is part of our ongoing effort to address this last issue.
It is known that the matter flow after merger is subject to a number of magnetohydrodynamics (MHD) instabilities, such as the Kelvin-Helmholtz (KH) instability and the magnetorotational instability (MRI) [69][70][71][72][73][74][75]. These inject turbulence at very small scale and can potentially impact the qualitative outcome of the merger [76][77][78][79]. However, even the best resolved GRMHD simulations to date [65,74,75] cannot capture the scale of the fastest growing mode of the MRI, unless artificially large initial magnetic fields are adopted to increase the cutoff length scales associated with some of these instabilities. Even in these cases, simulations are far from being able to capture the dynamics of the turbulent cascade all the way to the viscous scale, at which neutrino viscosity and drag damps the turbulent eddies [80], as would be required for a DNS simulation.
In Ref. [81] we proposed the general-relativistic large-eddy simulation (GRLES) method as an alternative to performing ultra-high resolution GRMHD simulations. In particular, we proposed to evolve the coarse-grained GRHD equations with a turbulent closure models design to capture the effect of turbulence operating at sub-grid scales. In parallel, a similar, but technically distinct, approach based on the Israel-Steward formalism was proposed by Shibata and collaborators [82]. More recently, a rigorous first-principle theory of relativistic turbulence that, among other things, strengthens the mathematical foundation for the GRLES method, has been proposed by Eyink and Drivas [83]. An extension of the method to GRMHD, taking into account also terms that we neglected in our initial formulation (more on this below) has been proposed in Refs. [84,85]. Rosofsky and Huerta [86] proposed to use machine learning to calibrate subgrid turbulence models for 2D MHD. Finally, a variant of the GRLES method has also been implemented into the SpEC code by the SXS collaboration to perform 2D axisymmetric simulations [66].
The GRLES or viscous approaches are the only way to perform long-term simulations of the postmerger evolution for multiple systems [56,61,62]. However, the results from these simulations inevitably depend on the adopted subgrid model. In earlier work we used turbulence models based on dimensional analysis and linear perturbation theory. Here, we calibrate a subgrid model using results from very high resolution GRMHD simulations performed by Kiuchi and collaborators [75], which were able to resolve all the unstable scales of the MRI for a binary system with an initially large magnetic field. We perform BNS merger simulations with microphysics and compare results obtained with the newly calibrated turbulence model with those obtained using the prescription we proposed in Ref. [81], which was used in several other works [58,61,[87][88][89][90], and to those obtained with traditional GRHD simulations having no subgrid model. We find that turbulence can impact the postmerger evolution of BNSs in a quantitative way, and we discuss the implications for the interpretation of synthetic GW, EM, and nucleosynthesis yields from BNS merger simulations.
The rest of this paper is organized as follows. In Section 2 we review the GRLES formalism and discuss the calibration of the subgrid model. In Section 3 we present our simulation results. Finally, Section 4 is dedicated to discussion and conclusions.

WhiskyTHC
All simulations are performed with the WhiskyTHC code [91][92][93][94]. WhiskyTHC separately evolves the proton and neutron number densities where J µ p,n = n p,n u µ are the proton and neutron four-currents, n p = Y e n is the proton number density, n n is the neutron number density, n = n p + n n is the baryon number density (including baryons in nuclei), u µ the fluid four-velocity, and Y e is the electron fraction of the material. R p = −R n is the net lepton number deposition rate due to the absorption and emission of neutrinos and anti-neutrinos, which is computed using the M0 scheme [50,58].
NS matter is treated as a perfect fluid with stress energy tensor where e is the energy density and p the pressure. We solve the equations for the balance of energy and momentum where Q is the net energy deposition rate due to the absorption and emission of neutrinos, also treated using the M0 scheme. The spacetime is evolved using the Z4c formulation of Einstein's equations [95,96] as implemented in the CTGamma code [97,98], which is part of the Einstein Toolkit [99,100]. CTGamma and WhiskyTHC are coupled using the method of lines. For this work we use the optimal strongly-stability preserving third-order Runge-Kutta scheme [101] as time integrator. Mesh adaptivity is handled using the Carpet mesh driver [102] which implements Berger-Oliger style adaptive mesh refinement (AMR) with subcycling in time and refluxing [103][104][105].

GRLES
According to the Valencia formalism for GRHD [106] the fluid for velocity is decomposed as the sum of a vector parallel and one orthogonal to the t = const hypersurface normal n µ (not to be confused with the neutron and proton number densities) as: where W is the Lorentz factor and v µ is the three velocity. Accordingly, the proton and neutron currents can be written as J µ n,p = n n,p W(n µ + v µ ) =: D n,p (n µ + v µ ) .
In a similar way, the stress energy tensor is decomposed as where are respectively the energy density, the linear momentum density, and the stress tensor in a frame having four-velocity n µ , and γ µν is the spatial metric. With these definitions in place, and neglecting the neutrino source terms to keep the notation simple, the GRHD equations read The GRLES methodology derives a set of equations for the large scale dynamics of the flow, in the sense precisely defined in Ref. [83], by applying a linear filtering operator X → X to derive a set of equations for the coarse grained quantities. For example, the cell averaging done in the context of a finite volume method can be considered as a type of filtering. The averaged equations read: Here, we have implicitly assumed that the metric quantities are unaffected by averaging, because they are already large scale quantities. This is the only approximation made when going from Eqs. (10)- (12) to Eqs. (13)- (15). Although these equations can be considere exact, they are obviously not closed, since not all terms can be expressed solely as a function of the evolved quantities D n,p , S i , and E. This is a manifestation of the nonlinearity of the equations. To close the equations it is necessary to provide a closure for some of the terms. The most obvious terms that need to be closed are the quadratic terms: The correlation terms τ ij and µ i are the subgrid-scale, or turbulent, stress and rest-mass diffusion. We remark that these terms are always present in any numerical discretization of the GRHD equations even if not explicitly included: this is the so-called implicit large-eddy simulation (ILES) approach.
Since the equation of state (EOS) is also non-linear, the filtered pressure is not equal to the EOS evaluated from the coarse-grained quantities, so an additional closure would also be needed when evaluating p, that is: Similarly, the three velocity v i is also a nonlinear function of the evolved coarse grained quantities, so we would need to include a closure also for v i . These terms are treated in full generality in Refs. [84,85], to which we refer for the details. Here, we neglect these corrections, e.g., we assume Π = 0, because we  expect them to be subdominant, since turbulence in the postmerger remnant is subsonic and subrelativistic, meaning that its character should be fully captured by τ ij . This assumption could in principle be verified using GRMHD simulation data. However, the simulations data that we use for calibration [75] is not publicly available, and we only have access to the value of the α parameter, which maps to τ rφ . Consequently, we cannot check the validity of this assumption.
We employ the relativistic extension of the classical turbulence closure of Smagorinsky [107], which we proposed in Ref. [81]: where D i is the covariant derivative associated with γ ij , and ν T , the turbulent viscosity. On the basis of dimensional analysis arguments it is natural to parametrize ν T in terms of a characteristic velocity, the sound speed c s , and a characteristic length scale of turbulence, the mixing length mix , as For MRI-driven turbulence, one can assume mix to be related to the length scale of the most unstable mode of the MRI [77] λ MRI ∼ 20 m Ω 6 rad ms −1 which is the scale at which turbulence is predominantly driven according to linear theory. Accordingly, in our previous work we explored the impact of turbulence by varying mix between 0 and 50 m respectively corresponding to no and very efficient turbulent mixing.
In the context of accretion disk theory, turbulent viscosity is typically parametrized in terms of a dimensionless constant α linked to mix through the relation mix = α c s Ω −1 , where Ω is the angular velocity of the fluid [108]. Recently, Kiuchi and collaborators [75] performed very high resolution GRMHD simulations of a NS merger with sufficiently high seed magnetic fields (10 15 G) to be able to resolve the MRI in the merger remnant and reported averaged α values for different rest-mass density shells. Combining their estimate of α with values of c s and Ω estimated from a simulation performed with mix = 0 we are able to estimate mix as a function of the rest mass density (Fig. 1). We find that the mixing length is well fitted by the expression where Our analysis reveals that mix is relatively small even for the highly-magnetized binary considered by Kiuchi et al. [75]. The peak value of mix estimated from the GRMHD simulations is remarkably close to the analytic prediction given by Eq. (20). We also find that the turbulence weakens at high densities inside the massive NS (MNS) product of the merger. This is expected, because the angular velocity deep inside the remnant grows with radius stabilizing the flow against the MRI [81,[109][110][111][112]. On the other hand, the drop of mix at low density is an artifact of our fitting procedure. Since Kiuchi et al. [75] do not provide the value of α for densities below 10 10 g cm −3 we perform a log-linear extrapolation of α to lower density which results in α becoming zero at the density ρ * . That said, the value of mix at those densities is inconsequential for our simulations, because the orbital period for the part of the disk with density ρ * is comparable to the total postmerger simulation time. Overall, we find that turbulence is strongest in the mantle of the MNS and in the inner part of the disk, at densities between a few times 10 9 g cm −3 to 10 13 g cm −3 .

Models and Simulation Setup
We consider a BNS system with component masses (at infinite separation) M A = M B = 1.35 M . The LS220 EOS [113] is used to describe the nuclear matter. The initial data is constructed with the Lorene code [114], while the evolution is performed with WhiskyTHC using the setup discussed in Ref. [58]. We simulate the same binary multiple times: once with the calibrated mix from Sec. 2.2, and then with fixed constant values for mix : 0, 5 m, 25 m, and 50 m. Additionally, each configuration is run twice: with and without the inclusion of neutrino reabsorption in the simulations. Neutrino cooling is instead always included. The resolution in the finest refinement level of the grid, which covers the NSs during the inspiral and the MNS after merger, is of 185 m. Finally, to quantify finite-resolution effects we also rerun the simulations with no neutrino reabsorption also at the lower resolution of 246 m. The results presented here are thus based on a total of 15 simulations for a total cost of about 3M CPU hours. The simulations with constant mix were already presented 1 in Refs. [58,81,88]. In Ref. [65] postmerger profiles from the mix = 0 no neutrino reabsorption binary was mapped into a high-resolution grid and simulated with the inclusion of a magnetic field. The simulations with calibrated turbulence model are new. For clarity, we 1 However, we rerun the mix = 0 simulation with neutrino reabsorption, which we now continue for a longer time after merger than in our previous work. only include the high-resolution simulations in the figures. If not otherwise specified, the figures refer to the simulations that included both neutrino emission and neutrino reabsorption. The low-resolution data follows the same qualitative trends, although there are quantitative differences.

Qualitative Dynamics
We refer to Refs. [58,115] for a detailed description of the qualitative evolution of the binary considered in this work. Here, we only mention that the simulations span the last ∼3−4 orbits prior to merger and continue for 20−25 ms afterwards. The merger produces a MNS remnant that collapses to black hole (BH) surrounded by a massive accretion torus, typically within the simulation time. Notable exceptions are the simulations with mix = 50 m for which collapse appears to be significantly delayed by viscosity, as we reported in Ref. [81].
The evolution of the maximum rest-mass density for the 5 binaries that did not include neutrino heating is shown in Fig. 2. As previously reported in Ref. [81], we find that the turbulence on the lifetime of the remnant is non monotonic. This is due to the complex interplay between angular momentum transport and suppression of angular momentum losses to GWs operated by the turbulence [81]. However, only the simulation with the largest mixing length ( mix = 50 m), corresponding to very efficient turbulent transport, shows truly significant differences in the contraction rate and in the lifetime of the remnant when compared to the baseline model with mix = 0. In the other cases the changes are quantitative rather than qualitative. This is not surprising: turbulent viscosity plays a small role in the inner core of the MNS because the velocity gradients are relatively small towards the center of the MNS [109][110][111]. These findings are consistent with the results of the GRMHD simulations of the same binary presented in Ref. [65]. There it was found that the inclusion or omission of the magnetic field (and hence of the MRI-induced turbulence) has a modest effect on the collapse time of the remnant, with difference of the same order as those found in our Fig. 2.    The effects of turbulent viscosity are more pronounced in the outer layers of the MNS and in the disk, where velocity gradients are larger. Moreover, mix is maximum in this region according to the calibrated turbulence model. The impact of turbulence on the structure of the disk is shown in Fig. 3, where we report the profiles of entropy, electron fraction, and density for the mix = 0 and mix = mix (ρ) runs (both with neutrino absorption included). The accretion disk is formed in the first milliseconds after the merger, as hot material is expelled from the collisional interface between the NSs. During this phase, turbulence dissipation enhances the thermalization of the flow resulting in the formation of a disk with larger initial entropy and electron fraction compared to that of the baseline model with mix = 0. At later times, turbulence has an opposite effect: turbulent stresses drive the mixing of the hot material in the inner disk with fresh low-entropy material from the mantle of the MNS lowering entropy and electron fraction in the inner part of the disk. Over longer timescales, turbulent angular momentum transport also drives flows of matter to larger radii. This manifests itself in the increase of the density in the midplane of the disk and the smoothing of the isodensity contours in the disk. We remark once again that the internal structure of the remnant is not visibly affected by the turbulence viscosity (with the exception of the mix = 50 m run). The apparent differences in the density isocontours in the MNS in Fig. 3 arise because the MNS has a strong quadrupolar deformation and the mix = 0 and mix = mix (ρ) simulations are dephased with respect to each other.

Gravitational Waves
The impact of turbulent viscosity on the GW signal is shown in Fig. 4. The most prominent feature of the postmerger spectrum is a peak at f GW 3 kHz associated with the rotational period of the quadrupolarly deformed MNS [78,[116][117][118]. We find this feature to be robust against turbulence: deviations in the peak frequency with mix are typically small and of the same order of the nominal uncertainty  Outflow rate for the baseline and the calibrated turbulence models. The thin lines denote simulations that did not include neutrino reabsorption. For clarity we smooth the data using a rolling average with amplitude 0.5 ms. Turbulent dissipation has a modest impact on the dynamical ejecta mass and is subdominant in comparison to neutrino heating. of the Fourier transform of the time domain data. Thus, our results confirm that the measurement of the postmerger peak frequency is a promising avenue to constrain the EOS of dense matter using 3rd generation detectors, e.g., Ref. [119]. Models with turbulent viscosities larger than those of our calibrated model, i.e., mix = 25 m and mix = 50 m, also show an overall decrease in the power of the GW signal, suggesting that turbulence might suppress GW emission, as also found in Refs. [81,120]. That said, significant reduction in the GW signal as found by Ref. [120] would require turbulent viscosities significantly larger than those estimated from GRMHD simulations [75].
Turbulent angular momentum transport instead has a significant impact on the secondary features of the postmerger GW spectrum. In particular, we find that turbulence can shift and amplify secondary peaks in the spectrum. Such features in the spectrum have been proposed as a possible signature of first order phase transitions [121]. However, our simulations show that turbulence might produce the same qualitative effect on the GW spectrum. Consequently, we caution that the search for new physics in the postmerger signal must account for the uncertainties related to the development of turbulence in the MNS. Follow up studies are necessary to precisely quantify them.

Outflows
Tidal interaction between the stars prior to merger and shocks after merger drive the ejection of neutron rich matter as the NSs coalesce [68,122]. In Ref. [88] we studied the impact of turbulent viscosity on the dynamical ejection of matter during BNS mergers. We found that turbulence can boost the ejection of matter in asymmetric binaries, but has only a small impact on the mass ejection from comparable mass binaries, such as the system considered in this study. Here, we confirm that these results hold also when using a calibrated turbulent viscosity model. In Fig. 5 we show the outflow rate of unbound matter (with u t ≤ −1) across a coordinate sphere with radius R = 200 G/c 2 M 295 km. The differences in the overall ejecta mass are not large, considering the large numerical uncertainties associated with the calculation of the ejecta [58]. However, there is a clear trend in the outflow data. On the one hand, turbulent dissipation does not effect the outflow rate at early time, when fast material accelerated when the remnant rebounds reaches the detection sphere. On the other hand, the mix > 0 runs have significantly larger outflow rates at later times, when the slower part of the dynamical ejecta crosses the detection sphere. This is visible for the mix = mix (ρ) simulations in Fig. 5. Simulations with other values of mix follow the same trend, but are not included to avoid cluttering the figure. This late time boost of the outflow rate is likely the result of the increased thermalization of the flow due to turbulent viscosity and the resulting larger pressure gradients. The effect of turbulent viscosity is, however, subdominant compared to the influence of neutrino reabsorption, which are found to play a key role in driving the "second wave" of the dynamical ejecta. Indeed, as shown in Fig. 5, this part of the outflow is almost completely absent when neutrino reabsorption is switched off in the simulations. This second component of the dynamical ejecta is predominantly constitute of neutrino-irradiated material at high latitudes above the merger remnant.
Both viscosity and neutrino reabsorption play an important role in the determination of the electron fraction Y e of the ejecta [48][49][50]55,58], as shown in Fig. 6. This quantity is of particular interest since it most directly determines the outcomes of the r-process nucleosynthesis for the thermodynamic conditions typical of NS merger ejecta [123]. Turbulent dissipation increases the temperatures in the outflows activating pair capture processes which drive Y e to higher values through the reaction e + + n → p +ν e . This effect is particularly pronounced when considering the tidal tail, which corresponds to the low-Y e peak in the outflow distribution. Neutrino absorption has an even larger effect via the reaction ν e + n → p + e − . Neutrinos generate relatively high Y e ejecta, especially at high latitudes [50]. They also effect the tidal tail, which is irradiated after colliding with the faster shock-accelerated outflows that are launched after the merger [58].
Turbulent angular momentum transport in the remnant accretion disk and neutrino heating are expected to power additional outflows on a timescale of a few seconds [47,56,57,61,62,[124][125][126][127][128][129][130][131]. These secular ejecta are expected to dominate the overall nucleosynthesis and electromagnetic signal from BNS mergers [58,132]. We considered the role of turbulent viscosity on the long term mass ejection in Ref. [61]. However, due to the large computational costs, we could not systematically vary mix in Ref. [61]. Unfortunately, the simulations we present here do not span a sufficiently long time to be able to study the secular ejecta, so we leave this investigation to future work.

Discussion
MHD instabilities active in BNS merger remnants drive turbulence at many different scales [74,75]. Turbulence can generate large scale magnetic fields with potentially dramatic consequences for the evolution of the MNS remnants and its EM emissions [65,71,79]. Due to the large Reynolds number in the flow, directly capturing all scales of the turbulent flow in the MNS is beyond the reach of numerical simulations for the foreseeable future.
In Ref. [81] we proposed a scheme to include subgrid scale turbulence effects into global simulations and study the associated uncertainties on the multimessenger signatures of BNS mergers. Our method extends the LES methodology to GR. We derived evolution equations for coarse the grained fluid number, momentum, and energy densities. These equations are exact, but are not closed, so a closure must be provided. This closure represents the effect of small scale (subgrid) turbulence on the evolution of the large scale quantities. Here, we proposed to use a closure based on the classical turbulent viscosity ansatz [107], which we calibrated against very-high-resolution GRMHD simulations from Ref. [75]. We performed BNS merger simulations with microphysics and neutrinos using the newly proposed turbulence model. We showed that our scheme is robust and gives sensible results. We compared simulations performed with the newly calibrated scheme with simulations performed either with no subgrid model, or with a simpler scheme in which turbulent viscosity is assumed to be a constant fixed on the basis of dimensional analysis arguments.
Our results show that subgrid turbulence has a quantitative impact on the evolution of the remnant. Turbulence can affect the lifetime of the remnant, although large differences in the postmerger evolution are only found for values of the turbulent viscosity that are significantly larger than those found in GRMHD simulations. The peak frequency of the postmerger GW spectrum is found to be insensitive to turbulence, but secondary features in the spectrum and the overall luminosity in the postmerger are affected by turbulent dissipation. Turbulent angular momentum transport and dissipation alter the structure and composition of the remnant disk and the amount and composition of the dynamical ejecta, potentially impacting r-process nucleosynthesis yields and EM counterparts. However, turbulence is found to be subdominant when compared to neutrino effects.
More simulations are needed to establish the degree to which systematic uncertainties due to turbulence will limit our ability to search for new physics, such as phase transitions, in multimessenger observations of BNS mergers. The method we propose here can be and has been extended to the full-GRMHD equations [84,85]. Including GRMHD in our simulations is crucial to capture also the large scale effects due magnetic fields, such as jet launching [37,65], which are presently not included. In parallel, local high-resolution simulations should be performed to develop and calibrate turbulence models. These are all objectives of our future work. Acknowledgments: It is a pleasure to thank S. Bernuzzi for discussions that motivated this work, A. Prakash for carefully proofreading the manuscript, and S. Hild for the ET-D noise curve data.