Numerical relativity simulations of neutron star merger
using conservative mesh refinement
We study equal and unequal-mass neutron star mergers by means of new numerical relativity simulations in which the general relativistic hydrodynamics solver employs an algorithm that guarantees mass conservation across the refinement levels of the computational mesh. We consider eight binary configurations with total mass , mass-ratios and , and four different equations of state (EOSs), and one configuration with a stiff EOS, and , which is the largest mass ratio simulated in numerical relativity to date. We focus on the post-merger dynamics and study the merger remnant, dynamical ejecta and the postmerger gravitational wave spectrum. Although most of the merger remnant are a hypermassive neutron star collapsing to a black hole+disk system on dynamical timescales, stiff EOSs can eventually produce a stable massive neutron star. During the merger process and on very short timescales, about of material become unbound with kinetic energies . Ejecta are mostly emitted around the orbital plane; and favored by large mass ratios and softer EOS. The postmerger wave spectrum is mainly characterized by the non-axisymmetric oscillations of the remnant neutron star. The stiff EOS configuration consisting of a and a neutron star, simulated here for the first time, shows a rather peculiar dynamics. During merger the companion star is very deformed; about of rest-mass becomes unbound from the tidal tail due to the torque generated by the two-core inner structure. The merger remnant is a stable neutron star surrounded by a massive accretion disk of rest-mass . This and similar configurations might be particularly interesting for electromagnetic counterparts. Comparing results obtained with and without the conservative mesh refinement algorithm, we find that post-merger simulations can be affected by systematic errors if mass conservation is not enforced in the mesh refinement strategy. However, mass conservation also depends on grid details and on the artificial atmosphere setup; the latter are particularly significant in the computation of the dynamical ejecta.
pacs:04.25.D-, 04.30.Db, 95.30.Sf, 95.30.Lz, 97.60.Jd 98.62.Mw
Binary neutron star (BNS) mergers are extreme events associated to a variety of observable phenomena in the gravitational and electromagnetic spectra, e.g. Eichler et al. (1989); Andersson et al. (2013); Rosswog (2015). BNS coalescence is primarily driven by the emission of gravitational waves (GWs). Indirect evidence for GWs has been indeed inferred by radio observation of double pulsars Hulse and Taylor (1975); Weisberg et al. (2010); Burgay et al. (2003); Lyne et al. (2004); Kramer et al. (2006), but a direct detection of GWs is still pending. The GW signal emitted during the last minutes of the coalescence and merger is in the band of ground-based laser interferometer network made of LIGO LIG () and Virgo Vir (). Within the next years, this network will start to operate at sensitivities where detections per year are expected Abadie et al. (2010); Aasi et al. (2013). Several electromagnetic counterparts are expected both during and following BNS mergers; joint observations of the gravitational and electromagnetic emissions will maximize the scientific returns Metzger and Berger (2012). Neutron star mergers are usually associated to short-gamma ray burst (and afterglows) Paczynski (1986); Eichler et al. (1989). Although the precise injection mechanism has not been clearly identified, BNSs remain the most plausible triggers of these powerful emissions. Dynamical ejecta from BNS are currently the most plausible site of origin of heavy nuclei () rapid neutron-capture process Lattimer and Schramm (1974); Rosswog et al. (1999); Goriely et al. (2011). The radioactive decay of some of these newly produced heavy elements is likely to lead to strong electromagnetic transients called kilonova (or macronova) events Li and Paczynski (1998); Tanvir et al. (2013); Metzger et al. (2010). Finally, a large amount of energy is released in neutrinos, produced by the merger remnant either via shocks Waxman (2004); Dermer and Holmes (2005) and neutron-rich outflows Bahcall and Meszaros (2000), or, at lower energies, in the hot dense regions of the hypermassive neutron star (HMNS) Dessart et al. (2008); Perego et al. (2014). However, the steep energy dependence of neutrinos of the interaction cross sections and their moderate energies ( MeV) make them hard to detect.
Modeling BNS mergers requires relativistic hydrodynamics simulations in dynamical spacetimes, i.e. the solution of the full set of Einstein’s field equations. General relativistic BNS simulations are typically performed in the framework of 3+1 numerical relativity using Cartesian-grids, finite volume methods, and explicit time evolutions, see Faber and Rasio (2012) for a review. A crucial ingredient in such numerical setups is the use of adaptive mesh refinement (AMR), in particular the methods of Berger and Oliger (1984), which were implemented for various applications in numerical relativity Choptuik (1989); Brügmann (1996); Schnetter et al. (2004); Evans et al. (2005); Brügmann et al. (2008). Nested Cartesian boxes with 2:1 grid spacing refinement and Berger-Oliger time stepping Berger and Oliger (1984) proved to be a robust and stable solution for the computation of black hole Brügmann (1999); Baker et al. (2006); Campanelli et al. (2006); Brügmann et al. (2008) and neutron star mergers Shibata and Uryu (2000); Shibata et al. (2005); Baiotti et al. (2008); Thierfelder et al. (2011a) as well as rotational collapse of neutron stars Baiotti et al. (2005); Reisswig et al. (2013a); Dietrich and Bernuzzi (2015) or massive stars Ott et al. (2004); Shibata et al. (2006).
One of the main problems in the simulation of hydrodynamical flows with finite volume methods and AMR techniques is to preserve global conservation of mass and other quantities, especially in the presence of shocks, contact discontinuities and large gradients. In a seminal paper, Berger and Colella have proposed a refluxing scheme which guarantees conservation across refinement level Berger and Colella (1989). Essentially, the refluxing scheme enforces the fluxes in and out across a coarse/fine cell boundary to be the same to round-off level. Algorithmically it consists in a correction step applied to the solution at certain grid points after each time step.
The importance of a conservative mesh refinement is illustrated in Fig. 1, for the simplest case of the 1D advection equation, with and discontinuous initial data, for , otherwise. We employ a grid composed of 3 fixed levels with grid points, centered around . The coarse level has grid spacing and the others are successively refined by factors two, . Time evolution is performed with a 4th order Runge-Kutta and the Berger-Oliger method; fluxes are computed with a linear reconstruction using the Van Leer MC2 limiter. Figure 1 shows the evolution of the mass of the system on each refinement level, , which without mesh refinement (unigrid) is conserved to round off precision. Using mesh refinement, one observes that every time the mass flows in/out a refinement level (see e.g. , top panel), a mass violation takes place (bottom panel). Notably, the mass either decreases or increases in a way that depends on the scheme’s truncation error and that in general is not predictable. Also, after the wave has left all inner refinement levels an error in the mass on is still present. Instead, using the Berger and Colella correction step, mass conservation is verified at the round-off error, exactly as in the unigrid case.
Conservative AMR schemes have been introduced in numerical relativity only very recently East et al. (2012); Reisswig et al. (2013a). They have been used to simulate eccentric mergers of both black holes - neutron star and double neutron star systems, including the merger remnant, post-merger disks and ejecta Stephens et al. (2011); East and Pretorius (2012). Also, they have been employed in massive star and core-collapse supernovae evolutions in general relativity Ott et al. (2013); Reisswig et al. (2013b); Abdikamalov et al. (2014). Recent studies of rotating neutron star collapse to black hole greatly benefit of the use of conservative AMR, and allowed an accurate calculation of the gravitational wave signal Reisswig et al. (2013a); Dietrich and Bernuzzi (2015) and a local comparison of the end state with black hole spacetimes Dietrich and Brügmann (2014).
In binary simulations one expects that conservative AMR can significantly improve numerical relativity simulations, especially simulations of the merger remnant. A direct comparison of the performance of a conservative mesh refinement algorithm in coalescing BNS systems is presently missing. In the context of spinning equal-mass quasi circular mergers, we have pointed out that the simulation of the hypermassive neutron star is sensitive to the mesh boxes size and their extension Bernuzzi et al. (2014a). The latter factors influence mass conservation (for a fixed resolution), and a conservative scheme is desirable. Another potentially important application of conservative AMR is the simulations of low-density material in postmerger accretion disks and dynamical ejecta. Ejecta have densities several orders of magnitude smaller than the typical neutron star maximum densities; thus, their calculation employing grid-based codes is very challenging. Dynamical ejecta in full general relativistic BNS merger simulations have been previously studied only in Hotokezaka et al. (2013a); Sekiguchi et al. (2015) in more detail. Those works do not employ a conservative AMR strategy, thus the accuracy of the result can be, in principle, seriously compromised.
The purpose of this paper is threefold.
First, we present our implementation of a conservative AMR algorithm and present a set of single star spacetime evolutions in which we assess the performances of the algorithm. We focus on the evolution of different single star spacetimes since such tests (i) received little attention in the literature; (ii) are computationally relatively cheap; (iii) are highly nontrivial and preparatory cases for the application of the code to BNs evolutions.
Second, we apply our upgraded code to the study of equal-mass (mass ratio ) and unequal-mass () BNS system described by various equations of state (EOS). We directly compare results obtained with and without the conservative AMR. We focus on the postmerger dynamics and investigate the physical properties of the remnant. In particular we study as a function of the EOS and the mass ratio the following properties: (i) the merger outcome; (ii) mass and kinetic energy of the dynamical ejecta; (iii) GW spectra.
Third, we consider for the first time the evolution of a BNS system with a stiff EOS and mass ratio (total mass ). This binary has the largest mass ratio simulated so far (see also Shibata and Taniguchi (2006)). The particular combination of EOS, , and total mass properties lead to a peculiar merger remnant composed of a stable massive neutron star surrounded by an extended, massive accretion disk. Also, the binary configuration favors mass ejection during merger. These kind of binary configurations are possible and might be particular relevant for electromagnetic counterparts. However, they have received little attention in numerical relativity, although some recent observations are in favor for a stiff EOS Hambaryan et al. (2011, 2014).
The article is structured as follows. After a brief review of the equations (Sec. II), we present our numerical strategy in Sec. III focusing on the novel implementation of the conservative mesh refinement. Section IV describes the main quantities employed for the analysis of our BNS simulations. In Sec. V we describe a variety of single star tests in which the performance of the conservative AMR is investigated for different combinations of the relevant parameters of the simulations (restriction and prolongation operators and artificial atmosphere parameters.) Section VI summarizes the BNS configurations and the grid setup used for evolutions. In Sec. VII we apply the new algorithm and evolve 16 BNS systems with mass ratios and and different EOS. In Sec. VIII we consider a BNS with , total mass , and the stiff equation of state MS1b. Finally, the conclusions are presented in Sec. IX. Throughout this article, geometrical units are employed unless otherwise stated. At some places units of are given explicitly for clarity.
Ii Summary of the Equations
Let us summarize briefly the most important equations employed in this work. We work with the 3+1 formalism (e.g. Gourgoulhon (2007)) and indicate with the 3-metric, and with and the lapse and shift vector.
General-relativistic hydrodynamics (GRHD) equations are solved in conservative form,
with being the vector of the Eulerian conservative variables defined in terms of the primitive variables as,
The primitive variables are the rest-mass density , the pressure , the specific internal energy , and the 3-velocity . Additionally, we define the Lorentz factor , the enthalpy , and the determinant of the 3-metric . On the right-hand-side of Eq. (1) one has the divergence of the fluxes and source terms depending on the metric, metrics first derivatives and fluid variables. We stress that only the first equation of (1) is a “strict” conservation law,
in the sense that the source term is zero and a conserved quantity can be associated: the rest-mass . We refer to Font (2007); Rezzolla and Zanotti (2013) for an extensive discussion of these equations.
The PDE system in (1) is closed by an equation of state (EOS) in the form . A simple EOS is the -law , or its barotropic version (polytropic EOS). Several barotropic – zero-temperature EOS developed to describe neutron star matter can be fit with piecewise polytropic models, and efficiently used in simulations. In our work we employ four segment fitting models following the construction of Read et al. (2009). Each segment is given by a certain rest-mass density interval ; the pressure is then calculated as where the polytropic constants are determined by demanding continuity of at the interfaces, . The parameters of our EOS are reported in Tab. 1; notice that we specify in our units. Thermal effects are simulated with an additive thermal contribution in the pressure in a -law form, , with , see Shibata et al. (2005); Thierfelder et al. (2011a); Bauswein et al. (2010).
The Einstein equations are written in 3+1 form, either as the BSSN Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1999) or the Z4c Bernuzzi and Hilditch (2010); Hilditch et al. (2013a) system. In the gauge sector, we use the 1+log-slicing condition Bona et al. (1995) for the lapse and the Gamma driver shift Alcubierre et al. (2003); van Meter et al. (2006). The fundamental role of this gauge in the numerical simulation of the gravitational collapse and singularity formation/evolution was investigated in different physical scenarios Thierfelder et al. (2011b); Staley et al. (2012); Hilditch et al. (2013b); Dietrich and Bernuzzi (2015).
Iii Numerical Method
In this work we use the numerical relativity methods implemented in the BAM code Thierfelder et al. (2011a); Brügmann et al. (2008, 2004); Brügmann (1999). Our new implementation of the conservative mesh refinement for hydrodynamics fields is based on the Berger-Colella method Berger and Colella (1989) and follows East et al. (2012); we describe it in detail in this section.
iii.1 Computational Grid
The computational grid is made of a hierarchy of cell-centered nested Cartesian grids. The hierarchy consists of levels of refinement labeled by . A refinement level has one or more Cartesian grids with constant grid spacing and points per direction. The grid spacing in each refinement level is refined according to . The grids are properly nested in such a way that the coordinate extent of any grid at level , , is completely covered by the grids at level . Some of the mesh refinement levels can be dynamically moved and adapted during the time evolution according to the technique of “moving boxes”, e.g. Yamamoto et al. (2008); Brügmann et al. (2008); Baiotti et al. (2008). BAM’s grid can be further extended in the wave zone using a multipatch “cubed-sphere” as described in Ronchi et al. (1996); Thornburg (2004); Pollney et al. (2011); Hilditch et al. (2013a). Every refinement level has buffer zones populated by interpolation; interpolation from the parent (coarse) to the child (fine) level is the prolongation (P) operation, the one from the fine to the coarse level is the restriction (R) operation. For metric variables these operations are performed with sixth order Lagrangian operators. Spatial interpolation of matter variables is discussed below.
The grid variables are evolved in time with the method of lines, using an explicit fourth order Runge-Kutta and employing the Berger-Oliger (BO) algorithm Berger and Oliger (1984). For efficiency, we typically use only six buffer zones and perform a linear interpolation in time to update the buffer zones during the Runge-Kutta step, see Brügmann et al. (2008) for more details. A Courant-Friedrich-Lewy factor of is employed in all runs, if not stated differently. Standard finite differencing 4th order stencils are employed for the spatial derivatives of the metric. GRHD is solved by means of a high-resolution-shock-capturing method Thierfelder et al. (2011a) based on primitive reconstruction and the Local-Lax-Friedrich’s (LLF) central scheme for the numerical fluxes. Primitive reconstruction is performed with the 5th order WENO scheme of Borges et al. (2008) as in Bernuzzi et al. (2012).
iii.2 Conservative mesh refinement
Let us review the main idea of the new conservative AMR algorithm implemented in the BAM code. Without loss of generality we restrict the presentation to the first equation of (1), and to the flux in the -direction only. Although the algorithm is applied to all the fluid variables, Eq. (3), the -equation, is the only one which is a strict conservation law. Directions different from the direction are treated in a similar way.
The discrete model equation reads,
where denotes the -component of the numerical flux across the
cell face (boundary of cell and
), , denotes the time level, and
the time step. Consider the model equation on two
sequential levels of refinement with , and on cells at the
boundary of refinement . Mass violation happens during a BO step, because:
(i) the buffer zones of level are set by prolongation (P) from level
; (ii) the prolongation carries a certain
truncation error, so the fluxes on differ from those on ;
(iii) after restriction (R) from level , the solution on level is
not consistent with the fluxes on .
The process is illustrated in Fig. 2.
After the time step , the changes of the variable on level due to the flux going through the cell face is given by
After level , level advances by two time steps and one has
In general, these two changes are different at truncation error level. Similarly the mass flows across the face are different, , and, after restriction, the mass conservation is violated in a way .
The original Berger-Colella algorithm corrects the solution at level after the refinement level has completed its time step and both levels are time-aligned Berger and Colella (1989). The correction (C) operation is , where is a flux correction stored on the cell face. First, is initialized with before advancing in time level . Then, during each time step of level , it receives and sums up the contributions (two contributions in our example). The C step guarantees consistency of the fluxes. East et al. East et al. (2012) proposed to store the mass correction rather than , and perform the correction as where is the cell volume East et al. (2012). This method is simpler and has the advantage of using grid variables defined on cell centers instead of faces. We follow this approach.
Our implementation is as follows:
We introduce a mask to label the cells involved in the C step. These are the innermost buffer points of level (red in Fig. 2) and the corresponding parent cells (blue in Fig. 2). The mask also stores the information about the box face, i.e. one of the possibilities . The mask has to be recomputed after each regridding step.
After each evolution step we store the mass change of the parent cells
and, similarly, after each sub-step, for level . Notice that the particular sign depends on the entry in the mask, e.g. -surfaces refer to a positive sign in (7).
When the parent and the child level are aligned in time, we sum up the contributions and correct the cell values with,
We observe that the effectiveness of the algorithm depends crucially on the specific RP operators. For hydrodynamics fields the R step is conservative if the operation is performed using local averages, which are second order accurate, . Similarly, a safe choice for the P step is linear interpolation using limiters in order to control oscillations. However, for neutron star spacetime simulations high-order operators may be important for accuracy and faster convergence. As indicated in Tab. 2, we have implemented several RP operators, including ENO 2nd order Harten (1989), Lagrangian and WENO 4th order Thierfelder et al. (2011a). In the next sections we will present results for various combinations of RPC operators.
iii.3 Atmosphere treatment
For the simulation of neutron star spacetimes, the vacuum region outside the stars requires special treatment. As described in Thierfelder et al. (2011a), we use a low-density static and barotropic atmosphere at a density level
During the recovery of the primitive variables from the conservative variables, a point is set to atmosphere if the density is below the threshold
The atmosphere treatment violates mass conservation and can potentially affect, or invalidate, the improvements related to the conservative AMR. In the following we will investigate this aspect in some detail experimenting with parameters and (see in particular Sec. V.0.3).
Iv Simulation analysis
In this section we describe the quantities employed for the analysis of our simulations. Notably, we introduce some diagnostic which are helpful to investigate the energetics/geometry of the ejected material.
The performances of the conservative AMR scheme are mainly tested using the baryonic, or rest-mass, mass integral,
which should remain constant during the evolution, compare Eq. (3). The rest-mass, and the other integrals discussed in this work, are calculated on each refinement level. We usually report results for a given level, which is the appropriate one for the particular quantity; e.g. the baryonic mass is reported on the level.
The merger remnant of several BNS configurations considered here is a hypermassive neutron star (HMNS) which collapses to black hole on a dynamical timescale. The lifetime is typically calculated from the moment of merger (see below) to the time an apparent horizon forms. The black hole is then characterized by its horizon mass and spin computed from the apparent horizon with average radius .
The rest-mass of the accretion disk that forms after collapse is computed as,
where the domain of integration excludes the spherical region inside the apparent horizon.
The ejected material is defined by the two conditions,
where is the first lower component of the fluid 4-velocity, and . The first condition in (13) assumes fluid elements follow geodesics and requires that the orbit is unbound. This is a simple criterion we use for continuity with previous work, e.g. East and Pretorius (2012); Hotokezaka et al. (2013a), and should at least capture the correct order of magnitude. The condition requires that the material has an outward pointing radial velocity; it has been used in East and Pretorius (2012) but not in Hotokezaka et al. (2013a). The total ejecta mass is computed as,
where the integral is computed on the region,
on which material is unbound according to (13).
In order to investigate the energetics and geometry of the ejecta we consider different sets of integrals in the -plane, in the -plane, and the full 3D-domain. The kinetic energy of the ejecta can be approximated as the difference between the total energy (excluding gravitational potential energy), and the rest-mass and the total internal energy Hotokezaka et al. (2013a),
where . Additionally, we compute the -weighted integral of ,
and the quantities,
and roughly estimate the geometric distribution of the ejecta. Similar integrals have been proposed in Hotokezaka et al. (2013a), but in that case they were employed in three dimensions. We will use the approximation of the kinetic energy in Sec. VII.1 and discuss the weighted velocities and the for the case-study in Sec. VIII.
We compute the entropy “indicator”,
where and are locally determined by the value of , see Tab. 1. In cases where the additional thermal contribution to the pressure is small , while in presence of shock heating .
Finally, gravitational waveforms are calculated via the curvature invariant and performing multipole decomposition on extraction spheres Brügmann et al. (2008). We work with the metric multipoles , which are reconstructed from the curvature multipoles using the frequency domain integration of Reisswig and Pollney (2011), using the initial circular gravitational wave frequency as a cutting frequency, see Tab. 4. All the waveforms are plotted against the retarded time,
where the extraction radius is .
V Single neutron star tests
|Single star test|
Our conservative AMR implementation has been tested and validated in full-general relativistic simulations of single star spacetimes. In this section, we present five different tests, namely
a static (nonrotating) neutron star with refinement levels inside the star (Sec. V.0.1);
a boosted, nonrotating neutron star crossing refinement levels (Sec. V.0.2);
a migration of an unstable spherical configuration to a stable one crossing refinement levels (Sec. V.0.3);
a uniformly rotating neutron star with refinement levels inside the star (Sec. V.0.4);
a neutron star close to the Kepler limit, which is perturbed and finally disrupted crossing refinement levels (Sec. V.0.5).
For each test we perform simulations for different combinations of the restriction (R), prolongation (P), and correction (C) step, as indicated in Tab. 2. The grid parameters are reported in Tab. 3. The most important quantity we are focusing on is the rest-mass.
For the simulations we use both the BSSN and the Z4c system. Although no major differences between BSSN and Z4c are observed for what concerns mass conservation, Z4c evolutions show overall smaller violations of Einstein constraints.
We investigate a spherical star with a gravitational mass of . The initial data are calculated with a polytropic EOS with and . The star is then evolved with the -law EOS. The grid is prepared such that the finest refinement level is fully contained in the star covering half diameter, and level ends at the star surface. This is shown in the top panel of Fig. 3, which collects the results. Although at the continuum the solution is trivial (static), in numerical simulations some dynamics is observed due to truncation errors. This is mostly triggered by the artificial atmosphere treatment close to the star surface, and by truncation errors on the refinement levels . Thus, differences in the RPC steps influence the overall dynamics of the system.
The middle and bottom panel of Fig. 3 show the relative error in the rest-mass and its time derivative. The conservative AMR (C step) improves the mass conservation of about orders of magnitude, independently on the particular RP choice. Additionally, we observe differences between the RP combinations. Even using C, the 4th order WENO and Lagrangian RP introduce spurious oscillations in the rest-mass derivative (see green and orange solid lines). In general, using the average R leads to the smallest errors.
These results refer to an atmosphere density . We have experimented with the a2e2 RP setup and an atmosphere density of . The result is shown as a black dashed line in the middle and bottom panel of Fig. 3. A lower atmosphere significantly improves the mass conservation. In this test, the error in the rest-mass derivative related to the C step is about , while the one related to the atmosphere treatment is about . Hence, optimal results can only be obtained with a proper combination of RPC and .
Initial data is prepared using the same star model as Sec. V.0.1, which is now boosted in the negative -axis direction111We have further tested our implementation by boosting the star in all the directions, and both applying bitant symmetry, i.e. evolving only , and simulating the full numerical domain.. The grid is prepared such that the star is initially entirely covered by the finest refinement level . During motion, the star crosses completely the two finest refinement levels, as shown in Fig. 4 (top panel).
As visible in Fig. 4 (middle and bottom panels) the C step improves mass conservation in most of the cases, but here its effectiveness depends more significantly on the RP choice than in the TOV test. In particular, the C step is not effective with WENO RP. The a2e2 and a2wz6 schemes perform best, indicating the importance of a conservative R. Similarly to the previous test, we test the role of the atmosphere parameters on the optimal a2e2 setup. Lowering the atmosphere by a factor 100 improves mass conservation by a factor 10 in this case (see dotted black line).
We investigate an unstable single neutron star configuration. Initially the central density is and the gravitational mass , see e.g. Thierfelder et al. (2011a). The star is in an unstable equilibrium, truncation errors trigger a migration to a stable configuration, which involves violent nonlinear oscillations on dynamical timescale. During these expansions and contractions, matter crosses the grid refinement levels. When matter reaches the grid outer boundary some rest-mass falls out of the grid, but typically mass conservation is mostly affected by the interaction between the star low-densities outer layers and the atmosphere. Results are summarized in Fig. 5 for ().
We observe the conservative AMR is effective up to times , that corresponds to bounces of the star core; up to the first bounce the rest-mass conservation improves of about two order of magnitude if the C step is used. At times matter densities reach outer regions, where the resolution is dropped by a factor of and interaction with atmosphere becomes significant.
Figure 6 summarizes our experiments with atmosphere parameters. Lowering by an order of magnitude leads to an improvement of the mass-conservation by approximately one order of magnitude for the beginning of the simulation, while for different and the same the error stays the same, as expected. Relative rest-mass violation can be minimized up to using and . One can notice that, if the C step is not applied and the atmosphere is small enough (), a dramatic mass violation happens as soon as matter crosses the first refinement boundary (), see dotted lines in the figure. The same does not happen with the C step. As time advances, rest-mass conservation is progressively corrupted in all the cases due to the drop in resolutions in the outer region reached by the low-density star outer layers bouncing back and forth.
Initial data is a stable uniformly rotating neutron star described by a polytropic EOS with and , and with , axes ratio , and gravitational mass , e.g. Dimmelmeier et al. (2006). The initial data are computed with the RNS code Stergioulas and Friedman (1995); Nozawa et al. (1998). The star is evolved with the -law EOS, for about 6 periods. Results are shown in Fig. 7.
As in the previous tests, the C step improves the results in many cases; the best RP setup is a2e2. The l4l4 and l4l4n RP perform equally good at late times. Surprisingly, the nonconservative w4w4n RP is here observed to give good results, and at the end of the simulation, it is comparable to a2e2.
Initial data is a rotating neutron star at the Kepler limit modeled by a polytropic EOS with and , and with , axes ratio , and gravitational mass . The star is evolved with the -law EOS with ; the lower polytropic exponent triggers the star expansion with matter crossing several refinement levels.
The left panels of Fig. 8 show how the matter expands along the -axis and the -axis over time, i.e. perpendicular and along the symmetry axis. The right panels of the figure show the mass conservation. The best RPC combinations are again a2wz6 and a2e2.
v.0.6 Summary of single star tests
Summarizing the results of the single star tests, we find the best mass conservation using the a2e2- scheme, i.e. the average restriction operation and a 2nd order ENO interpolation for the prolongation. The a2e2 simulations show, on average, the smallest and no artificial oscillations in . The latter are present in at least one test for all other setups than a2e2. Additionally, the TOV, TOV, and TOV tests suggest that also the artificial atmosphere treatment leads to mass violation. The stability of the simulation improves with higher atmosphere values, but the mass conservation improves for lower atmosphere thresholds. An optimal setup is necessarily a compromise between these two effects. The largest violations of rest-mass conservations are observed in the lowest resolved regions; where the violation becomes independent on the C step and the atmosphere values (i.e. it is mostly due to resolution).
Vi BNS configurations & Grid setup
For this work we have prepared several BNS irrotational configurations in quasiequilibrium and circular orbits; all the configurations are reported in Tab. 4. Initial data are calculated with the LORENE Gourgoulhon et al. () code.
Our BNS sample spans the EOS sample of Tab. 1, and, for each EOS, two mass ratios222We define the mass-ratio to be always . are considered for a fixed total binary mass of . All EOS support maximal neutron star masses in agreement with recent observations Demorest et al. (2010); Antoniadis et al. (2013), and the adiabatic speed of sound is for a density range up to the maximum density supported by a stable TOV star. The compactnesses of the stars lie within . The tidal coupling constant spans (see (23) for the definition). Notice that stiffer EOS have larger , and, for the same EOS and , a larger implies a larger . The initial GW frequency of all binaries is .
Additionally, we computed a and configuration with the MS1b EOS (MS1b-100150). MS1b-100150 has a highly deformable EOS and . The choice of parameters (EOS, , ) of this configuration could be considered as “extreme” given the double pulsar population, e.g. Kiziltan et al. (2013). However, the double pulsars sample is rather small to be a significant statistics and the MS1b-100150 parameters are possible.
Some of our configurations have already been investigated in full general relativity in Hotokezaka et al. (2013a, b). Thus, the choice of initial data allows us to compare results with the literature. We will also compare with results of Bauswein et al. (2013) employing smooth particles hydrodynamics and conformal flatness, although the evolution method differs and initial data are not prepared in the same way as here.
For the BNS evolutions we use the Z4c scheme Bernuzzi and Hilditch (2010); Hilditch et al. (2013a), and constraint preserving boundary conditions Ruiz et al. (2011); Hilditch et al. (2013a). For all our runs a grid consisting of refinement levels is used, levels with are dynamically moved. The grid spacing and outer boundary position depends on the employed model and is reported in Tab. 4. For refinement level we employ the spherical patches, as described in Sec. III; but we do not evolve matter on them. Indicating with the inradius on refinement level , GRHD is evolved up to with resolution , and up to with resolutions for .
For all binary evolutions we run both the a2e2 and a2e2n RPC schemes. The a2e2 scheme is chosen because of its robustness and best performances in our previous tests; a2e2n is considered in order to assess the effect of the C step in the AMR strategy. We have not considered other combinations due to the computational overhead that they would imply. We set and for all simulations.
|H4-135135-R2c||1804||0.129||1.54||HMNSBH||5130 (25)||0.146||1.75||0.214||2.57||0.6||0.3 (0.5)||10.8||2.48||0.62|
|H4-135135-R2n||1803||0.130||1.55||HMNSBH||4470 (22)||0.145||1.73||0.216||2.58||0.6||0.3 (0.6)||8.5||2.54||0.65|
|ALF2-135135-R2c||2148||0.142||1.71||HMNSBH||3760 (19)||0.168||2.01||0.235||2.81||3.5||0.4 (0.7)||17.8||2.43||0.62|
|ALF2-135135-R2n||2145||0.142||1.71||HMNSBH||3770 (19)||0.165||1.98||0.230||2.75||2.0||0.4 (0.7)||21.1||2.44||0.63|
|SLy-135135-R2c||2504||0.168||2.01||HMNSBH||2159 (11)||0.206||2.46||0.292||3.49||12.2||4.0 (7.1)||8.4||2.48||0.64|
|SLy-135135-R2n||2495||0.168||2.01||HMNSBH||2577 (13)||0.207||2.48||0.290||3.47||14.2||5.9 (10.5)||9.6||2.49||0.64|
|SLy-125145-R2c1||2353||0.162||1.93||HMNSBH||3020 (15)||0.184||2.20||0.286||3.42||6.5||2.8 (5.1)||17.9||2.40||0.58|
|SLy-125145-R2n1||2350||0.161||1.93||HMNSBH||2870 (14)||0.187||2.24||0.283||3.39||4.5||1.7 (3.0)||14.5||2.46||0.61|
|SLy-125145-R2c2||2350||0.161||1.92||HMNSBH||3310 (16)||0.186||2.23||0.285||3.41||6.2||2.1 (3.7)||18.4||2.40||0.58|
|SLy-125145-R2n2||2348||0.160||1.91||HMNSBH||2180 (11)||0.184||2.20||0.283||3.39||5.4||2.5 (4.5)||11.1||2.49||0.62|
Vii BNS mergers with and
We first discuss our results for configurations with mass ratios and a total binary mass . We focus on the effect of mass ratio and the EOS on the merger dynamics, ejecta and gravitational waves. Also, we show the use of conservative AMR significantly improves the simulation of the merger remnant. Several results are reported in Tab. 5, and collected in Fig. 9, Fig. 10, Fig. 11 and Fig. 12, to which we refer during the discussion.
vii.1 Effect of EOS and on merger dynamics
The initial configurations, prepared in quasicircular orbits at the same GW frequency , evolve for about 3 to 5 orbits before merger, depending on the EOS and mass ratio . Here, the moment of merger is defined as the time corresponding to the peak of the multipole of the GW amplitude (see below). Stiffer EOSs give shorter inspiral (less revolutions) and lower dimensionless GW frequency at merger, see in Tab. 5. Unequal-mass systems are characterized by slightly shorter inspiral than equal-mass ones and smaller merger frequencies of about . These properties can be understood considering the values of the main () tidal polarizability parameter (tidal coupling constant hereafter) Damour and Nagar (2010),
where are the dimensionless Love numbers of the individual stars Hinderer (2008); Damour and Nagar (2009); Binnington and Poisson (2009); Hinderer et al. (2010), in our sample. The results agrees with the analysis of Bernuzzi et al. (2014b). Essentially, for the same mass, stars with stiff EOS have larger radii than those with soft EOS, and attractive tidal interactions are stronger for larger values of ; thus, stiffer EOS binaries merge at lower frequencies. Notice that: (i) configuration have slightly larger values of than ; (ii) in our sample of configurations, EOS effects are typically larger than mass-ratio effects. The late-inspiral dynamics and GWs have been subject of recent work, e.g. Bernuzzi et al. (2012); Bernuzzi et al. (2014c) and we do not discuss them any further here. In the following we focus on the postmerger dynamics.
The postmerger dynamics has a rich phenomenology depending on the main binary properties: total mass, mass-ratio, EOS and stars’ spin (see e.g. Shibata and Taniguchi (2006); Shibata et al. (2011); Hotokezaka et al. (2013a, b); Bernuzzi et al. (2014a); Kastaun and Galeazzi (2014) for recent work). In the case of irrotational binaries and , equal-mass mergers result in a massive differentially rotating compact object, which oscillates violently (see the evolution in Fig. 10 right after merger). The compact object’s angular momentum is redistributed from the inner region to outer ones by torque and nonlinear hydrodynamical interaction. The stability of the object depends on the mass, EOS and dissipative processes (see below). Following the literature Baumgarte et al. (2000), we define this object as a hypermassive neutron star (HMNS), in case its rest-mass is larger than the maximum rest-mass of a stable uniformly rotating star described by the same EOS, or a supramassive neutron star (SMNS), in case its rest-mass is smaller. If the object does not exceed the rest mass of a stable TOV-solution, we simply refer to it as massive neutron star (MNS). These definitions apply to equilibrium configurations, in particular to cold EOS and axisymmetry, hence, although of common use, they cannot be rigorously applied to the merger remnants. In most cases HMNS are objects that are dynamically unstable and collapse to a black hole on timescales of ms; whereas SMNSs are objects that appear stable on those timescales, but may eventually collapse later on due dissipative processes, e.g. loss of angular momentum radiated via GWs. Snapshots of the density distribution and velocities in the orbital plane are presented in Fig. 9; the simulation time is close to the moment of merger.
Three of our configurations, H4-135135, ALF2-135135, and SLy-135135, merge in a HMNS which collapses to a black hole (BH) within ms from the merger moment. All these EOSs support approximately the same maximum mass regarding single spherical stars, but the stiffer the EOS, the longer is . This fact can be understood by the following considerations. The range for the tidal coupling constant is , where soft (stiff) EOS binaries correspond to small (large) values in this range. Stiff EOS binaries are gravitationally less bound systems than soft EOS binaries: their binding energy at merger is larger (less negative) as well as the angular momentum. As a result, the HMNS has more angular momentum support at formation. However, the initial angular momentum is not the only factor that determines the lifetime of the HMNS. At formation, the HMNS density in the star core increases, the pressure response depends on the (effective) adiabatic index of the fluid which is different for each EOS. As a result, the HMNS nonlinear oscillations and the efficiency of the angular momentum redistribution depend on the EOS Hotokezaka et al. (2013b). Stiffer EOSs have larger pressure support against gravity, especially at high densities. Finally, in a more realistic situation than the one simulated here (and on longer timescales), thermal support, angular momentum transport driven by magnetic fields333We notice the largest simulations with present techniques and resolutions have not properly resolved magnetic field amplification effects Kiuchi et al. (2014). and cooling mechanisms (neutrinos) are expected to play a role. The lifetimes of the HMNS are stated in Tab. 5 and our results agree with Hotokezaka et al. (2013b) within ms.
The merger of MS1-135135, differently from the other configurations, produces a differentially rotating object that is stable over the whole simulation time, i.e. ms after merger. Non-rotating stars described by the MS1 EOS can support a maximum rest-mass of . According to the previous definition, we classify the merger remnant for the MS1 models as a MNS. Considering the physics simulated here, we expect that the merger remnant will stabilize via GW emission reaching a uniformly rotating and cold configuration on the characteristic timescale, ms.
The unequal-mass configurations H4-125145 and ALF2-125145 have a different merger remnant than the corresponding configurations. In these cases we find an object stable over ms, but since the mass is still larger than the supported mass of the uniform rotating model, remnants are HMNSs. We expect these configurations will collapse within . References Hotokezaka et al. (2013a, b) found that similar configurations with a slightly different thermal component form BHs. A similar dynamics as for the MS1-135135 is observed in the merger of MS1-125145, where a stable MNS is produced. The SLy-125145 forms, as in the case, a black hole, but, following the general trend, the HMNS lifetime is longer.
Due to the unequal mass ratio, the merger remnant is typically more deformed than the corresponding and strongly non-axisymmetric at formation, see Fig. 9. Unequal-mass binaries have more stable merger remnants than corresponding equal-mass ones (e.g. larger ). The HMNS/MNS are characterized by slightly larger radii than the ones, and a different central density, Fig. 10. Additionally, the mass-ratio has an effect on the ejecta as we shall see below.
At formation, all the merger remnants show violent oscillations, visible in the evolution of in Fig. 10. The softer the EOS, the larger are the oscillations, see in particular the SLy panels in the figure. This property reflects the pressure response of the EOS for density jumps around (Cf. above and also Hotokezaka et al. (2013b)). The oscillations have a quasi-radial character, and relax either within few radial periods or before the onset of collapse.
In cases with BH formation, the BH masses are of order , and the dimensionless BH spin is of the order for all the configurations. The evolution of the BH parameters is presented in Fig. 11 (top and middle panels). These results suggest that, in this scenario, the BH formation and properties are mostly determined by the total mass of the system and depend only weakly on other details. However, uncertainties on these numbers are of the order of , and it is difficult to draw precise conclusions.
The final BH is surrounded by an accretion disk of rest-mass , see Tab. 5 and the bottom panel of Figure 11. The disk geometry is essentially axisymmetric for all the configurations. During the evolution, the maximum density inside the disk decreases from to . At the moments the BH masses and spins reach their plateaus (late times in our simulations), the dense regions of the disk extend up to distances km. Lower density, gravitationally bound regions larger than extend up to km. The accretion rate is of the order .
vii.2 Assessment of conservative AMR
Before continuing the analysis of the physical properties of the merger remnant we discuss here the accuracy improvements due to the numerical algorithm described in Sec. III.2. Figure 10 reports result obtained with (solid lines) and without (dashed lines) the C step in the AMR algorithm. The information of the figure is complemented with the entries of Tab. 5. For all the configurations the C step is crucial for the simulation accuracy after merger.
Let us first discuss rest-mass conservation. As pointed out in the introduction, the rest-mass can in general increase or decrease. In our BNS simulations we identify two main and competitive causes for the violation of conservation: (i) when fluid crosses refinement boundaries rest-mass tends to increase, (ii) the artificial atmosphere treatment tends to decrease the rest-mass. Clearly, the C step can improve only violations of type (i).
For most of the configurations the use of the C step leads to an improvement of a factor of , except for the MS1-135135-R2 configuration where an improvement by more than a factor of is observed. The only simulation were no significant improvement is observed is SLy-125145-R2, where the violation is from merger to the end of the run.
Overall, the data show some dependence on the EOS. Without C step the mass conservation is in general better for softer EOS; this is probably related to the smaller star deformations. On the contrary, with C step, slightly larger errors are observed for softer EOS. We suggest that this is caused by the influence of numerical viscosity, that, in these runs, is more significant than in the runs without C step due to better overall conservation. Notice that the performance of the conservative AMR algorithm is always better than (or at most comparable to) the corresponding simulations without C step. In Bernuzzi et al. (2014a) we have employed larger grid boxes without C step in an attempt to optimize the performances of the nonconservative AMR for the remnant simulation. In Appendix A we present some experiments along this line showing that conservative AMR is, in general, a better strategy.
Mass-violations influence the behavior and lifetime of the merger remnant, as evident from Fig. 10. We observe systematic shifts in the collapse time of several HMNS although there are no qualitative differences due to the sufficiently high grid resolutions of our runs. For H4-125145 the mass violation in the outer layers in the H4-125145-R2n run determines a slightly different evolution of the MNS and a lower .
We observe maximal differences of a factor of 3 in the ejecta mass measured on the coarsest level () between the runs with and without the C step. Figure 10 (bottom panel for each EOS) shows that the differences is larger either shortly after merger time or at later times: no clear trend is identifiable. Thus, low density ejecta remain challenging to simulate even with conservative AMR (as long as nested boxes are used as opposed to local AMR tracking the ejecta). In particular, the artificial atmosphere influence is probably significant: (i) during inspiral we observe some spurious ejecta due to atmosphere fluctuation, and (ii) at late times, when ejecta have expanded into larger radii (coarser resolutions) we expect an effect as the one discussed for the TOV test in Sec. V.0.3.
Differences in the black hole and disk remnant are also observed, see Tab. 5 and Fig. 11. If the C step is not employed the estimated disk mass changes up to . In all configurations the final black hole mass and spin is overestimated when no C step is applied, which is probably related to the increase of the rest-mass visible in the upper panels (for each EOS panel) of Fig. 10 (dashed lines).
Finally, we mention that the GWs calculation during the inspiral is basically not influenced by the use of the C step. This is due to the fact that we have not attempted to refine the grid inside the star during that phase. During orbital motion the stars stay compact and there is no need of further improving mass conservation. GWs in the post merger reflect the slightly different dynamics, but the characteristic frequencies (see below) are essentially unaffected.
In the following we will discuss exclusively the results obtained with conservative AMR scheme.
In this section we discuss the EOS and mass-ratio effect on the dynamical ejecta. A detailed analysis of the dynamical formation of the ejecta will be presented in Sec. VIII.
Figure 10 shows the evolution of the ejecta mass for the various configurations; Tab. 5 reports the maximum value. Ejecta peaks happen during and shortly after the merger moment; the ejecta rest-masses at this time are about , and in some cases reach .
The amount of ejected material depends on the EOS and on the mass ratio. If larger ejecta are observed for softer EOS. For a given EOS (but except for SLy EOS), configurations have larger ejecta than ones. Similarly, the kinetic energy estimate computed according to Eq. (16) is larger for softer EOS and larger than for stiff ones. Our results for MS1, H4, and ALF2 configurations essentially agree with Hotokezaka et al. (2013a); Bauswein et al. (2013).
We stress that ejecta computations are challenging. At present, mass conservation and artificial atmosphere are the main factors limiting the accuracy. This is evident in the case of SLy configurations. The results in Sec. VII.2 suggest that the evolution of this soft EOS is less accurate than the others (see also discussion in Hotokezaka et al. (2013a)). We believe this is the reason why is larger for SLy-135135 than for SLY-125145. The poor mass conservation in SLY-125145 certainly affects the ejecta calculation. Notice also that a similar setup as SLy-135135 has been evolved in Bauswein et al. (2013); there, the ejecta mass was estimated to lie in the range between .
vii.4 Gravitational waves
The multipoles of the GWs are shown in Fig. 12 for all the configurations. For each EOS, each panel shows the real part of the wave (top) and the instantaneous GW frequency (bottom). The vertical line in each panel marks the moment of merger, defined as the peak of the amplitude .
The emission from the orbital motion is the characteristic chirping signal, in which frequency and amplitude monotonically increase. At these separations, the dynamics is strongly affected by tidal interactions (parametrized by ), and the GWs phase carry information about the EOS. A detailed and accurate semi-analytical modeling of the inspiral up to merger has been given in Bernuzzi et al. (2014c). The chirp signal ends at the amplitude peak.
After the merger moment, the amplitude instantaneously drops down since the two stars merge in a single body which has, for one instant, a quasispherical geometry Thierfelder et al. (2011a) (see also the frequency spikes). The postmerger signal is mainly characterized by the nonlinear oscillations of the merger remnant. As discussed above and elsewhere, e.g. Stergioulas et al. (2011), the merger remnant can be approximated by a compact star oscillating nonlinearly at the proper frequencies. The -mode with frequency is the most efficient emitter of GWs, and it is strongly excited at formation. Thus, the GW emitted by the HMNS/MNS is dominated by this frequency. Looking at the frequency in Fig. 12, large oscillations are present right after the merger moment and correspond to the very nonlinear phase described in Sec. VII.1; softer EOS show larger oscillations. During early stages of the HMNS/MNS evolution, different modes are excited, see e.g. the spectrogram in Bernuzzi et al. (2014a). Nonlinearity results in mode couplings, the main ones being the combination between the quasiradial mode and the Stergioulas et al. (2011). In cases where a MNS is formed (MS1 EOS), the frequency oscillations relax quickly; the power in the channels decreases, and the frequency essentially settles on the mode. In cases where a HMNS is formed, the frequency monotonically increases as a result of the star contraction prior to collapse.
Let us finally discuss the GW spectra shown in Fig. 13. The figure includes, for each configuration, the spectrum of the entire signal as a thick line and the spectrum considering only the signal for as a thin line. Some of the relevant frequencies are marked with bullets: the frequency at the waveform amplitude peak (triangles), a frequency related to a secondary postmerger peak (circles), and the frequency corresponding to the main postmerger peak (diamonds). Recently, there has been intense research about the identification and characterization of this postmerger GW spectrum frequencies Bauswein and Janka (2012); Stergioulas et al. (2011); Hotokezaka et al. (2013a, b); Bauswein et al. (2014); Takami et al. (2014); Bauswein and Stergioulas (2015). For most of the configurations, the frequency is clearly identifiable. Note however the double peak for the MS1 models.
The origin of the secondary peak is not well understood. appears mostly related to the very late inspiral phase: several peaks are not present, or strongly suppressed, if the PSD is computed using only times . However, for configuration SLy135135 and H4135135 one can notice a clear secondary peak also in the PSD of the signal at times . We observe that the peaks generated by signals at times are suppressed for unequal-mass configurations (). Our values of are in good agreement with the frequencies called in Takami et al. (2015). Our PSD analysis might be compatible with the interpretation of Bauswein and Stergioulas (2015) according to which the peak of the spectrum close to is due to two different effects: the nonlinear mode coupling (that can be extracted clearly using the signal only), and motion of spiral arms during the last stage of the merger process (but mostly at times ) at a frequency called there .
Viii The MS1b-100150 configuration