Long-term evolution of a supernova remnant hosting a double neutron star binary

An ultra-stripped supernova (USSN) is a type of core-collapse SN explosion proposed to be a candidate formation site of a double neutron star (DNS) binary. We investigate the dynamical evolution of an ultra-stripped supernova remnant (USSNR), which should host a DNS at its center. By accounting for the mass-loss history of the progenitor binary using a model developed by a previous study, we construct the large-scale structure of the {circumstellar medium (CSM)} up to a radius $\sim 100\,{\rm pc}$, and simulate the explosion and subsequent evolution of a USSN surrounded by such a CSM environment. We find that the CSM encompasses an extended region characterized by a hot plasma with a temperature $\sim 10^8\,$K located around the termination shock of the wind from the progenitor binary ($\sim 10\,$pc), and the USSNR blastwave is drastically weakened while penetrating through this hot plasma. Radio continuum emission from a young USSNR is sufficiently bright to be detectable if it inhabits our Galaxy but faint compared to the observed Galactic SNRs, and thereafter declines in luminosity through adiabatic cooling. Within our parameter space, USSNRs typically exhibit a low radio luminosity and surface brightness compared to the known Galactic SNRs. Due to the small event rate of USSNe and their relatively short observable lifespan, we calculate that USSNRs account for only $\sim 0.1$-$1$ % of the total SNR population. This is consistent with the fact that no SNR hosting a DNS binary has been discovered in the Milky Way so far.


INTRODUCTION
A double neutron star (DNS) binary is believed to be the fossil object from a binary system of two massive stars which have both exploded as core-collapse supernovae (SNe) in the past (e.g., Podsiadlowski et al. 2005). Observations of Galactic radio pulsars have revealed that some DNS binaries are in an orbit tight enough to merge within the cosmic age (Burgay et al. 2003). Indeed, previous observations for the short gamma-ray burst GRB 130603B have implied the association between the gamma-ray emission and kilonova in the DNS merger (Tanvir et al. 2013;Hotokezaka et al. 2013). Furthermore, recent gravitational wave detectors and rapid follow-up electromagnetic observations have succeeded in probing the coalescence of a DNS, confirming the link of these objects to the origin of short gamma-ray bursts and the nucleosynthesis of r-process elements (e.g., Abbott et al. 2017a,b;Tanaka et al. 2017).
The formation of a DNS requires that the binary system is not disrupted by the evolution history of the massive stars all the way through their core-collapses. One of the plausible scenarios of DNS formation invokes an ultra-stripped supernova (USSN, Tauris et al. 2017;Yoshida et al. 2017). In a close binary consisting of two massive stars, the primary star first explodes as a SN. After a phase as a high-mass X-ray binary, the outer layer of the secondary star is stripped away in two steps: (1) the ejection of its hydrogen-rich envelope through a phase of common envelope (CE) interaction, and (2) the stripping of the helium layer through Roche lobe overflow (RLO). These binary interactions lead to the formation of an helium star ( 2M ), which eventually explodes as a USSN. Indeed, some of the rapidly evolving transients such as SN 2005ek (Drout et al. 2013), iPTF14gqr (De et al. 2018), and SN 2019dge (Yao et al. 2020) are suggested to be possible candidates for USSNe (Moriya et al. 2017). In addition, it has been proposed that during the operation period of the Zwicky Transient Facility (ZTF, Graham et al. 2019;Bellm et al. 2019), roughly 10 USSNe within 300 Mpc will be detected per a year (Hijikawa et al. 2019). Hence, it is expected that future surveys and follow-up observations of transients will enable us to examine in detail the validity of the USSN scenario as the formation mechanism of DNS binaries.
Another way to experimentally test the USSN scenario is to search for supernova remnants (SNRs) hosting a DNS binary. After the explosion, the ejecta of the USSN sweeps up the surrounding CSM while expanding into the interstellar space. Intriguingly, this kind of system can be potentially detected as a SNR hosting a DNS binary, which we will refer to as an ultra-stripped supernova remnant (USSNR) hereafter. While the current SNR surveys have not identified any of these remnants so far, we note that the observable characteristics of a USSNR have not been discussed and quantified in the literature. It is hence essential to investigate the dynamical evolution and emission properties of USSNRs using a dedicated simulation model to shed light on how they can be identified. Tauris et al. (2013) developed a progenitor evolution model for the USSN, and showed that the masstransfer rate through RLO can be enhanced up toṀ ∼ 10 −5 M yr −1 in the last 0.1 Myr prior to the core collapse. Because the mass-transfer rate is orders of magnitude larger than the Eddington accretion rate onto the neutron star, a large fraction of the stripped gas escapes the binary system and distributes around the progenitor as CSM. Assuming a wind velocity v w ∼ 1000 km s −1 , the gas which has been expelled from the binary system in the RLO phase can reach a distance of ∼ 100 pc from the progenitor, implying that the evolution of the USSNR is heavily influenced by the CSM created by the RLO mass loading process. However, detailed models for the mass-loss history driven by binary interaction are in most cases not incorporated in the simulations of SNR dynamics, which is particularly critical for understanding the properties of USSNRs.
In this study, we investigate the characteristics of a USSNR using a grid of one-dimensional hydrodynamic simulations. By employing the binary evolution model presented in Tauris et al. (2013), we first construct the large-scale structure of the CSM surrounding the USSN progenitor. We next calculate the hydrodynamics of the USSN ejecta interacting with the composed CSM and the resulted synchrotron radiation. Our simulations reveal that the blastwave of USSNRs has a difficulty in penetrating the hot plasma, which had been shaped by the preceding mass loss from the progenitor binary. Radio emission from a young USSNR is predicted to be bright enough to be detected if it inhabits our Galaxy, while its luminosity starts to decrease at t 10 3 years, making the USSNR observable for a relatively short time period. Besides, the low surface brightness of a USSNR predicted by our models at its typical diameters (D ∼ O(10 pc)) can serve as a key to the identification of these remnants in the future.
This paper is organized as follows. In Section 2, we review the USSN scenario as a formation theory of DNS, and describe the progenitor models used in our simulations. In Section 3, we discuss the formation sequence of the CSM, followed by a description of the procedures for constructing our CSM models. In Section 4, we examine the hydrodynamic evolution of a USSNR and show the properties of the expected radio signals, including the light curve and surface brightness. Their implications are discussed in Section 5, and our results are summarized in Section 6.
2. PROGENITOR MODEL Tauris et al. (2013) investigated the binary stellar evolution of a 2.9M He star with a neutron star companion, having an initial orbital period of 0.1 day. They found that the He star reduces its own mass down to 1.5M through RLO, and suggested that the He star explodes as a USSN which can be a candidate for some rapidly evolving transients. Here we overview the stellar evolution of the progenitor of a USSN, which is crucial for understanding the formation of the CSM adopted in this study. Figure 1 shows the time evolution of the mass (M ), radius of the Roche lobe (R), escape velocity (V esc ), and mass-transfer rate (Ṁ ) of the USSN progenitor presented in Tauris et al. (2013). Here, the escape velocity is defined as V esc = 2GM/R, where G is the gravitational constant. When the progenitor is in the state illustrated by the blue line, its outer layer is stripped away by the companion neutron star through RLO. Until the core collapse, the He star experiences RLO three times; the first phase is at 1.78 Myr t 1.84 Myr during which the core has exhausted its He-burning fuel (A). The second is at t ∼ 1.851 Myr when the core C-burning has ended (B), and the third is at t 1.854 Myr in which the off-center O-burning is about to onset (C and D). The CSM around the USSN progenitor is hence expected to be shaped by these three phases of mass loss activities. We note that the progenitor spends most of its lifetime in the state shown by the orange line prior to A, and that the increase of the mass-loss rate is realized in the last 0.1 Myr before core collapse.
The progenitor does not experience RLO in the detached phases, during which we conservatively assume a mass loss rate of 10 −7 M yr −1 . This mimics the stellar wind from the progenitor, but the mass and kinetic energy released by this wind are smaller than those carried by the gas stripped away through the RLO. Thus, we can assume that the stellar wind from the progenitor has an insignificant influence on the overall wind hydrodynamics, and that the consistency with the stellar evolution model is maintained. The model developed by Tauris et al. (2013) covers the lifetime of the He star only until 10 years prior to its core collapse. To trace the evolution up to the moment of the explosion, we use the final values of M, R, V esc , andṀ from the model for the last 10 years.
The gas transferred from the progenitor first flows toward the neutron star with an accretion rate orders of magnitude larger than the Eddington accretion rate (Tauris et al. 2013). The neutron star cannot feed up anymore and thus drives the accreted gas outward by mechanisms such as propeller effect (Tauris et al. 2017). However, resolving the detail of this outflow dynamics is beyond the scope of this work. For simplicity, we assume that the material which has been stripped away from the He star launches outward spherically at the radius of the Roche lobe R with a velocity V esc and massloss rateṀ . Then, the mass density at the Roche lobe radius (Ṁ /4πR 2 V esc ) can be estimated. Given the density and velocity of the gas at the Roche lobe radius as an inner boundary condition, we can solve the hydrodynamics of the gas launched from the progenitor binary to model the CSM formation around the progenitor. Combined with a parametric survey described in the following sections, this strategy allows us to demonstrate the long-term evolution properties of a USSNR with the mass loss history of the progenitor taken into account.

