# Carbon monoxide cooling in radiative transfer modelling of supernovae

Collin McLeod<sup>1,★</sup> D. John Hillier<sup>1</sup> and Luc Dessart<sup>2</sup>

Accepted 2024 June 14. Received 2024 June 14; in original form 2024 February 2

## 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⁻<fnref n="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 ~4.5 µ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.

**Key words:** molecular processes – radiative transfer – (stars:) supernovae: general.

## 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 (~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).<fnref n="1" /> 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 γ-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)<fnref n="2" /> 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 (> 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 ⁵⁶Ni and ⁵⁶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 (≤ 0.5 eV) energy separations, and transitions between them emit photons in the infrared range.

<footnote n="1"><sup>1</sup>Department of Physics and Astronomy & Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O'Hara Street, Pittsburgh, PA 15260, USA</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, CNRS-Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote>* E-mail: clm224@pitt.edu</footnote>

<footnote n="1">¹In some SNe, the late-time spectra may be powered by shock heating from ongoing interaction with circumstellar medium (CSM).</footnote>

<footnote n="2">²Some authors use ‘CO’ when referring to the carbon/oxygen-rich regions or the carbon/oxygen core of the progenitor. To avoid ambiguity, we use ‘CO’ solely to refer to carbon monoxide, and ‘C/O’ to refer to carbon and oxygen.</footnote>
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<fnref n="3" /> (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} 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} 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} 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} 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

<footnote n="3">³In the absence of interaction with CSM (Dessart & Hillier 2022)</footnote>
(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

$$ \rho \frac{D n_i}{D t} = \sum_j (n_j R_{ji} - n_i R_{ij}), \quad (1) $$

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 OI 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:

$$ \sum_{j=1}^{N} \frac{D N_j}{D t} = 0 \quad (2) $$

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

⁴Fluid instabilities during the explosion produce large bubbles or 'fingers' of metal-rich material from the core (which later expand as a result of heating from radioactive ⁵⁶Ni and ⁵⁶Co). These ⁵⁶Ni-rich regions break spherical symmetry and extend out to large velocities but may not introduce significant microscopic mixing between layers (Fryxell, Mueller & Arnett 1991; Kifonidis et al. 2003; Wongwathanarat, Müller & Janka 2015). Observed line profiles in many objects suggest minimal mixing between ejecta shells: for example, microscopic mixing of calcium into the O-rich shell would produce weaker oxygen lines and stronger calcium lines than are typically observed (Dessart & Hillier 2020b; Dessart et al. 2021). Direct imaging of SN1987A also indicates little mixing between C/O- and Si/O-rich shells (Abellán et al. 2017). For further discussions of mixing (which is not the focus of this work), see Dessart & Hillier (2020a), Ono et al. (2023), and Dessart & Hillier (2020b).

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:

$$ r = \alpha n_{R_1} n_{R_2} \left( \frac{T}{300 K} \right)^{\beta} \exp(-\gamma/T), \quad (3) $$

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:

(i) 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 (O1), this is the $2s^22p^{43}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.

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

(iii) Potential atomic product states are populated proportional to their statistical weights. For example, the reaction $C_2 + O^+ \rightarrow C_2^+ + O$ populates all three J states of the OI ground term: it populates the $J=2$ state ($g=5$) at $5/9$ of the total reaction rate, the $J=1$ state ($g=3$) at $3/9$ the total rate, and the $J=0$ state ($g=1$) at $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.<fnref n="5" />

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 $C_2 + O \rightarrow C + CO$. The species $C_2$, $O_2$, $CO_2$, and $C_2O$ 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

<footnote n="5">⁵In the CO-dominated region of our BSG model, radiative association dominates the formation of CO. The cooling rate (enthalpy change) from this reaction is ~ $5 \times 10^{-12}$ erg s⁻¹ cm⁻³, compared to a CO line cooling rate of ~ $1.6 \times 10^{-7}$ erg s⁻¹ cm⁻³ and a net cooling rate of ~ $3.7 \times 10^{-9}$ erg s⁻¹ cm⁻³. This reaction exceeds by at least a factor of 100 the cooling rate from any other molecular reaction.</footnote>
<table><caption>**Table 1.** Compton (non-thermal) reaction processes included in CMFGEN.</caption><thead><tr><th>ID</th><th>Reaction</th><th></th></tr></thead><tbody><tr><td>1</td><td>CO + e<sup>+</sup><sub>non-therm</sub> → C + O + e<sup>-</sup></td><td></td></tr><tr><td>2</td><td>CO + e<sup>+</sup><sub>non-therm</sub> → C<sup>+</sup> + O + e<sup>-</sup> + e<sup>-</sup></td><td></td></tr><tr><td>3</td><td>CO + e<sup>+</sup><sub>non-therm</sub> → C + O<sup>+</sup> + e<sup>-</sup> + e<sup>-</sup></td><td></td></tr></tbody></table>

*Notes.* Data for cross-sections come from Itikawa (2015), by way of the LXCAT database (Pancheshnyi et al. 2012). Data were extrapolated using a power law to match the sampling of the non-thermal electron spectrum in CMFGEN.

we develop our method. Future work will include a much larger network with additional molecules.<fnref n="6" />

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) ⁵⁶Ni and ⁵⁶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

