© The Authors 2023

1 Introduction

When massive stars (M_ZAMS ≳ 8 M_⊙) reach the end of their lives, they explode as core-collapse supernovae (CCSNe). These events enrich the cosmos in elements produced in both hydrostatic and explosive nucleosynthesis, and leave behind exotic remnants such as black holes and neutron stars. Likely part of the SNe ejecta condenses into different dust species after a few years or decades and may define a major source of dust in the Universe (see Sarangi et al. 2018, for a review). Understanding galactic chemical evolution of elements as well as dust then relies on understanding this element ejection.

An open question in SN research is the stellar origin of Type Ibc SNe, collectively named stripped-envelope SNe (SESNe). This is a subtype of CCSNe which lacks the spectral signatures of either hydrogen (Type Ib) or both hydrogen and helium (Type Ic). The indication is then that the outer hydrogen envelope, or both the hydrogen and helium envelopes, have been stripped away during the star’s evolution. The mechanism for this is, however, not well understood. Two main scenarios have been proposed. The progenitors of Type Ibc SNe may be massive single stars, M_ZAMS ≳ 25 M_⊙, that lose their outer layers due to strong stellar winds (Conti 1975). The progenitors may also be lower-mass stars, M_ZAMS ≲ 20–25 M_⊙, that have their envelopes stripped through interaction with a binary companion (Nomoto et al. 1995). Over the years, results from light curve modeling (e.g., Ensman & Woosley 1988; Woosley et al. 1994; Bersten et al. 2014; Dessart et al. 2020; Ertl et al. 2020), progenitor studies (Smartt 2015; Eldridge & Maund 2016; Yoon et al. 2017), and rate constraints (Smith et al. 2011) have indicated that the second channel is likely significant. The emergence of the transitional IIb type as a common SN class (rate of ~ 10% of all CCSNe, e.g., Shivvers et al. 2017) has further strengthened this. The posed question has therefore morphed toward trying to answer what contribution the single-star channel makes, and what different outcomes are possible for massive Wolf-Rayet stars at the end of their evolution.

To determine the detailed properties and origin of an SN, spectral studies are necessary. In particular, once the photosphere has receded after a few months, observations reveal the inner region of the exploded star, including its nucleosynthesis. Studies in this so-called nebular phase are, however, challenging. Although sustained by radioactivity, the SN continuously fades, limiting the epochs over which it can be observed and the signal-to-noise ratio (S/N) of the collected spectra. The physics of the spectral formation is complex, involving fast differential expansion, a radioactive environment with strong nonlocal thermal equilibrium (NLTE) effects, and energy transfer from MeV particles down to eV or sub-eV gas (see Jerkstrand 2017, for a recent review).

To complicate things further, at some point the SN will start to form molecules, and later, dust. Direct evidence for this comes from observations of molecular and dust emission at IR wavelengths, and also indirectly by dust obscuration of the optical atomic emission. Ground-based spectral data generally extend up to ~2.5μm, which includes the first overtone emission of CO, but no other molecular emission. For Type Ibc SNe, a handful of such observations have been published (Gerardy et al. 2002; Hunter et al. 2009; Rho et al. 2021). Observations of the emission from other molecules requires MIR spectroscopy of the 2-16 μm region where the bulk of molecular emission lies. MIR observations of SNe are typically not possible from the ground, with the single exception of the very close-by Type II SN 1987A (see e.g., Aitken et al. 1988; Roche et al. 1991). The Spitzer Space Telescope was used to study a few nearby (< 10 Mpc) Type II SNe (see e.g., Kotak et al. 2006; Kotak 2009). In both the case of SN1987A and in these Type II SNe observed with Spitzer, molecular emissions from CO and SiO was identified. For SESNe, there are so far no MIR spectral observations, although the James Webb Space Telescope holds promise to change this.

Regarding modeling, there are studies focused on the formation processes of molecules (e.g, Petuchowski et al. 1989; Lepp et al. 1990; Liu et al. 1992; Liu & Dalgarno 1995, Liu & Dalgarno 1996; Gearhart et al. 1999), or both molecules and dust (Clayton et al. 1999, 2001; Cherchneff & Dwek 2009, 2010; Biscaro & Cherchneff 2014, 2016; Sluder et al. 2018). These studies have provided valuable insights into the chemistry of SNe and produced quantitative results for molecule and dust masses for different SN types. These calculations have, however, been limited to using parameterized temperatures and have not performed radiative transfer to make predictions for the molecular emission. In radiative transfer modeling predicting SN spectra, on the other hand, emission and cooling by molecules have been ignored or treated in a parameterized fashion (e.g., Dessart et al. 2013; Jerkstrand et al. 2012, 2014).

The coupling of molecular chemistry with spectral synthesis is therefore so far a missing ingredient on the modeling side. Both from the observed strength of the molecular emission (e.g., Meikle et al. 1989; Kotak 2009), and the only two studies calculating molecular cooling so far (Liu & Dalgarno 1995; Liljegren et al. 2020) molecules appear capable of cooling the ejecta significantly, perhaps by several thousand degrees. This would strongly affect the optical spectrum which for many SNe is the only spectral region observed. Because there is an intricate feedback loop between molecule formation and temperature (Liljegren et al. 2020), the cooling and formation processes need to be modeled self-consistently and simultaneously.

To address this need, we here present the first SN spectral synthesis models where detailed prescriptions of the molecular formation and cooling processes are coupled to the radiative transfer. This is a continuation of the work presented in Liljegren et al. (2020), where the methodology foundation was outlined. Here we extend on this work by including molecular formation in multizone models, adding silicon and sulfur-bearing molecules in addition to carbon-bearing ones, and calculating the molecular destruction by Compton electrons using self-consistent methods. In contrast to L20, where a highly simplified test model was presented to investigate different relevant processes, the models here are meant to represent realistic SN spectra, including the molecular emission in the IR region. Such simulations can consequently be used both to interpret available optical and near-infrared (NIR) observations and future mid-infrared (MIR) observations with, for example, JWST and ELT.

To explore the impact and diagnostic potential of molecule formation in Type Ibc SNe, we here start with a case study of a Type Ic model in the 100–600 days interval. With this model, we investigate the type of molecules that form, at what epochs, and in which quantities. We study the influence of molecule formation on the temperature and the resulting spectra throughout the UVOIR range. These results can be used to predict what is possible to observe with upcoming instrumentation, and what specific spectral regions and phases are of most interest to collect data for.

This paper is structured as follows: Sect. 2 gives an overview of the modeling methods used. Section 3 presents the molecular mass results. Section 4 discusses the cooling effects of the molecules. Section 5 presents the Compton electron destruction rates, and Sect. 6 is dedicated to a detailed description of the resulting IR spectra for the standard model. In Sect. 7, an investigation on how line shapes depend on physical conditions is presented. In Sect. 8, we make a parameter space investigation, varying the density of the standard model. Section 9 discusses how the optical spectra are influenced by these changes in model density, and in Sect. 10 we compare our results with observations. In Sect. 11, we investigate the observational potential with JWST. Summary and conclusions are presented in Sect. 12.

2 Modeling methods

We have implemented a treatment of the formation and NLTE cooling of molecules in the SN spectral synthesis code SUMO (Jerkstrand et al. 2011, 2012). The physical modeling approach is described in this section, while a detailed description of the model molecules is available in Appendix A. Appendix D contains the molecular reactions rates used.

2.1 SUMO

SUMO is an NLTE spectral synthesis code, which employs detailed physical treatments of all key steps in the spectral formation (Jerkstrand et al. 2011, 2012; Jerkstrand 2017). It takes a hydrodynamic supernova ejecta model as input, with density and composition as functions of velocity, and calculates the physical conditions and emergent spectrum by solving for the temperature, the statistical equilibrium of ion abundance and level populations, and the radiation field. SUMO is specialized to epochs post the diffusion phase and can provide both snapshot steady-state solutions at specified phases as well as fully time-dependent solutions (Pognan et al. 2022). Currently, the atomic treatment considers 22 elements between hydrogen and nickel, and a few r-process elements, each with multiple ionic species.

2.2 Molecules

Molecular species are divided into two groups. In the first group, full inclusion, the species is included both in the formation processes and in the temperature and radiative transfer calculations (through NLTE solutions for its levels). In the second group, partial inclusion, it is included only in the formation process calculations. This division is done to avoid unnecessary, expensive NLTE calculations for rare molecules (or dipole-free ones, e.g., S2) which have no impact on thermal conditions. The included molecules and their categories are listed in Table 1.

The molecular abundances are calculated using a chemical network approach. For full inclusion group species, the impact on the temperature and spectrum are then solved for by treating the molecules on the same footing as the atoms, with complete coupling to temperature and radiation field. We treat transfer in the rovibrational lines in the Sobolev limit. The molecular electronic UV excitations are, however, not included - at nebular times the temperatures are too low for these to play any thermal role. Four molecular species; carbon monoxide (CO), silicon monoxide (SiO), silicon monosulfide (SiS), and sulfur monoxide (SO), receive this full treatment. We choose these because they are the most common and abundant species formed in the SN ejecta.

For the molecules in the partial inclusion group, we calculate the molecular abundances in the same chemical network calculation as previously mentioned molecules. In other parts of SUMO, however, the partially included molecules act as passive components, that is to say they do not influence either the temperature or spectrum, only acting as storage sinks for the component atoms.

2.2.1 Chemical network

To describe molecule formation and destruction, and calculate the molecular abundances, we use a chemical kinetic network approach. We here make an overview of the methodology. For an in-depth description, details and discussion see Liljegren et al. (2020).

The conditions under which molecules form in supernovae are harsh; as the ejecta moves at a few percent of the speed of light the material quickly moves into a low-density regime, and both the population of Compton electrons and the strong radiation field can ionize or dissociate the molecules that do form.

Whereas in Liljegren et al. (2020) the chemistry was modeled in the steady state approximation, we here employ a time-dependent formalism (both for temperature and the molecular and atomic abundances), solving (e.g., Cherchneff et al. 1992) $$\frac{\text{d}{n}_i}{\text{d}t}+\frac{3n_i}{t}=C_i-n_i{D}_i,$$(1)

where C_i is the total volumetric creation rate of particles of type i (cm⁻³ s⁻¹) and D_i is the total destruction rate per particle (s⁻¹). We also treat the temperature equation in its time-dependent form (Kozma & Fransson 1998; Jerkstrand 2011; Pognan et al. 2022).

By comparing solutions in steady state versus time-dependent mode, we could assess the differences due to these treatments. We found only quite small differences in general. The only significant difference was the CO abundance curves at the latest epochs. As the O⁺ charge transfer reaction switches off at a critical temperature, the destruction rate plummets by a factor 10–100. In steady state, this leads to an immediate increase of the CO abundance by the same factor. However, the formation time scale is several days at the low prevailing densities at 500–600 days, so in time-dependent model it takes 100–200 days to reach the new equilibrium with a factor 10–100 higher abundance.

The importance of time-dependent terms in a particular rate equation depends on the total formation and destruction rates for that species. Parts of a network may be affected but not others. For the main molecular species here, we find that, apart from the late CO phase described above, these total rates are always high compared to the rates of change of density and radioactive power; thus the steady state mostly holds over the modeled time period.

