© The Authors 2025

1. Introduction

Kilonovae (KNe) are transient events formed after the merger of binary neutron stars (Li & Paczyński 1998; Metzger et al. 2010; Rosswog 2015; Tanaka 2016; Fernández & Metzger 2016; Hotokezaka et al. 2018; Metzger 2019). Their emission initially peaks in the UV-optical-IR range and arises from heated ejecta material. The heat is generated from radioactive decay of heavy elements produced through r-process nucleosynthesis. The ejecta mass is expected to be in the range of 0.01 − 0.1 M_⊙, with a total kinetic energy of ∼10⁵¹ erg. If the merger results in a long-lived rapidly spinning and highly magnetized neutron star (millisecond magnetar), then a large amount of energy (typically ∼5 × 10⁵² erg up to ∼2 × 10⁵³ erg – limited primarily by the maximum possible rotation that the neutron star can attain before being torn apart, which is known as the mass-shedding limit; see, e.g. Giacomazzo & Perna 2013; Beniamini & Lu 2021) can be injected into the outflow. This can significantly enhance the brightness of the resulting KN (Yu et al. 2013; Metzger & Piro 2014; Kisaka et al. 2016; Ai et al. 2022; Wang et al. 2024).

In association with the gravitational-wave-detected neutron star merger, GW170817 (Abbott et al. 2017), a KN with ejecta of much lower kinetic energy, of the order of ∼10⁵¹ erg, has been observed with optical follow-up studies Coulter et al. (2017), Soares-Santos et al. (2017), Arcavi et al. (2017), Lipunov et al. (2017), Valenti et al. (2017), Troja et al. (2017), Kilpatrick et al. (2017), Smartt et al. (2017), Drout et al. (2017), Evans et al. (2017), suggesting a long-lived magnetar was not formed in this event. While supermassive neutron star (NS) formation might have been thought to be common based on considerations of total mass alone (e.g. Margalit & Metzger 2019), the occurrence rate of magnetar-boosted KNe (and the underlying long-lived central engines) is in practice strongly constrained by inferences based on observations. Relevant examples are (i) the analysis in Beniamini & Lu (2021), which is based on the typical durations and ejected energies of KN central engines as compared with the plateau-like features in short gamma-ray burst (sGRB) gamma-ray and X-ray observations; (ii) the long-term radio follow-up of low-redshift sGRBs (e.g. Horesh et al. 2016; Schroeder et al. 2020; Ricci et al. 2021), which puts limits on the kinetic energy of the outflow that are well below values expected for magnetar-boosted KNe (excluding the maximum ejecta energy of ≳5 × 10⁵² erg for M_ej ≲ 0.05 M_⊙, and constraining the fraction of short GRBs that produce stable remnants to ≲0.5); and (iii) the lack of bright transients in optical all-sky surveys such as Zwicky Transient Facility (ZTF; Bellm et al. 2019), which would naturally arise from magnetar-boosted KNe (Wang et al. 2024). The fraction of magnetar-boosted KNe is important for constraining the NS equation of state, because the latter is a major factor in determining the outcome of a NS merger (Shibata & Hotokezaka 2019).

While the optical transient from a KN substantially fades after a week or two, the ejecta continues to expand into the surrounding medium. As a result, a forward shock is driven into the surrounding material. This shock enhances the magnetic fields and accelerates electrons in the shocked medium, which then radiate their energy via synchrotron and synchrotron self-Compton. This can give rise to emission that is sometimes known as the ‘kilonova afterglow’ and can observationally manifest as approximately decades-long radio flares (Nakar & Piran 2011; Takami et al. 2014; Margalit & Piran 2015; Hotokezaka & Piran 2015; Kathirgamaraju et al. 2019). In the case of GW170817, the KN afterglow has not yet been observed even 7 years after the merger (e.g. Hajela et al. 2019; Balasubramanian et al. 2022).

Since the KN ejecta blast waves are expected to be mildly relativistic at most, the physics is essentially the same as in supernova remnants (SNRs). In analogy with the latter, this phase might also be dubbed the kilonova remnant (KNR; Beniamini & Lu 2021). Indeed, at times later than the ‘deceleration peak’ (also known as the Sedov time), the evolution depends only on the energy of the blastwave, and as such –considering the similarity between the energies of supernovae (SNe) and KNe–, it becomes challenging to distinguish between the two remnant types. The problem is compounded by the fact that the occurrence rate of SNe is about 300 times greater than that of KNe, and so SNRs are the major contaminant preventing potential identification of KNRs. Nevertheless, we show in the present work that there are important spectral, temporal, and environmental properties that can be used to distinguish KNRs from SNRs in observations.

In this paper, we perform an analytic study of the expected observable properties of KNRs and highlight their main differences from SNRs. In addition, we also focus on the class of highly energetic magnetar-boosted KN remnants (EKNRs), which, though expected to be rare, can be seen to much greater distances; therefore even a null detection of KNRs could strongly constrain their intrinsic occurrence rate. This paper is organised as follows. We start by presenting the dynamics of the ejecta of KNe in Sect. 2. We use Newtonian dynamics, which is relevant for our choice of fiducial parameters. However, we add a brief section (Appendix A) with fully relativistic expressions. In Sect. 3, we obtain expressions for the flux from emission originating from non-thermal electrons. We describe the combined spectrum due to synchrotron, bremsstrahlung, and synchrotron self-Compton emission and emphasise the differences between SNRs and KNRs. We predict the number of expected KNRs in radio surveys in Sect. 4 and show that we are likely to detect a candidate object with current (or future) facilities such as VLASS (Very Large Area Sky Survey) and VAST (Variables and Slow Transients). We are even more likely to detect an extreme KN, if they exist, due to the strong dependency of the observed flux on ejecta energy. As of yet, no such candidate object has been reported in the literature, which allows us to put moderate constraints on their abundance. In Sect. 5, we discuss how to extract the Sedov time from its observed light curve. For EKNRs, the Sedov time is of the order of 10–20 years, and therefore there is a good chance to obtain the Sedov time of extreme KNs with dedicated observations over a period of a few years. We also discuss how information about the host galaxies can be used to break degeneracies between KNRs and SNRs. We conclude in Sect. 6.

2. Ejecta dynamics

We compute the evolution of an ejecta (associated with either a KN or a SN) assuming it expands into a uniform surrounding medium in a spherically symmetric fashion. We assume Newtonian physics for the calculations shown here (see Appendix A for a fully relativistic treatment). The ejecta’s kinetic energy and velocity are connected via

$$\begin{array}{c}\hfill E=\frac{1}{2}(M_{\mathrm{ej}}+M_{\mathrm{ISM}})v_{\mathrm{ej}}^2,\end{array}$$(1)

where M_ej is the ejecta mass and M_ISM is the mass of displaced surrounding gas with $M_{\mathrm{ISM}}=\frac{4\pi }{3}{m}_{\mathrm{p}}{n}_{\mathrm{ISM}}{R}_{\mathrm{ej}}^3$, m_p is the mass of proton, n_ISM is the density of the surrounding gas and R_ej is the size of expanding ejecta. The kinetic and the internal energy of the ejecta can be comparable. Therefore, if we were to use the total energy of ejecta in Eq. (1), it would roughly cancel the factor 1/2 on the right-hand side. We assume a uniform density profile without any spatial gradient. We normalise the physical parameters by defining, $E_{51}=\frac{E}{{10}^{51}\mathrm{erg}},M_{\mathrm{ej},\odot }=\frac{M_{\mathrm{ej}}}{M_{\odot }}$, and $n_0=\frac{n_{\mathrm{ISM}}}{1{\mathrm{cm}}^{-3}}$. During the initial phase, the ejecta expands freely into the medium (we call this the ‘coasting phase’ regime). This lasts until M_ej ∼ M_ISM, when the ejecta starts decelerating and the evolution enters the Sedov-Taylor (ST) phase. The deceleration, or ST radius, R_ST can be obtained by solving

$$\begin{array}{c}\hfill \frac{E}{2}=\frac{4\pi }{3}{n}_{\mathrm{ISM}}{R}_{\mathrm{ST}}^3{m}_{\mathrm{p}}{\mathrm{c}}^2\frac{{\beta }_{\mathrm{ej}}^2}{2},\end{array}$$(2)

where v_ej = β_ejc. This simplifies to $R_{\mathrm{ST}}=0.52n_0^{-1/3}{M}_{\mathrm{ej},\odot }^{1/3}\times {10}^{19}\mathrm{cm}$. The corresponding timescale is given by

$$\begin{array}{c}\hfill t_{\mathrm{ST}}=165n_0^{-1/3}{E}_{51}^{-1/2}{M}_{\mathrm{ej},\odot }^{5/6}\mathrm{yr}.\end{array}$$(3)

The velocity in each regime is

$$\begin{array}{cc}\hfill {\beta }_{\mathrm{ej}}& =\frac{v_{\mathrm{ej}}}{\mathrm{c}}=0.03{\left(\frac{E_{51}}{M_{\mathrm{ej},\odot }}\right)}^{1/2}(\mathrm{coasting}\mathrm{phase})\hfill \\ \hfill {\beta }_{\mathrm{ej}}& =0.52E_{51}^{1/5}{n}_0^{-1/5}{t}_{\mathrm{yr}}^{-3/5}(\mathrm{ST}\mathrm{phase}).\hfill \end{array}$$(4)

