What if GW190425 did not produce a black hole promptly?

ABSTRACT

It is widely believed that the binary neutron star merger GW190425 produced a black hole promptly upon merger. Motivated by the potential association with the fast radio burst FRB 20190425A, which took place 2.5 h after the merger, we revisit the question of the outcome of GW190425 by means of numerical relativity simulations. We show that current laboratory and astrophysical constraints on the equation of state of dense matter do not rule out the formation of a long-lived remnant. However, the formation of a stable remnant would have produced a bright kilonova, in tension with upper limits by ZTF at the location and time of FRB 20190425A. Moreover, the ejecta would have been optically thick to radio emission for days to months, preventing a putative FRB from propagating out. The predicted dispersion measure is also several orders of magnitude larger than that observed for FRB 20190425A. Our results indicate that FRB 20190425A and GW190425 are not associated. However, we cannot completely rule out the formation of a long-lived remnant, due to the incomplete coverage of the relevant sky regions. More observations of GW190425-like events, including potential upper limit, have the potential to constrain nuclear physics. To this aim, it is important that follow-up observational campaigns of gravitational wave events are informed by the properties of the source, such as their chirp mass, and we urge the LIGO-Virgo-KAGRA collaboration to promptly release them publicly.

gravitational waves, neutron star mergers, stars: neutron

1 INTRODUCTION

GW190425 (Abbott et al. 2020, 2023) was the second neutron star (NS) merger detected by the Advanced LIGO (Aasi et al. 2015) and Virgo (Acernese et al. 2015) detectors. It had a chirp mass |$\mathcal {M} = 1.44 \pm 0.02\ {\rm M_{\odot }}$| and a total mass |$M = 3.4_{-0.1}^{+0.3}\ {\rm M_{\odot }}$|⁠, both more than 5σ larger than the mean for Galactic binary NS systems. It is not excluded that GW190425 was the merger of a light NS and a low-mass black hole (BH; Han et al. 2020), although this might be in mild tension with the lack of a kilonova associated with this event (Coughlin et al. 2019; Antier et al. 2020; Kyutoku et al. 2020; Kilpatrick et al. 2021; Rastinejad et al. 2022). Assuming that both objects were slowly spinning restricts the mass ratio and the total mass to q = 1 − 1.25 and M = 3.2 − 3.4 M_⊙ (Abbott et al. 2020). These values are compatible with the hypothesis that there were two NSs in GW190425. Assuming that both objects in the binary were NSs, Foley et al. (2020) argued that the total mass of GW190425 can be restricted even further to |$M = 3.33_{-0.06}^{+0.10}\ {\rm M_{\odot }}$| and the mass ratio |$q = M_1/M_2 = 1.27^{+0.40}_{-0.27}$|⁠.

It is commonly believed that GW190425 produced a BH promptly after merger. Indeed, Abbott et al. (2020) estimate that the remnant underwent prompt collapse with 96 per cent probability with the method proposed by Agathos et al. (2020). Furthermore, numerical relativity simulations targeted to GW190425 (Camilletti et al. 2022; Dudi et al. 2022) found prompt BH formation in all considered cases. Perhaps unexpectedly (Barbieri et al. 2021; Raaijmakers et al. 2021), simulations, particularly those with realistic microphysics (Camilletti et al. 2022), predicted small ejecta masses and faint kilonova for most configurations. As such, the observational upper limits turned out to be only weakly constraining on the binary parameters. It should also be remarked that not all of the plausible sky locations for GW190425, which spanned ∼10⁴ deg² of the celestial sphere, have been observed.

A possible short γ-ray burst (SGRB) coincident in time with GW190425 was identified (Pozanenko et al. 2020), but further analysis revealed that its likely origin is a magnetar flare in a nearby galaxy (NGC 253; Minaev & Pozanenko 2020). Searches for SGRB afterglow emission in GW190425 also did not identify any candidate (Boersma et al. 2021).

Moroianu et al. (2023) discussed the possible association between GW190425 and the fast radio burst (FRB) 20190425A. FRBs are millisecond duration pulses of bright radio emission that are located at cosmological distances (Lorimer et al. 2007; Thornton et al. 2013; Petroff, Hessels & Lorimer 2022). Although a substantial sub-population of repeating FRBs has been detected (see e.g. Spitler et al. 2016; Fonseca et al. 2020; Bhardwaj et al. 2021; Niu et al. 2022), cataclysmic origins for the larger sample of FRBs that do not repeat has not been ruled out. Among many FRB source models within the realm of such cataclysmic events, BNS mergers remain as a potential source of radio emission at different stages of evolution of the central remnant (see Platts et al. 2019). FRB 20190425A occurred 2.5 h after the merger in the gravitational wave (GW) sky localization area. Moroianu et al. (2023) estimate the probability of a chance coincidence to be 0.0052 (corresponding to 2.8 σ). A possible interpretation of this signal can be made in the context of the ‘blitzar’ mechanism of FRB (Falcke & Rezzolla 2014; Most, Nathanail & Rezzolla 2018). Accordingly, coherent radio emission can be launched by the snapping of the field lines of a NS when it suddenly collapses to BH. Such interpretation would imply not only that GW190425 did not promptly produce a BH, but also that the remnant survived for 2.5 h before collapsing to BH, which would have profound implications for the physics of dense matter (Fattoyev et al. 2020; Godzieba, Radice & Bernuzzi 2021; Lim et al. 2021; Tews et al. 2021). On the other hand, Bhardwaj et al. (2023) argued that the radio waves produced by the collapse of the remnant in GW190425 would not be able to propagate through the merger ejecta. This would exclude an association between GW190425 and FRB 20190425A, unless the merger produced a very small amount of ejecta.

