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

M. A. Cordiner  Astrochemistry Laboratory, NASA Goddard Space Flight Center,
8800 Greenbelt Road, Greenbelt, MD 20771, USA.  Department of Physics,
Catholic University of America, Washington, DC 20064, USA.  I. M. Coulson
East Asian Observatory, 660 N. A’ohoku Place, Hilo, HI 96720, USA.  E. Garcia-
Berrios  Astrochemistry Laboratory, NASA Goddard Space Flight Center, 8800
Greenbelt Road, Greenbelt, MD 20771, USA.  Department of Physics, Catholic
University of America, Washington, DC 20064, USA.  C. Qi  Harvard-Smithsonian
Center for Astrophysics, 60 Garden Street, MS 42, Cambridge, MA 02138, USA.
F. Lique  Université de Rennes 1, Campus de Beaulieu, 263 avenue du Général
Leclerc, 35042 Rennes Cedex, France.  M. Zołtowski  Université de Rennes 1,
Campus de Beaulieu, 263 avenue du Général Leclerc, 35042 Rennes Cedex, France.
LOMC - UMR 6294, CNRS-Université du Havre, 25 rue Philippe Lebon, BP1123, F-76
063 Le Havre cedex, France.  M. de Val-Borro  Astrochemistry Laboratory, NASA
Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA.
Department of Physics, Catholic University of America, Washington, DC 20064,
USA.  Y.-J. Kuan  National Taiwan Normal University, Ting-Chou Road, Taipei
11677, Taiwan, ROC.  Institute of Astronomy and Astrophysics, Academia Sinica,
Taipei 106, Taiwan, ROC.  W.-H. Ip  Graduate Institute of Space Science,
National Central University, Taoyuan City, Taiwan.  S. Mairs  East Asian
Observatory, 660 N. A’ohoku Place, Hilo, HI 96720, USA.  N. X. Roth
Astrochemistry Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt
Road, Greenbelt, MD 20771, USA.  Department of Physics, Catholic University of
America, Washington, DC 20064, USA.  S. B. Charnley  Astrochemistry
Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt,
MD 20771, USA.  S. N. Milam  Astrochemistry Laboratory, NASA Goddard Space
Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA.  W.-L Tseng
National Taiwan Normal University, Ting-Chou Road, Taipei 11677, Taiwan, ROC.
Y.-L Chuang  National Taiwan Normal University, Ting-Chou Road, Taipei 11677,
Taiwan, ROC.

###### 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/H2O and
N${}_{2}^{+}$/H2O 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 (JCMT) and
Submillimeter Array (SMA) to elucidate the outgassing behaviour of C/2016 R2.
Results are analyzed using a new, time-dependent, three dimensional radiative
transfer code (SUBLIME), 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 $(3.8-7.6)\times 10^{28}$ s-1
between 2018-01-13 and 2018-02-01, with a mean value of $(5.3\pm 0.6)\times
10^{28}$ s-1 at $r_{H}=2.8$–2.9 au. The emission is concentrated in a near-
sunward jet, with an outflow velocity $0.51\pm 0.01$ km s-1, compared to
$0.25\pm 0.01$ km s-1 in the ambient (and night-side) coma. Evidence was also
found for an extended source of CO emission, possibly due to icy grain
sublimation around $1.2\times 10^{5}$ 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, CO2, N2 and
CH3OH, and a predominantly frozen (or less abundant), polar phase containing
more H2O, CH4, H2CO and HCN.

Comets: individual (C/2016 R2 (PanSTARRS)) — Submillimeter astronomy — High
resolution spectroscopy — Radio interferometry — Radiative transfer
simulations

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

##  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${}_{2}^{+}$ (Cochran & McKay, 2018), consistent with very strong outgassing
of CO and an above average N2/CO ratio in the nucleus. Both CO and N2
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 H2O) implies that C/2016 R2 formed very
cold, and was maintained at temperatures $\lesssim 20$ K for the duration of
its lifetime.

