Late-time non-thermal emission from mildly relativistic tidal ejecta of compact objects merger

ABSTRACT

Mergers of compact objects [binary neutron stars or neutron star-black hole (NSBH)] with a substantial mass ratio (⁠|$q\gt 1.5$|⁠) are expected to produce a mildly relativistic ejecta within |$\sim 20^\circ$| from the equatorial plane. We present a semi-analytic approach to calculate the expected synchrotron emission observed from various viewing angles, along with the corresponding radio maps, that are produced by a collisionless shock driven by such ejecta into the interstellar medium. This method reproduces well (up to |$\sim 30~{{\ \rm per\ cent}}$| deviations) the observed emission produced by 2D numerical calculations of the full relativistic hydrodynamics. We consider a toroidal ejecta with an opening angle of |$15^\circ \le \theta _ \text{open}\le 30^\circ$| and broken power-law mass distribution, |$M(\gt \gamma \beta)\propto (\gamma \beta)^{-s}$| with |$s=s_{\rm KN}$| at |$\gamma \beta \lt \gamma _0\beta _0$| and |$s=s_{\rm ft}$| at |$\gamma \beta \gt \gamma _0\beta _0$| (where |$\gamma$| is the Lorentz factor). The parameter values are chosen to characterize merger calculation results – a ‘shallow’ mass distribution, |$1\lt s_{\rm KN}\lt 3$|⁠, for the bulk of the ejecta (at |$\gamma \beta \approx 0.2$|⁠), and a steep, |$s_{\rm ft}\gt 5$|⁠, ‘fast tail’ mass distribution. While the peak flux is dimmer by a factor of |$\sim$|2–3, and the peak time remains roughly the same (within 20 per cent), for various viewing angles compared to isotropic equivalent ejecta (⁠|$\theta _\text{open}=90^\circ$|⁠) considered in preceding papers, the radio maps are significantly different from the spherical case. The semi-analytic method can provide information on the ejecta geometry and viewing angle from future radio map observations and, consequently, constrain the ejection mechanism. For NSBH mergers with a significant mass ejection (⁠|$\sim 0.1\,{\rm M}_\odot$|⁠), this late non-thermal signal can be observed to distances of |$\lesssim 200$| Mpc for typical parameter values.

1 INTRODUCTION

Binary mergers involving compact objects with a substantial mass ratio (⁠|$q\gt 1.5$|⁠) or ones in which a prompt black hole (BH) is produced, such as neutron star-black hole (NSBH; Rosswog 2005; Just et al. 2015; Kawaguchi et al. 2015; Foucart et al. 2017; Kyutoku et al. 2018; Foucart 2020) or binary neutron star (BNS; Radice et al. 2018; Bernuzzi et al. 2020; Nedora et al. 2021b) systems, yield expanding ejecta that attains mildly relativistic velocities, |$\gamma \beta \gt 0.1$|⁠, through tidal interaction between the BH/NS and the NS (see Bernuzzi 2020; Kyutoku, Shibata & Taniguchi 2021; Chen et al. 2024, for review). Such ejecta is expected to produce non-thermal electromagnetic (EM) emission on a time-scale of |$\sim$|10 yr (Nakar & Piran 2011; Sadeh, Guttman & Waxman 2023). BNS and NSBH mergers are likely the sources of (short) gamma-ray bursts (see Mészáros 2002; Piran 2004; Nakar 2007, for reviews). Consequently, these events are expected to produce highly relativistic jets that may dominate the non-thermal emission from the tidally driven ejecta, depending on the jet’s energy, opening angle, and the angle at which it is observed. Typically, because of the jet’s lower mass and higher Lorentz factor, the late-time emission is expected to be dominated by the tidal ejecta. To interpret future observations accurately, it will be necessary to disentangle these two components.

