General-Relativistic Large-Eddy Simulations of Binary Neutron Star Mergers

TL;DR

This paper introduces a new general-relativistic large-eddy simulation method to model turbulence in binary neutron star mergers, revealing its impact on remnant structure, emissions, and gravitational-wave signals, with implications for astrophysical observations.

Contribution

It develops the first large-eddy simulation formalism for turbulence in general relativity applied to neutron star mergers, enhancing understanding of turbulence effects on merger outcomes.

Findings

Turbulence influences remnant structure and lifetime.

Gravitational-wave signals are only weakly affected by turbulence.

Post-merger GW models are validated by the simulations.

Abstract

The flow inside remnants of binary neutron star (NS) mergers is expected to be turbulent, because of magnetohydrodynamics instability activated at scales too small to be resolved in simulations. To study the large-scale impact of these instabilities, we develop a new formalism, based on the large-eddy simulation technique, for the modeling of subgrid-scale turbulent transport in general relativity. We apply it, for the first time, to the simulation of the late-inspiral and merger of two NSs. We find that turbulence can significantly affect the structure and survival time of the merger remnant, as well as its gravitational-wave (GW) and neutrino emissions. The former will be relevant for GW observation of merging neutron stars. The latter will affect the composition of the outflow driven by the merger and might influence its nucleosynthetic yields. The accretion rate after black-hole formation is also affected. Nevertheless, we find that, for the most likely values of the turbulence mixing efficiency, these effects are relatively small and the GW signal will be affected only weakly by the turbulence. Thus, our simulations provide a first validation of all existing post-merger GW models.

Figures

Equations

$$?[c]T_{\mu }{}_{\nu }?=\rho h{u}_{\mu }{u}_{\nu }+p?[c]g_{\mu }{}_{\nu }?,$$

$$?[c]T_{\mu }{}_{\nu }?=E{n}_{\mu }{n}_{\nu }+S_{\mu }{n}_{\nu }+S_{\nu }{n}_{\mu }+S_{\mu \nu },$$

$$E=T_{\mu \nu }{n}^{\mu }{n}^{\nu }=\rho h{W}^2-p,$$

$$S_{\mu }=-{\gamma }_{\mu \alpha }{n}_{\beta }{T}^{\alpha \beta }=\rho h{W}^2{v}_{\mu },$$

$$S_{\mu \nu }={\gamma }_{\mu \alpha }{\gamma }_{\mu \beta }{T}^{\alpha \beta }=S_{\mu }{v}_{\nu }+p{\gamma }_{\mu \nu },$$

$$\displaystyle\begin{split}\partial_{t}\big{(}\sqrt{\gamma}&S_{i}\big{)}+\partial_{j}\Big{[}\alpha\sqrt{\gamma}\big{(}?[c]S_{i}^{j}?+S_{i}n^{j}\big{)}\Big{]}=\\ &\alpha\sqrt{\gamma}\Big{(}\frac{1}{2}S^{jk}\partial_{i}\gamma_{jk}+\frac{1}{\alpha}S_{k}\partial_{i}\beta^{k}-E\partial_{i}\log\alpha\Big{)}\,,\end{split}$$

$$\displaystyle\begin{split}\partial_{t}\big{(}\sqrt{\gamma}E\big{)}+&\partial_{j}\Big{[}\alpha\sqrt{\gamma}\big{(}S^{j}+En^{j}\big{)}\Big{]}=\\ &\alpha\sqrt{\gamma}\Big{(}K_{ij}S^{ij}-S^{i}\partial_{i}\log\alpha\Big{)}\,,\end{split}$$

$$\displaystyle\begin{split}\partial_{t}\big{(}\sqrt{\gamma}&\overline{S_{i}}\big{)}+\partial_{j}\Big{[}\alpha\sqrt{\gamma}\big{(}\overline{?[c]S_{i}^{j}?}+\overline{S_{i}}n^{j}\big{)}\Big{]}=\\ &\alpha\sqrt{\gamma}\Big{(}\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\Big{)}\,,\end{split}$$

