Merger of Compact Stars in the Two-families Scenario

The first detection of gravitational waves from the merger of two compact stars⁸ in 2017 August (Abbott et al. 2017) represents a breakthrough for the astrophysics of compact objects and the physics of dense nuclear matter. The gravitational wave signal associated with the process of inspiral of two stars encodes information on the average tidal deformability $\tilde{{\rm{\Lambda }}}$ of the binary, which is strongly dependent on the stiffness of the equation of state (EoS) of dense matter. Several recent studies have exploited the limits on $\tilde{{\rm{\Lambda }}}$ to obtain constraints on the radii of compact stars (Abbott et al. 2018; Burgio et al. 2018; Fattoyev et al. 2018; Lim & Holt 2018; Most et al. 2018). The electromagnetic counterparts of GW170817, i.e., the short gamma-ray burst (sGRB) GRB170817A and the kilonova AT2017gfo, allow one to further constrain the EoS by indicating that the most probable merger's remnant is a hypermassive compact star and requiring that the mass released and powering the kilonova is of the order of 0.05 M_⊙ (Bauswein et al. 2017; Margalit & Metzger 2017; Annala et al. 2018; Radice et al. 2018b; Rezzolla et al. 2018; Ruiz et al. 2018). The main conclusion of those studies can be summarized by stating that the EoS of dense matter cannot be very stiff: EoSs leading to radii larger than about 13.5 km are basically ruled out. Moreover, the maximum mass of the nonrotating configuration should be less than approximately $\sim 2.2{M}_{\odot }$, although still with a large error bar of the order of $0.2{M}_{\odot }$. This is an interesting result: indications of radii larger than approximately 13.5 km would indeed point to a stiff EoS in which only nucleons are present, while for smaller radii, nonnucleonic degrees of freedom could appear in compact stars.

The relevant degrees of freedom for the composition of a compact star can be inferred from the value of the radius of a star having a mass of about $1.5{M}_{\odot }$. First, if the radius is large, 13 km ≲ R_1.5 ≲ 13.5 km, the density at the center is low enough that nonnucleonic degrees of freedom appear only in the most massive stars (Lonardoni et al. 2015). Second, if 11.5 km ≲ R_1.5 ≲ 13 km, resonances (hyperons and/or deltas) will appear even for stars having masses of the order of $\sim 1.5{M}_{\odot }$ (Maslov et al. 2015). Also, deconfined quarks can appear at the center of the star, which becomes a so-called hybrid star (HybS; Nandi & Char 2018). Finally, for R_1.5 ≪ 11.5 km, only disconnected solutions of the Tolman–Oppenheimer–Volkov equation exist (Alford et al. 2013; Burgio et al. 2018). These solutions are characterized by the presence of a mass window, inside which two different configurations are possible: one made only of hadronic matter and the other made partially or totally of deconfined quark matter. These configurations correspond either to the twin-stars scenario (Alford et al. 2013; Most et al. 2018; Paschalidis et al. 2018) or to the two-families scenario (Drago et al. 2014a, 2016b; Drago & Pagliara 2016; Wiktorowicz et al. 2017).

The main differences between these two scenarios concern the mass–radius relation and the microphysics embedded in the EoS. In the twin-stars scenario, the stars containing quarks are both the most massive and the ones with the smallest radius, which reaches its minimum at the maximum mass configuration. On the other hand, in the two-families scenario, the stars having the smallest radius belong to the hadronic branch (hadronic stars (HSs)), and they have a mass of the order of $(1.5\mbox{--}1.6){M}_{\odot }$. Instead, the second family, which corresponds to strange quark stars (QSs), is characterized by not-too-small radii (of the order of or larger than about 12 km) and a maximum mass ${M}_{\mathrm{TOV}}^{Q}$ that can exceed $2{M}_{\odot }$. Concerning the EoS, the twin-stars scenario assumes the existence of a strong first-order phase transition from hadronic matter to quark matter (often modeled via a Maxwell construction), whereas in the two-families scenario, one assumes the validity of the Bodmer–Witten hypothesis (Bodmer 1971; Witten 1984), which implies the existence of a global minimum of the energy per baryon at a finite density and for a finite value of strangeness.

The investigation of these scenarios, as well as their consequences for the interpretation of GW170817, has been discussed in previous papers (Burgio et al. 2018; Gomes et al. 2019; Nandi & Char 2018; Paschalidis et al. 2018). In particular, in Burgio et al. (2018), it was shown that the observational constraints on $\tilde{{\rm{\Lambda }}}$ would imply that the radius of the $1.4{M}_{\odot }$ configuration, R_1.4, is larger than about 12 km within the one-family scenario, while when considering two families of compact stars, R_1.4 could be significantly smaller, down to about 10 km.

An important point we need to stress here is that within the two-families scenario, QSs and HSs coexist. A way to classify compact objects is to look at their masses and/or radii: in a few cases, the observational indications allow one to establish to which of the two families a specific object should belong (Char et al. 2019). Once the object has been identified as a QS or an HS, one can test the model by checking whether all of the other properties are consistent with that classification. For example, the X-bursts associated with magnetar crustal modes that have been observed in a few cases are difficult to explain if those objects are QSs (see Watts & Reddy 2007), and therefore in our scheme they must be HSs. This prediction can be tested if the radius of those stars will be measured in the future. Concerning glitches, it is not yet clear whether QSs can produce such phenomena. A possibility is that at least part of the QS is in the color crystalline phase, which is rigid as well as superfluid. These are the two basic ingredients needed to model glitches; see Anglani et al. (2014) and Mannarelli et al. (2014).

An important criticism against the coexistence of QSs and HSs was put forward in Madsen (1988), where it was argued that if QSs exist, then all compact stars are QSs. This argument is based on the expectation of an abundant galactic strangelets pollution produced during the merger of two QSs. Such a flux of strangelets would trigger the conversion of all compact stars into QSs. Actually, this argument was criticized in Wiktorowicz et al. (2017), where it was shown that (i) the produced strangelets have a very large mass (order of 10⁴⁰ baryons), and thus the chance of capture by compact or main-sequence stars is very small; and (ii) smaller strangelets most likely evaporate into nucleons due to the dynamics of the merger and the strong reheating of matter after the merger; see discussion in Section 7 and N. Bucciantini et al. (2019, in preparation).

In this paper, we extend a previous preliminary work on the predictions of HS–HS mergers within the two-families scenario (Drago & Pagliara 2018). After an overview of the phenomenology of mergers within the two-families scenario, we present numerical simulations of the process of the merger of two HSs by using the Einstein Toolkit, an open-source, modular code for numerical relativity (Loffler et al. 2012). We consider two EoSs: a soft one, which contains hyperons and deltas and leads to HSs belonging to the first family, and, as a reference, a stiffer and purely nucleonic one. The main outcome of the simulations is a precise estimate of the threshold mass above which a prompt collapse occurs for the merger of two HSs. We also provide estimates of the mass ejected, which is an important quantity for the phenomenology of kilonovae. Also, we discuss in which region of the star quark matter nucleation can take place. Finally, we estimate the various possible merger processes within the two-families scenario by using the population synthesis code Startrack (Belczynski et al. 2002, 2008; Wiktorowicz et al. 2017).

The two-families scenario is based on the idea that in the hadronic EoS, delta resonances and hyperons do appear at large densities. To model the EoS, we adopt the relativistic mean field model SFHo of Steiner et al. (2013), with the inclusion of delta resonances and hyperons; see Drago et al. (2014b). Here we use the EoS SFHo-HD already investigated in Burgio et al. (2018) that corresponds to a coupling of delta resonances with the σ meson of x_σΔ = 1.15 (see Drago et al. 2014b; Burgio et al. 2018). A detailed motivation for this choice of the hadronic EoS will be presented in the Appendix.

For quark matter, we use the simple bag-like parameterizations of Alford et al. (2005) and Weissenborn et al. (2011), where an effective bag constant B_eff has been introduced together with a parameter a₄ that encodes pQCD corrections. By setting ${B}_{\mathrm{eff}}^{1/4}=137.5\,\mathrm{MeV}$, a₄ = 0.75, and the mass of the strange quark ${m}_{s}=100\,\mathrm{MeV}$, we obtain ${M}_{\mathrm{TOV}}^{Q}\sim 2.1{M}_{\odot }$.

2.1. Transition from Hadronic to Quark Matter