NSBH mergers can exhibit disruptive or non-disruptive behaviors. In the former, the NS is disrupted by the tidal field of the BH, triggering mass ejection, formation of a disc around the BH, and the emission of EM radiation. A disruptive NSBH merger occurs when the disruption radius, |$R_\text{dis}$|⁠, at which the tidal forces of the BH are larger than the gravitational forces that hold the NS as a compact object, is larger than the radius of the innermost stable circular orbit, |$R_\text{ISCO}$|⁠. Conversely, in non-disruptive mergers, the NS plunges as a whole into the BH, and the phenomenon is detectable solely through gravitational wave (GW) measurements. Numerical relativity (NR) simulations of BNS mergers suggest that the mass ratio between the BNS significantly affects the ejecta geometry (Bernuzzi et al. 2020, 2024). BNS mergers of comparable mass (⁠|$q\sim 1$|⁠) eject mass dynamically due to shocks driven by the collision and also due to tidal forces. The ejecta geometry for such mergers is quasi-spherical (Radice et al. 2018; Nedora et al. 2021b). Additional spiral-wave wind component (Nedora et al. 2021b; Radice & Bernuzzi 2024) should be considered in case of long-lived neutron star remnant (absent for NSBH systems). In cases where the mass ratio is considerable, |$q\gt 1.5$|⁠, the ejected mass is concentrated within |$20^\circ$| from the equatorial plane since most of the mass is ejected by the tidal forces. In BNS mergers in which a prompt BH is produced, the dynamical ejecta is also dominated by the tidal component. Correspondingly, the angular spread of the ejecta beyond the equatorial plane is within |$\sim 20^\circ$| (Radice et al. 2018; Nedora et al. 2021b). In all cases the NR simulations suggest dynamical ejecta of |$\sim$||$10^{-3}{\!-\!}10^{-2}\,{\rm M}_\odot$|⁠, for various equation of states (EoSs) and mass ratios (see Radice, Bernuzzi & Perego 2020), with velocity profile that is well approximated by the one suggested in Sadeh et al. (2023) and confirmed in Zappa et al. (2023). NSBH mergers with comparable mass have not been observed so far and, thus, have not been considered in NR simulations. NSBH merger simulations of large mass ratio, |$q\gt 3$|⁠, provide some constraints over the dynamical ejecta of such events (Rosswog 2005; Foucart et al. 2013; Just et al. 2015; Kyutoku et al. 2015; Foucart et al. 2017; Kyutoku et al. 2018; Chen et al. 2024). In cases of tidal disruption, they imply the dynamical ejecta has only a tidal component of |$\sim$||$10^{-2}{\!-\!}10^{-1}\,{\rm M}_\odot$|⁠, for various initial parameters and EoSs. This ejecta shares a similar geometry to BNS dynamical ejecta with a high-mass ratio according to numerical relativity simulations (Bernuzzi et al. 2020, 2024); furthermore, it is consistent with the velocity profile suggested in Sadeh et al. (2023) for BNS mergers (Kyutoku et al. 2015; Chen et al. 2024). The amount of ejected mass, its geometry, and its velocity distribution depend on the parameters of the merging system (whether it is NSBH or BNS), such as the mass and the spin of the BH (in the case of NSBH), the mass of the NS, the nuclear EoS and the mass ratio, q (Rosswog 2005; Shibata & Taniguchi 2008; Lovelace et al. 2013; Kyutoku et al. 2015; Foucart et al. 2017; Hayashi et al. 2022; Chen et al. 2024). For example, for a mass ratio of |$q\gt 8$|⁠, no considerable mass ejection is expected since the NS is not tidally disrupted (Foucart et al. 2017). Measurements of the EM emission generated by the expanding ejecta will enable us to constrain its properties and, hence, the above-mentioned parameters. The opening angle of the ejecta is approximately conserved (Kyutoku et al. 2015) because the velocity direction does not change appreciably once the hydrodynamic merger interaction becomes negligible. The azimuthal opening angle of the ejecta varies from |$180^\circ$| to |$360^\circ$| for different initial conditions (Kyutoku et al. 2015; Kawaguchi et al. 2016; Bernuzzi et al. 2020). Margalit & Piran (2015) consider the effect of a non-spherical (and non-relativistic) ejecta structure on the synchrotron light curve. They provide a rough approximation of the peak time scale that is delayed with respect to the isotropic equivalent case and reach the conclusion that the peak flux scale remains roughly the same as in the isotropic equivalent case. However, they do not provide full numerical calculations and base their conclusions on an analytic approximation from Nakar & Piran (2011). It should be noted here that more material is expected to be ejected from the remnant disc (Fernández et al. 2017). Due to the lower velocities and time delay compared to the dynamical ejecta, no mixing between the two components is expected (Fernández et al. 2015). Furthermore, since the non-thermal signal is highly sensitive to the shocked plasma Lorentz factor (Sadeh et al. 2023), the lower velocity of the disc outflow (Fernández et al. 2017, 2019) would make its potential non-thermal signal much less bright.

The rate of BNS mergers was estimated by the Galactic BNS systems to be |$\sim 10\left(\frac{D}{200\text{Mpc}}\right)^3\text{yr}^{-1}$| (Phinney 1991; Pol, McLaughlin & Lorimer 2019, 2020). Such estimation is not possible for the rate of NSBH mergers since our Galaxy has no known NSBH binary systems. Due to the low number of detected BNS and NSBH mergers via GW, the merger rate obtained from these observations is extremely uncertain. For BNS mergers it is between 0.1 to |$14\left(\frac{D}{200\text{Mpc}}\right)^3\text{yr}^{-1}$| and for NSBH mergers it is between 0.1 and |$1.5\left(\frac{D}{200\text{Mpc}}\right)^3\text{yr}^{-1}$| (LIGO Scientific Collaboration and KAGRA Collaboration et al. 2023; The LIGO Scientific Collaboration, the Virgo Collaboration & the KAGRA Collaboration 2024). The fraction of EM bright events out of the GW detectable NSBH mergers events is currently unknown due to the low statistics- only three confirmed events so far in the third and fourth LIGO-Virgo-KAGRA run (Abbott et al. 2021; The LIGO Scientific Collaboration et al. 2024). However, first attempts have been made (Fragione 2021; Biscoveanu, Landry & Vitale 2023; The LIGO Scientific Collaboration et al. 2024), suggesting a fraction of |$\sim 10{\!-\!}30~{{\ \rm per\ cent}}$|⁠. Martineau et al. (2024) estimated the ejected mass in the observed event, GW230529, to be negligible, |$\sim 10^{-5}{\!-\!}10^{-3}\,{\rm M}_\odot$|⁠, due to the low BH spin |$\chi _\text{BH}\lesssim 0.1$|⁠. While the mass ratio in NSBH mergers is likely to be above |$q=1.5$| considered the observed mergers so far, for the observed BNS mergers GW170817 and GW190425, the mass ratio estimate varies between |$q\approx 1$| to |$q\approx 2.5$| for different spin priors.

