General Relativistic Simulations of High-Mass Binary Neutron Star Mergers: rapid formation of low-mass stellar black holes

Almost a hundred compact binary mergers have been detected via gravitational waves by the LIGO-Virgo-KAGRA collaboration in the past few years providing us with a significant amount of new information on black holes and neutron stars. In addition to observations, numerical simulations using newly developed modern codes in the field of gravitational wave physics will guide us to understand the nature of single and binary degenerate systems and highly energetic astrophysical processes. We here presented a set of new fully general relativistic hydrodynamic simulations of high-mass binary neutron star systems using the open-source Einstein Toolkit and LORENE codes. We considered systems with total baryonic masses ranging from 2 . 8 𝑀 ⊙ to 4 . 0 𝑀 ⊙ and used the SLy equation of state. We analyzed the gravitational wave signal for all models and reported potential indicators of systems undergoing rapid collapse into a black hole that could be observed by future detectors like the Einstein Telescope and the Cosmic Explorer. The properties of the post-merger black hole, the disk and ejecta masses, and their dependence on the binary parameters were also extracted. We also compared our numerical results with recent analytical fits presented in the literature and provided parameter-dependent semi-analytical relations between the total mass and mass ratio of the systems and the resulting black hole masses and spins


INTRODUCTION
Since the first binary black hole merger observations detected by the LIGO-Virgo collaboration (Abbott et al. 2016), a new window has been opened in astrophysics in the gravitational wave (GW) field (Barack et al. 2019).In 2017, for the first time, the GW signal from a binary neutron star (BNS) merger, named GW170817, was also detected (Abbott et al. 2017a) and it was accompanied by electromagnetic (EM) counterpart observations from gamma rays to radio (Abbott et al. 2017b).Thus, the multi-messenger era has begun.Two years after the first detection of a BNS merger, GW190425, a new BNS system, was detected by the LIGO-Virgo collaboration (Abbott & et al. 2020).For the latter, there is however no detected EM counterpart.GW190425 observation is significant due to being the heaviest BNS system ever observed, with a total mass (3.4 +0.3  −0.1  ⊙ assumption of high-spin prior) much larger than the one measured for the galactic BNS systems to date (Farrow et al. 2019;Zhang et al. 2019).
The possible outcomes of a BNS merger can be either a prompt collapse to a black hole (BH), or the formation of a short-lived ★ E-mail: kutay.arinc.cokluk@ege.edu.trhypermassive neutron star (HMNS) or long-lived supermassive neutron star (SMNS) that eventually collapses to a BH, or a stable NS (Piro et al. 2017).If the total gravitational mass of the system, , is higher than a threshold mass,  thres (e.g., (Hotokezaka et al. 2011;Bauswein et al. 2013;Köppel et al. 2019;Barack et al. 2019;Kashyap et al. 2022;Perego et al. 2022)), then the remnant of the merger will promptly collapse to a BH.Although from the galactic population of BNSs the masses of NSs in BNS systems are expected to lie in a range 1.3 − 1.4 ⊙ , studies (Paschalidis & Ruiz 2019;Margalit & Metzger 2019) suggest, respectively, that up to 25% and 32% of BNS mergers might produce a rapid collapse to BH.For GW190425, it is estimated that the probability of the binary promptly collapsing into a BH after the merger is 96% , with the low-spin prior, or 97% with the high-spin prior (Abbott & et al. 2020).
BNS simulations in recent years (see, e.g., Shibata et al. (2003Shibata et al. ( , 2005)); Shibata & Taniguchi (2006); Kiuchi et al. (2009); Rezzolla et al. (2010); Hotokezaka et al. (2013); Kastaun et al. (2013); East et al. (2016); Endrizzi et al. (2016); Dietrich et al. (2017a,b); Endrizzi et al. (2018); Paschalidis & Ruiz (2019); Most et al. (2019); East et al. (2019); Ruiz et al. (2019Ruiz et al. ( , 2020)); Tootle et al. (2021); Papenfort et al. (2022); Sun et al. (2022)) have revealed that the amount of disk mass surrounding BH and the amount of ejected matter from the system are strongly dependent on the mass of the system, mass ratio, equation of Table 1.We list the BNS systems studied in this paper.The first column lists the model names, e.g.M32q10 refers to a BNS system with a total baryonic mass of 3.2 ⊙ and mass ratio  = 1.0.While  1,2 represent the baryonic mass of the components,  1,2 are the gravitational mass of each neutron star at infinite separation. 0 ,  0 and  0 are, respectively, the total initial ADM mass, the total initial angular momentum and the initial orbital frequency of the system.The last three columns show the compactness parameters of each NS and the combined dimensionless tidal deformability Λ, all computed at infinite separation.state (EOS), spin, and magnetic field of neutron stars.In addition, the effects of the BNS mass, mass ratio, EOS, spin, and magnetic field constitute a highly degenerate parameter space, which makes precise arguments difficult to make, as long as simulations that include all of the above are still lacking.Because nearly all matter is swallowed by the remnant BH, it is expected that equal-mass BNS systems with > thres will have a disk with negligible mass and a low amount of ejected matter.High-mass mergers, such as the ones discussed in this study, are important to study because they can provide insight into high-mass neutron stars and the remnants produced by their mergers.We can also learn about the internal structure of massive neutron stars (e.g., their EOS) by analyzing these types of mergers.We conducted simulations of high-mass binary systems that underwent rapid collapse in this study, and we investigated the effects of mass ratio and total mass on gravitational wave emission and system dynamics.The paper is organized as follows: we review the numerical setup and initial data for our models in Section 2. The numerical results are presented in Section 3. Finally, we discuss our findings in Section 4. We use a system of units in which  =  =  ⊙ = 1, unless specified otherwise.

