Type Ia supernova ejecta-donor interaction: explosion model comparison

The close-binary system initially consists of a MS donor star ($M_{\mathrm{2}}^{\mathrm{i}}=2.5\,M_{\sun}$) and a WD accretor ($M_{\mathrm{WD}}^{\mathrm{i}}=1\,M_{\sun}$), with an orbital period $P_{\mathrm{orb}}^{\mathrm{i}}=2\,\mathrm{days}$. This is a typical SN Ia progenitor system according to the population synthesis studies of Chen et al. (2014) and Claeys et al. (2014). To create an initial donor star model we use an one-dimensional (1D) stellar evolution code, MESA (version 4906; Paxton et al., 2011; Paxton et al., 2013, 2015, 2018). The WD is modelled as a point mass with retention efficiencies of H and He burning from Yaron et al. (2005) and Kato & Hachisu (2004), respectively. Orbital evolution, including gravitational-wave radiation and associated angular-momentum loss, is modelled as in Chen et al. (2014).

The binary system is evolved until the WD accretes enough material to reach the Chandrasehkar-mass limit ($\sim 1.4\,M_{\sun}$). We show in Fig. 1 the evolution of mass transfer rate and the HR diagram of the donor star. At the time of explosion, the binary separation is $A=6.44\,R_{\sun}$, while the donor radius and mass are $R_{2}=1.68\times 10^{11}\,\mathrm{cm}$ and $M_{\mathrm{2}}=1.17\,M_{\sun}$, respectively. The density structure of the donor star at this moment is shown in Fig. 2.

The W7 carbon deflagration model (Nomoto et al., 1984; Maeda et al., 2010) ²²2In particular, the W7 model of Maeda et al. (2010) is used in this work. and the N100 delayed-detonation model (Seitenzahl et al., 2013) are employed in this work. With both, a one-dimensional, angle-averaged shell model of the explosion $20\,s$ after ignition is interpolated onto a two-dimensional grid at the maximum refinement level of our hydrodynamical model to reduce the severity of numerical grid artifacts in the expanding gas front. Both explosions eject $1.4\,M_{\sun}$ but the N100 model has a higher kinetic energy in its ejecta ($1.40\times 10^{51}\,\mathrm{erg}$ compared to $1.23\times 10^{51}\,\mathrm{erg}$ in W7). Fig. 3 compares the density as a function of expansion velocity in N100 and W7. The N100 model also has a higher-velocity tail in the outer region of ejecta compared to W7 (see Fig. 3).

The simulation setup is similar to what has been described in Pan et al. (2010, 2012a) but upgraded to FLASH 4 version. To initialise the hydrodynamic simulation, the MS donor-star model is interpolated onto a two-dimensional grid from the one-dimensional shell model provided by MESA. The system is then relaxed for 10${}^{4}\,\mathrm{s}$ which artificially damps the momentum of the star so it settles into hydrostatic equilibrium. Next, the W7 or N100 SN Ia explosion is placed at a distance $A$ from the centre of the donor star on the $z$-axis. The simulation then runs for a further 10${}^{4}\,\mathrm{s}$ until the ejecta-donor interaction finishes and the mass of the donor star stabilises. Because the donor fills its Roche lobe at the moment of supernova, the binary separation ($A$) and the donor radius ($R_{2}$) follow Equation 1 (Eggleton, 1983),

where $q$ is the mass ratio of binary system at the moment of explosion.

In our simulations we employ 12 levels of adaptive-mesh refinement (AMR) in the SN and donor-star regions, and a maximum of 10 levels of refinement outside these regions. We refer to this setup as the 10/12 level. This is equivalent to a uniform grid spacing of $6.1\times 10^{7}\,\mathrm{cm}$ ($\sim 3.6\times 10^{-4}\,R_{2}$) inside the star and supernova regions, and a maximum of $9.8\times 10^{6}\,\mathrm{cm}$ ($\sim 5.8\times 10^{-3}\,R_{2}$) outside of these regions. The cylindrical coordinates $r,z$ lie in the ranges $r=(0,2\times 10^{12})$ and $z=(-2\times 10^{12},2\times 10^{12})$, which are about 12 times larger than the donor-star radius in each direction.

