Binary Neutron Star Merger Simulations with a Calibrated Turbulence Model

Abstract

Magnetohydrodynamic (MHD) turbulence in neutron star (NS) merger remnants can impact their evolution and multi-messenger 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.

1. 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 multi-messenger 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].

Multi-messenger observations of BNS mergers are starting to constrain the poorly-known properties of matter at extreme densities [11,12,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36] and the physical processes powering short $\gamma$-ray bursts (SGRBs) [37,38,39,40,41,42]. They are also beginning to reveal the role played by compact binary mergers in the chemical enrichment of the galaxy with r-process elements [8,13,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62]. The key to the solution of some of the most pressing open problems in nuclear and high-energy astrophysics—such as the origin of heavy elements, the nature of neutron stars (NSs), and the origin of SGRBs—is encoded in these and future observations. However, theory is essential to turn observations into answers.

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 multi-messenger 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 direct Navier Stokes (DNS) simulation.

In [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 model design to capture the effect of turbulence operating at subgrid scales. In parallel, a similar, but technically distinct, approach based on the Israel–Stewart formalism was proposed by Shibata and collaborators [82]. More recently, a rigorous first-principles theory of relativistic turbulence that, among other things, strengthens the mathematical foundation for the GRLES method was 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) was proposed in [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 was also implemented in 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 the results obtained with the newly calibrated turbulence model with those obtained using the prescription we proposed in [81], which was used in several other works [58,61,87,88,89,90], and with 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 the discussion and conclusions.

2. Methods

2.1. WhiskyTHC

All simulations were performed with the WhiskyTHC code [91,92,93,94]. WhiskyTHC separately evolves the proton and neutron number densities:

$${\nabla }_{\mu }\left(J_{p,n}^{\mu }\right)=R_{p,n},$$

(1)

where $J_{p,n}^{\mu }=n_{p,n}{u}^{\mu }$ 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^{\mu }$ is 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:

$$T_{\mu \nu }=(e+p)u_{\mu }{u}_{\nu }+p{g}_{\mu \nu },$$

(2)

where e is the energy density and p the pressure. We solve the equations for the balance of energy and momentum:

$${\nabla }_{\nu }{T}^{\mu \nu }=Q{u}^{\mu },$$

(3)

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 the 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].

2.2. 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=\mathrm{const}$ hypersurface normal $n^{\mu }$ (not to be confused with the neutron and proton number densities) as:

$$u^{\mu }=(-u_{\mu }{n}^{\mu })(n^{\mu }+v^{\mu })=:W(n^{\mu }+v^{\mu }),$$

(4)

where W is the Lorentz factor and $v^{\mu }$ is the three velocity. Accordingly, the proton and neutron currents can be written as:

$$J_{n,p}^{\mu }=n_{n,p}W(n^{\mu }+v^{\mu })=:D_{n,p}(n^{\mu }+v^{\mu }).$$

(5)

In a similar way, the stress energy tensor is decomposed as:

$$T_{\mu \nu }=E{n}_{\mu }{n}_{\nu }+S_{\mu }{n}_{\nu }+S_{\nu }{n}_{\mu }+S_{\mu \nu },$$

(6)

where:

$$\begin{array}{cc}\hfill & E=T_{\mu \nu }{n}^{\mu }{n}^{\nu }=(e+p)W^2-p,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & S_{\mu }=-{\gamma }_{\mu \alpha }{n}_{\beta }{T}^{\alpha \beta }=(e+p)W^2{v}_{\mu },\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & S_{\mu \nu }={\gamma }_{\mu \alpha }{\gamma }_{\nu \beta }{T}^{\alpha \beta }=S_{\mu }{v}_{\nu }+p{\gamma }_{\mu \nu },\hfill \end{array}$$

(9)

are respectively the energy density, the linear momentum density, and the stress tensor in a frame having four-velocity $n^{\mu }$, and ${\gamma }_{\mu \nu }$ is the spatial metric.

With these definitions in place and neglecting the neutrino source terms to keep the notation simple, the GRHD equations read:

