Non-thermal emission from mildly relativistic dynamical ejecta of neutron star mergers: spectrum and sky image

Binary neutron star mergers are expected to produce fast dynamical ejecta, with mildly relativistic velocities extending to 𝛽 = 𝑣 / 𝑐 > 0 . 6. In a preceding paper, we derived an analytic description of the time-dependent radio to X-ray synchrotron ﬂux produced by collisionless shocks driven by such fast ejecta into the interstellar medium, for spherical ejecta with broken power-law mass (or energy) distributions, 𝑀 ( > 𝛾𝛽 ) ∝ ( 𝛾𝛽 ) − 𝑠 with 𝑠 = 𝑠 KN at 𝛾𝛽 < 𝛾 0 𝛽 0 and 𝑠 = 𝑠 ft at 𝛾𝛽 > 𝛾 0 𝛽 0 (where 𝛾 is the Lorentz factor). Here, we extend our analysis and provide analytic expressions for the self-absorption frequency, the cooling frequency, and the observed angular size of the emitting region (which appears as a ring in the sky). For parameter values characteristic of merger calculation results - a "shallow" mass distribution, 1 < 𝑠 KN < 3, for the bulk of the ejecta (at 𝛾𝛽 ≈ 0 . 2), and a steep, 𝑠 ft > 5, "fast tail" mass distribution – the analytic results reproduce well (to tens of percent accuracy) the results of detailed numeric calculations, a signiﬁcant improvement over earlier order-of-magnitude estimates (based on extrapolations of results valid for 𝛾𝛽 ≪ 1).