To model the formation of quark matter, we adopt the scheme based on quantum nucleation developed in many papers; see, e.g., Heiselberg et al. (1993), Iida & Sato (1998), Berezhiani et al. (2003), Bombaci et al. (2004), and Niebergal et al. (2010). The basic idea is that the process of formation of the first drop of quark matter preserves the flavor composition, and it is therefore impossible to nucleate strange quark matter if hyperons or kaons are not already present in the hadronic phase. In Figure 1 we show the fraction Y_i of baryons as a function of the baryon density n_b. Note that, as n_b increases, more and more hyperons are produced, making the conversion to strange quark matter more likely. Let us define the strangeness fraction ${Y}_{S}=({n}_{{\rm{\Lambda }}}+2({n}_{{{\rm{\Xi }}}^{0}}+{n}_{{{\rm{\Xi }}}^{-}}))/{n}_{b}$. A possible way to define a threshold for the conversion of hadronic matter into strange quark matter is to request that d_s, the average distance of strange quarks inside the confined phase, is of the order of or smaller than the average distance of nucleons in nuclear matter, ${d}_{n}={n}_{0}^{-1/3}$, where n₀ is the nuclear matter saturation density, n₀ = 0.16 fm⁻³. The reason is that strange quarks should be close enough to be able to mutually interact in order to help the nucleation of a first drop of deconfined quark matter. This way of defining a critical density is simple and intuitively corresponds to a necessary condition, but it is less quantitative than the one based on discussing the minimal size of the stable drop of quark matter used in many papers, including the ones cited above. One should also notice that at the high temperature reached during the merger, thermal nucleation could be possible, and it could make the whole process of quark deconfinement significantly faster (e.g., in Di Toro et al. 2006, the temperature above which thermal nucleation becomes dominant has been estimated to be around 60 MeV).

Zoom In Zoom Out Reset image size

In Figure 1 (upper panel), we show the fraction of the strange and not-strange baryons, together with the fraction of strangeness, at zero temperature. The previous condition on the minimal density of strangeness is satisfied above a baryon density of the order of (0.9–1) fm⁻³, for which Y_S ∼ 0.2−0.3 (see the two points on the dashed line in the figure) and which corresponds to about (6–7) n₀. Note that it is a density significantly larger than the threshold density of the formation of hyperons, which slightly exceeds 0.6 fm⁻³, and from this viewpoint, our conditions agree with the analysis of, e.g., Bombaci et al. (2004, 2009), indicating that the first droplet of quark matter is nucleated at densities larger than the threshold density of hyperon formation. It is important to notice that the condition for quark nucleation at T = 0 is reached at densities smaller than those corresponding to the maximum mass mechanically stable, ${M}_{\mathrm{TOV}}^{H}$. Following Bombaci et al. (2004), we therefore introduce a separate notation ${M}_{\max }^{H}$ for the maximum mass of an HS, stable with respect to quark nucleation.

A specific feature of the two-families scenario and the quark deconfinement mechanism is the possibility of also forming (hot) HybSs as a result of the partial turbulent conversion of hadronic matter into quark matter. These stellar objects live for a time of the order of 10 s, which is necessary for hadronic matter to convert completely into quark matter within the diffusive regime (Drago & Pagliara 2015), and they are characterized by a two-phase structure with a hot quark matter core and a hadronic matter layer (initially cooler) that are in mechanical equilibrium but out of thermal and chemical equilibrium. It is possible to construct a model for the EoS describing these transient HybSs by using Coll's condition (Coll 1976) to define the density separating the low-density hadronic region from the high-density region already made of quarks; see Drago & Pagliara (2015) for details. The reason to discuss these short-living configurations here is that the post-merger remnant, for the cases of HS–HS and HS–QS mergers, can indeed form such HyBs in which the conversion is still proceeding (we label it Coll-hyb). Notice that since QSs are generated in this scenario by the conversion of HSs, they have a minimum mass (${M}_{\min }^{Q}$) that is fixed by the maximum mass of the HSs: the conservation of baryon number during the conversion implies that the total baryon number of the newly born QS is equal to the total baryon number of the progenitor HS. Since the QS is more bound than its progenitor HS, its gravitational mass is of the order of $0.1{M}_{\odot }$ smaller than the gravitational mass of the progenitor, and its minimum value is $\sim (1.35\mbox{--}1.45){M}_{\odot }$, depending on the value of the maximum mass of HSs; see Wiktorowicz et al. (2017). In turn, this implies that, at fixed gravitational mass, QSs are always larger than HSs; see red and black solid lines in Figure 2.

Zoom In Zoom Out Reset image size

Another distinctive feature of the two-families scenario concerns the value of the (square) speed of sound, c²_s, in dense matter: in the hadronic EoS, its value remains significantly below ∼0.4 due to the several softening channels associated with the appearance of hyperons and deltas (see red line in the lower panel of Figure 1). For the quark matter EoS, within the model adopted here, at increasing density, it approaches from below the asymptotic limit of the noninteracting massless gas one-third (the conformal limit); see green line in the lower panel of Figure 1. Nevertheless, it is possible to fulfill the 2 M_☉ limit by assuming that such massive stars are QSs.⁹ On the other hand, within the one-family scenario, the 2 M_☉ limit requires values of the speed of sound that are significantly larger than the asymptotic limit, as was observed in Alford et al. (2015), Bedaque & Steiner (2015) and Chamel et al. (2013). For instance, by using the SFHo EoS, one can obtain a maximum mass of $2.06{M}_{\odot }$ (Steiner et al. 2013) with a corresponding value of the central density of ∼1.2 fm⁻³ and ${c}_{s}^{2}\sim 0.8$; see the black line in the lower panel of Figure 1. Similarly, within the twin-star scenario, the quark matter EoS is assumed to have a large value of the speed of sound, ${c}_{s}^{2}\geqslant 0.8$ in Paschalidis et al. (2018), whereas Christian et al. (2018), Montaña et al. (2019), and Alford et al. (2019) assumed ${c}_{s}^{2}=1$. Since at asymptotically large densities, one should anyway recover the limit of ${c}_{s}^{2}=1/3$, one would need to clarify which are the physical mechanisms responsible for an initial increase of the speed of sound to values close to the causal limit and then its decrease to the conformal limit.

The fate of the merger of two compact stars, in particular the properties of the post-merger remnant and its emissions in gravitational and electromagnetic waves, is determined in general by two quantities: the total mass above which a prompt collapse to a black hole (BH) takes place, M_threshold, and the maximum mass of the supramassive configuration, M_supra. Both quantities strongly depend on the EoS. The fact that the outcome of a merger is (mainly) determined by those two quantities is true both in the one-family and in the two-families scenario, but in the latter, the situation is more complicated, because those quantities need to be determined for both families (and also for the short-living Coll-hyb configuration).

Let us first present the results for the structure of static and rotating stars when adopting the EoSs previously introduced. Several works have been dedicated to the study of the dependence of M_threshold on the stiffness of the EoS: a very interesting result, presented in Bauswein et al. (2013a, 2016) and Bauswein & Stergioulas (2017) and based on explicit numerical simulations, is that the ratio $k={M}_{\mathrm{threshold}}/{M}_{\mathrm{TOV}}$ scales linearly with the ratio ${M}_{\mathrm{TOV}}/{R}_{\mathrm{TOV}}$, i.e., with the compactness of the maximum mass configuration. Depending on the EoS, k varies between 1.3 and 1.6. In Section 5 we will discuss how to estimate this parameter from direct numerical simulations of HS–HS mergers. Concerning M_supra, many studies have found that ${M}_{\mathrm{supra}}\sim 1.2{M}_{\mathrm{TOV}}$ (Lasota et al. 1996; Breu & Rezzolla 2016) for the case of stars with a crust (i.e., HSs and HyBs), whereas for QSs, ${M}_{\mathrm{supra}}\sim 1.4{M}_{\mathrm{TOV}}$ (Gourgoulhon et al. 1999; Stergioulas 2003).