$$\begin{array}{cc}\hfill {\partial }_t\left(\sqrt{\gamma }{D}_{n,p}\right)+{\partial }_j[\alpha \sqrt{\gamma }\left(v^j+n^j)D_{n,p}\right]& =0,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\partial }_t\left(\sqrt{\gamma }{S}_i\right)+{\partial }_j\left[\alpha \sqrt{\gamma }\left({S_i}^j+S_i{n}^j\right)\right]& =\alpha \sqrt{\gamma }\left(\frac{1}{2}{S}^{jk}{\partial }_i{\gamma }_{jk}+\frac{1}{\alpha }{S}_k{\partial }_i{\beta }^k-E{\partial }_{i}log\alpha \right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\partial }_t\left(\sqrt{\gamma }E\right)+{\partial }_j\left[\alpha \sqrt{\gamma }\left(S^j+E{n}^j\right)\right]& =\alpha \sqrt{\gamma }\left(K_{ij}{S}^{ij}-S^i{\partial }_{i}log\alpha \right).\hfill \end{array}$$

(12)

The GRLES methodology derives a set of equations for the large-scale dynamics of the flow, in the sense precisely defined in [83], by applying a linear filtering operator $X\mapsto \overline{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:

$$\begin{array}{cc}\hfill {\partial }_t\left(\sqrt{\gamma }\overline{D_{n,p}}\right)+{\partial }_j[\alpha \sqrt{\gamma }\left(\overline{D_{n,p}{v}^j}+\overline{D}{n}^j)\right]& =0,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\partial }_t\left(\sqrt{\gamma }\overline{S_i}\right)+{\partial }_j\left[\alpha \sqrt{\gamma }\left(\overline{{S_i}^j}+\overline{S_i}{n}^j\right)\right]& =\alpha \sqrt{\gamma }\left(\frac{1}{2}\overline{S^{jk}}{\partial }_i{\gamma }_{jk}+\frac{1}{\alpha }\overline{S_k}{\partial }_i{\beta }^k-\overline{E}{\partial }_{i}log\alpha \right),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\partial }_t\left(\sqrt{\gamma }\overline{E}\right)+{\partial }_j\left[\alpha \sqrt{\gamma }\left(\overline{S^j}+\overline{E}{n}^j\right)\right]& =\alpha \sqrt{\gamma }\left(K_{ij}\overline{S^{ij}}-\overline{S^i}{\partial }_{i}log\alpha \right).\hfill \end{array}$$

(15)