We include several different thermal and nonthermal processes, an overview of such reactions relevant to molecular formation is given in Table 2. The rates for these processes depend on the number density of the reactants, and the efficiency of the reactions, described by a temperature-dependent rate coefficient k(T) (except for destruction by Compton electrons, which is discussed in Sect. 2.2.3).

The rate coefficients for the thermal and recombination reactions are here expressed in a modified Arrhenius form as: $$k\left(T\right)=\alpha \times {\left(\frac{T}{300\text{\hspace{0.17em}K}}\right)}^{\beta }\times \exp \left(-\gamma /T\right),$$(2)

where α, β, and γ are parameters specific to each reaction. The unit of k(T) depends on the number of reactants, being s⁻¹ for unimolecular reactions and cm³ s⁻¹ for bimolecular reactions. The reaction database in this work is an extension of the one used in Liljegren et al. (2020), now including also sulfur and silicon-based species and reactions. The primary sources are the UMIST Database for Astrochemistry¹ (McElroy et al. 2013), the online database KIDA² (Wakelam et al. 2012) and the chemical kinetic database NIST³ (Manion et al. 2015). Some reactions, for example the radiative association C + O → CO (Gustafsson & Nyman 2015), were updated with newer rates. Our reaction database was also supplemented by reaction rates from Cherchneff & Dwek (2009) and Sluder et al. (2018). Readers can refer to Appendix D for a complete reference list.

2.2.2 Ionization energies

For the full inclusion in the SUMO framework, the ionization energies of the included molecular species are needed. These were available for CO, CO⁺, SiO, SiS and SO, but had to be calculated for SiO⁺, SiS⁺ and SO⁺. For details see Appendix C.

2.2.3 Dissociation by Compton electrons

Gamma rays, produced by the radioactive decay in SNe, will repeatedly Compton scatter on free and bound electrons to create a population of high-energy Compton electrons (referred to as ${\text{e}}_{\text{C}}^{-}$ in this work). Collisions between these Compton electrons and molecules is a significant and sometimes dominant destruction process for molecules. For a diatomic molecule AB the dominant channel is typically ionization: $$\text{AB\hspace{0.17em}}+{\text{e}}_{\text{C}}^{-}\to {\text{AB}}^{+}+{\text{e}}^{-}+{\text{e}}_{\text{C}}^{-},$$(3)

followed by the weaker dissociation channel: $$\text{AB\hspace{0.17em}}+{\text{e}}_{\text{C}}^{-}\to \text{A\hspace{0.17em}}+\text{\hspace{0.17em}B}+{\text{e}}_{\text{C}}^{-}.$$(4)

Other reactions (e.g., $\text{AB\hspace{0.17em}}+{\text{e}}_{\text{C}}^{-}\to \text{A}+{\text{B}}^{+}+{\text{e}}^{-}+{\text{e}}_{\text{C}}^{-}$) are less important and disregarded in this work. In previous works, for example Liu & Dalgarno (1995), Cherchneff & Dwek (2009), Sluder et al. (2018), Liljegren et al. (2020), such rates were estimated based on work by Liu & Dalgarno (1995), where the destruction of CO was investigated and presented as a relationship between the (zone) power L, the number of particles N_tot and the mean energy per ion pair W (or molecule pair in the case of dissociation): $$k_C=\frac{L}{N_{\text{tot}}W}.$$(5)

For a fixed composition (49.5% 01,49.5% C I, 1% CO I), Liu & Dalgarno (1995) calculated W = 34 eV for CO ionization. Here we perform direct calculations of the Compton destruction rates for the reactions (3), for the most abundant molecular species. This is possible as SUMO contains a Compton electron distribution solver (see Kozma & Fransson 1992, for details), which combined with electron impact cross-sections for the relevant molecules, can yield the destruction rates for these channels. For reaction (3), the cross sections used are listed in Appendix B. For reaction (4), cross-section are available for CO, but not for the other molecules (to the authors’ knowledge). We estimate the rate of reaction (4) for SiO, SiS and SO by scaling the calculated rates for reaction (3) with the same factor that differentiates the CO rates for reactions (4) and (3) (which is ~0.2). In Sect. 5, we compare our calculated rates with those estimated by Liu & Dalgarno (1995).

2.2.4 Photoionization

In this work, photoionization is included for the atoms and atomic ions as usual in SUMO, however, the process is ignored for the molecular species. Previous work indicates that this is a subdominant process for molecules (Cherchneff & Dwek 2009). We also found this to be the case, when experimenting with a hydrogenic approximation for the relevant cross-sections of the molecules. We hope to include molecular photoionization in subsequent work, however, some necessary cross-section data is still not available.

2.3 Molecular spectra and cooling

Molecular emission has been observed to be significant in the near-to-mid IR region in Type II SNe, and may be the case for Type Ibc SNe as well, if molecules form in a considerable amount. Additionally, it is known from previous work (Liu & Dalgarno 1995; Liljegren et al. 2020) that molecules can add significant cooling channels, if abundant enough, lowering the temperature by several thousand degrees. We model the NLTE (i.e., not assuming local thermodynamical equilibrium, LTE) populations and cooling for four molecular species: CO, SiO, SiS and SO.

The emission in the IR region comes from ro-vibrational lines, that is from transitions between rotational levels in one vibrational state to another. For the investigated physical conditions, the rotational level populations within a vibrational level remain close to thermal equilibrium (see discussion in Liu et al. 1992); the radiative lifetimes of such levels are long compared to the time-scales of pure rotational collisions with either neutral atoms or electrons. Vibrational populations are, however, not in LTE, especially not at the later times (see discussion in Liu & Dalgarno 1994). We therefore calculate the ro-vibrational population by solving the equations of statistical equilibrium, and ensure relative LTE within a vibrational level by putting the collision strengths for the pure rotational transitions to large values.

For each level with vibrational number v and rotational number j, we have (somewhat simplified, ignoring exchange with other ions through ionizations and/or recombinations): $$n_{v,j}{\displaystyle \sum _{v\prime ,j\prime }{P}_{v,j\to v\prime ,j\prime }}={\displaystyle \sum _{v\prime ,j\prime }{n}_{v\prime ,j\prime }{P}_{v\prime ,j\prime \to v,j},}$$(6)

where the left-hand side represents the outgoing rate from the level n_v,j, and the right-hand side is the incoming rate from all levels $n_{{\nu }^{\prime },j^{\prime }}{}_{\cdot }{P}_{\nu ,j{\to }_{{\nu }^{\prime },j^{\prime }}}$ is the transition probability, which depends on the Einstein coefficients A and the Sobolev escape probability β for spontanous emission, the radiative photoexci-tation/deexcitation rates R, and the collision-induced transitions C as: $$\begin{array}{cc}{P}_{v,j\to v\prime ,j\prime }=A{\beta }_{v,j\to v\prime ,j\prime }+R_{v,j\to v\prime ,j\prime }+C_{v,j\to v\prime ,j\prime }& ,v>v\prime \end{array}$$(7) $$\begin{array}{cc}=R_{v,j\to v\prime ,j\prime }+C_{v,j\to v\prime ,j\prime }& ,v\prime >v\end{array}.$$(8)

Equation (7) represents the de-excitation probability, and Eq. (8) represents the excitation probability, for depopulating level v, j.

For these calculations, we adopt the energy levels and A-values from Li et al. (2015) for CO, from Barton et al. (2013) for SiO and from Upadhyay et al. (2018) for SiS. ₄his data is available online through the Exomol collaboration⁴. The line list for SO was received through private communication with Ryan Brady and Sergei N. Yurchenko, based on in-prep work to be published (Brady et al. 2022).

The collision-induced transitions C are here due to collisions between molecules and thermal electrons, and therefore depend on the electron density n_e as $$C=q\left(T\right)\times n_e,$$(9)

where q(T) is the transition rate coefficient, which depends on the electron temperature. We use data from Ristić et al. (2007) for CO. To the authors’ knowledge, equivalent data is not available for any of the other molecular species investigated here. For our NLTE molecules, SiO, SiS and SO, we therefore used the CO data as well, for lack of a better alternative. For details on how we extrapolate rate coefficients for transitions without available data, and details about the model molecules, see Appendix A.

We included levels up to the vibrational quantum number v = 7 and rotational quantum number j = 80. Through experimentation we found that higher levels had no significant impact on the temperature or spectra. We include ro-vibrational radiative transitions up to Δv = 5.

2.4 The supernova model

To represent a typical Type Ic SN we define a standard model based on a M_ZAMS = 25 M_⊙ progenitor core-collapse supernova from Woosley & Heger (2007), with the hydrogen and helium layers removed. The progenitor has a CO core of 6.8 M_⊙, of which 1.8 M_⊙ forms a compact object and 5.0 M_⊙ is ejected in the supernova. As in Jerkstrand et al. (2014, 2015), we divide the ejecta into five main compositional zones - Fe/He, Si/S, O/Si/S, O/Ne/Mg and O/C, further specified in Table 3 This will be referred to as the standard model throughout this work. By choice of the inner zone mass the ⁵⁶Ni mass in the model is 0.06 M_⊙. The zones are macroscopically (but not microscopically) mixed throughout a region between 0 and 3500 km s⁻¹, following the treatment in Jerkstrand et al. (2015). We assume a distance of 10 Mpc for the model.

There are no robustly determined major differences between Type Ic, Ib and IIb ejecta. The thin hydrogen layer in Type IIb SNe has no influence on any of the inner ejecta properties or spectra (Jerkstrand et al. 2015). Type Ib/IIb SNe seem to have similar ejecta masses and velocities as Ic SNe (e.g., Prentice et al. 2019). As long as the helium in these SNe does not mix into the metal zones and quench the chemistry by He+ reactions, we therefore expect to see little differences in models with our without a He layer, and the model here should be representative for Ib/IIb SNe as well, apart from He and H lines.

Further following the method in Jerkstrand et al. (2015), we assume a uniform density for the three oxygen zone (O/Si/S, O/Ne/Mg, and O/C zones), with a lower density for the two innermost zones (Fe/He and Si/S). This differentiation is based on expansion of the ⁵⁶Ni-rich zones during the first days due to trapped radioactive heating. The density structure of the model is set to obey $${\rho }_{rel}=\left(1-10-\chi -\chi -\chi \right),$$(10)

where the numbers/symbols on the right hand side represent the five zones so that the density of for example zone 4 is ρ₄ = χ×ρ₁. From previous work on nebular spectra of stripped-envelope SNe values of χ ~ 30–200 have been indicated (Jerkstrand et al. 2015), although the number is not very well constrained, see for example discussion in Dessart et al. (2021). For the standard model we use χ=100, making the density of the O/Si/S, O/Ne/Mg and O/C zones 100 times larger than that of the Fe/He zone. This clumping structure can equivalently be described by the density of the zone relative to the constant density in a model where all the ejecta is distributed uniformly throughout the expansion volume. Relative to the density in such uniform model, χ = 100 corresponds to an overdensity of 2.8 for the O-zones, an under-density of 0.026 for the Fe/He zone, and an underdensity of 0.23 for the Si/S zone. In Sect. 8 we make a small parameter space exploration for different χ values. The lowest χ model there corresponds to an O-zone overdensity of 1.8, and the largest to an overdensity of 14.

From the density structure the filling factors, which is the measure of clumping used in SUMO, can be calculated. The filling factor of a given zone is defined as the fraction of the expansion volume that the zone occupies, and the sum of all filling factors f_i must equal one $${\displaystyle \sum _1^{N_{zone}}{f}_i=1}.$$(11)

