Carbon monoxide cooling in radiative transfer modelling of supernovae

ABSTRACT

Carbon monoxide (CO) emission has been observed in a number of core-collapse supernovae (SNe) and is known to be an important coolant at late times. We have implemented a chemical reaction network in the radiative-transfer code cmfgen to investigate the formation of CO and its impact on SN ejecta. We calculate two 1D SN models with and without CO: a BSG explosion model at one nebular epoch and a full time-sequence (50–300 d) for a red supergiant explosion. In both models, CO forms at nebular times in the dense, inner regions at velocities <2000 km s|$^{-1}$| where line emission from CO can dominate the cooling and reduce the local temperature by as much as a factor of 2, weakening emission lines and causing the optical light curve to fade faster. That energy is instead emitted in CO bands, primarily the fundamental band at |$\sim 4.5 \mathrm{\mu m}$|⁠, which accounts for up to 20 per cent of the total luminosity at late times. However, the non-monotonic nature of the CO cooling function can cause numerical difficulties and introduce multiple temperature solutions. This issue is compounded by the sensitivity of the CO abundance to a few reaction rates, many of which have large uncertainties or disparate values across literature sources. Our results also suggest that, in many SNe, CO-level populations are far from their local thermodynamic equilibrium values. Unfortunately, accurate collisional data, necessary to compute non-local thermodynamic equilibrium populations, are limited to a few transitions.

1 INTRODUCTION

Late-time spectra of supernovae (SNe) facilitate the analysis of progenitor and explosion properties by probing the innermost, metal-rich ejecta. In the nebular phase (⁠|$\sim$|months after the explosion), SN spectra are dominated by emission lines from metals like oxygen, calcium, and iron in the interior of the ejecta, heated by radioactive decay (Fransson & Chevalier 1989; Swartz, Harkness & Wheeler 1989; Jerkstrand 2017).¹ These lines encode information about the structure and composition of the progenitor, albeit in a complex way.

Line formation in SNe ejecta involves a diverse set of processes, including energy deposition from |$\gamma$|-rays, fast-moving non-thermal electrons, and non-local thermodynamic equilibrium (NLTE). Modelling line formation in SNe thus requires a complex radiative-transfer approach that can take into account many processes simultaneously. Existing codes for modelling SN spectra include Monte Carlo codes like sumo (Jerkstrand 2011) and jekyll (Ergon et al. 2018), as well as full radiative-transfer solvers like cmfgen (Hillier & Dessart 2012).

SN ejecta have been posited as a potential site for the formation of molecules and dust (Gehrz 1989; Draine 2011). Observations of SN1987A have detected the presence of emission lines from multiple molecules, including carbon monoxide (CO)² and silicon monoxide (SiO), as well as infrared emission from substantial amounts of dust (Spyromilio et al. 1988; Bouchet & Danziger 1993; Liu & Dalgarno 1994). Any dust and molecules present in a stellar atmosphere or nearby environment at the time of explosion would likely be destroyed by the high (⁠|$\gt 10000$| K) temperatures and kinetic energies experienced during shock breakout; the presence of molecules at nebular times indicates they must have formed in situ as the SN ejecta expanded and cooled, and that they are being actively heated by radioactive decay of |$^{56}$|Ni and |$^{56}$|Co (Cherchneff & Sarangi 2011). Even if molecules are not directly observed, the possibility of core-collapse SNe as significant dust producers (Watson et al. 2015; Matsuura 2017; Sarangi, Matsuura & Micelotta 2018) requires that molecules like CO and SiO were produced as precursors to dust formation.

Accurate spectral modelling of SNe must include contributions from molecular species in order to recreate observed molecular lines and, more importantly, to take their cooling into account. Molecules possess rotational and vibrational degrees of freedom which atoms do not, giving them more possible energy states. These rovibrational energy states have small (⁠|$\lesssim 0.5$| eV) energy separations, and transitions between them emit photons in the infrared range. Transitions between these closely spaced energy states produce spectral lines that overlap (at velocities |$\gtrsim$| tens of km s|$^{-1}$|⁠) and are observed as a broad-band (Herzberg & Mrozowski 1951). Molecules are very efficient at converting thermal energy into photons, i.e. cooling the gas, as the low-lying energy states are populated by collisional excitation followed by radiative decay. If molecules form in substantial quantities in SN ejecta, they can alter the ejecta’s temperature evolution by accelerating its cooling and affecting the appearance of spectral lines which form in the same region as the molecules (for example, if CO forms in the C/O-rich region, O emission lines will be weakened).

Nebular spectra of SNe are dominated by forbidden emission lines from the inner, metal-rich ejecta, which originate in the core of the progenitor³ (Woosley 1988; Fransson & Chevalier 1989; Jerkstrand 2017; Dessart et al. 2021). The emission features in this phase provide direct clues to the progenitor and explosion properties. For example, the oxygen mass in the ejecta (which can be estimated from the flux of the [O i] |$\lambda \lambda$| 6300, 6364 emission doublet) is correlated with the helium or carbon/oxygen core mass of the progenitor (Woosley & Weaver 1995; Maguire et al. 2012; Dessart et al. 2021). If CO forms in the oxygen-rich region, however, the [O i] |$\lambda \lambda$| 6300, 6364 emission feature will be weakened, lowering estimates of the progenitor core mass derived from it (Liljegren et al. 2022). Properly accounting for the cooling effect of CO can mitigate potential systematic errors.

The nature of CO formation in SN ejecta has been investigated using radiative transfer and chemical evolution models covering a range of complexity (though few works have calculated radiation and chemistry simultaneously). Spyromilio et al. (1988) estimated a CO mass of |$\sim 10^{-4}\,{\rm M}_{\odot }$| in SN1987A using a uniform spherical model, assuming the CO-level populations followed LTE and that the ejecta was optically thin. Liu, Dalgarno & Lepp (1992) extended this treatment to NLTE populations and optically thick CO lines, deducing higher CO masses (⁠|$\sim 10^{-3}\,{\rm M}_{\odot}$|⁠).

Chemical models have been used to understand the formation of CO in SN ejecta, most with specific application to SN 1987A. Petuchowski et al. (1989), assuming steady-state chemistry and a handful of reactions, concluded that radiative association dominates formation of CO, with destruction dominated by photonic processes. Lepp, Dalgarno & McCray (1990) found charge exchange with helium to be among the fastest reactions, but later observations and chemical models (Liu et al. 1992) concluded that very little helium must be present in the oxygen-carbon zone in order to reproduce observed CO abundances (⁠|$\sim 10^{-3}\,{\rm M}_{\odot}$|⁠). While these models all assumed steady-state chemistry, later work (Cherchneff & Sarangi 2011) suggested that time-dependent chemistry is necessary. Chemical mixing (or the lack thereof) also plays an important role in molecular formation, as the presence of |$^{56}$|Ni in the oxygen-carbon zone can quickly dissociate molecules (directly via gamma rays or indirectly via heating the gas), limiting molecular formation to only regions with a very low abundance of radioactive isotopes (Ono et al. 2023).

