Delayed Detonation Thermonuclear Supernovae with an Extended Dark Matter Component

DM is believed to be the dominant component of galaxy clusters (Clowe et al. 2006), large-scale structures of the universe (Davis et al. 1985), as well as a host of other astrophysical objects (Oppenheimer et al. 2001). Yet the searches for the DM particles have not yielded any convincing candidates. The Large Hadron Collider searches for DM signals are still ongoing. The possibly positive signals from the Xenon1T DM detector (Aprile et al. 2020) may provide hints to the existence of DM. Nevertheless, astrophysical searches for the DM will remain an important channel to complement terrestrial experiments.

Given that DM is the major component in the universe, stellar objects comprising normal matter (NM) and DM may form. For instance, a significant amount of DM could be captured by the gravitational potentials of massive stars located near the galactic center, where the DM density is high (Sulistiyowati et al. 2014). Studies of self-interacting and non-self-interacting DM-admixed neutron stars suggested that compact stars can be a probe to DM (de Lavallaz & Fairbairn 2010). Studies on the equilibrium structures of non-self-interacting DM-admixed neutron stars provided hints to search for DM from the neutron star diversities (Leung et al. 2011; Ellis et al. 2018). Ongoing studies have been focusing on different aspects of compact objects, including bosonic DM stars (Eby et al. 2016; Lee et al. 2021, in preparation), formation of neutron stars through the accretion-induced collapse of DM-admixed white dwarfs (WDs; Leung et al. 2019; Zha et al. 2019), tidal properties of DM-admixed neutron stars (Leung et al. 2021, in preparation), thermalization of the WD core to produce a thermonuclear runaway (Bramante 2015; Acevedo & Bramante 2019; Janish et al. 2019; Steigerwald et al. 2019), and thermonuclear supernovae with point-mass DM cores (Leung et al. 2015a). Other than compact objects, effects of DM admixture on stellar formation (Arun et al. 2019) and evolution (Scott et al. 2008; Freese et al. 2016; Lopes & Lopes 2019) are actively investigated. In particular, the evolution of AGB stars is used to constrain the properties of axion-like particles (Dolan et al. 2021). DM-admixed astrophysical objects have become a new channel to search for astrophysical DM.

Acevedo & Bramante (2019), Bramante (2015), and Steigerwald et al. (2019) pointed out that accreting DM to a WD could lead to a kinetic energy transfer from the DM to the NM through particle scattering, which helps to heat up the carbon–oxygen core toward the temperature necessary for triggering a thermonuclear runaway. On the other hand, Leung et al. (2013b) showed that the admixture of non-self-interacting DM could reduce the Chandrashekhar mass of a WD so that a thermonuclear flame is generated when the mass of the WD is below the traditional Chandrashekhar limit.

Leung et al. (2015a) have demonstrated that DM-admixed thermonuclear supernovae could produce subluminous light curves. If the heavy fermionic DM consists of particles with a mass of around 1 GeV, the DM component is small enough to be approximated as a stationary point gravitational source located at the center of the WD. Recent constraints on the DM particle mass derived from lab experiments (Cirelli et al. 2021) and observed neutron star masses (Ivanytskyi et al. 2020) have not ruled not the parameter space for sub-GeV DM particles. Some extensions of the Standard Model of particles physics have also suggested the existence of sub-GeV DM particles (Dutta et al. 2019). Therefore, we are motivated to study the effects of admixing sub-GeV DM to a thermonuclear supernova. As we shall see, admixing sub-GeV DM would create an extended DM component that we could not treat as a point mass, and so we need to describe the explosive hydrodynamics using a two-fluid formalism.

WDs are believed to be the progenitors for thermonuclear supernovae (Hillebrandt et al. 2013; Maoz et al. 2014; Livio & Mazzali 2018). However, it is not clear which kind of WDs and explosion mechanisms are responsible for such energetic events (Ruiter 2019). It is believed that most of the progenitor WDs consist of carbon and oxygen (Livio & Mazzali 2018). A mixture of carbon, oxygen, and neon seems to be also possible (Denissenkov et al. 2014, 2015; Marquardt et al. 2015; Bravo et al. 2016; Willcox et al. 2016; Brooks et al. 2017; Augustine et al. 2019; Ruiter 2019). The progenitor models for thermonuclear supernovae are classified as either the single-degenerate scenario or the double-degenerate scenario (Maoz et al. 2014; Livio & Mazzali 2018; Ruiter 2019).

The single-degenerate scenario has been studied for several decades. Early simulations investigated thermonuclear supernovae that are driven by a pure detonation (Arnett 1969) and a pure deflagration (Nomoto et al. 1976) initiated near the center of the WD. While there were overproductions of iron-group elements in the former, the latter was too slow to generate a healthy explosion. Results from either pure detonation or pure laminar deflagration do not fully reproduce the observed features in a typical thermonuclear supernova. Moreover, it has been argued that a pure detonation would be quenched at high density, casting doubts on the possibilities of a prompt detonation (Kriminski et al. 1998).