Additionally, the volume and filling factors of two zones can be related using the density structure as $$\frac{M_j}{f_j\times {\rho }_j}=\frac{M_{j+1}}{{\rho }_{j+1}\times f_{j+1}},j\in \left(1,N_{\text{zone}}-1\right),$$(12) $$\begin{array}{l}\to f_i-\left(\frac{{\rho }_{j+1}}{{\rho }_j}\frac{M_j}{M_{j+1}}\right)f_{j+1}=\hfill \\ f_i-\left(\frac{{\chi }_{j+1}}{{\chi }_j}\frac{M_j}{M_{j+1}}\right)f_{j+1}=0.\hfill \end{array}$$(13)

This forms a system of linear equations that is solved for f_i. The f_i values for the standard model are available in Table 3.

3 Molecular formation results

We model the supernova between 100 and 600 days after explosion, at time intervals of 50 days. Earlier epochs than 100 days are difficult and time-consuming to converge with full NLTE atomic and molecular physics. Such epochs also still may have diffusion effects taking place, and relatively high opacities makes emergent diagnostics complex to interpret. Later than 600 days virtually no observations of stripped-envelope SNe exist, and issues with dust formation (which we do not model) become pressing.

Figure 1 shows the time evolution of the amount of molecules formed in the standard model. The five most abundant molecules to form are carbon monoxide (CO), silicon monoxide (SiO), silicon monosulfide (SiS), sulfur monoxide (SO), and disulfur (S₂). Of these, CO, SiO, and SO increase with time and their highest masses are reached at the last epoch of 600 days. SiS and S₂, on the other hand, have their highest masses early on and decline at later times.

For all molecules, the abundances are relatively constant between 100 and 300 days; the model gives ~10⁻³ M_⊙ for SiO, ~10⁻⁴ M_⊙ for S₂, CO and SiS, and ~10⁻⁷–10⁻⁶ M_⊙ for SO. After some rapid changes around 330 days, the abundances then again settle, now with significantly more SO but less SiS and S₂. At the last epoch an acceleration in the growth of CO is seen. The final masses at 600 days are 2 × 10⁻³ M_⊙ of CO, 3 × 10⁻³ M_⊙ of SiO, 7 × 10⁻⁴ M_⊙ of SO, and 3 × 10⁻⁵ M_⊙ of SiS and S₂. The fraction of atoms locked in molecules for the relevant species at this final time from a few per-mile for S (0.2%), Si (0.4%) and O (0.5%) to a percent for C (1%). This is a smaller fraction found than by for example Sarangi & Cherchneff (2013), however, the amount of molecules formed partly depends on the assumed density of a model and will therefore vary for different types of SNe (previously investigated in Biscaro & Cherchneff (2014) and discussed here in Sect. 8).

Other molecules are much less abundant; about 10⁻⁷ M_⊙ of O₂ and 10⁻⁸ M_⊙ of C₂ forms, both constant in time, and other included species form less than 10⁻⁸ M_⊙.

The reactions governing the formation and destruction balance of each molecule, and their evolution over the simulated times, are discussed below.

3.1 Carbon monoxide

Initially the total mass of CO is ~10⁻⁴ M_⊙, increasing to ~2×10⁻² M_⊙ after 450 days. CO predominantly forms in the O/C zone, with a small amount (always less than 10% of the total) forming in the O/Ne/Mg zone (Fig. 2).

As seen in Fig. E.2, CO is primarily formed by radiative association, $$\text{C\hspace{0.17em}}+\text{\hspace{0.17em}O}\to \text{CO\hspace{0.17em}}+\text{\hspace{0.17em}h}v.$$(14)

Other prominent creation channels are by recombination $${\text{CO}}^{+}+{\text{e}}^{-}\to \text{CO\hspace{0.17em}}+\text{\hspace{0.17em}h}v,$$(15)

through neutral exchange reaction with C₂ $${\text{C}}_2+\text{O}\to \text{C\hspace{0.17em}}+\text{\hspace{0.17em}CO},$$(16)

and through charge exchange between O and CO⁺ $${\text{CO}}^{+}+\text{O}\to \text{CO\hspace{0.17em}}+{\text{O}}^{+}.$$(17)

Up until 450 days a major destruction pathway is through the charge exchange with O⁺, $$\text{CO\hspace{0.17em}}+{\text{O}}^{+}\to {\text{CO}}^{+}+\text{O.}$$(18)

This reaction is endothermic and turns off when the gas gets too cool, consequently resulting in a growth of the CO abundance after 550 days once this destruction channel is quenched. This is discussed further in Liljegren et al. (2020).

Other important destruction channels are collisions with Compton electrons, $$\text{CO\hspace{0.17em}}+{\text{e}}_{\text{C}}^{-}\to {\text{CO}}^{+}+{\text{e}}^{-}\text{or\hspace{0.17em}C\hspace{0.17em}}+\text{\hspace{0.17em}O\hspace{0.17em}}+{\text{e}}_{\text{C}}^{-},$$(19)

and through the reverse reaction of (16). These results broadly agree with the findings in Liljegren et al. (2020), where CO formation was investigated in a Type IIP supernova using a smaller network model.

Fig. 1 Total mass formed for different molecular species as function of time. Solid lines are neutrals and dashed lines ionized molecules.

3.2 Silicon monoxide

The amount of SiO that forms in the standard model is ~10⁻³ M_⊙ throughout the simulated time interval. As seen in Fig. 2, SiO is almost exclusively formed in the O/Si/S zone. The SiO mass in the O/Ne/Mg zone does grow with time, ending at almost 10⁻⁴ M_⊙ at 600 days, but this still a small fraction of the total mass.

The main production channel is at all epochs (Fig. E.1) through radiative association, $$\text{Si\hspace{0.17em}}+\text{\hspace{0.17em}O}\to \text{SiO\hspace{0.17em}}+\text{\hspace{0.17em}h}v.$$(20)

Other important channels are through neutral exchange with O₂, $${\text{O}}_2+\text{Si}\to \text{SiO\hspace{0.17em}}+\text{\hspace{0.17em}O,}$$(21)

and through recombination, $${\text{SiO}}^{+}+{\text{e}}^{-}\to \text{SiO\hspace{0.17em}}++\text{\hspace{0.17em}h}v.$$(22)

The major destruction reactions are Compton electron destruction and the reverse reaction of Eq. (21). While the SiO mass is mostly stable over the simulated times, it is still connected to the behavior of the sulfur-bearing species in the O/Si/S zone (discussed in Sect. 3.3); SO and SiS.

3.3 Silicon monosulfide, sulfur monoxide and disulfur

There is an intricate relationship between the sulfur-bearing molecular species SiS, SO, and S₂, which makes their behavior complex and temperature-dependent. This is illustrated in Fig. 3, which contains the part of the chemical network that is relevant to this connection. Each of the species can form through radiative association and through recombination. There are strong neutral exchange reactions which connect these three molecules to each other. SO is connected to S₂ by $$\text{SO\hspace{0.17em}}+\text{\hspace{0.17em}S}\rightleftharpoons {\text{S}}_2+\text{O,}$$(23)

and, in turn, S₂ to SiS through $${\text{S}}_2+\text{Si}\rightleftharpoons \text{SiS\hspace{0.17em}}+\text{\hspace{0.17em}S.}$$(24)

Initially significant amounts of S₂ and SiS form, (1 – 3) × 10⁻⁴ M_⊙, while SO is much less abundant (around 5 × 10⁻⁷ M_⊙). After 300 days there is a swift decline of the molecular masses of S₂ and SiS, simultaneously as the SO mass increases to ~2 × 10⁻⁴ M_⊙, and then slowly continues to increase for the rest of the evolution. This is due to reaction (23), marked red in Fig. 3, being endothermic (in the S₂ forming direction) and becoming inefficient at low temperatures. A temperature drop from 3000 K to 2000 K (the relevant temperatures at around 300 days) will result in a decrease of an order of magnitude for this reaction rate. For the reverse reaction of (23), marked in green in Fig. 3, data for the rate is only available for lower temperatures. This pathway between SO and S2 is then effectively quenched, while the pathway between S2 and SO turns on, when the gas cools at later times. The overall effect is that the net flow from S₂ to SO becomes larger. The consequence is then a higher SO abundance at later phases, as the other involved reaction rates are not very sensitive to temperature, and a lower S₂ and SiS abundance, as more of it is turned into SO.

The molecule SO is connected to O₂ by $${\text{O}}_2+\text{S}\rightleftharpoons \text{SO\hspace{0.17em}}+\text{\hspace{0.17em}O.}$$(25)

This then further connects the behavior of the sulfur species to the SiO abundances through reaction (21) and its reverse reaction. The rates, and therefore the importance, of these reactions are however lower than the reactions (23) and (24) previously discussed. An important destruction process for all three molecules is collisions with Compton electrons.

Initially, SiS and S₂ are formed primarily in the O/Si/S zone. After 300 days, when the production of SO becomes much more efficient, these species are instead mostly formed in the Si/S zone: the lack of oxygen in this zone prevents the previously discussed effects involving O and O₂ that quenches SiS and S₂ in the O/Si/S zone. SO forms almost exclusively in the O/Si/S zone at all simulated times.

Fig. 2 Masses of the four molecular species included in the cooling and radiative transfer calculations as function of time. The different lines show the different masses in each zone.

Fig. 3 Graphical representation of the part of the chemical network that connects the sulfur-bearing molecular species.

3.4 Comparison to previous results

These are the first model calculations for molecule formation in a Type Ic supernova. It is of interest to compare the results with earlier models for Type II SNe. Because both ejecta structures and methodologies differ between works, it is not straightforward to identify the specific driving factors for differences, but a discussion can be had on this. We focus the comparison on CO and SiO.

In the Type II model of Liljegren et al. (2020), the CO abundance had an initial rapid growth phase, then settled at a quasi-constant abundance after 200–300 days. It was demonstrated that both the turn-off time and the settling value were sensitive to the zone density and the amount of radioactive power absorbed. For the standard model, the settling value was a few times 10⁻³ M_⊙.

The density in the standard model here is about a factor of 10 lower, for a fixed time, due to the higher expansion velocities in a Type Ic SN. This means that our first epoch of 100 days corresponds, density-wise, to a bit more than 200 days (ρ ∝ t⁻³)in a Type II SN and we are already on the “settled” plateau here by 100 days. The CO masses here are lower, 10⁻⁴–10⁻³ M_⊙ (Fig. 1) – the lower densities lead to less efficient CO cooling which, as analyzed in L20, in turn, leads to less CO. The formation curves here are qualitatively similar to the “no cooling” variants in L20.

In the Type II model of a M_ZAMS = 20 M_⊙ star in Cherchneff & Dwek (2009), CO forms rapidly and then settles on a quasi-constant value after about 200 days. The settling value, a few times 10⁻¹ M_⊙, is quite a bit higher than in the L20 model, due to some combination of different ejecta models, chemical network, and temperature evolution. Only at significantly later times (≳500d in the standard Type Ic model here, which would correspond to ≳1000 days in a Type II model) are large CO masses formed in our models. In the M_ZAMS = 15 M_⊙ Type II model of Sarangi & Cherchneff (2013), the CO forms over a longer time-scale, settling at a similar final value as in Cherchneff & Dwek (2009) but only after about 600–700 days.