Followup observations at infrared and millimeter wavelengths confirmed the
presence of an exceptionally CO-rich coma, with a low production rate of H2O
(usually the dominant cometary volatile) and extremely unusual abundance
ratios for other species compared with the typical cometary population
(Wierzchos & Womack, 2018; Biver et al., 2018; McKay et al., 2019). Due to
their very different sublimation temperatures, the CO to H2O ratio in cometary
comae is known to vary as a function of temperature of the nucleus, so some
enhancement in CO/H2O is expected at the relatively large heliocentric
distance ($r_{H}$) at which this comet was observed. However, the coma
composition of C/2016 R2 clearly differs from other comets observed at similar
$r_{H}$. Based on Table 7 of Biver et al. (2018), the abundance ratios
CO:H2O:CH3OH:HCN (normalized to CH3OH) were 94:0.3:1:0.004 in C/2016 R2 at
$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 H2O 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 H2O and HCN depletions. Cometary
and interstellar ice observations have identified the presence of separate
‘polar’ and ‘apolar’ ice phases, dominated by H2O and [CO + CO2], 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 + CO2 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, leads to ambiguity in its
interpretation, and estimates for the CO production rate (in 2018 January)
range from $Q({\rm CO})=(4.6\pm 0.4)\times 10^{26}$ s-1 (Wierzchos & Womack,
2018) to $(10.6\pm 0.5)\times 10^{26}$ 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 behaviour 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 12CO emission in this comet, we searched for the 13CO
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:CH3OH:H2CO
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) dataset is analyzed using a new,
time-dependent, three dimensional radiative transfer code (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 provides new insights into the coma properties and intrinsic nature
of C/2016 R2 (PanSTARRS), one of the most unusual comets of our time.

##  2 Observations

###  2.1 James Clerk Maxwell Telescope (JCMT)

Observations of comet C/2016 R2 (PanSTARRS) were conducted during the period
2018-01-13 to 2018-02-01 using the 15 m JCMT atop Mauna Kea. The comet’s
position on the sky was tracked using Jet Propulsion Laboratory (JPL) Horizons
ephemeris #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 mm-wave receiver
(operating at 212 to 274 GHz) and the HARP 16-element sub-mm 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 record a
spectrum at an offset point in the coma. The $30^{\prime\prime}$ spacing
between receptors means that the target field is under-sampled during stare
observations, as a result of the $14^{\prime\prime}$ telescope beam FWHM at
346 GHz ($20^{\prime\prime}$ at 230 GHz). Jiggle observations, on the other
hand, are designed to obtain spatially complete sampling of the target field.
We used the HARP4 jiggle pattern to obtain spectra for each point on a
$7^{\prime\prime}\times 7^{\prime\prime}$ grid, covering a $2^{\prime}\times
2^{\prime}$ map area. Position switching (offset by $300^{\prime\prime}$ 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
sub-optimal 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 ACSIS digital autocorrelation spectrometer. For
the 13CO and C18O spectral windows (observed simultaneously), 61 kHz
resolution was used, and for CH3OH, a relatively coarse resolution of 448 kHz
(over a 1000 MHz bandpass) was used to cover multiple transitions of CH3OH 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 $\sim 2^{\prime\prime}$ 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^{\prime\prime}$. Comparison of the actual telescope tracking
position (apparent RA and dec.) with the ephemeris showed tracking errors to
be less than $0.1^{\prime\prime}$ 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 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.

Table 1: Log of JCMT spectral line observations

Date | Species | Transition | Obs. Mode | $\tau_{0}$ | Int. | $\int{T_{MB}dv}$ | $\bar{v}$ | $r_{H}$ | $\Delta$ | $\phi$
---|---|---|---|---|---|---|---|---|---|---
|  | $-$ |  |  |  | (s) | (mK km s-1) | (km s-1) | (au) | (au) | (∘)
2018-01-13.215 | CO | 3$-$ | 2 | HARP stare | 0.10 | 600 | 870(92) | $-0.21(0.04)$ | 2.88 | 2.13 | 14.8
2018-01-13.253 | HCN | 4$-$ | 3 | HARP stare | 0.10 | 1200 | $<78$ | — | 2.88 | 2.13 | 14.9
2018-01-13.338 | CH3OH | 7$-$ | 6 | HARP stare | 0.11 | 900 | 73(20)a |  $-0.08(0.13)$a | 2.88 | 2.13 | 14.9
2018-01-13.401 | HCO+ | 4$-$ | 3 | HARP stare | 0.11 | 900 | $<93$ | — | 2.88 | 2.13 | 14.9
2018-01-13.437 | H2CO |  $5_{1,5}$$-$ | $4_{1,4}$ | HARP stare | 0.11 | 900 | $<110$ | — | 2.88 | 2.13 | 14.9
2018-01-14.206 | CO | 2$-$ | 1 | RxA3m | 0.17 | 600 | 542(71) | $-0.18(0.08)$ | 2.87 | 2.14 | 15.1
2018-01-14.248 | CO | 3$-$ | 2 | HARP stare | 0.17 | 120 | 1078(154) | $-0.24(0.10)$ | 2.87 | 2.14 | 15.1
2018-01-14.260 | CO | 3$-$ | 2 | HARP stare | 0.17 | 600 | 926(103) | $-0.22(0.05)$ | 2.87 | 2.14 | 15.1
2018-01-14.273 | CO | 3$-$ | 2 | HARP stare | 0.17 | 120 | 831(131) | $-0.21(0.06)$ | 2.87 | 2.14 | 15.2
2018-01-14.333 |  13CO | 3$-$ | 2 | HARP stare | 0.16 | 900 | $<229$ | — | 2.87 | 2.14 | 15.2
2018-01-14.333 | C18O | 3$-$ | 2 | HARP stare | 0.16 | 900 | $<282$ | — | 2.87 | 2.14 | 15.2
2018-01-14.367 | CO | 3$-$ | 2 | HARP jiggle | 0.15 | 1780 | 834(114) | $-0.25(0.04)$ | 2.87 | 2.14 | 15.2
2018-01-14.384 | CO | 3$-$ | 2 | HARP stare | 0.15 | 120 | 881(140) | $-0.27(0.12)$ | 2.87 | 2.14 | 15.2
2018-01-15.197 | CO | 3$-$ | 2 | HARP stare | 0.09 | 120 | 1084(130) | $-0.21(0.06)$ | 2.87 | 2.14 | 15.4
2018-01-15.213 | CO | 3$-$ | 2 | HARP jiggle | 0.08 | 1780 | 1145(131) | $-0.23(0.02)$ | 2.87 | 2.14 | 15.4
2018-01-15.231 | CO | 3$-$ | 2 | HARP stare | 0.08 | 120 | 1028(117) | $-0.23(0.05)$ | 2.87 | 2.14 | 15.4
2018-01-15.253 |  13CO | 3$-$ | 2 | HARP stare | 0.08 | 900 | $<84$ | — | 2.87 | 2.14 | 15.4
2018-01-15.253 | C18O | 3$-$ | 2 | HARP stare | 0.08 | 900 | $<112$ | — | 2.87 | 2.14 | 15.4
2018-01-15.271 | CO | 3$-$ | 2 | HARP stare | 0.08 | 120 | 1015(114) | $-0.19(0.05)$ | 2.87 | 2.14 | 15.4
2018-01-20.247 | CO | 3$-$ | 2 | HARP stare | 0.06 | 600 | 745(76) | $-0.23(0.03)$ | 2.85 | 2.18 | 16.8
2018-01-25.322 | CO | 3$-$ | 2 | HARP stare | 0.14 | 600 | 881(96) | $-0.23(0.04)$ | 2.83 | 2.23 | 17.9
2018-01-31.338 | CO | 3$-$ | 2 | HARP stare | 0.08 | 600 | 910(93) | $-0.23(0.03)$ | 2.80 | 2.29 | 19.1
2018-02-01.278 | CO | 3$-$ | 2 | HARP stare | 0.14 | 600 | 768(84) | $-0.18(0.04)$ | 2.80 | 2.30 | 19.3

a Integrated over the three strongest CH3OH transitions (see Section 4.6.)

###  2.2 Submillimeter Array (SMA)

Interferometric observations of C/2016 R2 were made on 2018-02-21 using the
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 observations, resulting in a spatial resolution
of $\approx 4.4^{\prime\prime}\times 7.7^{\prime\prime}$ 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 SWARM correlator provides 8 GHz bandwidth per sideband,
divided into four equal sized chunks with uniform spectral resolution of 140
kHz (0.18 km s-1 at 230 GHz). Position and Doppler-tracking of the comet over
a 5 hour observing period was performed using JPL Horizons orbital solution
#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
MIR software package111http://www.cfa.harvard.edu/$\sim$cqi/mircook.html. The
calibrated SMA data were then exported in uvfits format for subsequent imaging
and analysis using the CASA software package version 5.1 (Jaeger, 2008).

Interferometric imaging and deconvolution was performed using the Hogbom clean
algorithm, with a $20^{\prime\prime}$ 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^{\prime\prime}$.

##  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 topocentric 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 night side.

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 and rotational temperature ($T_{rot}$) were probably
close to steady state during the time period of our observations, although
variations in both quantities could have cancelled each other to some degree —
the CO $J=3-2$ line intensity varies by a factor of 1.5 between
$T_{rot}=15$–25 K (a plausible range for this comet based on the findings of
Biver et al. (2018); see also Section 4.2).

Figure 1: Time sequence of CO $J=3-2$ spectra observed using JCMT (in the
topocentric 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).
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. Horizontal
line shows the error-weighted mean. Error bars include the statistical
uncertainty, with an additional 10% calibration uncertainty added in
quadrature.

###  3.2 CO Single-dish Mapping with JCMT HARP

HARP jiggle-map spectral cubes from 2018-01-14 and 2018-01-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\sigma{n_{c}^{-0.5}}$, where
$\sigma$ is the average RMS noise of the data cube (equal to $T_{MB}=72$ mK km
s-1) 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^{\prime\prime}\times 120^{\prime\prime}$ HARP data cube, was $164\pm 19$
mK km s-1 on 2018-01-14 and $202\pm 22$ mK km s-1 on 2018-01-15. These values
are sufficiently similar (given the uncertainties) to justify combining the
two datasets. The CO spatial distributions were also apparently identical
(within the uncertainties) on the two dates, with consistent radial intensity
profiles.

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 is detected (at the $3\sigma$ level) up to a radial distance
$82,000$ km west of image-center, which is in a direction similar to the sky-
projected sunward vector ($103^{\circ}$ clockwise from North), whereas the
$3\sigma$ contour extends only $58,000$ km to the east. The 2nd contour (at
$6\sigma$) also shows an excursion in the sunward direction into a pixel with
an intensity $2.2\sigma$ 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.

Figure 3: Spectrally integrated CO $J=3-2$ map for comet C/2016 R2
(PanSTARRS), from the average of HARP jiggle observations on 2018-01-14 and
2018-01-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 lower right. Contours are in units of $3\sigma$
and the axes are aligned with the equatorial (RA/dec.) grid.

###  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.5^{\prime\prime}$ 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 N-S direction, and
is less well resolved in the E-W 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\sigma$ 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.

Figure 4: Spectrally integrated CO $J=2-1$ emission map for comet C/2016 R2,
obtained using the SMA on 2018-02-21. The FWHM (and orientation) of the
elliptical Gaussian restoring beam is indicated by the ellipse (upper right),
and the sky-projected solar and orbital trail vectors are shown lower right.
Contours are in units of $3\sigma$ and the axes are aligned with the
equatorial (RA/dec.) grid, with the origin at the SMA phase tracking center.
Negative contours are shown with a dashed line style.

##  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 (SUBlimating gases in LIME) for simulating the rotational
emission lines from cometary coma molecules in three dimensions (two spatial
and one spectral), where LIME refers to the LIne Modeling Engine by (Brinch &
Hogerheijde, 2010).

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). Our radiative transfer code is
based on the open-source version 1.9.3 of LIME (available at
https://github.com/lime-rt/lime/releases). LIME, however, employs the steady
state approximation.

We therefore modified substantially the LIME 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 Bockelee-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 on an unstructured 3D grid
(Delaunay, 1934), which can be 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 $\sim$ km-sized nucleus
to the $\sim 10^{6}$ km-scale outer coma) can thus be sampled much more
efficiently than with a uniform grid spacing. The final level populations
$N_{i}$ as a function of radius are mapped onto the Delaunay grid for
raytracing. 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 article, 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, subject to a
logarithmic density distribution with radius between $r_{c}=500$ m (the
assumed radius of the nucleus) and $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 raytracing, 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^{\prime\prime}$ and 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 sub-pixel scales towards 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\times 30$
rays) for each pixel within the central $4\times 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_{1}J_{2}}$) is therefore required to correctly model the CO emission
from the comet. Previous studies (Wierzchos & Womack, 2018; Biver et al.,
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.222Ndengué 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.

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 Kłos & 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 rotational state characterized by quantum number $J_{1}$), and
the second is the ‘collider’ ($J_{2}$). Collisional rate coefficients
$k_{J_{1}J_{2}\to J_{1}^{\prime}j_{2}^{\prime}}$ 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}^{\prime}J_{2}^{\prime}}(T)=\left(\frac{8}{\pi\mu
k^{3}_{B}T^{3}}\right)^{\frac{1}{2}}\int\limits_{0}^{\infty}\sigma_{J_{1}J_{2}\to
J_{1}^{\prime}J_{2}^{\prime}}(E_{c})E_{c}e^{-\frac{E_{c}}{k_{B}T}}dE_{c}$ (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
collisional energy. Rate coefficients used in the radiative transfer model
($k_{J_{1}\to J_{1}^{\prime}}(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:

$\begin{split}k_{J_{1}\to
J_{1}^{\prime}}(T)=N_{J_{2}=0}\sum_{J_{2}^{\prime}}k_{J_{1},J_{2}=0\to
J_{1}^{\prime},J_{2}^{\prime}}\\\
+N_{J_{2}=1}\sum_{J_{2}^{\prime}}k_{J_{1},J_{2}=1\to
J_{1}^{\prime},J_{2}^{\prime}}+...\end{split}$ (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}}}$
(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 ($\sim 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-H2 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-01-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_{ij}n$
(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 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: C1 and C2, corresponding to solid angle regions $\Omega_{1}$,
$\Omega_{2}$ with independent production rates ($Q_{1}$, $Q_{2}$) and outflow
velocities ($v_{1}$, $v_{2}$). 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 2-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 parameterisations 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 C1 is defined by a conical region about the subsolar
point, with half-opening angle $\theta$, whereas C2 is the remaining, ambient
coma. We used the MPFIT nonlinear least-squares routine (Markwardt, 2012) to
find the optimal set of parameters ($Q_{1}$, $Q_{2}$, $v_{1}$, $v_{2}$,
$\theta$). The model images were convolved with the (Gaussian) JCMT beam
pattern and normalized by their respective beam efficiency factors before
comparing them with observations. Statistical ($1\sigma$) 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}$ s-1, $Q_{2}=(2.2\pm 0.3)\times 10^{28}$
s-1, $v_{1}=0.50\pm 0.01$ km s-1, $v_{2}=0.29\pm 0.02$ km s-1,
$\theta=27^{\circ}\pm 5^{\circ}$, and $T_{kin}=18.7\pm 2.3$ K. The retrieved
kinetic temperature is consistent with the values $18.6\pm 2.6$ to $24.2\pm
7.9$ derived by Biver et al. (2018) using multiple lines of CH3OH on 2018
January 23-24. The total CO gas production rate from our model is $Q_{t}({\rm
CO})=(3.9\pm 0.7)\times 10^{28}$ s-1, and the ratio of production rates per
unit solid angle between the jet and ambient coma is
$R_{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_{Q}$ is a more physically meaningful quantity than the
simple ratio of production rates ($Q_{1}/Q_{2}$), revealing the degree to
which the comet’s activity is enhanced due to heating by the Sun in the
vicinity of the sub-solar 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 fit was not significantly
improved. Additional models were run allowing different $T_{kin}$ values for
the two coma components, with best-fitting results $T_{kin_{1}}=19.3\pm 3.7$ K
and $T_{kin_{2}}=17.9\pm 4.8$ K, but again, the overall quality of fit was not
improved. The consistency of $T_{kin_{1}}$ and $T_{kin_{2}}$ (within their
respective errors) provides further justification for adopting a uniform
$T_{kin}$ value throughout the rest of this study.

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 red
curves

###  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_{t}({\rm
CO})=(5.3\pm 0.2)\times 10^{28}$ s-1, $R_{Q}=7.4\pm 0.4$, $v_{1}=0.51\pm 0.01$
km s-1, $v_{2}=0.25\pm 0.01$ km s-1, $\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-01-14. Allowing for an additional 10% calibration error, the total
uncertainty on the CO production rate is $0.6\times 10^{28}$ s-1.

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. Dot-dashed grey curve shows the best-fitting 1D
model, assuming a spherically-symmetric coma.

Individual fits were also performed to the time series of 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 14th. On February 1st, 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 behaviour over time.

Table 2: Best-fit model results for individual JCMT (and SMA) spectra of
C/2016 R2

Date | $Q_{t}$ | $R_{Q}$ | $v_{1}$ | $v_{2}$ | $\theta$ | $\phi_{jet}$
---|---|---|---|---|---|---
| ($10^{28}$ s-1) |  | (km s-1) | (km s-1) | (deg.) | (deg.)
2018-01-13 | 4.5 (0.7) | 6.8 (4.5) | 0.51 (0.01) | 0.23 (0.02) | 60 (5) | 18 (6)
2018-01-14 | 4.0 (0.9) | 10.8 (2.9) | 0.50 (0.01) | 0.29 (0.02) | 32 (7) | 7 (10)
2018-01-15 | 7.4 (1.6) | 9.6 (3.1) | 0.52 (0.02) | 0.31 (0.03) | 77 (13) | 37 (19)
2018-01-20 | 3.9 (1.8) | 17.2 (10.8) | 0.50 (0.01) | 0.24 (0.01) | 42 (19) | 34 (11)
2018-01-25 | 5.0 (0.9) | 10.4 (1.7) | 0.53 (0.01) | 0.24 (0.02) | 58 (6) | 20 (7)
2018-01-31 | 5.7 (0.8) | 7.3 (0.7) | 0.53 (0.01) | 0.26 (0.01) | 64 (3) | 12 (4)
2018-02-01 | 4.4 (0.8) | 4.9 (0.9) | 0.49 (0.01) | 0.33 (0.02) | 46 (8) | 8 (12)
JCMT avg. | 5.3 (0.6) | 7.4 (0.4) | 0.51 (0.01) | 0.25 (0.01) | 62 (2) | 24 (2)
2018-02-21a | 6.7 (0.9) | 9.2 (1.5) | 0.64 (0.02) | 0.20 (0.02) | 74 (7) | 21b

Notes — $1\sigma$ statistical uncertainties on each value are given in
parentheses. Uncertainties on $Q_{t}$ include a 10% intensity calibration
error, added in quadrature with the statistical error. aSMA best-fitting
parent visibility model. bJet 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.

###  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 inetnsity 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.

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 CO2 photolysis, with $Q({\rm CO_{2}})/Q({\rm CO})=6.5$
(center panel), and CO from the nucleus, with additional CO shell at
$r_{c}\sim 1.4\times 10^{5}$ km (right panel). Contour spacings are the same
as in Figure 3.  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 text). An improved fit is
obtained at large radii when including a CO source from CO2 photolysis (with a
large $Q({\rm CO_{2}})/Q({\rm CO})$ ratio of 6.5), or with the inclusion of an
extended CO shell at $r_{s}=1.2\times 10^{5}$, with a density enhancement
factor of $f=2.0$.

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),
which implies the presence of an excess of CO emission at large radii, that
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})/\sigma)^{2}/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 model and observations is due to statistical chance.
With $\chi_{R}^{2}=2.1$ 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 ($\Gamma_{\rm CO}$) would not reproduce the observed
emission excess — at $r_{H}=2.85$ au, $\Gamma_{\rm CO}=9.2\times 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\times 10^{6}$ km. CO photodissociation is therefore
negligible within the 80,000 km JCMT field of view, so reducing $\Gamma_{\rm
CO}$ by a factor of a few (its range of uncertainty) has no noticeable impact
on the modeled CO density profile.

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_{ij}$) into this level are closely
balanced by radiative transitions out of the level, so on the larger distance
scales $\gtrsim 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_{kin}=16.4$ K (the $1\sigma$ lower limit
derived in Section 4.2), the best-fitting $\chi_{R}^{2}$ value is 1.6, with
$P=0.17$, which still does not represent a very good fit. Lower values of
$T_{kin}$ also appear less likely given the range of possible $T_{kin}$ values
(16–32 K) observed by Biver et al. (2018), so we seek alternative explanations
for the shape of the $\langle{T_{MB}}\rangle_{az}$ profile.