In this paper, we revisit the question of the fate of GW190425. We perform numerical relativity NS merger simulations targeted to GW190425 using the ‘Big Apple’ (BA) relativistic mean-field theory equation of state (EOS) (Fattoyev et al. 2020). This EOS satisfies many extant laboratory and astronomical constraints, yet predicts the formation of dynamically stable remnant massive neutron star (RMNS) remnants for GW190425, which could conceivably be supported over a time-scale of hours after the merger. However, we show that such scenario produces a bright kilonova that would have been detected at the distance and location of FRB 20190425A. Moreover, our simulations also show that the column density of the ejecta emitted during the formation and the evolution of long-lived remnants would have been too dense to be traversable by a putative FRB pulse. Our simulations thus exclude an association between FRB 20190425A and GW190425. They also disfavour the BA EOS, with the caveat that, since not all the sky region associated with GW190425 was observed, GW190425 might have produced a bright kilonova that escaped detection. All in all, our work makes the case for continued effort to follow up GW190425-like events electromagnetically.

The rest of this paper is organized as follows. In Section 2 we present the merger simulation results. The associated kilonova signal is discussed in Section 3, while the properties of a potential FRB are discussed in Section 4. Finally, Section 5 is dedicated to discussion and conclusions.

2 MERGER SIMULATIONS

For our simulations we use the BA EOS (Fattoyev et al. 2020), which predicts maximum non-rotating NS mass of 2.6 M_⊙ and radius of a non-rotating 1.4 M_⊙ NS of R_1.4 = 12.96 km. We construct irrotational initial data with the lorene pseudo-spectral code (Gourgoulhon et al. 2001) at the initial separation of 45 km, corresponding to the last ∼4 orbits of the binary. We fix the chirp mass of the system to |$\mathcal {M} = 1.44\ {\rm M_{\odot }}$| and consider four mass ratios q = 1, 1.17, 1.33, and 1.67. The properties of the initial data are summarized in Table 1.

Table 1.

Mass of the each binary component M₁ and M₂, total binary baryonic mass M_b, reduced tidal parameter |$\tilde{\Lambda }$|⁠, disc and ejecta masses M_disc and M_ej at the end of the simulations, and root-mean-squared opening angle |$\sqrt{\langle (90^\circ - \theta)^2 \rangle _{\rm ej}}$| and velocity |$\sqrt{\langle v^2 \rangle _{\rm ej}}$| of the ejecta for each binary.

Model	M1	M2	Mb	|$\tilde{\Lambda }$|	Mdisc	Mej	|$\sqrt{\langle (90^\circ - \theta)^2 \rangle _{\rm ej}}$|	|$\sqrt{\langle v^2 \rangle _{\rm ej}}$|
	[M⊙]	[M⊙]	[M⊙]		[M⊙]	[M⊙]	(deg)	[c]
q1.00-LR	1.654	1.654	3.700	286	0.28	0.0120	34	0.23
q1.00-SR	1.654	1.654	3.700	286	0.18	0.0071	34	0.22
q1.17-LR	1.795	1.527	3.721	283	0.26	0.0054	33	0.23
q1.17-SR	1.795	1.527	3.721	283	0.20	0.0049	32	0.26
q1.33-LR	1.914	1.437	3.766	278	0.25	0.0128	25	0.22
q1.33-SR	1.914	1.437	3.766	278	0.25	0.0136	29	0.22
q1.67-SR	2.149	1.289	3.904	217	0.12	0.0070	13	0.22

We evolve the initial data using the general-relativistic hydrodynamics code whiskythc (Radice & Rezzolla 2012; Radice, Rezzolla & Galeazzi 2014b,a, 2015), which is based on the Einstein Toolkit (Loffler et al. 2012; Etienne et al. 2021). For the simulations discussed here, we use the carpet adaptive mesh-refinement (AMR) driver (Schnetter, Hawley & Hawke 2004; Reisswig et al. 2013a), which implements the Berger-Oilger scheme with refluxing (Berger & Oliger 1984; Berger & Colella 1989). All simulations are performed at two resolutions: LR (corresponding to Δx_min = 0.167GM_⊙/c² ≃ 246 m) and SR (Δx_min = 0.125GM_⊙/c² ≃ 185 m). However, we do not report the results of the q1.67-LR simulation, since it failed shortly after BH formation. We evolve the space--time geometry using the CTGamma code (Pollney et al. 2011; Reisswig et al. 2013b), which solves the Z4c formulation of Einstein’s equations (Bernuzzi & Hilditch 2010; Hilditch et al. 2013). The standard first-order moving puncture gauge is employed (van Meter et al. 2006). The zero-temperature BA EOS is augmented with a ideal-gas thermal component with adiabatic index Γ_th = 1.7, following e.g. Shibata, Taniguchi & Uryu (2005), Bauswein, Janka & Oechslin (2010), Endrizzi et al. (2018), and Figura et al. (2020).

The BA EOS not only predicts a 2.6 M_⊙ non-rotating maximum NS but also agrees with a host of both nuclear and astrophysical observables. The parameters of this model are tuned to reproduce the ground-state properties of finite nuclei, such as binding energies and charge radii. The BA EOS also successfully predicts the bulk parameters of nuclear matter in agreement with bulk parameters constrained from both the giant monopole resonances and neutron skins. The soft density dependence of symmetry energy with L = 39.8 MeV allows the model to predict the 1.4 M_⊙ NS with radius of 12.96 km and dimension-less tidal polarizability of 717.3, consistent with astrophysical observations (Abbott et al. 2018; Miller et al. 2019; Riley et al. 2019; Breschi et al. 2021). The required stiffening of the BA EOS at high densities allows the model to predict a maximum non-rotating NS mass of 2.6 M_⊙, which is in mild tension with the existing constraints obtained from energetic heavy-ion collisions (Danielewicz, Lacey & Lynch 2002), but it is still viable within 2σ. In fact, even more extreme models can be produced achieving maximum non-rotating NS masses in excess of 2.8 M_⊙ as shown by Mueller & Serot (1996) and Fattoyev et al. (2020).