For SiO, both Cherchneff & Dwek (2009) and Sarangi & Cherchneff (2013) obtain an abundance peak quite early on (150–200 days) and then a monotonic decrease, having almost completely vanished by 500 days. The Si and O instead get locked up in O₂, SO, and silicate dust precursors. In our Type Ic model the SiO abundance is declining already from the first epoch of 100 days. It then starts to build up again by about 350 days, together with SO, as SiS rapidly declines. The chemistry is complex, with SiS, SO, and O₂ all playing roles for the SiO abundance – and in addition, dust precursors which are included in the Sarangi & Cherchneff (2013) model but not in ours. Obtaining observations of the time evolution of SiO, SiS, SO, and silicate dust emission with JWST would provide important tests and diagnostics of the silicate chemistry of the supernova ejecta.

Fig. 4 Temperature evolution for the standard model, with molecules included (dashed lines), and excluded (solid lines), for the different zones.

4 Molecular cooling effects

Figure 4 shows the time evolution of the temperature in the different zones, for the standard model with molecular cooling included (dashed lines) and excluded (solid lines). It is immediately clear that our standard model suggests that temperatures in Type Ic SNe are, with a few exceptions, not significantly affected by molecular cooling. The strongest effects occur in the O/Si/S zone and the O/C zone, which are the ones with the highest molecular masses.

In the O/Si/S zone, there are molecular effects on the temperature early on, between 100 and 300 days. The contributions from the most important cooling species (both molecular and atomic) are shown in the left panel of Fig. 5. The initial molecular cooling is primarily done by SiO, which accounts for up to 20% of the cooling, and to a lesser degree SiS (up to 5%). With molecular cooling, this zone has a temperatures up to 200 K lower compared to the molecule-free model, at early phases. The effect then lessens with time, and after 300 days there are negligible differences between the models with and without molecular cooling.

A similar pattern can be seen in the O/C zone (right-most panel in Fig. 5); at the early stages, there is significant CO cooling (up to 35% of the total cooling), which leads to a ~400 K lower temperature in this zone. However, by 300 days the molecular cooling contribution are less than 10%. From around 450 days, CO starts to become important again and a temperature gap arises, reaching a few 100 K at 600 days. The trends in both of these zones are similar; initially, there are prominent molecular cooling effects, which then lessen with time. As the abundances of CO and SiO, the most important molecular coolants, are essentially constant before for the first few hundred days (Fig. 1), a decrease of these species is not the cause of them becoming less relevant coolants with time. Rather, this behavior can be explained by looking at the evolution of cooling per particle, which is shown in Fig. 6. The figure shows that the decline in cooling per particle for the molecules is steeper than those for the atomic species. Consequently, if the number of particles of each species is approximately constant, molecules will become less important when compared to atomic coolants.

For the other zones, the molecular influence on the temperature is minor for most of the modeled epochs. In the Si/S zone, where the only significant molecule formed is SiS, at 10⁻⁶–10⁻⁵ M_⊙, the temperature effect is always small, ≲50 K. In the O/Ne/Mg zone, as seen in Fig. 2, there is an increasing amount of molecules forming with time, mainly CO and SiO. The masses are still relatively small, however, and only at the last modeled phase do we see any temperature effect. At this time (600 days), CO causes about 10% of the cooling, which leads to a temperature decrease of about 90 K.

The Fe/He zone has negligible molecular formation and is consequently not impacted by molecular cooling.

Our results broadly support the parameterized ansatz in Jerkstrand et al. (2015), where it was assumed that molecules have neglegible thermal effect in the Fe/He, Si/S and O/Ne/Mg zones in a stripped-envelope supernova. However, in some of those models, strong molecular cooling of the O/Si/S and O/C zones were assumed. The results here indicate that this can generally not be assumed for all epochs.

While a temperature difference of <1000 K may not seem very large, the atomic emission lines have exponential dependencies on temperature and can thus be significantly affected by also quite small relative temperature changes. In Sect. 9 we explore the impact on optical atomic emission by the molecular cooling and the sensitivity to the gas density.

5 Destruction by Compton electrons

As discussed in previous sections, collisions with nonthermal electrons is an important destruction mechanism for all the investigated molecules. In this work, we directly calculate the Compton ionization rate of reaction (3), and from this infer the Compton dissociation rate of reaction (4), as described in Sect. 2.2.3. In previous works these rates have often been estimated using the Eq. (5) prescription, for which Liu & Dalgarno (1995) reported W values for different Compton destruction processes of CO. Typically these W value are assumed for molecules other than CO as well (see e.g., Liu & Dalgarno 1996; Cherchneff & Dwek 2009; Sluder et al. 2018; Liljegren et al. 2020).

Figure 7 shows the ratios between our calculated rates for ionization by Compton electrons and the rate estimated by using Eq. (5) with the W value Liu & Dalgarno (1995) estimated for CO for this process (W = 34 eV). For SiO, SiS and SO we compare the two rates using results from the O/Si/S zone, and for CO we use results from the O/C zone, as these are the zones where the different molecules are most abundant.

For CO the two estimates agree to within a factor 2 for all the investigated phases in the standard model, with the W = 34 eV estimate consistently yielding a lower rate by about 35%. We can conclude that using the Liu & Dalgarno (1995) formula gives a quite satisfactory estimate of the Compton electron destruction rate for CO. For SiS and SiO our calculated rates are around a factor 3.5 larger compared to what Eq. (5) gives with W = 34 eV. The destruction rate of SO, similar to CO, is also higher in our calculation, but not by as much.

The differences between our calculations of the Compton destruction rate and that of Liu & Dalgarno (1995) can be explained by a few key differences in methods and data used. First, they calculate their rates for a fixed composition of 49.5% Ο I, 49.5% C I and 1% CO I, while our rate is calculated during a SUMO run, which uses the SUMO composition and ionization solutions. Second, for CO we use a different, more recent cross-section (Itikawa 2015) than that of Liu and Dalgarno (Liu & Victor 1994), which likely accounts for some of the difference. Third, using the Liu and Dalgarno CO estimate for other molecules will result in uncertainties. As seen in Appendix B, the cross-sections do differ between the molecules, which will yield differences in the rates.

Calculating backward, we can infer what the corresponding W values for CO, SiO, SiS and SO, would be for our standard model, these are listed in Table 4). Equation (5) does not, however, exactly reproduce the behavior of our calculated Compton rates even for these values as there is a slight downward slope in the ratio of the two for any used W. Consequently, using Eq. (5) with the updated W values still yields an error of up to 10% when compared to the calculated Compton destruction rates. While our reported W values are derived for the standard model, results from the parameter investigation described in Sect. 8 indicate that using these yields gives satisfactory estimates of the Compton ionization rates for the included molecules. We find at most deviations of around 20% between Compton destruction rates using our estimated W values and rates calculated with SUMO.

Fig. 5 Fraction of the cooling due to different species, for the three oxygen zones. Species are plotted if they contribute more than 5% of the cooling at any of the modeled epochs. We note that the same species can cool at one phase, and heat at another phase (i.e., have a negative cooling fraction).

Fig. 6 Cooling per particle, for the species plotted in Fig. 5. Only epochs when the species cool (rather than heat) are included.

Fig. 7 Ratio between the Compton destruction rate calculated with SUMO, for our standard model, and the one estimated using the Liu & Dalgarno (1995) prescription, which was calculated for CO.

Fig. 8 Model spectra in the MIR, where each line is the spectra at different times in days indicated by the right label, for the standard model. The fundamental (dv = 1) and first overtone (dv = 2) molecular bands are indicated. Unblended molecular emission has a light gray background, blends from two bands mid-gray, and blends from three bands dark gray. The flux here is on a log-scale.

6 Mid-infrared spectra

Figure 8 shows the mid-infrared spectra of the standard model, at different times. For all simulated phases there are strong molecular emissions in several spectral regions. The details of the spectral signatures from the different molecular species are discussed below. For a discussion of the atomic MIR emission in SNe, see Jerkstrand et al. (2012). Possible future observations of the model IR spectra are further discussed in Sect. 11.

6.1 CO

Carbon monoxide has a fundamental band emission that extends between ~4.3 and ~6.8 μm, with a peak at around 4.6 μm. This is the most luminous molecular emission at all the investigated times, and, except for the very late phases, additionally the most luminous emission in the MIR spectral range. The CO fundamental band should consequently be the easiest MIR emission to observe. This emission has been observed for a few Type II supernovae, most notably SN 1987A (see e.g., Meikle et al. 1989; Rank et al. 1988; Wooden et al. 1993), and with the Spitzer Space Telescope for SN 2004di (Kotak et al. 2005), SN 2004et (Kotak et al. 2009), and SN 2005af (Kotak et al. 2006). The Spitzer observations consist of the red tail of the band, as the spectral window starts at 5 μm. Any photometric measurements from this region will likely also be dominated by CO emission, which is discussed in for instance Catchpole et al. (1988), Andrews et al. (2011), Ergon et al. (2014).

The rapid growth in CO mass seen in the model after 550 days (Fig. 1) leads to a plateauing of the CO luminosity decline between 500–600 days, with a possible regrowth at yet later times. Observational monitoring with 50–100 days cadence would be of importance to test such predictions and constrain the ejecta properties, as the transition time depends on the ejecta density.

It should be noted that the CO fundamental band overlaps with the overtones of SiO and SO (and marginally with SiS), as seen in Fig. 9. In the model here, only the SiO first overtone is predicted to be visible of these, with both SO and SiS overtones being very weak compared to the CO emission. For a discussion about the SiO first overtone see Sect. 6.2.

The CO first overtone ranges from ~2.2 μm to ~3 μm with a peak at around 2.3 μm. While weak compared to the fundamental band, it contributes around 50% of the flux in this wavelength region, and is visible especially in the early phases (100–200 days). The CO overtone is the most observed molecular feature for all types of supernovae, as its wavelength is positioned so that it possible to observe with ground-based telescopes. For Type II SNe the first overtone has been observed in the following: SN 1987A (Catchpole et al. 1988), SN 1995ad (Spyromilio & Leibundgut 1996), SN 1998S (Gerardy et al. 2000), SN 1998de (Spyromilio et al. 2001), SN 1999em (Spyromilio et al. 2001), SN 2002hh (Pozzo et al. 2006), SN 2004di (Kotak et al. 2005), and SN 2011dh (Ergon et al. 2015). For Type Ibc there are another handful of observations available: SN 2000ew (Gerardy et al. 2002), SN 2007gr (Hunter et al. 2009), SN 2013eg (Drout et al. 2016), and SN 2020oi (Rho et al. 2021). It should be noted that CO is not universally observed in type Ibc SNe; for instance Pandey et al. (2021) reports the absence of CO first overtone features up to 319 days in SN 2012au. In Sect. 10 we compare the luminosity in the CO first overtone emission in our model with the observed luminosity for Type Ibc SNe.

Fig. 9 Blend between the fundamental bands of CO (blue) and the first overtone of SiO (green), SO (pink) and SiS (green). Only SiO will have significant contribution from the first overtone bands; both SO and SiS first overtone emission is weak at all simulated times. The flux here is on a linear scale.

Fig. 10 Blend between the fundamental bands of SiO (orange) and SO (pink), at different times.