$$\displaystyle\begin{split}\partial_{t}\big{(}\sqrt{\gamma}\overline{E}\big{)}+&\partial_{j}\Big{[}\alpha\sqrt{\gamma}\big{(}\overline{S^{j}}+\overline{E}n^{j}\big{)}\Big{]}=\\ &\alpha\sqrt{\gamma}\Big{(}K_{ij}\overline{S^{ij}}-\overline{S^{i}}\partial_{i}\log\alpha\Big{)}\,.\end{split}$$

$$\overline{S_i{v}_j}=\overline{S_i}\overline{v_j}+{\tau }_{ij}.$$

$$\tau_{ij}=-2\nu_{T}\rho hW^{2}\left[\frac{1}{2}\big{(}\nabla_{i}\overline{v_{j}}+\nabla_{j}\overline{v_{i}}\big{)}-\frac{1}{3}\nabla_{k}\overline{v^{k}}\gamma_{ij}\right]\,,$$

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

$${\lambda }_{\mathrm{MRI}}\sim 3\mathrm{m}{\left(\frac{\Omega }{4\mathrm{rad}{\mathrm{ms}}^{-1}}\right)}^{-1}\left(\frac{B}{10^{14}\mathrm{G}}\right).$$

$$\Omega =\frac{{\Omega }_c}{1+{\left(\frac{\varpi }{R_e}\right)}^2},$$

Full text

\tensordelimiter

? \tensorformatlrc

1. Introduction

The typical outcome of the merger of two neutron stars is expected to be the formation of an hypermassive neutron star (HMNS): a massive NS temporarily supported against gravitational collapse by its fast differential rotation, although prompt black hole (BH) formation might occur for large masses and/or soft equations of state (Baiotti & Rezzolla, 2016, and references therein). Its survival time and, in general, its properties, are important for the multimessenger signature of NS mergers and for their nucleosynthetic yields. The HMNS has a magnetar-level B-field, and it is a bright source of neutrinos (Sekiguchi et al., 2011; Kiuchi et al., 2014). These could drive baryon-rich winds (Dessart et al., 2009; Siegel et al., 2014). The presence of an HMNS could significantly boost the neutrino annihilation rates at high-latitudes (Richers et al., 2015; Perego et al., 2017) and perhaps contribute to the launching of a relativistic jet and a short $\gamma$-ray burst (SGRB) (Nakar, 2007). Neutrinos could also affect the yield and electromagnetic signature of the r-process nucleosynthesis in the binary ejecta (Wanajo et al., 2014; Metzger & Fernández, 2014). Long-lived massive NSs created in mergers might power the X-ray tails observed in some SGRB (Rowlinson et al., 2013; Lasky et al., 2014; Gao et al., 2016). Finally, gravitational waves from the HMNS could be used to constrain its EOS (Bauswein & Janka, 2012; Takami et al., 2014; Bernuzzi et al., 2015a; Radice et al., 2016a).

Despite the rapid recent progress of general-relativistic magnetohydrodynamics (GRMHD) simulations (Rezzolla et al., 2011; Kiuchi et al., 2014; Ruiz et al., 2016), the impact of magnetoturbulence on the structure and survival time of the HMNS is highly uncertain. The magnetorotational instability (MRI) (Balbus & Hawley, 1991) is expected to operate inside the HMNS, drive the redistribution of angular momentum and affect its lifetime and properties (Duez et al., 2006; Siegel et al., 2013). Unfortunately, the fastest growing mode of the MRI in these systems is inaccessible even to the highest-resolution simulations (Kiuchi et al., 2015).

