Shock-powered radio precursors of neutron star mergers from accelerating relativistic binary winds

The spin-down Poynting luminosity of an isolated rotating magnetized neutron star can be expressed as (e.g., Bucciantini et al. 2006)

where $\Omega$ is the angular rotation rate, $\mu\equiv B_{\rm d}R_{\rm ns}^{3}$ is the dipole moment of the neutron star of radius $R_{\rm ns}$, and $B_{\rm d}$ is the surface dipole field strength. Here $f_{\Phi}$ is the fraction of the magnetic flux threading the NS surface which is open to infinity, normalized to its value $\tilde{f}_{\Phi}\approx R_{\rm ns}/2R_{\rm L}$ for an isolated dipole wind with aligned magnetic and rotation axes in the limit $R_{\rm L}\gg R_{\rm ns}$, where $R_{\rm L}\equiv c/\Omega$ is the light cylinder radius.

One effect of the close orbiting binary companion of radius $R_{\rm c}\sim R_{\rm ns}$ in a binary of semi-major axis $a\ll R_{\rm L}$ may be to open additional closed field lines (e.g., Palenzuela et al. 2013) for finite ranges of companion’s resistivity (Lai, 2012), such that the amount of magnetic flux is enhanced over that of an isolated pulsar by an amount :

Compared to the case of an isolated pulsar, the open flux fraction has been reduced by a factor of $\sim 2R_{\rm c}/(4\pi a)$ to account for the azimuthal angle subtended by the binary companion, and increased by a factor of $\approx R_{\rm L}/a=c/(2\Omega_{\rm orb}a)$ because field lines crossing the equatorial plane exterior to the binary separation $a$ (instead of the light cylinder radius) are now open, where $\Omega_{\rm orb}=(GM_{\rm tot}/a^{3})^{1/2}$ is the orbital frequency of the binary of total mass $M_{\rm tot}$ (see Fernández & Metzger 2016 for further discussion). This corresponds to the “U/u” case simulated by Palenzuela et al. (2013), which as already mentioned is the most likely magnetic field configuration characterizing merging NS binary progenitors. A similar configuration should characterize BH-NS mergers because the BH is not magnetized.

Substituting Eq. (2) into (1), the power of the binary wind,

is an extremely sensitive function of the binary separation, $\dot{E}\propto a^{-7}$, where in the final line we have assumed an equal mass binary $M_{\rm tot}=2M_{\rm ns}$ with $M_{\rm ns}=1.4M_{\odot}$, and $R_{\rm ns}=12$ km. Lai (2012) derive a similar expression similar to Eq. (3) as the maximum power which can be dissipated in the closed “binary circuit” before the generation of a strong toroidal magnetic field from the poloidal current leads to magnetic field line inflation and reconnection. Here we hypothesize that an order-unity fraction of this maximum power emerges as a quasi-steady magnetized wind from the binary equatorial plane.

Even if the companion does not itself open the magnetic flux bundles, a separate source of electromagnetic power arises due to the magnetic dipole radiation generated by the orbital acceleration of the magnetized primary, which is independent of the properties of the companion. This gives a contribution to the power of (Ioka & Taniguchi, 2000; Carrasco & Shibata, 2020)

similar in its normalization and scaling with binary separation as Eq. 3.

The power emitted by the “binary wind” (Eqs. 3,4) eclipses that of the single more highly-magnetized pulsar ($\dot{E}$ from Eq. 1 with $f_{\Phi}=\tilde{f}_{\Phi}$ and where now $\Omega=2\pi/P$ for the pulsar spin-period $P$) for semi-major axes smaller than a critical value,

The binary is driven to coalescence by gravitational waves, such that its semi-major axis shrinks with time $t$ according to,

where $a(t=0)\equiv a_{0}$ is the initial binary separation, $t_{\rm m,0}\equiv t_{\rm m}(a_{0})$ is the time to merge starting from $a_{0}$, and

where in the first line we have assumed equal bodies of mass $M_{\rm ns}=M_{\rm tot}/2$ and in the second line we have taken $M_{\rm ns}=1.4M_{\odot}$. It will prove convenient to introduce the time until merger as an alternative variable,

in terms of which