$$r = n_{CO} \int_{0}^{\infty} \sigma(v) v \frac{dn_{e}}{dv} dv, \qquad (4)$$

where $n_{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{dn_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 (¹²C¹⁶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 ≈5000 K, the first excited electronic state of CO, $a^3\Pi$, which lies ~6 eV above the ground state, would have a Boltzmann factor of ~10⁻<fnref n="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 \to v'$, 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):

$$\frac{C(T|v \to v_1)}{C(T|v \to v_2)} = \frac{2v_1 + 1}{2v_2 + 1} \qquad (5)$$

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.).<fnref n="7" /> 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:

$$n_{\text{carbon}} = \sum_{i=0}^{l_{\text{max}}} n_{\text{C}^{i+}} + n_{\text{CO}} + 2n_{\text{C}_2} + n_{\text{CO}_2} + 2n_{\text{C}_2\text{O}}. \qquad (6)$$

## 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 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 M$_{\odot}$ progenitors, evolved with the stellar evolution code MESA (Paxton et al. 2018) and then exploded using the radiation hydrodynamics

<footnote n="6">⁶Our results from this study indicate that the computation time for chemical reactions is insignificant compared to the radiative transfer calculation, so a larger reaction network should not significantly impact the total program runtime in future work.</footnote>

<footnote n="7">⁷Via radioactive decay, some elements can change: ⁵⁶Ni decays to ⁵⁶Co over time. This process, however, has a fixed half-life and does not depend on the local conditions in the ejecta, while the population of CO does.</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Mass fractions of important elements at 330 d post-explosion in the BSG model, which is microscopically mixed using a boxcar. This model has a total ejecta mass of 13.22 M☉, of which 7.31 M☉ is hydrogen and 0.0843 M☉$^{56}$Ni. The interior of this model is dominated by helium and oxygen, with a smaller proportion of carbon.</figcaption></figure>

code V1D (Livre 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.

(i) BSG: This model is based on a 15 M☉ 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 M☉) 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.

(ii) RSG: This model is based on a 15M☉ 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 (> 9M☉) hydrogen envelope at the time of explosion. The explosion model produces 11.87 M☉ of ejecta, of which 0.036 M☉ 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 ~ 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

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Mass fractions of important elements at 46 d post-explosion in the RSG model, which is microscopically mixed using a boxcar. The innermost region is dominated by helium, oxygen, and silicon, followed by a shell of primarily helium, oxygen, and carbon. The outermost ~ 9M☉ is dominated by hydrogen.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Density profiles of the BSG and RSG models at the same epoch (300 d). Both models have a power-law density distribution at high velocities, but they differ at low velocities. Near the inner parts of the ejecta, the RSG model is more stratified, with a dense C/O shell. The primary difference between the models is that the RSG progenitor has a much larger radius at the time of the explosion.</figcaption></figure>

composition profiles for the RSG explosion model.

Both models are computed from a 15 M☉ 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}$ M☉ of CO has formed by 330 d, while our BSG model (at 330 d) contains $6.1 \times 10^{-5}$ M☉ 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 ~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 ~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 (≥ 300 d post-explosion) times and low (~1000K) 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 $C + O \rightarrow CO + \gamma$ typically dominates the formation of CO, though the charge exchange reaction $CO^+ + O \rightarrow CO + O^+$ can become important when the density of neutral oxygen is high. Charge exchange with oxygen: $CO + O^+ \rightarrow CO^+ + O$ typically dominates the destruction of CO, but dissociation by non-thermal electrons ($CO + e_{non-therm}^- \rightarrow C/C^+ + O/O^+$) can also contribute significantly, especially at lower densities. CO ionization balance is dictated by charge exchange with O and $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 (> 25 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⁻¹ with the remaining CO spread over the broader oxygen/carbon-rich region extending from 2500 to 4500 km s⁻¹. 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⁻¹.

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⁻¹, where the oxygen ionization fraction is lower than in the hotter region below 1500 km s⁻¹. Ionization of

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Colourmap of CO formation in our RSG model. Each pixel is coloured according to the total mass of CO present within a range of 100 km s⁻¹. Most of the CO is concentrated in the dense O/Si-rich region at ~1500 km s⁻¹. A secondary region of CO formation is apparent between ~3000 and 5000 km s⁻¹, which is more extended and has lower density. The two regions are separated by a ~1000 km s⁻¹ wide region which is almost devoid of CO. The mass of CO exterior to 6000 km s⁻¹ is negligible.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Total mass of CO present in the RSG model (black) as a function of time since explosion. Also shown are ionization fractions for carbon (orange) and oxygen (blue) in the CO region. We calculate the ionization fractions as an average weighted by the final CO density: $\frac{\int(N^+/N_{tot})\times\rho_{CO}dV}{\int\rho_{CO}dV}$ (where $N_{tot}$ is the total number density of oxygen [carbon] and $N^+$ the number density of singly ionized oxygen [carbon]) to illustrate the importance of the ionization structure on CO formation. The amount of CO present is negligible before 60 d. Most of the CO is formed in a short window between ~110 and ~130 d post-explosion, as oxygen and carbon start to recombine. The amount of CO reaches a maximum at ~150 d, before roughly plateauing as carbon remains ionized for some time. After 220 d, the mass of CO begins to rise again as neutral carbon becomes the dominant ionization stage.</figcaption></figure>

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}M_{\odot}$.
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Temperature structures in our RSG models at 200 d post-explosion. The two separate regions of CO formation are visible here as lower temperature regions: a narrow one near 1400 km s⁻¹, and a broader region between 2500 and 4000 km s⁻¹. The presence of CO in these regions has accelerated the cooling of the ejecta significantly by 200 d.</figcaption></figure>

## 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α 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 ~4.5–5.5 µ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 ~2 × 10⁻⁴ M☉ 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,

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Optical spectra computed from BSG models at 330 d without CO (NOCO), with CO-level populations in LTE (LTECO), and with CO-level populations allowed to depart from LTE (NLTECO). The strongest emission lines are identified. The CO models are cooler in the inner region (< 3500 km s⁻¹), reducing the emissivity of lines formed there, with the most significant effect on Ca II lines in the optical. Also shown are transmission curves for Johnson–Cousins optical filters (Bessell 2005).</figcaption></figure>
<figure><img src="image_8.png" /></figure>

**Figure 8.** Infrared spectra for the same BSG models as Fig. 7. Emission from CO primarily emerges in the first overtone (~2.5 µm, Δν = 2) and fundamental bands (~5 µm, Δν = 1). Emission from CO is reduced when level populations are allowed to depart from LTE, reducing the net cooling effect from CO. Also shown are some infrared filter transmission functions: J, H, and K bands used by the European Southern Observatory's Very Large Telescope (Wamsteker 1981), and the JWST's NIRCam F277W and F444W bands (Rieke, Kelly & Horner 2005), which cover the main regions of CO emission. The plotted curves for F277W and F444W include all mirror and detector effects, while the J-, H-, and K-band curves account for only the filter transmission.

<figure><img src="image_9.png" /></figure>

**Figure 9.** Light curves derived from our RSG model with and without CO formation included. Effective wavelengths for each filter are indicated in the legend. Magnitudes are given as AB magnitudes observed at a nominal distance of 10 Mpc. In optical and NIR filters, the CO model decays more rapidly at late times as CO cooling begins to dominate. In Mid-IR filters which overlap the CO emission bands (like F444W), the CO model is much brighter, with its brightness still increasing at 200 d.

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

<figure><img src="image_10.png" /></figure>

**Figure 10.** NIR RSG model spectrum at 210 d with a NIR spectrum of Type IIP SN2017eaw taken at 215 d after explosion. The SN2017eaw spectrum is taken from Rho et al. (2018). The inset plot shows the CO first overtone emission band. The model overestimates the CO flux, especially in the higher excitation bandheads, despite containing less CO mass than was estimated for 2017eaw. Some discrepancy is expected from our assumption of LTE-level populations in CO. The model also significantly overestimates the strength of several Mg I lines in this region, suggesting that the model's magnesium abundance may be much higher than SN2017eaw. The strength of the Mg I lines is slightly increased when CO is included, as a result of lower ionization in the CO-dominated regions.

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
<figure><img src="image_11.png" /><figcaption>Figure 11. Heating and cooling balance at a single region in our BSG model. The equilibrium temperature is where total cooling (green) equals total heating (red, from radioactive decay). The addition of CO cooling causes the cooling curve to rise at temperatures near 2000 K. The CO cooling drops off again below 2000 K, with the result that the total cooling curve intersects the heating curve at up to three distinct temperatures, labelled T<sub>A</sub>, T<sub>B</sub>, and T<sub>C</sub>.</figcaption></figure>

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 ~30 when the temperature decreases from 4000 to 2000 K. However, CO emissivity per molecule is an increasing function of temperature, roughly proportional to e<sup>-hv/k<sub>B</sub>T</sup>. 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:

$$ \int_{\text{CO lines}} N_u h \nu Z_{\text{line}} = \int_{\text{CO lines}} N_u h \nu \left( 1 - \frac{\bar{J}_{\text{line}}}{S_{\text{line}}} \right), \quad (7) $$

where N<sub>u</sub> is the upper state population, Z<sub>line</sub> the net radiative bracket, J<sub>line</sub> the frequency-averaged mean intensity, and S<sub>line</sub> 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 T<sub>A</sub>, T<sub>B</sub>, and T<sub>C</sub> 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 (~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 T<sub>C</sub> 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 ~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 T<sub>C</sub>, 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 T<sub>A</sub> as a temperature solution regardless of the initial conditions. Conversely, a 50 per cent decrease in CO abundance moves it to scenario 3, leaving T<sub>C</sub> as the only solution. An equivalent change in the radioactive decay heating rate can also change the solution scenario. Thus a factor of ~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 ~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 ~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
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Comparison of CO line appearance with and without the assumption of LTE. The electron densities in our BSG model are below the critical density for the CO *v*: 0 → 1 transition, resulting in excited-state populations well below the corresponding LTE populations. This causes a factor ~3 reduction in flux in the fundamental and first overtone bands. In both cases, the first overtone lines remain optically thin (τ ~ 0.5), while the fundamental band is optically thick (τ ~ 5). The 'simple model' shown here assumes a spherical emission region of LTE, optically thin CO gas with total mass, max velocity, and excitation temperature chosen to match the CMFGEN model properties, as used in Spyromilio et al. (1988). It closely reproduces the first overtone band (in LTE), but overestimates even the LTE fundamental band emission by a factor of ~4. Also shown is the spectrum of SN1987A at 260 d, scaled to match the change in bolometric magnitude between 260 and 330 d (Arnett et al. 1989; Meikle et al. 1989). Our LTE model matches the flux in the first bandhead, but not the overall structure in the first overtone, and overestimates the fundamental band flux. Our NLTE model also overestimates the fundamental band flux despite underestimating the flux in the first overtone.</figcaption></figure>

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 ≥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 $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 $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 ~ 1 × 10⁻⁴ M☉ at 330 d, observational estimates vary (Section 6.1), mostly ≳ 10⁻⁴ M☉

(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$)<fnref n="8" /> and the fundamental/first overtone flux ratio is < 20. Observations of SN1987A at 260 d post-explosion from Meikle et al. (1989) show a fundamental/first overtone flux ratio of ~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 ~ $10^{-4} 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} 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} 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} 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-excitation 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} 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 (~50) of fundamental band to first overtone band flux, which is in tension with the observed ratio (~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} M_\odot$ of CO in our BSG model (total ejecta mass 13.2 $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 $O^+$, which we find contributes significantly even after 200 d) and we find a similar impact from CO cooling, which reduces the