6.2 SiO

The second strongest molecular emission in the model is from the fundamental band of SiO, which ranges from about 7 μm to 10 μm, with a peak at about 7.9 μm. The SiO fundamental band is initially almost as luminous as its CO counterpart, but gets noticeably weaker compared to the CO fundamental band at later stages, when the CO mass increases after 450 days. The fundamental band of SiO has been observed for SN 1987A (Aitken et al. 1988; Roche et al. 1991) and additionally in two Type 11 SNe; SN 2004et (Kotak 2009) and SN 2005af (Kotak et al. 2006). It has not been observed in any Type Ibc SN.

As seen in Fig. 10, which shows three different time instances, there is a blend between the SiO fundamental band and the fundamental band of SO. Initially, when the SO abundance is low, all the emission in this range comes from SiO. However, when the SO molecular mass increases the emission eventually becomes comparable to that of SiO. While blended, there is still about 1 μm separating the peaks of the SiO and SO fundamental bands, so the two contributions can potentially be distinguished.

The SiO overtone lies between ~3.95 μm and ~4.8 μm with a peak at ~4.1 μm. This has partial overlap with the blue wing of the strong CO fundamental band, and consequently only part of the SiO emission, between 3.95 and ~4.2 μm, is readily visible (Fig. 9). With time this emission rapidly dims as the gas cools; identification prospects are therefore best at the early epochs. The SiO first overtone was tentatively identified in SN 1987A by Meikle et al. (1989), however there was no clear sign of it in the Kuiper observations (Wooden et al. 1993). Thus there is yet no unambiguous detection of this emission in any supernova.

6.3 SO

The fundamental band emission of SO extends from ~8.5 μm to ~10.7 μm, with a peak at around ~8.7 μm. This emission becomes visible once the mass of SO starts to become significant at about 300 days. As mentioned in Sect. 6.2, there is an overlap between the SiO fundamental band and the SO fundamental band, shown for three different time instances in Fig. 10. As seen, at later times the emission from SO is visible and distinct from that of SiO due to the separation of around 1 μm. These results indicate that if SO forms in the amount predicted by the model, it should be possible to observe and separate from SiO emission. The best time window for this is 300–600 days.

The first overtone of SO ranges from ~4.4 μm to ~5.1 μm, which overlaps completely with the strong emission of the fundamental band of CO. Therefore, the SO first overtone will likely not be possible to observe, and in our model it is always much weaker than the CO band. Rovibrational SO emission has not yet been identified in any supernova. There are, however, observations of pure rotation transitions in the radio range for SN 1987A at an age of ~10 000 days, which has been attributed to SO (Matsuura et al. 2017).

Fig. 11 Blend of the SiS fundamental band (green) and [Ne II] 12.81 μm (brown), at three different times. At times >300 days the SiS band is not visible anymore, swamped by the neon line.

6.4 SiS

SiS has a fundamental band between ~12.5 μm and ~15.6 μm with a peak at ~13.2 μm, shown at three different times in Fig. 11. Emission here is distinct up until 300 days, at which time the SiS mass decreases significantly. This SiS band is weaker than the fundamental bands of CO and SiO at all phases. The model suggests, nevertheless, that this emission could be possible to observe as there is no blending with other molecular emission, and blending with atomic emission (mainly Ne II emission at 12.8 μm and two Co II lines at 14.7 and 15.4 μm) should be possible to distinguish.

The first overtone band of SiS is positioned between ~6.2 μm and ~7.8 μm, with a peak at ~6.8 μm. Due to blending with strong atomic lines of Ni II at 6.6 μm, Ar II at 6.9 μm and Ni III at 7.3 μm, this emission is not distinct in model spectra at any phases and would be hard to detect.

SiS has so far not been identified in any supernova. Identification has been made difficult by the challenges to observe far out into the MIR, where the fundamental band is positioned, and by this band being relatively weak compared to the SiO and CO bands.

7 Model properties and line shapes

It may be possible to derive SN characteristics from the line shapes of molecular bands. The line shape will depend on the excitation structure of the molecule and this, in turn, depends on different properties, most importantly the temperature and the density. A higher total density typically also means a higher electron density, which leads to higher collisional excitation rates for the molecules, as the collision-induced transition rate depends on n_e (Eq. (9)). A higher temperature also leads to a more broadly excited molecule, and emission over a broader wavelength range.

In the standard model, however, temperature and density have a complex relationship; higher density leads to more molecules being formed, which in turn can lead to a more efficient molecular cooling and a lower temperature. Additionally, both the temperature and the density will impact how close to LTE a specific species is. To disentangle the effects of temperature and density on the shape of the bands we here calculate NLTE molecular spectra for a grid of different temperatures and densities. One parameter is the total density n_tot. We then set ratios of electron density n_e = 0.1 × n_tot and molecule density n_mol = 10⁻³ × n_tot, which corresponds roughly to the standard model results (which is 11% ionized, and CO makes up 6 × 10⁻⁴ of the total density at 100 days). It should be noted that these results are not from full SUMO simulations, as the temperature and molecular formation are not calculated self-consistently, and we assume that the lines do not interact (self-absorption is included by the Sobolev formalism but we ignore line-to-line absorption). These types of calculations will, however, give an indication of how the shape of the bands depends on density and temperature when varied independently.

Figure 12 shows the CO first overtone, at 100 days, for a range of different densities and temperatures. For reference, at this time instance the total number density of the O/C zone, where the majority of CO is formed, is 8.8 × 10⁹ cm⁻³ and the temperature is 5110 K. The black lines of each panel are the resulting bands for the specified temperature and density, and the gray lines are the band shape of the lowest temperature and lowest density plotted as a reference. All lines have been normalized to have peak flux unity.

By comparing the resulting black lines with the gray reference lines we can deduce that increasing temperature and increasing density have similar effects. The line shapes become broader and higher-order band heads (e.g., 3-1, 4-2) become more prominent as higher vibrational bands are more populated. There is then a degeneracy, and changes to the line shape can be driven by either of these factors. While this broadening is larger for the investigated range of densities, the density range is 4 orders of magnitude in comparison to a factor 2 range in temperature. Increasing the density also shifts the maximum of the band to longer wavelengths. This is the tendency when increasing the temperature as well, however the effect is smaller for the investigated parameter space. Thus, band head position may be linked to density constraints.

The pink curve in the figures represents the observed CO first overtone in the Type Ic supernova SN2007gr, from Hunter et al. (2009), at 137 days after explosion (which are the observations with best signal-to-noise available). For comparison to the model results, the observations have similarly been normalized, and are over-plotted on the density and temperature combinations that give model spectra agreeing reasonably well with the observations. Several parameter combinations reproduce the line shape of these observations to similar degree. This demonstrates the degeneracy between temperature and density, as the observations can be fit with either a model with low density at high temperature, or with higher density but lower temperature. Information from these line shapes thus probably need to be complemented with other diagnostics. For other molecules the trend is similar. While the details might differ, due to the inherent differences between the various molecules, increasing the temperature or the density leads to a broader line in both cases. Typically, very high densities will also shifts the line maximum to longer wavelengths, as for the CO dv = 2 line.

It should be noted that these calculations assume that the lines in the band are optically thin; this likely does not hold for any of the fundamental bands of the included molecules at early times (≲200 days). Additionally, some of the lines in the CO dv = 1 band become optically thick again after 500 days, with the large increase in CO mass. This makes the line shape of these bands more complex and requires full SUMO calculations to predict more accurately.

Fig. 12 Line shapes of the CO dy = 2 line, for a grid of temperatures and densities, at 100 days. The black line is for the specified density and temperature in each panel. The gray line is for the lowest density and temperature and is used as a reference. The pink lines are observations from Hunter et al. (2009), plotted in panels where the model results are a good fit to the observations. In these simplified simulations we use nCO = 10−3ntot and ne = 0.1ntot.

Fig. 13 Evolution of the masses of CO, SiO, SiS and SO, for models with different χ values (higher χ means higher O-zone density).

8 Parameter space investigation

From previous studies, it is known that ejecta structure can significantly affect the amount of molecules formed (e.g., Gearhart et al. 1999; Cherchneff & Dwek 2009; Liljegren et al. 2020). We perform here a small parameter study aimed to investigate the sensitivity to density for the Type Ibc case. This serves as a catch-all for varying either expansion velocity or the clumping of the gas. We do this by varying χ in Eq. (10); higher χ means denser O-zones and less dense Fe/He and Si/S zones. We here test χ = [50, 200, 400, 800] (the standard model has χ = 100.) These values correspond to O-zone overdensities of 1.8–14. For each χ model SUMO is rerun completely, meaning also variations in gamma-deposition per zone, etc.

8.1 Molecular masses and cooling

Figure 13 shows how the molecular masses are affected by density variations. Overall, the trends and patterns seen in the standard model (discussed in detail in Sect. 3) are retained when changing the density. The timings of those trends are, however, shifted due to the change in the physical conditions. With a few exceptions, the molecular abundances increase when the zone densities increase (the Si/S zone density decreases when the O-zone densities increase).

For CO, the steep increase of CO mass, which occurs at 600 days in the standard model due to the reduced efficiency of the destruction reaction (18) at lower temperatures, occurs at significantly earlier times in higher density models - down to 200 days in the χ = 800 case. The increase in CO abundance leads to a stronger cooling, which lowers the temperature compared to the standard model and shuts off reaction (18) earlier. This lowering in temperature in the higher density models is shown in the lower right panel of Fig. 14, where the O/C zone temperature evolution for the different models is plotted. As seen, the higher density models will have lower temperatures at all stages. Consistently, the lower density model (χ = 50) has a lower CO abundance, a later increase in CO mass, and a later decrease in temperature, following the same logic. At 600 days the CO mass is ~10⁻³ M_⊙ for the lowest density model, and ~5 × 10⁻² M_⊙ for the highest density model.

SiO shows a very similar trend for all models at different phases, with an increased SiO molecular mass with increased density. The SiO mass in the highest density χ = 800 model is larger at 600 days, and the largest mass of any of the investigated models.

For SiS and SO, which have strongly connected chemistry (Sect. 3), the overall trends are also the same as in the standard model. The higher density models form more molecules at earlier times, which in turn cools the surroundings more efficiently. The behavior of these two species then depends on the temperature-sensitive reaction (23), as previously discussed. As the temperature is lower in the high density models there is consequently a decline in SiS mass and associated rise in SO mass at earlier phases: 300 days for the χ = 600 model, and 200 days for the χ = 800 model, compared to the 330 days for the standard model.

The temperature evolution for the O/Si/S zone is shown in the upper right panel of Fig. 14. The increased molecular abundance in this zone influences the temperature; a higher density leads more than 1000 K lower temperature at some of the earlier phases, between 200 and 400 days. This is primarily due to the increase in SiO, which provides ~80% of the cooling at 200 days for the χ = 800 model. There are smaller cooling contribution from SiS and SO, as well, for the respective models where they are abundant.

Similarly, the increase in molecules also affects the temperature in the O/Ne/Mg zone. As there is very little sulfur in this zone, the most abundantly formed molecules are CO, which is the main molecular coolant, and SiO. The difference in temperature is up to ~1500 K in the highest-density model, at 100 days. After 300 days we see little difference in the temperature between the investigated models. As seen in Fig. 14, there is no significant cooling either due to more molecules or the increased densities in the Si/S zone, for any of the investigated models.