INTRODUCTION
Mergers of binary neutron stars (BNS) are expected to produce mildly relativistic ( > 0.1) ejecta that propagate into the interstellar medium (ISM) (see Fernández & Metzger 2016;Radice et al. 2020, for reviews).The structure, velocity, and geometry of these ejecta are currently explored by numerical relativity (NR) simulations of BNS mergers.Those simulations suggest that a small fraction (10 −7 ⊙ to 10 −4 ⊙ ) of the ejecta mass (the "fast tail") may reach relativistic velocities ( > 0.6 ), while the bulk of the ejecta (10 −3 ⊙ to 5 × 10 −2 ⊙ ) is expected to attain velocities of ∼ 0.1 − 0.3 (Radice et al. 2018;Nedora et al. 2021a,b;Fujibayashi et al. 2023;Hajela et al. 2022;Radice et al. 2022;Rosswog & Korobkin 2024).Radioactive heating of the bulk of the ejecta is expected to produce IR-UV "Kilonova" (KN) emission (Li & Paczyński 1998), while the fast tail is predicted to generate a non-thermal radio to X-ray flux on a time scale of years by synchrotron emission from collisionless shocks driven by the fast ejecta into the interstellar medium (Nakar & Piran 2011;Sadeh et al. 2023).Ghosh et al. (2024) found no evidence of late radio emission in a follow-up campaign of five short gamma-ray bursts, which are considered to originate in compact objects mergers.The expected peak flux of such a late-time component following BNS mergers at the distances examined in Ghosh et al. (2024), >500Mpc, is ∼ 0.5 Jy (Sadeh et al. 2023), well below the Very Long Array (VLA) sensitivity.
The mass and velocity profile of the fast tail depends strongly on the parameters of the binary system and on the EoS.Depending on parameter values, the ejecta mass at > 0.6 varies, for example, between 10 −7 ⊙ and 10 −4 ⊙ , and the maximal velocity varies ★ E-mail: gilad.sade@weizmann.ac.il from ≈ 0.6 to > 3.This implies that measurements of the non-thermal emission driven by the fast ejecta may provide stringent constraints on model parameters.Two points should, however, be noted in this context.First, a reliable determination of constraints would require a significant reduction of the current large numerical uncertainty in the numerically calculated fast tail parameters, as reflected by the large variations of results between different simulations (e.g.Radice et al. 2018;Dean et al. 2021;Nedora et al. 2021a).
Second, if BNS mergers are the sources of (short) gamma-ray bursts (see Mészáros 2002;Piran 2004;Nakar 2007, for reviews), then they are expected to launch highly relativistic jets that, depending on the jet energy, opening and observing angles, may dominate the non-thermal emission driven by the "dynamical" (shock and tidally driven) ejecta obtained in NR BNS merger calculations.In the case of the observed non-thermal emission of GW170817 (Hallinan et al. 2017;Troja et al. 2017;Lyman et al. 2018), the swift decay (Troja et al. 2019;Fong et al. 2019;Lamb et al. 2019;Hajela et al. 2019;Troja et al. 2020;Nakar 2020;Makhathini et al. 2021;Balasubramanian et al. 2021Balasubramanian et al. , 2022) ) and the superluminal motion of the radio centroid suggest a relativistic jet as the driver of the observed non-thermal emission (noting that a ∼ 5 jet may account for the observations, hence a much higher Lornetz factor cannot be directly inferred, Mooley et al. 2018Mooley et al. , 2022)).Generally, due to the jet's lower mass and higher Lorentz factor, the late-time emission is expected to be dominated by the dynamical ejecta.The interpretation of future observations will require disentangling the two components.
In a preceding paper, Sadeh et al. (2023) (hereafter Paper 1), we considered the non-thermal radio to X-ray synchrotron emission produced by collisionless shocks driven by mildly relativistic ejecta expanding into an interstellar medium (ISM) with a uniform (number) density .Similar to other earlier papers (Nakar & Piran 2011;Kathirgamaraju et al. 2019;Nedora et al. 2023b), we did not provide an accurate description of synchrotron self-absorption, which is expected to affect the emission at low radio frequencies, and of the cooling of electrons on time-scales shorter than the expansion timescale, that is expected to be important for electrons emitting at hard X-ray bands (and beyond).The analytic expressions we provided hold, therefore, in the frequency range < < , where is the self-absorption frequency (below which absorption is significant), and is the synchrotron emission frequency of electrons with cooling time comparable to the expansion time.Note that the frequency of synchrotron radiation emitted by the lowest energy electrons, , is typically much smaller than (e.g.Eq. ( 15) of Paper 1).In this paper, we extend our analysis and provide analytic expressions for , , and for the observed angular size of the emitting region (which appears as a ring in the sky), which are calibrated using the results of numeric calculations.Observational identification of (in radio) and (in X-rays), and measurements of the angular size of the source (which may be possible with the Very Long Baseline Interferometry (VLBI) for sources at distances of ∼ 100Mpc, as the size is expected to reach ≈ 10 pc on years time scale) will provide essential constraints on model parameters.We provide approximations that are valid at times earlier than the peak time, peak at which the non-thermal emission at < < peaks, which is of the order of a few years.
Similar to Paper 1, we consider spherical ejecta with a broken power-law dependence of mass on momentum, with parameter values characteristic of the results of numerical calculations of the BNS ejecta; 0.3 < 0 < 0.5, 5 < ft < 12, 0.5 < KN < 3, and 10 −6 < 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).The effects of deviations from spherical symmetry, which may be significant for mergers of objects with a large, > 1.5, mass ratio, will be discussed in a follow-up paper under preparation.We show in appendix A that for the case of a broken power-law dependence of ejecta energy on momentum, accurate analytic approximations are obtained by substituting (similar to Paper 1) (3) The non-thermal flux is derived assuming that fractions and of the post-shock internal energy density are carried by nonthermal electrons and magnetic fields, respectively, and assuming that the electrons are accelerated to a power-law energy distribution1 , / ∝ − , where is the electron Lorenz factor (at the plasma rest frame) and 2 ≤ ≤ 2.5.This phenomenological description of the post-shock plasma conditions is supported by a wide range of observations and plasma calculations, for both relativistic and non-relativistic shocks (see Blandford & Eichler 1987;Waxman 2006;Bykov & Treumann 2011;Sironi et al. 2013;Pohl et al. 2020;Ligorini et al. 2021;Kobzar et al. 2021).
This paper is organized as follows.In § 2, we provide the analytic derivation of the observed emission ring radius, of the self-absorption frequency, and of the cooling frequency.In § 3 (and appendix B), we describe the numerical calculation scheme.The accuracy of the analytic formulae is determined by comparisons to the results of numeric calculations in § 4. A comparison to earlier results is given in § 5.A concise summary is given in § 6.

Image radius
As the ejected plasma propagates into the ISM, a forward-reverse shock structure is formed; the forward shock propagates into the cold medium ahead, while the reverse shock propagates into the ejecta and decelerates it.Radiation emitted from an angle with respect to the line of sight and observed at time was emitted when the forward shock radius was (Paper 1) where ej,RS ( ej,RS ) is the initial velocity (Lorentz factor) of the ejecta shell that was reached by the reverse shock (at the time the forward shock reached ).The perpendicular distance of the emitting plasma from the line of sight is given by ⊥ ( ) ≡ ( , ) sin .For a spherical explosion, the outer edge of the observed image is a circle, and its radius is obtained by solving ⊥ ( ) = 0.The solution is = , with cos = ej,RS , sin and In Paper 1 we showed that ej,RS ej,RS = where is the "relativistic mass" with momentum > 1, and is approximately the time at which the reverse shock crosses .Using these results, we have where the 1.05 factor is obtained by fitting to the results of numerical calculations.