A possible way to model the impact of turbulent transport of angular momentum in the HMNS would be to use an effective viscosity (Duez et al., 2004). This approach is made difficult by the fact that the Navier-Stokes equations describing relativistic viscous flows are known to exhibit a number of unphysical pathologies (Hiscock & Lindblom, 1985; Kostadt & Liu, 2000). There are more complex fluid models that do not have these shortcomings (Andersson & Comer, 2007; Rezzolla & Zanotti, 2013). However, they are also not entirely without problems (Majorana & Motta, 1985; Hiscock & Lindblom, 1988), are difficult to implement (e.g. Takamoto & Inutsuka, 2011), and their non-linear properties are poorly understood. More importantly, they contain a large number of transport coefficients that have no classical counterpart. These have no clear physical meaning and are not even in principle measurable (Geroch, 1995; Lindblom, 1996).

Here, we propose an alternative approach. Our starting point is the observation that turbulence models do not have to be restricted to the class of equations describing fluids with physical viscosity or heat transfer. Instead, we develop an effective model based on a general-relativistic (GR) extension of the Newtonian large-eddy simulation (LES) framework (e.g. Miesch et al., 2015). Our model, while recovering the Navier-Stokes equations in the Newtonian limit, does not correspond to or have the same limitations as any relativistic theory of viscous flows.

In this Letter, after a brief description of the general-relativistic large-eddy simulation (GRLES) formulation, we present, for the first time, simulations in full-GR of merging NS with a realistic, tabulated, nuclear EOS, neutrino cooling, and parametrized turbulent transport. We show that turbulence could influence the structure of the HMNS, as well as its GW and neutrino emissions. On the other hand, for the most realistic values of the turbulent viscosity, these effects appear to be small and our simulations provide an important confirmation of a number of previous results where turbulent transport was not included.

2. Formulation

Our starting point is the stress energy tensor of a perfect fluid

[TABLE]

where $\rho$, $h$, $u_{\mu}$ and $g_{\mu\nu}$ are, respectively, density, specific enthalpy, four-velocity, and the spacetime metric.

In numerical relativity, spacetime is decomposed in space-like slices with normal $n^{\mu}$. We decompose $?[c]T_{\mu}{}_{\nu}?$ with respect to $n^{\mu}$ as

[TABLE]

where

[TABLE]

and $\gamma_{\mu\nu}$, $v^{\mu}$, $p$, and $W$ are, respectively, the spatial metric, the three-velocity, the pressure, and the Lorentz factor.

The equations of energy and momentum conservation are

[TABLE]

where $\alpha$, $\beta^{i}$, $K_{ij}$ are, respectively, the lapse function, shift vector, three-metric, and extrinsic curvature. $\sqrt{\gamma}$ is the spatial volume element. These equations are then closed with an EOS and equations describing the conservation of the baryon and lepton numbers.

Equations (6) and (7) contain modes at all scales, but, in numerical simulations, only modes resolved with sufficiently many grid zones can develop. In essence, any simulation deals only with a coarse-grained version of the hydrodynamics equations. In the LES framework, this observation is made rigorous with the introduction of a linear filtering operation $u\mapsto\overline{u}$ that removes features at scales smaller than a given $\Delta$. Here, we adopt for the filtering operator the cell-averaging of the finite-volume discretization of the equations. We leave the investigation of more advanced filters for future work. If we filter Eqs. (6) and (7) we obtain

[TABLE]

Note that Eqs. (8) and (9) are exact, but are not closed. The reason is that $\overline{S_{i}v_{j}}$ cannot be expressed only in terms of other coarse-grained quantities. A closure is needed:

[TABLE]

$\tau_{ij}$ is the so-called subgrid-scale turbulent tensor. Similar terms appear in the coarse graining of the baryon and lepton number conservation equations, but, for simplicity, we will neglect them here. Simulations usually assume $\tau_{ij}=0$. Here, instead, we will use $\tau_{ij}$ to model small-scale turbulence in merger simulations. To do so, in analogy with the classical Newtonian closure of Smagorinsky (1963), we choose the ansatz

[TABLE]