The maximum baryonic mass for a uniformly rotating NS predicted by the BA EOS is 3.834 M_⊙. As can be seen in Table 1, all binaries with the exception of q = 1.67, have baryonic mass below this value. As such, they are expected to form supermassive remnants that are stable on secular time-scales (Cook, Shapiro & Teukolsky 1992; Baumgarte, Shapiro & Shibata 2000; Radice et al. 2018a; Radice & Bernuzzi 2023). This expectation is confirmed by the simulations. In Fig. 1 we plot the maximum density for all SR simulations. After a period of violent oscillations, the merger remnants for the q = 1, q = 1.17, and q = 1.33 settle to a quasi-steady configuration comprising a massive RMNS surrounded by a thick accretion disc. In contrast, the q = 1.67 binary undergoes prompt BH formation upon merger.

Figure 1.

Maximum density for all SR binaries as a function of time. With the exception of the q1.67-SR binary, all models produce stable remnants.

Fig. 2 shows the dynamics of the material in the orbital plane for the q1.00-SR binary, which is representative of all the stable remnant cases. This should be contrasted with Fig. 3 which shows the dynamics for the q1.67-SR. In this case, the secondary NS is tidally disrupted during its last orbit. The primary NS collapses to form a BH as a result of the accretion from the tidal stream. A similar dynamics has also been reported in Bernuzzi et al. (2020) for a lower mass system with a different EOS.

Figure 2.

Rest-mass density on the orbital plane for the q1.00-SR binary. The merger results in the formation of a long-lived remnant.

Figure 3.

Rest-mass density on the orbital plane for the q1.67-SR binary. The black contour encloses the region with α < 0.2, which is expected to be located within the apparent horizon. During the merger the secondary star is disrupted and the primary collapses to BH. The outcome of the merger is a BH surrounded by a thick accretion disc.

All binaries produce massive ejecta, driven by tidal torques and shocks during the merger, including the prompt collapse q = 1.67 binary. We estimate the mass ejection by integrating the flux of gravitationally unbound matter (with u_t < −1, u being the fluid four-velocity) crossing a coordinate sphere of radius |$200\, G/c^2\, M_\odot \simeq 295\ {\rm km}$|⁠. The results are reported in Table 1 and Fig. 4 (left panel). The dynamical ejecta mass ranges from 0.49 × 10⁻² M_⊙ to 1.36 × 10⁻² M_⊙. Because dynamical ejecta are driven by a combination of shock heating, which weakens as the mass ratio increases, and tidal torques, which has the opposite mass ratio dependence, the dynamical ejecta mass is not monotonic with q (Dietrich et al. 2017). Finite-resolution errors are larger for the q = 1 binary, since shocks are more challenging to capture numerically. That said, we caution the reader that our simulations might be overestimating the ejecta mass, since we are neglecting finite-temperature and neutrino effects (Nedora et al. 2022). However, the resulting uncertainty is subdominant compared to that in the mass of the secular ejecta from the merger disc.

Figure 4.

Dynamical ejecta mass (left panel) and remnant disc mass (right panel) as a function of time for all the SR models. All systems produce massive outflows and result in the formation of very massive discs.

All binaries also result in the formation of accretion discs. To quantify this, we compute the total mass of material with rest-mass density ρ < 10¹³ g cm⁻³, which we assume to be belonging to the disc, following Shibata et al. (2017) and Radice et al. (2018b). This is motivated by the fact that this density roughly corresponds to the transition to a quasi-Keplerian rotational profile in the remnant (Hanauske et al. 2017). For the q = 1.67 binary, we also exclude material with lapse function below 0.2, which would be located inside the BH (Bernuzzi et al. 2020). The results are reported in Table 1 and Fig. 4 (right panel). We find that all binaries form very massive discs, up to 0.28 M_⊙. It is expected that |$20{-}50{{\ \rm per\ cent}}$| of the discs will evaporate due to neutrino-driven winds, magnetic torques, and nuclear recombination (Metzger & Fernández 2014; Fahlman & Fernández 2018; Radice et al. 2018a; Nedora et al. 2019; Fujibayashi et al. 2020; Just et al. 2023). As such, this secular component of the ejecta is expected to be dominant. We remark that the disc is still rapidly accreting onto the BH for the q1.67-SR binary at the time when we stop the simulation. This means that the disc mass we report for q1.67-SR should be taken as an upper bound. That said, because this binary underwent prompt collapse, it is not a suitable candidate for FRB 20190425A. For this reason, we will focus on the q = 1, q = 1.17, and q = 1.33 binaries for the rest of the discussion.

3 KILONOVA LIGHT CURVES