MODELS AND NUMERICAL SETUPS
We consider a set of irrotational equal and unequal mass BNS systems with a total baryonic mass between 2.8 − 4.0 ⊙ .We computed the initial data using the pseudo spectral elliptic solver LORENE (Gourgoulhon et al. 2001;Gourgoulhon et al. 2016) assuming irrotational NSs on a quasi-circular orbit.For all models, the initial coordinate separation is 40 , corresponding to ≈ 3 orbits before the merger.We report the initial parameters of our models in Table 1.The first column refers to the names of the models, e.g.M32q10 means that the initial total baryonic mass of the BNS is 3.2 ⊙ and that the mass ratio is 1.0 (the mass ratio is defined such that  =  1 / 2 > 1).From the second to fifth columns we show, respectively, neutron stars' initial baryonic masses ( 1,2 ) and gravitational ( 1,2 ) masses at infi-nite separation. 0 ,  0 and  0 are the values of the initial ADM mass, total angular momentum, and orbital frequency.The compactness parameters of each NS are given as  1,2 where  1,2 =   1,2 / 1,2 2 and  1 and  2 are computed considering the NS at infinite separation.The last column refers to the reduced tidal parameter, Λ (Favata 2014): where the quadrupolar tidal parameter of the individual stars (Flanagan & Hinderer 2008;Damour & Nagar 2010) is ,  = 1, 2, and (2)  is the dimensionless Love number (Damour 1983;Hinderer 2008;Damour & Nagar 2009;Binnington & Poisson 2009).The notation (1 → 2) indicates a second term identical to the first, but 1 and 2 are exchanged.
This work employs the SLy EOS (Douchin & Haensel 2001) to describe the NS matter.To build the initial data, we used the EOS table provided by Parma Gravity Group1 .During the evolution, we instead used the seven piece polytropic version (SLyPP) described in De Pietri et al. (2016) and added a thermal component given by Γ thermal = 1.8.
To solve the general relativistic hydrodynamic (GRHD) equations, we used the publicly available Einstein Toolkit (ET) (Etienne et al. 2021;Löffler et al. 2012) code, based on Cactus Computational Toolkit 2 .In particular, we used the nineteenth release version of ET named "Katherine Johnson" (ET_2021_11) (Brandt et al. 2021).
We evolve the space-time metric using the BSSN formalism (Nakamura et al. 1987;Shibata & Nakamura 1995;Baumgarte & Shapiro 1998;Alcubierre et al. 2000Alcubierre et al. , 2003) ) as implemented in the McLachlan thorn3 .We use the GRHydro code (Baiotti et al. 2005;Hawke et al. 2005;Mösta et al. 2013) to solve the GRHD equations and employ a fifth-order WENO-Z reconstruction method (Borges et al. 2008).The time integration is handled with a fourth-order Runge-Kutta method (Runge 1895; Kutta 1901) with a Courant factor of 0.4.We also set an artificial atmosphere of    = 1.01 × 10 −11 .We also would like to state that the value of    /  for the less massive companion, which is from model M32q20 and has 1.07 ⊙ baryonic mass, is 9.08 × 10 −9 .
We employ an adaptive mesh refinement (AMR) approach provided by the Carpet driver (Schnetter et al. 2004) and we force the finest grids to follow each NS during the inspiral.The grid hierarchy consists of seven nested refinement levels with a 2 : 1 refinement factor.The resolution of each simulation is characterised by the resolution of the innermost grid, which is  = 0.16 ≈ 237 m, while the radius of the outer boundary of the numerical domain is 1024 ≈ 1514 km.We apply reflection symmetry along the z=0 plane.After merger, if the final remnant's mass exceeds the EOS's threshold mass, the AHFinderDirect thorn (Thornburg 2003) is used to detect the formation of an apparent horizon and extract the properties of the BH.
Note that results presented in this paper are extracted from the simulations having the standard resolution (SR) of  = 0.16.For two models, we also run a high resolution (HR) and low resolution (LR) simulations with  = 0.12 and  = 0.20 to estimate our numerical error.The accuracy of our results is discussed in Appendix A.