where $\nabla$ is the covariant derivative compatible with $\gamma_{ij}$. The quantity $\nu_{T}$ has a dimension of a viscosity. On dimensional grounds, we are led to assume

[TABLE]

where $c_{s}$ is the local sound speed, and $\ell_{\rm mix}$, often called the mixing length, is a characteristic length over which turbulence operates. Note that $\nu_{T}$ is not a physical viscosity; indeed, its definition depends on the numerical grid and on the Eulerian observer $n^{\mu}$. This is expected, because the notion of resolved and unresolved scales is observer dependent in relativity. $\nu_{T}$ should be calibrated on the basis of highly-resolved simulations and/or using self-similarity methods (e.g. Germano et al., 1991). We leave this task for future work. For now, we will treat $\ell_{\rm mix}$ as a free parameter. Assuming MRI turbulence, it is then natural to set $\ell_{\rm mix}\sim\lambda_{\rm MRI}$, where (Duez et al., 2006)

[TABLE]

Equations (8), (9),  (11), and (12), together with the EOS, and the continuity equations are what we refer to as the GRLES equations. We verified, by repeating the analysis of Hiscock & Lindblom (1985) and numerically, that the GRLES equations are not affected by the same pathologies as the relativistic Navier-Stokes equations.

3. Implementation

We implement the GRLES equations into the general-relativistic hydrodynamics (GRHD) code WhiskyTHC (Radice & Rezzolla, 2012; Radice et al., 2014a, b). With our current choice of the filtering operator, this only amounts to the inclusion of $\tau_{ij}$ in the equations. We treat the viscous fluxes in a flux-conservative way and we self-consistently include the turbulent stress-tensor in the energy and momentum source terms, as well as in the calculation of the spacetime geometry.

For the simulations presented here, we use the microphysical EOS of Lattimer & Swesty (1991) with nuclear compressibility parameter $K=220\ {\rm MeV}$. Neutrino cooling is treated with the scheme presented in (Radice et al., 2016b). For $\ell_{\rm mix}$, we consider the values 0 (our reference run), 5, 25, and 50 meters. Over this range, 5 meters is the most likely value for $\ell_{\rm mix}$ given Eq. (13), while 50 meters might be unphysically large, in the light of the lack of convergence observed in the 17-meter resolution simulation of Kiuchi et al. (2015).

As a first example, we consider the evolution of an equilibrium configuration constructed with the RNS code (Stergioulas & Friedman, 1995). The initial configuration has gravitational mass $M\simeq 2.45\ M_{\odot}$ and angular momentum $J/M^{2}\simeq 0.66\ G/c$. We use the differential rotation law of Komatsu et al. (1989), which, in the Newtonian limits reduces to

[TABLE]

where $\varpi$ is the cylindrical radius, $\Omega_{c}$ is the angular velocity at the center, and $R_{e}$ is the stellar equatorial radius. The resolution for this test is $\simeq 370\ {\rm m}$.

We plot the maximum density as a function of time in Fig. 1. As expected on the basis of previous work (Duez et al., 2004), the inclusion of turbulent viscosity results in the transport of angular momentum leading to gravitational collapse. This test shows that WhiskyTHC is able to capture the effect of turbulent viscosity even at low resolution.

4. Binary Neutron Star Mergers

We consider the last ${\sim}4$ orbits and merger of two $1.35$-$M_{\odot}$ NSs. We already evolved this binary in Bernuzzi et al. (2015b), where a description of the properties of the initial data is also given. For the evolution, we use the high-resolution setup of Bernuzzi et al. (2015b), with the improvements discussed in Radice et al. (2016a). We perform simulations with resolutions, on the finest refinement level of ${\sim}185\ {\rm m}$ and ${\sim}246\ {\rm m}$. We present results from the high-resolution simulations. In the low-resolution simulations, there are quantitative, but not qualitative differences.