Based on the available ejecta properties extracted from the simulations, we compute synthetic, broad-band kilonova light curves during the first few days after the merger. We assume the kilonova emission to be powered by the radioactive decay of freshly synthesized r-process elements in the ejecta, while we neglect possibly spin-down energy injection which would have resulted in a merger nova that can be significantly brighter than kilonovae (Yu, Zhang & Gao 2013; Metzger & Piro 2014; Kisaka, Ioka & Nakar 2016; Murase et al. 2018; Ai, Zhang & Zhu 2022; Sarin et al. 2022; Wang, Beniamini & Giannios 2024). Thus, we consider our prediction as a robust lower limit for the electromagnetic counterpart. We employ xkn, a multicomponent, anisotropic kilonova model, presented in Ricigliano et al. (2023) and based on the framework first described and implemented in Martin et al. (2015) and Perego, Radice & Bernuzzi (2017). In contrast to the original version, xkn computes the light curves starting from a solution of the diffusion equation for photon radiation emitted at the photosphere, obtained using a semi-analytical approach first proposed by Wollaeger et al. (2018). The model employs a self-similar homologous density profile, which is expected to well describe the ejecta profile at the time of the kilonova emission (see e.g. Neuweiler et al. 2023). Ricigliano et al. (2023) extended this solution by including composition and time-dependent radioactive heating rates, thermalization efficiencies and opacities, and by accounting for the contribution to the emission coming from the external, optically thin layers of the ejecta. xkn requires as input of each component the total ejecta mass, the average expansion velocity, and the grey opacity, as well as the angular distribution of these quantities. Moreover, the model uses heating rates that are parametrized in terms of the three quantities that determine the nucleosynthesis in the homologously expanding ejecta, namely the electron fraction, the specific entropy, and the expansion time-scale (Hoffman, Woosley & Qian 1997). This approach was also adopted in Wu et al. (2022) (see e.g. Lippuner & Roberts 2015; Setzer et al. 2023, for similar fits.)

We consider two axisymmetric ejecta components for our kilonova model: a faster, external dynamical component, and a slower, inner wind one. The angular profiles of the ejecta (mass fraction with angle, electron fraction, and entropy) are only weakly dependent on the EOS, but are strongly impacted by the neutrino treatment used in the simulations (Radice et al. 2016; Perego, Radice & Bernuzzi 2017). Given that our simulations do not implement neutrino transport, we take values based on the results of neutrino-radiation simulations from the literature. Specifically, regarding the dynamical ejecta, the values of its total mass, M_ej, and expansion velocity, |$\sqrt{\langle v^2 \rangle }_{\rm ej}$|⁠, are the same as that of the simulated dynamical ejecta, as reported in Table 1. Since the opening angle of the expanding ejecta, |$\sqrt{\langle \theta ^2 \rangle }_{\rm ej}$|⁠, is between 25 and 34 degrees, we consider a broad angular distribution in the mass, characterized by a sin θ dependence on the polar angle θ, measured from the orbital axis of the binary. As a consequence of the more intense neutrino irradiation inside the polar regions, we expect an electron fraction Y_e ≳ 0.25 in this region, so we set the grey photon opacity to κ_dyn = 1 cm² g⁻¹ for θ ≲ 45°. In the region across the equator, for θ ≳ 45°, we expect lower Y_e and set κ_dyn = 15 cm² g⁻¹. Such a value is expected if a non-negligible amount of lanthanides are synthesized inside the ejecta (see e.g. Tanaka et al. 2020, for typical grey opacity values of r-process enriched material). We assume the average velocity of the ejecta to be constant across different polar angles, and we set its average specific entropy to s = 10 k_B baryon⁻¹, while we compute the ejecta expansion time-scale as |$\tau _{\rm exp}= 1~{\rm ms} \left(c/\sqrt{\langle v^2 \rangle _{\rm ej}} \right)$|⁠. For the second ejecta component, we cannot rely on detailed ejecta properties extracted from our simulations, since the latter are not long enough to model the disc evolution. We consider two possibilities: a spherical secular wind, or an aspherical spiral wave wind. In both cases, we assume the ejecta mass to vary between 20 per cent and 40 per cent of the disc mass and we adopt again a sin θ matter distribution in the case of the aspherical spiral wave wind. In both instances, we set a uniform entropy s = 20 k_B baryon⁻¹, and an isotropic average velocity, with lower values relative to the dynamical ejecta, i.e. |$\sqrt{\langle v^2 \rangle _{\rm ej}}=0.06~c$| for the secular wind (as expected from the nuclear recombination of matter inside the disc and consistent with AT2017gfo modelling, see e.g. Villar et al. 2017; Perego, Radice & Bernuzzi 2017), and |$\sqrt{\langle v^2 \rangle _{\rm ej}}=0.15~c$| for the spiral wave wind (as obtained by detailed spiral wave wind simulations in Nedora et al. 2019, 2022). Due to the presence of an RMNS, the photon opacity of this ejecta component is expected to be smaller than that of the equatorial part of the dynamical ejecta (see e.g. Lippuner et al. 2017). Therefore, in the case of the secular wind, we set the grey opacity to κ_wind = 5 cm² g⁻¹, while for the spiral wave wind we assume κ_wind = 1 cm² g⁻¹ for θ ≲ 45°, and κ_wind = 5 cm² g⁻¹ for θ ≳ 45°. Finally, we consider a source luminosity distance D_L varying between 141.8 and 142.7 Mpc, correspondent to the value and uncertainty in the redshift of the FRB host galaxy UGC10667 (Abazajian et al. 2009), together with an inclination angle θ_view between 46° and 70°, as inferred by Bhardwaj et al. (2023).