Thus, from Eq. (3) we have $\dot{E}\propto a^{-7}\propto\Delta t^{-7/4}$ for times

The open magnetic flux lines will carry rest-mass as well as energy, likely in the form of electron/positron pairs. For ordinary pulsar winds, the wind mass-loss rate, $\dot{M}$, is typically thought to scale with the Goldreich-Julian (GJ) value,

where $\mu_{\pm}$ is the pair multiplicity, and

are the GJ current and charge density at the light cylinder, respectively.

The dynamics of the binary-induced winds considered here, as well as the external electromagnetic and photon environment, are drastically different than isolated pulsar winds (e.g., Wada et al. 2020). For example, the outflow may be loaded by plasma generated via $\gamma-\gamma$ annihilation due to high energy radiation from reconnection or other forms of dissipation in the magnetosphere (e.g., Metzger & Zivancev 2016). Nevertheless, performing the naïve exercise of replacing the neutron star rotation rate with that of the binary orbit in Eq. (11), i.e. $\Omega\rightarrow\Omega_{\rm orb}$, then one finds $\dot{M}_{\rm GJ}\propto a^{-3}\propto\Delta t^{-3/4}$.

Throughout this paper we parameterize the dependence of the wind mass-loss rate on binary separation as

where $m$ is a free parameter which takes on the value $m=3$ if $\dot{M}\propto\dot{M}_{\rm GJ}$ but is highly uncertain as it depends on the mechanism of wind pair-loading.

We define a wind “mass loading” parameter $\eta$ according to

A wind that fully converts its energy to kinetic form can achieve a bulk Lorentz factor $\Gamma\simeq\eta$ by large radii, but in general $\Gamma\ll\eta$ in cases where the bulk remains in Poynting flux or is otherwise lost (e.g., as radiation). Substituting Eqs. (3), (11) into Eq. (14), we obtain

for the case with $\dot{M}=\dot{M}_{\rm GJ}$. For young, high-voltage pulsars such as the Crab pulsar, observations and theory show that $\mu_{\pm}\lesssim 10^{5}$ (e.g., Kennel & Coroniti 1984; Timokhin & Harding 2019), which applied to our case near the end of the inspiral (e.g., $a\lesssim 4R_{\rm ns}$) would imply $\eta\gtrsim 10^{5}$ for $B_{\rm d}=10^{12}$ G. However, this may substantially under-estimate the mass-loading in the NS merger case (e.g., due to $\gamma-\gamma$ pair creation in the inner magnetosphere), in which case $\eta$ would be substantially less.

For $\dot{M}=\dot{M}_{\rm GJ}$ and $\mu_{\pm}={\rm constant}$ ($m=3$), Eq. (15) predicts that $\eta\propto a^{-4}$ will increase approaching the merger. Even in the more general case $m\neq 3$ (Eq. 14), $\eta$ increases with time as long as $m<7$. In this case the binary-driven wind will grow in both power and speed, inevitably giving rise to self-interaction and the generation of internal shocks within the wind.

The strength of the internal shocks depends not just on the evolution of $\eta$, but also on the four-velocity attained by the wind at large radii $r$ well outside its launching point near the binary orbit. In models of axisymmetric pulsar winds (Goldreich & Julian, 1970), the bulk Lorentz factor achieved near the fast magnetosonic surface outside the light cylinder is $\Gamma=\eta^{1/3}$, while the asymptotic magnetization of the wind (ratio of Poynting flux to kinetic energy flux) is $\sigma=\eta^{2/3}\gg\Gamma$, i.e. the acceleration is highly inefficient and most of the energy remains trapped in magnetic form. However, observations of pulsar wind nebulae reveal that—at least by the location of the wind termination shock—the wind has nearly fully converted its magnetic energy into bulk kinetic energy, i.e. $\Gamma\simeq\eta$; $\sigma\ll 1$ (Kennel & Coroniti, 1984). Several ideas have been proposed to explain this anomalous acceleration (the “$\sigma$ problem”; Lyubarsky & Kirk 2001; Sironi & Spitkovsky 2011; Porth et al. 2013; Zrake & Arons 2017; Cerutti et al. 2020), but none are universally agreed upon. A similar inference of rapid and efficient acceleration of pulsar winds has been invoked to account for the gamma-ray flares of the Crab pulsar (Aharonian et al., 2012) and of gamma-ray binaries (Khangulyan et al., 2012).