CSM FORMATION
In this section, we describe our procedure for modeling the formation of the CSM surrounding the USSN progenitor. First, we construct the initial profile of the interstellar medium (ISM) in Section 3.1. We then explain our methodology for simulating the hydrodynamics of the mass-loss material in Section 3.2, and the properties of the composed CSM in Section 3.3.

Initial setup
The progenitor experiences a hydrogen-rich envelope ejection driven by the CE interaction before the stripping of the helium gas through RLO. The distribution of this expelled hydrogen-rich gas is important because it interacts with the helium gas released through RLO later on. Although some recent multi-dimensional simulations have succeeded in completely ejecting the hydrogen-rich envelope of a red supergiant through the CE interaction under some assumptions and realizations (Law-Smith et al. 2020;Lau et al. 2021, but see also Vigna-Gómez et al. 2021), the distribution of the material ejected by the CE interaction is still not completely understood. Figure 2 shows three models we adopt for the initial density profile of the CE material. We consider a situation where the ejected gas with a mass M CE = 10M is distributed within a radius R CE which smoothly connects with the ISM. Given that the characteristic timescale of the CE interaction is around thousands of years (Ivanova et al. 2013), the gas ejected with a speed ∼ 100 km s −1 can reach a radius R CE ∼ 10 18 cm. Since there is a variety in the ISM properties such as density and temperature (e.g., Berkhuijsen & Fletcher 2008;Draine 2011), we consider two ISM phases; a warm phase (ρ ism = 10 −24 g cm −3 , T ism = 10 4 K) and a hot phase (ρ ism = 10 −26 g cm −3 , T ism = 10 6 K). We remark that the thermal pressure in these two initial profiles are equal to each other. In addition, we prepare a reference model 'UNIFORM', in which a static and uniform ISM resides throughout the simulation domain with a density ρ ism = 10 −24 g cm −3 , to evaluate the effect of the CE ejection activity. The specific profiles of the initial density for each model are described in Table 1. The derivation of the exact value of ρ CE is explicated in Appendix A. We consider a static ISM profile (v = 0). The initial velocity profile of the CE component does not have an important role in the hydrodynamics of the CSM formation because the expected V CE is negligibly lower than the velocity of the wind from the progenitor binary. To verify this we conducted simulations in which the initial velocity of the CE component is assumed to be 100 km s −1 and confirmed that the outcome is not changed. We assume the temperature T = T ism and a solar metallicity throughout the entire profiles at this stage. A comparison of the results among these models enables us to evaluate how much the properties of the   CE ejection affect the CSM formation and the subsequent SNR evolution.