To determine the robustness of numerical results, Pan et al. (2010) performed a convergence test using different maximum AMR levels from 6/8 to 10/12 in their 2D impact simulations with FLASH. They found that the differences in the results, e.g. to the unbound mass and kick velocity, given by these different resolution levels are not significant, and that level 10/12 is sufficient to study ejecta-donor-star interaction. In addition, the convergence tests of Pan et al. (2010) were performed with a helium-star companion which is more compact than our main-sequence donor. The turbulence generated from ejecta-companion impact is weaker in the main-sequence channel and we expect better convergence with main-sequence companions. Because the FLASH code is used in this work, albeit a newer version, and our initial setup and assumptions are similar to those of Pan et al. (2010), we use the recommended 10/12 level in all our simulations.

We neglect orbital and rotational motion in our simulations because the SN Ia ejecta velocity ($\sim 10^{4}\,\mathrm{km\,s^{-1}}$) far exceeds the orbital velocity (a few $\times 10^{2}\,\mathrm{km\,s^{-1}}$), and the two stars are tidally locked at the time of explosion so the donor surface also rotates at the orbital velocity.

Fig. 4 shows the evolution of the ejecta-donor interaction in our N100 model. The exploding star is located vertically above the figure, out of shot, and the blast wave moves downwards towards the donor. The SN ejecta expands freely for a few minutes then hits the surface of the donor star. The impact strips some H-rich material from and increases the pressure on the surface of the donor. A high-pressure shock is driven into the stellar interior and a bow shock develops in the outer regions which diverts SN ejecta material around the star. As the diverted material further expands around the donor star, the interstellar medium (ISM) is compressed in the downstream region (below the star in the figure). The two opposing fronts collide behind the star at a few hundred seconds (Fig. 4a), forcing ISM material back towards the star, whose impact can be observed as a small notch at its rear. The shock continues to pass through the donor star, slowing as it moves through the higher density core (Fig. 4b). The opening hole is created behind the donor star.

Instabilities develop on the donor surface closest to the SN, disrupting its initially-spherical shape. At about $t=1700\,\mathrm{s}$ the shock passes through the entire star and pushes stellar material out of its rear as it exits (Fig. 4c). The symmetries of the donor star are completely destroyed as the shock passes through it. The stripped donor material mixes with the SN ejecta material in the opening hole behind the star (Fig. 4d). The outer layers of the star expand in the shape of an inverted teardrop, surrounded by Rayleigh-Taylor instabilities. Instabilities are also observed at the edge closest to the SN, travelling outwards from the centre of the donor but shielded from incoming material by a bow shock. The donor star expands because of the shock heating, and is totally out of hydrostatic and thermal equilibrium. By $t=6000\,\mathrm{s}$, the donor starts returns to hydrostatic equilibrium as the interaction finishes and a relatively spherical remnant remains (Fig. 4f). At the same time, an expanded low-density envelope has forms in the outer region of the star. In addition, the centre of the donor star moves away from its original position because it receives a velocity kick because of momentum imparted by the SN ejecta during the interaction.

A comparison of the interaction between the W7 and N100 explosion models at equal times after the supernova explosion is shown in Fig. 5. Both explosion models lead to similar interaction characteristics: a turbulent front behind the frontal bow shock, eddies to the sides of the rear of the star, and the teardrop-shaped plume ejected after the shock passes through the star in its entirety. The clearest difference is that the N100 explosion has faster ejecta, hence the interaction evolves more quickly than with W7 and strips more material from the donor star.

Fig. 6 shows the unbound donor mass as a function of time using the two explosion models in our simulations. To determine whether the donor material is unbound, we calculate its total energy by summing the kinetic, potential and internal energies, $E_{\mathrm{tot}}=E_{\mathrm{kin}}+E_{\mathrm{pot}}+E_{\mathrm{int}}$, and if $E_{\mathrm{tot}}>0$ the material is unbound. As shown in Fig. 6, the unbound donor mass stabilises about $5000\,\mathrm{s}$ after the explosion. The time-evolution of unbound mass is similar in both models, beginning with a steep rise as the shock front transfers momentum and thermal energy to the donor star, followed by a slight decrease during thermal relaxation which returns some material to a bound state. At the end of our simulations, we find that $0.293\,M_{\sun}$ and $0.368\,M_{\sun}$ of H-rich material, representing $25$ and $31\,\mathrm{per~{}cent}$ of the donor mass at the time of explosion, is stripped from the MS donor surface by the SN explosion using the W7 and N100 models respectively.