The radius of the ejecta is then R_ej ≈ β_ejct in the coasting phase and $R_{\mathrm{ej}}\approx \frac{5}{2}{\beta }_{\mathrm{ej}}ct$ in the ST phase. We consider three classes of objects with fiducial physical parameters as stated in parentheses: a SNR (E₅₁ = 1, M_ej, ⊙ = 1, n₀ = 1), a KNR (E₅₁ = 1, M_ej, ⊙ = 0.1, n₀ = 0.1), and an EKNR (E₅₁ = 10, M_ej, ⊙ = 0.1, n₀ = 0.1). We note that this implicitly assumes that, compared to SNRs, both KNRs and EKNRs occur in environments that exhibit less star formation and are therefore less dense. While this is true on average, it means that KNRs and EKNRs appear less bright than they would in a similar environment to a SNR, which is a topic we return to in Sect. 5.3.

In these calculations, we have assumed that the external medium is undisturbed until the ejecta collides with it. However, relativistic jets from merging binaries can disturb the external medium before the ejecta. In that case, the early phase of the KNR emission could be suppressed due to the blast-wave passing through a region that was evacuated by the passage of the jet (Margalit & Piran 2020). The timescale of this suppression is related to t_ST by $t_{\mathrm{col}}=2.04{\left(\frac{E_{\mathrm{j}}}{E}\right)}^{1/3}{t}_{\mathrm{ST}}$. For a jet with an isotropic equivalent energy of 10⁵¹ erg (Duran & Giannios 2015) and an opening angle of θ_j ∼ 0.2 radians, the energy of the jet $E_{\mathrm{j}}=\frac{{\theta }_j^2}{2}{E}_{\mathrm{iso}}\sim {10}^{49}$ erg. For this case, we find that the suppression timescale is significantly shorter that the Sedov timescale. This means that propagation through the evacuated region typically only changes the KNR flux at early times, when the KNR is generally very dim. This conclusion is further enhanced in the case of EKNRs.

3. Synchrotron, synchrotron self-Compton, and Bremsstrahlung emission

The expanding ejecta shocks the surrounding medium, compressing it and accelerating the electrons. We assume that a fraction ϵ_e of the shock energy goes into accelerating electrons and a fraction ξ_e of electrons are accelerated to relativistic energies.

The relativistic electron Lorentz factor distribution is taken to be of the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_e}{\mathrm{d}{\gamma }_e}({\gamma }_e)=\frac{N_0}{{\gamma }_m}{\left(\frac{{\gamma }_e}{{\gamma }_m}\right)}^{-p},\end{array}$$(5)

where γ_m is the minimal Lorentz factor. Integrating over γ_e and equating the electron number density to n_ISM, N₀ = (p − 1)n₀. The expression for γ_m in the bulk comoving frame is given by

$$\begin{array}{c}\hfill {\gamma }_m-1=\frac{p-2}{p-1}\frac{{ϵ}_e}{{\xi }_e}\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}}(\mathrm{\Gamma }-1),\end{array}$$(6)

where it is assumed that p > 2 with a fiducial value of p = 2.5 and Γ is the Lorentz factor of the shocked fluid. In the non-relativistic regime of expanding ejecta and defining $\overline{{ϵ}_e}={ϵ}_e\frac{p-2}{p-1}$, $\overline{{ϵ}_e}=0.1{\overline{ϵ}}_{e,-1}$, the expression simplifies to, ${\gamma }_m-1=91.85\frac{{\overline{ϵ}}_{e,-1}}{{\xi }_e}{\beta }_{\mathrm{ej}}^2$, or

$$\begin{array}{cc}\hfill {\gamma }_m-1& =8.2\times {10}^{-2}\frac{{\overline{ϵ}}_{e,-1}}{{\xi }_e}\left(\frac{E_{51}}{M_{\mathrm{ej},\odot }}\right)(\mathrm{coasting}\mathrm{phase})\hfill \\ \hfill {\gamma }_m-1& =45.9\frac{{\overline{ϵ}}_{e,-1}}{{\xi }_e}{E}_{51}^{2/5}{n}_0^{-2/5}{t}_{\mathrm{yr}}^{-6/5}(\mathrm{ST}\mathrm{phase}).\hfill \end{array}$$(7)

For a fixed ξ_e and at sufficiently late times, Eq. (7) would lead to γ_m − 1 ≪ 1, which will result in electrons being non-relativistic, at which point they cannot emit synchrotron photons at ν≳ GHz. The reason for this is straightforward: at such times, there is simply not enough energy to share between those electrons, while still allowing all of them to maintain relativistic velocities. This phase is known as the deep Newtonian phase. A natural resolution, which is consistent with giant flare outflows and GRB afterglows (Granot et al. 2006; Sironi & Giannios 2013), is to take γ_m to be the maximum between the value given in Eq. (7) and $\sqrt{2}$. The fraction of relativistic electrons then needs to be modified according to ξ_e = ξ_e, 0min[1, (β_ej/β_dn)²] (Beniamini et al. 2022). In the analysis below, we assume ξ_e, 0 = 1. The expression for β_dn is then ${\beta }_{\mathrm{dn}}={\left(2^{3/2}\frac{p-1}{p-2}\frac{10}{{\overline{ϵ}}_{e,-1}}\frac{m_{\mathrm{e}}}{m_{\mathrm{p}}}\right)}^{1/2}$. Using p = 2.5 and ${\overline{ϵ}}_{e,-1}=1$, β_dn = 0.21. For our fiducial class of objects, we find that only the EKNRs transition to the deep Newtonian regime, while SNRs and KNRs are already in this regime in the coasting phase. As the ejecta transitions to the deep Newtonian regime in the ST era, requiring β_ej = β_dn in Eq. (4), we obtain the transition timescale as,

$$\begin{array}{c}\hfill t_{\mathrm{dn}}=7.43{\left({\overline{ϵ}}_{e,-1}^{1/2}{E}_{51}^{1/5}{n}_0^{-1/5}\right)}^{5/3}\mathrm{yr}.\end{array}$$(8)

For E₅₁ = 10, n₀ = 0.1, and ${\overline{ϵ}}_{e,-1}=1$, t_dn ≈ 30 yr. The expression for ξ_e is given by

$$\begin{array}{cc}\hfill {\xi }_e& \sim 2.25\times {10}^{-2}\frac{E_{51}}{M_{\mathrm{ej},\odot }}{\overline{ϵ}}_{e,-1}(\mathrm{coastingphase})\hfill \\ \hfill {\xi }_e& \sim 6.3E_{51}^{2/5}{\overline{ϵ}}_{e,-1}{n}_0^{-2/5}{t}_{\mathrm{yr}}^{-6/5}(\mathrm{ST}\mathrm{phase}).\hfill \end{array}$$(9)

3.1. Synchrotron emission

The energetic electrons in the presence of magnetic field emit synchrotron emission (Chevalier 1982, 1998; Frail et al. 2000). Assuming a fraction ϵ_B of shocked energy goes into amplifying the magnetic field, the comoving magnetic energy density in the non-relativistic regime is $\frac{B^2}{8\pi }=2m_{\mathrm{p}}{\beta }_{\mathrm{ej}}^2{\mathrm{c}}^2{n}_{\mathrm{ISM}}{ϵ}_B$. Substituting the expression for β_ej (Eq. (4)), the expression for the magnetic field becomes

$$\begin{array}{cc}\hfill B& =9.1\times {10}^{-4}{n}_0^{1/2}{ϵ}_{B,-2}^{1/2}{E}_{51}^{1/2}{M}_{\mathrm{ej},\odot }^{-1/2}\mathrm{G},(\mathrm{coasting}\mathrm{phase})\hfill \\ \hfill B& =0.014{ϵ}_{B,-2}^{1/2}{E}_{51}^{1/5}{n}_0^{3/10}{t}_{\mathrm{yr}}^{-3/5}\mathrm{G}(\mathrm{ST}\mathrm{phase}),\hfill \end{array}$$(10)

where ϵ_B ≡ 0.01ϵ_B, −2. The characteristic synchrotron frequency ν_m is given by ν_m = γ_m²ν_B, where ${\nu }_B=\frac{\mathrm{e}B}{2\pi m_{\mathrm{e}}\mathrm{c}}=2.8\times {10}^{6}B$ Hz is the cyclotron frequency. For our fiducial cases of SNR and KNR, the solutions are in the deep Newtonian regime, while for EKNR, γ_m ∼ 10 in the coasting phase (Eq. (7)) and the deep Newtonian regime is typically obtained in the ST phase. In the deep Newtonian regime, the expression for ν_m is given by

$$\begin{array}{cc}\hfill {\nu }_m& \sim {10}^4{n}_0^{1/2}{ϵ}_{B,-2}^{1/2}{E}_{51}^{1/2}{M}_{\mathrm{ej},\odot }^{-1/2}\mathrm{Hz},(\mathrm{coasting}\mathrm{phase})\hfill \\ \hfill {\nu }_m& \sim 1.6\times {10}^5{ϵ}_{B,-2}^{1/2}{E}_{51}^{1/5}{n}_0^{3/10}{t}_{\mathrm{yr}}^{-3/5}\mathrm{Hz}(\mathrm{ST}\mathrm{phase}),\hfill \end{array}$$(11)