A deflagration wave is subject to several hydrodynamical instabilities (Hillebrandt et al. 2013; Röpke & Sim 2018). Turbulence would induce convolutions on the deflagration surface, which would allow the deflagration to propagate at an effective speed greater than the laminar value (Nomoto et al. 1984; Röpke 2017). Early studies of the pure turbulent deflagration model used the mixing-length theory described in Nomoto et al. (1976, 1984) (the W7 model). A fast deflagration wave with a propagating speed around 20%–30% of the local speed of sound was assumed. These models could produce healthy explosions. However. they produced too much ⁵⁸Ni when compared with the constraints from the galactic chemical evolution (Iwamoto et al. 1999). ³ Multidimensional studies on the pure turbulent deflagration model also underproduced the amount of ⁵⁶Ni and explosion energy (Reinecke et al. 1998a; Calder et al. 2004; Fink et al. 2014; Röpke & Sim 2018), with most of the unburned material left behind at the center of the WD. This is inconsistent with spectroscopic observations (Ma et al. 2013). The carbon–oxygen fuel should only exist at the outer ejecta (Gamezo et al. 2003). Although later studies showed that some good agreements with observations in terms of the nickel mass and the ejecta velocity could be obtained (Reinecke et al. 2002a; Röpke & Hillebrandt 2005), another major shortcoming of the pure turbulent deflagration model is the underproduction of the silicon-group elements at the outer layers of the ejecta (Röpke et al. 2007), which are strongly mixed with the iron-group elements (Ma et al. 2013).

A transition to detonation from a deflagration is observed in several terrestrial experiments with closed boundaries (Poludnenko et al. 2019; Brooker et al. 2020a). It is therefore natural to hypothesize that a transition to detonation could also happen in open boundaries such as the interior of a WD. Khokhlov (1991) proposed the DDT model, which gives a comparable amount of iron-group and silicon-group elements to the observed thermonuclear supernovae. The DDT model has also been shown to fit the observed supernova light curves in the V and R bands reasonably well (Hoeflich et al. 1995). To trigger the DDT, the Zel'dovich gradient mechanism (Zel'dovich et al. 1972) is used. Besides, there are also other proposed mechanisms (Charignon & Chièze 2013). In the Zel'dovich mechanism, the local eddy motions near the flame front would help flatten the temperature gradient to trigger the DDT. (Niemeyer 1999; Branch 2017; Poludnenko et al. 2019; Brooker et al. 2020a.) Whether the turbulence could robustly help to attain the required criteria for a DDT is surely a concern (Niemeyer 1999; Lisewski et al. 2000; Woosley 2007). In principle, turbulent motions could mix the cold fuel and hot ash to obtain the required preconditioned field (Pan et al. 2008; Hillebrandt et al. 2013). However, numerical simulations have shown that the required turbulent strength (Pan et al. 2008) would be too large, so that the implied velocity fluctuations (Roepke 2007) is only marginally sufficient for triggering the DDT. Yet some semi-analytical studies (Shu et al. 2020) and hydrodynamical simulations using a sub-grid scale turbulent model have confirmed that a WD can meet the criteria for a DDT (Ciaraldi-Schoolmann et al. 2013). The required strength of the turbulence was later found to be less than that previously believed (Poludnenko et al. 2019), and a DDT is possible in an environment with weak turbulence (Brooker et al. 2020b). Furthermore, the DDT model has been shown to produce spectra that match well with some of the observed supernovae (Blondin et al. 2011, 2013; Sim et al. 2013).

In Leung et al. (2015a), the DM-admixed thermonuclear supernova is studied by considering the pure turbulent deflagration model without a DDT. In this work, we extend the study of Leung et al. (2015a) by (1) focusing on the robust DDT model and (2) adding the dynamics of the DM component in the picture by using the fully nonlinear two-fluid simulations. The plan of the paper is as follows: Section 2 describes the construction of the supernova progenitor, the tools for simulating the explosions, and the post-processing methods. Section 3 summarizes the simulation results. Section 4 discusses the parameter dependences of our models and the implications of our results, and Section 5 concludes our study.

We construct a series of DMWDs as the supernova progenitors, by solving the Newtonian hydrostatic equilibrium equations (see, for example, Sandin & Ciarcelluti 2009; Ciarcelluti & Sandin 2011; Leung et al. 2019): ⁴

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{i}}{{dr}}=-\displaystyle \frac{G({m}_{1}+{m}_{2}){\rho }_{i}}{{r}^{2}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{i}}{{dr}}=4\pi {r}^{2}{\rho }_{i}.\end{eqnarray} \tag{ 2 }$$