GW Extraction and Analysis
We computed the GWs from −7  to 10 , where  = 0  corresponds to the maximum strain amplitude and refers to the merger time, and the GWs include approximately the last three orbits before merger.
We extract the GW signal at a coordinate radius of ≈ 444  from the BNS center of mass, calculating the Newman-Penrose scalar  4 (Newman & Penrose 1962) (Eq.2), provided by ET module WeylScal4 , and decomposed in spin-weighted spherical harmonics using the ET module Multipole: Since  4 is the second time derivative of the GW strain, it is required to integrate it twice in time to obtain ℎ  .For the time integration and further GW signal analysis, we used a publicly available Python module, Kuibit (Bozzola 2021), and computed ℎ + and ℎ  with Eq. 3 via the FFI method described in (Reisswig & Pollney 2011): For all of our GW analysis, the complex combination of the extracted GW amplitude polarizations, ℎ = ℎ + − ℎ × , is used, and only the dominant mode, i.e.  =  = 2, is considered.Note that the ℎ 22 mode we extracted from the simulation is the amplitude of the (2,2) mode without the spherical harmonics term.Therefore it does not take into account the effect of the viewing angle.We also plot the waveform as a function of retarded time (Eq.4), where the areal radius is  = √︁ ()/4, and () is the surface of the sphere of coordinate radius (): The instantaneous frequency of GW is computed as where  = (ℎ × /ℎ + ) is the phase of the GW signal.The frequency at merger is  merger =   ( = 0).The GW strain (upper) and the phase velocity (bottom) for equal (left) and unequal (right) mass models are shown in Figure 1.Gravitational waveforms include the last part of inspiral, approximately three orbits before merger, and the coalescence and ringdown stages for all models.Models with larger total mass show, as expected, higher GW amplitudes at merger, while models with higher mass ratios have smaller amplitudes.
Because of the prompt formation of a BH, all GW signals go to zero less than ≈ 2 ms after the merger.The frequency evolution for all models is shown in the bottom panels.Before the merger stage, the frequencies are nearly identical for all models and increase monotonically with GW amplitude in time.Except for high mass ratio models, the frequency increases after merger due to increased compactness and faster rotation of the remnant, as explained in Endrizzi et al. (2016) and Kastaun et al. (2016).They behave differently for the latter in that, while model M32q20 oscillates around 1  , the frequencies of the other models rapidly increase shortly before 1 .
We then compute the amplitude spectral density (ASD) as , where h+ and h× are the Fourier transforms of ℎ 22 + and ℎ 22 × computed from the beginning of the simulations up to 10 after merger.Figure 2 shows the amplitude spectral density of the GW signals for equal (left panel) and unequal (right panel) models placed at a distance of 50   and the sensitivity curves of current (aLIGO+, KAGRA) and future planned (Einstein Telescope, Voyager, Cosmic Explorer) detectors.The initial frequency peak at around 1000  is equal to double of the initial orbital frequency of each model, and it just indicates the beginning of the inspiral stages in our simulations.From Figure 2 we can see that all systems can be observed in their merger phase by all detectors but KAGRA (which would be able to see only the earlier part of the inspiral, not simulated here).
Besides, we observe plateaus between 4000−7000  for all models, but M38q10 and M40q10 which are those that form a BH earlier.
For the unequal-mass models, the frequency of plateaus changes with the mass ratio and is more evident for larger mass ratios.Moreover, we realise that frequencies of plateaus show a turning point.While a more noticeable plateau is observed in models with higher mass ratios, the clarity decreases as we move towards models with equal mass.Moreover, plateaus are observed at higher frequencies and higher ASD values as the mass ratios decrease, while the ASD value decreases after  = 1.0 − 1.2.This suggests that there is a turning point in this mass ratio range.In Kiuchi et al. (2010) it is suggested that plateaus at high frequencies are related to both the formations and evolution of black holes and their surrounding disks.Therefore, we compare the ASD of model M30q10, which shows delayed collapse into a BH surrounded by torus after a short-lived HMNS phase, and prompt collapse models with torus (e.g.M32q10) and without torus (e.g.M34q10, M36q10).We realise that those plateaus are only related to the case of prompt collapse models regardless of having a torus or not.
We calculated the quasi-normal-mode frequencies of the corresponding BH mass and spin using Eq.6 from Nakamura et al. (1987).
According to Kiuchi et al. (2010), the peak frequency and width of the plateau show a correlation with the disk mass.Thus, these plateaus provide us with information about the formation and evolution of matter around the central object.If the frequency increases, we observe a decrease in the amplitude spectral density, as expected, due to the QNM ring-down of the formed BH.In Fig. 2, the black star marks show the QNM ring-down frequencies for each model.