for slow-cooling electrons, or ν_m < ν_c, sync (Sari et al. 1998; Granot & Sari 2002), where ν_c, sync is the cooling frequency due to synchrotron losses alone (Appendix B). The power in the comoving frame due to one electron at the frequency ν_m is $P_{{\nu }_m}=3.8\times {10}^{-22}\left(\frac{B}{1\mathrm{G}}\right)\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}$. The number of shocked electrons (which is a Lorentz invariant) within a dynamical time is $\sim \frac{4\pi }{3}{R}^3{\xi }_e{n}_{\mathrm{ISM}}$. Therefore, we multiply this number with the power per electron to obtain the maximum optically thin synchrotron luminosity, $L_{{\nu }_{m,0}}\sim {\xi }_e{n}_{\mathrm{ISM}}\left(\frac{4\pi }{3}\right)R^3{P}_{{\nu }_{m,0}}$. All together,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L_{{\nu }_{m,0}}& \sim 9\times {10}^{23}{\overline{ϵ}}_{e,-1}{ϵ}_{B,-2}^{1/2}{E}_{51}^3{M}_{\mathrm{ej},\odot }^{-3}{n}_0^{3/2}{t}_{\mathrm{yr}}^3{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}(\mathrm{coasting}\mathrm{phase}),\hfill \\ \hfill L_{{\nu }_{m,0}}& \sim 2\times {10}^{32}{\overline{ϵ}}_{e,-1}{ϵ}_{B,-2}^{1/2}{E}_{51}^{6/5}{n}_0^{3/10}{t}_{\mathrm{yr}}^{-3/5}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}(\mathrm{ST}\mathrm{phase}).\hfill \end{array}\end{array}$$(12)

The equation above applies in the deep Newtonian regime. Furthermore, for the physical situation being discussed, we verified that ν_m < ν_sa < ν_c, sync, where ν_sa is the synchrotron self-absorption frequency (Appendix C). We caution that in this case, the luminosity at ν_m is not given by the expression given in Eq. (12). Nonetheless, as we are interested in luminosity or flux at ν > ν_sa > ν_m, we can still use these expressions to compute flux by extrapolation, as done below. We use the subscript ‘0’ in the expressions for L_νm, 0 and F_νm, 0 to clarify this difference. The spectral flux at a given observation frequency ν with ν_m < ν_sa < ν < ν_c, sync is given by $F_{\nu }=F_{{\nu }_{m,0}}{\left(\frac{\nu }{{\nu }_m}\right)}^{-(p-1)/2}$, which after substitution becomes

$$\begin{array}{c}\hfill \begin{array}{c}\hfill F_{\nu }=5.4\times {10}^{-7}\frac{{\overline{ϵ}}_{e,-1}{ϵ}_{B,-2}^{7/8}{E}_{51}^{27/8}{M}_{\mathrm{ej},\odot }^{-27/8}{n}_0^{15/8}{t}_{\mathrm{yr}}^3}{{\mathrm{D}}_{10}^2}{\left(\frac{\nu }{3\mathrm{GHz}}\right)}^{-3/4}\mathrm{mJy}\\ \hfill (\mathrm{coasting}\mathrm{phase}),\\ \hfill F_{\nu }=8.85\times {10}^2\frac{{\overline{ϵ}}_{e,-1}{ϵ}_{B,-2}^{7/8}{E}_{51}^{27/20}{n}_0^{21/40}{t}_{\mathrm{yr}}^{-21/20}}{{\mathrm{D}}_{10}^2}{\left(\frac{\nu }{3\mathrm{GHz}}\right)}^{-3/4}\mathrm{mJy}\\ \hfill (\mathrm{ST}\mathrm{phase}),\end{array}\end{array}$$(13)

where we have assumed Euclidean geometry, as appropriate for distances of D ≲ 10 Mpc with $D_{10}=\frac{D}{10\mathrm{Mpc}}$. We normalised the expression for flux with respect to the central frequency of VLASS (Gordon et al. 2023). In Fig. 1, we plot the observed flux for the three classes of objects considered in this work. The peak of the observed flux occurs at t_ST with a peak flux of ∼1 mJy (D = 10 Mpc) for both the SNR and KNR. EKNRs are significantly brighter and also have shorter t_ST. This will significantly impact the detection prospects for such extreme objects; we discuss this topic in subsequent sections. We have also included the corrections by Margalit & Piran (2020) in this figure, which are briefly discussed in Sect. 2. We note that the calculations discussed above are valid at t > t_col, while at t < t_col, there is an extremely sharp decrease in flux (∼t¹⁰) due to the decrease in external density caused by the previous propagation of the jet.

Fig. 1. Flux at 3 GHz as a function of age and ejecta size with D = 10 Mpc for an SNR (E51 = 1, Mej, ⊙ = 1, n0 = 1), a KNR (E51 = 1, Mej, ⊙ = 0.1, n0 = 0.1), and a EKNR (E51 = 10, Mej, ⊙ = 0.1, n0 = 0.1). We have chosen ${\overline{ϵ}}_{e,-1}=1$ and ϵB, −2 = 1. We have included the corrections by Margalit & Piran (2020) in the left panel of the figure to highlight the qualitative features of propagation through jet-evacuated material. We ignore this aspect in the rest of the calculation and the remaining figures. The red and black curves overlap in the right panel, which is a coincidence because the difference in typical density compensates for that in ejecta mass.

3.1.1. Flux at t_ST

As objects have the highest detection probability when they are brightest, it is useful to have the expression for the peak flux, which occurs at t_ST (Fig. 1). Using the expression for t_ST (Eq. (3)) in the coasting phase solution for the observed flux (Eq. (13)), we have

$$\begin{array}{c}\hfill F_{3\mathrm{GHz}}^{\mathrm{ST}}=1.1\times {10}^2{\xi }_e\frac{E_{51}^{7/8}{M}_{\mathrm{ej},\odot }^{1/8}{\left({ϵ}_{B,-2}{n}_0\right)}^{7/8}}{D_{10}^2}\mathrm{mJy},\end{array}$$(14)

where we we keep ξ_e explicit. Using the expression for ξ_e in the deep Newtonian regime (Eq. (9)),

$$\begin{array}{c}\hfill F_{3\mathrm{GHz}}^{\mathrm{ST}}=2.4\frac{{\overline{ϵ}}_{e,-1}{E}_{51}^{15/8}{\left(\frac{{ϵ}_{B,-2}{n}_0}{M_{\mathrm{ej},\odot }}\right)}^{7/8}}{D_{10}^2}\mathrm{mJy}.\end{array}$$(15)

The peak flux is a strong function of ejecta energy, which explains the dramatic boost to observed flux in the case of a EKNR. The peak flux is a degenerate combination of ϵ_B, −2, n₀, and M_ej, ⊙. This is visible in Fig. 1 for the SNR and KNR cases, which both have n₀/M_ej, ⊙ = 1.

3.1.2. Peak flux at X-ray frequency

We also compute the flux at t_ST in the X-ray band. The luminosity at 1 keV (or ν = 2.4 × 10¹⁷ Hz) with ν > ν_c is given by $L_{\nu }=L_{{\nu }_{m,0}}{\left(\frac{{\nu }_c}{{\nu }_m}\right)}^{-(p-1)/2}{\left(\frac{\nu }{{\nu }_c}\right)}^{-p/2}$. At t = t_ST,

$$\begin{array}{c}\hfill \nu F_{\nu }{|}_{1\mathrm{keV}}^{\mathrm{ST}}=0.62\times {10}^{-13}{\xi }_e\frac{E_{51}^{5/8}{ϵ}_{B,-2}^{1/8}{n}_0^{11/24}{M}_{\mathrm{ej},\odot }^{1/24}}{D_{10\mathrm{Mpc}}^2(1+Y)}\mathrm{e}rgc{m}^{-2}{s}^{-1}.\end{array}$$(16)

Using the expression for ξ_e in the deep Newtonian regime (Eq. (9)),

$$\begin{array}{c}\hfill \nu F_{\nu }{|}_{1\mathrm{keV}}^{\mathrm{ST}}=1.4\times {10}^{-15}\frac{{\overline{ϵ}}_{e,-1}{E}_{51}^{13/8}{ϵ}_{B,-2}^{1/8}{n}_0^{11/24}{M}_{\mathrm{ej},\odot }^{-23/24}}{D_{10Mpc}^2(1+Y)}\mathrm{erg}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}.\end{array}$$(17)

In the equation above, and in the full results shown in Fig. 2, we account for inverse-Compton losses affecting ν_c (Sect. 3.2). For our fiducial cases, Y ≲ 1, and so synchrotron self-Compton losses are a minor correction. The flux limit for eROSITA in the soft X-ray band (0.5–2 keV) is ∼10⁻¹⁵ − 10⁻¹⁴ erg cm⁻² s⁻¹ (Merloni et al. 2012). The ongoing Einstein probe (Yuan et al. 2022) has a similar energy band with slightly worse flux limit, which depends on the length of exposure. Therefore, we expect to see objects in X-rays up to distances of a few megaparsecs, which can complement radio observations. If synchrotron is the dominant process in X-rays (as indeed expected; see Sect. 3.2), we expect a correlation between radio and X-ray flux, as both arise from the same population of electrons.

Fig. 2. νLν as a function of ν at tST for a SNR (black), a KNR (red), and a EKNR (blue).

3.1.3. Synchrotron emission cutoff

We compute the synchrotron cutoff by requiring that the diffusion length scale of electrons, while gyrating in the magnetic field, be less than the size of the remnant. The diffusion coefficient in the Bohm limit is given by D = cr_g/3, where r_g is the gyration radius, r_g = γm_ec²/eB, which can be simplified to $r_g=1.7\times {10}^3\frac{\gamma }{B}$ cm. Imposing the condition as stated above, we have Dt < β_ej²c²t², or $t>\frac{D}{{\beta }_{\mathrm{ej}}^2{\mathrm{c}}^2}$. We then use the synchrotron cooling timescale (Appendix B) in the expression to obtain

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{s},\mathrm{co}}=2\times {10}^8{\beta }_{\mathrm{ej}}{B}^{-1/2}.\end{array}$$(18)

