A Dynamical Model for IRAS 00500+6713: The Remnant of a Type Iax Supernova SN 1181 Hosting a Double Degenerate Merger Product WD J005311

Figure 1 shows 0.5–5.0 keV X-ray images obtained with XMM EPIC-pn and Chandra ACIS. One can see a central source (≲20'') and a faint diffuse emission around it (≲100''). We exclude data in the lower energies, considering less accurate detector calibrations, and in the higher energies because of higher background rates. The central source is found to be slightly extended with the spatial resolution of Chandra (≈05). The diffuse emission extending to ∼100'' is not visible with Chandra because of the high detector background rate, which is ≈10 times higher than the expected diffuse flux. In order to evaluate the spatial distribution of the source, we extract radial profiles around the central source. The central position for extraction is set to (R.A., decl.)₂₀₀₀ = (132967, 675007) by referring to the SIMBAD catalog. ¹³ We take into account spatially dependent exposures (depending on off-axis angles), bad pixels, and gaps between the CCDs. For XMM EPIC, we obtain radial flux profiles for each instrument for each observation. For Chandra ACIS, we merge all the observations and extract radial profiles considering almost the same pointing directions and observation dates close to each other. In Figure 2, we show flux profiles extracted from one XMM observation in 2019 (ObsID 0841640101) and the merged Chandra data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare the observed distributions with a phenomenological model. The instrumental point spread functions (PSFs), i.e., the radial profiles of an ideal point source on the CCDs, are calculated for the source position with the SAS tool psfgen and CIAO tool simulate_psf. We assume a simple emission model for the central and diffuse sources: a sphere with a uniform and isotropic emissivity over the whole volume. The model for the radial flux profile (F(r), where r is the radius) is described as

$$\begin{eqnarray}&&F(r)=P(r,{\theta }_{\mathrm{in}})+D(r,{\theta }_{\mathrm{out}})+B,\end{eqnarray} \tag{ 1 }$$

with P(r, θ_in), D(r, θ_out), and B is the fluxes of the central source, diffuse source, and background (assumed to be uniform), respectively. The parameters θ_in and θ_out indicate the spatial extents of the central and diffuse sources, respectively. P(r, θ_in) and D(r, θ_out) are spherical emission convolved with the instrumental PSFs. We mainly focus on the evaluation of θ_in and θ_out.

We use the Markov Chain Monte Carlo algorithm ¹⁴ to determine the best-fit model parameters and their confidence ranges. To perform a precise evaluation, the profiles extracted from all the instruments in all the observations (21 in total) are modeled simultaneously via linked θ_in and θ_out for the XMM data. The normalizations of P(r, θ_in) and D(r, θ_out), and B for individual instruments are treated as free parameters. Radial ranges of 0–300 and 0''–10'' are used for the modeling of the XMM and Chandra data, respectively. As a result, we obtain θ_in < 15 and ${\theta }_{\mathrm{out}}\,=\,{131}_{-2-2}^{+2}$'' ¹⁵ with the XMM data. We obtain the corresponding $-2\mathrm{log}L/\mathrm{dof}\,=\,2597/2026$, where L and dof are the likelihood value and degree of freedom, respectively. This likelihood value is marginally acceptable, considering that we are using a simple phenomenological model. The central source is consistent with a pointlike source for XMM. The spatial extent of the diffuse source is similar to that in IR (Gvaramadze et al. 2019), as pointed out by Oskinova et al. (2020). On the other hand, we are able to constrain the radius of the central source as ${\theta }_{\mathrm{in}}\,=\,{1.52}_{-0.08-0.67}^{+0.08}$'' ¹⁶ with Chandra. The corresponding $-2\mathrm{log}L/\mathrm{dof}\,=\,29.0/9$, is also marginally acceptable. The best-fit models are overlaid in Figure 2. We also search for a possible time variation of the spatial distributions by separately modeling the XMM flux profiles in 2019 and 2021. We find results consistent with each other, θ_in < 15 and θ_out = 132'' ± 2'' in 2019 and θ_in < 15 and θ_out = 129'' ± 3'' in 2021, which means that there is almost no time variation in the sizes of the X-ray emitting regions. Throughout this paper, we adopt the value ${\theta }_{\mathrm{out}}={131}_{-2-2}^{+2}$'' obtained using all the observations simultaneously.

