# Non-thermal emission from mildly relativistic tidal ejecta of compact
objects merger

Gilad Sadeh
Dept. of Particle Phys. & Astrophys., Weizmann Institute of Science, Rehovot
76100, Israel  E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Mergers of compact objects (binary neutron stars, BNS, or neutron star-black
hole, NSBH) with a substantial mass ratio ($q>1.5$) are expected to produce a
mildly relativistic ejecta within $\sim 20^{\circ}$ from the orbital 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 is supported by 2D numerical
calculations of the full relativistic hydrodynamics with an agreement in the
observed emission of $\sim 30\%$. We consider a toroid ejecta with an opening
angle of $15^{\circ}\leq\theta_{\text{open}}\leq 30^{\circ}$ and broken power-
law mass distribution, $M(>\gamma\beta)\propto(\gamma\beta)^{-s}$ with
$s=s_{\rm KN}$ at $\gamma\beta<\gamma_{0}\beta_{0}$ and $s=s_{\rm ft}$ at
$\gamma\beta>\gamma_{0}\beta_{0}$ (where $\gamma$ is the Lorentz factor).
While the peak flux is dimmer by a factor of $\sim$2-3, and the peak time
remains roughly the same (within $20\%$), 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. Due to the larger ejected mass
in NSBH merger, this late non-thermal signal can be observed to distances of
$\lesssim 200$Mpc for typical parameter values.

###### keywords:

(transients:) neutron star mergers – (transients:) black hole-neutron star
mergers – gravitational waves – radio continuum: transients – relativistic
processes

††pubyear: 2024††pagerange: Non-thermal emission from mildly relativistic
tidal ejecta of compact objects merger–B

##  1 Introduction