<footnote n="8">⁸We calculate the optical depth in CO bands using direct integration in the observer's frame, accounting for the overlap of multiple rovibrational transitions.</footnote>
<table><caption>**Table 2.** CO masses computed at ~300 d post-explosion for models of SN1987A and SN1987A-like models.</caption><thead><tr><th rowspan="2">Study</th><th colspan="2">CO mass in SN chemical evolution models</th></tr><tr><th>CO mass (M<sub>☉</sub>)</th><th>Epoch (days after explosion)</th></tr></thead><tbody><tr><td>Petuchowski et al. (1989)</td><td>2.5 × 10<sup>-5</sup></td><td>300 d</td></tr><tr><td>Lepp et al. (1990)</td><td>4.1 × 10<sup>-7</sup>–2.8 × 10<sup>-5</sup></td><td>349 d</td></tr><tr><td>Liu & Dalgarno (1995)</td><td>3.0 × 10<sup>-3</sup></td><td>300 d</td></tr><tr><td>Cherchneff & Sarangi (2011)</td><td>2 – 4 × 10<sup>-3</sup></td><td>300 d</td></tr><tr><td>Sluder, Milosavljević & Montgomery (2018)</td><td>~ 10<sup>-2</sup></td><td>300 d</td></tr><tr><td>Liljegren et al. (2020)</td><td>3 × 10<sup>-3</sup></td><td>300 d</td></tr><tr><td>Ono et al. (2023)</td><td>3 × 10<sup>-2</sup>–2 × 10<sup>-1</sup></td><td>300 d</td></tr><tr><td>This work</td><td>~ 10<sup>-4</sup></td><td>330 d</td></tr></tbody></table>