Recently, Liljegren, Jerkstrand & Grumer (2020) included CO formation (assuming steady-state chemical equilibrium) and emission in the Monte-Carlo spectral synthesis code sumo, investigating its effects in Type II (SN1987A) and Type Ibc SNe (Liljegren et al. 2022), finding that molecules form in significant quantities (⁠|$\ge 10^{-3}\,{\rm M}_{\odot}$|⁠) in the ejecta of both SN types and accelerate the cooling.

We undertake a similar investigation by including CO in the radiative-transfer code cmfgen. While sumo, jekyll, and many other codes compute the radiative transfer using a Monte Carlo approach, cmfgen solves for the moments of the radiation field numerically. cmfgen also has the advantage of broad applicability, having been used to model Type Ia (Blondin et al. 2013; Dessart et al. 2014), Ib/c (Dessart et al. 2020), and Type II SNe (Dessart & Hillier 2005; Jacobson-Galán et al. 2022) as well as stars, including O (Martins, Schaerer & Hillier 2005; Bestenlehner et al. 2014), B (Thompson et al. 2008; Keszthelyi et al. 2021), and Wolf–Rayet (De Marco et al. 2000; Massey et al. 2014) stars. While early work using Monte-Carlo codes required simplifying assumptions about the ionization structure or excited state populations (Mazzali & Lucy 1993), more recent work (e.g. using jekyll) treats NLTE fully and includes a very similar set of physical processes to cmfgen. It is useful to compare the results of different approaches applied to the same complex problem.

With the launch of the JWST, high-quality infrared observations of core-collapse SNe should become available in the coming years. While many chemical models have been produced to explain the formation of CO in SN ejecta, very few works have explicitly connected chemical evolution and spectral observations by including molecules in radiative transfer calculations. Moreover, CO has been detected in a variety of SN types, including Type IIP (Kotak et al. 2005, 2009; Rho et al. 2018), IIL (Yuan et al. 2016), Ic (Banerjee et al. 2016; Drout et al. 2016; Rho et al. 2021; Stritzinger et al. 2023), and IIn (Gerardy et al. 2000; Fassia et al. 2001), suggestive that formation of molecules may be commonplace in SNe. Thus it is crucial to compare chemical reaction models with observations of a wide variety of SN types.

Using two simple SN models, we investigate the in situ formation of CO in SN ejecta via a self-consistent set of chemical reactions, and consider the impact of ejecta cooling due to rovibrational emission from CO. The SN models used in this work feature extensive microscopic mixing of composition and are intended as a first demonstration of the impacts of molecules in cmfgen. More physically realistic models (featuring only macroscopic mixing) require longer computation times (Dessart & Hillier 2020b), but will be presented in future work. Section 2 describes the function of cmfgen and the alterations made to it which allow for CO formation and emission, Section 3 outlines the details of the SN models used in our investigation, and Section 4 discusses our results. Section 5 details the behaviour of the CO cooling function which warrants particular attention, Section 6 compares our results to other chemical models of SNe and to observations of CO in Type II SNe, and Section 7 contains our conclusions.

2 CMFGEN

cmfgen is a non-LTE radiative transfer code originally designed to model hot stars with winds (Hillier 1990; Hillier & Lanz 2001). It uses a Newton–Raphson technique to solve for the moments of the radiative transfer equation simultaneously with the rate equations and radiative equilibrium equation in the comoving frame. cmfgen has been adapted to simulate SNe, treating time-dependent effects (Hillier & Dessart 2012), deposition of radioactive decay energy, and non-thermal processes (Li, Hillier & Dessart 2012).

cmfgen is a one-dimensional code, which necessitates the use of some simplifying assumptions. We assume that the SN ejecta is spherically symmetric and expanding homologously. cmfgen includes multiple approaches for chemical mixing: standard boxcar smoothing (Dessart & Hillier 2010) or a shuffled-shell approach (Dessart & Hillier 2020b) which allows for macroscopic mixing but not microscopic mixing.⁴ In this work, we use the simpler boxcar approach to reduce complexity.

Each ion included in a model adds a number of rate equations to the total system. Every included energy level is associated with one equation of the form

where |$\rho$| is the mass density, |$n_i$| the number density of the state i, and |$R_{ij}$| the rate of all processes from state i to state j. |$R_{ij}$| includes photoionization, recombination, line absorption/emission, etc. (see Hillier & Dessart 2012). In general, these processes couple energy states within a particular ion and neighbouring ionization states. When molecular reactions are included (see Section 2.1.1), j can refer additionally to states in different species (e.g. the equation for the O i ground state will include terms for rates involving CO energy levels).

To reduce the total number of equations which must be solved, cmfgen makes use of superlevels, in which multiple energy levels are grouped together. Such a group of (usually closely spaced) energy levels constitutes one superlevel, which is associated with a single equation. The relative populations of levels grouped within a superlevel are assumed to be in LTE. For a more complete discussion of superlevels, see Hillier & Miller (1998). Each atomic species (i.e. hydrogen, helium, carbon, etc.) also gives rise to a constraint equation of the form of equation (2) which specifies that the total amount of a particular species remains fixed:

for |$N_j$| the population of ionization state j, referred to as the number conservation equation. This constraint must be modified when molecules are included; see Section 2.1.3.

2.1 Carbon monoxide

We add CO to cmfgen, and treat it similarly to an atom or ion. Unlike atoms, CO can form in situ in the SN ejecta via chemical reactions.

2.1.1 Chemical reactions

We allow for the formation of CO in the ejecta by way of a set of chemical reactions, with data obtained from the UMIST database (McElroy et al. 2013). A full list of reactions is included in Table A1. The change to a species from each chemical reaction is included in the relevant kinetic equation. These equations are then iteratively solved with a Newton–Raphson approach in the full matrix solution. In a full time-sequence of SN models, we calculate the time-dependent chemical equilibrium (using equation 1) to determine the CO population. However, we also compute some time-independent models in which we assume steady-state chemical equilibrium (modifying equation 1 to eliminate the time-dependent term). The rate of each reaction is calculated using an equation of the form:

where |$n_{R1}$| and |$n_{R2}$| are the number densities of the reactants, T is the local electron temperature, and |$\alpha$|⁠, |$\beta$|⁠, and |$\gamma$| are constants specific to each reaction. When solving equation (1), the rate r replaces either |$n_i R_{ij}$| or |$n_j R_{ji}$|⁠, depending on whether the state |$n_i$| is a reactant or product. Reverse reaction rates are unavailable for many of the reactions we consider, and we do not attempt to compute reverse reaction rates when they are not available. Data on the dependence of reaction rates on the electronic/rovibrational state are likewise missing. We make the following assumptions regarding reactant/product energy states and reaction rates:

• Only ground term states in atomic species participate in chemical reactions (as reactants or products). Generally, this is the ground LS term, which can consist of multiple energy states. For example, for neutral oxygen (O i), this is the |$2\rm {s}^22\rm {p}^4$||$^3\rm {P}$| term, which has 3 J states. Other than determining participation in chemical reactions, details of the individual states (energy, angular momentum, etc.) have no impact on the rate of the reaction, which depends only on the population of the state. Higher energy states are ignored when calculating molecular reaction rates.

• All molecular energy states can participate in chemical reactions, and the rate has no dependence on energy states of molecular reactants/products.

• Potential atomic product states are populated proportional to their statistical weights. For example, the reaction |$\mathrm{C}_2 \,\, + \,\, \mathrm{O}^+ \,\, \rightarrow \,\, \mathrm{C}_2^+ \,\, + \,\, \mathrm{O}$| populates all three J states of the OI ground term: it populates the |$J=2$| state (⁠|$g=5$|⁠) at |$\frac{5}{9}$| of the total reaction rate, the |$J=1$| state (⁠|$g=3$|⁠) at |$\frac{3}{9}$| the total rate, and the |$J=0$| state (⁠|$g=1$|⁠) at |$\frac{1}{9}$| the total rate.

We neglect potential cooling due to these chemical reactions, as the rates are very small compared to other heating/cooling processes.⁵

Creating a complete reaction network necessitates the inclusion of additional molecular species, as CO may be formed via multistep processes that include reactions such as |$\mathrm{C}_2 \,\, + \,\, \mathrm{O} \,\, \rightarrow \,\, \mathrm{C} \,\, + \,\, \mathrm{CO}$|⁠. The species |$\mathrm{C}_2$|⁠, |$\mathrm{O}_2$|⁠, |$\mathrm{CO}_2$|⁠, and |$\mathrm{C}_2\mathrm{O}$| are included in cmfgen. They are, however, excluded from radiative transfer calculations. All molecules other than CO are treated using only a single energy state, no transitions, and no photoionization pathways. We currently limit the reaction network to carbon and oxygen – keeping the reaction network as simple as possible and its computational cost low while we develop our method. Future work will include a much larger network with additional molecules.⁶

Molecules can also be destroyed by collisions with high-energy electrons. In SNe, some electrons are accelerated via inverse-Compton scattering with gamma rays produced by radioactive decay of (primarily) |$^{56}$|Ni and |$^{56}$|Co. These fast-moving electrons are often referred to as Compton electrons, and can collide with molecules and dissociate them into constituent atoms/ions – we include three pathways (Table 1).

We calculate the rate of destruction of CO according to

where |$n_{\mathrm{CO}}$| is the total number density of CO, |$\sigma (v)$| is the energy-dependent cross-section for interaction, v is the electron velocity, and |$\frac{d n_e}{dv}$| is the non-thermal electron spectrum (distribution of non-thermal electrons as a function of velocity). cmfgen explicitly calculates the spectrum of non-thermal electrons present in the ejecta using the Spencer–Fano equation (Li et al. 2012). The cross-sections for interaction were taken from Itikawa (2015), fitted with a power law, and interpolated to the same grid as the electron spectrum. These reactions have no dependence on the temperature nor on the number density of thermal electrons.

2.1.2 CO levels and lines

We obtain rovibrational level and line lists for CO (⁠|$^{12}\mathrm{C}^{16}\mathrm{O}$| isotopologue) from Li et al. (2015). We include vibrational levels up to |$v=10$| and rotational levels up to |$J=75$|⁠. These limits were selected based on available collisional data for CO – data from Ogloblina et al. (2020) cover |$v=0$| up to v = |$10$|⁠. The resulting data consist of 836 rovibrational levels and 9075 transitions, including level energies, quantum numbers, and statistical weights, along with transition wavelengths, oscillator strengths, and Einstein A values. We include only the ground electronic state for CO (⁠|$X^1 \Sigma ^+$|⁠), assuming that excited electronic state populations are negligible under the conditions in which CO forms. In LTE, at |$\approx$|5000 K, the first excited electronic state of CO, |$a^3\Pi$|⁠, which lies |$\sim 6$| eV above the ground state, would have a Boltzmann factor of |$\sim 10^{-7}$|⁠, i.e. a negligible population (Essenhigh 2005). These data are included in cmfgen in the same format as atomic energy levels. When calculating the radiative energy loss due to CO emission, we account for all radiative transfer and optical depth effects. See Section 5 for a discussion of how we characterize the cooling effect of CO line emission.

In cmfgen, CO is treated as any other atom, except without a fixed abundance (see Section 2.1.3). It is treated using the full line-blanketing approach. We make no changes to the radiative transfer in cmfgen aside from the inclusion of CO lines.

We include two approaches to the CO-level populations. The total abundance of CO is determined by the chemical reactions as described in Section 2.1.1, and is independent of the level population approach. In approach 1, we assume that all CO-level populations are in LTE, and they are treated in the equations using a single superlevel (see Hillier & Miller 1998, for discussion of superlevels). In approach 2, we allow CO-level populations to deviate from LTE, using electron-collision cross-sections from Ogloblina et al. (2020) to compute the populations. The cross-sections from Ogloblina et al. (2020), which cover vibrational transitions of the form |$v|0 \rightarrow v^{\prime }$|⁠, were integrated with a Maxwell–Boltzmann distribution to obtain thermal collision rates for vibrational transitions, then scaled using the relation given in Chandra & Sharma (2001):

to cover all possible vibrational transitions. These collision rates, together with radiative de-excitation rates (quantified with Einstein coefficients) determine the relative vibrational-state populations. The level populations are calculated using a single superlevel for each vibrational state, assuming that all rotational states within a given vibrational state have the same departure coefficient. The validity of the full LTE assumption is discussed in Sections 4 and 6.

2.1.3 Other changes to cmfgen

In the SN ejecta, the total number density of a given element is fixed, i.e. all oxygen atoms remain oxygen atoms over time (though some may be neutral O i, some ionized O ii, etc.).⁷ This is quantified in the number conservation equation. The quantity of molecular species, however, is not fixed – CO does form over time, and can be destroyed. We modify the oxygen and carbon conservation equations to account for molecules. The new conservation equation specifies that the total quantity of carbon (oxygen) is fixed, with some in the form of atomic carbon (oxygen) ions, and some in the form of molecules:

3 MODEL SETUP