Binary mergers involving compact objects with a substantial mass ratio
($q>1.5$) or ones in which a prompt black hole (BH) is produced, such as NSBH
or BNS systems, yield expanding ejecta that attains mildly relativistic
velocities (for NSBH Rosswog, 2005; Just et al., 2015; Kawaguchi et al., 2015;
Kyutoku et al., 2018; Foucart et al., 2017; Foucart, 2020), (for BNS Radice et
al., 2018; Bernuzzi et al., 2020; Nedora et al., 2021b), $\gamma\beta>0.1$,
through tidal interaction between the BH/NS and the NS (see Bernuzzi, 2020;
Kyutoku et al., 2021, for review). Such ejecta is expected to produce non-
thermal electromagnetic (EM) emission on a time scale of $\sim$10 years (Nakar
& Piran, 2011; Sadeh et al., 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 disk around the BH, and the emission of EM radiation.
A disruptive NSBH merger occurs when the disruption radius, $R_{\text{dis}}$,
in 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). 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 should be considered in case of long-
lived neutron star remnant (Nedora et al., 2021b; Radice & Bernuzzi, 2023). In
cases where the mass ratio is considerable, $q>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}M_{\odot}$ (for various equation of states (EoSs) and mass
ratios, see Radice et al., 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>3$, provide some constraints over the
dynamical ejecta of such events (Rosswog, 2005; Just et al., 2015; Kyutoku et
al., 2015, 2018; Foucart et al., 2013, 2017; Chen et al., 2024). In cases of
tidal disruption, they imply the dynamical ejecta (notice that the spiral wind
component (Nedora et al., 2021b; Radice & Bernuzzi, 2023) is absent for NSBH
systems) has only a tidal component of $\sim$ $10^{-2}-10^{-1}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); 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>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 abovementioned 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).

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
et al., 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 it is between $0.1$ to
$14\left(\frac{D}{200\text{Mpc}}\right)^{3}\text{yr}^{-1}$ and for NSBH merger
it is between $0.1$ to
$1.5\left(\frac{D}{200\text{Mpc}}\right)^{3}\text{yr}^{-1}$ (LIGO Scientific
Collaboration et al., 2023; The LIGO Scientific Collaboration et al., 2024).
The fraction of EM bright events out of the detectable (via GW) 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 et al., 2023; The LIGO Scientific
Collaboration et al., 2024), suggesting a fraction of $\sim 10-30\%$.
Martineau et al. (2024) estimated the ejected mass in the observed event,
GW230529, to be negligible, $\sim 10^{-3}M_{\odot}-10^{-5}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 estimation
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 of an opening angle of between $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
§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, 2018; Foucart et al., 2017;
Bernuzzi et al., 2020; Nedora et al., 2021b). 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 under the
assumption of a power-law energy distribution for electrons, given by
$dn_{e}/d\gamma_{e}\propto\gamma_{e}^{-p}$. Here, $\gamma_{e}$ represents the
electron Lorentz factor within the plasma’s rest frame, while $p$ is within
the range of $2\leq p\leq 2.5$. Concurrently, 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 and non-relativistic shocks
(see Blandford & Eichler, 1987; Keshet & Waxman, 2005; Waxman, 2006; Keshet,
2006; Bykov & Treumann, 2011; Sironi et al., 2013; Pohl et al., 2020; Ligorini
et al., 2021; Kobzar 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 § 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 § 3. Then, we
describe the 2D relativistic-hydrodynamics numerical calculation and the post-
processing we use to produce light curves in § 4. The accuracy of the semi-
analytic method is estimated by comparison to the detailed numeric
calculations in § 5, while in § 6, we provide a few examples for light curves
and radio maps produced by the semi-analytic approach. Finally, our
conclusions are summarized in § 7.

##  2 Calculation framework

###  2.1 Geometry description

Figure 1: A schematic illustration showing a 2D slice of the 3D toroid ejecta
(described in § 2), and the forward-reverse shock structure that is formed and
described in § 3.1. The symmetry axis corresponds to the azimuthal symmetry of
toroids. $\theta_{\text{open}}$ is the angle between the equatorial plane and
the initial ejecta opening angle. $\theta_{\text{v}}$ is the angle between the
symmetry axis and the line of sight, such that $\theta_{\text{v}}=\pi/2$
corresponds to the equatorial plane.

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

$M(>\gamma\beta)=M_{0}\begin{cases}\left(\frac{\gamma\beta}{\gamma_{0}\beta_{0}}\right)^{-s_{\text{ft}}}&\text{for}\quad\gamma_{0}\beta_{0}<\gamma\beta,\\\
\left(\frac{\gamma\beta}{\gamma_{0}\beta_{0}}\right)^{-s_{\text{KN}}}&\text{for}\quad
0.1<\gamma\beta<\gamma_{0}\beta_{0},\end{cases}$ (1)

with parameter values characteristic of the results of numerical calculations
of the BNS ejecta; $0.3<\beta_{0}<0.5$, $5<s_{\text{ft}}<12$,
$0.5<s_{\text{KN}}<3$, and $10^{-6}<M_{0}<10^{-4}$. 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>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 merger with comparable masses (Kyutoku et al., 2015; Foucart et al., 2017;
Radice et al., 2018; Bernuzzi et al., 2020; Nedora et al., 2021b). 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 Eq. 1, and various
opening angles that are consistent with the results of NR simulations
($15^{\circ}-30^{\circ}$), 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

$\displaystyle\cos\xi$
$\displaystyle=\cos\theta_{\text{v}}\cos\theta-\sin{\theta_{\text{v}}}\sin\theta\sin\phi,$
(2) $\displaystyle\sin\varphi$
$\displaystyle=\frac{\cos\theta-\cos\theta_{\text{v}}\cos\xi}{\sin\xi\sin\theta_{\text{v}}},$

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

$\displaystyle R$ $\displaystyle=\sqrt{r^{2}+z^{2}},$ (3)
$\displaystyle\tan\theta$
$\displaystyle=\frac{\sqrt{(r\cos{\phi})^{2}+(r\sin{\phi}\cos{\theta_{\text{v}}}-z\sin{\theta_{\text{v}}})^{2}}}{r\sin{\phi}\sin{\theta_{\text{v}}}+z\cos{\theta_{\text{v}}}}.$

##  3 Semi-analytic calculation

A semi-analytic approach for calculating the expected synchrotron emission
from a toroidal ejecta with a velocity profile given by Eq. (1) is presented
below. This approach is valid for 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 will be shared following
a reasonable request to the corresponding author, offering a computation
method for the approximated dynamics, anticipated flux, and the radio image on
the sky.

###  3.1 Dynamics

The mass of a toroid ejecta, $M$, with an opening angle of
$\theta_{\text{open}}$ from the equatorial plane corresponds to
$\sin{\theta_{\text{open}}}M_{\text{iso}}$, where $M_{\text{iso}}$ is the
isotropic equivalent mass. The ejecta propagates through uniform ISM density,
forming a forward-reverse shock structure; see Fig. 1. As long as the reverse
shock did not cross the fast tail (the steep, fast part of the ejecta), 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: due to the steep decline of the ejecta profile
($s_{\text{ft}}\geq 5$) the injected energy (and momentum) to the shocked
medium increases faster with decreasing velocity (as the ejecta mass increases
faster with decreasing velocity). The swept ISM mass is dominated by the
freshly shocked material since the shock is propagating radially and, due to
momentum conservation, the momentum of the freshly shocked material is
dictated by added momentum (to the shocked medium) from the ejecta at the
position of the reverse shock, which is in the radial direction. 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. 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
modeling the dynamics of spherical mildly-relativistic ejecta with a velocity
profile given by Eq. (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
of $e(R,t),\rho^{\prime}(R,t)$ and $\gamma(R,t)$ (internal energy density,
proper mass density, and Lorentz factor correspondingly as a function of
radius and time in the rest frame of the ISM) behind the shocks front, with
different values for the forward and reverse shocks. We approximate the
shocked plasma as two separate uniform layers behind the shocks with the flow
profiles described above. 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,
$dn_{e}/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 $\theta_{\text{open}}$ of the ejecta.
Mathematically, this is expressed as:

$j^{\prime}_{\nu}=\begin{cases}j^{\prime}_{\nu}\left(e,\rho,\gamma\right)&\text{for}\quad\frac{\pi}{2}-\theta_{\text{open}}<\theta<\frac{\pi}{2}+\theta_{\text{open}},\\\
0&\text{else},\end{cases}$ (4)

where $j^{\prime}_{\nu}\left(e,\rho,\gamma\right)$ is the emissivity as a
function of the flow profiles, and the prime (′) notation indicates proper
frame quantities. Subsequently, the emissivity and intensity are computed
following the methodology outlined in Sadeh et al. (2024), where we only
considered frequencies above the self-absorption frequency and below the
cooling frequency, $\nu_{a}<\nu<\nu_{c}$. To calculate the flux, we consider
the fact that the azimuthal observer symmetry breaks in the case of
$\theta_{\text{v}}\neq 0$, so we calculate each azimuthal angle separately by

$\displaystyle F_{\nu}$ $\displaystyle=\frac{1}{D^{2}}\int\int
rI_{\nu}drd\phi.$ (5)

For the light curve temporal cutoff, we use the peak time estimation from
Sadeh et al. (2023),

$t_{\text{peak}}=550g(\beta_{0})\left(\frac{M_{0,-4}}{n_{-2}}\right)^{\frac{1}{3}}\text{days},\quad
g(\beta_{0})=\frac{1.5-\sqrt{0.25+2\beta_{0}^{2}}}{\gamma_{0}^{\frac{1}{3}}\beta_{0}},$
(6)

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}M_{\odot}$. For
$0^{\circ}\leq\theta_{\text{v}}\lesssim 20^{\circ}$ a delay is expected
between the peak time in the isotropic equivalent case given by Eq. (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)

$\displaystyle F_{\nu,\rm peak}$
$\displaystyle\approx\sin\theta_{\text{open}}F^{\text{iso}}_{\nu,\rm peak},$
(7) $\displaystyle\approx 4D_{26.5}^{-2}\varepsilon_{e,-1}^{p-1}f_{\rm
ft}(p)\varepsilon_{B,-2}^{\frac{p+1}{4}}n_{-2}^{\frac{p+1}{4}}M_{0,-4}\nu_{9.5}^{\frac{1-p}{2}}g_{q}\mu\text{Jy},$

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}$,