Figure 9: Fractional energy level populations as a function of radius for the
lowest 7 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.

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 photodissociation 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-1, $Q_{p}=7.0\times 10^{28}$ s-1) is shown in
Figure 8, and systematically under-fits 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 $Q_{t}=5.5\times 10^{28}$ s-1, and
an additional CO source in the coma from CO2 photolysis (CO2 \+
$h\nu\longrightarrow$ CO + O at $\Gamma_{\rm CO_{2}}=1.5\times 10^{-7}$ s-1;
Huebner & Mukherjee 2015). Although this composite model may appear plausible
at first glance, and accurately reproduces the CO emission profile, it
requires a CO2 production rate 6.5 times larger than that of CO. Considering
McKay et al. (2019) derived an upper limit of $Q({\rm CO_{2}})\lesssim
1.5\times 10^{28}$ s-1 based on Spitzer observations from 2018-02-12 and
2018-02-21 (_i.e._ 3.7 times less than the parent CO source), it seems
unlikely that the comet could have been producing so much CO2 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 $\sim 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,
$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 $\sim 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 km s-1, a factor of $\sim$ two density
enhancement at $r_{s}=1.2\times 10^{5}$ km would be consistent with a
corresponding drop in $Q({\rm CO})$ around 67 hours earlier (on January 12th),
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 14th and 15th 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}(\rm CO)=(6.1\pm 0.8)\times 10^{28}$ s-1 to $(5.4\pm 1.1)\times 10^{28}$
s-1. 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-1) 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, re-gridding 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-1, $v_{2}=0.20\pm 0.02$ km s-1 and
$\theta=74^{\circ}\pm 7^{\circ}$; the jet outflow axis was held fixed towards
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 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.