Fig. 5 summarizes the results of the kilonova emission, with the uncertainty in the amount of disc material expelled, the structure of the disc wind component, the source distance, and inclination angle being represented by the width of the bands. The colour curves are compared to the observational upper limits in the r and g bands taken from https://treasuremap.space/ (Wyatt et al. 2020), imposed by the Zwicky Transient Facility (ZTF; Graham et al. 2019), the Gravitational-wave Optical Transient Observer (GOTO; O’Brien 2018) and the Searches After Gravitational waves Using ARizona Observatories (SAGUARO/CSS; Christensen et al. 2012; Lundquist et al. 2019; Paterson et al. 2021).¹ In the latter case, all images were taken without a filter and were calibrated to Gaia G-band using Gaia DR2 (Brown et al. 2018). We compare the resulting upper limits for the magnitudes computed both in the g and r band. While limits from GOTO and CSS do not place any significant constraint on the model, we find that the kilonova emission, especially in its redder component, appears too bright to be compatible with the ZTF limits before and around 1 d, regardless of the different simulated binary merger. This is due to the relatively high amount of material expelled by the disc during the RMNS lifetime, which exceeds by about one order of magnitude the amount of dynamical ejecta. We note that, in the context of the kilonova model, the considered sources of uncertainty only place up to one magnitude of variation on the overall emission. In order to reconcile the model with the ZTF limits, we would need to invoke considerably larger uncertainties. The latter can be potentially found in the modelling of the nuclear heating and of the ejecta opacity. Different nuclear inputs for the nucleosynthesis calculations lead to variations of up to ∼1 order of magnitude on the heating rate. The opacities are even more uncertain, due to the incompleteness of atomic data (see Ricigliano et al. 2023, for a more detailed discussion of the uncertainties of our kilonova model). However, here we absorb the ignorance on the ejecta opacity by considering different ejecta setups. We note that the heating rate uncertainty can be translated in an additional ∼0.25 mag variation in our light curves. Even so, our results would barely meet the ZTF constraints and, as such, they tend to disfavour the formation of a long lived remnant at the time of GW190425 and sky location of FRB 20190425A. In absence of an association with FRB 20190425A, we can only exclude the formation of a long-lived remnant if the distance to the source was D_L ≲ 150 Mpc and if the merger happened in an area of the sky that has been surveyed down to ∼21 magnitudes. On the other hand, the distance inferred from the GW alone is |$159^{+69}_{-71}\ {\rm Mpc}$| (Abbott et al. 2020), and therefore the current data are not strongly constraining.

Figure 5.

Synthetic kilonova colour curves in the r and g bands, for a source at a distance D_L ≃ 142 Mpc and inclination angle θ_view ≃ 56°, as obtained using as input the outcome of the q = 1, q = 1.17, and q = 1.33 simulations, together with observational upper limits from ZTF, GOTO, and SAGUARO/CSS. The SAGUARO/CSS upper limit was taken without a filter and we compare it with the magnitudes computed both for the g and r filters. The width of each band encodes the estimated uncertainties in the ejecta modelling and source properties (see the main text for more details). The upper limits in the r band from ZTF appears to rule out the formation of a stable RMNS in GW190425, assuming to know the sky position as the one of FRB 20190425A.

4 FAST RADIO BURSTS

A fast-rotating pulsar/magnetar that can be left as a merger remnant may power not only FRB emission but also non-thermal emission from pulsar/magnetar wind nebulae. A significant fraction of the spin-down and/or magnetic energy can be extracted by an outgoing relativistic wind and carried by non-thermal electrons and positrons via magnetic energy dissipation, leading to synchrotron and inverse-Compton emission. This scenario is widely considered to explain multiwavelength spectra of Galactic pulsar wind nebulae (Gaensler & Slane 2006, for a review). Quasi-steady synchrotron emission from embryonic nebulae (Murase, Kashiyama & Mészáros 2016), can explain persistent radio counterparts of FRB 121102 and FRB 190520 (Chatterjee et al. 2017; Niu et al. 2022). Such synchrotron nebular emission from a pulsar/magnetar has also been studied in the context of merger remnants and the accretion-induced collapse of white dwarfs (Kashiyama & Murase 2017; Murase et al. 2018; Yamasaki, Totani & Kiuchi 2018).

The fact that FRB 20190425A occurred 2.5 h after GW190425 motivates us to consider FRB emission from a merger remnant in the post-merger phase. However, radio waves can be significantly depleted via various absorption processes, including free–free absorption in the ejecta and synchrotron absorption in the nebula (Murase, Meszaros & Fox 2017; Bhardwaj et al. 2023). The dense ejecta can also make a significant contribution to the observed dispersion measure (DM), which can be used as an important constraint on the model.

Assuming UGC 10667 at redshift z = 0.03136 as the host for FRB 190425, the DM inferred from the FRB can be written in terms of five primary components, |${\rm DM}_{\rm obs} \simeq {\rm DM}_{\rm MW} + {\rm DM}_{\rm halo} + {\rm DM}_{\rm IGM} + {\rm DM}_{\rm host}+{\rm DM}_{\rm src} = 128\, {\rm pc\, cm^{-3}}$|⁠, where DM_MW/halo is the Milky Way interstellar medium (ISM)/halo contribution, DM_IGM is the intergalactic medium (IGM) contribution and DM_host is the FRB host galaxy contribution including its halo. Here DM_src = DM_ej/w + DM_nb includes contributions from both the merger ejecta/wind and magnetized nebula. We mainly focus on the merger ejecta, since their contribution is dominant. However, we give an estimate for the optical depth of the pair nebula below. From theoretical and data-driven estimates of DM_MW, DM_halo, and DM_IGM, using models for the average electron number density along the FRB line-of-sight (see e.g. Cordes & Lazio 2002, 2003; Ade et al. 2016; Prochaska & Zheng 2019; Platts, Prochaska & Law 2020), we infer |${\rm DM}_{\rm host}+{\rm DM}_{\rm src} \lesssim 30\, {\rm pc\, cm^{-3}}$|⁠. Assuming a negligible host galaxy contribution to DM_obs, we numerically compute the number density of free electrons in the ejecta/wind as |$n_{\rm e,ej/w} \approx 3M_{\rm ej/w}/(4\pi R_{\rm ej/w}^3 \mu _e m_{\rm H})$|⁠, to ascertain DM_src ≈ n_{e, ej}R_ej + n_{e, w}R_w. Here M_ej/w is the ejecta/wind mass, R_ej/w is the ejecta/wind radius, and 1/μ_e corresponds to the number of free electrons per baryon for an electron temperature of |$\mathcal {T}_e \sim 10^4\, {\rm K}$| (Hamilton & Sarazin 1984).