We have computed both static and Keplerian configurations by using the RNS code of Stergioulas & Friedman (1995). Results are presented in Figure 2. For HSs, QSs, and neutron stars (NSs), we confirm the standard results previously described concerning the relation between ${M}_{\mathrm{TOV}}$ and M_supra. We have also computed ${M}_{\mathrm{supra}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$ for our Coll-hyb configurations, which turns out to be of the order of 2.6 M_⊙, respecting the general relation between ${M}_{\mathrm{TOV}}$ and M_supra of compact stars with a crust. As discussed above, while Coll-hyb configurations are chemically and thermally out of equilibrium, they represent a necessary intermediate stage of the evolution of the post-merger remnant in the two-families scenario. Their stability is therefore important both for HS–HS mergers and for HS–QS mergers. We have not computed the maximum mass of a differentially rotating Coll-hyb configuration ${M}_{\mathrm{threshold}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$ through direct numerical simulations, but we can safely assume that ${M}_{\mathrm{threshold}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}\geqslant {M}_{\mathrm{supra}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$.

For the sake of the following discussions, let us summarize the possible values of the masses discussed above. First, let us introduce the quantity ${M}_{\mathrm{tot}}$, which is the sum of the gravitational masses of the two objects undergoing the merger:

$$\begin{eqnarray}&&{M}_{\mathrm{tot}}={M}_{1}+{M}_{2}.\end{eqnarray} \tag{ 1 }$$

When discussing the critical masses, one has to be careful to distinguish between ${M}_{\mathrm{tot}}$ and the gravitational mass that remains after mass ejection (and gravitational wave emission) from the system, ${M}_{\mathrm{remnant}}$. The relation between these two masses is of the form

$$\begin{eqnarray}&&{M}_{\mathrm{remnant}}=(1-\alpha ){M}_{\mathrm{tot}},\end{eqnarray} \tag{ 2 }$$

where α is of the order of a few percent. Plausible values of the masses discussed above are

• 1.  
  HSs: ${M}_{\mathrm{TOV}}^{H}\sim {M}_{\max }^{H}\sim 1.6{M}_{\odot }$ ${M}_{\mathrm{supra}}^{H}\sim 1.9{M}_{\odot }$ ${M}_{\mathrm{threshold}}^{H}\sim 2.5{M}_{\odot }$;
• 2.  
  QSs: ${M}_{\mathrm{TOV}}^{Q}\sim 2.1{M}_{\odot }$ ${M}_{\mathrm{supra}}^{Q}\sim 3{M}_{\odot }$ ${M}_{\mathrm{threshold}}^{Q}\sim {M}_{\mathrm{supra}}^{Q}\sim 3{M}_{\odot };$ ¹⁰
• 3.  
  Coll-HyBs: ${M}_{\mathrm{TOV}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}\sim {M}_{\mathrm{TOV}}^{Q}\sim 2.1{M}_{\odot }$ ${M}_{\mathrm{supra}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}\sim 2.6{M}_{\odot }$ ${M}_{\mathrm{remnant}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}=(1-\alpha$) ${M}_{\mathrm{threshold}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}\sim 1.5\,{M}_{\mathrm{TOV}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}\sim 3.1{M}_{\odot }$.

Concerning HSs, within the two-families scenario, there is an uncertainty on ${M}_{\mathrm{TOV}}^{H}$ (whose value determines ${M}_{\mathrm{supra}}^{H}$ and ${M}_{\mathrm{threshold}}^{H}$) and ${M}_{\max }^{H}$ due to the difficulty in estimating the critical value of the central density at which the formation of quark matter is triggered. A possible range of values for ${M}_{\max }^{H}$ is ${M}_{\max }^{H}=(1.5-1.6){M}_{\odot }$, and ${M}_{\mathrm{TOV}}^{H}$ takes only slightly higher values.

The value of ${M}_{\mathrm{TOV}}^{Q}$ has a larger uncertainty stemming from the unknown value of the maximum mass of compact stars, which, as mentioned in the Introduction, is estimated to be $(2.2\pm 0.2){M}_{\odot }$. Interestingly, this same range is the one allowed by the microphysics of nucleation of quark matter in hadronic matter. For significantly larger values of ${M}_{\mathrm{TOV}}^{Q}$, the process of nucleation of quark matter would no longer be energetically convenient (Drago et al. 2019).

Finally, the uncertainty on ${M}_{\mathrm{TOV}}^{Q}$ translates into a similar uncertainty on the values of ${M}_{\mathrm{TOV}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$. Once a value for ${M}_{\mathrm{TOV}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$ is chosen, the value of ${M}_{\mathrm{supra}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$ can been explicitly computed, while the value of ${M}_{\mathrm{threshold}}^{\mathrm{Coll} \mbox{-} \mathrm{hyb}}$ has only been estimated by using the results of Weih et al. (2018), indicating that the maximum mass of differentially rotating stars is a factor of ∼1.5 larger than ${M}_{\mathrm{TOV}}$.

Recently, a reanalysis of the observational values of masses and radii within the two-families scenario was performed in Char et al. (2019). It was shown that all existing constraints can be satisfied within the model. Interestingly, the ranges of the masses of HSs and QSs turn out to be in agreement with the distribution of the masses of the two subpopulations discussed in Antoniadis et al. (2016), Alsing et al. (2018), and Tauris et al. (2017). The two distributions are centered at ${m}_{1}=1.39{M}_{\odot }$ with a dispersion ${\sigma }_{1}=0.06{M}_{\odot }$ and at ${m}_{2}=1.81{M}_{\odot }$ with a dispersion ${\sigma }_{2}=0.18{M}_{\odot }$, respectively (Antoniadis et al. 2016). Our values for ${M}_{\max }^{H}$ and ${M}_{\min }^{Q}$ are very close to the upper and lower limits for m₁ and m₂, respectively. It is therefore tempting to interpret the first subpopulation as composed of HSs and the second of QSs.

Within the two-families scenario, there are three possible merger events: HS–HS, HS–QS, and QS–QS. Notice that HS–HS is obviously the only possible type of merger in a one-family scenario and that QS–QS has been discussed in the past in a few papers by various authors without making any assumption on the number of families of compact stars (Haensel et al. 1991; Bauswein et al. 2009, 2010). The HS–QS merger is instead possible only within the two-families scenario. These different possibilities are displayed in Figure 3 in which the two axes correspond to the two parameters determining the main properties of the binary, i.e., the chirp mass M_chirp and the mass asymmetry $q={m}_{2}/{m}_{1}$ (for the sake of discussion, we have made a specific choice of the maximum and minimum masses of HSs and QSs; see the caption). It is clear how rich the phenomenology of mergers can be within the two-families scenario; for instance, within the small window labeled "All-comb.," there could be events of HS–HS, HS–QS, and QS–QS merger. On the other hand, there are regions of the diagram in which we predict only one specific type of merger to be possible (for instance, only HS–HS mergers to the left of the green line). As a consequence, for the same value of M_chirp, there could be a prompt collapse to a BH if the merger is of the HS–HS type, while a hypermassive or supramassive configuration could form if it is an HS–QS merger; see Drago & Pagliara (2018). Note that the classification presented in Figure 3 depends only on the maximum and minimum masses of the HSs and QSs, since it does not suggest the outcome of the mergers (this question will be discussed in Section 4.2) but only indicates which types of mergers are possible within the two-families scheme. Finally, the probabilities of the various merger processes will be addressed in Section 6 by using a population synthesis code.

Zoom In Zoom Out Reset image size

4.1. Quark Deconfinement in the Merger of Compact Stars

4.2. Possible Outcomes of a Merger

We considered binary NS merger processes using numerical relativity in the case of pure HSs and for two EoSs, namely, the SFHo and SFHo-HD EoS, presented in Section 2 and we studied the behavior of equal-mass models, as summarized in Table 1, for both of them.

5.1. Numerical Methods and Initial Data

The numerical methods used to perform these simulations are the same as those of De Pietri et al. (2016, 2018) and Maione et al. (2016, 2017). Here we report the general simulation setup and parameters and refer to those previous articles for more details. In particular, the resolution used in this work is dx = 0.1875 CU = 277 m (see De Pietri et al. 2016 for a discussion of the convergence properties of the code).

The simulations were performed using the Einstein Toolkit (Loffler et al. 2012), an open-source modular code for numerical relativity based on Cactus (Allen et al. 2011). The evolved variables were discretized on a Cartesian grid with six levels of fixed mesh refinement, each using twice the resolution of its parent level. The outermost face of the grid was set at $720{M}_{\odot }$ (1040 km) from the center. We solved the BSSN-NOK formulation of Einstein's equations (Nakamura et al. 1987; Shibata & Nakamura 1995; Baumgarte & Shapiro 1998; Alcubierre et al. 2000, 2003) implemented in the McLachlan module (Brown et al. 2009) and the general relativistic hydrodynamics equations with high-resolution shock-capturing methods implemented in the publicly available module GRHydro (Baiotti et al. 2005; Mösta et al. 2014). In particular, we used a finite-volume algorithm with the HLLE Riemann solver (Harten et al. 1983; Einfeldt 1988) and the WENO reconstruction method (Liu et al. 1994; Jiang & Shu 1996). The combination of the BSSN-NOK formalism for Einstein's equations and the WENO reconstruction method was found in De Pietri et al. (2016) to be the best setup within the Einstein Toolkit, even at low resolutions. For the time evolution, we used the method of lines with fourth-order Runge–Kutta (Runge 1895; Kutta 1901). For numerical reasons, the system is evolved on an external-matter atmosphere set to ${\rho }_{\mathrm{atm}}=6.1\times {10}^{5}$ g cm${}^{-3}$ (as in Lehner et al. 2016) that it is just slightly larger than the one used in Sekiguchi et al. (2015) and Radice et al. (2018a), but it is still sufficiently low to avoid the atmosphere's inertial effects on the ejected mass on a timescale of ≃10 ms after the onset of the merger (Sekiguchi et al. 2015). Initial data were generated with the LORENE code (Gourgoulhon et al. 2001) as irrotational binaries in the conformal thin sandwich approximation. In this work, we analyze different combinations of equal-mass setups for binary systems using the SFHo-HD (Drago et al. 2014b; Burgio et al. 2018) and SFHo EoSs (Steiner et al. 2013). For both models, the initial distance was set to 44.3 km, like in Maione et al. (2016) and Feo et al. (2017). Matter description is performed using a 10 piece piecewise polytropic approximation for an EoS of the type $P=P(\rho ,\epsilon )$ supplemented by an additional thermal component described by ${{\rm{\Gamma }}}_{\mathrm{th}}=1.8$. The parameterization of the EoS uses the following prescriptions:

$$\begin{eqnarray}&&\epsilon ={\epsilon }_{0}(\rho )+{\epsilon }_{\mathrm{th}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&p={p}_{0}(\rho )+({{\rm{\Gamma }}}_{\mathrm{th}}-1)\rho {\epsilon }_{\mathrm{th}},\end{eqnarray} \tag{ 4 }$$

where _th is an arbitrary function of the thermodynamical state that has the property of being zero at T = 0, and ${\epsilon }_{0}(\rho )$ and ${p}_{0}(\rho )$ are the internal energy and pressure at T = 0. In particular, for a piecewise polytropic approximation with N pieces, we can write

$$\begin{eqnarray}&&{p}_{0}(\rho )={K}_{i}{\rho }^{{{\rm{\Gamma }}}_{i}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\epsilon }_{0}(\rho )={\epsilon }_{i}+\displaystyle \frac{{K}_{i}}{{{\rm{\Gamma }}}_{i}-1}{\rho }^{{{\rm{\Gamma }}}_{i}-1},\end{eqnarray} \tag{ 6 }$$

and all of the coefficients are set once the polytropic index ${{\rm{\Gamma }}}_{i}$ (i = 0, ... , $N-1$), the transition density ${\rho }_{i}$ (i = 1, ... , $N-1$), and K₀ are chosen. One should note that here _th plays the role of temperature, and it may be converted (using its interpretation as a perfect fluid thermal component) to a temperature scale as $T=({{\rm{\Gamma }}}_{\mathrm{th}}-1){m}_{b}{\epsilon }_{\mathrm{th}}$, where we assumed a conventional value of ${m}_{b}=940\,\mathrm{MeV}$/c² for the mass of a free baryon.

5.2. Gravitational Waves Extraction

During the simulations, the gravitational wave signal is extracted using the Newman–Penrose scalar Ψ₄ (Newman & Penrose 1962; Baker et al. 2002) using the module WeylScalar4. The scalar is linked to the gravitational wave strain by the following relation, which is valid only at spatial infinity:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{4}={\ddot{h}}_{+}-i{\ddot{h}}_{\times },\end{eqnarray} \tag{ 7 }$$

where ${h}_{+}$ and h_× are the two polarizations of the gravitational wave strain h. The signal is then decomposed into multipoles using spin-weighted spherical harmonics of weight $(-2)$ (Thorne 1980). This procedure is made by the module MULTIPOLE using the following relation:

$$\begin{eqnarray}&&{\psi }_{4}{\left(t,r,\theta ,\phi )=\displaystyle \sum _{l=2}^{\infty }\displaystyle \sum _{m=-l}^{l}{\psi }_{4}^{{lm}}(t,r\right)}_{-2}{Y}_{{lm}}(\theta ,\phi ).\end{eqnarray} \tag{ 8 }$$

In this work, we focus only on the dominant $l=m=2$ mode; therefore, we will refer to ${h}_{\mathrm{2,2}}$ as h for the rest of this paper. In order to extract the gravitational wave strain from Ψ₄ and minimize the errors due to the extraction, one has to extrapolate the signal extracted within the simulation at a finite distance from the source to infinity, in order to satisfy the previous equation. Then the extrapolated Ψ₄ is integrated twice in time, employing an appropriate technique in order to reduce the amplitude oscillations caused by the high-frequency noise aliased in the low-frequency signal and amplified by the integration process (Reisswig & Pollney 2011; Nakano et al. 2015). The procedure adopted in this work is extensively discussed in Maione et al. (2016). First, Ψ₄ is extrapolated to spatial infinity using the second-order perturbative correction of Nakano et al. (2015),

$$\begin{eqnarray}\begin{array}{rcl}R{\psi }_{4}^{{lm}}({t}_{\mathrm{ret}}){| }_{r=\infty } & = & \left(1-\displaystyle \frac{2M}{R}\right)\left(r{\ddot{\bar{h}}}^{{lm}}({t}_{\mathrm{ret}})\right.\\ & & -\,\displaystyle \frac{(l-1)(l+2)}{2R}{\dot{\bar{h}}}^{{lm}}({t}_{\mathrm{ret}})\\ & & +\left.\,\displaystyle \frac{(l-1)(l+2)({l}^{2}+l-4)}{8{R}^{2}}{\bar{h}}^{{lm}}({t}_{\mathrm{ret}})\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where the gravitational wave strains (${\bar{h}}^{{lm}}$) at finite radius (in our case, at a fixed coordinate radius R = 1033 km) are computed by integrating the Newman–Penrose scalar twice in time with a simple trapezoid rule, starting from zero coordinate time, and fixing only the two physically meaningful integration constants Q₀ and Q₁ by subtracting a linear fit of itself from the signal,

$$\begin{eqnarray}&&{\bar{h}}_{{lm}}^{(0)}={\int }_{0}^{t}{dt}^{\prime} {\int }_{0}^{t^{\prime} }{dt}^{\prime\prime} {\psi }_{4}^{{lm}}(t^{\prime\prime} ,r),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\bar{h}}_{{lm}}={\bar{h}}_{{lm}}^{(0)}-{Q}_{1}t-{Q}_{0}.\end{eqnarray} \tag{ 11 }$$

After the integration, we apply a digital high-pass Butterworth filter, which is designed to have a maximum amplitude reduction of 0.01 dB at the initial gravitational wave frequency ${f}_{{t}_{0}}$ (assumed to be two times the initial orbital angular velocity, as reported by the LORENE code) and an amplitude reduction of 80 dB at frequency $0.1{f}_{{t}_{0}}$. All of the gravitational wave–related information is reported in a function of the retarded time,

$$\begin{eqnarray}&&{t}_{\mathrm{ret}}=t-{R}^{* },\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{R}^{* }=R+2{M}_{\mathrm{ADM}}\mathrm{log}\left(\displaystyle \frac{R}{2{M}_{\mathrm{ADM}}}-1\right).\end{eqnarray} \tag{ 13 }$$

The physical observable we are focusing on in this work, which is the gravitational wave amplitude spectral density $| \tilde{h}(f)| {f}^{1/2}$, is calculated from the gravitational wave strain obtained with the aforementioned procedure using

$$\begin{eqnarray}&&| \tilde{h}(f)| =\sqrt{\displaystyle \frac{| {\tilde{h}}_{+}(f){| }^{2}+| {\tilde{h}}_{\times }{| }^{2}}{2}},\end{eqnarray} \tag{ 14 }$$

where $\tilde{h}(f)$ is the Fourier transform of the complex gravitational wave strain,

$$\begin{eqnarray}&&\tilde{h}(f)={\int }_{{t}_{i}}^{{t}_{f}}h(t){e}^{-2\pi {ift}}{dt},\end{eqnarray} \tag{ 15 }$$

and the associated gravitational wave energy in the post-merger phase in the time interval from 1 ms after the start of the simulation up to 15 ms after the merger time. In this respect, we have defined the merger time as the time at which the l = 2, m = 2 gravitational wave strain ${h}^{{lm}}$ has a maximum. All times are reported by fixing t = 0 as the merger time.

5.3. Estimate of the Ejected and Disk Mass

5.4. Results of the Numerical Simulations

In this work, we simulate equal-mass binary systems for HS–HS and NS–NS configurations during about four orbits before merger and for about 20 ms after merger.

All of the simulated models show a very similar behavior. To show the main properties of the evolution, snapshots for model SFHo-HD 118vs118 are shown in Figures 4 and 5. In particular, Figure 4 shows the density of the expelled matter and its temperature at different times, while in Figure 5 its velocity and the localization of unbound mass are shown at the same times. The computed gravitational wave signal of our simulation is shown in Figures 6 and 7, where we display the power spectrum extracted using a detector at R = 1034 km for models described by the SFHo-HD EoS and the SFHo EoS, respectively. Associated with each simulation, we report the values of M_ej and M_disk in Table 2. Their values are also plotted in Figure 12 as a function of ${M}_{\mathrm{tot}}$. A striking feature of the gravitational wave spectrum is that models with ${M}_{\mathrm{tot}}\geqslant 2.48$ and ${M}_{\mathrm{tot}}\geqslant 2.84$ for SFHo-HD and SFHo, respectively, do not show any post-merger f₂ peak. This is related to the fact that, in these cases, we have a direct collapse to a BH (less than 1 ms after the merger). The fact that a prompt collapse occurs is clearly shown in the plot of the maximum density as a function of time for models described by the SFHo-HD (Figure 8) and SFHo (Figure 9) EoSs. In the cases in which a prompt collapse does not occur, the main structure of the spectrum is the same as the one discussed in Takami et al. (2015), Bauswein & Stergioulas (2015), and Maione et al. (2017), and references therein. Importantly, the post-merger spectrum is characterized by a main f₂ peak whose frequency depends on the EoS. Moreover, its frequency increases with the mass of the star, while its tidal parameter decreases. One should also notice that secondary peaks are present, but their relative importance decreases as the total mass of the system is reduced.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Values of Ejected Mass (M_ej), Disk Mass (M_disk), Energy Emitted in Gravitational Waves in the Post-merger Phase (${E}_{\mathrm{gw}}^{\mathrm{POST}}$), Frequency of the Main Peak of the Gravitational Wave Signal (f₂), and Time to Collapse to BH (t_BH) for All of the Models Considered in This Work

Model	Mej	Mdisk	${E}_{\mathrm{gw}}^{\mathrm{POST}}$	f2	tBH
	(mM⊙)	(mM⊙)	(mM⊙)	(kHz)	(ms)
SFHo-HD 118vs118	12.993	12.92	25.42	3.71	3.82
SFHo-HD 120vs120	9.435	13.81	22.42	4.00	3.16
SFHo-HD 122vs122	4.290	8.34	6.06	⋯	1.91
SFHo-HD 124vs124	3.011	2.89	0.66	⋯	1.00
SFHo-HD 126vs126	0.737	2.45	0.20	⋯	0.79
SFHo-HD 128vs128	0.055	0.74	0.04	⋯	0.70
SFHo-HD 130vs130	0.043	0.71	0.01	⋯	0.59
SFHo 118vs118	1.968	76.66	42.16	2.88	⋯
SFHo 120vs120	2.085	71.72	43.87	2.90	⋯
SFHo 122vs122	1.730	91.81	42.00	2.90	⋯
SFHo 124vs124	1.824	65.58	52.98	2.96	⋯
SFHo 126vs126	2.375	60.86	58.33	2.98	⋯
SFHo 128vs128	3.145	112.24	50.33	3.05	⋯
SFHo 130vs130	4.523	73.82	59.33	3.06	⋯
SFHo 132vs132	6.007	88.87	67.29	3.18	25.75
SFHo 134vs134	9.511	49.27	65.09	3.25	13.55
SFHo 136vs136	16.244	30.71	58.76	3.40	9.42
SFHo 138vs138	10.367	16.09	46.06	3.55	5.06
SFHo 140vs140	4.170	6.45	22.39	⋯	2.13
SFHo 142vs142	2.247	2.01	2.02	⋯	0.98

Download table as:  ASCIITypeset image

The presence of a well-defined main peak and its frequency f₂ could lead to a clear determination of the kind of merger and the properties of the EoS. In fact, that would be a precise signature that the system did not directly collapse to a BH and that its spectrum encodes information on the properties (mainly the stiffness) of dense matter. For example, the merger of two compact stars of mass 1.18 M_⊙ shows a main peak at frequency ${f}_{2}=3.71$ or $2.88\,{\rm{k}}{\rm{H}}{\rm{z}}$ if the EoS describing the matter is SFHo-HD or SFHo, respectively. It is interesting to remark that the value of this frequency is strictly related to the value of ${R}_{1.6}$ as discussed in Bauswein et al. (2016). Therefore, our result does not depend on the details of the EoS but on the request of having very small radii in the hadronic branch. Unfortunately, the amount of gravitational wave energy available for the mode detection (see Table 2) is going to be at best of the order of 0.1 M_⊙; therefore, it is unlikely that such a peak will be detected by the present generation of gravitational wave detectors, but they will be a distinct feature to be observed in third-generation detectors.

Let us now discuss the results for the ejected mass and the amount of mass that is left in the disk. In Figure 10 we show an example of the outflow for SFHo-HD 118vs118. By following the flow of the unbound matter during the evolution, we can estimate the ejected mass by integrating (at each radius) on time the unbound matter flow. In this model, the total flows of unbound mass at coordinate radii of 96, 220, 443, and 590 km are (5.0, 10.7, 11.24, 11.45) ${\rm{m}}{M}_{\odot }$, respectively. Indeed, this procedure for computing the total ejected mass is very sensitive to the extraction radius, and this is a common property of all of the simulated models. In particular, we compared the results obtained by using this method with those obtained by considering the maximum of total unbound mass present on the whole computational domain at a given time. We note that the difference between these two estimates of M_ej is always less than about 20%, and the two results get closer by increasing the extraction radius. From the lower panel of Figure 10 one should note that most of the dynamically unbound mass is generated between 100 and 200 km from the center of the remnant, and it is related to three main outbursts of ejected matter associated with the three bouncings of the maximum density. This is more clear looking at Figure 11 where we display the radial density profile (as a function of the coordinate radius) of the unbound matter at different times. There, one should note that three shocks are created within a 100 km radius from the center of the star (t = 0.14, 1.15, and 2.15 ms), the first being of lower amplitude. Then they move out, amplify, and spread while the system evolves. This is analogous to what is observed for the flux of matter at different radii (see Figure 10, upper panel), where one can see that at R = 100 km, one has three clearly separated peaks that join and spread as the unbound matter flows outside the computational domain.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The same overall dynamics is associated with the mass ejection for all of the considered models (except models that feature a direct collapse to a BH and for which the mass ejection is suppressed). All of the results for M_ej and M_disk as a function of ${M}_{\mathrm{tot}}$ are summarized in the lower and upper panels, respectively, of Figure 12. The explicit values are also reported in Table 2. A striking feature concerning the mass ejected is the presence of a maximum, which is located at a value of ${M}_{\mathrm{tot}}$ slightly smaller than the value of M_threshold. This is one of the main results of the present work. In particular, this maximum is located at $2.72{M}_{\odot }$ for the SFHo EoS, and it corresponds to an ejected mass of ≃16 mM_⊙. For the SFHo-HD case, we have not determined the maximum of the ejected mass, which should be located below or close to $\sim 2.36{M}_{\odot }$ (the simulation with the lowest ${M}_{\mathrm{tot}}$). For this specific value, we have found an ejected mass of ≃13 mM_⊙, very close to the maximum for the SFHo EoS. The plots of the ejected mass show a very steep decay of its value as the total mass of the binary increases, and the ejection is almost completely suppressed in the case of a prompt collapse to a BH.

Zoom In Zoom Out Reset image size

Another clear difference between the two EoSs concerns the value of M_disk: in the case of SFHo-HD, it is always $\lesssim 0.01{M}_{\odot }$, whereas for SFHo, it can be an order of magnitude larger. These are potentially very interesting results for the phenomenology of the kilonova. As we will discuss in Section 8, the luminosity of the different components of the kilonova (red, blue, and purple kilanova; Perego et al. 2017) strongly depends on the specific mechanisms for the ejection of mass during the different stages of the merger. A relevant fraction of the ejected mass can, in principle, come from the disk; therefore, a very small value of M_disk corresponds to a strong suppression of certain components of the kilonova signal.

Let us now compare our results with the ones obtained by similar numerical simulations present in the literature. In particular, there are a few simulations using the SFHo EoS for masses similar to our model SFHo 136vs136. The outflow dynamics for this model is shown in Figure 13. The results we obtain are in general agreement with the results obtained in Sekiguchi et al. (2015, 2016) for the SFHo EoS. In those works, it was found that the associated ejected mass is of the order of 10 mM_⊙ for the mass of the two stars of 1.35 M_⊙ and of the same order of magnitude for unequal-mass systems of approximately the same total mass (Shibata et al. 2017). In the same work, it was also found that the expected disk mass (for the same systems) is of the order of 50–120 mM_⊙ and indeed of the same order of magnitude as the values found in the present work. Computations of the ejected mass for the SFHo EoS were also presented in Bauswein et al. (2013b) using SPH dynamics, Lehner et al. (2016), and, using the WhiskyTHC code, Bovard et al. (2017) and Radice et al. (2018a), where a comparison among all of these numerical results is shown in their Table 3. One can notice that the variability among the various estimates is rather large. Concerning our work, in particular, a few potentially relevant effects have not been incorporated, such as a treatment of the neutrino transport and a fully consistent description of the thermal component of the EoS, and we use a piecewise polytropic approximation. However, our results are in very good agreement with Sekiguchi et al. (2015), where full beta equilibrium, thermal evolution, and approximate neutrino cooling and absorption were considered.

Zoom In Zoom Out Reset image size

Another important effect that one should consider is related to the dependence of M_ej and M_disk on the mass asymmetry q. This issue has been discussed in Rezzolla et al. (2010), Giacomazzo et al. (2013), and, more recently, Kiuchi et al. (2019). It is expected that ${M}_{\mathrm{disk}}$ is always larger for unequal-mass binaries because of the more efficient tidal interactions and angular momentum transfer during the merger. The same consideration seems not to apply in general for M_ej, and further studies are needed, in particular when considering models close to the threshold mass, because the state after the merger may rapidly change from producing a direct collapse to forming a short-lived remnant (Kiuchi et al. 2019).

5.5. Estimate of the Threshold Mass

5.6. Trigger of the Phase Conversion

For ${M}_{\mathrm{tot}}\lt {M}_{\mathrm{threshold}}^{H}$, the two-families scenario predicts that the phase conversion of hadronic matter to quark matter necessarily occurs. Simulating the process of conversion itself after the merger is numerically very challenging (one needs to couple the code used, for instance, in Herzog & Ropke 2011 and Pagliara et al. 2013 to Cactus). We can, however, estimate when the conversion should start. As explained previously, the conversion is triggered only when a significant amount of strangeness is produced during the evolution of the remnant, i.e., ${Y}_{S}\gtrsim 0.2$. In Figure 8 (left panel), we display the stripe in the mass density–temperature diagram for which the value of Y_S is large enough for the conversion to start. In the right panel, we show the temporal evolution of the maximum baryon density for all of the runs of the HS–HS merger. Let us first discuss the case of the run SFHo-HD 118vs118. One can see that after the first two oscillations of the remnant, the condition for quark matter nucleation is fulfilled; see also Figures 14 and 15, which allow one to visualize the regions inside the remnant where the density and temperature reach the condition for the beginning of the conversion process. The subsequent temporal evolution of the remnant would be a dramatic change of the structure of the star: within a few ms, a big part of the star is converted into quark matter, the radius of the star increases by a few km, and the heat released by the conversion heats up the star. Also, a significant amount of differential rotation can develop due to the very fast change of the total moment of inertia of the star (Pili et al. 2016). Finally, the stiffening of the EoS would stabilize the remnant with respect to the collapse. We have a similar evolution for the SFHo-HD 12 model, with the only difference being that the conditions for the phase transition are met during the second bounce instead of the third one.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A comment concerning the results for the mass ejected in these cases is in order. The oscillations of the maximum density are associated with the shock waves, which, in turn, are responsible for most of the dynamically ejected mass. Thus, even if we are not implementing the conversion to quark matter (which could modify the temporal behavior of the maximum density) in our simulations, we can assume that the results of the mass ejected by HS–HS mergers before the trigger of the phase transition (before the third peak in the SFHo-HD 118vs118 model) are fairly reliable.

Let us now discuss the cases with larger masses. For the model SFHo-HD 124vs124, which leads to a prompt collapse, the nucleation condition is met only when the collapse has already started and the formation of the quark phase cannot halt it. The intermediate case of SFHo-HD 122vs122 instead shows a first oscillation during which the maximum density reaches the threshold for nucleation, but the second oscillation already leads to the collapse of the star. In this case, it is possible that a significant part of the star is converted to quark matter, but it is quite likely that the consequent stiffening of the EoS is not strong enough to prevent the collapse.

The crucial information that we now want to provide is an estimate of the rate of events that we expect for the different types of mergers. For this purpose, we studied the evolution toward mergers of double compact objects in the two-families scenario. Specifically, we used the Startrack population synthesis code (Belczynski et al. 2002, 2008) with further updates described in Wiktorowicz et al. (2017, and references therein). As GW170817 is our main point of interest, we concentrated on systems having parameters similar to the observed ones and adapted the simulation parameters to reproduce an environment similar to that of the host galaxy (NGC 4993).

Recently, Belczynski et al. (2018) performed a study of merger rates of two NSs in the NGC 4993 environment. Here we also include QSs as possible components of the progenitor binary. Specifically, following Belczynski et al. (2018), we assumed the metallicity of the environment in which GW170817 was formed to be Z = 0.01, which corresponds to about 50% of the solar metallicity. The results were scaled to the observable volume of a LIGO detector for events similar to GW170817, which may be approximately represented by a sphere with a radius of $D\approx 100$ Mpc. On the basis of an Illustris simulation (Vogelsberger et al. 2014), this volume contains about $1.1\times {10}^{14}{M}_{\odot }$ of stellar mass in elliptical galaxies. The star formation history was chosen to be uniform between 3 and 7 Gyr ago (Troja et al. 2017), which is different from Belczynski et al. (2018), who used a burst-like star formation history occurring 1, 5, or 10 Gyr ago.

As far as the binary evolution is concerned, all mass transfer was assumed to be conservative; i.e., all mass lost by the donor is transferred to the accretor. Such an approach increases the rates for double compact object mergers (Chruslinska et al. 2018). On the other hand, we kept the Maxwellian distribution of natal kicks to $\sigma =265$ km s⁻¹, which conforms to the observational estimates from Galactic pulsar proper motions (Hobbs et al. 2005), although lower natal kicks would increase the rates (Chruslinska et al. 2018).

Similar to Wiktorowicz et al. (2017), we adopted the two-families scenario by assuming that an HS converts into a QS when its gravitational mass reaches a maximal value, ${M}_{\max ,\mathrm{ns}}={M}_{\max }^{H}$. During the transition, the gravitational mass of the compact object changes instantly (compared to the typical timescales of population synthesis analysis), and we update the orbital parameters accordingly.

Calculated merger rates are presented in Table 3. Two models adopt the two-families scenario with deconfinement taking place at HS masses of 1.5 and $1.6{M}_{\odot }$ (models ${M}_{max}^{H}=1.5$ and $1.6\,{M}_{\odot }$, respectively). In the third model (one-family), which we provide for reference, all compact objects are HSs, and the deconfinement process never occurs. We note that in all models, it was assumed that all compact objects with masses above $M\gt 2.5{M}_{\odot }$ formed prior to the merger, collapse to a BH, and are not included in our study.

Table 3.  Merger Rates from Population Synthesis $[\times {10}^{-3}\,{\mathrm{yr}}^{-1}$]

	All Mergers			GW170817-like
				$0.7\lt q\lt 1.0$			$0.7\lt q\lt 0.85$
Model	HS–HS	HS–QS	QS–QS	HS–HS	HS–QS	QS–QS	HS–HS	HS–QS	QS–QS
${M}_{\max }^{H}=1.5\,{M}_{\odot }$	9.1	3.1	0.2	6.4	0.4	0.01	0.03	0.2	⋯
${M}_{\max }^{H}=1.6\,{M}_{\odot }$	9.2	3.2	0.02	6.5	0.3	⋯	0.1	0.2	⋯
One-family	12.8	⋯	⋯	6.6	⋯	⋯	0.3	⋯	⋯

Note. For the two-families scenario, we have adopted two possible values of ${M}_{\max }^{{\rm{H}}}$, and we report the results for all of the merger combinations: HS–HS, HS–QS, and QS–QS mergers. Results for the one-family scenario are also reported for comparison. Notice that for values of q ≲ 0.85, an HS–QS merger is more probable than an HS–HS merger.

Download table as:  ASCIITypeset image

In our analysis, we compare merger rates for all systems and those corresponding to events similar to the GW170817 event. The latter were chosen as events with chirp masses in the range ${M}_{\mathrm{chirp}}=(1.188\pm 0.1){M}_{\odot }$ (the $\pm 0.1{M}_{\odot }$ was chosen arbitrarily, but the specific adopted value of the range has a negligible effect on the conclusions) and a mass ratio of $q\gt 0.7$ ($q={m}_{2}/{m}_{1}\gt 0.7$, where ${m}_{2}\leqslant {m}_{1}$).

First, when considering all mergers, it is noticeable that the rate of HS–QS is not strongly suppressed in respect to the rate of HS–HS; they constitute roughly one-fourth of the total rate.

Second, concerning the GW170817-like events, different types of mergers are allowed depending on the value of the mass ratio q. In Table 3 we compare the merger rates for the entire range of mass ratios ($q\gt 0.7$) and, separately, for lower mass ratios ($0.7\lt q\lt 0.85$). The latter was motivated by an analysis of the GW170817 event (Abbott et al. 2018) indicating that a mass ratio $q\sim 0.85$ is highly compatible with the data. Results for both ranges of mass ratios are significantly different. The entire range ($q\gt 0.7$) is dominated by high mass ratio systems that are mostly HS–HS, because for $q\approx 1$ and ${M}_{\mathrm{chirp}}\approx 1.188$, the components' masses can barely enter the QS mass range ${M}_{min}^{Q}\gt 1.37$ or $1.46{M}_{\odot }$ (for ${M}_{max}^{H}=1.5$ and $1.6\,{M}_{\odot }$, respectively). On the other hand, for low mass ratios, one of the stars usually becomes more massive than ${M}_{\min }^{Q}$; thus, according to the two-families scenarios, it converts into a QS. As a result, although the mergers within the entire mass ratio range ($q\gt 0.7$) show a preference for HS–HS mergers, the ones with a lower mass ratio ($0.7\lt q\lt 0.85$) are typically HS–QS mergers. However, the former give much higher merger rates (∼6.8 kyr⁻¹) than the latter (∼0.2–0.3 kyr⁻¹).

It is also important to note that the rates of QS–QS mergers are in general very much suppressed, reaching at most a few percent of the total expected events involving at least one QS in the merger, a result in agreement with the findings of Wiktorowicz et al. (2017). In particular, the probability that GW170817 was due to a QS–QS merger is totally negligible.

Finally, we note that all of the merger rates obtained from population studies are significantly below the one estimated from the GW170817 event (by more than 2 orders of magnitude), as previously shown by Belczynski et al. (2018). It may be a result of, e.g., low observational statistics (just one event), poorly understood binary evolution phases (e.g., common envelope survival of low mass ratio binaries), or the existence of an unknown evolutionary scenario.

In the following, we discuss the interpretations of GW170817 in terms of the different types of mergers and their probability.

6.1. HS–HS Binary

6.2. HS–QS Binary

6.3. QS–QS Binary

In the case of GW170817, the total mass is ${M}_{\mathrm{tot}}^{170817}={2.74}_{-0.01}^{+0.04}{M}_{\odot }$, and, as explained before, within the two-families scenario, this event can be interpreted only as due to the merger of an HS with a QS. Indeed, ${M}_{\mathrm{threshold}}^{H}\lt {M}_{\mathrm{tot}}^{170817}$, and therefore an HS–HS binary would promptly collapse to a BH. However, a direct collapse is excluded by the observation of an sGRB ∼2 s after the merger. Also, as discussed previously, it cannot be interpreted as a merger of two QSs.

A first question concerns the probability of detecting such an HS–QS merger event. From the binary evolution point of view, the population synthesis analysis has shown that the delay time, which is ∼4.6 Gyr, is consistent with the lack of recent star formation in the host galaxy.

Concerning the phenomenology of GW170817, without detailed numerical simulations of the merger of an HS with a QS, it is difficult to firmly establish whether we can explain all the features of the observed signals in our scenario. We can, however, verify whether our model satisfies the upper and lower limits on $\tilde{{\rm{\Lambda }}}$. The lower limit on $\tilde{{\rm{\Lambda }}}$ in particular, empirically found in Radice et al. (2018b), stems from the request for a description of the kilonova signal AT 2017gfo (this constraint was recently criticized in Kiuchi et al. 2019).

In Figure 16, we show the tidal deformabilities of the two components of the binary and $\tilde{{\rm{\Lambda }}}$ for the chirp mass of GW170817 for the one-family and two-families scenarios (for the cases of HS–HS, HS–QS, and QS–QS mergers¹⁴ ). Both the upper limit on $\tilde{{\rm{\Lambda }}}$, as obtained by the gravitational wave's signal, and the lower limit, as inferred from the kilonova analysis (Radice et al. 2018b), can be fulfilled within the two-families scenario in the case of an HS–QS system (see also Burgio et al. 2018). Notice that the small value of Λ for the HS component is due mainly to the softening associated with the appearance of delta resonances.

Zoom In Zoom Out Reset image size

We observe that the case of HS–HS is not compatible with the constraints on $\tilde{{\rm{\Lambda }}}$ because it provides a too-small value for it (Radice et al. 2018b; Kiuchi et al. 2019). Thus, also from the constraints on $\tilde{{\rm{\Lambda }}}$, we can rule out the possibility that GW170817 was due to an HS–HS merger. We need to discuss now how the electromagnetic counterparts of GW170817 would be produced in our scheme.

Concerning the sGRB, we can assume that it has been produced by the disk accreting onto the BH in agreement with the classification discussed in Section 4. Indeed, for GRB 170817A, there is no clear indication of the existence of an extended emission after the prompt signal. Concerning the kilonova, a crucial issue is connected with the fate of the quark matter ejected from the QS during the merger. There are a few studies indicating that quark matter evaporates into nucleons and that this process can be so efficient that only strangelets with baryon numbers much larger than 10⁴⁰ can survive (Alcock & Farhi 1985; Madsen et al. 1986; N. Bucciantini et al. 2019, in preparation). Most of the matter ejected from the QS will therefore rapidly evaporate into nucleons and can contribute to the kilonova signal. If this idea is correct, the kilonova produced after an HS–QS merger would not show significant differences, concerning the type of matter ejected, with respect to a kilonova produced after an HS–HS merger. Actually, the process of matter ejection can be very efficient, and it opens up the possibility of explaining AT 2017gfo as an HS–QS merger: in the case of AT 2017gfo, the radii of the two compact objects are both rather small, the system is asymmetric, and the threshold mass is large, as discussed in Section 3. These are exactly the requests posed in Kiuchi et al. (2019), and, as discussed in that paper, there are no examples of EoSs based on a single-family scenario able to satisfy all of those requests.

Finally, we note that while the scenario of a delayed collapse (a few tens of ms after the merger) is adopted by most of the literature on GW170817, there are some studies indicating that the merger remnant could be a supramassive or even a stable star (Ai et al. 2018; Li et al. 2018b) and that it could have been active for a much longer timescale. In van Putten & Della Valle (2018), a post-merger gravitational wave signal has been reported with an associated gravitational wave energy lower than the sensitivity estimates of the LIGO-VIRGO collaboration. This signal would have lasted a few seconds, with an initial frequency of about 700 Hz. Another possible signal of a long-lived remnant comes from the analysis of Piro et al. (2019), in which was found a low-significance temporal feature in the X-ray spectrum ∼150 days after the merger that is consistent with a sudden reactivation of the central compact star. If these analyses are correct, they can have very important implications for the EoS of dense matter that we will address in a forthcoming study.

We have explored the phenomenological consequences of the existence of two families of compact stars on the observations of mergers. A specific feature of the two-families scenario is the possibility that at fixed values of the chirp mass and mass asymmetry, the merging binary can be composed by two HSs, two QSs, or an HS and a QS. Each of these possibilities has its own specific signatures.

Let us first consider systems similar to the one associated with GW170817. Our population synthesis analysis basically excludes the possibility that such binary systems are made of two QSs (notice that this result depends only on the adopted value for ${M}_{\max }^{H}$, which is rather tightly determined in our scenario). They can be formed either by two HSs or by an HS and a QS.

If we look to the whole range of values for the mass ratio q of the two stars, the rate of HS–HS merger events is a factor of (16–20) larger than the rate of HS–QS mergers, while for a small mass ratio ($0.7\leqslant q\leqslant 0.85$), HS–QS systems are about two to six times more probable than HS–HS systems. While in the HS–HS case, we would predict a prompt collapse (no sGRB and a very faint kilonova), in the HS–QS case, a hypermassive remnant is more likely to form (with associated electromagnetic signals similar to GRB 170817A and AT 2017gfo). Notice that within the one-family scenario, the total mass associated with GW170817 can be regarded as a lower limit for the threshold mass, since there are clear indications that a direct collapse to a BH did not take place. In the two-families scenario, the threshold mass instead depends on the type of merger, and it is dramatically different for an HS–HS and an HS–QS merger. In particular, for the HS–HS case, the threshold mass, as found in Section 5.5, is of the order of $2.5{M}_{\odot }$ and thus significantly smaller than the mass associated with GW170817.

The smoking gun of our scenario would be just a single detection of a source of gravitational waves with a total mass smaller than the one of GW170817 but lacking a significant electromagnetic counterpart (indeed, it would be interpreted as an HS–HS merger). This possible signature is valid for values of the total mass in the range $(2.48\mbox{--}2.74){M}_{\odot }$, i.e., between the threshold mass for the merger of two HSs and the mass of the system associated with GW170817.

Another promising way to test the two-families scenario is related to the observation of binaries with low values of ${M}_{\mathrm{tot}}$. Let us concentrate on ${M}_{\mathrm{tot}}=2.4{M}_{\odot }$ with two stars of $1.2{M}_{\odot }$. A first difference between the two-families and the one-family scenarios concerns the value of $\tilde{{\rm{\Lambda }}}$, which can be significantly smaller (almost a factor of 2 for our specific choice of the hadronic EoS) for a merger of HS–HS with respect to the merger of two NSs (see the values for SFHo-HD and SFHo in Table 1). This is due to the appearance of delta resonances in hadronic matter at densities of about twice the saturation density. Accurate future measurements of the gravitational wave's signal associated with the inspiral phase should put interesting upper limits on $\tilde{{\rm{\Lambda }}}$ and therefore test the possibility of the merger of very compact stellar objects, i.e., HSs in our scheme.

Additionally, if observations of the post-merger signal will be feasible in future experiments (most likely the third-generation detectors), one can test the two-families scenario through the measurement of the frequency of the f₂ mode, which, for an HS–HS merger, is of the order of 1 kHz larger than the case of an NS–NS merger during the first milliseconds of the life of the remnant. The subsequent formation of quark matter, which in our scenario corresponds to a stiffening of the EoS, should then decrease the value of f₂ (see Bauswein et al. 2016 for a preliminary analysis of this signature). It is interesting to note that in recent studies discussing a phase transition to quark matter in the post-merger, the appearance of quark matter reduces the lifetime of the remnant and leads to a shift of f₂ to frequencies larger than those of the "progenitor" made only of nucleonic matter (Bauswein et al. 2019; Most et al. 2019).

Let us discuss the possible signatures associated with the kilonova. Its signal depends strongly on the amount of mass dynamically ejected by the merger and in the disk. One can notice from Table 2 that for an HS–HS merger, both of these masses are of the order of $0.01{M}_{\odot }$. Instead, in the case of an NS–NS merger (in the one-family scenario), the dynamically ejected mass is significantly smaller than the mass of the disk. Qualitatively, we would then expect a significant difference between the resulting kilonova signals: within the two-families scenario (in the case of an HS–HS merger), the intermediate-opacity purple component is strongly suppressed with respect to the components associated with the dynamically ejected matter (Perego et al. 2017).

Finally, it is interesting to notice that the suggestion of the two-families scenario originated from the requirement of explaining the observations of stars having large masses and the possible existence of stars with small radii. The same problem presents itself again when trying to model the kilonova AT 2017gfo: as suggested in Siegel (2019), the presence of a strong shock-heated component in the ejecta would favor radii smaller than (11–12) km. Such small radii can be easily accommodated in the two-families scenario but are basically ruled out within the one-family scenario.

We would like to thank the thousands of volunteers who, by their participation in the Universe@Home project,¹⁵ performed the calculations necessary for the population synthesis studies. We benefited from the availability of public software that enabled us to conduct the 3D general relativistic simulations, namely, "LORENE" and the "Einstein Toolkit." We express our gratitude to the many people that contributed to their realization. This work would not have been possible without the CINECA-INFN agreement that provides access to resources on GALILEO and MARCONI at CINECA. We acknowledge PRACE for awarding us access to MARCONI at CINECA, Italy, under grant Pra14_3593. G.W. was partly supported by the Presidents International Fellowship Initiative (PIFI) of the Chinese Academy of Sciences under grant No. 2018PM0017 and by the Strategic Priority Research Program of the Chinese Academy of Science Multi-waveband Gravitational Wave Universe (grant No. XDB23040000). Support from INFN "Iniziativa Specifica NEUMATT" is kindly acknowledged.

In this Appendix, we will summarize the results of the papers published during the last 5 yr leading to a precise definition of the two-families scenario. The aim of this section is to justify the selection of the two EoSs we use for the hadronic and quark families and to explain why we consider those equations as representative of the entire class of EoSs that can be used in that scenario. One important point is that our approach is mainly phenomenological, based on the attempt to interpret the existing data (although controversial) and provide predictions directly related to the interpretation of those data.

The main motivation for proposing the two-families scenario in 2014 (Drago et al. 2014a) was the need to reconcile the existence of very massive stars and the possible existence of very compact objects. In particular, several studies on low-mass X-ray binaries have suggested values for the radii significantly smaller than ∼12 km (see Figure 17 where we display the observational constraints indicated by a couple of very recent analyses; for other examples, see Drago et al. 2016b). While the values of the masses of those objects are not known (but masses smaller than about $1.5{M}_{\odot }$ are favored; d'Etivaux et al. 2019), it is unlikely that all of those objects have masses much larger than the canonical value of $1.4{M}_{\odot }$. Therefore, these very recent analyses (as many others in the recent past) suggest that ${R}_{1.4}$ could be significantly smaller than 12 km. If we assume that those estimates for the radii are correct and interpret those stars as HSs, we need to tune the parameters of our hadronic model in order to explain such small radii. Typically, small radii imply large central densities and make the opening of softening channels, such as the production of delta resonances (Schürhoff et al. 2010) and/or hyperons (Chatterjee & Vidaña 2016), very likely. It is therefore very simple and direct to relate small radii to the appearance of new degrees of freedom: since, in the two-families scenario, the hadronic branch need not support very massive stars, one need not introduce ad hoc repulsive mechanisms that from one side would avoid the production of, e.g., hyperons but on the other side would not allow for small values of ${R}_{1.4}$ (an example of that mechanism is discussed in, e.g., Lonardoni et al. 2015).

Zoom In Zoom Out Reset image size

An important point of our scenario is that the nucleation of quark matter takes place only when a sizable amount of strangeness is present in the HS. An early appearance of hyperons would therefore imply an early transition to quark matter, implying that stars having a mass larger than about $1.4{M}_{\odot }$ would convert to QSs. This scenario would most likely be ruled out by the phenomenology associated with, for instance, magnetar X-bursts. The early formation of delta resonances solves this problem because it shifts the formation of hyperons to larger densities, as found in Drago et al. (2014b) and confirmed in, e.g., Li et al. (2018a) and Li & Sedrakian (2019). Another issue associated with the late appearance of hyperons concerns the cooling: as was shown in Prakash et al. (1992), if, for a certain stellar configuration, the direct nucleon URCA mechanism is forbidden, then the direct URCA process for ${{\rm{\Delta }}}^{-}$ is also forbidden. Thus, we expect stars containing only nucleons and Δs to behave as slow coolers. On the other hand, stars close to the end of the hadronic branch (with central densities above $\sim 5{n}_{0}$) contain enough Λs that the direct URCA is active and would instead cool rapidly. Interestingly, the phenomenology of the cooling of compact stars indicates that, at least for some cases, fast cooling is needed to explain the data (Brown et al. 2018).

Let us first discuss the delta resonances. The couplings between them and mesons (within the relativistic mean field approach that we use here) because of SU(2) symmetry must be of the same order of magnitude as the couplings between nucleons and mesons. In particular, in Drago et al. (2014b), it was realized that the properties of the EoS are determined by the relative magnitude between the couplings of deltas and scalar mesons and of deltas and vector mesons. Thus, we can fix the couplings with vector mesons to be equal to the couplings between nucleons and vector mesons and regard the coupling with scalar mesons x_σΔ as a free parameter. Actually, electron–nucleus scattering data constrain x_σΔ to vary in the range 1–1.15. With the corresponding EoSs, we obtain the dashed mass–radius curves shown in Figure 17 for the extremes of the allowed range for x_σΔ.

Let us now discuss the effect of the appearance of hyperons. As in the case of delta resonances, their effects on the EoS in our framework are determined by the couplings with mesons. Moreover, one should extend the discussion to SU(3) and consider the full scalar and vector mesons nonets. There are clearly too many free parameters and only a few experimental data: one can use SU(3) symmetry as a guideline (as assumed here), but violations are possible; see Weissenborn et al. (2012). Interestingly, the value of the threshold mass does not strongly depend on the exact density at which the hyperons appear. This can be tested by artificially switching off the hyperons when computing the EoS. By fixing ${x}_{\sigma {\rm{\Delta }}}=1.15$ as in the numerical simulations and using the fitting formula of Köppel et al. (2019), we obtain ${M}_{{\rm{t}}{\rm{h}}{\rm{r}}{\rm{e}}{\rm{s}}{\rm{h}}{\rm{o}}{\rm{l}}{\rm{d}}}=2.52$ and $2.5{M}_{\odot }$ for the cases with and without hyperons, respectively. Instead, by setting ${x}_{\sigma {\rm{\Delta }}}=1$ (see black lines in Figure 17), we obtain ${M}_{{\rm{t}}{\rm{h}}{\rm{r}}{\rm{e}}{\rm{s}}{\rm{h}}{\rm{o}}{\rm{l}}{\rm{d}}}=2.6$ and $2.68{M}_{\odot }$ for the cases with deltas and hyperons and only deltas, respectively. Those values for the threshold mass are clearly larger with respect to the case of SFHo-HD but still smaller than the threshold mass predicted by the SFHo EoS and the mass associated with GW170817.

In this respect, our prediction for the cases of prompt collapses would be only mildly dependent on the threshold for the appearance of hyperons, and it would be mainly determined by the effect of introducing the delta resonances in the hadronic EoS (with a specific value of the coupling).

The formation of hyperons is crucial to fix the value of ${M}_{\max }^{H}$, which cannot be very small, as discussed above, but also cannot be too large. For instance, a value of ${M}_{\max }^{H}\sim 1.8{M}_{\odot }$ would lead to a minimum mass of QSs, ${M}_{\min }^{Q}\sim {M}_{\max }^{H}\,-0.1{M}_{\odot }\sim 1.7{M}_{\odot }$. Therefore, only stars more massive than $1.7{M}_{\odot }$ could have large radii.

On the other hand, there are some indications of stars with a mass close to $1.5{M}_{\odot }$ and a large radius: of the order of 13 km for PSR J0437–4715 (see Gonzalez-Caniulef et al. 2019). This stellar object would be interpreted in our scenario as a QS.

To summarize, the range of masses in which HSs and QSs coexist should be centered around 1.5–1.6 M_⊙ in order to interpret the phenomenology.