The background rates B determined from the XMM data, e.g., ≈0.0023 counts s⁻¹ arcmin⁻² (0.5–5.0 keV) for MOS1 in 2021, are consistent with or slightly higher than the pure particle background rate (Kuntz & Snowden 2008). This is reasonable considering the contribution of the sky background. The diffuse plus background flux D(r, θ_out) + B for the Chandra data are also consistent with the sum of the diffuse flux determined from the XMM data and the particle background rate (Bartalucci et al. 2014; Suzuki et al. 2021). Using the observed constraints on ${\theta }_{\mathrm{in}}={1.52}_{-0.08-0.67}^{+0.08}$'' with 1σ uncertainties, the minimum and maximum values of θ_in are 077 and 160, respectively. Adopting the Gaia distance of 2.3 kpc, the minimum and maximum physical radii are 0.0087 and 0.018 pc. These values are used to compare with our model. The spread of the inner source r_in = 0.0087–0.018 pc is much larger than the photospheric radius ∼1 × 10¹⁰ cm (Gvaramadze et al. 2019). This indicates that the inner shocked region has certainly spread, likely due to the wind from the central WD. However, since the wind-crossing time r_in/v_w ≲ 1 yr is much shorter than the SN age, this region may have been formed very recently.

We here extract spectra and estimate the temperatures, metal abundances, ionization timescale, and EMs of the central and diffuse sources. We use the XMM data, which have more than twice greater statistics than the Chandra data. Spectral extraction regions for the EPIC data are selected based on the imaging analysis results (Figure 2) and are shown in Figure 1. A circle with a 20'' radius is used for the central source, and an annulus with inner and outer radii of 30'' and 100'', respectively, is assumed for the diffuse source. As the background, we use the diffuse source and an annulus region with inner and outer radii of 140'' and 200'', respectively. The contamination of the diffuse emission to the background region is <1% of the total flux, according to the results of our imaging analysis. The central position for extraction is the same as that used for the extraction of the radial profiles. Since the RGS spectrum has much poorer statistics (<1% of the EPIC data), we only use it as a consistency check for the spectral analysis of the bright central source. We merge all the spectra in each year and for each instrument: the first- and second-order spectra are also combined for the RGS data. The spatial variation of the detector background is negligible, considering the on-axis position and sizes of our analysis regions (Kuntz & Snowden 2008).

As the spectral model, we consider optically thin thermal plasmas in two different assumptions of metalicity by referring to Oskinova et al. (2020): near solar and pure metal. Under each assumption, we try several models with different free parameters and finally adopt the simplest model that explains the data sufficiently well. Note that we do not pre-assume the abundance patterns reported in other wavelengths. We find an absorbed ionizing plasma model (tbabs × nei in XSPEC) explains both central and diffuse sources. The treatment of the model parameters for the central source is as follows: in the near-solar case, the hydrogen column density (N_H), electron temperature (k_B T_e), ionization timescale (τ), EM, and abundances of O, Mg, and Fe (=Ni) are treated as free parameters. The other metal abundances are fixed to solar. Here, k_B denotes the Boltzmann constant. In the pure-metal case, the same model but without H, He, N, Ar, and Ca is used. ¹⁷ For the diffuse source, the free parameters are the same as above but the abundance of O is fixed to solar and Ne is set free. The spectra obtained with all the instruments in 2019 and 2021 (six in total) are modeled simultaneously. We use the energy ranges of 0.5–5.0 keV and 0.5–2.0 keV for the spectral fit of the central and diffuse sources, respectively, considering their statistics.