A more detailed calculation for the cutoff Lorentz factor is given by Zirakashvili & Aharonian (2007), which in the Bohm regime leads to a pre-factor of 3.6 × 10⁷ instead of 2 × 10⁸. Using the value from Zirakashvili & Aharonian (2007), the corresponding cutoff frequency is ${\nu }_{\mathrm{s},\mathrm{co}}={\gamma }_{\mathrm{s},\mathrm{co}}^2{\nu }_B=3.6\times {10}^{21}{\beta }_{\mathrm{ej}}^2\mathrm{Hz}$. Plugging in β_ej,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill h{\nu }_{\mathrm{s},\mathrm{co}}& \sim {10}^4\left(\frac{E_{51}}{M_{\mathrm{ej},\odot }}\right)\mathrm{eV}(\mathrm{coasting}\mathrm{phase}),\hfill \\ \hfill h{\nu }_{\mathrm{s},\mathrm{co}}& \sim 4\times {10}^6{E}_{51}^{2/5}{n}_0^{-2/5}{t}_{\mathrm{yr}}^{-6/5}\mathrm{eV}(\mathrm{STphase}).\hfill \end{array}\end{array}$$(19)

As β_ej is constant in the coasting phase, the corresponding expression gives a rough estimate at t_ST (at later times ν_s, co only declines). For SNRs, the expected cutoff is of the order of ∼10 keV, which becomes ∼100 keV – 1 MeV for our choice of KN parameters. Therefore, this feature, along with others (which we discuss later), can be used to separate SNRs from KNRs.

3.2. Synchrotron self-Compton emission

In addition to synchrotron, electrons also cool via Compton losses by the synchrotron photons, known as synchrotron self-Compton (SSC) loss. Therefore, Eq. (B.1) is modified as

$$\begin{array}{c}\hfill {\gamma }_c{m}_{\mathrm{e}}{\mathrm{c}}^2=\frac{4}{3}{\sigma }_{\mathrm{T}}{\gamma }_c^2(\frac{B^2}{8\pi }+U_{\mathrm{sync}})\mathrm{c}t,\end{array}$$(20)

where U_sync is the energy density of the synchrotron photon spectrum. As a result, the cooling Lorentz factor is modified as ${\gamma }_c\to \frac{{\gamma }_{\mathrm{c},\mathrm{sync}}}{1+Y}$, where Y is the Compton parameter whose expression is given by Beniamini et al. (2015), $Y=\frac{{ϵ}_e}{{ϵ}_B(3-p)(1+Y)}{\left(\frac{{\nu }_m}{{\nu }_c}\right)}^{(p-2)/2}$, where ν_c is the cooling frequency with the Compton correction included, that is, ν_c = ν_c, sync/(1 + Y)², where ν_c, sync is the cooling frequency without SSC losses. Plugging in the expressions with p = 2.5, we have

$$\begin{array}{c}\hfill Y{(1+Y)}^{0.5}=0.2{\overline{ϵ}}_{e,-1}{ϵ}_{B,-2}^{-3/4}{E}_{51}^{1/4}{n}_0^{1/3}{M}_{\mathrm{ej},\odot }^{-1/12}.\end{array}$$(21)

Using fiducial parameters, we obtain Y(1 + Y)^0.5 ∼ 0.18, 0.24, and 0.43 for SNRs, KNRs, and EKNRs, respectively. Therefore, Compton cooling has a minor effect on the scenarios that we consider. However, if ϵ_B, −2 is much smaller than 1, then this effect is moderately important. Nevertheless, we include this effect in our calculation of the synchrotron spectrum as shown in Fig. 2. The Compton-scattered frequency corresponding to ν_c is given by ν_c,  IC ∼ γ_c²ν_c ∼ γ_c⁴ν_B. The value of γ_c at t_ST is of the order of 10⁵. Therefore, ν_c, IC ∼ 10²⁴ Hz¹. We solve for the SSC flux as a function of frequency, as explained in Appendix D.

In Fig. 2, we plot νL_ν for synchrotron, SSC, and bremsstrahlung at t_ST for our fiducial cases. We find that bremsstrahlung (Appendix E) is subdominant compared to the other channels over the entire frequency domain. For SNRs, with larger n₀, bremsstrahlung may dominate over synchrotron in the X-ray band because it is highly sensitive to the value of the gas density. For our fiducial parameters of KNR, the synchrotron emission is further boosted compared to bremsstrahlung. Therefore, we expect the X-ray flux to be dominated by synchrotron, which may be a diagnostic tool to separate it from SNRs, especially with the synchrotron cut-off feature (Sect. 3.1.3). The SSC process completely dominates the emission in the gamma-ray band. For KNRs, the overall transition from synchrotron to SSC happens around 1-10 MeV. In this context, we briefly discuss γ-ray detectability with the Fermi gamma-ray burst monitor (Fermi-GBM). The sensitivity of this instrument is of the order of 0.5 cm⁻² s⁻¹ between 10 keV and 25 MeV. We note that this sensitivity is given in terms of photon flux, which is the number of photons per unit area per second. For the EKNR case in Fig. 2, the flux will be of the order of ∼10⁻⁷ cm⁻² s⁻¹ at 100 keV for an object located at D = 10 Mpc. Therefore, it would be difficult to detect features in the gamma-ray part of the spectrum. We also show the relevant energy range of ZTF in Fig. 2. Using its magnitude limit, we find that we can see objects up to a few megaparsecs. In subsequent sections, we concentrate on the radio part of spectrum and discuss how the light curves (Sect. 5.1) may help us distinguish KNRs from SNRs.

4. Number of detectable remnants

We now consider the detection prospects of KNRs with ongoing and future telescope surveys. The very large area sky survey (VLASS) is observing the entire sky north of −40° (33 885 deg²) in the 2 GHz ≲ ν ≲ 4 GHz band with an angular resolution of the order of 3 arcsec. The current sensitivity is of the order of 1 mJy (Gordon et al. 2021), which will eventually reach the level of 0.1 mJy. VLASS will cover its footprint over three epochs separated by 32 months. For the current flux threshold of ∼1 mJy, the expected distance up to which a KN can be seen is ∼10 Mpc (Fig. 1). For a generic flux threshold F_ths, the maximum distance is given by D^max = (L_3 GHz/4πF_ths)^1/2, which can be rewritten as $D_{10}^{\max }={\left(\frac{F_{3\mathrm{GHz},10}}{F_{\mathrm{ths}}}\right)}^{1/2}$, where F_{3 GHz, 10} is the flux at 10 Mpc. Objects of a given age t are detected at a rate of

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dt}}=\mathcal{R}(t)\frac{4\pi }{3}{D}^3,\end{array}$$(22)

where ℛ(t) is the rate of formation of KNe per volume per time, which is ∼300 Gpc⁻³ yr⁻¹ (Abbott et al. 2023). Substituting the numbers, we have, $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dt}}=1.2\times {10}^{-3}{(D_{10}^{\max })}^3{\mathrm{yr}}^{-1}$. For L_ν ∝ t^α, $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dt}}\propto t^{3\alpha /2}$ or equivalently $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogt}}\propto t^{1+3\alpha /2}$. While Eq. (22) is denoted in terms of the observer’s time, the luminosity is a function of the age of the system. In practice, we have many systems with different ages and we are averaging over their properties at a given observer time. Under our assumptions, this does not change our results. For the coasting phase regime, α = 3, while for the ST regime, α = −21/20 (Eq. (13)). Substituting, we have

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogt}}\propto t^{11/2}(\mathrm{coasting}\mathrm{phase}),\hfill \\ \hfill & \frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogt}}\propto t^{-1/2}(\mathrm{ST}\mathrm{phas}).\hfill \end{array}\end{array}$$(23)

For the fiducial KNR parameters, $D_{10}^{\max }=1.5$.

We expect the quantity $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogt}}$ to be dominant at the Sedov time. At t_ST, we can use Eq. (15) to obtain the general expression,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dt}}=0.004\left(\frac{\mathcal{R}(t)}{300{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}}\right){\overline{ϵ}}_{e,-1}^{\frac{3}{2}}{E}_{51}^{\frac{45}{16}}{\left(\frac{{ϵ}_{B,-2}{n}_0}{M_{\mathrm{ej},\odot }}\right)}^{\frac{21}{16}}{\left(\frac{F_{\mathrm{ths}}}{1\mathrm{mJy}}\right)}^{-\frac{3}{2}}{\mathrm{yr}}^{-1}.\end{array}$$(24)

For our fiducial KNR case, t_ST ∼ 50 yr or $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}\log \mathrm{t}}\sim 0.2$, meaning that there is a roughly 20% chance that such an object exists within the current VLASS dataset. In the future, with a reduction of the VLASS flux threshold to 0.1 mJy (Gordon et al. 2023), these numbers will be boosted by a factor of about 30. In Fig. 3, we plot $\frac{\mathrm{d}N}{\mathrm{dlog}t}$ for SNRs, KNRs, and EKNRs with our choice of fiducial parameters. We also consider a few cases with a distribution over density and ejecta mass of EKNRs, which we describe below. We take the SN occurrence rate to be ∼10⁵ Gpc⁻³ yr⁻¹ (Li et al. 2011), which is approximately 300 times higher than the KN rate. A consistency check for the number of detectable SNRs given our assumptions is provided in Appendix F.