This paper explores synchrotron emission spanning the radio-to-X-ray range, originating from a collisionless shock driven by mildly relativistic toroidal ejecta with an opening angle of |$15^\circ$| to |$30^\circ$| and a broken power-law dependence of mass on momentum, with parameter values characteristic of the results of numerical calculations of the BNS ejecta (see Section 2). The ejecta expands into a uniform interstellar medium (ISM) characterized by a number density n. We derive a semi-analytic calculation method for the case of axisymmetric distribution of mass (as in Gompertz et al. 2023) exhibiting a power-law momentum-mass dependence (Sadeh et al. 2023), where the parameter values are consistent with the results of numerical calculations of the tidal ejecta (Kyutoku et al. 2015; Foucart et al. 2017; Kyutoku et al. 2018; Bernuzzi et al. 2020; Nedora et al. 2021b; Bernuzzi et al. 2024; Chen et al. 2024). This semi-analytic method is valid for as long as the reverse shock did not cross the fast tail part of the ejecta, during this phase the shocked plasma is not expected to go through a significant lateral expansion, so we assume angular independent distribution of shocked plasma parameters. We derive the non-thermal emission spectra under the assumption of a power-law energy distribution for electrons, given by |${\rm d}n_e/{\rm d}\gamma _e\propto \gamma _e^{-p}$|⁠. Here, |$\gamma _e$| represents the electron Lorentz factor within the plasma’s rest frame, where p is within the range of |$2\le p\le 2.5$|⁠. Concurrently, the fractions |$\varepsilon _e$| and |$\varepsilon _B$| correspond to the post-shock internal energy density proportions attributed to non-thermal electrons and magnetic fields, respectively. This phenomenological description, capturing the post-shock plasma conditions, finds support across a diverse spectrum of observations and plasma calculations encompassing both relativistic (see Keshet & Waxman 2005; Keshet 2006; Waxman 2006; Bykov & Treumann 2011; Sironi, Spitkovsky & Arons 2013; Kobzar et al. 2021) and non-relativistic shocks (see Blandford & Eichler 1987; Pohl, Hoshino & Niemiec 2020; Ligorini et al. 2021). We verified the accuracy of the semi-analytic results for the non-thermal emission presented in this paper through a comprehensive comparison with full 2D numerical computations across a wide range of parameter configurations.

The paper is organized as follows. In Section 2, we describe the framework of the calculations presented later on. We derive the semi-analytic calculation for the synchrotron light curve and corresponding radio map in Section 3. Then, we describe the 2D relativistic-hydrodynamics numerical calculation and the post-processing we use to produce light curves in Section 4. The accuracy of the semi-analytic method is estimated by comparison with the detailed numeric calculations in Section 5, while in Section 6, we provide a few examples for light curves and radio maps produced by the semi-analytic approach. Finally, our conclusions are summarized in Section 7.

2 CALCULATION FRAMEWORK

2.1 Ejecta geometry description

In Sadeh et al. (2023) and Sadeh, Linder & Waxman (2024), we considered a spherical ejecta with the following mass profile:

with parameter values characteristic of the results of numerical calculations of the BNS ejecta; |$0.3\lt \beta _0\lt 0.5$|⁠, |$5\lt s_\text{ft}\lt 12$|⁠, |$0.5\lt s_\text{KN}\lt 3$|⁠, and |$10^{-6}\,{\rm M}_\odot \lt {\rm M}_0\lt 10^{-4}\,{\rm M}_\odot$|⁠. This analytic form provides a good approximation for the variety of ejecta profiles obtained in NR simulations of BNS mergers (Zappa et al. 2023). A merger of compact objects with a mass ratio of |$q\gt 1.5$|⁠, or one in which a prompt BH is produced, is expected to eject mass mainly within an angle of |$\sim 20^\circ$| from the equatorial plane and follow the same velocity profile as the dynamical ejecta of mergers with comparable masses (Kyutoku et al. 2015; Foucart et al. 2017; Radice et al. 2018; Bernuzzi et al. 2020; Nedora et al. 2021b; Bernuzzi et al. 2024). The ejecta azimuthal structure varies from full azimuthal symmetry to a crescent shape with an azimuthal opening of |$\phi \approx 180^\circ$| (Kyutoku et al. 2015). The geometric structure we consider to describe such ejecta is a toroid (sphere truncated by a double-sided cone), which limits us to the case of full azimuthal symmetry, with a mass profile that is given by equation (1), and various opening angles, |$15^\circ {\!-\!}30^\circ$|⁠, that are consistent with the results of NR simulations, at some initial lab time |$t_0$| (rest frame of the ISM). The angle between the equatorial plane and the ejecta ‘edge’ is defined as |$\theta _\text{open}$|⁠, while |$\theta _\text{v}$| is the angle between the symmetry axis and the line of sight. In Fig. 1, we provide a sketch illustrating the ejecta geometry. As the ejecta expands, a forward shock is driven into the ISM, and a reverse shock is driven into the ejecta.

2.2 Coordinate systems

The natural coordinate system for the ejecta initial conditions is a spherical coordinate system |$(R,\phi ,\theta)$|⁠, where R is the distance from the origin, |$\phi$| is the azimuthal angle, and |$\theta$| is the polar angle, where the pole is aligned with the ejecta symmetry axis (see Fig. 1). To calculate the synchrotron emission, it is convenient to define a coordinate system that is oriented with the line of sight. We define another spherical coordinate system, |$(R,\varphi ,\xi)$|⁠, where |$\varphi$| is the azimuthal angle (chosen such that |$\varphi =0$| corresponds to the equatorial plane in case of |$\theta _\text{v}=\pi /2$|⁠), |$\xi$| is the polar angle, and the pole is aligned with the line of sight. The coordinate transformation is

similar to Govreen-Segal & Nakar (2023), although one difference is the choice of |$\phi$| such that for |$\phi =\pm \pi /2$| we have |$\xi =\theta _\text{v}\pm \theta$|⁠.

A cylindrical coordinate system, |$(r,\varphi ,z)$|⁠, such that the z-axis is aligned with the line of sight, is useful for the semi-analytic calculation presented below. For such a coordinate system, the coordinate transformation is

3 SEMI-ANALYTIC CALCULATION

A semi-analytic approach for calculating the expected synchrotron emission from a toroidal ejecta with a velocity profile given by equation (1) is presented below. This approach is valid as long as the reverse shock propagates through the ejecta steep, fast tail. After the reverse shock crosses the fast tail, a full 2D calculation is needed for a reliable calculation of the synchrotron emission. The semi-analytic analysis code offers a computation method for the approximated dynamics, anticipated flux, and the radio image on the sky.