Self-absorption frequency
We estimate the self-absorption frequency at time as the frequency for which = Δ = 1, where and Δ are the typical absorption coefficient and the typical path length traversed by photons through the shocked plasma, when the forward shock reached the radius , dominating the emission of radiation observed at time .
The contribution to the observed flux from emission at shock radius is dominated by plasma located at an angle where is the velocity of the shocked plasma right behind the forward shock at radius (Paper 1).Noting that the post-shock density is approximately 4 2 times the pre-shock density, we estimate the thickness Δ of the emitting layer as Using Eqs. 12 and 4, and approximating ≈ ej,RS , we have and To derive the self-absorption frequency, we first approximate the absorption coefficient Using the results from Paper 1 for the post-shock magnetic field strength, ( max , min are the maximal, minimal electron Lorentz factors respectively), and the relation ′ / ≈ 1.5 − 2 2 + 0.25 the frequency in the rest frame of the emitting plasma and the observed frequency, the absorption coefficient is given by (using the results of Rybicki & Lightman 1979, for power-law electron distribution) where Following the same power-law approximation used in Paper 1, appendix A, the optical depth at (observed) frequency is given by and can be approximated as (by Eq. ( 7)) The self-absorption frequency, defined by where the ft subscript implies (following the notation of Paper 1) that the result is valid as long as the reverse shock propagates through the fast tail.The 1.2 pre-factor and the − 3 5.5+ ft power-law index are obtained from fitting to the results of numeric calculations.

Cooling frequency
In Paper 1, we estimated the cooling frequency as a function of observed time.It is defined as the (observer frame) synchrotron frequency of electrons for which the (plasma frame) synchrotron loss time, 2 /( (4/3) ), is equal to the time measured at the plasma frame, / / 1.5 − 2 2 + 0.25 .Here, we numerically calibrate the values of the pre-factor and temporal power-law index of the result given in Paper 1 to better fit the results of numerical calculations, obtaining The dependence on { , , } is the same as that obtained in Paper 1.

Hydrodynamics
For an accurate description of the dynamics, we employ the hydrodynamic code described in Paper 1, which solves the full special relativistic 1D hydrodynamics equations.Our numerical 1D Lagrangian code utilizes a predictor-corrector scheme, incorporates artificial viscosity, and sets time steps using Courant-Friedrichs-Lewy (CFL) conditions.The code was validated by rigorous comparisons of its results to well-established solutions of benchmark problems (Sedov 1946;Martí & Müller 2003;Guzmán et al. 2012).The equation of state we use is = ( ˆ − 1) , where the adiabatic index ˆ ( / ) varies within the range of 4/3 to 5/3, following Synge (1957).The initial ejecta density and velocity profiles are given in Eq. ( 1), and the initial, = 0 , radius of a mass shell (> ) is determined by 0 = 0 .The expanding ejecta is initially (at = 0 ) embedded (at > 0 ) in a static, uniform, cold (zero pressure) gas with number density .The ejecta comprises 700 numeric cells out of a total of 5000 cells.Convergence is verified by doubling the resolution in the radial direction and reducing the time steps by a factor of 2, yielding indistinguishable results.