Here, the subscript i = 1(2) denotes the DM (NM) quantities, and ρ, p, m, and r are the density, pressure, enclosed mass, and radial coordinate, respectively. We adopt the ideal degenerate Fermi gas equation of state (EOS), assuming spin $\tfrac{1}{2}$ for the DM component (Teukolsky & Shapiro 2008). We vary the DM particle mass from 0.1–0.3 GeV in 0.05 GeV steps. For each value of the DM particle mass, we vary the DM central density from 10⁸–10¹⁰ g cm⁻³. We show how we obtain the supernova progenitors by analyzing this series of models.

Since the DM only affects the dynamics of the NM by gravity, in this study we identify the primary effects of admixing DM to thermonuclear supernovae by considering the spherically symmetric Newtonian two-fluid hydrodynamics in the Eulerian framework. We show in Appendix C that the dynamics of DM are indeed necessary for this scenario. Multidimensional two-fluid simulations in this scenario are computationally demanding because during the late phase of the explosion, a large simulation box is required when the NM is extended to a much larger radius, while a sufficiently fine grid resolution is needed to resolve the central DM component. In a typical one-dimensional simulation, the Lagrangian formalism is often used. We do not consider this formalism because a further interpolation of the local gravitational force by both fluid components becomes necessary when they do not overlap with each other. The two fluids may also have different masses and radii, which makes the definition of the mass coordinate difficult.

Our code is constructed based on the two-fluid hydrodynamics code developed by Ka-Wing Wong (2011). Unless otherwise noted, we use geometric units with G = c = M_⊙ = 1. Different from the work by Leung et al. (2015a), we have adopted the finite-volume approach where the Euler equations in the spherical coordinate are given as (Mignone 2014)

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{{\boldsymbol{U}}}_{i}=-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}({r}^{2}{{\boldsymbol{F}}}_{i})+{{\boldsymbol{S}}}_{i}+{{\boldsymbol{G}}}_{i},\end{eqnarray} \tag{ 3 }$$

with the state vector U and the flux vector F given as

$$\begin{eqnarray}&&{{\boldsymbol{U}}}_{i}=\left(\begin{array}{c}{\rho }_{i}\\ {\rho }_{i}{v}_{i}\\ {E}_{i}\\ {{\rm{\Psi }}}_{i}\end{array}\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{i}=\left(\begin{array}{c}{\rho }_{i}{v}_{i}\\ {\rho }_{i}{v}_{i}^{2}+{p}_{i}\\ ({E}_{i}+{p}_{i}){v}_{i}\\ {\rho }_{i}{{\rm{\Psi }}}_{i}\end{array}\right).\end{eqnarray} \tag{ 5 }$$

Here, _i is the internal energy, v_i is the radial velocity, and ${E}_{i}={\rho }_{i}({\epsilon }_{i}+\tfrac{1}{2}{v}_{i}^{2})$ is the total energy density. Ψ_i is any passive scalar quantity. For example, Ψ_i could be the isotope mass fraction X(¹² C) of the element ¹²C. The geometric source term S and the extra source term are given as

$$\begin{eqnarray}&&{{\boldsymbol{S}}}_{i}=\left(\begin{array}{c}0\\ 2{p}_{i}/r\\ 0\\ 0\end{array}\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{i}=\left(\begin{array}{c}0\\ -{\rho }_{i}\displaystyle \frac{d\phi }{{dr}}\\ -{\rho }_{i}{v}_{i}\displaystyle \frac{d\phi }{{dr}}\\ 0\end{array}\right),\end{eqnarray} \tag{ 7 }$$

where ϕ is the gravitational potential, governed by the two-fluid Poisson equation:

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\phi =4\pi ({\rho }_{1}+{\rho }_{2}).\end{eqnarray} \tag{ 8 }$$

We adopt the monotonicity-preserving 5th-order scheme (Suresh & Huynh 1997) to reconstruct the primitive variables at the cell interfaces. We use the Lax–Friedrich Riemanns solver (Jiang & Shu 1996) to compute the interface numerical fluxes, and we discretize the temporal evolution using the method of lines where the strong stability-preserving 3rd-order Runge–Kutta method is implemented (Gottlieb et al. 2011).

We adopt a modified monopole solver (Couch et al. 2013) for the gravitational potential because the error due to computing the potential at the cell center while using the mass density at that point could be reduced. We implement a moving-grid algorithm (Roepke 2005) to follow the expansion of the ejecta until it becomes homologous. The moving-grid equations contain a grid velocity v_f so that the fluxes and source term vectors are modified as follows (Leung & Nomoto 2018):

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{i}\to {{\boldsymbol{F}}}_{i}-{{\boldsymbol{U}}}_{i}{v}_{f},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{2}\to {{\boldsymbol{G}}}_{2}-{{\boldsymbol{U}}}_{2}\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}({r}^{2}{{\boldsymbol{v}}}_{f}).\end{eqnarray} \tag{ 10 }$$

To capture all exothermic reactions before the expanding material reaches the boundaries so that the moving-grid algorithm is triggered, we set the initial grid box size as 2.7 × 10⁹ cm, which is more than 10 times the typical size of the progenitor.

To reduce the computational time, we implement a reduced nuclear reaction network, which was used by Leung & Nomoto (2018, 2020). It was developed based on the original works by Townsley et al. (2007) and Calder et al. (2007). The nuclear network includes seven isotopes: ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, and ⁵⁶Ni. The nuclear network separates nuclear burning into three stages:

• 1.  
  Combustion of ¹²C to ²⁴Mg;
• 2.  
  Burning of ¹⁶O and ²⁴Mg to ²⁸Si;
• 3.  
  Nuclear statistical equilibrium (NSE).

When the deflagration first sweeps through the fuel, step 1 is triggered. The burning proceeds to step 2, which is called the nuclear quasi-statistical equilibrium (NQSE) when the hydrodynamical time step δ t is larger than the quasi-equilibrium time step ${\tau }_{\mathrm{NQSE}}=\exp (182.06/{T}_{9}-46.064)$. Here, T₉ is the temperature in 10⁹ K. Finally the burning proceeds to step 3 if the hydrodynamical time step is larger than the NSE timescale ${\tau }_{\mathrm{NSE}}=\exp (196.02/{T}_{9}-41.645)$. These timescales are provided in, for example, the work by Townsley et al. (2007) and Calder et al. (2007). For any incomplete burning (δ t < τ_NQSE or τ_NSE), we will use linear interpolation to compute the resulting compositions and energy released. We will discuss the accuracy of the linear interpolation in Appendix A.

To capture the propagation of a subsonic deflagration, which has a width of ∼10⁻⁴ cm (Branch 2017), we adopt the level-set method (Sethian 1999) as our flame capturing scheme. This method has been used by several authors in their studies of multidimensional thermonuclear supernovae (Reinecke et al. 1998b, 2002a, 2002b). To initiate the deflagration, we plant a level-set enclosing ∼0.02 M_⊙ of NM for all models. The spherical symmetry does not allow us to compute the turbulence production explicitly. Therefore, we follow the method in the literature (Nomoto et al. 1984; Hoeflich et al. 1995; Woosley 1997; Blondin et al. 2012, 2015) to parameterize the deflagration speed in terms of the local speed of sound c_s . The typical value would be a fraction, possibly a few percent. After various numerical experiments on the pure NM explosions, we choose the fraction to be 0.06, and the DDT is manually triggered once the deflagration front reaches a density of 1.7 × 10⁷ g cm⁻³ (Iwamoto et al. 1999), so that the production of ⁵⁶Ni lies within the acceptable range of 0.4–0.7 M_⊙. During the onset of the DDT, the speed of the level-set is immediately raised to that of the detonation. The detonation for a density lower than ∼2 × 10⁷ g cm⁻³ is a Chapman–Jouguet (CJ) detonation with the speed of sound, while that for a density higher than ∼2 × 10⁷ g cm⁻³ is a pathological one. ⁵ The latter one is solved by the method described in Sharpe (1999) and Leung & Nomoto (2020).

Recent studies on thermonuclear supernovae show that neutrino production could be a key to distinguish different explosion models (Wright et al. 2016, 2017). It is therefore likely that the neutrino signals will also play an important role in DM-admixed models. We use an open-source neutrino emission subroutine ⁶ , which calculates the pair, photo-, bremsstrahlung, and recombination neutrinos with formulae derived in Itoh et al. (1996). The plasmon neutrinos are computed using the method by Kantor & Gusakov (2007). To compute the differential production rate for high-energy neutrinos (1–5 MeV), we adopt the method by Misiaszek (2006) to calculate the pair-neutrino spectrum. In addition, the plasma-neutrino spectrum is calculated using the method by Odrzywołek (2007). These methods have also been applied in some previous works (Leung et al. 2015b, 2020). The differential production rate is then integrated to obtain the total number of thermoneutrinos produced. Roepke (2005) found that the supernova ejecta would approach a homologous expansion in ∼5–10 s. Furthermore, most of the exothermic nuclear reactions would come to an end after ∼1 s, and the distribution of the isotope mass fractions does not change too much between 5 and 10 s. Moreover, we performed a number of numerical experiments for the pure NM explosions, and we found that due to the significant increment in the grid size, the mass and energy conservations could become worse if we keep increasing the simulation time. Therefore, we terminated the simulation at 5 s after the runaway has started. The resulting density, temperature, velocity, and isotope profiles are then mapped from our one-dimensional Eulerian form into a one-dimensional Lagrangian form, which is then inputted to the SuperNova Explosion Code (SNEC code; Morozova et al. 2015) to calculate the light curves.

We use the Helmholtz EOS (Timmes & Swesty 2000) for the NM, which is an equal mixture of ¹²C and ¹⁶O initially. We assume an initial constant temperature of T = 10⁸ K. We adopt the central density for the NM to be 3 × 10⁹ g cm⁻³. This value gives rise to a thermonuclear explosion that is close to the average thermonuclear supernova. Similar values have also been adopted by a number of previous studies (Maeda et al. 2010; Seitenzahl et al. 2011, 2013; Nomoto & Leung 2017; Kobayashi et al. 2020). We vary the DM central density from 10⁸–10¹⁰ g cm⁻³ to obtain the masses and radii for the NM and the DM components. We find that a higher DM admixture decreases the NM mass, which is consistent with the results in the previous study on DMWD using the Baym–Pethick–Sutherland (BPS) EOS (Leung et al. 2013a). Figures 1 and 2 show that for a smaller DM particle mass, the DM radii are more extended, and there could exist configurations where the DM and the NM components are comparable in radii and masses. These results are in contrast with the results by Leung et al. (2013a), for which heavier DM would reside as a point object. Figure 2 shows that for a given NM mass, a vertical line that crosses the corresponding DM mass curve gives the minimum amount of DM to be admixed so that the WD could attain the conditions necessary for generating a thermonuclear runaway. For the DM particle mass of 0.1 GeV, there is a sharp transition of the DM mass such that it increases drastically near the DM central density of ∼8.75 × 10⁸ g cm⁻³, while the changes in NM mass are relatively small. In that case, there exist configurations with ∼1 M_⊙ of the NM mass and ∼0.5 M_⊙ of the DM mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We see that the effects of admixing 0.1 GeV DM on the stellar parameters, such as the masses and radii, are particularly interesting. We compute five configurations using this DM particle mass as the supernova progenitors: DM-1, DM0, DM1, DM2, and DM3. The model with no DM admixture is named model NM. They are chosen with an increasing ratio of the DM to NM masses from 0.05 to ∼0.5. The parameters for the progenitors are shown in Table 1, while the density profiles for the NM and DM of each model are shown in Figure 3. The NM radius of Model DM2 is nearly the same as its DM radius, while model DM3 has an extended DM component that covers the entire NM. Note that in all of the above models, the NM is more massive than the DM. We will show in Appendix D that admixing DM of the order of O(10⁻¹)M_⊙ is indeed possible. We will investigate the thermonuclear explosions of these configurations. In particular, we are interested in how the extended DM component will affect the explosion dynamics and the supernova observables.

Zoom In Zoom Out Reset image size

The simulation results are summarized in Table 2, where we list the energy, ⁵⁶Ni generation, and time of the DDT. In general, the presence of the DM delays the transition time and reduces energy production. The DM gravitational potential and the interactions between the DM and NM prohibit the expansion of the NM during the deflagration phase. Therefore, the NM needs a longer time to reach the density necessary for a DDT. On the other hand, the ⁵⁶Ni masses synthesized are similar in magnitude when compared with the pure NM explosion. The ⁵⁶Ni mass first decreases with more DM admixture to a minimum of the model DM2, then it rebounds when we proceed to the model DM3. We plot the total ⁵⁶Ni mass variations in Figure 4. In general, models having more DM admixture would give lower ⁵⁶Ni production after the DDT. The steep density gradient would reduce the mass content at the outer envelope of the NM. The amount of ⁵⁶Ni for the DM-admixed models is compensated by a higher production before the DDT. The prolonged deflagration phase gives more time for the NM to synthesize ⁵⁶Ni. We also include the total energy variations in Figure 4. The energy production after the DDT is smaller for models with more DM admixture, which is consistent with the variations of the ⁵⁶Ni after the DDT.

Zoom In Zoom Out Reset image size

We show the maximum and central temperatures in Figure 5. We find that the presence of the DM tends to keep the NM hotter. This can be explained by the suppressed expansion during the deflagration phase, and thus the NM cannot dissipate its internal energy by the pressure work done. Energy conservation suggests that the loss in the internal energy would be converted not only to kinetic energy but also to gravitational energy. Therefore, for models with more DM admixture, the expansion velocity of the NM is lower. This can be seen in Figures 6 and 7, where the temperature, density, and velocity profiles for the NM of different models at the moment of the DDT are shown. This also explains why these models have later occurrences of the DDT. The first peak in the maximum temperature corresponds to the DDT, and it is consistent with the transition time t_DDT we have recorded in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We plot the central densities of the DM and NM in Figure 8. The presence of DM tends to make the NM expand slower, which is similar to what we have found for the central and maximum temperatures. The sudden increases in the DM and NM central densities lead to the second peaks of the central and maximum temperature. This is because a converging acoustic wave is generated during the DDT. The acoustic wave travels toward the center of the NM, compressing and heating the material.

Zoom In Zoom Out Reset image size

We also tabulate the DM and the NM energies in Table 2 and show their variations in Figure 9. The energy E_i of the ith component of the fluid is defined by

$$\begin{eqnarray}&&{E}_{i}=\sum _{\mathrm{All}\ \mathrm{Grids}}{\rho }_{i}({\epsilon }_{i}+{\phi }_{j}+\displaystyle \frac{1}{2}({v}_{i}^{2}+{\phi }_{i})),\end{eqnarray} \tag{ 11 }$$

where j ≠ i. The extended DM component has negative energy at the end of the simulations, for all of our considered models, which means that they remain bounded. Figure 9 shows that the NM and DM energies would eventually approach constants toward the ends of the simulations. The increase in the DM energy after all major exothermic nuclear reactions are completed is because the contribution from the NM gravitational potential is weaker when the NM expands to a larger size. It hints at the decoupling between some of the most energetic NM and all of the bound DM. The most energetic NM is governed by the explosion timescale, which is much shorter than the dynamical timescale of the DM. They can no longer transfer energy to unbind the remaining DM through gravitational interaction.

Zoom In Zoom Out Reset image size

We tabulate the time-integrated total thermoneutrino production and the neutrino energy loss in Table 3. We find that models having more DM admixture tend to have a higher neutrino energy loss and more neutrino production in the 1–5 MeV channels. We plot the neutrino energy loss rate and the differential production rate for the 1 and 2 MeV channels in Figures 10–12. The sudden peaks in all neutrino observables at the transition time t_DDT are called neutrino bursts and they have been observed in several previous studies (Wright et al. 2016, 2017). The neutrino burst signals are weaker for models that are admixed with DM, but the neutrino production for these models is compensated by the prolonged deflagration phase. This observation is consistent with what we have found for the production of the ⁵⁶Ni and energy. We sum up the contributions of the thermoneutrinos from the 1–5 MeV channels and compute the ratio of the total neutrino production of the DM-admixed models to that of the pure NM one. The total neutrino production for the DM-admixed models could at most increase by a factor of 2.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is an exception to the above discussion, which is the model DM-1. The neutrino observables of this model are weaker when compared with those of the pure NM one. As we have discussed above that the DM gravitational potential traps the NM, which helps to suppress the thermal expansion of the hot NM. However, the DM mass in DM-1 is only 0.067 M_⊙. The gravitational potential of the DM is not deep enough to trap the NM (see Figure 6). On the other hand, the initial NM mass is substantially decreased by 0.132 M_⊙. As a result, there is a much lower amount of NM contributing to the production of thermoneutrinos at a similar temperature, and hence the neutrino signal is substantially reduced.

After we terminate the simulation, we map the isotope, density, internal energy, electron fraction, and velocity profiles into a uniform Lagrangian grid. We then input them into the SNEC code. Since the hydrodynamic simulation shows that the DM remains bound while the most energetic NM component is decoupled and ejected, the NM has a typical size much larger than that of the DM at the end of the explosion. To approximate the gravitational effects from the DM, we include the DM gravitational field as a point source g = − GM_DM/r² in the momentum equation of the SNEC code for all of the DM-admixed models.

The luminosity, radius, and temperature of the photosphere are plotted in Figures 13 and 14. We find that the light curves of the DM-admixed models are similar in peak luminosities when compared with that of the pure NM one. In particular, the peak luminosity of the model DM3 is only two times lower. This could be understood, as the amounts of ⁵⁶Ni generated are similar. The peak luminosity decreases as more DM is admixed, even though the model DM3 generates slightly more ⁵⁶Ni than the model DM2. Models that have more DM-admixed produce light curves that are dimmer and slower in rising and decline. These light curves are therefore flatter and broader. Note that in general, there exists a point of inflection for a normal thermonuclear supernova light curve. We find that such a point disappears in the range of 0–100 days for those models that have a more massive DM component, such as DM2 and DM3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

DM-admixed models tend to give a higher luminosity, say after 60 days of the explosions. The increased luminosity can be explained by the radii of their photospheres. For the pure NM model, the photosphere shrinks to the center at around 50 days, earlier than all DM-admixed models. The radii of the photospheres even increase to a second peak for the models DM2 and DM3, contributing to their higher luminosities and effective temperatures.

Admixing DM of 0.1 GeV to a WD could give rise to some interesting equilibrium structures. The DM component no longer acts as a stationary point mass, but it has an extended radius comparable to that of the NM. The subsequent supernova depends on the competition between the NM and DM components. If the DM is massive enough, the extra DM–NM gravitational interaction would have a significant impact on the expansion of the NM during the deflagration phase. It traps and keeps the NM hotter and denser, hence producing more fractions of ⁵⁶Ni. Though the model DM3 has less NM content compared to that of the model DM2, more ⁵⁶Ni is produced in the model DM3 than DM2. What if the DM content increases even further? Would it be possible to synthesize more ⁵⁶Ni than that of the pure NM model? We note that increasing the DM mass would reduce the NM mass. Hence, there should be an upper bound of the amount of ⁵⁶Ni being produced. Furthermore, the DM content may become so massive that no NM could eventually escape from the deep gravitational potential. The explosions would then behave like a failed detonation (Kasen & Plewa 2007; Plewa 2007; Jordan et al. 2012) even though the DDT occurs. It is an interesting future work to model the thermonuclear flame reignition of such an explosion.

We show that the DM is left behind as a compact dark star with a mass ranging from ∼0.07–0.5 M_⊙ in all of our considered models. This suggests an alternative way to search for astrophysical DM through thermonuclear supernovae—to look for any dark compact remnant in a thermonuclear supernova, through the microlensing effect for example. We note that recent development in gravitational-wave astronomy has opened up a new window to search for astronomical compact objects, especially sub-solar-mass black holes (Shandera et al. 2018; Nitz & Wang 2021). Our results show that if sub-solar-mass dark gravitational sources are observed, they could be compact dark stars remaining after DM-admixed thermonuclear supernovae. ⁷

Models with a more massive DM component tend to produce more thermoneutrinos. It would be interesting to estimate the number of expected neutrino events for the models considered in this work. It was shown that the number of events N scale as (Cavanna et al. 2004)

$$\begin{eqnarray}&&N={N}_{t}\times {\sigma }_{\nu }\times \displaystyle \frac{{\dot{N}}_{\nu }}{4\pi {d}^{2}}\times \tau ,\end{eqnarray} \tag{ 12 }$$

where N_t is the number of targets, σ_ν is the neutrino interaction cross section, ${\dot{N}}_{\nu }$ is the number of neutrinos produced per second, d is the distance of the supernova event from Earth, and τ is the time elapsed. There are previous results on estimating the event counts for some neutrino detectors proposed or under construction (Odrzywolek & Plewa 2011). They include the Hyper-Kamiokande (Proto-Collaboration et al. 2018), Memphys (Autiero et al. 2007; Patzak & the MEMPHYS Collaboration 2012), Glacier (Autiero et al. 2007; Rubbia 2009), and Low Energy Neutrino Astronomy (LENA; Autiero et al. 2007; Wurm et al. 2015). Here, we use Equation (12) to estimate the ratios in Table 3. The total number of neutrino events is presented in Table 4. In general, the neutrino production is distinguishable from that of the pure NM model only for the models DM2 and DM3. We have not done post-processing on the weak interaction neutrinos due to the limited number of isotopes. This together with a more detailed analysis in nucleosynthesis using a full nuclear network will be an interesting research direction in the future.

Although the amount of ⁵⁶Ni generated for the model DM3 is more than that of model DM2, the corresponding peak luminosity is lower. Moreover, the light curves for models that have a massive DM component (DM1–DM3) are flatter and broader, and the corresponding points of inflection disappear. It is interesting to investigate factors governing the peak luminosity and the shape of the light curve. In particular, the major difference between the models DM2 and DM3 is the amount of the DM remnant. To explore the effect of the DM gravity, we add one more light curve for each of the models DM0–DM3 in Figure 15. They are labeled as the model "NoDM," and the DM gravity is not included in the momentum equation of the SNEC code. We find that the presence of the DM gravitational force significantly alters the peak luminosity and the shape of the light curve. After switching off the DM gravity, the resulting light curves become narrower, brighter, and their points of inflection re-appear, especially for the models DM2 and DM3. Thus, we conclude that the effects of the DM gravity cannot be neglected in the light-curve calculation.

Zoom In Zoom Out Reset image size

Although the NM in all of the models have positive energies at the end of the simulation, some NM at the tail of the ejecta may remain bound, especially in the zone that is overlapping with the DM remnant. Nevertheless, we are ultimately interested in the supernova observables, and the amount of the free NM is important only in calculating the light curves using an analytic approach (Dado & Dar 2015). In our studies, we have already added the effect from the DM consistently in the SNEC's momentum equation, and so the SNEC code would keep track of the falling NM, which makes estimating the amount of unbound NM mass not necessary. Moreover, the fate of the NM in the post-explosion phase is governed by the ⁵⁶Ni and ⁵⁶Co decays and the γ-ray deposition. The kinetic energy transfer from the photons would also help to increase the kinetic energy of the NM so that it can escape the DM gravitational potential. However, modeling the radiative transfer of such a model will require a more realistic two-fluid radiative transfer code, involving complex gamma-ray heating in the bound core, which is beyond the scope of this study.

Peculiar thermonuclear supernovae had been observed in the past few decades. Depending on their luminosities and spectra, they can be further classified as the super-luminous (e.g., SN 1991T Phillips et al. 1992), subluminous (e.g., SN 1991bg Filippenko et al. 1992), and subluminous with strong mixing (also known as Type Iax supernovae). Among many observed Type Iax supernovae, SN 2002cx is one of the classical examples (Li et al. 2003). In particular, some of these events show a slower decline rate than that of a normal supernova. We extract the R-band light curves ⁸ of the following supernova events: SN 1999ac (Candia et al. 2003; Ganeshalingam et al. 2010; Silverman et al. 2012), SN 2002cx, SN 2002es (Ganeshalingam et al. 2012), SN 2005hk (Chornock et al. 2006; Phillips et al. 2007; Silverman et al. 2012; Stahl et al. 2019), and SN 2012Z (Silverman et al. 2012; Stahl et al. 2019), from The Open Supernova Catalog (White et al. 2015), and we show them together with the R-band light curves of our models in Figure 16. We find that some of the events, for instance, SN 1999ac and SN 2002cx are well fitted by the model DM1 before 15 days post-R-band maximum, ⁹ suggesting that other than some non-Chandrasekhar models, these two events may also be explained by thermonuclear supernovae having DM admixtures.

Zoom In Zoom Out Reset image size

Recent studies using observational data from the Palomar Transient Factory (PTF; Law et al. 2009) also found some peculiar thermonuclear supernovae that have a lower luminosity, longer rise time, and a slower decline rate (White et al. 2015). Such features resemble the drop of the peak luminosity and the flattened light curve in our DM-admixed models. We plot the peak magnitudes of the R-band light curve ¹⁰ against the decline rate Δm₁₅, for the PTF-observed supernovae: PTF 09ego, PTF 09eiy, PTF 10xk, PTF 11hyh, iPTF 13an, PTF 10ujn, and PTF 10acdh, together with those of our models in Figure 17. Our models show a wide range of peak magnitudes ranging from −17.8 to −18.8 and decline rates from 0.28–1.26. We find that the data of PTF 09eiy are closest to the best-fit line for the series of the DM-admixed models in Figure 17, followed by PTF 09ego, PTF 10acdh, and PTF 10ujn. These events seem to be possible candidates for DM-admixed thermonuclear supernovae.

Zoom In Zoom Out Reset image size

We show that DM-admixed models give rise to orthogonal deviations from the Phillips relation. The prevailing view on thermonuclear supernovae is that they are a group of diverse objects corresponding to multiple explosive models and progenitors, e.g., the single-degenerate and Chandrashekhar scenario. Even the single-degenerate scenario fails to produce the entire observed width–luminosity relation (Sim et al. 2013; Ohlmann et al. 2014). Therefore, not only the observed normal supernovae but also peculiar supernovae that deviate from the mainstream of the width–luminosity relation are now believed to be caused by other explosion mechanisms, such as the double-degenerate and sub-Chandrashekhar scenarios (Taubenberger 2017). We note that there is a recent simulation study on double-detonation explosion models and non-LTE light-curve calculations (Shen et al. 2021). Their theoretical B-band data that are generated by progenitors of different masses and compositions could span the observed B-band Phillips relation. This indicates that non-peculiar supernovae can correspond to explosive models in the sub-Chandrasekhar mass regime, and it also implies that other progenitors or explosion mechanisms are needed to account for peculiar supernovae. We note that several PTF entities presented in Figure 17 deviate significantly from the Phillips relation, being slightly dimmer but declining very slowly. We show in Appendix B that most canonical models, including single-degenerate and double-degenerate models, cannot account for this observed feature. DM-admixed models may provide an alternative path to explain some of these peculiar events and to span the parameter space that is not reachable by canonical models. However, we do not intend to claim that all peculiar supernovae are caused by a single-degenerate progenitor admixed with DM, particularly since we have only done a preliminary analysis through the bolometric light curve here. Instead, we propose the DM-admixed model as an alternative to span the possibilities for some peculiar supernovae.

Besides the 1D modeling for the explosive dynamics, several limitations of this work should be noted. We constructed a series of thermonuclear supernova models based on a model without DM admixture. To match the diversified thermonuclear supernova data, a wide range of normal models, coupled with different amounts of DM admixture, are needed to span the parameter space. Given that the effects of the DM observed in this work are generic, we expect that when DM is admixed with other ordinary thermonuclear supernova models, the trend will be similar.

Another caveat is that the Type Iax supernova spectra show a strong mixing feature and weak intermediate-mass element lines. Our DM-admixed models show that the delayed DDT can strongly suppress the strength of the detonation. This suggests that less intermediate-mass elements are produced. However, a firm conclusion will require a combination of radiative transfer coupled with an extended network of nuclear reactions to capture the explosive nucleosynthesis, which will be an interesting project in the future.

Last but not least, we note that our calculations do not include a full radiation transport, and the multicolor light curve is unclear. It would be interesting to investigate the post-explosion light curve using a multiband radiative transfer model. We have shown that including the gravitational force from the DM in the radiative transfer phase is necessary for more accurate treatment.