3.1 Dynamics

The toroidal ejecta has mass, |$M=M_{\text{iso}}\sin {\theta _\text{open}}$|⁠, where |$M_{\text{iso}}$| is the isotropic equivalent mass and |$\theta _\text{open}$| is the ejecta opening angle. The ejecta propagates through uniform ISM density, forming a forward-reverse shock structure, as shown in Fig. 1. As long as the reverse shock did not cross the steep fast tail, the total internal energy (and thus, the emissivity) of the system is dominated by the freshly shocked material within the initial opening angle because of the following argument: Given the steep ejecta profile with a power-law index of |$s_\text{ft} \ge 5$|⁠, the energy and momentum injected into the shocked medium increase rapidly with decreasing velocity. This is because the mass of the ejecta, which scales as |$M(\gt \gamma \beta) \propto (\gamma \beta)^{-s_\text{ft}}$|⁠, grows more quickly as velocity decreases. Consequently, the freshly shocked material, which is continuously energized by this injection, dominates the internal energy of the system. As the shock front propagates radially through the interstellar medium (ISM), the swept-up ISM mass is primarily composed of this freshly shocked material. The unshocked ejecta, confined within an opening angle |$\theta _\text{open}$|⁠, also moves radially, injecting additional energy and momentum into the shocked plasma. Due to conservation of momentum, the momentum of the freshly shocked material is dictated by this added momentum from the ejecta at the position of the reverse shock. We conclude the shocked plasma lateral expansion should have a limited effect on the observed emission before the reverse shock crosses the steep, fast tail. For larger values of |$s_\text{ft}$|⁠, this effect becomes more pronounced, resulting in a less substantial lateral expansion. Thus, we approximate the ejecta deceleration dynamics using the dynamics of stratified spherical ejecta with a corresponding mass of |$M_{\text{iso}}=M/(\sin {\theta _\text{open}})$|⁠.

3.2 Semi-analytic spherical dynamics

In Sadeh et al. (2023), we developed a semi-analytical calculation method for modelling the dynamics of spherical mildly relativistic ejecta with a velocity profile given by equation (1), propagating into a cold and uniform ISM particularly applicable when the reverse shock propagates through the ejecta fast tail. Within the framework of this calculation, we approximated both shocked ISM and shocked ejecta as two uniform separate layers between the forward-reverse shock structure. This method is consistent with the full 1D relativistic hydrodynamics simulations (Sadeh et al. 2023). The detailed derivation of this approach is presented in Appendix  A.

3.3 Intensity and flux

Using the approximated dynamics mentioned above, we obtain the flow profiles behind the shocks front: the internal energy density |$e(R,t)$|⁠, the proper mass density |$\rho ^{\prime }(R,t)$|⁠, and the Lorentz factor |$\gamma (R,t)$|⁠, all as a function of the shock radius and time in the rest frame of the ISM. We approximate the shocked plasma as two separate uniform layers behind the shocks with different flow profiles for the two different layers. The thickness of these layers is determined by conservation of mass both for forward and reverse shocks. We employ a phenomenological approach to describe synchrotron emissivity: the fraction of energy carried by electrons is denoted by |$\varepsilon _e$|⁠, and the fraction of energy carried by magnetic fields is denoted by |$\varepsilon _B$|⁠. We assume a power-law electron distribution within the shocked layers, |${\rm d}n_e/{\rm d}\gamma _e\propto \gamma _e^{-p}$|⁠. Following the assumption that the emission is dominated by shocked plasma within the initial opening angle, we restrict the emissivity to within |$\theta _\text{open}$|⁠. Mathematically, this is expressed as

where |$j^{\prime }_\nu \left(e,\rho ,\gamma \right)$| is the emissivity as a function of the flow profiles, and the prime (⁠|$^{\prime }$|⁠) notation indicates proper frame quantities. Subsequently, the emissivity and intensity are computed by the methodology outlined in Sadeh et al. (2024): |$j_\nu$| is defined by following Rybicki & Lightman (1979), and |$I_\nu$| is calculated by integrating over the |$j_\nu$| while taking into account arrival time effects and relativistic effects. We only consider frequencies above the self-absorption frequency and below the cooling frequency, |$\nu _a\lt \nu \lt \nu _c$|⁠, since the bulk of the radio to X-ray observations are expected to be in this part of the spectrum (Sadeh et al. 2024). To calculate the flux, we consider the fact that the azimuthal observer symmetry breaks in the case of |$\theta _\text{v}\ne 0$|⁠, so |$I_\nu$| varies for different |$\varphi$| values. We integrate over the contribution from each azimuthal angle by

where R is the shock radius. For the light-curve temporal cut-off, we use the peak time estimation from Sadeh et al. (2023),

where |$M_0$| here is the isotropic equivalent mass of the fast tail (⁠|$M_{\text{iso}}=M/(\sin {\theta _\text{open}})$|⁠), |$n=10^{-2}n_{-2}\,{\rm cm}^{-3}$|⁠, and |$M_0=10^{-4}M_{0,-4}\,{\rm M}_\odot$|⁠. For |$0^\circ \le \theta _\text{v}\lesssim 20^\circ$|⁠, a delay is expected between the peak time in the isotropic equivalent case given by equation (6) and the toroidal ejecta peak time due to lateral expansion of the shocked plasma after the reverse shock crosses the fast tail. This is due to relativistic beaming as |$\sin ^{-1}(1/\gamma _0)\sim 70^\circ$| for typical values of |$\gamma _0\approx 1.05$|⁠. The peak flux in the isotropic equivalent case is proportional to |$M_0$| (Sadeh et al. 2023). Since |$\theta _\text{open}\sim 15^\circ {\!-\!}30^\circ$|⁠, we expect the peak flux in the toroidal case to be lower by a factor of |$\sin {15^\circ }-\sin {30^\circ }\sim 2-3$| compared to the isotropic equivalent case. Thus, the peak flux is approximately (Sadeh et al. 2023)