Dust can possibly start forming once the temperature falls below ~2000 K (e.g., Tielens et al. 2005; Sarangi & Cherchneff 2013, 2015; Sluder et al. 2018), which is indicated in the figure. Some commonly found dust species, such as amorphous carbon or corundum dust (Al₂O₃ monomers) can condense at higher temperatures than for example fosterite dust (Mg₂SiO₄ monomers) or iron dust (Gail 2010; Lattimer & Grossman 1978). The indicated start of dust formation then represents a lower limit on the time at which dust can potentially form, in the different zones. Fig. 14 shows that, depending on the O-zone densities, such a transition is predicted to begin in the range 250–380 days.

Fig. 14 Thermal evolution in the Si/S, O/Si/S, O/Ne/Mg O/C zones, for models with different χ values. Assuming dust formation to begin at 2000 K, the models predicts an on-set of this sometime in the time interval 250–380 days. At high χ the first dust formed would be carbonaceous (in the O/C zone), whereas at low χ-values it would be silicate-based (in the Si/S and O/Si/S zones).

Fig. 15 Spectra of models with different densities, at three different times. The fundamental and first overtone emissions of CO, SiO, SiS and SO are marked.

8.2 Changes to the IR spectra

The difference in density, molecular formation, and temperature will also influence the IR spectra, which is shown in Fig. 15 where the emission is plotted for models with different χ, at three time instances. It is clear that molecular contribution is important in this spectral domain; for all the models with molecules there is significant molecular emission for the three time instances, for all densities.

Increasing the density of the models typically leads to a stronger molecular emission, however, the effects are complex. At 100 days the emission from the higher-density models are very similar to that of the lower-density models, due to optical depth effects as many lines are optically thick at this epoch. Consequently, observations of these times could be used to constrain other quantities such as zone volumes.

At 330 and 600 days, the differences are much larger, and the emission of the CO and SiO fundamental bands are more than an order of magnitude stronger at the highest O-zone density (χ = 800) compared to the standard model. Observations at these epochs could therefore be potentially used to deduce the O-zone densities.

Some of the complex behaviour for the molecular masses also imprints on the spectra; at 330 days the χ = 800 model has very weak SiS fundamental band emission because the SiS mass has decreased significantly in this model at this time (as discussed in Sect. 8.1). At 500 days the SiS emission is not visible for any of the models.

For the higher-density models, the CO, SiO and SiS first overtone bands become much stronger, even though at first glance this can seem counter-intuitive as the temperatures of these models are lower. However, several different processes influence how strong the molecular emission for a species is, and they can interact in a complex manner, which makes the impact of for example changing the density in the models complex to deduce. A higher excited population of molecules will typically lead to a stronger first overtone emission. The electron collision excitation rates, the main excitation mechanism for molecules, increase with temperature (which is lower at higher densities) and the electron number density n_e (which is higher at higher density). Additionally, the Sobolev self-absorption becomes more prominent at higher densities, which can lead to higher excitation of the molecules. The net effect for the mentioned species is a higher excitation, and, thus, more emission in the first overtone band.

For SiS, the band at 7 μm is visible at 100 days for both the χ = 400 and χ = 800 models, which is not the case for the standard model. The SiO first overtone band is much stronger at 330 days. The CO first overtone band is also more prominent; at 100 days this emission dominates the region, with a peak flux similar to the strong CO fundamental band.

Fig. 16 Optical spectra at three different times for different densities as labeled. By comparing to Fig. F.1, one can conclude that the [O I] 6364, 6364 reduction for higher χ is driven by molecular cooling at 330 and 600 days, whereas at 100 days is it driven by other density-dependent processes.

9 Optical spectra

While IR observations can enable direct detection of molecular emission, the presence of molecules may also influence the optical spectra. This can happen by a few different mechanisms. 1) Atoms get locked up into molecules. 2) Cooling by the molecules can decrease the zone temperature. 3) The molecular opacity can change the radiation field. Here we focus on the first two effects. We include rovibrational line opacity in the transfer, however this likely have quite minor effects on the radiation field.

In our standard model, the molecules do not form in large enough amounts to act as significant sinks for any atomic species to affect their emission by the first mechanism. The cooling by the molecules can be significant, however, lowering the temperature by several hundred degrees at specific times. The amount of molecular cooling varies with time and is very zone-dependent, its impact is thus quite complex. In the O/C zone, there is an effect at early and late times, but not at intermediate times. Fig. 16 shows the optical spectra from 4000 Å to 10 000 Å at three times (100, 330 and 600 days) for the standard model (dark blue line), the standard model without any molecules included (green dashed line), and three different density models (χ= 50, black line, χ = 400, pink line, χ= 800, yellow line).

Comparing the standard model with molecules (blue solid line) and without (green dashed line), there are some differences.

At 100 days, some cooling due to molecules occurs in the O/Si/S and C/O zones, which slightly impacts the [O I] 6300, 6364 Å line and the Ca II and [Fe I] blends at 8400–8800 Å. At 330 days there is no longer significant molecular cooling in any of the zones, and the spectra of the standard model are more or less identical to the spectra of the model with no molecules. At 600 days, the O/C zone is cooled by the increased amount of CO. This has a large impact on [O I] 6300, 6364 Å (almost 50% weaker), which is commonly used as diagnostic for the core mass. For the standard model, we can conclude that the optical spectra are only slightly affected by molecule formation and cooling, except for the mentioned [O I] 6300, 6364 Å feature which can become significantly quenched at late times.

An import aspect of how results may vary for other CO core masses is that for this relatively massive one (6.8 M_⊙), the molecule-forming zones (Si/O and O/C) constitute only 43% of the total O-zone mass (with 57% being in the O/Ne/Mg zone). At lower CO core masses, this number is typically higher, for example 68% in a 4 M_⊙ CO core (Jerkstrand et al. 2015). As such, the spectral impact of molecule formation is almost certainly larger for smaller CO cores. However, in the context of trying to determine whether some Ibc SNe come from relatively massive Wolf-Rayet progenitors, the results here suggest that optical models ignoring molecular chemistry are not likely to severely underestimate the CO core mass – optical emission is not predicted to be strongly quenched by molecular IR cooling – at least not in the early nebular phase.

Figure 16 also shows the optical spectra of models with different densities; χ = 50 (black line), χ = 400 (pink lines) and χ = 800 (yellow line). The variation in the spectra are due to two effects; the denser models tend to produce more molecules, which in turn leads to a larger cooling, and higher density in itself can change the emission even if the molecular impact stays the same. To disentangle these effects we plot in Fig. F.1 spectra for models without molecules for different χ.

From Fig. 16, the single model with lower O-zone densities (χ = 50) has a spectrum is very similar to the standard model at all times. The higher-density models, however, differ significantly. At 100 days both the χ = 400 and the χ = 800 models have less distinct O-zone emission lines. Inspection of Fig. F.1 shows that this is so also with molecules switched off – it is therefore an effect not driven by the chemistry.

At 330 days, however, the [O I] 6300, 6364 feature is not sensitive to density in the molecule-free models, but is drastically weakened (factor 2) in the chemical models. Mg I] 4571 is, on the other hand, still density-sensitive also without molecular effects.

At 600 days, higher density gives stronger [O I] 6300, 6364 in molecule-free models, but weaker [O I] 6300, 6364 in models with molecule cooling. At any χ, molecular cooling here strongly reduces the strength of the line. Other lines are however little affected.

Fig. 17 Total luminosity in the first overtone of CO, for models with different χ values (lines), compared to observations of stripped-envelope supernovae (data taken from Gerardy et al. 2002, Hunter et al. 2009, Ergon et al. 2015 and Rho et al. 2021). The observed values are scaled with the difference in 56Ni mass between the model and the SN. For SN 2000ew there is no observational estimate of this so the scaling factor is set to 1, with an error bar of factor 5 attached. The error bars indicate the uncertainties in the 56Ni mass when such are reported. There is no reported error for SN 2020oi, so it was arbitrarily set to be a factor 2.

10 Comparison with molecular emission observations

All spectral observations of molecules in Type Ibc SNe have been of the first CO overtone emission (see Table 5), which is observable from the ground. While it should be emphasized that we are not attempting to reproduce any specific SN with our single model here, we nevertheless compare the model luminosity in the CO first overtone to the available observations to gauge if the model is in the right ballpark for this feature. The four observed SNe (three Type Ic and one Type IIb) to which we are comparing the model results are listed in Table 5.

Figure 17 shows the observations as scatter points and the results for models with different χ values as the lines. The measured luminosities are additionally scaled with the ⁵⁶Ni mass of the model divided by the known ⁵⁶Ni mass of the SN. This is to take differences in total powering level into account. For SN 2000ew no ⁵⁶Ni mass has been estimated, we then put this factor to 1.

As can be seen, the models predict an increase in the degree of O-zone clumping (χ) to give a higher luminosity in the CO first overtone. The growth flattens out as one exceeds χ ~ 600. The strong increase in emission seen in the χ = 800 model between 200 and 250 days is due to the rapid CO mass growth that occurs in these models at that time.

Comparing the modeling results with the observations we can conclude that the modeling results fall into the same range as the observed CO first overtone luminosities. At face value, there is a good match between the observations and the standard model, but we do not ascribe any particular significance to this until more detailed modeling where also other properties of each SN are taken into account. Rather, the takeaway point from the figure is that the modeling does not produce results in any glaring conflict with current observations of the CO overtone in stripped-envelope SNe, and this is encouraging for the methodology and continued work.

11 Future possibilities for SNe IR spectra observations

To get more information about different types of molecules observations throughout the infrared region are needed. While these SNe, in their nebular evolution, will be the brightest early (around 100-200 days) observations from two or more epochs would be desirable to investigate the time evolution of the formed molecules. A suggestion based on the model results from this work would be to observe at ~100 days and -350 days.

At ~100 days, the CO, SiO, and SiS fundamental bands, as well as the CO and SiO first overtones, are predicted to be observable. As the molecular fundamental bands are insensitive to densities at these times (discussed in Sect. 6), observations could be used to investigate other properties, such as zone masses or volumes. It is a promising epoch to search for the first SiS detection in an SN.

At -350 days, SiS becomes weak in the model, but SO has become abundant and its fundamental band contributes significantly to the blend with the SiO fundamental band at 7–10 μm. At these epochs, there are large differences in the molecular emission between models with different densities, which therefore can be constrained. It is a promising epoch for the first SO detection in an SN.

The IR regions of the spectrum discussed in Sect. 6 will fall within the wavelength coverage of JWST. The fundamental and first overtone of CO and the first overtone of SiO are in the range of JWST/NIRSpec , while the bands further to the red will be observable with the Medium Resolution Spectroscopy (MRS) instrument of JWST/MIRI.

To investigate the feasibility of observing the different molecular lines, at different phases with JWST, the JWST Exposure Time Calculator⁵ (Pontoppidan et al. 2016) was used, with the model spectra presented in Fig. 8 as template. With the ETC tool, the exposure time needed for the relevant filters to achieve a given S/N for the line peaks can be estimated. We placed the SN at a distance of 10 Mpc for these calculations.