The spectra and best-fit models are presented in Figure 3. Table 2 shows the best-fit parameters. Both models explain well both central and diffuse sources with acceptable fit statistics. The EMs (hereafter EM_in and EM_out for the central and diffuse sources, respectively) derived for the pure-metal and near-solar plasmas differ by 2 orders of magnitude because of the difference in the total ion emissivity. The abundance ratios are C/O = 0.014 ± 0.003, Ne/O = 0.005 ± 0.001, Mg/O = 0.0010 ± 0.0004, Fe/O = 0.005 ± 0.002 for the central source, and C/O = 0.43 (fixed), Ne/O = 0.41 ± 0.07, Mg/O = 0.04 ± 0.01, Fe/O < 0.003 for the diffuse source (all in number). The mass fractions in the pure-metal case are presented in Table 3. Oskinova et al. (2020) applied a two-temperature, CIE thermal plasma model to these sources. Their assumption of the abundance pattern is partly based on the IR observations (Gvaramadze et al. 2019) and thus different from ours. Although the mass fractions in this work differ from previous studies (Gvaramadze et al. 2019; Oskinova et al. 2020; Lykou et al. 2023), the tendency is similar: C/O < 1 and Ne/O < 1 for the central source and C/O∼1 and Ne/O∼1 for the diffuse source. Compared to Oskinova et al. (2020), the electron temperatures we obtain here lie between their higher and lower temperatures. The absorption column densities are similar.

Zoom In Zoom Out Reset image size

Table 3. Mass Fractions of the Central and Diffuse Sources Derived in the Pure-metal Case

Element	Central	Diffuse
C	${0.008}_{-0.008}^{+0.072}$	0.18
O	${0.97}_{-0.24}^{+0.03}$	0.56
Ne	${0.005}_{-0.005}^{+0.045}$	0.16 ± 0.03
Mg	${0.004}_{-0.004}^{+0.037}$	0.09 ± 0.02
Fe	${0.010}_{-0.010}^{+0.092}$	0.004 ± 0.004

Download table as:  ASCIITypeset image

Note that we use the data with approximately twice better statistics than those used in Oskinova et al. (2020). Oskinova et al. (2020) suggested that an additional nonthermal component might contribute to the central region. If we add an additional nonthermal component described with a power law with the fixed index of 2.4 following Oskinova et al. (2020), we indeed obtain an improved fit with C-statistic/dof = 4294/5408. However, since the fit statistics before adding a nonthermal component are already acceptable, we cannot claim that the central source requires a nonthermal component. We also search for possible spectral changes as a function of time but find no such evidence.

As we have shown in Section 2, the angular radius of the outer X-ray nebula is ${\theta }_{\mathrm{out}}={131}_{-2-2}^{+2}{^{\prime\prime}}$ based on the XMM-Newton observation. In our model, this corresponds to the forward shock radius r_out,f at t = 840 yr:

$$\begin{eqnarray}&&{r}_{\mathrm{out},{\rm{f}}}={1.46}_{-0.04}^{+0.02}\,\mathrm{pc}.\end{eqnarray} \tag{ 10 }$$

Note that the estimated error of r_out,f is smaller than the errors of the other parameters, thus we neglect the error hereafter. Applying this condition to the self-similar solution described above, we obtain the condition

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{E}_{\mathrm{ej}}}{{10}^{48}\,\mathrm{erg}}\right)}^{-1/2}{\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{0.5\,{M}_{\odot }}\right)}^{1/6}{\left(\displaystyle \frac{{\mu }_{\mathrm{ISM}}{n}_{\mathrm{ISM}}}{0.1\,{\mathrm{cm}}^{-3}}\right)}^{1/3}=0.86.\end{eqnarray} \tag{ 11 }$$

The EM of the outer X-ray nebula can be described as

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{EM}}_{\mathrm{out}}={n}_{{\rm{e}}}{n}_{\mathrm{ion}}V\\ \,=\displaystyle \frac{4\pi }{{m}_{{\rm{u}}}^{2}}\left(\displaystyle \frac{1}{{\mu }_{{\rm{r}}}^{2}}{\int }_{{r}_{\mathrm{out},{\rm{r}}}}^{{r}_{\mathrm{out},{\rm{c}}}}{r}^{2}\rho {\left(r\right)}^{2}{dr}+\displaystyle \frac{1}{{\mu }_{{\rm{f}}}^{2}}{\int }_{{r}_{\mathrm{out},{\rm{c}}}}^{{r}_{\mathrm{out},{\rm{f}}}}{r}^{2}\rho {\left(r\right)}^{2}{dr}\right).\end{array}\end{eqnarray} \tag{ 12 }$$