Here, we implicitly assume that the metric quantities are unaffected by averaging, because they are already large-scale quantities. This is the only approximation made when going from Equations (10)–(12) to Equations (13)–(15). Although these equations can be considered exact, they are obviously not closed, since not all terms can be expressed solely as a function of the evolved quantities $\overline{D_{n,p}},\overline{S_i},$ and $\overline{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:

$$\begin{array}{ccc}\hfill \overline{S_{ij}}=\overline{S_i}\overline{v_j}+\overline{p}{\delta }_{ij}+{\tau }_{ij},& \hfill & \hfill \overline{D{v}^i}=\overline{D}\overline{v^i}+{\mu }^i,\end{array}$$

(16)

The correlation terms ${\tau }_{ij}$ and ${\mu }^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 $\overline{p}$, that is:

$$\overline{p}=p(\overline{D_{n,p}},\overline{S_i},\overline{E})+\Pi$$

(17)

Similarly, the three-velocity $\overline{v^i}$ is also a nonlinear function of the evolved coarse-grained quantities, so we would need to include a closure also for $\overline{v^i}$. These terms were treated in full generality in [84,85], to which we refer for the details. Here, we neglect these corrections, e.g., we assume $\Pi =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 ${\tau }_{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 $\alpha$ parameter, which maps to ${\tau }_{r\varphi }$. 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 [81]:

$$\begin{array}{ccc}\hfill {\tau }_{ij}=-2{\nu }_T(e+p)W^2\left[\frac{1}{2}\left(D_i\overline{v_j}+D_j\overline{v_i}\right)-\frac{1}{3}{D}_k\overline{v^k}{\gamma }_{ij}\right],& \hfill & \hfill {\mu }^i=0,\end{array}$$

(18)

where $D_i$ is the covariant derivative associated with ${\gamma }_{ij}$ and ${\nu }_T$ the turbulent viscosity. On the basis of dimensional analysis arguments, it is natural to parametrize ${\nu }_T$ in terms of a characteristic velocity, the sound speed $c_s$, and a characteristic length scale of turbulence, the mixing length ${\ell }_{\mathrm{mix}}$, as: 

$${\nu }_T={\ell }_{\mathrm{mix}}{c}_s.$$

(19)

For MRI-driven turbulence, one can assume ${\ell }_{\mathrm{mix}}$ to be related to the length scale of the most unstable mode of the MRI [77]:

$${\lambda }_{\mathrm{MRI}}\sim 20\mathrm{m}{\left(\frac{\Omega }{6\mathrm{rad}{\mathrm{ms}}^{-1}}\right)}^{-1}\left(\frac{B}{{10}^{15}\mathrm{G}}\right),$$

(20)

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 ${\ell }_{\mathrm{mix}}$ between 0 and $50\mathrm{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 $\alpha$ linked to ${\ell }_{\mathrm{mix}}$ through the relation ${\ell }_{\mathrm{mix}}=\alpha c_s{\Omega }^{-1}$, where $\Omega$ is the angular velocity of the fluid [108]. Recently, Kiuchi and collaborators [75] performed very-high-resolution GRMHD simulations of an NS merger with sufficiently high seed magnetic fields $\left({10}^{15}\mathrm{G}\right)$ to be able to resolve the MRI in the merger remnant and reported averaged $\alpha$ values for different rest-mass density shells. Combining their estimate of $\alpha$ with the values of $c_s$ and $\Omega$ estimated from a simulation performed with ${\ell }_{\mathrm{mix}}=0$, we are able to estimate ${\ell }_{\mathrm{mix}}$ as a function of the rest-mass density (Figure 1). We find that the mixing length is well fitted by the expression:

$${\ell }_{\mathrm{mix}}=\left\{\begin{array}{cc}a\xi exp\left(-{\left|b\xi \right|}^{5/2}\right)\left[\mathrm{m}\right],\hfill & \mathrm{if} \xi >0,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(21)

where:

$$\xi ={log}_{10}\left(\frac{m_b(n_p+n_m)}{{\rho }^{*}}\right),$$

(22)

where $m_b$ is the atomic mass unit in grams, $a=22.31984$, $b=-0.4253832$, and ${\rho }^{*}=1.966769\times {10}^9\mathrm{g}{\mathrm{cm}}^{-3}$. This fit and the constant values of ${\ell }_{\mathrm{mix}}$ used in our previous studies are shown in Figure 1.

Our analysis reveals that ${\ell }_{\mathrm{mix}}$ is relatively small even for the highly-magnetized binary considered by Kiuchi et al. [75]. The peak value of ${\ell }_{\mathrm{mix}}$ estimated from the GRMHD simulations is remarkably close to the analytic prediction given by Equation (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 the radius stabilizing the flow against the MRI [81,109,110,111,112]. On the other hand, the drop of ${\ell }_{\mathrm{mix}}$ at low density is an artifact of our fitting procedure. Since Kiuchi et al. [75] did not provide the value of $\alpha$ for densities below ${10}^{10}\mathrm{g}{\mathrm{cm}}^{-3}$, we perform a log-linear extrapolation of $\alpha$ to lower density, which results in $\alpha$ becoming zero at the density ${\rho }^{*}$. That said, the value of ${\ell }_{\mathrm{mix}}$ at those densities is inconsequential for our simulations, because the orbital period for the part of the disk with density ${\rho }^{*}$ 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\mathrm{g}{\mathrm{cm}}^{-3}$ to ${10}^{13}\mathrm{g}{\mathrm{cm}}^{-3}$.

2.3. Models and Simulation Setup

We considered a BNS system with component masses (at infinite separation) $M_A=M_B=1.35M_{\odot }$. The LS220 EOS [113] was used to describe the nuclear matter. The initial data were constructed with the Lorene code [114], while the evolution was performed with WhiskyTHC using the setup discussed in [58]. We simulated the same binary multiple times: once with the calibrated ${\ell }_{\mathrm{mix}}$ from Section 2.2 and then with fixed constant values for ${\ell }_{\mathrm{mix}}$: 0, 5 m, 25 m, and 50 m. Additionally, each configuration was run twice: with and without the inclusion of neutrino reabsorption in the simulations. Neutrino cooling was instead always included. The resolution in the finest refinement level of the grid, which covered the NSs during the inspiral and the MNS after merger, was $185\mathrm{m}$. Finally, to quantify finite-resolution effects, we also reran the simulations with no neutrino reabsorption also at the lower resolution of $246\mathrm{m}$. The results presented here were thus based on a total of 15 simulations for a total cost of about 3M CPU hours. The simulations with constant ${\ell }_{\mathrm{mix}}$ were already presented (however, we reran the ${\ell }_{\mathrm{mix}}=0$ simulation with neutrino reabsorption, which we now continued for a longer time after merger than in our previous work) in [58,81,88]. In [65], postmerger profiles from the ${\ell }_{\mathrm{mix}}=0$ no neutrino reabsorption binary were mapped into a high-resolution grid and simulated with the inclusion of a magnetic field. The simulations with the calibrated turbulence model were new. For clarity, we only included 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 followed the same qualitative trends, although there were quantitative differences.

3. Results

3.1. Qualitative Dynamics

We refer to [58,115] for a detailed description of the qualitative evolution of the binary considered in this work. Here, we only mention that the simulations spanned the last $\sim 3-4$ orbits prior to merger and continued for $20-25\mathrm{ms}$ afterwards. The merger produced an MNS remnant that collapsed to a black hole (BH) surrounded by a massive accretion torus, typically within the simulation time. Notable exceptions were the simulations with ${\ell }_{\mathrm{mix}}=50\mathrm{m}$ for which the collapse appeared to be significantly delayed by viscosity, as we reported in [81].

The evolution of the maximum rest-mass density for the five binaries that did not include neutrino heating is shown in Figure 2. As previously reported in [81], we found that the turbulence on the lifetime of the remnant was non-monotonic. This was 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 (${\ell }_{\mathrm{mix}}=50\mathrm{m}$), corresponding to very efficient turbulent transport, showed truly significant differences in the contraction rate and in the lifetime of the remnant when compared to the baseline model with ${\ell }_{\mathrm{mix}}=0$. In the other cases, the changes were quantitative rather than qualitative. This was not surprising: turbulent viscosity played a small role in the inner core of the MNS because the velocity gradients were relatively small towards the center of the MNS [109,110,111]. These findings were consistent with the results of the GRMHD simulations of the same binary presented in [65]. There, it was found that the inclusion or omission of the magnetic field (and hence of the MRI-induced turbulence) had a modest effect on the collapse time of the remnant, with the difference of the same order as those found in our Figure 2.

The effects of turbulent viscosity were more pronounced in the outer layers of the MNS and in the disk, where velocity gradients were larger. Moreover, ${\ell }_{\mathrm{mix}}$ was maximum in this region according to the calibrated turbulence model. The impact of turbulence on the structure of the disk is shown in Figure 3, where we report the profiles of entropy, electron fraction, and density for the ${\ell }_{\mathrm{mix}}=0$ and ${\ell }_{\mathrm{mix}}={\ell }_{\mathrm{mix}}\left(\rho \right)$ runs (both with neutrino absorption included). The accretion disk was formed in the first milliseconds after the merger, as hot material was expelled from the collisional interface between the NSs. During this phase, turbulence dissipation enhanced 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 ${\ell }_{\mathrm{mix}}=0$. At later times, turbulence had the opposite effect: turbulent stresses drove the mixing of the hot material in the inner disk with fresh low-entropy material from the mantle of the MNS, lowering the entropy and electron fraction in the inner part of the disk. Over longer timescales, turbulent angular momentum transport also drove flows of matter to larger radii. This manifested 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 was not visibly affected by the turbulence viscosity (with the exception of the ${\ell }_{\mathrm{mix}}=50\mathrm{m}$ run). The apparent differences in the density isocontours in the MNS in Figure 3 arose because the MNS had a strong quadrupolar deformation, and the ${\ell }_{\mathrm{mix}}=0$ and ${\ell }_{\mathrm{mix}}={\ell }_{\mathrm{mix}}\left(\rho \right)$ simulations were dephased with respect to each other.

3.2. Gravitational Waves

The impact of turbulent viscosity on the GW signal is shown in Figure 4. The most prominent feature of the postmerger spectrum was a peak at $f_{\mathrm{GW}}\simeq 3\mathrm{kHz}$ associated with the rotational period of the quadrupolarly-deformed MNS [78,116,117,118]. We found this feature to be robust against turbulence: deviations in the peak frequency with ${\ell }_{\mathrm{mix}}$ were typically small and of the same order of the nominal uncertainty of the Fourier transform of the time domain data. Thus, our results confirmed that the measurement of the postmerger peak frequency was a promising avenue to constrain the EOS of dense matter using third-generation detectors, e.g., see [119]. Models with turbulent viscosities larger than those of our calibrated model, i.e., ${\ell }_{\mathrm{mix}}=25\mathrm{m}$ and ${\ell }_{\mathrm{mix}}=50\mathrm{m}$, also showed an overall decrease in the power of the GW signal, suggesting that turbulence might suppress GW emission, as also found in [81,120]. That said, significant reduction in the GW signal as found by [120] would require turbulent viscosities significantly larger than those estimated from GRMHD simulations [75].

Turbulent angular momentum transport instead had a significant impact on the secondary features of the postmerger GW spectrum. In particular, we found that turbulence could shift and amplify secondary peaks in the spectrum that were formed at early time after merger [121]. Such features in the spectrum were analogous to those due to first order phase transitions [122]. 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 quantify them precisely.

3.3. 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,123]. In [88], we studied the impact of turbulent viscosity on the dynamical ejection of matter during BNS mergers. We found that turbulence could boost the ejection of matter in asymmetric binaries, but had only a small impact on the mass ejection from comparable mass binaries, such as the system considered in this study. Here, we confirmed that these results held also when using a calibrated turbulent viscosity model. In Figure 5, we show the outflow rate of unbound matter (with $u_t\le -1$) across a coordinate sphere with radius $R=200G/c^2{M}_{\odot }\simeq 295\mathrm{km}$. The differences in the overall ejecta mass were not large, considering the large numerical uncertainties associated with the calculation of the ejecta [58]. However, there was a clear trend in the outflow data. On the one hand, turbulent dissipation did not effect the outflow rate at early time, when fast material accelerated when the remnant rebounded reached the detection sphere. On the other hand, the ${\ell }_{\mathrm{mix}}>0$ runs had significantly larger outflow rates at later times, when the slower part of the dynamical ejecta crossed the detection sphere. This is visible for the ${\ell }_{\mathrm{mix}}={\ell }_{\mathrm{mix}}\left(\rho \right)$ simulations in Figure 5. Simulations with other values of ${\ell }_{\mathrm{mix}}$ followed the same trend, but are not included to avoid cluttering the figure. This late time boost of the outflow rate was 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 was, however, subdominant compared to the influence of neutrino reabsorption, which was found to play a key role in driving the “second wave” of the dynamical ejecta. Indeed, as shown in Figure 5, this part of the outflow was almost completely absent when neutrino reabsorption was switched off in the simulations. This second component of the dynamical ejecta was predominantly constituted of neutrino-irradiated material at high latitudes above the merger remnant.

Both viscosity and neutrino reabsorption played an important role in the determination of the electron fraction $Y_e$ of the ejecta [48,49,50,55,58], as shown in Figure 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 [124]. Turbulent dissipation increased the temperatures in the outflows, activating pair capture processes, which drove $Y_e$ to higher values through the reaction $e^{+}+n\to p+{\overline{\nu }}_e$. This effect was particularly pronounced when considering the tidal tail, which corresponded to the low-$Y_e$ peak in the outflow distribution. Neutrino absorption had an even larger effect via the reaction ${\nu }_e+n\to p+e^{-}$. Neutrinos generated relatively high $Y_e$ ejecta, especially at high latitudes [50]. They also affected the tidal tail, which was irradiated after colliding with the faster shock-accelerated outflows that were 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 [44,47,56,57,61,62,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 [61]. However, due to the large computational costs, we could not systematically vary ${\ell }_{\mathrm{mix}}$ in [61]. Unfortunately, the simulations we presented here did not span a sufficiently long time to be able to study the secular ejecta, so we leave this investigation to future work.

4. 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 their 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 [81], we proposed a scheme to include subgrid-scale turbulence effects into global simulations and studied the associated uncertainties on the multi-messenger signatures of BNS mergers. Our method extended the LES methodology to GR. We derived evolution equations for coarse-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 [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 was 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 multi-messenger observations of BNS mergers. The method we propose here could be and has been extended to the full-GRMHD equations [84,85]. Including GRMHD in our simulations will be 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.