Wind hydrodynamics
We solve the one-dimensional equations of ideal gas hydrodynamics where the internal energy is taken away by radiative cooling in spherical coordinates. The governing equations are described as follows: where ρ is the mass density, v is the velocity, p is the pressure, e is the specific internal energy, h = e + p/ρ is the specific enthalpy, n i and n e are the number density of ions and electrons. Λ(T ) represents the radiative cooling function, for which we employ the power-law formalism introduced by Chevalier & Fransson (1994). The energy loss by radiative cooling is calculated only in the optically thin region where τ ≤ 1 which is sufficient for tracing the evolution of the blastwave (see also Section 5.4). These governing equations are closed with the equation of state, p = (γ − 1)ρe, where γ = 5/3 is the adiabatic index. The equations are solved by a Roe Riemann solver with the second entropy fix by Harten and Hyman to treat the contact discontinuity and the shock wave (Harten & Hyman 1983). The numerical accuracy of the code used in this study is verified in Appendix B.
We divide the simulation domain from 10 16 cm to 3 × 10 21 cm into 2047 zones in a logarithmic scale. Inside 10 16 cm as an inner boundary condition, we inject the blowing He-rich gas whose time evolution is described in Section 2.
We trace the distribution of the chemical abundances by advection, assuming that no mixing of the chemical composition occurs. The abundance distribution is required in order to accurately estimate the number density of ions and electrons in the radiative cooling term. Figure 3 shows the snapshot of the density structure of the CSM at the moment of core collapse of the progenitor. The model 'WARM' and 'UNIFORM' share an identical CSM structure in the entire simulation domain. It is also the case for the model 'HOT' within ∼ 3 pc, but its outer configuration deviates from the other two models. The distribution of the density within ∼ 3 pc reflects the mass loss history. Namely, the dense CSM being distributed around r ∼ 0.01 pc and 0.1 pc are originated from the mass loss at points D and C in Figure 1. Yet, a segment resides around 10 pc in which the density is roughly constant with some fluctuations. This nonsmooth segment is created by the collision between the wind launched at point B and the reverse shock generated by the gas ejected at point A before. The ISM wall is located at a radius of 20 pc in the model 'WARM' and 'UNIFORM' and 30 pc in the model 'HOT', respectively.

Composed CSM
We will briefly elaborate on the importance of the CE component on the ISM profile. The reference model The dashed black line shows the distribution realized for the steady wind with its mass-loss rateṀ = 10 −7 M yr −1 . The distributions of the gas pointed out by cursive alphabets represent that they are from the mass loss activity referred in Figure 1.
'UNIFORM' without the CE component allows us to investigate the contribution of the CE component on the hydrodynamics of the wind. The results obtained from this reference model are found to be almost identical to the outcome from 'WARM', being nearly indistinguishable just in Figure 3. This can be interpreted as follows. The radius of the ISM wall is roughly determined by the balance between the ram pressure of the wind and the thermal pressure of the swept-up material (Weaver et al. 1977), which is computed as ∼ 20 pc in our simulations. The enclosed mass of the initial ISM profile at r ∼ 20 pc is ∼ 400M , indicating that the mass of the CE component can be regarded to be negligibly small. Hence, the composed CSM has similar characteristics between 'WARM' and 'UNIFORM'. We confirmed that even when considering a uniformly distributed hot ISM (ρ ism = 10 −26 g cm −3 , T ism = 10 6 K), the consequent CSM structure does not differ from the model 'HOT' significantly other than slight quantitative modifications. This implies that as long as the CE ejection before the USSN is considered within a range of typical time and energy scales, it does not play an important role in the formation of the CSM around the USSN progenitor. Figure 4 shows the temperature structure of the CSM at the moment of core collapse. Similar to the density structure, the models 'WARM' and 'UNIFORM' have the same temperature structure over the entire region. The model 'HOT' also possesses the identical distribution with the other models within 3 pc, but the quantitatively different structure is formed outside 3 pc. A hot plasma with ∼ 10 8 K is located in the vicinity of the ISM wall in all models. The location of the inner edge of this hot plasma coincides with the radius of the termination shock of the wind driven by RLO. The geometrical thickness of the plasma is ∼ 10 pc. The existence of this hot plasma region plays a critical role in weakening the SNR blastwave as it propagates through the region as discussed later.

SNR EVOLUTION
In this section, we investigate the evolution of a USSNR interacting with the CSM constructed in the previous section. In Section 4.1, we show the method to simulate the dynamics of the ejecta and the expected synchrotron emission, and the results are presented in Section 4.2. As was confirmed in the previous section, the solution derived from the model without a CE component ('UNIFORM') converges to that of 'WARM'. We will therefore examine results from the models 'WARM' and 'HOT' hereafter.

Ejecta dynamics
The initial profile of the USSN progenitor is taken from Moriya et al. (2017), who evolved the model of the He star previously presented by Tauris et al. (2013) further until core collapse. Then we attach the CSM composed in Section 3 to the progenitor while retaining the distribution of the density, velocity, temperature, and chemical abundance.
We next examine the hydrodynamics of the SN explosion to obtain the SN ejecta structure. We excise the remnant mass M rem = 1.35M from the inner region of the progenitor, and inject an explosion energy E exp = 10 50 erg to the rest of the material in the progenitor (M ej ∼ 0.15M ) as a thermal energy following the method developed by Morozova et al. (2015). The explosion energy is chosen based on light curve models (Moriya et al. 2017), which is also consistent with that proposed by state-of-the-art simulations (Suwa et al. 2015;Müller et al. 2018). The profile is resolved into more than 4000 meshes with a logarithmic spacing, and the hydrodynamics of the ejecta is calculated by the same method as described in Section 3.2, except that a reflective condition is employed at the inner boundary. As a result, we obtain the time evolution of the blastwave velocity and the trajectory of Lagrangian particles, which are used to compute the energy distribution of relativistic electrons and the amplified magnetic field (see the next section).
As the SNR evolves into the Sedov phase, its reverse shock begins to propagate towards the inner region and heats up the ejecta (e.g., Truelove & McKee 1999). Since the simulation domain is resolved under a logarithmic mesh, the high temperature in the inner region can cause small timesteps, making it difficult for the simulation to progress. To solve this numerical difficulty, we excise the Eulerian meshes in the innermost region within 10 18 cm when the blastwave radius has reached 10 19 cm. This does not affect the consistency of the simulations since the total gas mass within 10 18 cm at the moment of the excision is negligibly small and hence dynamically unimportant. This allows us to trace the long-term evolution of the USSNR within a reasonable simulation time. The computations are terminated at 10 5 years since the explosion.