We compute the time evolution of DM, in which for the ejecta we assume a spherically expanding distribution of mass until the NS collapse time |$t\simeq 2.5\, {\rm h}$|⁠. For the magnetar, we set the NS initial spin period |$P_i \sim 1\, {\rm ms}$|⁠, assuming that the remnant will achieve solid body rotation with a spin close to the mass-shedding limit (Radice et al. 2018a; Radice & Bernuzzi 2023). The surface dipole field is taken to be |$B_{\rm dip}\sim 2\times 10^{14}\, {\rm G}$|⁠, so that the spindown time-scale is consistent with the 2.5 h delay between GW190425 and FRB 20190425A. The interior magnetic field is |$B_{\rm int}\sim 10^{16}\, {\rm G}$| and is predominantly toroidal. The radii of the merger ejecta and magnetar-powered nebula are then calculated following the merger nova model presented in Murase et al. (2018). Due to a faster radius evolution for the ejecta and wind with an increase in dipolar field strength B_dip (and for fixed |$B_{\rm int}\sim 10^{16}\, {\rm G}$|⁠), DM_src ≈ n_{e, ej}R_ej + n_{e, w}R_w decreases as the corresponding electron number densities are smaller.

The unshaded area in Fig. 6 shows the time evolution of DM_src prior to the collapse of the remnant NS to a BH. We consider separately the contribution of dynamical ejecta and disc wind. As expected, the dominant contribution to DM_src arises from the dense merger ejecta. We estimate the dynamical ejecta mass (M_ej) and disc wind mass (M_w) from the results obtained in Table 1 and Fig. 4. For the dynamical ejecta, we consider a fiducial mass |$M_{\rm ej}=0.01\, M_{\odot }$|⁠, expansion velocity |$v_{\rm ej}=0.25\, c$|⁠, energy |$\mathcal {E}_{\rm ej}=3\times 10^{50}\, {\rm erg}$| and initial electron fraction Y_{e, ej} ∼ 0.15 corresponding to mean atomic mass number |$\bar{A}\simeq 166$|⁠, average charge |$\bar{Z}\simeq 65$| and |$1/\mu _e\sim 3.6/\bar{A}$| (at 2.5 h). The corresponding parameters at 2.5 h for the disc wind are |$M_{\rm w}=0.1\, M_{\odot }$|⁠, |$v_{\rm w}=0.1\, c$|⁠, |$\mathcal {E}_{\rm w}=10^{51}\, {\rm erg}$|⁠, and Y_{e, w} ∼ 0.3, resulting in |$\bar{A}\simeq 81$|⁠, |$\bar{Z}\simeq 34$| and |$1/\mu _e\sim 2.4/\bar{A}$|⁠. The average atomic mass number and charges have been obtained by performing single-trajectory nucleosynthesis calculation with the SkyNet nuclear reaction network code (Lippuner & Roberts 2017). The average ionization fraction is obtained assuming local thermodynamic equilibrium.

Figure 6.

Time evolution of the DM contribution from dynamical ejecta and disc wind are shown. The representative ejecta masses are obtained from Table 1 and Fig. 4. The dashed lines for both components correspond to mass variation by a factor of two. The horizontal grey-shaded area bound by the solid black line shows the allowed region where |${\rm DM}_{\rm src} \lesssim 30\, {\rm pc\, cm^{-3}}$| (assuming DM_host ≈ 0), as inferred from the radio observations. The vertical grey-shaded area (⁠|$t \gtrsim 2.5\, {\rm h}$|⁠) denotes the region disallowed by the collapse of remnant NS to a BH.

We stress that our constraint from the observed DM is conservative. It should be noted that the assumption that the ionization fraction is given by the Saha equation at 2.5 h underestimates the ejecta DM and opacity. In particular, the ejecta is expected to be hotter and potentially more ionized at early times. X-ray and gamma-ray photons may diffuse out from the nebula and increase the opacity in the radio band (Metzger & Piro 2014; Murase et al. 2018; Wang, Beniamini & Giannios 2024). For our ejecta parameters, we find that DM_src significantly overestimates the inferred near-source contribution |${\rm DM}_{\rm src} \lesssim 30\, {\rm pc\, cm^{-3}}$|⁠, even when assuming DM_host ≈ 0.

DM can also be larger due to free electron–positron pairs that can be created in the magnetized nebula or NS magnetosphere through electromagnetic cascades. For a spin-down luminosity L_sd ∼ 3 × 10⁴⁷ erg s⁻¹ and nebular radius R_nb ∼ 10¹⁴ cm at 2.5 h, the compactness parameter is estimated to be l_e ∼ L_sdσ_T/(R_nbm_ec³) ∼ 8 × 10⁴. The resulting Thompson optical depth is |$\tau _T\approx n_{\pm ,\rm nb}R_{\rm nb}\sigma _T \sim {({l_e {\rm PY}})}^{1/2}\sim 90$|⁠, where PY ∼ 0.1 is the pair yield (Lightman & Zdziarski 1987), implying that the dispersion measure from the nebula, |${\rm DM}_{\rm nb}\approx n_{\pm ,\rm nb}R_{\rm nb}$|⁠, would also exceed DM_obs. Although we assume that the FRB emission is produced near the RMNS surface at |$R \sim 10^{6}\, {\rm cm}$|⁠, this does not affect our constraints on |${\rm DM}_{\rm src} \lesssim 30\, {\rm pc\, cm^{-3}}$|⁠. This is expected as most of the DM is accumulated either in the merger ejecta with |$R_{\rm ej} \sim 10^{13-14}\, {\rm cm}$| (baryon dominated) or in the magnetized nebula (pair dominated) at |$R_{\rm ej} \sim 10^{11-13}\, {\rm cm}$|⁠, both with a significantly larger radius compared to the remnant.