Figure 10: CO $J=2-1$ spectrum observed with the SMA (extracted at the CO
emission peak), with best-fitting 2-component SUBLIME model overlaid using a
red curve.

From this best-fitting ‘base’ model, additional models were constructed
assuming (1) CO solely as a daughter species, with $\Gamma_{p}=7.3\times
10^{-5}$ s-1 derived from the fit to the JCMT HARP data (Section 4.4), (2) CO
as a parent with additional outer-coma source from CO2 photolysis, and (3) CO
as a parent with additional extended shell at $r_{s}=1.2\times 10^{5}$ km,
with $f=2.0$. The resulting 3D model images were integrated in the spectral
domain then multiplied by the (FWHM $=55^{\prime\prime}$) 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 observations). A power-
law curve ($ax^{b}$) was fit through each set of model results, (with $x$ the
baseline length, and $a$ and $b$ 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_{t}(\rm CO)=(6.7\pm 0.6)\times 10^{28}$ s-1 (the
uncertainty increases to $0.9\times 10^{28}$ s-1 after inclusion of a 10%
amplitude calibration error). This value is $1.6\sigma$ from the mean CO
production rate obtained using JCMT, corresponding to a barely-significant
increase in $Q_{t}$ with time.

Figure 11: CO $J=2-1$ visibility amplitudes _vs._ baseline length for C/2016
R2 observed using the SMA, including $1\sigma$ statistical error bars. Model
visibility curves are overlaid for four different CO distributions. The CO
parent, parent + CO2 photolysis, and parent + shell model curves all lie on
top of each other, whereas the CO daughter curve differs significantly, under-
fitting the observations at large baseline (small angular scales).

Three of the model curves (parent, parent + CO2 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 the 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 daughter
species in this comet, although significant production at large radii from CO2
photolysis is still possible.

###  4.6 Other Molecules from JCMT: CH3OH, HCN, H2CO, 13CO and HCO+

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,
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_{K}=7_{0}-6_{0}\,E$, $7_{-1}-6_{-1}\,E$ and $7_{0}-6_{0}\,A^{+}$. 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-1), with a spectrally-
integrated intensity (in the range $-1.0$ to $0.5$ km s-1 shown by the blue
shaded region) of $\int{T_{MB}dv}=73\pm 20$ mK km s-1, corresponding to a
$3.6\sigma$ detection. Fitting the same three transitions simultaneously using
SUBLIME (allowing only the CH3OH abundance to vary) gave $Q({\rm
CH_{3}OH})=(6.7\pm 2.2)\times 10^{26}$ s-1, which corresponds to a CH3OH/CO
abundance ratio of $1.3\pm 0.4$% at the nucleus. For this model, we adopted
CH3OH–CO collisional transition rates based on the CH3OH–H2 rates from Rabli &
Flower (2010), with solar pumping rates from Roth et al. (2021a).
Uncertainties of a factor of five in the CH3OH–CO collision rates lead to at
most a 7% error on the CH3OH abundance.

Figure 12: Average (in velocity space) of the three strongest CH3OH
transitions in our JCMT observations ($7_{0}-6_{0}\,E$, $7_{-1}-6_{-1}\,E$ and
$7_{0}-6_{0}\,A^{+}$). A tentative feature is present around 0 km s-1 (the
comet’s rest velocity); blue shaded region shows the velocity range over which
significant CO emission was detected.

The HCN, H2CO and 13CO molecules were not detected in our data, as shown by
the spectra in Figure 13. Three-sigma upper limits were derived by comparing
the spectrally-integrated noise level (from $-0.8$ to 0.5 km s-1) 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 H2CO we used the Roth et al.
(2021a) rates, and for 13CO we used the same rates as for CO. Production rate
$3\sigma$ upper limits are $Q({\rm HCN})<8.3\times 10^{25}$ s-1, $Q({\rm
H_{2}CO})<2.3\times 10^{26}$ s-1 and $Q({\rm{}^{13}CO})<2.7\times 10^{27}$
s-1, corresponding to abundance ratios HCN/CO $<0.16$%, H2CO/CO $<0.45$% and
CO/13CO $>19$. The 12C/13C ratio in CO is therefore consistent with the
typical Solar System value of 89, as well as the value of $86\pm 9$ measured
in comet 67P (Altwegg et al., 2019).

We also obtained a non-detection of HCO+ $J=4-3$, with $\int{T_{MB}dv}<164$ mK
km s-1 (integrated over the velocity range $\pm 1.5$ km s-1). 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 SMT. However, the
CO production rate in C/2016 R2 was also less than 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).