$f_{\rm ft}(p)=960\times 8.8^{-p}(2.5-0.7p)(p-1)^{2-p}(l_{p})^{p-1},$ (8)

and $F^{\text{iso}}_{\nu,\rm peak}$ is the peak flux in the equivalent
isotropic case.

##  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$ to $4/3$) depending on the internal
energy density

$\hat{\gamma}=\frac{4+\left(1+1.3\frac{e}{\rho^{\prime}c^{2}}\right)^{-1}}{3},$
(9)

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<e/\rho^{\prime}c^{2}<2$ (corresponding
to $1<\gamma<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$. 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 Eq. (1). The initial, $t=t_{0}$, radius of a mass shell
$M(>\gamma\beta)$ is determined by $r_{0}=\beta ct_{0}$. The expanding ejecta
is initially (at $t=t_{0}$) embedded (at $r>ct_{0}$) in a static, uniform,
cold (negligible pressure) gas.

###  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}<\nu<\nu_{c}$),
in the plasma 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,

$\Delta
F_{\nu}\left(t_{\text{obs}}=t-\frac{R\cos{\xi}}{c}\right)=\frac{j^{\prime}_{\nu}(t)}{\gamma(t)^{2}(1-\beta(t)\cos{\xi_{v}})^{2}}\cdot\frac{\Delta
V}{D^{2}},$ (10)

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_{v}$ is the angle between the cell velocity
direction and the line of sight. $\xi$ and $\xi_{v}$ pose different values for
different observation angles following Eq. (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 5.1, we verify our estimation regarding the dominant contribution to the
internal energy by the freshly shocked material within the initial opening
angle of the ejecta. This verification holds true as long as the reverse shock
hasn’t yet crossed the ejecta fast tail. In 5.2, we compare the light curves
produced by the full 2D numerical calculation to our semi-analytic scheme to
verify its validity. Finally, in 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 approve validation.

###  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\%$ of the internal energy, excluding rest
mass energy, is included. A few examples are given in Fig. 2. We indeed find
that the contribution to the internal energy is dominated by the material
within the initial opening angle (up to $\sim 30\%$ corrections) as long as
the reverse shock doesn’t cross the ejecta fast tail.

Figure 2: The opening angle (relative to the ejecta initial opening angle)
within which $90\%$ of the internal energy, excluding rest mass energy, is
included for various toroid ejecta parameters;
$\\{s_{\text{ft}},s_{\text{KN}},\gamma_{0},n,M_{R},\theta_{\text{open}}\\}=\\{6,1.5,1.2,3\times
10^{-2}\text{cm}^{-3},10^{-8}M_{\odot},30^{\circ}\\}$
$/\\{7,1.5,1.15,3\text{cm}^{-3},10^{-6}M_{\odot},20^{\circ}\\}$
$/\\{8,1.5,1.1,3\text{cm}^{-3},10^{-6}M_{\odot},15^{\circ}\\}$
$/\\{9,1.5,1.15,3\times 10^{-1}\text{cm}^{-3},10^{-7}M_{\odot},25^{\circ}\\}$.

###  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}<\nu<\nu_{c}$, since the radio observations are expected to be
in this part. In Fig. 3, we compare the light curve obtained semi-analytically
(for as long as the reverse shock doesn’t cross the fast tail, § 3) with the
full 2D numerical calculation (§ 4). The agreement is within $10$’s of percent
for a wide range of relevant values of $\\{\beta_{0},s_{\rm
ft},\theta_{\text{v}},\theta_{\text{open}}\\}$. As expected, a $\sim 20\%$
delay is observed for $\theta_{\text{v}}=0^{\circ}$.

Figure 3: A comparison between the flux calculated semi-analytically in (solid
lines) and the flux calculated in the full numerical calculation (dashed
lines), for various parameters
$\\{s_{\text{ft}},s_{\text{KN}},\gamma_{0},n,M_{R},\varepsilon_{e},\varepsilon_{B},p,\theta_{\text{open}}\\}=\\{6,1.5,1.2,3\times
10^{-2}\text{cm}^{-3},10^{-8}M_{\odot},10^{-1},10^{-1},2.1,30^{\circ}\\}$
$/\\{7,1.5,1.15,3\text{cm}^{-3},10^{-6}M_{\odot},10^{-1},10^{-1},2.2,20^{\circ}\\}$
$/\\{8,1.5,1.1,3\text{cm}^{-3},10^{-6}M_{\odot},10^{-1},10^{-1},2.3,15^{\circ}\\}$
$/\\{9,1.5,1.15,3\times
10^{-1}\text{cm}^{-3},10^{-7}M_{\odot},10^{-1},10^{-1},2.4,25^{\circ}\\}$. We
used a distance of $100$Mpc and a frequency of $1$GHz.

###  5.3 Sky image

To further validate our semi-analytic calculation, in Fig. 4, 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.

Figure 4: The fraction of flux that is emitted at different times by different
annuli, $df_{\nu}/dr$, for $p=2.2,s_{\text{ft}}=7,\gamma_{0}=1.15$. In solid
lines: full numerical 1D calculation. In dashed lines: our semi-analytic
scheme. $r$, the radius of the annulus (its transverse distance to the line of
sight) is normalized to $R_{R}=ct_{R}$ (approximately the radius at which the
post-shock plasma momentum drops to $\gamma\beta=1$), where
$t_{R}\equiv\left(\frac{M_{R}}{16\pi nm_{p}c^{3}}\right)^{\frac{1}{3}}$. The
$+$’s are the analytic estimate of the image radius, given in Sadeh et al.
(2024).

##  6 Results

In Figs. 5-6, 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,
$\\{s_{\text{ft}},s_{\text{KN}},\beta_{0},n,M_{R},\varepsilon_{e},\varepsilon_{B},\theta_{\text{open}}\\}=\\{7,1.5,0.3,10^{-2}\text{cm}^{-3},10^{-4}M_{\odot},10^{-1},10^{-2},20^{\circ}\\}$.

Figure 5: The observed sky radio map for a distance of $100$Mpc, frequency of
$1$GHz at $t=7110$days, for parameters expected for NSBH merger,
$\\{s_{\text{ft}},s_{\text{KN}},\beta_{0},n,M_{0},\varepsilon_{e},\varepsilon_{B},p,\theta_{\text{open}}\\}=\\{7,1.5,0.3,10^{-2}\text{cm}^{-3},8\times
10^{-2}M_{\odot},10^{-1},10^{-2},2.2,20^{\circ}\\}$, and for different
observation angles
$\theta_{\text{v}}=0^{\circ},30^{\circ},60^{\circ},90^{\circ}$, calculated by
the semi-analytic method described in §3. In blue: the analytic estimation
from Sadeh et al. (2024) for observed image radius (the image appears as a
narrow ring in the spherical case). Figure 6: The observed light curves for a
distance of $100$Mpc and a frequency of $1$GHz for parameters expected for
NSBH merger,
$\\{s_{\text{ft}},s_{\text{KN}},\beta_{0},n,M_{0},\varepsilon_{e},\varepsilon_{B},p,\theta_{\text{open}}\\}=\\{7,1.5,0.3,10^{-2}\text{cm}^{-3},8\times
10^{-2}M_{\odot},10^{-1},10^{-2},2.2,20^{\circ}\\}$, and for different
observation angles
$\theta_{\text{v}}=0^{\circ},30^{\circ},60^{\circ},90^{\circ}$, calculated by
the semi-analytic method described in §3. The black line represents the
analytical model from Sadeh et al. (2023) for a spherical ejecta with the same
isotropic equivalent mass. The black ’+’ is given by our estimations for the
peak time from Eq. (6) and peak flux from Eq. (7).

##  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
Eq. (1). We provided a semi-analytic calculation method for the synchrotron
emission, assuming emission only from the initial opening angle of the ejecta
(§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 Eqs. (6) and Eq. (7). Then, by calculating the full 2D hydrodynamics (§4),
we validated our semi-analytic approach in §5 for a wide range of ejecta
parameter values characteristic of merger calculation results; an equatorial
ejecta opening angle of $15^{\circ}<\theta<30^{\circ}$, a steep ($s_{\rm
ft}>5$, $\gamma_{0}\beta_{0}\geq 0.3$) mass distribution for mildly
relativistic velocities and a moderate mass distribution for the bulk of the
ejecta ($1<s_{\rm KN}<3$, $0.1<\gamma_{0}\beta_{0}\leq 0.3$). We also used
various values of electrons power-law distribution ($2.1<p<2.4$). 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>1.5$) mass ratio. We found that the
assumption that the freshly shocked material within the initial opening angle
of the ejecta has a dominant contribution to the internal energy during the
phase in which the reverse shock propagates through the fast tail is in
agreement with the full 2D numerical calculation (Fig. 2). Additionally, we
found that, in agreement with Eq. (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 (Eq. (6)) for spherical ejecta (with the same isotropic
equivalent mass) is within $20\%$ deviation from the peak time results of the
full numerical calculation for toroidal ejecta presented in this paper (Fig.
3). The radio map is significantly different between the toroidal and the
spherical case in which the image is a symmetric ring (Fig. 5), for
$\theta_{\text{v}}\geq 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), while 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. 5), providing the
possibility to constrain the viewing angle in future radio map observations.
For $0^{\circ}\leq\theta_{\text{v}}\leq 20^{\circ}$, the image is observed as
a narrow ring (Fig. 5), similar to the spherical case. The semi-analytic
scheme we presented in this work can be used to fit model parameters 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.

Different mass ratios between the compact objects correspond to different
ejecta geometries (Bernuzzi et al., 2020). Consequently, bounds over the
ejecta geometry can constrain this mass ratio. For example, in GW170817, the
estimation for the mass ratio varies between $1<q<2.6$ (for low spin
assumption, it is $1<q<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 $0.1M_{\odot}$
in disturbing merger systems (Bernuzzi et al., 2020), implying an isotropic
equivalent mass of $M\sim 0.3M_{\odot}$ (for
$\theta_{\text{open}}=20^{\circ}$). Following Eq. (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.01M_{\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. 6), 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. 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. 6.
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

  * Abbott et al. (2021) Abbott R., et al., 2021, The Astrophysical Journal Letters, 915, L5
  * Ayache et al. (2022) Ayache E. H., van Eerten H. J., Eardley R. W., 2022, Monthly Notices of the Royal Astronomical Society, 510, 1315
  * Bernuzzi (2020) Bernuzzi S., 2020, General Relativity and Gravitation, 52, 108
  * Bernuzzi et al. (2020) Bernuzzi S., et al., 2020, Monthly Notices of the Royal Astronomical Society, 497, 1488
  * Biscoveanu et al. (2023) Biscoveanu S., Landry P., Vitale S., 2023, Monthly Notices of the Royal Astronomical Society, 518, 5298
  * Blandford & Eichler (1987) Blandford R., Eichler D., 1987, Physics Reports, 154, 1
  * Blandford & McKee (1976) Blandford R. D., McKee C. F., 1976, The Physics of Fluids, 19, 1130
  * Bykov & Treumann (2011) Bykov A. M., Treumann R. A., 2011, The Astronomy and Astrophysics Review, 19, 42
  * Chen et al. (2024) Chen S., Wang L., Hayashi K., Kawaguchi K., Kiuchi K., Shibata M., 2024, Black hole-neutron star mergers with massive neutron stars in numerical relativity, https://ui.adsabs.harvard.edu/abs/2024arXiv240418714C
  * Foucart (2020) Foucart F., 2020, Frontiers in Astronomy and Space Sciences, 7
  * Foucart et al. (2013) Foucart F., et al., 2013, Physical Review D, 87, 084006
  * Foucart et al. (2017) Foucart F., et al., 2017, Classical and Quantum Gravity, 34, 044002
  * Fragione (2021) Fragione G., 2021, The Astrophysical Journal, 923, L2
  * Gompertz et al. (2023) Gompertz B. P., Nicholl M., Smith J. C., Harisankar S., Pratten G., Schmidt P., Smith G. P., 2023, Monthly Notices of the Royal Astronomical Society
  * Govreen-Segal & Nakar (2023) Govreen-Segal T., Nakar E., 2023, Monthly Notices of the Royal Astronomical Society, 524, 403
  * Hayashi et al. (2022) Hayashi K., Fujibayashi S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2022, Physical Review D, 106, 023008
  * Just et al. (2015) Just O., Bauswein A., Pulpillo R. A., Goriely S., Janka H.-T., 2015, Monthly Notices of the Royal Astronomical Society, 448, 541
  * Kawaguchi et al. (2015) Kawaguchi K., Kyutoku K., Nakano H., Okawa H., Shibata M., Taniguchi K., 2015, Physical Review D, 92, 024014
  * Kawaguchi et al. (2016) Kawaguchi K., Kyutoku K., Shibata M., Tanaka M., 2016, The Astrophysical Journal, 825, 52
  * Keshet (2006) Keshet U., 2006, Physical Review Letters, 97, 221104
  * Keshet & Waxman (2005) Keshet U., Waxman E., 2005, Physical Review Letters, 94
  * Klein (2023) Klein Y. Y., 2023, Construction of a Multidimensional Parallel Adaptive Mesh Refinement Special Relativistic Hydrodynamics Code for Astrophysical Applications, doi:10.48550/arXiv.2310.02331, https://ui.adsabs.harvard.edu/abs/2023arXiv231002331K
  * Kobzar et al. (2021) Kobzar O., Niemiec J., Amano T., Hoshino M., Matsukiyo S., Matsumoto Y., Pohl M., 2021, The Astrophysical Journal, 919, 97
  * Kyutoku et al. (2015) Kyutoku K., Ioka K., Okawa H., Shibata M., Taniguchi K., 2015, Physical Review D - Particles, Fields, Gravitation and Cosmology, 92
  * Kyutoku et al. (2018) Kyutoku K., Kiuchi K., Sekiguchi Y., Shibata M., Taniguchi K., 2018, Physical Review D, 97, 023009
  * Kyutoku et al. (2021) Kyutoku K., Shibata M., Taniguchi K., 2021, Living Reviews in Relativity, 24, 5
  * LIGO Scientific Collaboration et al. (2023) LIGO Scientific Collaboration and KAGRA Collaboration V. C., et al., 2023, Physical Review X, 13, 011048
  * Ligorini et al. (2021) Ligorini A., et al., 2021, Monthly Notices of the Royal Astronomical Society, 502, 5065
  * Lovelace et al. (2013) Lovelace G., Duez M. D., Foucart F., Kidder L. E., Pfeiffer H. P., Scheel M. A., Szilágyi B., 2013, Classical and Quantum Gravity, 30, 135004
  * Margalit & Piran (2015) Margalit B., Piran T., 2015, Monthly Notices of the Royal Astronomical Society, 452, 3419
  * Martineau et al. (2024) Martineau T., Foucart F., Scheel M., Duez M., Kidder L., Pfeiffer H., 2024, Black Hole-Neutron Star Binaries near Neutron Star Disruption Limit in the Mass Regime of Event GW230529, doi:10.48550/arXiv.2405.06819, https://ui.adsabs.harvard.edu/abs/2024arXiv240506819M
  * Meliani et al. (2007) Meliani Z., Keppens R., Casse F., Giannios D., 2007, Monthly Notices of the Royal Astronomical Society, 376, 1189
  * Mignone & McKinney (2007) Mignone A., McKinney J. C., 2007, Monthly Notices of the Royal Astronomical Society, 378, 1118
  * Mészáros (2002) Mészáros P., 2002, Annual Review of Astronomy and Astrophysics, 40, 137
  * Nakar (2007) Nakar E., 2007, Physics Reports, 442, 166
  * Nakar & Piran (2011) Nakar E., Piran T., 2011, Nature, 478, 82
  * Nedora et al. (2021a) Nedora V., Radice D., Bernuzzi S., Perego A., Daszuta B., Endrizzi A., Prakash A., Schianchi F., 2021a, Monthly Notices of the Royal Astronomical Society, 506, 5908
  * Nedora et al. (2021b) Nedora V., et al., 2021b, The Astrophysical Journal, 906, 98
  * Phinney (1991) Phinney E. S., 1991, The Astrophysical Journal, 380, L17
  * Piran (2004) Piran T., 2004, Reviews of Modern Physics, 76, 1143
  * Pohl et al. (2020) Pohl M., Hoshino M., Niemiec J., 2020, Progress in Particle and Nuclear Physics, 111, 103751
  * Pol et al. (2019) Pol N., McLaughlin M., Lorimer D. R., 2019, The Astrophysical Journal, 870, 71
  * Pol et al. (2020) Pol N., McLaughlin M., Lorimer D. R., 2020, Research Notes of the American Astronomical Society, 4, 22
  * Radice & Bernuzzi (2023) Radice D., Bernuzzi S., 2023, Secular Outflows from Long-Lived Neutron Star Merger Remnants, https://ui.adsabs.harvard.edu/abs/2023arXiv231009934R
  * Radice et al. (2018) Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., Roberts L. F., 2018, The Astrophysical Journal, 869, 130
  * Radice et al. (2020) Radice D., Bernuzzi S., Perego A., 2020, Annual Review of Nuclear and Particle Science, 70, 95
  * Rosswog (2005) Rosswog S., 2005, The Astrophysical Journal, 634, 1202
  * Rybicki & Lightman (1979) Rybicki G. B., Lightman A. P., 1979, Radiative processes in astrophysics. https://ui.adsabs.harvard.edu/abs/1979rpa..book.....R
  * Sadeh et al. (2023) Sadeh G., Guttman O., Waxman E., 2023, Monthly Notices of the Royal Astronomical Society, 518, 2102
  * Sadeh et al. (2024) Sadeh G., Linder N., Waxman E., 2024, Monthly Notices of the Royal Astronomical Society
  * Shibata & Taniguchi (2008) Shibata M., Taniguchi K., 2008, Physical Review D, 77, 084015
  * Sironi et al. (2013) Sironi L., Spitkovsky A., Arons J., 2013, The Astrophysical Journal, 771, 54
  * Synge (1957) Synge J. L., 1957, The Relativistic Gas, North-Holland Publ. Comp.
  * The LIGO Scientific Collaboration et al. (2024) The LIGO Scientific Collaboration the Virgo Collaboration the KAGRA Collaboration 2024, Observation of Gravitational Waves from the Coalescence of a $2.5-4.5~M_\odot$ Compact Object and a Neutron Star, %****␣main.bbl␣Line␣300␣****https://ui.adsabs.harvard.edu/abs/2024arXiv240404248T
  * Waxman (2006) Waxman E., 2006, Plasma Physics and Controlled Fusion, 48, B137
  * Zappa et al. (2023) Zappa F., Bernuzzi S., Radice D., Perego A., 2023, Monthly Notices of the Royal Astronomical Society, 520, 1481
  * Zhang & MacFadyen (2006) Zhang W., MacFadyen A. I., 2006, The Astrophysical Journal Supplement Series, 164, 255

##  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. 7). The deceleration imposed by the reverse
shock results in a reduction of kinetic energy in the ejecta and this energy
is subsequently transmitted to the heated ionized plasma within the shocked
ISM.

Figure 7: A schematic illustration showing the four regions of the forward-
reverse shock structure described in appendix A: (1) un-shocked ISM, (2)
shocked ISM, (3) shocked ejecta, (4) un-shocked ejecta. The pressure in the
unshocked regions is negligible compared to that in the shocked regions,
$p_{1},p_{4}\ll p_{2},p_{3}$. The pressure and velocity are approximated in
our analytic calculations as uniform within the shocked region, and the
density is approximated as uniform within the two regions separated by the
contact discontinuity.

We consider an ejecta with a mass profile given by
$M(>\gamma_{4}\beta_{4})=M_{R}(\gamma_{4}\beta_{4})^{-s}$ (numeric subscripts
relate to the region presented in Fig. 7) and assume that the ISM density is
uniform. The initial energy of the ejecta (excluding the rest-mass energy) is
given by

$E(>\gamma_{4}\beta_{4})=\int\frac{dM(>\gamma_{4}\beta_{4})}{d\gamma_{4}}(\gamma_{4}-1)c^{2}d\gamma_{4}.$
(11)

As mentioned in $\S$ 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),

$\displaystyle e$ $\displaystyle=(\gamma-1)\rho^{\prime}c^{2},$ (12)
$\displaystyle\rho^{\prime}$
$\displaystyle=\frac{\hat{\gamma}\gamma+1}{\hat{\gamma}-1}\rho,$

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(>\gamma_{4}\beta_{4})$, to
the mass densities in the rest frames by

$\displaystyle\rho_{1}$ $\displaystyle=\frac{3M_{\text{ISM}}}{4\pi R^{3}},$
(13) $\displaystyle\rho_{4}$
$\displaystyle=\frac{\gamma_{4}sM(>\gamma_{4}\beta_{4})}{4\pi(R_{\text{cd}}-\Delta_{\text{ej}})^{3}},$

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

$\displaystyle\frac{4\pi}{3}\rho_{1}R^{3}$
$\displaystyle=\frac{4\pi}{3}\rho_{1}\gamma_{2}\frac{\hat{\gamma}_{2}\gamma_{2}+1}{\hat{\gamma}_{2}-1}\left(R^{3}-(R-\Delta)^{3}\right),$
(14) $\displaystyle\Delta$
$\displaystyle=R\left(1-\frac{\left(\gamma_{2}\frac{\hat{\gamma}_{2}\gamma_{2}+1}{\hat{\gamma}_{2}-1}-1\right)^{\frac{1}{3}}}{\left(\gamma_{2}\frac{\hat{\gamma}_{2}\gamma_{2}+1}{\hat{\gamma}_{2}-1}\right)^{\frac{1}{3}}}\right).$

The thickness of the shocked ejecta layer is given by

$\displaystyle M(>\gamma_{4}\beta_{4})$
$\displaystyle=\frac{4\pi}{3}\rho_{4}\gamma_{2}\frac{\hat{\gamma}_{3}\gamma_{3}^{*}+1}{\hat{\gamma}_{3}-1}\left(R_{\text{cd}}^{3}-(R_{\text{cd}}-\Delta_{\text{ej}})^{3}\right),$
(15) $\displaystyle\Delta_{\text{ej}}$
$\displaystyle=R_{\text{cd}}\left(1-\frac{\left(\frac{\gamma_{2}\gamma_{4}s}{3}\cdot\frac{\hat{\gamma}_{3}\gamma_{3}^{*}+1}{\hat{\gamma}_{3}-1}\right)^{\frac{1}{3}}}{\left(1+\frac{\gamma_{2}\gamma_{4}s}{3}\cdot\frac{\hat{\gamma}_{3}\gamma_{3}^{*}+1}{\hat{\gamma}_{3}-1}\right)^{\frac{1}{3}}}\right),$

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

$(R_{\text{cd}}-\Delta_{\text{ej}})^{3}=R^{3}\left(1-\frac{\hat{\gamma}_{2}-1}{\gamma_{2}(\hat{\gamma}_{2}\gamma_{2}+1)}\right)\left(\frac{\left(\frac{\gamma_{2}\gamma_{4}s}{3}\cdot\frac{\hat{\gamma}_{3}\gamma_{3}^{*}+1}{\hat{\gamma}_{3}-1}\right)^{\frac{1}{3}}}{\left(1+\frac{\gamma_{2}\gamma_{4}s}{3}\cdot\frac{\hat{\gamma}_{3}\gamma_{3}^{*}+1}{\hat{\gamma}_{3}-1}\right)^{\frac{1}{3}}}\right).$
(16)

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

$\int T^{00}dV=\int((e+p)\gamma^{2}-p)dV.$ (17)

This can be expressed as

$\displaystyle E(>\gamma_{4}\beta_{4})$
$\displaystyle=E_{\text{k},2}+E_{\text{th},2}+E_{\text{k},3}+E_{\text{th},3},$
(18)

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

$\displaystyle E_{\text{k},2}$
$\displaystyle=(\gamma_{2}-1)M_{\text{ISM}}c^{2},$ (19) $\displaystyle
E_{\text{th},2}$
$\displaystyle=(\gamma_{2}-1)\left(\gamma_{2}\hat{\gamma}_{2}-\frac{\hat{\gamma}_{2}-1}{\gamma_{2}}\right)M_{\text{ISM}}c^{2},$
$\displaystyle E_{\text{k},3}$
$\displaystyle=(\gamma_{2}-1)M(>\gamma_{4}\beta_{4})c^{2},$ $\displaystyle
E_{\text{th},3}$
$\displaystyle=(\gamma_{3}^{*}-1)\left(\gamma_{2}\hat{\gamma}_{3}-\frac{\hat{\gamma}_{3}-1}{\gamma_{2}}\right)M(>\gamma_{4}\beta_{4})c^{2}.$

Eq. (18) enables us to derive an expression for $M_{\text{ISM}}$

$M_{\text{ISM}}=\frac{\frac{E(>\gamma_{4}\beta_{4})}{c^{2}}-M(>\gamma_{4}\beta_{4})\left(\gamma_{2}-1+(\gamma_{3}^{*}-1)\left(\gamma_{2}\hat{\gamma}_{3}+\frac{1-\hat{\gamma}_{3}}{\gamma_{2}}\right)\right)}{(\gamma_{2}-1)\left(1+\gamma_{2}\hat{\gamma}_{2}+\frac{1-\hat{\gamma}_{2}}{\gamma_{2}}\right)}.$
(20)

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

$(\hat{\gamma}_{2}\gamma_{2}+1)(\gamma_{2}-1)\rho_{1}c^{2}=(\hat{\gamma}_{3}\gamma_{3}^{*}+1)(\gamma_{3}^{*}-1)\rho_{4}c^{2}.$
(21)

We numerically solve Eq. (21) by using the expressions in Eqs. (13), (16) and
(20) 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 Eqs. (12), (13), and (20). The forward shock Lorentz factor,
$\Gamma$ is given by (Blandford & McKee, 1976)

$\Gamma^{2}=\frac{(\gamma_{2}+1)(\hat{\gamma}_{2}(\gamma_{2}-1)+1)^{2}}{\hat{\gamma}_{2}(2-\hat{\gamma}_{2})(\gamma_{2}-1)+2}.$
(22)

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

$\frac{dR}{dt}=\sqrt{1-\frac{1}{\Gamma^{2}}}.$ (23)

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)

$h=\frac{K_{3}(1/\Theta)}{K_{2}(1/\Theta)},$ (24)

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

$h=1+\frac{e}{\rho^{\prime}c^{2}}+\frac{p}{\rho^{\prime}c^{2}},$ (25)

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

$\hat{\gamma}=\frac{h-1}{h-1-\Theta}=\frac{\frac{K_{3}(1/\Theta)}{K_{2}(1/\Theta)}-1}{\frac{K_{3}(1/\Theta)}{K_{2}(1/\Theta)}-1-\Theta},$
(26)

and

$\Theta=(\hat{\gamma}-1)\frac{e}{\rho^{\prime}c^{2}}.$ (27)

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. 8, we compare
this result with the following approximation

$\hat{\gamma}=\frac{4+\left(1+Y\cdot\frac{e}{\rho^{\prime}c^{2}}\right)^{-1}}{3},$
(28)

for different values of $Y=1/1.1/1.3$. We also compare the result to the
approximation suggested in Ayache et al. (2022),

$\hat{\gamma}_{\text{Ayacha}}=\frac{4+\left(1+\frac{e}{\rho^{\prime}c^{2}}\right)^{-2}}{3}.$
(29)

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

$\displaystyle Y=\begin{cases}1,\quad\text{for
}9<\frac{e}{\rho^{\prime}c^{2}},\\\ 1.1,\quad\text{for
}2<\frac{e}{\rho^{\prime}c^{2}}<9,\\\ 1.3,\quad\text{for
}0<\frac{e}{\rho^{\prime}c^{2}}<2.\end{cases}$ (30)

Figure 8: Comparison between the numerically calculated adiabatic index
following Synge (1957) to our analytic approximation in Eq. (28), and to the
analytic approximation suggested in Ayache et al. (2022).