For FRB pulses that are produced at inner regions to be observed, both of the ejecta and nebula should be optically thin to scattering or absorption effects that can suppress the radio signal. The optical depth for free–free absorption in the merger ejecta is given by (Murase, Meszaros & Fox 2017)

$$\begin{eqnarray} \tau _{\rm ff}\approx 2.1\times 10^{-25}\, \mathcal {T}_{e,4}^{-1.35} \nu _{9}^{-2.1} \int {\rm d}r\, n_{\rm e,ej}\, n_{\rm i,ej}\, \bar{Z}^2\, , \end{eqnarray}$$

(1)

and should be sufficiently small. Here n_{i, ej} is the number density of ions in the ejecta and |$\bar{Z}$| is the average charge of the nuclei. Similarly, the optical depth for synchrotron absorption in the nebula,

$$\begin{eqnarray} \tau _{\rm sa} = R_{\rm nb} \int {\rm d}\mathcal {E}_e\, \sigma _{\rm sa}(\nu , \mathcal {E}_e) \frac{{\rm d}n_{\mathcal {E}_e}}{{\rm d}\mathcal {E}_e}\, , \end{eqnarray}$$

(2)

should also not exceed unity for a given energy distribution of electrons |${\rm d}n_{\mathcal {E}_e}/{\rm d}\mathcal {E}_e$| (see Murase, Kashiyama & Mészáros 2016 for more details) and an energy-dependent SSA cross section σ_sa (Ghisellini & Svensson 1991).

For dynamical ejecta with |$M_{\rm ej} \sim 0.01\, M_{\odot }$|⁠, the optical depths for free–free absorption τ_ff ∼ 2.3 × 10¹⁰ and synchrotron absorption τ_sa ∼ 1.6 × 10¹¹ both significantly exceed unity at |$t \simeq 2.5\, {\rm h}$|⁠, when the remnant NS undergoes collapse to a BH. Similarly, for the disc wind with |$M_{\rm w} \sim 0.1\, M_{\odot }$|⁠, the corresponding optical depths are τ_ff ∼ 8.1 × 10¹³ and τ_sa ∼ 7.1 × 10¹⁰, respectively. These demonstrate that both dynamical ejecta and disc winds are too dense at |$t \lesssim 2.5\, {\rm h}$| which results in significant attenuation of the radio signal, thereby excluding an association between GW190425 and FRB 190425.

Although our results qualitatively agree with those of Bhardwaj et al. (2023), it should be noted that the values of DM_src derived in our model are significantly smaller. While they obtain |$M_{\rm ej} \lesssim 10^{-8}\, M_{\odot }$| to satisfy the observed DM constraint, we infer |$M_{\rm ej/w} \lesssim 10^{-7}\, M_{\odot }$| which is less restrictive. This is due to the fact that Bhardwaj et al. (2023) assume the full ionization of the ejecta nuclei in their analysis, which significantly overestimates the number of free electrons, and therefore their integrated column density. In our model, we find that |$1/\mu _e\sim (2-3)\bar{A}^{-1}$| for both the dynamical ejecta and disc wind, which is a more realistic assumption. In addition, Bhardwaj et al. (2023) obtain |$M_{\rm ej} \lesssim 7\times 10^{-15}\, M_{\odot }$| for the radio signal to be not attenuated due to free–free absorption. From our model, we have a less stringent limit |$M_{\rm ej/w} \lesssim 10^{-8}\, M_{\odot }$| for the optical depth of free–free absorption to be less than unity.

5 CONCLUSIONS

Motivated by the potential association between GW190425 and FRB 20190425A (Moroianu et al. 2023), we have considered a scenario in which the NS merger remnant survived over secular time-scale. We have performed numerical relativity simulations targeted to GW190425 with the BA EOS (Fattoyev et al. 2020), which satisfies most presently available astronomical and laboratory constraints on the dense matter EOS and predicts a maximum mass for a non-rotating NS of 2.6 M_⊙. The chirp mass for all simulations has been fixed to |$\mathcal {M} = 1.44\ {\rm M_{\odot }}$|⁠, and we have explored four mass ratios: q = 1, 1.17, 1.33, and 1.67. With the exception of the q = 1.67 binary, all other binaries produce secularly stable remnants.

We have found that all binaries produce massive outflows with M_ej ≲ 0.01 M_⊙ through a combination of tidal torques and shocks at merger. Even the q = 1.67 binary, which undergoes prompt BH formation, produces a massive outflow as a result of the tidal disruption of the secondary star. The dynamics for this latter binary are similar to that discussed in Bernuzzi et al. (2020) for a lower mass system with a softer EOS. In addition to the dynamical ejecta, all systems produce remnants with massive accretion discs (M_disc = 0.1 M_⊙ − 0.3 M_⊙). It is expected that |$20{-}40{{\ \rm per\ cent}}$| of the disc will also be unbound on secular time-scale, further contributing to the mass ejection (Shibata & Hotokezaka 2019).

