Evolution of Main-sequence-like Surviving Companions in Type Ia Supernova Remnants

Recent theoretical and numerical studies of Type Ia supernovae (SNe Ia) explosions within the single-degenerate scenario suggest that the nondegenerate companions could survive during the supernova impact and could be detectable in nearby supernova remnants. However, observational efforts show less promising evidence of the existence of surviving companions from the standard single-degenerate channels. The spin-up/spin-down models are possible mechanisms to explain the nondetection of surviving companions. In these models, the spin-up phase could increase the critical mass for explosion, leading to a super-Chandrasekhar-mass explosion, and the spin-down phase could lead to extra mass loss and angular momentum redistribution. Since the spin-down timescale for the delayed explosion of a rotating white dwarf is unclear, in this paper we explore a vast parameter space of main-sequence-like surviving companions via two-dimensional hydrodynamic simulations of supernova impact and the subsequent stellar evolution of surviving companions. Tight universal relations to describe the mass-stripping effect, supernova kick, and depth of supernova heating are provided. Our results suggest that the not-yet-detected surviving companions from observations of nearby SN Ia remnants might favor low-mass companions, short binary separation, or stronger supernova explosion energies than the standard single-degenerate channels.

Observational evidence suggests that SNe Ia must come from thermonuclear explosions of carbon-oxygen white dwarfs in binary systems (Hillebrandt et al. 2013). The binary companion could be either another degenerate white dwarf, the double-degenerate scenario (DDS; Iben & Tutukov 1984;Webbink 1984), or a nondegenerate companion, such as mainsequence (MS) stars, helium stars, or red giants (Ruiz-Lapuente 2019), the single-degenerate scenario (SDS; Whelan & Iben 1973;Nomoto 1982). One major theoretical prediction to distinguish between the DDS and the SDS is that a surviving companion is expected in the SDS through several multidimensional hydrodynamics calculations on the interaction of the supernova ejecta with a nondegenerate companion (Pan et al. 2010(Pan et al. , 2012a(Pan et al. , 2012b(Pan et al. , 2013(Pan et al. , 2014Liu et al. 2013;Bauer et al. 2019;Zeng et al. 2020;. Direct observation of a surviving companion could link the progenitor channel to the explosion environment, which would place meaningful constraints on the intrinsic properties of SNe Ia. Numerous observational attempts to search for surviving companions in Type Ia supernova remnants (SNRs) have been made in the past decade, including searches in galactic Ia SNRs, e.g., Tycho's SNR (Ruiz-Lapuente et al. 2004;González Hernández et al. 2009;Kerzendorf et al. 2013), Kepler's SNR (Ruiz-Lapuente et al. 2018), and SN1006 (González Hernández et al. 2012;Kerzendorf et al. 2018), as well as extragalactic Ia SNRs in the Large Magellanic Cloud (Schaefer & Pagnotta 2012;Pagnotta & Schaefer 2015;Li et al. 2017Li et al. , 2019Litke et al. 2017). A few abnormal stars have been found, but no conclusive evidence for the existence of surviving companions (Ruiz-Lapuente 2019).
Theorists have been puzzling about the nondetection of surviving companions for more than a decade. The prediction of surviving companions is based on the assumption that the mass transfer from the nondegenerate companion is through Roche lobe overflow (RLOF) or winds in close binary with a binary separation, a ∼ 3 R * , in the majority of SDS channels, where R * is the companion radius (Pan et al. 2012a;Liu et al. 2013;Bauer et al. 2019).
In order to avoid high mass-transfer rates during the binary evolution, which might end up in a common envelope instead of a supernova, the optically thick wind model (Hachisu et al. 1996) and the common envelope wind model (Meng & Podsiadlowski 2017) are successful models to explain the observed Type Ia supernova (SN Ia) rates and allow more massive companions in the initial configuration. For instance, in the MS channel in Hachisu et al. (2008), the allowed companion mass at the onset of the SN Ia explosion is about 1.5-2 M e . However, in Pan et al. (2012b), the predicted surviving companions should still be detectable in nearby Ia SNRs (L ∼ 10-300 L e ).
Di Stefano et al. (2011) proposed that if the mass transfer from a binary companion carries angular momentum, the accreting white dwarf could be spun up and increase the Chandrasekhar mass for explosion. If the rotating and accreting white dwarf could reach the new Chandrasekhar mass, a super-Chandrasekhar-mass explosion could happen; otherwise, it takes extra time to spin down. The explosion will be delayed after an unknown spin-down period. During this spin-down phase, the binary configuration could be significantly changed depending on the spin-down timescale. For instance, extra mass loss of the companion star and redistribution of orbital angular momentum could happen (Di Stefano et al. 2011;Di Stefano & Kilic 2012). Therefore the prediction of surviving companions could be very different, as well (Wang et al. 2010;Di Stefano & Kilic 2012;Meng & Luo 2021). It should be noted that this spin-up/spin-down model still has large uncertainty. Thus, we use a parameterized method to mimic this effect.
In this paper, we relax the assumption of RLOF and conduct a parameter study on the evolution of surviving companions in the MS channel. The paper is organized as follows. In Section 2, we describe our simulation codes and numerical setup. In Section 3, we present the results of our simulations of surviving companions. Finally, we summarize our results and conclude in Section 4.

Numerical Methods
Our simulation setup is essentially similar to what has been described in Pan et al. (2010Pan et al. ( , 2012aPan et al. ( , 2012b, but with several improvements and updates. We divide the simulations into three phases. The first phase uses the stellar evolution code MESA to construct the nondegenerate binary companions at the onset of the explosion. The second phase involves twodimensional (2D) hydrodynamics simulations of the supernova ejecta and the companion interactions using the FLASH code. The third phase reads the outcome from the hydrodynamics simulations as input parameters to perform long-term surviving companion evolutions using MESA. In the following subsections, we describe the detailed set-up of each phase.

Main-sequence-like Companions
We use MESA r10398 (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019 to construct four zero-age MS models with masses M i = 0.8, 1, 1.5, and 2 M e and with metalicity z = 0.02. Table 1 summarizes the mass and radius information of our four considered companions. Note that realistic companions should be evolved during the mass-transfer stage, as considered in Pan et al. (2012b), However, the detailed binary evolution will make the companion's structure more complex, and it is unclear how the spin-down process will make a difference in the companion's structure. Therefore, the dependence of companion mass on the surviving companion evolution will be less meaningful than zero-age MS stars. For simplicity, we consider zero-age MS companions in this study for a parameter study.

Supernova Impact Simulations
Once the companion models are created, we map the density, temperature, and compositions into axisymmetric, 2D grids in a hydrodynamics code, FLASH (version 4; Fryxell et al. 2000;Dubey et al. 2008). We consider only hydrogen ( 1 H), helium ( 4 He), and carbon ( 12 C) from the companion model in FLASH. All elements heavier than carbon are approximated as carbon. We use the Helmholtz equation of state (Timmes & Arnett 1999;Timmes & Swesty 2000) and the new multiple Poisson solver (Couch et al. 2013) for self-gravity with the maximum angular moment for spherical harmonics l 80 max = .
Since the companion star will receive a kick from the supernova ejecta, we consider a simulation box that includes r < 15 R * and −20 R * < z < 10 R * in axisymmetric 2D cylindrical coordinates with 10 levels of refinement, yielding a finest cell spacing of 1/546 R * , where R * is the companion radius (see Table 1). The simulations with initial companion mass M i = 0.8 M e have a twice bigger simulation box in each dimension due to a higher kick velocity with this companion, but we keep the same spatial resolution by increasing one additional level of refinement. The "outflow" boundary condition is applied in both outer boundaries. The companion star is initialized at grid origin and relaxed in the 2D grids for ∼10 dynamical timescales, During the relaxation phase, we artificially damp the velocity by a factor of 0.97 in each time step. After relaxation, a W7-like explosion (Pan et al. 2012a) is imposed on the positive z-axis at z = a × R * , where a is the binary separation. We use nickel ( 56 Ni) as a tracer to represent the supernova ejecta. In the W7 model (Nomoto et al. 1984), the white dwarf mass M W7 is 1.378 M e , the total explosion energy E W7 is 1.233 × 10 51 erg, and the averaged ejecta speed V W7 is 8.527 × 10 8 cm s −1 . To further consider super-or sub-Chandrasekhar-mass explosion models, we vary the explosion energy from 0.1 to 10 times of the original value in the W7 model. The corresponding averaged ejecta speeds are scaled as , where E SN is the explosion energy. We ignore the orbital motion and companion's rotation in this study for simplicity.

Long-term Evolutions of Surviving Companions
Since the time step in FLASH simulations is limited by the local sound speed, we could not afford long-term surviving companion evolution within FLASH, and we do not consider nuclear burning in our simulations. Therefore, we terminate the FLASH simulations after ∼70 dynamical timescales. Following the method described in Bauer et al. (2019), we measure the stripped mass and entropy changes (Δs) in our hydrodynamics simulations and assume that the specific entropy will be Note. M i is the mass of the companion star right before the supernova explosion, R * is the companion radius, and conserved during the hydrostatic equilibrium process. We have tested the entropy changes at different times during the FLASH simulation. No significant changes will appear in the subsequent MESA simulation.
We relax the mass of the preexplosion companion models to the final bound mass in the postimpact FLASH simulations within MESA. Once the stellar mass is relaxed, we then apply a user-defined heating rate,  T s t  = D D , in MESA (Bauer et al. 2019), where Δt is the heating timescale and Δt = 3 days. After the heating process, we let MESA continue running for up to 1 million years to describe the evolution of postimpact surviving companions. A convergence test with different heating timescales has been performed. We find no significant difference in the long-term evolution after 10Δt.

Results
All the simulation parameters we considered in this study are summarized in Table 2. For each companion model in Table 1, we vary the binary separation from the distance that just allowed Roche lobe overflow (denoted R RLOF ) to 5 R * or 6 R * to investigate the separation effects on the evolution of surviving companions. In addition, we simulate nine extra models with different explosion energies to explore the energy effects. Figure 1 shows the density evolution of the model 2M3R (see Table 2) after a SN Ia explosion at different times. The overall features are consistent with what has been reported in Pan et al. (2012a). Important features that could provide evidence in searches for surviving companions can be summarized as follows.

Overview of the Supernova Impact Simulations
(1) The collision of supernova ejecta with the binary companion at around one shock crossing time (t a v cross SN = ) might emit strong X-ray or UV signals (Kasen 2010;Burke et al. 2021), corresponding to the left panel in Figure 1. (2) The ejecta ram pressure (heat) will strip (ablate) a certain amount of mass of the companion star (Wheeler et al. 1975;Marietta et al. 2000;Pan et al. 2012a). (3) During the collision, the linear momentum of the supernova ejecta could transfer to the surviving companion, resulting in a kick velocity normal to the original orbital velocity (Wheeler et al. 1975;Marietta et al. 2000;Pan et al. 2012a).   Note. M f is the final mass of the companion star after the supernova, a is the binary separation between the companion star and the white dwarf, E SN is the supernova explosion energy in units of W7 energies (E W7 = 1.233 × 10 51 erg), v orb is the orbital speed of the companion right before the supernova explosion, v kick is the kick is the final linear velocity of the surviving companion, D heat is the heating depth which describes the location of maximum heating in the mass coordinate (see Section 3.1), E heat is the amount of energy deposition in the surviving companion, and M Ni is the final bound mass of 56 Ni onto the companion.
impact than massive companions (M i = 1.5 and 2.0 M e ) with the same a/R * ratio. In addition, massive companions have higher gravitational binding energies ( ). Therefore, the fractions of final unbound mass with lighter companions are significantly higher than massive companions.
The model 0.8M3R has the closest binary separation to the exploding white dwarf and the lowest gravitational binding energy among all considered companion models. We find that this model is almost completely destroyed by the supernova explosion. The final bound mass is less than 0.3 M e at the end of the simulation, and the kick velocity has reached >160 km s −1 . Although the model 0.8M3R has a twice larger simulation box, the companion reaches the bounding box at t = 10 hr due to its high kick velocity. Therefore, we label the final bound mass (kick velocity) as an upper (lower) limit in Table 2.
The relations between the final unbound mass and binary separation can be fitted by a power-law relation, and so can the kick velocity (see Figure 4). The fitted power-law relations are summarized as below: ⎞ ⎠ Note that, in principle, there are two physical mechanisms to unbound the companion mass, including stripping and ablation (Wheeler et al. 1975;Pan et al. 2012a). In the previous literature, some authors use the term "stripping mass" to describe the total unbound mass that includes both stripped mass and ablated mass. In this paper, we use the term "final unbound mass" to avoid confusion.
The shaded gray band in Figure 4 represents one standard deviation of the fitted parameters. One should note that, in the previous literature, the power-law relations are evaluated based on the binary separation (a) only. The scaling factor, * M R i , could   Table 2). The scattering between evolved MS models and zeroage main sequence models in this study could be explained by the changes of stellar structure during the binary evolutions, but the overall fitting is quite successful, suggesting that this fitted formula could be universal, in first order of magnitude, for all MS channels. Note that for three-dimensional (3D) hydrodynamics simulations that include the initial orbital motion, only the final linear velocity (including orbital and kick velocity) is recorded. Thus, the kick velocities of these simulations are not trivial to extract, and, therefore, we do not include the data points from Pan et al. (2012b) for the kick velocity panel in Figure 4.
The postimpact surviving companion receives heating from the supernova ejecta, as well. Pan et al. (2012b) suggest that the postimpact evolution of surviving companions depends on both the amount of the supernova heating and the depth of energy deposition. The amount of supernova heating is measured by the total heating from the entropy evolution history in MESA and is summarized in Table 2. In general, the models with close binary separation receive a stronger supernova energy due to a larger solid angle, but they also suffer a more violent massstripping effect. The net effect makes no clear relations between the supernova heating and binary separation.
However, the location, or the depth, of supernova heating is sensitive to the binary separation (see Figure 5). Here, we define the heating depth ) m T local max is the mass with local temperature maximum due to the supernova heating, and M f is the final companion mass after supernova impact.
A linear relationship between the heating depth and the binary separation can be fitted as below: Note that the heating depth is insensitive to the companion structure of our considered MS companions. Furthermore, it is possible that part of the supernova ejecta could be bound on the surface of the surviving companion and therefore contaminate its spectrum (Pan et al. 2012a). This supernova could be another example of a "smoking gun" evidence for searches of surviving companions in SNRs. In Table 2, we summarize the final bound nickel mass on the surviving companions. We notice that some cases with strong supernova impact have bound nickel masses close to rounding errors (∼10 −14 M e ). Thus, we only list an upper limit of the bound nickel mass in these cases. It is found that only lowexplosion-energy cases (models 1M3RE0.1 and 1.5M3RE0.1) or large-binary-separation cases (0.8M6R) have a bound nickel mass 10 −5 M e , which is lower than the previous estimation of the nickel contamination from 3D simulations in Pan et al. (2012a). The lower-bound nickel mass in our high-resolution 2D simulations might be due to the stronger turbulent convections around the companion envelope during the supernova impact, compared to 3D simulations.

Explosion Energy Dependence
The amount of supernova explosion energy is another important factor that should affect the conditions of a surviving companion but has not been fully discussed in the past (except in Pakmor et al. 2008). The energy difference between subluminous SNe Ia and overluminous SNe Ia could be more  Table 2 with different binary separations. The black line shows a linear fit in log-log scale using numpy.polyfit. The gray shaded region indicates a 2σ confidence interval. Right: similar to the left panel but for the kick velocity. corresponds to the energy density that hit on the companion star. The black line is the linear fit using numpy.polyfit. We label the 2σ interval using a gray shaded region. than an order of magnitude (Taubenberger 2017;Jha et al. 2019). Figure 6 shows the influence of supernova explosion energies on the final unbound mass and kick velocities with models 1M and 1.5M. Similar to what has been reported in Pakmor et al. (2008), the final unbound mass is linearly proportional to the explosion energies regardless of companion mass, and can be described by Equation (4): is the percentage of unbound mass normalized to the corresponding model with a = 3 R * and E E

SN W7
= , and ΔM = M i − M f is the mass difference. In Pakmor et al. (2008), only one companion model is considered. In Equation (4), we additionally show that the percentage of the unbound mass is inversely proportional to the companion radius. This could be understood by the fact that the energy received on the companion is * R a 2 µ and we have chosen a = 3 R * .
The right panel in Figure 6 shows the fraction of kick velocity as functions of explosion energy, where ( ) * v R kick 3 is the kick velocity normalized to the corresponding model with a = 3 R * and E E SN W7 = . The relation could be fitted by a power law in Equation (5): The power-law index 0.448 is less than 0.5, suggesting that a small fraction of ejecta energy is converted to thermal energy by supernova heating.

Evolution of Surviving Companions
As described in Section 2.3, we measure the mass losses and entropy changes after the supernova impact in FLASH and use them as parametrized heating and mass loss in MESA. Figure 7 shows the profiles of the entropy changes for our considered models. For higher mass loss models, i.e., models 0.8M3R, 1MRL, and 1M3R, the supernova energy could penetrate the companion and significantly heat the companion up to a depth of 0.4m/M f , where m/M f is the fraction of the enclosed mass in the mass coordinate. For more compact companions, such as models 1.5M and 2M, the supernova impact contributes heating mainly around the companion surface (∼0.8m/m * ). The depth of supernova heating determines the local thermal timescale of a surviving companion, and therefore affects the timescale of the postimpact surviving companion's evolution (Pan et al. 2012b). Figure 8 shows the time evolution of the bolometric luminosity, effective temperature, and photosphere radius of our considered models. The overall features are similar to what have been reported in Pan et al. (2012b) with evolved MS companions. We find that the luminosity ranges and effective temperatures in our cases with binary separation that just allowed RLOF match with the ranges in Pan et al. (2012b), suggesting that the different numerical treatment of supernova impact in MESA in this study is consistent with the method described in Pan et al. (2012b). The sharp increment at around 10 −3 yr is due to the sudden heating from the supernova ejecta in MESA with a heating timescale Δt = 3 days. A second brightening happens at ∼10 1 -10 2 yr for models 1M, 1.5M, and 2M, when the deposited energy diffuses away by expanding the photosphere. This timescale is described by the local thermal timescale in the stellar envelope (Pan et al. 2012b). The deeper the energy deposition, the longer the local thermal timescale. In addition, models with larger binary separations tend to have shallower energy deposition (see Figure 7). Therefore, a more massive companion (less mass stripping) with large binary separation is expected to be overluminous in young SNRs. Once the deposited energy is released, the surviving companion starts to contract due to gravity.   Table 2). Since the final unbound mass of model 0.8M3R is still increasing at the end of the simulation (see Figure 2), we plot it with a dashed line. Figure 9 shows the corresponding evolutionary tracks of our considered models in Hertzsprung-Russell diagrams. We find that the evolutionary tracks can be categorized in two types based on their compactness (similar to the behavior in the final unbound mass). Soft companions (models 0.8M and 1.0M) have higher final unbound mass and deeper energy deposition, resulting in less luminous surviving companions with L < 10 L e . On the other hand, stiff companions (models 1.5M and 2M) have less final unbound mass and shallower energy deposition, leading to a luminous surviving companion with L > 10 3 L e for hundreds of years after the explosion time (t 0 ). The bottom panels in Figure 9 present the evolutionary tracks of surviving companions with different explosion energies. Regardless of compactness, stronger explosions lead to deeper energy deposition and a larger final unbound mass, and therefore produce a dimmer surviving companion. For instance, model 1.5M3R has a peak luminosity  L L 500 max~a t around 35 yr after the supernova explosion. However, if the supernova explosion energy is 5 E W7 , the peak luminosity could drop to  L L 23 max < in Figure 9. The current nondetection of surviving companions in nearby SNRs might be explained by soft surviving companions, shorter binary separations due to delayed explosions, or super-Chandrasekhar-mass explosions with a higher explosion energy. In the latter case, the supernova impact will be stronger than a typical binary separation by RLOF, which might show stronger UV excess in the early lightcurve (Kasen 2010).  Table 2. Different colors represent models with different initial binary separation. Different rows indicate different companion models.

Summary and Conclusions
We have performed 2D hydrodynamics simulations of supernova impacts and the subsequent one-dimensional stellar evolution of surviving companions with a wide range of parameters of mainsequence-like companions. It should be noted that the simulations presented are 2D, which ignores both spin and orbital motions, and with nonevolved binary companions. Pan et al. (2012a) suggest that the spin and orbital motions could make some difference during the supernova impact and therefore might affect the later-on surviving companion evolutions (but less important, as suggested by . Despite this, we find three universal relations (Equations (1), (2), and (3)) to describe the final unbound mass, kick velocity, and supernova heating depth as functions of binary separation, regardless of the companion mass or evolution stage in the MS. Stiff surviving Figure 9. Evolutionary tracks on Hertzsprung-Russell diagrams for our models described in Table 2. Different colors represent different initial binary separations (top two rows) or different explosion energies (bottom rows). We label the times since explosion with different symbols. t 0 indicates the time of the supernova explosion. Dotted lines present the phase of supernova heating in MESA. The same as in Figure 7, we label the 0.8M3R, 1M3RE2, and 1.5M3RE10 models with dashed lines. companions with larger progenitor masses (models 1.5M and 2M) have less mass stripping during the supernova impact and receive shallower supernova heating in the envelope. The deposited energy in the envelope will be diffused away within hundreds of years, allowing the photosphere to expand. This type of surviving companion will be overluminous (  -L L 10 10 max 2 4 > ) and have a higher chance to be detected in nearby SNRs. In contrast, soft surviving companions with lower progenitor masses (models 0.8M and 1M) result in larger final unbound masses and deeper supernova heating. These low-mass surviving companions could remain  L L 10 max < , which could explain the nondetection of surviving companions in nearby SNRs. In addition, a super-Chandrasekharmass explosion with a higher explosion energy or a shorter binary separation during the spin-up/spin-down phase could produce lowluminosity surviving companions. In these cases, the surviving companions will have a higher kick velocity, and the early UV or X-ray light from the collision between the supernova ejecta and companion will be stronger than the previous estimate in Kasen (2010) for MS-like companions.
We thank the anonymous referee for his/her valuable comments and suggestions that helped to improve this paper. S.J.R. and K.C.P. acknowledge valuable discussions with You-Hua Chu and Hsiang-Yi Karen Yang on the searches for surviving companions. S.J.R. also thanks Dr. Stephen Justham for his useful comments and suggestions during her virtual visit in IMPRS. This work is supported by the Ministry of Science and Technology of Taiwan through grant Nos. MOST 107-2112-M-007-032-MY3 and MOST 110-2112-M-007-019, by the Center for Informatics and Computation in Astronomy (CICA) at National Tsing Hua University through a grant from the Ministry of Education of Taiwan. FLASH was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. The simulations and data analysis have been carried out at the Taiwania supercomputer in the National Center for High-Performance Computing (NCHC) in Taiwan, and on the CICA cluster at National Tsing Hua University. Analysis and visualization of simulation data were completed using the analysis toolkit yt.