Remnant Properties
We consider the matter whose time-component of the four-velocity is   < −1 as ejected matter, and we indicate its mass as    .We also define the disk mass as the baryonic mass of regions with rest mass density above the artificial atmosphere,    = 6.2 × 10 6 / 3 , and that is not ejected (i.e., we exclude  ejec from the disk mass).
Table 2.We report the values of numerical results from the standard resolution (SR) simulations and the value inside brackets (third row) shows the error on the last digit as the absolute semi-difference between low resolution (LR) and high resolution (HR) simulations, except for the disk and ejecta masses where we provide the relative errors.Details on the estimates of such errors can be found in the Appendix A. The second column shows the instantaneous frequency at merger,  merger , corresponding to the maximum GW amplitude.The third column is the elapsed time,   , from the merger to the formation of BH.   is the mass of the remnant.If the remnant is a BH, it is a quasi-local measure mass, otherwise, it is the ADM mass of the remnant NS.The fifth to last columns represent, respectively, the dimensionless spin of the BH,   , the disk mass,  disk , the amount of ejected matter,    and total radiated energy in GW,   .Besides, to consider the effects of the artificial atmosphere's density to disk mass and ejected matter, we simulated a model of M32q10 with three different atmosphere densities.Details are given in Appendix A.
All simulations run in this study result in the rapid collapse into a BH, except M2810, M30q11 (hyper-massive NS) and M30q10 (delayed collapse into a BH).In Table 2 we report the formation time of the BH after the merger.As seen from the values of   , all promptly collapse models form a BH before 0.9 .The results suggest that for massive binary systems, the rapid formation of a BH occurs sooner because there is not enough time for the redistribution of angular momentum to support the remnant, resulting in an earlier collapse into a BH.As suggested in Shibata et al. (2003) and Shibata & Taniguchi (2006), nearly all the matter tends to fall into the BH in the prompt collapse scenario.Moreover, it is realised that the elapsed time for the formation of a BH changes in a parabolic way with the increase of the mass ratio.When the mass ratio is between  = 1.0 − 1.4, the required time to form a BH increases, but after  = 1.4, it starts to decrease.
We report the final mass and dimensionless spins of BH found by AHFinder after the merger in Table 2.The final black hole swallows 94 − 99% of the initial binary mass.Similarly, we estimate that approximately 70 − 80% of the initial total angular momentum is transferred to the final BH with dimensionless spin between 0.66 − 0.79 as also suggested in Kiuchi et al. (2009).We also estimated that the energy carried away by gravitational radiation should be between 0.01 − 0.07 ⊙  2 .
Figure 3 shows the dependence of  merger ,  BH ,  BH ,  GW ,  disk and  eject on the different initial total mass of the system and mass ratios for the eighteen simulations performed in this study.As can be seen in Figure 3, the relationships between some of the parameters are clearly visible.
We plot in Fig. 3 E GW , M eject , M disk ,  BH , M BH , t BH , and f merger in function of the initial total baryonic mass of the BNS (  ) and mass ratio () (Table 2).The dashed black curve is fit to a secondorder polynomial using statistical analysis given by Equation 7. The uncertainties in the equations are expressed as percentages in paren-  2 from total initial mass (left panels) and mass ratios (right panel).Dashed lines on both plots are fits calculated by considering only results from promptly collapsed models.
theses.We believe these relations can be used as a prediction for future simulations using similar equations of state.
In Fig. 3, we show how the values extracted from our simulations change with initial binary neutron star system mass (left panels) and mass ratio (right panels).From the left panels in the figure, it can be seen that the dimensionless spin (mass) of the BH decreases (increases) with increasing BNS mass.When increasing the BNS mass, also the disk mass and the amount of ejected matter from the system decrease.Energy radiated via GW is also increasing for higher mass binary neutron star systems.From these plots, we can estimate that, in the prompt collapse case, the BH swallows 98 − 99% of the initial mass of the BNS systems and that the remaining 1 − 2% of the initial mass is radiated away from the system as   .On the other hand, for systems with different mass ratios (right panel of Fig. 3) the final mass of the BH decreases with the increase of mass ratio while the BH spin decreases.Similarly, we observe more massive disks and a higher amount of ejected matter, which are also consistent with the decrease in   for higher .