Extra mass stripping in our N100 model is caused by the larger kinetic energy of the N100 explosion, $1.40\times 10^{51}\,\mathrm{erg}$ (Seitenzahl et al., 2013), compared with W7 which has $1.23\times 10^{51}\,\mathrm{erg}$ (Nomoto et al., 1984). Previous studies have found that more powerful ejecta generally strip more donor mass for a given donor model and binary separation, and the total stripped mass is linearly proportional to the SN kinetic energy (Pakmor et al., 2008; Pan et al., 2012a; Liu et al., 2013a). The stripped mass in our W7 model increases to $0.334\,M_{\sun}$ if we scale it with a higher kinetic energy of $1.40\times 10^{51}\,\mathrm{erg}$. On the other hand, as shown in Fig. 3, the outer SN ejecta of the N100 model is more powerful than that of W7, hence N100 is also expected to strip more.

As described in Section 1, the interaction between SN Ia ejecta and a donor star has been investigated by previous studies with both analytical methods, and two- and three-dimensional hydrodynamical simulations (e.g., Wheeler et al., 1975; Marietta et al., 2000; Pakmor et al., 2008; Pan et al., 2012a; Liu et al., 2012; Boehner et al., 2017). Most studies found that $\gtrsim 0.1\,M_{\sun}$ can be stripped from a main-sequence donor star, from a red giant the whole envelope is stripped. In particular, Marietta et al. (2000) found that $15\,\mathrm{per~{}cent}$ of the donor star is stripped in their HCV model. The HCV model has a mass $M_{2}=1.017\,M_{\sun}$, radius $R_{2}=6.8\times 10^{10}\,\mathrm{cm}$ and ratio of binary separation to radius $A/R_{2}=3.0$ (Marietta et al., 2000). Compared to our main-sequence donor model, the HCV model is more compact and has a larger $A/R_{2}$ (Table 1). The stripped mass in HCV increases to $\sim 0.2\,M_{\sun}$ (Pan et al., 2012a, see their Fig. 12) if we reduce its $A/R_{2}$ to 2.67 similarly to our model. Liu et al. (2012) found between 6 and 24 per cent depending on the MS donor model (see also Pan et al., 2012a). More recently, Boehner et al. (2017) predicted $21-26\,\mathrm{per~{}cent}$ of a main-sequence donor is stripped in their 2D simulations using FLASH. Our results are comparable with these previous studies. As shown in Figs. 2 and 1, our donor star is near the end of its main sequence at the time of SN explosion, hence it has a dense core with an extended envelope. This explains why our stripped mass is a greater than in Boehner et al. (2017), in which the donor stars are normal MS stars and they have a more compact envelope than our donor model. However, our stripped mass greatly exceeds the $\sim 1$–$5\,\mathrm{per~{}cent}$ predictions of Pakmor et al. (2008, their Table 2). The differences are mainly caused by the detailed structure of the donor star models and binary separations at the moment of SN Ia explosion (Liu et al., 2012; Pan et al., 2012a; Boehner et al., 2017).