In this paper, we assume the binary-driven wind manages to solve its “$\sigma$ problem” and hence to achieve a bulk Lorentz factor close to the maximal value, i.e.

by the radii $r\lesssim r_{\rm sh}$ at which internal shocks occur. The wind self-interaction will then be mediated by strong shocks of moderate magnetization, which can be approximately modeled as hydrodynamical (i.e. neglecting magnetic fields). If the winds instead were to maintain $\sigma\gg 1$ to large radii, then the resulting strongly magnetized shocks could be considerably weaker. However, as we show below (eq. 20), under a reasonable assumptions for acceleration of the wind, we naturally expect $\sigma\sim 1$ by the radius of internal shocks. Nonetheless, the shock can still be strong for $\sigma\gg 1$ if the post-shock region is weakly magnetized—likely as a result of shock-driven magnetic reconnection, which can be faithfully probed with resistive magnetohydrodynamical simulations. We leave a detailed study of the magnetized shock case to future work, though we comment on how it would effect our predictions in $\S\ref{sec:discussion}$.

The accelerating binary wind will collide with the earlier, slower phases of the wind ejecta, generating magnetized ultra-relativistic shocks. In this section we provide analytic estimates of this interaction and the resulting radio emission in particular, tractable limits, which we later test with numerical simulations (§4).

Given the rapid rise of the outflow power $\dot{E}\propto a^{-7}$, most of the total wind energy,

is released near the end of the inspiral, where $a_{\rm f}$ is the final separation (e.g., as defined by the physical collision between the NSs) and $\dot{E}_{\rm f}\equiv\dot{E}[t_{\rm f}]$. This energy will emerge in a “fast shell” of width $t_{\rm f}\equiv t_{\rm m}[a_{\rm f}]\sim 1$ ms and Lorentz factor $\Gamma_{\rm f}\equiv\Gamma[t_{\rm f}]$. The shell will initially expand ballistically into the pre-existing wind released by the earlier phases of the inspiral, before sweeping up enough mass to decelerate (depending on the time evolution of $\dot{M}$). In what follows, it is convenient to express the energy-loss rate and Lorentz factor of the binary wind as

The fast shell released at the end of the merger initially expands with $\Gamma_{\rm f}\gg 1$, reaching a radius $r_{\rm sh}\simeq ct$ by a time $t$ after the merger. Throughout the remainder of this section we ignore deceleration of the fast shell; the conditions for this “freely coasting” evolution to be valid are estimated below.

A question immediately arises: Did the inspiral wind material met by the fast shell at radius $r_{\rm sh}$ expand freely prior to that point, or had it earlier collided/shocked with the even slower material ahead of it? Let us assume the former and check the conditions under which this assumption is violated. If the upstream inspiral wind has not decelerated appreciably since injection, then the slow wind released $\Delta t$ prior to the merger will reach $r_{\rm sh}$ by a time $t$ after the merger given by

where the second equality makes use of Eq. (18). This expression also assumes that $\Gamma[\Delta t]$ is much smaller than the fast shell $\Gamma_{\rm f}$ and hence is not accurate for small $\Delta t$. For $m\gtrsim 5$, the wind met by the fast shell at $r_{\rm sh}$ is “pristine” and corresponds to that ejected at the time $\Delta t$ satisfying Eq. (19). By contrast, for $m\lesssim 5$ we obtain the seemingly unphysical result that $\Delta t$ increases with increasing $t$, which instead shows that the inspiral wind will hit earlier ejecta prior to interacting with the fast shell (and hence Eq. 19 cannot be used).

As discussed near the end of the previous section, we have assumed that the wind has largely converted its Poynting flux into bulk kinetic energy by the radius of internal shock interaction. We now check whether such efficient acceleration is possible under a specific assumption regarding the acceleration mechanism of the wind. Following Granot et al. (2011), the Lorentz factor of the impulsively ejected fast shell will grow with radius from the wind launching point $r_{0}\simeq 2a$ (near the fast magnetosonic point, where $\Gamma\simeq\eta^{1/3}$; Goldreich & Julian 1970) as