We find that turbulent viscosity has a much less obvious impact on the evolution of the HMNS than what could have been anticipated on the basis of the idealized model in Sec. 3. In the first few milliseconds after merger, turbulent transport results in a decrease of the compactness, as can be seen from the maximum density evolution (Fig. 2; left panel). Over longer timescales, the behavior is non-linear. The $\ell_{\rm mix}$-25-m HMNS is the most compact remnant and collapses to a BH ${\sim}17\ {\rm ms}$ after merger. The $\ell_{\rm mix}$-5-m remnant is only slightly less compact than that of the reference simulation $\ell_{\rm mix}$-0-m. BH formation occurs at ${\sim}20\ {\rm ms}$ and ${\sim}22\ {\rm ms}$ after merger for the $\ell_{\rm mix}$-5-m and $\ell_{\rm mix}$-0-m binaries, respectively. The $\ell_{\rm mix}$-50-m HMNS is the least compact and does not collapse to a BH within our simulation time. For the models that collapse within our simulation time, we observe the formation of a massive (${\sim}0.1\ M_{\odot}$) accretion disk (Fig. 2; right panel). As could have been anticipated, the accretion rate is larger for simulations with larger $\ell_{\rm mix}$.

The reason for the different evolutions of the remnant can be understood from the analysis of its internal structure (Fig. 3). The rotational profile established in the HMNS after the initial, very dynamical, phase is qualitatively different from that of Eq. (14), as also pointed out by Shibata et al. (2005); Kastaun et al. (2016); Hanauske et al. (2016); Ciolfi et al. (2017). Consistently with these previous studies, we find in the $\ell_{\rm mix}$-0-m simulation an HMNS composed of a slowly rotating core and a rotationally supported massive envelope. As the mixing length increases, the structure of the HMNS is altered due to interplay between three competing effects. First, angular momentum redistribution spins up the core, reducing its compactness. Second, the loss of angular momentum from the massive envelope results in a compression the HMNS. Third, as more kinetic energy is converted into thermal energy by turbulent dissipation, the inner core becomes hotter and expands because of the increased pressure. The interplay between these effects is complicated by the fact that the angular momentum of the HMNS is not conserved, but is radiated in GWs at a rate proportional to that of the gravitational binding energy (Bernuzzi et al., 2015b). For this reason, as the HMNS contracts, it becomes more bound and at the same time it looses angular momentum support.

The first and third effect are dominant at early times, so the effect of turbulent viscosity is to monotonically reduce the compactness in the first few milliseconds after merger. Later, all three effects become important. At this stage, energy and angular momentum losses to GW play an important role. In the case of the $\ell_{\rm mix}$-5-m binary, the envelope remains centrifugally supported (Fig. 3; upper-left panel), so the compactness is slightly decreased compared to the reference run without turbulence dissipation. For the $\ell_{\rm mix}$-25-m binary, the effect of turbulent transport is qualitatively similar to the $\ell_{\rm mix}$-5-m binary at early times. Later, its envelope contracts causing the growth of the central density (Fig. 2; left panel) and early BH formation. Finally, in the case of the $\ell_{\rm mix}$-50-m run the thermal effect prevails; the hot spots formed in the contact layer at the time of merger sink to the center and enhance the core temperature to ${\sim}70\ {\rm MeV}$. The increased thermal support in the layers with ${\sim}5\cdot 10^{14}\cdot{\rm g}\cdot{\rm cm}^{-3}$ inflates the HMNS. The reduced compactness, in turn, results in a decrease of the angular momentum loss due to GW and prevents its collapse within the simulation time.

The changes in the HMNS structure are reflected in its multimessenger emissions. The total energy radiated in GWs (Fig. 4; left panel) is closely related to the rate of increase of the HMNS compactness. For this reason, at early times, the amplitude of the signal slightly decreases with $\ell_{\rm mix}$, while, over longer timescales, the behavior is non-monotonic. The characteristic GW frequency after merger (Fig. 4; right panel) is, instead, only weakly affected, with the exception of a slight growth before BH formation, which is a commonly observed feature (e.g. Radice et al., 2016a).