Fig. 3. Abundance of SNR, KNR and EKNRs. Top: Number of systems observed by an all-sky radio survey per logarithmic unit of age, $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlog}t}$, for SNR, KNR, and EKNRs. Bottom: $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlog}L}$ for the same cases as in the upper panel. For the SN case, we show $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlog}L}$ for the coasting regime in grey in order to highlight its non-monotonic nature. We assume a minimum flux threshold of Fths = 1 mJy, an observation frequency of 3 GHz, and fiducial ejecta parameters as in Fig. 1. EKNRs evolve on a short-enough timescale for their evolution to be observable (as explored in more detail in Fig. 5).

While Eq. (23) is written in terms of time, for an object at a known distance, the luminosity is directly inferred from observations. Therefore, it is useful to convert the expression into one in terms of luminosity using the relation, $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}L}=\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}t}{\left(\frac{\mathrm{d}L}{\mathrm{d}t}\right)}^{-1}$. Using the dependence of luminosity on time, we have $L\left(\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}L}\right)=\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}t}\frac{t}{\alpha }$. In terms of luminosity,

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogL}}\propto L^{11/6}(\mathrm{coasting}\mathrm{phase}),\hfill \\ \hfill & \frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogL}}\propto L^{1/2}(\mathrm{ST}\mathrm{phase}).\hfill \end{array}\end{array}$$(25)

In the bottom panel of Fig. 3, we plot the same as in the upper panel but in terms of luminosity. We only consider the ST regime for this plot because $\frac{\mathrm{d}N}{\mathrm{dlog}L}$ is a non-monotonic function of L and the number of objects in the ST regime dominates over the number of sources in the coasting phase (see the thin grey line in the figure).

If we are able to spatially resolve the remnants, we can directly track their evolution in size as a function of time (see e.g. Barniol Duran et al. 2016). This provides additional information and can inform us as to whether the object is in the coasting or ST phase. In Fig. 1 (right panel), we plot the observed flux as a function of the size of the ejecta. We find that SNRs and KNRs are degenerate, with the Sedov radii of all three cases being similar. Nevertheless, observation of the ejecta evolution over time can help distinguish remnants from potential contaminants. The related observable is the surface brightness S_ν given by $S_{\nu }=\frac{L_{\nu }}{4\pi R_{\mathrm{ej}}^2}$. The observed distribution of S_ν is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{S}_{\nu }}{\mathrm{dt}}=\mathcal{R}(t)\frac{4\pi }{3}{\left(\frac{L_{\nu }}{4\pi F_{\mathrm{ths}}}\right)}^{3/2}\left(\frac{L_{\nu }}{4\pi R_{\mathrm{ej}}^2}\right).\end{array}$$(26)

As the size of a KNR varies as t and t^2/5 in blast-wave and ST-regime, respectively,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{\mathrm{d}{S}_{\nu }}{\mathrm{dlogt}}\propto t^{13/2}(\mathrm{coasting}\mathrm{phase}),\\ \hfill \frac{\mathrm{d}{S}_{\nu }}{\mathrm{dlogt}}\propto t^{-23/10}(\mathrm{ST}\mathrm{phase}).\end{array}\end{array}$$(27)

Given the flux sensitivity of a survey, we can see objects up to the maximum distance $D_{10}^{\max }$, as introduced above. The physical size of a resolvable object is given by $R_{\mathrm{rsl}}=4.5\times {10}^{20}\left(\frac{{\theta }_{\mathrm{rsl}}}{3\prime }\right)D_{10}^{\max }\mathrm{cm}$, where θ_rsl is the angular resolution of the survey, which we have normalised to 3 arcsec, for applicability to VLASS. As the Sedov radius of SNRs and KNRs are of the order of 10¹⁹ cm, in order to resolve a KNR, we have the relation

$$\begin{array}{c}\hfill \left(\frac{{\theta }_{\mathrm{rsl}}}{3″}\right)D_{10}^{\max }\lesssim 0.02.\end{array}$$(28)

We can rewrite this equation in terms of F_ths as

$$\begin{array}{c}\hfill \left(\frac{{\theta }_{\mathrm{rsl}}}{3″}\right)F_{\mathrm{ths}}^{-1/2}\lesssim 0.02F_{3\mathrm{GHz},10}^{-1/2},\end{array}$$(29)

where we can use Eq. (15) for the expression of F_{3 GHz, 10} at t_ST. Assuming a conservative value of θ_rsl = 1″, we can resolve KNRs up to ∼600 kpc (Eq. (28)). We find the number of objects within 600 kpc to be $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogt}}\sim {10}^{-5}$ at t_ST, even when assuming ℛ(t) to be ten times higher than the assumed value in Eq. (22). Therefore, we do not expect to find a spatially resolvable KNR with a survey like VLASS.

4.1. Extended distribution over density and ejecta mass

Up to this point, we have considered single values of the external density and ejecta mass for KNe. Here, we consider a distribution of external density, and we discuss the same for ejecta mass below. We assume the external density of EKNRs to be

$$\begin{array}{c}\hfill \frac{\mathrm{d}P}{\mathrm{d}{n}_0}=A{n}_0^{-a}.\end{array}$$(30)

The prefactor A is fixed by the normalisation ${\displaystyle {\int }_{n_{0,\min }}^{{\mathrm{n}}_{0,\max }}}\frac{\mathrm{d}P}{\mathrm{d}{n}_0}\mathrm{d}{n}_0=1$ or A = (∫_n0, min^n_0, maxn₀^−adn₀)⁻¹. The rate of detectable objects is

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{R}(t)\frac{4\pi }{3}{\displaystyle {\int }_{n_{0,\min }}^{{\mathrm{n}}_{0,\max }}}D{(n_0)}^3\frac{\mathrm{d}P}{\mathrm{d}{n}_0}\mathrm{d}{n}_0,\end{array}$$(31)

where D is a function of other ejecta parameters as well as the flux threshold, but we have suppressed them for brevity. We compute this quantity numerically to obtain the results shown in Fig. 4. Here, we give a simple estimate of how $\frac{\mathrm{d}N}{\mathrm{d}t}$ changes at t_ST depending on the value of a. As L_ν ∝ n₀^7/8 at t_ST (Eq. (15)) and the observable distance $D^{\max }\propto L_{\nu }^{1/2}$, the rate of detectable objects is

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t}\sim {\displaystyle \int n_0^{\frac{21}{16}-a}}\mathrm{d}{n}_0\sim n_0^{\frac{37}{16}-a}.\end{array}$$(32)

Fig. 4. Dependence of EKNR abundance on surrounding density and mass distribution. Top: $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlog}t}$ for EKNRs taking a distribution of external densities n0 with $n_{0,\max }=1{\mathrm{cm}}^{-3}$ and varying n0, min and a (Eq. (30)) shown by different line styles. In green, we consider a distribution over both density and ejecta mass but not over ejecta energy. For reference, we show our fiducial case of EKNR as in Fig. 3 in solid blue lines. Bottom: Areal density for the same cases as in the top panel. We also show the KNR case from Fig. 3, which is a factor of about 102 − 103 smaller than the EKNR cases. We have used the results of Anderson et al. 2020, Murphy et al. 2021, Mooley et al. 2022 as upper limits, because these authors do not detect a KNR.

As t_ST ∼ n₀^−1/3, the number of detected objects at t_ST is approximately N ∼ n₀^{95/48 − a}. We do not know the exact value of a at present. To motivate a guess for a, we assume n_0, min = 10⁻³, n_0, max = 1 cm⁻³ and require that the number of objects at n = 1 cm⁻³ is of the order of 1 percent of the population at n₀ = 10⁻³ (this is guided by observations of short GRB afterglows by O’Connor et al. 2020, as well as by the inferred gravitational wave delay time distribution between formation and merger of binary neutron star mergers; see Beniamini & Piran 2024). This rough estimate gives a = 5/3. Incorporating this value into Eq. (32), the number of detectable objects goes as N ∝ n₀^15/48. Therefore, a small number of objects with large external density can naturally dominate the number of detectable objects. The total number of detected objects depends on n_0, min, n_0, max, and a through the normalisation parameter A. We consider a few cases where we fix n_0, max = 1 and vary a and n_0, min such that the number of object at n = 1 cm⁻³ is of the order of 1 percent of the population at n_0, min. We plot these cases in Fig. 4. We see that the number of detected objects depends sensitively on our choice of the density distribution.

We also consider a distribution of ejecta masses with,

$$\begin{array}{c}\hfill \frac{\mathrm{d}P}{\mathrm{d}{M}_{\mathrm{ej}}}=B{M}_{\mathrm{ej}}^{-b},\end{array}$$(33)

where B is the normalisation parameter. We take the lower and upper limit of ejecta mass to be 0.01 and 0.1 M_⊙, respectively. Without proper knowledge of the value of b, we choose b = 1, which corresponds to equal logarithmic contributions. Similar steps as above lead to

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{R}(t)\frac{4\pi }{3}{\displaystyle {\int }_{{\mathrm{n}}_{0,\min }}^{{\mathrm{n}}_{0,\max }}}{\displaystyle {\int }_{{\mathrm{M}}_{\mathrm{ej},min}}^{{\mathrm{M}}_{\mathrm{ej},\max }}}D{(n_0,M_{\mathrm{ej}})}^3\frac{\mathrm{d}P}{\mathrm{d}{n}_0}\frac{\mathrm{d}P}{\mathrm{d}{M}_{\mathrm{ej}}}\mathrm{d}{n}_0\mathrm{d}{M}_{\mathrm{ej}},\end{array}$$(34)

