Numerical relativity simulations of the neutron star merger GW190425: microphysics and mass ratio effects

ABSTRACT

GW190425 was the second gravitational wave (GW) signal compatible with a binary neutron star (BNS) merger detected by the Advanced LIGO and Advanced Virgo detectors. Since no electromagnetic counterpart was identified, whether the associated kilonova was too dim or the localization area too broad is still an open question. We simulate 28 BNS mergers with the chirp mass of GW190425 and mass ratio 1 ≤ q ≤ 1.67, using numerical-relativity simulations with finite-temperature, composition dependent equations of state (EOS) and neutrino radiation. The energy emitted in GWs is |$\lesssim 0.083\mathrm{\, M_\odot }c^2$| with peak luminosity of 1.1–|$2.4\times ~10^{58}/(1+q)^2\, {\rm {erg \, s^{-1}}}$|⁠. Dynamical ejecta and disc mass range between 5 × 10⁻⁶–10⁻³ and 10⁻⁵–|$0.1 \mathrm{\, M_\odot }$|⁠, respectively. Asymmetric mergers, especially with stiff EOSs, unbind more matter and form heavier discs compared to equal mass binaries. The angular momentum of the disc is 8–|$10\mathrm{\, M_\odot }~GM_{\rm {disc}}/c$| over three orders of magnitude in M_disc. While the nucleosynthesis shows no peculiarity, the simulated kilonovae are relatively dim compared with GW170817. For distances compatible with GW190425, AB magnitudes are always dimmer than ∼20 mag for the B, r, and K bands, with brighter kilonovae associated to more asymmetric binaries and stiffer EOSs. We suggest that, even assuming a good coverage of GW190425’s sky location, the kilonova could hardly have been detected by present wide-field surveys and no firm constraints on the binary parameters or EOS can be argued from the lack of the detection.

gravitational waves, hydrodynamics, nuclear reactions, nucleosynthesis, abundances, methods: numerical

1 INTRODUCTION

The advent of the network of terrestrial gravitational wave (GW) detectors formed by Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015), recently joined also by KAGRA (Aso et al. 2013; Akutsu et al. 2019), has opened the era of GW astronomy. At the end of the third observing run, the GW emission resulting from the late inspiral or from the merger of two black holes (BHs), one BH and a neutron star (NS), or two NSs have all been observed (Abbott et al. 2019b, 2021a, c).

So far, two GW signals compatible with the inspiral of a BNS system were reported: GW170817 (Abbott et al. 2017a) and GW190425 (Abbott et al. 2020). GW170817 was interpreted as the merger of a BNS system with a chirp mass |$\mathcal {M}_{\rm chirp}=(1.186 \pm 0.001) \mathrm{\, M_\odot }$|⁠. The masses of the individual stars were |$M_A = (1.46^{+0.12}_{-0.10}) \mathrm{\, M_\odot }$| and |$M_B = (1.27^{+0.09}_{-0.09}) \mathrm{\, M_\odot }$|⁠, at 90 per cent credible level, resulting in a total mass in the range |$2.72\!-\!2.76 \mathrm{\, M_\odot }$| (Abbott et al. 2017a, 2019a). The total mass of such a system is thus well within the expected range of Galactic BNS systems, as resulting from electromagnetic (EM) observations of pulsars in BNS systems (see e.g. Özel & Freire 2016). The NS nature of the colliding objects was further corroborated by the detection of several EM counterparts originated from a galaxy located at 40Mpc from us, including a short gamma-ray burst and its afterglow, a kilonova, and possibly the non-thermal emission produced by the high speed tail of the dynamical ejecta expelled in the merger (see e.g. Radice 2020; Margutti & Chornock 2021, and references therein). The possibility of detecting GW170817 counterparts crucially depended on the availability of three detectors, which drastically reduced the sky localization area to 16 deg² (Abbott et al. 2017a, b, 2021c).

GW190425 represented a significantly different event with respect to GW170817 in many aspects (Abbott et al. 2020, 2021a). The rest-frame chirp mass was |$(1.44 \, \pm \, 0.02) \,\, {\rm M_{\odot }}$|⁠, while the NS mass ranges were |$M_A = (2.0^{+0.6}_{-0.3}) \mathrm{\, M_\odot }$| and |$M_B = (1.4^{+0.3}_{-0.3}) \mathrm{\, M_\odot }$|⁠, at 90 per cent credible level, resulting in a total mass in the range |$3.3\!-\!3.7 \mathrm{\, M_\odot }$|⁠. Such a high total mass qualifies GW190425 as a possible outlier in the Galactic BNS system distribution (Abbott et al. 2020, 2021b).

During the passage of the GW signal, the Livingston LIGO detector was offline and Virgo was unable to contribute to the measure because of the small signal-to-noise ratio (2.5) resulting from the large inferred distance (D ≈ 70–250 Mpc). The effective presence of only one GW detector did not allow a good sky localization (∼10⁴ deg²).

Despite an intense followup campaign within the first days after the gravitational wave (GW) detection, no firm identification of EM counterparts was possible so far (see e.g. Coughlin et al. 2019; Steeghs et al. 2019). In particular, the GROWTH and GRANDMA collaborations performed dedicated follow-up campaigns. GROWTH made use of the Zwicky Transient Facility (ZTF) and the Palomar Gattini-IR telescopes. The ZTF system covered 21 per cent of the probability integrated skymap and achieved a depth of 21 AB magnitudes in the g- and r-bands, while Palomar Gattini-IR covered 19 per cent of the probability integrated skymap in J-band to a depth of 15.5 mag (Coughlin et al. 2019). With 9 of its 21 heterogeneous telescopes, the GRANDMA network imaged 70 galaxies covering ≲ 2 per cent of the probability integrated skymap, attaining a depth of 17-23 AB magnitudes depending on the telescope (Antier et al. 2020). In absence of an optical or infrared counterpart, Apertif-WSRT searched for afterglow radio emission in a 9.5 deg² region of the high probability skymap (Boersma et al. 2021). Despite the reduced fraction of the covered skymap, the apparent lack of EM counterparts and the unusually high total mass of the binary leave open questions both on the origin of the system and on the remnant properties.

Numerical modelling of BNS mergers is a necessary step to properly interpret results, address open questions, and extract the largest amount of information from available data, even from the potential lack of detections. In particular, simulations of the inspiral, merger, and early merger aftermath allow to extract the GW signal, the properties of the so-called dynamical ejecta, and the properties of the merger remnant (see e.g. Baiotti & Rezzolla 2017; Shibata & Hotokezaka 2019; Radice, Bernuzzi & Perego 2020; Bernuzzi 2020, for recent reviews). GW170817 was the privileged target of several simulation campaigns in numerical relativity (see e.g. Nedora et al. 2021b). Recently, an independent study on GW190425 in numerical relativity has been proposed in Dudi et al. (2021; hereafter Dudi et al.). The authors set up 36 BNS simulations targeted to GW190425 considering four mass ratios and three nuclear EOSs at different resolutions. They used cold EOSs with a density dependent composition fixed by neutrino-less beta-equilibrium conditions, and with thermal effects included by an effective Γ-law. Dudi et al. compute kilonova light curves employing a wavelength-dependent radiative transfer code (Kawaguchi, Shibata & Tanaka 2020), for which the post-merger ejecta composition is fixed for all components. They concluded that, assuming an effective coverage of the event localization region in the GROWTH follow-up campaign, the lack of kilonova detection suggests that GW190425 is incompatible with a face-on, unequal BNS merger with more than 20 per cent of mass difference between the two NSs. In all other cases (soft EOSs, edge-on and more distant mergers, or more symmetric binaries) the lack of detection is still compatible with a fainter kilonova signal.

Several other works focused on GW190425 have recently appeared. For example Han et al. (2020) and Kyutoku et al. (2020) investigated the possibility that GW190425 originated from a BH-NS merger by studying the corresponding GW and kilonova signal, respectively.

In Raaijmakers et al. (2021) and Barbieri et al. (2021) kilonova light curves for GW190425 were computed under the assumption that the originating event was a BH-NS or a BNS merger.¹ In both cases, the properties of the ejecta powering the kilonova signal were computed using fitting formulae derived from broad simulation samples, while the kilonova signals were computed using models with different levels of sophistication. In Barbieri et al. (2021), the BNS fitting formulae were taken from Radice et al. (2018b) and from the appendix of Barbieri et al. (2021). The NS masses were chosen to be compatible with the GW190425 chirp mass, while the two employed NS EOSs were compatible with present nuclear and astrophysical constraints. Additionally, using the same model, they also computed light curves directly using GW190425 posteriors (Abbott et al. 2020). They concluded that a light BH in GW190425 would have produced a brighter kilonova emission compared to BNS case, allowing to distinguish the nature of the binary. However also in the BNS case, the merger could have produced kilonovae bright enough to have been possibly detected by ZTF, especially for stiff EOSs and for more asymmetric systems.

In Raaijmakers et al. (2021), only the posteriors from GW190425 (Abbott et al. 2020) and the EOS obtained from GW170817 analysis (Abbott et al. 2018) were used as input for the BNS fitting formulae from Krüger & Foucart (2020) and Foucart et al. (2017). Based on the obtained ejecta and disc properties, kilonova light curves were computed using the semi-analytic model from Hotokezaka & Nakar (2019). The latter adopts the radioactive heating rate fit from Korobkin et al. (2012) and assumes a spherical symmetry for the ejecta geometry. Additionally, tests using the same kilonova model but fitting formulae from Radice et al. (2018b), Barbieri et al. (2021), and Dietrich, Hinderer & Samajdar (2021) were also performed.

Despite these works, several open questions regarding GW190425 still remain. For example, how robust are the results obtained in numerical relativity for GW190425-like events? And, in particular, what is the impact of input physics that was so far neglected in GW190425-targeted simulations, including finite temperature, composition dependent EOS, and neutrino radiation? What are the detailed properties of the dynamical ejecta expelled in these events and how do they depend on the binary properties and on the NS EOS? Is there a characteristic nucleosynthesis signature in these ejecta? Based on these results, what can we infer from the missing detection of electromagnetic counterparts for GW190425?

To answer these questions, we setup 28 simulations in numerical relativity targeted to GW190425 with finite temperature, composition dependent NS EOSs, and with neutrino radiation. We investigate the binary evolution up to the first ≈10 ms after merger. We extract both remnant and dynamical ejecta properties, to give credible answers to some of the above questions. In particular, we use the detailed outcome of our simulations to compute nucleosynthesis yields and to set up kilonova models. We found that, for a distance compatible with GW190425, only in the case of a very stiff EOS and a very asymmetric binary the resulting kilonova could have been bright enough to be observed by the ZTF facility. This suggests that the possible lack of kilonova counterpart for GW190425 provides much weaker constraints than previously thought.

The paper is structured as follows: after a brief recap of the numerical setup and of the simulations properties in Section 2, we resume the qualitative behaviour of the merger dynamics in Section 3.1 and analyse the GW energetics in Section 3.2. The quantitative description of the remnant is reported in Section 3.3, while we discuss the main properties of the dynamical ejecta in Section 3.4. In Sections 4.1 and 4.2 we describe the output from the nucleosynthesis process and its related kilonova signal. We compare our results with the one discussed in the literature in Section 5. We summarize our results in the conclusions in Section 6.

2 METHODS AND MODELS

2.1 Binary merger calculations

We evolve BNS systems in full general relativity (GR) through 3 + 1 numerical relativity simulations encompassing the latest orbits, the merger, and the early post-merger phase.

The spacetime metric is evolved with the Z4c formulation of Einstein’s equations (Bernuzzi & Hilditch 2010; Hilditch et al. 2013) using the ctgamma code (Pollney et al. 2011; Reisswig et al. 2013b), developed within the EinsteinToolkit framework (Loffler et al. 2012; Brandt et al. 2021).

We use the WhiskyTHC code (Radice & Rezzolla 2012; Radice, Rezzolla & Galeazzi 2014), implemented within the Cactus (Goodale et al. 2003; Schnetter et al. 2007) framework to solve the GR hydrodynamic equations. WhiskyTHC evolves the proton and neutron number density equations, in addition to the relativistic version of the momentum and energy conservation equations, written in conservative form.

To properly resolve the NS structure and merger dynamics, and at the same time track the evolution of the ejecta on a large enough domain, we employ a mesh refinement (Schnetter, Hawley & Hawke 2004; Reisswig et al. 2013a) consisting in seven nested grids characterized by a 1:2 linear scaling between consecutive grids, with the most refined level covering the two NSs during the inspiral and the central remnant after merger. We characterize each simulation by the resolution of the innermost grid, h, and in particular h ≈ 246 m for low resolution (LR) and |$h \approx 185 \,\, \rm m$| for standard resolution (SR) runs. Once the symmetry along the z = 0 plane is taken into account, the simulated space is a cube of side |$3024 \,\, \rm km$|⁠. For further details on the numerical setup we refer to Radice et al. (2018b).

Thanks to the use of a puncture gauge, the spacetime evolution can handle the formation of a singularity within the computational domain (Thierfelder et al. 2011; Dietrich & Bernuzzi 2015). The apparent horizon (AH) can possibly be detected by the AHFinderDirect thorn (Thornburg 2004) of the EinsteinToolkit, from which the BH properties can be extracted.

In all simulations we include compositional and energy changes due to the emission and absorption of neutrinos of all flavours. In particular, a grey leakage scheme (Ruffert, Janka & Schäfer 1996; Galeazzi et al. 2013; Neilsen et al. 2014) is used to model the net neutrino emission rates both from optically thick regions, where neutrinos are expected to form a diffusing gas in thermal and weak equilibrium with matter, and optically thin regions. Neutrinos are then transported by an M0 scheme (Radice et al. 2018b) through optically thin regions, where the reabsorption of streaming electron flavours (anti)neutrinos can happen in addition to local emission.

We use four finite-temperature, composition-dependent EOSs compatible with current astrophysical (Cromartie et al. 2019; Miller et al. 2019; Riley et al. 2019) and nuclear (Capano et al. 2020; Jiang et al. 2020) constraints: BLh (Bombaci & Logoteta 2018; Logoteta, Perego & Bombaci 2021), HS(DD2) (Hempel & Schaffner-Bielich 2010; Typel et al. 2010), SFHo (Steiner, Hempel & Fischer 2013) and SRO(SLy4) (Douchin & Haensel 2001; Schneider, Roberts & Ott 2017). In the following, we will refer to the second and fourth ones simply as DD2 and SLy4. All these EOSs include neutrons, protons, nuclei, electrons, positrons, and photons as relevant thermodynamics degrees of freedom, and assume baryon matter in nuclear statistical equilibrium.