with |$D=10^{26.5}D_{26.5}$| cm, |$\nu =10^{9.5}\nu _{9.5}$| Hz, |$\varepsilon _{e}=10^{-1}\varepsilon _{e,-1}$||$\varepsilon _{B}=10^{-2}\varepsilon _{B,-2}$|⁠,

and |$F^\text{iso}_{\nu ,\rm peak}$| is the peak flux in the equivalent isotropic case. |$f_{\rm ft}(p)$| deviates from unity by less than 30  per cent for |$2\lt p\lt 2.5$| (Sadeh et al. 2023).

4 NUMERICAL CALCULATION

4.1 Numerical set-up

We employed numerical solutions of the 2D special relativistic hydrodynamics equations to validate our model through the RELDAFNA code (Klein 2023). RELDAFNA is an Eulerian code employing the Godunov scheme, incorporating adaptive mesh refinement (AMR) and second-order accuracy in time and space integration. Its efficient parallelization allows for high-resolution calculations, even for complex multiscale problems. RELDAFNA’s accuracy was confirmed by comparing it with similar codes (Zhang & MacFadyen 2006; Meliani et al. 2007) using standard test problems. Utilizing RELDAFNA, we solved the equations governing a relativistic outflow of an ideal fluid with a smoothly varying adiabatic index (between |$5/3$| and |$4/3$|⁠) depending on the internal energy density

which is in very good agreement with the EoS of mildly relativistic fluid provided by Synge (1957). The 1.3 factor is added for better agreement with the exact solution in the range of |$0\lt e/\rho ^{\prime } c^2\lt 2$| (corresponding to |$1\lt \gamma \lt 3$|⁠); see Appendix  B. Our simulations were performed in a 2D cylindrical coordinate system |$(r,z)$| assuming axial symmetry, and we adopted a Courant–Friedrichs–Lewy (CFL) number of 0.5, which provided efficient convergence (physical justified as long as it is |$\le 1$|⁠). The initial grid prioritizes resolving the ejecta by placing a higher concentration of cells within it, with a coarser spacing outside. Typically, the simulation begins with |$500\times 500$| cells, focusing most cells within the ejecta. AMR dynamically adjusts the resolution throughout the simulation, optimizing cell distribution along both the r and z axes to effectively capture regions of significant variation in pressure, density, or Lorentz factor. The initial pressure in all of the cells is set to |$P=10^{-10}\rho _{\text{ISM}} c^2$| (where |$\rho _{\text{ISM}}$| is the ISM mass density), and the computational domain boundaries are set at |$r=10^{18}$|cm and |$z=10^{18}$|cm. We systematically refined the initial spatial grid and temporal steps to assess convergence. We also iteratively increased the number of times the regridding scheme is allowed to multiply the number of cells within a time-step. This process continues until no significant changes are observed in the hydrodynamics and emergent flux. Reflective boundary conditions were implemented at both the symmetry axis (⁠|$r=0$|⁠) and the symmetry plane (⁠|$z=0$|⁠) to account for symmetry considerations. Our numerical calculations adopt initial conditions comprising toroid ejecta, as illustrated in Fig. 1. The ejecta is characterized by a flow profile described in equation (1). The initial radius (at |$t=t_0$|⁠) of a mass shell |$M(\gt \gamma \beta)$| is determined by |$r_0=\beta ct_0$|⁠. The expanding ejecta is initially embedded at |$r\gt c t_0$| in a static, uniform, cold gas with negligible pressure. See Fig. 2 for an example of a density map snapshot calculated with RELDAFNA.

4.2 Numerical light curve

From the 2D grid obtained via RELDAFNA, we construct a 3D grid by subdividing cells into segments that represent different azimuthal angles such that the number of segments is appropriate for achieving flux convergence. By post-processing the hydrodynamic profiles, we calculate specific emissivity, |$j^{\prime }_\nu (\vec{r},t)$| (following Rybicki & Lightman 1979; Sadeh et al. 2024, Appendix B), for frequency below the cooling frequency and above the synchrotron and self-absorption frequencies, |$\nu _a,\nu _m\lt \nu \lt \nu _c$|⁠. The emissivity is calculated in the plasma proper frame for each cell within the 3D grid at a discrete set of lab times. Each fluid element is assumed to exist for a period of |$\Delta t$| interval between time-steps. This calculation adopts a phenomenological approach for the magnetic field and electron energy distribution, assuming fractions |$\varepsilon _B$| and |$\varepsilon _e$| of internal energy density for magnetic fields and electrons, respectively. Furthermore, we assume a power-law energy distribution for electrons, |$dn_e/d\gamma _e\propto \gamma _e^{-p}$|⁠. The flux contribution from each cell for every observed time is computed as follows,

where R is the distance from the origin, |$\Delta V$| is the cell’s volume, D is the luminosity distance, t is the time in the lab frame, |$t_\text{obs}$| is the observer time, |$\xi$| is the angle between the cell and the line of sight and |$\xi _\beta$| is the angle between the cell velocity direction and the line of sight. |$\xi$| pose different values for different observation angles following equation (2). Similarly, |$\xi _\beta$| shares the same property by considering the following transformation, |$\xi \rightarrow \xi _\beta ,\theta \rightarrow \theta _\beta$| (where |$\theta _\beta$| is the angle between the velocity direction and the symmetry axis), to equation (2). The total observed flux per unit frequency and observed time |$t_\text{obs}$| is obtained by summing the different contributions |$\Delta F_\nu$| from all of the fluid elements at all the discretized times.