At 100 days the first overtone of CO and the blend between the fundamental band of CO and the first overtone of SiO are covered by JWST/NIRSpec with the filters G235M and G395M, respectively, and bright enough that short exposure times (~1 min per filter) will yield S/N > 50. The SiO and SiS fundamental bands are covered by JWST/MIRI MRS filters and require exposure times of ~2 min for the SiO line and ~5 min for the SiS line to get an S/N > 10. Consequently, at 100 days the discussed molecular lines could with ease be observed with JWST for an SN similar to the standard model in this work.

At 360 days the SN is significantly dimmer and while longer exposure times are needed observations at these times are still attainable. The strong CO fundamental band still needs short exposure times only (~1 min yields an S/N of ~20). However, the SiO and SO fundamental band blend require a longer exposure time of ~1 h to get an S/N > 10. The SiO emission may also be quenched by silicate dust depletation at these times (Sarangi & Cherchneff 2013). The SiS fundamental band flux predicted by our models is at this time not bright enough to be observable.

Observations beyond 400 days become challenging; the fundamental band of CO is still observable with JWST/NIRSpec at 600 days at achievable exposure times, however, bands of the other molecules would require a very long exposure. Observations at these epochs would realistically require an SN significantly closer than 10 Mpc.

It is therefore realistically achievable to observe two epochs of an SN at ~10Mpc at ~100 and ~350 days. Such observations would be of great value to extend our understanding of Type Ibc SNe, both in terms of the evolution of molecules and early stages of the dust production, and for pinning down the composition and the ejecta mass of the SN. This would give insight into the possible progenitors of stripped SNe.

12 Summary and conclusion

We have presented the inclusion of molecular chemistry and cooling into the NLTE spectral synthesis code SUMO for supernova modeling. The treatment involves the calculation of Compton destruction rates for molecules, a large chemical kinetic network, and NLTE ionization, excitation, and cooling by CO, SiO, SiS and SO. We carry out the calculations with fully time-dependent atomic and molecular kinetics, and temperature.

We have applied this new framework to investigate how molecules form and affect physical conditions, and optical and IR emission in a Type Ic supernova (ejecta mass 5 M_⊙) over the 100–600 days time interval. Our main results are as follows:

• Most molecules have destruction governed by Compton electrons. From our SUMO-calculated destruction rates we provide updated values for the work-per-ion-par quantity W suitable for SESNe: 23 eV for CO, 10 eV for SiO, 9 eV for SiS and 24 eV for SO.

• The five most abundant molecules predicted to form in stripped-envelope SNe are CO, SiO, SO, SiS, and S2. The CO, SiO, and SiS abundances increase over time and are most abundant at our last modeled epoch of 600 days. SO and S₂, on the other hand, are most abundant at our first modeled epoch of 100 days. The maximum molecule masses produced are 2 × 10⁻³ M_⊙ of CO, 3 × 10⁻³ M_⊙ of SiO, 7 × 10⁻⁴ M_⊙ of SO, 3 × 10⁻⁴ M⊙ of S₂, and 1 × 10⁻⁴ M_⊙ of SiS.

• CO is mainly formed in the O/C zone, while SiO and SO primarily form in the O/Si/S zone. The sulfur species, SiS and S₂, initially mostly forms in the O/Si/S zone, however, after 300 days they are most abundant in the Si/S zone. The molecules, therefore, provide tracers of three different nuclear burning layers: helium burning (CO), neon burning (SiO, SO, early SiS and S2), and explosive oxygen burning (late SiS and S₂).

• Molecular cooling occurs primarily in the O/C zone and in the O/Si/S zone, by CO and SiO, respectively. In our standard model, a modest temperature effect of a few hundred K occurs in these zones.

• The impact on optical spectra, such as on the [O I] 6300, 6364 feature commonly used to attempt O-mass determinations, depends on the O-zone density. After some transition time, this line is significantly quenched by molecule cooling. This occurs around 550 days in our default density model but already at 200 days in the highest density model.

• Temperatures low enough to allow for the onset of dust formation (-2000 K) are reached at 350–450 days for zones where silicate dust can form, and 250–600 days for a zone where carbon dust can form. Formation starts earlier in the range for higher density models.

• The molecules are predicted to give strong ro-vibrational emission in the infrared, primarily from the fundamental bands of CO at 4.6 μm, SiO at 7.9 μm, and SiS at 13.2 μm. The SO fundamental band, with a peak at 8.7 μm is blended with the stronger SiO fundamental band. However, after 350 days, when there is a large increase in SO mass, the models indicate it should be possible to separate the two bands and detect SO.

• The first overtone emissions are predicted to be significantly weaker than the fundamental bands, and clearly distinguishable only for SiO (4.1 μm) and CO (2.3 μm). The SO overtone at 4.6 μm coincides with the much stronger CO fundamental band, and the SiS overtone at 6.8 μm coincides with much stronger atomic emission by [Ni II] 6.8 μm and [Ar II] 6.9 μm – both of these are therefore predicted to be difficult to detect. We show that the luminosity of the CO first overtone is sensitive to the density of the O/C zone.

• We investigate the observability of the predicted molecular IR emission with JWST. For typical distances (10 Mpc), most molecular emissions can be observed with good S/N with minutes or hours of exposure times. We identify an optimal strategy to observe the SN at 2 or 3 epochs in the 100–400 days window. We demonstrate how such observations can give information on molecular masses and zone properties, and give the first detections of SiS and SO in supernovae.

The models presented in this study offer valuable insights into the molecular formation process in SESNe, and we predict the observability of different molecules at various stages of a SESN. By combining these results with future observations obtained from state-of-the-art optical, near-IR, and MIR instruments (e.g., JWST), we have the potential to resolve degeneracies and gain a comprehensive understanding of the evolutionary trajectory of type Ibc SNe. Furthermore, these combined efforts may offer clues regarding the origin of SESNe.

Acknowledgements

S.L. and A.J. acknowledge funding from the Swedish National Space Board (SNSB/Rymdstyrelsen), grant 2017-R 95/17. A.J. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (ERC Starting Grant 803189 - SUPERSPEC). The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing, KTH Royal Institute of Technology, partially funded by the Swedish Research Council through grant agreement no. 2018-05973. P.S.B. acknowledges support from the Swedish Research Council through an individual project grant with contract no. 2020-03404. R.B. and S.N.Y. acknowledge funding by the UK Space Agency and Science and Technology Funding Council (STFC) grant ST/R000476/1 and ERC Advanced Grant Project No. 883830. We thank Vincenzo Laporta for clarifications regarding calculations for electron collisions with CO.

Appendix A Molecular structure

Fig. A.1 The possible de-excitation channels for the molecules.

Fig. A.2 The possible excitation channels for the molecules.

Currently, neither theoretical nor experimental data is available for all the collision rates needed to model molecule statistical ro-vibrational populations in SNe environments. We here detail the assumptions and extrapolations used in this work to fill in the gaps in the data, and describe the setup of the model molecule.

Appendix A.1 General setup

For a specific level in a molecule, the possible de-excitation channels are spontaneous emission, stimulated emission, and collisional de-excitation and possible excitation channels are collisional excitation and photoexcitation. The molecules included in the radiative transfer are all heteronuclear diatomic molecules, and follow the selection rules for electric dipole allowed rovibrational radiative transitions, which are Δv = ±1, ±2,,…, and, Δj = ±1. We assume that the dominant channel for v-changing collisions is rotationally elastic, such that the rotation number j is not changed in such an interaction (Laporta et al. 2012). Finally, we assume that the rotational levels within a vibrational level are in relative thermal equilibrium, so that they are distributed according to a Boltzmann distribution with respect to each other. Consequently, if we know the total population of a specific vibrational level and the temperature, we can calculate the population of the rovibrational levels. A schematic of the possible excitation and de-excitation channels can be seen in Figs. A.1 and A.2.

Appendix A.2 Collisional vibrational transitions

For collisions with electrons, cross-sections σ(E) for vibrational transitions to the ground level, such as σ(ν = 1 → ν = 0), σ(v = 2 → v = 0) etc, are available for CO (Ristić et al. 2007). From these cross-sections, the rate coefficients as functions of temperature, q(v → ν = 0)(T), can then be calculated assuming a Maxwellian distribution for the thermal electrons. Rate coefficients for excitation can consequently be obtained through detailed balance relations.

To account for other such transitions in CO where no data are available, scaling laws have to be applied. There are several different implementations of such scaling laws available for collisional transitions, and while physical motivations are typically given for a specific formulation, there is no general consensus for which scheme is most accurate. This issue is discussed in detail in Thi et al. (2013), where they compared three different scaling laws: from Scoville et al. (1980), from Elitzur (1983), and from Chandra & Sharma (2001).

Scoville et al. (1980), basing their reasoning upon the Bethe approximation, suggests that the collisional transition probabilities are proportional to the spontaneous radiative transition probabilities A, such that q(v → v′)/q(1 → 0) ∝ A(v → v′)/A(1→ 0). Therefore, $$q\left(v\to v\prime \right)=\frac{A\left(v\to v\prime \right)}{A\left(1\to 0\right)}q\left(1\to 0\right).$$(A.1)

Elitzur (1983) base their expression on the initial work by Procaccia & Levine (1975), in which vibrational collisions were investigated and analyzed. They suggest that: $$\begin{array}{l}q\left(v\to v\prime \right)=\left(v-v\prime \right)\times q\left(1\to 0\right)\hfill \\ \times \exp \left(-\left(v-v\prime -1\right)\frac{1.5E_v/kT}{1+1.5E_v/kT}\right),\hfill \end{array}$$(A.2)

where E_v is the energy level of v with respect to the ground energy level, Τ is the temperature, and k is the Boltzmann constant.

Based on the Landau-Teller relationship for transitions between adjacent states, Chandra & Sharma (2001) proposed the following relationship: $$q\left(v\to v\prime \right)=\frac{v\left(2v\prime +1\right)}{2v-1}\times q\left(1\to 0\right).$$(A.3)

Another scaling law of interest is from Liu & Dalgarno (1995), where they, similar to this work, investigate molecules in SNe environments. Their physical motivation is in part similar to that of Scoville et al. (1980), however, instead of scaling with the Α-value, they use the absorption oscillator strength f_vv′ for a radiative transition from level v to v′, namely $$q\left(v\to v\prime \right)=\frac{f\left(v\to v\prime \right)}{f\left(1\to 0\right)}q\left(1\to 0\right).$$(A.4)

While the individual estimates for a specific collision rate can vary widely depending on which rule is used, we found through experimentation that the final level populations and, in particular, the final spectra are not very sensitive to this choice. This indicates that the excitation channels described with the scaling laws are secondary to those where data are available. These tests were done for both high-density and low-density models.

Throughout this work, we then chose to use the Chandra & Sharma (2001) scaling law as this is the most modern such scheme, and it is the one deemed most appropriate and used in Thi et al. (2013).

Fig. A.3 Schematic of the vibrational collisional de-excitation through collisions with electrons. The left part shows the collisional transition rate coefficients available from experiments and theoretical work. The right part shows how we divide these into rovibrational rates.

Appendix A.3 Collisional rovibrational transitions

From the data and scaling laws we obtain the collisional transition rate coefficient q(v → v′), for vibrational transitions v → v′. For the purpose of this work, however, we are interested in the rovibrational transitions, and, therefore, need to divide the available vibrational transition rate coefficients q(v → v′) further into rovibrational transition rate coefficients q(v, j → v′, j′), as illustrated in Fig. A.3. q is likely independent, or only very weakly dependent on the j-number in the case of CO (Laporta et al. 2012).