Figure 13: HCN, H2CO and 13CO spectra observed using JCMT, showing no evidence
for detections of these molecules.

##  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({\rm CO})=(5.3\pm 0.6)\times 10^{28}$ s-1 (averaged over the
JCMT observing period 2018-01-13 to 2018-02-01). This was subject to a
marginally significant rise to $(6.7\pm 0.9)\times 10^{28}$ s-1 on 2018-02-21,
observed using the SMA. These results confirm C/2016 R2 as having among the
highest CO production rates ever observed in a comet — only a factor of 4 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\pm
0.4)\times 10^{28}$ s-1 between 2017-12-22 and 2018-01-16 observed by
Wierzchos & Womack (2018) using the SMT, and similar to the $(5.5\pm
0.9)\times 10^{28}$ s-1 observed on 2018-02-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_{H}=2.98$–2.73 au). The observations by Biver et al.
(2018) using the IRAM 30 m telescope on 2018-01-24 gave $Q({\rm CO})=(10.6\pm
0.5)\times 10^{28}$ s-1, which is approximately double our JCMT value. As
shown in Appendix C, the 1D versions of ours and Biver et al.’s models 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’ ($\gamma$; angle measured
from the Earth-comet vector), with outflow velocities $v=0.56$ km s-1 between
$\gamma=0$–60∘, $0.50$ km s-1 between $\gamma=60$–120∘, and no outflow for
$\gamma>120$. Our model, on the other hand, gives $v_{1}=0.51$ km s-1 in a
sunward-facing conical jet and $v_{2}=0.25$ km s-1 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 anti-
sunward sides of the nucleus. The large reduction in the outflow velocity of
our model across the entirety of the night side (most of which overlaps with
the $\gamma>120^{\circ}$ region of no outflow in Biver et al. (2018)’s model)
necessitates a lower total production rate in our model compared with that of
Biver et al. (2018).

Although our 2-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 contraints 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 2-component
model for C/2016 R2 leads to a $\sim 30$% lower CO production rate, as a
result of reduced outgassing on the anti-sunward 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 range 29–48% 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 H2O outgassing behaviour
was similar to that of CO and CH3OH (see Biver et al. 2018 Figure 14), the
McKay et al. (2019) H2O production rate may also need to be adjusted. However,
the correction factor is lower for OH (the H2O daughter fragment observed by
McKay et al. 2019), due to the OH kinetic energy gained in the photolysis
reaction H2O + $h\nu$ $\longrightarrow$ OH + H. Based on the outflow
velocities from our best-fitting two-component CO parent model (Section 4.3),
the vectorial model used by McKay et al. (2019) implies OH outflow velocities
of 1.05 km s-1 in the sunward jet and 0.78 km s-1 in the ambient coma (M.
Knight, private communication 2022). This geometry would correspond to a 22%
reduction in the H2O production rate, assuming $\theta=62^{\circ}$.

Abundance ratios from a given study tend to be less susceptible to model-
dependent and instrumental uncertainties. Our CH3OH/CO abundance ratio of
$1.3\pm 0.4$% matches Biver et al. (2018)’s value of $1.04\pm 0.08$%, and our
HCN/CO and H2CO/CO upper limits are also consistent with their values of
$(3.8\pm 1.0)\times 10^{-3}$% and $0.043\pm 0.006$%, respectively. The
CH3OH/CO ratio is significantly higher than McKay et al. (2019)’s upper limit
of 0.38%, which could be indicative of significant temporal variability in the
CH3OH outgassing rate during 2018 January. Alternatively, the higher CH3OH
abundance observed using radio spectroscopy could have been due to additional
CH3OH 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)

After correcting the McKay et al. (2019) H2O production rate for coma
asymmetry (Section 5.1), we find a CO/H2O ratio $\sim 194$, which is $\sim
3200$ times greater than the average value observed Oort cloud comets (Dello
Russo et al., 2016), and $\sim 42$ 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/H2O 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}({\rm CO})=24$ K _vs._
$T_{sub}({\rm H_{2}O})=152$ K), so a reduction in H2O outgassing is expected
for comets at heliocentric distances $\gtrsim 3$ au, where the ice temperature
falls below $T_{sub}({\rm H_{2}O})$. The CO/H2O 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 (Mumma & Charnley, 2011; A’Hearn et al., 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 H2O 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 H2O production
rates could be enhanced by several orders of magnitude by the presence of an
insulating dust mantle $\sim 5$–10 cm thick that hinders H2O 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 H2O sublimation (Lisse et al., 2021). In that
case the comet could maintain a lower H2O outgassing rate (relative to CO) for
a longer duration as it approached the sun. Either case may not require an
extremely anomalous CO/H2O 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 CH3OH and HCN (99 K and 95 K, respectively),
these two molecules were enriched in the coma (relative to H2O) by very
different amounts: the CH3OH/H2O ratio was 163 times the Oort cloud comet
average, whereas HCN/H2O was only 5.3 times the average. Evidently, such
enrichment patterns cannot be produced by reduced H2O outgassing alone, or by
a simple temperature-dependence of the CH3OH and HCN sublimation rates.
Similarly, the coma CH4/H2O ratio is a factor of 206 higher than average, but
this is several hundred times less than the CO/H2O 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}({\rm
H_{2}O})$ (_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 (H2O-dominated) ices at temperatures
well above their sublimation points, and may be only partially released until
the H2O sublimation temperature is reached ($\sim 150$ K). Release of trapped
hypervolatiles can also occur as H2O 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, CO2 and CH3OH in comet 67P (Läuter et al., 2019; Biver
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; Ioppolo et al., 2011; Garrod & Pauly, 2011). Observations of young
stellar objects also indicate mixing of CO, CO2 and CH3OH ices in an apolar
phase, distinct from the polar, H2O-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, A. et al., 2016; Gasc et al., 2017)
discovered spatial and temporal variations in the coma H2O/CO2 ratios,
implying non-uniform 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 H2O/CO2 ratio
intrinsic to the nucleus are also possible. It is therefore plausible that
CH3OH in C/2016 R2 exists as a component within an apolar ice matrix dominated
by CO and CO2, and some of this CH3OH is released into the coma when CO
sublimates (as observed). Meanwhile, a significant CH3OH component also likely
remains frozen as part of the comet’s polar, H2O-rich ices (Qasim et al.,
2018). HCN, on the other hand, could be primarily associated with the (still
frozen) H2O ice rather than the (sublimating) CO component, and would then be
only partially outgassed at $T_{sub}({\rm HCN})$ or $T_{AC}$. Future studies
of cometary CO, CH3OH and HCN spatial distributions could help test this
hypothesis.

To complete this picture, the C2H6 and NH3 upper limits from McKay et al.
(2019) are consistent with moderate-to-no enrichment (with respect to H2O),
and could therefore be associated primarily with the H2O-dominated ice, while
the strongly-enriched N2 and CO2 are associated more with CO. Observations and
laboratory studies show that interstellar CH4 ice tends to be more associated
with H2O than CO (Öberg et al., 2008; Qasim et al., 2020), while the location
of H2CO ice is less well constrained. The CH4 and H2CO enrichment factors in
C/2016 R2 (206 and 43, respectively) could thus be explained by release of
trapped volatiles (in H2O ice) above their respective sublimation temperatures
(31 K and 64 K) (Collings et al., 2004). The increasing similarity to typical
abundances (with respect to H2O) in the sequence CH4 $\to$ H2CO $\to$ HCN is
consistent with the decreasing volatility (increasing $T_{sub}$) of these
three species, such that they each behave progressively more like H2O. We
therefore postulate the existence of a rapidly sublimating (apolar) component
of ice in C/2016 R2 rich in CO, CO2, N2 and CH3OH, and a second (polar)
component containing more CH4, H2CO and HCN mixed with H2O. If the N2/CO ratio
is higher in the CO-rich ice phase than the H2O-rich ice, this could explain
why the N2/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 Gunnarsson et al. (2008)’s analysis of
the 29P coma, and we find a similar ratio of day-to-night hemisphere $Q({\rm
CO})$ and $v_{out}$ values in C/2016 R2.