5 VALIDATION OF THE SEMI-ANALYTIC METHOD

In Section 5.1, we verify our estimation regarding the dominant contribution of the freshly shocked material within the initial opening angle of the ejecta to the internal energy. This verification holds true as long as the reverse shock has not yet crossed the ejecta fast tail. In Section 5.2, we compare the light curves produced by the full 2D numerical calculation with our semi-analytic scheme to verify its validity. Finally, in Section 5.3, we compare the flux angular distribution obtained by our semi-analytic method to the full 1D numerical calculation of a spherical ejecta case to test its validity.

5.1 Hydrodynamics

In our semi-analytic method, we assume that the internal energy is dominated by material within the initial opening angle of the ejecta. We ran numeric calculations with various initial conditions to test our assumption. The opening angle as a function of lab time is defined as the angle from the equatorial plane within which 90  per cent of the internal energy, excluding rest mass energy, is included. A few examples are given in Fig. 3. We consider various parameter options since all of the following parameters |$\lbrace s_\text{ft},n,M_R,\rbrace$| affect the ejecta dynamics, where |$M_R\equiv M(\gamma \beta \gt 1)= M_0(\gamma _0\beta _0)^{s_{\rm ft}}$|⁠, is the ‘relativistic mass’ with momentum |$\gamma \beta \gt 1$|⁠. We consider only typical values for the opening angle, |$\theta _\text{open}=20^\circ$|⁠, and the velocity cut-off, |$\beta _0=0.3-0.4$|⁠, that are consistent with the ejected material in merger simulations. We indeed find that the contribution to the internal energy is dominated by the material within the initial opening angle (up to |$\sim 30~{{\ \rm per\ cent}}$| corrections) as long as the reverse shock does not cross the ejecta fast tail. An evident trend in Fig. 3 is that larger values of |$s_\text{ft}$| lead to a less pronounced lateral expansion. This is because, at larger |$s_\text{ft}$|⁠, the freshly shocked material makes up a larger fraction of the shocked mass, as explained in Section 3.1.

5.2 Light curves

We consider only frequencies that are between the self-absorption frequency/peak frequency to the cooling frequency, |$\nu _a,\nu _m\lt \nu \lt \nu _c$|⁠, since the radio observations are expected to be in this part. In Fig. 4, we compare the light curve obtained semi-analytically (for as long as the reverse shock does not cross the fast tail, Section 3) with the full 2D numerical calculation (Section 4). The agreement is within 10’s of per cent for a range of relevant values of |$\lbrace n,M_R,\beta _0,s_{\rm ft},\theta _\text{v},\theta _\text{open},p\rbrace$|⁠, along with |$\varepsilon _e=10^{-1}$| and |$\varepsilon _B=10^{-2}$|⁠, which modify the light curve based on a known analytical dependence, and |$s_\text{KN}=1.5$|⁠, which alters the light curve after the peak. As expected, a |$\sim 20~{{\ \rm per\ cent}}$| delay is observed for |$\theta _\text{v}=0^\circ$|⁠.

5.3 Sky image

To further validate our semi-analytic calculation, in Fig. 5, we show an example of the flux emitted from different annuli obtained numerically and semi-analytically for spherical ejecta (we considered |$\theta _\text{open}=90^\circ$|⁠), along with the analytic estimate of the image radius from (Sadeh et al. 2024). Due to relativistic beaming and time arrival effects, the image is a relatively narrow ring. This ring radius increases with time since the blast wave expands and the relativistic beaming effect decreases. The emission from large viewing angles, which was initially negligible due to the relativistic beaming, becomes brighter as the velocity declines.

6 RESULTS

The mass ejection in NSBH mergers simulations typically reaches |$\sim 0.1\,{\rm M}_\odot$| in disruptive merger systems (Foucart et al. 2017; Kyutoku et al. 2021; Chen et al. 2024), implying an isotropic equivalent mass of |$M\sim 0.3\,{\rm M}_\odot$| (for a typical opening angle of |$\theta _\text{open}=20^\circ$|⁠). This is consistent with values of |$M_R=10^{-4}\,{\rm M}_\odot ,s_\text{ft}=7$|⁠. In Figs 6 –7, we provide several examples of the obtained radio maps and the corresponding light curves from the semi-analytic calculation for several values of |$\theta _\text{v}$| and typical parameters expected for NSBH merger, |$\lbrace s_\text{ft},s_\text{KN},\beta _0,n,M_R,\varepsilon _e,\varepsilon _B,\theta _\text{open}\rbrace =\lbrace 7,1.5,0.3,10^{-2}\,\text{cm}^{-3},10^{-4}\,{\rm M}_\odot ,10^{-1},10^{-2},20^\circ \rbrace$|⁠. Notice that at the peak, the shocked plasma velocity is |$\sim \beta _0=0.3$|⁠, and the radiation emitted from such plasma is beamed into an angle of |$\sim \sin ^{-1}(\frac{1}{\gamma _0})\approx 70^\circ$|⁠. This is the reason for the small difference at the peak flux between viewing angles of |$0^\circ {\!-\!}30^\circ$|⁠, in which a larger fraction of the radiation is beamed away, to viewing angles of |$60^\circ{\!-\!}90^\circ$|⁠. We add our analytic estimation for the light curve and image radius for a spherical ejecta case with the same equivalent isotropic mass for comparison.

We find that the toroidal ejecta case is less bright by a factor of |$\sim 2{\!-\!}3$| than the isotropic equivalent case (Fig. 7), and expected to peak on a time-scale of 10 yr. Furthermore, the peak flux’s dependence on the viewing angle is minimal for this toroidal ejecta. Oppositely, the radio map is significantly different between the different viewing angles (Fig. 6), providing an opportunity to constrain the viewing angle in future radio map observations.