The BLh EOS Logoteta et al. (2021) is an extension of the zero-temperature BL EOS Bombaci & Logoteta (2018) that includes finite-temperature effects and arbitrary particle composition. It was obtained within the finite-temperature version of the Brueckner–Bethe–Goldstone quantum many-body theory in the Brueckner–Hartree–Fock approximation. The underlying two-body and three-body interactions have been derived in chiral perturbation theory taking into account the effect of nucleon–nucleon and nucleon–nucleon-nucleon interactions.

DD2 and SFHo were computed in the framework of relativistic mean field theories. The two EOSs differ because of the different parameterizations and coupling constants for modelling the mean-field nuclear interactions. The transition to inhomogeneous nuclear matter was done using an excluded volume approach.

The SLy4 used in the present work is the finite temperature extension of the Skyrme effective nuclear interactions introduced in Douchin & Haensel (2001). The SLy4 EOS reproduces well empirical saturation properties of nuclear matter as well as several observables deduced from the mass of finite nuclei.

In Fig. 1 we show the mass-radius, the mass-central density and the mass-quadrupolar tidal polarizability relations computed for the equilibrium NS sequences for the different EOSs used in this work. The quadrupolar tidal polarizability is computed as Λ = (2/3)k₂C⁻⁵, where k₂ is the dimensionless quadrupolar Love numbers (Damour 1983; Hinderer 2008), and C the stellar compactness C = GM/(c²R). The SLy EOS produces the most compact NSs, while NSs modelled with the DD2 EOS have the largest radii around 13 km for |$1 \mathrm{\, M_\odot }\lesssim M \lesssim 2.1 \mathrm{\, M_\odot }$|⁠.

Figure 1.

TOV sequences for the NS EOSs used in this work. Left-hand panel: Gravitational mass versus radius. Central panel: Gravitational mass versus central density normalized to the nuclear saturation density, ρ₀ = 2.67 × 10¹⁴ g cm⁻³. Right-hand panel: gravitational mass versus tidal polarizability Λ. The different markers refer to the different mass ratios of the binaries evolved in the simulations.

Initial data for our simulations are constructed using the pseudo spectral elliptic solver Lorene (Gourgoulhon et al. 2001), using the EOS slice at the lowest available temperature and assuming neutrino-less beta-equilibrium. All simulations are initialized as irrotational binaries on quasicirular orbits of coordinate radius 45 km. The residual initial eccentricity, estimated following Kyutoku, Shibata & Taniguchi (2014), is between 0.02 and 0.06 for all models.

We set up and analyse a total of 28 simulations, 14 at SR and 14 at LR. In Table 1 we report a summary of all initial parameters characterizing our simulations, in particular: the values of the individual stellar masses M_{A, B} with M_A > M_B, the total gravitational mass M, the mass ratio q ≡ M_A/M_B > 1, the total ADM mass, and angular momentum of the system M_ADM and J_ADM, the stellar compactness C_i for i = A, B, the the tidal deformability of the binary, |$\tilde{\Lambda }$|⁠, defined as

$$\begin{eqnarray} \tilde{\Lambda }= \frac{16}{13} \frac{(M_A + 12M_B) M_A^4}{M^5} \, \Lambda _A + A \leftrightarrow B \, , \end{eqnarray}$$

(1)

and the coefficients |$k_2^{\rm L}$| as defined in equation (4) of Zappa et al. (2018), namely,

$$\begin{eqnarray} \kappa _2^L = 6 \left[ \frac{(3 M_B + M_A) M_A^4}{M^5} \Lambda _A + A \leftrightarrow B \right] \, , \end{eqnarray}$$

(2)

where the notation (A↔B) indicates a second term identical to the first except that the indices A and B are exchanged. We also report the GW initial frequency f_GW(0) measured in Hertz. All BNS parameters are compatible with the ones inferred from the GW signal GW190425 (Abbott et al. 2020) using both the low- and high-spin priors, except for the ones characterized by q = 1.33 and q = 1.67, which are compatible only with high-spin prior.

Table 1.

NS initial properties grouped by EOS. From left to right: EOS, maximum Tolman–Oppenheimer–Volkoff (TOV) mass |$M_{\rm TOV}^{\rm max}$|⁠, maximum TOV compactness |$C_{\rm TOV}^{\rm max}$|⁠, NS masses M_A, M_B, total gravitational mass M, BNS mass ratio q ≡ M_A/M_B, compactness of the two NSs C_A, C_B, tidal deformability of the BNS |$\tilde{\Lambda }$| defined in equation (1), the coefficient |$k_2^{\rm L}$| defined in equation (4) of Zappa et al. (2018), equation (2), the initial GW frequency f_GW(0), the total ADM mass of the system M_ADM and the initial ADM angular momentum J_ADM.

EOS	|$M_\text{TOV}^\text{max}$|	|$C_\text{TOV}^\text{max}$|	MA	MB	M	q	CA	CB	|$\tilde{\Lambda }$|	|$\kappa _2^{\rm L}$|	fGW(0)	MADM	JADM
	|$[\mathrm{\, M_\odot }]$|		|$[\mathrm{\, M_\odot }]$|	|$[\mathrm{\, M_\odot }]$|	|$[\mathrm{\, M_\odot }]$|						[Hz]	|$[\mathrm{\, M_\odot }]$|	|$[\mathrm{\, M_\odot }^2]$|
BLh	2.103	0.299	1.654	1.654	3.308	1.0	0.201	0.201	129.525	194.3	608	3.272	10.23
BLh	2.103	0.299	1.750	1.557	3.307	1.12	0.215	0.187	133.008	198.6	603	3.271	10.19
BLh	2.103	0.299	1.795	1.527	3.322	1.18	0.222	0.183	131.172	195.0	609	3.284	10.23
BLh	2.103	0.299	1.914	1.437	3.351	1.33	0.242	0.172	134.612	196.8	611	3.313	10.24
DD2	2.420	0.300	1.654	1.654	3.308	1.0	0.184	0.184	257.963	386.9	608	3.270	10.23
DD2	2.420	0.300	1.795	1.527	3.322	1.18	0.200	0.170	256.534	382.8	609	3.285	10.24
DD2	2.420	0.300	1.914	1.437	3.351	1.33	0.214	0.160	254.057	375.1	611	3.312	10.24
DD2	2.420	0.300	2.149	1.289	3.438	1.67	0.244	0.144	247.763	354.8	616	3.400	10.25
SFHo	2.059	0.294	1.654	1.654	3.308	1.0	0.209	0.209	101.708	152.6	608	3.275	10.25
SFHo	2.059	0.294	1.795	1.527	3.322	1.18	0.230	0.191	102.689	152.7	609	3.290	10.26
SFHo	2.059	0.294	1.914	1.437	3.351	1.33	0.251	0.179	104.653	153.0	611	3.320	10.28
SLy4	2.055	0.303	1.654	1.654	3.308	1.0	0.212	0.212	89.251	133.9	608	3.271	10.23
SLy4	2.055	0.303	1.795	1.527	3.322	1.18	0.234	0.194	90.538	134.6	609	3.285	10.24
SLy4	2.055	0.303	1.914	1.437	3.351	1.33	0.256	0.181	93.140	136.0	611	3.314	10.25

To better characterize the binaries used in this work and their properties in relation to the different EOSs, in Fig. 1 we also highlight the properties of the NSs initially forming the binaries evolved by our simulations. Note that the initial conditions span a broad range of central densities, from 2.2ρ₀ to 6.0ρ₀ (in terms of the nuclear saturation density ρ₀ = 2.67 × 10¹⁴ g cm⁻³) depending on the EOS and mass ratio. For the more asymmetric binaries, the central density of the heaviest NS is roughly 1.5 times larger than the one of the lightest NS, while in the equal mass case the two identical NSs have a central density ∼1.2 times larger than the one of the lightest NS in our sample. The single star tidal polarizability varies between two orders of magnitudes and, again, to asymmetric BNS corresponds two NSs with rather different tidal polarizability: a more compact and less deformable NS along with a larger and more deformable one. Interestingly, |$\tilde{\Lambda }$| varies only by a few per cents within the same EOS, while it changes by almost a factor of three between the SLy4 and the DD2 EOS.

2.2 GWs and remnant properties

We analyse the GW signal of the BNS mergers as extracted at a coordinate radius of |$\approx 591 \,\, \rm km$| from the BNS centre of mass for all the simulations in the present work. We simulate the last 3–4 orbits before merger. The latter is defined as the moment in retarded time at which the amplitude of the l = m = 2 mode of the GW waveform reaches its maximum. The short inspiral phase and the prompt collapse of the remnant to a BH do not permit to test in detail inspiral-merger-post-merger waveform models. Instead, we focus on the characterization of the GW emission during the inspiral, merger and post-merger phases through integrated and peak quantities.

In particular, we define the rescaled total energy radiated in GWs, |$e_{\rm GW}^{\rm tot}$|⁠, and the rescaled angular momentum of the remnant, j_rem, as

$$\begin{eqnarray} e_{\rm GW}^{\rm tot} = \frac{(M - M_{\rm ADM})c^2 + E_{\rm GW}^{\rm rad}}{\nu M c^2} \, , \end{eqnarray}$$

(3)

and

$$\begin{eqnarray} j_{\rm rem} = \frac{J_{\rm ADM} - J_{\rm GW}^{\rm rad}}{\nu G M^2/c} \, , \end{eqnarray}$$

(4)

where |$E_{\rm GW}^{\rm rad}$| and |$J_{\rm GW}^{\rm rad}$| are the energy and angular momentum radiated in GWs during the whole simulation, and ν is the symmetric mass-ratio, ν = M_AM_B/M².

Our remnants are characterized by the presence of a central BH surrounded by an accretion disc. We extract the properties of both from our simulations.

In particular, we define the disc as the portion of the remnant outside the apparent horizon whose rest mass density is smaller than 10¹³ g cm⁻³, (see e.g. Shibata et al. 2017). Moreover, we express the mass of the BH as

$$\begin{eqnarray} M_{\rm BH}^2 = M_{\rm irr}^2 + \left(\frac{c J_{\rm BH}}{2 G M_{\rm irr}} \right)^2 \, , \end{eqnarray}$$

(5)

where M_BH and J_BH are the gravitational mass and spin of the BH, respectively, while M_irr is the irreducible BH mass

$$\begin{eqnarray} M_{\rm irr} = \frac{c^2}{G} \sqrt{\frac{A_{\rm H}}{16 \pi }} \, , \end{eqnarray}$$

(6)

with A_H the AH area. For a Kerr-BH, the irreducible mass is a non-decreasing quantity and it coincides with the gravitational mass for non rotating BHs. In analogy with the Kerr solution, we define the dimensionless spin parameter as |$a_{\rm BH} \equiv (c J_{\rm BH}) / (G M_{\rm BH}^2)$|⁠.

The AH finder is able to give an estimate of such quantities by locating the AH of the singularity, albeit it is not guaranteed that it does locate the AH with sufficient accuracy. This issue can clearly have an impact on the estimated BH properties. We compare the gravitational mass provided by the AH finder with the expected BH mass

$$\begin{eqnarray} M_{\rm BH}^{\text{exp}} = M_{{\rm ADM}} - M_{\rm disc} - E_{\rm GW}^{\rm rad}/c^2 \, , \end{eqnarray}$$

(7)

where |$E_{\rm GW}^{\rm rad}$| is the total energy radiated in GWs. In the above expression, we have neglected the ejecta mass and for the disc we have considered only the rest-mass energy. Similarly, for the spin parameter we compute the expected value as:

$$\begin{eqnarray} a_{\rm BH}^{\rm exp} = \frac{c J_{\rm BH}^{\rm exp}}{G \left(M_{\rm BH}^{\rm exp} \right)^2} = \frac{c (J_{{\rm ADM}} - J_{\rm GW}^{\rm rad} - J_{\rm disc})}{G \left(M_{\rm BH}^{\rm exp} \right)^2} \, , \end{eqnarray}$$

(8)

where |$J_{\rm GW}^{\rm rad}$| is the angular momentum radiated in GWs and J_disc is the angular momentum of the surrounding disc.

2.3 Ejecta and nucleosynthesis calculations

From each simulation we consider the dynamical ejecta as the matter that becomes unbound within the end of the simulation on the basis of the geodesic criterion, i.e. when |u_t| ≥ c, where u_t is the time-component of the four-velocity. The properties of the ejecta are determined as matter crosses a spherical detector of coordinate radius |$r_{\rm E} = 200 G\mathrm{\, M_\odot }/c^2 \approx 294\ \textrm {km}$|⁠, discretized in N_θ = 51 polar and N_ϕ = 93 azimuthal uniform angular bins. For the unbound matter, the speed reached at infinity is computed as |$v_\infty = c \sqrt{1 - (c/u_t)^2}$|⁠.

The distribution of nuclei within the expanding ejecta is computed using the same approach and the same input data as the ones reported in Perego et al. (2022), that we briefly summarize in the following. We note that a similar approach was already used in Radice et al. (2016, 2018b), Nedora et al. (2021b), but with different input data. To obtain time-dependent yield abundances we employ SkyNet (Lippuner & Roberts 2017), a publicly available nuclear network which computes the nucleosynthesis depending on the evolution of a given Lagrangian fluid element. We evolve several trajectories with different initial parameters, with the aim of modelling the long-term expansion of the unbound matter measured in the simulations at the detector. All the trajectories start in nuclear statistical equilibrium (NSE) from an initial temperature of T₀ = 6.0 GK. The corresponding initial density, |$\rho _0 \equiv \rho (s,Y_e,T=6\, \textrm {GK})$|⁠, is determined by the NSE solver implemented in SkyNet depending on the initial values of the electron fraction Y_e and of the specific entropy s. The subsequent evolution of the density is set by the expansion time-scale τ, first as an exponentially decaying phase and then as a homologous expansion:

$$\begin{eqnarray} \rho (t) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho _0\, e^{-t/\tau } & \textrm {if}\ t\le 3\tau \, ,\\ \rho _0\bigg (\dfrac{3\tau }{et}\bigg)^3 & \textrm {if}\ t\gt 3\tau \, . \end{array}\right. \end{eqnarray}$$

(9)

Parametric nucleosynthesis calculations are repeated for a set of fluid elements characterized by different values of s, τ and Y_e, ranging on a 26 × 18 × 25 regular grid that spans the typical conditions characterizing the ejecta in compact binary mergers, i.e. 1.5 ≤ s [k_B baryon⁻¹] ≤ 300, 0.5 ≤ τ [ms] ≤ 200 and 0.01 ≤ Y_e ≤ 0.48, approximately logarithmic in the two former parameters while linear in the latter. To compute the nucleosynthetic yields in the ejecta we take the convolution of the output given by SkyNet with the distribution of the ejecta properties extracted from the numerical simulation at r_E. While s and Y_e are directly extracted from the numerical simulation, τ is computed following the procedure described in Radice et al. (2016, 2018b).

2.4 Kilonova light curves calculations