The first term in the right-hand side of Equation (12) represents the contribution from the reverse shock or the shocked ejecta; μ_r is the mean molecular weight and ρ(r) is the density. The second term represents the contribution from the forward shock, or the shocked ISM, where μ_f is the mean molecular weight. Using the shock structure obtained from the self-similar solution, the integrals of Equation (12) can be calculated as

$$\begin{eqnarray}&&{\int }_{{r}_{\mathrm{out},{\rm{c}}}}^{{r}_{\mathrm{out},{\rm{f}}}}{r}^{2}{\rho }^{2}{dr}=2.15\times {10}^{6}\,{{\rm{g}}}^{2}\ {\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{\mu }_{\mathrm{ISM}}{n}_{\mathrm{ISM}}}{0.1\,{\mathrm{cm}}^{-3}}\right)}^{2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\int }_{{r}_{\mathrm{out},{\rm{r}}}}^{{r}_{\mathrm{out},{\rm{c}}}}{r}^{2}{\rho }^{2}{dr}=8.96\times {10}^{6}\,{{\rm{g}}}^{2}\ {\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{\mu }_{\mathrm{ISM}}{n}_{\mathrm{ISM}}}{0.1\,{\mathrm{cm}}^{-3}}\right)}^{2}.\end{eqnarray} \tag{ 14 }$$

Note that they do not depend on either E_ej or M_ej for a given r_out,f. The mean molecular weights, μ_{r and} μ_f, are evaluated as follows. When the composition is pure metal, which is suggested by our analysis of the observed X-ray spectrum in Section 2, the mean molecular weights can be expressed as follows: in the outer regions

$$\begin{eqnarray}&&{\mu }_{\mathrm{out},\mathrm{pm}}^{2}={\mu }_{\mathrm{out},{\rm{e}},\mathrm{pm}}{\mu }_{\mathrm{out},\mathrm{ion},\mathrm{pm}},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{\mathrm{out},{\rm{e}},\mathrm{pm}}}=\displaystyle \sum _{i}\displaystyle \frac{\langle n{\rangle }_{i}{X}_{i}}{{A}_{i}}=0.47,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{\mathrm{out},\mathrm{ion},\mathrm{pm}}}=\displaystyle \sum _{i}\displaystyle \frac{{X}_{i}}{{A}_{i}}=0.060.\end{eqnarray} \tag{ 17 }$$

Here, the subscript "pm" represents the physical quantities for the case of pure-metal abundance, X_i denotes the mass fraction of element i, A_i its mass number, and 〈n〉_i is the average number of free electrons originating from element i. We adopt these values shown in Table 3. The number of free electrons is obtained from AtomDB 3.0.9 (Foster et al. 2012) using the abundance ratios and the ionization degree for each element. We have used k_B T_e and τ obtained from our X-ray analysis to estimate the degree of ionization. On the other hand, the mean molecular weight for the solar abundance is as follows:

$$\begin{eqnarray}&&{\mu }_{\mathrm{out},\mathrm{solar}}^{2}={\mu }_{\mathrm{out},{\rm{e}},\mathrm{solar}}{\mu }_{\mathrm{out},\mathrm{ion},\mathrm{solar}},\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{\mathrm{out},{\rm{e}},\mathrm{solar}}}\approx X\times \displaystyle \frac{1}{1}+Y\times \displaystyle \frac{2}{4}=0.85,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{\mathrm{out},\mathrm{ion},\mathrm{solar}}}\approx X\times \displaystyle \frac{1}{1}+Y\times \displaystyle \frac{1}{4}=0.78.\end{eqnarray} \tag{ 20 }$$

Here, X and Y are the mass fractions of hydrogen and helium, respectively.