We calculate two different SN models: a red supergiant (RSG) model (Hillier & Dessart 2019) and a blue supergiant (BSG) model (Dessart & Hillier 2019), each with multiple approaches. Both models are based on a progenitor with a zero-age main sequence (ZAMS) mass of 15 |${\rm M}_{\odot }$|⁠. First each model is calculated without allowing for the formation of molecules (NOCO). Then each model is recalculated allowing for the formation of CO with level populations and line emission treated in LTE (LTECO). The BSG model is also recalculated allowing CO vibrational state populations to deviate from LTE (NLTECO). Both models are based on simulations of 15 |${\rm M}_{\odot }$| progenitors, evolved with the stellar evolution code mesa (Paxton et al. 2018) and then exploded using the radiation hydrodynamics code v1d (Livne 1993). Both models feature a mixed composition, which is performed using a boxcar method following the initial explosion simulation (following the method of Dessart & Hillier 2010). Mixing allows CO to form throughout the ejecta, wherever carbon and oxygen are present.

• BSG: This model is based on a 15 |${\rm M}_{\odot }$| BSG progenitor which produces a Type II-peculiar SN (observational examples of this type include SN1987A, SN2000cb, and SN2009E) (Arnett et al. 1989; Kleiser et al. 2011; Pastorello et al. 2012). Full details of the model can be found in Dessart & Hillier (2019), from which we use the ‘a4’ model. This model is evolved at sub-solar metallicity and scaled to match the |$^{56}$|Ni mass (0.08 |${\rm M}_{\odot }$|⁠) observed in SN1987A (see Spyromilio et al. 1988; Arnett et al. 1989). After the initial explosion, all species are macroscopically mixed throughout the ejecta (see Dessart & Hillier 2019), creating a strongly mixed composition. At 1 d after explosion, the v1d model is remapped to cmfgen assuming homologous expansion (which is a good approximation at that time). This model is evolved to 330 d without the inclusion of molecules. At that epoch, we remove all time-dependent terms (which are small) and add molecules, assuming the ejecta is now in chemical equilibrium. In our molecular analysis, we consider only this single snapshot of the evolution. Fig. 1 shows the composition profiles in this model.