Compared to our model, the main-sequence donor models of Pakmor et al. (2008) are more compact and have longer binary separations (i.e., they have a larger $A/R_{2}$) so in their simulations less mass is stripped. For instance, model rp3_24a in Pakmor et al. (2008) has a binary separation of $4.39\times 10^{11}\,\mathrm{cm}$ which is similar to our $4.48\times 10^{11}\,\mathrm{cm}$. However, our companion radius exceeds theirs ($7.21\times 10^{10}\,\mathrm{cm}$) by a factor of $2.3$. Additionally, model rp3_20a in Pakmor et al. (2008) has the same mass as ours (i.e., $M_{2}$=$1.17\,M_{\odot}$), but their $A/R_{2}=4.5$ exceeds ours ($A/R_{2}=2.67$) by a factor of 1.7. If we assume that their rp3_20a model also has $A/R_{2}=2.67$, their stripped mass increases to $0.19\,M_{\odot}$ based on their Equation 4. The main differences between the main-sequence donor models of Pakmor et al. (2008) and ours are the different treatments of pre-SN mass-loss of the MS stars and the binary separations adopted at the moment of SN explosion, which have been discussed in detail in our previous work (Liu et al., 2012, their Section 4.1). Pakmor et al. (2008) start the mass accretion phase with less-evolved companion stars and they constructed their companion star models by rapidly removing mass while evolving a single main-sequence star to mimic pre-SN mass-transfer in a binary system. As a consequence, their models have a very large $A/R_{2}$, and they are too compact to fill their Roche-lobe at the moment of SN explosion. In contrast, our main-sequence donor model is constructed by consistently following mass-transfer in a binary system. Our main-sequence star has a more extended structure and fills its Roche-lobe when the SN explodes.

The donor star receives a kick as a result of momentum transfer from the ejecta during the interaction. The time evolution of kick velocity behaves similarly to the stripped mass (Fig. 6). At the end of simulations ($t=10^{4}\,\mathrm{s}$), the kick velocities are $73.8\,\mathrm{km\,s^{-1}}$ and $86.0\,\mathrm{km\,s^{-1}}$ in our W7 and N100 models respectively. The faster kick in the N100 model is probably because its ejecta have a higher kinetic energy than in W7. The movement of the donor away from its original position is seen in Fig. 4. The final spatial velocity of the surviving donor star is only slightly altered by the kick, with most of the velocity arising from its pre-explosion orbital motion.

Hydrogen-emission lines will be seen in late-time spectra of SNe Ia, about a few hundred days after the explosion, if sufficient hydrogen-rich material ($>0.001-0.058\,M_{\sun}$, which depends on the distance of the observed SNe Ia) is stripped from the donor surface during the interaction of the explosion with the star (Mattila et al., 2005; Lundqvist et al., 2015; Botyánszki et al., 2018; Dessart et al., 2020). Therefore, searching for hydrogen-emission lines in late-time spectra of SNe Ia should provide strong evidence in favour, or otherwise, of a SD progenitor system to observed Ia supernovae. To date, no hydrogen-emission lines have been detected in late-time spectra of normal SNe Ia (Leonard, 2007; Lundqvist et al., 2013, 2015; Maguire et al., 2016; Shappee et al., 2016; Tucker et al., 2020). However, hydrogen features have been detected in late-time spectra of two fast-declining, sub-luminous SNe Ia, SN 2018fhw (Kollmeier et al., 2019) and SN 2018cqj (Prieto et al., 2020). All this late-time observations limit the hydrogen-rich mass in type-Ia SN progenitor systems to $\leq 0.001-0.058\,M_{\sun}$.

This observational upper limit to the stripped mass is less than our simulations predict, $\geq 0.30\,M_{\sun}$, and our results are similar to previous studies (Marietta et al., 2000; Pakmor et al., 2008; Pan et al., 2012a; Liu et al., 2012; Liu & Stancliffe, 2017; Boehner et al., 2017). This suggests that the non-detection of hydrogen features in late-time spectra of SNe Ia is seriously challenges the SD scenario if we assume that our current hydrodynamical and radiative transfer models of late-time hydrogen emission are accurate.

The distribution of velocities of H/He-rich material stripped from the donor star hints at expected H- and He-emission line features in late-time spectra of SNe Ia (Mattila et al., 2005; Botyánszki et al., 2018; Dessart et al., 2020). In Fig. 7, we show post-impact velocity distributions of stripped material from the surface of our donor star in our two explosion models. The typical velocities of stripped material in our W7 and N100 models are $350\,\mathrm{km\,s^{-1}}$ and $500\,\mathrm{km\,s^{-1}}$, respectively, which are consistent with previous works (Marietta et al., 2000; Pan et al., 2012a; Liu et al., 2012).