7 CONCLUSIONS

In this work, we study the non-thermal emission of a toroidal (see Fig. 1), mildly relativistic ejecta with broken power-law mass distribution, given by equation (1). We provided a semi-analytic calculation method for the synchrotron emission, assuming emission only from the initial opening angle of the ejecta (Section 3) at times in which the reverse shock propagates through the steep, fast tail. Furthermore, we provided analytic estimations for the peak time and flux in equations (6) and (7). Then, by calculating the full 2D hydrodynamics (Section 4), we validated our semi-analytic approach in Section 5 for a wide range of ejecta parameter values characteristic of merger calculation results; an equatorial ejecta opening angle of |$15^\circ \lt \theta \lt 30^\circ$|⁠, a steep (⁠|$s_{\rm ft}\gt 5$|⁠, |$\gamma _0\beta _0\ge 0.3$|⁠) mass distribution for mildly relativistic velocities and a moderate mass distribution for the bulk of the ejecta (⁠|$1\lt s_{\rm KN}\lt 3$|⁠, |$0.1\lt \gamma _0\beta _0\le 0.3$|⁠). We also used various values of electrons power-law index, |$2.1\lt p\lt 2.4$|⁠, of the electron energy distribution. This is the first semi-analytic method supported by full 2D relativistic numeric calculations for the expected late non-thermal emission from mergers of compact objects with a substantial (⁠|$q\gt 1.5$|⁠) mass ratio. We found that the full 2D numerical calculation is consistent with our assumption (Fig. 3), i.e. that during reverse shock propagation through the fast tail, the freshly shocked material within the initial opening angle of the ejecta has a dominant contribution to the total internal energy. Additionally, we found that, in agreement with equation (7), the peak flux for toroidal ejecta is smaller by a factor of |$\sim 2{\!-\!}3$| (varies for different viewing angles) than spherical ejecta with the same isotropic equivalent mass, while the peak time estimation from (equation 6) for spherical ejecta (with the same isotropic equivalent mass) is within 20 per cent deviation from the peak time results of the full numerical calculation for toroidal ejecta presented in this paper (Fig. 4). The radio map is significantly different between the toroidal and the spherical case in which the image is a symmetric ring (Fig. 6). For |$\theta _\text{v}\ge 20^\circ$| the width of the observed image (for |$\theta _\text{v}=90^\circ$| the width is parallel to the equatorial plane) is comparable with the ring diameter in the spherical case (Sadeh et al. 2024). Concurrently, the length (perpendicular to the width) is considerably smaller, providing a clear indication of equatorial ejecta. The different viewing angles affect the observed image orientation (see Fig. 6), providing the possibility to constrain the viewing angle in future radio map observations. For |$0^\circ \le \theta _\text{v}\le 20^\circ$|⁠, the image is observed as a narrow ring (Fig. 6), similar to the spherical case. The semi-analytic scheme we presented in this work can be used to fit all of the model parameters |$\lbrace s_\text{ft},\gamma _0,n,M_R,\varepsilon _e,\varepsilon _B,p,\theta _\text{open},\theta _\text{v}\rbrace$| to future observations of light curves and radio maps by an iterative process that scans the possible parameter space, e.g. Markov chain Monte Carlo. Degeneracies between the different parameters are expected to arise during the fitting process, for example, see equations (6)–(7). There are two distinctions: the radio map shape is expected to provide stringent bounds over |$\theta _\text{open}$| and |$\theta _\text{v}$|⁠, while the synchrotron spectrum provides p.

Different mass ratios between the compact objects correspond to different ejecta geometries (Bernuzzi et al. 2020, 2024). Consequently, bounds over the ejecta geometry can constrain this mass ratio. For example, in GW170817, the estimation for the mass ratio varies between |$1\lt q\lt 2.6$| (for low spin assumption, it is |$1\lt q\lt 1.4$|⁠), while the estimation for the viewing angle from gravitational waves is |$32^\circ \pm 10^\circ$|⁠. Observing the radio map of the late non-thermal emission, in this case, can potentially resolve the large uncertainty in the mass ratio.

The mass ejection in NSBH mergers simulations typically reaches |$\sim 0.1\,{\rm M}_\odot$| in disruptive merger systems (Foucart et al. 2017; Kyutoku et al. 2021; Chen et al. 2024), implying an isotropic equivalent mass of |$M\sim 0.3\,{\rm M}_\odot$| (for |$\theta _\text{open}=20^\circ$|⁠). Following equation (7), we expect these events to be brighter than the late non-thermal emission expected to follow BNS mergers, which typically have an isotropic equivalent mass of |$M\sim 0.01\,{\rm M}_\odot$| (Radice et al. 2018; Nedora et al. 2021a), by a factor of |$\sim 5{\!-\!}10$|⁠. Thus, although the toroidal ejecta case is less bright by a factor of |$\sim 2-3$| than the isotropic equivalent case (see Fig. 7), the late non-thermal emission of NSBH merger is expected to be observed in radio at distances of |$\lesssim 200$| Mpc by the Very Long Array (VLA) and expected to peak on a time scale of ten years for typical parameter values (⁠|$M_0=10^{-1}\,{\rm M}_\odot$|⁠, |$\beta _0=0.3$|⁠, |$n=10^{-2}$| cm|$^{-3}$|⁠). In conclusion, the peak flux’s dependence on the viewing angle is minimal for this highly non-spherical geometry. Instead, among the ejecta parameters, the peak flux primarily relies on the total mass in the fast tail, as illustrated in Fig. 7. This suggests that our formula for the peak flux (Sadeh et al. 2023) remains a valid estimate, within a factor of a |$\sim$|few, for the non-thermal peak flux from various 3D ejecta structures produced in full NR calculations of compact object mergers, provided that |$M_0$| is replaced with the total mass in the fast tail.