The EM obtained from the X-ray analysis assumes a one-zone uniform shocked gas, while our EM model is a two-zone with different abundances. Though this makes a straightforward comparison difficult, the abundance of the mixture of the two regions is expected to be close to the ISM composition in the above self-similar solution because the mass ratio of the shocked ejecta to the shocked ISM is determined to be 0.28 (Chevalier 1982). Thus, we compare the calculated EM (Equation (12), adopting μ_r = μ_out,pm, μ_f = μ_out,solar) with the value obtained by the X-ray analysis assuming the near-solar abundance (EM_out = 3.7 ± 1.0 × 10⁵⁴ cm⁻³, Table 4). By substituting Equations (13)–(20) into Equation (12), we can determine the ISM number density as

$$\begin{eqnarray}&&{n}_{\mathrm{ISM}}=0.11\,{\mathrm{cm}}^{-3}\,{\left(\displaystyle \frac{{\mathrm{EM}}_{\mathrm{out}}}{3.7\times {10}^{54}\,{\mathrm{cm}}^{-3}}\right)}^{1/2}.\end{eqnarray} \tag{ 21 }$$

Here, we have used μ_ISM in Equation (7). Substituting this into Equation (11), we obtain the condition

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{E}_{\mathrm{ej}}}{{10}^{48}\,\mathrm{erg}}\right)}^{-1/2}{\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{0.5\,{M}_{\odot }}\right)}^{1/6}=0.97{\left(\displaystyle \frac{{\mathrm{EM}}_{\mathrm{out}}}{3.7\times {10}^{54}\,{\mathrm{cm}}^{-3}}\right)}^{-1/6}.\end{eqnarray} \tag{ 22 }$$

Figure 6 shows the obtained relation between E_ej and M_ej. Given that IRAS 00500+6713 hosts a WD remnant, the ejecta mass would not be larger than M_ej ≲ 1 M_⊙. In this case, the ejecta kinetic energy is constrained to be E_ej ≲ 10⁴⁸ erg only from the outer X-ray nebula. These results are broadly consistent with the previous work (Oskinova et al. 2020).

Zoom In Zoom Out Reset image size

As we show in Section 2, the angular radius of the inner X-ray nebula is constrained as $0\buildrel{\prime\prime}\over{.} 77\lt {\theta }_{\mathrm{in}}\lt 1\buildrel{\prime\prime}\over{.} 6$ based on the Chandra observation in 2021. In our model, this corresponds to the thin-shell radius, i.e., the radius of the wind termination shock r_sh at t = 840 yr:

$$\begin{eqnarray}&&0.0087\,\mathrm{pc}\lt {r}_{\mathrm{sh}}\lt 0.018\,\mathrm{pc}.\end{eqnarray} \tag{ 26 }$$

Figure 7 shows the time evolution of the inner nebula radius calculated for several parameter sets. One can see that r_sh at t = 840 yr becomes larger for a smaller M_ej, a larger ${\dot{M}}_{{\rm{w}}}$, and/or smaller t_w. This is simply because the wind keeps pushing the fallback ejecta since its launch. Therefore, too small M_ej or too small t_w leads to a too-expanded shell, which is inconsistent with the Chandra observation. This can give a constraint on the parameter set $({\dot{M}}_{{\rm{w}}},{t}_{{\rm{w}}},{M}_{\mathrm{ej}})$.

Zoom In Zoom Out Reset image size

By modeling the expansion of the inner nebula, we can also estimate the total mass in the shocked region. Figure 8 shows the time evolution of the shocked mass for a case with $({\dot{M}}_{{\rm{w}}},{t}_{{\rm{w}}},{M}_{\mathrm{ej}})=(1\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1},820\,\mathrm{yr},0.3\,{M}_{\odot })$. The solid and dashed lines indicate the contribution from the shocked wind swept by the reverse shock and the shocked-SN ejecta swept by the forward shock, respectively. One can see that the latter is 2–3 orders of magnitude larger than the former. Importantly, however, the shocked-SN ejecta component will not contribute to the EM of the inner X-ray nebula; as we show in the Appendix, the electron temperature is kept below 3–4 eV in the shocked-SN ejecta. Thus, we only consider the shocked wind as the source of the inner X-ray nebula. In this case, the EM of the inner nebula can be estimated as