Gunnarsson et al. (2002) also discovered an extended shell of CO emission
surrounding 29P, at a cometocentric distance $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 followup 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 CO2 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 CO2 production rate than found by McKay et al. (2019).

A time-variable, diffuse, ring-shaped feature (consistent with excess gas
emission) also appeared in Spitzer IRAC CO + CO2 images of C/2016 R2 on
2018-02-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, 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 CO2 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 2-component coma
model.

###  5.4 Uncertainties in the CO Excitation

Figure 14: Modeled CO rotational temperature derived using the 7 lowest energy
level levels, from our time-dependent SUBLIME model (solid blue curve) and the
steady-state version of the model (dotted grey 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^{\prime\prime}$ JCMT beam FWHM at the Geocentric distance of C/2016 R2
(2.14 au), at 346 GHz.

Our new CO–CO collisional rate coefficients are considered to be accurate to
within about a factor of 2 (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_{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) 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 (Bockelee-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 behaviour. 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 born in mind that this issue is less severe for molecules with larger
dipole moments (and therefore, smaller $\tau_{R}$ values) such as H2O, HCN,
CH3OH and H2CO, 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({\rm H_{2}O})=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^{\prime\prime}$ JCMT beam) is overestimated in the steady-state
model by 46%, whereas at $Q({\rm H_{2}O})=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-01-13 to 2018-02-21, we confirm the presence
of extremely strong, asymmetric CO outgassing from the hypervolatile-rich
comet C/2016 R2 (PanSTARRS), with $Q({\rm CO})$ in the range $(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$–$90^{\circ}$ (and an average of
$\theta=62\pm 2^{\circ}$, offset by $7\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 2.0 times 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 a value for $Q({\rm
H_{2}O})$ from McKay et al. (2019), corrected for asymmetric outgassing, our
CO/H2O ratio is $\sim 42$ times larger than seen in any comet to-date
(including the distant Centaur 29P), which, combined with previously-noted
chemical peculiarities, suggests 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, that prevented the initiation of a more
conventional, H2O-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, CO2, N2 as
well as CH3OH ices, and (2) a polar phase containing larger abundances of CH4,
H2CO and HCN mixed in with H2O ice. More observations of coma chemistry in
distant comets (at $r_{H}\gtrsim 2.5$ au, for which H2O 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 H2O sublimation
point.

This work was supported by the National Science Foundation under Grant No.
AST-2009253. The work of SNM, NXR, MAC, EGB and SBC was also supported by
NASA’s Planetary Science Division Internal Scientist Funding Program through
the Fundamental Laboratory Research work package (FLaRe). SNM and SBC 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.

##  Appendix A Equations of Radiative Transfer and Molecular Excitation

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

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

where $j_{\nu}$ and $\alpha_{\nu}$ 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 follows:

$j_{\nu}=\frac{h\nu}{4\pi}N_{j}A_{ij}\psi_{\nu},$ (A2)

$\alpha_{\nu}=\frac{h\nu}{4\pi}(N_{i}B_{ij}-N_{j}B_{ji})\psi_{\nu}$ (A3)

where $N_{i}$, $N_{j}$ are the number of gas particles per unit volume in
levels $i$ and $j$, respectively, and $\psi_{\nu}$ 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):

$\begin{split}\frac{dN_{i}}{dt}=-N_{i}\left[\sum_{j<i}A_{ij}+\sum_{j\neq
i}(B_{ij}J_{\nu}+k_{ij}n+G_{ij})\right]\\\ +\sum_{j>i}N_{j}A_{ji}+\sum_{j\neq
i}N_{j}(B_{ji}J_{\nu}+k_{ji}n+G_{ji}).\end{split}$ (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_{\nu}$, 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_{\nu}$ is achieved. In practice however, for species
other than H2O 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._
$\tau_{\nu}\ll 1$). In that case, the stimulated emission and absorption terms
($B_{ji}J_{\nu}$) 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}^{\prime}}=k_{ij}$ (in cm-3 s-1) 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 3: CO–CO state-to-state collision rate coefficients $k_{J_{1}\to J_{1}^{\prime}}$ as a function of $T_{kin}$ |  | $T_{kin}$
---|---|---
$J_{1}$ | $J_{1}^{\prime}$ | 5 K | 10 K | 20 K | 30 K
1 | 0 | 2.69e-11 | 3.61e-11 | 4.01e-11 | 4.13e-11
2 | 0 | 2.73e-11 | 3.03e-11 | 2.97e-11 | 2.81e-11
2 | 1 | 5.56e-11 | 6.23e-11 | 6.69e-11 | 6.70e-11
3 | 0 | 1.81e-11 | 2.02e-11 | 1.96e-11 | 1.81e-11
3 | 1 | 4.73e-11 | 5.14e-11 | 5.40e-11 | 5.30e-11
3 | 2 | 8.19e-11 | 7.88e-11 | 7.67e-11 | 7.50e-11
4 | 0 | 1.42e-11 | 1.49e-11 | 1.43e-11 | 1.30e-11
4 | 1 | 3.55e-11 | 3.74e-11 | 3.73e-11 | 3.49e-11
4 | 2 | 6.01e-11 | 6.22e-11 | 6.28e-11 | 6.10e-11
4 | 3 | 5.78e-11 | 6.66e-11 | 7.01e-11 | 6.98e-11
5 | 0 | 9.94e-12 | 9.69e-12 | 8.43e-12 | 7.43e-12
5 | 1 | 2.82e-11 | 2.95e-11 | 2.77e-11 | 2.54e-11
5 | 2 | 4.17e-11 | 4.30e-11 | 4.13e-11 | 3.84e-11
5 | 3 | 5.98e-11 | 6.46e-11 | 6.36e-11 | 6.06e-11
5 | 4 | 5.59e-11 | 6.40e-11 | 6.74e-11 | 6.64e-11

##  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 (private
communication). Their model has been used to analyze mm/sub-mm 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({\rm CO})=5\times 10^{28}$ s-1, a constant kinetic
temperature $T_{kin}=20$ K and outflow velocity $v_{out}=0.5$ km s-1, at a
heliocentric distance $r_{H}=2.8$ au. The results of the Biver et al. 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).

As a test of the SUBLIME raytracing 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$ K km s-1, compared with 0.36 K km s-1 from Biver et al’s
model, and for $J=3-2$, we have $\int{T_{R}dv}=0.64$ K km s-1, compared with
0.63 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 style) the level populations derived
from the same SUBLIME model, but using CO–H2 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–H2, with a mean
ratio between the new and H2-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–H2 rates, however, is relatively small.

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({\rm CO})=5\times 10^{28}$ s-1, $T_{kin}=20$ K, $v_{out}=0.5$ km s-1 and
$r_{H}=2.8$ au. Solid curves are using the time-dependent version of SUBLIME,
with CO–CO collision rates from Section 4.1; dotted curves are assuming the
CO–CO collision rates are the same as the CO–H2 rates from Yang et al. (2010);
dashed curves are using the model of Biver et al. (2018).