We therefore assume that the collisional rate coefficient q is independent of the j-number, such that q(v, j = 0 → v′, j = 0) = q(v, j → v′, j) = q(v → v), for any j. This further ensures that the collisional transition rates between vibrational levels are conserved with respect to the vibrational transition rate coefficient from previous work, as the sum of all individual transition rates then is equal to the total given rate, such that $$\begin{array}{l}{n}_v\times q\left(v\to v\prime \right)={\displaystyle \sum _j{n}_{v,j}\times q\left(v,j\to v\prime ,j\right)}\hfill \\ ={\displaystyle \sum _j{n}_{v,j}\times q\left(v\to v\prime \right),}\hfill \end{array}$$(A.5)

for a vibrational state v with rotational levels j.

Data for molecule-electron collisional excitation rates, are only available for CO, to the authors’ knowledge. For lack of any better alternative, the collisional transitions rate coefficients derived from both data and scaling laws for CO, were also used for SiO, SiS and SO.

Appendix A.4 Rotational levels

The typical energies of pure rotational transitions are at least an order of magnitude smaller than transitions between rovibrational levels, and the rotational levels within a vibration level will be in relative LTE in the investigated environment, as mentioned. This means that the internal distribution of rotational levels, within a vibration level, can be described by the Boltzmann distribution. In practice, the relative LTE is enforced by assuming that 1) all rotational levels within a vibrational level are connected through collisional rate coefficients C_v,j→v,j′, and 2) these collisional rates are large enough to ensure relative LTE. We arbitrarily picked C_v,j→v,j′=1000 s⁻¹, which is about 1000 times larger than the largest collision-induced transitions between vibrational levels. Through testing, we have deduced that this value always enforced relative LTE in the rotational levels within a vibrational level, for the parameter space used in this work.

Appendix B Compton electron collision cross-sections

Fig. B.1 Overview of the cross-sections of included molecular species, with the original data and the fits, as described in the text.

We use the form defined in Arnaud & Rothenflug (1985) for the Compton electron collision cross-sections σ_Compton[cm²] in SUMO as $${\sigma }_{Compton}=\frac{1}{u\times i^2}\left(a\left(1-\frac{1}{u}\right)+b{\left(1+\frac{1}{u}\right)}^2+c\log \left(u\right)+\frac{d\log u}{u}\right),$$(Β.1)

where u is the impact electron energy e normalized by the ionization energy i such that u = e/i, and a, b, c, and d are constants specific for the cross-section for each species. We fit this function to the cross-sections given in Itikawa (2015) for CO, Joshipura et al. (2007) for SiO and SiS, and Joshipura & Gangopadhyay (2008) for SO. The values for a, b, c, and d for the molecular species included in the radiative transfer calculations are given in Table B.1 These fits and the original data are plotted in Fig. B.1.

Appendix C Ionization energies

The ionization energies of CO, SiO, SiS, SO, CO⁺, SiO⁺, SiS⁺ and SO⁺ are needed for the calculation of ionization rates. Literature values have been used for CO, CO⁺, SiO, SiS, and SO. For SiO⁺, SiS⁺ and SO⁺ the values have been calculated here using the MOLPRO package Werner et al. (2012, 2020). Vertical ionization energies were obtained by performing full valence MRCI calculations using the aug-cc-pV5Z basis sets. Cluster corrected, Davidson relaxed, energies are reported.

We first calculated the equilibrium bond lengths for the SiO⁺, SiS⁺ and SO⁺ doublet spin states. These bond lengths were then used to calculate the energies for the corresponding doubly ionized species. The lower energies, at the given bond lengths, of the singlet and triplet states for Si0⁺⁺, SiS⁺⁺ and SO⁺⁺ were used to obtain the ionization energies, i.e. the singlet state for SO⁺⁺ and the triplet states for SiO⁺⁺ and SiS⁺⁺. In this way the vertical ionization energies were found to be 20.6, 17.4 and 19.4 eV for SiO⁺, SiS⁺ and SO⁺ respectively. The ionization energies used in SUMO are listed in Table C.1.

Appendix D Reactions

The following tables contain the parameters used to calculate the reaction rates used in this work, divided into thermal collision reactions in Tables D.1, D.2 and D.3, and recombination reactions in Table D.4. When possible we refer to the original work, in other cases the database from which the value is from is indicated.

Using the Eq. (2), the reaction rate k(T) for a specific reaction can be retained. The units of both α and k(T) depend on the number of reactants: s⁻¹ for unimolecular reactions and cm³s⁻¹ for bimolecular reactions. The units of γ is K, and β is a dimensionless number.

Appendix E Reaction rate plots

The reaction rates the most abundant molecules in the standard model shown in Figs. E.1 and E.2. The rates of the sulfur and silicon species are from the O/Si/S zone, and the carbonaceous species are from the O/C zone, which are the zones where they are predominantly formed. It should be noted that these rates depend on several different properties; the rate coefficients, which can be temperature dependent, and the amount of available reacting species, which can change with time. These plots consequently do not necessarily reflect the changes of the physical conditions, but indicate which reactions are important for the formation and destruction of a specific species. For a discussion about the interpretation of these plots, see Liljegren et al. (2020).

Fig. E.1 Log of the reaction rates, divided by the density in the O/Si/S zone, with time, for the most abundant molecules in this zone. This indicates which rates are important for the formation and destruction of different molecules, with time. In these plots XI means the neutral atom or molecule, XII represent the once-ionized ion of an atom or molecule.

Fig. E.2 The reaction rates divided with the density in the O/C zone, for formation and destruction processes of CO. In these plots XI means the neutral atom or molecule, XII represent the once-ionized ion of an atom or molecule.

Appendix F Spectra with different χ, without molecules

The differences in the spectra between models with varying densities can be either due to more or less molecules forming, which influences the cooling or due to higher ddensities themselves, which influences the emission. This is shown in Figs. F.1 and F.2, where models with different densities, but with no molecules included, are plotted.

Fig. F.1 The optical spectra for models with different densities and with no molecules included.

Fig. F.2 The IR spectra for models with different densities when no molecules are included.

References

All Tables

All Figures

	Fig. 1 Total mass formed for different molecular species as function of time. Solid lines are neutrals and dashed lines ionized molecules.
In the text

	Fig. 2 Masses of the four molecular species included in the cooling and radiative transfer calculations as function of time. The different lines show the different masses in each zone.
In the text

	Fig. 3 Graphical representation of the part of the chemical network that connects the sulfur-bearing molecular species.
In the text

	Fig. 4 Temperature evolution for the standard model, with molecules included (dashed lines), and excluded (solid lines), for the different zones.
In the text

	Fig. 5 Fraction of the cooling due to different species, for the three oxygen zones. Species are plotted if they contribute more than 5% of the cooling at any of the modeled epochs. We note that the same species can cool at one phase, and heat at another phase (i.e., have a negative cooling fraction).
In the text

	Fig. 6 Cooling per particle, for the species plotted in Fig. 5. Only epochs when the species cool (rather than heat) are included.
In the text

	Fig. 7 Ratio between the Compton destruction rate calculated with SUMO, for our standard model, and the one estimated using the Liu & Dalgarno (1995) prescription, which was calculated for CO.
In the text

	Fig. 8 Model spectra in the MIR, where each line is the spectra at different times in days indicated by the right label, for the standard model. The fundamental (dv = 1) and first overtone (dv = 2) molecular bands are indicated. Unblended molecular emission has a light gray background, blends from two bands mid-gray, and blends from three bands dark gray. The flux here is on a log-scale.
In the text

	Fig. 9 Blend between the fundamental bands of CO (blue) and the first overtone of SiO (green), SO (pink) and SiS (green). Only SiO will have significant contribution from the first overtone bands; both SO and SiS first overtone emission is weak at all simulated times. The flux here is on a linear scale.
In the text

	Fig. 10 Blend between the fundamental bands of SiO (orange) and SO (pink), at different times.
In the text

	Fig. 11 Blend of the SiS fundamental band (green) and [Ne II] 12.81 μm (brown), at three different times. At times >300 days the SiS band is not visible anymore, swamped by the neon line.
In the text

Fig. 12 Line shapes of the CO dy = 2 line, for a grid of temperatures and densities, at 100 days. The black line is for the specified density and temperature in each panel. The gray line is for the lowest density and temperature and is used as a reference. The pink lines are observations from Hunter et al. (2009), plotted in panels where the model results are a good fit to the observations. In these simplified simulations we use nCO = 10−3ntot and ne = 0.1ntot.

In the text

	Fig. 13 Evolution of the masses of CO, SiO, SiS and SO, for models with different χ values (higher χ means higher O-zone density).
In the text

	Fig. 14 Thermal evolution in the Si/S, O/Si/S, O/Ne/Mg O/C zones, for models with different χ values. Assuming dust formation to begin at 2000 K, the models predicts an on-set of this sometime in the time interval 250–380 days. At high χ the first dust formed would be carbonaceous (in the O/C zone), whereas at low χ-values it would be silicate-based (in the Si/S and O/Si/S zones).
In the text

	Fig. 15 Spectra of models with different densities, at three different times. The fundamental and first overtone emissions of CO, SiO, SiS and SO are marked.
In the text

	Fig. 16 Optical spectra at three different times for different densities as labeled. By comparing to Fig. F.1, one can conclude that the [O I] 6364, 6364 reduction for higher χ is driven by molecular cooling at 330 and 600 days, whereas at 100 days is it driven by other density-dependent processes.
In the text

Fig. 17 Total luminosity in the first overtone of CO, for models with different χ values (lines), compared to observations of stripped-envelope supernovae (data taken from Gerardy et al. 2002, Hunter et al. 2009, Ergon et al. 2015 and Rho et al. 2021). The observed values are scaled with the difference in 56Ni mass between the model and the SN. For SN 2000ew there is no observational estimate of this so the scaling factor is set to 1, with an error bar of factor 5 attached. The error bars indicate the uncertainties in the 56Ni mass when such are reported. There is no reported error for SN 2020oi, so it was arbitrarily set to be a factor 2.

In the text

	Fig. A.1 The possible de-excitation channels for the molecules.
In the text

	Fig. A.2 The possible excitation channels for the molecules.
In the text

	Fig. A.3 Schematic of the vibrational collisional de-excitation through collisions with electrons. The left part shows the collisional transition rate coefficients available from experiments and theoretical work. The right part shows how we divide these into rovibrational rates.
In the text

	Fig. B.1 Overview of the cross-sections of included molecular species, with the original data and the fits, as described in the text.
In the text

	Fig. E.1 Log of the reaction rates, divided by the density in the O/Si/S zone, with time, for the most abundant molecules in this zone. This indicates which rates are important for the formation and destruction of different molecules, with time. In these plots XI means the neutral atom or molecule, XII represent the once-ionized ion of an atom or molecule.
In the text

	Fig. E.2 The reaction rates divided with the density in the O/C zone, for formation and destruction processes of CO. In these plots XI means the neutral atom or molecule, XII represent the once-ionized ion of an atom or molecule.
In the text

	Fig. F.1 The optical spectra for models with different densities and with no molecules included.
In the text

	Fig. F.2 The IR spectra for models with different densities when no molecules are included.
In the text