ACKNOWLEDGEMENTS

We thank Y. Klein and N. Wygoda for their authorization to use the RELDAFNA code and for their contribution. We also thank Noya Linder for her helpful contribution. Finally, we thank Eli Waxman for helpful discussions and insightful comments.

DATA AVAILABILITY

The data underlying this article will be shared following a reasonable request to the corresponding author.

REFERENCES

APPENDIX A: SEMI-ANALYTIC CALCULATION DESCRIPTION

As the ejected material expands, it initiates a forward shock that propagates into the ISM, while simultaneously, a reverse shock is generated within the ejected material itself (see Fig. A1). The deceleration imposed by the reverse shock reduces kinetic energy in the ejecta, which is subsequently transmitted to the heated ionized plasma within the shocked ISM.

We consider an ejecta with a mass profile given by |$M(\gt \gamma _4\beta _4)=M_R(\gamma _4\beta _4)^{-s}$| (numeric subscripts relate to the region presented in Fig. A1) and assume that the ISM density is uniform. The initial energy of the ejecta (excluding the rest-mass energy) is given by

As mentioned in Section 3 we approximate the shocked layers behind the shocks with uniform flow profiles (⁠|$e,\rho$| and |$\gamma$|⁠), given by the shock jump conditions (Blandford & McKee 1976),

where |$\rho ^{\prime }$| and |$\rho$| are the mass density of the shocked material and the unshocked medium correspondingly, both in their own rest frame. |$\hat{\gamma }$| is the adiabatic index relating the internal energy density, e, and the pressure, p (⁠|$p=(\hat{\gamma }-1)e$|⁠). We use an adiabatic index |$\hat{\gamma }(e/\rho ^{\prime })$| varying from |$4/3$| to |$5/3$| following the analysis of Synge (1957) for a plasma of protons and electrons. This is inaccurate for the shocked ejecta plasma, which is composed of a wide range of nuclei. However, the dependence of the results on the exact form of |$\hat{\gamma }(e/\rho ^{\prime })$| for the shocked ejecta plasma is weak (as verified by the numeric calculations). We can relate the shocked ISM mass, |$M_\text{ISM}$|⁠, and the shocked ejecta mass, |$M(\gt \gamma _4\beta _4)$|⁠, to the mass densities in the rest frames by

where R is the forward shock radius, |$R_\text{cd}$| is the contact discontinuity radius, and |$\Delta _\text{ej}$| the thickness of the shocked ejecta layer. To find the shocked ISM thickness, |$\Delta$|⁠, we consider mass conservation

The thickness of the shocked ejecta layer is given by

where |$\gamma _3^{*}=\gamma _2\gamma _4\left(1-\sqrt{1-\frac{1}{\gamma _2^2}-\frac{1}{\gamma _4^2}+\frac{1}{\gamma _2^2\gamma _4^2}}\right)$| is the Lorentz factor of region 3 in the frame of the unshocked ejecta. Since |$R-\Delta =R_\text{cd}$| we can write

The total energy conservation equation, including pressure, is derived from the conservation of the energy-momentum tensor, |$\partial _\mu T^{\mu 0}=0$|⁠, such that the conserved quantity is

This can be expressed as

where |$E_{\text{k},i}/E_{\text{th},i}$| is the kinetic/thermal energy in region i. They are given by

equation (A8) enables us to derive an expression for |$M_\text{ISM}$|

Finally, the equal pressure between regions 2 and 3, which are separated by the contact discontinuity, leads to

We numerically solve equation (A11) by using the expressions in equations (A3), (A6), and (A10) to find |$\gamma _2(\gamma _4)$| for a given |$\gamma _4$|⁠. Then we find the shock radius, |$R(\gamma _4)$|⁠, and the flowfields, |$\rho _3(\gamma _4),\rho _2(\gamma _4),e_3(\gamma _4),e_2(\gamma _4)$| by considering equations (A2), (A3), and (A10). The forward shock Lorentz factor, |$\Gamma$|⁠, is given by (Blandford & McKee 1976)

Thus, the radius as a function of time, |$R(t)$|⁠, can be found by numerically solving the following equation:

Using the above results, we finally define the shock radius, emitting layers thickness, and flow fields as a function of shock radius and/or lab time. All calculations above are done in terms of |$R_R,t_R$|⁠, and |$M_R$| (Sadeh et al. 2023) for better accuracy and efficiency. We multiply |$R(t)$| by 0.95 and |$\gamma _2\beta _2$| by 0.9 for better agreement with full numerical calculations.

APPENDIX B: EQUATION OF STATE APPROXIMATION

The theory of relativistic perfect gases gives us an expression for the specific enthalpy as a function of |$\Theta \equiv \frac{p}{\rho ^{\prime } c^2}$|⁠, that holds for a gas that is composed of the same particles and in the limit of a small free path when compared to the sound wavelength (Mignone & McKinney 2007). It has the following form (Synge 1957)

where |$K_3$| and |$K_2$| are, respectively, the modified Bessel functions of the second kind of order 3 and 2. The specific enthalpy is defined by

plugging the relation |$p=(\hat{\gamma }-1)e$| we obtain

and

These equations are solved numerically to find the adiabatic index of the fluid as a function of |$\frac{e}{\rho ^{\prime } c^2}$|⁠. In Fig. B1, we compare this result with the following approximation:

for different values of |$Y=1/1.1/1.3$|⁠. We also compare the result to the approximation suggested in Ayache, van Eerten & Eardley (2022),

We find that our approximation fits well with the numerical results with the following Y values,