Radiation
To calculate the emission of radiation, noting the azimuthal symmetry, we use cylindrical coordinates, and , where the axis is aligned with the line of sight, such that the radiation emitted by a cell at lab arrives to the observer at obs = lab − . ( We divide space into × ( , ∼ few 1000 each) cells.Using the flow fields ( , , ) obtained in the hydrodynamic calculations, the synchrotron emission is calculated in the following manner (full details are given in appendix B): We assume that fractions and of the internal energy density are held by magnetic fields and electrons and that the electrons are accelerated to a power-law distribution / ∝ − .The emissivity, , and the absorption coefficient, , in each location and time, are defined in the plasma frame, and then Lorentz boosted to the lab frame.
To consider the effect of synchrotron cooling of high energy electrons, we numerically solve (as described in appendix B) the energy evolution of the electrons, determine the (spatially and temporally) dependent Lorentz factor above which the electron distribution is strongly suppressed by cooling, and calculate the radiation emitted by a power-law energy distribution of electrons (properly shifted down in energy due to adiabatic expansion) truncated at .
The contribution to the observed intensity of each cell is given by where Δ = Δ and is the optical depth along the path to the observer (In calculating obs + we take into account the evolution of the shocked plasma during the photons' propagation through it, with properly evaluated = obs + along the path).To accelerate numeric convergence, the factor Δ obs + in Eq. ( 24) is replaced with 1 − −Δ ( obs + ) .
The flux and fraction of the flux originating from different annuli are defined as where is the distance to the observer.The convergence of our numeric radiation calculations is verified by doubling the resolution in both and axes, yielding indistinguishable results.

COMPARISON OF ANALYTIC AND NUMERIC RESULTS
In comparing our numeric and analytic results, we have explored a wide range of values of the dimensionless parameters determining the hydrodynamic behavior, 5 < ft < 10, 0.3 < 0 < 0.5, and the range 2 < < 2.4 of the electron spectral index.The dependence of the dynamics on the dimensional parameters, { 0 , , }, follows directly from dimensional considerations, while the dependence on and is analytically obtained straightforwardly.Since we provide results for times earlier than the peak time, the time at which the reverse shock reaches the shell with initial momentum 0 0 , our results are independent of KN .

Sky image
In Fig. 1, we show an example of the flux emitted from different annuli obtained numerically along with the analytic estimate of the image radius, Eq. ( 9).Due to relativistic beaming and time arrival effects, the image is a relatively narrow ring.In Fig. 2, we compare the image radius obtained numerically (defined by the radial position of the peak of / ) with the analytic estimate, Eq. ( 9).The agreement is to within a few percent for a wide range of relevant values of { 0 , ft }.

Spectrum
Fig. 3 presents an example of spectra obtained by our numeric calculations, compared with the analytic estimates of the self-absorption and cooling frequencies, Eqs. ( 21) and ( 22) respectively.In Figs. 4  and 5, we compare the cooling and self-absorption frequencies obtained numerically with those given by the analytic approximations, Eqs. ( 21) and ( 22).We define the self-absorption[cooling] frequency in the numerical calculation as the minimal frequency at which the spectrum frequency dependence reaches log( )/ log( ) ≤ 0[(1 − 2 )/4].The values 0 and (1 − 2 )/4 are chosen as they are intermediate values between those of the different spectral regimes: The agreement is within 10's of percent for a wide range of relevant values of { 0 , ft }.
The ratio between our (numerically verified) analytic result and that of Nakar & Piran (2011) For typical values of 0 = 0.3, the order of magnitude estimate of Nakar & Piran ( 2011) is accurate to within a factor of ∼ 5.

CONCLUSIONS
Analytic expressions were derived in § 2 for the image radius, the self-absorption frequency, , and the cooling frequency, , of the non-thermal emission from a collisionless shock driven into the ISM by mildly relativistic spherical ejecta with broken power-law mass or energy distributions, given by Eqs.(1) or (2), at times earlier than the peak time of the emission at < < (Eqs.( 9), ( 21) and ( 22)).The analytic model results were compared in § 4 to the results of 1D numeric calculations for a wide range of ejecta parameter values characteristic of merger calculation results, 5 < ft < 10, 0.3 < 0 < 0.5, and for the range 2 < < 2.4 of the electron spectral index (The dependence of the dynamics on the dimensional parameters, { 0 , , }, follows directly from dimensional considerations, while the dependence on and is analytically obtained straightforwardly).We showed that the analytic model expressions reproduce the results of numeric calculations with tens of percent accuracy; see Figs. 1-5.This is a significant improvement over earlier order-of-magnitude  estimates, based on extrapolations of results valid for ≪ 1 to ≈ 1.The results presented in this paper, combined with the analytic results derived in Paper 1 for the peak time, peak flux, and temporal power-law indices of the flux rise and decline, will enable to constrain the parameters of the model, { ft , KN , 0 , , 0 0 , , , }, using future observational data.We note that an analytic description is essential for model parameter inference from data since the model depends non-trivially on several parameters, { ft , KN , 0 0 , }, and numerical calculations of the model predictions for each set of parameter values requires significant computational resources and time.
The spectral flux peaks at ∼ , typically in the range of 10 − 100 MHz, and suppressed beyond , typically at ∼ 10 10 GHz.For sources at a distance of ∼ 100 Mpc, the VLA may measure the self-absorption frequency, while the cooling frequency may be determined using Chandra X-ray data.

Figure 1 .Figure 2 .
Figure 1.Left panel: The fraction of flux that is emitted at different times by different annuli, / , (obtained numerically) for = 2.2, ft = 7, 0 = 1.15., the radius of the annulus (its transverse distance to the line of sight) is normalized to = (approximately the radius at which the post-shock plasma momentum drops to= 1), where is defined in Eq. (8).The +'s are the analytic estimate of the image radius, given by Eq. (9).Right panel: A corresponding intensity map at observed time = 300days ( and / are related through Eq. 25).