Accretion-induced prompt black hole formation in asymmetric neutron star mergers, dynamical ejecta, and kilonova signals

ABSTRACT

1 INTRODUCTION

Binary neutron stars (BNS) mergers are key astrophysical laboratories to explore the fundamental interactions in dynamical and strong gravity. This was clearly demonstrated by the observation of GW170817 (Abbott et al. 2017a, 2019a,b) and its related counterparts Abbott et al. (2017b). After GW170817, a second event (GW190425) compatible with a BNS source was reported, indicating a merger rate of 250–2810 Gpc³ per year (Abbott et al. 2020). The interpretation of current and future observations rely on quantitative simulations of astrophysically relevant binaries in the framework of numerical relativity (NR). In particular, observational signatures are strongly dependent on the possible neutron star (NS) masses and the still uncertain equation of state (EOS). The latter determine the properties of the final compact object, the eventual accretion merger remnant, and the observed gravitational and electromagnetic spectra (see Radice, Bernuzzi & Perego 2020 and reference therein for a recent review by some of us).

The possible NS mass range is ∼0.9–3 M_⊙, where the lower bound is inferred from the formation scenario (gravitational collapse) and from current observations (e.g. Rawls et al. 2011; Ozel et al. (2012). The upper bound is inferred from a stability argument (Buchdahl limit) and from precise measurements of ∼2-M_⊙ NSs in compact binaries containing a millisecond pulsar and a white dwarf (Demorest et al. 2010; Antoniadis et al. 2013; Cromartie et al. 2019). Coalescing circularized BNS were long expected to have nearly equal-mass NS with individual masses around M_A ∼ 1.35–1.4 M_⊙ and individual spin periods above the millisecond (Lattimer 2012; Kiziltan et al. 2013; Swiggum et al. 2015). For example, the source of GW170817 has a total mass of M ≃ 2.73–2.77 M_⊙ and a mass ratio q ∼ 1, with the largest uncertainties coming from the spin prior utilized in the analysis (Abbott et al. 2019a). This expectation was however challenged by GW190425, which is associated with the heaviest BNS source known to date with M ≃ 3.2–3.7 M_⊙ (Abbott et al. 2020). Spins distributions in GW170817 and GW190425 are both compatible with zero (Abbott et al. 2019a, 2020).

The mass ratio distribution in BNS (here conventionally defined as the ratio between the most massive primary and the secondary NS, i.e. q ≡ M_A/M_B ≥ 1) is very uncertain. BNS population from pulsar observations indicate mass ratios 1– ≤ q ≲ 1.4 (Lattimer 2012; Kiziltan et al. 2013; Swiggum et al. 2015). The mass ratio of GW170817 could be as high as q ∼ 1.37 (q ∼ 1.89) for low (high) spin priors. Similarly, the mass ratio of GW190425 can be as high as q ∼ 1.25 (q ∼ 2.5). Given the expected mass values and the recent observations, it is accepted that BNS mass ratios can reach ‘extreme’ mass ratio q ≲ 2. While these values are not as extreme as those that can be reached in black hole (BH) binaries, significant differences for the remnant and radiation signals are expected for BNS with q ∼ 1 and q ∼ 2.

NR simulations with microphysical EOS performed so far focused on comparable-mass cases and mass ratios q ≲ 1.4 (Sekiguchi et al. 2011a,b, 2015,2016; Neilsen et al. 2014; Palenzuela et al. 2015; Bernuzzi et al. 2016; Lehner et al. 2016; Radice et al. 2016, 2017, 2018a,b,c,d). The highest mass ratios of q = 1.5 and 2 have been simulated with a very stiff piecewise polytropic EOS (Dietrich et al. 2015, 2017), which is currently disfavoured by the GW170817 observation. Mergers of BNS with total mass M ∼ 2.7–2.8 M_⊙ and moderate mass ratios up to q ≲ 1.4 with an EOS supporting ∼2 M_⊙ are likely to produce remnants that are at least temporarily stable against gravitational collapse to BH, as opposed to remnants that collapse immediately to BH (prompt BH formation). However, the conditions for prompt BH formation at high-q have not been studied in detail to date. For a given total mass, moderate mass ratios can extend the remnant lifetime with respect to equal-mass BNS because of the less violent fusion of the NS cores and a partial tidal disruption that distributes angular momentum at larger radii in the remnant (Bauswein, Baumgarte & Janka 2013a). However, large-mass asymmetries q ≳ 1.6 can favour BH formation due to the larger mass of the primary NS.

Tidal disruption in asymmetric BNS can significantly affect the properties of the dynamical ejecta, favouring a redder kilonova peaking at late times (see e.g. Rosswog et al. 2018; Wollaeger et al. 2018). Moderate mass ratios up to q ∼ 1.3–1.4 are found to produce more massive discs than q = 1 BNS (Shibata, Taniguchi & Uryu 2003; Shibata & Taniguchi 2006; Kiuchi et al. 2009; Rezzolla et al. 2010; Dietrich et al. 2017). But BH formation can significantly alter the remnant disc properties, both in terms of compactness and composition (Perego, Bernuzzi & Radice 2019). In turn, this can impact the secular (viscous) ejecta component and the kilonova, with bright emissions generally favoured by the presence of a long-lived NS remnant (e.g. Radice et al. 2018d; Nedora et al. 2019). Moreover, assuming the Blandford–Znajek mechanism to be the mechanism launching the relativistic jet that produces a gamma-ray burst, a more massive disc is also expected to power a more energetic jet through a more intense accretion process (Shapiro 2017). Many of these aspects are currently not well quantified and they require NR simulations of high-mass-ratio BNS with microphysics.

In this work, we perform 32 new NR simulations with microphysical EOS, fixed chirp mass |$\mathcal {M}_{\rm c}=1.188\,{\rm M}_\odot$|⁠, and mass ratios up q = 1.8 for four microphysical EOS, including a new microscopic EOS BLh (Section 2). The simulations show that for a sufficiently high value of the mass ratio (and in an EOS-dependent way), the remnant promptly collapses to BH as consequence of the accretion of the companion on the massive primary NS (Section 4). These prompt-collapse dynamics is not well described by current NR fitting formulas. By analysing the gravitational waveforms, we further verify current quasiuniversal NR relations for the merger and postmerger gravitational waveforms in the high-q limit (Section 5). We find an overall agreement of the merger relations and characteristic postmerger GW frequencies. But the accurate modelling of postmerger waveforms with high q will require more simulations and improved methods than those currently employed. We discuss in detail the differences in the dynamical ejecta between the q = 1 and the high-q mergers in terms of overall ejecta masses, morphology, and composition (Section 6). High mass ratio and large chirp mass leading to prompt BH formation maximize the dynamical tidal ejecta mass, which is expelled with a peculiar geometry. The r-process nucleosynthesis in these neutron-rich ejecta result in bright (more luminous than the q = 1 case), redder, and temporally extended kilonovae (Section 7).

We employ SI units in most of the paper except for masses, reported in solar masses (M_⊙), lengths in km, and densities reported in g cm⁻³. Nuclear density is indicated as ρ₀ ≈ 2.3 × 10¹⁴ g cm⁻³. If units are not reported, we then use geometric units c = G = 1 in a context where those are more appropriate (e.g. Section 5 and appendices).

2 EQUATIONS OF STATE

In this work we consider four finite-temperature, composition-dependent EOS: the LS220 EOS (Lattimer & Swesty 1991), the SFHo EOS (Steiner, Hempel & Fischer 2013), the SLy4-SOR EOS (Schneider, Roberts & Ott (2017), and the BLh EOS (Bombaci & Logoteta 2018). All these EOS include neutrons (n), protons (p), nuclei, electron, positrons, and photons as relevant thermodynamics degrees of freedom. Cold, neutrino-less β-equilibrated matter described by these microphysical EOS predicts NS maximum masses and radii within the range allowed by current astrophysical constraints, including the recent GW constraint on tidal deformability (Abbott et al. 2017a, 2018, 2019a; De et al. 2018). All four models have symmetry energies at saturation density within experimental bounds. However, LS220 has a significantly steeper density dependence of its symmetry energy than the other models (see e.g. Lattimer & Lim 2013; Danielewicz & Lee 2014), and it could possibly underestimate the symmetry energy below saturation density.

The LS220 EOS is based on a non-relativistic Skyrme interaction with the modulus of the nuclear bulk incompressibility set to 220 MeV. Non-homogeneous nuclear matter is modelled by a compressible liquid-drop model including surface effects, and considers an ideal, classical gas formed by α particles and heavy nuclei. The latter are treated within the single nucleus approximation (SNA). The transition between homogeneous and non-homogeneous matter is performed through a Gibbs construction.

The SFHo EOS combines a relativistic mean field approach for the homogeneous nuclear matter to an ideal, classical gas treatment of a statistical ensemble of several thousands of nuclei in nuclear statistical equilibrium (NSE) for the inhomogeneous nuclear matter. The transition between the two phases is achieved by an excluded volume mechanism.

The SLy4 Skyrme parametrization was introduced in Douchin & Haensel (2001) for cold nuclear and NS matter. In this work, we employ its extension to finite temperature presented in Schneider et al. (2017) using an improved version of the LS220 model that includes non-local isospin asymmetric terms and a better treatment of nuclear surface properties, and treats the size of heavy nuclei more consistently. The transition between the uniform and non-uniform phase is achieved by a first-order transition, i.e. choosing the phase with lower free energy.

A main novelty of this work is the use of the BLh EOS, a new finite temperature EOS derived in the framework of non-relativistic Brueckner–Hartree–Fock (BHF) approach (Logoteta et al., in preparation). The corresponding cold, β-equilibrated EOS was first presented in Bombaci & Logoteta (2018) and applied to BNS mergers in Endrizzi et al. (2018). For the homogeneous nuclear phase, this EOS employs a purely microphysical approach based on a specific nuclear interaction. Consistently with Bombaci & Logoteta (2018), the interactions between nucleons are described through a potential derived perturbatively in the chiral effective field theory (Machleidt & Entem 2011). Specifically, the local potential reported in Piarulli et al. (2016) and calculated up to next-to-leading order (N3LO) was used as the two-body interaction. This potential takes into account the possible excitation of a Δ-resonance in the intermediate states of the nucleon–nucleon interaction. The above potential was then supplemented by a three-nucleon force calculated up to N2LO and including again the contributions from the Δ-excitation. The parameters of the three-nucleon force were determined to reproduce the properties of symmetric nuclear matter at saturation density (Logoteta, Bombaci & Kievsky 2016). For the non-homogeneous nuclear phase, there is no straightforward extension of these microphysical methods to subsaturation densities. Thus, the low-density part (n ≤ 0.05 fm⁻³) of the SFHo EOS has been smoothly connected to the high-density BLh EOS. This necessary extension has been tested with different finite-temperature, composition-dependent tabulated EOS (Hempel & Schaffner-Bielich 2010). They all use (1) relativistic mean field approaches for the homogeneous phase, (2) an ideal, classical gas of a statistical ensemble of several thousands of nuclei in NSE for the non-homogeneous nuclear phase, and (3) an excluded volume mechanism to model the transition. No appreciable differences were found in all relevant quantities at subnuclear densities between different high-density treatments.

The LS220 and SLy4-SRO EOS are based on Skyrme effective nuclear interactions. In these models, thermal effects are introduced starting from a zero-temperature internal energy functional that contains an explicit nuclear density dependence. The interaction part of this functional is split into a quadratic term in the nuclear density (playing the role of a two-body nucleon–nucleon interaction) plus a term proportional to some power of the nuclear density. The latter term mimics the effect of many-body nuclear forces. The temperature dependence of the effective nuclear interaction is encoded in the effective mass dependence of the kinetic energy as well as in the single-particle potentials. The latter are calculated by the variation of the internal energy with respect to the neutron and proton densities. Assuming indeed a constant entropy, smaller effective masses translate into larger kinetic energies and thus higher matter temperatures. The LS220 EOS assumes that the effective nucleon mass is the bare nucleon mass at all densities while for the SLy4-SRO we have |$m_N^*/m_N = 0.695$| at saturation density, with |$m_N^*$| and m_N the effective and the bare nucleon masses, respectively.

In the relativistic Lagrangian underlaying the SHFo EOS, nuclear interactions are described by σ-, ω- and ρ-meson exchanges. The resulting Euler–Lagrange equations are then solved in mean field approximation. In this approach, thermal effects are included by introducing Fermi–Dirac distributions at finite temperatures for the various nuclear species. Mesons and nucleon fields, and consequently all thermodynamical quantities, acquire automatically a temperature dependence through the self-consistent solution of the mean field equations.

Differently from the other models considered in this work, thermal effects enter in a quite different way in the BLh EOS. The calculation of the Free energy in the BHF approach (Bombaci, Kuo & Lombardo 1993) requires first the determination of an effective in-medium nuclear interaction, starting from the bare nuclear potential. This effective interaction (G-matrix) is obtained by solving the Bethe–Goldstone integral equation that describes the nucleon–nucleon scattering in the nuclear medium and properly takes into account the Pauli principle. Finally, the nucleon single-particle potentials U_i(k, T) (i = n, p) are obtained through the integration of the on-shell G-matrix. U_i(k, T) is a sort of mean field felt by a nucleon of momentum k due to the presence of the surrounding nucleons. The determination of U_i(k, T) allows for the calculation of the Free energy from which all the other thermodynamical quantities can be derived. The procedure described above is complicated by the non-linear and non-local dependence of U_i(k, T) in Bethe–Goldstone equation. We finally note that this scheme provides many-body correlations that are beyond the mean field approximation. Such correlations are not present in the other EOS models considered in this paper.

2.1 EOS constraints and NS equilibrium models

Fundamental differences in the EOS models translate into different NS structures. For cold, non-rotating NSs, the considered EOS can support maximum masses in the range |$M_\text{max}^\text{TOV}\sim 2.06\!-\!2.10\,{\rm M}_\odot$|⁠, while the predicted radii of a 1.4-M_⊙ NS lay in the range R_1.4 ∼ 11.78–12.74 km. More specifically, LS220, SFHo, SLy4-SRO, and BLh EOS have |$M_\text{max}^\text{TOV}$| of 2.04, 2.06, 2.06, and 2.10 M_⊙, and R_1.4 of 12.8, 12.0, 11.9, and 12.5 km, respectively. The predicted maximum NS masses and the 1.4-M_⊙ NS radii are all compatible at 1σ level with the recent detection of an extremely massive millisecond pulsar (Cromartie et al. 2019) and with results obtained by the NICER collaboration (Miller et al. 2019; Riley et al. 2019), although being systematically on the lower side. Note that EOS allowing NS radii R_1.4 ≫ 13 km are currently disfavoured by both GW BNS and X-ray pulsar observations (Abbott et al. 2019a; Miller et al. 2019; Riley et al. 2019).

Finite-temperature effects introduce additional pressure support. On the one hand, for the typical central entropies expected for nuclear matter during a BNS merger (s ≲ 2k_B/baryon), this additional support is not sufficient to significantly alter the maximum TOV mass (Kaplan et al. 2014) or the central baryon density due to the large degree of degeneracy of matter above saturation density. On the other hand, thermal effects can provide a more significant impact for matter at lower densities, increasing the NS radius. In Fig. 1, we report equilibrium sequences in the mass–radius and mass–central density plane obtained for the EOS used in this work, considering both a cold (continuous lines) and an isentropic (dashed lines) EOS with s = 2 k_B baryon⁻¹. Due to thermal effects, R_1.4 increases by |$15.6 {{\ \rm per\ cent}}$| for the LS220 EOS and |$36.4 {{\ \rm per\ cent}}$| for the SLy4 EOS. while for the BLh and the SFHo EOS, the variation is |${\sim}21\!-\!22 {{\ \rm per\ cent}}$|⁠. The different relative impacts on the NS radius clearly correlate with the different values of the nucleon effective mass.

Figure 1.

Open in new tabDownload slide

Equilibrium NS sequences obtained for the different EOS used in this work. Solid lines correspond to irrotational, T = 0 configurations; dashed lines to irrotational, isentropic (s = 2 k_B baryon⁻¹) configurations; and dotted lines represent T = 0, rigidly rotating NS spinning at the mass shedding limit (Ω = Ω_K). Left-hand panel: gravitational mass versus radius (equatorial for rotating NS). Right-hand panel: gravitational mass versus central density. Markers along the cold, non-rotating sequences indicate the NSs used in this work: squares, triangles, and stars refer to q = 1, 1.67, and 1.8 binaries, respectively.

Rotational support also increases the maximum NS mass. For example, in the limiting case of rigid rotation at the Keplerian limit, the maximum NS mass is increased by |${\sim}20{{\ \rm per\ cent}}$| for all EOS models, as visible in Fig. 1, dotted lines. Since this affects the whole star, the NS radius is typically increased by |${\sim}40{{\ \rm per\ cent}}$|⁠, but at the same time, the central density is decreased by a similar amount, if one compares non-rotating and Keplerian NSs of identical masses. These properties emphasize the importance of using the full EOS (i.e. including thermal effects) in merger simulations. Thermal (and composition; see Kaplan et al. 2014) effects are indeed key to quantify the prompt-collapse dynamics, mass-shedding in the remnant, and disc properties.

3 SIMULATIONS

3.1 Methods

We construct initial data for irrotational binaries in quasicircular orbit solving the constraint equations of 3+1 general relativity in the presence of a helical Killing vector and under the assumption of a conformally flat metric (Gourgoulhon et al. 2001). The equations are solved with the pseudo-spectral multidomain approach implemented in the lorene library.¹ The EOS used for the initial data are constructed from the minimum temperature slice of the EOS table employed for the evolution, assuming neutrino-less β-equilibrium. Initial data have a residual eccentricity of ∼0.01, which is radiated away before merger (e.g. Thierfelder, Bernuzzi & Brügmann 2011b).

The initial data are then evolved with the 3+1 Z4c free-evolution scheme for Einstein’s equations (Bernuzzi & Hilditch 2010; Hilditch et al. 2013) coupled to general relativistic hydrodynamics. For the latter, we use the whiskythc code (Radice & Rezzolla 2012; Radice, Rezzolla & Galeazzi 2014a,b) that implements the approximate neutrino transport scheme developed in Radice et al. (2016, 2018d) and the general-relativistic large eddy simulations method (GRLES) for turbulent viscosity (Radice 2017). The interactions between the fluid and neutrinos are treated with a leakage scheme in the optically thick regions (Ruffert, Janka & Schäfer 1996; Rosswog & Liebendoerfer 2003; Neilsen et al. 2014) while free-streaming neutrinos are evolved according to the M0 scheme (Radice et al. 2018d). The latter is a computationally efficient scheme that incorporates an approximate treatment of gravitational and Doppler effects, which is well adapted to the geometry of BNS mergers and free of the radiation shock artefact that plagues the M1 scheme (Foucart et al. 2018). The turbulent viscosity in the GRLES is parametrized as σ_T = ℓ_mixc_s, where c_s is the sound speed and ℓ_mix is a free parameter sets the intensity of the turbulence. For the simulations of this work, σ_T is prescribed as a function of the rest-mass (baryon) density using the results of the high-resolution general relativistic magnetohydrodynamics simulations results of an NS merger of Kiuchi et al. (2018). A detailed description of the model can be found in Radice (2020); simulations with this model are also presented in Perego et al. (2019) and Nedora et al. (2019).

We remark that the GRLES method introduces parabolic terms for which there is no maximum characteristic velocity. However, parabolic equations still have an effective, wavelength-dependent, propagation speed (Weymann 1967; Kostadt & Liu 2000). Accordingly, only disturbances that have spatial scale that are small compared to the mixing-length parameter are propagated with an effective velocity larger than the speed of light. These modes are absent in our simulations because the mixing length is always set to be smaller than the minimum grid scale in the simulation. Indeed, the GRLES method becomes invalid precisely when the mixing length becomes comparable with the grid scale, which would correspond to turbulent motion on a scale that is resolved in the simulations and should be included directly, not through a subgrid model. This is also the reason why it is possible to integrate the GRLES equations using an explicit time-integration scheme. Moreover, it is possible to show that in the long-wavelength, low-frequency limit that is relevant for us, the solution of the parabolic model is always arbitrarily close to those of an associated hyperbolic model obtained with the introduction of relaxation terms (Nagy, Ortiz & Reula 1994). In other words, the parabolic and the (significantly more complex) hyperbolic models of turbulent viscosity should give the same results in our context.

whiskythc is implemented within the cactus (Goodale et al. 2003; Schnetter et al. 2007) framework and coupled to an adaptive mesh refinement driver and a metric solver. The space–time solver is implemented in the ctgamma code (Pollney et al. 2011; Reisswig et al. 2013a), which is part of the Einstein Toolkit (Loffler et al. 2012). We use fourth-order finite-differencing for the metric’s spatial derivatives method of lines for the time evolution of both metric and fluid. We adopt the optimal strongly-stability preserving third-order Runge–Kutta scheme (Gottlieb & Ketcheson 2009) as the time integrator. The time-step is set according to the speed-of-light Courant–Friedrich–Lewy (CFL) condition with CFL factor 0.15. While numerical stability requires the CFL to be less than 0.25, the smaller value of 0.15 is necessary to guarantee the positivity of the density when using the positivity-preserving limiter implemented in whiskythc.

The computational domain is a cube of 3024 km in diameter whose center is at the center of mass of the binary. Our code uses Berger–Oliger conservative AMR (Berger & Oliger 1984) with subcycling in time and refluxing (Berger & Colella 1989; Reisswig et al. 2013b) as provided by the Carpet module of the Einstein Toolkit (Schnetter, Hawley & Hawke 2004). We set up an AMR grid structure with seven refinement levels. The finest refinement level covers both NSs during the inspiral and the remnant after the merger and has a typical resolution of h ≃ 246 (grid setup named LR), 185 (SR), or 123 m (HR).

BH formation is indicated by the appearance of an apparent horizon (AH) that is computed with the module AHFinderDirect (Thornburg 2004). With the gauge conditions employed in the simulations, the BH is formed and simulated as a puncture (Thierfelder et al. 2011a; Dietrich & Bernuzzi 2015). In a first series of simulations, the AH finder could not find an AH in the simulations using the GRLES scheme. We have thus repeated them and found that in those cases was necessary to (i) increase the number of more guess spheres for the finder, and (ii) switch off the GRLES scheme in regions with α < 0.1 in order to compute the AH robustly. The latter is analogous to the ‘hydro-excision’ implemented in many codes to facilitate AH location. This way, we obtained horizon data for all the LR and most of the SR simulations. We could not rerun the HR simulations for which the AH was not found initially for lack of computational resources. Finally, the employed grid structure is not optimal to follow the dynamics of the BH+disc remnant; thus, simulations are stopped ∼5–10 ms after BH formation.

3.2 BNS models

We consider 10 binaries with fixed chirp mass |${\cal M}_{\rm c}\simeq 1.188\,{\rm M}_\odot$| and simulate them at different resolutions. The chirp mass is |${\cal M}_{\rm c}=M\nu ^{3/5}$|⁠, where ν = M_AM_B/M² = q/(1 + q)². The main properties of the BNS initial data are summarized in Table 1. We simulated the equal-mass case and mass ratio q = 1.67 for all the EOS with the GRLES scheme. The highest mass ratio simulated are q = 1.8 for the BLh and SLy EOS. A subset of models were simulated also without turbulent viscosity to directly assess its impact on the merger dynamics and on the ejecta properties. The initial separation between the NS is set to |$45\, {\rm km}$|⁠, corresponding to approximately four to six orbits to merger. Note that similar equal-mass LS220 and SFHo BNS were already presented in Perego et al. 2019, but the mass here is slightly larger. An equal-mass SLy4 without turbulent viscosity was instead presented in Endrizzi et al. (2020).

4 MERGER DYNAMICS AND REMNANT

Starting at a GW frequency of ∼570 Hz, the binaries revolve for approximately four to six orbits before reaching the moment of merger. The latter is defined as the peak amplitude of the (2, 2) GW mode and marks the end of the chirp signal. A summary of the merger dynamics for all the runs is given in Fig. 2, which shows the maximum mass density (fluid frame) and the minimum of the lapse function, α. BH formation is indicated by the lapse dropping below α ≲ 0.3. Note that the 1+log slicing of the spherical puncture has lapse function at the horizon α_AH ≃ 0.376 (Hannam et al. 2007, 2008), but punctures formed in our simulations have dimensionless spins ∼0.7 for which α_AH ≃ 0.3 (see Appendix  A). At the same time, the lapse decreases below α_AH, the maximum density increases beyond 6ρ₀, and it is then unresolved on the grid due to the gauge conditions (Thierfelder et al. 2011a).

Figure 2.

Open in new tabDownload slide

Evolution of the maximum density (normalized to nuclear saturation density) and of the minimum lapse in all simulations. The horizontal line indicates the value of the lapse at the horizon for a puncture with the indicated dimensionless spins (cf. Table 2). Data refer to simulations with turbulent viscosity and resolution SR.

The remnants of BLh q = 1.8, LS220 q = 1.67, SFho q = 1.67, and SLy4 q = 1, 1.67, 1.8 collapse to BH within ∼2–3 ms from merger. We call prompt BH collapse mergers those in which the NS cores collision has no bounce but instead the remnant immediately collapse at formation; see Table 2. This usually happens within 1–2 ms from the moment of merger and can be identified by the maximum density monotonically increasing to the collapse; see Fig. 2. Note that this definition of prompt collapse implies negligible shocked dynamical ejecta because the bulk of this mass emission comes precisely from the (first) core bounce (Radice et al. 2018d). The BH masses of the prompt-collapsed remnants are M_BH ≃ 2.49, 2.44, 2.45, 2.47, and 2.52 M_⊙ for BLh q = 1.8, LS220 q = 1.67, SFho q = 1.67, and SLy4 q = 1.67, 1.8, respectively. The BH spins are a_BH ≃ 0.66, 0.7, 0.68, 0.69, and 0.66. The remnants of LS220, SFHo, and SLy q = 1 also form BH within the simulated time and they have M_BH = 2.41, 2.41, and 2.38 and spins a_BH = 0.5, 0.75, and 0.76, respectively. The SLy4 q = 1 merger was simulated in Endrizzi et al. (2020) without viscosity, and in that case, the remnant survives for ∼10 ms. The earlier collapse here is a consequence of the angular momentum redistribution by the subgrid model for turbulent viscosity (Radice 2017). Overall, these results for the BH spins consistently indicate an upper limit on the BH rotation of a_BH ≲ 0.8, also when including q ∼ 2 BNS (Kiuchi et al. 2010; Bernuzzi et al. 2014, 2016; Dietrich et al. 2017). In the following, we first discuss the details of the BH formation highlighting the effect of high mass ratio and the main differences with respect to the (well-studied) equal-mass cases. Then, we discuss the properties of the remnant discs.

For comparable masses, the NS cores enter in contact before reaching the moment of merger (Thierfelder et al. 2011b) and the last two 2–3 GW cycles before the amplitude’s peak are emitted by the cores collision and remnant formation. At high mass ratios, a new effect is the tidal disruption of the companion and its accretion on to the primary NS. This has been reported also in previous simulations with a stiff polytropic EOS (Dietrich et al. 2017), and we confirm it here for softer and microphysical EOS. As a representative example, we show in Fig. 3 the cases of BLh q = 1 versus q = 1.8. The accreting material has initially low temperatures but as soon as the accretion becomes more massive and faster the temperature raises. At approximately the time of the snapshot, the accreting material shocks against the primary NS core and there the temperature raises up to ∼100 MeV. As a consequence of this shock, some material becomes unbound, although the exact amount of ejecta cannot be confidently measured in the simulations (see Section 6).

Figure 3.

Open in new tabDownload slide

Snapshots of premerger dynamics for BLh q = 1.8 (top panel) and 1.0 (bottom panel) simulations. Shown is the rest-mass density in the orbital plane at ∼9 ms corresponding to the third orbit from the beginning of the simulations and two orbits to the moment of merger. The companion in the q = 1.8 BNS is tidally disrupted and a significant accretion on to the primary is taking place. Accretion starts approximately after one orbits from the beginning of the simulations.

The new aspect highlighted by our simulations is the dynamics of prompt collapse for high-mass-ratio BNS. In a q ∼ 1.5–2 binary, the tidal disruption and accretion of the companion NS on to the massive primary NS can drive the remnant unstable and causes a prompt collapse to BH. The process is shown in Fig. 4 (top panels) in a 3D volume rendering of the rest-mass density for the representative case of the BLh EOS. The BLh q = 1.8 has a rather massive primary NS with M_A = 1.856 M_⊙ as compared to the maximum TOV mass for the BLh EOS (⁠|$M_\text{max}^\text{TOV}=2.103\,{\rm M}_\odot$|⁠), and a companion NS of small compactness (M_B = 1.020 and C_B ≃ 0.12). The companion NS is almost completely destroyed by tidal effects and its accretion results in the prompt formation of a BH surrounded by a massive accretion disc. In contrast, the lower mass-ratio and equal-mass binaries with the same chirp mass produce a less compact remnant and none of them collapse to the end of the simulated time (middle and bottom panels). Note that the equal-mass BLh was evolved beyond the 80-ms postmerger. We checked with a sequence of simulations at intermediate-mass ratios that the behaviour is continuous in the mass-ratio parameter (See Appendix  B).

Figure 4.

Open in new tabDownload slide

3D volume rendering of the rest-mass density ρ in g cm⁻³ expressed in logarithmic scale for the BLh models. Each column represents a different time inside the simulation: merger time (left-hand panel), early postmerger (∼2 ms, middle panel), and later stages (∼10 ms, right-hand panel). In each row, we show a different mass ratio q = M_A/M_B: q = 1.8 (top panel), 1.67 (middle panel), and 1.0 (bottom panel). The BH AH is shown as a bright green isosurface of the lapse function α = α_AH.

Figure 5.

Open in new tabDownload slide

Masses of simulated BNS as function of the TOV maximum compactness as predicted by the zero-temperature EOS. Prompt collapsing remnant are marked with additional black circles. The solid line is the mass threshold given by equation (2) (q-independent): Binaries above the line are predicted to collapse. The dashed and solid lines are the thredshold model modified by the mass-ratio correction in equation (4). The plot refers to simulations with turbulent viscosity only.

Let us now discuss disc formation, evolution, and properties. Following a common convention, we define disc the baryon material either outside the BH’s AH or the one with densities ρ ≲ 10¹³ g cm⁻³ around an NS remnant. The baryonic mass of the discs are computed as volume integrals of the conserved rest-mass density |$D=\sqrt{\gamma }~W\rho$| from 3D snapshots of the simulations in postprocessing (γ is the three-metric’s determinant and W the Lorentz factor). Estimates for the disc masses are reported in Table 2. The disc mass is reported as measured at the time when it is maximum during the simulation. For the remnant collapsing to BH, this can be interpreted as the mass at BH formation, since the disc mass can only decrease with time due to accretion. For NS remnants, the disc (remnant at lower densitites) can also increase its mass over time as it acquires matter expelled from the higher densities shells.

Examples of the disc mass evolutions for different remnants are shown in Fig. 6. Note that we show the BLh q = 1.8 and LS220 q = 1 simulations at resolution SR but without turbulent viscosity and the q = 1.67 with viscosity but at LR because these are the longer data sets available to us (see below for a discussion about turbulence).

Figure 6.

Open in new tabDownload slide

Evolution of the disc rest mass in representative remnants. The disc is defined as the remnant outside the AH if the remnant has collapsed to a BH or as the portion of the remnant whose rest-mass density satisfies ρ < 10¹³ g cm⁻³ otherwise.

In the case of comparable-mass BNS, the accretion disc is formed during and after the merger. As time evolves, if the remnant does not collapse, it continuously sheds mass and angular momentum increasing the mass of the disc and generating outflows (Radice et al. 2018a; Nedora et al. 2019). This is why in Fig. 6, the accretion disc mass is increasing with time for these binaries. These processes terminate with BH formation, which is accompanied by the rapid accretion of a substantial fraction of the disc. An important consequence is that, in the case of comparable-mass ratio binaries, prompt BH formation results in very small accretion disc masses (Radice et al. 2018b; Kiuchi et al. 2019) because the mechanism primarily responsible for the formation of the disc is shut off immediately in these cases.

In high-mass-ratio BNS mergers, the companion star is tidally disrupted (Fig. 4). In these cases, the bulk of the accretion disc is constituted by the tidal tail, which is, for the most part, still gravitationally bound to the remnant. This tail is launched prior to merger. So massive accretion discs are possible even if prompt BH formation occurs (see also Kiuchi et al. 2019). In general, high-q binaries are found to generate more massive discs than binaries with the same chirp mass but lower q (Shibata et al. 2003; Shibata & Taniguchi 2006; Kiuchi et al. 2009; Rezzolla et al. 2010; Dietrich et al. 2017). The postmerger evolution of these discs is also very different. While in the massive NS case the central object pushes material into the disc and drives outflows, in the case of high mass-ratio binaries forming BHs, the fallback of the tidal tail perturbs the disc and drives rapid accretion on to the BH as evinced by the rapid decrease of the disc masses with time shown in Fig. 6.

These different formation mechanisms are imprinted in the structure and composition of the discs formed in comparable and very unequal mass binaries, as shown in Fig. 7. In the case of equal-mass binaries, the disk is composed of material squeezed out of the collisional interface between the NSs (Radice et al. 2018d); see also Fig. 4. This matter is heated to temperatures of tens of MeV before being pushed out of the central part of the remnant, so its electron fraction is reset by pair processes (Perego et al. 2019); see Fig. 7. Due to the absence of strong compression and shocks, the discs formed in high-mass-ratio binaries are initially colder and more neutron-rich (Fig. 7). Since high-q BNS mergers launch tidal tails to large radii, comparable-mass-ratio binaries create discs that are initially more compact and have higher Y_es. Besides the mass ratio, the structure of the disc is also strongly dependent on the nature of the remnant. In the case of BH remnants, the discs are typically more compact and thin than those around massive NS remnants. In the latter case, because of the additional pressure support, the discs reach higher densities ∼10¹³ g cm⁻³ and become partly optically thick to neutrinos (Perego et al. 2019; Endrizzi et al. 2020).

Figure 7.

Open in new tabDownload slide

Volume rendering of density at the end of the BLh q = 1.8 (left-hand panel) and 1 (right-hand panel) simulations. In each panel, on the left-hand side, we show the electron fraction distribution, while on the right-hand side, we plot the rest-mass density; for both quantities, we show only matter with densities ρ ≥ 3 × 10⁸gcm⁻³. The spatial scale is the same in the two panels.

5 GRAVITATIONAL WAVES

In this section, we analyse the GW signals computed from the simulations. The latter are too short for a quantitative comparison with inspiral-merger models (e.g. Akcay et al. 2019). Hence, we focus on the merger and postmerger signal. The new simulations allow us to verify (and extend) the quasiuniversal relations characterizing the merger and study the postmerger waveform (Bernuzzi, Dietrich & Nagar 2015b; Breschi et al. 2019).

Fig. 8 compares the new NR data of this paper (Table 3) with the fits of simulations of the CoRe collaboration for q ≤ 1.5 proposed in Breschi et al. (2019). The fits are consistent with the new data with q > 1.5 within the the uncertainties, indicating the robustness of the model (and especially the ansatz equation 8). The fits for the binding energy and angular momentum at the merger were not presented in Breschi et al. (2019) and are thus are given here in Appendix  C.

Figure 8.

Open in new tabDownload slide

Comparison between GW data at moment of merger extracted from the simulations introduced in this work (and listed in Table 3) and fits calibrated using NR simulations with q ≤ 1.5. The fits for energy and angular momentum are reported in Appendix  C, while frequency and amplitude fits are from Breschi et al. (2019). The quantities |$e_{\rm b}^{\rm mrg}$| and j^mrg are defined from the GW data as |$e_{\rm b}^{\rm mrg} = E_{\rm b}^{\rm mrg}/(\nu M)$| and j^mrg = J^mrg/(νM²). The results extracted from high-mass-ratio simulations are consistent with the current fits in the limit of the uncertainties.

Regarding the postmerger waveform, Fig. 9 (top panel) shows a comparison between the waveforms from the BLh BNS for the three mass ratio considered here. The figure clearly indicates that for similar (though not identical) initial frequencies, the moment of merger occurs earlier for unequal-mass simulations due to tidal disruption of the high-q binaries, where the companion has larger radius than the primary NS (see also Figs 3 and 4). As expected, the dependence of the waveform on q is smooth as shown explicitly in Appendix  B. Note that, in general, the postmerger amplitude is smaller for high q than for equal mass due to a less violent shock between the two NS cores and either a less compact remnant or the formation of a BH that quickly rings down to a stationary state.

Figure 9.

Open in new tabDownload slide

Example of GW signals. Amplitude and real part of the (2, 2) mode strains (top panel) and respective frequency-domain spectra (bottom panel) for the three simulations with BLh EOS at standard resolution. The vertical dashed lines mark the merger times, both in time- and frequency-domain, while the dotted lines in the bottom panel show the initial point of the simulations. The frequency-domain spectra are compared with NRPM postmerger model evaluated with the same physical values, except for the prompt-collapse case (BLh, q = 1.8).

The only new unequal-mass simulation with a long postmerger GW signal is the BLh with q = 1.67. For this case, the value of the characteristic postmerger frequency f₂ is properly captured from NR fits presented in Breschi et al. (2019): From the simulation, we get f₂ ≈ 3.31 kHz, while for the same binary, the NR fit predicts |$f_2^\text{NRPM}\approx 3.01$| kHz, which is within the uncertainty of the fits (⁠|${\sim }12{{\ \rm per\ cent}}$|⁠). This result is in line with the interpretation of Bernuzzi et al. (2015b) and Radice et al. (2017): The postmerger f₂ frequency is mostly determined by |$\kappa _2^\text{T}$| and the merger physics. Fig. 9 (bottom panel) shows the comparison between the spectrum of this NR simulation and the respective spectrum generate with the NRPM model of Breschi et al. (2019). While the NRPM model captures well the characteristic frequencies, it does not reproduce the morphology of these high-mass-ratio waveforms due to imperfect modelling of the characteristic amplitudes and damping times. This fact further stresses the need of new simulations to improve postmerger models and/or of more agnostic approaches to kiloHertz GW modeling (cf. Breschi et al. (2019).

Figure 10.

Open in new tabDownload slide

Amplitude at merger and maximum postmerger amplitude of different modes of the GW strain computed with BLh EOS. The upper panel shows the equal-mass binary, while the lower one is a q = 1.67 case. For unequal-mass coalescences, odd-parity higher order modes are boosted, e.g. (2, 1) and (3, 3) modes.

6 DYNAMICAL EJECTA

Mass ejecta is calculated on coordinate spheres at r ≃ 300 km, assuming stationary space–time and flow, and flagging the unbound mass according to the geodesic criterion. A particle on geodesics is unbound if the four-velocity component u_t ≤ −1 and thus it reaches infinity with velocity |$v_\infty \simeq (1-u^2_t)^{1/2}$|⁠. This geodesic criterion neglects the fluid’s pressure, thus potentially underestimating the mass, but it is considered appropriate for the dynamical ejecta that are moving on ballistic trajectories (e.g. Kastaun & Galeazzi 2015).

We compute the mass histograms of the ejecta main properties and show them in Fig. 11. In the case of equal-mass BNS (top panels), dynamical ejecta is distributed all over the solid angle and composed of both the tidal and the shocked component (e.g. Hotokezaka et al. 2013; Bauswein, Goriely & Janka 2013b; Sekiguchi et al. 2015; Radice et al. 2018d). The velocity of the material peaks at v ∼ 0.1–0.2c and has high-speed tails extending up to v ∼ 0.6–0.8c. The largest tail velocities are reached by the softest EOS and in the polar regions, where baryon pollution is minimal, as a consequence of the NS cores’ bounce (fig. 3 of Radice et al. 2018d). Note, however, masses' ejecta ≪10⁻⁵ and velocities ≳0.9c cannot be trusted and can suffer from large numerical errors due to the atmosphere treatment and imperfect mass conservation (Appendix  B). The ejecta’s composition is characterized by a wide range of Y_e; for the LS220 and the BLh EOS, two peaks at ∼0.1–0.15 and at ∼0.4 are clearly visible, which roughly correspond to the shocked and tidal components, although the former has also a significant amount of material with low Y_e material. Note the SFHo model peaks instead at ∼0.25. Comparing the two equal-mass LS220 BNS, we find a small effects of turbulent viscosity: The viscous ejecta have a more prominent peak at lower Y_e and a slightly reduced tidal component, possibly due to the difference in the early-postmerger dynamics around the moment of core bounce.

Figure 11.

Open in new tabDownload slide

Distributions of the ejecta mass in the polar angle (left-hand panel), velocity (middle-hand panel), and electron fraction (right-hand panel) of the dynamical ejecta. Each row refers to a different mass ratio: from the top to bottom, q = 1, 1.67, and 1.8. Note that the angle θ = 0^○ identifies the orbital plane, while θ = 90^○ is the pole above the remnant. Data refer to resolution SR; data from simulations without turbulent viscosity are also shown.

The dynamical ejecta of asymmetric BNS with q ≳ 1.67 (middle panels of Fig. 11) are quantitatively different from symmetric BNS. The ejected material is distributed more narrowly about the orbital plane and decreases almost monotonically to the polar latitudes. The dependence on the azimuthal angle is also very different from the equal-mass cases. Because the matter is almost entirely expelled by tidal torques, the ejecta is distributed over a fraction of the azimuthal angle around its ejection angle and has a crescent shape; see Fig. 12. This is similar to what observed in BH–NS binaries (Kyutoku et al. 2015; Kawaguchi et al. 2016). Hence, the ejecta for high-q BNS are not formed isotropically. Most of the unbound mass has low Y_e ≲ 0.1, although several q = 1.67 BNS have a second peak at ≳0.4. Thus, while the tidal component is the dominant for asymmetric BNS, a small shocked component persists. The velocity distributions have comparable peak values, indicating that the tidal component has velocities comparable to those of shocked component (cf. fig. 6 of Dietrich et al. 2017). Note that the fast tails are suppressed for increasing mass ratio. This is because of the less violent merger and bounce experienced by these binaries. These features are even more extreme for the q = 1.8 case (bottom panels of Fig. 11). The above results appear consistent with those reported in Sekiguchi et al. (2015) and Lehner et al. (2016), although different EOS and more moderate mass ratio were used there.

Figure 12.

Open in new tabDownload slide

Distribution of the cumulative dynamical mass ejecta in the azimuthal angle for the BLh BNS. The ejecta of the high-q binaries expands with a crescent-like geometry similar to what found in simulations of BH–NS binaries. The mass is normalized to the total ejecta mass.

The mass-averaged properties of the dynamical ejecta computed from the histograms are reported in Table 4. We show results for all the resolutions available in order to convey an idea of the uncertainties. The latter are difficult to precisely quantify since strict convergence is not observed in the data. However, the results are robust for a large variation of the grid resolution with mass variations at the |${\sim }20{{\ \rm per\ cent}}$| level between SR and HR and less than a factor of 2 between LR and SR. Note there is a factor of 2 (1.5) between the spacing of LR and SR (SR and HR) grids. The following discussion mostly refers to highest resolutions available, as the LR is not always sufficient to properly resolve the composition.

The large mass asymmetry can boost the mass ejecta by up to a factor of 4 with respect to the equal-mass cases. The average electron fraction of the dynamical ejecta from asymmetric BNS is ∼0.11, a factor of 2 smaller than for the respective equal-mass BNS. The mass distribution is concentrated around the equatorial plane. The rms of the polar angle is ∼5–15^○ for asymmetric BNS with q = 1.8–1.67, while it is ∼35^○ for symmetric BNS. Overall, these results show that while the tidal component of the dynamical ejecta is dominant with respect to the shocked ejecta in high-mass-ratio binaries, a delayed collapse can produce unbound mass with electron fractions that can extend to Y_e ∼ 0.4. The rms of the azimuthal angle is reduced from ≲106^○ of symmetric BNS to less than half, ≳50^○, for asymmetric BNS. We recall that the rms of a uniform distribution with support on the segment 2α ∈ (0, 2π] is |$\langle \phi \rangle =\sqrt{3}/3(\pi -\alpha)$|⁠, thus giving 〈ϕ〉 ≃ 104^○ if the support is the full interval (360^○) and 〈ϕ〉 ≃ 54^○ if the support is half of the interval (180^○). A similar argument holds also for the polar angle support around the equator, π/2 − α ≤ θ ≤ π/2 + α, for which |$\langle \theta \rangle =(\sqrt{3}/3) \alpha$|⁠. This is correct as far as the ejecta is emitted uniformly over a small portion around the equator (a good approximation in the case of high-mass-ratio BNS).

As an example, we discuss mass histograms for the shocked and tidal components separately for the BLh q = 1.67; see Fig. 13. The tidal component is confined within an angle of θ ≲ 10^○ from the orbital plane; most of the mass has Y_e ∼ 0.05 with the largest electron fractions Y_e ∼ 0.15 reached at those latitudes. The velocities are uniformly distributed v ∼ 0.1c. The shocked component, instead, has mass mostly distributed at angles θ ∼ 25^○ but it extends to polar latitudes. The ejecta has electron fraction Y_e ∼ 0.17–0.25 for θ ≲ 25^○ and Y_e ∼ 0.30–0.35 for θ > 60^○. The velocity of the bulk ejecta at the orbital latitudes is v ≲ 0.25c, minimal at around θ ∼ 27^○, and has a peak v ≲ 0.3c at polar latitudes. In general, the shocked component is slightly delayed with respect to the tidal component because it is generated when the NS cores’ bounce (Radice et al. 2018d).

Figure 13.

Open in new tabDownload slide

Dynamical ejecta cumulative profiles of shocked (left-hand panel) and tidal (right-hand panel) components for BLh q = 1.67 (SR) as a function of the polar angle. The different lines are the mass, the velocity and the electron fraction. Tidal and shocked components are conventionally separated by flagging fluid elements with specific entropy smaller or larger than 10k_B per baryon, respectively. Note that the angle θ = 0^○ identifies the orbital plane, while θ = 90^○ is the pole.

Table 4 also highlights a dependence on resolution, especially for high-mass-ratio BNS. This is expected since resolving NS with different sizes is more challenging than with equal sizes for the box-in-box AMR. In particular, the LR resolutions does not seem sufficient to deliver quantitatively robust results for all the cases, especially at high q and with viscosity. Note, for example, that ejecta mass decreases with resolutions indicating numerical dissipation plays a role enhancing the ejecta. Moreover, Y_e raises very rapidly from the NS surface; in the case the latter is not well resolved, the tidal ejecta might be spuriously composed of material from the interior, as observed in the BLh q = 1.8 LR simulation.

We finally comment on the effect of viscosity on the dynamical ejecta. Radice et al. (2018c) pointed out that the dynamical ejecta in asymmetric BNS can be enhanced by the thermalization of mass accretion streams between the secondary and the primary NS. This viscous component of the dynamical ejecta is characterized by large asymptotic velocities and have masses that depend on the efficiency of the viscous mechanism. Fig. 14 shows the ejecta mass for the BLh q = 1.8 and LS220 q = 1.67 BNS. The viscous dynamical ejecta is not present because the shocked component ejecta is negligible. Actually, the turbulent viscosity here can reduce the tidal dynamical ejecta as a consequence of the different angular momentum distribution due to turbulence. Note the effect is significant and robust with respect to the variation of the grid resolution. The effect of viscosity is much reduced in the LS220 q = 1.67 BNS and practically negligible considering the numerical uncertainties (only the SR is shown for clarity). This might be related to the differences in the EOS at low density (Section 2). The simulations of Radice et al. (2018c) employed the GRLES scheme as those presented here, but using ℓ_mix = const and varying systematically the constant for the turbulent parameter. We cannot currently exclude that the specific subgrid model ℓ_mix(ρ) built from Kiuchi et al. (2018) determines a different effect with respect to the ℓ_mix = const model. A detail investigation of the viscous dynamical ejecta with the subgrid model ℓ_mix(ρ) for intermediate values of q will be presented elsewhere.

Figure 14.

Open in new tabDownload slide

Dynamical mass ejecta in viscous and non-viscous simulations. The viscous dynamical ejecta reported in Radice et al. (2018c) is here not present because the shocked component ejecta is negligible. Furthermore, the angular momentum distribution introduced by the subgrid model ℓ_mix(ρ) appears to be dependent on the EOS model. Note the turbulent viscosity subgrid model employed here is different from the ℓ_mix employed in previous simulations.

7 SYNTHETIC KILONOVA LIGHT CURVES

We compute synthetic kilonova light curves for each of the BNS mergers presented in this work. We use a semi-analytical multicomponent, anisotropic kilonova model that takes into account the angular distribution of the ejecta properties as well as the presence of different kinds of ejecta (Perego, Radice & Bernuzzi 2017; Radice et al. 2018a,d; Barbieri et al. 2020). The latter differ by the mechanisms that cause the ejection and the time-scales over which they operate. Within this framework, the homologously expanding ejecta is discretized in velocity space and the photon diffusion time is estimated by time-scale arguments. Radiation is assumed to be in local thermodynamical equilibrium up to the relevant photosphere and photon emission is modelled as a superposition of blackbody spectra. The different ejecta components comprise the dynamical ejecta discussed in Section 6 and possibly winds expelled by the remnant disc on longer time-scales (0.1–1 s) by means of neutrino irradiations and turbulent viscosity of magnetic origin. The kilonova emission produced by each component depends mainly on three quantities that characterize the ejecta, namely the amount of mass, M_ej, its average expansion velocity, 〈v_ej〉, and an (effective) grey photon opacity, κ_ej. In all our kilonova models, we locate the merging BNS at a distance of 40 Mpc and we consider a reference viewing angle of π/6 with respect to the rotational axis of the binary. If not otherwise specified, the model parameters and input physics are assumed to be as in the best-fitting model named BF to AT2017gfo of Perego et al. (2017).

We first examine the kilonova emission obtained by considering only the dynamical ejecta discussed in Section 6. In the case of q = 1 mergers, matter is expelled over the entire solid angle and we follow the model presented in Perego et al. (2017) and Radice et al. (2018a,d). We assume the ejecta to be axisymmetric and the photon diffusion to proceed mostly radially. In these cases, we discretize the polar angle in 30 slices over the whole solid angle. We use azimuthal averages of the angular distribution of the ejected mass, electron fraction, and mean expansion velocity directly extracted from the latest stages of our NR simulations. While the ejecta mass and mean velocity are directly input into the kilonova model, the electron fraction is used to assign the ejecta opacity according to κ_dyn = 1 cm² g⁻¹ for 〈Y_e〉 > 0.25 and 20 cm² g⁻¹ otherwise (cf. Kasen, Badnell & Barnes 2013; Tanaka et al. 2020; Fontes et al. 2020). Alternatively, for the q = 1.67 and 1.8 cases, the dynamical ejecta is confined (in very good approximation) within a crescent across the equatorial plane (see Section 6) and we imply the model described in Barbieri et al. (2020) (see also Kawaguchi et al. 2016) in which the photon emission is the combination of radial and lateral emissions from an optically thick disc. In this case, we use the total ejecta mass, M_ej, and mean velocity, 〈v_ej〉, obtained by our NR simulations to initialize a vertically homogeneous, radially expanding disc. For the grey opacity, we assume always κ_dyn = 20 cm² g⁻¹, since in these cases, 〈Y_e〉 < 0.25 (often 〈Y_e〉 < 0.10). For the disc half-opening angle in the polar direction, we use |$\theta _{\rm disc} = \sqrt{3} \langle \theta _{\rm ej} \rangle$|⁠, while for the azimuthal disc opening, we set |$\phi _{\rm disc} = 2 \sqrt{3} \langle \phi _{\rm ej} \rangle$| (See Section 6). The crescent shape breaks the axisymmetry of the emission. In our calculations, we always assume the dynamical ejecta to be emitted toward the observer. For small polar opening angles, this assumption is not very relevant, since the radial emission is subdominant. In the case of larger discs (as in the BLh q = 1.67 case), the radial emission can be relevant and our model assumptions can be more questionable.

In Fig. 15, we present light curves in three different photometric bands (g, z, and K_s) to span the relevant wavelength interval from visible to near-infrared radiation, for the three different models obtained with the BLh and SLy4 EOS. We first notice that, even in the case of prompt collapse, BNS mergers can power bright kilonovae (Kawaguchi, Shibata & Tanaka 2020; Kyutoku et al. 2020). In particular, in the high-q models, the light curves from the dynamical ejecta are possibly brighter, with wider light curves peaking at later times compared with the q = 1 mergers. This is due to the crescent-like configuration of the expanding ejecta (cf. Tanaka et al. 2014; Korobkin et al. 2020). On the one hand, when matter is emitted over a large portion of the solid angle (as it usually happens for q ∼ 1), the hotter ejecta is buried inside the optically thick region and high-energy photons have to diffuse and thermalize before being emitted in the kilonova. On the other hand, thanks to the disc-like geometry of the crescent, the innermost, hotter portion of the disc provides a significant contribution to the kilonova emission at any time, explaining the brighter and more substained emission. These effects are visible in all bands, but the increase in magnitude moving from q = 1 to higher qs is more pronounced in the infrared band as a consequence of the lower electron fraction (and thus of the higher opacity) of the dynamical ejecta in the crescent. This effect is even amplified by the larger amount of dynamical ejecta observed in the high-q models (with the only exception of the SFHo models).

Figure 15.

Open in new tabDownload slide

Kilonova light curves from the simulated dynamical ejecta (one component) for the BLh and SLy BNS. The ejecta properties are taken from the simulation at the highest grid resolution. The light curves are computed with the axisymmetric models of Perego et al. (2017) for q = 1 and of Barbieri et al. (2020) for q = 1.67 and 1.8. Binaries are always assumed to be located at a distance of 40 Mpc and to be observed under a viewing angle of 30^○ with respect to the BNS rotational axis. The bump observed in the K_s band for the BLh q = 1.67 model results from the radial emission from the crescent pointing towards the observer.

The peak times of all kilonova models are shown as a function of the mass ratio in Fig. 16. In addition to the models presented in Table 1, we include here also a few more LR models computed with the BLh EOS (see Appendix  B) to better explore the dependence on q. The kilonova peak times of mergers undergoing accretion-induced prompt collapse are significantly delayed with respect to the q = 1 cases. For the BLh merger, the emission in g, z, and K_s bands peaks between few hours and within a day, respectively, if q = 1, and between a day and a week if q = 1.8. The near-infrared frequencies are those that vary most as a function of the mass ratio. The SLy4 light curves show a similarly behaviour, although less data points are available. Less variation in the peak times is observed in the LS220 and SFHo mergers between q = 1 and 1.67, but note that in those cases, the dynamical ejecta mass also vary less with the mass ratio.

Figure 16.

Open in new tabDownload slide

Peak time of the one component kilonova models as a function of the mass ratio for all the simulated BNS. Data from BLh simulations at intermediate mass ratios are also included (see Appendix  B). Note the logscale and that straight lines connecting the points are due to the limited number of simulations available and do not represent accurately the functional behaviour in the mass ratio.

We tested that the features described above do not depend on the specific velocity profile for the homologously expanding ejecta, in which most of the mass resides in the innermost part of the disc. Indeed, a flat distribution in the expanding velocity as employed in Kawaguchi et al. (2016) provides very similar results. This is due to a compensation effect between the larger amount of decaying material and the denser (thus, optically thicker) vertical profile of the disc in our models. These features are robust also with respect to the uncertainties on the ejecta properties of numerical origin. Considering the ejecta properties extracted from simulations at different resolutions gives some quantitative changes that mostly affect the light curves’ luminosity. Here is worth to remark that a factor of 2 uncertainty in the ejecta mass can translate into an order of magnitude in luminosity. Moreover, current light-curve models suffer from larger systematic uncertainties in nuclear (e.g. mass models, fission fragments and β-decay rates) and atomic (e.g. detailed wavelength-dependent opacities for r-process element) physics (Eichler et al. 2015; Rosswog et al. 2017; Gaigalas et al. 2019).

The models presented in Fig. 15 do not contain potentially relevant contributions to the total ejecta coming from disc winds. Thus, the resulting light curves could be considered as lower limits for the kilonova emission. To estimate the potential impact of the disc wind emission on our results, in Fig. 17, we also present light curves obtained by considering a three-component kilonova model for the same three photometric filters and models of Fig. 15. The dynamical ejecta profiles are NR-informed as previously discussed. For the disc winds, we consider both a neutrino-driven and a viscosity-driven wind. Since wind ejection is expected over a wide portion of the solid angle, we model the related kilonova emission using again the framework described in Perego et al. (2017) and Radice et al. (2018a,d). For the neutrino-driven wind component, the amount of ejecta is assumed to be 5 per cent (1 per cent) of the disc mass if the remnant is a long-lived (short-lived or promptly collapsing) massive NS. Due to the effects of neutrino irradiation, the effective grey photon opacity is set κ_w = 1 cm² g⁻¹, while the wind expands within a π/4 angle around the polar axis with an average speed of 〈v_w〉 = 0.08c. For the viscous wind component, the amount of ejecta is always assumed to be 20 per cent of the disc mass, expanding with an average speed of 〈v_v〉 = 0.06c, while the grey opacity is set to κ_v = 5 cm² g⁻¹. To compute the masses of the wind ejecta, we consider the disc masses presented in Section 4.

Figure 17.

Open in new tabDownload slide

Kilonova light curves as in Fig. 15, but employing a three-component model for the BLh and SLy BNS. The dynamical ejecta component is taken as in Fig. 15. The other two components are assumed from a neutrino-driven and a viscosity-driven wind. The neutrino-wind mass is assumed 5 per cent (1 per cent) of the disc mass if the remnant is a long-lived (short-lived or promptly collapsing) massive NS; the effective grey photon opacity is set κ_w = 1 cm² g⁻¹, while the wind expands within a π/4 angle around the polar axis with an average speed of 〈v_w〉 = 0.08c. The viscos-wind mass is assumed 20 per cent of the disc mass, expanding with an average speed of 〈v_v〉 = 0.06c, and with a grey opacity set to κ_v = 5 cm² g⁻¹. The observational data of AT2017gfo are shown as black markers for comparison (see the discussion in the text).

Since the disc ejecta is usually more relevant than the dynamical one (see, e.g. Radice et al. 2018d), the large differences between q = 1 and high-q models in the kilonova light curves observed in the one-component models reduce for the multicomponent cases. Nevertheless, since BNS mergers with higher mass ratios tend to produce also more massive discs, also these possibly more complete models confirm that BNS mergers undergoing prompt merger can power bright kilonovae and high-qs can possibly produce kilonovae that are brighter and charaterized by wider peaks in all relevant bands, compared to more symmetric mergers mergers that have the same chirp mass. More specifically, in the case of high-q binary models for which the dynamical ejecta has a relatively large mass (up to 10⁻² M_⊙) and is highly anisotropic (e.g. BLh and SLy4 q = 1.8), the emission from the crescent is significant at all time and possibly dominant for mergers forming discs of not too large masses (M_disc ≲ 0.1 M_⊙). The opposite scenario is realized in symmetric binaries: In all q = 1 models, irrespectively of the EOS, the low-mass, widely distributed dynamical ejecta has a visible impact on the light curves only at very early times and in the blue portion of the kilonova spectrum. At later times, and especially at red and infrared frequencies, the emission is dominated by disc winds.

The observations of AT2017gfo (Villar et al. 2017 and refs20200623 therein) are also included in Fig. 17 and can be qualitatively compared to the light curves from the simulations (note that the simulated BNS have chirp mass consistent with GW170817). The light curves from high-q mergers are generically flatter and more extended in time than those of AT2017gfo. Assuming these particular light-curve models, the observation of AT2017gfo would exclude high-q and stiff EOS with |$\tilde{\Lambda }\gtrsim 600$| (long-lived NS remnants) consistently with the low-spin prior GW analysis (Abbott et al. 2019a,b). The plots also highlight that the light curves in different bands favour a different mass ratio, thus anticipating systematics (and degeneracies) between the multicomponent light curves and the binary parameters. We finally remark that the kilonova model employed here avoids the solution of the challenging radiative transfer problem in multidimension (e.g. Wollaeger et al. 2018; Bulla 2019; Kawaguchi et al. 2020) and approximates the time- and frequency-dependent r-process opacities with constant, grey opacities. This procedure likely introduces systematic uncertainties that are not easy to quantify. Direct comparisons between simpler analytical models and the outcome of radiative transfer calculations indicate that the former tend to predict lower luminosities and later peaks, especially for κ ≳ 100 cm² g⁻¹ (Wollaeger et al. 2018). The usage of input parameters gauged on AT2017gfo and of opacities ≲25 cm² g⁻¹ possibly limits these uncertainties. We conservatively estimate a residual uncertainty of ±0.5 mag at peak. Even including these uncertainties, the qualitative differences between AT2017gfo and the light curves obtained for high-q and stiff EOSs still hold.

8 CONCLUSIONS

In this paper, we explored in a systematic way the dynamics, the ejecta, and the expected kilonova light curves of highly asymmetric BNS mergers by means of detailed simulations in NR. The latter employed different finite-temperature, composition-dependent EOS, and numerical resolutions. The prompt-collapse dynamics discussed here for high-q BNS have an underlying mechanisms different from the equal-mass prompt collapse: In the former case, the collapse is driven by the accretion of the companion on to the massive primary star. For binaries with increasing mass ratio and fixed chirp mass, the companion NS undergoes a progressively more significant tidal disruption. Thus, in these BNS sequences, accretion-induced prompt collapse should be always present after a critical mass ratio in connection to the maximum NS mass. For example, for the BLh EOS, the critical mass ratio should fall in the interval 1.54 < q_thr < 1.67.

The remnant BH in these high-mass-ratio mergers is surrounded by a massive accretion disc in contrast to comparable-mass prompt-collapse merger that have no significant disc left outside the BH. The accretion discs of high-mass-ratio mergers are primarily constituted of tidally ejected material, hence they are initially cold and neutron-rich. The simulations show that fallback of the tidal tail perturbs the disc and affect its accretion. The long-term disc and fallback dynamics are relevant to understand the complete kilonova emission and also for GRB afterglow (extended) emission (Rosswog 2007; Metzger et al. 2010; Desai, Metzger & Foucart 2019). This study is left for the future.

Perhaps the most relevant astrophysical consequence of our work is the possibility of having massive dynamical ejecta from these accretion-induced prompt collapsing remnants. The ejecta mass can reach M_ej ∼ 0.007–0.01 M_⊙ and are mostly emitted within 10^○–20^○ about the orbital plane and in a portion of 100^○–180^○ in the azimuthal angle. The ejecta is neutron-rich with Y_e ≲ 0.1 and with velocities v ≲ 0.1c. The related kilonova light curves are predicted to be usually significantly brighter than the equal-mass case (at fixed chirp mas) in all the bands as a consequence of the crescent-like geometry of the expanding dynamical ejecta. The light curves peak at later times and are powered by the sustained emission of the innermost, hotter portion of the crescent especially in the infrared bands.

We suggest that the confident detection (or confident non-detection) of an electromagnetic counterpart for a high-mass binary can directly inform us about the binary mass ratio. Because the latter is currently poorly constrained by GW analysis, the kilonova counterpart can deliver significant complementary information. Multimessenger analysis of high-mass events are thus particularly relevant. They will require a precise NR characterization of the ejecta in terms of the binary parameters that is not currently available, as well as improved nuclear and atomic physics input or suitably parametrized models for the light curves.

Our results can help interpreting GW190425 in the scenario that the GW was produced by an asymmetric binary with q ≳ 1.6. (Note the chirp mass for GW190425 is even larger than the one simulated here, while large mass ratios are excluded for GW170817 if spins are small). Using the methods developed in Agathos et al. (2020), Abbott et al. (2020) estimated that the probability for the remnant to prompt collapsed to BH is |${\sim }97{{\ \rm per\ cent}}$|⁠. The NR fitting models used in Agathos et al. (2020) refer to equal masses and thus are to be considered conservative for q ≳ 1.6. Hence, if GW190425 is interpreted as a such asymmetric BNS merger, the BH remnant scenario is further strengthened by our results. Moreover, a bright and temporally extended red kilonova could have been expected as a counterpart if GW190425 was produced by a high-q merger (cf. Foley et al. 2020). The kilonova signal in this case could be similar to the one produced in BH–NS binaries (Radice et al. 2018a; Kyutoku et al. 2020).

ACKNOWLEDGEMENTS

SB, MB, BD, NO, and FZ acknowledge the support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. Numerical relativity simulations were performed on the supercomputer SuperMUC at the LRZ Munich (Gauss project pn56zo), supercomputer Marconi at CINECA (ISCRA-B project number HP10BMHFQQ), supercomputers Bridges, Comet, and Stampede (NSF XSEDE allocation TG-PHY160025), NSF/NCSA Blue Waters (NSF AWD-1811236), and ARA cluster at Jena FSU. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231. Data postprocessing was performed on the Virgo ‘Tullio’ server at Torino, supported by INFN. The authors gratefully acknowledge the Gauss Centre for Supercomputing eV (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de).

DATA AVAILABILITY

All of our GW waveforms and ejecta data will be publicly available as part of the CoRe data base at http://www.computational-relativity.org/.

Footnotes

REFERENCES