##  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-01-14 and 2018-01-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 16: Azimuthally averaged CO $3-2$ spectra, based on JCMT HARP jiggle
maps (combined data from 2018-01-14 and 2018-01-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.

## References

  * A’Hearn et al. (2011) A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396, doi: 10.1126/science.1204054
  * Altwegg et al. (2019) Altwegg, K., Balsiger, H., & Fuselier, S. A. 2019, ARA&A, 57, 113, doi: 10.1146/annurev-astro-091918-104409
  * Bensch & Bergin (2004) Bensch, F., & Bergin, E. A. 2004, ApJ, 615, 531, doi: 10.1086/424439
  * Biver (1997) Biver, N. 1997, PhD Thesis, Univ. Paris 7-Didero
  * Biver et al. (1997) Biver, N., Bockelée-Morvan, D., Colom, P., et al. 1997, Science, 275, 1915, doi: 10.1126/science.275.5308.1915
  * Biver et al. (1999) Biver, N., Bockelée-Morvan, D., Crovisier, J., et al. 1999, AJ, 118, 1850, doi: 10.1086/301033
  * Biver et al. (2016) Biver, N., Moreno, R., Bockelée-Morvan, D., et al. 2016, A&A, 589, A78, doi: 10.1051/0004-6361/201528041
  * Biver et al. (2018) Biver, N., Bockelée-Morvan, D., Paubert, G., et al. 2018, A&A, 619, A127, doi: 10.1051/0004-6361/201833449
  * Biver et al. (2019) Biver, N., Bockelée-Morvan, D., Hofstadter, M., et al. 2019, A&A, 630, A19, doi: 10.1051/0004-6361/201834960
  * Bockelee-Morvan (1987) Bockelee-Morvan, D. 1987, A&A, 181, 169
  * Bockelée-Morvan et al. (2004a) Bockelée-Morvan, D., Crovisier, J., Mumma, M. J., & Weaver, H. A. 2004a, Comets II, ed. G. W. Kronk, 391
  * Bockelée-Morvan et al. (2004b) Bockelée-Morvan, D., Biver, N., Colom, P., et al. 2004b, Icarus, 167, 113, doi: 10.1016/j.icarus.2003.04.001
  * Bockelée-Morvan et al. (2010) Bockelée-Morvan, D., Hartogh, P., Crovisier, J., et al. 2010, A&A, 518, L149, doi: 10.1051/0004-6361/201014655
  * Bockelée-Morvan et al. (2012) Bockelée-Morvan, D., Biver, N., Swinyard, B., et al. 2012, A&A, 544, L15, doi: 10.1051/0004-6361/201219744
  * Bøgelund & Hogerheijde (2017) Bøgelund, E. G., & Hogerheijde, M. R. 2017, A&A, 604, A131, doi: 10.1051/0004-6361/201629197
  * Boissier et al. (2007) Boissier, J., Bockelée-Morvan, D., Biver, N., et al. 2007, A&A, 475, 1131, doi: 10.1051/0004-6361:20078380
  * Boogert et al. (2015) Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541, doi: 10.1146/annurev-astro-082214-122348
  * Brinch & Hogerheijde (2010) Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25, doi: 10.1051/0004-6361/201015333
  * Buckle et al. (2009) Buckle, J. V., Hills, R. E., Smith, H., et al. 2009, MNRAS, 399, 1026, doi: 10.1111/j.1365-2966.2009.15347.x
  * Capaccioni et al. (2015) Capaccioni, F., Coradini, A., Filacchione, G., et al. 2015, Science, 347, aaa0628, doi: 10.1126/science.aaa0628
  * Cochran et al. (2012) Cochran, A. L., Barker, E. S., & Gray, C. L. 2012, Icarus, 218, 144, doi: 10.1016/j.icarus.2011.12.010
  * Cochran & McKay (2018) Cochran, A. L., & McKay, A. J. 2018, ApJ, 854, L10, doi: 10.3847/2041-8213/aaab57
  * Collings et al. (2004) Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133, doi: 10.1111/j.1365-2966.2004.08272.x
  * Cordiner et al. (2014) Cordiner, M. A., Remijan, A. J., Boissier, J., et al. 2014, ApJL, 792, L2, doi: 10.1088/2041-8205/792/1/L2
  * Cordiner et al. (2017) Cordiner, M. A., Biver, N., Crovisier, J., et al. 2017, ApJ, 837, 177, doi: 10.3847/1538-4357/aa6211
  * Cordiner et al. (2019) Cordiner, M. A., Palmer, M. Y., de Val-Borro, M., et al. 2019, ApJL, 870, L26, doi: 10.3847/2041-8213/aafb05
  * Cordiner et al. (2020) Cordiner, M. A., Milam, S. N., Biver, N., et al. 2020, Nature Astronomy, 4, 861, doi: 10.1038/s41550-020-1087-2
  * Cottin & Fray (2008) Cottin, H., & Fray, N. 2008, Space Sci. Rev., 138, 179, doi: 10.1007/s11214-008-9399-z
  * Coulson et al. (2017) Coulson, I. M., Cordiner, M. A., Kuan, Y.-J., et al. 2017, AJ, 153, 169, doi: 10.3847/1538-3881/aa6440
  * Coulson et al. (2020) Coulson, I. M., Liu, F.-C., Cordiner, M. A., et al. 2020, AJ, 160, 182, doi: 10.3847/1538-3881/abafc0
  * Crifo et al. (1999) Crifo, J. F., Rodionov, A. V., & Bockelée-Morvan, D. 1999, Icarus, 138, 85, doi: 10.1006/icar.1999.6051
  * Crovisier (1987) Crovisier, J. 1987, A&A, 68, 223
  * Crovisier & Encrenaz (1983) Crovisier, J., & Encrenaz, T. 1983, A&A, 126, 170
  * Crovisier et al. (1997) Crovisier, J., Leech, K., Bockelee-Morvan, D., et al. 1997, Science, 275, 1904, doi: 10.1126/science.275.5308.1904
  * Davidsson et al. (2022) Davidsson, B. J. R., Samarasinha, N. H., Farnocchia, D., & Gutiérrez, P. J. 2022, MNRAS, 509, 3065, doi: 10.1093/mnras/stab3191
  * Davidsson et al. (2021) Davidsson, B. J. R., Birch, S., Blake, G. A., et al. 2021, Icarus, 354, 114004, doi: 10.1016/j.icarus.2020.114004
  * de Val-Borro et al. (2018) de Val-Borro, M., Milam, S. N., Cordiner, M. A., et al. 2018, The Astronomer’s Telegram, 11254
  * Delaunay (1934) Delaunay, B. 1934, Bulletin de l’Académie des Sciences de l’URSS, 6, 793
  * Dello Russo et al. (2016) Dello Russo, N., Kawakita, H., Vervack, R. J., & Weaver, H. A. 2016, Icarus, 278, 301, doi: 10.1016/j.icarus.2016.05.039
  * Festou et al. (2001) Festou, M. C., Gunnarsson, M., Rickman, H., Winnberg, A., & Tancredi, G. 2001, Icarus, 150, 140, doi: 10.1006/icar.2000.6553
  * Fougere et al. (2012) Fougere, N., Combi, M. R., Tenishev, V., et al. 2012, Icarus, 221, 174, doi: 10.1016/j.icarus.2012.07.019
  * Fougere et al. (2016) Fougere, N., Altwegg, K., Berthelier, J. J., et al. 2016, A&A, 588, A134, doi: 10.1051/0004-6361/201527889
  * Fuchs et al. (2009) Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, A&A, 505, 629, doi: 10.1051/0004-6361/200810784
  * Garcia-Berrios et al. (2020) Garcia-Berrios, E., Cordiner, M. A., Charnley, S. B., de Val-Borro, M., & Biver, N. 2020, in European Planetary Science Congress, EPSC2020–426
  * Garrod & Pauly (2011) Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15, doi: 10.1088/0004-637x/735/1/15
  * Gasc et al. (2017) Gasc, S., Altwegg, K., Balsiger, H., et al. 2017, MNRAS, 469, S108, doi: 10.1093/mnras/stx1412
  * Gordon et al. (2021) Gordon, I., Rothman, L., Hargreaves, R., et al. 2021, Journal of Quantitative Spectroscopy and Radiative Transfer, 107949, doi: https://doi.org/10.1016/j.jqsrt.2021.107949
  * Green (1975) Green, S. 1975, J. Chem. Phys., 62, 2271
  * Gunnarsson (2003) Gunnarsson, M. 2003, A&A, 398, 353, doi: 10.1051/0004-6361:20021610
  * Gunnarsson et al. (2008) Gunnarsson, M., Bockelée-Morvan, D., Biver, N., Crovisier, J., & Rickman, H. 2008, A&A, 484, 537, doi: 10.1051/0004-6361:20078069
  * Gunnarsson et al. (2002) Gunnarsson, M., Rickman, H., Festou, M. C., Winnberg, A., & Tancredi, G. 2002, Icarus, 157, 309, doi: 10.1006/icar.2002.6839
  * Gunnarsson et al. (2003) Gunnarsson, M., Bockelée-Morvan, D., Winnberg, A., et al. 2003, A&A, 402, 383, doi: 10.1051/0004-6361:20030178
  * Haser (1957) Haser, L. 1957, Bulletin de la Societe Royale des Sciences de Liege, 43, 740
  * Hindmarsh (2019) Hindmarsh, A. C. 2019, ODEPACK: Ordinary differential equation solver library.  http://ascl.net/1905.021
  * Hogerheijde et al. (2009) Hogerheijde, M. R., Qi, C., de Pater, I., et al. 2009, AJ, 137, 4837, doi: 10.1088/0004-6256/137/6/4837
  * Huebner & Mukherjee (2015) Huebner, W. F., & Mukherjee, J. 2015, Planet. Space Sci., 106, 11, doi: 10.1016/j.pss.2014.11.022
  * Hutson & Green (1994) Hutson, J. M., & Green, S. 1994
  * Ioppolo et al. (2011) Ioppolo, S., van Boheemen, Y., Cuppen, H. M., van Dishoeck, E. F., & Linnartz, H. 2011, MNRAS, 413, 2281, doi: 10.1111/j.1365-2966.2011.18306.x
  * Ip (1986) Ip, W. H. 1986, ApJ, 300, 456, doi: 10.1086/163821
  * Jaeger (2008) Jaeger, S. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 394, Astronomical Data Analysis Software and Systems XVII, ed. R. W. Argyle, P. S. Bunclark, & J. R. Lewis, 623
  * Kawakita et al. (2014) Kawakita, H., Dello Russo, N., Vervack, Ron, J., et al. 2014, ApJ, 788, 110, doi: 10.1088/0004-637X/788/2/110
  * Kłos & Lique (2018) Kłos, J., & Lique, F. 2018, in Cold Chemistry: Molecular Scattering and Reactivity Near Absolute Zero (The Royal Society of Chemistry), 46–91, doi: 10.1039/9781782626800-00046
  * Läuter et al. (2019) Läuter, M., Kramer, T., Rubin, M., & Altwegg, K. 2019, MNRAS, 483, 852, doi: 10.1093/mnras/sty3103
  * Lisse et al. (2021) Lisse, C. M., Young, L. A., Cruikshank, D. P., et al. 2021, Icarus, 356, 114072, doi: 10.1016/j.icarus.2020.114072
  * Loomis et al. (2018) Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, AJ, 155, 182, doi: 10.3847/1538-3881/aab604
  * Luspay-Kuti et al. (2015) Luspay-Kuti, A., Hässig, M., Fuselier, S. A., et al. 2015, A&A, 583, A4, doi: 10.1051/0004-6361/201526205
  * Marboeuf & Schmitt (2014) Marboeuf, U., & Schmitt, B. 2014, Icarus, 242, 225, doi: 10.1016/j.icarus.2014.07.001
  * Markwardt (2012) Markwardt, C. 2012, MPFIT: Robust non-linear least squares curve fitting.  http://ascl.net/1208.019
  * McKay et al. (2019) McKay, A. J., DiSanti, M. A., Kelley, M. S. P., et al. 2019, AJ, 158, 128, doi: 10.3847/1538-3881/ab32e4
  * Migliorini, A. et al. (2016) Migliorini, A., Piccioni, G., Capaccioni, F., et al. 2016, A&A, 589, A45, doi: 10.1051/0004-6361/201527661
  * Milam et al. (2004) Milam, S. N., Savage, C., Ziurys, L. M., & Wyckoff, S. 2004, ApJ, 615, 1054, doi: 10.1086/424701
  * Mumma & Charnley (2011) Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471, doi: 10.1146/annurev-astro-081309-130811
  * Mumma et al. (2011) Mumma, M. J., Bonev, B. P., Villanueva, G. L., et al. 2011, ApJL, 734, L7, doi: 10.1088/2041-8205/734/1/L7
  * Ndengué et al. (2015) Ndengué, S. A., Dawes, R., & Gatti, F. 2015, Journal of Physical Chemistry A, 119, 7712, doi: 10.1021/acs.jpca.5b01022
  * Öberg et al. (2008) Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2008, ApJ, 678, 1032, doi: 10.1086/533432
  * Ootsubo et al. (2012) Ootsubo, T., Kawakita, H., Hamada, S., et al. 2012, ApJ, 752, 15, doi: 10.1088/0004-637X/752/1/15
  * Paganini et al. (2010) Paganini, L., Villanueva, G. L., Lara, L. M., et al. 2010, ApJ, 715, 1258, doi: 10.1088/0004-637X/715/2/1258
  * Penteado et al. (2015) Penteado, E. M., Boogert, A. C. A., Pontoppidan, K. M., et al. 2015, MNRAS, 454, 531, doi: 10.1093/mnras/stv1987
  * Qasim et al. (2018) Qasim, D., Chuang, K. J., Fedoseev, G., et al. 2018, A&A, 612, A83, doi: 10.1051/0004-6361/201732355
  * Qasim et al. (2020) Qasim, D., Fedoseev, G., Chuang, K. J., et al. 2020, Nature Astronomy, 4, 781, doi: 10.1038/s41550-020-1054-y
  * Rabli & Flower (2010) Rabli, D., & Flower, D. R. 2010, MNRAS, 406, 95, doi: 10.1111/j.1365-2966.2010.16671.x
  * Roth et al. (2021a) Roth, N. X., Milam, S. N., Cordiner, M. A., et al. 2021a, ApJ, 921, 14, doi: 10.3847/1538-4357/ac0441
  * Roth et al. (2021b) —. 2021b, PSJ, 2, 55, doi: 10.3847/PSJ/abdd3a
  * Sun et al. (2020) Sun, Z.-F., van Hemert, M. C., Loreau, J., et al. 2020, Science, 369, 307, doi: 10.1126/science.aan2729
  * Surin et al. (2007) Surin, L. A., Fourzikov, D. N., Giesen, T. F., et al. 2007, Journal of Physical Chemistry A, 111, 12238, doi: 10.1021/jp0743471
  * Tenishev et al. (2008) Tenishev, V., Combi, M., & Davidsson, B. 2008, ApJ, 685, 659, doi: 10.1086/590376
  * van der Tak et al. (2007) 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, doi: 10.1051/0004-6361:20066820
  * Vissers et al. (2003) Vissers, G. W. M., Wormer, P. E. S., & van der Avoird, A. 2003, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 5, 4767, doi: 10.1039/B309468E
  * Wierzchos & Womack (2018) Wierzchos, K., & Womack, M. 2018, AJ, 156, 34, doi: 10.3847/1538-3881/aac6bc
  * Womack et al. (2017) Womack, M., Sarid, G., & Wierzchos, K. 2017, PASP, 129, 031001, doi: 10.1088/1538-3873/129/973/031001
  * Yang et al. (2010) Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062, doi: 10.1088/0004-637X/718/2/1062
  * Zakharov et al. (2007) Zakharov, V., Bockelée-Morvan, D., Biver, N., Crovisier, J., & Lecacheux, A. 2007, A&A, 473, 303, doi: 10.1051/0004-6361:20066715