*Notes.* Differences among the models can be attributed to different assumptions about the ejecta composition, reaction networks, and reaction data.

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⁻¹), in agreement with observations of CO. In our RSG model, CO becomes important (i.e. contributes a significant fraction of the total cooling) at ~100 d after the explosion, and by ~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 ~4.5 µm. In this band, the ejecta remain optically thick well past 300 d. The CO first overtone band at ~2.5 µ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.

## DATA AVAILABILITY

CMFGEN is available on Github or through the CMFGEN website at http://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.
REFERENCES

Abellán F. J. et al., 2017, ApJ, 842, L24

Arnett W. D., Bahcall J. N., Kirshner R. P., Woosley S. E., 1989, ARA&A,
27, 629

Banerjee D. P. K. et al., 2016, The Astronomer's Telegram, 8976, 1

Bessell M. S., 2005, ARA&A, 43, 293

Bestenlehner J. M. et al., 2014, A&A, 570, A38

Blondin S., Dessart L., Hillier D. J., Khokhlov A. M., 2013, MNRAS, 429,
2127

Bouchet P., Danziger I. J., 1993, A&A, 273, 451

Chandra S., Sharma A. K., 2001, A&A, 376, 356

Cherchneff I., Sarangi A., 2011, Proc. IAU Symp. 7, The Molecular Universe,
Kluwer, Dordrecht. p. 228