Due to the large mass ejection, our models produce bright kilonovae. At the distance inferred for FRB 20190425A, the kilonova is expected to peak at 20 mag in both g and r bands within 1 d after merger, for the q = 1, 1.17, and 1.33 binaries, which form long-lived remnants. There are several stringent upper limits from ZTF, GOTO, and CSS at the sky position of FRB 20190425A in the first two days. The three r-band ZTF upper limits of 21 mag in the first day post-merger are particularly constraining and disfavour a scenario in which FRB 20190425A is triggered by the delayed collapse of the merger remnant. The formation of a stable remnant, and hence the BA EOS, is disfavoured even if we discard the association with FRB 20190425A, due to the non-detection of any kilonova despite the extensive search (Coughlin et al. 2019). However, it is important to keep in mind that not all the plausible sky locations for GW190425 have been observed and our kilonova model possibly suffers from systematic uncertainties difficult to quantify with accuracy. Therefore it is not possible to definitively exclude the formation of a long-lived remnant for GW190425 from the lack of a robust kilonova identification. The kilonova would also have escaped detection if the merger was at more than ∼150 Mpc, which is not unlikely given that the inferred distance for GW190425 from GW data alone was |$159^{+69}_{-71}\ {\rm Mpc}$| (Abbott et al. 2020).

The presence of a dense shell of ejecta between us and the compact object formed in GW190425 also affects the propagation of a possible radio burst generated by the remnant. The requirement that the medium should be transparent to radio waves at the relevant frequencies translates to a very stringent constraint on the ejecta mass. The observed DM of the FRB poses an additional constraint. We have revised the estimates of Bhardwaj et al. (2023) using the data from our simulations. We have also taken into account the fact that heavy elements have a very large ionization barrier, so they are not fully ionized, but only singly or doubly ionized at the relevant times. Accordingly, we find significantly smaller DMs and absorption opacities. While Bhardwaj et al. (2023) placed an upper limit on the ejecta mass of 7 × 10⁻¹⁵ M_⊙, the constraint that we obtain are much weaker: M_ej ≲ 10⁻⁸ M_⊙. However, the dynamical ejecta mass from the simulations is at least five orders of magnitude larger. As discussed previously, the secular disc wind will contribute even more material, likely increasing the overall mass of the outflows to ∼10⁻¹ M_⊙. As such, our conclusions agree with those of Bhardwaj et al. (2023): FRB 20190425A and GW190425 are not associated.

The DM contribution of the ejecta evolves with time as DM_src∝t⁻², so an FRB with a DM comparable to that of FRB 20190425A, could only have been produced with a delay of |${\sim } 18\, {\rm d}$|⁠, taking into account just the dynamical ejecta. With the inclusion of the disc wind contribution to DM, this time increases to |${\sim } 173\, {\rm d}$|⁠. Similarly, for the material to become optically thin, the dynamical ejecta (disc wind) needs to evolve for |${\sim } 12.4\, (63)\, {\rm d}$|⁠. As such, FRBs from possible stable NS merger remnants should only be expected with a delay of a few months from the merger. This picture could be modified for events observed close to on-axis, because the SGRB could clear the region at high-latitude of ejecta, thereby lowering the opacity and DM (Zhang 2014; Wang, Zhang & Zhu 2023). On the other hand, this argument would not apply to GW190425, because the inferred observation angle assuming an association to FRB 20190425A would be larger than 40° with high confidence (Bhardwaj et al. 2023).

This study is a cautionary tale that underscores our incomplete understanding of NS mergers and their dynamics. Contrary to the widespread notion that GW190425 was a prompt collapse event, we have shown that the formation of a long-lived remnant cannot be excluded on the basis of our current knowledge of the NS EOS. We can confidently exclude that GW190425 was associated with FRB 20190425A, but we cannot exclude the formation of a long-lived remnant, due to the incomplete coverage of the relevant sky regions and the faintness of the potential optical/infrared counterpart. Future observations of GW190425-like events, or deeper upper limits, can constrain models, like the BA EOS, with ‘large’ maximum NS mass. This would have a profound impact on our understanding of the properties of matter at the most extreme densities realized inside NSs. Perhaps, the most important lesson from this work is that the allocation of resources for electromagnetic follow up of NS mergers should be informed by the properties of the source, such as their chirp mass, rather than by theoretical prejudiced metrics. We urge the LIGO-Virgo-KAGRA collaboration to publicly release the chirp mass of future NS merger events promptly.

ACKNOWLEDGEMENTS

The authors thank Mattia Bulla, Antonella Palmese, Bing Zhang, and Eleonora Loffredo for useful discussions. DR 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 (NSF) under Grants No. PHY-2011725, PHY-2020275, PHY-2116686, and AST-2108467. GR acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID 279384907-SFB 1245. GR acknowledges support by the State of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006). MB acknowledges support from the Eberly Postdoctoral Research Fellowship at the Pennsylvania State University. The work of AP is partially funded by the European Union under NextGenerationEU. PRIN 2022 Prot. n. 2022KX2Z3B. KM is supported by the NSF Grant No. AST-1908689, No. AST-2108466 and No. AST-2108467, and KAKENHI No. 20H01901 and No. 20H05852. 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. Computations for this research were also performed on the Pennsylvania State University’s Institute for Computational and Data Sciences’ Roar supercomputer.

DATA AVAILABILITY

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

Footnotes

1

When our manuscript was in the final editing stages, a new analysis of ATLAS and Pan-STARRS data was announced (Smartt et al. 2023). We do not include these data, but it would not alter our conclusions.

Author notes

Alfred P. Sloan fellow