In order to compute kilonova light curves from the outcome of our simulations, we employ the multi-component anisotropic framework presented in Perego, Radice & Bernuzzi (2017). In this framework, axial symmetry and symmetry with respect to the BNS orbital plane are assumed, while the polar angle θ is discretized in N_θ = 30 angular bins equally spaced in cos θ. The kilonova emission is then computed in a ray-by-ray fashion by summing up the photon fluxes coming from each angular slice, properly projected along the line of sight of an observer located at a polar angle θ_view. Inside each slice, a 1D kilonova model is used. The latter depends on the mass and (root mean square) speed of the ejecta, as well as on an effective grey opacity κ. Inside each ray, several ejecta components are considered, resulting from the expulsion of matter operated by different mechanisms, acting on different time-scales and providing distinct ejecta properties. The total luminosity is found by summing over the contributions of the different ejecta components, assuming that the energy emitted by the innermost ones is quickly reprocessed and emitted by the outermost component.²

Differently from the model originally implemented in Perego et al. (2017) and later employed, for example, in Radice et al. (2018a, b), Breschi et al. (2021), Barbieri et al. (2019, 2020, 2021), here we adopt a new semi-analytical 1D kilonova model for each angular slice that we present in the following. The model assumes a spherically symmetric and optically thick outflow with a constant average grey opacity. The outflow expands with an homologous expansion law, i.e. the density of each fluid element decreases as t⁻³ while its expansion speed stays constant, starting from a few hours after merger. The kilonova emission is calculated as the combination of two contributions, one emitted at the photosphere and one coming from the optically thin layers above it.

The contribution coming from the photosphere is computed starting from the semi-analytic formula for the luminosity originally proposed by Wollaeger et al. (2018) and derived from a solution of the radiative transfer equation in the diffusion approximation (Pinto & Eastman 2000). This formula was further validated in Wu et al. (2022), where it showed a very reasonable agreement with results provided by the radiation hydrodynamics code snec. While the original model assumes that the whole ejecta are in optically thick conditions, an increasing fraction of it resides outside of the photosphere, becoming optically thin to thermal radiation. For this reason, the outcome of this computation is rescaled by a factor M_thick/M_ej, where M_thick is the mass of the optically thick part of the ejecta, defined as the region enclosed by the photosphere. The photospheric radius R_ph(t) is found analytically by imposing the condition τ_γ(R_ph) = 2/3, where τ_γ is the optical depth of the material, and by using the homologous density profile as in Wollaeger et al. (2018):

$$\begin{eqnarray} \rho (t,x)=\rho _0\left(\frac{t_0}{t}\right)^3\left(1-x^2\right)^3, \end{eqnarray}$$

(10)

where ρ₀ is the density at the initial time t₀ and x = v/v_max is the dimensionless radial variable. The photospheric temperature T_ph(t) is computed from the photospheric luminosity and radius using the Stefan–Boltzmann law. A temperature floor of 2000 K for T_ph(t) is applied in order to account for electron-ion recombination in the expanding ejecta. When T_ph(t) reaches the temperature floor, R_ph(t) is redefined using again the Stefan–Boltzmann law. Furthermore a Planckian blackbody spectrum is assumed at the photosphere.

The contribution to the luminosity from the thin part of the ejecta is computed by partitioning the latter into equal mass shells and by assuming that each shell with temperature T emits its radioactive decay energy assuming local thermodynamics equilibrium. To characterize the temperature of the thin part of the ejecta, we adopt a temperature profile similar to the one derived in Wollaeger et al. (2018) under the assumption of radiation dominated, homologous expansion: T(t, x) = T₀(x)(t_tr(x)/t), where T₀(x) is the temperature of the photosphere as it transits through the shell centred in x at the time t_tr(x). The bolometric luminosity contribution from the thin region is computed by multiplying the mass of each shell by the specific heating rate.

For the nuclear heating rates powering the kilonova emission, we employ the analytic fitting formula first presented in Wu et al. (2022) and based on the results from the nucleosynthesis calculations reported in Perego et al. (2022): |$\dot{\epsilon }_{\mathrm{r}}(t)=At^{-\alpha }$|⁠, where A and α are fit parameters. The latter are interpolated from tabulated values on the same (Y_e, s, τ) grid used for the nucleosynthesis calculations (see Section 2.3).

A constant thermalization efficiency ϵ_th = 0.5 is employed for the thick region of the ejecta, while we construct a thermalization efficiency profile for the thin part starting from the analytic fitting formula proposed in Barnes et al. (2016). The expression for the thermalization efficiency profile reads,

$$\begin{eqnarray} \epsilon _{\mathrm{th}}(t,x)= 0.36\left[\exp (-aX)+\frac{\ln (1+2bX^d)}{2bX^d}\right] \, , \end{eqnarray}$$

(11)

where a, b and d are the fitting parameters reported in Barnes et al. (2016) and interpolated from tabulated values on a grid spanning the intervals 1 × 10⁻³M_⊙ < M_ej < 5 × 10⁻²M_⊙ and 0.1c < v_ej < 0.3c. In the original formulation of Barnes et al. (2016), obtained assuming ρ(t) = ρ₀(t/t₀)³, X(t, x) = t. Due to the use in our model of the density profile equation (10), we adopt X(t, x) = t/(1 − x²), instead. In this work, we consider two ejecta components: a dynamical ejecta and a disc ejecta component, both symmetric with respect to the equatorial plane and to the polar axis. Following the same procedure described in Section 2.3, we directly extract from the simulations the profiles of the properties of the dynamical component, namely the distributions of the ejecta mass, of the root mean square velocity at infinity, of the average electron fraction, average entropy, and average density at the extraction radius, as a function of the polar angle θ, averaged over the azimuthal angle ϕ. The opacity κ is computed by interpolating the results of the atomic calculations performed in Tanaka et al. (2020) for a wide range of the electron fraction 0.01 ≤ Y_e ≤ 0.50. Additionally, inspired by disc simulations of Wu et al. (2016), Lippuner et al. (2017), Siegel & Metzger (2017), Fernández et al. (2019), Fahlman & Fernández (2022), we assume that a fraction between ∼20 and ∼40 per cent of the disc mass inferred from our simulations (see Section 3.3) is ejected in the form of a viscosity-driven wind. We model the mass of this disc wind as uniformly distributed in θ, as we do not expect preferential latitudes for the ejection. Moreover, for the disc ejecta we assume a root mean square velocity of 0.06c, a uniform opacity of 5 cm² g⁻¹, an average entropy of 20 k_B baryon⁻¹ and an expansion time-scale of 30 ms. We stress that our kilonova model relies on a large number of assumptions and simplifications which limit its accuracy. However, for the parameters that are not directly fixed by our simulations, we chose representative values in broad agreement with what obtained by fitting AT2017gfo data with the original kilonova model (Perego et al. 2017).

3 RESULTS

3.1 Merger dynamics

All simulations in our sample follow a qualitative common evolution pattern with quantitative differences, mainly due to the different tidal deformability provided by the EOSs and BNS mass ratios. All simulations result in the prompt collapse of the central part of the remnant into a BH. In this context, we say that a BNS simulation has resulted in a prompt collapse if the minimum of the lapse function inside the computational domain decreases monotonically immediately after merger without showing core bounces.

We define the moment of formation of the BH as the time at which the lapse function drops below 0.2. In all simulations presented here the BH forms within a fraction of a ms after the merger (t_BH < 0.47 ms, see Table 2).

Table 2.

For each simulation the table reports the rescaled angular momentum of the remnant, j_rem; the rescaled total energy radiated in GWs, |$e_{\rm GW}^{\rm tot}$|⁠; the BH expected mass (spin), |$M_{\rm BH}^{\rm exp}$| (⁠|$a_{\rm BH}^{\rm exp}$|⁠) as defined in equation (7) (equation 8); the BH mass (spin) as detected from the AH finder, M_BH (a_BH), together with the related average on a sample time, 〈M_HB〉 (〈a_BH〉). We report values from the SR simulations and the error inside brackets estimated as the absolute semidifference between the SR and LR values. Uncertainties refers to the least significant digit(s).

EOS	q	AH finder	tBH − tmrg	jrem	|$e_{\rm GW}^{\rm tot}$|	Lpeak	|$M_{\rm BH}^{\rm exp}$|	MBH	〈MBH〉	|$a_{\rm BH}^{\rm exp}$|	aBH	〈aBH〉
			(ms)			1055 [erg s−1]	[M⊙]	[M⊙]	[M⊙]
BLh	1.0	|$\checkmark$|	0.185	2.994	0.099	8.23	3.2259	3.2349	3.245	0.788	0.7860	0.801
			(2)	(8)	(1)	(13)	(2)	(2)	(2)	(2)	(1)	(2)
BLh	1.12	|$\checkmark$|	0.209	3.012	0.097	7.75	3.2250	3.2330	3.245	0.789	0.7865	0.802
			(2)	(8)	(1)	(22)	(5)	(<10−1)	(2)	(2)	(3)	(2)
BLh	1.18	|$\checkmark$|	0.209	3.020	0.098	7.19	3.2411	3.2458	3.259	0.789	0.7866	0.803
			(30)	(6)	(1)	(9)	(18)	(4)	(2)	(2)	(1)	(3)
BLh	1.33	|$\checkmark$|	0.221	3.067	0.090	5.53	3.2559	3.2573	3.273	0.780	0.7779	0.796
			(8)	(6)	(1)	(8)	(2)	(6)	(1)	(5)	(<10−1)	(3)
DD2	1.0	✗	0.422	3.122	0.092	5.46	3.2210	– –	0.826	–	–
			(10)	(9)	(2)	(18)
DD2	1.18	✗	0.445	3.117	0.091	4.96	3.2298	–	–	0.820	–	–
			(6)	(6)	(1)	(12)
DD2	1.33	✗	0.469	3.149	0.0877	4.06	3.2315	–	–	0.780	–	–
			(41)	(2)	(2)	(3)
DD2	1.67	✗	0.374	3.204	0.077	2.89	–	–	–	–	–	–
			(2)	(3)	(3)	(4)
SFHo	1.0	|$\checkmark$|	0.138	2.953	0.102	9.98	3.223	3.25	3.26	0.778	0.774	0.79
			(2)	(14)	(1)	(22)	(1)			(1)
SFHo	1.18	|$\checkmark$|	0.138	2.976	0.097	8.86	3.240	3.27	3.28	0.776	0.775	0.79
			(18)	(8)	(1)	(17)	(1)			(2)
SFHo	1.33	|$\checkmark$|	0.126	3.066	0.0872	7.32	3.268	3.29	3.29	0.783	0.770	0.79
			(8)	(17)	(4)	(16)
SLy4	1.0	✗	0.138	3.031	0.105	10.90	3.2167	–	–	0.801	–	–
			(18)	(6)	(1)	(32)	(1)			(2)
SLy4	1.18	✗	0.114	3.010	0.103	9.67	3.2323	–	–	0.791	–	–
			(14)	(12)	(1)	(23)	(6)		(3)
SLy4	1.33	✗	0.114	3.043	0.097	7.97	–	–	–	–	–	–
			(2)	(9)	(1)	(7)

Tidal forces deform the NSs during the inspiral, especially the lighter and less compact one. This effect is more pronounced for BNS with stiffer EOSs, providing, for the same gravitational mass, a less compact NS. The subsequent merger dynamics is able to unbind matter from the tidal tails on a few dynamical time-scales. The neutron-rich matter ballistically expelled during this phase from the tidal tails has low entropy and can have large enough velocity to escape the potential barrier, contributing to the dynamical ejecta. The otherwise gravitationally bound matter forms a disc with toroidal shape around the forming BH. BNS models characterized by a stiffer EOS expel more matter, such that more dynamical ejecta and larger discs are found, as discussed in detail below.

During the few fractions of ms that precede BH formation, a small amount of very high-entropy matter coming from the NS contact interface is expelled, see Fig. 2. This extremely shocked matter is characterized by higher entropy and electron fraction than the ones that characterize matter expelled by tidal forces. This small component with entropy of |$90\!-\!120~\rm {k_B}~\rm {baryon^{-1}}$| is responsible of the bimodal distribution of the entropy shown in Fig. 3. Its unbound component contributes to the dynamical ejecta, while the bound mass contributes to the disc formation, spanning in both cases a broader polar angle than the bound and unbound matter of tidal origin.

Figure 2.

Snapshot of the rest mass density (left) and the entropy per baryon (right) taken at ∼0.3 ms after BH formation across the orbital plane for the equal mass BNS merger SR simulation with the SFHo EOS. Matter inside the dashed contour with entropy 90–120 k_B baryon⁻¹ and densities <10⁸ g cm⁻³ comes from the rotationally non-symmetric central object, expelled from the contact surface of the two stars. Since equal mass binaries eject few |$10^{-5} \mathrm{\, M_\odot }$|⁠, this shocked matter have a prominent role in the median properties of the ejecta.

Figure 3.

Histograms of the dynamical ejecta. From the first to the last column: velocity at infinity v_∞, electron fraction Y_e, entropy per baryon s and polar angle θ_ej. Each row represents a different EOS. From the first to the last line: BLh, DD2, SFHo, SLy4. As a representative case, we represent the median and the average values of all quantities for the q = 1.33 cases as vertical solid and dashed lines, respectively. The high Y_e tail in the BLh, q = 1.33 case is not robust due to the finite size of the EOS tables not extending above Y_e = 0.6.

The resulting disc, ejecta, and the central BH will be the focus of Sections 3.3 and 3.4.

3.2 Gravitational-wave luminosity

In the left columns of Table 2, we report GW data (i.e. j_rem, |$e_{\rm GW}^{\rm tot}$|⁠, and L_peak) as extracted from our GW190425-like BNS simulations.

We first test the quasi-universal relation between |$e_{\rm GW}^{\rm tot}$| and j_rem given in Zappa et al. (2018): |$e_{\rm fit}^{\rm tot}(j_{\rm rem}) = c_2 j_{\rm rem}^2 + c_1 j_{\rm rem} + c_0$|⁠, with c₀ = 0.95, c₁ = −0.44 and c₂ = 0.053.³ These coefficients were fitted over a data set containing more than 200 BNS merger simulations performed with the bam (Brügmann et al. 2008) and thc codes. The BNS simulations were grouped in four categories according to the fate of the remnant: prompt collapse to a BH, short-lived hypermassive NS, supramassive NS, and stable NS. This simple quadratic polynomial in j_rem was very effective in relating the angular momentum of the remnant with the total radiated energy in the whole data set, despite the different fates of the remnants, nuclear EOSs, and intrinsic BNS parameters. Moreover, the ranges j_rem ∈ [2.944, 3.204] and |$e_{\rm GW}^{\rm tot} \in [0.077, 0.105]$| are compatible with the respective ranges presented in Zappa et al. (2018) for the case of BNS resulting in a prompt collapse. We notice that the absolute value of the relative error |$\left| e_{\rm fit}^{\rm tot}- e_{\rm GW}^{\rm tot} \right| / e_{\rm GW}^{\rm tot}$||$\lesssim \mathcal {O}(0.1)$| is in accordance with the residuals plotted in fig. 4 of Zappa et al. (2018). Additionally, we remark that |$e_{\rm GW}^{\rm tot} \lt e_{\rm GW}^{\rm fit}$|⁠, also in accordance with the behaviour of the prompt-collapse simulations in Zappa et al. (2018). To further test the quality of the fit results with respect to the uncertainties of numerical origin we compute the ratio between the residuals and the estimated total error due to resolution uncertainties, |$\sqrt{ {\delta e_{\rm GW}^{\rm tot}}^2 + \delta {e_{\rm fit}^{\rm tot}}^2 }$|⁠, where |$\delta e_{\rm fit}^{\rm tot} = \sqrt{ 4 c_2^2 j_{\rm rem}^2 + c_1^2} ~ \delta j_{\rm rem}$|⁠. The uncertainties of numerical origin, δj_rem and |$\delta e_{\rm GW}^{\rm tot}$|⁠, are computed as the absolute value of the semidifference between SR and LR results.