Additionally, Marietta et al. (2000) suggests that the high-velocity tail of stripped material is an important diagnostic that allows discrimination between various donor models (see also Boehner et al., 2017). Here, we check whether the two explosion models investigated in this work can be discriminated similarly. Fig. 7 shows that although a high-velocity, $\gtrsim 1000\,\mathrm{km\,s^{-1}}$, tail is seen in both our explosions models, there is little difference between their velocity distributions. This suggests that it may be difficult to discriminate between the W7 and N100 models with the features caused by the high-velocity tail for a given main-sequence donor star model.

In the SD SNIa scenario, the donor star survives the explosion after ejecta-donor interaction. However, surviving donor stars are generally not expected in the double-degenerate scenario, although Shen et al. (2018) have recently suggested that a WD donor may survive in their D6 double-generate model. Therefore, searching for surviving donor stars in SN remnants (SNRs) is a promising way to distinguish between the SD and DD scenarios of SNe Ia (Ruiz-Lapuente et al., 2004; Ruiz-Lapuente, 2019; Kerzendorf et al., 2009, 2013; Schaefer & Pagnotta, 2012; Pan et al., 2014; Li et al., 2019).

To provide the predicted observational features of surviving donor stars for comparison with observations of SNRs, it is important to address long-term post-impact evolution of surviving donor stars based on the results of multi-dimensional hydrodynamical impact simulations. After the explosion, the surviving donor star becomes a runaway star. It escapes from the binary system with a spatial velocity caused mostly by its pre-explosion orbital motion. In Figure 8, we present 1D-averaged density (left panel) and entropy (right panel) profiles of our surviving main-sequence donor star at the end of our simulations with the W7 and N100 explosions. For comparison, the initial profiles of the donor star at the beginning of the simulations are also shown. The central density of the donor star decreases by about $50\,\mathrm{per~{}cent}$ because of the shock that passes through the star. In addition, the donor star is impacted and heated during the interaction, and expands by a factor of about two by the end of our simulations.

While the donor subsequently thermally relaxes, it is overluminous compared to its pre-explosion equilibrium state (Podsiadlowski, 2003; Pan et al., 2012b, 2013; Shappee et al., 2013; Liu & Zeng, 2020; Liu et al., 2021, 2022). For instance, Pan et al. (2012a) and Shappee et al. (2013) suggest that the post-impact main-sequence companion stars are significantly more luminous (up to about $10$–$1000\,L_{\odot}$) within a few thousand years, depending on the main-sequence companion model. Using the approach of Liu & Zeng (2020), the 2D surviving donor model computed from our impact simulations in this work can be mapped into a 1D stellar evolution code such as MESA to follow their subsequent post-impact evolution for more than a few hundred years to obtain expected observational features (Pan et al., 2013; Liu et al., 2021, 2022). This is important to understand how different explosion models affect the observational features of the surviving main-sequence companion stars, and is useful in the search for such events in SN Ia remnants. This investigation will be performed in a forthcoming study.

Our simulations are performed in 2D without considering the orbital motion and rotation of the binary system. This is expected to only slightly affect our results (Liu et al., 2012; Pan et al., 2012a) because these two velocities are slower than the typical ejecta velocity of $10^{4}\,\mathrm{km\,s^{-1}}$ by more than an order of magnitude. Pan et al. (2012a) have performed the ejecta-donor interaction in 2D and 3D for a given main-sequence donor model with the FLASH code. They find that the results of the ejecta-donor interaction are quite consistent in 2D and 3D if the orbital motion and spin of the star are not included. However, the orbital motion and spin can lead to slightly different results in 3D (Pan et al., 2012a, their Table 2). Future 3D simulations are still needed to make more strict theoretical predictions such as post-SN rotation features of the surviving donor star. For instance, exploring the post-impact spin features of the surviving donor star requires the impact simulation to be in 3D.

In addition, the rotation of an accreting WD is not considered when we construct the donor model at the moment of SN explosion in our 1D binary evolution calculation. The accretion from the donor star spins up the WD, leading to a delay to spin down the WD before it explodes (i.e., the so-called “spin up/spin down” model, Justham 2011; Di Stefano et al. 2011). Because the spin-down timescale is quite uncertain, the donor star in the MS donor scenario could be any type of star from MS (short delay) to WD (long delay) at the time of explosion. In such a case, the amount of stripped donor mass and donor properties will very likely differ considerably from our results with a MS donor.