thus reaching its maximal value ($\Gamma\simeq\eta;\sigma\lesssim 1)$ by the radius

Assuming efficient acceleration, the radius of the shock generated by the fast shell corresponding to an observer time $\sim t_{\rm f}$ is given by (eq. 19),

Thus,

where in the final line we have taken typical values $a\sim 12$ km, $t_{\rm f}\sim 1$ ms. Thus, we generally expect that the shell will have accelerated to nearly its maximal Lorentz factor ($\Gamma\sim\eta$) but also that the residual magnetization of the wind material will not be too low $\sigma\sim 1$.

We now check the condition for the fast shell to appreciably decelerate by sweeping up pristine (previously unshocked) inspiral wind. The rest mass swept up by the fast shell as it meets material released at times $\lesssim\Delta t$ prior to the merger grows as $M_{\rm rest}\sim\dot{M}\Delta t\propto\Delta t^{(4-m)/4}$ and hence is dominated by the matter already present in the fast shell for $m>4$. On the other hand, what matters for deceleration is the swept-up comoving inertial mass $M_{\rm th}\propto M_{\rm rest}\gamma_{\rm th}$ which accounts also for the thermal Lorentz factor $\gamma_{\rm th}\gg 1$ of relativistically hot gas behind the shock. The latter should scale with the Lorentz factor of the forward shock $\Gamma_{\rm sh}$, i.e. $\gamma_{\rm th}\propto\Gamma_{\rm sh}\simeq\Gamma_{\rm f}/2\Gamma\propto\Delta t^{(7-m)/4}$, where we have made use of Eq. (19). Thus, $M_{\rm th}\propto\Delta t^{(11-m)/2}$ will grow with increasing $\Delta t$ for $m<11/2=5.5$, leading to a violation of our assumption of constant velocity of the fast shell and hence also a violation of the assumptions entering Eq. (19). The regimes of wind interaction are summarized in Table 1.

We now estimate the properties of the maser emission in the regime of freely-coasting evolution of the fast shell ($5.5<m<7$). The upstream rest-frame pair density met by the fast shell, as it meets inspiral wind released at time $\Delta t$, is given by

where in the final line we have used Eq. (19).

Observer time is connected to lab-frame time according to the standard expression,

such that Eq. (24) becomes,

For moderate upstream magnetization $\sigma\lesssim 1$, the spectral energy distribution of the maser peaks at a few times the plasma frequency of the upstream medium $\nu_{\rm p}=(4\pi n_{\pm}e^{2}/m_{\rm e})^{1/2}$ (Plotnikov & Sironi, 2019), which boosted into the observer frame corresponds to a peak emission frequency

where

where we have used Eq. (17) and take $t_{\rm f}=t_{\rm m}[a_{f}]\approx 1.2$ ms (Eq. 7). The spectral peak of the radio emission will thus start at high frequencies at $t_{\rm f}$ and drift lower after that time.

The luminosity of the maser emission is given by

where $L_{\rm e}$ is the electron kinetic luminosity entering the forward shock seen by the external observer, $f_{\xi}\sim 10^{-3}-10^{-2}$ is the maser radiative efficiency defined relative to the kinetic energy of the electrons (Plotnikov & Sironi, 2019), and $f_{\rm b}<1$ is geometric beaming factor due to the fraction of the total solid angle subtended by the binary outflow (e.g., the wedge near the binary equatorial plane; Fig. 1).

The shock luminosity in the observer frame is given by

Even though the fast shell is not significantly decelerated for $m>5.5$, the total energy dissipated by the shock, $\int_{t_{\rm f}}L_{\rm e}{\rm d}t_{\rm obs}\sim\dot{E}_{\rm f}t_{\rm f}$, is of the same order of magnitude as the total energy of the fast shell.

We consider the wind launched from the binary during the final stages of the inspiral, as the semi-major axis shrinks from some specified initial distance $a_{\rm 0}$ to a final size $a_{\rm f}$. We take $a_{\rm f}=2R_{\rm ns}=24$ km in our simulations, appropriate to NS-NS mergers. For yet smaller separations $a\lesssim a_{\rm f}$, additional complications arise (e.g., deviations from point-mass inspiral picture, tidal mass-loss from the star(s)—baryon loading in the wind, non-dipolar contributions to the magnetic field, etc.), and therefore we ignore the very final phase.