The typical values are ≲ 1, indicating that the numerical error accounts for a significant fraction of the observed discrepancy.

Finally we emphasize that the rescaled GW peak luminosity, |$(q/\nu)^2 \, L_{\rm peak}$|⁠, and |$\kappa _2^L$| coefficient span the same range of the prompt collapse BNSs reported in fig. 2 of Zappa et al. (2018), i.e. [1.11, 2.36] × 10⁵⁸ erg s⁻¹ and [134, 387], respectively. We recall that |$\kappa _2^L$| is the coefficient that parametrizes the leading effect of tides on the GW emission from a BNS merger in the post-Newtonian expansion, equation (2).

3.3 Remnant properties

Remnants in our simulations are characterized by a light accretion disc surrounding a spinning BH formed |$\lesssim 0.5~\rm ms$| after the merger. In the following we present the properties of both as extracted from our simulations.

3.3.1 Accretion disc

During the last few orbits, the disc starts to form because of the tidal interaction between the two stars. In high-mass binaries resulting in prompt BH formation, the tidal interaction that occurs before and at merger is the major source of the disc. A few ms after merger the disc mass and angular momentum reach a quasi-steady phase, and slowly decrease until the end of the simulation.

In Fig. 4, we report the mass (filled markers) and angular momentum (unfilled markers) of the discs once they have reached their quasi-steady phase (i.e. ∼5–7 ms after merger), computed as the integral of mass and angular momentum densities⁴ extracted from our simulations. The masses (angular momenta) span a broad range from |$\sim 10^{-5} \mathrm{\, M_\odot }$| to |$0.1 \mathrm{\, M_\odot }$||$(10^{-4} \!-\! 1 \,\, \mathrm{\, M_\odot }^2)$| depending on the BNS parameters. Both the disc mass and angular momentum increase as a function of the mass ratio q. We find that the increase is more pronounced for stiffer EOSs, where the tidal interaction is more efficient due to the larger |$\tilde{\Lambda }$|⁠. For example, considering the trend for fixed q = 1.33, the DD2 simulation (⁠|$\tilde{\Lambda }= 254$|⁠) leads to the formation of a disc twice more massive than the one formed in the BLh simulation (⁠|$\tilde{\Lambda }= 135$|⁠) and roughly six times more massive than those in the SFHo (⁠|$\tilde{\Lambda }= 105$|⁠) and SLy4 (⁠|$\tilde{\Lambda }= 93$|⁠) simulations. The errors on the disc mass, estimated when both resolutions are available as the absolute semidifference between the SR and LR are in the range 25–40 per cent for very light discs and get smaller (∼3 per cent) as the disc mass increases above |$10^{-3}\mathrm{\, M_\odot }$|⁠. Resolution effects are higher for the BLh simulation with q = 1.18, for which the disc mass of the LR simulation is ∼14 times larger than the SR one. Despite efforts, we did not find the origin of such difference.

Figure 4.

Disc mass (filled markers) and angular momentum (empty markers) at 4–7 ms after merger for SR simulations. Mass and angular momentum increase with the mass ratio. The trends suggest a link between mass and angular momentum since |$c J_{\rm disc}/ G \sim (8 - 10) \ \mathrm{ M}_{\odot } \, M_{\rm disc}$|⁠. Errors are estimated as |SR − LR| when the LR is available.

Fig. 4 suggests a correlation between the mass and the angular momentum of the disc, i.e. |$J_{\rm disc} \sim (8 - 10) \mathrm{ M}_{\odot } \, G M_{\rm disc} /c$|⁠, possibly independent from the EOS and mass ratio. Stated differently, the mean specific angular momentum of the disc is (roughly) constant: |$J_{\rm disc}/ M_{\rm disc} \sim (8 - 10) \mathrm{ M}_{\odot } \, G /c$|⁠.

To provide a possible explanation, we consider the radial density distributions, σ(r) = ∫dϕdz ρ(r, ϕ, z), as obtained from our numerical simulations, and we approximate it with a Gaussian peak smoothly connected to a radial power law,

