# A SUBLIME 3D Model for Cometary Coma Emission: The Hypervolatile-rich Comet C/2016 R2 (PanSTARRS)

M. A. Cordiner, I. M. Coulson, E. Garcia-Berrios, C. Qi, François Lique, Michal Zoltowski, M. de Val-Borro, Y-J Kuan, W-H Ip, S. Mairs, et al.

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

## ► To cite this version:

M. A. Cordiner, I. M. Coulson, E. Garcia-Berrios, C. Qi, François Lique, et al., A SUBLIME 3D Model for Cometary Coma Emission: The Hypervolatile-rich Comet C/2016 R2 (PanSTARRS). The Astrophysical Journal, 2022, 929 (1), pp.38. 10.3847/1538-4357/ac5893. [hal-03660671v2](https://hal.science/hal-03660671v2)

HAL Id: hal-03660671

https://hal.science/hal-03660671v2

Submitted on 6 May 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# A SUBLIME 3D Model for Cometary Coma Emission: The Hypervolatile-rich Comet C/2016 R2 (PanSTARRS)

M. A. Cordiner<sup>1,2 ID</sup>, I. M. Coulson<sup>3 ID</sup>, E. Garcia-Berrios<sup>1,2,4 ID</sup>, C. Qi<sup>5 ID</sup>, F. Lique<sup>6 ID</sup>, M. Zołtowski<sup>6,7</sup>, M. de Val-Borro<sup>1,2 ID</sup>, Y.-J. Kuan<sup>8,9 ID</sup>, W.-H. Ip<sup>10</sup>, S. Mairs<sup>3 ID</sup>, N. X. Roth<sup>1,2 ID</sup>, S. B. Charnley<sup>1 ID</sup>, S. N. Milam<sup>1 ID</sup>, W.-L Tseng<sup>8 ID</sup>, and Y.-L Chuang<sup>8 ID</sup>

<sup>1</sup> Astrochemistry Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA; martin.cordiner@nasa.gov

<sup>2</sup> Department of Physics, Catholic University of America, Washington, DC 20064, USA

<sup>3</sup> East Asian Observatory, 660 N. A'ohoku Place, Hilo, HI 96720, USA

<sup>4</sup> Department of Astronomy, University of Illinois, 1002 W. Green St., Urbana, IL 61801, USA

<sup>5</sup> Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, MS 42, Cambridge, MA 02138, USA

<sup>6</sup> Université de Rennes 1, Campus de Beaulieu, 263 avenue du Général Leclerc, F-35042 Rennes Cedex, France

<sup>7</sup> LOMC—UMR 6294, CNRS-Université du Havre, 25 rue Philippe Lebon, BP1123, F-76 063 Le Havre cedex, France

<sup>8</sup> National Taiwan Normal University, Ting-Chou Road, Taipei 11677, Taiwan, ROC

<sup>9</sup> Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 106, Taiwan, ROC

<sup>10</sup> Graduate Institute of Space Science, National Central University, Taoyuan City, Taiwan

Received 2021 November 3; revised 2022 February 23; accepted 2022 February 23; published 2022 April 11

## Abstract

The coma of comet C/2016 R2 (PanSTARRS) is one of the most chemically peculiar ever observed, in particular due to its extremely high CO/H₂O and N₂⁺/H₂O ratios, and unusual trace volatile abundances. However, the complex shape of its CO emission lines, as well as uncertainties in the coma structure and excitation, has lead to ambiguities in the total CO production rate. We performed high-resolution, spatially, spectrally, and temporally resolved CO observations using the James Clerk Maxwell Telescope and Submillimeter Array to elucidate the outgassing behavior of C/2016 R2. Results are analyzed using a new, time-dependent, three-dimensional radiative transfer code (SUBlimating gases in LIME; SUBLIME, based on the open-source version of the LIne Modeling Engine), incorporating for the first time, accurate state-to-state collisional rate coefficients for the CO–CO system. The total CO production rate was found to be in the range of (3.8–7.6) × 10²⁸ s⁻¹ between 2018 January 13 and February 1 (at r<sub>H</sub> = 2.8–2.9 au), with a mean value of (5.3 ± 0.6) × 10²⁸ s⁻¹. The emission is concentrated in a near-sunward jet, with a half-opening angle of ~62° and an outflow velocity of 0.51 ± 0.01 km s⁻¹, compared to 0.25 ± 0.01 km s⁻¹ in the ambient (and nightside) coma. Evidence was also found for an extended source of CO emission, possibly due to icy grain sublimation around 1.2 × 10⁵ km from the nucleus. Based on the coma molecular abundances, we propose that the nucleus ices of C/2016 R2 can be divided into a rapidly sublimating apolar phase, rich in CO, CO₂, N₂, and CH₃OH, and a predominantly frozen (or less abundant), polar phase containing more H₂O, CH₄, H₂CO, and HCN.

*Unified Astronomy Thesaurus concepts: Comets (280); Long period comets (933); Submillimeter astronomy (1647); Comet volatiles (2162); High resolution spectroscopy (2096); Radio interferometry (1346); Radiative transfer simulations (1967); De-excitation rates (2066)*

## 1. Introduction

Comets are composed of ice, dust, and debris accreted during the epoch of planet formation. Having spent most of their lives frozen, at large distances from the Sun, cometary nuclei contain some of our solar system's most pristine (thermally unprocessed) material. From studies of their gaseous atmospheres (comae), the properties of the nucleus can be inferred, thereby providing unique insights into the conditions prevalent at the dawn of the solar system.

The long-period comet C/2016 R2 (PanSTARRS) is one of the most chemically peculiar comets ever observed. Early observations at a heliocentric distance of around 3 au revealed a visually blue coma with a highly unusual optical spectrum, dominated by CO⁺ emission, in addition to a rare detection of N₂⁺ (Cochran & McKay 2018), consistent with very strong outgassing of CO and an above average N₂/CO ratio in the

nucleus. Both CO and N₂ sublimate at low temperatures compared to other cometary ices (e.g., Womack et al. 2017), and are therefore considered hypervolatile. Their presence in very high abundances (relative to H₂O) implies that C/2016 R2 formed very cold, and was maintained at temperatures ≲20 K for the duration of its lifetime.

Follow-up observations at infrared and millimeter wavelengths confirmed the presence of an exceptionally CO-rich coma, with a low production rate of H₂O (usually the dominant cometary volatile) and extremely unusual abundance ratios for other species compared with the typical cometary population (Biver et al. 2018; Wierzchos & Womack 2018; McKay et al. 2019). Due to their very different sublimation temperatures, the CO-to-H₂O ratio in cometary comae is known to vary as a function of temperature of the nucleus, so some enhancement in CO/H₂O is expected at the relatively large heliocentric distance (r<sub>H</sub>) at which this comet was observed. However, the coma composition of C/2016 R2 clearly differs from other comets observed at similar r<sub>H</sub>. Based on Table 7 of Biver et al. (2018), the abundance ratios CO:H₂O:CH₃OH:HCN (normalized to CH₃OH) were 94:0.3:1:0.004 in C/2016 R2 at

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
$r_H = 2.8$ au, 14:43:1:0.08 in C/1995 O1 at 2.8 au, and 24:27:1:0.10 in C/2006 W3 at 3.2–3.3 au (the H₂O abundances in R2 and W3 are from McKay et al. 2019 and Bockelée-Morvan et al. 2010, respectively), highlighting the CO richness of the coma, as well as its strong H₂O and HCN depletions. Cometary and interstellar ice observations have identified the presence of separate polar and apolar ice phases, dominated by H₂O and [CO + CO₂], respectively (Mumma et al. 2011; Boogert et al. 2015; Luspay-Kuti et al. 2015). Studies of coma abundances in C/2016 R2 therefore provide a rare opportunity to investigate volatiles outgassed primarily from the apolar (CO + CO₂ dominated) phase, which may provide unique insights into the origin, storage, and outgassing mechanisms of the less abundant ices in cometary nuclei.

Complexity of the C/2016 R2 blueshifted CO rotational line profile, as well as uncertainties in the CO excitation calculation, lead to ambiguity in its interpretation, and estimates for the CO production rate (in 2018 January) range from $Q(\text{CO}) = (4.6 \pm 0.4) \times 10^{26} \text{ s}^{-1}$ (Wierzchos & Womack 2018) to $(10.6 \pm 0.5) \times 10^{26} \text{ s}^{-1}$ (Biver et al. 2018). Limited CO mapping by both studies found marginal evidence for extended CO emission beyond that expected based purely on nucleus-driven outgassing, but this possibility has not yet been investigated in detail. Uncertainties therefore remain regarding (1) the intrinsic CO production rate of this comet, (2) the processes by which CO is released into the gas phase, and (3) the detailed coma morphology.

During 2018 January and February, we undertook a program of time-resolved high-resolution spectroscopy, spatial-spectral mapping, and interferometry using the James Clerk Maxwell Telescope (JCMT) and Submillimeter Array (SMA) to elucidate the outgassing behavior of C/2016 R2. The CO $J=3-2$ and $J=2-1$ lines were observed as a probe of the coma kinetic temperature, and based on the strength of the ¹²CO emission in this comet, we searched for the ¹³CO isotopologue as a tracer of any unusual isotopic processing in this comet's natal carbon. We also sought to confirm the comet's peculiar CO:HCN:CH₃OH:H₂CO ratios through observations of submillimeter rotational lines from these species, and HCO⁺ was observed as a probe of ion chemistry in the outer coma.

The resulting (spectral-spatial-temporal) data set is analyzed using a new, time-dependent, three-dimensional radiative transfer code (SUBlimating gases in LIME; SUBLIME), which is an evolution of the steady-state model used previously by Paganini et al. (2010), Bøgelund & Hogerheijde (2017), de Val-Borro et al. (2018), and Roth et al. (2021a), and includes excitation via collisions with CO and electrons, as well as radiative processes. In Section 5.4, we briefly discuss the conditions for which a time-dependent solution of the molecular excitation produces more accurate results than the steady-state treatment.

The emission line profiles from C/2016 R2 are simulated using a two-component (conical jet + ambient coma) outgassing model, allowing the molecular production rates and outflow velocities to vary as a function of coma position. Due to a lack of known collision cross-sections, previous models for cometary CO emission have relied on approximations for the CO collisional excitation rates. We therefore performed quantum scattering calculations (using the coupled-states approximation; Kłos & Lique 2018) to determine, for the first time, accurate state-to-state collisional (de-)excitation rate coefficients for the CO-CO system. The parameters retrieved

from our SUBLIME modeling provide new insights into the coma properties and intrinsic nature of C/2016 R2 (Pan-STARRS), one of the most unusual comets of our time.

## 2. Observations

### 2.1. James Clerk Maxwell Telescope

Observations of comet C/2016 R2 (PanSTARRS) were conducted during the period 2018 January 13 to 2018 February 1 using the 15 m JCMT atop Maunakea. The comet's position on the sky was tracked using Jet Propulsion Laboratory (JPL) Horizons ephemeris number 14. During this time, the heliocentric distance decreased from $r_H=2.9$ to 2.8 au, the geocentric distance increased from $\Delta = 2.1$ to 2.3 au, and the Sun-comet-observer (phase) angle increased from $\phi = 15^\circ$ to $19^\circ$.

Spectral line observations were carried out using the RxA3m millimeter-wave receiver (operating at 212 to 274 GHz) and the Heterodyne Array Receiver Program (HARP) 16-element submillimeter focal-plane receiver array (operating at 325 to 375 GHz; Buckle et al. 2009). HARP observations of CO $J=3-2$ were made in two modes: stare and jiggle. During stare observations the pointing receptor, H05, tracked the coordinates of the target while the other receptors (in the square array) each recorded a spectrum at an offset point in the coma. The 30" spacing between the receptors means that the target field is undersampled during stare observations, as a result of the 14" telescope beam FWHM at 346 GHz (20" at 230 GHz). Jiggle observations, on the other hand, are designed to obtain a spatially complete sampling of the target field. We used the HARP4 jiggle pattern to obtain spectra for each point on a 7" × 7" grid, covering a 2' × 2' map area. Position switching (offset by 300" in azimuth) was performed for the purpose of subtracting spectral contributions from the terrestrial atmosphere and telescope optics. The performance of the individual HARP receptors was monitored throughout the observations, and any suboptimal receptors were flagged and masked during data reduction.

Most of the spectral data were obtained at a resolution of 31 kHz, across a 250 MHz bandpass, using the Auto Correlation Spectral Imaging System (ACSIS) digital autocorrelation spectrometer. For the ¹³CO and C¹⁸O spectral windows (observed simultaneously), 61 kHz resolution was used, and for CH₃OH, a relatively coarse resolution of 448 kHz (over a 1000 MHz bandpass) was used to cover multiple transitions of CH₃OH in the $J=7-6$ band at 338 GHz.

Standard observatory amplitude and pointing calibrations were performed at regular (at most, hourly) intervals, demonstrating nominal performance. The all-sky pointing performance of JCMT is ~2" in each of the (azimuth, elevation) coordinates (Coulson et al. 2020), but telescope tracking accuracy over the course of an hour for a given source is usually better than 1". Comparison of the actual telescope tracking position (apparent R.A. and decl.) with the ephemeris showed tracking errors to be less than 0.1" in each coordinate on each night. The observed antenna temperatures ($T_A^*$) were corrected for sky opacity, forward scattering, and spillover. Beam efficiency ($\eta$) was corrected using the standard values of $\eta = 0.60$ for RxA3m and 0.64 for HARP. The resulting main beam brightness temperature scale ($T_{MB}$) is expected to be accurate to within $\pm 10\%$.
<table><caption>**Table 1**
Log of JCMT Spectral Line Observations</caption><thead><tr><th>Date</th><th>Species</th><th>Transition</th><th>Obs. Mode</th><th>τ<sub>0</sub></th><th>Int. (s)</th><th>∫T<sub>MB</sub>dv (mK km s<sup>-1</sup>)</th><th>v&#773; (km s<sup>-1</sup>)</th><th>r<sub>H</sub> (au)</th><th>Δ (au)</th><th>φ (°)</th></tr></thead><tbody><tr><td>2018-01-13.215</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.10</td><td>600</td><td>870(92)</td><td>-0.21(0.04)</td><td>2.88</td><td>2.13</td><td>14.8</td></tr><tr><td>2018-01-13.253</td><td>HCN</td><td>4-3</td><td>HARP stare</td><td>0.10</td><td>1200</td><td>&lt;78</td><td>...</td><td>2.88</td><td>2.13</td><td>14.9</td></tr><tr><td>2018-01-13.338</td><td>CH<sub>3</sub>OH</td><td>7-6</td><td>HARP stare</td><td>0.11</td><td>900</td><td>73(20)<sup>a</sup></td><td>-0.08(0.13)<sup>a</sup></td><td>2.88</td><td>2.13</td><td>14.9</td></tr><tr><td>2018-01-13.401</td><td>HCO<sup>+</sup></td><td>4-3</td><td>HARP stare</td><td>0.11</td><td>900</td><td>&lt;93</td><td>...</td><td>2.88</td><td>2.13</td><td>14.9</td></tr><tr><td>2018-01-13.437</td><td>H<sub>2</sub>CO</td><td>5<sub>1.5</sub>-4<sub>1.4</sub></td><td>HARP stare</td><td>0.11</td><td>900</td><td>&lt;110</td><td>...</td><td>2.88</td><td>2.13</td><td>14.9</td></tr><tr><td>2018-01-14.206</td><td>CO</td><td>2-1</td><td>RxA3m</td><td>0.17</td><td>600</td><td>542(71)</td><td>-0.18(0.08)</td><td>2.87</td><td>2.14</td><td>15.1</td></tr><tr><td>2018-01-14.248</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.17</td><td>120</td><td>1078(154)</td><td>-0.24(0.10)</td><td>2.87</td><td>2.14</td><td>15.1</td></tr><tr><td>2018-01-14.260</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.17</td><td>600</td><td>926(103)</td><td>-0.22(0.05)</td><td>2.87</td><td>2.14</td><td>15.1</td></tr><tr><td>2018-01-14.273</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.17</td><td>120</td><td>831(131)</td><td>-0.21(0.06)</td><td>2.87</td><td>2.14</td><td>15.2</td></tr><tr><td>2018-01-14.333</td><td><sup>13</sup>CO</td><td>3-2</td><td>HARP stare</td><td>0.16</td><td>900</td><td>&lt;229</td><td>...</td><td>2.87</td><td>2.14</td><td>15.2</td></tr><tr><td>2018-01-14.333</td><td>C<sup>18</sup>O</td><td>3-2</td><td>HARP stare</td><td>0.16</td><td>900</td><td>&lt;282</td><td>...</td><td>2.87</td><td>2.14</td><td>15.2</td></tr><tr><td>2018-01-14.367</td><td>CO</td><td>3-2</td><td>HARP jiggle</td><td>0.15</td><td>1780</td><td>834(114)</td><td>-0.25(0.04)</td><td>2.87</td><td>2.14</td><td>15.2</td></tr><tr><td>2018-01-14.384</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.15</td><td>120</td><td>881(140)</td><td>-0.27(0.12)</td><td>2.87</td><td>2.14</td><td>15.2</td></tr><tr><td>2018-01-15.197</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.09</td><td>120</td><td>1084(130)</td><td>-0.21(0.06)</td><td>2.87</td><td>2.14</td><td>15.4</td></tr><tr><td>2018-01-15.213</td><td>CO</td><td>3-2</td><td>HARP jiggle</td><td>0.08</td><td>1780</td><td>1145(131)</td><td>-0.23(0.02)</td><td>2.87</td><td>2.14</td><td>15.4</td></tr><tr><td>2018-01-15.231</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.08</td><td>120</td><td>1028(117)</td><td>-0.23(0.05)</td><td>2.87</td><td>2.14</td><td>15.4</td></tr><tr><td>2018-01-15.253</td><td><sup>13</sup>CO</td><td>3-2</td><td>HARP stare</td><td>0.08</td><td>900</td><td>&lt;84</td><td>...</td><td>2.87</td><td>2.14</td><td>15.4</td></tr><tr><td>2018-01-15.253</td><td>C<sup>18</sup>O</td><td>3-2</td><td>HARP stare</td><td>0.08</td><td>900</td><td>&lt;112</td><td>...</td><td>2.87</td><td>2.14</td><td>15.4</td></tr><tr><td>2018-01-15.271</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.08</td><td>120</td><td>1015(114)</td><td>-0.19(0.05)</td><td>2.87</td><td>2.14</td><td>15.4</td></tr><tr><td>2018-01-20.247</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.06</td><td>600</td><td>745(76)</td><td>-0.23(0.03)</td><td>2.85</td><td>2.18</td><td>16.8</td></tr><tr><td>2018-01-25.322</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.14</td><td>600</td><td>881(96)</td><td>-0.23(0.04)</td><td>2.83</td><td>2.23</td><td>17.9</td></tr><tr><td>2018-01-31.338</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.08</td><td>600</td><td>910(93)</td><td>-0.23(0.03)</td><td>2.80</td><td>2.29</td><td>19.1</td></tr><tr><td>2018-02-01.278</td><td>CO</td><td>3-2</td><td>HARP stare</td><td>0.14</td><td>600</td><td>768(84)</td><td>-0.18(0.04)</td><td>2.80</td><td>2.30</td><td>19.3</td></tr></tbody></table>

**Note.**

Table 1 shows a summary of the JCMT spectral line observations, including the zenith opacity at 225 GHz ($\tau_0$) and on-source integration time per observation (Int.). Integrated line intensities ($\int T_{MB}dv$) and centroid velocities ($\bar{v}$) are also given (with $1\sigma$ uncertainties in parentheses). The uncertainties on $\int T_{MB}dv$ include the 10% intensity calibration error, added in quadrature with the statistical error.

## 2.2. Submillimeter Array

Interferometric observations of C/2016 R2 were made on 2018 February 21 using the Submillimeter Array (SMA), when the comet was at $r_H = 2.7$ au; $\Delta = 2.5$ au; $\phi = 21^\circ$. Five of the SMA (6 m) antennas were online and in the subcompact (SUB) configuration at the time of the observations, resulting in a spatial resolution of $\approx 4'' \times 7''$ at 230 GHz. The two SMA receivers were both tuned to 1.3 mm to cover the CO $J = 2-1$ transition at 230.538 GHz. The SMA Wideband Astronomical ROACH2 Machine (SWARM) correlator provides 8 GHz bandwidth per sideband, divided into four equal-sized chunks with a uniform spectral resolution of 140 kHz (0.18 km s<sup>-1</sup> at 230 GHz). Position and Doppler-tracking of the comet over a 5 hr observing period was performed using JPL Horizons orbital solution number 14, and compensated for by the correlator in real time.

Calibration of visibility phases and amplitudes was performed using periodic observations of quasars 0336 + 323 and 3c111, at 15 minute intervals. Measurements of Uranus were used to obtain an absolute scale for calibration of the flux densities. All data were phase- and amplitude-calibrated using the Millimeter Interferometer Reduction (MIR) software package.<fnref n="11" /> The calibrated SMA data were then exported in

uvfits format for subsequent imaging and analysis using the Common Astronomy Software Applications (CASA) software package version 5.1 (Jaeger 2008).

Interferometric imaging and deconvolution was performed using the Hogbom clean algorithm, with a 20" mask centered on the peak of CO emission and a threshold of twice the rms noise per channel ($\sigma = 0.19$ Jy). Natural visibility weighting was used and the pixel size was set to 0.5".

# 3. Results

## 3.1. JCMT CO Spectral Line Time Series

A time series of the observed JCMT CO $J = 3-2$ spectra is shown in Figure 1. Individual spectra have been Doppler-corrected for the comet's radial velocity with respect to the observer, and were baseline-subtracted using low-order polynomial fits to the emission-free regions. Spectra in this figure were selected based on on-source integration times of at least 600 s, to exclude the noisier data obtained from shorter-integration spectra on some of the observing dates.

The comet's radial velocity in the cometocentric rest frame increased steadily throughout the period of observations, but the double-peaked line profile, consisting of a strong, narrow blueshifted peak and a weaker redshifted peak remained apparently constant (within the noise). Given the relatively small phase angle ($\phi$) of our observations, such blueshifted emission is explained as a result of enhanced outgassing from the side of the nucleus facing the Sun (and Earth), compared with the nightside.

To search for temporal variability in the comet's activity, the spectrally integrated emission line area ($\int T_{MB}dv$) is plotted as a function of time in Figure 2. No significant time variability was detected, implying that the CO outgassing rate

<footnote>a Integrated over the three strongest CH<sub>3</sub>OH transitions (see Section 4.6).</footnote>

<footnote n="11">¹¹ https://lweb.cfa.harvard.edu/~cqi/mircook.html</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Time sequence of CO *J* = 3 − 2 spectra observed using JCMT (in the cometocentric rest frame), shown with additive baseline offsets. The observing date for each spectrum is given in the format YYYYMMDD. For clarity, only spectra with an on-source observing time of at least 600 s are shown. Best-fitting spectral models are overlaid with red curves (see Section 4.3).</figcaption></figure>

and rotational temperature (*T*<sub>rot</sub>) were probably close to steady
state during the time period of our observations, although
variations in both quantities could have canceled each other
out to some degree—the CO *J* = 3 − 2 line intensity varies by
a factor of 1.5 between *T*<sub>rot</sub> = 15−25 K (a plausible range for
this comet based on the findings of Biver et al. (2018); see
also Section 4.2).

3.2. CO Single-dish Mapping with JCMT HARP

HARP jiggle-map spectral cubes from 2018 January 14 and
2018 January 15 were Doppler-corrected and averaged together
(with rejection of masked pixels). The spectrally integrated
(moment 0) map of CO J = 3 − 2 emission is shown in
Figure 3, with contour levels plotted in units of 3σn<sub>c</sub><sup>−0.5</sup>, where

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Spectrally integrated CO *J* = 3 − 2 line intensities observed using the JCMT as a function of time, based on the spectra shown in Figure 1. The horizontal line shows the error-weighted mean. Error bars include the statistical uncertainty, with an additional 10% calibration uncertainty added in quadrature.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Spectrally integrated CO *J* = 3 − 2 map for comet C/2016 R2 (PanSTARRS), from the average of HARP jiggle observations on 2018 January 14 and 2018 February 15. The FWHM of the circular Gaussian JCMT beam is indicated upper right, and the sky-projected solar and orbital trail (negative of the comet’s velocity) vectors are shown in the lower right. Contours are in units of 3σ and the axes are aligned with the equatorial (R.A./decl.) grid.</figcaption></figure>

$\sigma$ is the average rms noise of the data cube (equal to $T_{\text{MB}}$ =
72 mK km s⁻¹) and $n_c$ is the number of spectral channels in the
moment 0 integral. The CO $J=3-2$ integrated line intensity,
averaged over the entire 120" × 120" HARP data cube,
was 164 ± 19 mK km s⁻¹ on 2018 January 14 and 202 ±
22 mK km s⁻¹ on 2018 January 15. These values are sufficiently
similar (given the uncertainties) to justify combining the two
data sets. The CO spatial distributions were also apparently
identical on the two dates, with consistent radial intensity
profiles (within the uncertainties).

The CO coma is diffuse and spatially extended compared
with the 22,000 km JCMT beam, with some weak evidence for
deviations from circular symmetry about the nucleus. Emission
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Spectrally integrated CO *J* = 2 − 1 emission map for comet C/2016 R2, obtained using the SMA on 2018 February 21. The FWHM (and orientation) of the elliptical Gaussian restoring beam is indicated (upper right), and the sky-projected solar and orbital trail vectors are shown in the lower right. Contours are in units of 3σ and the axes are aligned with the equatorial (R.A./decl.) grid, with the origin at the SMA phase tracking center. Negative contours are shown with a dashed line style.</figcaption></figure>

is detected (at the 3σ level) up to a radial distance 82,000 km west of the image center, which is in a direction similar to the sky-projected sunward vector (103° clockwise from north), whereas the 3σ contour extends only 58,000 km to the east. The second contour (at 6σ) also shows an excursion in the sunward direction into a pixel with an intensity 2.2σ larger than the mean of the other pixels at the same cometocentric distance. These features provide tentative evidence for preferential outgassing on the illuminated (sunward) side of the nucleus. The (normalized, azimuthally averaged) CO spectral line profile is plotted as a function of distance from the center of the image in Appendix D (Figure 16), and reveals a consistent excess in the blueshifted emission out to cometocentric distances of at least $r_c \sim 80,000$ km.

### 3.3. CO Interferometric Mapping with the SMA

The CO *J* = 2 − 1 intensity map observed using SMA is shown in Figure 4, integrated over the velocity width of the detected emission. The intensity reaches a peak 2.15 west of the phase center, which may be partly a result of asymmetrical outgassing in the sunward direction, or errors in the position of the comet nucleus compared with the JPL Horizons ephemeris orbital solution.

The coma shows an extended morphology in an approximately north–south direction, and is less well resolved in the east–west direction, where significant large-scale flux appears to have been resolved out by the interferometer, resulting in negative side lobes (regions with dashed contours) apparent on either side of the comet. Although the orientation of the spatially extended emission (defined by the outermost 3σ contour) matches closely the direction of the comet’s (sky-projected) orbital trail, it also aligns with an axis of strong artifacts in the interferometric point-spread function, so the reality of this asymmetric, extended feature remains questionable.

## 4. Radiative Transfer Modeling

As a result of near-spherical expansion, cometary comae span an extremely broad range of densities over a short distance. Consequently, their gases are often subject to a range of excitation conditions within a single telescope beam, governed by a balance of microscopic collisional and radiative processes (Bockelée-Morvan et al. 2004a), and are generally not in local thermodynamic equilibrium (LTE). To interpret cometary rotational spectra therefore requires detailed excitation and radiative transfer modeling. Here we introduce a new code called SUBLIME for simulating the rotational emission lines from cometary coma molecules in three dimensions (two spatial and one spectral). The SUBLIME code is based on the open-source version 1.9.3 of the LIME (LIne Modeling Engine) code by Brinch & Hogerheijde (2010).<fnref n="12" />

The basic equations of radiative transfer and excitation used in our model are described in Appendix A. Some recently published coma radiative transfer models have invoked the steady-state approximation (e.g., Bøgelund & Hogerheijde 2017; de Val-Borro et al. 2018; Cordiner et al. 2019), setting $dN_i/dt = 0$ in Equation (A4), which allows the energy-level populations to be solved independently (in parallel) at a large number of discrete positions within the region of interest. This approximation facilitates the treatment of complex (3D) coma morphologies, but comes at the expense of discarding the effects of the outflow dynamics, which can be important for molecules such as CO, with small dipole moments and hence slow rotational transitions relative to the dynamical timescale (see Section 5.4).

LIME employs the steady-state approximation. We therefore substantially modified the code to enable the more physically accurate, time-dependent solution of Equation (A4), using the CVODE solver (Hindmarsh 2019) to calculate the molecular excitation along radial trajectories of the outflowing coma gases. The time-dependent solution has previously been implemented in the models of Bockelée-Morvan (1987) and Biver et al. (1999), and allows the temporal evolution of molecular excitation in the rapidly expanding coma to be properly accounted for. The time-dependent version of our code, as used in the present study, ignores the impact of opacity/photon trapping on the molecular excitation, which is negligible for the molecules considered in this study (see Appendix A).

SUBLIME calculates the molecular excitation along radial vectors, which are then interpolated onto an unstructured 3D grid (Delaunay 1934), configured with a density of grid points proportional to the gas density. The broad range of density and size scales in the coma (covering many orders of magnitude, from the ~kilometer-sized nucleus to the ~10⁶ km-scale outer coma) can thus be sampled much more efficiently than with a uniform grid spacing. To simulate coma asymmetries and jets using this method, the spatial domain is divided into multiple solid-angle regions ($\Omega_i$), each with its own outflow velocity ($v_i$), kinetic temperature ($T_i$), and molecular production rate ($Q_i$). A separate CVODE calculation is performed for each solid-angle region. In the present paper, we use two regions, corresponding to (1) the ambient coma and (2) a conical jet with its apex at the center of the nucleus.

For all models, we used a Delaunay grid with 10,000 points, between $r_c=500$ m (the assumed radius of the nucleus) and an

<footnote n="12">¹² Available at https://github.com/lime-rt/lime/releases.</footnote>
outer radius of $10^6$ km. Tests showed that this grid density was sufficient to produce reliable spectral line models at the resolution of our observations; adding more grid points did not significantly change the results. The precise grid point locations are selected pseudo-randomly for each model run (Brinch & Hogerheijde 2010), but we fixed the random number generator seed so that an identical grid was produced every time, providing numerical stability in order to facilitate reliable (repeatable) parameter retrievals. During ray tracing, we employed the LIME traceray.smooth algorithm, which interpolates the level populations between grid points, thus reducing grid-related artifacts in the output image. A pixel size of 0.5 and a channel spacing of 25 m s$^{-1}$ were chosen for the model images to sufficiently sample the spatial and spectral-resolution elements of the JCMT and SMA observations. To accurately capture the rapid (nonlinear) flux increase on subpixel scales toward the center of the image, due to the strongly increasing coma density with decreasing $r_c$, we employed cartesian supersampling (on a regularly spaced grid of 30 × 30 rays) for each pixel within the central 4 × 4 pixel region of each image.

## 4.1. CO-CO Collision Rate Coefficients

As the dominant coma gas (McKay et al. 2019), CO is the primary collision partner in our model, and is therefore the main species responsible for the redistribution of thermal energy among the rotational states (J) of the observed gases. Knowledge of the CO-CO collision rate coefficients ($k_{J_1J_2}$) is therefore required to correctly model the CO emission from the comet. Previous studies (Biver et al. 2018; Wierzchos & Womack 2018) made gross approximations for these rates, so their results remain uncertain. Here we employ quantum calculations to model the CO-CO collisions, allowing us to accurately determine the CO excitation for the first time in a CO-dominated cometary coma.<fnref n="13" />

To describe the interaction between colliding CO molecules, we used the 4D potential energy surface (PES) with rigid CO molecules calculated by Vissers et al. (2003). The PES was calculated using the coupled-cluster single, double, and perturbative triples (CCSD(T)) method with augmented triple zeta basis (aug-cc-pVTZ). The accuracy of the PES was benchmarked with respect to experimental studies (Surin et al. 2007; Sun et al. 2020). The scattering calculations were performed with the Moslcat code under the assumption of distinguishable particles (Hutson & Green 1994). A series of tests was performed, revealing that the most accurate close-coupling (CC; Green 1975) approach would not be feasible in terms of computer memory and processing time. We therefore explored the possibility of using the coupled-states (CS) approximation (see Klos & Lique 2018 for a review of these methods). The differences between CC and CS were found, on average, to be less than a factor of 1.5–2 (and never higher than a factor of 3).

Assuming that we can distinguish the CO molecules, the first one is the target (with the rotational state characterized by quantum number $J_1$), and the second is the collider ($J_2$). Collisional rate coefficients $k_{J_1J_2 \to J_1'J_2'}$ were then computed by averaging the cross-sections over the Boltzmann distribution of

collisional energies (Equation (1)),

$$k_{J_1 J_2 \to J_1' J_2'}(T) = \left(\frac{8}{\pi \mu k_B^3 T^3}\right)^{\frac{1}{2}} \int_0^\infty \sigma_{J_1 J_2 \to J_1' J_2'} (E_c) E_c e^{-\frac{E_c}{k_B T}} dE_c, \quad (1)$$

where $\mu$ is the reduced mass of the system, $k_B$ is a Boltzmann constant, $\sigma$ is the cross-section of a given transition, and $E_c$ is the collisional energy. Rate coefficients used in the radiative transfer model ($k_{J_1 \to J_1'}(T_{kin})$; see Appendix B, Table 3) were calculated by averaging over a thermal rotational distribution for the initial excitation state of the collider, and summed over its final state as follows:

$$k_{J_1 \to J_1'}(T) = N_{J_2=0} \sum_{J_2'} k_{J_1, J_2=0 \to J_1', J_2'} + N_{J_2=1} \sum_{J_2'} k_{J_1, J_2=1 \to J_1', J_2'} + \dots \quad (2)$$

where summation is over all possible final states of the collider. $N_{J_2}$ is the (Boltzmann) energy-level population of the collider, given by

$$N_{J_2} = \frac{(2J_2 + 1) e^{-\frac{E_{J_2}}{k_B T}}}{\sum_{J_2} (2J_2 + 1) e^{-\frac{E_{J_2}}{k_B T}}}, \quad (3)$$

where $E_{J_2}$ is its energy level. Calculations were performed up to a collision energy of 300 cm$^{-1}$ and maximum rotational quantum number $J_1 = J_2 = 5$, leading to rate coefficients for kinetic (and rotational) temperatures valid up to 30 K. Such a restricted set of rate coefficients is appropriate for modeling the CO excitation in C/2016 R2 due to the relatively low coma kinetic temperature (~20 K; Biver et al. 2018); at this temperature, energy levels above $J=5$ comprise less than 1% of the total CO population in the collisionally dominated zone. Our model for C/2016 R2 includes levels up to $J=40$, but for $J>5$, the CO-CO collision rates are assumed to be the same as for para-H$_2$ colliding with CO (Yang et al. 2010). A more detailed explanation of the collision rate calculation method, including calculation of rate coefficients for kinetic temperatures up to 100 K, is currently in preparation.

## 4.2. Modeling the JCMT CO J = 2 − 1 and J = 3 − 2 Data

On 2018 January 14, the CO $J=2-1$ and $J=3-2$ lines were observed in close succession (within 1.3 hr of each other). Given the lack of significant temporal variability in the $J=3-2$ line strength around this date (Figure 1), differences between the strengths of these two lines can be assumed to originate as a result of (1) their different intrinsic line strengths and (2) differences in the excitation of the upper-state energy level. According to Equation (A4), the level populations depend on the collision rates $k_{ijn}$ (where n is the local density), which in turn depend on the CO production rate, coma outflow velocity, and kinetic temperature. Modeling the two CO lines simultaneously therefore provides a diagnostic of the coma temperature, while the line profile provides information on the coma outflow velocity along the line of sight.

There is insufficient information from these relatively low signal-to-noise, 1D spectra to infer the entire coma physical structure, so we adopt a modified Haser (1957) model. The

<footnote n="13">¹³ Ndengué et al. (2015) previously studied the energy transfer in CO-CO collisions, but their results were only partly converged and contain a systematic error for some (but not all) transition rates.</footnote>
assumption of isotropically expanding gas, with a constant production rate and outflow velocity (as described by the 1D Haser model), is routinely used for analysis of cometary spectra and images observed across the range of wavelengths (e.g., Bockelée-Morvan et al. 2004a; Cochran et al. 2012; Cordiner et al. 2014). In the case of strongly asymmetric outgassing, however, as observed for C/2016 R2, the assumption of spherical symmetry is no longer applicable, so a more flexible 3D model is required.

In the present study, we adopt the simplest physically reasonable model capable of fitting the observed data. We consider two different outflow components, C₁ and C₂, corresponding to solid-angle regions Ω₁, Ω₂ with independent production rates (Q₁, Q₂) and outflow velocities (v₁, v₂). This assumes that the nucleus can be divided into two different activity regimes: (1) an ambient outflow from the majority of the (thermally activated) sublimating area of the nucleus, and (2) enhanced gas production from a (set of) spatially confined vent(s) or jet(s), in particular, on the sunward-facing side of the nucleus. Similar two-component models (involving a sunward jet and ambient/isotropic coma, or a hemispherically asymmetric outflow) have been invoked previously to explain asymmetries in high-resolution spectral line profiles from several comets, including 29P/SW1 (Festou et al. 2001), O1/Hale-Bopp (Gunnarsson 2003), 19P/Borelly (Bockelée-Morvan et al. 2004b), 2I/Borisov (Cordiner et al. 2020), 46P/Wirtanen (Roth et al. 2021b), and C/2015 ER61 (Roth et al. 2021a). The increased production rate and outflow velocity measured on the sunward side of these comets is consistent with the results of fluid dynamic and Monte Carlo coma models (Crifo et al. 1999; Fougere et al. 2016), and arises as a result of elevated temperatures in both the coma and the nucleus, due to increased solar insolation. While more complex coma parameterizations can be envisaged, the simplest model capable of fitting the data should be the best constrained, with the additional benefit of being more efficient to configure and run.

We construct a SUBLIME model to simultaneously fit the CO *J* = 2 − 1 and *J* = 3 − 2 line profiles, assuming a constant coma kinetic temperature (after Biver et al. 2018). Component C₁ is defined by a conical region about the subsolar point of the nucleus, with half-opening angle θ, whereas C₂ is the remaining, ambient coma. We used the MPFIT nonlinear least-squares routine (Markwardt 2012) to find the optimal set of parameters (Q₁, Q₂, v₁, v₂, θ). The model images were convolved with the (Gaussian) JCMT beam pattern and normalized by their respective beam efficiency factors before comparing them with the observations. Statistical (1σ) error estimates are obtained for each parameter from the diagonal elements of the MPFIT covariance matrix.

The best-fitting model spectra are shown in Figure 5, corresponding to $Q_1 = (1.8 \pm 0.4) \times 10^{28} \text{ s}^{-1}$, $Q_2 = (2.2 \pm 0.3) \times 10^{28} \text{ s}^{-1}$, $v_1 = 0.50 \pm 0.01 \text{ km s}^{-1}$, $v_2 = 0.29 \pm 0.02 \text{ km s}^{-1}$, $\theta = 27^\circ \pm 5^\circ$, and $T_{\text{kin}} = 18.7 \pm 2.3 \text{ K}$. The retrieved kinetic temperature is consistent with the range of values 18.6 ± 2.6 to 24.2 ± 7.9 derived by Biver et al. (2018) using multiple lines of CH₃OH on 2018 January 23–24. The total CO gas production rate from our model is $Q_\text{\textit{CO}} = (3.9 \pm 0.7) \times 10^{28} \text{ s}^{-1}$, and the ratio of production rates per unit solid angle between the jet and ambient coma is $R_\text{\textit{Q}} = (Q_1/\Omega_1)/(Q_2/\Omega_2) = 14.5 \pm 3.9$. Given the strong dependence of the production rate for both components on the size of their solid-angle regions, $R_\text{\textit{Q}}$ is a more physically

<figure><img src="image_8.png" /><figcaption>**Figure 5.** JCMT spectra of CO *J* = 3 − 2 (top) and *J* = 2 − 1 (bottom), in the cometocentric rest frame. The best-fitting SUBLIME model is overlaid with the red curves.</figcaption></figure>

meaningful quantity than the simple ratio of production rates (Q₁/Q₂), and reveals the degree to which the comet’s activity is enhanced due to heating by the Sun in the vicinity of the subsolar point. The overall quality of fit is good considering the noise, and reproduces well the asymmetry of the *J* = 3 − 2 line. We also performed fits allowing the angle of the jet axis to vary with respect to the Sun–comet vector, but the quality of the fit was not significantly improved. Additional models were run allowing different $T_{\text{kin}}$ values for the two coma components, with best-fitting results $T_{\text{kin1}} = 19.3 \pm 3.7 \text{ K}$ and $T_{\text{kin2}} = 17.9 \pm 4.8 \text{ K}$, but again, the overall quality of the fit was not improved. The consistency of $T_{\text{kin1}}$ and $T_{\text{kin2}}$ (within their respective errors) provides further justification for adopting a uniform $T_{\text{kin}}$ value throughout the rest of this study.

### 4.3. Modeling the Average and Time-series JCMT CO
*J* = 3 − 2 Spectra

Given the relative constancy of the comet’s CO emission over time, a useful estimate for the total (time-averaged) CO production rate over the course of our observations was obtained by taking the average of all the JCMT CO (3 − 2) spectral data listed in Table 1. After Doppler-correcting each spectrum to the cometocentric rest frame and weighting them by $1/\sigma^2$, the resulting average spectrum was modeled using the same procedure described above. The best-fitting model is shown in Figure 6, and corresponds to $Q_\text{\textit{CO}} = (5.3 \pm 0.2) \times 10^{28} \text{ s}^{-1}$, $R_\text{\textit{Q}} = 7.4 \pm 0.4$, $v_1 = 0.51 \pm 0.01 \text{ km s}^{-1}$,
<table><caption>**Table 2**
Best-fit Model Results for Individual JCMT (and SMA) Spectra of C/2016 R2</caption><thead><tr><th>Date</th><th>Q<sub>r</sub><br>(10<sup>28</sup> s<sup>-1</sup>)</th><th>R<sub>Q</sub></th><th>v<sub>1</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>2</sub><br>(km s<sup>-1</sup>)</th><th>θ<br>(deg.)</th><th>φ<sub>jet</sub><br>(deg.)</th></tr></thead><tbody><tr><td>2018-01-13</td><td>4.5 (0.7)</td><td>6.8 (4.5)</td><td>0.51 (0.01)</td><td>0.23 (0.02)</td><td>60 (5)</td><td>18 (6)</td></tr><tr><td>2018-01-14</td><td>4.0 (0.9)</td><td>10.8 (2.9)</td><td>0.50 (0.01)</td><td>0.29 (0.02)</td><td>32 (7)</td><td>7 (10)</td></tr><tr><td>2018-01-15</td><td>7.4 (1.6)</td><td>9.6 (3.1)</td><td>0.52 (0.02)</td><td>0.31 (0.03)</td><td>77 (13)</td><td>37 (19)</td></tr><tr><td>2018-01-20</td><td>3.9 (1.8)</td><td>17.2 (10.8)</td><td>0.50 (0.01)</td><td>0.24 (0.01)</td><td>42 (19)</td><td>34 (11)</td></tr><tr><td>2018-01-25</td><td>5.0 (0.9)</td><td>10.4 (1.7)</td><td>0.53 (0.01)</td><td>0.24 (0.02)</td><td>58 (6)</td><td>20 (7)</td></tr><tr><td>2018-01-31</td><td>5.7 (0.8)</td><td>7.3 (0.7)</td><td>0.53 (0.01)</td><td>0.26 (0.01)</td><td>64 (3)</td><td>12 (4)</td></tr><tr><td>2018-02-01</td><td>4.4 (0.8)</td><td>4.9 (0.9)</td><td>0.49 (0.01)</td><td>0.33 (0.02)</td><td>46 (8)</td><td>8 (12)</td></tr><tr><td>JCMT avg.</td><td>5.3 (0.6)</td><td>7.4 (0.4)</td><td>0.51 (0.01)</td><td>0.25 (0.01)</td><td>62 (2)</td><td>24 (2)</td></tr><tr><td>2018-02-21<sup>a</sup></td><td>6.7 (0.9)</td><td>9.2 (1.5)</td><td>0.64 (0.02)</td><td>0.20 (0.02)</td><td>74 (7)</td><td>21<sup>b</sup></td></tr></tbody></table>

**Notes.** $1\sigma$ statistical uncertainties on each value are given in parentheses. Uncertainties on $Q_r$ include a 10% intensity calibration error, added in quadrature with the statistical error.

<sup>a</sup> SMA best-fitting parent visibility model.

<sup>b</sup> The jet axis set to direction of the comet–Sun vector due to a lack of constraints caused by the lower spectral resolution of the SMA data.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Weighted average CO $J = 3-2$ spectrum from the JCMT HARP pointing receptor (in the cometocentric rest frame). The best-fitting SUBLIME 3D model is overlaid using a red curve. The dotted-dashed gray curve shows the best-fitting 1D model, assuming a spherically symmetric coma.</figcaption></figure>

$v_2 = 0.25 \pm 0.01$ km s$^{-1}$, and $\theta = 62^\circ \pm 2^\circ$. In this model the direction of the jet's axis with respect to the line of sight ($\phi_{jet}$) was also allowed to vary, and the best fit was for $\phi_{jet} = 24^\circ \pm 2^\circ$ ($7^\circ \pm 2^\circ$ away from the mean Sun–comet vector). These values are within $2\sigma$ of the parameters derived for the simultaneous fit to the $J = 2-1$ and $J = 3-2$ data on 2018 January 14. Allowing for an additional 10% calibration error, the total uncertainty on the CO production rate is $0.6 \times 10^{28}$ s$^{-1}$.

Individual fits were also performed to the time series of the spectra presented in Figure 1. A good fit was obtained in each case, and the results are given in Table 2. While the majority of these individual fit results are consistent (within $1-2\sigma$) with the average JCMT spectrum, notable deviations include the significantly narrower jet opening angle of $\theta = 32^\circ \pm 7^\circ$ on January 14. On February 1, the ambient coma showed an enhanced outgassing velocity (with $v_2 = 0.33 \pm 0.02$ km s$^{-1}$), which was also accompanied by a reduction in $R_Q$ to $4.9 \pm 0.9$. Such variability in the detailed outflow morphology on specific epochs implies some inhomogeneity of the nucleus (and/or its

heating rate), leading to modest changes in the CO outgassing behavior over time.

### 4.4. Modeling the CO Spatial Distribution

To investigate the spatial morphology of the CO coma, we generated 3D (spectral–spatial) models for comparison with the JCMT HARP image cube. First, a SUBLIME model fit was performed for the CO $3-2$ spectrum extracted from the central map pixel (coinciding with the intensity peak in Figure 3), by optimizing the model parameters $Q_1$, $Q_2$, $v_1$, $v_2$, $\theta$, and $\phi_{jet}$. The model orientation (in the plane of the sky) was fixed so that the jet azimuth angle matched the direction of the sky-projected comet–Sun vector. This best-fitting CO parent model was then convolved by the JCMT beam shape, multiplied by the main beam efficiency and spectrally integrated to obtain the 2D map in the left panel of Figure 7.

Comparison of the contours in this figure with the observed CO map (Figure 3) shows that while the central intensity peak is accurately reproduced, the outermost contour of the observations lies, on average, well outside that of the model. This is more clearly demonstrated by comparing the azimuthally averaged radial profiles ($\langle T_{MB} \rangle_{az}$) of the modeled and observed CO emission, which are plotted as a function of distance from the nucleus in Figure 8. The best-fitting parent model (green dotted line) falls off more rapidly with distance than the observed data (filled black circles). This implies the presence of an excess of CO emission at large radii, which cannot be explained solely by CO released from the nucleus at a constant outflow velocity. The JCMT spectral baselines were well behaved over the course of our observations, and care was taken to ensure robust baseline removal through low-order polynomial subtraction, giving confidence regarding the reality of the observed extended CO emission (see also Figure 16).

We quantify the statistical significance of the extended emission using the reduced chi-square statistic $\chi_R^2 = \sum_i (y_i - y_m)_i / d$, where $y_i$ and $y_m$ are the observed and modeled intensities, respectively, and $d$ is the number of degrees of freedom (equal to the number of independent data points minus the number of free model parameters). The model $\chi_R^2$ value is given in the legend of Figure 8, along with the associated probability ($P$) that the difference between the model and observations is due to statistical chance. With $\chi_R^2 = 2.1$
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Spectrally integrated model CO *J* = 3 − 2 maps (for comparison with Figure 3), assuming CO solely as a parent molecule (left panel), CO produced from the nucleus plus CO₂ photolysis, with *Q*(CO₂)/*Q*(CO) = 6.5 (center panel), and CO from the nucleus, with an additional CO shell at *r*<sub>c</sub> ~ 1.4 × 10⁵ km (right panel). Contour spacings are the same as in Figure 3.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Azimuthal average (in 1-pixel bins) of the JCMT CO *J* = 3 − 2 emission map (Figure 3) as a function of distance from the central pixel (black circles). Best-fitting SUBLIME model results are overlaid, including Haser-type parent and daughter models (see the text). An improved fit is obtained at large radii when including a CO source from CO₂ photolysis (with a large *Q* (CO₂)/*Q*(CO) ratio of 6.5) or with the inclusion of an extended CO shell at *r*<sub>s</sub> = 1.2 × 10⁵, with a density enhancement factor of *f* = 1.8.</figcaption></figure>

and *P*=0.07, the CO parent model does not represent a good fit to the observations.

We consider three possibilities to account for the CO excess:
(1) a CO *J* = 3 level population that increases more rapidly
with distance than predicted, (2) an additional Haser-type
extended CO source, and (3) a step-like enhancement in the
coma density at a large distance from the nucleus. A reduction
in the CO photodissociation rate (ΓCO) would not reproduce the
observed emission excess—at rH = 2.85 au, ΓCO = 9.2 × 10−8
s−1 (Huebner & Mukherjee 2015), so at an outflow velocity of
0.2 km s−1, the CO scale length is 2.2 × 106 km. CO
photodissociation is therefore negligible within the 80,000 km
JCMT field of view, so reducing ΓCO by a factor of a few (its
range of uncertainty) has no noticeable impact on the modeled
CO density profile.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Fractional energy-level populations as a function of radius for the lowest seven CO rotational levels, using our SUBLIME time-dependent excitation model (solid curves) and the LIME steady-state solver (dotted curves). For clarity, only the populations from the higher-density jet component of the model (C1) is shown. The steady-state model’s failure to account for coma dynamics leads to the onset of non-LTE effects too close to the nucleus.</figcaption></figure>

For option 1, we need to consider the CO rotational level
populations. As shown by Figure 9, the population of the *J* = 3
level (red curve) remains relatively constant as a function of
distance from the nucleus. This is because the effective
pumping rates (*G*<sub>ij</sub>) into this level are closely balanced by
radiative transitions out of the level, so on the larger distance
scales, ≳30,000 km at which the excess CO emission becomes
most visible, the *J* = 3 level population is already close to the
value attained at fluorescence equilibrium, making excitation
effects an unlikely explanation for the observed extended
emission.

On smaller distance scales comparable with the size of the
(central) JCMT beam, a lower rotational temperature could
reduce the CO 3 − 2 intensity, leading to a shallower radial
profile in better agreement with the observations. However,
assuming T<sub>kin</sub> = 16.4 K (the 1σ lower limit derived in
Section 4.2), the best-fitting χ<sup>2</sup> value is 1.6 with *P* = 0.17,
which still does not represent a very good fit. Lower values of
T<sub>kin</sub> also appear less likely given the range of possible T<sub>kin</sub>
values (16–32 K) observed by Biver et al. (2018), so we seek
alternative explanations for the shape of the ⟨T<sub>MB</sub>⟩<sub>az</sub> profile.

The presence of an additional CO source (option 2) is worth
considering, in light of previous evidence for extended CO
sources in cometary comae (Cottin & Fray 2008). However, the JCMT map is inconsistent with CO produced solely as a (Haser-type) photochemical daughter species. Our best-fitting SUBLIME model was modified to produce CO from photo-dissociation of an (unknown) parent molecule, with the production rate ($Q_p$) and parent photodissociation rate $\Gamma_p$ free to vary (where $\Gamma_p = v_i/L_p$; $v_i$ is the outflow velocity in conical region $i$ and $L_p$ is the parent scale length). The best-fitting, azimuthally averaged CO daughter model (with $\Gamma_p = 7.3 \times 10^{-5}$ s⁻¹, $Q_p = 7.0 \times 10^{28}$ s⁻¹) is shown in Figure 8, and systematically underfits the observed CO emission at large radii (with $\chi_R^2 = 1.6$).

A better fit to the observations can be obtained using a composite CO density profile, where CO is included as both a parent and a daughter species. Figures 7 and 8 show the best-fitting results from such a model, with CO originating from the nucleus at a total production rate of $Q_t = 5.5 \times 10^{28}$ s⁻¹ and an additional CO source in the coma from CO₂ photolysis (CO₂ + $h\nu$ → CO + O at $\Gamma_{CO_2} = 1.5 \times 10^{-7}$ s⁻¹; Huebner & Mukherjee 2015). Although this composite model may appear plausible at first glance, and accurately reproduces the CO emission profile, it requires a CO₂ production rate 6.5 times larger than that of CO. Considering that McKay et al. (2019) derived an upper limit of $Q(\text{CO}_2) \lesssim 1.5 \times 10^{28}$ s⁻¹ based on Spitzer observations from 2018 February 12 and 2018 February 21 (i.e., 3.7 times less than the parent CO source), it seems unlikely that the comet could have been producing so much CO₂ at the time of our JCMT observations only a few weeks prior. Indeed, the large production rate required for any distributed molecular source to adequately fit the excess CO emission at large radii renders this an unlikely scenario in general, regardless of the assumed CO parent.

The possibility of non-Haser-type extended sources should also be considered, which could produce different radial density profiles, potentially resulting in a good fit to the observations. For example, Gunnarsson et al. (2002) developed a model for the sublimation and fragmentation of CO-rich icy grains to explain the extended CO distribution observed in comet 29P. Such detailed physical modeling is beyond the scope of our present study, however, but could be usefully investigated in a future article.

The third explanation we consider for the CO excess is an increase in coma density at large radii, which could result from temporal modulation of the gas production and/or outflow velocity, or a rapid onset of icy grain sublimation far from the nucleus. Slowing of the gas expansion rate in the outer coma (for example, due to sublimation of slow-moving icy grains—see, e.g., Ip 1986; Fougere et al. 2012), would also lead to an increase in gas density. From Appendix D, Figure 16, the CO line profile does not show any obvious evidence for coma deceleration with increasing distance, although we note that the weakness of the extended CO emission shell (representing ~1/6 of the parent source emission at $r_c = 60, 000$ km) would make its velocity signature difficult to detect given the noise.

To model such an outer CO shell, we implemented a multiplicative increase in the coma density $n$ at radius $r_s$, smoothed by an exponential function such that $n(r) = 1 + (f-1)/(1 + e^{(r_s-r)/w_s})$. The CO initial abundance, density enhancement factor ($f$), step radius ($r_s$), and smoothing width ($w_s$) were optimized to obtain the best fit to $\langle T_{MB} \rangle_{az}(r)$, and the resulting azimuthally averaged model emission profile and moment 0 map are shown in Figures 8 and 7, respectively.

The best-fitting model parameters are $f = 1.8$, $r_s = 1.2 \times 10^5$ km, and $w_s = 10^3$ km, corresponding to $\chi_R^2 = 0.33$ and $P = 0.65$. The small value for $w_s$ implies that the density enhancement occurs abruptly, although larger values of $w_s$ (up to ~$10^4$ km) also produce radial CO profiles consistent with the observations, within errors. The distance over which the implied density enhancement occurs is therefore not well constrained by our data. This is primarily because the radius of the shell (projected in the plane of the sky) lies just beyond the spatial extent of our JCMT map.

For an outflow velocity of $0.5 \text{ km s}^{-1}$, a factor of approximately two density enhancement at $r_s = 1.2 \times 10^5$ km would be consistent with a corresponding drop in $Q(\text{CO})$ around 67 hr earlier (on January 12), but this seems perhaps unlikely given how uniform the comet's activity was over the weeks following that date (Figure 1). The individual HARP maps from January 14 and 15 are also consistent with each other (within the $2\sigma$ noise level), and do not show any evidence for outward-moving CO density structure(s) over this period, so an onset of icy grain sublimation at large radii is our favored explanation.

Accounting for the extended CO shell leads to a marginally significant reduction in the nucleus production rate retrieved from the JCMT map, from $Q_t(\text{CO}) = (6.1 \pm 0.8) \times 10^{28}$ s⁻¹ to $(5.4 \pm 1.1) \times 10^{28}$ s⁻¹. Adding the same extended CO component to our model for the time-averaged JCMT CO spectrum (Figure 6) does not significantly alter the quality of the spectral fit, but the resulting CO production rate from the nucleus ($(4.3 \pm 0.1) \times 10^{28}$ s⁻¹) is 17% lower.

### 4.5. SMA CO Visibility Analysis

Radio interferometry is a powerful technique for analyzing the radial distributions of cometary gases to provide insight into their physical and chemical origins in the coma (Boissier et al. 2007; Cordiner et al. 2014; Roth et al. 2021b). This is due to the interferometer's ability to simultaneously sample emission from a range of spatial scales at high accuracy, from the near-nucleus environment to the outer coma. Interferometric maps (such as the SMA CO map in Figure 4), however, suffer from significant artifacts due primarily to the sparsely filled telescope aperture, as well as Fourier image processing, regridding, and deconvolution artifacts. Cometary coma images suffer in particular from a lack of information on the largest spatial scales missed by the interferometer. Consequently, the preferred method for robustly analyzing interferometric data of such extended sources is by directly modeling the calibrated visibilities recorded by the telescope (i.e., the cross-correlation amplitudes between all antennas, as a function of baseline length).

Visibility models were generated for the SMA CO $J = 2-1$ observations of C/2016 R2, using our two-component SUBLIME model to test the same four scenarios as in Section 4.4. First, a fit was performed to the spectral line profile extracted at the emission peak of the SMA CO map (Figure 4). The best-fitting spectral line model is shown in Figure 10, corresponding to $R_Q = 9.2 \pm 1.5$, $v_1 = 0.64 \pm 0.02$ km s⁻¹, $v_2 = 0.20 \pm 0.02$ km s⁻¹, and $\theta = 74^\circ \pm 7^\circ$; the jet outflow axis was held fixed toward the Sun due to a lack of constraints (resulting from the lower spectral resolution of this data). The different component outflow velocities compared to those derived from the JCMT observations 3–5 weeks earlier (Section 4.3) imply an increase in the jet outflow velocity
<figure><img src="image_13.png" /><figcaption>**Figure 10.** CO *J* = 2 − 1 spectrum observed with the SMA (extracted at the CO emission peak), with the best-fitting two-component SUBLIME model overlaid using a red curve.</figcaption></figure>

accompanied by a slowing of the ambient coma. This could be
due to a stronger manifestation of intrinsic coma asymmetry at
the higher spatial resolution of the SMA data compared with
the JCMT. Changes in CO outflow velocity could also occur as
a result of temporal activity variations as the comet moved
closer to the Sun.

From this best-fitting base model, additional models were
constructed assuming (1) CO solely as a daughter species, with
Γ<sub>p</sub> = 7.3 × 10<sup>-5</sup> s<sup>-1</sup> derived from the fit to the JCMT HARP
data (Section 4.4); (2) CO as a parent with an additional outer-
coma source from CO<sub>2</sub> photolysis; and (3) CO as a parent with
an additional extended shell at r<sub>s</sub> = 1.2 × 10<sup>5</sup> km, with f = 1.8.
The resulting 3D model images were integrated in the spectral
domain and were then multiplied by the (FWHM = 55′′) SMA
primary beam pattern before sampling in the Fourier domain
using the vis_sample code (Loomis et al. 2018). Visibility
amplitude sampling was performed using the same set of uv
distances (baselines) as for the SMA observations (based on the
time-averaged antenna positions during the comet observa-
tions). A power-law curve (ax<sup>b</sup>) was fit through each set of
model results (with x as the baseline length, and a and b as free
parameters) and plotted along with the observed, time-averaged
visibility amplitudes in Figure 11. For clarity, the model
visibility curves were scaled vertically in order to pass through
the shortest-baseline point. The CO production rate of the best-
fitting (scaled) parent model is Q<sub>r</sub>(CO) = (6.7 ± 0.6) × 10<sup>28</sup> s<sup>-1</sup>
(the uncertainty increases to 0.9 × 10<sup>28</sup> s<sup>-1</sup> after the inclusion
of a 10% amplitude calibration error). This value is 1.6σ from
the mean CO production rate obtained using JCMT, corresp-
onding to a barely significant increase in Q<sub>r</sub> with time.

Three of the model curves (parent, parent + CO₂ photolysis, and parent + extended shell) all represent an equally good fit to the observations, falling precisely on top of each other in Figure 11. This implies that SMA would have been blind to the additional CO component observed at large radii in the JCMT maps—i.e., the extended CO component was smooth enough on large angular scales not to be detected on even the shortest SMA baselines. In contrast, the CO daughter model does not fit the observed visibilities well, with insufficient flux on small angular scales (long baselines). Consequently, the SMA observations rule out the possibility of CO being solely a

<figure><img src="image_14.png" /><figcaption>**Figure 11.** CO *J* = 2 − 1 visibility amplitudes vs. baseline length for C/2016 R2 observed using SMA, including 1σ statistical error bars. Model visibility curves are overlaid for four different CO distributions. The CO parent, parent + CO₂ photolysis, and parent + shell model curves all lie on top of each other. whereas the CO daughter curve differs significantly, underfitting the observations at large baseline (small angular scales).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Average (in velocity space) of the three strongest CH₃OH transitions in our JCMT observations (7₀ − 6₀ E, 7₋₁ − 6₋₁ E, and 7₀ − 6₀ A⁺). A tentative feature is present around 0 km s⁻¹ (the comet’s rest velocity); the blue shaded region shows the velocity range over which significant CO emission was detected.</figcaption></figure>

daughter species in this comet, although significant production
at large radii from CO₂ photolysis is still possible.

4.6. *Other Molecules from JCMT: CH<sub>3</sub>OH, HCN, H<sub>2</sub>CO, <sup>13</sup>CO, and HCO<sup>+</sup>*

Our JCMT CH3OH spectrum covered at least 12 lines of the
K = 7 - 6 band around 338 GHz (for details of the observed
transitions, see Cordiner et al. 2017, their Table 1), but no
individual CH3OH lines were clearly detected in our data. We
constructed a preliminary model for the CH3OH spectrum
based on the retrieved coma physical parameters from CO
(Section 4.3), and identified three transitions that clearly stood
out as stronger than the rest: J<sub>K</sub> = 7<sub>0</sub> - 6<sub>0</sub> E, 7<sub>-1</sub> - 6<sub>-1</sub> E, and
7<sub>0</sub> - 6<sub>0</sub> A<sup>+</sup>. The observed spectra for these three transitions
were then averaged together in velocity space, producing the
spectrum in Figure 12. A tentative feature is present around the
comet's rest velocity (0 km s<sup>-1</sup>), with a spectrally integrated
<figure><img src="image_16.png" /><figcaption>**Figure 13.** HCN, H₂CO, and ¹³CO spectra observed using JCMT, showing no evidence for detections of these molecules.</figcaption></figure>

intensity (in the range of -1.0 to 0.5 km s⁻¹ shown by the blue shaded region) of ∫T<sub>MB</sub>dv = 73 ± 20 mK km s⁻¹, corresponding to a 3.6σ detection. Fitting the same three transitions simultaneously using SUBLIME (allowing only the CH₃OH abundance to vary) gave Q(CH₃OH) = (6.7 ± 2.2) × 10²⁶ s⁻¹, which corresponds to a CH₃OH/CO abundance ratio of 1.3 ± 0.4% at the nucleus. For this model, we adopted CH₃OH–CO collisional transition rates based on the CH₃OH–H₂ rates from Rabli & Flower (2010), with solar pumping rates from Roth et al. (2021a). Uncertainties of a factor of five in the CH₃OH–CO collision rates lead to at most a 7% error on the CH₃OH abundance.

The HCN, H₂CO, and ¹³CO molecules were not detected in our data, as shown by the spectra in Figure 13. Upper limits of 3σ were derived by comparing the spectrally integrated noise level (from -0.8 to 0.5 km s⁻¹) to SUBLIME model line intensities for each molecule (again, using the same coma physical parameters derived in Section 4.3). For HCN we used the same collision and pumping rates as Cordiner et al. (2020), for H₂CO we used the Roth et al. (2021a) rates, and for ¹³CO we used the same rates as for CO. The production rate 3σ upper limits are Q(HCN) < 8.3 × 10²⁵ s⁻¹, Q(H₂CO) < 2.3 × 10²⁶ s⁻¹, and Q(¹³CO) < 2.7 × 10²⁷ s⁻¹, corresponding to abundance ratios

HCN/CO < 0.16%, H₂CO/CO < 0.45%, and CO/¹³CO > 19. The ¹²C/¹³C ratio in CO is therefore consistent with the typical solar system value of 89, as well as the value of 86 ± 9 measured in comet 67P (Altwegg et al. 2019).

We also obtained a nondetection of HCO⁺ J = 4 – 3, with ∫T<sub>MB</sub>dv < 164 mK km s⁻¹ (integrated over the velocity range ±1.5 km s⁻¹). This is about an order of magnitude less than the integrated HCO⁺ J = 3 – 2 line brightness observed in comet Hale–Bopp by Milam et al. (2004) using the Sub-Millimeter Telescope (SMT). However, the CO production rate in C/2016 R2 was also less than that in Hale–Bopp by at least an order of magnitude at the time of observation (Biver et al. 1997), so our result does not imply that C/2016 R2 had an unusually low HCO⁺ production rate (considering the importance of CO in coma HCO⁺ synthesis).

**5. Discussion**

**5.1. Comparison with Previous Observations of this Comet**

Our primary result is the determination of a revised CO
outgassing rate for C/2016 R2: Q(CO) = (5.3 ± 0.6) × 10²⁸
s⁻¹ (averaged over the JCMT observing period 2018 January
13 to 2018 February 1). This was subject to a marginally
significant rise to (6.7 ± 0.9) × 10²⁸ s⁻¹ on 2018 February 21,
observed using SMA. These results confirm C/2016 R2 as
having among the highest CO production rates ever observed in
a comet—only a factor of four less than C/1995 O1 (Hale-
Bopp) at similar heliocentric distances. Our CO production
rates are slightly higher than the average value of
(4.6 ± 0.4) × 10²⁸ s⁻¹ between 2017 December 22 and 2018
January 16 observed by Wierzchos & Womack (2018) using
SMT, and similar to the (5.5 ± 0.9) × 10²⁸ s⁻¹ observed on
2018 February 23 by McKay et al. (2019), consistent with a
slow, steady rise in CO activity as the comet moved closer to
the Sun (between r<sub>H</sub> = 2.98–2.73 au). The observations by
Biver et al. (2018) using the Institute for Radio Astronomy in
the Millimeter Range (IRAM) 30 m telescope on 2018 January
24 gave Q(CO) = (10.6 ± 0.5) × 10²⁸ s⁻¹, which is approxi-
mately double our JCMT value. As shown in Appendix C, the
1D versions of ours and the models of Biver et al. produce
near-identical results. The discrepancy in CO production rates
is therefore primarily attributable to the unusual complexity of
the CO spectral line profile combined with differences between
our 3D radiative transfer modeling strategies.

To derive CO production rates, Biver et al. (2018) adopted a
spherical coma model divided into three regions of colatitude
(γ; angle measured from the Earth–comet vector), with outflow
velocities v = 0.56 km s⁻¹ between γ = 0–60°, 0.50 km s⁻¹
between γ = 60°–120°, and no outflow for γ > 120. Our
model, on the other hand, gives v₁ = 0.51 km s⁻¹ in a sunward-
facing conical jet and v₂ = 0.25 km s⁻¹ in the remaining
(ambient) coma. Our chosen model geometry is physically
justified based on fluid dynamic and Monte Carlo coma models
(e.g., Crifo et al. 1999; Fougere et al. 2016), which show
significantly different outflow velocities on the sunward and
antisunward sides of the nucleus. The large reduction in the
outflow velocity of our model across the entirety of the
nightside (most of which overlaps with the γ > 120° region of
no outflow in the model of Biver et al. 2018) necessitates a
lower total production rate in our model compared with that of
Biver et al. (2018).
Although our two-component outflow model is sufficient to fit the data in the present study, we concede that it likely represents a simplification of the true physical situation. Continuously variable Q and v as a function of (3D) coma position may be more correct, and could account for smoothly varying temperatures and mixing ratios within the nucleus, while avoiding discontinuities in the model parameters at the (conical) jet boundaries. The reality for C/2016 R2 may therefore involve a continuously variable Q (as in the more complex model presented in Figure 3 of Biver et al. 2018), but with v also allowed to vary more substantially/continuously, as may be expected based on physical models (e.g., Tenishev et al. 2008). Further constraints on the 3D outflow velocity distribution may be obtained from theoretical models, to help break the degeneracy that can occur between Q and v (particularly in the parts of the coma around $\gamma \sim 90^{\circ}$), in order to obtain the most accurate production rates for an asymmetric coma.

Compared to a completely spherically symmetric outflow model (as used by Wierzchos & Womack 2018 and McKay et al. 2019), our best-fitting two-component model for C/2016 R2 leads to a ~30% lower CO production rate, as a result of reduced outgassing on the antisunward side (see Figure 6). The exact correction factors are, however, dependent on the opening angle of the jet, and are larger if the jet is narrower. Based on the range of opening angles ($\theta = 32-77^{\circ}$) obtained on our different JCMT observing dates (Table 2), correction factors in the 29%-48% range are possible, so the CO production rate reported by McKay et al. (2019) should be adjusted to (2.9–3.9) $\times 10^{28}$ s$^{-1}$. Assuming the comet’s H$_2$O outgassing behavior was similar to that of CO and CH$_3$OH (see Biver et al. 2018, their Figure 14), the McKay et al. (2019) H$_2$O production rate may also need to be adjusted. However, the correction factor is much lower for OH (the H$_2$O daughter fragment observed by McKay et al. 2019), due to the OH kinetic energy gained in the photolysis reaction H$_2$O + hν → OH + H, which causes asymmetries in the gas velocity distribution to be smoothed out.

Abundance ratios from a given study tend to be less susceptible to model-dependent and instrumental uncertainties. Our CH$_3$OH/CO abundance ratio of 1.3% ± 0.4% matches the value of Biver et al. (2018) of 1.04% ± 0.08%, and our HCN/CO and H$_2$CO/CO upper limits are also consistent with their values of (3.8 ± 1.0) $\times 10^{-3}$% and 0.043% ± 0.006%, respectively. The CH$_3$OH/CO ratio is significantly higher than the upper limit of McKay et al. (2019) of 0.38%, which could be indicative of significant temporal variability in the CH$_3$OH outgassing rate during 2018 January. Alternatively, the higher CH$_3$OH abundance observed using radio spectroscopy could have been due to additional CH$_3$OH production in the extended coma (from icy grain sublimation; e.g., Coulson et al. 2017), which was not detected on the smaller angular scales probed by the infrared observations of McKay et al. (2019).

## 5.2. Molecular Abundances and the True Nature of C/2016 R2 (PanSTARRS)

Using the H$_2$O production rate from McKay et al. (2019), our JCMT results are consistent with a CO/H$_2$O ratio ~170, which is ~2800 times greater than the average value observed for Oort cloud comets (Dello Russo et al. 2016), and ~37 times greater than the highest previously observed in a comet (29P; Ootsubo et al. 2012), thus confirming the extremely CO-rich

nature of C/2016 R2’s coma. The CO/H$_2$O abundance in cometary comae is observed to vary strongly as a function of heliocentric distance (see Wierzchos & Womack 2018; McKay et al. 2019) due to the very different sublimation temperatures of these gases ($T_{sub}(CO) = 24$ K versus $T_{sub}(H_2O) = 152$ K), so a reduction in H$_2$O outgassing is expected for comets at heliocentric distances $\ge 3$ au, where the ice temperature falls below $T_{sub}(H_2O)$. The CO/H$_2$O ratio in the coma of C/2016 R2 is several hundred to several thousand times greater than that found in other Oort cloud comets at similar heliocentric distances (Crovisier et al. 1997; Ootsubo et al. 2012; Kawakita et al. 2014), so this comet appears anomalous compared with all those previously observed. It should be kept in mind, however, that the coma abundances are not necessarily representative of those in the nucleus ices.

Cometary nuclei are heterogeneous, containing mixtures of ices in different phases and compositions (A’Hearn et al. 2011; Mumma & Charnley 2011; Altwegg et al. 2019). Chemical differentiation as a function distance below the surface, or the presence of a volatile-depleted outer crust (e.g., Capaccioni et al. 2015), could lead to H$_2$O being insulated from solar heating while CO continues to sublimate. For example, the 3D numerical nucleus model of Marboeuf & Schmitt (2014) demonstrates that the ratio of CO-to-H$_2$O production rates could be enhanced by several orders of magnitude by the presence of an insulating dust mantle ~5–10 cm thick that hinders H$_2$O sublimation. Alternatively, a moderate overabundance of CO ice close to the surface could lead to increased cooling by CO sublimation, helping keep the nucleus at a low enough temperature to inhibit H$_2$O sublimation (Lisse et al. 2021). In that case the comet could maintain a lower H$_2$O outgassing rate (relative to CO) for a longer duration as it approached the Sun. Either case may not require an extremely anomalous CO/H$_2$O abundance in the bulk nucleus.

Peculiarities in the coma abundance ratios for several other molecules were reported by Biver et al. (2018) and McKay et al. (2019). Despite similar sublimation temperatures for CH$_3$OH and HCN (99 K and 95 K, respectively), these two molecules were enriched in the coma (relative to H$_2$O) by very different amounts: the CH$_3$OH/H$_2$O ratio was 163 times the Oort cloud comet average, whereas HCN/H$_2$O was only 5.9 times the average. Evidently, such enrichment patterns cannot be produced by reduced H$_2$O outgassing alone, or by a simple temperature dependence of the CH$_3$OH and HCN sublimation rates. Similarly, the coma CH$_4$/H$_2$O ratio is a factor of 206 higher than average, but this is several hundred times less than the CO/H$_2$O enrichment, despite similar sublimation temperatures for these two molecules (31 K and 24 K, respectively). It is tempting to take such unusual coma abundance patterns to be directly representative of a peculiar nucleus ice composition, but before doing so, it is worth considering the possible role played by ice heterogeneity and molecular trapping at temperatures below $T_{sub}(H_2O)$ (i.e., with the comet not yet fully activated).

As shown by laboratory ice sublimation experiments (Collings et al. 2004), volatile gases can be trapped in mixed (H$_2$O-dominated) ices at temperatures well above their sublimation points, and may be only partially released until the H$_2$O sublimation temperature is reached (~150 K). The release of trapped hypervolatiles can also occur as H$_2$O ice undergoes a phase change from amorphous to crystalline (at $T_{AC} \sim 130$ K). Such trapping and release processes are inevitable in mixed
cometary ices, and could explain some of the observed abundance patterns in C/2016 R2. Closely related outgassing behaviors were observed for CO, CO₂, and CH₃OH in comet 67P (Biver et al. 2019; Läuter et al. 2019), consistent with our understanding regarding a common origin for these species in carbon- and oxygen-rich interstellar ices (Fuchs et al. 2009; Garrod & Pauly 2011; Ioppolo et al. 2011). Observations of young stellar objects also indicate mixing of CO, CO₂, and CH₃OH ices in an apolar phase, distinct from the polar, H₂O-dominated ice (Boogert et al. 2015; Penteado et al. 2015). In situ observations of comets 103P (A'Hearn et al. 2011) and 67P (e.g., Migliorini et al. 2016; Gasc et al. 2017) discovered spatial and temporal variations in the coma H₂O/CO₂ ratios, implying nonuniform mixing ratios for these two volatiles within the nucleus. As shown by Davidsson et al. (2021, 2022) for 67P, this could be caused by physical/chemical evolution of the comet's surface layers due to anisotropic illumination of the comet, but primordial variations in the H₂O/CO₂ ratio intrinsic to the nucleus are also possible. It is therefore plausible that CH₃OH in C/2016 R2 exists as a component within an apolar ice matrix dominated by CO and CO₂, and some of this CH₃OH is released into the coma when CO sublimates (as observed). Meanwhile, a significant CH₃OH component also likely remains frozen as part of the comet's polar, H₂O-rich ices (Qasim et al. 2018). HCN, on the other hand, could be primarily associated with the (still frozen) H₂O ice rather than the (sublimating) CO component, and would then be only partially outgassed at *T*<sub>sub</sub>(HCN) or *T*<sub>AC</sub>. Future studies of cometary CO, CH₃OH, and HCN spatial distributions could help test this hypothesis.

To complete this picture, the C₂H₆ and NH₃ upper limits from McKay et al. (2019) are consistent with moderate-to-no enrichment (with respect to H₂O) and could therefore be associated primarily with the H₂O-dominated ice, while the strongly enriched N₂ and CO₂ are associated more with CO. Observations and laboratory studies show that interstellar CH₄ ice tends to be more associated with H₂O than CO (Öberg et al. 2008; Qasim et al. 2020), while the location of H₂CO ice is less well constrained. The CH₄ and H₂CO enrichment factors in C/2016 R2 (206 and 44, respectively) could thus be explained by the release of trapped volatiles (in H₂O ice) above their respective sublimation temperatures (31 and 64 K; Collings et al. 2004). The increasing similarity to typical abundances (with respect to H₂O) in the sequence CH₄ → H₂CO → HCN is consistent with the decreasing volatility (increasing *T*<sub>sub</sub>) of these three species, such that they each behave progressively more like H₂O. We therefore postulate the existence of a rapidly sublimating (apolar) component of ice in C/2016 R2 rich in CO, CO₂, N₂, and CH₃OH, and a second (polar) component containing more CH₄, H₂CO, and HCN mixed with H₂O. If the N₂/CO ratio is higher in the CO-rich ice phase than in the H₂O-rich ice, this could explain why the N₂/CO ratio is lower in fully activated comets than in C/2016 R2.

## 5.3. Coma Morphology

The highly asymmetric, blueshifted CO line profile of C/2016 R2 is similar to that of the large Centaur 29P/SW 1 (Festou et al. 2001; Gunnarsson et al. 2002), as well as to the CO line profile observed in C/1995 O1 at $r_H \gtrsim 8$ au (Gunnarsson et al. 2003). Our interpretation of the line shape in terms of enhanced CO production and outflow velocity on

the sunward side of the nucleus is consistent with the analysis of the 29P coma by Gunnarsson et al. (2008), and we find a similar ratio of day-to-night hemisphere *Q*(CO) and $v_{\text{out}}$ values in C/2016 R2.

Gunnarsson et al. (2002) also discovered an extended shell of CO emission surrounding 29P, at a cometocentric distance of $r_c \sim 1.4 \times 10^5$ km, based on mapping observations of the J = 2 − 1 line in 1998 (although the shell was no longer apparent in 2003 follow-up observations, demonstrating an intermittent nature; Gunnarsson et al. 2008). The spatial properties of this shell are remarkably similar to those of the extended CO emission structure we detected at a similar cometocentric distance in our C/2016 R2 JCMT maps (Section 4.4). Gunnarsson (2003) interpreted the extended CO structure in 29P as arising from sublimation of a population of icy grains long-lived enough to reach 29P's outer coma. Such an explanation also appears plausible for C/2016 R2, although more detailed modeling would be required to confirm this possibility, and to test the other possible origins for the shell considered in Section 4.4. We do not consider CO₂ to be a likely source for the majority of the extended CO in C/2016 R2 as it would require at least an order of magnitude larger CO₂ production rate than that found by McKay et al. (2019).

A time-variable, diffuse, ring-shaped feature (consistent with excess gas emission) also appeared in the Spitzer IRAC CO + CO₂ images of C/2016 R2 on 2018 February 21—a feature that was not apparent 9 days earlier (see Figure 5 of McKay et al. 2019). The approximate diameter of this ring was also $\sim 1.4 \times 10^5$ km (M. Kelley, 2021, private communication), so it could plausibly be related to the shell-like feature observed using JCMT. Detailed modeling of the Spitzer images will be required to determine the physical properties of this ring, and to confirm whether icy grain sublimation or coma deceleration could be responsible.

On both dates, the Spitzer images also show enhanced gas emission within an angular wedge (of opening angle $\sim 80^\circ$), oriented toward the (sky-projected) comet–Sun vector, with a morphology consistent with a jet or fan emanating from the nucleus. Although this feature may be attributable to a combination of both CO and CO₂ emission, its apparent qualitative similarity to the jet properties derived from modeling our JCMT and SMA data provides further evidence for preferential sunward outgassing from a confined region of the nucleus. This provides additional validation of our two-component coma model.

## 5.4. Uncertainties in the CO Excitation

Our new CO–CO collisional rate coefficients are considered to be accurate to within about a factor of two (Section 4.1). As shown in Figure 15 (Appendix C), differences in the rates on that order lead only to relatively small discrepancies in the CO rotational level populations, the evolution of which is controlled to a large extent by radiative processes, especially in the outer coma. As a result, varying the collision rates by $\pm$ a factor of two leads to changes of only $\pm 1\%$ in our retrieved CO production rates.

An additional source of error in the CO excitation calculation arises from the fact that the collision rates are tabulated as a function of $T_{\text{rot}}$, which provides only an average measure of the true distribution of energy-level populations. In a non-LTE regime, individual-level populations may deviate from the Boltzmann distribution, in which case the (thermally averaged)
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Modeled CO rotational temperature derived using the seven lowest energy levels, from our time-dependent SUBLIME model (solid blue curve) and the steady-state version of the model (dotted gray curve), based on the data shown in Figure 9. The rotational temperature evolves from thermal equilibrium at the gas kinetic temperature (18.7 K) in the collisionally dominated inner coma, to fluorescence equilibrium (at 75 K) in the solar radiation-dominated outer coma. The vertical dashed line shows the radial extent of the 14" JCMT beam FWHM at the geocentric distance of C/2016 R2 (2.14 au), at 346 GHz.</figcaption></figure>

state-to-state collision rate coefficients are no longer accurate. Perhaps more importantly, the rate coefficients are dependent on $T_{rot}$ as well as $T_{kin}$, yet for the purpose of choosing rates in our model, we assume $T_{rot} = T_{kin}$, which is only strictly accurate in the collisionally dominated inner coma. Both these issues have negligible impact on the results of our study, however, since the collision rates $k_{ij}$ are found not to vary strongly as a function of $T_{rot}$, changing by on average only 9% between 10–30 K.

For deriving accurate CO abundances, rather than having the most precise set of collision rate coefficients, it is more important to use the most physically accurate, time-dependent solution to the equation of statistical equilibrium (Bockeleer-Morvan 1987; Biver 1997), as opposed to the steady-state approximation (Hogerheijde et al. 2009; de Val-Borro et al. 2018). Figure 9 shows a comparison between the CO level populations as a function of radius using these two different methods, while Figure 14 shows the corresponding rotational temperature behavior. As described by Garcia-Berrios et al. (2020), species with small dipole moments such as CO have long timescales $\tau_R$ with respect to spontaneous (and induced) rovibrational radiative transitions. If $\tau_R$ is greater than the dynamical timescale in the outflow (i.e., the time for the coma gas to move a given distance $r_d$), then an accurate calculation of the excitation over distances $r_d$ requires the gas motion to be considered, and the steady-state approximation is no longer valid.

An observational consequence for molecules with large $\tau_R$ is that the non-LTE effects occur further out in the coma, because the time-dependent solution (solid curves in Figures 9 and 14) is delayed with respect to the steady-state solution (dashed curves). For the $J=3$ level populations, the discrepancy is minimal (amounting to <1% difference in the $J=3-2$ line integrated intensity). However, for the $J=2-1$ line, the steady-state solution overestimates the line intensity by 19%. For larger beam sizes more sensitive to the non-LTE region of

the coma (between the thermal and fluorescence equilibrium extremes highlighted in Figure 14), the discrepancy can be even larger. Consequently, the use of a time-dependent excitation model is recommended for correct analysis of single-dish cometary CO data. It should be borne in mind that this issue is less severe for molecules with larger dipole moments (and therefore, smaller $\tau_R$ values) such as H$_2$O, HCN, CH$_3$OH, and H$_2$CO, although in low-activity comets the premature departure from LTE of the steady-state solution can still lead to some large discrepancies. For example, in a comet at $r_H = \Delta = 1$ au, with $Q(H_2O) = 10^{27}$ s$^{-1}$, $T_{kin} = 50$ K, and $v_{out} = 0.8$ km s$^{-1}$, the HCN $J=4-3$ line intensity (for a 14" JCMT beam) is overestimated in the steady-state model by 46%, whereas at $Q(H_2O) = 10^{29}$ s$^{-1}$ the discrepancy is reduced to only 4%.

## 6. Conclusions

Based on high-resolution spectral-spatial observations using the JCMT and SMA telescopes during the period 2018 January 13 to 2018 February 21, we confirm the presence of extremely strong, asymmetric CO outgassing from the hypervolatile-rich comet C/2016 R2 (PanSTARRS), with $Q(CO)$ in the range of $(3.8-7.6) \times 10^{28}$ s$^{-1}$. The observational data were analyzed using a new, time-dependent, 3D radiative transfer code, adopting a two-component model for the expanding coma, for the first time using state-to-state CO-CO collision rates based on quantum scattering calculations using the coupled-states method. We determined the presence of a (near-)sunward CO jet with a (time-variable) half-opening angle of $\theta = 25^\circ-90^\circ$ (and an average of $\theta = 62^\circ \pm 2^\circ$, offset by $7^\circ \pm 2^\circ$ from the Sun-comet vector). The average jet outflow velocity determined from our JCMT data was $0.51 \pm 0.01$ km s$^{-1}$, while the ambient coma outflow velocity was found to be $0.25 \pm 0.01$ km s$^{-1}$. The total amount of CO produced by the jet was, on average, a factor of two more than that of the ambient coma.

On 2018 January 14–15, we found evidence for extended CO emission that cannot be easily explained by standard nucleus outgassing or excitation effects. The extended emission is therefore interpreted as a possible result of modulation in the CO outgassing rate, deceleration in the outer coma, or sublimation of long-lived icy grains. Subtraction of such an extended CO component from the best-fitting (time-averaged) JCMT 3–2 model results in a 17% reduction in the CO production rate derived for the nucleus. Using the $Q(H_2O)$ value from McKay et al. (2019), our CO/H$_2$O ratio is ~37 times larger than seen in any comet to date (including the distant Centaur 29P), which, combined with previously noted chemical peculiarities, suggests that C/2016 R2 is among the most unusual comets ever observed.

However, the heterogeneous nature of cometary ices, combined with knowledge that molecular outgassing rates from mixed ices do not necessarily correlate with their sublimation temperatures, means that we cannot yet rule out a bulk composition for C/2016 R2 more similar to the general population of Oort cloud comets than previously inferred. We hypothesize that the ice temperature of C/2016 R2 may have been suppressed by sublimative cooling, or the presence of an unusually thick insulating crust, which prevented the initiation of a more conventional, H$_2$O-dominated outgassing regime, leading to sublimation rates more heavily influenced by trapping and binding of individual molecules within the bulk
ices. We propose that the observed abundance patterns can be explained by the existence of two ice phases in C/2016 R2, similar to those observed in young stellar objects: (1) an apolar phase, rich in CO, CO₂, N₂, as well as CH₃OH ices; and (2) a polar phase containing larger abundances of CH₄, H₂CO, and HCN mixed in with H₂O ice. More observations of coma chemistry in distant comets (at $r_H \gtrsim 2.5$ au, for which H₂O sublimation is not yet fully activated) will be crucial to better understand this comet's peculiar nature, and to constrain the physical and chemical processes that govern the formation, storage, and release of cometary volatiles below the H₂O sublimation point.

This work was supported by the National Science Foundation under grant No. AST-2009253. The work of S.N.M., N.X. R., M.A.C., E.G.B., and S.B.C. was also supported by NASA's Planetary Science Division Internal Scientist Funding Program through the Fundamental Laboratory Research work package (FLaRe). S.N.M. and S.B.C. were supported by the NASA Astrobiology Institute through the Goddard Center for Astrobiology. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China, and the Chinese Academy of Sciences (grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. We gratefully acknowledge the assistance of R. Loomis and S. Andrews with the use of the vis_sample code. We thank N. Biver for providing comparison radiative transfer model calculations. Thanks also to A. McKay for providing comments on the manuscript regarding the nucleus and coma composition.

Software: CASA (v5.1; Jaeger 2008), SMA-MIR (https://github.com/qi-molecules/sma-mir), LIME code (Brinch & Hogerheijde 2010), CVODE solver (Hindmarsh 2019), Mosaic (Hutson & Green 1994), vis_sample (Loomis et al. 2018), MPFIT (Markwardt 2012).

## Appendix A
Equations of Radiative Transfer and Molecular Excitation

The intensity of radiation ($I_ν$) propagating through the coma at a frequency $ν$ is calculated by integrating the equation of radiative transfer as a function of distance along the line of sight (s),

$$ \frac{dI_ν}{ds} = j_ν - \alpha_ν I_ν, \qquad (A1) $$

where $j_ν$ and $\alpha_ν$ are the gas emission and absorption coefficients, respectively. These are derived from the Einstein coefficients of the gas ($A_{ij}$, $B_{ij}$, and $B_{ji}$, for a transition between the upper energy level i and lower level j), as

$$ j_ν = \frac{hν}{4π} N_j A_{ij} ψ_ν, \qquad (A2) $$

$$ \alpha_ν = \frac{hν}{4π} (N_j B_{ij} - N_i B_{ji}) ψ_ν, \qquad (A3) $$

where $N_i$, $N_j$ are the number of gas particles per unit volume in levels i and j, respectively, and $ψ_ν$ is a (normalized) line-broadening function for the spectral line of interest (typically a Gaussian for individual, thermally broadened lines).

The number of molecules (per unit volume) in energy level i is obtained as a function of time (t) in the outflowing coma gas by solving the following differential equation (e.g., Crovisier 1987):

$$ \frac{dN_i}{dt} = -N_i \left[ \sum_{j<i} A_{ij} + \sum_{j \neq i} (B_{ij}J_ν + k_{ij}n + G_{ij}) \right] + \sum_{j>i} N_j A_{ji} + \sum_{j \neq i} N_j (B_{ji}J_ν + k_{ji}n + G_{ji}). \qquad (A4) $$

In this equation, $k_{ij}$ are the rates at which transitions occur between rotational levels i and j (in the ground vibrational state) due to collisions between the gas particles, n is the number density of colliders (in this case, the CO gas density), and $G_{ij}$ are the effective transition rates due to fluorescence/vibrational pumping by solar radiation, summed over the relevant rovibrational bands (see, e.g., Crovisier & Encrenaz 1983; Bensch & Bergin 2004). The energy-level populations at a given location in the coma depend on the local radiation field, $J_ν$, which is calculated by summing the incident radiant energy received at that point from all solid angles. In general, this means that Equation (A4) needs to be solved iteratively until convergence of $J_ν$ is achieved. In practice however, for species other than H₂O the optical depth for photons leaving the less dense parts of the coma where non-LTE effects are important tends to be low (i.e., $τ_ν \ll 1$). In that case, the stimulated emission and absorption terms ($B_{ji}J_ν$) are small, and can be neglected. For more optically thick gases, the escape probability method can be used as a quick (and easy to implement) approximation for the effects of photon trapping (e.g., Bockelee-Morvan 1987).

Solar radiation-induced fluorescence (pumping) is responsible for modifying the rotational level populations. Effective pumping rates ($G_{ij}$) for CO were calculated using the method of Crovisier & Encrenaz (1983), incorporating the latest infrared transition data from the HITRAN catalog (Gordon et al. 2021). The effective pumping rates were summed over all rovibrational transitions involving the ground vibrational state of CO. Excitation of the gases of interest due to collisions with coma electrons has also been implemented in SUBLIME using the method of Biver (1997) and Zakharov et al. (2007), but for the present study focusing on CO emission from a CO-dominated coma, the electron-collision rates are found to be small enough that they can be neglected.

## Appendix B
State-to-state Collision Rates for the CO-CO System

Table 3 shows collisional (de-excitation) rate coefficients $k_{J_1 \to J_1'} = k_{ij}$ (in cm⁻³ s⁻¹) for gas-phase CO molecules undergoing transitions $J_1 \to J_2$, as a function of kinetic temperature ($T_{kin}$). The reverse (excitation) rates are calculated from the principle of detailed balance according to Equation (6) of van der Tak et al. (2007). For further details of the rate coefficient calculations, see Section 4.1.
<table><caption>**Table 3**
CO–CO State-to-state Collision Rate Coefficients $k_{J_1 \to J_1'}$ as a Function of $T_{\text{kin}}$</caption><thead><tr><th rowspan="2">J<sub>1</sub></th><th rowspan="2">J<sub>1</sub>'</th><th colspan="4">T<sub>kin</sub></th></tr><tr><th>5 K</th><th>10 K</th><th>20 K</th><th>30 K</th></tr></thead><tbody><tr><td>1</td><td>0</td><td>2.69e-11</td><td>3.61e-11</td><td>4.01e-11</td><td>4.13e-11</td></tr><tr><td>2</td><td>0</td><td>2.73e-11</td><td>3.03e-11</td><td>2.97e-11</td><td>2.81e-11</td></tr><tr><td>2</td><td>1</td><td>5.56e-11</td><td>6.23e-11</td><td>6.69e-11</td><td>6.70e-11</td></tr><tr><td>3</td><td>0</td><td>1.81e-11</td><td>2.02e-11</td><td>1.96e-11</td><td>1.81e-11</td></tr><tr><td>3</td><td>1</td><td>4.73e-11</td><td>5.14e-11</td><td>5.40e-11</td><td>5.30e-11</td></tr><tr><td>3</td><td>2</td><td>8.19e-11</td><td>7.88e-11</td><td>7.67e-11</td><td>7.50e-11</td></tr><tr><td>4</td><td>0</td><td>1.42e-11</td><td>1.49e-11</td><td>1.43e-11</td><td>1.30e-11</td></tr><tr><td>4</td><td>1</td><td>3.55e-11</td><td>3.74e-11</td><td>3.73e-11</td><td>3.49e-11</td></tr><tr><td>4</td><td>2</td><td>6.01e-11</td><td>6.22e-11</td><td>6.28e-11</td><td>6.10e-11</td></tr><tr><td>4</td><td>3</td><td>5.78e-11</td><td>6.66e-11</td><td>7.01e-11</td><td>6.98e-11</td></tr><tr><td>5</td><td>0</td><td>9.94e-12</td><td>9.69e-12</td><td>8.43e-12</td><td>7.43e-12</td></tr><tr><td>5</td><td>1</td><td>2.82e-11</td><td>2.95e-11</td><td>2.77e-11</td><td>2.54e-11</td></tr><tr><td>5</td><td>2</td><td>4.17e-11</td><td>4.30e-11</td><td>4.13e-11</td><td>3.84e-11</td></tr><tr><td>5</td><td>3</td><td>5.98e-11</td><td>6.46e-11</td><td>6.36e-11</td><td>6.06e-11</td></tr><tr><td>5</td><td>4</td><td>5.59e-11</td><td>6.40e-11</td><td>6.74e-11</td><td>6.64e-11</td></tr></tbody></table>

## Appendix C
Benchmarking the SUBLIME CO Model

To confirm the accuracy of our coma radiative transfer and excitation model, we compared results with the similar, well-tested model of N. Biver (2021, private communication). Their model has been used to analyze millimeter/submillimeter rotational spectra of numerous comets over the last few decades (e.g., Biver et al. 1999; Bockelée-Morvan et al. 2012; Biver et al. 2016, 2019), and is based on the formalism presented by Crovisier (1987) and Bockelee-Morvan (1987), described in more detail by Biver (1997) and Bockelée-Morvan et al. (2004a). Figure 15 shows a comparison between the computed CO energy-level populations from their model and our SUBLIME model, adopting a spherically symmetric (1D) outflow geometry with $Q(\text{CO}) = 5 \times 10^{28} \text{ s}^{-1}$, a constant kinetic temperature of $T_{\text{kin}} = 20 \text{ K}$, and an outflow velocity of $v_{\text{out}} = 0.5 \text{ km s}^{-1}$, at a heliocentric distance of $r_H = 2.8 \text{ au}$. The results of the Biver et al. (2021, private communication) model are shown with dashed lines, whereas our model is shown with solid lines.

The model populations are in very close agreement considering the complexity of the calculation and the different assumptions regarding the relevant molecular parameters. Outside of the collisionally dominated (LTE) zone, the models diverge slightly, primarily as a result of the different treatments of CO–CO collision rates—we are using quantum-mechanically derived state-specific rate coefficients, whereas their model uses thermalizing rate coefficients based on an (assumed) uniform collisional cross-section. Slight differences are also evident at fluorescence equilibrium (largest nucleocentric distances, where the populations reach a steady state with respect to the solar radiation field), presumably due to small differences in the rovibrational Einstein A coefficients used to derive the pumping rates (our model uses the latest HITRAN data).

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Fractional CO rotational energy-level populations for $J = 0-5$ as a function of radius, based on a non-LTE, spherically symmetric coma model with $Q(\text{CO}) = 5 \times 10^{28} \text{ s}^{-1}$, $T_{\text{kin}} = 20 \text{ K}$, $v_{\text{out}} = 0.5 \text{ km s}^{-1}$, and $r_H = 2.8 \text{ au}$. Solid curves are using the time-dependent version of SUBLIME, with CO–CO collision rates from Section 4.1; dotted curves are assuming that the CO–CO collision rates are the same as the CO–H$_2$ rates from Yang et al. (2010); and dashed curves are using the model of Biver et al. (2018).</figcaption></figure>

As a test of the SUBLIME ray-tracing algorithm, we compared integrated line fluxes from our output model spectral images for two different CO lines, convolved to the JCMT spatial resolution. For the $J=2-1$ line, we obtained $\int T_R dv = 0.37 \text{ K km s}^{-1}$, compared with $0.36 \text{ K km s}^{-1}$ from the model of Biver et al. (2019), and for $J=3-2$ we have $\int T_R dv = 0.64 \text{ K km s}^{-1}$, compared with $0.63 \text{ K km s}^{-1}$, again, demonstrating very good agreement between our models, at a level much less than the observational uncertainties.

Figure 15 also shows (with a dotted line) the level populations derived from the same SUBLIME model, but using CO–H$_2$ rates from Yang et al. (2010) to approximate the CO–CO collision rates, instead of our new state-specific CO–CO rates calculated in Section 4.1 (for $J<6$). At 20 K, the new rate coefficients differ by up to a factor of 7.1 from those of CO–H$_2$, with a mean ratio between the new and H$_2$-derived rates of 2.3. The resulting discrepancy in the final results between coma models using the new rate coefficients as opposed to adopting CO–H$_2$ rates, however, is relatively small.

## Appendix D
Azimuthally Averaged JCMT HARP Spectra

Azimuthal averages (about the central pixel) of the JCMT HARP CO 3–2 spectral-spatial data cube are shown in Figure 16. These are based on the average of the HARP jiggle-map observations from 2018 January 14 and 2018 January 15. Each spectrum has been scaled (normalized to the same peak value as the spectrum from the central pixel) to cancel out the rapid decay in the overall line intensity with radius due to the falling coma density. Within the noise, there is no obvious evolution in the spectral line profile with distance from the comet.
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Azimuthally averaged CO 3–2 spectra, based on JCMT HARP jiggle maps (combined data from 2018 January 14 and 2018 January 15). These are the mean spectra within successive ring-shaped annuli (1 pixel thick), centered on the nucleus, and show the spectral profile as a function of radial offset from the comet (from top to bottom). The numbers to the right of each spectrum indicate the radius of each annulus in units of $10^4$ km. The spectra have each been scaled for display purposes (normalized to the same peak value as the central, 0.0 km spectrum) with scale factors denoted by the numbers in parentheses.</figcaption></figure>

## ORCID iDs

M. A. Cordiner https://orcid.org/0000-0001-8233-2436

I. M. Coulson https://orcid.org/0000-0002-7316-4626

E. Garcia-Berrios https://orcid.org/0000-0002-1069-2931

C. Qi https://orcid.org/0000-0001-8642-1786

F. Lique https://orcid.org/0000-0002-0664-2536

M. de Val-Borro https://orcid.org/0000-0002-0455-9384

Y.-J. Kuan https://orcid.org/0000-0002-4336-0730

S. Mairs https://orcid.org/0000-0002-6956-0730

N. X. Roth https://orcid.org/0000-0002-6006-9574

S. B. Charnley https://orcid.org/0000-0001-6752-5109

S. N. Milam https://orcid.org/0000-0001-7694-4129

W.-L. Tseng https://orcid.org/0000-0001-7322-3801

Y.-L. Chuang https://orcid.org/0000-0001-5653-3584

## References

A'Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, *Sci*, 332, 1396

Altwegg, K., Balsiger, H., & Fuselier, S. A. 2019, *ARA&A*, 57, 113

Bensch, F., & Bergin, E. A. 2004, *ApJ*, 615, 531

Biver, N. 1997, PhD Thesis, Univ. Paris 7-Didero

Biver, N., Bockelée-Morvan, D., Colom, P., et al. 1997, *Sci*, 275, 1915

Biver, N., Bockelée-Morvan, D., Crovisier, J., et al. 1999, *AJ*, 118, 1850

Biver, N., Bockelée-Morvan, D., Hofstadter, M., et al. 2019, *A&A*, 630, A19

Biver, N., Bockelée-Morvan, D., Paubert, G., et al. 2018, *A&A*, 619, A127

Biver, N., Moreno, R., Bockelée-Morvan, D., et al. 2016, *A&A*, 589, A78

Bockelee-Morvan, D. 1987, *A&A*, 181, 169

Bockelée-Morvan, D., Biver, N., Colom, P., et al. 2004b, *Icar*, 167, 113

Bockelée-Morvan, D., Biver, N., Swinyard, B., et al. 2012, *A&A*, 544, L15

Bockelée-Morvan, D., Crovisier, J., Mumma, M. J., & Weaver, H. A. 2004a, in *Comets II*, ed. G. W. Kronk (Tucson, AZ: Univ. Arizona Press), 391

Bockelée-Morvan, D., Hartogh, P., Crovisier, J., et al. 2010, *A&A*, 518, L149

Bøgelund, E. G., & Hogerheijde, M. R. 2017, *A&A*, 604, A131

Boissier, J., Bockelée-Morvan, D., Biver, N., et al. 2007, *A&A*, 475, 1131

Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, *ARA&A*, 53, 541

Brinch, C., & Hogerheijde, M. R. 2010, *A&A*, 523, A25

Buckle, J. V., Hills, R. E., Smith, H., et al. 2009, *MNRAS*, 399, 1026

Capaccioni, F., Coradini, A., Filacchione, G., et al. 2015, *Sci*, 347, aaa0628

Cochran, A. L., Barker, E. S., & Gray, C. L. 2012, *Icar*, 218, 144

Cochran, A. L., & McKay, A. J. 2018, *ApJL*, 854, L10

Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, *MNRAS*, 354, 1133

Cordiner, M. A., Biver, N., Crovisier, J., et al. 2017, *ApJ*, 837, 177

Cordiner, M. A., Milam, S. N., Biver, N., et al. 2020, *NatAs*, 4, 861

Cordiner, M. A., Palmer, M. Y., de Val-Borro, M., et al. 2019, *ApJL*, 870, L26

Cordiner, M. A., Remijan, A. J., Boissier, J., et al. 2014, *ApJL*, 792, L2

Cottin, H., & Fray, N. 2008, *SSRv*, 138, 179

Coulson, I. M., Cordiner, M. A., Kuan, Y.-J., et al. 2017, *AJ*, 153, 169

Coulson, I. M., Liu, F.-C., Cordiner, M. A., et al. 2020, *AJ*, 160, 182

Crifo, J. F., Rodionov, A. V., & Bockelée-Morvan, D. 1999, *Icar*, 138, 85

Crovisier, J. 1987, *A&AS*, 68, 223

Crovisier, J., & Encrenaz, T. 1983, *A&A*, 126, 170

Crovisier, J., Leech, K., Bockelee-Morvan, D., et al. 1997, *Sci*, 275, 1904

Davidsson, B. J. R., Birch, S., Blake, G. A., et al. 2021, *Icar*, 354, 114004

Davidsson, B. J. R., Samarasinha, N. H., Farnocchia, D., & Gutiérrez, P. J. 2022, *MNRAS*, 509, 3065

de Val-Borro, M., Milam, S. N., Cordiner, M. A., et al. 2018, *ATel*, 11254, 1

Delaunay, B. 1934, Bulletin de l'Académie des Sciences de l'URSS, 6, 793

Dello Russo, N., Kawakita, H., Vervack, R. J., & Weaver, H. A. 2016, *Icar*, 278, 301

Festou, M. C., Gunnarsson, M., Rickman, H., Winnberg, A., & Tancredi, G. 2001, *Icar*, 150, 140

Fougere, N., Altwegg, K., Berthelier, J. J., et al. 2016, *A&A*, 588, A134

Fougere, N., Combi, M. R., Tenishev, V., et al. 2012, *Icar*, 221, 174

Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, *A&A*, 505, 629

Garcia-Berrios, E., Cordiner, M. A., Charnley, S. B., de Val-Borro, M., & Biver, N. 2020, in 14th Europlanet Science Congress 2020, EPSC2020-426, https://www.epsc2020.eu/

Garrod, R. T., & Pauly, T. 2011, *ApJ*, 735, 15

Gasc, S., Altwegg, K., Balsiger, H., et al. 2017, *MNRAS*, 469, S108

Gordon, I., Rothman, L., Hargreaves, R., et al. 2021, *JQSRT*, 277, 107949

Green, S. 1975, *JChPh*, 62, 2271

Gunnarsson, M. 2003, *A&A*, 398, 353

Gunnarsson, M., Bockelée-Morvan, D., Biver, N., Crovisier, J., & Rickman, H. 2008, *A&A*, 484, 537

Gunnarsson, M., Bockelée-Morvan, D., Winnberg, A., et al. 2003, *A&A*, 402, 383

Gunnarsson, M., Rickman, H., Festou, M. C., Winnberg, A., & Tancredi, G. 2002, *Icar*, 157, 309

Haser, L. 1957, BSRSL, 43, 740

Hindmarsh, A. C. 2019, ODEPACK: Ordinary differential equation solver library, Astrophysics Source Code Library, ascl:1905.021

Hogerheijde, M. R., Qi, C., de Pater, I., et al. 2009, *AJ*, 137, 4837

Huebner, W. F., & Mukherjee, J. 2015, *P&SS*, 106, 11

Hutson, J. M., & Green, S. 1994, MOLSCAT computer code, v.14, https://www.theochem.ru.nl/molscat/Section1.html

Ioppolo, S., van Boheemen, Y., Cuppen, H. M., van Dishoeck, E. F., & Linnartz, H. 2011, *MNRAS*, 413, 2281

Ip, W. H. 1986, *ApJ*, 300, 456

Jaeger, S. 2008, in ASP Conf. Ser. 394, Astronomical Data Analysis Software and Systems XVII, ed. R. W. Argyle, P. S. Bunclark, & J. R. Lewis (San Francisco, CA: ASP), 623

Kawakita, H., Dello Russo, N., Vervack, Ron, J., et al. 2014, *ApJ*, 788, 110
Klos, J., & Lique, F. 2018, in Cold Chemistry: Molecular Scattering and Reactivity Near Absolute Zero, ed. O. Dulieu & A. Osterwalder (London: Royal Society of Chemistry), 46

Läuter, M., Kramer, T., Rubin, M., & Altwegg, K. 2019, MNRAS, 483, 852

Lisse, C. M., Young, L. A., Cruikshank, D. P., et al. 2021, Icar, 356, 114072

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, AJ, 155, 182

Luspay-Kuti, A., Hässig, M., Fuselier, S. A., et al. 2015, A&A, 583, A4

Marboeuf, U., & Schmitt, B. 2014, Icar, 242, 225

Markwardt, C. 2012, MPFIT: Robust nonlinear least squares curve fitting, Astrophysics Source Code Library, ascl:1208.019

McKay, A. J., DiSanti, M. A., Kelley, M. S. P., et al. 2019, AJ, 158, 128

Migliorini, A., Piccioni, G., Capaccioni, F., et al. 2016, A&A, 589, A45

Milam, S. N., Savage, C., Ziurys, L. M., & Wyckoff, S. 2004, ApJ, 615, 1054

Mumma, M. J., Bonev, B. P., Villanueva, G. L., et al. 2011, ApJL, 734, L7

Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471

Ndengué, S. A., Dawes, R., & Gatti, F. 2015, JPCA, 119, 7712

Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2008, ApJ, 678, 1032

Ootsubo, T., Kawakita, H., Hamada, S., et al. 2012, ApJ, 752, 15

Paganini, L., Villanueva, G. L., Lara, L. M., et al. 2010, ApJ, 715, 1258

Penteado, E. M., Boogert, A. C. A., Pontoppidan, K. M., et al. 2015, MNRAS, 454, 531

Qasim, D., Chuang, K. J., Fedoseev, G., et al. 2018, A&A, 612, A83

Qasim, D., Fedoseev, G., Chuang, K. J., et al. 2020, NatAs, 4, 781

Rabli, D., & Flower, D. R. 2010, MNRAS, 406, 95

Roth, N. X., Milam, S. N., Cordiner, M. A., et al. 2021a, ApJ, 921, 14

Roth, N. X., Milam, S. N., Cordiner, M. A., et al. 2021b, PSJ, 2, 55

Sun, Z.-F., van Hemert, M. C., Loreau, J., et al. 2020, Sci, 369, 307

Surin, L. A., Fourzikov, D. N., Giesen, T. F., et al. 2007, JPCA, 111, 12238

Tenishev, V., Combi, M., & Davidsson, B. 2008, ApJ, 685, 659

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

Vissers, G. W. M., Wormer, P. E. S., & van der Avoird, A. 2003, PCCP, 5, 4767

Wierzchos, K., & Womack, M. 2018, AJ, 156, 34

Womack, M., Sarid, G., & Wierzchos, K. 2017, PASP, 129, 031001

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062

Zakharov, V., Bockelée-Morvan, D., Biver, N., Crovisier, J., & Lecacheux, A. 2007, A&A, 473, 303