The binary separation at the beginning of the simulation, $a_{\rm 0}$, must be chosen large enough so as to generate a sufficient upstream medium for the subsequent wind (the bulk of the energy of which is released near the end of the inspiral) to interact with and, potentially, be decelerated by. On the other hand, larger values of $a_{\rm 0}$ require a longer simulation time prior to and after the merger to follow the full gravitational wave inspiral and shock interaction, respectively. In our numerical simulations, we take $a_{\rm 0}/a_{\rm f}=3-10=6-20~{}R_{\rm ns}$ as a compromise between these considerations. This is still a factor $\gtrsim$ 2 smaller than the total range of binary separation over which the binary wind dominates that of the single pulsar (Eq. 5). Defining the start of the calculation as $t=0$, the merger takes place at $t=t_{\rm m}[a_{0}]\equiv t_{\rm m,0}\approx 0.1-12$ s (for $a_{0}/a_{\rm f}=3-10$), and the physical part of the wind injected from the inner boundary of the simulation is assumed to terminate just prior to the merger at $t\simeq t_{\rm m,0}-t_{\rm f}\approx t_{\rm m,0}$, where $t_{\rm f}\approx 1.2$ ms.

Since the shock interaction takes place at radial distances $\gg a_{\rm f}$, we fix the inner boundary of our simulation grid at $r_{\rm in}=a_{\rm f}$, rather than following the binary contraction self-consistently. In general, the injected wind power, $\dot{E}_{\rm in}$, can be divided into kinetic and magnetic components,

where $\sigma_{\rm in}\equiv\dot{E}_{\rm mag}/\dot{E}_{\rm kin}$ is the magnetization, and

where $B_{\rm in}(t)$ and $\rho_{\rm in}(t)$ are the lab-frame magnetic field strength and rest-frame mass density injected at the inner boundary, and $\Gamma_{\rm in}$ is the bulk Lorentz factor of the injected wind.

As discussed in §2.1, we assume that the wind is able to accelerate over a short radial baseline, effectively immediately converting its energy entirely into kinetic form. Hence we take $\sigma_{\rm in}=0$ and perform purely hydrodynamical simulations of the shock interaction (however, note that a moderate value of $\sigma=0.1$ is assumed in calculating the properties of the shock maser emission; §4.1.2). Future work will perform MHD simulations exploring the more general case of highly magnetized winds $\sigma_{\rm in}\gg 1$.

As summarized in Fig. 2 and Table 2, the time-dependent wind properties are specified via boundary conditions at the inner domain of the simulation grid. Following Eq. (18), the kinetic luminosity of the injected wind $\dot{E}=\dot{E}_{\rm kin}$ evolves according to

where the initial wind power $\dot{E}_{0}\approx\dot{E}_{\rm in}[t_{\rm f}](a_{\rm 0}/a_{\rm f})^{-7}$ is chosen such that the final wind power $\dot{E}_{\rm in}[t_{\rm f}]\equiv\dot{E}_{\rm f}\approx(10^{3}-10^{7})\dot{E}_{0}$ following Eq. (3), for different values of $m$. Likewise, the mass-loss rate of the wind injected from the inner boundary evolves as (Eq. 18)

where $\dot{M}_{\rm 0}=\dot{E}_{\rm 0}/\Gamma_{\rm 0}c^{2}$ (Eq. 14) and $\Gamma_{0}=1.5$. For numerical stability, we gradually taper the wind power and mass-loss rate just after the merger ($t\gtrsim t_{\rm m,0}$), instead of abruptly shutting the engine off. This is achieved in Eqs. (38, 39) by multiplying the wind power by a tapering function,

We have checked that the precise form of the tapering does not significantly affect the overall shock dynamics relevant to the radio emission.