$$\begin{eqnarray} \sigma (r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b \exp \left(- \dfrac{(r-r_{\rm peak})^2}{2 s^2} \right) & 0 \le r \le r^{*} \\ \sigma _0 \left(\dfrac{r}{r^{*}} \right)^{-\alpha } & r \gt r^{*} \end{array}\right. \end{eqnarray}$$

(12)

where b, r_peak, s, and α are fitted against the actual radial density distribution in our simulations, while σ₀ and r* are fixed requiring σ(r) to be differentiable in r*. The parameter values and the quality of the fit are described in Appendix  A. Additionally, we assume a Keplerian angular velocity profile, |$\omega _{\rm kep}(r) = \sqrt{G M_{\rm BH}/r^3}$|⁠, inside the disc. The mass and angular momentum of the resulting Keplerian disc are

$$\begin{eqnarray} M_{\rm disc}^{\rm kep} = \int _{0}^{\infty } r \sigma (r) {\rm d}r, J_{\rm disc}^{\rm kep} = \int _{0}^{\infty } r^3 \sigma (r) \, \omega _{\rm kep}(r) {\rm d}r. \end{eqnarray}$$

(13)

In Fig. 5, we show the result of the fit for σ(r) (blue dashed line) on the numerical one (blue dots) for the simulation with the BLh EOS and q = 1.33. We also show the radial angular momentum density from the numerical simulation (purple points) and the corresponding Keplerian analogue computed from equation (13) with the fitted σ(r) (purple dashed line). We found that |$J_{\rm disc}^{\rm kep} \lesssim J_{\rm disc}$|⁠, usually within 30 per cent over more than two orders of magnitudes in J_disc. We excluded the discs of equal mass BNS from this analysis since they are very light and 40–100 per cent of their mass is inside the innermost stable circular orbit (ISCO) predicted according to the BH properties. Such discs will be accreted by the BH on the viscous timescale.

Figure 5.

Disc’s radial density (blue points, left y-axis) and radial angular momentum density (purple points, right y-axis) for the BNS with BLh EOS and q = 1.33. The blue dashed line is σ(r) fitted on the numerical data, while the purple dashed line is the corresponding Keplerian angular momentum density. The vertical dashed line is the boundary between the Gaussian and the power-law r* in equation (12). The vertical solid line is R_ISCO.

Given equations (12)–(13), the ratio between |$J_{\rm disc}^{\rm kep}$| and |$M_{\rm disc}^{\rm kep}$| can be written as (see Appendix  A for a derivation)

$$\begin{eqnarray} \frac{J^{\rm kep}_{\rm disc}}{M^{\rm kep}_{\rm disc}} = \left(\eta \, \frac{\alpha - 2}{\alpha - 5/2} \sqrt{\frac{M_{\rm BH}}{\mathrm{\, M_\odot }} \frac{2 r^{*}}{R^{\rm Sch}_{\odot }}} \right)~\frac{G \mathrm{\, M_\odot }}{c} \, , \end{eqnarray}$$

(14)

where η is defined as in equation (A9) and varies between 0.78 and 0.90 with average 0.83 in our numerical simulations, |$R_{\odot }^{\rm Sch}$| is the Schwarzschild radius of the Sun, r* is such that 21 km ≲ r* ≲ 40 km, while |$M_{\rm BH} \approx 3.21-3.26 \mathrm{\, M_\odot }$| (see Section 3.3.2). The parameter which is subject to more significant variation is α ∈ [4.0, 13.9] whose average is 7.5 (see Appendix  A for the values of α and r*). Inserting these ranges of values in equation (14), one obtains |$J^{\rm kep}_{\rm disc}/M^{\rm kep}_{\rm disc} \sim 6-9 \mathrm{\, M_\odot }$| with average of |$7.3 \mathrm{\, M_\odot }$|⁠, in agreement within ≈83 per cent with the average |$\langle J_{\rm disc} / M_{\rm disc} \rangle = 8.8 \mathrm{\, M_\odot }$| obtained by our numerical simulations.

3.3.2 Black hole

In Fig. 6 we report the BH irreducible and gravitational masses, and the dimensionless spin parameter as a function of time after the BH formation for the BLh simulation at SR with q = 1.33. We see that all the three quantities increase abruptly as the AH finder detects the apparent horizon. The horizontal dashed lines indicate the expected values |$M^{\rm exp}_{\rm BH}$| and |$a^{\rm exp}_{\rm BH}$|⁠, while the vertical dashed line indicates the time at which the irreducible mass reaches its maximum value (a few ms after the BH formation). Although M_irr is expected to remain constant or to increase, we find that after having reached the maximum it starts to slowly decrease. We attribute this behaviour to numerical and discretization errors in tracing the AH location. While the AH shrinks, M_BH and a_BH continue to increase without reaching saturation. Matter accretion from the disc is not sufficient to explain this growth. The rise of M_BH after the maximum of M_irr is due to the continuous increase of the BH spin, which is an artefact of our simulations. Due to these uncertainties, we decide to focus on the gravitational mass and spin parameter of the BH at the moment when the irreducible mass is maximum.

Figure 6.

Evolution of the normalized BH irreducible mass M_irr/M, gravitational mass M_BH/M and dimensionless spin parameter a_BH for a SR simulation based on the BLh EOS with q = 1.33. Horizontal dashed lines represent the expected values for the gravitational mass |$(M_{\rm ADM} - E_{\rm GW}^{\rm rad} - M_{\rm disc})/M$| and the spin parameter |$(J_{\rm ADM} - J_{\rm GW}^{\rm rad} - J_{\rm disc}) / (M_{\rm BH}^{\rm exp})^2$|⁠. Vertical dashed lines indicate the time at which the irreducible mass starts to decrease and the corresponding value on the plotted line.

In Table 2 we report the gravitational mass M_BH and the spin parameter a_BH of the BH computed on the basis of the latter definition. To give more conservative values of the BH properties, we report also the time averages of the BH mass, 〈M_BH〉, and spin parameter, 〈a_BH〉, over the first |$7~\rm ms$| after the time at which M_irr is maximum. We report the available data obtained by SR simulations and we estimate the uncertainties (when available) as the semidifference with respect to the data from the corresponding LR simulations when available. In the case of simulations employing the BLh or SFHo EOS, the AH is resolved by the AH finder and the BH properties can be analysed with appropriate accuracy. More quantitatively, M_BH and a_BH differ from the respective expected values less than 1 per cent. On the other hand, the AH finder was unable to detect the AH for the simulations employing the DD2 or SLy4 EOS. In these cases we decided not to report the corresponding values in Table 2.

Regarding the dependence of the BH properties on the initial binary parameters, the final outcome depends mostly on two effects. On one hand, energy and angular momentum are extracted from the central object via the ejection of matter and the formation of a remnant disc. On the other hand, GWs carry energy and angular momentum away. Both these effects reduce at the same time M_BH and J_BH. Since |$J_{\rm disc} \approx 10~\mathrm{\, M_\odot }~G/c~M_{\rm disc}$|⁠, the formation of a massive disc is particularly efficient in reducing the BH angular momentum, and ultimately also the spin parameter since the variation of |$a_{\rm BH}^{\rm exp}$| due to the disc formation only becomes |$\left. \delta a_{\rm BH}^{\rm exp} \right|_{\rm disc} \approx (2a_{\rm BH}^{\rm exp}-10 \mathrm{\, M_\odot }/M_{\rm BH}^{\rm exp}) \delta M_{\rm disc}/M_{\rm BH}^{\rm exp} \sim - 0.468 ~ \delta M_{\rm disc}/(\mathrm{\, M_\odot })$|⁠. As visible in Fig. 7, (quasi) equal mass binary simulations employing the DD2 EOS have the largest spin parameters, since their symmetric character produces a smaller disc mass, while their larger |$\kappa _2^{\rm L}$| implies a lower GW emission. However, very asymmetric binaries employing the same EOS produce massive discs reducing efficiently both M_BH and a_BH. A similar, but less significant effect, is also observed for simulations employing the BLh and SFHo EOSs. For simulations employing the SLy EOS (whose discs are usually the lightest), a_BH decreases with q, while M_BH/M stays roughly constant. Focusing on the (quasi-)equal mass simulations using the BLh, SFHo, or SLy4 EOS, the removal of mass and angular momentum through the disc formation becomes subdominant, while the dominant process is the GW emission. More symmetric binaries modelled with the SLy4 EOS (corresponding to lower values of |$\kappa _2^{\rm L}$|⁠), have indeed the smallest BH masses.

Figure 7.

M_BH/M and dimensionless spin parameter a_BH distribution for the SR simulations of this work. Filled markers represent the values computed by the AH finder, while empty markers represent the expected ones. Errors are computed as the absolute semidifference between SR and LR when available. For the filled markers errors are smaller than the symbol size.

3.4 Dynamical ejecta

In Table 3, we present the properties of the dynamical ejecta as extracted from our simulations, namely the mass of the ejecta, M_ej; the standard deviation (SD) of the polar (θ ∈ [0°, 180°]) and azimuthal (ϕ ∈ [0°, 360°], see Appendix  C for more details on its calculation) angular distributions, |$\theta _{\rm ej}^{\rm SD}$| and |$\phi _{\rm ej}^{\rm SD}$|⁠, respectively; the median of the distribution of the velocity at infinity, |$v_\infty ^{\rm med}$|⁠, of the electron fraction, |$Y_{e}^{\rm med}$|⁠, and of the entropy per baryon, |$s_{\rm ej}^{\rm med}$|⁠. The last column refers to the fraction of shocked ejecta X_s, defined as the fraction of the ejecta whose entropy is larger than |$10 \, k_{\rm B}~{\rm baryon^{-1}}$|⁠. We report the values for both SR and LR simulations accompanied by the 15–75 percentile range around the median computed from the respective mass-weighted histogram. We do not report the ejecta properties when |$M_{\rm ej}\lt 10^{-5} \mathrm{\, M_\odot }$|⁠, since the properties of such a small amount of ejected matter cannot be trusted due to numerical uncertainties.

Table 3.

Dynamical ejecta properties for each simulation. M_ej is the total mass of the ejecta; |$\theta _{\rm ej}^{\rm SD}$| and |$\phi _{\rm ej}^{\rm SD}$| are the mass-weighted standard deviation of the polar and azimuthal angle, respectively; |$v_\infty ^{\rm med}$|⁠, |$Y_{e}^{\rm med}$|⁠, and |$s_{\rm ej}^{\rm med}$| are the median values of the electron fraction, speed and entropy distributions. The last column is the ratio |$X_s\equiv M_{\rm ej}^{\rm shocked} / M_{\rm ej}$|⁠, where the shocked and tidal ejecta are defined as the components with entropy respectively above and below the threshold of 10 k_B baryon⁻¹. The subscript and superscript numbers indicate the 15 and 75 percentile around the median of the respective quantity.

EOS	q	Resolution	Mej	|$\theta _{\rm ej}^{\rm SD}$|	|$\phi _{\rm ej}^{\rm SD}$|	|$v_\infty ^{\rm med}$|	|$Y_{e}^{\rm med}$|	|$s_{\rm ej}^{\rm med}$|	Xs
			[10−4 M⊙]			[c]		[kB baryon−1]
BLh	1.0	SR	0.002	–	–	–	–	–	–
		LR	0.023	–	–	–	–	–	–
BLh	1.12	SR	0.039	–
		LR	0.090	–	–	–	–	–	–
BLh	1.18	SR	0.164	21.3	82.0	|$0.24 ^{{+ 0.08}} _{{- 0.12}}$|	|$0.21 ^{{+ 0.07}} _{{- 0.08}}$|	|$18.1 ^{{+ 39.4}} _{{- 11.6}}$|	0.78
		LR	0.182	23.3	89.8	|$0.21 ^{{+ 0.07}} _{{- 0.10}}$|	|$0.25 ^{{+ 0.04}} _{{- 0.07}}$|	|$41.2 ^{{+ 55.4}} _{{- 31.5}}$|	0.94
BLh	1.33	SR	0.508	18.2	74.0	|$0.27 ^{{+ 0.10}} _{{- 0.14}}$|	|$0.17 ^{{+ 0.9}} _{{- 0.5}}$|	|$9.71 ^{{+ 17.4}} _{{- 4.21}}$|	0.61
		LR	0.959	20.7	78.6	|$0.29 ^{{+ 0.10}} _{{- 0.15}}$|	|$0.16 ^{{+ 0.14}} _{{- 0.5}}$|	|$12.3 ^{{+ 22.0}} _{{- 6.87}}$|	0.63
DD2	1.0	SR	0.586	26.3	95.1	|$0.28 ^{{+ 0.09}} _{{- 0.12}}$|	|$0.27 ^{{+ 0.04}} _{{- 0.06}}$|	|$33.2 ^{{+ 38.8}} _{{- 18.3}}$|	1.00
		LR	0.416	23.8	92.1	|$0.32 ^{{+ 0.06}} _{{- 0.08}}$|	|$0.29 ^{{+ 0.03}} _{{- 0.05}}$|	|$47.1 ^{{+ 42.4}} _{{- 31.4}}$|	1.00
DD2	1.18	SR	7.16	21.4	122	|$0.27 ^{{+ 0.10}} _{{- 0.14}}$|	|$0.17 ^{{+ 0.05}} _{{- 0.06}}$|	|$10.28 ^{{+ 7.18}} _{{- 4.12}}$|	0.57
		LR	9.67	18.1	87.3	|$0.27 ^{{+ 0.11}} _{{- 0.15}}$|	|$0.19 ^{{+ 0.06}} _{{- 0.08}}$|	|$9.36 ^{{+ 5.42}} _{{- 3.80}}$|	0.63
DD2	1.33	SR	4.00	17.3	76.6	|$0.23 ^{{+ 0.08}} _{{- 0.11}}$|	|$0.15 ^{{+ 0.05}} _{{- 0.05}}$|	|$9.38 ^{{+ 3.64}} _{{- 3.66}}$|	0.65
		LR	3.94	21.7	80.7	|$0.19 ^{{+ 0.10}} _{{- 0.11}}$|	|$0.13 ^{{+ 0.8}} _{{- 0.05}}$|	|$9.34 ^{{+ 5.15}} _{{- 3.29}}$|	0.52
DD2	1.67	SR	4.05	11.1	103	|$0.20 ^{{+ 0.14}} _{{- 0.14}}$|	|$0.10 ^{{+ 0.03}} _{{- 0.07}}$|	|$5.66 ^{{+ 4.27}} _{{- 1.87}}$|	0.29
		LR	6.20	13.0	95.8	|$0.13 ^{{+ 0.13}} _{{- 0.8}}$|	|$0.06 ^{{+ 0.08}} _{{- 0.03}}$|	|$6.15 ^{{+ 3.70}} _{{- 3.33}}$|	0.37
SFHo	1.0	SR	0.023	–	–	–	–	–	–
		LR	0.033	–	–	–	–	–	–
SFHo	1.18	SR	0.071	–	–	–	–	–	–
		LR	0.151	24.5	90.6	|$0.22 ^{{+ 0.07}} _{{- 0.10}}$|	|$0.26 ^{{+ 0.03}} _{{- 0.04}}$|	|$72.3 ^{{+ 51.3}} _{{- 53.1}}$|	0.97
SFHo	1.33	SR	0.603	12.7	68.8	|$0.26 ^{{+ 0.10}} _{{- 0.13}}$|	|$0.13 ^{{+ 0.04}} _{{- 0.06}}$|	|$7.55 ^{{+ 4.97}} _{{- 3.30}}$|	0.37
		LR	1.87	13.1	85.0	|$0.32 ^{{+ 0.10}} _{{- 0.16}}$|	|$0.13 ^{{+ 0.05}} _{{- 0.05}}$|	|$6.45 ^{{+ 5.08}} _{{- 2.50}}$|	0.32
SLy4	1.0	SR	0.030	–	–	–	–	–	–
		LR	0.024	–	–	–	–	–	–
SLy4	1.18	SR	0.055	–	–	–	–	–	–
		LR	0.114	21.4	79.5	|$0.22 ^{{+ 0.10}} _{{- 0.10}}$|	|$0.24 ^{{+ 0.05}} _{{- 0.06}}$|	|$38.1 ^{{+ 97.5}} _{{- 31.4}}$|	0.79
SLy4	1.33	SR	2.29	9.0	71.5	|$0.40 ^{{+ 0.12}} _{{- 0.20}}$|	|$0.10 ^{{+ 0.03}} _{{- 0.02}}$|	|$5.48 ^{{+ 1.82}} _{{- 3.15}}$|	0.22
		LR	1.12	14.6	70.8	|$0.30 ^{{+ 0.10}} _{{- 0.14}}$|	|$0.12 ^{{+ 0.09}} _{{- 0.5}}$|	|$7.40 ^{{+ 8.42}} _{{- 4.44}}$|	0.49

Additionally, in Fig. 3, we present mass histograms of the v_∞, Y_e, s_ej, and θ_ej distributions for simulations at SR for which |$M_{\rm ej}\ge 10^{-5} \mathrm{\, M_\odot }$|⁠. The vertical solid (dashed) lines represent the medians (average) of the ejecta properties for the q = 1.33 cases, taken as representative case. While the difference between mean and median is small or even negligible for the velocity and the electron fraction, a significant difference is clear in the entropy distribution.

The ejecta mass ranges from values smaller than |$10^{-5}\mathrm{\, M_\odot }$| up to |$\sim 6 \times 10^{-4} \mathrm{\, M_\odot }$|⁠, increasing with the mass ratio q and the stiffness of the EOS, as visible in Fig. 8. For asymmetric systems (q ≠ 1) and stiffer EOSs, the tidal interaction is more efficient in deforming the secondary NS and the resulting merger dynamics is more effective in expelling matter from its tidal tails (see e.g. Hotokezaka et al. 2013; Bauswein, Goriely & Janka 2013; Sekiguchi et al. 2015; Rosswog 2015; Lehner et al. 2016; Dietrich et al. 2017; Bernuzzi et al. 2020). Simulations employing the DD2 EOS exhibit a deviation from this trend at higher mass ratios (⁠|$q = 1.33, \, 1.67$|⁠), for which the value of the ejecta mass saturates or even tends to decrease, similarly to what found in Dudi et al. (2021; see section 5). We speculate that the ejection process at high q’s is more sensitive to usually subdominant effects, including the detailed behaviour of the NS radius and of |$\tilde{\Lambda }$|⁠, see Fig. 1 and Table 1. For the latter quantity, for high-q BNSs, models employing the DD2 show a decreasing |$\tilde{\Lambda }$| (see Table 1). It suggest that for asymmetric enough BNS (q ≳ 1.2 in our case), if an additional increase of the asymmetry is not accompanied by and increase of |$\tilde{\Lambda }$|⁠, the ejecta mass can saturate or even decrease. More simulations at higher resolutions are needed to confirm the robustness of this trend.

Figure 8.

Dynamical ejecta mass as a function of the mass ratio q of the binary. Different symbols denote numerical simulations with different EOS. Simulations with |$M_{\rm ej}\lt 10^{-6} \,\, \mathrm{\, M_\odot }$| have been excluded, while only ejecta with |$M_{\rm ej}\gt 10^{-5} \,\, \mathrm{\, M_\odot }$| is trusted due to numerical uncertainties. Errors are computed as the absolute difference between SR and LR values.

The SD of the geometrical angles gives an indication of the spatial distribution of the ejected matter. We find that the ejecta spread over the whole space, but it is mostly concentrated close to the equator, with an opening angle |$2 \theta _{\rm ej}^{\rm SD}$| that varies across the range 18°–54°, depending on the binary properties and where higher values correspond to more symmetric binaries. This can be understood since the tidal interaction tends to distribute matter along the orbital plane. The SD of the azimuthal angle |$\phi _{\rm ej}^{\rm SD}$| is related to the rotational symmetry of the dynamical ejecta around the orbital axis. For a mass distribution uniform in ϕ and centred in 180° with symmetric support on 2α ∈ [0, 360°], we expect a SD of |$\phi _{\rm ej}^{\rm SD}= (\sqrt{3}/3)~\alpha \approx 52^{\circ } (\alpha /90^{\circ })$|⁠. The values of |$\phi _{\rm ej}^{\rm SD}$| obtained in our simulations range within 65°–96° and are compatible with a uniform distribution centred in 180° with support on ∼225°–360° respectively, where higher values correspond to equal-mass systems. This indicates that the dynamical ejecta expelled by symmetric binaries is distributed over the whole azimuthal angle, while the anisotropy increases with q (see e.g. Bovard et al. 2017; Radice et al. 2018b; Bernuzzi et al. 2020).

The distribution of the radial velocity at infinity has |$v_\infty ^{\rm med}$| ranging from |$\sim 0.2 \, c$| to |$\sim 0.4 \, c$|⁠, with fast tails reaching |$\sim 0.6-0.9 \, c$| for the highest mass ratios.

The median of the electron fraction distribution is always smaller than 0.3 and is lower for higher mass ratios: tidal interaction ejects cold neutron rich material only marginally subject to composition reprocessing from positron and neutrino captures (e.g. Wanajo et al. 2014; Sekiguchi et al. 2015; Perego et al. 2017; Martin et al. 2018).

Finally, the entropy per baryon has a distribution with a marked peak at relatively low entropy, between ∼5 k_B baryon⁻¹ and ∼20 k_B baryon⁻¹, and a slow decrease towards higher entropy, with medians that in the SR cases range between ∼5 k_B baryon⁻¹ and ∼18 k_B baryon⁻¹ (with the only exception of the q = 1 simulation employing the DD2 EOS, and more often ≲ 10 k_B baryon⁻¹). All the entropy distributions show a second peak around s_ej ∼ 120 k_B baryon⁻¹ whose relative importance decreasing with q and with the stiffness of the EOS, ranging approximately between 10⁻² and 10⁻³. This high-entropy component reflects the presence of a shocked fraction of the ejecta coming from the collisional interface of the two NSs (see Section 3.1 and Fig. 2). We expect this component to be present also in BNS mergers characterized by lower total masses (and often not resulting in a prompt collapse), in which the total amount of ejected matter is typically larger than what found in our simulations.

The compositional properties of the dynamical ejecta show distributions comparable to what studied in Most et al. (2021) for the case of an irrotational binary, with similar fast-tail, high ye and high entropy components.

In the analysis outlined above, we have found that many properties of the ejected matter correlate with q and with the EOS stiffness. We now explicitly explore correlations among the different ejecta properties. In Fig. 9, we show M_ej, |$Y_{e}^{\rm med}$| and |$\theta _{\rm ej}^{\rm SD}$| as a function of |$s_{\rm ej}^{\rm med}$| for each BNS simulation producing more than |$10^{-5}\mathrm{\, M_\odot }$| of dynamical ejecta. We recall that lower |$s_{\rm ej}^{\rm med}$| correspond to higher values of q. In the left-hand panel we observe that M_ej is larger for lower values of |$s_{\rm ej}^{\rm med}$| and it is usually greater for stiffer EOSs. In the two middle panels, we observe that both |$\theta _{\rm ej}^{\rm SD}$| and |$Y_{e}^{\rm med}$| increase almost linearly with the logarithm of the median of the entropy distribution. This confirms that the tidal interaction tends to distribute cold, low-entropy ejecta along the orbital plane. Only for simulations in which the shock-heated component is relevant (i.e. symmetric or nearly symmetric BNSs), the angular distribution of the ejecta departs significantly from the orbital plane, indicating that shocked matter spreads more over the solid angle. Similar results were found also for unequal-mass binaries that do not collapse promptly into a black hole. (see e.g. Bauswein et al. 2013; Lehner et al. 2016; Dietrich et al. 2017; Radice et al. 2018b; Bernuzzi et al. 2020; Nedora et al. 2021a). In the right-hand panel, we study the correlations between the median of the entropy and the median of the velocity at infinity. In our simulations |$v_\infty ^{\rm med}$| decrease with |$s_{\rm ej}^{\rm med}$|⁠, indicating that higher mass ratios result in faster ejecta, contrary to what usually found in relation to systems characterized by smaller total masses. This could be indeed a peculiar property of very massive BNSs.

Figure 9.

Correlation of the ejecta mass M_ej, standard deviation of the polar angle |$\theta _{\rm ej}^{\rm SD}$|⁠, median of the electron fraction |$Y_{e}^{\rm med}$| and median of the velocity at infinity |$v_\infty ^{\rm med}$| with the median of the entropy |$s_{\rm ej}^{\rm med}$|⁠. Uncertainties are estimated as the absolute difference between SR and LR simulations, while SR values are used to represent the points. The simulations with higher mass ratios have higher values of the ejected mass.

4 NUCLEOSYNTHESIS AND KILONOVA

4.1 Nucleosynthesis

Using the procedure outlined in Section 2.3, we compute nucleosynthesis yields for the dynamical ejecta of all our GW190425 targeted simulations.

In Fig. 10, we present nucleosynthesis yields for a subset of representative simulations at t = 30 years after merger, superimposed to the Solar residual r-process abundances taken from Prantzos et al. (2020) as a useful point of reference. To guide the comparison between the different models, the Solar residuals are scaled in order to reproduce the abundance of the simulation with q = 1.33 and the DD2 EOS at A = 130.

Figure 10.

Nucleosynthesis pattern at t = 30 years after the merger as a function of the mass number A. Left: Comparison between relative abundances from simulations employing the DD2 EOS. Right: Comparison between relative abundances from numerical relativity (NR) simulations with mass ratio q = 1.33. Black dots represent the Solar r-process abundances, taken from Prantzos et al. (2020). To guide the comparison, the Solar residuals are scaled in order to reproduce at A = 130 the abundance of the simulation with q = 1.33 and the DD2 EOS.

Unequal-mass merger simulations employing the DD2 EOS (left-hand panel) robustly produce elements between the second and the third r-process peak, without showing any substantial difference between the various mass ratios. Relative abundances are comparable to the Solar residuals with a significant excess in the third peak height with respect to the height of the second peak, and a significant production of translead nuclei. On the other hand, A ≲ 120 nuclei are systematically underproduced. A weak dependence on the value of the mass ratio is visible, with more asymmetric mergers producing on average a larger amount of heavy nuclei. These behaviours are expected given the prompt collapse of the central remnant into a BH, the tidal character of the ejection mechanism and the consequent absence of a significant high-Y_e tail in the dynamical ejecta above a critical value Y_e ≳ 0.22 (e.g. Lippuner & Roberts 2015; Radice et al. 2016), that is associated with the production of less than 10 per cent of the mass fraction of heavy nuclei above the second peak through an incomplete r-process.

The situation changes significantly when considering the DD2 equal-mass case (blue line). In fact, the relative abundances of heavy r-process nuclei (A ≳ 130 and even more for A ≳ 140) are less significant with respect to the unequal mass cases, while around the first peak the q = 1 pattern is the largest and the closest one to the Solar abundances. This is consistent with the fact that, despite having a small total mass, the bulk of the ejecta Y_e distribution for the equal-mass case lies within the interval 0.20–0.40 (see Fig. 3).

The right-hand panel of Fig. 10 shows, instead, the comparison between simulations characterized by the same mass ratio, namely q = 1.33, but different EOSs. Since the mass ratio differs significantly from 1, the nucleosynthesis outcome is in all cases similar to what described for unequal-mass merger simulations in the comparison between the DD2 simulations. All the curves are quite close to each other except around the first peak, where the spread between the various distributions becomes more evident and sensitive to the nuclear EOS, with the largest (smallest) relative values for the abundances obtained for the BLh (SLy4) EOS. Usually (and especially for equal or nearly equal mergers that do not promptly collapse to a BH), the synthesis of light r-process elements within BNS ejecta should be favoured by soft EOSs, since the higher temperatures achieved in the shock-heated ejecta component leptonise matter in a more efficient way. However, we notice that for A ≲ 120 the relative production of light r-process elements does not follow exactly this trend. This is because, for such asymmetric binaries promptly collapsing to BHs, the dynamical ejection of matter is usually dominated by the cold, neutron-rich tidal component. However a small, but non-negligible fraction of the dynamical ejecta comes from the contact surface of the colliding NSs and is characterized by relatively high entropies (see the X_s column in Table 3). The corresponding larger peak temperatures produce a tail in the Y_e distribution above ≈0.22. These ejecta are likely present in all BNS mergers, but their relatively low amount make them more relevant only in the case of mergers characterized by a very small dynamical ejecta mass. Moreover, these ejecta can more likely escape in the case of stiffer EOSs, characterized by larger radii and less deep gravitational well.

We conclude that the nucleosynthesis patterns show a mild variability, depending on the mass ratios and EOSs. However, they are comparable with the ones obtained by BNS merger simulations of lighter binary systems and do not show peculiar behaviours (see e.g. Wanajo et al. 2014; Just et al. 2015; Bovard et al. 2017; Radice et al. 2018b; Nedora et al. 2021b). Nevertheless, we point out that the nucleosynthesis yields obtained exhibit different features with respect to the Solar residuals, for example in the position and shape of the second and third r-process peaks. The fine structure of the abundance pattern in this region is indeed affected by the particular choice of the nuclear input data made for the nucleosynthesis calculations, like for example the nuclear mass model, the different fission channels considered (spontaneous, neutron-induced, β-delayed etc.) or the fission fragment distribution employed (see e.g. de Jesús Mendoza-Temis et al. 2015; Eichler et al. 2015; Goriely 2015). However, since we do not expect dynamical ejecta from high-mass BNS mergers to represent the dominant contribution to the r-process enrichment in the Universe, possible discrepancies with the solar pattern are not an issue. In addition, one should also remember that, even for high mass BNS mergers, the nucleosynthesis from the disc ejecta is expected to dominate the dynamical ejecta one.

4.2 Kilonovae

Using the model described in Section 2.4, we compute synthetic kilonova light curves for each of the SR models presented in this work for which the mass of the dynamical ejecta is larger than |$10^{-5} \mathrm{\, M_\odot }$|⁠. In Fig. 11, we present the evolution of the AB magnitudes in three representative bands (B-, r-, and K-band), for two EOSs (the stiff DD2 and the soft SLy4), and two mass ratios (q = 1.18 and q = 1.33). In general, kilonova magnitudes depend both on the distance and on the viewing angle. Regarding the former, the wide range of distances compatible with GW190425 (D = 70–250 Mpc) implies a possible uncertainty of ∼3 magnitudes, with lower magnitudes corresponding to shorter distances. On the other hand, the inclination angle is almost unconstrained by the GW190425 signal. Due to the degeneracy between viewing angle and distance, viewing angles close to the polar axis (θ_view ∼ 0°) are more compatible with larger distances, while shorter distances would imply edge-on configurations (θ_view ∼ 90°). In Fig. 11, we set D = 130 Mpc while we explore all possible viewing angles, θ_view ∈ [0°, 90°]. The amount of ejecta and their composition are the most relevant parameters in shaping kilonova light curves. In general, since GW190425-like events are expected to eject a relatively small amount of mass, the resulting kilonovae are predicted to be relatively dim and fast-evolving, compared for example with GW170817-like events. More specifically, in Fig. 11 we observe that the kilonova associated to the simulation employing the DD2 EOS and with q = 1.33 is brighter and lasts longer with respect to both the simulation employing the same EOS but with q = 1.18, and the simulation with the same mass ratio but employing the SLy4 EOS, for all bands. This mostly reflects the difference in the amount of ejecta between the different models, see Sections 3.3 and 3.4, with greater mass ejection resulting in brighter peak luminosities due to the stronger availability of nuclear fuel required for the kilonova emission.

Figure 11.

AB magnitudes in the blue, red and IR bands of CTIO telescope as a function of time. We report the results for the DD2 and SLy4 EOSs and for a binary mass ratio of q = 1.18 and q = 1.33 at standard resolution. The uncertainty in the source inclination angle (varying between 0°–90°) is represented using solid lines for θ = 0° and dotted lines for θ = 90°, with intermediate values enclosed by the above lines. The source distance is set to 130 Mpc.In each panel, the darker and lighter areas refer to two different scenarios in which 20 per cent and 40 per cent of the disc mass is expelled, respectively.

Differences in the viewing angle affect the light curves at times shorter than a couple of days, while our results are insensitive to the specific viewing angle at later times. This can be explained by considering that the slower and significantly more massive disc wind component, eventually powering the kilonova at late times (t ≳ 1 day), is assumed to be isotropic in our model. Conversely, within the first days after merger, the dynamical ejecta component plays a relevant role. The angular distribution of its mass and composition are thus reflected in the band magnitude evolution. In particular, we obtain brighter light curves in the visual bands at angles closer to the pole (θ ∼ 0°), where matter with a higher initial Y_e (and thus lower opacity) can be found. Conversely, the emission in the IR band is typically brighter close to the equatorial plane (θ ∼ 90°), where the most neutron-rich (and thus more opaque) matter is concentrated, with respect to higher latitudes. Since for each of our SR models the disc wind ejecta component is determinant in generating the kilonova emission, we test our results sensitivity with respect to its mass. In particular, we notice that the increase in the fraction of ejected disc mass from a plausible 20 per cent to an optimistic 40 per cent results in an overall gain in brightness of ∼1 magnitude for all bands at late times, when the disc ejecta component becomes dominant. We also test the sensitivity of light curves on the disc ejecta mass and composition angular distributions. We consider a density distribution ρ_wind(θ) ∝ sin θ as alternative to the isotropic case and an opacity distribution shaped as a step function with k = 1 cm² g⁻¹ for θ < 45° and k = 10 cm² g⁻¹ for θ > 45°. While such modifications on the opacity can vary the final bolometric light curves up to a factor of a few, the different mass distribution results in a model dependence on the viewing angle also at late times. More specifically, since the wind density gradually increases towards the equator, the magnitudes decrease accordingly for all bands, and we obtain the brightest emission for θ_view ∼ 90°, ∼1 magnitude below the polar one. Despite the non-negligible dependences, these tests place our uncertainty in the luminosity due to the disc parameters well below the one due to the source distance and viewing angle.

For simulations with q = 1.33, providing a prominent tidal low-Y_e ejecta component, the infrared K-band lasts several days and nearly always dominates over bluer bands, due to the prevailing presence of lanthanides-rich material synthesized through a strong r-process both in the dynamical and in the disc wind ejecta. On the other hand, in the case of the simulation with q = 1.18 and the SLy4 EOS, the considerably lower ejecta mass with a broader Y_e distribution results in lower material opacities and slightly brighter blue band light curves at early times.

Due to the evolution of the photospheric temperature, the B-band magnitude is the first to peak, within the very first few hours, promptly followed by the r-band magnitude, dominating within the first half-day after merger, while the infrared band peaks much later in time, possibly on a time-scale of days. While the precise peak times and magnitudes vary depending on the specific simulation, the presence of common trends in the light curve behaviour allow us to identify characteristic time-scales for each band in which the latter typically dominates over or is comparable to the others. In Fig. 12, we present the values of the AB magnitudes in the same three bands as in Fig. 11 at three corresponding characteristic times for each available simulation, namely at 0.3, 1.1, and 3.2 d for the B, r, and K band, respectively. Since we want now to address the possible detectability of GW190425, two possible ranges for the source distance and inclination angle are considered in order to account for the large degeneracy in the estimation of these parameters for GW190425 (see also Dudi et al. 2021, for a similar choice). Regardless of the specific band, magnitudes tend to decrease with the increase of the mass ratio, leading to emissions up to ∼8 magnitudes brighter, moving from equal-mass to strongly asymmetric mergers. Likewise, the stiffest EOS corresponds to luminosities which can be as bright as ∼6 magnitudes below the same results obtained using softer EOSs. Exceptions to these trends can be directly traced back to already emerged distinctive mass ejections. For example, the simulation employing the BLh EOS and a mass ratio of q = 1.12 returns brighter red and infrared luminosities with respect to the simulation employing the same EOS but with q = 1.18: this is due to the fact that in the first instance the computed disc mass is greater, leading to a more massive disc wind (which dominates over the dynamical component). Based on our analysis, from Fig. 12 it is clear that almost none of our models can be fully ruled out by the ZTF upper limits to the kilonova of GW190425 (shown as a dashed horizontal line), meaning that current data cannot help further constraining the model parameters. This leaves open the question as to whether the detection of events like GW190425 can shed light on the source properties, and hints to the necessity of determining the sky localization with high accuracy for these events, to employ deeper observations in order to resolve such EM counterparts.

Figure 12.

AB magnitudes in the blue, red, and IR bands of CTIO telescope at fixed characteristic times as a function of the binary mass ratio q. The kilonova is obtained assuming an ejection of |$20{{\%}}$| of the disc mass. Results are colour-coded to indicate different EOSs. Only standard resolution simulations are shown. Two cases for the source distance and inclination angle are reported, with the error bars representing the uncertainty in the source distance. The dashed horizontal line represents the upper limit for GW190425 obtained with the ZTF by the GROWTH collaboration for the r and g-band (Coughlin et al. 2019).

5 DISCUSSION

In this section, we compare the results of our work with recent publications about the modelling of GW190425 and of its EM counterparts, in particular with results reported in Dudi et al. (2021), Raaijmakers et al. (2021), and Barbieri et al. (2021).

During the preparation of this work, Dudi et al. published an independent study on GW190425 in NR. They used the bam code, a NR code which was shown to produce results consistent with WhiskyTHC (see e.g. Dietrich et al. 2018). They considered four mass ratios, ranging from 1 to 1.43, and for each of them they employed three cold, beta-equilibrated EOSs: the piecewise-polytropic EOS MPA1 (Read et al. 2009), a piecewise-polytropic representation of the tabulated DD2 EOS at the lowest available temperature, and the softer APR4 EOS (Akmal, Pandharipande & Ravenhall 1998). Each model was run at three different resolutions, with our SR being intermediate between their worst and middle resolution.

Similarly to what we found in our simulations, all the BNS models presented by Dudi et al. result in a prompt collapse. Regarding the properties of the remnant, the two works predict a comparable range for M_BH/M, while we notice that the dimensionless spin parameter obtained by Dudi et al. is systematically lower than the one obtained by our simulations by several per cents, corresponding to Δa_BH ∼ 0.05, when comparing simulations characterized by similar mass ratios and EOSs. Both analyses agree in predicting more massive discs when considering more asymmetric binaries and stiffer EOSs. In particular, the disc results for the DD2 EOS share the same trend with respect to q, both on a qualitative and quantitative level.

Moving to the comparison of the dynamical ejecta, we first notice that the amount of matter obtained for the MPA1 and APR4 EOSs by Dudi et al. increases as the binary becomes more asymmetric, similarly to what observed in our BLh, SFHo and SLy4 simulations. Similarly, the amount of ejecta from the DD2 simulations first increases then decreases with q in both analyses. However, while in the former cases the amount of ejecta are comparable among them, the values obtained for the DD2 EOS differ significantly, with the ejecta reported in Dudi et al. larger by ∼ one order of magnitude. According to the reported values, uncertainties due to different resolutions seem to account only for a fraction of this discrepancy and higher resolution seems to result in smaller ejecta masses. A potentially relevant source of discrepancy could be the different microphysical input. In addition to a more accurate temperature treatment, the presence of neutrino radiation can influence the dynamical ejecta, since simulations accounting for neutrino emission show systematically smaller dynamical ejecta masses (see e.g. Nedora et al. 2022), due to the emission of neutrinos occurring during the ejection process.

The different amount of ejecta obtained employing the DD2 EOS is directly reflected in the kilonova light curves, where for a similar mass ratio the r-band magnitudes reported in Dudi et al. are systematically brighter. In particular, while for edge-on views the results are in good agreement, for a viewing angle close to the polar axis we find up to ∼5 magnitudes of difference between light curves corresponding to the same binary configurations. On the one hand, this may reflect the substantially different mass and composition distributions resulting from the NR models. On the other hand, we also stress that the models employed for the light curves computation are significantly different: as opposed to our semi-analytic model described in Section 2.4, Dudi et al. employ a more advanced wavelength-dependent radiative transfer approach (Kawaguchi et al. 2020), for which the post-merger ejecta composition is fixed for all components. Additionally, our kilonova model decomposes the solid angle in radial slices. While this approach is reasonable for ejecta expelled over the entire solid angle, it could be inadequate for ejecta expelled only close to the equator for which it tends to underestimate magnitudes up to a few since it neglects possible lateral effects (Kawaguchi et al. 2016; Kawaguchi, Shibata & Tanaka 2018; Barbieri et al. 2019; Bernuzzi et al. 2020). Keeping in mind the above differences for the GW190425 event and working under the assumption that the location of the source was covered by ZTF, Dudi et al. disfavoured a higher number of models with respect to this work, i.e. the ones employing DD2 or MPA1 EOSs with a high mass ratio and a source configuration similar to that used in the top panels of Fig. 12. On the contrary, our results imply that only the model employing the DD2 EOS with the highest mass ratio and a source distance close to D ∼ 70 Mpc (corresponding to a edge-on view) would be disfavoured (as visible in the bottom panels of Fig. 12).

Raaijmakers et al. (2021) studied the expected photometric light curves of BNS mergers with masses in the range compatible with the posteriors of GW190425. We recall that, due to the spherical symmetry of the employed kilonova model, it was not possible to investigate the light curve dependence on the viewing angle, even if selected tests with the multidimensional possis code were performed (Bulla 2019). By fixing the source distance to 130 Mpc, we find that the spread in the magnitudes generated by the different NR models considered in this work is comparable to the comprehensive results displayed in Raaijmakers et al. (2021), which span ∼4 magnitudes at times shorter than ∼1 d. In the same time period, our light curves are generally dimmer with respect to those computed in Raaijmakers et al. (2021), with an average difference of ∼3 magnitudes. A plausible source of this systematic discrepancy lies in the different ways in which the ejecta and disc masses were computed. In our case, they are the outcome of BNS merger simulations, while in Raaijmakers et al. (2021) they are estimated on the basis of the fitting formulae for the mass of the dynamical ejecta and of the disc proposed in Krüger & Foucart (2020, equations 4 and 6), and for the average dynamical ejecta speed proposed in Foucart et al. (2017). These formulae take as input parameters the compactness and the masses of the binary components. We compare the outcome of these fitting formulae with our numerical results in Appendix  B. We found significant differences in the ejected mass and in the expansion speed, and less severe disagreement for the disc mass, which is consistent with the numerical data when errors are taken in consideration. In particular, the mass of the ejecta predicted by the fitting formulae is ∼10–100 higher than in our simulations. Our comparison reveals how NR fitting formulae can become inaccurate when used far from their calibration regime.

Finally, we compare the light curves computed in this work with those obtained in Barbieri et al. (2021) for BNS systems, and, as in the case of Raaijmakers et al. (2021), we find typically lower peak luminosities. Since also Barbieri et al. (2021) used fitting formulae to predict the ejecta properties (see Appendix  B for a more detailed discussion), we argue that disc and ejecta masses larger by one or even two orders of magnitudes can account for the observed differences. In addition, our results employing the DD2 EOS are significantly more sensitive to the binary configuration, as peak luminosities in the r-band and at IR frequencies vary by ≲ 7 magnitudes for a mass ratio varying between 1 ≤ q ≲ 1.7, while in Barbieri et al. (2021) the same bands exhibit a variation of ∼3.5 magnitudes for a mass ratio between 1 ≲ q ≲ 2. Also in this case, at least a part of these differences is possibly due to disc later irradiation, which is expected to occur in very asymmetric system, which was taken into account by Barbieri et al. (2021).

Both in Raaijmakers et al. (2021) and Barbieri et al. (2021), the overall brighter kilonovae allow the identification of some binary configurations potentially detectable by the ZTF within the first few days from merger, in addition to a major portion of the BHNS configurations considered in those works. In particular, in Barbieri et al. (2021) several configurations employing the DD2 EOS and the APR4 EOS can be ruled out by the GW190425 EM follow-up. Conversely, here almost all the our BNS simulations employing the DD2 EOS and the totality of those employing softer EOSs produce kilonovae which are not detectable by ZTF in a GW190425-like event at a comparable distance.

6 CONCLUSIONS

In this work, we investigated in detail the outcome of BNS merger simulations targeted to GW190425 with detailed microphysics. We set up 28 simulations with finite temperature, composition dependent NS EOSs, and neutrino radiation. For each simulation we extracted remnant and dynamical ejecta properties, and we computed in post-processing nucleosynthesis yields and kilonova light curves.

Using four EOSs compatible with present constraints and considering a broad range of mass ratios, we aimed at giving an accurate description of GW190425-like BNS mergers and answering a number of questions, including: what can we expect from future detection of this kind of events in terms of remnant, dynamical ejecta, nucleosynthesis signature and kilonova light curves? Despite the wide sky localization of GW190425, can the lack of an EM counterpart give constraints on the EOS and/or the binary parameters?

We found that such BNS mergers, characterized by an unusual high total mass of |$3.4 \mathrm{\, M_\odot }$| and a chirp mass of |$1.44 \mathrm{\, M_\odot }$|⁠, prompt collapse to a light black hole of |$\sim 3.2 \mathrm{\, M_\odot }$| with a dimensionless spin parameter that ranges from 0.73 to 0.83, surrounded by a light disc formed by tidal interactions. Asymmetric BNS mergers with stiffer EOS have more massive remnant disc, ranging from |$10^{-5} \mathrm{\, M_\odot }$| for equal mass binaries with soft EOS, to |$0.1 \mathrm{\, M_\odot }$| for the most asymmetric BNS in our sample.

During the late inspiral and merger, previous to the collapse, the simulated binaries expel a small amount of matter in the form of dynamical ejecta. The high compactness is responsible for less deformable NSs while the prompt collapse inhibits the production of shock-heated ejecta. This explains the lower values of ejected mass compared to what previously found for BNS whose chirp mass is closer to what is observed in the Galactic BNS population and in GW170817. Since tidal interactions are the main cause of dynamical ejection, we found that asymmetric BNS mergers with a stiff EOS are able to unbind up to |$\sim 10^{-3} \mathrm{\, M_\odot }$| of ejecta, while equal mass BNS with a soft EOS only eject |$\lesssim 5 \times 10^{-6} \mathrm{\, M_\odot }$| of matter. Also the properties mostly depend on the mass ratio and on the EOS of the BNS merger. Dynamical ejecta spread all over the space but it is mainly concentrated along the orbital plane in an opening angle which goes from 54° for symmetric BNS to 18° for the more asymmetric BNS in our sample. We also discuss the distributions of electron fraction, velocity at infinity and entropy of the dynamical ejecta and their trends with the binary parameters.

In all the considered simulations, the resulting r-process nucleosynthesis pattern does not show peculiar behaviours and reflects directly the properties of the matter outflow. For ejecta dominated by cold, neutron-rich matter, we noticed a remarkably robust production of heavy elements between the second and the third r-process peaks, as opposed to the less significant one of lighter elements. The latter is however more sensitive to the binary parameters. In fact, around the first peak the nucleosynthesis pattern changes depending on the EOS considered (even if not with a clear trend) and increases with decreasing mass ratio, but always on a lower level with respect to the Solar residuals.

For the kilonova, we found that narrow-band light curves in the B- and r-bands peak within the first few hours after the merger with a rapid subsequent decline, while the emission at IR frequencies lasts several days. Assuming a distance of 70–130 Mpc or 130–250 Mpc, compatible with what was inferred for GW190425, and combined with a edge-on or face-on inclination, respectively, the peak magnitude in every band is not brighter than ∼20 magnitudes, as opposed to the case of kilonovae resulting from BNS more compatible with the Galactic BNS population or with GW170817. As such, we conclude that it could be difficult to observe such a transient at the distances inferred for GW190425 with present wide-field surveys, unless a good sky localization allows for deeper and localized searches. This can be traced back to the low mass of the dynamical ejecta and of the disc remnant. Only a BNS with a particularly stiff EOS, a high mass ratio and a source distance around ∼70 Mpc would have been detected by the ZTF facility according to our findings. This would favour a BH-NS merger in the case of a kilonova detection resulting from a compact binary merger similar to GW190425 by ZTF.

Future follow-up campaigns will be joined by Vera Rubin (LSST) observatory. In spite of the relatively small field of view (∼10 deg²) compared to ZTF, the short read-out time, the all-sky reference and a sensitivity of 24.7–27.5 AB magnitudes in the r-band will permit Vera Rubin to be a powerful resource to detect faint kilonovae (Andreoni et al. 2022). Vera Rubin is potentially able to detect kilonova signals from some of the simulated BNS mergers. For a kilonova at a distance of 130–250 Mpc, a kilonova signal would be detectable for BNS mergers with q > 1.33 and, in the case of a very stiff EOS (as DD2) for the BNS with q = 1.18. In addition, for smaller distances, i.e. 70–130 Mpc, also kilonovae resulting from slightly asymmetric BNS mergers could be observable. Finally, for a distance comparable to the one of GW170817, all the simulated kilonovae could be potentially detected. However, despite the increased sensitivity, Vera Rubin’s field of view will cover efficiently up to 200 deg², far less than the confidence region of GW190425. Thus, a better sky localization will be crucial.

We compared our results with recent works that aim to predict the remnant and ejecta properties, as well as the kilonova light curves of GW190425. We find overall similar qualitative trends, but with some quantitative differences. In the case of Dudi et al., who explored a comparable set of simulations in numerical relativity, trends in the ejecta masses and disc masses are very similar, with a better quantitative agreement for the latter than for the former. We speculate that these differences could be due to the different microphysical setups (both polytropic EOSs and the lack of neutrino radiation tend to overestimate the dynamical ejecta) as well as resolution effects. All these uncertainties could be even amplified in this case due to the small amount of ejecta, that makes their identification and tracking inside the computational domain more challenging. Raaijmakers et al. (2021) and Barbieri et al. (2021) computed kilonova light curves for GW190425-like events and they found kilonova transients systematically brighter than ours. A plausible source of discrepancy could be the use of existing fitting formulae to predict the dynamical ejecta and the disc mass. Indeed the peculiarity of GW190425 slip to the predictions given by the formulae presented in previous works (Foucart et al. 2017; Radice et al. 2018b; Barbieri et al. 2019; Nedora et al. 2022) that we took into exam. Fitted on large sample of numerical simulations of BNS mergers with parameters however different from the ones of GW190425, they usually predict an enhancement of the dynamical ejecta and of the disc mass with respect to our simulations, with observable consequences on the kilonova. This result underlines the difficulty in providing fitting formulae for the ejecta properties valid over a broad range of binary parameters and even outside of the fitting range. This could indeed strongly affect their effectiveness.

The detection of GW190425 demonstrated that, in addition to the sample of BNS mergers whose properties are close to the ones observed in the current population of Galactic BNS systems, there could be a population of GW-loud events characterized by larger chirp masses. Their modelling is less developed and their properties (including the smaller ejecta and disc masses) are possibly more challenging to study. Our work represents a step forward in the direction of better characterizing such systems. Considering the GW190425 follow-up campaign, we conclude that, even assuming that the sky coverage was enough and the binary was a BNS system, no strong constraints on the BNS parameters nor on the EOS can be inferred by the lack of EM signal. Only the corner case of very stiff EOS and extreme mass ratios could be possibly excluded. Future observations of EM counterparts by wide-field surveys, such as ZTF or Paolmar Gattini-IR telescope, for such a population outsider will be non-trivial, unless the merger distance decreases to ≲ 40 Mpc. However, large uncertainties still remain. We mostly quantified errors due to finite resolutions, but we expect possibly larger uncertainties due to systematics and modelling limitations. Further works in the modelling of both BNS mergers and their EM counterparts is required to properly assess these limitations.

ACKNOWLEDGEMENTS

We thank Andrea Endrizzi for initial work on the project. The Authors acknowledge the INFN (Instituto Nazionale di Fisica Nucleare) and Virgo for the usage of computing and storage resources through the tullio cluster in Torino. AP acknowledge PRACE (Partnership for Advanced Computing in Europe) for awarding him access to Joliot-Curie at GENCI@CEA (Grand Équipement National de Calcul Intensif). He also acknowledges the usage of computer resources under a CINECA-INFN (Consorzio Interuniversitario per il Calcolo Automatico) agreement (allocation INF20_teongrav and INF21_teongrav). S.B. acknowledges funding from the EU H2020 under ERC (European Research Council). I added Starting Grant, no. BinGraSp-714626, and from the Deutsche Forschungsgemeinschaft, DFG, project MEMI number BE 6301/2-1. D.R. acknowledges funding from the U.S. Department of Energy Office of Science, Division of Nuclear Physics under Award Number(s) DE-SC0021177 and from the National Science Foundation under Grants No. PHY-2011725, PHY-2020275, PHY-2116686, and AST-2108467. FMG acknowledges funding from the Fondazione Caritro, program Bando post-doc 2021, project number 11745. NR simulations were performed on Joliot-Curie at GENCI@CEA (PRACE-ra5202), SuperMUC-LRZ (Gauss project pn56zo), Marconi-CINECA (ISCRA-B project HP10BMHFQQ, INF20_teongrav and INF21_teongrav allocation); Bridges, Comet, Stampede2 (NSF XSEDE allocation TG-PHY160025), NSF/NCSA Blue Waters (NSF AWD-1811236), supercomputers. 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 U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

DATA AVAILABILITY

Data generated for this study will be made available upon reasonable request to the corresponding authors.

Footnotes

1

In both works, the focus was broader than GW190425 kilonova characterization, but this event was extensively studied as realistic test case.

2

The location of the components is determined by the location of the photospheres.

3

We notice that, despite referring to the same fit, the fitting values reported in this work have one more figure than the ones originally reported by Zappa et al. (2018).

4

This approach assumes that the metric is axisymmetric.

5

Multiple δ* that satisfy this condition can exist, so we also add the condition that the mode of the distribution lies in the interval π − π/4 ≤ ϕ ≤ π + π/4.

APPENDIX A: DETAILS OF THE KEPLERIAN MODEL

To deduce equation (14) we define

$$\begin{eqnarray} M_{\rm disc}^{\rm kep} \equiv M_{\rm disc}^{\rm G} + M_{\rm disc}^{\alpha }~, && J_{\rm disc}^{\rm kep} \equiv J_{\rm disc}^{\rm G} + J_{\rm disc}^{\alpha }~, \end{eqnarray}$$

(A1)

where the superscript G and α indicate the Gaussian and power-law parts of the Keplerian disc in equations (12) and (13),

$$\begin{eqnarray} M_{\rm disc}^{\rm G} \equiv \int _{R_{\rm ISCO}}^{r^{*}} r \sigma (r) {\rm d}r~, && M_{\rm disc}^{\rm \alpha } \equiv \int _{r^{*}}^{r_{\rm max}} r \sigma (r) {\rm d}r~, \end{eqnarray}$$

(A2)

and similar for the angular momentum. We can solve the integration,

$$\begin{eqnarray} M_{\rm disc}^{\rm G} &=& bs^2 \bigg[ \sqrt{\frac{\pi }{2}} \frac{r_{\rm peak}}{s} {\rm erf}\left(\frac{r - r_{\rm peak}}{\sqrt{2} s} \right)\\ &&-\, \exp \left(- \frac{(r - r_{\rm peak})^2}{2 s^2} \right) \bigg]_{R_{\rm ISCO}}^{r*}, \end{eqnarray}$$

(A3)

$$\begin{eqnarray} M_{\rm disc}^{\alpha } = \frac{\sigma _0}{\alpha -2 } \left(1 - \frac{(r^{*})^{\alpha - 2}}{r_{\rm max}^{\alpha - 2}} \right)(r^{*})^2~, \end{eqnarray}$$

(A4)

$$\begin{eqnarray} J_{\rm disc}^{\rm G} = \sqrt{\frac{1}{2} G M_{\rm BH} r_{\rm peak}^3 (bs)^2} \sum _{k=0}^{\infty } \binom{3/2}{k} \left(\frac{ \sqrt{2} s}{r_{\rm peak}} \right)^k \end{eqnarray}$$

(A5)

$$\begin{eqnarray} \hphantom{J^G=}\times \left. \Gamma \left(\frac{k+1}{2},~\frac{(r - r_{\rm peak})^2}{2 s^2} \right) \right|_{R_{\rm ISCO}}^{r^{*}}~, \end{eqnarray}$$

(A6)

$$\begin{eqnarray} J_{\rm disc}^{\alpha } = \frac{\sigma _0 \sqrt{G M_{\rm BH}}}{\alpha - 5/2} \left(1 - \frac{(r^{*})^{\alpha - 5/2}}{r_{\rm max}^{\alpha - 5/2}} \right)(r^{*})^{5/2}~, \end{eqnarray}$$

(A7)

where |${\rm erf}(x) \equiv (2/\sqrt{\pi }) \int _0^x e^{-t}~{\rm d}t$| is the error function and |$\Gamma (a,x) \equiv \int _{x}^{\infty } t^{a-1} e^{-t} \mathrm{ d}t~$| the incomplete gamma function. One can write

$$\begin{eqnarray} \frac{J^{\rm kep}_{\rm disc}}{M^{\rm kep}_{\rm disc}} = \eta \frac{J^{\alpha }_{\rm disc}}{M^{\alpha }_{\rm disc}}~, \end{eqnarray}$$

(A8)

where

$$\begin{eqnarray} \eta = \frac{1 + J^{\rm G}_{\rm disc} / J^{\alpha }_{\rm disc}}{1 + M^{\rm G}_{\rm disc} / M^{\alpha }_{\rm disc}} ~. \end{eqnarray}$$

(A9)

Assuming r* ≪ r^max (with an error ≲ 1 per cent) we arrive at

$$\begin{eqnarray} \frac{J^{\rm kep}_{\rm disc}}{M^{\rm kep}_{\rm disc}} = \eta \frac{\alpha - 2}{\alpha - 5/2} \sqrt{G M_{\rm BH} r^{*}}~. \end{eqnarray}$$

(A10)

As showed in Fig. 5, the model tends to underestimate the radial angular momentum density, especially for r < r*. To better quantify this difference, in Fig. A1 we compare the angular momentum of the discs from our simulations at SR with the corresponding Keplerian analogue, equation (13). With the exception of DD2 EOS with q = 1.67, the discrepancy is <30 per cent.

Figure A1.

Top: Comparison between the disc angular momentum outside the ISCO from numerical simulations, J_disc, and the one obtained by constructing a Keplerian disc whose radial density profile was fitted over the numerical results using equation (12), |$J_{\rm disc}^{\rm kep}$|⁠. Bottom: Relative difference between the two values. Unfilled markers represent discs for which the Keplerian mass differs from the numerical one by more than 20 per cent.

In Fig. A2, we show the power-law exponent α, obtained by fitting equation (13) over the numerical data as a function of |$M^{\rm num}_{\rm disc}$|⁠. Unfilled markers represent discs for which the mass of the Keplerian disc differs from the actual one by more than 0.2. The exponent α changes considerably within our sample, from 4 up to 14, and more massive discs (⁠|$M_{\rm disc} \gt 10^{-2} \mathrm{\, M_\odot }$|⁠) have a shallower decline, characterized by 4.0 ≲ α ≲ 5.4.

Figure A2.

Power-law exponent, α, for each disc in our numerical simulation sample, as a function of the disc mass, M_disc. Unfilled markers represent discs for which the mass inside the Keplerian disc differs from the numerical one by more than 0.2. Massive discs have a shallower decline corresponding to smaller values of α′s.

The relevant parameters for the radial distributions of simulations at SR are summarized in Fig. A3. The radius of the ISCO R_ISCO (crosses), of the density peak r_peak (up-triangles), of the junction between the Gaussian and the power decay r* (stars) and of the half density peak |$r_{\sigma _{\rm max}/2}$| span a small range, indicating similar radial density distributions despite the mass spans almost three orders of magnitude. R_ISCO is found at 13–16 km from the centre, while the density peak is around 17–29 km.

Figure A3.

Fitted values of R_ISCO, r_peak, and r* as defined in equation (12) for the discs reported in A1 except the simulation with error on the mass above 0.2. Solid lines represent the radius spanned by the Gaussian, while dashed lines represent the power decay branch of σ(r) up to the radius |$r_{\sigma _{\rm max}/2}$| at which the value of the density is half of its maximum.

APPENDIX B: COMPARISON WITH THE FITTING FORMULAE USED TO COMPUTE GW190425 KILONOVA LIGHT CURVES

In this appendix, we test the fitting formulae for the ejecta and disc properties used in Raaijmakers et al. (2021) and Barbieri et al. (2021) in the parameter range of GW190425 to predict the associated kilonova light curves. Some of these formulae were originally proposed in Foucart et al. (2017), Krüger & Foucart (2020), and Radice et al. (2018b; see also Dietrich & Ujevic 2017). Additionally, we include in the comparison fitting formulae from Nedora et al. (2022) in the form of their equation (6), i.e. a second-order polynomial in the mass ratio and tidal deformability. In particular, we use coefficients fitted on the data set RefM0Set & M0/M1Set, i.e. on a set of simulations including neutrino emission and absorption, and microphysical EOSs.

We stress that we examine the different formulae in an unexplored parameter region since the binary systems within the calibration data set are overall lighter and involve more deformable objects than those in our simulations.

In Fig. B1, we compare the disc (top) and ejecta (bottom) masses predicted by the various fitting formulae with the ones obtained by our simulations. The uncertainties in the fitted values are 50 per cent of the estimated value, summed to a floor value of |$5 \times 10^{-4} \mathrm{\, M_\odot }$| for the disc mass and |$5 \times 10^{-5} \mathrm{\, M_\odot }$| for the ejecta mass. The bisector is the ‘agreement line’, while the dashed lines represent the 35 per cent deviation from the exact prediction. For the mass of the dynamical ejecta only simulations with |$M_{\rm ej}\gt 10^{-5}\mathrm{\, M_\odot }$| have been taken into account.

Figure B1.

Top: Comparison of the disc masses obtained from our numerical simulations and from the fitting formulae used in Raaijmakers et al. (2021; originally from Krüger & Foucart 2020) and in Barbieri et al. (2021). Bottom: Comparison of the dynamical ejecta masses obtained from our numerical simulations and from the fitting formulae used in Raaijmakers et al. (2021; originally from Krüger & Foucart 2020) and in Barbieri et al. (2021; originally from Radice et al. 2018b). Fitting formulae from Nedora et al. (2022) are also reported. The error bars on the vertical (horizontal) axis are estimated as the 50 per cent of the predicted value (absolute difference between the SR and LR values). For the BNS in our sample with |$M_{\rm disk}^{\rm num} \lesssim 10^{-3} \mathrm{\, M_\odot }$| (⁠|$M_{\rm disk}^{\rm num} \lesssim 10^{-4} \mathrm{\, M_\odot }$|⁠), the formulae from Krüger & Foucart (2020; Nedora et al. 2022) result in non-physical values for the disc mass.

In most of the cases, the fitting formulae significantly overestimate both the mass of the disc and the mass of the dynamical ejecta, and sometimes even predict opposite trends with respect to the binary parameters. Only in the case of the disc masses predicted by Krüger & Foucart (2020; used in Raaijmakers et al. 2021) and of the ejecta masses by Radice et al. (2018b; used in Barbieri et al. 2021) there is a partial agreement, at least within the estimated uncertainties.

The estimates of Nedora et al. (2022) is rather insensitive to the detailed binary parameters, giving rather similar ejecta mass and disc mass for each binary configuration.

Another physical input needed in kilonova light curves calculations is the velocity at which ejected matter is expelled from the binary system. In Fig. B2, we show the mass-weighted average asymptotic velocity of the dynamical ejecta obtained from our numerical simulations and from the fitting formulae presented in Radice et al. (2018b), Foucart et al. (2017), and Nedora et al. (2022). Only simulations with |$M_{\rm ej}\gt 10^{-5}\mathrm{\, M_\odot }$| have been taken into account. We assume a conservative uncertainty of the 30 per cent on the values obtained from the fitting formulae. We observe that the formulae from Radice et al. (2018b) and Nedora et al. (2022) work reasonably well for outflow speed with |$\langle v_\infty ^{\rm num}\rangle$| in the range |$0.24-0.30\, c$|⁠, while they underestimate the average velocity in the simulation with the fastest ejecta. The fitting formula from Foucart et al. (2017) used in Raaijmakers et al. (2021) to make predictions on the kilonova from the GW190425 event, but originally tailored for the dynamical ejecta of BHNS systems, predicts a very similar average velocity for all the binaries, that is systematically smaller than the outcome of the simulations. This is because the expression assumes that the average velocity of the ejecta is given by a constant value of ∼0.15 plus a linear correction in the mass ratio, which is tiny in the case of BNS systems (q ∼ 1–2).

Figure B2.

Comparison of the mass-weighted average velocity of the dynamical ejecta as obtained in our simulations and from the fitting formulae employed in the kilonova calculations of Raaijmakers et al. (2021) and Barbieri et al. (2021), taken from Foucart et al. (2017) and Radice et al. (2018b), respectively. Results from the fitting formulae from Nedora et al. (2022) are also reported. The (symmetric) uncertainties on the vertical axis are conservatively estimated as the 30 per cent of the values obtained from the fitting formulae. Error bars on the horizontal axis are estimated as the difference between the values inferred from the SR and LR simulations.

APPENDIX C: STANDARD DEVIATION OF THE AZIMUTHAL ANGLE

The azimuthal angle of the dynamical ejecta distribution ϕ_ej has a 2π-rotational symmetry. So its mass weighted SD |$\phi _{\rm ej}^{\rm SD}$| depends on an arbitrary chosen reference. For each angular bin ϕ_i of normalized weight w_i of the ejecta distribution we define the periodic shift S_δ(ϕ_i) as:

$$\begin{eqnarray} S_{\delta }(\phi _i) := \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\phi _i + \delta &{\rm if} \,\, \phi _i \lt 2\pi - \delta ~,\\ \phi _i + \delta - 2\pi &{\rm if} \,\, \phi _i \ge 2\pi - \delta ~.\\ \end{array}\right. \end{eqnarray}$$

(C1)

Let us indicate with S_δ(ϕ_ej) the distribution obtained after the shift of awl the ϕ_i. The average 〈ϕ_ej〉_δ ≡ 〈S_δ(ϕ_ej)〉 is then

$$\begin{eqnarray} \langle \phi _{\rm ej}\rangle _{\delta } = \langle \phi _{\rm ej}\rangle _0 + \delta - 2\pi W_{\delta } \, , \end{eqnarray}$$

(C2)

where W_δ is the total weight of the bins ϕ_i ≥ 2π − δ,

$$\begin{eqnarray} W_\delta = \sum _{\phi _i \ge 2\pi - \delta } w_i \le 1 \, . \end{eqnarray}$$

(C3)

We choose δ = δ* such that 〈ϕ_ej〉_δ is centred in the half of the interval, i.e in π:⁵

$$\begin{eqnarray} \delta ^{*} - 2\pi W_{\delta ^{*}} = \pi - \langle \phi _{\rm ej}\rangle _0\, . \end{eqnarray}$$

(C4)

The root mean square (RMS) of ϕ_ej after the shift S_δ is

$$\begin{eqnarray} \begin{split} {\rm RMS}_{\delta }(\phi _{\rm ej}) =\bigg [ &{\rm RMS}_0(\phi _{\rm ej})^2 + 2 \delta \langle \phi _{\rm ej}\rangle _0 + \delta ^2 + \\ &+ 4\pi \left((\pi - \delta) W_{\delta } - \overline{\langle \phi _{\rm ej}\rangle }_{\delta } \right) \bigg ]^{1/2} \, , \end{split} \end{eqnarray}$$

(C5)

where RMS₀(ϕ_ej) and 〈ϕ_ej〉₀ are the unshifted RMS and average of ϕ and |$\overline{\langle \phi _{\rm ej}\rangle }_{\delta }$| is the average of the bins ϕ_i ≥ 2π − δ,

$$\begin{eqnarray} \overline{\langle \phi _{\rm ej}\rangle }_{\delta } = \sum _{\phi _i \ge 2\pi - \delta } w_i \phi _i \, . \end{eqnarray}$$

(C6)

Finally, the SD with respect to the new average 〈ϕ_ej〉_δ is

$$\begin{eqnarray} \begin{split} \sigma _{\delta }(\phi _{\rm ej}) &= \sqrt{\sum _i w_i (S_{\delta }(\phi _i) - \langle \phi _{\rm ej}\rangle _{\delta })^2} \\ & = \sqrt{{\rm RMS}_{\delta }(\phi _{\rm ej})^2 - \langle \phi _{\rm ej}\rangle _{\delta }^2} \, . \end{split} \end{eqnarray}$$

(C7)