$$\begin{eqnarray}&&{\mathrm{EM}}_{\mathrm{in}}={n}_{{\rm{e}}}{n}_{\mathrm{ion}}{V}_{\mathrm{sh},{\rm{w}}}=\displaystyle \frac{{M}_{\mathrm{sh},{\rm{w}}}^{2}}{{\mu }_{{\rm{e}}}{\mu }_{\mathrm{ion}}{m}_{{\rm{u}}}^{2}{V}_{\mathrm{sh},{\rm{w}}}},\end{eqnarray} \tag{ 27 }$$

where ${M}_{\mathrm{sh},{\rm{w}}}={\dot{M}}_{{\rm{w}}}(t-{t}_{{\rm{w}}})$ is the total mass of the shocked wind, and

$$\begin{eqnarray}&&{V}_{\mathrm{sh},{\rm{w}}}\approx \displaystyle \frac{4}{3}\pi \left({r}_{\mathrm{sh}}^{3}-{r}_{\mathrm{in},{\rm{r}}}^{3}\right)=\displaystyle \frac{4}{3}\pi {r}_{\mathrm{sh}}^{3}\left[1-{\left(\displaystyle \frac{2}{1+\gamma }\right)}^{3}\right]\end{eqnarray} \tag{ 28 }$$

is the volume of the shocked-wind region. Here, we estimate the positions of the reverse shock fronts as

$$\begin{eqnarray}&&{r}_{\mathrm{in},{\rm{r}}}\approx (1-{ \mathcal C }){r}_{\mathrm{sh}}=\displaystyle \frac{2}{1+\gamma }{r}_{\mathrm{sh}}.\end{eqnarray} \tag{ 29 }$$

With an approximation that the shock region is plane-parallel and adiabatic, the compression ratio is ${ \mathcal C }=(\gamma -1)/(\gamma +1)=4$ for γ = 5/3 in the strong shock limit. Since we consider the inner X-ray region to be formed by the interaction between the C-O-rich wind from the central WD and the unshocked ejecta, we consider the abundance in this region to be pure metal, and adopt the value of EM_in obtained from the analysis under the pure-metal assumption. The mean molecular weights in the shocked-wind region are calculated for pure-metal abundance with ionization inferred from the X-ray spectroscopy in Section 2:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{\mathrm{in},{\rm{e}},\mathrm{pm}}}=\displaystyle \sum _{i}\displaystyle \frac{\langle n{\rangle }_{i}{X}_{i}}{{A}_{i}}=0.49,\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{\mathrm{in},\mathrm{ion},\mathrm{pm}}}=\displaystyle \sum _{i}\displaystyle \frac{{X}_{i}}{{A}_{i}}\sim 0.062.\end{eqnarray} \tag{ 31 }$$

Zoom In Zoom Out Reset image size

By comparing the calculated EM_in of the shocked wind (Equation (27)) with the observed value for the pure-metal model (Tables 2 and 4), we can obtain an additional constraint on our model parameters. Here, as shown in Table 4, the value of EM_in is constrained from the X-ray analysis as follows:

$$\begin{eqnarray}&&5.0\times {10}^{52}\,{\mathrm{cm}}^{-3}\lt {\mathrm{EM}}_{\mathrm{in}}\lt 7.6\times {10}^{52}\,{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 32 }$$

Figure 9 shows the time evolution of EM_in calculated for several parameter sets $({\dot{M}}_{{\rm{w}}},{t}_{{\rm{w}}},{M}_{\mathrm{ej}})$ =(2 × 10⁻⁶ M_⊙ yr⁻¹, 820 yr, 0.6 M_⊙), (2 × 10⁻⁶ M_⊙ yr⁻¹, 820 yr, 0.3 M_⊙), and (4 × 10⁻⁶ M_⊙ yr⁻¹, 820 yr, 0.6 M_⊙). As expected, a larger ${\dot{M}}_{{\rm{w}}}$ gives a larger shocked mass and thus a larger EM_in. In addition, a larger M_ej also gives a larger EM_in; this is because the inner nebula expansion becomes slower for a larger M_ej and thus the volume of the inner nebula V_sh,w becomes smaller (see Equation (27)).