Particle acceleration and magnetic field amplification
Once the gas is heated by the forward shock, the diffusive shock acceleration (DSA) imparts relativistic energies to the injected charge particles and induces amplification of the turbulent magnetic field (e.g., Fermi 1949;Drury 1983). The region shocked by the blastwave serves as a site of synchrotron emission from SNRs (Reynolds 2008;Dubner & Giacani 2015). In this study, we define the blastwave as the discontinuity which satisfies the following two conditions: (1) the pressure jump is the largest in the simulation domain, and (2) the Mach number is greater than 3. The latter is justified because strong shocks have a potential to drive DSA, whilst weak shocks are less capable of efficient particle acceleration, confirmed by the observations for radio relics in galaxy clusters (e.g., Botteon et al. 2020, and references therein).
We first consider a Lagrangian mesh a s through which the blastwave passes at time t s . As the shock sweeps through the mesh, the charged particles are accelerated to relativistic energies, coupled with an amplification of the turbulent magnetic field. We model the energy densities of the accelerated relativistic electrons (u e ) and the magnetic field (u B ) in the Lagrangian mesh a s as follows: where e and B are the acceleration and amplification efficiencies, ρ sh is the mass density in the Lagrangian mesh a s , V b is the velocity of the blastwave, and v u is the velocity of the unshocked gas upstream of the shock, respectively. These parametrizations are conventionally used in the modeling of radio SNe (e.g., Chevalier et al. 2006;Chevalier & Fransson 2006;Matsuoka et al. 2019). These equations apply to the mesh only when The energy distribution of the accelerated electrons, N (a s , E), is described by a power-law distribution as follows: where E and p are the energy and the spectral index of the electrons, respectively. The coefficient C is determined by performing a normalization of the energy density: where E min = 2m e c 2 is used (see Section 5.6 for a discussion on the uncertainty related to E min ).
As the system evolves, the ejecta expands and the blastwave propagates to the next Lagrangian mesh. Meanwhile, the relativistic electrons lose their energies by both synchrotron and adiabatic cooling, and the magnetic field also decays with the adiabatic expansion. We consider a Lagrangian mesh (a) which had been heated by the shock at mass coordinate a s and time t s , and assume that the relativistic particles are confined within the mesh and the magnetic field is frozen in the plasma. We calculate the cooling processes of the accelerated particles and the time evolution of the energy distribution following previous studies (e.g., Reynolds 1998;Orlando et al. 2011;Ferrand et al. 2014). An electron's energy E declines to E through synchrotron and adiabatic cooling, which can be written as follows: where c is the speed of light, q is the elementary charge, and m e is the mass of electron, respectively. The energy distribution of the electrons evolves following number conservation, i.e.
As for the strength of the magnetic field, we consider a magnetic flux conservation in each Lagrangian mesh.

Synchrotron emission
Given the energy distribution of electrons and the strength of the magnetic field, the intensity of the synchrotron emission I ν can be calculated by integrating the radiative transfer equation written as follows: where j ν,syn , α ν,syn , and α ν,ff are the synchrotron emissivity, synchrotron self-absorption and free-free absorption coefficient, respectively (Rybicki & Lightman 1979), and R b is the blastwave radius. We also calculate the surface brightness Σ(θ) which is often used as a diagnostic observable for SNRs. Σ-D diagrams which show the relation between the surface brightness and the diameter of SNRs are commonly used for determining the distance to the objects (see e.g., Poveda & Woltjer 1968;Pavlović et al. 2013, and references therein). Since the surface brightness Σ(θ) is independent of the distance to the SNR, it can be a useful quantity for investigating the intrinsic nature of the USSNR compared to the rest of the SNR population. Σ(θ)δθ, the power per unit surface area and unit frequency emitted from a ring with sky projection angles θ to θ + δθ, can be evaluated by integrating the total power of the synchrotron emission per unit volume along the line of sight as follows: where d, ν = 4πj ν,syn , δA(θ) = δ(πd 2 θ 2 ), and ∆Ω(θ) are the distance to the SNR, the total power of the synchrotron emission per unit volume, the area of the ring with projection angle θ, and the total solid angle of the SNR. The angle-averaged surface brightness can then be estimated, which allows us to examine the position of USSNRs on the Σ-D diagram.