We solve this expression numerically and plot it in Fig. 4. We note that for extreme KNe with M_ej, ⊙ = 0.01, the bulk velocity is mildly relativistic and we use the dynamical expressions provided in Appendix A. Still, some of the aspects can be understood from scaling relations obtained in the non-relativistic regime. While the flux at t_ST increases with decreasing M_ej (Eq. (15)), the Sedov time itself is also a strong function of M_ej (Eq. (3)). Therefore, the two effects compensate for each other in the quantity $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlog}t}$ and the ejecta mass does not have a huge effect on the number of detectable objects at t_ST. Reducing M_ej further has a much more drastic effect, with the ejecta moving at relativistic speeds, and causes a huge increase in the expected number of detected objects with extremely short t_ST.

From the above calculations, we find that we should be in a position to detect a KNR in the near future with a VLASS-like survey with an achievable goal of F_ths ∼ 0.1 mJy. Even more interestingly, if EKNRs constitute a significant population of all KNRs, we expect to find 10–100 such objects. If we do not find any in the survey, we should be able to put strong constraints on the population of EKNRs. As no such candidate objects are known in the literature to our knowledge, we take the conservative approach and obtain an upper limit on the abundance of EKNRs (Sect. 4.2). We discuss the potential signatures that can be used to differentiate these objects from SNRs in Sect. 5.

We also plot the areal density, or the number of objects per square degree, for the objects considered in Fig. 4. For the number of objects, we consider the peak of the distribution of $\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{dlogt}}$. The number of objects is proportional to (D^max)³, which is turn is proportional to F^−3/2 for an object of given luminosity. This gives rise to the qualitative behaviour seen in the bottom panel of Fig. 4.

4.2. Constraints on extreme KN abundance

In order to constrain the abundance of EKNRs, we show presently available constraints from transient surveys in the bottom panel of Fig. 4. These are dedicated surveys that do not report a single KNR. Therefore, we use their results to put upper limits on the abundance of EKNRs. Presently, the strongest constraints come from the observations of Murphy et al. (2021), but they strongly depend on the assumed density and mass distribution of EKNRs, as can be seen in Fig. 4. For our fiducial case (E₅₁ = 10, M_ej, ⊙ = 0.1, n₀ = 0.1), we find f_EKN ≲ 0.3, where f_EKN is the fraction of EKNRs with respect to the overall KNR population. As this constraint essentially depends on the value of $t\frac{{\mathrm{dN}}_{\mathrm{obs}}}{\mathrm{dt}}$ at t_ST, we derive a general expression using Eq. (24) and (3):

$$\begin{array}{c}\hfill f_{\mathrm{EKN}}\lesssim 0.3{\left(\frac{\mathcal{R}(t)}{{300\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}}\right)}^{-1}{\overline{ϵ}}_{e,-1}^{-\frac{3}{2}}{ϵ}_{B,-2}^{-\frac{21}{16}}{E}_{52}^{-\frac{37}{16}}{n}_{-1}^{-\frac{47}{48}}{M}_{\mathrm{ej},\odot ,-1}^{\frac{23}{48}}{\left(\frac{F_{\mathrm{ths}}}{1\mathrm{mJy}}\right)}^{\frac{3}{2}}.\end{array}$$(35)

4.3. Future constraints

In Fig. 4, we also plot future constraints expected from VAST wide (Murphy et al. 2021) and VLASS (Mooley et al. 2016; Gordon et al. 2021). The expected areal densities from these surveys are 4 × 10⁻⁷ at 2.5 mJy (1.4 GHz) and 10⁻⁵ at 1 mJy (3 GHz), respectively. These surveys are expected to find EKNRs if they constitute a significant fraction of the total KNR population. If no candidate objects are found in these surveys, we should be able to put strong constraints on the occurrence rate of EKNRs irrespective of their assumed density and mass distribution.

5. Detection prospect of kilonova remnants

5.1. Estimating t_ST from the observed light curve

We turn next to the detection prospects for KNRs from the evolution of their light curves. The observed flux F_3 GHz ∝ t^α, where α is a constant, while observing deep within either the coasting or ST phases. The proportionality constant sets the overall magnitude of flux and depends on ejecta parameters and distance. By taking the derivative of the logarithm of the flux, the constant term drops out and we have

$$\begin{array}{c}\hfill \frac{\mathrm{dlog}{F}_{3\mathrm{GHz}}}{\mathrm{d}t}=\frac{\alpha }{t}+\frac{\mathrm{d}\alpha }{\mathrm{d}t}\log \mathrm{t}.\end{array}$$(36)

At t ≪ t_ST or t ≫ t_ST, the second term drops out with α = 3 or −21/20, respectively. In Fig. 5, we plot the discretised version of $\frac{|\mathrm{dlog}(F_{3\mathrm{GHz}})|}{\mathrm{d}t}$, assuming observation intervals of Δt = 5 years for our fiducial cases. Far from their respective Sedov time, it is not possible to differentiate the three objects, because the light curves approach the same (featureless) power-law scaling. However, most of the objects are expected to be detected in the Sedov-Taylor regime, because the objects spend most of their time in that domain.

Fig. 5. Detection prospects for KNRs from their temporal evolution information. Top: $\frac{\mathrm{\Delta }F}{F\mathrm{\Delta }t}$ for our fiducial cases with Δt = 5 yr. $\frac{\mathrm{dlog}{F}_{3\mathrm{GHz}}}{\mathrm{d}t}$ approaches 3/t and ∼1/t in the blast wave and ST regimes, respectively. Bottom: $(\frac{\mathrm{\Delta }F}{F})\times (\frac{t_{\mathrm{ST}}}{\mathrm{\Delta }t})$ for the three cases with Δt = 5 yr. We note that tST is not available to us directly and has to be extracted from the light-curve data. As $(\frac{\mathrm{\Delta }F}{F})\times (\frac{t_{\mathrm{ST}}}{\mathrm{\Delta }t})$ is similar for all cases around the Sedov time, we can estimate tST from observed quantities, such as $\frac{\mathrm{\Delta }F}{F}$ and Δt (Eq. (37)). As an example, we show the case for an EKNR for which $\frac{\mathrm{\Delta }F}{F\mathrm{\Delta }t}\sim 0.02{\mathrm{yr}}^{-1}$ at around the time when it starts to enter deep into the Sedov regime (upper panel).

Assuming we detect an object around t_ST, we want to determine this timescale from the observed $\frac{\mathrm{\Delta }F}{F}$ at intervals of Δt. From Fig. 5, we see that, at around t_ST, $\frac{\mathrm{\Delta }F}{F}$ is higher for EKNRs, which is expected given that t_ST is smaller than for SNRs, and thus the flux changes faster. As the shape of the light curve is the same in all cases, in dimensionless units they must all be identical. As a visual demonstration, we plot $(\frac{\mathrm{\Delta }F}{F})\times (\frac{t_{\mathrm{ST}}}{\mathrm{\Delta }t})$ for our fiducial cases in the bottom panel of Fig. 5. We find that $(\frac{\mathrm{\Delta }F}{F})\times (\frac{t_{\mathrm{ST}}}{\mathrm{\Delta }t})$ line up horizontally around their respective t_ST. As an example, $(\frac{\mathrm{\Delta }F}{F})\times (\frac{t_{\mathrm{ST}}}{\mathrm{\Delta }t})\sim 0.5$ when the light curves are about to enter deep into the ST regime. Therefore, if we detect an object at this epoch, we can obtain an estimate for t_ST as,

$$\begin{array}{c}\hfill t_{\mathrm{ST}}\sim 0.5\mathrm{\Delta }t{\left(\frac{\mathrm{\Delta }F}{F}\right)}_{\mathrm{pk}}^{-1},\end{array}$$(37)

where ${\left(\frac{\mathrm{\Delta }F}{F}\right)}_{\mathrm{pk}}$ corresponds to the value close to the light curve peak. For EKNRs, ${\left(\frac{\mathrm{\Delta }F}{F\mathrm{\Delta }t}\right)}_{\mathrm{pk}}\sim 0.02{\mathrm{yr}}^{-1}$. By incorporating this value into Eq. 37, we have t_ST ∼ 25 years. We expect to detect a greater number of KNRs with higher external density (Sect. 4.1), which have even lower t_ST. If we find such an object at around its t_ST, 5–10 years of detailed observations with finer resolution in Δt could help to obtain a precise light curve, which would be helpful in distinguishing KNRs from SNRs and other contaminants. As an example, with finer resolution, we can reconstruct the feature that we see around t_ST. With coarse resolution, we may be unable to find this feature in the data as t_ST ∼ 20 yr for EKNRs.

We should point out that observing such an object can be challenging, as ΔF/F ∼ 0.1 over a 5 yr time period can be confused with spatial variation of the surrounding gas density, scintillation, or calibration issues, and so on. As an example, Mooley et al. (2016) used the criteria of ΔF/F ≳ 0.3 at S/N ≳ 5. We note that the radio afterglow of GW170817 does not show a deviation of more than 30 percent from a power-law spectrum Mooley et al. (2018). In addition, the calibration error can also be of the order of 5 percent (Ofek et al. 2011). Therefore, one needs high-S/N events so that the observational uncertainty is < 10 percent. However, the fractional change we are discussing here is smooth, whereas fluctuations due to scintillation or density inhomogeneities will likely be stochastic. Furthermore, large sources at close distances –that is, with R ≳ 10¹⁹ cm at D < 100 Mpc, as in the remnants discussed here– will be resolved by the scintillation screen (see Kumar et al. 2024 for details) and as such are very weakly affected by scintillation.