We vary the power-law index $m\in[3,7]$ in order to explore a diversity of wind mass-loading behavior (§2.2). The values of $\dot{E}_{0}$ and $\dot{M}_{0}$ are chosen such that the initial outflow Lorentz factor $\Gamma_{\rm in}(t=0)=\Gamma_{\rm 0}+1.01\simeq 2.5$ for all models. For our chosen set-up, different values of $m$ and initial inspiral radius $a_{\rm 0}/a_{\rm f}$ result in different values for the peak $\Gamma_{\rm in}[t_{\rm f}]\simeq\Gamma_{\rm f}\sim 3-115$ of the outflow achieved near the end of the inspiral (Fig. 2; Table 2). These relatively modest peak Lorentz factors were chosen for purposes of numerical stability ($\S\ref{sec:code}$), despite being much smaller than are likely achieved in the actual wind (e.g., Eq. 15 and surrounding discussion). However, insofar as we are able to capture the shock dynamics in the ultra-relativistic regime, in our final analysis we can simply rescale the shock properties and radio emission to the physical values of $\dot{E}_{\rm f}$ and $\Gamma_{\rm f}$ ($\S\ref{sec:discussion}$).

The medium exterior to the binary, into which the inspiral wind enters at the start of the simulation ($t=0$), is taken to be that of a steady wind of power $\dot{E}_{0}$, mass-loss rate $\dot{M}_{\rm 0}$ and Lorentz factor $\Gamma_{0}$. Physically, this wind could represent the ordinary spin-down powered pulsar wind prior to the binary orbit entering the light cylinder, albeit for larger values of $a_{0}\approx a_{\rm bin}$ (Eq. 5) than adopted in our numerical simulations. We have checked that the pre-existing wind does not affect the shock interaction during the times we simulate. All of the wind material is ejected with a low thermal pressure ($p_{\rm gas}/\rho c^{2}=10^{-3}$, where $p_{\rm gas}$ is the gas pressure as defined in the upcoming section).

After the binary wind is shut-off at $t_{\rm m,0}-t_{\rm f}$, we continue to run the simulations to later times, $t_{\rm sim}\gg\Gamma_{\rm f}^{2}t_{\rm f}(a_{\rm 0}/a_{\rm f})^{4/3}$, to capture the interaction between the fast shell released at the end of the inspiral with the slowest material released at the beginning (Table 2 compiles the simulation run time). For numerical stability, the same small residual wind power $\dot{E}_{0}$ and mass-loss rate $\dot{M}_{0}$ of the pre-inspiral wind are also injected at late times after the merger $\gg t_{\rm m,0}$; however, this slow late-time wind does not immediately catch up to the fast shell and has no appreciable impact on the forward shock properties at larger radii of greatest interest.

Our numerical simulations are performed using Mara3, an open-source higher-order Godunov code first described in Zrake & MacFadyen (2012). The code has been extended with an Eulerian-Lagrangian moving mesh scheme that evolves the mesh vertices along with the outflowing plasma, similar to the scheme described in Duffell & MacFadyen (2013). This approach automatically concentrates numerical resolution in regions of high compression naturally arising around the forward shock. As mentioned above, we neglect the dynamical importance of the magnetic fields, assuming that the wind converts its Poynting flux to bulk kinetic energy at smaller radii than we simulate.

The code solves the special relativistic hydrodynamical equations in flux-conservative form,

where the vector of conserved quantities $\mathbf{U}=(D,S,\tau)^{T}$ consists of the lab-frame mass density $D=\Gamma\rho$, radial momentum density $S=\rho h\Gamma^{2}v/c^{2}$, and energy density (excluding rest-mass) $\tau=\rho h\Gamma^{2}-p_{\rm gas}-Dc^{2}$. $\Gamma=(1-v^{2}/c^{2})^{-1}$ is the gas Lorentz factor, and $h=c^{2}+\varepsilon+p_{\rm gas}/\rho$ is specific enthalpy. $\varepsilon$, $p_{\rm gas}$, $\rho$ denote the specific internal energy, gas pressure, and proper mass density, respectively. We adopt a gamma-law equation of state $p_{\rm gas}=(\gamma_{\rm ad}-1)\varepsilon\rho c^{2}$, with adiabatic index $\gamma_{\rm ad}=4/3$. In making this choice, we are neglecting the effects of radiative cooling of the post shock, i.e. we take the evolution to be adiabatic. The validity of this assumption will depend on the detailed composition of the binary wind (e.g., the presence of ions, which cool less effectively than pairs) and the synchrotron cooling timescale of these particles relative to the expansion rate. A more thorough exploration of the dynamical effects of cooling is left to future work (see also a discussion in §5.3).