Characteristics of a USSNR
Firstly, we discuss the hydrodynamics of the interaction between the USSN ejecta and the CSM. In Figure 5, the time evolutions of the density and velocity profile are shown. Here we mention on the dependence of the density profile on the ISM state. We can see that the model 'HOT' has a larger radius of the ISM wall than the model 'WARM', even though these two models have initially the same pressure. This suggests that the ISM density is important for dictating the location of the ISM wall; the lower ISM density (hot ISM) allows the exploding SNR gas to further expand. This feature is critical for quantifying the surface brightness of the USSNRs (see Figure 9 and Figure 10).
From the density distributions, we can see that the ejecta keeps expanding until t ∼ 10000 years but starts decelerating around the ISM wall. The system can expand further for another ∼ 3 and 10 pc at most from the location of the ISM wall in the model 'WARM' and 'HOT', respectively. This can be observed in the panel of the velocity profile; the system experiences fast expansion at t 3000 years, while after the collision with the ISM wall it only possesses several hundreds km s −1 of the outward velocity. This implies that the diameter of the USSNR is highly constrained by the location of the ISM cavity wall, which in turn depends on the pre-SN mass loss activity of the progenitor. This picture can be applied to all core-collapse SNRs in general, for which the diameters of SNRs are associated with the pre-SN mass loss activity of their progenitors (e.g., Yasuda et al. 2021a,b). Figure 6 shows the time evolution of the Mach number and the blastwave velocity. Within the first 300 years, these two quantities both in the model 'WARM' and 'HOT' behave similarly each other since in this phase the identical CSM structure is traced. We can see two epochs at which the blastwave accelerates at r ∼ 0.01 pc and r ∼ 0.1 pc respectively, where the CSM density drops by orders of magnitude. Correspondingly, the Mach number also increases by more than an order of magnitudes at the same time. Overall, the velocity stays at about 10 9 cm s −1 , leaving the USSNR active for the first 300 years. Furthermore, at 5 years t 50 years when the swept CSM mass begins to exceed the ejecta mass, the velocity of the blastwave decays proportional roughly to t −1/3 . This agrees with the expected time dependence of the velocity in the Sedov phase for a CSM density profile proportional to r −2 (Book 1994). The gradual increase of the Mach number during that phase can be also observed, due to the decrease of the upstream temperature (see Figure 4).
After t ∼ 300 years, the blastwave decelerates down to ∼ 10 8 cm s −1 , and then simply disappears out, as well as   the Mach number decreases rapidly down to O(1). This phenomenon can be observed both in 'WARM' and in 'HOT' though there are some quantitative differences between these two models. This is caused by the hot plasma at r ∼ 5 pc shown in Figure 4; as the blastwave plunges into the plasma where the sound speed is high, the Mach number of the blastwave quickly decreases down to unity. It is implied that such a weak shock cannot support an efficient DSA. Additionally, the density jump at r ∼ 3 pc can also give rise to the deceleration of the blastwave. In conclusion, this result indicates that the blastwave in a USSNR dies out by propagating into a region of hot plasma at 10 3 years. Figure 7 shows the long-term 1 GHz radio light curves from the models shown in Table 2. The observed flux density F ν shown in the right y-axis is normalized by a distance d = 10 kpc. The peak luminosity of the light curve is determined by synchrotron self-absorption with their shapes slightly modified by free-free absorption (see also Matsuoka et al. 2019). Note that for a USSN candidate iPTF14gqr non-detections of radio signals at the frequency 6 GHz and 22 GHz within 10 days have been reported, placing upper limits (De et al. 2018). In such a very early phase, free-free absorption completely damps the centimeter radio emissions, much more for 1 GHz (Matsuoka & Maeda 2020). Since more electrons are accelerated and magnetic field is more intensively amplified in the models which assume larger efficiencies for DSA, brighter radio emission from USSNRs can be expected in the model with ( e , B ) = (10 −2 , 10 −1 ) than  those with ( e , B ) = (10 −3 , 10 −2 ). Besides, a harder spectral index increases the number of more energetic electrons in the shocked region, which also results in the luminous radio signals. This behavior can be confirmed by comparing the luminosity between the models with p = 2.1 and those with p = 2.5, 3.0. We note that there are no qualitative difference in the light curve behaviors between the two CSM models over the entire timespan up to 10 5 years. Actually as 'HOT' has a more extended structure than 'WARM' as seen in Figure 5, a difference between these two models is expected in their surface brightness as we will discuss later.
We first look at the behaviors of the young USSNR at ages less than 1000 years, and compare them with SNe well-observed at the frequency ∼ 1 GHz even 1 year after their explosions such as SN 1993J (Martí-Vidal et al. 2011), SN 1995N (Chandra et al. 2009), and SN 2006jd (Chandra et al. 2012, and one of the youngest Galactic SNR Cas A (DeLaney & Rudnick 2003, the point plotted at t ∼ 300 years). As seen in Figure 7, our models show that young USSNRs at an age t ∼ 10 years and t ∼ 300 years produce fainter radio signal than those from the bright SNe and Cas A, respectively. The relatively weak emissions can be partially attributed to the shock velocity which is by a factor of a few lower than what is inferred for these objects (see, e.g., Fransson & Björnsson 1998). Another possible reason is that at t ∼ 100 years the blastwave is propagating at r ∼ 1 pc where the dense CSM formed by the mass loss driven by the RLO is absent. Then the density of the CSM swept by the blastwave is considerably small there, making the DSA less efficient. However, we note that the expected flux density of the radio emission from the USSNR at d = 10 pc keeps greater than 0.1 mJy within an age t 1000 years, which is bright enough to be detected by the present radio surveys such as Very Large Array Sky Survey (VLASS, Lacy et al. 2020), if it inhabits inside our galaxy.
Next we discuss the properties of the light curves of USSNR at larger ages (1000 years t 10 5 years). At t ∼ 1000 years, the radio emission brightens by a factor to an order of magnitude compared to t ∼ 300 years, even though the synchrotron emission in this phase is optically thin to self-absorption. This enhancement stems from the interaction between the SN ejecta and the relatively dense CSM located at ∼ 3-10 pc; a larger amount of the gas injection into the shocked region leads to a larger number of the synchrotron emitting electrons, resulting in a higher radio luminosity. In addition, the compression of the gas around the blastwave by the collision with the dense CSM brings about the further amplification of the magnetic field through the conservation of the magnetic flux (see Figure 8). This can also be a cause of the brightening of the radio luminosity. We note that this brightening is one of the characteristics of a USSNR associated with the time dependent mass loss driven by RLO, since a CSM with a simple powerlaw distribution cannot reproduce such a rise in radio luminosity in the optically-thin regime. Yet, the subsequent radio signals are fainter than those observed from the Galactic SNRs enumerated in Table 3. The stalled blastwave at t ∼ 300 years can no longer execute efficient DSA any further. Even so, it is worth mentioning that SNRs discovered so far are biased towards bright  Table 3 (black points with error bars), estimated by the distances to each objects. The right y-axis stands for the observed flux densities with which the source with the luminosity shown in the left y-axis is observed at a distance d = 10 kpc. The red dotted line indicates the detection limit of VLASS (Lacy et al. 2020).
objects. Deep surveys such as VLASS will have potential to uncover the population of the SNRs as faint as the aged USSNRs. After the death of the blastwave, DSA will no longer be triggered, and the non-thermal emissions are forced to decline through adiabatic cooling. The timing of dominance by adiabatic cooling is roughly 1000 years, and is more-or-less determined by the location of the hot plasma (Figure 4). The hot plasma is formed by the interaction between the He-rich wind blown from the progenitor binary and the H-rich gas originated from the CE ejection or the uniform ISM. Our result implies that the location of the hot plasma in the CSM is key to determining the lifetime of the blastwave and hence the observable lifespan of the USSNR.
We also observe oscillations of the light curves at t 10 4 years. This is an one-dimensional artifact due to the reflective condition at the inner boundary of the simulation domain. As the reverse shock of the USSNR brings along an inward gas flow back to the explosion center, it rebounds back to the outer interacting region. Then the material around the shocked region is compressed, inducing an amplification of the magnetic field through flux conservation. A repeating occurrence of this inward and outward motion results in the oscillation of the radio luminosity in our models for the aged USSNR. In practice, multi-dimensional dynamics should suppress the motion of the gas described above due to a broken spherical symmetry. Even so, it can be noted that the global evolution of the radio luminosity of the aged USSNR roughly follows an adiabatic evolution when averaged over a longer timescale. Figure 9 shows the time evolution of the surface brightness as a function of the sky projection angle. The model 'HOT' has fainter surface brightness and larger projection angles at which the surface brightness becomes maximum (θ max ) than those in the model 'WARM', because the model 'HOT' has a more extended CSM density structure than the model 'WARM' (see Figure 3). Yet the qualitative behavior of the surface brightness as a function of the sky projection angle is similar between these two models. θ max is mainly dic-tated by the location of the ISM wall, which prevents the gas in the shocked region from expanding any further outward (see Figure 5). As mentioned before, the hot plasma and the ISM cavity wall are shaped by the wind colliding with the CE and/or the ISM, which ultimately determines the detectability of the USSNR.
The evolution of the relation between the surface brightness and diameter of the USSNR can be assessed by the Σ-D diagram shown in Figure 10. For the same reason as the relation between Σ and θ max (Figure 9), the model 'HOT' has a fainter surface brightness and larger diameter than the model 'WARM'. This results in the lower right position of the evolutionary path of the model 'HOT' in the Σ − D diagram. The magnitude of the surface brightness strongly depends on the parameters relevant to the DSA (i.e., p, e , and B ). The surface brightness of the model appears to be relatively faint compared to those of the Galactic SNRs in the models such as 'H SNR', 'I SNR', 'S SN', and 'S SNR', in which the expected flux density of the radio emission from the aged USSNRs are approximately 0.1 mJy. This poses a challenge to detection and is consistent with the current non-detection of the SNR hosting a DNS binary in our Galaxy. On the other hand, in all of our models the USSNR diameter is in the order of 10 pc, which is also typical of the observed Galactic SNRs (Pavlović et al. 2013). We suggest that a faint surface brightness combined with a diameter D ∼ 10 pc can be a characteristics of a USSNR, and might be useful diagnostics for searching SNRs hosting a DNS binary.
At last we comment on the role of the ISM state on the radiative characteristics of the USSNRs. Comparisons of the solid and dashed lines in Figure 9 and Figure 10 demonstrate that the surface brightness of 'HOT' is fainter than that of 'WARM' at the same age. This is attributed to the fact that 'HOT' has a larger diameter and sky-projected angular size than 'WARM' and that the luminosities of these two models are similar to each other. Our simulations of CSM formation (Section 3) assume that the models 'WARM' and 'HOT' have the same thermal pressure but different densities in the initial profiles; we have shown that the model with a lower initial density leads to a larger diameter of the USSNR. Thus we conclude that the ISM density plays a role in determining the physical scale of the USSNR, which also affects the surface brightness.  the SNR population. One is the observable lifespan of the SNR, t snr , defined here as the timescale in which the radio emission from the SNR can be detected. The other one is t sn , the time interval between subsequent SNe or the inverse of the SN rate in a galaxy. The number of active SNRs can then be estimated as t snr /t sn . As for USSNRs, Hijikawa et al. (2019) predicted the event rate of USSNe as 510.88 gal −1 Myr −1 in their feasible population synthesis model, leading to t sn ∼ 2 × 10 3 years 1 . For the SNR lifetime, t snr ∼ 100-10 5 years can be implied from our models depending on the DSA efficiencies and the spectral index of accelerated electrons. Hence, the expected number of active USSNRs can be derived as 0.002 − 20.
These estimations involve uncertainties from observational conditions (e.g., sensitivity) as well as the DSA parameters. Models with high DSA efficiencies or hard power-law index for the accelerated electrons (e.g., 'H SN' and 'H SNR') probably over-estimate the observable lifespan; typical shock acceleration efficiency constrained by SNR observations are usually found to be lower than those inferred from the observations of radio SNe (Lee et al. 2012). Moreover, it has been suggested that the spectral index of the accelerated particles in young SNRs can be modified and steepened by non-linear effects associated with magnetic field amplification in an efficient DSA and Alfvénic drift effect (Vink et al. 2006;Zirakashvili & Ptuskin 2008;, whereas in mature SNRs it tends to follow the prediction by the standard DSA (Reynoso & Walsh 2015). The former is more appropriate for the situation considered in the present work, since our simulations indicate that the blastwave dies out at a young age in our CSM model. From these arguments, we can refer 'I SNR' as our fiducial models for the evolution of a USSNR, which predicts an observable lifespan t snr ∼ 10 4 years. Then we can further constrain the expected number of the observable USSNRs to be ∼ 2. Since the detected number of the Galactic SNRs reaches ∼ 400 (Green 2019), the most probable fraction of USS-NRs is then at most ∼ 0.5 % of all active SNRs. We note however that the quantification of the observable lifespan of the USSNRs involves uncertainties and depends on the sensitivity of the detectors as well.
The expected number of active USSNRs in a galaxy, ∼ 2, poses a severe challenge on the search of USS-NRs. Radio observation facilities capable of deep surveys such as the Square Kilometre Arrays (SKA) are requisite to solving this difficulty. A Galactic SNR survey with a sensitivity ∼ 0.1 mJy is one of the solutions to search for USSNRs, as well as for eliminating the possible bias against faint SNRs. Another possible solution is to extend the search to other galaxies in the local group. SNRs producing radio emission brighter than ∼ 10 21 (d/Mpc) 2 (F lim /µJy) erg s −1 Hz −1 in the local group can be detected by making use of the deepest observation projects, where F lim is the maximum sensitivity of SKA (Braun et al. 2019). This sensitivity enables us to detect the radio emission from USSNRs (proved by the model 'H SN' and 'H SNR' in Figure 7). Assuming that the galaxies in the local group have the same proportion of USSNRs to all kinds of SNRs (∼ 0.5%), this attempt might offer an opportunity to discover USSNRs.

General implications for stripped-envelope SNRs
We have shown that the blastwave of a USSNR suddenly loses its punch by being blunted by the hot plasma. The lifetime of the blastwave is limited to 10 3 years, and the diameter is roughly a few 10 pc. The evolution of USSNRs is different from that elucidated classically. Generally, after the Sedov phase at ∼ 10 4 years, radiative cooling from the swept ISM drains the internal energy away from the system, leading to a fast deceleration of the blastwave. Through the pressure-driven snowplow phase and momentum-driven snowplow phase, the SNR merge with the surrounding ISM at t ∼ 5×10 5 years (Cioffi et al. 1988). On the other hand, the evolution of USSNRs is heavily influenced by the non-uniform CSM density distribution and the presence of a hot plasma in the vicinity of the ISM wall, both of which are attributed to the wind driven by the binary interaction. The binary interaction is a key physical process in the evolutionary behaviors of USSNRs that deviate from the classical picture of SNR evolution.
Besides USSNe, it is widely believed that strippedenvelope SNe (Type IIb, Ib, and Ic SNe) are explosions of a massive star involved in binary interaction (e.g., Yoon et al. 2010;Ouchi & Maeda 2017;Fang et al. 2019). It can be speculated that the evolution of SNRs originated from stripped-envelope SNe also deviates from the classical theory. Considering that some fraction of the observed SNe are classified as stripped-envelope SNe (Type IIb, Ib, and Ic SNe, Eldridge et al. 2013), it is natural that some of the confirmed SNRs in our Galaxy also come from a stripped-envelope SN origin. Previously, in terms of hydrodynamics, the effect of the wind bubble and its multi-dimensional behaviors on the subsequent SNR evolutions have been investigated by making use of simple models for stellar mass loading (Tenorio-Tagle et al. 1990, 1991Dwarkadas 2005Dwarkadas , 2007, but models for mass loss history based on detailed binary evolution calculations have not been incorporated. We thus suggest that such stripped-envelope SNRs should be modeled with the mass loss history of the progenitor binary taken into account for their surrounding CSM environments (e.g., Yasuda et al. 2021a,b).

Radio emission from the hot plasma region
The velocity of the RLO wind is high, reaching ∼ 1000 km s −1 . It is therefore possible that in the formation process of the hot plasma driven by the RLO wind, electron acceleration and magnetic field amplification can happen through the DSA mechanism. Such effects can contribute to the radio luminosity and surface brightness of the subsequent USSNRs, and thus an evaluation of this process is required. Our simulations show that the velocity of the RLO wind shock is V sh,RLO ∼ 200 km s −1 . If we consider a hot ISM state (T ism ∼ 10 6 K), the Mach number of the shockwave launched by the RLO wind is in the order of unity. This indicates that the contribution from the hot plasma to the total flux of the radio emission from USSNRs is negligible in a hot ISM. On the other hand, in a warm ISM (T ism ∼ 10 4 K), the Mach number is large enough to sustain DSA. The hot plasma can then be a potential emitter of synchrotron radiation. Assuming that the region is optically thin for synchrotron radiation, the radio luminosity is written as L ν ∼ 4π 2 R 3 j ν,syn , where R is the position of the RLO wind shock. Based on the formulae introduced in this study, the luminosity can be roughly estimated as L ν ∼ 10 21 R 3 . Comparing this magnitude to the models, we can see that this contribution from the hot plasma in a warm ISM is negligibly small with respect to the predicted luminosities of young USSNRs (t 1000 years), but can be comparable to or even brighter than those of older USSNRs (t 1000 years), especially for the 'I SNR', 'S SN', and 'S SNR' models. In addition, we note that at later ages, the hot plasma can experience a compression from the expanding remnant, and the radio emission contribution from the plasma can be boosted further by this compression. Although the primary purpose of our study is on the modeling of USSNRs, the above discussion further advocates the importance of taking into account the CSM environment formed by the pre-SN mass loss activity of the progenitor in the USSNR emission model.

Treatment of radiative cooling
Apart from the models presented so far, we have also performed extra simulations in which radiative cooling occurs in regions with a broader range of optical depths with τ < c/v to approximate the contribution of photon diffusion to the energy loss. While this approach overestimates the energy loss from radiative cooling, it is helpful nonetheless for assessing the robustness of our results. In these models with an enhanced energy loss, we found that the blastwave velocity is decreased by a few percent. This confirms that the impact of radiative cooling on the overall dynamics is small enough that it plays an insignificant role in the modeling of USSNRs.

Effects of non-linear diffusive shock acceleration
We have employed the simplified treatment of particle acceleration and magnetic field amplification. In our study, non-linear effects in DSA are not considered, and the contribution of the pressure from cosmic-rays and its feedback to the hydrodynamics are not included. These effects can soften the energy distribution of accelerated electrons and could decrease the luminosity of non-thermal emission, including X-rays and gamma-rays (e.g., Vink et al. 2006;. Our estimate of the USSNR population can thus be altered by including such effects (see Section 5.1). On the other hand, however, the dynamics of the USSNR blastwave is mainly determined by the distribution of the CSM. The lifetime of the blastwave is mainly limited by its interaction with the hot plasma in the vicinity of the ISM wall formed by the pre-SN mass loss. Thus, improving the treatment of the microphysics in shock acceleration plays a secondary role in the observable lifespan of a USSNR.

Parametrizations of e and E min
There are two major simplifications in the parametrization for particle acceleration adopted in our study. First, some particle-in-cell simulations imply that the decrease of the Mach number (or the blastwave velocity) leads to a drop of the acceleration efficiency of protons (Caprioli & Spitkovsky 2014;Ha et al. 2018). This suggests the possibility that the acceleration efficiency of electrons also declines with a decreasing Mach number, while our study fixes e at a constant value with time. Second, it is believed that electrons with momentum greater than ∼ √ m e m p V b follow a power-law distribution even below the relativistic regime. However, we have fixed the minimum energy of the power-law distri-bution at E min in Equation (7) (see also Sironi & Giannios 2013). Hence, a decrease of the blastwave velocity leads to an increase of the number of electrons with a momentum p mom within √ m e m p V b p mom E min /c. This effect is not included in our models. In summary, our study is over-estimating the radio luminosities and the actual brightness of the USSNRs could be fainter if the above two factors are accounted for. However, the blastwave velocity in our calculations is in the order of ∼ 10 9 cm s −1 and the Mach number is sufficiently high in the young phase before the collision with the hot plasma. In the late phase (t 1000 years), the blastwave dies away rapidly. Therefore, the system considered in our study is not prone to the situation described above. Moreover, even if we include the two effects mentioned above in our modeling, the resulted radio luminosities should be fainter than those reported in Section 4.2, so that our conclusions on the characteristics and populations of USSNRs would not be affected qualitatively. Furthermore, we have examined two values for e shown in Table 2, and believe that the effect of the microphysics noted here can be investigated within this parameter space.

Asphericity
An aspherical configuration of the CE component and its effect on the wind hydrodynamics can be important as well. For instance it has been suggested that the material released by the CE ejection tends to distribute along the equatorial plane (Iaconi et al. 2019). Thus if the CE component resides in the vicinity of the SN progenitor it could affect the subsequent wind hydrodynamics. The gas ejected through the CE interaction should concentrate on the equatorial plane of the binary, while in the polar direction a static ISM should dominate. Then, the propagation of the wind driven by the RLO in the direction of the equatorial plane and the polar axis are regulated by the interaction between the ISM with and without the CE component, respectively. Our simulations in Section 3 show that the effect of the presence of the CE component is not significant regardless of the state of the ISM. From this point of view, by assuming a spherically blown wind from the progenitor binary, we can qualitatively speculate that the effect of possible non-spherical CE distributions would not be important.
Besides, an anisotropy of the conformation of the wind can be expected to shape the non-spherical geometry of the CSM as proposed in the literature of Type IIn SNe (Patat et al. 2011;Katsuda et al. 2016;Kumar et al. 2019). It is worth investigating the multi-dimensional structures of the composed CSM taking into account the anisotropy of the circumstellar environment and the wind outflows. These aspherical configurations of the CSM can alter the properties of the radiation from the SNe or SNRs, which will be examined in detail in a future work (see also e.g., Kurfürst & Krtička 2019;Suzuki et al. 2019).

SUMMARY
In this paper, we have investigated the characteristics of a SNR hosting a DNS binary, which we have termed a USSNR, using a grid of numerical models. A USSN has been proposed to be a transient event preceding the formation of a DNS binary. Before the USSN, the He star envelope is stripped away by the companion neutron star and escapes the binary system. By employing the mass-transfer history presented by Tauris et al. (2013), we simulated the hydrodynamics of the wind expelled from the progenitor binary, and constructed the largescale CSM structure around the USSN progenitor up to ∼ 100 pc. A hot plasma is formed in the vicinity of the ISM wall, which is found to play a critical role in governing the lifetime of the blastwave of the USSNR.
We also examined the dynamical and radiative evolution of a USSNR by considering a progenitor surrounded by the CSM composed by our simulation. We found that within the first ∼ 1000 years the blastwave traces the inner part of the CSM, producing a radio emission bright enough to be detected if the USSNR inhabits inside our Galaxy, though it is still fainter than those from typical SNRs. Once the blastwave collides with the hot plasma, it stalls rapidly and the radio luminosity also starts to decrease steadily. This dynamical behavior does not depend much on the strength of the CE ejection before the release of the helium gas from the progenitor binary. The surface brightness of the USSNR tends to be fainter than those of typical SNRs, while the diameter settles at D ∼ O(10 pc) similarly to the Galactic SNRs. Therefore, the USSNRs populate in the lower portion on the Σ-D diagram compared to the observed Galactic SNRs, and this can serve as a useful diagnostics for the search of a USSNR. We also confirmed that the initial ISM profile with a lower density allows the USSNR to expand further, leading to a lower surface brightness and a larger diameter. Furthermore, we evaluated the observable lifespan of a USSNR to be ∼ 10 4 years, defined as the time interval from the explosion to the point when the radio luminosity has declined beyond the detection limit of the present radio surveys. Combining the short observable lifespan of the USSNRs with the small event rate of USSNe, we conclude that the expected number of active USSNRs is less than one out of the observed 10 2−3 SNRs, which is consistent with the current non-detection of a SNR hosting a DNS. ACKNOWLEDGMENTS The authors thank the anonymous referee for his or her fruitful suggestions, and Haruo Yasuda and Norita Kawanaka for their comments which have deepened the discussion in our study. T.M. acknowledges the support from the Iwadare Scholarship Foundation in the fiscal year 2020 and from Japanese Society for the Promotion of Science ( In the simulation of the CSM formation, the value of ρ CE must be specified to determine the initial density profile. We consider a CE component with a total mass M CE = 10M ejected into a static uniform ISM. The required condition is where R ∞ = 3×10 21 cm is the outermost radius of the simulation domain. For the case ρ(r) = ρ CE exp(−r/R CE )+ρ ism , this can be analytically integrated, so that ρ CE M CE 8πR 3 CE = 7.96 × 10 −22 gcm −3 (A2) can be derived.

B. TESTS FOR THE NUMERICAL CODE
The numerical simulation code for the hydrodynamics employed in this study is verified in this section. Figure 11 shows the result of the shock tube problem with an adiabatic index γ = 5/3. At t = 0, a static (v = 0) gas is put into the simulation box, with a step function profile for its density and pressure centered at x = 0 as follows; ρ L = 1.0, v L = 0.0, p L = 1.0, ρ R = 0.125, v R = 0.0, p R = 0.1 (the subscripts L and R denote x < 0 and x ≥ 0, respectively). The numerical solution successfully reproduces the profiles given by the exact solution. Furthermore, Figure 12 displays the Sedov solution at t = 1.0 second in which an explosion energy E sedov = 1 erg is deposited into a uniform medium with ρ sedov = 1.0 × 10 −24 g cm −3 (Sedov 1959). The results are again in a good agreement with the analytical solutions for the density, velocity, and pressure profiles, as well as a good match of the shock radius given by R = 1.15(E sedov t 2 /ρ sedov ) 0.2 . These two experiments assure us a good accuracy of our numerical code.