• RSG: This model is based on a 15|${\rm M}_{\odot }$| RSG progenitor which explodes to produce a typical Type IIP SN (examples include SN1999em, SN2004et, and SN2017eaw (Leonard et al. 2002; Kotak et al. 2009; Rho et al. 2018). Specific details of the model can be found in Hillier & Dessart (2019), from which we use the ‘x2p0’ model. This model is calculated from a progenitor which retains a significant (⁠|$\gt 9{\rm M}_{\odot }$|⁠) hydrogen envelope at the time of explosion. The explosion model produces 11.87 |${\rm M}_{\odot }$| of ejecta, of which 0.036 |${\rm M}_{\odot }$| is |$^{56}$|Ni, which provides the heating at late times. After explosion, the composition is mixed using a boxcar algorithm, with weaker mixing than in the BSG model. We do not include CO until 50 d after the explosion. At early times, the high temperatures in the oxygen-rich region should result in very little molecular formation. 50 d was chosen so that CO begins to form before |$\sim 100$| d, per the earliest detections of CO in SN1987A (Spyromilio et al. 1988). The model is evolved until 300 d post-explosion, including all time-dependent terms and the effect of CO cooling on its evolution. Fig. 2 shows composition profiles for the RSG explosion model.

Both models are computed from a 15 |${\rm M}_{\odot }$| progenitor but with different explosion properties, notably the density structure and progenitor radius at the time of the explosion (Fig. 3).

4 RESULTS AND DISCUSSION

4.1 CO formation

In our RSG model, |$1.76 \times 10^{-4}\,{\rm M}_{\odot}$| of CO has formed by 330 d, while our BSG model (at 330 d) contains |$6.1 \times 10^{-5} {\rm M}_{\odot }$| of CO. In both models, CO forms preferentially in conditions of high density and low temperature. It is formed in the inner regions of the ejecta (<2000 km s|$^{-1}$| in our BSG model, <5000 km s|$^{-1}$| in our RSG model), with the majority of CO forming after |$\sim$|100 d post-explosion in our RSG model. The formation of CO depends strongly on temperature: in regions where the temperature exceeds |$4000{-}4500$| K (depending on the local density), negligible quantities of CO form, as it is quickly destroyed by non-thermal electrons and charge exchange with oxygen. In lower temperature regions, CO can form in significant quantities, though the local number density of CO never exceeds |$\sim$|1 per cent of the carbon number density in any region within our models. A temperature below 4000K is not sufficient for CO to form in large quantities; relatively high densities are also required (though we do not explore enough of the parameter space to make a quantitative estimate of a density threshold). We do not observe any scenario in which most of the carbon or oxygen becomes locked up in molecules, even at late (⁠|$\ge 300$| d post-explosion) times and low (⁠|$\sim 1000$|K) temperatures. The equilibrium quantity of CO increases as density increases or temperature decreases.

The equilibrium population of CO is dictated by a few reactions. The radiative association reaction |$\mathrm{C} \,\, + \mathrm{O} \,\, \rightarrow \,\, \mathrm{CO}\,\, + \,\, \gamma$| typically dominates the formation of CO, though the charge exchange reaction |$\mathrm{CO}^+ \,\, + \,\, \mathrm{O} \,\, \rightarrow \,\, \mathrm{CO}\,\, + \,\, \mathrm{O}^+$| can become important when the density of neutral oxygen is high. Charge exchange with oxygen: |$\mathrm{CO}\,\, + \,\, \mathrm{O}^+ \,\, \rightarrow \,\, \mathrm{CO}^+ \,\, + \,\, \mathrm{O}$| typically dominates the destruction of CO, but dissociation by non-thermal electrons (⁠|$\mathrm{CO} + {\rm e}^{-}_{\rm non-therm} \rightarrow \mathrm{C}/\mathrm{C}^+ + \mathrm{O}/\mathrm{O}^+$|⁠) can also contribute significantly, especially at lower densities. CO ionization balance is dictated by charge exchange with O and |$\mathrm{O}^+$|⁠, as the reaction rates exceed by orders of magnitude the CO photoionization/recombination rate.

Our prescription for chemical mixing included in both models likely impacts where and how much CO forms. The chosen boxcar approach allows CO to form throughout the ejecta, following the distribution of carbon and oxygen. The shuffled-shell approach, while more realistic, requires higher spatial resolution and thus longer computation times. Our initial results using boxcar mixing likely overestimate the spatial extent of CO formation, but possibly underestimate the local density of CO that could form in dense, carbon/oxygen-dominated shells. Boxcar mixing also means that CO is distributed throughout the ejecta with all other elements, so its cooling effect impacts not only carbon and oxygen but also hydrogen and other species. With a stratified shell model, we would expect CO to affect only the species in the carbon/oxygen-rich shell, leaving the outer, hydrogen-rich layers unaffected.

4.1.1 Formation over time

In our RSG model, CO forms in two distinct locations, with a region of high temperature and low CO density separating them. In both regions, cooling due to CO lines constitutes a significant portion (⁠|$\gt 25~{{\ \rm per\ cent}}$|⁠) of the total cooling. Most (around three-fourth of the total mass) of the CO forms in the oxygen/silicon-rich region at velocities <1500 km s|$^{-1}$| with the remaining CO spread over the broader oxygen/carbon-rich region extending from 2500 to 4500 km s|$^{-1}$|⁠. The first region is not as cool but is denser, illustrating the conditions that lead to CO formation: high density, low temperature, or both. The CO band emission in the spectrum is dominated primarily by emission from the lower velocity region, where most of the CO mass is formed. The separating region has intermediate density but lower carbon abundance and higher temperature, reducing the amount of CO that forms. By nebular times, the two CO regions have cooled to a significantly lower temperature than the separating region. Fig. 4 shows the evolution of CO in our RSG model as a function of velocity and time. In our BSG model, all the CO is formed in the central regions, and the CO density quickly becomes negligible at velocities >1500 km s|$^{-1}$|⁠.

The amount of CO grows slowly between 50 and 100 d, after which its growth begins to accelerate (Fig. 5). The growth of CO then levels off at around 130 d post-explosion. The earliest formation of CO (before 100 d) occurs primarily in the cooler formation region between 3000 and 4000 km s|$^{-1}$|⁠, where the oxygen ionization fraction is lower than in the hotter region below 1500 km s|$^{-1}$|⁠. Ionization of oxygen and carbon dictates when CO forms, as reactions involving neutral carbon and oxygen dominate CO formation. The initial burst of CO formation occurs as oxygen recombines, corresponding to the plateau drop-off of the SN light curve when the photosphere has moved through the hydrogen envelope. The predominant formation pathway at this epoch is radiative association of neutral carbon and oxygen. Carbon remains ionized for longer, causing the amount of CO to plateau between 130 and 220 d. After 220 d, neutral carbon starts to become the dominant ionization stage, allowing for more CO formation. At this stage, the growth of CO is also due to a pair of formation reactions which involve neutral oxygen as a reactant (reactions 13 and 27 in Table A1). At 300 d after explosion, the total mass of CO in our RSG model is around |$10^{-4} {\rm M}_{\odot }$|⁠.

4.2 Spectral features

The addition of line emission from CO adds a new efficient cooling channel, causing the nebular-phase equilibrium temperature in the CO region to be lower than without CO, by as much as a factor of 2 (Fig. 6). While the SN is still in the photospheric phase, CO emission does not significantly affect the spectrum, as very little CO forms before the end of the plateau phase. As the SN evolves to the nebular phase, the impact of CO cooling becomes more pronounced: the lower temperature in the CO region reduces emissivity in the lines formed in that region, notably oxygen and calcium lines. In our mixed models, the reduced CO-region temperature depresses emission from all elements, but in more realistic unmixed models, it would be expected to impact only lines formed in the oxygen/carbon-rich region.

Optical spectra calculated from our BSG models at the nebular phase with and without CO are shown in Fig. 7. The impact of CO can be seen in all the lines, but is most pronounced in the strongest lines, including H|$\alpha$| and Ca ii: the effect on these features is highly dependent on our imposed mixing. The vast majority of CO emission emerges in the first overtone and fundamental bands, with most of the emission in the fundamental band at |$\sim 4.5{-} 5.5\,{\mu \rm m}$| (Fig. 8). Even with less than 1 per cent of carbon atoms in the form of CO, emission from the CO fundamental band accounts for as much as 20 per cent of the total luminosity of the CO model. The addition of cooling from CO lines means the ejecta cool faster than without CO, resulting in a faster declining light curve at optical wavelengths (Fig. 9). The radiation that would have been emitted in the optical is instead emitted in CO rovibrational bands in the infrared.

Broad-band filters which overlap the CO emission in the infrared (Fig. 8) show significantly more flux at late times; brightness in the F444W NIRCam filter actually rises after 120 d as CO fundamental band emission becomes more important. The K band and some JWST bands overlap with CO emission bands in the infrared, and should be sensitive to CO emission in SNe. Mid-infrared photometric observations using Spitzer and WISE have shown evidence for CO emission in multiple SNe (Gandhi et al. 2013; Ergon et al. 2014), and JWST’s greater sensitivity should improve constraints on CO emission. Most existing observational constraints on CO emission come from the first overtone, as in, for example, SN2017eaw (Rho et al. 2018). JWST, however, can observe both the first overtone and fundamental bands of CO emission with high resolution; its upcoming observational programs could substantially improve our understanding of CO emission in SNe. The relative strength of the first overtone and fundamental bands depends on the optical depth. Constraining both bands simultaneously can improve the accuracy of observational estimates of CO mass (see Section 6.1).

Fig. 10 shows our RSG model spectrum at 210 d compared to a similar-epoch near-infrared spectrum of SN2017eaw (from Rho et al. 2018). Note that the CO bands in the emergent spectrum are dominated by emission from the lower velocity formation region in our RSG model. Rho et al. (2018) estimated a CO mass of |$\sim 2 \, \times \, 10^{-4} \, {\rm M}_{\odot }$| from the first overtone emission at this epoch, roughly three times the CO mass present in our RSG model at 210 d. Despite this discrepancy, our model produces more emission in the first overtone than in the observations. The relative bandhead strength in our model also differs from the observations, with the model overestimating the relative strength of the longer wavelength bandheads (which come from higher energy vibrational states). Our models (with and without CO) significantly overestimate the strength of several Mg i lines present in this wavelength region, which could be due to errors in ionization or magnesium abundance. Magnesium lines in the infrared are sensitive to a number of progenitor properties, any of which could contribute to the discrepancy between our model and SN2017eaw. See Dessart & Hillier (2020b) for a discussion of how properties including the progenitor hydrogen envelope mass, the extent of chemical mixing, or the progenitor ZAMS mass affect the ionization in the ejecta and the appearance of Mg i lines.

4.3 Clumping and dust

Clumping in SN ejecta consists of small regions of enhanced density (clumps) separated by much lower density material in the interclump region. The total amount of CO formed in the ejecta depends strongly on density – a clumpier ejecta profile would likely produce more CO and cool faster than uniform ejecta. Dessart, Hillier & Wilk (2018) investigated the effects of clumping in a similar cmfgen model of a Type II-peculiar SN, finding that the inclusion of clumping resulted in a faster evolving light curve, redder colours, and an earlier transition to the nebular phase. The inclusion of both CO cooling and strongly clumped ejecta would likely allow for earlier formation of CO and exacerbate both the effects of clumping described in Dessart et al. (2018) and the effects of CO cooling inferred in this work. In fact, results from multidimensional simulations suggest the ejecta should be most strongly clumped (i.e. higher densities within clumps) in the inner regions of the ejecta where CO is most important as a result of local expansion driven by ‘nickel bubbles’ (Mueller, Fryxell & Arnett 1991; Wongwathanarat et al. 2015).

The presence of CO has implications for dust formation by cooling the ejecta and potentially affecting the abundance of carbon atoms available to form dust. The inclusion of CO speeds the rate of cooling in the densest ejecta, allowing for the formation of dust earlier. Our model also predicts CO abundances significantly smaller than the total carbon or oxygen abundance (<1 per cent), leaving both species available to form larger molecular chains, i.e. dust.

5 CO COOLING BEHAVIOUR

5.1 CO cooling function

The CO cooling function has a non-monotonic dependence on temperature, due to the interaction of processes that depend on temperature in different ways. CO molecules form preferentially at low temperatures, with the equilibrium CO abundance increasing by a factor of |$\sim 30$| when the temperature decreases from 4000 to 2000 K. However, CO emissivity per molecule is an increasing function of temperature, roughly proportional to |${\rm e}^{-h \nu /k_B T}$|⁠. The total CO cooling function, the product of these two processes, has a noticeable peak and decreases at higher and lower temperatures (Fig. 11). The exact location of the peak also depends on the local density (density, temperature, and composition being the primary drivers of molecule formation), with the peak shifting to higher temperatures as density increases. We compute the total CO cooling rate by converging a model at fixed temperature, then summing the net radiative bracket times the photon energy (i.e. the net emission in the line) over all CO lines:

where |$N_u$| is the upper state population, |$Z_{line}$| the net radiative bracket, |$\bar{J}_{line}$| the frequency-averaged mean intensity, and |$S_{line}$| the line source function. This represents the total cooling by CO line emission, taking into account non-local processes such as line overlap and radiative transfer effects.

The non-monotonic CO cooling function causes difficulties for finding the temperature solution in our models. Some of the complications that can arise due to the shape of the CO cooling curve can be seen in Fig. 11.

Without CO cooling included, the total cooling curve monotonically increases with temperature. When CO is included, it dominates the cooling behaviour at temperatures around 2000–3000 K, forming a local maximum near 2000  K and a local minimum near 3200 K. The new total cooling curve now crosses the heating curve at up to three points, each of which is a potentially valid temperature solution, labelled |${\rm T}_{\rm {A}}$|⁠, |${T}_{\rm {B}}$|⁠, and |${T}_{\rm {C}}$| in Fig. 11. We can reject solution B on the grounds of stability, but solutions A and C are both stable. Solution C lies closest to the temperature solution if CO cooling is not included, and could be considered the most physically realistic solution, as it is close to Spyromilio et al. (1988)’s estimate of the CO excitation temperature (⁠|$\sim$|3000 K) in SN 1987A.

Multiple temperature solutions (hereafter scenario 1) can arise in models that include molecular cooling (Fig. 11). Solutions A and C in Fig. 11 are both physically valid; recall that the radiative-transfer problem comes with no guarantee of uniqueness. We expect that SN ejecta will evolve from high temperature to low temperature, so the initial temperature should be higher than |$\rm {T}_{\rm {C}}$| and the gas should cool to that temperature without issue. However, the CO cooling phenomenon can cause more serious problems for our Newton–Raphson solution technique if the quantity of CO is large.

At some regions in our models, the density of CO is high enough that the local minimum in cooling at |$\sim$|3000 K is still greater than the heating rate, leaving solution A as the only temperature solution (scenario 2). Depending on our initial estimate of temperature, our iterative Newton–Raphson solver can become trapped in this local minimum and fail to converge to a temperature solution without manual intervention. Our initial temperature estimates, which typically come from similar models that do not include molecular cooling, lie at temperatures above |$\mathrm{T}_\mathrm{C}$|⁠, and early iterations tend to fall into the local minimum. Automating a correction for this effect, which would attempt to locate solution A, is non-trivial, as the location of the CO cooling peak depends on the local density, oxygen, and carbon abundances, and the local radiation field. On the other hand, it may be the case that the density of CO is low enough that the local maximum near 2000 K lies below the heating curve (scenario 3). In this scenario, only solution C is valid, and the local minimum is not a concern. In this case, the CO cooling effect is small, and the final temperature solution is similar to that without CO.

5.2 CO abundance and the cooling function

The separation among the scenarios discussed above is not particularly large. Because the total cooling curve does not vary strongly with temperature between 2000 and 4000 K, small changes in CO abundance can result in large changes in the inferred equilibrium temperature. If the CO abundance is increased by 50  per cent, the cooling curve in Fig. 11 would shift to scenario 2, leaving only |$\rm {T}_A$| as a temperature solution regardless of the initial conditions. Conversely, a 50  per cent decrease in CO abundance moves it to scenario 3, leaving |$\rm {T}_C$| as the only solution. An equivalent change in the radioactive decay heating rate can also change the solution scenario. Thus a factor of |$\sim$|2 change in the CO abundance, which is within the range of uncertainty in molecular reaction rates, separates the regime in which CO cooling is unimportant and the equilibrium temperature is |$\sim$|4000 K, to a regime in which CO cooling dominates and the equilibrium temperature is closer to 1500 K.

The equilibrium CO abundance depends sensitively on the rates of the most significant reactions that create and destroy it, notably radiative association (reaction 37 in Table A1, formation) and charge exchange with oxygen (reaction 9 in Table A1, destruction). If one or both reaction rates are off by |$\sim$|50  per cent, it can alter not only the final temperature solution but the number of valid temperature solutions. This is a significant concern, considering that most reaction rates available through UMIST (McElroy et al. 2013) are only accurate to within 25–50 per cent. Recent experimental results by Meng, Qin & Liu (2022) suggest the rate of radiative association may be as much as an order of magnitude higher than the UMIST rate we use in this work. Since this reaction typically dominates CO formation, an increased rate would significantly increase the amount of CO and associated cooling in our models. Using the radiative association reaction rate calculated by Meng et al. (2022) increases the CO abundance by a factor of 10 and reduces the equilibrium temperature in the CO-dominated region of our BSG model by |$\ge$|500 K compared to the UMIST rate (McElroy et al. 2013). Obtaining highly accurate reaction rates is thus of paramount importance, as deviations of a factor of 2 or less can change the temperature solution by thousands of Kelvin or introduce additional valid solutions.

Our assumptions of steady state and LTE level populations likely contribute to this issue. In a full time-dependent model, we expect the ejecta to evolve from high temperature to low temperature, and to settle into |$\rm {T}_C$| without issue if the CO abundance is not too high. If CO-level populations are allowed to deviate from LTE, then the net cooling from CO is reduced, potentially alleviating the problem of multiple solutions and leaving only the |$\rm {T}_C$| solution. Some of our other assumptions, however, may be causing us to underestimate the importance of CO cooling. For example, unmixed models would be expected to produce higher CO densities in the oxygen/carbon shell, increasing the importance of CO cooling in that region and exacerbating this issue. Including strong clumping in the ejecta would also be expected to result in increased CO abundance and potentially lower temperature solutions. Alternatively, a change in our assumptions about how molecular reaction rates depend on the energy states of reactants (see Section 2.1.1) could have a sizeable impact on the CO abundance and subsequently the equilibrium temperature if a substantial portion of carbon or oxygen atoms are not in the ground state (as is likely the case in ionized gas).

While our BSG model suggests a total CO mass |$\sim 1 \times 10^{-4}\,{\rm M}_{\odot}$| at 330 d, observational estimates vary (Section 6.1), mostly |$\gtrsim 10^{-4}\,{\rm M}_{\odot}$| (Spyromilio et al. 1988; Liu & Dalgarno 1995), slightly above our model estimates (though not inconsistent, as the observational estimates come with large uncertainties and we utilize simple proof-of-concept models).

6 CONTEXT

This work is intended primarily to demonstrate the developments made to treat molecules in cmfgen, not reproduce observations, so detailed comparison to observations will be revisited in future works that make use of more realistic models. However, we make some comparisons both to observations of SN1987A and to other models of molecular formation to ensure that our results are not unreasonable.

6.1 CO observations in SN1987A

The first detection of CO in an SN was in SN1987A, which exploded from a BSG progenitor in the Large Magellanic Cloud in 1987 February. Spectra of SN1987A taken as early as fall 1987 and spring 1988 showed evidence for CO fundamental band and first overtone emission, as well as SiO emission, with subsequent analysis deducing that the molecules must have formed within the ejecta (Arnett et al. 1989). Spyromilio et al. (1988) estimated the total mass of CO in the ejecta using a simple model which assumed a uniform, optically thin sphere of CO emitting in LTE. We compare a similar model to the LTE first overtone emission in our BSG model at 330 d. The simple model, computed with total CO mass, maximum velocity, and excitation temperature chosen to match the properties of our BSG model, closely replicates the appearance of the first overtone band, suggesting that the first overtone region is optically thin in our model (see Fig. 12). It does not accurately reproduce the fundamental band emission, however, which is optically thick in our BSG model. When the assumption of LTE in our CO model is removed, the simple model fails to fit the appearance of either band – at this phase, the CO-level populations are far from their LTE values.

Liu et al. (1992) extended this simple model by also considering models without the assumptions of LTE populations or optical thinness. If CO populations deviate from LTE, the relative strength of shorter wavelength bandheads, which correspond to radiative de-excitations from lower energy vibrational states, should increase compared to longer wavelength features (existing observations of SN1987A cover only the shortest wavelength bandheads in the first overtone band). In addition, higher optical depth in CO band regions can reduce the total emissivity at fixed CO mass, as photons are trapped until the transition is collisionally de-excited. This introduces a degeneracy when modelling only the first overtone band, as the feature could be produced by a large mass of CO with a high optical depth or a low mass with low optical depth. Using both the first overtone and the fundamental bands can potentially break this degeneracy, as the relative strength of the two should depend on the optical depth. In the optically thin limit, the total flux in the CO fundamental band should exceed the flux in the first overtone by about a factor of 20.

In our models, regardless of the assumption of LTE, the fundamental band region is always optically thick (⁠|$\tau \sim 2-5$|⁠)⁸ and the fundamental/first overtone flux ratio is |$\lt 20$|⁠. Observations of SN1987A at 260 d post-explosion from Meikle et al. (1989) show a fundamental/first overtone flux ratio of |$\sim 3$|⁠, suggesting that the optical depth in the fundamental band region is very high.

When using the same assumptions as in Spyromilio et al. (1988), Liu et al. (1992) deduced significantly (5 to 10 times) higher CO masses at the same epochs, though the reason for the discrepancy is unclear. Our BSG model produces |$\sim 10^{-4}\,{\rm M}_{\odot}$| of CO at 330 d (regardless of the assumption of LTE), which is consistent with the Liu et al. (1992) estimate of the CO mass in SN1987A at 349 d when assuming LTE in an optically thin gas (and roughly consistent with the Spyromilio et al. (1988) estimate at the same epoch), but well below their estimate of |$1.1 \times 10^{-3}\,{\rm M}_{\odot}$| of CO at the same epoch when the assumption of LTE is removed. The total mass of CO present in our BSG model is the same regardless of our assumption of LTE, but the strength of the CO bands changes dramatically (Fig. 12). Applying the method of Spyromilio et al. (1988) to the CO first overtone flux in our BSG LTECO model gives a mass estimate of |$6.1 \times 10^{-5} {\rm M}_{\odot }$|⁠, matching the mass of CO present in our model. Fitting our NLTE model flux with the same method suggests a CO mass of |$2.0 \times 10^{-5} {\rm M}_{\odot }$|⁠. Mass measurements made using this method are likely underestimated if the level populations are out of LTE.

When out of LTE, collisional processes determine CO-level populations. Unfortunately, existing collisional data for CO rovibrational levels only cover vibrational transitions from |$v=0$|⁠, requiring the use of theoretical scaling relations to cover all possible transitions (Chandra & Sharma 2001; Ogloblina et al. 2020). Using this data, our NLTE models fail to reproduce the first overtone band appearance in either SN19897A or SN2017eaw, with our models consistently overestimating the strength of higher order bandheads. This suggests that our treatment of NLTE-level populations requires more accurate and comprehensive collisional data. Other features of the chosen model, such as progenitor mass, nature of chemical mixing, or composition may also contribute to the appearance of the CO bands.

Alternatively, the first overtone band appearance in SN1987A (Fig. 12, Meikle et al. 1989) can be fit by an LTE model with lower temperature, as in Spyromilio et al. (1988) and Liu et al. (1992). At lower excitation temperatures, higher vibrational levels will be less populated, reducing the strength of higher order bandheads – i.e. de-excitations from higher vibrational states, which emit at longer wavelengths – and increasing the relative strength of the lowest order bandheads (at shorter wavelengths). A simple model with an excitation temperature around 2000 K and |$2 \times 10 ^{-4} {\rm M}_{\odot }$| of CO fits the appearance of the first overtone at 260 d quite well (Liu et al. 1992). This temperature is also close to the very low-temperature solution we obtain in the CO-dominated region of our BSG LTECO model (Section 5). However, such a model also predicts an even higher ratio (⁠|$\sim$|50) of fundamental band to first overtone band flux, which is in tension with the observed ratio (⁠|$\sim$|3) unless the fundamental band optical depth in 1987A is very high.

6.2 Chemical models

Several studies have investigated CO formation in SN ejecta using chemical evolution models, including Sluder et al. (2018), Petuchowski et al. (1989), and Liu & Dalgarno (1995). Table 2 shows a non-exhaustive list of CO masses computed using a variety of chemical evolution models for SN environments analogous to SN1987A. Recent estimates of the CO mass tend to be somewhat higher than early estimates as a result of improved reaction networks and reaction rate data. The earliest models in the table, Petuchowski et al. (1989) and Lepp et al. (1990), assume fully mixed ejecta and compute significantly lower CO masses due to charge exchange with He|$^+$|⁠, which can quickly destroy CO molecules if present in the same region. More recent works (including this one) neglect this reaction as a result of chemically stratified ejecta in which He and C/O reside in distinct shells, producing significantly more CO (Lepp et al. 1990; Gearhart, Wheeler & Swartz 1999).

We calculate somewhat lower total CO masses than contemporaneous works, which can be attributed partially to differences in chemical reaction networks, but is more likely a result of the strong mixing enforced in our BSG model, which results in lower C and O mass fractions in the inner regions where CO dominates as compared to an unmixed or macroscopically mixed model. Our NLTE BSG model predicts lower flux in the first overtone band than is observed in SN1987A, suggesting that our model could be underestimating the CO production. On the other hand, the same NLTE model still overestimates the strength of the CO fundamental band. Even |$10^{-4}\,{\rm M}_{\odot}$| of CO in our BSG model (total ejecta mass 13.2 |${\rm M}_{\odot }$|⁠) is significant, with 20 per cent of the total flux emitted in CO bands in the LTE model, while the NLTE model sees the CO emission reduced by about a factor of 10 relative to the LTE model. This mass of CO means that, in the C/O-rich region of the ejecta, less than 1 per cent of C atoms are trapped in CO molecules, which allows for the later production of carbonaceous dust.

The most directly comparable studies to ours are those performed using sumo, described in Liljegren et al. (2020) and Liljegren et al. (2022), which use a similar reaction network and include non-thermal and NLTE processes. We calculate lower CO masses than Liljegren et al. (2020), most likely as a result of the model structure: they utilize a one-zone model consisting of only carbon and oxygen, while our BSG model features carbon and oxygen mixed throughout the ejecta with other elements. However, our individual reaction rates agree well with their results (with the exception of CO charge exchange with |$\mathrm{O}^+$|⁠, which we find contributes significantly even after 200 d) and we find a similar impact from CO cooling, which reduces the temperature in the CO-dominated region (our models see a smaller cooling effect as a result of less CO). However Liljegren et al. (2020) do not report the same temperature solution issues described in Section 5 that result from CO cooling in our models. We also confirm the need for more accurate reaction rates, as the CO mass and equilibrium temperature depend strongly on the assumed rates: Liljegren et al. (2020) found that uncertainty in CO abundance scales roughly linearly with uncertainty in reaction rates. Since most rates in UMIST (McElroy et al. 2013) are accurate to about 50  per cent, our models can only constrain the CO abundance within 50  per cent, which covers all three temperature solution scenarios (Section 5).

7 CONCLUSIONS

We have included a molecular reaction network in the radiative-transfer code cmfgen to study the formation of CO and its impact on the gas and radiation properties in SN ejecta. We account for the formation and destruction of CO by a number of chemical reactions and its impact on the ejecta via cooling by ro-vibrational emission. We investigate CO formation and cooling in an SN II-pec (BSG) model at 330 d post-explosion and an SN IIP (RSG) model from 50 to 300 d post-explosion, both of which feature microscopic mixing throughout the ejecta. We confirm that CO is an important coolant at late times in oxygen- and carbon-rich material.

In both models, CO forms in the densest and coolest regions of the ejecta, limited to low velocities (<2000 km s|$^{-1}$|⁠), in agreement with observations of CO. In our RSG model, CO becomes important (i.e. contributes a significant fraction of the total cooling) at |$\sim$|100 d after the explosion, and by |$\sim$|200 d after explosion, CO line cooling can constitute as much as 50 per cent of the total cooling in some regions. The primary effect of CO cooling is to reduce the temperature by as much as two thousand Kelvin in CO-dominated regions. A reduced temperature weakens emission lines which are produced in the CO-dominated regions (due to the mixing in our models, we likely underestimate the impact on oxygen lines but overestimate the impact on lines from other species) and shifts that luminosity mostly to the CO fundamental band at |$\sim 4.5\,\mu$|m. In this band, the ejecta remain optically thick well past 300 d. The CO first overtone band at |$\sim 2.5\,\mu {\rm m}$| is nearly optically thin at 300 d – any estimates of CO mass made using these bands should take optical depth effects into account.

In our BSG model, CO dominates the total cooling in the inner regions. Both LTE and NLTE treatments of CO-level populations significantly overestimate the ratio of fundamental band flux to first overtone band flux compared to observations of SN 1987A, potentially indicative of high optical depth in the CO fundamental band in SN1987A. Allowing CO-level populations to deviate from LTE is necessary to reproduce the first overtone band structures in SN1987A (II-pec) and 2017eaw (IIP) (Meikle et al. 1989; Rho et al. 2018). Our LTE models overestimate the strength of higher order bandheads compared to observations of both SNe, suggesting the higher vibrational levels are underpopulated relative to LTE. Unfortunately, collisional excitation data for CO is limited, covering only a subset of possible transitions (Ogloblina et al. 2020).

CO cooling presents numerical issues when determining the ejecta temperature structure. CO forms more efficiently in low-temperature environments, causing the total CO population to be a decreasing function of temperature. However, its emission per molecule is an increasing function of temperature. As a result, the CO cooling curve shows a strong peak between 2000 and 3000 K. In dense regions of the ejecta where CO cooling is important, the total cooling function becomes non-monotonic with a local minimum near 3000–4000 K. The presence of a local minimum can trap numerical solvers or result in multiple temperature solutions.

The total amount of CO produced depends sensitively on the rates of a few formation and destruction mechanisms, predominantly radiative association (formation) and charge exchange with oxygen (destruction). Literature values for the rates of these reactions come with large uncertainties (>25 per cent) and different sources report different rates: the rate Dalgarno, Du & You (1990) measured for radiative association differs by an order of magnitude from the rate measured by Meng et al. (2022). Accurate rates are crucial for correctly estimating the mass of CO and the impact of CO cooling, and different rates can result in very different temperature solutions as a consequence of the non-monotonic CO cooling function (Section 5).

Using the new capabilities of cmfgen, we will investigate the effects of CO cooling in SNe in more detail in future work. We plan to include CO in full time-sequence models of Type II-Plateau, Type II-peculiar, and Type Ibc SNe and explore macroscopic mixing and clumping. In addition, other potentially important molecules like SiO and carbon monosulfide (CS), both of which have been observed in SN1987A (Meikle et al. 1993), will be added to cmfgen using the same approach detailed in this work.

Acknowledgement

CM and DJH gratefully acknowledge partial support through NASA astro-physical theory grant 80NSSC20K0524. DJH also acknowledges partial support from the STScI through grant HST-AR-16131.001-A.

DATA AVAILABILITY

cmfgen is available on Github or through the cmfgen website at https://sites.pitt.edu/~hillier/web/CMFGEN.htm. The molecular version (including all molecular data utilized in this work) can be found on Github at https://github.com/collinmcleod/cmfgen_mol. Specific models utilized in this paper are available upon request.

Footnotes

REFERENCES

APPENDIX A: REACTION DATA