5.2. Chance coincidence of KNR with other radio sources

As we expect the candidate KNRs to be point sources at distances of ∼10 Mpc (see Sect. 4), these can be confused with other radio sources. We compute the probability of a chance coincidence of a candidate source with another radio object. The source count at S > 40 μJy is given by (Richards 2000),

$$\begin{array}{c}\hfill \frac{\mathrm{dn}}{\mathrm{d}S}=8.3S^{-2.4}{\mathrm{Jy}}^{-1}{\mathrm{Sr}}^{-1}.\end{array}$$(38)

The source count of VLASS at a threshold of ∼1 mJy (Gordon et al. 2021) matches very well with a six-degree polynomial fit (Hopkins et al. 2003), which differs somewhat from the simple power-law fit at S > 0.1 mJy. For simplicity, we use the expression above for our estimation. The probability of a chance coincidence is given by (Bloom et al. 2002; Eftekhari & Berger 2017)

$$\begin{array}{c}\hfill P_{\mathrm{ch}}=1-{\mathrm{e}}^{-\pi \sigma {\theta }^2}\approx \pi \sigma {\theta }^2,\end{array}$$(39)

where σ is the number of sources per solid angle and θ can be the angular size of a source such as a galaxy. From Eq. (38), we obtain n(S) as a function of S. For a full sky survey, such as VLASS, we obtain σ ∼ 2 × 10⁻⁶ arcsec⁻² at 1 mJy. In order to identify a source as a potential remnant candidate, we need to obtain its distance, and from that we should be able to identify the host galaxy. For a galaxy of ∼10 kpc in size at D ∼ 10 Mpc, the angular size is of the order of arcminutes. Putting everything together with appropriate normalisation, the probability of a chance coincidence is given by

$$\begin{array}{c}\hfill P_{\mathrm{ch}}\sim 0.02{\left(\frac{S}{1\mathrm{mJy}}\right)}^{-1.4}{\left(\frac{\theta }{\mathrm{arcmin}}\right)}^2.\end{array}$$(40)

For a galaxy of 10 kpc in size, we find that the probability of a chance coincidence can be of the order of a few percent. This might be an issue where the overall number of candidate objects is a few, but they may be distinguishable using spectral (Sect. 3) and light-curve information.

The EKNRs are expected to be much brighter with flux of ∼100 mJy (Fig. 1) at a distance of 10 Mpc. In that case, the probability of a chance coincidence is given by

$$\begin{array}{c}\hfill P_{\mathrm{ch}}\sim 3.2\times {10}^{-5}{\left(\frac{S}{100\mathrm{mJy}}\right)}^{-1.4}{\left(\frac{\theta }{\mathrm{arcmin}}\right)}^2.\end{array}$$(41)

In these expressions, we have considered all objects that are seen in the radio sky. However, we can use the flux variability information to reduce the chance coincidence even further. The authors in Becker et al. (2010), Bannister et al. (2011) and Croft et al. (2010) show that the number of variable sources on a timescale of 7–22 yr is of the order of a few percent. Using VLA and three epochs spanning 15 yr, Becker et al. (2010) show that the number of extragalactic sources with a variability of greater than 2 is less than 10 percent within the flux-density range of 1–100 mJy. Therefore, using the variability information we can increase the confidence of our detection even further.

5.3. Distinguishing between KNRs and SNRs with flux information

In the above sections, we discuss the importance of temporal measures such as the Sedov time in order to distinguish a KNR from a SNR. The measurement of flux at a given epoch may still be used for that purpose provided we know the density of the galaxy in which the KNR is located. The reader may notice that the peak flux for both SNRs and KNRs, assuming fiducial parameters, is the same in Fig. 1, which follows from Eq. (15). This may appear misleading, as the assumed external density encountered by KNRs is ten times smaller than that for SNRs. For a given density, the peak flux of a KNR is expected to be many times higher than that of a SNR, which becomes even more drastic in the case of EKNRs. Also, the spectral flux at different frequency bands (Fig. 2) will also be higher in the case of KNRs. In addition, the SNRs are expected to show distinctive line features (Woosley & Weaver 1986; Badenes et al. 2007) as well as morphologies (Lopez et al. 2009) in X-ray, which allows us to distinguish Type Ia from core-collapse supernovae. If we find a candidate object in an elliptical galaxy, we can almost certainly rule out a core-collapse supernova as the progenitor. Then, the non-detection of expected line features or morphology of Type Ia supernovae may help us in distinguishing a candidate KNR. Furthermore, as explained above, the peak luminosity for a KNR will be significantly larger. If we detect an object with L_3 GHz ≳ 10²⁷ erg s⁻¹ Hz⁻¹ associated with a remnant in an elliptical galaxy, this would be a very strong indication that the remnant is a KNR (or EKNR).

6. Discussion and conclusion

We have solved for the evolution of kilonova remnants (KNRs) and compare their light curves and spectral properties with supernova remnants (SNRs). While the physics is the same in both cases, we expect differences in quantitative properties due to different fiducial parameters, such as lower ejecta mass and lower external density in the KNR case. Spectral signatures such as the synchrotron cutoff, the correlation between radio and X-ray flux, and the transition from the synchrotron- to the SSC-dominated regime are features that can be used to distinguish KNRs from SNRs. As an example, around the peak of their flux, the synchrotron cutoff is of the order of 10 keV, 100 keV, and MeV for SNRs, KNRs, and EKNRs, respectively. We also note that despite the fact that most remnants are much older than the time at which their light curve peaks, due to the strong decrease in the luminosity with age during the decelerating phase, the observable sample is dominated by events observed close to their peak flux. However, as dN/dlog(t) ∼ t^−1/2, one in three remnants will be approximately ten times older than t_ST, which is not insignificant.

For fiducial parameters (E₅₁ = 1, M_ej, ⊙ = 0.1, n₀ = 0.1), we expect to detect KNRs in the near future using existing telescopes, such as VLASS. Due to the strong dependence of the observed flux on ejecta energy, extreme KNe (E₅₁ = 10) are much easier to detect if they constitute a sizeable portion of the total KN population. As there have been no reported detections of EKNRs in the literature, we are able to put moderate constraints on their abundance relative to regular KNRs, with f_EKN ≲ 0.3 assuming fiducial parameters. This constraint depends sensitively on the ejecta mass and external density distribution, as well as the neutron star merger rate, which is somewhat uncertain. We should be able to detect such objects or put strong constraints on their occurrence rate with ongoing surveys, such as VLASS and VAST, or future surveys with DSA-2000 and SKA (Hallinan et al. 2019; Braun et al. 2019).

The Sedov time for a EKNR is ∼20 years, which is at least an order of magnitude less than for a SNR. Depending on the external density, the Sedov time can be even shorter. If the remnants are detected deep in the Sedov-Taylor regime, we will not be able to easily distinguish KNRs from SNRs. However, if we do detect a remnant around its Sedov time, it may be possible to obtain a good estimate of t_ST with a few data points taken with a time interval of several years. A combination of estimated t_ST and corresponding spectral properties can provide strong evidence for the existence of a KNR. Also, if we find a candidate object in an elliptical galaxy, we can distinguish a KNR from a type Ia SNR using its brightness information or using line features and morphology in X-rays.

From the only robust KN identification to date, which is associated with the GW-detected NS merger GW170817, we can infer the ejecta mass and velocity for a given electron mass fraction Y_e. Different choices of Y_e can result in different interpretations of the properties of the ejecta (Metzger 2019). If we are able to infer ejecta mass and velocity from a remnant, we should be able to break the degeneracy and can independently infer Y_e; this will have implications in terms of r-process nucleosynthesis and help us understand how much radioactive material was synthesised during the formation of the KN.

Acknowledgments

SKA is supported by ARCO fellowship. PB’s work was supported by a grant (no. 2020747) from the United States-Israel Binational Science Foundation (BSF), Jerusalem, Israel, by a grant (no. 1649/23) from the Israel Science Foundation and by a grant (no. 80NSSC 24K0770) from the NASA astrophysics theory program. KH’s work was supported by JSPS KAKENHI Grant-in-Aid for Scientific Research JP23H01169, JP23H04900, and JST FOREST Program JPMJFR2136.

References

Appendix A: Relativistic equations for dynamical properties of kilonovae

In this section we provide fully relativistic equations for the evolution of ejecta. For this work, we find that relativistic corrections are mildly important for some extreme parameter cases. For such cases, we directly solve these equations numerically. The dynamical equation for the ejecta can be written as

$$\begin{array}{c}\hfill E=(\mathrm{\Gamma }-1)M_{\mathrm{ej}}{\mathrm{c}}^2+\frac{4\pi }{3}{n}_{\mathrm{ISM}}{R}_{\mathrm{ej}}^3{m}_{\mathrm{p}}{\mathrm{c}}^2\mathrm{\Gamma }(\mathrm{\Gamma }-1),\end{array}$$(A.1)

where E is the kinetic energy, Γ is the bulk Lorentz factor and R_ej is the size of ejecta in observer frame. The ejecta size is related to the observed time as R_ej = (1 + β_ej)Γ²β_ejct. The rest of the symbols have the same interpretation as in Sect. 2. In the coasting phase, $\mathrm{\Gamma }=1.0+\left(\frac{E_{51}}{M_{\mathrm{ej},\odot }}\right)\times 5\times {10}^{-4}$.