Dynamics of Systems During Evolution
In Fig. 4, we present the evolution of the rest mass density for the equal mass BNS models (M28q10, M30q10, M32q10 and M34q10) in the equatorial plane (upper panel) at times −3, 0, 2, 10 .The columns indicate times, while the rows refer to the models.All models are inspiralling and losing matter from the tidal tails due to tidal emission 3  before the merger (first column in Fig. 4).At  = 0.00  (second column in Fig. 4), the cores of the two neutron stars merge.Due to the high orbital speed (i.e., high angular momentum), matter is released in a spiral flow.After 2  (third column in Fig. 4), we show that M32q10 and M34q10 collapse into a BH which rapidly accretes the surrounding matter.However, model M30q10 shows the formation of a short-lived HMNS right after the merger and the formation of a BH 2.6  after the merger.Similarly, the remnant of M28q10 is a high-mass neutron star with a thick torus.The comparison of the four models in the figure shows that the disk structure changes with the initial binary mass.The mass of the disk is crucially decreasing for high-mass BNS systems.The bottom panel of Fig. 4 shows the four models on the meridional plane.
Besides, in Fig. 5, we present the evolution of the rest mass density    , as well as the mass ratio and initial binary mass.We showed that, as expected, the disk mass and the amount of dynamically ejected matter either increase or decrease, respectively, with an increase in mass ratio and total binary mass.
We reported our estimates for the disk mass and amount of ejected matter and compared them to the literature's current fits.Our models' estimates of mass ejection agree with the fit, but the amount of disk mass is less than the fit's estimate.
Moreover, we calculated gravitational waves and their spectra.We compared the amplitude spectral density of our models to the KA-GRA, aLIGO+, Voyager, Einstein Telescope, and Cosmic Explorer sensitivity curves.We propose a possible feature in the amplitude spectral density that could indicate rapid collapse into a BH following the merger of two neutron stars, and that could be observed using future-planned detectors such as the Einstein Telescope and Cosmic Explorer.
Additionally, from the evolution of the rest-mass density, we discussed the possible impact of initial binary mass and mass ratio on the geometry and mass of the accretion disk that will be formed around the post-merger NS or BH.
Table A2.We reported the maximum relative error on the baryon mass for the all models simulated in this paper.The symbol (*) indicates that the model does not have enough matter to collapse into a BH after the merger.Therefore, we reported the maximum baryon mass error during the whole simulation.