Dalgarno A., Du M. L., You J. H., 1990, ApJ, 349, 675

De Marco O., Schmutz W., Crowther P. A., Hillier D. J., Dessart L., de Koter
A., Schweickhardt J., 2000, A&A, 358, 187

Dessart L., Hillier D. J., 2005, A&A, 437, 667

Dessart L., Hillier D. J., 2010, MNRAS, 405, 2141

Dessart L., Hillier D. J., 2019, A&A, 622, A70

Dessart L., Hillier D. J., 2020a, A&A, 642, A33

Dessart L., Hillier D. J., 2020b, A&A, 643, L13

Dessart L., Hillier D. J., 2022, A&A, 660, L9

Dessart L., Hillier D. J., Blondin S., Khokhlov A., 2014, MNRAS, 441, 3249

Dessart L., Hillier D. J., Wilk K. D., 2018, A&A, 619, A30

Dessart L., Yoon S.-C., Aguilera-Dena D. R., Langer N., 2020, A&A, 642, A106

Dessart L., Hillier D. J., Sukhbold T., Woosley S. E., Janka H.-T., 2021, A&A,
652, A64

Draine B. T., Spergel David, 2011, Physics of the Interstellar and Intergalactic Medium, Princeton University Press, Princeton, New Jersey, 440, https://ui.adsabs.harvard.edu/abs/2011piim.book.....D

Drout M. R. et al., 2016, ApJ, 821, 57

Ergon M. et al., 2014, A&A, 562, A17

Ergon M., Fransson C., Jerkstrand A., Kozma C., Kromer M., Spricer K.,
2018, A&A, 620, A156

Essenhigh K. A., 2005, PhD thesis, The Ohio State University, https://etd.ohiolink.edu/acprod/odb_etd/etd/r/1501/10?clear=10&p10_accession_num = osu1124292824

Fassia A. et al., 2001, MNRAS, 325, 907

Fransson C., Chevalier R. A., 1989, ApJ, 343, 323

Fryxell B., Mueller E., Arnett D., 1991, ApJ, 367, 619

Gandhi P. et al., 2013, ApJ, 767, 166

Gearhart R. A., Wheeler J. C., Swartz D. A., 1999, ApJ, 510, 944

Gehrz R., 1989, Proc. IAU Symp. Vol. 135 Interstellar Dust. Kluwer,
Dordrecht. p.445

Gerardy C. L., Fesen R. A., Höflich P., Wheeler J. C., 2000, AJ, 119, 2968

Herzberg G., Mrozowski S., 1951, Am. J. Phys., 19, 390

Hillier D. J., 1990, A&A, 231, 116

Hillier D. J., Dessart L., 2012, MNRAS, 424, 252

Hillier D. J., Dessart L., 2019, A&A, 631, A8