Zoom In Zoom Out Reset image size

In addition to measuring angular radius and EMs, our X-ray spectral analysis has determined relevant quantities of the X-ray-emitting plasma, such as ionization timescale and electron temperature (see Table 2). The ionization timescale is estimated as τ ≈ n_e × Δt, where n_e is the post-shock electron number density and Δt is the dynamical time of the plasma, and it provides an indicator of the degree of ionization.

In the allowed parameter space obtained by our analysis, the ionization timescale in the inner shocked region can be estimated as τ_in ≈ ρ_sh,w/(μ_e m_u) × Δt_w ∼ 10^10–11 s cm⁻³. Here, ρ_sh,w is the density of the inner shocked region obtained by ρ_sh,w = M_sh,w/V_sh,w. The ionization timescale in the outer shocked region has a value of τ_out ≈ 4n_ISM × 840 yr ∼ 10¹⁰ s cm⁻³. Both timescales are broadly consistent with the values obtained from the X-ray spectroscopy. In order to achieve CIE in the nebulae with the estimated temperatures, τ ≳ 10¹¹ s cm⁻³ is required for oxygen and ≳10¹² s cm⁻³ for heavier elements (e.g., Smith & Hughes 2010). Hence, we infer that both the inner and outer nebulae are likely not in CIE. To estimate the electron temperature accurately, a detailed plasma calculation coupled with shock dynamics is required (e.g., Hamilton et al. 1983; Masai 1994; Tsuna et al. 2021a), which is beyond the scope of this paper. Still, if the shock downstream is adiabatic and the electrons and ions are far from CIE, which are good approximations, particularly for the shocked wind or the inner nebula, the electron temperature can be estimated as ${k}_{{\rm{B}}}{T}_{{\rm{e}},\mathrm{in}}\approx 511\,\mathrm{keV}\times {({v}_{{\rm{w}}}/c)}^{2}\sim 1.3\,\mathrm{keV}$, which is also consistent with X-ray spectroscopy. Here, c denotes the speed of light.

IR observations of IRAS 00500+6713 revealed a dust-rich ring with an angular radius of θ_d ∼ 1' or a physical size of r_d ∼ 0.66 pc. The ring is surrounded by a diffuse IR halo (Gvaramadze et al. 2019; Lykou et al. 2023), whose expansion velocity was assumed to be similar to that of the shocked ejecta v_Sii ∼ 1, 100 km s⁻¹ based on the [S ii] line signatures observed for ${\theta }_{\mathrm{SII}}\sim 90^{\prime\prime}$ or physical size of r_Sii ∼ 0.99 pc (Ritter et al. 2021; Lykou et al. 2023; Fesen et al. 2023).

In our model, we expect that there are relatively low-temperature layers in the unshocked region, r_sh < r < r_out,r, where the reverse shock radius is calculated as r_out,r = 1.04 pc. Note that the ratio between the forward shock radius and the reverse shock radius is r_out,f/r_out,r = 1.39, according to the self-similar solution with n = 6 and s = 0 (Chevalier 1982). Therefore, we have r_sh < r_d, r_SII < r_out,r, which is consistent with the observed sizes of the dust-rich ring and the diffuse IR halo. Since the inferred expansion velocity of the outer edge of the unshocked region is ∼r_out,r/840 yr ∼ 1200 km s⁻¹, which is comparable to v_Sii = 1100 km s⁻¹, the region with the [S ii] line signature should be close to the reverse shock front. On the other hand, the outer edge of the dust-rich ring, which is at a slightly smaller radius, may correspond to the sublimation front of the dust that is irradiated by the radiation from the reverse shock. Although a dynamical calculation including the formation and destruction of dust is beyond the scope of this paper, it is important, particularly in light of the recently identified radially aligned filamentary structure around this region (Fesen et al. 2023).