Figure 1 .
Figure 1.The left and right figures represent, respectively, equal and unequal mass models.Upper panels: extracted gravitational waveforms of the dominant mode (l,m=2,2) computed at  = 444.The thin lines show the GW signal ℎ 22 = ℎ +,22 , while the thick outer lines represent the GW amplitudes, |ℎ| = √︁ |ℎ + | 2 + | ℎ × | 2 and it does not take into account the effect of the viewing angle.Bottom panels: instantaneous frequency, i.e. the time derivative of the GW phase of each model.Note that the instantaneous frequency is smoothed by convoluting with a blackman window function in 0.1 .

Figure 2 .
Figure 2. Amplitude spectral densities (ASD) of the GW signals placed at a distance of 50 .The sensitivity curves of the current and the future-planned detectors (KAGRA, aLIGO+, Voyager, Cosmic Explorer and Einstein Telescope) are shown as shaded thick lines and they are the ones included in the public version of the Kuibit Python library(Bozzola 2021).Vertical dotted lines correspond to the merger frequencies estimated from the instantaneous frequency at the merger.The quasi-normal-mode frequencies of the final black holes are marked with black star markers.

Figure 3 .
Figure 3. Dependence of the numerical results given in Table2from total initial mass (left panels) and mass ratios (right panel).Dashed lines on both plots are fits calculated by considering only results from promptly collapsed models.

Figure 4 .
Figure 4. Snapshots of the rest mass density evolution for the equal mass models.Upper panels: the rest-mass density is plotted on the equatorial plane 3  before the merger, at merger  = 0, 2 and 10  after the merger.Each row represents a different model: top to bottom, M28q10, M30q10, M32q10 and M34q10.Bottom panels: comparison of the rest mass density snapshots in the equatorial and meridional planes at the end of the simulations ( = 10.00 ).

Figure 5 .
Figure 5. Snapshots of the rest mass density evolution for the unequal mass models.Upper panels: the rest-mass density is plotted on the equatorial plane 3  before the merger, at merger  = 0, 2 and 10  after the merger.Each row represents a different model: top to bottom, M32q10, M32q13, M32q16 and M32q20.Bottom panels: comparison of the rest mass density snapshots in the equatorial and meridional planes at the end of the simulation ( = 10.00 ).

Table A3 .
Disk and ejecta mass for simulation of the model M32q10 with high, standard and low NS atmosphere density values.