Hillier D. J., Lanz T., 2001, ASP Conf. Ser. Vol. 247, Spectroscopic Challenges of Photoionized Plasmas, Astron. Soc. Pac., San Francisco. p. 343

Hillier D. J., Miller D. L., 1998, ApJ, 496, 407

Itikawa Y., 2015, J. Phys. Chem. Ref. Data, 44, 013105

Jacobson-Galán W. V. et al., 2022, ApJ, 924, 15

Jerkstrand A., 2011, PhD thesis, Accessed 6/2/2023, https://ui.adsabs.harvard.d.edu/abs/2011PhDT.......90J

Jerkstrand A., Alsabti A., Murdin P., 2017, Handbook of Supernovae, Springer International Publishing AG, New York, p. 842

Keszthelyi Z., Meynet G., Martins F., de Koter A., David-Uraz A., 2021, MNRAS, 504, 2474

Kifonidis K., Plewa T., Janka H. T., Müller E., 2003, A&A, 408, 621

Kleiser I. et al., 2011, Bulletin of the American Astronomical Society, AAS Meeting #217, 43, 337

Kotak R., Meikle P., van Dyk S. D., Höflich P. A., Mattila S., 2005, ApJ, 628, L123

Kotak R. et al., 2009, ApJ, 704, 306

Leonard D. C. et al., 2002, Publ. Astron. Soc. Pac., 114, 35

Lepp S., Dalgarno A., McCray R., 1990, ApJ, 358, 262

Li C., Hillier D. J., Dessart L., 2012, MNRAS, 426, 1671

Li G., Gordon I. E., Rothman L. S., Tan Y., Hu S.-M., Kassi S., Campargue A., Medvedev E. S., 2015, ApJS, 216, 15

Liljegren S., Jerkstrand A., Grumer J., 2020, A&A, 642, A135

Liljegren S., Jerkstrand A., Barklem P. S., Nyman G., Brady R., Yurchenko S. N., 2022, A&A, 674,32

Liu W., Dalgarno A., 1994, ApJ, 428, 769

Liu W., Dalgarno A., 1995, ApJ, 454, 472

Liu W., Dalgarno A., Lepp S., 1992, ApJ, 396, 679

Livne E., 1993, ApJ, 412, 634

Maguire K. et al., 2012, MNRAS, 420, 3451

Martins F., Schaerer D., Hillier D. J., 2005, A&A, 436, 1049

Massey P., Neugent K. F., Morrell N., Hillier D. J., 2014, ApJ, 788, 83

Matsuura M., 2017, in Alsabti A, Murdin P, eds, Dust and Molecular Formation in Supernovae. Springer International Publishing AG, New York, p.2158, https://ui.adsabs.harvard.edu/abs/2017hsn..book.2125M

Mazzali P. A., Lucy L. B., 1993, A&A, 279, 447

McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36

Meikle W. P. S., Allen D. A., Spyromilio J., Varani G. F., 1989, MNRAS,
238, 193

Meikle W. P. S., Spyromilio J., Allen D. A., Varani G. F., Cumming R. J., 1993, MNRAS, 261, 535

Meng H., Qin Z., Liu L., 2022, ApJ, 935, 148

Mueller E., Fryxell B., Arnett D., 1991, A&A, 251, 505

Ogloblina P., Tejero-del Caz A., Guerra V., Alves L. L., 2020, Plasma Sources Sci. Technol., 29, 015002

Ono M., Nozawa T., Nagataki S., Kozyreva A., Orlando S., Miceli M., Chen K.-J., 2023, ApJS, 271, 33

Pancheshnyi S., Biagi S., Bordage M. C., Hagelaar G. J. M., Morgan W. L., Phelps A. V., Pitchford L. C., 2012, Chem. Phys., 398, 148

Pastorello A. et al., 2012, A&A, 537, A141

Paxton B. et al., 2018, ApJS, 234, 34

Petuchowski S. J., Dwek E., Allen J. E., Jr, Nuth J. A., III, 1989, ApJ, 342, 406

Rho J., Geballe T. R., Banerjee D. P. K., Dessart L., Evans A., Joshi V., 2018, ApJ, 864, L20

Rho J. et al., 2021, ApJ, 908, 232

Rieke M. J., Kelly D., Horner S., 2005, Proc. SPIE, 5904, 1

Sarangi A., Matsuura M., Micelotta E. R., 2018, Space Sci. Rev., 214, 63

Sluder A., Milosavljević M., Montgomery M. H., 2018, MNRAS, 480, 5580

Spyromilio J., Meikle W. P. S., Learner R. C. M., Allen D. A., 1988, Nature, 334, 327

Stritzinger M. D. et al., 2023, A&A, 686, 30