The neutrino emission (Fig. 5) is also influenced by the turbulent dissipation and the consequently higher temperatures in the HMNS. The luminosity of neutrinos of all flavours increases with the mixing length parameter up to $\ell_{\rm mix}=25\ {\rm m}$. The luminosity of the $\ell_{\rm mix}$-50-m simulation is, however, smaller than that of the $\ell_{\rm mix}$-25-m simulation. This is possibly because, in the $\ell_{\rm mix}$-50-m HMNS, the maximum of the temperature occurs at the center, while for the other models it is off-centered (see also Kastaun et al., 2016).

5. Discussion and Conclusions

We have developed a new framework for the modeling of turbulence in full-GR simulations. Our approach is based on a relativistic extension of the large-eddy simulation technique, which represents the state-of-the-art for turbulence modeling in classical hydrodynamics (Miesch et al., 2015). Our method can naturally exploit turbulent closures developed in Newtonian physics, is simple to implement, robust, and stable.

As a first application, we have employed a turbulent viscosity closure to study the effect of angular momentum transport and dissipation in NS mergers. We have performed, for the first time, general-relativistic large-eddy simulations of merging NS with microphysical nuclear EOS and neutrino cooling. We have found that turbulence can modify the structure and collapse time of the merger remnant. These, in turn, are reflected in the GW and neutrino emissions from the HMNS. The accretion rate after BH formation is also affected.

The total energy radiated in GW is the most affected quantity, since it closely tracks the contraction of the HMNS on its way to the final collapse to BH. In the presence of very efficient turbulent transport, the effective viscosity might mask changes in the compactness of the HMNS that would otherwise be attributable to changes in the high-density component of the EOS (Radice et al., 2016a). This effect is, however, only modest for more conservative choices of the turbulent mixing-length parameter. In these cases, turbulence would not significantly affect the prospect of detecting phase transitions in the core of the HMNS using GW observations. However, a definitive statement will have to wait until sufficiently resolved GRMHD simulations are available to estimate $\ell_{\rm mix}$.

We have also found that the post-merger GW frequency is only weakly affected by the effective turbulent viscosity. Thus, our results provide an important validation of the several proposed methods relying on its measure to constrain the EOS of dense nuclear matter (Bauswein & Janka, 2012; Takami et al., 2014; Bernuzzi et al., 2015a). Our results also reaffirm the observation by Bernuzzi et al. (2015a) that the post-merger GW peak-frequency is set at the time of merger. Afterwards, the frequency stays close to constant and is largely insensitive to the evolution of the HMNS, with the exception of the signature of BH formation.

Finally, our results show that the neutrino signal is also influenced by the turbulent dissipation of kinetic energy into heat. The increased temperatures and luminosities, especially for the anti-electron neutrinos, will influence the proton fraction in the outflows and might have an effect on the resulting nucleosynthetic yields (Wanajo et al., 2014; Metzger & Fernández, 2014; Foucart et al., 2016). Our results strongly suggest that turbulent dissipation will have to be included in the next generation of neutrino-radiation-hydrodynamics models of the outflows from merging NSs.

Here, we presented a first application of the GRLES method. In the future, on the one hand, we will extend the present study to more binary configurations and EOS. On the other hand, work is already underway to develop closures tuned with highly-resolved GRMHD simulations of HMNSs. Finally, we will extend GRLES to GRMHD and couple it with a subgrid-scale dynamo model such as the one of Giacomazzo et al. (2015).

It is a pleasure to thank S. Bernuzzi for the many stimulating discussions on binary neutron star mergers. I also thank S. Hild for the ET-D noise curve data and A. Burrows, P. Mösta, L. Rezzolla, and C. D. Ott for discussions. I gratefully acknowledge support from the Schmidt Fellowship and the Sherman Fairchild Foundation. The simulations were performed on Stampede NSF XSEDE (TG-PHY160025), and employed computational resources provided by the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology.