The Impact of Effective Matter Mixing Based on Three-dimensional Hydrodynamical Models on the Molecule Formation in the Ejecta of SN 1987A

Masaomi Ono; Takaya Nozawa; Shigehiro Nagataki; Alexandra Kozyreva; Salvatore Orlando; Marco Miceli; Ke-Jung Chen

doi:10.3847/1538-4365/ad1a08

1. Introduction

More than 30 yr have passed since the discovery of supernova 1987A (SN 1987A) in the LMC on 1987 February 23; it has entered a young supernova remnant (SNR) phase. Therefore, SN 1987A provides a unique opportunity to study the evolution of a core-collapse supernova (CCSN) from the explosion to an early phase of the SNR. Inside the ejecta, as the temperature goes down mainly due to the expansion, molecules and dust could form during this transition phase to a young SNR. Since SN 1987A is the last nearby supernova (SN) and the inner ejecta are spatially resolved (as mentioned below, spatially resolved distribution of the emission from molecules and dust has already been observed), SN 1987A is also a unique target for the study of molecule and dust formation in CCSNe.

Actually, in infrared spectra of early observations of SN 1987A, carbon monoxide (CO; Catchpole et al. 1987; McGregor et al. 1987; Oliva et al. 1987; Spyromilio et al. 1988; Meikle et al. 1989; Wooden et al. 1993) and silicon monoxide (SiO; Aitken et al. 1988; Roche et al. 1989, 1991, 1993; Bouchet & Danziger 1993; Wooden et al. 1993) have been detected as the emission of rotational–vibrational (hereafter, rovibrational) transitions. The light curves of the CO fundamental (∼4.6 μm; Bouchet & Danziger 1993) and the first overtone (∼2.3 μm; Bouchet & Danziger 1993; Meikle et al. 1993) bands and the SiO fundamental (∼8 μm) band (Bouchet & Danziger 1993) have also been delineated. From the observations above, the masses of CO and SiO have been derived approximately as 10⁻⁵–10⁻⁴ M_⊙ depending on the age and the model (Oliva et al. 1987; Spyromilio et al. 1988; Meikle et al. 1989) and 4 ± 2 × 10⁻⁶ M_⊙ (Roche et al. 1991, 1993), respectively. The mass of CO has been reinterpreted as 10⁻³ M_⊙ with a non-LTE (NLTE) and an optically thick assumption (Liu et al. 1992), and later approximately as (0.3–4) × 10⁻² M_⊙ (Liu & Dalgarno 1995). Similarly, the mass of SiO was also reinterpreted as 10⁻⁴ M_⊙ (Liu & Dalgarno 1994). CO has also been detected in other SNe such as Type II SN 1995ad (Spyromilio & Leibundgut 1996), Type IIn SN 1998S (Gerardy et al. 2000; Fassia et al. 2001), Type Ic SN 2000ew (Gerardy et al. 2002), Type II-P SN 2004et (Kotak et al. 2009, SiO was also detected), and Type Ib/c SN 2013eg (Drout et al. 2016). CO has been found even in the young SNR Cassiopeia A (Cas A) with the age of >300 yr, in which CO may have survived in the ejecta without significant destruction (Rho et al. 2009, 2012). For a more complete list of the detections of molecules in SNe and SNRs earlier than 2017 including the ones by the emission from rotational transitions, see Table 1 in Matsuura (2017).

There have been theoretical studies of chemistry in the ejecta of CCSNe. Motivated by the early detections of CO and SiO in SN 1987A, the formation of molecules in the core of SN 1987A has been studied with a chemical reaction model (network) mainly focusing on CO (Petuchowski et al. 1989; Lepp et al. 1990; Rawlings & Williams 1990). Later, Liu & Dalgarno (1995) modeled CO vibrational transitions based on thermal–chemical calculations with the optically thick regime (already mentioned above), in which time-dependent chemical equations were solved by simultaneously solving the thermal balance with the effects of heating due to the energy deposition from the decay of ⁵⁶Co and cooling via CO vibrational transitions; the observed light curves of the fundamental (Δv = 1; here, v denotes the vibrational level) and the first overtone (Δv = 2) bands were successfully reproduced (Liu & Dalgarno 1995), although an arbitrary time-dependent factor to control the energy deposition efficiency is necessary; even now, there has been no other study that reproduces the CO light curves. There have been similar studies on the chemistry targeted at SiO (Liu & Victor 1994; Liu & Dalgarno 1996). Apart from SN 1987A, the chemistry, i.e., the formation and/or destruction of molecules, clusters, and dust, in the ejecta of primordial (Population III) SNe has been investigated with comprehensive (large) chemical reaction networks in a series of works (Cherchneff & Lilly 2008; Cherchneff & Dwek 2009, 2010). The chemistry of molecules and dust in the ejecta of Cas A (Type IIb SN; Biscaro & Cherchneff 2014, 2016) and Type II-P SNe (Sarangi & Cherchneff 2013, 2015) has also been explored with similar approaches based on chemical reaction networks. Recently, Sluder et al. (2018) revisited the chemistry in the ejecta of SN 1987A with a large chemical reaction network including the nucleation of clusters, and comprehensively applied a dust formation theory. Liljegren et al. (2020) studied the formation of CO with the cooling via rovibrational transitions of CO.

It should be noted that, in all the studies mentioned in the previous paragraph, one-zone or multizone ejecta models have been used; some models include the effects of matter mixing (see the next paragraph) with an ad hoc assumption, i.e., a uniform composition throughout the core/a specific zone. There has been no study on molecule formation in CCSN ejecta that takes into account the effect of matter mixing based on three-dimensional (3D) hydrodynamical models. Additionally, as pointed out in several studies, chemical reactions involved with He⁺, Ne⁺, Ar⁺ could play a role in the destruction of molecules (Lepp et al. 1990; Liu & Dalgarno 1996; Cherchneff & Dwek 2009), which means that how He, Ne, and Ar are mixed in the core and ionized by high-energy electrons produced via the decay of ⁵⁶Co could potentially be crucial for the survival of molecules.

Matter mixing has been evoked by several early observations of SN 1987A, i.e., the early detections of hard X-ray emission (Dotani et al. 1987; Sunyaev et al. 1987), direct gamma-ray lines from the decay of ⁵⁶Co (Matz et al. 1988; Varani et al. 1990), the fine structure in H_α line (Hanuschik et al. 1988), and [Fe ii] lines with high-velocity tails (∼4000 km s⁻¹; Haas et al. 1990); some mixing of radioactive ⁵⁶Ni to be conveyed into high-velocity outer layers, which consist of helium and hydrogen, is necessary. Motivated by the failures of early numerical studies based on spherical CCSN explosions (Arnett et al. 1989a; Hachisu et al. 1990, 1992; Fryxell et al. 1991; Herant & Benz 1991, 1992; Mueller et al. 1991), which have mainly focused on the impact of fluid instabilities such as RT instability, the impact of nonspherical CCSN explosions (Yamada & Sato 1991; Nagataki et al. 1998; Kifonidis et al. 2000, 2003, 2006; Nagataki 2000; Hungerford et al. 2003, 2005; Joggerst et al. 2009, 2010a, 2010b; Gawryszczak et al. 2010; Hammer et al. 2010; Ellinger et al. 2012; Ono et al. 2013; Mao et al. 2015; Wongwathanarat et al. 2015, 2017; Gabler et al. 2021) on the matter mixing has been studied based on multidimensional hydrodynamical simulations. It was shown that a bipolar-like explosion works to produce a large amount of ⁴⁴Ti, which has the advantage to explain the bright light curve of SN 1987A in the late phase (Nagataki et al. 1997; Nagataki 2000). From the point of view of the mechanisms of CCSN explosions, multidimensionality in the explosion due to neutrino-driven convection, standing accretion shock instability (SASI; Blondin et al. 2003), and magnetorotational effects are essential for a successful explosion (for reviews of the mechanisms of CCSN explosions, see Janka 2012; Janka et al. 2012; Kotake et al. 2012a, 2012b; Burrows 2013; Müller 2016). Gabler et al. (2021) investigated the so-called "Ni bubble" effect (e.g., Basko 1994), i.e., the acceleration and inflation of the ejecta by the heating due to the radioactive decay of ⁵⁶Ni (⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe), with 3D hydrodynamical models with an approximated model for the energy deposition. 3D hydrodynamical models with 3D matter mixing have also been applied to other CCSN objects, Cas A (Wongwathanarat et al. 2017), Crab (Stockinger et al. 2020), and SN 2007Y (van Baal et al. 2023). Among the studies above, Wongwathanarat et al. (2015) pointed out the importance of the density structure of the progenitor stars on the matter mixing. Actually, there have been debates (Arnett et al. 1989b) on the progenitor star of SN 1987A, Sanduleak −69° 202 (Sk −69° 202), which is known as a blue supergiant (BSG; Walborn et al. 1987; West et al. 1987), since it had been difficult to explain observational features of Sk −69° 202, e.g., the effective temperature, luminosity, anomalous abundance in the nebula, and how to form the well-known triple-ring structure, by a single-star evolution model without fine-tuning parameters. Recently, progenitor star models, which satisfy all the observational features of Sk −69° 202, based on binary mergers have been developed (Menon & Heger 2017; Urushibata et al. 2018), following the binary merger scenarios (Podsiadlowski et al. 1990, 1992; Morris & Podsiadlowski 2007, 2009). A series of theoretical studies on the optical light curves of SN 1987A (Utrobin et al. 2015, 2019, 2021; Menon et al. 2019) has shown that binary merger progenitor models better explain the observed light curves than single-star progenitor models. On the basis of part of the studies above, we performed (Ono et al. 2020) 3D hydrodynamical simulations of globally asymmetric (bipolar-like) SN explosions with progenitor models including a binary merger progenitor model (Urushibata et al. 2018). As a result, an explosion model with the binary merger model succeeds to explain the observed features of [Fe ii] lines (Haas et al. 1990) compared with other explosion models with single-star progenitor models at least in the paper. Then, to link such aspherical explosions to the consequential observables in early SNR phases, further evolutions were followed by 3D magnetohydrodynamical (MHD) simulations (Orlando et al. 2020); it was found that the explosion model with the binary merger progenitor model also well explains the observed X-ray spectra, morphology, and light curves (e.g., Frank et al. 2016; Miceli et al. 2019). Additionally, as discussed in Orlando et al. (2020), thanks to the modeled bipolar-like explosion, the model is also able to reproduce the Doppler shift and broadening of the observed ⁴⁴Ti lines with NuSTAR (Boggs et al. 2015) and the distributions of CO and SiO inferred from the observations by Atacama Large Millimeter/submillimeter Array (ALMA; Abellán et al. 2017; see the next paragraph).

Recent breakthrough observations of SN 1987A by ALMA (Abellán et al. 2017) have revealed the 3D distribution of emission lines from rotational transitions of CO and SiO molecules in the ejecta. The distribution is not spherically symmetric and is rather complex; in particular, the distribution of CO has a torus-like structure perpendicular to the observed equatorial ring in the nebula of SN 1987A. It is noted that the approximated distributions of CO and SiO with the 3D MHD model above (see Figure 10 in Orlando et al. 2020) are consistent with the sizes of the observed molecular-rich structures and the torus-like distribution of CO and its orientation to the equatorial ring; in the 3D model, the bipolar-like explosion axis is perpendicular to the CO torus. The consistency with the ALMA observations further supports the progenitor model, asymmetric explosion, and matter mixing demonstrated in the 3D MHD model (Orlando et al. 2020). Further observations by ALMA have also delineated the distributions of the emission from both dust and the two molecules above (Cigan et al. 2019). Interestingly, from the latter ALMA observations, a blob (hotspot) of dust emission has been found, and it has been interpreted that the required extra heating to maintain the blob may be attributed to the emission from the putative neutron star of SN 1987A (NS 1987A), which has not directly been detected yet regardless of more than 30 yr of searches (e.g., Alp et al. 2018). Recent spectral fittings of X-ray observations of SN 1987A by Chandra, NuSTAR, and XMM-Newton (Greco et al. 2021, 2022) have suggested a nonthermal component, which probably stems from a pulsar wind nebula activity, which has further supported the existence of NS 1987A (for an alternative interpretation of the X-ray data, see Alp et al. 2021). There have also been attempts to constrain the properties of NS 1987A and to predict its evolution (Page et al. 2020; Dohi et al. 2023). Observations mentioned above have further motivated us to investigate the impacts of aspherical explosions and matter mixing on the formation of molecules in the ejecta of SN 1987A. Moreover, SN 1987A is a target of the newly launched James Webb Space Telescope (JWST) to figure out features of shocked dust grains (approved JWST proposal for Cycle 1, Matsuura et al. 2021). Recently, the first observational results of SN 1987A by JWST NIRSpec have been reported (Larsson et al. 2023). The ejecta distribution traced by [Fe i] line looks like a broken dipole. It is remarkable that the morphology and orientation are roughly consistent with the modeled bipolar-like explosions in our previous studies (Nagataki 2000; Ono et al. 2020; Orlando et al. 2020).

In this study, we investigate the impact of effective matter mixing (before the molecule formation in the ejecta) on the formation of diatomic molecules approximately up to 10,000 days after the explosion, for the first time, based on 3D CCSN explosion models for SN 1987A (Ono et al. 2020). We construct a chemical reaction network to follow the formation of molecules in the ejecta. The thermal evolution of the ejecta is crucial for molecule formation. Then, heating by the radioactive decay of ⁵⁶Ni⁹ (and/or ⁵⁶Co) and cooling by CO rovibrational transitions are somehow phenomenologically introduced (in the calculations, adiabatic cooling is also taken into account). By comparing with early observations of the light curves of CO vibrational bands, the parameters involved with the phenomenological approaches for the heating and cooling are calibrated to some extent. A perfect reproduction of the CO light curves is out of the scope of this paper because of the lack of realistic treatments for radiative transfer, excitation, ionization, and energy depositions by the decay of ⁵⁶Ni. Instead, we focus on the qualitative impact of the matter mixing during the SN shock propagation inside the progenitor star, which determines the initial composition of the seed atoms. Before fully applying our methodology to the 3D hydrodynamical models (Ono et al. 2020; Orlando et al. 2020), we utilize one-zone models and 1D radial profiles derived from 3D SN explosion models (Ono et al. 2020) as the initial conditions in this paper. By comparing with the model results of spherical explosion cases and the ones with angle-specified 1D profiles, the qualitative impact of effective matter mixing on the molecule formation in the ejecta is investigated.

This paper is organized as follows. In Section 2, the methodology is described. In Section 3, the strategy and models for investigating the impact of matter mixing on the molecule formation in the ejecta of SN 1987A are delineated. Section 4 is dedicated to the results. In Section 5, several issues are discussed on the basis of the results. In Section 6, the study in this paper is summarized.

2. Method

In this study, the dimensionalities of the models are restricted up to 1D. To derive the initial conditions of the ejecta, i.e., the density, temperature, and composition of the seed atoms, 3D CCSN explosion models in one of our previous studies (Ono et al. 2020) are utilized. By averaging the last snapshots (∼1 day after the explosion) of the 3D models, the inner regions of the ejecta, which consist of seed atoms, are approximated, first, as one-zones for understanding overall trends and, then, as 1D radial profiles to investigate the impact of matter mixing. For the former case, the temperature and density evolutions are basically obtained by power laws. For the latter case, the evolutions are obtained by performing 1D (spherical) hydrodynamical simulations with the 1D profiles as the initial conditions. For both cases, the effects of extra heating and cooling on the thermal evolution of the ejecta are taken into account. Then, rate equations for chemical reactions are solved for the ejecta models. Here, the method is described in detail.

2.1. Chemical Reaction Network

Here, the chemical reaction network adopted in this paper is described. To follow the formation and destruction of (diatomic) molecules, rate equations are solved. Species taken into account in this paper are listed in Table 1: 11 atoms; 24 diatomic molecules; electrons; and 39 singly ionized ions; i.e., in total, 75 species are included. In this paper, we limit our discussion to the formation of molecules; the nucleation of clusters of dust, and grain growth are beyond the scope of this study. Since more than triatomic species are more related to the nucleation of clusters, we take into account diatomic molecules as molecular species¹⁰ for simplicity. All diatomic molecules considered in the previous theoretical studies on molecule formation in CCSNe (Cherchneff & Dwek 2009; Sarangi & Cherchneff 2013, 2015; Sluder et al. 2018) are included except for aluminum oxide (AlO) that cannot form from the seed atoms corresponding to the 19 nuclei taken into account in the 3D hydrodynamical models (Ono et al. 2020, in which aluminum was not included). As for ion species, considering the low ionization fraction in the environment in this paper (see Section 2.1.3), singly ionized species¹¹ corresponding to the adopted atoms and diatomic molecules are taken into account if the data for the rates of the reactions involved with those species are available.

Table 1. Included Species

Atoms	H, He, C, N, O, Ne, Mg, Si, S, Ar, Fe
Diatomic molecules	H₂, CH, C₂, CN, CO, CS, NH, N₂, NO,
	OH, O₂, MgO, MgS, SiH, SiC, SiN, SiO,
	Si₂, SiS, SO, S₂, FeO, FeS, Fe₂

Ions	e⁻, H⁻, C⁻, O⁻, H⁺, He⁺, C⁺, N⁺, O⁺,
	Ne⁺, Mg⁺, Si⁺, S⁺, Ar⁺, Fe⁺, ${{\rm{H}}}_{2}^{+}$ , HeH⁺,
	CH⁺, ${{\rm{C}}}_{2}^{+}$ , CN⁺, CO⁺, CS⁺, NH⁺, ${{\rm{N}}}_{2}^{+}$ , NO⁺,
	OH⁺, ${{\rm{O}}}_{2}^{+}$ , MgO⁺, MgS⁺, SiH⁺, SiC⁺, SiN⁺
	SiO⁺, Si₂ ⁺, SiS⁺, SO⁺, ${{\rm{S}}}_{2}^{+}$ , FeO⁺, FeS⁺, Fe₂ ⁺

Download table as: ASCII Typeset image

The rate equations to be solved are as follows

$\begin{eqnarray}&&\displaystyle \frac{{{dc}}_{i}}{{dt}}=\sum {F}_{i}-\,\sum {D}_{i},\end{eqnarray} \tag{ 1 }$

where c_i is the number density of species i. F_i and D_i are formation and destruction rates of species i, respectively. Here, the former (later) summation is taken for the reactions, which include species i in the products (reactants). In the case of two-body reactions, F_i for j + k → i + ⋯ reaction and D_i for i + j → k + ⋯ reaction are given by

$\begin{eqnarray}&&{F}_{i}={k}_{i}(T)\,{c}_{j}\,{c}_{k},\ \ \ {D}_{i}={k}_{i}(T)\,{c}_{i}\,{c}_{j},\end{eqnarray} \tag{ 2 }$

respectively, where k_i(T) is the temperature (T) dependent rate coefficient. c_i, c_j, and c_k are the number density of the reactants, i, j, and k, respectively. For the case of three-body reactions, F_i for j + k + l → i + ⋯ reaction and D_i for i + j + k → l + ⋯ reaction are given by

$\begin{eqnarray}&&{F}_{i}={k}_{i}(T)\,{c}_{j}\,{c}_{k}\,{c}_{l},\ \ \ {D}_{i}={k}_{i}(T)\,{c}_{i}\,{c}_{j}\,{c}_{k},\end{eqnarray} \tag{ 3 }$

respectively. The rate coefficient, k_i(T), is conventionally expressed by the so-called Arrhenius form,

$\begin{eqnarray}&&{k}_{i}(T)={A}_{i}{\left(\displaystyle \frac{T}{300\ {\rm{K}}}\right)}^{{\alpha }_{i}}\exp (-{\beta }_{i}/T),\end{eqnarray} \tag{ 4 }$

where A_i, α_i, and β_i are specific coefficients for each reaction. Given the species listed in Table 1, all the reactions for the types listed in Table 2 (other than TF, CM, and UV) in which the coefficient data are available in the UMIST database¹² for astrochemistry (McElroy et al. 2013) or in the literature (Woodall et al. 2007; Cherchneff & Dwek 2009, 2010; Sluder et al. 2018) are included in the network. The adopted coefficient values for each reaction and the corresponding reference are listed in Table A1 in Appendix A. Other than the reactions above, thermal fragmentation reactions of diatomic molecules via thermal collisions with arbitrary particles (TF reactions), ionization, or destruction reactions by Compton electrons that stem from the decay of radioactive ⁵⁶Ni and ⁵⁶Co (CM reactions), and dissociation by UV photons (UV reactions) are also taken into account (see Sections 2.1.1 and 2.1.2).

Table 2. Types of Chemical Reactions

Code	Reaction Type
`3B`	Three-body reaction
`CE`	Charge exchange reaction
`IN`	Ion–neutral reaction
`MN`	Mutual neutralization reaction
`NN`	Neutral–neutral reaction
`AD`	Associative electron detachment
`DR`	Dissociative recombination
`RA`	Radiative association
`REA`	Radiative electron attachment
`RR`	Radiative recombination
`CD`	Collisional dissociation

`TF` ^a	Thermal fragmentation reaction
`CM` ^b	Ionization/destruction by Compton electrons
`UV` ^b	Dissociation by UV photons

Notes.

^aReactions described in Section 2.1.1.^bReactions described in Section 2.1.2.

Download table as: ASCII Typeset image

The rate equations to be solved are essentially the same as the nucleosynthesis calculations performed in our previous studies (e.g., Ono et al. 2009, 2012); in this paper, the equations are iteratively solved by the same full-implicit method in the time integration but with Gauss elimination with lower-upper decomposition as a matrix solver, which is different from the one (a sparse matrix solver) used in our previous studies above, for solving the matrix robustly.

2.1.1. Thermal Fragmentation Reactions

For thermal fragmentation reactions (TF reactions in Table 2),¹³ we adopt a kinetic approach, i.e., it is assumed that particles whose kinetic energy is greater than the bond-dissociation energy of the molecule i can destruct the molecule by thermal collisions. Therefore, D_i is given by

$\begin{eqnarray}&&{D}_{i}={c}_{i}\displaystyle \sum _{j}\displaystyle \frac{1}{1+{\delta }_{{ij}}}{c}_{j}\,{S}_{{ij}}\,{G}_{{ij}}(T),\end{eqnarray} \tag{ 5 }$

where S_ij is the cross section for the collision between the molecule i and a particle j; S_ij is assumed to be

$\begin{eqnarray}&&{S}_{{ij}}=\pi {\left({a}_{0,i}+{a}_{0,j}\right)}^{2},\end{eqnarray} \tag{ 6 }$

where a_0,i is the radius of the species i. The values of the radii of species are listed in Tables 3 and 4. G_ij(T) corresponds to the mean relative velocity between the particles i and j in which it is assumed that the relative velocities follow the Maxwell–Boltzmann distribution. Then, G_ij(T) is given by

$\begin{eqnarray}\begin{array}{rcl}{G}_{{ij}}(T) & = & \displaystyle \frac{4}{{\left(2\pi {\mu }_{{ij}}\right)}^{1/2}{\left({k}_{{\rm{B}}}T\right)}^{3/2}}\\ & & \times \,{\displaystyle \int }_{{E}_{{\rm{b}},i}}^{\infty }E\exp \left(-\displaystyle \frac{E}{{k}_{{\rm{B}}}T}\right){dE},\end{array}\end{eqnarray} \tag{ 7 }$

where μ_ij is the reduced mass of the two-body system of particles i and j, k_B is the Boltzmann constant, E denotes the kinetic energy of the relative motion, and E_b,i is the bond-dissociation energy of the molecule i. E_b,i is given by a fitting form as below

$\begin{eqnarray}&&{E}_{{\rm{b}},i}(T)=\displaystyle \sum _{l=1}^{4}{\epsilon }_{i,l}{\left(\displaystyle \frac{T}{3000\,{\rm{K}}}\right)}^{l-1},\end{eqnarray} \tag{ 8 }$

where _i,l are the fitting coefficients listed in Table 4. In Table 4, the bond-dissociation energies at representative temperatures, T = 2000 and 4000 K, are also listed. It is noted that the factor, 1/(1 + δ_ij), in Equation (5) is introduced to avoid double counting of the rate in the case of i = j. Here, δ_ij is the Kronecker delta.

Table 3. Atomic Radii

Atom (i)	a_0,i
	(Å)
H	0.53
He	0.31
C	0.67
N	0.56
O	0.48
Ne	0.38
Mg	1.45
Si	1.11
S	0.88
Ar	0.71
Fe	1.56

Note. The values are taken from Clementi et al. (1967).

Download table as: ASCII Typeset image

Table 4. Radii and Values Related to the Bond-dissociation Energies of the Diatomic Molecules

Molecule (i)	a_0,i ^a	_i,1 ^b	_i,2 ^b	_i,3 ^b	_i,4 ^b	E_b,i ^c (T = 2000 K)	E_b,i ^c (T = 4000 K) ^c
	(Å)					(eV)	(eV)
H₂	0.668	4.4737	4.6595 (−1)	–2.0450 (−1)	2.8627 (−2)	4.7019	4.7992
CH	0.766	3.4906	4.4935 (−1)	–3.3945 (−1)	8.7186 (−2)	3.6651	3.6929
C₂	0.844	6.1532	1.8501 (−1)	–6.7833 (−2)	8.7819 (−3)	6.2490	6.3001
$\mathrm{CN}$	0.781	7.7735	4.7013 (−1)	–3.5102 (−1)	8.4160 (−2)	7.9558	7.9758
CO	0.744	11.1200	3.9010 (−1)	–2.0736 (−1)	6.1460 (−2)	11.3070	11.4180
CS	0.994	7.3673	2.9418 (−1)	–1.1130 (−1)	3.3346 (−2)	7.5238	7.6407
NH	0.687	3.2134	4.5679 (−1)	–2.4998 (−1)	5.5345 (−2)	3.4232	3.5092
N₂	0.706	9.7592	4.3238 (−1)	–2.8255 (−1)	9.4590 (−2)	9.9499	10.0580
NO	0.659	6.5114	3.8560 (−1)	–2.4091 (−1)	7.6626 (−2)	6.6841	6.7789
OH	0.638	4.3952	4.5843 (−1)	–2.2279 (−1)	4.8458 (−2)	4.6162	4.7253
O₂	0.605	5.1285	3.9141 (−1)	–2.5362 (−1)	5.8071 (−2)	5.2939	5.3372
MgO	1.467	3.5241	–1.9805 (−1)	2.2559 (−2)	4.1024 (−2)	3.4143	3.3974
MgS	1.551	2.9102	–2.4538 (−1)	1.2113 (−1)	1.2405 (−2)	2.8041	2.8277
SiH	1.149	2.9886	3.3233 (−1)	–1.2256 (−1)	2.0183 (−2)	3.1617	3.2617
SiC	1.186	4.6197	1.3931 (−2)	6.0514 (−2)	–8.0209 (−3)	4.6535	4.7268
$\mathrm{SiN}$	1.156	5.6667	3.9114 (−1)	–3.6496 (−1)	1.0959 (−1)	5.7977	5.7992
SiO	1.139	8.2581	2.8743 (−1)	–1.0153 (−1)	2.9247 (−2)	8.4132	8.5301
Si₂	1.399	3.2074	2.1689 (−2)	2.2195 (−2)	2.6906 (−3)	3.2325	3.2822
SiS	1.270	6.4142	2.1404 (−1)	–7.8795 (−3)	–1.3752 (−3)	6.5530	6.6823
SO	0.925	5.3678	3.6931 (−1)	–2.5105 (−1)	6.5632 (−2)	5.5219	5.5695
S₂	1.109	4.3814	3.1299 (−1)	–2.0727 (−1)	5.5606 (−2)	4.5145	4.5621
FeO	1.575	4.2540	3.5731 (−1)	–2.0600 (−1)	8.8265 (−2)	4.4268	4.5734
FeS	1.648	3.3050	3.2418 (−1)	–1.6523 (−1)	7.8517 (−2)	3.4710	3.6296
Fe₂ ^d	1.965	4.3814	3.1299 (−1)	–2.0727 (−1)	5.5606 (−2)	4.5145	4.5621

Notes.

^aApproximated radius of the molecule, i, calculated as ${a}_{0,i}={({a}_{0,j}^{3}+{a}_{0,k}^{3})}^{1/3}$ ${a}_{0,i}={({a}_{0,j}^{3}+{a}_{0,k}^{3})}^{1/3}$ with atomic radii in Table 3, where a_0,j and a_0,k are the atomic radii of the constituent atoms.^bCoefficients of the fitting formula for the bond-dissociation energies in Equation (8). The bond-dissociation energies are estimated from the changes in the enthalpies of formation (Δ_f H^°) in the reaction systems (https://janaf.nist.gov/).^cBond-dissociation energies at T = 2000 and 4000 K obtained by Equation (8).^dThe values for the bond-dissociation energy of Fe₂ are assumed to be same as S₂.

Download table as: ASCII Typeset image

2.1.2. Collisional Ionization and Destruction by Compton Electrons; Dissociation by UV Photons

Compton electrons produced by the decay of ⁵⁶Ni and ⁵⁶Co could collisionally ionize atoms and molecules and could destruct (dissociate) molecules (CM reactions in Table 2). Subsequent UV photons caused by the deposited energy could also destruct molecules (UV reactions in Table 2). Such ionization and destruction are taken into account in a similar way to ones in the literature (e.g., Cherchneff & Dwek 2009; Liljegren et al. 2020).

First, the local energy generation rates by the decay of ⁵⁶Ni and ⁵⁶Co, _Ni and _Co, respectively, are considered as described in Joggerst et al. (2009). _Ni is given by

$\begin{eqnarray}&&{\epsilon }_{\mathrm{Ni}}={\lambda }_{\mathrm{Ni}}\,{e}^{-{\lambda }_{\mathrm{Ni}}t}\,{q}_{\mathrm{Ni}}\,{X}_{\mathrm{Ni}}\ \ \ \mathrm{erg}\,{{\rm{g}}}^{-1}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 9 }$

where λ_Ni, t, q_Ni, and X_Ni are the decay rate of ⁵⁶Ni, the time after the explosion, the q-value of the decay from ⁵⁶Ni to ⁵⁶Co, and the local mass fraction of ⁵⁶Ni, respectively. The values of λ_Ni, q_Ni are 1.315 × 10⁻⁶ s⁻¹ and 2.96 × 10¹⁶ erg g⁻¹, respectively. _Co is given by

$\begin{eqnarray}&&\begin{array}{l}{\epsilon }_{\mathrm{Co}}=\displaystyle \frac{{\lambda }_{\mathrm{Ni}}}{{\lambda }_{\mathrm{Co}}-{\lambda }_{\mathrm{Ni}}}{\lambda }_{\mathrm{Co}}\left({e}^{-{\lambda }_{\mathrm{Ni}}t}-{e}^{-{\lambda }_{\mathrm{Co}}t}\right){q}_{\mathrm{Co}}\,{X}_{\mathrm{Ni}}\\ \,\,\,\,\,\,\mathrm{erg}\,{{\rm{g}}}^{-1}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 10 }$

where λ_Co and q_Co are the decay rate of ⁵⁶Co and the q-value of the decay from ⁵⁶Co to ⁵⁶Fe, respectively. The values of λ_Co, q_Co are 1.042 × 10⁻⁷ s⁻¹ and 6.4 × 10¹⁶ erg g⁻¹, respectively. Then, the local energy deposition rate is given as follows

$\begin{eqnarray}&&\epsilon =\left({\epsilon }_{\mathrm{Ni}}+{\epsilon }_{\mathrm{Co}}\right){D}_{\gamma }\ \ \ \mathrm{erg}\,{{\rm{g}}}^{-1}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 11 }$

where D_γ is the fraction of the energy deposition, which depends on the effective optical depth, τ_γ, to gamma rays that emerged from the decay as follows

$\begin{eqnarray}&&{D}_{\gamma }=1-\exp (-{\tau }_{\gamma }).\end{eqnarray} \tag{ 12 }$

The optical depth is obtained differently for one-zone and 1D cases. For the one-zone case, the optical depth is obtained as in Liu & Dalgarno (1995), i.e.,

$\begin{eqnarray}&&{\tau }_{\gamma }={\tau }_{0}{\left(t/{t}_{0}\right)}^{-2},\end{eqnarray} \tag{ 13 }$

where τ₀ is the optical depth at a reference time t₀. τ₀ and t₀ are 31.1 and 100 days, respectively. For the 1D case, the optical depth depends on the radius, r, and it is derived as

$\begin{eqnarray}&&{\tau }_{\gamma }(r)={\int }_{r}^{\infty }\displaystyle \frac{1}{\,{l}_{\gamma }\,}\,{dr},\end{eqnarray} \tag{ 14 }$

where l_γ is the effective mean free path of the gamma rays. The upper limit of the integration is practically taken to be the radius of the upper computational domain. l_γ is taken from Shigeyama & Nomoto (1990) as

$\begin{eqnarray}&&{l}_{\gamma }=1.67\times {10}^{13}\ {Y}_{{\rm{e}}}^{-1}{\left(\displaystyle \frac{\rho }{{10}^{-12}\ {\rm{g}}\,{\mathrm{cm}}^{-3}}\right)}^{-1}\ \ \ \mathrm{cm},\end{eqnarray} \tag{ 15 }$

where Y_e and ρ are the electron fraction (the total number of free and bound electrons per nucleon) and the density, respectively. l_γ corresponds to the effective distance that the emerged 1.28 MeV gamma rays travel before being degraded to 60 keV and absorbed (Shigeyama & Nomoto 1990).

It is noted that, as described above, nonlocal energy deposition is not taken into account. In reality, emitted energy can nonlocally be deposited over a region approximately corresponding to the length scale of the effective mean free path in Equation (15) (practically over several tracer particles; see Section 2.4), which may cause smoother gas temperature profiles across the optically thick-to-thin regions and may allow energy deposition outside the regions where ⁵⁶Ni exists. Therefore, the potential impact of nonlocal energy deposition is discussed with a simple model in Appendix B. The methodology to count the effects of the energy deposition (without nonlocal energy deposition) on thermal quantities of the gas is described in Sections 2.3 and 2.4.

Additionally, as described in Equation (14), the optical depth only in the radial direction, i.e., the escape of the emitted energies only in the radial direction, is considered. Optical depths in nonradial directions (e.g., Gabler et al. 2021) for counting nonradial energy escape, which may be effective if the distribution of ⁵⁶Ni is globally asymmetric, are not also taken into account since the hydrodynamical models presented in this study are assumed to be spherically symmetric, and the effect cannot be treated appropriately.

Given the energy deposition rate in Equation (11), the rate coefficients, k_C, of the ionization and destruction by Compton electrons are given as follows

$\begin{eqnarray}&&{k}_{{\rm{C}}}=\displaystyle \frac{\epsilon \,\rho }{W\,{n}_{\mathrm{tot}}}\ \ \ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 16 }$

where W is the mean energy per ion pair, and the values depend on the target and the reaction type. n_tot is the total number density of the particles. The values of W could also depend on the ionization fraction (see Figure 2 in Liu & Dalgarno 1995). However, we take the values as constant for simplicity. The constant values of W are taken from Cherchneff & Dwek (2009), Sluder et al. (2018). For the ionization and destruction reactions by Compton electrons, the reaction types of the forms described below are taken into account. The ionization of an atom X is as follows:

$\begin{eqnarray}&&{\rm{X}}+{{\rm{e}}}_{{\rm{C}}}^{-}\to {{\rm{X}}}^{+}+{{\rm{e}}}^{-}+{{\rm{e}}}_{{\rm{C}}}^{-},\end{eqnarray} \tag{ 17 }$

where ${{\rm{e}}}_{{\rm{C}}}^{-}$ denotes the Compton electron. The ionization of a diatomic molecule AB (A and B denote the constituent atoms) is as follows:

$\begin{eqnarray}&&\mathrm{AB}+{{\rm{e}}}_{{\rm{C}}}^{-}\to {\mathrm{AB}}^{+}+{{\rm{e}}}^{-}+{{\rm{e}}}_{{\rm{C}}}^{-}.\end{eqnarray} \tag{ 18 }$

The following are for the dissociative reactions:

$\begin{eqnarray}&&\begin{array}{l}\mathrm{AB}+{{\rm{e}}}_{{\rm{C}}}^{-}\to {\rm{A}}+{\rm{B}}+{{\rm{e}}}_{{\rm{C}}}^{-},\\ \mathrm{AB}+{{\rm{e}}}_{{\rm{C}}}^{-}\to {{\rm{A}}}^{+}+{\rm{B}}+{{\rm{e}}}^{-}+{{\rm{e}}}_{{\rm{C}}}^{-},\\ \mathrm{AB}+{{\rm{e}}}_{{\rm{C}}}^{-}\to {\rm{A}}+{{\rm{B}}}^{+}+{{\rm{e}}}^{-}+{{\rm{e}}}_{{\rm{C}}}^{-}.\end{array}\end{eqnarray} \tag{ 19 }$

The adopted values of the mean energy per ion pair W for the ionization and destruction reactions described above are listed in Table 5. Similarly, secondary UV photons destruct a diatomic molecule AB to produce A and B. The rate coefficient k_UV is given by

$\begin{eqnarray}&&{k}_{\,\mathrm{UV}}=\displaystyle \frac{\alpha \,\epsilon \,\rho }{{E}_{\mathrm{UV}}\,{n}_{\mathrm{tot}}}\ \ \ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 20 }$

where E_UV is the representative energy of UV photons. The value of E_UV is set to be the energy of a fiducial photon with the wavelength of 1302 Å as in Cherchneff & Dwek (2009). α is the fraction of the deposited energy to be converted to UV photons. For simplicity, the values of α at 200 and 1000 days are set to be 0.25 and 0.4, respectively, according to a discussion in Cherchneff & Dwek (2009). For an epoch between 200 and 1000 days, the value is linearly interpolated, and on the outside, the values are extrapolated, but arbitrarily, the maximum and minimum values are set to be 0.5 and 0.1, respectively. Both the Compton electrons and UV photons are not included as one of the explicit species in the chemical reaction network. For the reactions mentioned above, the destruction rate D_i in Equation (1) is obtained by regarding the reaction as a transition (one-body reaction); for example, the destruction rate D_i in Equation (1) for the case of the ionization of an atom X by Compton electrons is given by

$\begin{eqnarray}&&{D}_{{\rm{X}}}={k}_{{\rm{C}}}\,{c}_{{\rm{X}}}.\end{eqnarray} \tag{ 21 }$

As described above, the derived ionization and destruction rates are simply proportional to the number of available Compton electrons per particle estimated with the constant mean energy W; there may be large uncertainties in the derived rate coefficients. Therefore, to see the impact of such uncertainties, a constant factor f_d is introduced as a model parameter to control the efficiency of the rates as in Sluder et al. (2018). Then, the ionization and destruction rates in Equations (16) and (20) are practically introduced as f_d k_C and f_d k_UV in the chemical reaction calculations.

Table 5. Reactions Involved with Compton Electrons ${{\rm{e}}}_{{\rm{C}}}^{-}$ and the Adopted Mean Energy per Ion-pair W

Reaction (Type)	W (eV)
CO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → C + O⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	768
CO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → C⁺ + O + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	274
CO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → C + O + ${{\rm{e}}}_{{\rm{C}}}^{-}$	125
CO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → CO⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	34

SiO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → Si + O⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	678
SiO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → Si⁺ + O + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	218
SiO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → Si + O + ${{\rm{e}}}_{{\rm{C}}}^{-}$	110
SiO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → SiO⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	30

H₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → H⁺ + H + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	829
H₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → H + H + ${{\rm{e}}}_{{\rm{C}}}^{-}$	77
H₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → ${{\rm{H}}}_{2}^{+}$ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	37.7

N₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → N⁺ + N + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	264
N₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → N + N + ${{\rm{e}}}_{{\rm{C}}}^{-}$	133.5
N₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → ${{\rm{N}}}_{2}^{+}$ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	36.3

O₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → O⁺ + O + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	768
O₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → O + O + ${{\rm{e}}}_{{\rm{C}}}^{-}$	125
O₂ + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → ${{\rm{O}}}_{2}^{+}$ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	34

For oxygen involving molecules other than above
AO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → A + O⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	768
AO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → A⁺ + O + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	274
AO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → A + O + ${{\rm{e}}}_{{\rm{C}}}^{-}$	125
AO + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → AO⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	34

For molecules other than above
AB + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → A + B⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	274
AB + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → A⁺ + B + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	274
AB + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → A + B + ${{\rm{e}}}_{{\rm{C}}}^{-}$	125
AB + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → AB⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	34

Ionization of atoms
X + ${{\rm{e}}}_{{\rm{C}}}^{-}$ → X⁺ + e⁻ + ${{\rm{e}}}_{{\rm{C}}}^{-}$	47

Note. Here, A and B denote constituent atoms other than oxygen. X denotes an arbitrary atom.

Download table as: ASCII Typeset image

2.1.3. Initialization of the Chemical Reaction Calculations

Generally, if the gas temperature is higher than 10⁴ K, the destruction by thermal fragmentation reactions dominates the formation reactions of molecules. The transition temperature when molecules start to form, however, is not so clear. Then, in this study, the chemical reaction calculations are started if the gas temperature drops to a bit higher temperature than above, i.e., 1.5 × 10⁴ K, to be safe. The initial abundances of the species for the chemical reaction calculations are determined as follows. As the initial conditions, only the 11 atoms, singly ionized atoms, and electrons have nonzero abundances initially. The abundances of the atoms are determined to be consistent with the compositions derived from the 3D hydrodynamical calculation results (Ono et al. 2020). In Ono et al. (2020), in total, 19 nuclei were taken into account. Then, the abundances of the nuclei are assigned to the corresponding atoms depending on the isotopes. For example, the abundances of ⁵²Fe, ⁵⁴Fe, and ⁵⁶Ni (⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe) are assigned as the iron atom Fe. Stable isotopes are simply summed up to count the corresponding atomic abundances. It is noted that the electron fraction Y_e in Equation (15) does not depend on ionization and chemical reactions. Hence, the electron fraction is obtained from the original abundances of the 19 nuclei. The change in the electron fraction due to the radioactive decay is neglected since the change does not affect much the optical depth in Equation (14).

Before starting the molecule formation, some fraction of atoms could be ionized. Therefore, the ionization before starting the chemical reaction calculation is determined manually by introducing the ionization fraction¹⁴ X_e, i.e., the number of free electrons per all the nuclei. In the ejecta core of SN 1987A, X_e could range from ∼10⁻¹–10⁻² to 10⁻³–10⁻⁴ corresponding to from 200 to 2000 days after the explosion depending on the position in the core layers (see Figure 8 in Kozma & Fransson 1998). Therefore, two parameters, the initial (∼1 day after the explosion) ionization fraction X_e,i and the ionization fraction at 2000 days X_e,f, are introduced for simplicity. Then, X_e before 2000 days is obtained by linear interpolation of the powers of 10 with the values of X_e,i and X_e,f. X_e after 2000 days is set to be the same value at 2000 days. As the reference values of X_e,i and X_e,f, 10⁻¹ and 10⁻³ are adopted, respectively.

With the X_e just before the chemical reaction calculation, the initial abundances of the species are determined to be consistent with it; a fraction of X_e of the abundance of each atom is assigned to the abundance of the corresponding singly ionized atom. Then, the same fraction of the abundance is added to the electron abundance.

2.2. Cooling by Emission from Rovibrational Transitions of CO

Emission from CO rovibrational transitions could cool the ejecta gas and eventually affect the chemistry in the ejecta. To estimate emission lines of CO rovibrational transitions, rate equations for vibrational levels are solved. In this study, the rotational levels are assumed to be thermally populated. Here, we basically adopt the method described in Liu & Dalgarno (1995). The rate equations for vibrational levels (vibrational quantum number, $v=0,1,\ldots ,{v}_{\max }$ ) are as follows

$\begin{eqnarray}&&\displaystyle \frac{{{dX}}_{i}}{{dt}}=\displaystyle \sum _{j}\left({R}_{{ji}}{X}_{j}-{R}_{{ij}}{X}_{i}\right),\end{eqnarray} \tag{ 22 }$

where X_i and R_ij are the population of the vibrationally excited state i and the rate coefficient for the transition from the vibrational level i to j, respectively. Here, the rate coefficient is expressed as

$\begin{eqnarray}\begin{array}{rcl}{R}_{{ij}} & = & \overline{{A}_{{ij}}}+\overline{{B}_{{ij}}\,{J}_{{ij}}}+{n}_{{\rm{e}}}\,{q}_{{ij}},i\gt j,\\ & = & \overline{{B}_{{ij}}\,{J}_{{ji}}}+{n}_{{\rm{e}}}\,{q}_{{ij}},i\lt j,\end{array}\end{eqnarray} \tag{ 23 }$

where $\overline{{A}_{{ij}}}$ , $\overline{{B}_{{ij}}\,{J}_{{ij}}}$ are derived (explained later below) from the Einstein's A-coefficient A_{κ
λ} for spontaneous emission of the rovibrational transition; similarly, those are derived for the B-coefficient B_{κ
λ} for stimulated emission, and the mean radiation field J_{κ
λ} corresponding to the rovibrational transition, respectively. Here, Greek characters in subscripts denote rovibrational levels. Hence, a rovibrational level is specified by the vibrational level v and the rotational level, $J=0,1,\ldots ,{J}_{\max };$ the state κ can be expressed more specifically as κ = (v, J). The term, n_e q_ij, is the contribution from electron impact excitation/deexcitation. n_e is the number density of free (thermal) electrons. q_ij is the rate coefficient of the excitation (i < j) or deexcitation (i > j) of the vibrational level i to j (for the description, see the last paragraph in this section and Appendix C). Einstein's coefficients are followed by the two relations below

$\begin{eqnarray}&&{A}_{\kappa \lambda }=\displaystyle \frac{2h{\nu }_{\kappa \lambda }^{3}}{{c}^{2}}{B}_{\kappa \lambda },\ \ {E}_{\kappa }\gt {E}_{\lambda },\end{eqnarray} \tag{ 24 }$

$\begin{eqnarray}&&{g}_{\kappa }{B}_{\kappa \lambda }={g}_{\lambda }{B}_{\lambda \kappa },\end{eqnarray} \tag{ 25 }$

where ν_{κ
λ} is the frequency of the corresponding photon. h and c are the Planck constant and the speed of light, respectively. g_κ is the statistical weight of the state κ. E_κ is the energy of the rovibrational state κ.

As mentioned above, rotational levels are assumed to be thermally populated; the population probability of a rotational level J among the same vibrational level v is obtained as

$\begin{eqnarray}&&{P}_{(v,J)}=\displaystyle \frac{{g}_{(v,J)}\exp (-{E}_{(v,J)}/{k}_{{\rm{B}}}T)}{{\sum }_{J^{\prime} =0}^{{J}_{\max }}\,{g}_{(v,J^{\prime} )}\exp (-{E}_{(v,J^{\prime} )}/{k}_{{\rm{B}}}T)},\end{eqnarray} \tag{ 26 }$

where g_(v,J) and E_(v,J) are the statistical weight and the energy of the rovibrational state (v, J), respectively. T and k_B are the temperature and the Boltzmann constant, respectively. Then, $\overline{{A}_{{ij}}}$ are derived from the summation of the rovibrational Einstein's A coefficients over the transitions of the vibrational level i to j with the weight of the initial state's population probability as

$\begin{eqnarray}&&\overline{{A}_{{ij}}}=\displaystyle \sum _{J=0}^{{J}_{\max }}\displaystyle \sum _{J^{\prime} =0}^{{J}_{\max }}{P}_{(i,J)}{A}_{(i,J)(j,J^{\prime} )}.\end{eqnarray} \tag{ 27 }$

The mean radiation field J_{κ
λ} is obtained with the so-called Sobolev approximation (Sobolev 1960), in which as shown below the line optical depth is derived with only local variables under the condition that the velocity gradient is steep enough, as follows

$\begin{eqnarray}\begin{array}{rcl}{J}_{\kappa \lambda } & = & (1-{\beta }_{\kappa \lambda })\,\displaystyle \frac{2h{\nu }_{\kappa \lambda }^{3}}{{c}^{2}}{\left(\displaystyle \frac{{X}_{\lambda }/{g}_{\lambda }}{{X}_{\kappa }/{g}_{\kappa }}-1\right)}^{-1}\\ & & +\ {\beta }_{\kappa \lambda }\,{B}_{{\rm{d}}}({\nu }_{\kappa \lambda },{T}_{\mathrm{ph}}),\ \ {E}_{\kappa }\gt {E}_{\lambda },\end{array}\end{eqnarray} \tag{ 28 }$

where β_{κ
λ} is the escape probability. X_κ (X_λ) and g_κ (g_λ) are the population of the rovibrational state κ (λ) and the statistical weight of the state κ (λ), respectively. The escape probability β_{κ
λ} is obtained by

$\begin{eqnarray}&&{\beta }_{\kappa \lambda }=\displaystyle \frac{1-\exp (-{\tau }_{\kappa \lambda })}{{\tau }_{\kappa \lambda }},\end{eqnarray} \tag{ 29 }$

where τ_{κ
λ} is the line optical depth expressed as

$\begin{eqnarray}&&{\tau }_{\kappa \lambda }=\displaystyle \frac{{A}_{\kappa \lambda }\,{c}^{3}{n}_{\kappa }}{8\pi {\nu }_{\kappa \lambda }^{3}}\left(\displaystyle \frac{{X}_{\lambda }/{g}_{\lambda }}{{X}_{\kappa }/{g}_{\kappa }}-1\right)t,\end{eqnarray} \tag{ 30 }$

where n_κ is the number density of the state κ, and t is the time after the explosion. The second term of Equation (28) is the contribution from the background photons from the photosphere, and B_d is the Planck function diluted by a factor α_d (dilution factor);

$\begin{eqnarray}&&{B}_{{\rm{d}}}({\nu }_{\kappa \lambda },{T}_{\mathrm{ph}})=\displaystyle \frac{2h{\nu }_{\kappa \lambda }^{3}}{{c}^{2}}\displaystyle \frac{{\alpha }_{{\rm{d}}}}{\exp (h{\nu }_{\kappa \lambda }/{k}_{{\rm{B}}}{T}_{\mathrm{ph}})-1},\end{eqnarray} \tag{ 31 }$

where T_ph is the temperature at the photosphere. The dilution factor α_d is calculated as

$\begin{eqnarray}&&{\alpha }_{{\rm{d}}}=\displaystyle \frac{1}{2}\left(1-\sqrt{1-{\left({v}_{r,\mathrm{ph}}/{v}_{r}\right)}^{2}}\right),\end{eqnarray} \tag{ 32 }$

where v_r,ph is the radial velocity at the photosphere, and v_r is the local radial velocity. If v_r < v_r,ph, i.e., inside the photosphere, the contribution from the background photons is turned off, i.e., B_d(ν_{κ
λ}, T_ph) = 0, by practically setting the dilution factor to zero. The properties of the photosphere are unknown. Then, throughout this paper, the results of one of 1D radiative transfer calculations are adopted. With the angle-averaged 1D profiles based on the 3D hydrodynamical calculation for SN 1987A (the model b18.3-high; Ono et al. 2020; see Section 2.4) as the initial condition, the radiative transfer is calculated¹⁵ (A. Kozyreva et al. 2024, in preparation) by using a 1D radiation hydrodynamical code (STELLA; see, e.g., Blinnikov et al. 1998, 2006; Kozyreva et al. 2020). Figure 1 shows the photospheric temperatures, T_ph, and the radial velocities at the photospheres, v_r,ph, for several bands (U, B, V, R, I). Since the emission from rovibrational transitions (Δv = 1, and Δv = 2 transitions correspond to ∼4.6 and ∼2.3 μm, respectively) is observed outside the I band, we adopt the values of the I band as the reference. The properties of the photospheres are obtained only up to ∼170 days after the explosion. Therefore, the values after 170 days are crudely extrapolated with a power law (∼t⁻¹, see the dashed lines in Figure 1).

**Figure 1.** The photospheric temperatures (left) and the radial velocities at the photospheres (right) as a function of time after the explosion for the U, B, V, R, and I bands calculated by a 1D radiation hydrodynamical code with the angle-averaged 1D profiles based on the 3D hydrodynamical model (the model b18.3-high; Ono et al. 2020) as the initial condition. Arithmetic averages (mean) and arbitrary extrapolations (extend; see the text) are also plotted.
Download figure:
Standard image High-resolution image

Then, $\overline{{B}_{{ij}}\,{J}_{{ij}}}$ in the first line in Equation (23) is obtained similarly as in Equation (27), i.e.,

$\begin{eqnarray}&&\overline{{B}_{{ij}}\,{J}_{{ij}}}=\displaystyle \sum _{J=0}^{{J}_{\max }}\displaystyle \sum _{J^{\prime} =0}^{{J}_{\max }}{P}_{(i,J)}{B}_{(i,J)(j,J^{\prime} )}\,{J}_{(i,J)(j,J^{\prime} )}.\end{eqnarray} \tag{ 33 }$

Once the number density of the state κ, n_κ, is obtained by solving Equation (22), the line emissivity for the rovibrational transition from the state κ to λ is eventually obtained by

$\begin{eqnarray}&&\begin{array}{l}{j}_{\kappa \lambda }=\displaystyle \frac{h{\nu }_{\kappa \lambda }}{4\pi }{n}_{\kappa }\,{A}_{\kappa \lambda }\,{\beta }_{\kappa \lambda },\ \ \mathrm{erg}\,{\mathrm{cm}}^{-3}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1},\\ \,\,\,{E}_{\kappa }\gt {E}_{\lambda }.\end{array}\end{eqnarray} \tag{ 34 }$

The rate equations described in Equation (22) are solved with the same methodology for the rate equations in Equation (1).

As for the data necessary for the rovibrational transitions, i.e., Einstein's coefficients, energy levels (transition frequencies), and statistical weights, are taken from Li et al. (2015).¹⁶ The rate coefficients of the electron impact excitation/deexcitation of vibrational levels, i.e., q_ij, are derived from the values of the cross sections in Poparić et al. (2008) by taking Maxwellian averaging (see Appendix C). As for the maximum quantum numbers (levels) of vibrational and rotational states to be considered in this study, we set ${v}_{\max }=6$ , and ${J}_{\max }=128$ , respectively. It is noted that, in a recent theoretical study on the formation of CO in a CCSN ejecta (Liljegren et al. 2020), vibrational levels up to v = 6 were taken into account. We confirm that taking into account higher levels does not affect the results much. As for the initial population, it is assumed that only the grand state, i.e., v = 0, is initially populated; it is confirmed that changing the initial population does not affect the main results of this study. Once CO forms, Equation (22) is solved to obtain the population of (ro)-vibrational levels and the emissivity as in Equation (34). The calculation for the CO (ro)-vibrational levels is performed up to 1000 days after the explosion. After this epoch, there have been no observations of CO vibrational bands in SN 1987A, and the CO cooling would play only a minor role.

2.3. Calculations with the One-zone Approximation

In order to see the chemical evolution in the ejecta, first, for simplicity, we apply our method to one-zone models, which regard the inner core of the ejecta as one-zone. The initial condition of the ejecta is obtained based on the 3D hydrodynamical simulation results for SN 1987A (Ono et al. 2020). Among the models, we adopt the best model b18.3-high, i.e., an asymmetric bipolar-like explosion with a BSG progenitor star formed through a binary merger. The initial condition is obtained by averaging the physical quantities in a core region of the 3D model. First, an angle-averaged 1D profile is derived based on the last snapshot (∼1 day after the explosion) of the 3D model. Then, the physical quantities from the inner boundary to the mass coordinate of about 6 M_⊙ are averaged with the density as a weight for quantities other than the density. The obtained initial physical values, i.e., the masses of the atoms, the total mass, the density, temperature, and radial velocity, are listed¹⁷ in Table 6. The time evolutions of the density ρ and the specific internal energy e are followed by power laws with the powers of −3 and −3(γ_ad − 1) as below, respectively, assuming an adiabatic expansion,

$\begin{eqnarray}&&\rho (t)={\rho }_{0}{\left(t/{t}_{0}\right)}^{-3},\end{eqnarray} \tag{ 35 }$

$\begin{eqnarray}&&e(t)={e}_{0}{\left(t/{t}_{0}\right)}^{-3({\gamma }_{\mathrm{ad}}-1)},\end{eqnarray} \tag{ 36 }$

where γ_ad is the adiabatic index. t₀ is the time of the previous time step, and ρ₀ and e₀ are the density and specific internal energy at t = t₀, respectively. The temperature is derived through the ideal equation of state (EoS),

$\begin{eqnarray}&&P=\displaystyle \frac{R}{\mu }\,\rho \,T=e({\gamma }_{\mathrm{ad}}-1)\,\rho ,\end{eqnarray} \tag{ 37 }$

where P is the pressure, R is the gas constant, and μ is the mean molecular weight. If the radiation is dominant in the system, γ_ad = 4/3. On the other hand, if the monoatomic gas (diatomic gas) is dominant, γ_ad = 5/3 (7/5). At the age of the initial conditions (∼1 day after the explosion), the density of the inner core (see Table 6) is still high, and the radiation is dominant. As the expansion proceeds, the ejecta gas becomes optically thin, and a transition in the adiabatic index γ_ad occurs at some point as shown in the 1D radiation hydrodynamical calculation results mentioned in Section 2.2 (see Figure 2).

**Figure 2.** The gas temperatures of Lagrange meshes as a function of time after the explosion. The values are obtained from the 1D radiation hydrodynamical calculation results with the angle-averaged 1D profiles based on the 3D hydrodynamical model (the model b18.3-high; Ono et al. 2020). Each line shows the time evolution for each Lagrange mesh. The colors denote the corresponding initial positions (radii). The two straight lines represent the power-law evolutions with the powers of −1 and −2 corresponding to γ_ad = 4/3, and γ_ad = 5/3, respectively.
Download figure:
Standard image High-resolution image

Table 6. Initial Conditions for One-zone Calculations Obtained Based on the 3D Hydrodynamical Model b18.3-high (Ono et al. 2020)

Description	Parameter	Value	Unit
Hydrogen mass	M_H	1.28	M_⊙
Helium mass	M_He	1.66	M_⊙
Carbon mass	M_C	1.17 (−1)^a	M_⊙
Nitrogen mass	M_N	5.59 (−3)	M_⊙
Oxygen mass	M_O	3.21 (−1)	M_⊙
Neon mass	M_Ne	7.66 (−3)	M_⊙
Magnesium mass	M_Mg	1.30 (−2)	M_⊙
Silicon mass	M_Si	1.57 (−1)	M_⊙
Sulfur mass	M_S	4.36 (−2)	M_⊙
Argon mass	M_Ar	6.30 (−3)	M_⊙
Iron mass	M_Fe	8.49 (−2)	M_⊙

Total mass	M_tot	3.70	M_⊙

Density	ρ	7.86 (−7)	g cm⁻³
Temperature	T	5.76 (5)	K
Radial velocity	v_r	1.43 (8)	km s⁻¹

Note.

^aThe numbers in the parentheses denote the powers of 10.

Download table as: ASCII Typeset image

Therefore, when the density drops to a certain critical value (ρ_break) for the first time, the adiabatic index is changed manually (in the calculations, it is automatically switched according to the density) from 4/3 to 5/3 if the timing is before the start of the molecule formation calculation. If the timing of the change is after the start of the molecule formation calculation, the adiabatic index is changed from 4/3 to an effective adiabatic index derived from the averaging of the ones for diatomic species and the others with the mass fractions as a weight. Such sudden change in the adiabatic index results in a jump in the gas temperature before and after the change with the same internal energy. In reality, after the system becomes optically thin, the energy of radiation escapes from the gas. In this study, some fraction of the internal energy of the gas is subtracted at the transition point in the adiabatic index to keep the same gas temperature before and after it. As seen in Figure 2, after the break, gas temperatures seem to remain at temperatures ≳10³ K; the trend may be attributed to not a physical reason but a computational one in the 1D radiation hydrodynamical calculations (A. Kozyreva et al. 2024, in preparation). Therefore, we adopt the single break in the specific internal energy (temperature) as a fiducial thermal evolution. Based on the power-law evolution in Equations (35) and (36), additional heating and cooling effects are taken into account as follows.

As mentioned in Section 2.1.2, the decay of radioactive ⁵⁶Ni and ⁵⁶Co could deposit some energies into the gas depending on the optical depth τ_γ. In reality, the deposition processes, i.e., how the emerged gamma rays produce Compton electrons and how the energies of Compton electrons are degraded by the interaction with gas particles, are complicated; the deposited energies could be consumed by the excitation of electrons, ionization, and heating of gas (Liu & Victor 1994). Due to the lack of accurate treatments for the time evolution of the electron energy distribution, the ionization, and the radiative transfer in this paper, we adopt a simple approach to the heating of gas via the decay of ⁵⁶Ni, i.e., some fraction of the deposited energy expressed in Equation (11) is locally used for the heating of the gas. The potential impact of nonlocal energy deposition is discussed in Appendix B as mentioned in Section 2.1.2. To count the effective efficiency of the heating of the gas, a constant factor, f_h, which denotes the fraction of the deposited energy to be used for the heating of the gas, is introduced as a parameter.

It is noted that in this paper, radiative cooling by line emission from ionized atomic species (Kozma & Fransson 1998) is not taken into account. In particular, in the metal-rich ejecta core, such radiative cooling could be faster than adiabatic cooling (see, e.g., Figure 5 in Kozma & Fransson 1998). Actually, as can be seen in Figure 2, some of the inner particles have steeper temperature slopes than that of the case of adiabatic cooling for the ideal gas (n = − 2; dashed line). Such inner particles tend to contain ⁵⁶Ni inside. Therefore, in such inner particles, the heating of gas due to the decay of ⁵⁶Ni may compete with the radiative cooling by ionized atomic species; neglecting the radiative cooling may partly be compensated by adopting lower f_h values. We would regard that the parameter, f_h, effectively covers the uncertainties due to neglecting the radiative cooling by atomic species, although the effects may not be represented by such a constant factor. More realistic treatments of gas heating and cooling are beyond the scope of this study.

Given the effective energy deposition rates, the change in the specific internal energy of gas per unit time is expressed as below

$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}={f}_{{\rm{h}}}\,\epsilon ={f}_{{\rm{h}}}\,({\epsilon }_{\mathrm{Ni}}+{\epsilon }_{\mathrm{Co}})\,{D}_{\gamma }.\end{eqnarray} \tag{ 38 }$

With the temperature evolution, the abundance of CO is obtained through the chemical reaction calculations described in Section 2.1. Once CO forms in the ejecta, it could be an important coolant through the emission via rovibrational transitions. As the emissivity is given in Equation (34), the cooling rate of gas can be calculated as follows

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{de}}{{dt}} & = & -\displaystyle \frac{h}{\rho }\displaystyle \sum _{{E}_{\kappa }\gt {E}_{\lambda }}\,{\nu }_{\kappa \lambda }\,{n}_{\kappa }\,{A}_{\kappa \lambda }{\beta }_{\kappa \lambda }\\ & & \times \left(\displaystyle \frac{{S}_{\kappa \lambda }-{B}_{{\rm{d}}}({\nu }_{\kappa \lambda },{T}_{\mathrm{ph}})}{{S}_{\kappa \lambda }}\right)\\ & & \ \ \ \mathrm{erg}\,{{\rm{g}}}^{-1}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 39 }$

where S_{κ
λ} is the so-called source function, which can also be seen in the first term of Equation (28), and is expressed below

$\begin{eqnarray}&&{S}_{\kappa \lambda }=\displaystyle \frac{2h{\nu }_{\kappa \lambda }^{3}}{{c}^{2}}{\left(\displaystyle \frac{{X}_{\lambda }/{g}_{\lambda }}{{X}_{\kappa }/{g}_{\kappa }}-1\right)}^{-1}.\end{eqnarray} \tag{ 40 }$

Here, if the excitations are mainly maintained by the radiation from the background photons, the radiated energies should not be subtracted from the internal energy of the gas. Then, a controlling factor, i.e., inside the parentheses in Equation (39), is introduced as in de Jong et al. (1975). The source function, S_{κ
λ}, and the diluted Planck function, B_d(ν_{κ
λ}, T_ph), represent the contributions from the local and background photons, respectively. If S_{κ
λ} < B_d(ν_{κ
λ}, T_ph) (the contribution of the background photons dominates the local one), the controlling factor is set to zero, i.e., no energy subtraction from the internal energy of the gas. In the case of S_{κ
λ} ≥ B_d(ν_{κ
λ}, T_ph), the energy subtraction from the internal energy of the gas is controlled by the factor of less than unity depending on the magnitudes of the source function and the diluted Planck function. In the case of B_d(ν_{κ
λ}, T_ph) = 0, i.e., no contribution from the background photons, the controlling factor becomes unity.

If the specific internal energy is updated by the heating and/or cooling in an operator-splitting manner as in Equations (38) and (39), through the equation of state in Equation (37), the temperature of the gas is also updated. Additionally, if the extra heating (cooling) is dominant, the local gas could be inflated (deflated). To mimic such an effect, by assuming a pressure balance between the local gas and the surroundings, the density of the gas with such extra heating and/or cooling, ρ, is obtained by equating the two cases of the right-hand side of Equation (37), i.e., with and without net extra heating/cooling as follows

$\begin{eqnarray}&&\rho ={\rho }_{\mathrm{ad}}({e}_{\mathrm{ad}}/e),\end{eqnarray} \tag{ 41 }$

where ρ_ad and e_ad denote the density and the specific internal energy, respectively, obtained from the evolution with no extra heating by the decay of ⁵⁶Ni and cooling by CO rovibrational transitions, i.e., only with the adiabatic cooling described in Equations (35) and (36) throughout the calculation. e is the specific internal energy obtained from the evolution with net extra heating and/or cooling. Once the extra cooling becomes dominant, the local gas could be deflated, and the density increases. Such an enhancement in density may cause additional emissions, which may lead to thermal instability. In this manner, such a potential impact is effectively taken into account in the calculations. The methodology described in Equation (41) is, however, postprocessing; the feedback of the inflation and deflation to the hydrodynamics is not taken into account. The feedback to the hydrodynamics may result in expanded bubble-like structures around ⁵⁶Ni-rich regions (as seen in Gabler et al. 2021) and/or denser clumps around CO-rich regions.

Because of the adiabatic cooling and/or cooling by CO rovibrational transitions, with the methodology above, the gas temperature drops to ∼10² K at some point. In this study, practically, we set a minimum gas temperature to be 10² K assuming that the ejecta keeps the minimum temperature thanks to some heating, e.g., by the radioactive decay of ⁴⁴Ti.¹⁸ It is noted that the ALMA observations of molecular lines of CO, SiO, and HCO⁺ in the ejecta of SN 1987A have shown the temperature of the ejecta gas ranges over 13–132 K (Kamenetzky et al. 2013; Matsuura et al. 2017; Cigan et al. 2019). Because of this minimum temperature, the minimum specific internal energy is also implicitly set through the EoS. Practically, by turning off the adiabatic cooling described in Equation (36) after both e and e_ad become the minimum value, e can be equal to e_ad at some point without effective heating due to the decay of ⁵⁶Ni, i.e., the density evolution becomes again consistent with the original power-law evolution after that. In this way, we avoid introducing the effect described in Equation (41) at such a low temperature environment.

As a summary, the chemical reaction calculation, the calculation of CO (ro)vibrational levels for CO line cooling, and the updating of the specific internal energy and gas temperature are recurrently performed to obtain the time evolution. The time step for the iteration is empirically determined not to cause a change of a large fraction of the internal energy per step. Additionally, for the chemical reaction calculation and the rate equation calculation for the CO (ro)-vibrational levels, subtime steps are independently introduced.

As will be discussed in Section 4.1 and Appendix D, the calculations with the method described above tend to result in higher CO line emissions than that corresponding to the peak fluxes of the observed CO light curves (Bouchet & Danziger 1993; Meikle et al. 1993) at an early phase (≲200 days). Considering the lack of realistic treatments for the radiative transfer, ionization, and energy deposition through the decay of ⁵⁶Ni in this study, and to avoid a significant cooling in the ejecta gas, a time-dependent reduction factor (function) f_red for the escape probability β_{κ
λ} found in Equations (28), (34), and (39) is introduced. The function is arbitrarily given by

$\begin{eqnarray}&&{f}_{\mathrm{red}}=\exp \left[-({t}_{{\rm{c}}}-t)/{t}_{{\rm{s}}}\right]\ \ \ t\lt {t}_{{\rm{c}}},\end{eqnarray} \tag{ 42 }$

where t_c is an arbitrary critical time up to when the factor is effective. t_s is a timescale to control the magnitude of the factor. Then, the escape probability β_{κ
λ} is practically multiplied by the reduction factor f_red to reduce the escape probability (to increase the line optical depth effectively). As seen in Equation (29), the β_{κ
λ} is dependent on the line optical depth τ_{κ
λ} in which the age t is introduced as a rough approximation for the reciprocal of the velocity gradient of the system (see, e.g., de Jong et al. 1975). Then, we consider that there is some uncertainty at least in the approximation. The reduction factor is not dependent on the transition levels, but the resultant escape probability, f_red β_{κ
λ}, is varied depending on the transitions. The values of t_c and t_s are empirically determined to avoid a significant cooling incompatible with the observed light curves. The impact of the reduction factor f_red through different values of the parameter t_s on the CO rovibrational emissions is discussed in Section 4.1 and Appendix D.

The model parameters related to the calculations with the one-zone approximation are listed in Table 7.

Table 7. Model Parameters for One-zone and 1D Calculations

Parameter	Explanation	Value/Range
f_h	Efficiency of gas heating by the decay of ⁵⁶Ni	10⁻⁴–10⁻²
f_d	Efficiency of the destruction/ionization via the decay of ⁵⁶Ni	10⁻²–1.0
t_c	Critical time in the reduction factor f_red in Equation (42) for the escape probability β_{κ λ}	1000 days (fixed)
t_s	Timescale in the reduction factor f_red in Equation (42) for the escape probability β_{κ λ}	200–∞ days
X_e,i	Initial ionization fraction	10⁻¹ (fixed)
X_e,f	Ionization fraction at 2000 days	10⁻³ (fixed)
ρ_break ^a	Critical density for the transition of the adiabatic index γ_ad	10⁻⁹ g cm⁻³ (fixed)

Note.

^aThis parameter is only for the one-zone calculations.

Download table as: ASCII Typeset image

2.4. Calculations with 1D Radial Profiles

To figure out the chemical evolution in the ejecta in a more realistic way, next, 1D (radial) profiles derived from 3D hydrodynamical models for SN 1987A (Ono et al. 2020) are utilized as the initial conditions. By implementing angle-averaging with the density as a weight for quantities other than the density on the results of the last snapshots at ∼1 day after the explosion, 1D profiles are obtained (as partly mentioned in Section 2.3). Among the 3D hydrodynamical models in Ono et al. (2020), the models b18.3-high and n16.3-high are selected. The former model is the one with an aspherical (bipolar-like) explosion with the binary merger progenitor model (Urushibata et al. 2018), which explains the observational features of the progenitor star of SN 1987A, Sk −69° 202; the explosion model also best explains the observed [Fe ii] line profiles (Haas et al. 1990) among the models in Ono et al. (2020), and its further evolution (Orlando et al. 2020) also well reproduces the observed X-ray morphology and light curves (e.g., Frank et al. 2016) and the spatial distribution of CO and SiO (Abellán et al. 2017) in SN 1987A. The latter model is the model with the explosion with the same asymmetry and explosion energy as the model b18.3-high, but the progenitor model (Nomoto & Hashimoto 1988; Shigeyama & Nomoto 1990) is based on a single-star evolution, and it was modeled somehow artificially to fit the observations of Sk −69° 202. As mentioned in Section 1, the density structure of the progenitor star could play an important role in the matter mixing (Wongwathanarat et al. 2015). Actually, as seen in Ono et al. (2020), because of the differences between the two progenitor models, the mixing of the radioactive ⁵⁶Ni into high-velocity outer laters changes and affects the velocity distribution of ⁵⁶Ni corresponding to the observed [Fe ii] line profiles. Therefore, we investigate the dependence of chemical evolution on the two progenitor models.

To further see the impact of the matter mixing, we also investigate the molecule formation for the two cases where the SN explosion itself is spherical, but the evolution is 3D, in which mixing via hydrodynamical instabilities such as RT instability could be involved, and the whole evolution is purely spherical. To obtain the 1D profiles as the initial conditions for the two cases, additional hydrodynamical calculations are performed as follows. For the former case, two 3D hydrodynamical simulations are newly performed up to ∼1 day after the explosion as the counterparts of the models b18.3-high and n16.3-high. With the same method in Ono et al. (2020), i.e., the explosions are artificially initiated by injecting thermal and kinetic energies asymmetrically, 3D hydrodynamical simulations are performed with the same progenitor models and injection energies, but the energy injections are spherical. Then, angle-averaged 1D profiles are obtained based on the last snapshots of the 3D hydrodynamical simulation results. In the latter case, two 1D hydrodynamical simulations are performed with the same method in Ono et al. (2020) and the same progenitor models and injection energies but with a 1D version of the same open-source hydrodynamical code, FLASH (Fryxell et al. 2000), with the spherical coordinate (for the details of the 1D code, see later in this section).

In Figure 3, the angle-averaged 1D profiles of the mass fractions of atoms¹⁹ and hydrodynamical quantities, the density, temperature, and radial velocity, based on the models b18.3-high and n16.3-high (Ono et al. 2020), are shown. Overall, elements heavier than carbon (left panels) are smoothly extended to the radius of ∼3 × 10¹³ cm (∼3.5 × 10¹³ cm) for the model b18.3-high (n16.3-high). Between the two models, there are significant differences in the composition. For example, the mass fraction of oxygen (O) in the model n16.3-high is higher than that of b18.3-high inside the radius of ∼10¹³ cm. On the other hand, the mass fraction of carbon (C) in the model b18.3-high is higher than that of n16.3-high at the inner region. Additionally, parts of hydrogen (H) and helium (He) originally in the outer layers are mixed into the inner regions.

As mentioned above, matter mixing sensitively depends on the density structure of the progenitor star (Wongwathanarat et al. 2015; Ono et al. 2020); the progenitor model of b18.3-high (BSG through a binary merger; Urushibata et al. 2018) results in an efficient matter mixing thanks to the compact helium core. Then, in the model b18.3-high, the mass fractions of hydrogen and helium are higher than that of n16.3-high at the inner region. For the same reason, the mass fraction of iron (Fe; mostly coming from ⁵⁶Ni at inner regions) is affected to be mixed into outer layers. The mass fraction of iron in the model b18.3-high has a plateau between the radii of (1.6–2.6) × 10¹³ cm. It is noted that, in the model b18.3-high, the abundance of nitrogen (N) in the hydrogen envelope is apparently higher than that of the model n16.3-high. This is because the progenitor model of the b18.3-high is the binary merger model; during the merger process of the companion, additional hydrogen burning including CNO cycle is triggered. Then, CNO cycle-processed materials (nitrogen, more specifically ¹⁴N, is one of the main products of the CNO cycle) are mixed into the envelope. As for the hydrodynamical values, i.e., the density, temperature, and radial velocity, there are no significant differences between the two models.

Figure 4 shows the 1D profiles of the spherical explosion case. The profiles of the composition are apparently smoother (with lesser radial bumps/irregularities) than that of the previous case. However, the mass fractions of elements heavier than carbon other than iron are roughly peaked around the radius of ∼10¹³ cm, and the one of iron is concentrated in the innermost region for both models. Since iron (practically ⁵⁶Ni) is limited to the radius of ∼10¹³ cm, the destruction and ionization of molecules at outer regions due to the decay of ⁵⁶Ni may be restricted compared with the previous case. It is noted that the amount of ⁵⁶Ni obtained by the 3D hydrodynamical simulations in the spherical explosion case is higher than that of the original 3D models, i.e., b18.3-high and n16.3-high. The amount of ⁵⁶Ni in the model b18.3-high (n16.3-high) is 8.64 × 10⁻² M_⊙ (9.67 × 10⁻² M_⊙). On the other hand, the amount of ⁵⁶Ni of the spherical explosion case corresponding to the model b18.3-high (n16.3-high) is 1.92 × 10⁻¹ M_⊙ (1.89 × 10⁻¹ M_⊙). Therefore, the amounts of ⁵⁶Ni in the spherical explosion case are higher roughly by a factor of 2 than those in the 3D models. The high amounts of ⁵⁶Ni in the spherical case are attributed to the fact that, in the spherical explosion models, larger regions are heated to enough high temperature (∼5 × 10⁹ K; e.g., Thielemann et al. 1990) for the synthesis of ⁵⁶Ni during the explosive nucleosynthesis. In the 3D (bipolar-like explosion) models, only limited regions (steradians) have enough high temperatures. On the other hand, in the spherical explosion models, all directions are heated evenly to have enough high temperatures. The differences in oxygen, hydrogen, and helium between the two models due to the different adopted progenitor models seen in the previous case, i.e., angle-averaged 1D profiles (Figure 3), are roughly kept in the spherical explosion case.

In Figure 5, the 1D profiles of the purely spherical case are shown. The composition profiles are apparently more stratified compared with the previous two cases; from the outside, hydrogen-, helium-, carbon plus oxygen-, and iron-rich layers are distinct. In particular, ⁵⁶Ni is confined to the innermost regions. The positions of the distinct layers are different between the two models depending on the size of the cores of the progenitor star and its density structure, which affects the acceleration and deceleration of the SN shock. Since the size and mass of the helium core of the progenitor model for the model n16.3-high are both larger than those for the model b18.3-high, the inner cores (layers) are pushed to the helium layers due to the stronger deceleration than that for the model b18.3-high. The peak of the density profile for the model n16.3-high (bottom right panel) is positioned outward compared with that for the model b18.3-high (top right panel). The amount of ⁵⁶Ni of the purely spherical case corresponding to the model b18.3-high (n16.3-high) is 1.82 × 10⁻¹ M_⊙ (1.95 × 10⁻¹ M_⊙). The amounts of ⁵⁶Ni in the purely spherical case are also higher than the 3D models because of the same reason mentioned above. The purely spherical case is useful to see the impact of matter mixing as the representative case of no mixing as mentioned.

It is noted that the amounts of ⁵⁶Ni presented above may have some uncertainties. The explosions in the models in this study, i.e., the 3D models (Ono et al. 2020) and the spherical cases, are initiated by a kind of thermal bomb (thermal and kinetic energies are artificially injected around the surface of the original iron core). Conventionally, in the spherical CCSN explosion models, the so-called "mass cut" (e.g., Woosley & Weaver 1995), i.e., the mass coordinate from where the matter inside should fall back to the compact object, has been introduced to compensate for the limitations in the theories and to obtain an appropriate ⁵⁶Ni mass consistent with the light-curve observations (∼0.07 M_⊙; e.g., Shigeyama et al. 1988; Woosley 1988). In the 3D explosion models (Ono et al. 2020), the mass cut is not introduced to adjust the ⁵⁶Ni mass since the obtained ⁵⁶Ni mass is not far from the value obtained from the observations, and how to determine the boundary between the ejecta and compact object is not simple for the case of globally asymmetric explosions (see Nagataki 2000). In the spherical (spherical explosion and purely spherical) models in this paper, we also do not consider the mass cut to have consistent hydrodynamical conditions, e.g., the explosion energy, other than the explosion asymmetry for simple comparison. Considering the adopted method in the explosion models in this study, more realistic explosion models may result in different ⁵⁶Ni masses. Additionally, the nucleosynthesis calculations may also introduce uncertainties in the amounts of ⁵⁶Ni as shown in our previous study (Mao et al. 2015).²⁰ We expect that the uncertainties in the amounts of ⁵⁶Ni due to the points above may be at most a factor of 2 as seen in the results in this paper unless the explosion energy and the progenitor model are significantly different from those in this study. As shown in the later sections, the distribution of ⁵⁶Ni (how ⁵⁶Ni is extensively distributed in outer layers) is more important for the formation of molecules. Uncertainty within a factor of 2 in the amount of ⁵⁶Ni is acceptable for the purpose of this study (to see the qualitative impact of the matter mixing on the formation of molecules).

In order to obtain the hydrodynamical evolution from ∼1 day to ∼10⁴ days, 1D hydrodynamical simulations are further performed with the angle-averaged 1D profiles as the initial conditions. For the calculations, the 1D version of the hydrodynamical code FLASH (Fryxell et al. 2000) mentioned above is used again. In the code, the directionally split Eulerian hydrodynamical unit with the piecewise parabolic method (Colella & Woodward 1984) is utilized. As for the Riemann solver, a hybrid solver with an HLLE solver inside the shocks is adopted. The coordinate system is the spherical coordinates. In the code, the so-called adaptive mesh refinement is implemented, but a uniform grid with 1024 meshes is used in the calculations. As in our previous studies on the matter mixing (Ono et al. 2013, 2020; Mao et al. 2015), as the SN shock approaches the outer computational boundary, the computational domains are expanded with a scaling factor of 1.2 to resolve the structure inside the SN shock and follow the long-term evolution. In the pre-shocked regions outside the progenitor star, a spherically symmetric BSG wind with a constant mass-loss rate of $\dot{{M}_{{\rm{w}}}}={10}^{-7}{M}_{\odot }$ yr⁻¹, and wind velocity of v_w = 500 km s⁻¹ (Morris & Podsiadlowski 2007) is assumed as the initial condition as in Ono et al. (2020), i.e., the density profile follows the power law of ρ(r) ∝ r⁻², where r is the radius from the center of the explosion. The following points are different from the 1D hydrodynamical simulation mentioned in the previous paragraph. During the simulations, the unit for the nucleosynthesis is turned off, since, at epochs ∼1 day after the explosion, explosive nucleosynthesis no longer occurs. In the equation of state, the contribution from the radiation is omitted, i.e., the EoS is the ideal gas EoS with the adiabatic index of γ_ad = 5/3. The energy deposition from the decay of ⁵⁶Ni is also turned off. Neglecting the energy deposition from the decay does not affect the overall dynamics, but the effect of the decay of ⁵⁶Ni on the thermal evolution of the ejecta and the chemistry is taken into account similar to the one-zone calculations as mentioned later.

The chemical evolution in the ejecta is calculated in post-processing. For the sake of this, first, based on the 1D hydrodynamical simulation results, the time evolution of tracer particles (imaginary Lagrange particles) is obtained with the procedures described below. In total, one hundred tracer particles are radially distributed in the initial computational domain of the hydrodynamical simulation to cover the inner metal-rich ejecta, which includes seed atoms. Initially, for each particle, the mass and composition are assigned depending on its initial position (the values are obtained with the same interpolation method as described later in this paragraph); the mass is assigned to be a constant value so that the particle covers the mass corresponding to a shell in which the particle initially is. As for the composition, the initial one is respected, i.e., without a change via chemical reactions, each particle's composition is assumed to be constant. Then, the tracer particles at t = tⁿ (here, integer n denotes the time step) are moved along the velocity field obtained from the hydrodynamical simulation to estimate the positions at tⁿ⁺¹ = tⁿ + Δt. To obtain the time evolution accurately, as a time integration method, the so-called predictor-corrector method is used, i.e., generally, the position of a particle at t = tⁿ⁺¹, x ⁿ⁺¹, is obtained with the velocity field v ( x ) as follows

$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{x}}^{\prime} & = & {{\boldsymbol{x}}}^{n}+{{\boldsymbol{v}}}^{n}({{\boldsymbol{x}}}^{n})\,{\rm{\Delta }}t,\\ {{\boldsymbol{x}}}^{n+1} & = & {{\boldsymbol{x}}}^{n}+\displaystyle \frac{{{\boldsymbol{v}}}^{n+1}({\boldsymbol{x}}^{\prime} )+{{\boldsymbol{v}}}^{n}({{\boldsymbol{x}}}^{n})}{2}\,{\rm{\Delta }}t,\end{array}\end{eqnarray} \tag{ 43 }$

where ${\boldsymbol{x}}^{\prime}$ is the temporal value (predictor) to obtain the correct value x ⁿ⁺¹ (corrector). Shortly, we plan to apply our method to 3D Eulerian hydrodynamical simulation results. In such a 3D case, practically, it is difficult to save the data for all the time snapshots. Therefore, we check among several time integration methods, e.g., the fourth order Runge–Kutta method and an implicit midpoint method, which method best reproduces the full-time step results with limited time slices (outputs every 10 or 100 time steps). Then, we conclude that the predictor-corrector method is the most appropriate for this purpose. Even though there is no such a time snapshot problem in the calculations limited in 1D in this paper, the predictor-corrector method is adopted. The physical quantities, e.g., the density, temperature, and velocity, for each particle are derived by time and spacial interpolations with the data for nearby time slices and grid points, respectively. A monotonic cubic interpolation method (Steffen 1990) is utilized for the interpolations of physical quantities.

The thermal and chemical evolutions are calculated for each tracer particle. The basic methodology is the same as the case of the one-zone approximation as described in Section 2.3. From here, the points different from the one-zone case are mainly delineated. Although the time evolution of gas temperatures (internal energies) can be obtained from the interpolations with the 1D hydrodynamical simulation results, the thermal evolution of each particle is obtained by the basic power laws plus contributions from the additional heating and/or cooling as in the one-zone case. This is because, in the 1D hydrodynamical simulations, the contribution from the radiation (in the ideal EoS, the constant adiabatic index of γ_ad = 5/3 is used), the possible energy depositions via the decay of ⁵⁶Ni, and the unknown CO line cooling are neglected. As for the time evolution of the density and radial velocity, the values obtained by the interpolations, i.e., the hydrodynamical simulation results, are respected except for the effect described in Equation (41) in the density. In the case of the 1D calculations, ρ_ad in Equation (41) is replaced with the density obtained by the interpolation. As in the one-zone case, it is assumed that the gas temperature (specific internal energy) has a break at some time point when the gas becomes optically thin. At such point, the adiabatic index γ_ad is manually changed. However, the transition point is differently determined from the one-zone case as described below. To determine the transition point, as a reference, the optical depth to Thomson scattering, τ_Th, is derived for each particle as below

$\begin{eqnarray}&&{\tau }_{\mathrm{Th}}(r)={\int }_{r}^{\infty }{\sigma }_{\mathrm{Th}}\,{n}_{{\rm{e}}}\,{dr},\end{eqnarray} \tag{ 44 }$

where σ_Th is the Thomson scattering cross section. Here, the number density of electrons n_e reflects the gas density and the ionization and recombination through the chemical reaction calculation. In deriving n_e, the gas density obtained from Equation (41) is used. As described in Section 2.1.3, before starting the chemical reaction calculation, the degree of the ionization is controlled by the ionization fraction X_e. After starting the chemical reaction calculation, n_e is derived from the resultant abundance of (thermal) electrons, although most of the particles may have a break before the chemical reaction calculation starts. For each particle, once the τ_Th becomes less than 2/3, regarding that the gas becomes optically thin, the adiabatic index γ_ad is automatically switched from 4/3 to 5/3 (before the start of the molecule formation calculation) or the effective adiabatic index derived from the averaging of ones for the diatomic and monoatomic species (after that). Other parts of the calculations are the same as the one-zone case; the chemical reaction calculations, the calculations of CO (ro)-vibrational levels for CO line cooling, and the updating of the internal energies and gas temperatures are recurrently carried out for each tracer particle. The model parameters related to the calculations with the 1D profiles are listed in Table 7.

3. Strategy and Models

In this section, the strategy to investigate the impact of the matter mixing on molecule formation in the ejecta with the method delineated in the previous section and the related models are described.

First, with one-zone calculations, the impacts of the parameters and overall trends are briefly explained. By comparing the results of the time evolutions of the masses of CO and SiO estimated in previous studies (e.g., Liu & Dalgarno 1995) and observed light curves of CO vibrational bands in SN 1987A (Meikle et al. 1989, 1993; Bouchet & Danziger 1993; Wooden et al. 1993), better parameter values for the one-zone approximation are derived. Additionally, important chemical reactions in particular for CO and SiO molecules are extracted according to the one-zone calculation results for the parameter values derived above. The results are presented in Section 4.1.

Second, with the angle-averaged 1D profiles for the model b18.3-high (Ono et al. 2020), similarly to the one-zone calculations, an acceptable set (acceptable sets) of parameter values for 1D calculations are derived by comparing the results with the previous studies and observations mentioned above. Then, the parameter values are fixed to be the acceptable values to investigate the impact of the matter mixing on the molecule formation for the latter discussion. The impact of the matter mixing in the ejecta is effectively discussed by comparing the model results for angle-averaged 1D profiles (full mixing) and profiles corresponding to the spherical explosion (but the further evolution is in 3D) case (partial mixing), and purely spherical case (no mixing). The results are shown in Section 4.2.

A caveat on the discussion with the initial profiles above is that the angle-averaged profiles may overestimate the mixing and the molecule formation results could be different between the molecule formation calculation with a single angle-averaged profile and ones for all radial profiles in 3D models. Therefore, similar molecule formation calculations for a few radial profiles along several specified directions (angle-specified profiles) are performed to obtain an insight into the impact of more realistic matter mixing based on full 3D hydrodynamical models. To do this, a 1D hydrodynamical simulation is performed with each angle-specified profile as the initial conditions to obtain the evolution from ∼1 day to ∼10⁴ days after the explosion with the same method described in Section 2.4. Therefore, implicitly, it is assumed that the profile is spherically distributed. Here, selected representative directions are described. As already mentioned in Section 2.3, among the models investigated in Ono et al. (2020), the asymmetric bipolar-like explosion model, i.e., the model b18.3-high, with the BSG progenitor model through a binary merger (b18.3; Urushibata et al. 2018), best reproduces the observed [Fe ii] line profiles (Haas et al. 1990) among the models. Additionally, with the same explosion energy and the explosion asymmetry, a different matter mixing was demonstrated based on the explosion model (n16.3-high) with another BSG progenitor model of a single-star evolution (n16.3; Nomoto & Hashimoto 1988; Shigeyama & Nomoto 1990). Then, radial profiles are extracted from the two explosion models, b18.3-high and n16.3-high. As for the directions specified, two directions along the bipolar-like explosion axis, i.e., +Z- and −Z-directions, are selected, where the strongest explosion is directed to the +Z-direction, and the explosion is weaker in the −Z-direction (see Figure 23 in Ono et al. 2020). As another representative direction, the +Y-direction on the equatorial plane, against which the bipolar explosion is asymmetric, is selected.

Figure 6, for a reference, shows the distributions of representative elements, ⁵⁶Ni, ²⁸Si, ¹⁶O, and ¹²C, at the end of the simulation in the 3D models b18.3-high and n16.3-high (Ono et al. 2020), which are shown with the arrows denoting the +Z-, −Z-, and +Y-directions. The distributions reflect the asymmetric bipolar-like explosions, and those are different between the two models with the different progenitor models. It is noted that, compared with the distributions just before the shock breakout (see Figure 23 in Ono et al. 2020), smaller-scale structures are lost after the shock breakout. This is because the forward shock (FS) is accelerated due to the steep pressure gradient across the stellar surface after the shock breakout, and the inner ejecta becomes left far behind the FS. Then, to keep the resolution inside the FS with the limited number of meshes, the resolution of the inner ejecta relatively becomes low after the shock breakout. Even though the smaller-scale structures are lost, the asymmetry caused by the bipolar-like explosions and the overall distributions of the seed atoms remain; such a numerical artifact should not change the main findings presented in later sections.

In Figures 7 and 8, the angle-specified radial profiles of the mass fractions of atoms and hydrodynamical quantities, the density, temperature, and radial velocity, at the end of the simulation (corresponding to initial conditions in this paper) are plotted for the models b18.3-high and n16.3-high, respectively. The positions (radii) of the FS depend on the directions (see the left panels). In the case of the model b18.3-high (Figure 7), the radii of the FS along the three directions are ∼2.2 × 10¹⁴ cm, ∼2.0 × 10¹⁴ cm, and ∼1.5 × 10¹⁴ cm in descending order corresponding to the directions, +Z, +Y, and −Z, respectively, depending on the initial explosion asymmetry. The fact that the FS radius in the +Y-direction is larger than that in the −Z-direction sounds a bit strange since, at the shock-revival phase, the explosion (energy injection) is weakest in the +Y-direction. It is difficult to specify the reason clearly, but it may be interpreted as follows. Because of the initial weakest explosion in the +Y-direction (due to the artificially imposed bipolar-like explosions), the density of the expanding ejecta in the +Y-direction tends to be higher than those along the bipolar-like explosion axis during the propagation of the SN shock inside the star (see the left panels in Figure 19 in Ono et al. 2020)²¹ ; the dense material in the +Y-direction maintains its outward velocity for a longer time than that in other directions, thanks to its high inertia. Consequently, the FS in the +Y-direction is influenced and deformed by the dense material. Additionally, Rayleigh–Taylor (RT) instability develops in the +Y-direction after the shock wave enters the hydrogen envelope, due to the steep density gradient between the dense material and the surroundings.²² This instability may partly play a role in exaggerating the expansion of the dense material and the deformation of the FS. Actually, extended finger-like structures can be recognized along the equatorial plane (see the left panel in Figure 23 in Ono et al. 2020). On the other hand, in the model n16.3-high, the positions of the FS are more identical to each other. The composition profiles are radially fluctuated compared with the 1D profiles described above (Figures 3–5). The positions of the outer boundaries of the metal-rich cores do not necessarily reflect the FS radii; for example, in the model b18.3-high (Figure 7), the outer boundaries along −Z- and +Y-directions are similar to each other in contrast to the FS radii; in the model n16.3-high, the outer boundary along the +Z-direction is more extended than those along the −Z- and +Y-directions approximately by a factor of 2. It is noted that ⁵⁶Ni in the +Y-direction is very concentrated in the innermost regions for both b18.3-high and n16.3-high models. This is because ⁵⁶Ni is not abundantly synthesized during the explosive nucleosynthesis due to the less energetic explosion in the equatorial plane. Again, although the amounts of ⁵⁶Ni could be uncertain at most by a factor of 2, the qualitative differences among the directions (how ⁵⁶Ni is extensively distributed) due to the bipolar-like explosions may not be affected by the bipolar-like explosions obtained by more realistic models, e.g., with a better fallback model and/or a larger nuclear reaction network, with the similar explosion energies and the progenitor models.

**Figure 7.** Initial radial profiles along the specific directions, +Z (top), −Z (middle), and +Y (bottom), for the model b18.3-high at the age of 19.0 hr after the explosion. Mass fractions of species as a function of the radius (left) and the density, temperature, and radial velocity as a function of radius (right).
Download figure:
Standard image High-resolution image

**Figure 8.** Same as Figure 7 but for the model n16.3-high and at the age of 22.3 hr after the explosion.
Download figure:
Standard image High-resolution image

By the way, in reality, depending on the explosion mechanism, e.g., the delayed neutrino-driven explosion aided by convection/SASI or the magnetorotationally driven explosion as mentioned in Section 1, the distribution of ⁵⁶Ni may be different from the idealized bipolar-like explosions in this paper, although some neutrino-driven explosion models seem to result in an asymmetric dipole-like morphology (e.g., Vartanyan et al. 2019). As shown in the later sections, the distribution of ⁵⁶Ni could affect the formation of molecules through several processes. The comparisons of theoretical models like those presented in this study with observations may provide insights to pin down the explosion mechanism in the future.

Here, in order to easily distinguish calculation models, we name the models for 1D calculations as follows. Each model is distinguished by specifying the adopted progenitor model (b18.3 or n16.3) and the type of the initial 1D profiles (angle-averaged/spherical explosion/purely spherical/angle-specified) after fixing the values of the parameters, f_h, f_d, and t_s. Then, a model name is expressed as "Progenitor model"-"Type of 1D profile," for example, b18.3-mean. Here, "b18.3" and "mean" denote the binary merger progenitor model (b18.3) and the angle-averaged profiles, respectively. Similarly, for the first key, "n16.3" denotes the single-star progenitor model (n16.3). For the second key, "sphel," "sphel-pure," "zp," "zn," and "yp" denote the 1D profiles for the spherical explosion case, the purely spherical case, along the +Z-, −Z-, and +Y-axes, respectively.

4. Results

4.1. One-zone Calculation Results

To find a set (sets) of values of the parameters, which provides (provide) reasonable results compared with previous studies and observed CO light curves, the impact of those parameters is investigated. Showing all the calculated results is rather lengthy, and the purpose of the paper is not the investigation of the impact of those parameters but the impact of matter mixing. Instead, only several calculation results are presented to describe the overall evolution and features. Among several parameters listed in Table 7, critical parameters, which affect the molecule formation results significantly, are f_h, f_d, and t_s. Other parameters, i.e., t_c, X_e,i, X_e,f, and ρ_break, are practically fixed to be a reasonable value. For simplicity, the value of t_c is set to be 1000 days up to when the calculation of CO (ro)-vibrational transitions is performed. The values of X_e,i and X_e,f, are set to be 0.1 and 10⁻³, respectively by reference to the values shown in Figure 8 in Kozma & Fransson (1998), in which the thermal and ionization evolutions in the ejecta of SN 1987A were modeled. We confirmed that changing the values of X_e,i and X_e,f affects negligibly the molecule formation results. The value of ρ_break changes the timing of the break in the gas temperature, and a higher (lower) ρ_break value causes an early (late) break and an early (late) molecule formation, since the gas temperature drops to ∼10⁴ K, at which molecules start to form, at an early (late) phase. The value of ρ_break is set to be 10⁻⁹ g cm⁻³, which is consistent with the treatment in the EoS of the hydrodynamical simulations (Ono et al. 2020). Then, by changing the values of the parameters, f_h, f_d, and t_s, in total, 100 cases are calculated: the cases investigated are the combinations of f_h = 10⁻⁴, 5 × 10⁻⁴, 10⁻³, 5 × 10⁻³, 10⁻², f_d = 10⁻², 5 × 10⁻², 10⁻¹, 5 × 10⁻¹, 1.0, t_s = 200, 300, 500, and ∞ (practically f_red = 1.0) days.

In Section 4.1.1, the evolution of physical quantities is described with several representative cases, and the most reasonable parameter set for one-zone calculations is selected. Section 4.1.2 is devoted to explaining important chemical reactions with the best parameter set as a representative.

4.1.1. Overall Trends of Physical Quantities and the Selection of Reasonable Parameter Sets

In the following subsections (Sections 4.1.1.1–4.1.1.3), the results of the representative three cases are presented as candidates of the reasonable parameter set, respectively. Among all the calculated cases, the three cases are arbitrarily selected according to the three criteria described in Section 4.1.1.4; then, among the three cases, the most reasonable one is selected.

4.1.1.1. The Case with f_h = 10⁻⁴, f_d = 5 × 10⁻², and t_s = 200 days

In the top panels in Figure 9, the results with the parameter values of f_h = 10⁻⁴, f_d = 5 × 10⁻², and t_s = 200 days are shown. The gas temperature evolution is plotted in the left panel (violet line). The temperature has a break at around 14 days, and soon after (∼20 days), the temperature goes down to ∼10⁴ K at which molecules start to form as seen in the amounts of CO and SiO in the same panel. In this case, the gas temperature quickly goes down to below ∼100 K at ∼200 days due to the CO line cooling after around 100 days.

By using this case as an example, the population of vibrational levels of CO and the evolution of the gas density are described in Figure 10.

**Figure 10.** One-zone calculation results for the case of f_h = 10⁻⁴, f_d = 5 × 10⁻², and t_s = 200 days. Top panel: the time evolution of the population of CO vibrational levels X_i and the mean escape probability $\overline{\beta }$ . Bottom panel: the time evolution of the gas density; for reference, a power-law evolution with the power of −3 is also shown. For the details, see the text.
Download figure:
Standard image High-resolution image

In the top panel, the time evolution of the population of vibrational levels of CO is shown. Initially, the population seems to be determined almost as an equilibrium state because of the high density and temperature at this phase. Throughout the evolution, the lower the vibrational level is, the higher the populational level is; the level v = 0 is dominated at all phases, but the higher levels are also populated in early phases until a few hundred days. The time evolution of the population of the levels is sensitively determined by the balance of the excitation (absorption) and deexcitation (emission) processes, which depends on the evolution of the gas density and temperature, the background radiation field, and the escape probabilities (line optical depths). The line optical depths depend on the gas density and the population of the levels; the evolution is complicated due to the dependencies. In the same panel, the mean escape probability $\overline{\beta }$ obtained by the summation of the effective escape probabilities f_red β_{κ
λ} (see Equations (29) and (42)) for all rovibrational lines with the weight of the population of the initial state is also shown. The value of the reduction factor f_red at 100 days is approximately 10⁻² for t_s = 200 days. Therefore, the escape probabilities at early phases are effectively reduced by the reduction factor f_red to restrict CO line emissions and the resultant cooling. As mentioned above, the gas temperature quickly goes down to ∼100 K at ∼200 days even with the effective reduction through the factor f_red.

It is found that, in one-zone calculations, overall, f_red = 200 days results in better fluxes of emission via CO rovibrational transitions. The impact of the reduction factor f_red through the parameter t_s is discussed in Appendix D by taking several one-zone calculations as a reference.

As seen in the time evolution of the gas density in the bottom panel of Figure 10, the value increases after 100 days due to the effect described in Equation (41), i.e., the density evolution deviates from the power-law evolution with the power of n = − 3 due to the deflation that stems from the cooling. The timing of the density increase corresponds to the rapid temperature decrease due to the cooling by CO rovibrational transitions around 100 days. It is noted that, after 200 days, the density evolution again gradually approaches the original power-law evolution. This is because, in the calculation, the minimum temperature (specific internal energy) is set. Then, after the gas temperature reaches the minimum temperature (100 K), without effective heating due to the decay of ⁵⁶Ni, e becomes equal to e_ad in Equation (41) at about 300 days, at which point the evolution of the gas density becomes consistent with the original power-law evolution by the assumption.

In the top right panel in Figure 9, the fluxes of CO vibrational bands derived from the calculation and observed ones (open and closed points; the former points are in the fundamental band, and the latter ones are in the first overtone band) are plotted. Only before about 100 days, the calculated flux of the first overtone band (Δv = 2) is slightly higher than that of the fundamental band (Δv = 1). The calculation underestimates the fluxes of the fundamental and the first overtone bands a bit after 100 days compared with the observed fluxes. In particular, it is difficult to explain the peaks around 200–250 days and the comparable flux levels between the two bands around the peaks seen in the observed fluxes. Such difficulties are also common for other calculation models in this paper. Therefore, we would expect that supra-thermal electrons supplied by the decay of ⁵⁶Ni may play a role in the excitation. Actually, the timings of the peaks of the observed fluxes are consistent with one at which the gas temperatures peak due to the heating by the decay of ⁵⁶Ni²³ (the feature of the peaks in the gas temperatures due to the decay of ⁵⁶Ni can be more recognized in the 1D calculation results seen, e.g., in the left panel in Figure 15). It is noted that, as mentioned in Section 1, the observed CO vibrational bands are well reproduced by a one-zone model with time-dependent thermal–chemical calculations (Liu & Dalgarno 1995), in which the thermal evolution was calculated with a more realistic method including a time-dependent heating efficiency of gas and the line cooling not only from CO rovibrational transitions but also from ones of ions with nonthermal excitations, although an arbitrary time-dependent reduction factor for the energy deposition was introduced for the reproduction (Liu & Dalgarno 1995). However, since the reproduction of the observed CO light curves is beyond the scope of this paper, we leave these issues related to CO rovibrational transitions for future works (see Section 5.2).

The time evolution of the amounts of CO and SiO and the seed atoms, carbon, oxygen, and silicon, is shown in the top left panel in Figure 9 with the estimations (including theoretical calculations) for CO and SiO by the previous studies (closed points for CO, and open points for SiO, Liu et al. 1992; Liu & Dalgarno 1994, 1995; Liljegren et al. 2020) and ones from the ALMA observations (vertical bars with errors, Matsuura et al. 2017). It is noted that some of the estimated values in the previous study (Liu et al. 1992) are based on the fitting of observed spectra, but in the fitting, some assumptions (e.g., steady states, LTE, and NLTE) and models are involved. Therefore, the estimated values in the previous studies may not necessarily be the actual amounts as can be seen from the fact that the estimated values for CO are not fully consistent among the studies. Nevertheless, those values are useful for reference. The results of the molecules other than CO and SiO are presented later in some of 1D calculation results (see Sections 4.2.3 and 4.2.4).

As seen in the top left panel, CO and SiO molecules start to form at about 30 days quickly after the gas temperature reaches ∼10⁴ K, and the amounts of both CO and SiO reach over 0.1 M_⊙ at 40 days. Since there were no observations (at least for SN 1987A) before 100 days, unfortunately, it is difficult to set a constraint on the amounts of molecules before 100 days. After 100 days, some fractions of CO and SiO are destructed until 300 days (the major formation and destruction reactions of CO and SiO are discussed later in Section 4.1.2). The amount of CO is consistent in particular with the values of LD95 (Liu & Dalgarno 1995), and the amount of SiO is also consistent with the estimations of LD94 (Liu & Dalgarno 1994) relatively well. After 300 days, CO and SiO are recovered again. At the final phase of ∼10⁴ days, the amount of CO is a bit less than the estimated range from ALMA; the amount of SiO is overestimated compared with the ALMA estimation. It is noted that, in the ejecta of SN 1987A, the existence of a large mass of dust (0.3 M_⊙ amorphous carbon and 0.5 M_⊙ silicate) has been suggested by the observations by Herschel (Matsuura et al. 2015). Thus, in reality, some fractions of molecules, in particular, SiO, must have been converted to dust as partly supported by dust formation theories in CCSNe (Nozawa et al. 2003; Nozawa & Kozasa 2013) and the fact that the mass of CO estimated by the ALMA observations (Matsuura et al. 2017) is higher than that of SiO. Therefore, some degree of overestimations of SiO compared with the ALMA estimation can be explained in this way. Then, the calculated amount of CO and SiO would be roughly consistent with the previous studies and the ALMA in this specific case.

As mentioned above, the gas temperature, however, quickly goes down to ∼100 K in spite of the underestimation of the fluxes of the CO vibrational bands. Actually, in the previous studies (Liu et al. 1992; Liu & Dalgarno 1995; Liljegren et al. 2020),²⁴ the temperatures of carbon monoxide have been calculated as >500 K before 1000 days. Since the adopted models and assumptions are different among the studies including this paper, such a rapid cooling would not be necessarily rejected. Though, we prefer to find a milder gas cooling case (gas temperature evolution) as presented below.

4.1.1.2. The Case with f_h = 5 × 10⁻⁴, f_d = 10⁻², and t_s = 200 days

In the middle panels in Figure 9, the results of the parameter values, f_h = 5 × 10⁻⁴, f_d = 10⁻², and t_s = 200 days, are shown. In this case, the gas temperature (the left panel) remains greater than 100 K before 1000 days in contrast to the previous case. After approximately 300 days, the gas temperature increases by the heating due to the decay of radioactive ⁵⁶Ni. As observed in the previous case, the flux of the first overtone band dominates that of the fundamental band before about 100 days. The monotonically increasing trend both for the fundamental and the first overtone bands before 100 days is different from the previous case with a dip around 70 days because of a bit different evolution of the gas temperature (density). Compared with the previous case, the fluxes (the right panel) of the CO fundamental and first overtone bands are slightly improved to be more consistent with the observations after 200 days, although at around 100 days the fluxes are a bit overestimated, and after 200 days the fluxes are still underestimated by more than 1 order of magnitude.

The amount of CO (the left panel) before 100 days is not so different from the previous case, but the amount after 100 days is apparently higher than the previous case and the estimations by the previous studies. This trend is partly due to the lower f_d value than that of the previous case. The amount of SiO also has a similar trend. However, as can be seen, the final amounts of CO and SiO are not so different from the ones in the previous case.

4.1.1.3. The Case with f_h = 10⁻³, f_d = 10⁻², and t_s = 200 days

In the bottom panels in Figure 9, the results of the parameter values of f_h = 10⁻³, f_d = 10⁻², and t_s = 200 days are shown. The behavior of the gas temperature (the left panel) is similar to the previous case, but the value remains higher than those of the previous two cases after 100 days and before 1000 days due to the higher heating efficiency via the f_h value. The early trend (before 100 days) of the fluxes of CO vibrational bands is similar to the previous case; the high gas temperature, however, slightly increases the fluxes of CO vibrational bands (the right panel) compared with the previous two cases.

The amount of CO (the left panel) before 100 days is not so different among the three cases including this case. However, the amount after 100 days is further shifted to higher values, and the deviation from the estimations by the previous studies becomes significant. The amount of SiO is also a similar trend to the one of CO. Although the deviation of the amounts of CO and SiO from those estimated in the previous studies with data analysis and theoretical calculations is significant, the final amounts are more consistent with the ALMA estimations. The amount of CO at the final phase is at this time within the error bars of the ALMA estimation, and the amount of SiO becomes close to the ALMA one.

4.1.1.4. The Criteria for Selecting the Best Parameter Set

In order to pick up reasonable parameter sets, the following three points are taken into account as the criteria to be satisfied: (1) the consistency of the amounts of CO and SiO with the previous estimations and the ALMA values, i.e., points and bars with errors shown in the left panels in Figure 9; (2) the consistency of the calculated fluxes of CO vibrational bands with the observed fundamental and first overtone bands; (3) the calculated gas temperatures should not be too low (∼100 K) before 1000 days. Among the three cases presented above and the other models calculated, we arbitrarily select the most acceptable model to be the second case (the case in Section 4.1.1.2) presented above according to the three points above. Hereafter, the impact of the parameters, f_h and f_d, and trends of the other models are briefly described with a few representative/extreme cases.

If the values of f_h are high (f_h = 5 × 10⁻³ or 10⁻²), i.e., the efficiency of the gas heating by the decay of ⁵⁶Ni is high, the gas temperature goes down to ∼10⁴ K at late phases (∼1000 days). Then, molecules start to form at such late phases (the left panel in Figure 11), which is inconsistent with the fact that CO vibrational bands were observed as early as ∼100 days. Additionally, at the late phases, the gas densities are low (∼10⁻¹⁵ g cm⁻³) compared with the ones at ∼100 days (∼10⁻¹² g cm⁻³; see the bottom panel in Figure 10); the amounts of molecules tend to be lower due to the low efficiency of two-body reactions at later phases compared with the cases of lower f_h values.

In the case of high f_d values, i.e., the destruction efficiency is high, CO and SiO are destructed from ∼100 to 300 days (the right panel in Figure 11) due to the ionizations by Compton electrons (CM reactions) and some reactions related to ionized atoms for this specific case (for the details of important chemical reactions for the formation and destruction of CO and SiO, see Section 4.1.2), although after 300 days the amounts are recovered to some extent. In the destruction processes, ionized atoms could also play a role through ion–neutral reactions as described in the next section.

It should be noted that, in the one-zone approximation, the chemical composition is forced to be uniform throughout the ejecta core; matter mixing is then overestimated compared to the real environment. Since both the gas heating and the ionization and destruction by Compton electrons (CM reactions) depend on the local amount of ⁵⁶Ni (see Equations (9), (10), (11), and (16)), these processes may also be overestimated with the effectively high local ⁵⁶Ni abundance. Therefore, in the case of the one-zone approximation, lower values of f_h and f_d may be selected as reasonable parameter values compared with the cases with more realistic matter mixing as discussed later with 1D calculation results.

4.1.2. Important Chemical Reactions for CO and SiO

In this section, with the most acceptable case (the case in Section 4.1.1.2; the middle panels in Figure 9), i.e., f_h = 5 × 10⁻⁴, f_d = 10⁻², and t_s = 200 days, as a representative, important chemical reactions for the formation and destruction of molecules are described focusing on CO and SiO. To pick up important reactions of which the contribution to the formation/destruction is high, the first and second terms in the right-hand side in Equation (1), more specifically F_i and D_i in Equation (2), are monitored. Molecules start to form as early as approximately 30 days after the explosion (see the middle left panel in Figure 9). The important reactions are picked up in chronological order for CO and SiO, separately (first CO, and later SiO). As seen in the evolution of the amounts of CO and SiO, the evolution can roughly be separated into four epochs, i.e., from the start of the formation of molecules to 40 days (the initial rapid rise), from 40 to 100 days (a plateau phase), from 100 to 300 days (a destruction dominant phase), and after 300 days (a recovery phase). In Figure 12, the contributing chemical reactions for the formation and destruction of CO in which CO is directly involved are shown as a function of time (after 30 days) for reference.

**Figure 12.** Contributing chemical reactions for the formation (left) and destruction (right) of CO in which carbon monoxide (CO) is directly involved as a function of time (after 30 days) for the one-zone calculation of the case with f_h = 5 × 10⁻⁴, f_d = 10⁻², and t_s = 200 days. The code inside the parentheses left side of each reaction denotes the corresponding reaction type listed in Table 2. Colors denote the relative contribution of each reaction, which is proportional to the reaction flow F_i (D_i) in Equation (1), normalized as the maximum F_i (D_i) among CO formation (destruction) reactions at each time to be unity. For neutral–neutral (`NN`) reactions, net contributions are counted by subtracting the reaction flow of the corresponding inverse reaction. Reactions are picked up if the relative contribution once becomes greater than 10⁻³.
Download figure:
Standard image High-resolution image

At the early phase from the start of the molecule formation (≲30 days), basically, the formation processes are dominant. The most dominant formation process of CO is initially the radiative association (RA) reaction below

$\begin{eqnarray}&&{\rm{C}}+{\rm{O}}\to \mathrm{CO}+\gamma .\end{eqnarray} \tag{ 45 }$

The secondary reactions are the neutral–neutral (NN) reactions,

$\begin{eqnarray}&&{\rm{O}}+\mathrm{CH}\to \mathrm{CO}+{\rm{H}},\end{eqnarray} \tag{ 46 }$

$\begin{eqnarray}&&{\rm{C}}+\mathrm{OH}\to \mathrm{CO}+{\rm{H}},\end{eqnarray} \tag{ 47 }$

where CH is dominantly formed by a radiative association reaction of the same form of Equation (45), and OH is dominantly produced by the NN reaction, O + H₂ → OH + H, where H₂ is produced by the associative electron detachment (AD) reaction, H⁻ + H → H₂ + e⁻; H⁻ can be produced by the radiative electron attachment (REA) reaction, H + e⁻ → H⁻ + γ. It is noted that, in most of the NN reactions adopted in this paper, there are corresponding inverse reactions, and the forward and inverse reactions compete with each other. Therefore, for NN reactions, the net contributions are taken into account in judging whether the reaction is important or not. Then, the contributions of the RA reaction in Equation (45) and the NN reaction in Equation (46) are gradually reduced, and the contribution of the NN reaction in Equation (47) instead increases (actually, the contribution of the NN reaction is initially negative, i.e., the inverse reaction is dominant) to be comparable with the RA reaction at ∼40 days. On the other hand, the main destruction processes before 40 days are the NN reactions below

$\begin{eqnarray}&&{\rm{H}}+\mathrm{CO}\to \mathrm{OH}+{\rm{C}},\end{eqnarray} \tag{ 48 }$

$\begin{eqnarray}&&\mathrm{Si}+\mathrm{CO}\to \mathrm{SiO}+{\rm{C}}.\end{eqnarray} \tag{ 49 }$

Initially, the NN reaction in Equation (48) dominates the reaction in Equation (49); the latter, however, becomes the primary destruction reaction after 30 days as the contribution of the former decreases to be negative after 35 days. The latter reaction is also the formation reaction of SiO. As seen in the middle left panel in Figure 9, by approximately 40 days, a large fraction of the seed carbon and oxygen atoms are consumed.

Then, after about 40 days, the dominant formation process of CO is replaced by the NN reaction in Equation (47) (until 80 days); the RA reaction in Equation (45) becomes no longer the primary in the formation processes. One of the secondary formation reactions is the NN reaction in Equation (46). As secondary reactions, the NN reactions below also partake in the formation processes

$\begin{eqnarray}&&{\rm{O}}+\mathrm{CN}\to \mathrm{CO}+{\rm{N}},\end{eqnarray} \tag{ 50 }$

$\begin{eqnarray}&&{\rm{O}}+\mathrm{SiC}\to \mathrm{CO}+\mathrm{Si},\end{eqnarray} \tag{ 51 }$

where SiC is mainly produced by the RA reaction, Si + C → SiC + γ, at an early phase; and later, it is produced by the NN reaction, C + SiH → SiC + H; SiH is produced by the NN reaction, Si + H₂ → SiH + H. The formation processes of CN will be mentioned later. On the other hand, the primary destruction process becomes the NN reaction in Equation (49) after 40 days. As one of the secondary destruction processes, the NN reaction below starts to contribute to the destruction

$\begin{eqnarray}&&{\rm{N}}+\mathrm{CO}\to \mathrm{NO}+{\rm{C}}.\end{eqnarray} \tag{ 52 }$

It is noted that CN seen in the reactants in Equation (50) is mainly produced by the NN reaction, C + NO → CN + O; NO can also be seen in the products in Equation (52). Therefore, through NN reactions, several molecules are complexly involved in the formation and destruction of CO, i.e., one of the products of the destruction process of CO is indirectly involved in the formation processes of CO. After 50 days, the contributions of the ion–neutral (IN) reactions below

$\begin{eqnarray}&&{{\rm{C}}}^{+}+\mathrm{SiO}\to \mathrm{CO}+{\mathrm{Si}}^{+},\end{eqnarray} \tag{ 53 }$

$\begin{eqnarray}&&{\mathrm{He}}^{+}+\mathrm{CO}\to {\rm{O}}+{{\rm{C}}}^{+}+\mathrm{He},\end{eqnarray} \tag{ 54 }$

become gradually larger and larger as one of the formation and destruction processes, respectively. The former reaction is also a destruction process of SiO. The ions C⁺ and He⁺ can be produced by ionization of the corresponding atoms by Compton electrons in Equation (17) (CM reactions). After 60 days, the former reaction becomes the secondary formation process, and the latter process becomes the primary destruction process. The contributions of CM reactions, i.e., the ionization and destruction by Compton electrons due to the decay of ⁵⁶Ni also gradually increase as destruction processes. In particular, the ionization process of CO by Compton electrons in Equation (18) and the destruction by UV photons (UV reaction) become the secondary destruction processes. The NN reaction below also partakes in the destruction processes

$\begin{eqnarray}&&{\rm{O}}+\mathrm{CO}\to {{\rm{O}}}_{2}+{\rm{C}}.\end{eqnarray} \tag{ 55 }$

As CO⁺ increases by the ionization in Equation (18), the charge exchange (CE) reaction below becomes one of the secondary formation processes

$\begin{eqnarray}&&{\rm{H}}+{\mathrm{CO}}^{+}\to \mathrm{CO}+{{\rm{H}}}^{+}.\end{eqnarray} \tag{ 56 }$

After 80 days, the IN reaction in Equation (53), i.e., the conversion of SiO to CO, overtakes the NN reaction in Equation (47). The situation does not change until about 100 days.

After 100 days, the destruction processes are dominant as seen in the middle left panel in Figure 9. The primary destruction process is still the IN reaction via He⁺ in Equation (54), and the ionization by Compton electrons in Equation (18) and the dissociation by UV photons are secondary (the situation on the destruction processes continues until a few thousand days). The NN reaction in Equation (47) again becomes the primary formation process by overtaking the IN reaction in Equation (53). After approximately 130 days, the NN reactions below also partake in the secondary formation processes as well as the NN reaction in Equation (46)

$\begin{eqnarray}&&{\rm{C}}+{{\rm{O}}}_{2}\to \mathrm{CO}+{\rm{O}},\end{eqnarray} \tag{ 57 }$

$\begin{eqnarray}&&{\rm{O}}+\mathrm{CS}\to \mathrm{CO}+{\rm{S}},\end{eqnarray} \tag{ 58 }$

$\begin{eqnarray}&&{\rm{C}}+\mathrm{SO}\to \mathrm{CO}+{\rm{S}}.\end{eqnarray} \tag{ 59 }$

Here, O₂, CS, and SO are mainly produced by the NN reactions, O + OH → O₂ + H, C + SO → CS + O, and S + OH → SO + H, respectively. It is interpreted that, after 100 days, those molecules are converted to CO through the NN reactions in Equations (57), (58), and (59). Among the secondary formation reactions, i.e., the IN reaction in Equation (53), the CE reaction in Equation (56), and the NN reactions in Equations (46), (57), (58), and (59), the significance of each reaction depends on time. However, overall, the situation continues until 300 days.

After 300 days, the formation processes overwhelm the destruction processes. The primary formation process is still the NN reaction in Equation (47), and the secondary formation reactions are the CE reaction in Equation (56) (until about 600 days) and the NN reactions in Equations (46), (57), (58), and (59). The primary destruction process is the IN reaction in Equation (54), and the secondary destruction processes are the dissociation by UV photons and the ionization and destruction by Compton electrons in Equations (18) and (19) until a few thousand days. After that, mostly, the NN reaction below contributes to the destruction, although the amount of CO does not vary so much at this later phase

$\begin{eqnarray}&&\mathrm{OH}+\mathrm{CO}\to \mathrm{CH}+{{\rm{O}}}_{2}.\end{eqnarray} \tag{ 60 }$

Hereafter, similarly to CO, the formation and destruction processes of SiO are described. In Figure 13, contributing chemical reactions for the formation and destruction of SiO in which SiO is directly involved are shown as a function of time for reference. As seen in the left panel, the contributing SiO formation reactions are limited compared with those for CO.

From the start of the molecule formation (∼30 days) to 40 days, the formation processes are dominant as seen in the middle left panel in Figure 9. Initially, the primary formation process is the RA reaction below

$\begin{eqnarray}&&\mathrm{Si}+{\rm{O}}\to \mathrm{SiO}+\gamma .\end{eqnarray} \tag{ 61 }$

The secondary formation reactions are the NN reaction in Equation (49) mentioned as a destruction process of CO and the NN reaction below

$\begin{eqnarray}&&\mathrm{Si}+\mathrm{OH}\to \mathrm{SiO}+{\rm{H}},\end{eqnarray} \tag{ 62 }$

where OH is mostly produced by the NN reaction, O + H₂ → OH + H, as mentioned above. The destruction processes of the initial phase are thermal fragmentation reactions (TF reactions) with hydrogen and helium described in Section 2.1.1. The secondary destruction process is the NN reaction below

$\begin{eqnarray}&&{\rm{N}}+\mathrm{SiO}\to \mathrm{SiN}+{\rm{O}}.\end{eqnarray} \tag{ 63 }$

After 35 days, the NN reaction in Equation (49) (CO destruction process) overtakes the RA reaction in Equation (61) to become the primary formation process, and the NN reaction below is added to the secondary formation reactions

$\begin{eqnarray}&&{\rm{O}}+\mathrm{SiH}\to \mathrm{SiO}+{\rm{H}},\end{eqnarray} \tag{ 64 }$

where SiH is primarily produced by the NN reaction, Si + H₂ → SiH + H, at this phase.

After 40 days, the formation and destruction processes compete with each other. The primary formation process is the NN reaction in Equation (49) for a while (until about 55 days), and the contribution of the RA reaction in Equation (61) becomes gradually small. The destruction process (NN reaction) in Equation (63) overtakes TF reactions of SiO with hydrogen and helium. After 50 days, the NN reaction below partakes in the secondary formation process

$\begin{eqnarray}&&{\rm{O}}+\mathrm{SiC}\to \mathrm{SiO}+{\rm{C}},\end{eqnarray} \tag{ 65 }$

where SiC is mainly formed by the RA reaction between silicon and carbon and the NN reaction, C + SiH → SiC + H, as mentioned above. As the primary destruction processes, the NN reactions below are added by overtaking the NN reaction in Equation (63)

$\begin{eqnarray}&&{\rm{H}}+\mathrm{SiO}\to \mathrm{SiH}+{\rm{O}},\end{eqnarray} \tag{ 66 }$

$\begin{eqnarray}&&{\rm{H}}+\mathrm{SiO}\to \mathrm{OH}+\mathrm{Si}.\end{eqnarray} \tag{ 67 }$

Simultaneously, the contributions from the IN reaction in Equation (53) (formation process of CO) and the CE and IN reactions below gradually increase as secondary destruction processes

$\begin{eqnarray}&&{{\rm{H}}}^{+}+\mathrm{SiO}\to {\mathrm{SiO}}^{+}+{\rm{H}},\end{eqnarray} \tag{ 68 }$

$\begin{eqnarray}&&{\mathrm{He}}^{+}+\mathrm{SiO}\to \mathrm{Si}+{{\rm{O}}}^{+}+\mathrm{He},\end{eqnarray} \tag{ 69 }$

$\begin{eqnarray}&&{\mathrm{He}}^{+}+\mathrm{SiO}\to {\mathrm{Si}}^{+}+{\rm{O}}+\mathrm{He}.\end{eqnarray} \tag{ 70 }$

Soon after that, the contribution of the NN reaction in Equation (67) as the destruction process becomes small (the inverse reaction in turn dominates the forward reaction). After 60 days, the NN reaction in Equation (62) overtakes the NN reaction in Equation (49) to be the primary formation process. The contributions of the NN reactions below also gradually increase as the secondary formation processes

$\begin{eqnarray}&&\mathrm{Si}+{{\rm{O}}}_{2}\to \mathrm{SiO}+{\rm{O}},\end{eqnarray} \tag{ 71 }$

$\begin{eqnarray}&&\mathrm{Si}+\mathrm{NO}\to \mathrm{SiO}+{\rm{N}},\end{eqnarray} \tag{ 72 }$

where O₂ is mainly produced by the NN reaction, O + OH → O₂ + H, and NO is formed by the NN reaction, O + SiN → NO + Si; SiN is produced by the NN reaction in Equation (63). After 80 days, the CE reaction in Equation (68) overtakes the NN reaction in Equation (66) to become the primary destruction process. Overall, the situation continues until 100 days.

After 100 days, the destruction processes are dominant as seen in the middle left panel in Figure 9. The primary formation process is the NN reaction in Equation (62) (until 140 days). The primary destruction process is the CE reaction in Equation (68); the NN reaction in Equation (66) and the IN reactions in Equations (53), (69), and (70) are followed as the secondary destruction processes. After 140 days, the NN reaction in Equation (71) overtakes the NN reaction in Equation (62) to be the primary formation process. The contribution of the formation processes, the RA reaction in Equation (61) and the NN reaction in Equation (65), recovers again to some extent. The contribution of the destruction reactions, the NN reaction in Equation (66) and the IN reaction in Equation (53), becomes smaller. After that, the situation overall continues until 300 days.

After 300 days, the formation processes become dominant until the reactions are settled (∼1000 days). The contributing formation and destruction processes, however, are overall the same as the situation at ≲300 days. After ∼1000 days, the contribution of the formation process, the NN reaction in Equation (62), becomes a bit small. The destruction process, the NN reaction in Equation (66), becomes the primary process by overtaking the CE reaction in Equation (68).

In this section, by taking the representative case, important reactions for CO and SiO were described. For some of the reactions, in particular, the NN and IN reactions, listed above, hydrogen and helium atoms and the corresponding ions, H⁺ and He⁺, are involved. Such reactions play a nonnegligible role in the formation and destruction processes of CO and SiO. In the one-zone approximation, as partly mentioned in the previous section, the mixing of hydrogen and helium into the ejecta core and the mixing of ⁵⁶Ni, which is important for the ionization and destruction of molecules, into outer layers would be overestimated. Therefore, the description in this section may not be universal, but as a reference, it may still be useful. Actually, most of the reactions related to the CO and SiO formation and destruction referred to in the later sections are covered above. The impact of matter mixing with more realistic 1D models on molecule formation is discussed in the next section.

4.2. 1D Calculation Results

In this section, 1D calculation results are presented according to the strategy described in Section 3. In Section 4.2.1, as in the one-zone calculations, with several cases, the evolution of physical quantities is described, and the most reasonable parameter set is selected. In Section 4.2.2, the results of the angle-averaged 1D profiles of the 3D models with the binary merger and the single-star progenitor models are presented including the description of the molecules other than CO and SiO. Section 4.2.3 is devoted to the comparison with the spherical cases. Finally, in Section 4.2.4, differences in the specified directions in the 3D models are presented.

4.2.1. The Fiducial Case and the Dependence on the Parameters, f_h, f_d, and t_s

As done in Section 4.1 for one-zone calculations, to find a reasonable set (reasonable sets) of parameters, f_h, f_d, and t_s, for 1D calculations, based on the angle-averaged 1D profiles of the model b18.3-high (Ono et al. 2020), in total, 100 cases are calculated with the combinations of the parameters, f_h = 10⁻⁴, 5 × 10⁻⁴, 10⁻³, 5 × 10⁻³, 10⁻², f_d = 10⁻², 5 × 10⁻², 10⁻¹, 5 × 10⁻¹, 1.0, t_s = 200, 300, 500, and ∞ (practically f_red = 1.0) days.

In the following subsections, the results of two representative cases are presented in Sections 4.2.1.1 and 4.2.1.2, respectively, as candidates for the best parameter set to be fixed for the latter discussion. Then, in Section 4.2.1.3, with a few additional cases, the parameter dependence is described, and among the two cases, the most reasonable parameter set is selected. For the selection of the best parameter set (the two candidates), the three criteria listed in Section 4.1.1.4 are taken into account. Section 4.2.1.4 is devoted to the description of the important chemical reactions for CO and SiO for the reasonable parameter set as a representative.

4.2.1.1. The Case with f_h = 10⁻², f_d = 10⁻², and t_s = 500 days

In the upper panels in Figure 14, the results with the parameters f_h = 10⁻², f_d = 10⁻², and t_s = 500 days are shown as one of the candidates for the reasonable parameter set.

As seen in the calculated fluxes of CO vibrational bands (the right panel), the peak flux level (∼10⁻⁸ erg s⁻¹ cm⁻²) of the fundamental band (Δv = 1) is comparable with the observed peak fluxes, but the calculation fails to reproduce the timing (200–250 days), and the dome-like feature of the observed peaks the same as in the one-zone calculations. The qualitative feature of the fluxes of the other bands is similar to that of the fundamental band, but the higher Δv is, the lower the flux; the flux of the first overtone band (Δv = 2) is rather underestimated compared with the observed one even in early phases before 100 days. In the one-zone calculation results presented in Section 4.1.1 (see the right panels in Figure 9), the flux of the first overtone band can be comparable or even higher compared with that of the fundamental band before 100 days. The discrepancy between the one-zone and 1D calculations may partly be attributed to the more effective mixing of ⁵⁶Ni in the one-zone calculations than that in the 1D calculations and the resultant higher electron densities due to more effective ionization by Compton electrons, which may enhance the excitation of higher CO vibrational levels. The increase of fluxes after 500 days is the contribution from some inner particles in which the formation of CO and CO rovibrational transitions start at such later phases.

The calculated flux levels depend on the parameter t_s as demonstrated in Appendix D for one-zone calculations. The higher (lower) t_s value, the higher (lower) the peak flux levels. Among the calculated cases, the parameter value of t_s = 500 days results in better flux levels not only for this case but also for the other 1D calculations (the better value is different from one for the one-zone calculations, i.e., t_s = 200 days). Then, hereafter, only the cases of the parameter of t_s = 500 days are focused (if not specifically mentioned, t_s = 500 days).

In Figure 15, the time evolution of the gas temperatures (left) and densities (right) are shown by taking this specific parameter set as one of the representatives. The gas temperatures of inner (tracer) particles are efficiently heated due to the decay of ⁵⁶Ni making peaks around 200 days. The particles whose initial positions (radii) are less than 3.5 × 10¹³ cm are more or less affected by the heating at some point. Some particles around 1.5 × 10¹³ cm undergo the cooling by CO rovibrational transitions around 100 days; the gas temperature evolution of such particles is a bit complicated since the heating and cooling are competing with each other. Some of the inner particles also go through the cooling at later phases (≲1000 days). The evolution of the gas densities fluctuates on top of the power-law evolution of the power of −3 due to the effect described in Equation (41). The behavior of the particles affected by the cooling is a bit complex reflecting the temperature evolution. The gas densities recover the original power-law evolution after the gas temperatures reach the minimum (100 K) by the assumption (see Section 2.3).

The time evolution of the total amounts (contributions from all the particles) of CO and SiO (the left panel) shows early (∼10 days) formation of CO and SiO and smoother increase of the amounts than those in the one-zone calculation models (see, e.g., the left panels in Figure 9). This is because of the variations in physical quantities among the particles as seen in Figure 15. Both the amounts of CO and SiO reach about 10⁻²–10⁻¹ M_⊙ by 100 days in contrast to the representative one-zone calculations (Figure 9), in which both the amounts exceed 10⁻¹ M_⊙ at the peak. Actually, the seed atoms for CO and SiO are not so consumed after 40 days compared with the representative one-zone cases (Figure 9). The calculated amounts of CO and SiO are overestimated roughly by 1 order of magnitude compared with the previous studies (points), and the decreasing trend of CO from 100 days to a few hundred days seen in LD95 (filled circles) is not well reproduced. After 100 days, in particular, SiO decreases by 1 order of magnitude by 600 days; CO slightly decreases during this phase. After that, the amount of SiO recovers by a factor of a few. The increasing and decreasing trends are overall milder than those of the representative one-zone cases. This is probably attributed to the wide range of physical quantities among the particles as mentioned and less effective matter mixing in the 1D calculations.

4.2.1.2. The Case with f_h = 5 × 10⁻³, f_d = 1.0, and t_s = 500 days

In the lower panels in Figure 14, the 1D calculation results with the parameters, f_h = 5 × 10⁻³, f_d = 1.0, and t_s = 500 days, are shown as another candidate for the reasonable parameter set.

The qualitative features of the fluxes of CO vibrational bands are similar to the previous case; the fluxes peak before 100 days in contrast to the observed fundamental and the first overtone bands; the higher Δv is, the lower the flux is throughout the evolution; fluxes increase to some extent after 500 days. The peak flux of the fundamental band is slightly lower than that of the previous case; instead, the flux is a bit higher and smoother at around 100 days.

The qualitative features of the amounts of CO and SiO (the left panel) are similar to the previous case, but quantitative differences can be recognized. The initial increase of CO and SiO is rather similar to the previous case. However, the amount of CO becomes greater than 10⁻¹ M_⊙ at the peak, and the amount of SiO nearly reaches 10⁻¹ M_⊙. After 100 days, both CO and SiO have a clear decreasing trend compared with the previous case; by 300 days, the amount of CO decreases by a factor of 5 in contrast to the previous case, and the amount of SiO is smaller than that of the previous case at this phase. After 300 days, the amount of CO only slightly recovers, and the amount of SiO recovers by a factor of more than a few. The final amount of CO is comparable with the previous case, and the amount of SiO is slightly greater than that of the previous case.

4.2.1.3. The Selection of the Best Parameter Set and the Parameter Dependence

To find the reasonable parameter set, the three criteria to be satisfied mentioned in Section 4.1.1.4 are taken into account for 1D calculations too. In the 1D calculations, the third point, which is related to the gas temperature, is not so problematic as long as the calculations are performed with the reduction factor f_red with t_s = 500 days. Actually, outer particles tend to reach such a low temperature at a few hundred days even without the cooling by CO rovibrational transitions; unless the majority of the particles have such an immediate low temperature, we do not make an issue out of the model. Then, the first and second points matter. In this regard, among the calculated 1D cases, the two cases presented above in this section would be candidates for the reasonable parameter set to be fixed for later discussion.

Before fixing the parameter set, the overall features for other calculated cases and the dependence on the values of the parameters, f_h and f_d, are described. To see the impact of the parameter f_d, as the counterparts of the two cases above, the results with the f_d value of the opposite end of the two cases, i.e., f_d = 10⁻², are shown in Figure 16, respectively. As seen in the results with the parameters, f_h = 10⁻², f_d = 10⁻², and t_s = 500 days (upper panel), there is no distinct decrease (destruction) of CO and SiO around 100–300 days compared with the counterpart case with f_d = 1.0 (upper left panel in Figure 14). This trend can simply be understood as the consequence of the inefficient ionization and destruction by Compton electrons (CM reactions) due to the low value of f_d = 10⁻² and several ion-involving reactions, e.g., the CE reaction in Equation (68).

**Figure 16.** 1D calculation results with the angle-averaged profiles for the model b18.3-high (Ono et al. 2020) with the parameters, f_h = 10⁻², f_d = 10⁻², and t_s = 500 days (upper panel), and f_h = 5 × 10⁻³, f_d = 10⁻², and t_s = 500 days (lower panel). The time evolution of the total amounts of CO and SiO, and the seed atoms, carbon, oxygen, and silicon, are shown. The points are the estimations for CO and SiO in the previous studies; for the details, see the caption of Figure 14.
Download figure:
Standard image High-resolution image

In the case of f_h = 5 × 10⁻³, f_d = 10⁻², and t_s = 500 days (lower panel), the destruction of SiO after 100 days can be recognized in contrast to the previous case; the destruction is, however, less evident compared with the counterpart case (lower left panel in Figure 14). The destruction of SiO after 100 days is not due to CM reactions but due to a few NN reactions. The decreasing trend of CO seen in the counterpart case after 100 days disappears in this case.

Figure 17 shows the results for additional two cases of f_d = 1.0, and f_d = 10⁻², with f_h = 10⁻⁴ to see the impact of the parameter f_h; the value of the parameter f_h is the low end in this study. Then, the gas heating is inefficient compared with the cases shown above. As expected, the gas temperatures reach about 10⁴ K at earlier phases, and molecule formation starts in a high-density environment compared with the cases with higher f_h values. In the case of f_h = 10⁻⁴, and f_d = 1.0 (left panel in Figure 17), apparently, the amounts of CO and SiO at 40–100 days are the highest compared with the previous cases, and the amounts are almost comparable with the initial amounts of the seeds. However, after 100 days, due to the efficient CM reactions (destruction and/or ionization by Compton electrons), the amounts are reduced to some extent, and the final amounts are not dramatically different from those of the previous cases, e.g., one with f_h = 5 × 10⁻³, and f_d = 1.0 (lower left panel in Figure 14), although the amount of SiO is the highest compared to that of the previous cases. On the other hand, in the case of f_h = 10⁻⁴, and f_d = 10⁻² (right panel in Figure 17), basically, the inefficient CM reactions with the low f_d value result in the highest amounts of CO and SiO throughout the evolution among at least all the presented 1D cases above.

As a summary of the impacts of the parameters f_h and f_d on the amounts of CO and SiO, the parameter f_h mainly affects the amounts in early phases before 100 days, and the parameter f_d affects the destruction of CO and SiO from 100 days to a few hundred days. The lower f_h value is, the higher the amounts of CO and SiO before 100 days. The higher f_d value is, the more efficient destruction of CO and SiO at 100–300 days is. These qualitative features are the same as observed in one-zone calculations (see Section 4.1).

In view of the first criterion (described in Section 4.1.1.4), the cases of the f_h value less than 5 × 10⁻³ would result in too high amounts of CO and SiO at 40–600 days compared with the previous estimations. Then, among all the calculated 1D cases (for the model b18.3-high), again, the first and second cases (the cases in Sections 4.2.1.1 and 4.2.1.2, respectively) would be candidates for the fiducial case to be fixed. Among the estimations of the amounts of CO in the previous studies, only LD95 (Liu & Dalgarno 1995) solved the time-dependent chemical evolution with a similar method for CO rovibrational transitions with an optically thick regime as described in Section 2.2. Therefore, we would arbitrarily prefer to reproduce the estimations by LD95, in particular, the trend with a valley (i.e., the initial rise, the middle destruction, and the final recovery). In this regard, the second case has a clearer decreasing trend of the amounts of both CO and SiO, although the amount of CO at 100 days is a bit higher than the value by LD95 compared with the first case. The estimations in the previous studies are not consistent with each other, and in the estimations, more or less, some theoretical models are involved. Therefore, a strict quantitative consistency of the amounts of CO and SiO with the previous estimations may not be indispensable as mentioned. In view of the second criterion (described in Section 4.1.1.4), the fluxes of CO vibrational bands for the two cases are qualitatively similar to each other, and it is difficult to find an apparently better case than those in the two cases among all the 1D calculations. Thus, we would select the second case presented in Section 4.2.1.2, i.e., f_h = 5 × 10⁻³, f_d = 1.0, and t_s = 500 days as the fiducial case to be fixed for later discussion on the impact of effective matter mixing. Hereafter, if the parameter values of f_h, f_d, and t_s are not specifically mentioned, the values are 5 × 10⁻³, 1.0, and 500 days, respectively.

4.2.1.4. Important Chemical Reactions for CO and SiO for the Fiducial 1D Calculation Model

As done in Section 4.1.2 for the fiducial one-zone calculation model (middle panels in Figure 9 in Section 4.1.1.2), important chemical reactions for the formation and destruction of CO and SiO are described briefly for the fiducial 1D calculation model (b18.3-mean), i.e., the calculation with the angle-averaged 1D profile based on the 3D model b18.3-high (Ono et al. 2020) with the parameters, f_h = 5 × 10⁻³, f_d = 1.0, and t_s = 500 days, comparing with the fiducial one-zone case.

Since the contributing chemical reactions depend on the time and each particle, it is difficult to describe the formation and destruction processes of molecules in detail; nonetheless, by integrating the contributions from all the particles with the particle's mass as a weight, a similar discussion is presented here. In Figure 18 (Figure 19), contributing chemical reactions for the formation and destruction of CO (SiO) are shown for reference. As for CO formation (left panel in Figure 18), the listed contributing reactions in the figure are the same as the corresponding ones for the fiducial one-zone calculation (left panel in Figure 12); differences, however, can be recognized in the relative contributions between the fiducial one-zone and 1D calculations. As for the CO destruction and SiO formation and destruction, the listed reactions are a bit different from those in the fiducial one-zone model, and the relative contributions are also different to some extent. Hereafter, important chemical reactions are described for each representative epoch in chronological order, focusing on the epochs when the changes in the quantities of CO and SiO are large; the formation (destruction) reactions are focused for formation (destruction) dominant phases.

**Figure 18.** Contributing chemical reactions for the formation (left) and destruction (right) of CO in which carbon monoxide (CO) is directly involved as a function of time (after 30 days) for the model b18.3-mean. The code inside the parentheses left side of each reaction denotes the corresponding reaction type listed in Table 2. Colors denote the relative contribution of each reaction (contributions from all the tracer particles are counted by taking the masses as a weight), which is proportional to the reaction flow F_i (D_i) in Equation (1), normalized as the maximum F_i (D_i) among CO formation (destruction) reactions at each time to be unity. For neutral–neutral (`NN`) reactions, net contributions are counted by subtracting the reaction flow of the corresponding inverse reaction. Reactions are picked up if the relative contribution once becomes greater than 10⁻³.
Download figure:
Standard image High-resolution image

**Figure 19.** Same as Figure 18 but for silicon monoxide (SiO).
Download figure:
Standard image High-resolution image

From the start of the molecule formation to about 100 days (initial increasing phase), the contributing formation processes of CO are the RA reaction in Equation (45), the NN reactions in Equations (46) and (47), the CE reaction in Equation (56), and the NN reaction in Equation (50); compared with the one-zone case (see the left panel in Figure 12), the contribution of CE reaction in Equation (56) becomes large, and the contributions of the NN reaction in Equation (51) and the IN reaction in Equation (53) become small. The reactions contributing to the formation of SiO are the RA reaction in Equation (61) and the NN reactions in Equations (49), (62), (72), and (71); compared with the one-zone case, the contribution of the NN reactions in Equations (72) and (71) becomes large, and the NN reactions in Equations (49) and (64) become small.

The amounts of CO and SiO decrease from 100 days to approximately 300 days. In this phase, the contributing destruction processes of CO are the ionization and destruction by Compton electrons (CM reactions) described in Equations (18) and (19), the dissociation by UV photons (UV reaction), and the NN reaction in Equation (49). It is noted that, in the one-zone case, the IN reaction with He⁺ in Equation (54) is the primary destruction process in contrast to the 1D case here. This feature in the one-zone case may be attributed to the efficient mixing of helium and ⁵⁶Ni compared with the 1D case. The IN reaction above does not contribute in the 1D case. The NN reaction in Equation (49) becomes effective in contrast to the one-zone case. The contributing destruction processes of SiO are the NN reaction in Equation (66) and the CE reaction in Equation (68), the ionization by Compton electrons (CM reaction), and the destruction by UV photons (UV reaction). The contribution of the NN reaction in Equation (66) becomes large compared with the one-zone case. The contributions of the IN reactions with He⁺ in Equations (69) and (70) are minor in the 1D case in contrast to the one-zone case.

After 300 days, the amount of CO slightly increases. On the other hand, the amount of SiO is distinctively recovered to some extent. The contributing formation processes of CO at this phase are the CE reaction in Equation (56), the NN reactions in Equations (47), (50), (57), and the NN reaction below

$\begin{eqnarray}&&{\rm{C}}+\mathrm{NO}\to \mathrm{CO}+{\rm{N}}.\end{eqnarray} \tag{ 73 }$

The NN reactions in Equations (58) and (59) also contribute as the secondaries. The contributions of the CE reaction in Equation (56) and the NN reactions in Equations (50) and (73) become large compared with the one-zone case. The major contributing formation processes of SiO at this phase are the NN reaction in Equations (62), (72), (71), and (49). The contributions of the NN reactions in Equations (49) and (72) become large compared with the one-zone case.

As a summary, compared with the fiducial one-zone model, the less efficient mixing of ⁵⁶Ni, hydrogen, and helium in the 1D model could decrease the contributions of the reactions involved with ionized hydrogen and helium atoms, i.e., the CE reactions in Equations (53) and (68) and the IN reactions in Equations (54), (69), and (70).

4.2.2. The Results including Molecules Other Than CO and SiO with the Angle-averaged 1D Profiles Based on the 3D Models; the Comparison between the Binary Merger and Single-star Progenitor Models

Having determined the fiducial values for the parameters (f_h = 5 × 10⁻³, f_d = 1.0, and t_s = 500 days), here, the results including molecules other than CO and SiO with the angle-averaged 1D profiles based on the 3D explosion models (the models b18.3-high and n16.3-high; Ono et al. 2020) with the binary merger (b18.3) and the single-star (n16.3) progenitor models, i.e., the results of the models b18.3-mean and n16.3-mean, are described in Sections 4.2.2.1 and 4.2.2.2, respectively.

4.2.2.1. The Results of the Model b18.3-mean

In the top panels in Figure 20, the results of the model b18.3-mean are shown. The qualitative features of the evolution of the gas temperatures (the right panel) are similar to the case of the parameters, f_h = 10⁻², f_d = 1.0, and t_s = 500 days (see the left panel in Figure 15), although some differences can be recognized due to the smaller f_h value; the inner particles are preferentially heated by the decay of ⁵⁶Ni making peaks around 200 days in the gas temperatures; some of the particles whose initial positions are around 1.5 × 10¹³ cm and a few inner particles experience both the heating and the cooling by CO rovibrational transitions. The highest peak temperature around 200 days (<10⁵ K) is lower than that of the case of f_h = 10⁻² (>10⁵ K) as expected. The time evolution of the fluxes of CO vibrational bands was already described in Section 4.2.1.2 (see the lower right panel in Figure 14; for the comparison in Section 4.2.3, the same panel is shown in the top panel in Figure 21).

**Figure 21.** Time evolution of the fluxes for CO vibrational bands, Δv = 1 (fundamental), Δv = 2 (first overtone), ..., Δv = 6, compared with the observed light curves for Δv = 1 (BD93; W+93; Bouchet & Danziger 1993; Wooden et al. 1993, respectively), and Δv = 2 (M+89, 93; BD93; Meikle et al. 1989, 1993; Bouchet & Danziger 1993, respectively), for the models b18.3-mean (top), b18.3-sphel (middle), and b18.3-sphel-pure (bottom). The top panel is the same as the lower right panel in Figure 14.
Download figure:
Standard image High-resolution image

In the left panel, the time evolution of the total amounts of the molecules including those other than CO and SiO is shown.²⁵ The time evolution of the amounts of CO and SiO and the important chemical reactions for their formation and destruction were already described in Sections 4.2.1.2 (see the lower left panel in Figure 14) and 4.2.1.4, respectively. As with CO and SiO, the molecules other than CO and SiO have qualitatively similar trends in the time evolution, i.e., the amounts increase until about 100 days; after that, those decrease to some extent; after a few hundred days, some of them increase again. Hereafter, the time evolution of the amounts and primary formation and destruction reactions are described for some of the molecules other than CO and SiO focusing on epochs when the changes in the amounts are remarkable; formation (destruction) reactions are focused during the formation (destruction) dominant phases.

At the end of the calculation, the molecules whose amounts are greater than 10⁻³ M_⊙ are CO, SiO, H₂, N₂, MgO, and FeO; there are some gaps in the amounts between the molecules above and the others. The sequence of the reactions for H₂ formation is the REA reaction, H + e⁻ → H⁻ + γ; the AD reaction, H⁻ + H → H₂ + e⁻ (as mentioned in Section 4.1.2). The destruction processes of H₂ after 60 days are basically the UV reaction, the ionization by Compton electrons (one of CM reactions), and the NN reaction of oxygen with H₂ (products: OH and helium). The primary formation processes of N₂ are the NN reaction of nitrogen with NO before about 200 days and the NN reaction of nitrogen with CN after 200 days. The destruction processes of N₂ after about 100 days are the CM reactions (ionization, dissociation, dissociative ionization) and the UV reaction. The formation and destruction processes of MgO and FeO are similar to each other. The formation process of MgO (FeO) before about 200 days is the NN reaction between magnesium (iron) and O₂, and the ones after that are three-body (3B) reactions with hydrogen and helium, e.g., Mg (Fe) + O + H → MgO (FeO) + H. The major destruction processes are the CM and UV reactions.

Among the molecules other than those mentioned above, the amounts of OH, NO, SO, O₂, and SiS become greater than 10⁻⁴ M_⊙ by 200 days, and those remain greater than 10⁻⁶ M_⊙ after the destruction dominant phase around a few hundred days. The primary formation reactions of OH are the RA reaction between oxygen and hydrogen atoms before 50 days and the NN reaction between oxygen and H₂ after 50 days. The destruction process of OH after 60 days is the NN reaction of OH with silicon (products: SiO + H). NO is primarily formed by the NN reactions of nitrogen with CO and OH before 200 days and the NN reaction of oxygen with NH after 200 days. The destruction processes of NO after 100 days are the NN reactions of NO with silicon (products: SiO and nitrogen) and with carbon (products: CO and nitrogen). The major formation processes of SO before 100 days are the RA reaction of sulfur and oxygen atoms and the NN reaction between oxygen and S₂. The formation process after 100 days is the NN reaction of sulfur with OH. The primary destruction reactions of SO after 100 days are the NN reactions of SO with nitrogen (products: NO and sulfur) and with carbon (products: CO and sulfur; CS and oxygen). The primary formation reactions of O₂ before 100 days are the RA reaction between two oxygen atoms and the NN reaction of oxygen and OH. The destruction processes of O₂ after 100 days are the NN reactions of O₂ with silicon (products: SiO and oxygen) and with carbon (products: CO and oxygen). SiS is primarily formed by the RA reaction between silicon and sulfur atoms (before about 100 days and after 200 days) and the NN reaction between silicon and S₂ (from 100 to 200 days). The destruction processes of SiS after 100 days are the NN reactions of SiS with nitrogen (products: SiN and sulfur) and with sulfur (products: S₂ and silicon). The TF reaction of SiS with hydrogen also contributes to the destruction.

It is noted that the products of some of the destructive NN reactions of the molecules above are CO and/or SiO. Therefore, such molecules are somehow converted to CO and/or SiO through the NN reactions.

The amount of S₂ becomes greater than 10⁻⁴ M_⊙ around 100 days, but the one becomes less than 10⁻⁶ M_⊙ due to destruction processes by 200 days. The primary formation processes of S₂ before 100 days are the NN reactions of sulfur with SiS and with FeS. The destruction processes of S₂ after 100 days are the NN reactions of S₂ with oxygen (products: SO and sulfur) and with carbon (products: CS and sulfur).

The amounts of the molecules, CS, CN, and NH, increase after about 100 days, and those become greater than 10⁻⁵ M_⊙ by the end of the calculation. The primary formation processes of CS after 100 days are the NN reaction of carbon with SO and the NN reaction of sulfur with CN (only after 600 days). The destruction process of CS after 130 days is the NN reaction of CS with oxygen (products: CO and sulfur). CN is primarily formed by the NN reactions of carbon with NO and with NH. The primary formation process of NH is the NN reaction of nitrogen with H₂.

Other molecules that are not mentioned above, CH, C₂, SiH, SiN, SiC, and Fe₂ except for MgS and FeS, are minor in the amounts throughout the evolution.

The amounts of MgS and FeS are qualitatively similar; the amount of MgS (FeS) once becomes greater than 10⁻⁵ (10⁻⁴) M_⊙ around 100 days, although later those become less than 10⁻⁶ (10⁻⁵) M_⊙. The primary formation process of MgS (FeS) is the NN reaction of magnesium (iron) with SO; the primary destruction process is the NN reaction of MgS (FeS) with sulfur.

4.2.2.2. The Results of the Model n16.3-mean

In this section, the results of the model n16.3-mean are presented by comparing with those of the model b18.3-mean mainly focusing on the differences. In the top panels in Figure 22, the results of the model n16.3-mean are shown.

**Figure 22.** Same as Figure 20 but for the models n16.3-mean (top panels), n16.3-sphel (middle panels), and n16.3-sphel-pure (bottom panels).
Download figure:
Standard image High-resolution image

As can be seen in the right panel, the gas temperatures of outer particles (≳5 × 10¹³ cm in the bottom panels in Figure 3) are overall higher than those of the model b18.3-mean (top right panel in Figure 20). Then, the timing when the gas temperatures of the outer particles go down to ∼10⁴ K, i.e., when the molecule formation starts, is delayed compared with one for the model b18.3-mean; actually, the molecule formation starts after 10 days (before 10 days for the case of b18.3-mean) as seen in the left panel. Overall, other qualitative features of the gas temperature evolution are similar to the case of the b18.3-mean; due to the gas heating by the decay of ⁵⁶Ni, in particular, the inner particles are heated up, and the gas temperatures peak at about 200 days; some particles initially positioned at (1.5–2.5) × 10¹³ cm go through the cooling via CO rovibrational transitions after about 60 days, and some inner particles also cooled at later phases (≲1000 days).

The fluxes of CO vibrational bands (the top panel in Figure 23) peak before 100 days, and the calculation fails to reproduce the peaks of the observed fundamental and first overtone bands same as other results shown above. The peak flux levels are similar to those of the b18.3-mean case but slightly lower around the observed peaks. The increases in the fluxes at later phases (after 300 days for the fundamental and first overtone bands) are attributed to the cooling of inner particles.

Related to the formation of molecules, the amounts of initial seed atoms are different from those of the model b18.3-mean as seen in the top left panel in Figure 22. Apparently, the amount of oxygen atoms (∼1 M_⊙) is higher than that of the model b18.3-mean (see the top left panel in Figure 20). On the other hand, the amounts of carbon and silicon atoms are lower than those of the model b18.3-mean. Then, the abundance ratios of oxygen atoms to carbon and silicon are higher than the case of the model b18.3-mean. Additionally, in the model n16.3-mean, the amount of carbon atoms is higher than that of silicon in contrast to the model b18.3-mean. The initial amount of sulfur is also different from the case of b18.3-mean. Such differences in the initial abundance ratios of seed atoms could potentially affect the formation of molecules.

The qualitative trends of the evolution of the amount of CO are similar to the case of b18.3-mean except for the initial rise. On the other hand, the amount of SiO is rather different between the two models. In the model n16.3-mean, the amount of SiO at around 40–200 days is lower compared with the model b18.3-mean roughly by a factor of 5. The amount of SiO at about 200–600 days is rather lower (at most 1 order of magnitude) than that of the model b18.3-mean. The final amount of SiO is also lower than that of the model b18.3-mean approximately by a factor of a few.

The contributing formation and destruction processes of SiO at about 40–200 days are not so different between the models b18.3-mean and n16.3-mean, although the order of the significance of those processes is a bit different between the two models depending on time. After 200 days, the primary destruction process of SiO is the CE reaction in Equation (68) for both the models b18.3-mean and n16.3-mean. The secondary destruction processes are the NN reaction in Equation (66), the ionization by Compton electrons (one of CM reactions), and the UV reaction; in the model n16.3-mean, among the secondary processes, the latter processes dominate the NN reaction (in the model b18.3-mean, the NN reaction dominates the processes related to Compton electrons). The difference in the amount of SiO between the two models may be attributed to the different initial abundance ratios of the seed atoms, in particular, the ratio of carbon to silicon, and/or the differences in gas densities and temperatures.

As for other molecules, for example, probably due to the differences in the initial abundance ratios of sulfur to oxygen, the model n16.3-mean leads to a higher amount of O₂ and a lower amount of SO than those of the model b18.3-mean. Similarly, in the model n16.3-mean, the amount of SiS is apparently lower than that of the model b18.3-mean due to the lower silicon and sulfur abundances.

Here, the primary formation processes of O₂ is the RA reaction between two oxygen atoms and the NN reaction of oxygen with OH as the same as the model b18.3-mean. SO is mainly formed by the RA reaction between sulfur and oxygen atoms and the NN reaction of sulfur with O₂. It is noted that, in the model b18.3-mean, the NN reaction of oxygen with S₂ is dominant instead of the NN reaction above. The discrepancy is probably due to the high initial abundance ratio of oxygen to sulfur in the model n16.3-mean. SiS is mainly produced by the RA reaction between silicon and sulfur atoms. The NN reaction of silicon with S₂, which is effective in the model b18.3-mean, however, does not contribute to the formation of SiS in the model n16.3-mean. This can also be attributed to the low initial abundance ratio of sulfur to the others.

After about 130 days, for both two models, O₂, SO, and SiS are significantly destructed (in the model n16.3-mean, SiS does not appear in the top left panel in Figure 22 before 2000 days). The main destruction processes of O₂ at this phase are the NN reactions of O₂ with silicon and with carbon. The CE reaction between H⁺ and O₂ also contributes to some extent in contrast to the model b18.3-mean. SO is primarily destructed by the NN reactions of carbon with SO (two cases of products CS and oxygen; CO and sulfur), the NN reaction of nitrogen with SO, and the CE reaction between H⁺ and SO. The contribution of the CE reaction above is limited in the model b18.3-mean. The destruction processes of SiS are a bit complicated. SiS is destructed mainly by the NN reactions of SiS with nitrogen (products: SiN and sulfur) and with hydrogen (products: SiH and sulfur). The thermal fragmentation reaction of SiS with helium also primarily contributes to the destruction. The CE reaction between H⁺ and SiS contributes to some extent in contrast to the model b18.3-mean. Among those destruction processes, the significance varies depending on time and the models. It is noted that, since the CE reactions mentioned above involve ions, the mixing of ⁵⁶Ni, which influences the ionization by Compton electrons, may affect the significance. On the other hand, the primary destruction processes of O₂, SO, and SiS are not always ion-related reactions like the CE reactions above but are rather NN reactions; the amounts would be determined by the complicated balance between the NN reactions depending on the gas temperatures.

Among the other molecules, the trends of H₂, OH, MgO, and FeO are qualitatively similar to those of the model b18.3-mean. The amount of MgO is, however, larger than that of FeO in the model n16.3-mean in contrast to the model b18.3-mean. The contributing formation and destruction processes of MgO and FeO are consistent with the model b18.3-mean (see Section 4.2.2.1). This feature is probably attributed to the initial high abundance of magnesium in the model n16.3-mean (see, e.g., the left panels in Figures 3 and 4).

A distinct difference between the models b18.3-mean and n16.3-mean is the amount of N₂; the amount of N₂ in the model n16.3-mean is significantly lower than that of the model b18.3-mean. The contributing formation and destruction reactions are the same as the model b18.3-mean. The feature above is likely attributed to the fact that the nitrogen abundance, in particular, at the outer layers in the progenitor model b18.3 is higher than that in n16.3. As mentioned, since the progenitor model b18.3 is based on a binary merger, the CNO cycle additionally triggered during the merger process increases the CNO cycle-processed materials including nitrogen (¹⁴N) into the envelope.

As a summary, it is difficult to clearly figure out what is the primary reason for the differences in the amounts of those molecules between the two models. Though, the initial abundance ratios of seed atoms, ⁵⁶Ni, hydrogen, and helium, which reflect the matter mixing, likely affect the results as presented above.

4.2.3. Comparison with the Spherical Cases

In this section, in order to see the impact of the effective matter mixing, by comparing the models b18.3-mean and n16.3-mean, the results based on the initial profiles for spherical cases (spherical explosion and purely spherical cases), i.e., the results for the models b18.3-sphel, b18.3-sphel-pure, n16.3-sphel, and n16.3-sphel-pure, are presented; the former two and the latter two are described in Sections 4.2.3.1 and 4.2.3.2, respectively.

4.2.3.1. The Cases with the Binary Merger Progenitor Model

In the middle panels in Figure 20, the results of the model b18.3-sphel (the spherical explosion case) are shown. The inner particles (initially less than about 1 × 10¹³ cm) are heated by the decay of ⁵⁶Ni (the right panel) since ⁵⁶Ni is initially distributed only inside ∼1 × 10¹³ cm in contrast to the model b18.3-mean, in which the particles initially inside ∼3 × 10¹³ cm are a matter of heating. As a result, the peak temperatures of the inner particles at around 200 days become higher than those of the model b18.3-mean due to the high local ⁵⁶Ni abundance. Some particles initially at (1–2) × 10¹³ cm and a few inner particles go through the cooling by CO rovibrational transitions.

In the middle panel in Figure 21, the time evolution of the fluxes of CO vibrational bands is shown. The calculated fluxes of CO vibrational bands are qualitatively similar to those of the b18.3-mean (the top panel), but quantitatively, the flux levels are slightly lower than those of the b18.3-mean.

The qualitative features of the evolution of the amounts of CO and SiO (the middle left panel in Figure 20) are roughly consistent with that of the model b18.3-mean; some differences, however, can also be recognized. The amounts of CO and SiO are higher than those of the model b18.3-mean throughout the evolution.

The initial rise of the amounts of molecules is steeper than that of the model b18.3-mean. At approximately 50 days, the amounts of both CO and SiO exceed 0.1 M_⊙. Then, later, the amounts are reduced by destruction processes; the destruction of some of the molecules, in particular SiO, is, however, milder than that of the model b18.3-mean. This is because ⁵⁶Ni is concentrated only in inner regions compared with the model b18.3-mean. In the inner particles heated by the decay of ⁵⁶Ni, the molecule formation is delayed. On the other hand, in the outer particles, with the relatively low local ⁵⁶Ni abundance, molecules start to form earlier in a higher-density environment, and the destruction processes caused by the decay of ⁵⁶Ni, i.e., the CE reaction in Equation (68) and the CM reactions (ionization, dissociation, and dissociative ionization by Compton electrons), are less effective compared with the model b18.3-mean. In total, less efficient mixing of ⁵⁶Ni in the model b18.3-sphel results in earlier effective formation and less effective destruction of molecules compared with the model b18.3-mean. Actually, it is confirmed that, in the model b18.3-sphel, the primary destruction processes of SiO after about 100 days are the NN reactions in Equation (66), and the CE reaction in Equation (68) is the secondary in contrast to the model b18.3-mean (in the model b18.3-mean, the CE reaction is the primary), which supports the statement above.

Additionally, the evolution of the amounts of other molecules such as SO and O₂ is distinctively different. The amounts of SO and O₂ at around 50–130 days are rather higher (at most, more than 1 order of magnitude) than those of the model b18.3-mean, and the final amounts are also affected. At about 50–100 days (SO and O₂ increasing phase), the primary formation processes of SO are the NN reactions of sulfur with O₂ and with OH. The primary formation process of O₂ is the NN reaction of oxygen with OH. For both SO and O₂, the major formation processes at this phase are consistent with the model b18.3-mean.

After 100 days, in the model b18.3-sphel, the amounts of SO and O₂ start to decrease in contrast to the model b18.3-mean. The major destruction reactions of SO are the NN reactions of SO with carbon (products: CO and sulfur; CS and oxygen). The primary destruction process of O₂ is the NN reaction of O₂ with silicon. Therefore, the destruction processes of SO and O₂ are also consistent with the model b18.3-mean. The abundances of SO and O₂ are probably determined by the balance of the NN reactions depending on the gas temperatures.

The behavior of H₂ is different from that of b18.3-mean. In the model b18.3-sphel, the early formation of H₂ is limited, and the ratios of the amount of H₂ to CO and SiO are apparently smaller than those in the model b18.3-mean. Since the primary formation process is the sequence, the REA reaction, H + e⁻ → H⁻ + γ; the AD reaction, H⁻ + H → H₂ + e⁻, the smaller amount of H₂ in the model b18.3-sphel may be attributed to the less effective ionization due to less effective mixing of ⁵⁶Ni.

Among the other molecules, N₂, MgO, and FeO have a different feature from that of the model b18.3-mean; after the formation of those molecules before 100 days, there is no distinct destruction of those species. Since the destruction processes of those are basically the CM reactions, the less effective destruction in the model b18.3-sphel may be attributed to the less effective mixing of ⁵⁶Ni, again.

In the bottom panels in Figure 20, the results of the model b18.3-sphel-pure (the purely spherical case) are shown. First, as seen in the right panel, only several inner particles are heated by the decay of ⁵⁶Ni, since the ⁵⁶Ni is initially distributed at the innermost region (≲0.5 × 10¹³ cm; see the top left panel in Figure 5). The peak temperatures of the heated particles are significantly higher than those of the models b18.3-mean and b18.3-sphel due to the high local ⁵⁶Ni abundance. Only several inner particles are affected by the cooling due to the CO rovibrational transitions. The gas temperature evolution of the other particles not affected by the heating and cooling follows the broken power law.

In the bottom panel in Figure 21, the fluxes of CO vibrational bands are shown. The fluxes at early phases (before 300 days) are distinctively lower than those of the models b18.3-mean and b18.3-sphel. It may be understood that, due to the inefficient ionization by Compton electrons (because of the less effective mixing of ⁵⁶Ni), the number density of (thermal) electrons for exciting the vibrational levels of CO is lower than that of efficient ionization cases (efficient mixing cases). Actually, this tendency, i.e., less efficient ionization results in lower CO vibrational fluxes, is also seen in between the models b18.3-mean and b18.3-sphel; the fluxes of CO vibrational transitions in the model b18.3-mean (the top panel in Figure 21), in which the mixing of ⁵⁶Ni and ionization by Compton electrons are efficient compared with the model b18.3-sphel, are overall higher than that of the model b18.3-sphel (the middle panel in the same figure). The sudden rise in the fluxes at around 300 days is apparently attributed to the particle that has a rapid cooling (see the bottom right panel in Figure 20); this feature would be artificial due to the insufficient number of particles around the regions where the particle with the rapid cooling is (the gas temperatures have huge gaps between adjacent particles). Although a higher resolution (number) of tracer particles may result in smoother light curves of CO vibrational bands and the radial distributions of molecules, it should not affect the total amounts of molecules much.

The evolution of the amounts of molecules has noticeable differences from the models b18.3-mean and b18.3-sphel. CO and SiO molecules are sharply formed after about 40 days when the gas temperatures of most of the particles go down to ∼10⁴ K. After the rapid formation of CO and SiO, plateaus of the amounts of them can be recognized. This phase corresponds to the timing when molecules start to form in the tracer particles inside the helium-rich shell (≲0.8 × 10¹³ cm in the top left panel in Figure 5), i.e., their gas temperatures cool down to ∼10⁴ K. The amount of SiO has a plateau with sharp edges compared with that of CO. After SiO enters the plateau phase, the amount of CO continues to increase rapidly for a while. Then, the increase becomes slow. The plateau of SiO with the sharp edges is probably due to the abundance ratios among the seed atoms, carbon, oxygen, and silicon, inside the helium-rich shell; the high abundance ratio of carbon to oxygen (actually, N(C) > N(O) in the helium-rich shell, where N(C) and N(O) denote the number density of the seed carbon and oxygen atoms, respectively) and the low abundance ratios of silicon to carbon and oxygen may inhibit the formation of SiO. After 70 days, the amounts of CO and SiO again start to increase. This phase corresponds to the timing when molecules start to form in the tracer particles inside oxygen- and silicon-rich shells (≲0.5 × 10¹³ cm in the top left panel in Figure 5). The high abundance ratios of oxygen and silicon allow SiO to form and contribute to increasing both the amounts of CO and SiO. After 100 days, the amounts of CO and SiO, in particular, SiO, decrease. The primary destruction processes of CO are the NN reaction in Equation (49) (SiO formation reaction) and the CM and UV reactions. SiO is primarily destructed by the CM and UV reactions. For the destruction by CM and UV reactions, several particles affected by the decay of ⁵⁶Ni (deviated from the broken power-low evolution) except for ones whose peak temperatures after 100 days are greater than 10⁵ K contribute to the destruction. After 500 days, SiO distinctly increases again, and the amount becomes comparable with that of CO. The contributing formation processes of SiO are the NN reactions in Equations (49) and (71).

It is noted that, in the innermost particles, due to the strong heating with the high local ⁵⁶Ni abundance, the gas temperatures go down to 10⁴ K only at rather late phases after a few thousand days when the gas densities are very low. Then, the contributions of the innermost particles to the amounts of molecules may be relatively small.

The properties of the formation of other molecules are also different from the models b18.3-mean and b18.3-sphel. In the models b18.3-mean and b18.3-sphel, the formation of SO and O₂ is recognized at about 40–200 days; in the model b18.3-sphel-pure, the two molecules are, however, not distinctly formed at least before 100 days. Instead, CS and SiC are formed just after 40 days, and the amount of CS changes remarkably at around 70 days. CS is mainly formed by the RA reaction between carbon and sulfur atoms, the NN reaction of sulfur with C₂, and the NN reaction of carbon with SO before 70 days. SiC is primarily formed by the NN reaction of carbon with SiO, the RA reaction between silicon and carbon atoms, and the NN reaction of silicon with C₂ before 70 days. The major destruction processes of CS after 70 days are the NN reaction of hydrogen with CS (products: CH and sulfur) and the NN reaction of oxygen with CS (products: CO and sulfur). SiC is mainly destructed by the NN reactions of oxygen with SiC (products: SiO and carbon; CO and silicon) after 70 days, although the decrease of the amount is not distinctively recognized. It is noted that, in both of the formation processes of CS and SiC, C₂ plays a major role. In the model b18.3-sphel-pure, C₂ is significantly formed after 40 days in contrast to the models b18.3-mean and b18.3-sphel; the amount exceeds 10⁻² M_⊙. C₂ is mainly formed by the RA reaction between two carbon atoms and the NN reactions of carbon with CN and CH. Similarly, for example, SiN is also markedly formed in contrast to the models b18.3-mean and b18.3-sphel, where SiN is primarily formed by the NN reaction of nitrogen with SiO.

On the other hand, as observed in the model b18.3-sphel, the formation of H₂ is further limited compared with that of the model b18.3-mean due to the less effective mixing of ⁵⁶Ni. Additionally, the formation of OH is also significantly reduced compared with the models b18.3-mean and b18.3-sphel. The primary formation processes of OH before 100 days are the RA reaction between oxygen and hydrogen atoms and the NN reaction of oxygen with H₂. Therefore, the formation of OH is also affected by the mixing of ⁵⁶Ni through the latter reaction above.

As a summary, the distinct features found in the model b18.3-sphel-pure are probably attributed to the initial stratified profiles of the seed atoms and the abundance ratios of each layer (see the top left panel in Figure 3) owing to the purely spherical hydrodynamical evolution and the less effective mixing of ⁵⁶Ni into outer layers.

4.2.3.2. The Cases with the Single-star Progenitor Model

In this subsection, the results of spherical cases based on the single-star progenitor model (n16.3), i.e., ones for the models n16.3-sphel and n16.3-sphel-pure, are presented mainly focusing on the differences from those based on the binary merger progenitor model (b18.3).

In the middle and bottom panels in Figure 22, the results of the models n16.3-sphel and n16.3-sphel-pure are shown, respectively. The qualitative features of the gas temperatures are similar to the models b18.3-sphel and b18.3-sphel-pure. Some recognizable differences are as follows. The temperatures of the innermost particles in the model n16.3-sphel are higher than those in the model b18.3-sphel. In the models n16.3-sphel and n16.3-sphel-pure, the number of significantly heated particles whose peak temperatures at around 200 days greater than 10⁵ K is higher than those of the corresponding counterpart models b18.3-sphel and b18.3-sphel-pure. This feature may partly be attributed to the differences in the initial distribution of ⁵⁶Ni at the innermost region between the spherical models with the progenitor models b18.3 and n16.3 (see the left panels in Figures 4 and 5).

In the middle and bottom panels in Figure 23, the fluxes of CO vibrational bands for the models n16.3-sphel and n16.3-sphel-pure are shown, respectively. As seen in the cases of the binary merger progenitor model (b18.3), the more spherical the evolution is, the lower the fluxes of CO vibrational bands are; the fluxes of CO vibrational bands in the models n16.3-mean, n16.3-sphel, and n16.3-sphel-pure decrease in this order. Again, this may be attributed to the inefficient ionization by Compton electrons produced by the decay of ⁵⁶Ni in the spherical cases. In the model n16.3-sphel-pure, there is no distinct increase of fluxes at later phases (after 300 days) in contrast to the model b18.3-sphel-pure. This is because no particle goes through rapid cooling at such late phases.

The evolution of the amounts of molecules in the models n16.3-sphel and n16.3-sphel-pure is different from the corresponding counterpart models based on the binary merger progenitor model (b18.3) in several points.

As seen in the model n16.3-mean, in both models n16.3-sphel and n16.3-sphel-pure, molecules start to form a bit later time (after 10 days) compared with the models based on the progenitor model b18.3 by reflecting the fact that the gas temperatures of outer particles are relatively higher than those of the corresponding counterpart models with the progenitor model b18.3. In the models n16.3-sphel and n16.3-sphel-pure, the amounts of CO and SiO are more stable after 40 days compared with the models b18.3-sphel and b18.3-sphel-pure. In the model b18.3-sphel-pure, in particular, the amount of SiO varies in time partly due to the stratified distributions of the seed atoms. The more stable feature in the amounts of CO and SiO in the spherical models based on the single-star progenitor model (n16.3) might partly be attributed to the less contribution from the significantly heated particles, in which molecule formation starts only at the later phases in the low-density environment due to the heating, (the number of those particles is larger than the counterpart models with the binary merger model as mentioned) and the relatively larger contribution from the particles that are not significantly affected by the decay of ⁵⁶Ni, in which the destruction processes induced by the decay of ⁵⁶Ni are not so significant.

As mentioned in Section 4.2.2.2, the initial abundance ratios of the seed atoms in the model n16.3-mean are different from those of the counterpart model b18.3-mean; the abundance ratios of oxygen to silicon and carbon are higher than that of the model b18.3-mean, and the amount of oxygen in n16.3-mean is initially as high as ∼1 M_⊙. This feature does not change in the spherical cases, i.e., the models n16.3-sphel and n16.3-sphel-pure. By reflecting this feature, in the models n16.3-sphel and n16.3-sphel-pure, the amount of O₂ is significantly higher than those of the models based on the binary merger progenitor (b18.3) and the model n16.3-mean. In the model n16.3-sphel, the amount of O₂ is as high as 0.6 M_⊙ and is higher than that of CO at approximately 40–200 days. Later, it becomes almost comparable to that of CO after the destruction from 100 days to a few hundred days. O₂ is primarily formed by the NN reaction of oxygen with OH and the NN reaction between two OH molecules. The contributing destruction processes after 100 days are the NN reactions of O₂ with silicon (products: SiO and oxygen), with sulfur (products: SO and oxygen), and with carbon (products: CO and oxygen); the UV reaction; and the ionization by Compton electrons (one of the CM reactions).

In the model n16.3-sphel-pure, the amount of O₂ is the highest among the molecules after 30 days, and it is as high as ∼1 M_⊙ at the last phase. The primary formation process of O₂ in the model n16.3-sphel-pure is the AD reaction, O⁻ + O → O₂ + e⁻, where O⁻ can be formed by the REA reaction, O + e⁻ → O⁻ + γ. The NN reaction of oxygen with MgO also contributes to the increase of O₂ after 100 days. Therefore, the formation processes of O₂ in the n16.3-sphel-pure are peculiar compared with the other models. In the model n16.3-sphel-pure, the seed oxygen is initially concentrated in the oxygen-rich shell (see the lower left panel in Figure 5), and the abundance ratios of oxygen to others, in particular, the ratio of oxygen to carbon, are higher than those of the other models. Additionally, around the inner edge of the oxygen-rich shell, ⁵⁶Ni is overlapped, where the number density of electrons is enhanced by the ionization by Compton electrons. Therefore, in such regions, the former formation sequence above would preferentially work compared with the other models.

As seen between the models n16.3-mean and b18.3-mean, the amount of MgO is higher than that of FeO in the models n16.3-sphel and n16.3-sphel-pure. On the other hand, the amount of FeO dominates that of MgO in the models based on the binary merger progenitor (b18.3) model at least in the last phase. As mentioned (in Section 4.2.2.2), this feature is probably attributed to the higher initial magnesium abundance in the models based on the single-star progenitor model (n16.3). Another distinct feature in the models n16.3-sphel and n16.3-sphel-pure is that the amount of MgO is significantly higher than that of FeO compared with the model n16.3-mean. The primary formation processes of MgO (FeO) are consistent with the models based on the binary merger model and the model n16.3-mean, i.e., the NN reaction of magnesium (iron) with O₂ and the 3B reaction among magnesium (iron), oxygen, and hydrogen. The feature of the high MgO abundance in the models n16.3-sphel and n16.3-sphel-pure is probably attributed to the higher significance of the former reaction compared with the other models due to the high O₂ abundance.

N₂ is also formed in the models n16.3-sphel and n16.3-sphel-pure, and the amount is smaller than that in the counterpart models b18.3-sphel and b18.3-sphel-pure. This is because of the initial smaller abundance ratios of the seed nitrogen atoms to the others as mentioned.

4.2.4. Differences among the Specified Directions

In Sections 4.2.1–4.2.3, the impact of matter mixing on the molecule formation was described by using the angle-averaged 1D profiles. In reality, radial profiles could, however, vary depending on the direction. Such angle-averaged profiles may overestimate the effects of matter mixing compared with the case of the direct application to the 3D models. Additionally, in the 3D models b18.3-high and n16.3-high (Ono et al. 2020), the explosion is globally asymmetric (bipolar-like) as mentioned. It is worth seeing the dependence of the molecule formation on the specified direction (see Kozyreva et al. 2021, for discussion on the light-curve dependence on different directions in the highly asymmetric explosion). In this section, the molecule formation results for the specified directions, +Z, −Z, and +Y, are presented. As a reminder, +Z- and −Z-axes are pointed to the stronger and weaker bipolar explosion directions, respectively; the +Y-axis is perpendicular to the bipolar explosion axis. It should be noted that the initial energy injections (depositions) in the −Z- and +Y-directions are lower than that in the +Z-direction by factors of 2.2 and 3.7 × 10², respectively (Ono et al. 2020). The different energy depositions significantly affect the amounts and the distribution (radial extension) of ⁵⁶Ni as can be seen in Figures 7 and 8; in the +Z- and −Z-directions, the amounts and the radial extension of ⁵⁶Ni are much higher than those in the +Y-direction, which could affect the formation of molecules.

In Sections 4.2.4.1 and 4.2.4.2, the time evolution of physical quantities for the models b18.3-zp, b18.3-zn, b18.3-yp, n16.3-zp, n16.3-zn, and n16.3-yp are presented (the former three and the latter three models are in Sections 4.2.4.1 and 4.2.4.2, respectively). Radial distributions of the molecules for the former three and the latter three models are described in Sections 4.2.4.3 and 4.2.4.4, respectively.

4.2.4.1. The Time Evolution of Quantities in the Models with the Binary Merger Progenitor Model

Figure 24 shows the results of the models b18.3-zp, b18.3-zn, and b18.3-yp. First, from the evolution of gas temperatures (right panels), the evolution in the +Y-direction (bottom one) is apparently different from the ones in the +Z- and −Z-directions. This feature is clearly attributed to the differences in the deposited energies during the explosion, i.e., in the +Y-direction, the deposited energy is much lower (by a factor of a few hundred) than those in the other two directions (bipolar-like explosion axis); the amount of ⁵⁶Ni synthesized in the +Y-direction is significantly lower than those in the other two directions. Hence, in the +Y-direction, only several inner particles are affected by the decay of ⁵⁶Ni by reflecting the fact that ⁵⁶Ni is initially located only in the innermost region (less than 2 × 10¹² cm; see the distribution of Fe corresponding to ⁵⁶Ni in the bottom left panel in Figure 7). On the other hand, ⁵⁶Ni is more extensively distributed (up to ∼2 × 10¹³ cm) in the +Z- and −Z-directions. Then, as seen in the top and middle right panels in Figure 24, the numbers of heated particles are apparently higher than that in the +Y-direction; the particles initially at less than 2 × 10¹³ cm are distinctively heated, and the temperatures peak at around 200 days. In the +Z- and −Z-directions, some of the particles initially located at approximately (1.5–2.5) × 10¹³ cm go through the cooling by CO rovibrational transitions. In the +Y-direction, the inner particles initially located inside 1.5 × 10¹³ cm are affected by the cooling. It is noted that, in the inner regions where ⁵⁶Ni is initially distributed, helium is also present (see the left panels in Figure 7); a part of helium in such a region may stem from the unburned alpha particles during the explosive nucleosynthesis.

The time evolution of the mass fractions of molecules (not the amounts) is shown in the left panels. For all three directions, CO and SiO rapidly increase around 40 days, but the mass fractions along the +Z- and −Z-directions are rather different from ones along the +Y-direction. In the +Z- and −Z-directions, the mass fractions of both CO and SiO at 50–100 days range over 10⁻³–10⁻². On the other hand, in the +Y-direction, the mass fractions CO and SiO at this phase are as much as 10⁻²–10⁻¹; actually, large fractions of the seed atoms, i.e., carbon, oxygen, and silicon, are consumed in the +Y-direction (see the bottom left panel). In the +Y-direction, the mass fraction of SiO is comparable with that of CO at the final phase. This trend (the mass fractions of CO and SiO in the +Y-directions are higher than those in the other two directions) should be attributed to the different energy depositions during the explosion again, which result in higher initial abundances of the seed atoms in the +Y-direction and higher ⁵⁶Ni abundances in the +Z- and −Z-directions (see Figure 7). In the +Y-direction, molecules start to form at earlier phases with a high density because the inefficient gas heating due to the lower local ⁵⁶Ni abundance makes the gas temperatures about 10⁴ K at earlier phases compared with the other two directions.

After 100 days, the mass fractions decrease due to destructive reactions. In the +Z- and −Z-directions, in particular, SiO is distinctively destructed compared with the +Y-direction. At this phase, CO is mainly destructed by the ionization, dissociation, and dissociative ionization due to Compton electrons (CM reactions) and the destruction by UV photons (UV reaction) in the +Z-and −Z-directions. On the other hand, the primary destruction process of CO is the NN reaction in Equation (49) in the +Y-direction. SiO is mainly destructed by the CE reaction in Equation (68) in the +Z- and −Z-directions; the NN reaction in Equation (66), the ionization and dissociation through CM reactions, and the UV reaction are followed as secondary destruction processes. In the +Z (−Z) -direction, the contributions from the CM and UV reactions are higher (lower) than that in the −Z (+Z) -direction. In the +Y-direction, the primary destruction processes of SiO are the NN reactions in Equations (66) and (67) in contrast to the +Z- and −Z-directions. The differences in the destruction processes between the bipolar-like directions (+Z and −Z) and the +Y-direction are because of the different local ⁵⁶Ni abundances due to the different initial energy depositions, again. The differences in the significance of the SiO destruction processes between +Z- and −Z-directions may partly be attributed to the distribution of ⁵⁶Ni; in the −Z-direction, the innermost oxygen- and silicon-rich region (6 × 10¹² cm) in the initial seed atom distributions (middle left panel in Figure 7), the mass fraction of ⁵⁶Ni is smaller than those of oxygen and silicon in contrast to the +Z-direction.

As for other molecules, from here, several features different among the directions are briefly described. For example, H₂ is also remarkably formed after 40 days in the +Z- and −Z-directions. On the other hand, in the +Y-direction, H₂ is formed at later phases after 100 days. In the +Z-and −Z-directions, H₂ is primarily formed by the AD reaction, H⁻ + H → H₂ + e⁻ (H⁻ can be formed by the REA reaction, H + e⁻ → H⁻ + γ) after 40 days. Therefore, the ionization by Compton electrons plays a role in the formation of H₂. This means that the higher local ⁵⁶Ni abundances that stem from the higher initial energy depositions in the +Z- and −Z-directions could also be important for the formation of H₂. Actually, in the +Y-direction, the main formation process of H₂ after 100 days is not the AD reaction above but the NN reaction of hydrogen with SiH. In the −Z-direction, the mass fraction of H₂ is higher than that of CO after approximately 40 days in contrast to the +Z-direction, which may be attributed to the fact that the mass fraction of hydrogen in the −Z-direction is initially higher than that in the +Z-direction in the outer layers (greater than 1 × 10¹³ cm) as seen in Figure 7.

SO and O₂ are also markedly produced at around 40–200 days in the +Y-direction compared with the other two directions. However, after 200 days, the mass fractions are not so dramatically different among the directions. In the +Z- and −Z-directions, the primary formation processes of SO before 200 days are the RA reaction between sulfur and oxygen atoms (before 40–50 days), the NN reactions of oxygen with S₂, and the NN reaction of sulfur with OH. In the +Y-direction, the NN reaction of sulfur with O₂ also contributes as one of the primary formation reactions. The primary formation processes of O₂ are the RA reaction between two oxygen atoms (before 40–50 days) and the NN reaction of oxygen with OH in all three directions. In the +Y-direction, the NN reaction of oxygen with SiO also contributes temporarily. The contributing destruction processes of SO at the destruction dominant phase are the NN reaction of carbon with SO (products: oxygen and CS; CO and sulfur), the NN reaction of iron with SO (products: oxygen and FeS), and the NN reaction of magnesium with SO (products: oxygen and MgS). In the −Z (+Y) -direction, the contribution of the NN reaction of magnesium (iron) is limited compared with the other directions. The primary destruction process of O₂ is the NN reaction of silicon with O₂ (products: oxygen and SiO) in all three directions. Although the destruction processes of SO and O₂ are depending on the direction, basically, those are all NN reactions, i.e., the destructions are not at least directly related to the decay of ⁵⁶Ni.

As seen in the final phase, FeO and MgO are also formed to some extent. In the +Z- and −Z-directions, overall, FeO dominates MgO throughout the evolution. On the other hand, in the +Y-direction, MgO dominates FeO after approximately 40 days. This feature may be attributed to the low iron (decayed from ⁵⁶Ni; ⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe) abundance in the +Y-direction. FeO is primarily formed by the NN reaction of iron with O₂ or the 3B reaction, Fe + O + H → FeO + H, depending on the time before about 100 days. After 100 days, FeO is primarily formed by the 3B reaction above. Similarly, MgO is mainly produced by the NN reaction of magnesium with O₂ or the 3B reaction, Mg + O + H → MgO + H before 100 days. After 100 days, MgO is primarily formed by the 3B reaction. It is noted that, in the +Z- and −Z-directions, in particular, FeO increases even after a few thousand days. It is interpreted that inner particles heated by the decay of ⁵⁶Ni with high local iron (from the decayed ⁵⁶Ni, again) abundance partake in the molecule formation at about a few thousand days; at this phase, only 3B reactions play a role in the formation of FeO and MgO.

4.2.4.2. The Time Evolution of Quantities in the Models with the Single-star Progenitor Model

In Figure 25, the results of the models based on the single-star progenitor model n16.3, i.e., n16.3-zp, n16.3-zn, and n16.3-yp, are shown. From the evolution of the gas temperatures (left panels), it is recognized that the number of heated particles is the largest in the +Z-direction, and it is the smallest in the +Y-direction. In the −Z-direction, the number is in between the former two directions. This feature should be attributed to the different initial energy depositions, again. Actually, as can be seen in Figure 8, in the +Z-, −Z-, and +Y-directions, ⁵⁶Ni is initially distributed extensively in this order. The number of heated particles reflects how extensively ⁵⁶Ni is initially distributed. The temperatures of the heated particles peak at around 200 days. In the +Z-direction, the peak temperatures of some of the heated particles are higher than 10⁶ K, which is higher than the highest peak temperature, e.g., in the model b18.3-zp. As seen in the top left panel in Figure 8, there is a region (at around 5 × 10¹² cm) where the mass fraction of Fe (practically corresponding to ⁵⁶Ni) is as high as ∼0.7; such a high mass fraction of ⁵⁶Ni cannot be seen in the model b18.3-zp. The particles that have peak temperatures higher than 10⁶ K come from such a high local ⁵⁶Ni abundance region. In the +Z-direction, some of the particles initially at about (1.5–2.5) × 10¹³ cm go through the cooling by CO rovibrational transitions. In the −Z-direction, some of the particles initially at (1–2) × 10¹³ cm are affected by the cooling. In the +Y-direction, only several inner particles are affected by the cooling.

As seen in the evolution of the mass fractions of molecules (left panels), the decreases in the values of CO and SiO are limited compared with the models based on the progenitor model b18.3. In contrast to the models based on b18.3, molecules start to form after 10 days. In the +Z- and −Z-directions, CO is only slightly destructed after 200 days, and SiO decreases by a factor of a few. In the +Y-direction, there is no recognizable change in the values of CO and SiO after about 50 days.

In the +Z- and −Z-directions, CO is primarily formed by the RA reaction in Equation (45) before approximately 30–40 days; after that, the NN reaction in Equation (47) primarily contributes until the value enters a plateau (50–60 days). The CE reaction in Equation (56) also contributes to some extent as the secondary reaction. In the +Y-direction, the primary formation processes are the same as the +Z- and −Z-directions except the CE reaction does not contribute. In the +Z- and −Z-directions, SiO is primarily produced by the RA reaction in Equation (61) before approximately 30 days; later, SiO is mainly produced by NN reactions in Equations (62), (71), and (49) until the values have a plateau. In the +Y-direction, the situation is similar, but the contribution from the RA reaction is more limited; the RA reaction forms SiO until 25 days as the primary formation process. The NN reaction in Equation (62) also does not contribute distinctively. In the +Z- and −Z-directions, after approximately 200 days, CO and SiO are destructed to some extent as mentioned above. The primary destruction processes of CO at this phase are the ionization and dissociation by Compton electrons (CM reactions), and the UV reaction is followed as the secondary process. In both the +Z- and −Z-directions, the primary destruction process of SiO after 200 days is the CE reaction in Equation (68); the NN reactions in Equations (66) and (67) also contribute as the secondary processes.

A distinct feature that cannot be seen in the models based on the binary merger progenitor model b18.3 is that O₂ is remarkably formed compared with the models with b18.3. This feature is basically attributed to the high abundance ratios of the seed oxygen atom to the others as mentioned. In all three directions, the mass fraction of O₂ peaks at around 70–80 days. It is noted that the mass fraction of O₂ is higher than those of CO and SiO at the peaks. After the peak, in the +Z (−Z) -direction, the value decreases by roughly 2 orders (1 order) of magnitude. In the +Y-direction, the value decreases only by a factor of a few. In all the three directions, initially, the RA reaction between two oxygen atoms is dominant as the primary formation process until 30–40 days. After that, the NN reactions between oxygen and OH and between two OH molecules overtake the RA reaction; the significance of the two NN reactions varies depending on the direction and time. It is noted that, in the models with the progenitor model b18.3, the dominant NN reaction for the formation of O₂ is only the NN reaction of oxygen with OH at the formation dominant phase. As for the destruction processes after the peak (∼100 days), in the +Z- and −Z-directions, the contributing destruction processes of O₂ are the NN reactions of O₂ with silicon and sulfur and the CE reaction between H⁺ and O₂. In the +Y-direction, only the NN reaction of O₂ with silicon contributes in contrast to the other two directions. The contribution of the CE reaction in the +Z and −Z may be due to the efficient ionization by CM reactions owing to the high ⁵⁶Ni abundance.

As for other molecules, for example, SO qualitatively has a similar trend with O₂. The primary formation process of SO is initially the RA reaction between sulfur and oxygen until 20–40 days in all the three directions. After that, the NN reaction of sulfur with O₂ dominates other formation reactions. The primary destruction process after 100 days is the NN reaction of hydrogen with SO (products: sulfur and OH). Later, in the +Z- and −Z-directions, the CE reaction between H⁺ and SO and/or the NN reaction of carbon with SO (products: oxygen and CS; sulfur and CO) are dominant as the primary destruction process at around the rapidly decreasing phase (after 150 and 200 days, respectively); in the +Y-direction, the NN reaction of hydrogen with SO above remains the primary. In the +Y-direction, after about 1000 days, the mass fraction of SO distinctively increases again. This is the contribution from the innermost particles heated by the decay of ⁵⁶Ni, which partake in the molecule formation at such a later phase when the gas temperatures drop to 10⁴ K. The primary formation process is the NN reaction of oxygen with S₂ at this phase.

The behavior of H₂ is qualitatively similar compared with the counterpart model with the progenitor model b18.3. After 40 days, in all three directions, H₂ starts to form rapidly. In the +Z- and −Z-directions, the mass fraction of H₂ peaks at around 60 days. On the other hand, in the +Y-direction, the increase is slow and continues until 200 days. The primary formation processes of H₂ depend on the directions. In the +Z-direction, after 40 days, H₂ is primarily formed by the AD reaction, H⁻ + H → H₂ + e⁻. In the −Z-direction, H₂ is mainly formed by the NN reaction of hydrogen with CH and the NN reaction between two OH molecules; the AD reaction above contributes as the secondary and overtakes the NN reactions after 60 days. In the +Y-direction, the primary formation process of H₂ is the NN reaction of hydrogen with CH and/or the NN reaction between two OH molecules before 60 days; later, the NN reaction of hydrogen with SiH becomes the primary process. In the +Y-direction, the AD reaction contributes to the formation after 40 days as one of the secondaries. Therefore, the more extensively ⁵⁶Ni is initially distributed (see Figure 8), the more significant the contribution of the AD reaction is. This is because the number density of one of the reactants of the AD reaction, i.e., H⁻, increases due to the ionization by the CM reactions through the REA reaction, H + e⁻ → H⁻ + γ as mentioned.

Similarly to the models with the progenitor model b18.3, FeO and MgO are also formed to some extent; the mass fraction of MgO, however, dominates that of FeO in all three directions in contrast to the models with b18.3. This is because the initial mass fractions of magnesium are overall higher than those in the models with b18.3 (see Figures 7 and 8). After a few thousand days, in all three directions, FeO and MgO (in particular, FeO in the +Z- and −Z-directions) increase, although the increase of MgO in the +Y-direction is not remarkable. As described in the models with the progenitor model b18.3, such a later increase of FeO and MgO is attributed to the fact that the particles heated by the decay of ⁵⁶Ni with the initial high local ⁵⁶Ni abundance partake in the formation of the molecules at such a later phase; the heated particles contain plenty of the seed iron to form FeO. The primary formation process of MgO approximately before 100 days is the NN reaction of magnesium with O₂. Later, the 3B reaction Mg + O + H → MgO + H becomes the primary formation process until the end of the simulation in the +Z- and −Z-directions. In the +Y-direction, after the NN reaction-dominated phase, the 3B reaction among magnesium, oxygen, and helium also contributes as one of the primary processes. In all the three directions, the primary formation process of FeO is the NN reaction of iron with O₂ before 100 days. After a few thousand days, FeO increases again in particular in the +Z- and −Z-directions. At this phase, in the +Z- and −Z-directions, the 3B reaction, Fe + O + H → FeO + H, is the primary formation process; in the +Y-direction, the 3B reaction among iron, oxygen, and hydrogen (or second oxygen) is the primary depending on time.

4.2.4.3. Radial Distributions of Molecules in the Models with the Binary Merger Progenitor Model

So far, the directional dependence of the time evolution of the total mass fractions of molecules was discussed. Hereafter, the dependence of the radial distributions of molecules is described in detail.

In Figure 26, the radial profiles of the mass fractions of molecules and several seed atoms for the models b18.3-zp (top panels), b18.3-zn (middle panels), and b18.3-yp (bottom panels) are shown for two different epochs (approximately 200 days and the end of the calculation in the left and right panels, respectively). As mentioned above, during approximately 100–300 days, in particular, in the +Z- and −Z-directions, the destruction processes of some molecules are dominant (see left panels in Figure 24). As a representative snapshot of the destruction phase, the ones at approximately 200 days are selected.

At about 200 days (left panels), in the +Z (−Z)-direction, it is apparent that molecules are not formed at the regions inside ∼4 × 10¹⁵ cm (3.5 × 10¹⁵ cm). On the other hand, in the +Y-direction, molecules such as CO, SiO, and H₂ are remarkably formed even in inner regions. This is because, in the +Z- and −Z-directions, at this epoch, the gas temperatures of the inner particles are too high for molecule formation due to the gas heating by the decay of ⁵⁶Ni (see the right panels in Figure 24); on the other hand, in the +Y-direction, the local ⁵⁶Ni abundance is overall low (gas heating is less), and the molecule formation is already developed because of the early starting due to the early decrease of the gas temperatures to ∼10⁴ K. Again, the different initial ⁵⁶Ni abundances are due to the different energy depositions during the bipolar-like explosion.

In all three directions, the outer layers where all the mass fractions shown are approximately less than 10⁻³ and flat with some exceptions can be recognized. The transition points (radii) between the outer layers and the inner regions are similar in the three directions, although, in the +Z-direction, the radius of the point is slightly larger than those in the other two directions. In the outer layers, CO and SiO are also formed, although the mass fractions are not so high.

As for other molecules, in the outer layers, N₂, MgO, and FeO are distinctly formed; the formation of SO is also recognized, although the region is spatially limited (0.4–1.0 × 10¹⁶ cm). In the +Y-direction, N₂, MgO, and FeO are formed even in the inner regions; the formation of Si₂ can also be recognized at a limited region (0.1–0.4 × 10¹⁶ cm).

After 200 days, the regions where molecules are formed gradually extend inside as inner particles heated by the decay of ⁵⁶Ni partake in the molecule formation when the gas temperatures go down to ∼10⁴ K.

At the end of the calculation (right panels), in the +Z- and −Z-directions, the formation of molecules, e.g., CO, SiO, H₂, MgO, and FeO, is recognized at inner regions, although the mass fractions are not so high compared with those in the +Y-directions.

As seen in the +Z-direction, the mass fractions of, e.g., CO, SiO, and H₂ are not smoothly distributed in the inner regions (less than 4 × 10¹⁷ cm). In the region of (2.0–2.8) × 10¹⁷ cm (a dip region), the mass fractions of the three molecules above are apparently lower than those of the surroundings. It is noted that the transition point (∼4 × 10¹⁷ cm) at which the mass fractions of CO, SiO, and H₂ increase compared with those in the outer layers corresponds to the initial position of 2.4 × 10¹³ cm (see the top left panel in Figure 7) at which the mass fractions of the seed atoms increase.

What corresponds to the upper bound of the dip region? As seen in the top right panel in Figure 24, at 200 days, there is a gap in the gas temperatures between particles with and without a distinct gas cooling. Actually, the upper bound of the dip region (i.e., 2.8 × 10¹⁷ cm) corresponds to the temperature gap. Inside the upper bound (corresponding to the upper region from the temperature gap in the top right panel in Figure 24), particles go through the gas heating by the decay of ⁵⁶Ni without a distinct gas cooling. Then, the molecule formation is delayed compared with the particles outside the upper bound; some fraction of formed molecules may also be destructed by the processes involved with Compton electrons (CM reactions) and/or UV photons (UV reactions) and/or ions. On the other hand, outside the upper bound, the heating of the gas (ionization and/or destruction of molecules) is less effective, and some of the particles go through the gas cooling by CO rovibrational transitions; the molecule formation may partly be enhanced by the earlier start of the formation and the higher density caused by the cooling through the effect described in Equation (41).

Inside the lower bound of the dip region (i.e., 2.0 × 10¹⁷ cm), the mass fractions of CO, SiO, and H₂ are apparently higher than those in the dip region. What makes the difference?

First, it is found that the lower bound roughly corresponds to the initial position of 1.1 × 10¹³ cm at which the abundance of ⁵⁶Ni is a bit sharply changed (see the top left panel in Figure 7). Hence, inside the lower bound, the local abundance of ⁵⁶Ni is higher than that in the dip region. The decay of ⁵⁶Ni generally plays a role in the gas heating, ionization of atoms, and ionization and destruction of molecules. Therefore, considering the latter effect, it sounds a bit strange that, inside the lower bound, the formation of CO, SiO, and H₂ is efficient compared with the dip region.

The reason is probably that, inside the lower bound, the decay of ⁵⁶Ni actually enhances the formation of the three molecules above through the sequences as follows. The decay of ⁵⁶Ni increases the number density of electrons; H⁻ is produced by the REA reaction, H + e⁻ → H⁻ + γ; then, H₂ is produced by the AD reaction, H⁻ + H → H₂ + e⁻; OH is formed by the NN reaction, O + H₂ → OH + H; finally, CO and SiO are efficiently produced through the NN reactions in Equations (47) and (62) with OH. Actually, it is confirmed that, inside the lower bound, the contributions of the NN reactions in Equations (47) and (62) on the formation of CO and SiO are higher than those in the dip region.

What is the condition that the sequences that enhance the formation of CO and SiO dominate the destructive processes through the CM and UV reactions and/or ones involved with ions? Inside the lower bound of the dip region, the tracer particles are heated by the decay of ⁵⁶Ni, and the start of the formation of molecules is delayed. In this case, before the start of the molecule formation, the number density of electrons increases due to the ionization of the seed atoms by Compton electrons through the reaction in Equation (17). Hence, the abundance ratios of electrons to the other species become higher than that of the cases where the formation of molecules starts at an early phase of the active decay of ⁵⁶Ni, which may effectively increase the contribution from the formation processes in Equations (47) and (62) through the sequences above with the high electron abundance.

In the −Z-direction, similarly to the +Z-direction, at approximately (1.7–2.3) × 10¹⁷ cm, a sharp dip of CO, SiO, and H₂ is recognized. In the +Y-direction, compared with the profiles at 200 days, there are not so large differences, which means that the molecule formation is overall settled before 200 days.

As for other molecules, MgO and FeO are also formed in the inner regions in +Z- and −Z-directions. In all three directions, SO is formed at the bottom of the outer layers, and the distributions are not so changed from those at 200 days. In the +Z- and −Z-directions, at the end of the calculation, the formation of CS is recognized just inside the region where SO is formed with some overlaps.

4.2.4.4. Radial Distributions of Molecules in the Models with the Single-star Progenitor Model

In Figure 27, the profiles for the models, n16.3-zp (top panels), n16.3-zn (middle panels), and n16.3-yp (bottom panels), are shown. At about 200 days (left panels), it is apparent that, in the +Z-, −Z-, and +Y-directions, the inner regions where molecules are remarkably formed with the mass fractions over ∼10⁻³ are extended to ∼9 × 10¹⁵ cm, 4 × 10¹⁵ cm, and 5 × 10¹⁵ cm, respectively. The radii roughly reflect the initial extensions of the regions where the seed atoms are abundant (see Figure 8). In the +Z- and −Z-directions, there are inner regions where molecules are not distinctly formed; the regions are extended to roughly 4 × 10¹⁵ cm and 2 × 10¹⁵ cm, respectively. On the other hand, in the +Y-direction, molecules are remarkably formed except for the innermost region.

Another distinct feature compared with the models b18.3-zp, b18.3-zn, and b18.3-yp is the formation of O₂. In the models with the progenitor model b18.3, the formation of O₂ is not evidently recognized. As seen in the remained seed atoms at 200 days in the +Z- and −Z-directions, the mass fractions of the remained oxygen atoms are apparently larger than those in the models with b18.3. The high mass fractions of O₂ compared with the models with b18.3 are basically attributed to the initial high oxygen abundance. Hereafter, for each direction, the properties of O₂ and other major molecules at about 200 days are described.

In the +Z-direction, the mass fraction of O₂ peaks at approximately 5.8 (4.7) × 10¹⁵ cm with the mass fraction of ∼10⁻² (10⁻³), and in the other regions within (3.5–8.2) × 10¹⁵ cm, the mass fraction is roughly the order of 10⁻⁵. In contrast to O₂, CO, SiO, and H₂ are extensively formed within the regions of roughly (3.5–9.0) × 10¹⁵ cm. In the regions where the mass fraction of O₂ is ∼10⁻⁵, O₂ is actually more abundant at an earlier time. Then, O₂ is destructed (converted to other molecules) by the NN reactions, Si + O₂ → SiO + O and C + O₂ → CO + O, approximately from 60 days. Whether O₂ survives or not may sensitively be determined by a bit complicated balance of NN reactions depending on the thermal history. It is noted that, in the evolution of the gas temperatures (the top right panel in Figure 25), e.g., at 200 days, a gap of the gas temperatures (approximately from 1.5 × 10³ K to 4.5 × 10³ K) is recognized. Some of the particles whose temperatures at 200 days are below the lower bound have complicated temperature histories affected by both the heating and cooling. The regions of approximately (4.2–5.5) × 10¹⁵ cm are affected by both the heating and cooling.

In the −Z-direction, the formation of O₂ is more significant; the distribution of O₂ is roughly consistent with those of CO, SiO, and H₂, and there are regions where the mass fraction of O₂ is higher than that of CO with the peak value of ∼10⁻¹. It is noted that the regions of (1.7–4.0) × 10¹⁵ cm correspond to the initial positions of (0.9–2.0) × 10¹³ cm (see the middle left panel in Figure 8); the radial fluctuations of the mass fractions of seed atoms, in particular, oxygen and silicon, seen in the initial profiles seem to be somehow reflected in the profiles of those of CO, SiO, and O₂ at 200 days.

In the +Y-direction, O₂ is significantly formed inside ∼3 × 10¹⁵ cm; the distribution is roughly consistent with that of SiO. The peak value of the mass fraction of O₂ is as high as 0.4 at approximately 1 × 10¹⁵ cm. In contrast to O₂ and SiO, the regions where the mass fractions of CO and H₂ are high (greater than 10⁻³) are more extended to 5 × 10¹⁵ cm. The regions of (1.0–3.0) × 10¹⁵ cm and (3.0–5.0) × 10¹⁵ cm roughly correspond to the initial positions of (0.6–1.5) ×10¹³ cm and (1.5–2.4) × 10¹³ cm, respectively (see the bottom left panel in the Figure 8). In the latter region in the initial profiles, the mass fractions of oxygen and silicon are low compared with those in inner regions (the mass fraction of silicon is the order of 10⁻⁵); instead, the mass fractions of hydrogen and helium are high. On the other hand, in the former regions, there is plenty of oxygen and silicon. The distinct different features across 3.0 × 10¹⁵ cm seen in the profiles at 200 days are attributed to the initial distributions of seed atoms.

As for other molecules, in the +Z-direction, MgO, FeO, and N₂ are also formed in the inner regions with similar distributions as CO and SiO at (4.5–9.0) × 10¹⁵ cm. In particular, the mass fraction of MgO is comparable to that of CO in the regions. The formation of SO is also recognized at around the positions of the peaks of O₂. In the regions inside about 6.0 × 10¹⁵ cm, FeO slightly increases inward. The increase is due to the seed iron originating from ⁵⁶Ni. In the −Z-direction, the formation of MgO, FeO, N₂, and SO are recognized. Same as the +Z-direction, MgO dominates FeO. MgO is abundantly formed in the regions of (2.0–4.0) × 10¹⁵ cm and peaks at around 3.5 × 10¹⁵ cm with the value of 5 × 10⁻². N₂ and SO are a bit more abundantly formed in the regions of (2.0–3.0) × 10¹⁵ cm. In the +Y-direction, first, in the regions of (3.0–5.0) × 10¹⁵ cm, the formation of C₂, SiN, CN, and CS is recognized. In the regions of (1.0-3.0) × 10¹⁵ cm, the formation of MgO and SO is evident. N₂ is formed in both regions above. In the former regions, the formation of the molecules mentioned is attributed to the initial abundances (abundance ratios) of the seed atoms, i.e., the relatively low oxygen and silicon abundances compared with the inner regions.

After 200 days, the particles heated by the decay of ⁵⁶Ni gradually partake in the molecule formation from the outer ones.

At the end of the calculation, in the +Z-direction, the distributions of the molecules in the regions of approximately (2.4–5.3) × 10¹⁷ cm are overall consistent with those in the regions of (3.5–9.0) × 10¹⁵ cm at 200 days, although the mass fractions of SiO and O₂ in the regions of (2.4–3.5) × 10¹⁷ cm are a bit different from those at 200 days. In the regions inside 2.4 × 10¹⁷ cm, the formation of CO, SiO, and H₂ is recognized; the distributions, however, are not smooth. In these regions, particles are affected by the decay of ⁵⁶Ni. As found in the model b18.3-zp, a dip region, i.e., (2.2–2.4) × 10¹⁷ cm, where the mass fractions of CO, SiO, and H₂ are small compared with the surroundings is recognized. In the regions inside the upper bound, tracer particles are affected by the heating due to the decay of ⁵⁶Ni without a distinct cooling through CO rovibrational transitions. As mentioned, the decay basically plays a role in the ionization and destruction of molecules; it, however, could enhance the formation of CO and SiO probably in certain cases if the molecule formation starts well after the decay becomes effective. The transition point whether the destruction or enhancement is dominant is at around 2.2 × 10¹⁷ cm in the model n16.3-zp as also seen in the model b18.3-zp (2.0 × 10¹⁷ cm for b18.3-zp; see the top right panel in Figure 24). In the regions inside the outer bound of the dip region (≲2.4 × 10¹⁷ cm), in particular, in the dip region, the formation of MgO and FeO is prominent. It is interpreted that the destruction of CO and SiO by the decay of ⁵⁶Ni increases the seed oxygen abundance; then, the oxygen is converted to MgO and FeO. The high abundance ratio of FeO to MgO in this region than that in the outer layers may be due to the abundant seed iron atoms originating from ⁵⁶Ni. Although the mass fraction is rather small (the order of 10⁻⁶), the formation of SiS can also be recognized at the regions inside 2 × 10¹⁷ cm.

In the −Z-direction, the position at 1.6 × 10¹⁵ cm at 200 days corresponds to 1.1 × 10¹⁷ cm at the end of the calculation. Therefore, the mass fractions in the regions of approximately (1.1–2.2) × 10¹⁷ cm are further varied from those at 200 days. In particular, the radial distribution of O₂ becomes complicated; SO and NO have qualitatively similar trends. It is noted that the regions above correspond to the particles affected by both the heating and cooling (as seen in the middle right panel in Figure 25). The complicated gas temperature (density) evolution may partly play a role in the complex distributions. In the regions inside 1.1 × 10¹⁷ cm, the formation of CO, SiO, H₂, MgO, and FeO is recognized. As in the models n16.3-zp, b18.3-zp, and b18.3-zn, a dip in the mass fractions of CO, SiO, and H₂ can be recognized because of the same reason.

In the +Y-direction, the same as the model b18.3-yp, the mass fractions of molecules are barely changed from those at 200 days, since the molecule formation is already settled before 200 days.

5. Discussion

In this section, several issues related to the results of this paper are discussed.

5.1. Implication to the 3D Distributions of CO and SiO

As noted in Section 1, ALMA observations have revealed 3D distributions of line emissions (rotational transition lines) from CO and SiO (Abellán et al. 2017). The distributions of both CO and SiO delineated by several isosurfaces seem to be rather lumpy and nonspherical. In particular, the distribution of CO shows a torus-like (ring-like) structure perpendicular to the equatorial ring in the nebula of SN 1987A. Additionally, SiO seems to be overall engulfed by CO with some exceptions (in some directions, SiO is extended farther than CO), and both CO and SiO may have a hole in inner regions. In our previous study (Ono et al. 2020), the directions of the bipolar-like explosions against the equatorial ring were derived to be consistent with the observed [Fe ii] lines (Haas et al. 1990). The farther evolutions obtained (Orlando et al. 2020) by the 3D MHD simulations (consistent with the derived directions above) result in dipole-like iron distributions roughly parallel to the equatorial ring in the nebula. In Orlando et al. (2020), as a reference, the distribution of CO (SiO) was naively estimated from the simulation results to be the square root of the product of the (mass) densities of carbon and oxygen (silicon and oxygen). In the estimation, qualitatively, a torus-like structure of CO is somehow recognized (see Figure 10 in Orlando et al. 2020), and it is also roughly perpendicular to the equatorial ring as seen in the ALMA observations.

In this paper, the molecule formation in the representative directions (the bipolar-like explosion directions and one perpendicular to the former) is calculated. For an implication to the 3D distributions, here, the calculated distributions of CO and SiO in these directions are a bit more discussed. In Figure 28, the radial distributions of the number densities of CO and SiO (at the end of the calculation) in the +Z (top panels), −Z (middle panels), and +Y (bottom panels) -directions for both the progenitor models, i.e., the ones for the models b18.3-zp, b18.3-zn, b18.3-yp, n16.3-zp, n16.3-zn, and n16.3-yp, are shown. For reference, similarly to the previous study (Orlando et al. 2020), the distributions of the approximated number densities of CO and SiO, i.e., $\sqrt{N({\rm{C}})\times N({\rm{O}})}$ (C × O) and $\sqrt{N(\mathrm{Si})\times N({\rm{O}})}$ (Si × O), are also plotted, where N(C), N(O), and N(Si) are the number densities of the seed carbon, oxygen, and silicon atoms for the case with no chemical evolution.

**Figure 28.** Radial distributions of the number densities of CO and SiO (solid lines) at the end of the calculation for the models b18.3-zp, b18.3-zn, b18.3-yp, n16.3-zp, n16.3-zn, and n16.3-yp. For reference, approximated number densities of CO and SiO (dashed lines, C × O and Si × O, respectively) are also plotted. Here, C × O (Si × O) denotes $\sqrt{N({\rm{C}})\times N({\rm{O}})}$ ( $\sqrt{N(\mathrm{Si})\times N({\rm{O}})}$ ), where N(C), N(O), and N(Si) are the number densities of carbon, oxygen, and silicon atoms, respectively, for the case with no chemical evolution, i.e., the case that each tracer particle keeps the initial mass fractions of the species. From top to bottom, the models for the directions, +Z, −Z, and +Y, are shown, respectively. The left panels (right panels) are for the models based on the binary merger progenitor model b18.3 (single-star progenitor model n16.3).
Download figure:
Standard image High-resolution image

**Figure 28.** Radial distributions of the number densities of CO and SiO (solid lines) at the end of the calculation for the models b18.3-zp, b18.3-zn, b18.3-yp, n16.3-zp, n16.3-zn, and n16.3-yp. For reference, approximated number densities of CO and SiO (dashed lines, C × O and Si × O, respectively) are also plotted. Here, C × O (Si × O) denotes $\sqrt{N({\rm{C}})\times N({\rm{O}})}$ ( $\sqrt{N(\mathrm{Si})\times N({\rm{O}})}$ ), where N(C), N(O), and N(Si) are the number densities of carbon, oxygen, and silicon atoms, respectively, for the case with no chemical evolution, i.e., the case that each tracer particle keeps the initial mass fractions of the species. From top to bottom, the models for the directions, +Z, −Z, and +Y, are shown, respectively. The left panels (right panels) are for the models based on the binary merger progenitor model b18.3 (single-star progenitor model n16.3).
Download figure:
Standard image High-resolution image

As can be seen, in all three directions, CO is more extended to outer regions than SiO for both the progenitor models, if the positions of the same number densities are compared. Compared with the approximated distributions, the distribution of CO in the outer regions is roughly consistent with the approximation (C × O). On the other hand, the number density of SiO in the outer regions is quantitatively small compared with the approximation (Si × O) and qualitatively different from the approximation. The outer tails of the approximated SiO seen in some models, e.g., in the model n16.3-yp (see ≳2 × 10¹⁷ cm in the bottom right panel), cannot be seen in the calculation. The most distinct difference between the calculation results and the approximations is that the approximations significantly overestimate the number densities of inner regions in the +Z- and −Z-directions; the calculated distributions in the inner regions have complicated structures. The difference is at least partly attributed to the decay of ⁵⁶Ni through the delay of the molecule formation by the heating, the destructive processes by Compton electrons and UV photons, and the sequences induced by the decay for increasing the amounts of CO and SiO as mentioned.

In the models based on the progenitor model b18.3 (left panels), in the +Z- and −Z-directions, at around 1 × 10¹⁷ cm, there are some dips for both CO and SiO. Additionally, the position at the peak of CO is larger than that of SiO except for the innermost region. In the +Y-direction, the maximum number density is more than 1 order of magnitude higher than those of the two directions above, and the number densities of both CO and SiO roughly decrease as the radius increases except for the innermost region.

In the models based on the progenitor model n16.3 (right panels), in the +Z- and −Z-directions, the distributions of CO and SiO are qualitatively similar, and the positions of the peaks of CO and SiO are consistent with each other. In the +Y-direction, the qualitative feature is similar to the case of b18.3.

Since the binary merger progenitor model better explains several observational features of SN 1987A (Ono et al. 2020; Orlando et al. 2020), hereafter, the models based on the progenitor model b18.3 are focused. According to the descriptions above, in a solid angle along the +Z and −Z (polar) -directions, the distributions of CO and SiO may show a dipole-like (polar cap) structure and a gap (hole) in the inner regions. Around the perpendicular plane to the polar directions, with the gap in the polar directions and the higher number densities than those in the polar directions, the distributions may or may not look like a torus. For both the dipole- and torus-like structures, CO is overall more extended to outer regions than SiO. Therefore, some of the observed qualitative features may be consistent with the current models. Such features, however, should be confirmed by the direct application to the 3D models (Ono et al. 2020; Orlando et al. 2020). Since the distributions of CO and SiO are affected by the explosion asymmetries and the matter mixing through the distribution of ⁵⁶Ni and the initial abundance ratios of the seed atoms as presented, the observed 3D distribution of CO and SiO (Abellán et al. 2017) and further observations in the future may provide insights on the explosion mechanism and/or matter mixing by comparing with such 3D theoretical models.

5.2. Issues of the Calculation Method in This Study

As presented in Section 4, the evolution of the gas temperatures of the particles is an important key for the formation of molecules. In this study, as described in Section 2, in order to take into account the heating by the decay of ⁵⁶Ni, it is assumed that some fractions of the released energies are only locally deposited depending on the optical depth for gamma rays for simplicity. Additionally, to control the efficiency of the gas heating, a constant parameter, f_h, is introduced. In reality, the energy depositions, however, could occur nonlocally (potential impact of nonlocal energy deposition is discussed in Appendix B as mentioned; the effective nonlocal energy deposition model shows that it may change the amounts of molecules by at most a few tens of percent), and the deposited energies could be consumed by not only gas heating but also ionization and excitation (e.g., Liu & Dalgarno 1995); the fraction of energy for each process (gas heating or ionization or excitation) depends on the ionization fractions (Liu & Dalgarno 1995). Therefore, the efficiency of the heating (the fraction of the deposited energies to be used for the heating), f_h, is not generally constant. In order to obtain a more realistic temperature evolution model, the energy depositions via the decay of ⁵⁶Ni should be treated by taking into account the ionization, excitation, and gas heating simultaneously with line transitions (emissions) from ionized atomic species as another cooling source, which is not taken into account in this study, as calculated, e.g., in the reference (Kozma & Fransson 1998). Recently, van Baal et al. (2023) reported models for spectral lines in 3D based on 3D hydrodynamical simulations incuding those effects above.

As for the CO rovibrational transitions, it is difficult to reproduce the features of the observed light curves of the fluxes of the fundamental and first overtone bands (Meikle et al. 1989, 1993; Bouchet & Danziger 1993; Wooden et al. 1993), i.e., the dome-like shapes with the peaks at around 200–250 days. First, the calculations in this study tend to result in fluxes higher than the observed peak fluxes before the timings of the observed peaks. Then, an arbitrary reduction factor, f_red, to the escape probability is introduced expecting some uncertainty in the expression of the line optical depth. There is, however, no concrete physical basis to validate such an arbitrary factor. Second, in all the calculated models, the fluxes at around 200 days are apparently underestimated compared with the observations, which means that the excitations of CO vibrational levels at around the observed peaks may be insufficient. As mentioned, naively, we would expect that supra-thermal electrons originating from the decay of ⁵⁶Ni may play a role in the excitations. Actually, the rough coincidence between the timing of the peaks of the CO vibrational bands and the peaks of the gas temperatures of the particles heated by the decay of ⁵⁶Ni may support the expectation. As mentioned in Liu & Dalgarno (1995), the optical depth of the gamma rays from the decay of ⁵⁶Co could be a key to reproducing the observed peaks. Additionally, in the Sobolev approximation adopted in this paper, the line optical depths (not for the gamma rays but for the CO rovibrational lines) are locally determined, and the trapped photons controlled by the escape probabilities, which depend on the optical depths, can only locally be absorbed. If the Sobolev approximation is broken, emitted photons might be trapped nonlocally, and the emergence of the photons might be delayed, which may also play a role in the light curves. It would be necessary to take into account the energy distribution of electrons in the excitation processes of CO vibrational levels and the potential nonlocal effects mentioned above. Once a reasonable model for the CO rovibrational transitions is obtained, the time-dependent spectra of the CO rovibrational line emissions can be estimated in the future.

The balance of the heating and cooling is important for the evolution of the gas temperatures and eventually the formation of molecules. Some improvements in the treatments on the energy depositions from the decay of ⁵⁶Ni and the CO rovibrational transitions are necessary to obtain more realistic thermal histories of the ejecta particles in the future.

It is noted that, in this study, the constant parameters, f_h, f_d, and t_s, are introduced, and the values are somehow calibrated by comparing with the observed fluxes of CO vibrational bands and the estimations in the previous studies. The reasonable parameter values are different between one-zone models and 1D calculations with the angle-averaged profiles based on the 3D model (b18.3-high). Considering the application to the 3D models without angle-averaging in the future, the selection of such parameters could matter because it is difficult to perform wide-range parameter surveys due to the numerical cost. The selection of the different parameter values may partly be to compensate for the different degrees of matter mixing (practically, the mixing of ⁵⁶Ni) as mentioned. The degree of matter mixing in one-zone models would be higher than those in 1D models because of the uniform composition in one-zone models. Since, in reality (in the 3D models), the degree of matter mixing may be milder than the angle-averaged 1D profiles, the reasonable parameter values for the direct application to the 3D models would be the extrapolation between those in the one-zone and 1D models or the boundary value (maybe for f_d). In the 3D models, the matter mixing is probably more realistic than that in the one-zone and 1D models. Those parameters may be determined as more reasonable values in the 3D application models.

However, such parameters could depend on the time and region because of the different physical conditions in reality. For the quantitative comparison of observations, reducing uncertainties, and disentangling degeneracies of the problem, again, more realistic treatments for the issues above are indispensable in the future.

5.3. Implication to the Observations by JWST

JWST²⁶ is an infrared observatory covering wavelengths of 0.6–28 μm²⁷ launched in 2021 December. The fundamental (∼4.6 μm) and first overtone (∼2.3 μm) bands of CO vibrational transitions and the fundamental (∼8 μm) and first overtone (∼4 μm) bands of SiO can potentially be covered by the telescope. Actually, there has been an attempt for theoretical modeling (Liljegren et al. 2023) of the spectra of the rovibrational line emissions of not only the CO and SiO but also SO and SiS for type Ib/c SNe (at a typical distance of 10 Mpc); it is concluded that most molecular emissions can be observed with a good signal-to-noise ratio in the 100–400 days time window. Unfortunately, such rovibrational transitions cannot be expected for SN 1987A at the current age, but, for general type II and type Ib/c SNe, spectral modelings of the rovibrational line emissions would be useful for comparison with future JWST observations, which will shed light on the formation of molecules in the ejecta of CCSNe by comparing theoretical models.

The spatial resolution of JWST is as small as ∼0 farcs 1. Therefore, nearby SNRs, e.g., SN 1987A (distance 51.4 kpc) and Cas A (distance ∼3.4 kpc), can spatially be resolved. Actually, as mentioned in Section 1, there has been an approved proposal for Cycle 1 JWST observations of SN 1987A (Matsuura et al. 2021); the observations will spatially resolve hot dust grains and may figure out how dust is destructed by the blast wave and the reverse shock. Resent JWST NIRCam imaging of SN 1987A and spectral decomposion have shown the potential to distinguish overlapping emisions from different components (Arendt et al. 2023). There has also been an approved proposal for observations of Cas A by JWST (Milisavljevic et al. 2021) for investigating the properties of dust, i.e., the dust survival against the passage of the reverse shock, grain size distribution, and clump size, and so on. The initial survey of Cas A with JWST NIRCam and MIRI imaging (Milisavljevic et al. 2024) has already provided some information to address those properties. It is noted that, recently, the results of the first observations of SN 1987A by JWST have been reported (Larsson et al. 2023) as mentioned in Section 1; the 3D distribution of the ejecta delineated by [Fe i] shows a broken dipole-like structure, which is roughly parallel to the equatorial ring. The derived 3D morphology of the ejecta and the orientation seem to be in good agreement with those of our 3D models (Ono et al. 2020; Orlando et al. 2020), which may further support our 3D models.

Moreover, in the case of SN 1987A, there have also been observations of spatially resolved dust in the ejecta of SN 1987A by ALMA (Cigan et al. 2019). Therefore, in light of such spatially resolved observations of SN 1987A and Cas A, it is indispensable to make models of the formation and destruction of dust based on 3D hydrodynamical models to be compared with JWST observations. As for SN 1987A, we plan to apply our methodology in this paper to the 3D hydrodynamical models (Ono et al. 2020; Orlando et al. 2020) more directly. Furthermore, we also plan to connect our molecule formation calculations in the ejecta to dust formation and destruction theories (e.g., Nozawa et al. 2010; Nozawa & Kozasa 2013) in the near future. As for Cas A, such formation and destruction models of molecules and dust can be applied to the recent 3D hydrodynamical models for Cas A (Orlando et al. 2021, 2022) in the future.

5.4. A Caveat: Possible Impacts of Realistic 3D Explosion Models

In this study, the impact of matter mixing on molecule formation in CCSN ejecta was discussed based on the 3D bipolar-like explosion models for SN 1987A (Ono et al. 2020), in which the explosions are initiated with the rather simplified model, i.e., instantaneous asymmetric kinetic and thermal energy injections, as mentioned. Then, based on the obtained matter mixing, the impact of effective matter mixing was investigated by using the one-zone models, the angle-averaged 1D profiles, and the extracted radial profiles in the specific directions from the 3D models. However, in reality, CCSN explosions should be driven via complicated physical processes as observed in the theories of CCSN explosion mechanisms, e.g., neutrino-driven and magnetorotationally driven explosions (see Section 1). The 3D models adopted in this study may fail to capture some features obtained by more realistic 3D CCSN explosion models (e.g., Mösta et al. 2014; Vartanyan et al. 2019, for magnetorotationally driven and neutrino-driven machanisms, respectively). Therefore, such more realistic 3D CCSN explosion models may result in different matter mixing, i.e., distribution of the seed atoms and ⁵⁶Ni and the abundance ratios, from the adopted 3D models and potentially change the results on molecule formation presented in this paper in particular the ones for the specific directions in the 3D models (Section 4.2.4) and even for the one-zone and 1D models with the angle-averaged 1D profiles (Section 4.1 and Sections 4.2.1, 4.2.2, and 4.2.3, respectively). It is worth investigating the impact of more realistic 3D CCSN explosion models on the formation of molecules in the ejecta in the future.

6. Summary

In this paper, the molecule formation in the ejecta of CCSNe is calculated by using a chemical reaction network based on one-zone and 1D ejecta evolution models based on the two 3D hydrodynamical models for SN 1987A (b18.3-high and n16.3-high; Ono et al. 2020), i.e., asymmetric bipolar-like explosions with the binary merger (b18.3; Urushibata et al. 2018) and single-star (n16.3; Nomoto & Hashimoto 1988; Shigeyama & Nomoto 1990) progenitor models. In Section 4.1, with one-zone models, the impact of the parameters related to the heating and cooling of gas is investigated (Section 4.1.1). Then, with a reasonable set of the parameters, important chemical reactions for the formation and destruction of CO and SiO molecules are described (Section 4.1.2). Similar to the one-zone models, in Section 4.2, the impact of the parameters is investigated for 1D models, and a set of parameters is derived as a fiducial case (Section 4.2.1). With the parameters of the fiducial case, by comparing the model results among the 1D models with the initial profiles obtained by the angle-averaging of the bipolar-like explosion models and ones for the spherical explosion case (the explosions are spherical, but the evolution is 3D), and ones for the purely spherical case, the impact of matter mixing on the molecule formation in the ejecta is effectively investigated (Section 4.2.3). Finally, from the model results of the calculations for specified directions in the 3D hydrodynamical models (+Z-, −Z-, and +Y-directions; the positive and negative directions of the bipolar-like explosion axis and a perpendicular direction to the former two), the dependence of the molecule formation on the directions is examined (Section 4.2.4). The main points and findings in this paper are summarized with some compliments as follows.

1.
The angle-averaged 1D profiles (Figure 3) derived from the 3D models of the asymmetric bipolar-like explosions (b18.3-high and n16.3-high; Ono et al. 2020) effectively reflect the mixing of ⁵⁶Ni into outer layers and the mixing of seed atoms. In the angle-averaged 1D profiles, ⁵⁶Ni is rather extended to outer layers compared with the spherical explosion (but the evolution is in 3D) case (Figure 4) and purely spherical case (Figure 5). In the spherical explosion case, ⁵⁶Ni is more concentrated in inner regions, and the seed atoms are smoothly distributed. In the purely spherical case, ⁵⁶Ni is confined only in the innermost regions, and the seed atoms are rather stratified. The one-zone model derived from the radial averaging of the angle-averaged 1D profiles for the model b18.3-high represents the case that the mixing is more efficient compared with the angle-averaged 1D profiles because of the uniform composition.
2.
The decay of radioactive ⁵⁶Ni and/or ⁵⁶Co (hereafter, simply, the decay of ⁵⁶Ni as before this section) could heat the gas and increase the gas temperature. If the efficiency of the gas heating by the decay of ⁵⁶Ni, i.e., f_h, is higher, the timing when the gas temperature goes down to ∼10⁴ K, at which molecules start to form, is more delayed compared with lower f_h cases. Compared with lower f_h cases, molecules start to form in a low-density environment, which generally reduces the amounts of molecules, e.g., CO and SiO.
3.
The decay of ⁵⁶Ni produces Compton electrons, which could ionize and destruct molecules through the processes described in Equations (18) and (19). The secondary UV photons due to the decay could also dissociate molecules. If the destruction efficiency, f_d, is higher, the ionization of molecules by Compton electrons and the destruction of those by Compton electrons and/or UV photons become more effective; the higher f_d, generally, the more effective the destruction of CO and SiO at around a few hundred days after the explosion (however, under a certain condition, formation processes could locally dominate the destructive ones instead; see the seventeenth point below).
4.
In the calculations in this study, the fluxes of CO vibrational bands tend to become higher than the observed peak fluxes (at 200–250 days) at an early phase before 100 days. In order to avoid such early high fluxes, an arbitrary reduction factor, i.e., f_red, had to be introduced. Though, it is difficult to reproduce the observed features of the CO bands, i.e., the peaks at 200–250 days and the dome-like shapes. At around the observed peaks, additional excitations by supra-thermal electrons produced by the decay of ⁵⁶Ni might play a role in reality.
5.
According to the fiducial one-zone model (see the sixth point below), initially, CO and SiO are primarily formed by the RA processes in Equations (45) and (61), respectively, until a few tens of days; after that, many NN reactions are involved with both the formation and destruction processes. If the atomic ions are produced by the ionization due to Compton electrons in Equation (17), the IN reactions in Equations (53), (54), (70), and (69) (the first is for CO formation and SiO destruction, the second is for the destruction of CO, and the third and fourth are the destruction of SiO) could play an important role in the formation and/or destruction of CO and SiO. For contributing reactions for the formation and destruction of CO and SiO in the fiducial one-zone model, see Figures 12 and 13, respectively. The contributing reactions, however, generally depend on the thermal history and the initial composition; actually, some differences between the fiducial one-zone and 1D calculations can be recognized (see the seventh point). The previous studies on the formation of CO and SiO in SN 1987A (Lepp et al. 1990; Liu et al. 1992; Liu & Dalgarno 1994, 1995, 1996; Sluder et al. 2018; Liljegren et al. 2020) are useful for reference. Actually, some of the chemical reactions described in this paper have also been mentioned in the studies above. For chemical reactions in general/other types of SNe, e.g., the reference (Cherchneff & Dwek 2009) is useful.
6.
To be somehow consistent with the estimates of the amounts of CO and SiO derived from the previous studies (Liu et al. 1992; Liu & Dalgarno 1994, 1995; Matsuura et al. 2017; Liljegren et al. 2020) and the light curves of the observed CO vibrational bands (Meikle et al. 1989, 1993; Bouchet & Danziger 1993; Wooden et al. 1993), for the one-zone and 1D models, the parameter sets of (f_h, f_d, t_s) = (5 × 10⁻⁴, 10⁻², 200 days) and (5 × 10⁻³, 1.0, 500 days) might be selected as the reasonable sets, respectively (for the former set, see the middle panels in Figure 9; for the latter, see, e.g., the lower panels in Figure 14 and the top panels in Figure 20). It is interpreted that, in the one-zone models, lower values of the efficiencies of the gas heating, the ionization and destruction by Compton electrons, and the dissociation by UV photons may be preferred to compensate for the more effective matter mixing compared with the 1D models.
7.
The contributing chemical reactions for the formation of CO and SiO in the fiducial 1D model (see the sixth point) are not significantly different from those in the fiducial one-zone model; some of the primary processes, and the relative significances could be different between the two models depending on time (see Figures 12, 13, 18, and 19). The less effective mixing of ⁵⁶Ni, hydrogen, and helium in the fiducial 1D model could decrease the contributions of the reactions involved with ionized hydrogen and helium, i.e., the CE reactions in Equations (53) and (68) and the IN reactions in Equations (54), (69), and (70) compared with the fiducial one-zone model.
8.
The abundance ratios of the initial seed atoms are different (see, e.g., Figure 3 and the top left panels in Figures 20 and 22) between the models based on the binary merger progenitor model (b18.3) and the single-star progenitor model (n16.3); in the models based on n16.3, the abundance ratios of oxygen to carbon and silicon are apparently higher compared with those in the models based on b18.3; in the models based on n16.3, the abundance ratio of carbon to silicon is also higher than those in the models based on b18.3. Such different initial abundance ratios of the seed atoms could affect the consequential amounts of molecules. Actually, in the models based on n16.3, the amount of O₂ is higher than that in the models based on b18.3.
9.
MgO and FeO are primarily formed by the NN reactions of magnesium and iron with O₂ initially and by the 3B reactions, Mg + O + H → MgO + H and Fe + O + H → FeO + H, at later phases, respectively. In the models based on the single-star progenitor model (n16.3), the amount of MgO is higher than that of FeO. On the other hand, the amount of FeO dominates that of MgO at the final phase in the models based on the binary merger progenitor model (b18.3). This feature is attributed to the initial high magnesium abundance in the models based on n16.3.
10.
The formation of N₂ is distinctly recognized regardless of the models. Moreover, the amounts of N₂ in the models based on the binary merger progenitor model (b18.3) are overall higher than those in the models based on the single-star progenitor model (n16.3). This feature is probably attributed to the higher initial abundance of nitrogen in the models with b18.3, in particular, in the envelope (see, e.g., Figure 3); in the binary merger progenitor model, the CNO cycle additionally triggered during the merger process enhances the nitrogen (¹⁴N, one of the main products of the CNO cycle) abundance in the envelope via mixing.
11.
The fluxes of CO vibrational bands can be affected by the mixing of ⁵⁶Ni. In higher mixing efficiency models (the efficiencies are higher in the asymmetric explosion models, spherical explosion models, and purely spherical models, in this order; see the first point), the ionization by Compton electrons due to the decay of ⁵⁶Ni increases the number density of thermal electrons to excite CO vibrational levels. Actually, it is found that the calculated fluxes of CO vibrational bands are overall higher in high mixing efficiency models compared with lower mixing efficiency (more spherical) models.
12.
In the purely spherical models (the models b18.3-sphel-pure and n16.3-sphel-pure), due to the rather stratified distributions of the initial seed atoms, the formation of some molecules that are not markedly formed in the other models, e.g., C₂, CS, SiC, SiN, are recognized. In particular, in both the models b18.3-sphel-pure and n16.3-sphel-pure, the amounts of C₂ are significantly higher than those in the other (nonspherical) models and become over 10⁻² M_⊙. In the model n16.3-sphel-pure, further, due to the high initial oxygen abundance, O₂ is significantly formed, and the amount is the highest among the molecules and reaches ∼1 M_⊙.
13.
Depending on the direction, the distributions of the seed atoms and ⁵⁶Ni are different in the last snapshots of the 3D hydrodynamical models (b18.3-high and n16.3-high). In particular, the explosions are bipolar-like and the distributions are distinctively different between the directions along the bipolar axis (+Z and −Z) and a direction (+Y) perpendicular to the axis (in the equatorial plane). In the bipolar axis, ⁵⁶Ni is extensively distributed to outer layers (see Figures 7 and 8). Then, in the bipolar axis, the gas temperatures of the inner particles heated by the decay of ⁵⁶Ni increase making peaks at around 200 days; such gas heating delays the start of the formation of molecules. On the other hand, in the equatorial plane, ⁵⁶Ni is concentrated only in the innermost region; the gas temperatures of most of the particles go down to ∼100 K within a few hundred days.
14.
In the directions along the bipolar axis (+Z- and −Z-directions) for both the progenitor models (b18.3 and n16.3), some fractions of CO and, in particular, SiO are destructed approximately from 100 days to a few hundred days. At this phase, the main destruction processes of CO are the ionization by Compton electrons in Equation (18), the destruction by UV photons, and the destruction by Compton electrons in Equation (19). SiO is mainly destructed by the CE reaction in Equation (68) and the NN reactions in Equations (66) and (67). On the other hand, in the equatorial plane, there is no distinct destruction by processes related to Compton electrons.
15.
In the +Z- and −Z-directions, FeO increases even after a few thousand days by the participation of the inner particles heated by the decay of ⁵⁶Ni, which contain plenty of seed iron, to the formation of FeO by a 3B reaction. Probably by reflecting the initial abundance ratios of the seed atoms, in the models based on the progenitor model n16.3, the amount of MgO dominates FeO in all three directions. On the other hand, in the models based on the binary merger progenitor model b18.3, in the +Z- and −Z-directions, the amount of FeO dominates MgO; in the +Y-direction, MgO dominates FeO.
16.
The formation of molecules depends on the radial position. Molecules start to form from the outer regions, since, from these regions, the gas temperatures go down to ∼10⁴ K. At 200 days, in the +Z- and −Z-directions, the formation of molecules is only recognized (from the outer regions) up to the outer part of the regions where, initially, the seed atoms are abundant (see, e.g., the top and middle left panels in Figures 7 and 26). On the other hand, in the +Y-direction, at 200 days, molecules such as CO and SiO are already remarkably formed even in the inner regions.
17.
The decay of ⁵⁶Ni is generally effective in the ionization and destruction of molecules (in decreasing molecules) as mentioned at the third point above. It is, however, found that, under a certain condition, if the start of the molecule formation lags behind the activation of the decay, the sequence of the formation processes of CO and SiO described below could locally dominate the destructive reactions due to the decay of ⁵⁶Ni (see the third point above); the ionization of atoms by Compton electrons in Equation (17) increases the number density of thermal electrons; H⁻ is produced by the REA reaction, H + e⁻ → H⁻ + γ; H₂ is formed by the AD reaction, H⁻ + H → H₂ + e⁻; OH is formed by the NN reaction, O + H₂ → OH + H; CO and SiO are eventually formed by the NN reactions in Equations (47) and (62), respectively.
18.
In the +Z- and −Z-directions, even at the end of the calculations (∼10,000 days), in the inner regions, the formation of molecules is limited compared with those in the +Y-direction, and the molecules, e.g., CO and SiO, are not so smoothly distributed (see, e.g., the top right panels in Figures 26 and 27). The latter points may be partly attributed to the complicated thermal histories affected by both the heating due to the decay of ⁵⁶Ni and the cooling via CO rovibrational transitions. Another factor for the feature above is that there may be a transition point (radial position) inside which the decay of ⁵⁶Ni increases the amounts of CO and SiO (see the third and seventeenth points); outside the point, it decreases the amounts by destructive processes related to Compton electrons.
19.
The survival of some of the molecules, e.g., O₂ and SO, seems to sensitively be determined by the balance of complicated NN reactions depending on the thermal histories unless the seed atoms for these molecules are sufficient. Actually, the distributions of O₂ are rather complicated in the +Z- and −Z-directions as seen in the models based on the single-star progenitor model (see the top and middle left panels in Figure 27).
20.
From the results of the models based on the binary merger progenitor model b18.3, the calculated distributions of CO and SiO may partly be consistent with the observed 3D features obtained by ALMA (Abellán et al. 2017), e.g., a torus-like structure of CO and holes for both CO and SiO. Such features, however, should be confirmed by the direct application to the 3D models (Ono et al. 2020; Orlando et al. 2020).

In this paper, we limit the application of our methodology only to the one-zone and 1D models. We plan to apply it to the 3D hydrodynamical models for SN 1987A (Ono et al. 2020; Orlando et al. 2020) in order to more directly compare the observed 3D distributions of CO and SiO (Abellán et al. 2017) in the near future. We also have the plan to connect the molecule formation calculations with theoretical models of dust formation (e.g., Nozawa et al. 2010; Nozawa & Kozasa 2013) to investigate the impact of matter mixing on the formation of dust. The observations of not only SN 1987A but also general CCSNe and SNRs (as, for instance, Cas A) by the recently launched JWST will shed light on the molecule and dust formation in CCSNe through the comparison with theoretical models. It would be better to improve our current methodology in particular on the energy depositions by the decay of ⁵⁶Ni and the excitation of CO rovibrational levels to obtain more reliable thermal histories of the ejecta and the CO line emission to be compared with observations in the future.

As for the adopted chemical reactions, the uncertainties of the rate coefficients may not be so small (for a sensitivity study of chemical reaction rates to the formation of CO, see Liljegren et al. 2020) as can be partly seen in Table A1 that some of the rate coefficients have exactly the same values. Therefore, it would be better to use more reliable rate coefficients as much as possible, if available in the future. Additionally, there may be potentially important chemical reactions that are not included in this study (partly due to the limitation of the adopted species) depending on the local initial abundance. Related to this issue and dust formation, theoretically, there have been two approaches to the formation of clusters of dust (nucleation) in CCSNe, i.e., in one approach, dust formation is calculated with a chemical reaction network extended to clusters (Sarangi & Cherchneff 2013, 2015), and in another approach, it is followed by a model of kinematic theories (Nozawa et al. 2010; Nozawa & Kozasa 2013); the differences between the two approaches may affect the consequential dust formation. It would be worth seeing the impact of the differences on dust formation in the two approaches in the future.

Table A1. Adopted Rate Coefficients of the Included Chemical Reactions Other Than Thermal Fragmentation Reactions and Ones Involved with Compton Electrons and UV Photons

Type ^a	Reaction	A_i ^b	α_i ^b	β_i ^b	References
		(cgs)		(K)
`3B`	H + H + M → H₂ + M	6.84 (−33)	−1	0	CD09; NIST
`3B`	H + C + M → CH + M	1 (−33)	0	0	CD09
`3B`	H + O + M → OH + M	4.36 (−32)	−1	0	CD09; NIST
`3B`	Si + H + M → SiH + M	1 (−33)	0	0	CD09
`3B`	C + C + M → C₂ + M	4.25 (−34)	1	0	SL18

Notes. The species "M" in the 3B reactions denotes arbitrary species. The numbers in the parentheses denote the powers of 10. The references represented by the codes in the reference column are as follows. Cherchneff & Dwek (2009; "CD09"); NIST Chemical Kinetic Database (https://kinetics.nist.gov/kinetics/index.jsp; "NIST"); Sluder et al. (2018; "SL18"); UMIST database for astrochemistry 2012 ("UM12"; McElroy et al. 2013); UMIST database for astrochemistry 2006 ("UM06"; Woodall et al. 2007); Cherchneff & Dwek (2010; "CD10"); Giesen et al. (2003; "GI03"); Zhao et al. (2004; "ZH04"); Willacy & Cherchneff (1998; "WC98"); Andreazza & Singh (1997; "AS97"); Franz et al. (2011; "FR11"); Babb & Dalgarno (1995; "BD95"); Andreazza et al. (2009; "AN09"); Andreazza et al. (1995; "AN95"); Andreazza & Marinho (2007; "AM07"); Andreazza & Marinho (2005; "AM05").

^aCode for a reaction type listed in Table 2.^bCoefficients for the Arrhenius form in Equation (4).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Acknowledgments

We acknowledge the anonymous referee for careful reading of the paper and for providing many valuable comments. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. The numerical computations were carried out complementarily on XC40 (YITP, Kyoto University), Cray XC50 (Center for Computational Astrophysics, National Astronomical Observatory of Japan), HOKUSAI (RIKEN). This work is supported by JSPS KAKENHI grant Nos. 19H00693 and JP21K03545. M.O. and S.N. thank the support from RIKEN Interdisciplinary Theoretical and Mathematical Sciences Program. M.O. and S.N. also acknowledge the support from the Pioneering Program of RIKEN for Evolution of Matter in the Universe (r-EMU). M.O. and K.C. thank the support from the Ministry of Science and Technology, Taiwan, under grant Nos. MOST 110-2112-M-001-068-MY3 and MOST 111-2112-M-001-038-MY3 and the Academia Sinica, Taiwan, under grant Nos. AS-CDA-111-M04 and AS-IA-112-M04. S.O. and M.M. acknowledge the financial contribution from the PRIN INAF 2019 grant "From massive stars to supernovae and supernova remnants: driving mass, energy and cosmic rays in our Galaxy."

Appendix A: Chemical Reaction Rates

In this appendix, the adopted chemical reactions that are expressed by the Arrhenius form as in Equation (4), the values of the rate coefficients, and the references of the values are listed in Table A1. As mentioned in the text, other than the reactions listed in Table A1, the thermal fragmentation reactions in Equation (5) in Section 2.1.1, the ionization and destruction by Compton electrons described in Equations (17), (18), and (19), and the dissociation by UV photons in Section 2.1.2 are also taken into account in this paper.

Appendix B: Potential Impact of Nonlocal Energy Deposition due to the Radioactive Decay of ⁵⁶Ni

In this section, the potential impact of nonlocal energy deposition due to the decay of ⁵⁶Ni, which is neglected in the calculations described in the main text (see Sections 2.3 and 2.4), is discussed. As mentioned in Section 2.1.2, in reality, the emitted energies from the decay can be deposited nonlocally over a region corresponding to the length scale of the effective mean free path in Equation (15). Here, a model of such nonlocal energy deposition and the results of the test calculation with the model are presented.

Suppose the emitted energies at the radius $r^{\prime}$ can be deposited at a radius r ( $r\geqslant r^{\prime}$ ) and the energy deposited at r from at $r^{\prime}$ is proportional to $\exp (-(r-r^{\prime} )/{l}_{\gamma ,r^{\prime} })$ , where ${l}_{\gamma ,r^{\prime} }$ is the mean free path in Equation (15) measured at the radius $r^{\prime}$ . It is assumed that nonlocal energy deposition occurs only in the outer radial direction, and the time lag between the emission at $r^{\prime}$ and the deposition at r is neglected. Then, the effective energy generation rate per unit length (ergs per grams per seconds per centimeter), ${q}_{r^{\prime} }(r)$ , at r contributed from at $r^{\prime}$ is expressed as follows

$\begin{eqnarray}{q}_{r^{\prime} }(r)=\left\{\begin{array}{ll}{\epsilon }_{r^{\prime} }\,{N}_{r^{\prime} }\,\exp \left(-\displaystyle \frac{r-r^{\prime} }{{l}_{\gamma ,r^{\prime} }}\right) & \mathrm{if}\ r\geqslant r^{\prime} \\ 0 & \mathrm{if}\ r\lt r^{\prime} ,\end{array}\right.\end{eqnarray} \tag{ B1 }$

where ${\epsilon }_{r^{\prime} }={\epsilon }_{\mathrm{Co}}+{\epsilon }_{\mathrm{Ni}}$ is the summation of the energy generation rates in Equations (9) and (10) at $r^{\prime}$ . ${N}_{r^{\prime} }$ is the normalization factor, and the value is determined so that the equation (energy available nonlocally should be the same as the energy locally generated) below is true

$\begin{eqnarray}&&{\int }_{r^{\prime} }^{\infty }{q}_{r^{\prime} }(r)\,{dr}={\epsilon }_{r^{\prime} }.\end{eqnarray} \tag{ B2 }$

Therefore, ${N}_{r^{\prime} }$ is given as

$\begin{eqnarray}&&{N}_{r^{\prime} }=\displaystyle \frac{1}{{\int }_{r^{\prime} }^{\infty }\exp \left(-\tfrac{r-r^{\prime} }{{l}_{\gamma ,r^{\prime} }}\right){dr}}.\end{eqnarray} \tag{ B3 }$

Finally, the effective energy deposition rate (ergs per grams per second), Q(r), at the radius r is given by integrating the contributions from the radii <r as follows

$\begin{eqnarray}&&Q(r)={\int }_{0}^{r}\displaystyle \frac{{\rho }_{r^{\prime} }}{\rho (r)}{\left(\displaystyle \frac{r^{\prime} }{r}\right)}^{2}{q}_{r^{\prime} }(r)\,{dr}^{\prime} ,\end{eqnarray} \tag{ B4 }$

where ${\rho }_{r^{\prime} }$ and ρ(r) are the density at radii $r^{\prime}$ and r, respectively. The factor of ${({\rho }_{r^{\prime} }/\rho (r))(r^{\prime} /r)}^{2}$ is introduced to conserve the energy flux since the densities and effective surface areas between the radii r and $r^{\prime}$ are different. It is noted that, during the integration in Equation (B4), ${q}_{r^{\prime} }(r)$ is practically set to be zero if $r-r^{\prime} \gt {l}_{\gamma ,r^{\prime} }$ to count the escape of emitted energies from the system to some extent. Hereafter, the results of the test calculation with the model are presented.

To see the impact of the effective nonlocal energy deposition, the fiducial 1D model, i.e., the model b18.3-mean (f_h = 5 × 10⁻³, f_d = 1.0, t_s = 500 days; the top panels in Figure 20 in Section 4.2.1.2), is selected for comparison. With the same values of the three parameters, f_h, f_d, and t_s, and the same angle-averaged 1D profile based on the 3D model b18.3-high (Ono et al. 2020) as the initial condition, the calculation is performed by replacing the effective energy deposition rate, , in Equation (11) with Q(r) in Equation (B4), which may potentially affect the ionization and dissociation of molecules by Compton electrons (see Section 2.1.2) and gas temperatures through Equation (38). It is reminded that, to take into account the efficiency of the heating of the gas (ionization and/or destruction of atoms and molecules), the parameter f_h (f_d) is introduced. The test calculation model is called b18.3-mean-nledep, hereafter.

In Figure 29, the calculation results are shown. First, the differences in the gas temperatures (right panels) are discussed by comparing them with the fiducial model (b18.3-mean). Focusing on the innermost particles, it is recognized that the heating of gas in the model b18.3-mean-nledep is milder compared with the model b18.3-mean; the peaks at about 200 days seen in the model b18.3-mean are not so evident in the model b18.3-mean-nledep. The highest gas temperature at 200 days in the model b18.3-mean-nledep (∼3 × 10⁴ K) is less than that in the model b18.3-mean (∼6 × 10⁴ K). It is understood that part of locally generated energy is deposited in outer regions in the model b18.3-mean-nledep. As seen in the model b18.3-mean, several particles initially at around 1.5 × 10¹³ cm are affected by the cooling through CO rovibrational transitions after 50 days. The innermost particles are also cooled at later phases (≲1000 days). A distinct feature seen only in the model b18.3-mean-nledep is the heating of outer particles (initially at ≳3.5 × 10¹³ cm) with the peaks at around 600 days. As seen in the top left panel in Figure 3, initially, ⁵⁶Ni is distributed inside ∼3.5 × 10¹³ cm. Therefore, the heating of the outer particles is due to the effective nonlocal energy deposition.

Hereafter, the the time evolution of the total amounts of molecules (left panels) is discussed focusing on the differences between the two models. As can be seen, the qualitative features are, actually, rather similar to each other for all the molecules, although quantitatively slight differences can be recognized. In the model b18.3-mean-nledep, the amounts of CO and SiO at around 100 days slightly increase by 20%–30% from those in the model b18.3-mean. After the destruction dominant phase (after approximately 130 days), the destruction of SiO in the model b18.3-mean-nledep is slightly milder than that in the model b18.3-mean. The amount of SiO at 400 days in the model b18.3-mean-nledep is higher than that in the model b18.3-mean by 35%. At the final phase, the amount of SiO in the model b18.3-mean-nledep is higher than that in the model b18.3-mean by 20%. The amount of CO at the final phase is almost the same (the difference is about 5%). As for other molecules, although slight qualitative and quantitative differences in the evolution of the amounts can be recognized, the quantitative differences between the two models seem to be within a few tens of percent. The differences seen between the models may be attributed to the slightly milder gas heating and dissociative processes by Compton electrons (CM reactions, i.e., ionization, dissociation, and dissociative ionization) and by UV photons (UV reactions, i.e., ionization, dissociation, and dissociative ionization) and/or destructive reactions with ionized atomic species, e.g., H⁺ and He⁺.

As a summary, the impact of the effective nonlocal energy deposition at least on the total amounts of the molecules may not be so large compared with other uncertainties, e.g., the parameter value of f_h. The model presented here is, however, rather simple and has not been justified by the comparison with more realistic models based on radiative transfer calculations. In reality, nonlocal energy deposition could be more significant, particularly for local quantities in 3D models. It would be desirable to improve the treatments in the future.

Appendix C: Electron Impact Excitation/Deexcitation Rates

The rate coefficients of the electron impact excitation/deexcitation of vibrational levels in Equation (23), q_ij, are derived from the cross sections in Poparić et al. (2008; see Figure 1 in the reference) by taking a Maxwellian averaging as follows

$\begin{eqnarray}\begin{array}{rcl}{q}_{{ij}}(T) & = & \displaystyle \frac{4}{{\left(2\pi \mu \right)}^{1/2}{\left({k}_{{\rm{B}}}T\right)}^{3/2}}\\ & & \times {\displaystyle \int }_{0}^{\infty }{\sigma }_{{ij}}(E)E\,\exp \left(-\displaystyle \frac{E}{{k}_{{\rm{B}}}T}\right){dE},\end{array}\end{eqnarray} \tag{ C1 }$

where σ_ij(E) is the cross section of the transition of the vibrational state i to j, and μ is the reduced mass of the two-body system of carbon monoxide and electron. The values are plotted in Figure 30. For the transitions involved with the vibrational level v = 6, the cross sections cannot be obtained from Poparić et al. (2008). Therefore, such rates are crudely approximated (assumed) not to violate the qualitative features, for example, among the transitions of the state i = 0 to j ( $j=1,\ldots ,{j}_{\max }$ ), q_0j > q_0k (k > j).

**Figure 30.** Rate coefficients of the electron impact excitation/deexcitation of vibrational levels, q_ij, in Equation (23) as a function of temperature obtained based on Poparić et al. (2008). For example, the coefficient corresponding to the transition of the state i = 0, to j = 1 is denoted as "01" in the legend.
Download figure:
Standard image High-resolution image

How the rates involved with the level v = 6 are obtained is described as follows. First, we assume that the rates, q_6i (i = 0, 1,...,5), are proportional to the rate, q₅₀ (have the same temperature dependence as the rate, q₅₀). As can been seen in Figure 30, the rates, q_i0 (i = 1,...,5), become small as i increases, i.e., the transition energy from v = i to 0 (energy difference between the levels i and 0, i.e., ΔE_i0) increases. Therefore, the rate, q₆₀, is obtained by assuming that the ratio of the rates q₆₀ to q₅₀ is equal to the ratio of the inverse transition energies 1/ΔE₆₀ to 1/ΔE₅₀. The rates q_6i (i = 1,...,5) are obtained by assuming the ratios of the rates q_6i (i = 1,...,5) to q₆₀. The ratio is obtained by an arithmetic averaging of the ratios of the peak rates q_peak,ji to q_peak,j0 (j = 1,...,5). For example, the ratio of the rates q₆₅ to q₆₀ is obtained by the arithmetic averaging of the ratios of the peak rates q_peak,j5 to q_peak,j0 (j = 1, 2, 3, 4). The rates, q_i6 (i = 0,...,5), are obtained by the detailed balance relation between the transitions from v = i to 6, and from v = 6 to i. Although the assumptions here are not so physically validated, since the contributions from rovibrational levels higher than v = 6 are minor, the adopted approximation above should not affect the major results in this study.

Appendix D: Impact of the Reduction Factor f_red for the Escape Probabilities

As noted in Section 2.3, without the reduction factor f_red, at an early phase (less than about 100 days after the explosion), the calculated fluxes of CO vibrational bands tend to be higher than those of the observed peak fluxes. In this section, how the time-dependent reduction factor f_red described in Equation (42) changes the CO line emissions is demonstrated by showing the results with different values of the parameter t_s with representative one-zone calculations.

In Figure 31, the reduction factor f_red is plotted as a function time for different t_s values of 200, 300, and 500 days after the explosion. For example, the values at 100 (300) days are approximately 1.1 × 10⁻², 5.0 × 10⁻², and 1.7 × 10⁻¹ (3.0 × 10⁻², 9.7 × 10⁻², and 2.5 × 10⁻¹) for t_s = 200, 300, and 500 days, respectively. In Figure 32, the one-zone calculation results with f_h = 5 × 10⁻⁴, and f_d = 10⁻² for different t_s values are shown. From the values of the mean escape probabilities $\overline{\beta }$ (green dashed lines) in the left panels, it can be seen that the lower t_s, the lower $\overline{\beta }$ becomes by the reduction factor f_red, in particular, at early phases.

**Figure 32.** One-zone calculation results with f_h = 5 × 10⁻⁴, and f_d = 10⁻² for the comparison of different t_s values. From top to bottom, the results for t_s = 200, 300, 500, and ∞ days are shown, respectively. Left panels: time evolution of the population of vibrational levels of CO, the gas temperature, and the mean escape probability $\overline{\beta }$ . Right panels: time evolution of the fluxes of CO vibrational bands, Δv = 1 (fundamental), Δv = 2 (first overtone), ..., Δv = 6, compared with the observed light curves for Δv = 1 (BD93; W+93; Bouchet & Danziger 1993; Wooden et al. 1993, respectively), and Δv = 2 (M+89, 93; BD93; Meikle et al. 1989, 1993; Bouchet & Danziger 1993, respectively). The top right panel is the same as the middle right panel in Figure 9.
Download figure:
Standard image High-resolution image

Overall, the fluxes of CO vibrational bands (see left panels) tend to become higher than the observed peak fluxes at around 200–250 days (≲10⁻⁸ erg s⁻² cm⁻²) as early as a few tens of days, comparied to the case of t_s = 200 days (the top panels). Since there has been no observation before 100 days, if the fluxes before 100 days were higher than the observed peak fluxes, no detection before 100 days sounds a bit strange. Such early high CO line emission cases also result in rapid gas cooling, and the gas temperatures (violet lines in the left panels) quickly go down to ∼100 K before 100 days except for the case of t_s = 200 days. It is noted that, as partly mentioned in Section 4.1.1.1, the previous studies (Liu et al. 1992; Liu & Dalgarno 1995; Liljegren et al. 2020) have estimated the gas temperature to be >500 K before 1000 days (some of the estimations are based on the fitting of the observed CO spectra). Such rapid cooling is inconsistent with the previous estimations, although it may not necessarily be rejected since the adopted models and assumptions in the previous studies are at least partly different from the ones in this study. The formation of molecules is sensitive to the gas temperature (and density), and the fluxes of the CO line emissions depend on the amount of CO and the gas temperature again. Then, the degree of the gas cooling could affect nonlinearly the fluxes. Though, as can be seen, the lower t_s, roughly the lesser CO line emissions (gas cooling), avoids higher fluxes before 100 days than the observed peak fluxes.

On the other hand, the calculated fluxes with higher t_s, i.e., lower reduction factor f_red, better explain the observed fluxes after 200 days.

As a summary, considering that rapid cooling seen in the calculation results with t_s > 200 days is inconsistent with the previous estimations of the gas temperature, and it may cause a peculiar impact on the formation of molecules at early phases (important phases for the formation of molecules), in this specific case, we prefer t_s = 200 days to avoid the overestimation compared with the observed peak fluxes and such rapid gas cooling, although the calculated fluxes after 200 days are not so consistent with the observations.

It should be noted that the preferred parameter values could be different in other cases, in particular, in the 1D calculations as presented in Section 4.2.

The Impact of Effective Matter Mixing Based on Three-dimensional Hydrodynamical Models on the Molecule Formation in the Ejecta of SN 1987A

Article metrics

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Method

2.1. Chemical Reaction Network

2.1.1. Thermal Fragmentation Reactions

2.1.2. Collisional Ionization and Destruction by Compton Electrons; Dissociation by UV Photons

2.1.3. Initialization of the Chemical Reaction Calculations

2.2. Cooling by Emission from Rovibrational Transitions of CO

2.3. Calculations with the One-zone Approximation

2.4. Calculations with 1D Radial Profiles

3. Strategy and Models

4. Results

4.1. One-zone Calculation Results

4.1.1. Overall Trends of Physical Quantities and the Selection of Reasonable Parameter Sets

4.1.1.1. The Case with fh = 10−4, fd = 5 × 10−2, and ts = 200 days

4.1.1.2. The Case with fh = 5 × 10−4, fd = 10−2, and ts = 200 days

4.1.1.3. The Case with fh = 10−3, fd = 10−2, and ts = 200 days

4.1.1.4. The Criteria for Selecting the Best Parameter Set

4.1.2. Important Chemical Reactions for CO and SiO

4.2. 1D Calculation Results

4.2.1. The Fiducial Case and the Dependence on the Parameters, fh, fd, and ts

4.2.1.1. The Case with fh = 10−2, fd = 10−2, and ts = 500 days

4.2.1.2. The Case with fh = 5 × 10−3, fd = 1.0, and ts = 500 days

4.2.1.3. The Selection of the Best Parameter Set and the Parameter Dependence

4.2.1.4. Important Chemical Reactions for CO and SiO for the Fiducial 1D Calculation Model

4.2.2. The Results including Molecules Other Than CO and SiO with the Angle-averaged 1D Profiles Based on the 3D Models; the Comparison between the Binary Merger and Single-star Progenitor Models

4.2.2.1. The Results of the Model b18.3-mean

4.2.2.2. The Results of the Model n16.3-mean

4.2.3. Comparison with the Spherical Cases

4.2.3.1. The Cases with the Binary Merger Progenitor Model

4.2.3.2. The Cases with the Single-star Progenitor Model

4.2.4. Differences among the Specified Directions

4.2.4.1. The Time Evolution of Quantities in the Models with the Binary Merger Progenitor Model

4.2.4.2. The Time Evolution of Quantities in the Models with the Single-star Progenitor Model

4.2.4.3. Radial Distributions of Molecules in the Models with the Binary Merger Progenitor Model

4.2.4.4. Radial Distributions of Molecules in the Models with the Single-star Progenitor Model

5. Discussion

5.1. Implication to the 3D Distributions of CO and SiO

5.2. Issues of the Calculation Method in This Study

5.3. Implication to the Observations by JWST

5.4. A Caveat: Possible Impacts of Realistic 3D Explosion Models

6. Summary

Acknowledgments

Appendix A: Chemical Reaction Rates

Appendix B: Potential Impact of Nonlocal Energy Deposition due to the Radioactive Decay of 56Ni

Appendix C: Electron Impact Excitation/Deexcitation Rates

Appendix D: Impact of the Reduction Factor fred for the Escape Probabilities

Footnotes

4.1.1.1. The Case with f_h = 10⁻⁴, f_d = 5 × 10⁻², and t_s = 200 days

4.1.1.2. The Case with f_h = 5 × 10⁻⁴, f_d = 10⁻², and t_s = 200 days

4.1.1.3. The Case with f_h = 10⁻³, f_d = 10⁻², and t_s = 200 days

4.2.1. The Fiducial Case and the Dependence on the Parameters, f_h, f_d, and t_s

4.2.1.1. The Case with f_h = 10⁻², f_d = 10⁻², and t_s = 500 days

4.2.1.2. The Case with f_h = 5 × 10⁻³, f_d = 1.0, and t_s = 500 days

Appendix B: Potential Impact of Nonlocal Energy Deposition due to the Radioactive Decay of ⁵⁶Ni

Appendix D: Impact of the Reduction Factor f_red for the Escape Probabilities