Assuming a fraction ϵ_e of the shocked energy is transferred to the electrons, by the same procedure as in Sect. 3, we obtain

$$\begin{array}{c}\hfill {\gamma }_m-1=\frac{p-2}{p-1}\frac{{ϵ}_e}{{\xi }_e}\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}}(\mathrm{\Gamma }-1).\end{array}$$(A.2)

The comoving magnetic energy density is given by

$$\begin{array}{c}\hfill B={(32\pi \mathrm{\Gamma }(\mathrm{\Gamma }-1)n_{\mathrm{ISM}}{m}_{\mathrm{p}}{\mathrm{c}}^2{ϵ}_B)}^{1/2}.\end{array}$$(A.3)

The synchrotron energy loss rate per electron in observer frame is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_e}{\mathrm{d}t}=\frac{4}{3}{\sigma }_{\mathrm{T}}\mathrm{c}{\mathrm{\Gamma }}^2{\gamma }_m^2{\beta }_m^2{B}^2,\end{array}$$(A.4)

where E_e is the electron’s energy. The characteristic frequency in observer frame is given by

$$\begin{array}{c}\hfill {\nu }_m=\mathrm{\Gamma }{\gamma }_m^2{\nu }_B,\end{array}$$(A.5)

The power in observer frame is then

$$\begin{array}{c}\hfill P_m=3.8\times {10}^{-22}\mathrm{\Gamma }\left(\frac{B}{1\mathrm{G}}\right){\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1},\end{array}$$(A.6)

assuming β_m → 1. The calculation of luminosity and flux follows the same way as in Sect. 3.1.

Appendix B: Cooling due to synchrotron

The Lorentz factor γ_c, sync for which the electrons cool within the timescale of expansion of the ejecta, i.e. the dynamical time, can be evaluated as, ${\gamma }_{\mathrm{c},\mathrm{sync}}{m}_{\mathrm{e}}{\mathrm{c}}^2=\frac{4}{3}{\sigma }_{\mathrm{T}}{\gamma }_{\mathrm{c},\mathrm{sync}}^2\frac{B^2}{8\pi }\mathrm{c}t$, or

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{c},\mathrm{sync}}=\frac{6\pi m_{\mathrm{e}}\mathrm{c}}{{\sigma }_{\mathrm{T}}{B}^{2}t}\end{array}$$(B.1)

for γ_c, sync > > 1 which can be simplified to, ${\gamma }_{\mathrm{c},\mathrm{sync}}=\frac{24.5}{B^2{t}_{\mathrm{yr}}}$. The cooling frequency is given by, ν_c, sync = γ_c, sync²ν_B. In the coasting phase, B is constant, hence, ν_c, sync goes like $t_{\mathrm{yr}}^{-2}$. Substituting the expression for B is ST regime, we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\nu }_{\mathrm{c},\mathrm{sync}}\sim 1.7\times {10}^{18}{ϵ}_{B,-2}^{-3/2}{E}_{51}^{-3/2}{n}_0^{-3/2}{M}_{\mathrm{ej},\odot }^{3/2}{t}_{\mathrm{yr}}^{-2}\mathrm{Hz}(\mathrm{coasting}\mathrm{phase})\\ \hfill {\nu }_{\mathrm{c},\mathrm{sync}}\sim 6.12\times {10}^{14}{ϵ}_{B,-2}^{-3/2}{E}_{51}^{-3/5}{n}_0^{-9/10}{t}_{\mathrm{yr}}^{-1/5}\mathrm{Hz}(\mathrm{ST}\mathrm{phase})\end{array}.\end{array}$$(B.2)

Appendix C: Synchrotron self-absorption frequency (ν_sa)

In the optically thin regime, the power per Hz at ν_m per steradian is given by

$$\begin{array}{c}\hfill P_{{\nu }_m}=\frac{3.8\times {10}^{-22}}{4\pi }\left(\frac{B}{1\mathrm{G}}\right){\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}{\mathrm{sr}}^{-1}.\end{array}$$(C.1)

Assuming ν_sa > ν_m, the power at ν_sa is given by, $P({\nu }_{\mathrm{sa}})=P_{{\nu }_m}{\left(\frac{{\nu }_{\mathrm{sa}}}{{\nu }_m}\right)}^{-(p-1)/2}$. The intensity at ν_sa is given by

$$\begin{array}{c}\hfill I_{\nu }=0.3\times {10}^{-3}B{\left(\frac{{\nu }_{\mathrm{sa}}}{{\nu }_m}\right)}^{-(p-1)/2}\left(\frac{R}{{10}^{19}cm}\right)n_0{\xi }_e{\mathrm{erg}\mathrm{cm}}^{-2}{\mathrm{Hz}}^{-1}{\mathrm{sr}}^{-1}.\end{array}$$(C.2)

The intensity of a blackbody is given by

$$\begin{array}{c}\hfill I_{\nu }=\frac{2k_{B}T}{2.7}\frac{{\nu }_{\mathrm{sa}}^2}{{\mathrm{c}}^2},\end{array}$$(C.3)

where k_BT ∼ γ_sam_ec² with ${\gamma }_{\mathrm{sa}}=\sqrt{{\nu }_{\mathrm{sa}}/{\nu }_B}$. Equating the two and using p = 2.5, the expression for ν_sa is given by

$$\begin{array}{c}\hfill {\nu }_{\mathrm{sa}}^{3.25}=15\times {10}^{31}{B}^{2.25}{R}_{19}{n}_0{\xi }_e,\end{array}$$(C.4)

or ν_sa ∼ 8 × 10⁹(B^2.25R₁₉n₀ξ_e)^1/3.25 Hz. Substituting the expressions for the physical parameters, we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\nu }_{\mathrm{sa}}\sim 3\times {10}^6{{\overline{ϵ}}_{e,-1}}^{4/13}{ϵ}_{B,-2}^{9/26}{E}_{51}^{21/26}{M}_{\mathrm{ej},\odot }^{-21/26}{n}_0^{17/26}{t}_{\mathrm{yr}}^{0.31}\mathrm{Hz}(\mathrm{coasting}\mathrm{phase})\\ \hfill {\nu }_{\mathrm{sa}}\sim 4\times {10}^8{{\overline{ϵ}}_{e,-1}}^{4/13}{ϵ}_{B,-2}^{9/26}{E}_{51}^{21/65}{n}_0^{43/130}{t}_{\mathrm{yr}}^{-43/65}\mathrm{Hz}(\mathrm{STphase})\end{array}.\end{array}$$(C.5)

Appendix D: Computation of SSC flux

The number flux of inverse Compton photons is given by (Blumenthal & Gould 1970; Sari & Esin 2001)

$$\begin{array}{c}\hfill f_{\nu ,\mathrm{IC}}={\sigma }_T{R}_{\mathrm{ej}}{\displaystyle {\int }_{{\gamma }_m}^{\infty }d{\gamma }_e\frac{\mathrm{d}{N}_e}{\mathrm{d}{\gamma }_e}}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{\nu }_{s}g(x)f_{{\nu }_s}(x),\end{array}$$(D.1)

where $x=\frac{\nu }{4{\gamma }_e^2{\nu }_s}$, g(x) = 1 + x − 2x² + 2xlog(x) in Thomson approximation and $f_{{\nu }_s}=\frac{L_{{\nu }_s}}{h{\nu }_s}$ where L_νs is obtained in Sect. 3.1. Technically, the flux on both sides of the equation is at the location of shocked gas but they can be transformed to the observed flux with the same factor of $\frac{R_{\mathrm{ej}}^2}{D^2}$ which cancels out. The factor g(x) has support only over 0 ≲ x ≲ 1 which effectively cuts off the ν_s integral for a given ν. Expanding the definition of the quantities, we have

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\gamma ,\mathrm{I}C}}{\mathrm{d}t\mathrm{d}\nu }={\sigma }_T{R}_{\mathrm{ej}}{\displaystyle {\int }_{{\gamma }_m}^{\infty }}d{\gamma }_e\frac{\mathrm{d}{N}_e}{\mathrm{d}{\gamma }_e}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{\nu }_{s}g(x)f_{{\nu }_s}(x)\frac{1}{4{\gamma }_e^2{\nu }_s}.\end{array}$$(D.2)

We solve this equation numerically and obtain the number flux which can be converted to L_ν by multiplying hν. We compare synchrotron, SSC and bremsstrahlung flux for our fiducial cases in Fig. 2.

Appendix E: Bremsstrahlung emission

The thermal energy of electron is related to the shock energy via electron kinetic energy. Therefore, the temperature of electron is given by

$$\begin{array}{c}\hfill {\overline{ϵ}}_e\frac{1}{2}{m}_{\mathrm{p}}{\beta }_{\mathrm{ej}}^2{\mathrm{c}}^2=\frac{3}{2}{k}_B{T}_{\mathrm{e}},\\ \hfill \end{array}$$(E.1)

which after substitution of β_ej becomes

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T_{\mathrm{e}}=3\times {10}^8{\overline{ϵ}}_{e,-1}\frac{E_{51}}{M_{\mathrm{ej},\odot }}\mathrm{K}(\mathrm{coasting}\mathrm{phase})\\ \hfill T_{\mathrm{e}}=0.88\times {10}^{11}{\overline{ϵ}}_{e,-1}{E}_{51}^{2/5}{n}_0^{-2/5}{t}_{\mathrm{yr}}^{-6/5}\mathrm{K}(\mathrm{ST}\mathrm{phase})\end{array}.\end{array}$$(E.2)