The vector of fluxes is given by

and the geometrical source term is $\mathbf{S}_{\rm geom}=(0,0,2p_{\rm gas}/r)^{T}$.

Eq. (41) is discretized in terms of the extrinsic conserved quantities, $\mathbf{Q}_{i}=\int_{i}\mathbf{U}{\rm d}V$ integrated over each control volume. The control volumes are spherical shells whose inner and outer boundaries are moved to enforce a zero-mass-flux condition. The discretized solution is updated in time using the method of lines and an explicit second-order Runge-Kutta integration, whose base scheme is given by

Here $\mathbf{\hat{F}}_{i\pm 1/2}$ are the Godunov fluxes through the inner and outer cell boundaries, $A_{i+1/2}=4\pi r_{i+1/2}^{2}$ is the area of the cell boundary at $i+1/2$, $V_{i}$ is the cell volume, and $\mathbf{S}_{i,\rm geom}$ is the geometrical source term sampled half way between $r_{i+1/2}$ and $r_{i-1/2}$. The Godunov fluxes are given by $\mathbf{\hat{F}}=\mathbf{F}^{*}-\mathbf{U}^{*}v^{*}$, where $v^{*}$ is the speed of the contact discontinuity, and $\mathbf{U}^{*}$ and $\mathbf{F}^{*}$ are the intrinsic conserved quantities and associated fluxes at the $i+1/2$ cell interface. The starred quantities are obtained by sampling the solution (which is self-similar in $r/t$) of the Riemann problem $(\mathbf{U}_{i},\mathbf{U}_{i+1})$ at $r/t=v^{*}$. This solution is approximated using the Harten–Lax–van Leer contact (HLLC) solver of Mignone & Bodo (2006). The cell interfaces are advanced in time according to the same Runge-Kutta time-integration scheme, where

New cells are inserted at the radial inner boundary $r=r_{\rm in}$ by adding an additional face at $r_{\rm in}$ when the innermost face has advanced past $r_{\rm in}+\delta r$, where $\delta r$ is chosen based on the desired grid resolution. The fluid state of the inserted zone is set according to the wind inner boundary condition described in Eqs. (36, 39). Due to the high compression factor arising at high-$\Gamma$ shocks, cells can become very narrow, requiring a prohibitively short time step to satisfy the Courant Friedrichs-Lewy condition. Cells that are compressed below a threshold value $\delta r_{\rm min}$ are joined with an adjacent cell in a conservative fashion, by erasing the interface between the two cells and combining their extrinsic conserved quantities. To minimize the size disparity of adjacent zones, the removed interface is always the one between the too-small cell, and the smaller of its two neighbors. We also employ a conservative cell-splitting technique to keep from under-resolving rarified regions where the cells grow larger than $\delta r_{\rm max}$. The parameters $\delta r_{\rm min}$ and $\delta r_{\rm max}$ are typically kept within an order of magnitude of one another.

The outer edge of the simulation domain is placed at $L_{\rm box}/c=10^{3}$ s ($L_{\rm box}\approx 3\times 10^{13}\,\mathrm{cm}\gtrsim 10^{7}a_{\rm f}$) and the initial grid has a resolution of 512 cells per decade. The 1D nature of our special relativistic hydrodynamical simulations enables us to study the evolution of the ejecta at high resolution for longer durations—from the initial acceleration until late time deceleration phase. The duration of each simulated model, $t_{\rm sim}$, are listed in Table 2.

Two-dimensional simulations of decelerating relativistic ejecta by Duffell & MacFadyen (2013) have shown that the onset of Rayleigh-Taylor (RT) instability can disrupt the contact discontinuity, and generate turbulence in the shocked gas. The mixing causes the shocked gas to acquire more entropy and internal energy. This results in the reverse shock being pushed to smaller radii, and delays its observed emission relative to equivalent 1D models which are artificially stabilized against RT. We do not expect this phenomenon to change our overall results, as the radio emission in our model primarily arises from the forward shock (§4.1.2), which is generally unaffected by the development of RT.