Swartz D. A., Harkness R. P., Wheeler J. C., 1989, Nature, 337, 439

Thompson H. M. A., Keenan F. P., Dufton P. L., Trundle C., Ryans R. S. I., Crowther P. A., 2008, MNRAS, 383, 729

Wamsteker W., 1981, A&A, 97, 329

Watson D., Christensen L., Knudsen K. K., Richard J., Gallazzi A., Michałowski M. J., 2015, Nature, 519, 327

Wongwathanarat A., Müller E., Janka H. T., 2015, A&A, 577, A48

Woosley S. E., 1988, ApJ, 330, 218

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Yuan F. et al., 2016, MNRAS, 461, 2003

APPENDIX A: REACTION DATA
<table><caption>**Table A1.** Basic reaction data for all molecular reactions included in CMFGEN.</caption><thead><tr><th>ID</th><th>Reaction</th><th>α</th><th>β</th><th>γ</th><th>T<sub>low</sub></th><th>T<sub>high</sub></th></tr></thead><tbody><tr><td>1</td><td>C<sub>2</sub> + O<sub>2</sub><sup>+</sup> → C<sub>2</sub><sup>+</sup> + O<sub>2</sub></td><td>4.10E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>2</td><td>C<sub>2</sub> + O<sup>+</sup> → C<sub>2</sub><sup>+</sup> + O</td><td>4.80E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>3</td><td>C + C<sub>2</sub><sup>+</sup> → C<sup>+</sup> + C<sub>2</sub></td><td>1.10E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>4</td><td>C + CO<sup>+</sup> → C<sup>+</sup> + CO</td><td>1.10E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>5</td><td>C + O<sub>2</sub><sup>+</sup> → C<sup>+</sup> + O<sub>2</sub></td><td>5.20E-11</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>6</td><td>C<sup>+</sup> + C<sub>2</sub>O → C + C<sub>2</sub>O<sup>+</sup></td><td>1.00E-9</td><td>-0.5</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>7</td><td>CO<sub>2</sub><sup>+</sup> + O<sub>2</sub> → CO<sub>2</sub> + O<sub>2</sub><sup>+</sup></td><td>5.30E-11</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>8</td><td>CO<sub>2</sub><sup>+</sup> + O → CO<sub>2</sub> + O<sup>+</sup></td><td>9.62E-11</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>9</td><td>CO + O<sup>+</sup> → CO<sup>+</sup> + O</td><td>4.90E-12</td><td>0.5</td><td>4.58E3</td><td>2000</td><td>10 000</td></tr><tr><td>10</td><td>CO<sup>+</sup> + C<sub>2</sub> → CO + C<sub>2</sub><sup>+</sup></td><td>8.40E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>11</td><td>CO<sup>+</sup> + CO<sub>2</sub> → CO + CO<sub>2</sub><sup>+</sup></td><td>1.00E-9</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>12</td><td>CO<sup>+</sup> + O<sub>2</sub> → CO + O<sub>2</sub><sup>+</sup></td><td>1.20E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>13</td><td>CO<sup>+</sup> + O → CO + O<sup>+</sup></td><td>1.40E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>14</td><td>O<sup>+</sup> + O<sub>2</sub> → O + O<sub>2</sub><sup>+</sup></td><td>1.90E-11</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>15</td><td>C<sub>2</sub> + O<sub>2</sub><sup>+</sup> → CO + CO<sup>+</sup></td><td>4.10E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>16</td><td>C<sub>2</sub> + O<sup>+</sup> → C + CO<sup>+</sup></td><td>4.80E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>17</td><td>C<sub>2</sub><sup>+</sup> + O<sub>2</sub> → CO + CO<sup>+</sup></td><td>8.00E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>18</td><td>C<sub>2</sub><sup>+</sup> + O → C + CO<sup>+</sup></td><td>3.10E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>19</td><td>C + O<sub>2</sub><sup>+</sup> → CO<sup>+</sup> + O</td><td>5.20E-11</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>20</td><td>C<sup>+</sup> + CO<sub>2</sub> → CO + CO<sup>+</sup></td><td>1.10E-9</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>21</td><td>C<sup>+</sup> + O<sub>2</sub> → CO + O<sup>+</sup></td><td>4.54E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>22</td><td>C<sup>+</sup> + O<sub>2</sub> → CO<sup>+</sup> + O</td><td>3.42E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>23</td><td>CO<sub>2</sub> + O<sup>+</sup> → CO + O<sub>2</sub><sup>+</sup></td><td>9.40E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>24</td><td>CO<sub>2</sub><sup>+</sup> + O → CO + O<sub>2</sub><sup>+</sup></td><td>1.64E-10</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>25</td><td>C + CO<sub>2</sub> → CO + CO</td><td>1.00E-15</td><td>0</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>26</td><td>C<sub>2</sub> + O<sub>2</sub> → CO + CO</td><td>1.50E-11</td><td>0</td><td>4.3E3</td><td>298</td><td>1300</td></tr><tr><td>27</td><td>C<sub>2</sub> + O → C + CO</td><td>2.00E-10</td><td>-1.2E-1</td><td>0</td><td>10</td><td>8000</td></tr><tr><td>28</td><td>C<sub>2</sub>O + O → CO + CO</td><td>8.59E-11</td><td>0</td><td>0</td><td>295</td><td>296</td></tr><tr><td>29</td><td>C + C<sub>2</sub>O → CO + C<sub>2</sub></td><td>2.00E-10</td><td>0</td><td>0</td><td>10</td><td>300</td></tr><tr><td>30</td><td>C + CO → C<sub>2</sub> + O</td><td>2.94E-11</td><td>5.0E-1</td><td>5.80250E4</td><td>1934</td><td>41 000</td></tr><tr><td>31</td><td>C + O<sub>2</sub> → CO + O</td><td>5.56E-11</td><td>4.10E-1</td><td>-2.69E1</td><td>10</td><td>8000</td></tr><tr><td>32</td><td>CO<sub>2</sub> + O → CO + O<sub>2</sub></td><td>2.46E-11</td><td>0</td><td>2.6567E4</td><td>300</td><td>6000</td></tr><tr><td>33</td><td>CO + O<sub>2</sub> → CO<sub>2</sub> + O</td><td>5.99E-12</td><td>0</td><td>2.4075E4</td><td>300</td><td>6000</td></tr><tr><td>34</td><td>CO + O → C + O<sub>2</sub></td><td>1.00E-16</td><td>0</td><td>0</td><td>1.0</td><td>1.1</td></tr><tr><td>35</td><td>C + C → C<sub>2</sub> + photon</td><td>4.36E-18</td><td>3.5E-1</td><td>1.613E2</td><td>10</td><td>41 000</td></tr><tr><td>36</td><td>C + C<sup>+</sup> → C<sub>2</sub><sup>+</sup> + photon</td><td>4.01E-18</td><td>1.7E-1</td><td>1.015E2</td><td>10</td><td>41 000</td></tr><tr><td>37</td><td>C + O → CO + photon</td><td>4.69E-19</td><td>1.52</td><td>-5.05E1</td><td>10</td><td>14700</td></tr><tr><td>38</td><td>C + O<sup>+</sup> → CO<sup>+</sup> + photon</td><td>5.00E-10</td><td>-3.7</td><td>8.0E2</td><td>2000</td><td>10 000</td></tr><tr><td>39</td><td>C<sup>+</sup> + O → CO<sup>+</sup> + photon</td><td>3.14E-18</td><td>-1.5E-1</td><td>6.80E1</td><td>10</td><td>13 900</td></tr><tr><td>40</td><td>O + O → O<sub>2</sub> + photon</td><td>4.90E-20</td><td>1.58</td><td>0</td><td>10</td><td>300</td></tr><tr><td>41</td><td>CO<sup>+</sup> + e<sup>-</sup> → C + O</td><td>2.00E-7</td><td>-4.8E-1</td><td>0</td><td>10</td><td>41 000</td></tr><tr><td>42</td><td>CO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → CO + O</td><td>3.80E-7</td><td>-5.0E-1</td><td>0</td><td>10</td><td>300</td></tr><tr><td>43</td><td>C<sub>2</sub><sup>+</sup> + e<sup>-</sup> → C + C</td><td>3.00E-7</td><td>-5.0E-1</td><td>0</td><td>10</td><td>300</td></tr><tr><td>44</td><td>C<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → CO + C</td><td>3.00E-7</td><td>-5.0E-1</td><td>0</td><td>10</td><td>300</td></tr><tr><td>45</td><td>O<sub>2</sub><sup>+</sup> + e<sup>-</sup> → O + O</td><td>1.95E-7</td><td>-7.0E-1</td><td>0</td><td>10</td><td>300</td></tr></tbody></table>

Notes. α has units of cm⁻³ s⁻¹, β is dimensionless, and T<sub>high</sub>, T<sub>low</sub>, and γ are in K. Rates are calculated using the specified constants for each reaction following equation (3). All data were obtained in 2021 from the UMIST Database for Astrochemistry (McElroy et al. 2013).

This paper has been typeset from a TeX/LATEX file prepared by the author.