© ESO 2018

1 Introduction

Comets are the most pristine remnants of the formation of the solar system 4.6 billion yr ago. Investigating the composition of cometary ices provides clues to the physical conditions and chemical processes at play in the primitive solar nebula. Comets may also have played a role in the delivery of water and organic material to the early Earth (see Hartogh et al. 2011, and references therein). The latest simulations of the evolution of the early solar system (Brasser & Morbidelli 2013; O’Brien et al. 2014) suggest a more complex scenario. On the one hand, ice-rich bodies formed beyond Jupiter may have been implanted in the outer asteroid belt, participating in the supply of water to the Earth, or, on the other hand, current comets coming from either the Oort Cloud or the scattered disk of the Kuiper belt may have formed in the same trans-Neptunian region, sampling the same diversity of formation conditions. Understanding the diversity in composition and isotopic ratios of the comet material is therefore essential for the assessment of such scenarios (Altwegg & Bockelée-Morvan 2003; Bockelée-Morvan et al. 2015).

Recent years have seen significant improvement in the sensitivity and spectral coverage of millimetre receivers enabling sensitive spectral surveys of cometary atmospheres and simultaneous observations of several molecules. We report here observations of a very peculiar comet, C/2016 R2 (PanSTARRS), with the 30-m telescope of the Institut de radioastronomie millimétrique (IRAM). This object is a long-period, dynamically old, Oort cloud comet that passed close to the Sun 21 600 yr ago and has an orbit inclination of 58°. It is very peculiar in the sense that in autumn 2017, while approaching the Sun (perihelion was on 9.6 May 2018 UT at 2.60 AU), it exhibited a deep blue coma and tail, due to the presence of strong CO⁺ lines in the optical spectrum and very little dust or other emission lines. In contrast, usual comets show a dust tail and a coma of neutral or yellowish colour, resulting from the scattering of solar radiation by dust and the emission of C₂ Swan bands, if any at this heliocentric distance, which are not seen here. It seems to belong to a category of comets of which we know only very few examples: C/1908 R1 (Morehouse) (de la Baume Pluvinel & Baldet 1911) or C/1961 R1 (Humason) (Greenstein 1962), none having been observed with modern astronomical facilities. As in those two comets that were however observed closer to the Sun, optical spectroscopy of comet C/2016 R2 revealed unusually strong N${}_2^{+}$ lines (Cochran & McKay 2018).

Here, we present millimetre observations of CO, CH₃ OH, H₂ CO, HCN, CO⁺ and upper limits on several other molecules obtained with the IRAM 30-m, complemented by observations of the OH radical with the Nançay radio telescope, and amateur observations of CO⁺, N${}_2^{+}$ and dust continuum. From the millimetre data, we derive constraints on the gas temperature and outgassing pattern. Production rates and their time evolution on 23 and 24 January, as well as abundances relative to CH₃ OH or CO are derived. We then compare the abundances in comet C/2016 R2 with those measured in eight other comets observed at a similar heliocentric distance.

2 Observations of comet C/2016 R2 (PanSTARRS)

2.1 Observations conducted at IRAM 30 m

Comet C/2016 R2 (PanSTARRS) was observed with the IRAM 30-m radio telescope on two consecutive evenings, 23.8 and 24.8 January 2018 UT, under good weather (mean precipitable water vapour of 1 mm). We targeted the CO(J = 2–1) line at 230538.000 MHz first, before the HCN(J = 3–2) line at 265886.434 MHz and the H₂ S(1₁₀−1₀₁) line at 168762.762 MHz (see CDMS catalogue; Müller et al. 2005, the line frequencies uncertainties (<1 kHz) being much lower than the finest sampling of the spectrometer used, 20 kHz), with other molecular lines, especially ofmethanol, in the band or the other side band. We used the EMIR 1 and 2 mm band receivers (Carter et al. 2012) with the fast Fourier transform spectrometer (FTS) set to a 200 kHz sampling covering ~ 2 ×8 GHz in two polarizations simultaneously.

A log ofthe observations is given in Table 1. We used the secondary mirror wobbler, with a 180′′ throw at a frequency of 0.5 Hz to cancel the sky background. The comet was tracked with JPL#14 orbital solution. Coarse mapping of the strong CO line shows that the peak intensity was shifted by 2–3′′ to the west northwest (WNW; Figs. 1 and 2). We adjusted the pointing offset to integrate closer to the peak of intensity. All pointing offsets mentioned henceforth in Tables 1 and 5 are relative to the position of this peak of intensity.

The comet immediately revealed itself as very different from most other comets observed so far. The CO(2–1) line is one order of magnitude stronger than expected (based on a correlation between visual magnitudes and CO production rates from Biver 2001). Meanwhile the narrow and blue-shifted line shape (Fig. 3) strongly resembles the one of distant comet 29P/Schwassmann-Wachmann 1 (Gunnarsson et al. 2008), whose activity is dominated by a continuous, asymmetric and large outgassing of CO (≈ 4 × 10²⁸ molec s⁻¹). In total, 15 methanol lines are clearly detected; two lines of formaldehyde are also detected and the HCN(3–2) line is marginal. The sum of the three lines of HNCO covered by the observations that are expected to be the strongest is close to 3σ. Lines of neutral species (Figs. 3–7) extend from −0.7 to +0.5 km s⁻¹ on the axis of Doppler velocities relative to the comet rest frame, which is due to the projection of the expansion velocity vector of the gas. In contrast, the two CO⁺ lines at 236062.553 and 235789.641 MHz show a very marginal feature of ~ 2 km s⁻¹ in width superimposed on a more pronounced broad red-shifted line of 30 km s⁻¹. No other species are detected.

Sample spectra are shown in Figs. 1–8, and the line intensities and derived production rates are given in Table A.1.

Fig. 1 Small map of the CO(2–1) emission in the coma of comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 24.85 January 2018. The velocity and intensity scales of the spectra are given in the top right. The arrow points towards the projected direction of the Sun (phase angle of 17°).

2.2 Observation of OH radical with the Nançay radio telescope

In order to assess the water production rate of the comet, the OH lines at 1665 and 1667 MHz were searched for with the Nançay radio telescope between 2 January and 31 March 2018 (observation circumstances for selected periods are presented in Table 2). The telescope tracked the comet around transit for about 1 h on average every two days.

Fig. 2 Colour-coded map of the CO(2–1) line integrated intensity of comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 24.85 January 2018. The (0,0) offset position corresponds to the JPL ephemeris position, with a ~ 2′′ pointing uncertainty. The arrow points towards the projected direction of the Sun (phase angle of 17°).

Fig. 3 Average on-nucleus (offset 1.7′′) spectrum of the CO(2–1) line at 230.538 GHz obtained with the IRAM 30 m telescope on 24.80 January 2018. The vertical scale is the main beam brightness temperature and the horizontalscale is the Doppler velocity in the comet rest frame. The simulated profile with QCO = 5 × 1028 molec s−1 in a 0–60° cone at vexp = 0.56 km s−1 and QCO = 5.5 × 1028 molec s−1 in a 60–120° cone at vexp = 0.50 km s−1 is superimposed in green dotted line.

2.3 Complementary optical observations

In support of these observations conducted with large facilities, amateur astronomers followed the activity of the comet (which exhibited a total visual magnitude between 10 and 11) with their own telescopes. We report here on results obtained by Etienne Bertrand (St Sordelin), Hubert Boussier (IAU station K21), and François Kugel (IAU station A77). The first two astronomers obtained moderate-resolution (λ∕δλ ~ 500) visible spectra (λ = 390–600 nm) of the comet, while the latter two also provided Afρ values of the comet (measured outside the range of CO⁺ lines). The data and logs of the observations are also accessible from the comet observation database hosted by the Paris Observatory(Noël 2018).

Hubert Boussier took spectra on 29.9 November and 18.9 and 22.9 December 2017. In this paper, we present an analysis of the spectrum from the latter of the three (Figs. 9 and 10) obtained between 22:01 and 22:42 UT with a 28-cm Schmidt-cassegrain telescope working at F/D = 7.0 equipped with a LISA spectrometer (Thizy & Cochard 2011) and Atik 314L CCD detector. The slit width and integration window was 2.4′′ × 14′′. Sky-glow and light pollution lines were removed using the signal in 2.4′′ ×33′′ windows at −38′′ and +72′′ along the slit. The effective spectrum resolution is 1 nm (5.1 Åper pixel). Etienne Bertrand obtained his spectrum (Fig. 10) on 22 February 2018 between 19:09 and 20:09 UT during a 1 h integration with a 20-cm Schmidt-cassegrain telescope working at F/D = 6.3 equipped with an Alpy600 spectrometer (Thizy & Cochard 2011) and Atik 414EX CCD detector. The slit width and integration window was 3.8′′ × 53′′. Sky-glow and light pollution lines are removed using the signal in 3.8′′ × 37′′ windows at ±177′′ along the slit. The spectral resolution is also λ∕δλ ~ 500, that is, on the order of 1 nm. The average apertures used for these observations (equivalent diameter of 8′′ and 14′′ for the 22 Dececember and 22 February observations, respectively) are similar to the IRAM 30 m beam, which is useful for comparing both types of observation.

The optical spectra are presented in Figs. 10 and 11. The correction for the system and sky transmission was done using a nearby A-type reference star and an out-of-atmosphere synthetic spectrum. Only a relative calibration was obtained (no absolute calibration of the fluxes in e.g. W m⁻²Å⁻¹). Due to poortransmission of the atmosphere and optical system; and limited sensitivity of the detector towards the blue end of the spectrum, especially for the LISA spectrometer, the noise and uncertainty both increase below 400–430 nm. The wavelength calibration was performed using a reference internal source.

The spectra are clearly dominated by the CO⁺ doublets contrarily to most optical spectra of comets, dominated by CN and C₂ Swan lines (e.g. Fig. 1 of Feldman et al. 2004). The N${}_2^{+}$ line at 391 nm is also detected in both spectra, but in the noisiest part of the spectrum; we assumed that there is no contribution from the nearby CN line at 388 nm, both because it was not detected by Cochran & McKay (2018) and because the abundance of HCN is very low (see following section).

Fig. 4 Average off nucleus spectra of the CO(2–1) line obtained with the IRAM 30 m telescope on 23.78–24.85 January 2018. Left panel: 12.1′′ tail-ward (Δ RA > 0), right panel: 7.9′′ sun-ward. Scales are as in Fig. 3. The simulated profiles with QCO = 5 × 1028 molec s−1 in a 0–60° cone at vexp = 0.56 km s−1 and QCO = 5.5 × 1028 molec s−1 in a 60–120° cone at vexp = 0.50 km s−1 are superimposed in green dotted line.

Fig. 5 Series of methanol lines around 242 GHz from comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 23.8 January 2018. The inferred rotational temperature is Trot (242GHz) = 18 ± 2 K. The vertical scale is the main beam brightness temperature and the horizontal scale is the Doppler velocity in the comet rest frame (lower axis) or the rest frequency (upper axis).

Fig. 6 Average on-nucleus spectra of the ortho H2CO(312–211) line at 225.698 GHz (left panel) and para H2CO(303–202) line at 218.222 GHz (right panel) in comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 23.81 and 24.78 January 2018, respectively. Scales as in Fig. 3.

Fig. 7 Average on-nucleus spectrum of the HCN (3–2) line at 265886.434 MHz obtained with the IRAM 30 m telescope between 23 and 24 January 2018. Scalesare as in Fig. 3. Relative intensities and positions of the hyperfine components relative to F = 3 − 2 line at 265886.434 MHz are indicated.

Fig. 8 Weighted average of the two strongest CO+ N = 2 − 1 lines observed in C/2016 R2 (PanSTARRS) on 24.8 January 2018. Each line has been divided by its statistical weight (9/15 and 5/15 for the F = 5∕2 − 3∕2 and F = 3∕2 − 1∕2 transitions, respectively) before averaging. Scales are as in Fig. 3. The line shows a narrow peak close to the zero velocity and a broader component mostly red-shifted indicative of antisunward acceleration of CO+ ions by the solar wind.

3 Analysis and results

3.1 IRAM data

Thanks to the simultaneous detection of two to six methanol lines and to the very high S/N obtained on the CO line, we were able to constrain the excitation conditions of the molecules in the coma and the outgassing pattern of the comet. This is important to derive accurate outgassing rates and relative abundances.

3.1.1 Gas temperature

Table 3 provides the derived rotational temperatures of methanol lines and implied gas temperature. A rotational diagram of the 165 GHz lines is presented in Fig. 12. All groups of lines (at 165, 213–220, 242, 252, and 254–266 GHz) provide rotational temperatures T_rot in very good agreement in the 18–24 K range. Due to the radiative decay of the rotational levels within the ground vibrational state, we expect colder rotational temperatures for the 242 GHz lines. A higher collisional rate would limit this decay and observations suggest that indeed the collision rate might be higher than modelled, possibly due to the outgassing in narrow jets resulting in higher local densities. Nevertheless, we adopted a gas temperature T_gas = 23 K to derive all the molecular production rates. We also note that if we use the “jet” part (velocity interval −0.8 to −0.2 km s⁻¹) of the 165 GHz lines observed at high spectral (40 kHz) resolution we find a temperature slightly higher by 3.2 K. For comparison, at similar heliocentric distances, derived values are T_gas = 30–40 K for Hale-Bopp (Biver et al. 2002), 18 K for C/2006 W3 (Christensen) (Bockelée-Morvan et al. 2010), and 16 K for C/2002 T7 (LINEAR).

We can also constrain the kinetic temperature T_gas from the CO(2–1) line profile. Assuming that the width of the narrow component of the CO lineis due to thermal broadening, its full width at half maximum (FWHM) of 0.237 ± 0.003 km s⁻¹ implies T_gas ≤ 34 ± 1 K. This is compatible with T_gas = 23 K derived from the rotational temperatures. Modelling of the CO line shape (following section) yields information on the actual dispersion of velocities.

Fig. 9 Visible spectrum of comet C/2016 R2 (PanSTARRS) obtained by H. Boussier on 22.93 December 2018 UT with a 0.28 m telescope + LISA spectrometer and Atik 414L. The wavelengths covered are from 385 to 655 nm (from left to right). Vertical lines are atmospheric lines, mostly due to light pollution, while the cometary lines (mostly CO+ doublets) are less extended vertically. Scattering of the light of the mercury vapour lamps (Hg-I lines at 404.7, 435.8, 546.1, and 577+579 nm) and high pressure sodium lamps (568 nm line and broad emission around 589 nm) are the dominating sky emissions.

Fig. 10 Black line: visible spectrum of comet C/2016 R2 (PanSTARRS) obtained by H. Boussier on 22.92 December 2017 UT (40 min integration) from MPC station K21. The spectrum was extracted from the central 10 (binned) pixels rows. Green filled lines represent the simulated CO+ spectrum with respective average g-factors. Blue filled lines show the N${}_2^{+}$ spectrum fora N${}_2^{+}$/CO+ column density ratio of 0.06. The vertical intensity scale has been normalized. The dotted line is the sum of N${}_2^{+}$ and CO+ synthetic spectra. The signal towards the red end (>500 nm) of the spectrum is likely under-corrected because the LISA spectrometer is not optimised to cover the full wavelength range and the comet was slightly out of focus in that part of the spectrum. As a consequence, the signal of the lines beyond ~510 nm is lower than expected and some lines are not detected. This part of the spectrum is not used in the analysis.

3.1.2 Outgassing pattern and expansion velocity

After determining the gas temperature and implied thermal broadening of the lines, we modelled the CO line profile to infer the gas expansion velocity and outgassing pattern. For T_gas = 20–23 K the CO line shape asymmetry is well explained assuming an expansion velocity v_exp = 0.56 ± 0.01 km s⁻¹ on the Earth side and considering that most of the outgassing is restricted to a cone with a half-opening angle of 60°. Figure 13 shows the CO production rate per unit solid angle (assuming symmetry around the Comet-Earth axis) as a function of the Earth-Comet-V gas or colatitude angle (ϕ) inferred from the line profile. For this profile inversion (i.e. conversion of each channel intensity into a production rate per solid angle, with the channel velocity v_i converted into a colatitude angle ϕ_i from the relationship v_i = −v_exp cos(ϕ_i)), we assumed that v_exp (0.56 km s⁻¹) does not depend on ϕ. In reality, we expect the expansion velocity to be lower at higher phase angle due to lower surface temperature. This is consistent with the smaller blue-shift (by 0.06 km s⁻¹) observed in spectra acquired on the anti-sunward side (Fig. 4).

In any case, line profile inversion with a lower velocity down to 0.3 km s⁻¹ still produces a cut-off of production around the colatitude angle ϕ = 120–130°. To approximately simulate the profile of gas production that was to be used to derive total production rates, we assumed constant production rates per solid angle q₁ and q₂ in the colatitude ranges ϕ₁ = 0–60° and ϕ₂ = 60–120°, respectively. The expansion velocity was set to 0.56 km s⁻¹ (ϕ₁ = 0–60°) and 0.50 km s⁻¹ (ϕ₂ = 60–120°), respectively. We adjusted q₁ and q₂ to account for the observed mean Doppler shift of the line. q₂ is about half q₁ and the total production rates in these two regions (Q_i = 2π∫_ϕq_idϕ) are in the ratio Q₁:Q₂ = 10:11. The profile of the total production rate per solid angle q is plotted in Figs. 13 and 14. The resulting line profiles are superimposed to the observed line profiles in Figs. 3 and 4.

Although the agreement between simulation and observation is relatively good for the central position (Fig. 3), small differences appear at offset positions (Fig. 4). This shows that the actual outgassing pattern deviates somewhat from the axisymmetric description done in this study. Indeed the line intensity spatial distribution (Figs. 1 and 2) points at some asymmetry in the outgassing.

A similar profile inversion was undertaken for the average of the four strongest methanol lines (J = 1–4) observed at 165 GHz. The result (Fig. 14) is similar to the one from CO, although with a lower S/N, suggesting that CH₃OH and CO have similar production patterns from the nucleus. However, the peak of the methanol lines appears at slightly more negative velocities, possibly linked to the uncertainty on the absolute line frequencies. Indeed the shift is ~ 50 kHz, which is on the order of the differences between available frequencies in 2016 and previously published values (Müller et al. 2005, and references therein). If we expect the same line shape for CH₃OH and CO, the inferred CH₃OH line frequencies are 165050.258 ± 0.016, 165061.153 ± 0.023, 165099.220 ± 0.017 and 165190.544 ± 0.010 MHz, for the 1₁–1₀E, 2₁ – 2₀E, 3₁ – 3₀E, and 4₁ – 4₀E lines, respectively.

Fig. 11 Black line: visible spectrum of comet C/2016 R2 (PanSTARRS) obtained by E. Bertrand on 22.80 February 2018 UT (1 h integration) from Vaux-sur-Mer (France). The spectrum was extracted from the central 50 pixels rows. Green filled lines show the simulated CO+ spectrum with respective average g-factors. Blue filled lines indicate the N${}_2^{+}$ spectrum fora N${}_2^{+}$/CO+ column density ratio of 0.08. The vertical intensity scale has been normalized. Dotted lines is the sum of N${}_2^{+}$ and CO+ synthetic spectra.

Fig. 12 Rotational diagram of the CH3OH lines at 165 GHz. The green open dots correspond to the expected values for Tgas =23 K, corresponding to Trot = 20.5 K. The red line is the fit to the observations (Trot = 19.2 ± 3.9 K).

Fig. 13 Inversion of the CO(2–1) line profile from Fig. 3, assuming a constant expansion velocity of 0.56 km s−1. The production rate per steradian is derived as a function of the colatitude angle (ϕ = 0° on the Earthside, 180° on the opposite), assuming symmetry along the Comet-Earth line. The dots correspond to the values derived for each spectral channel and the connected line with error bars is the running average taking into account thermal broadening. The grey shaded region corresponds to the production rate profile used to determine the total production rates, with vexp =0.56 and 0.50 km s−1 in the dark and light grey regions, respectively (see text).

Fig. 14 Inversion of the average of four CH3OH lines (J = 1–4) at 165 GHz. Scales are as in Fig. 13.

3.1.3 Molecular production rates

In summary, to analyse all spectra, we used a gas temperature of 23 K and the outgassing pattern described above. A simpler description of the outgassing pattern assuming hemispheric outgassing at v_exp = 0.56 km s⁻¹ yields production rates that are only 5% higher. Production rates are provided in Tables A.1 and 5 and plotted in Fig. 16. They are calculated as described in Biver et al. (2016) and previous papers. Results from the maps are provided in Table 5 which gives average intensities at several radial offsets from the peak of intensity determined from the maps (Figs. 1 and 2) for the strongest lines and the derived production rates. There is only a marginal trend of increasing production rate with offset for CO (approximately +20%), not significant for the other molecules, which could also be related to instrumental effects (error beam and distortion of the beam shape at higher elevations). For the three molecules, observations are compatible with most of the production coming from the nucleus.

3.2 Search for OH

In order to constrain the water production rate of the comet, the 18-cm OH lines were observed with the Nançay radio telescope. From 2 January to 26 March, 2018, 48 individual daily spectra of about 1 h integration were secured¹. The observing and reduction procedures were as explained in Crovisier et al. (2002). During the observing period, the OH maser inversion varied from − 0.31 to − 0.15. Neither the individual spectra nor their average show any significant signal. For the 2 January– 9 February average (corresponding to a mean maser inversion of − 0.31), the 3-sigma upper limit on the 1667 MHz line area, integrated over a 3 km s⁻¹ width, is 12 mJy km s⁻¹. According to the model of Crovisier et al. (2002), this corresponds to an upper limit of 1.1 × 10²⁸ molec s⁻¹ for the OH production rate, that is, 1.2 × 10²⁸ molec s⁻¹ for the water production rate. Observations and resulting upper limits are summarised in Table 2.

3.3 N${}_2^{+}$/CO⁺ from visible spectra

Following the work presented by Cochran & McKay (2018), we estimate here the N${}_2^{+}$/CO⁺ column density ratio from the visible spectra presented in Sect 2.3. We used line frequencies from Lofthus and Krupenie (1977), Kepa et al. (2004), and Magnani and A’Hearn (1986). Fluorescence efficiencies (L/N or g-factors) of N${}_2^{+}$ are from Lutz et al. (1993) and those for CO⁺ are from Magnani and A’Hearn (1986). We use those pertaining to the heliocentric velocity range −10 to +25 km s⁻¹ (see Sect. 4.3).

Synthetic fluorescence spectra were computed assuming a width of ~ 0.4 nm for each v′ −v″ band and an instrumental resolution of 1 nm. A N${}_2^{+}$/CO⁺ abundance ratio of 6–8% was used to obtain the simulation plotted in Figs. 10 and 11. Band areas were computed to estimate the relative intensities (Table 4). Due to non-perfect relative calibration across the visible band and also variable uncertainties due to removal of atmospheric lines (mostly light pollution lines of Hg-I and high pressure sodium (Na-I) broad emission above 540 nm, especially for the spectrum of 22.9 December, 2017, in Fig. 9), we took the weighted average of the column densities derived from each detected line. The band observed around 427 nm is the sum of the N${}_2^{+}$(0–1) band at 427.4 nm and the CO⁺(2–0)Π_3∕2 band at 427.2 nm (band head; Magnani and A’Hearn 1986). In order to evaluate the possible contribution of the N${}_2^{+}$(0–1) band, we computed the signal of the CO⁺(2–0)Π_3∕2 using the average of the column density derived from the nearby CO⁺(2–0)Π_1∕2 band at 424.9 nm and the column density derived from all CO⁺ bands. The 424.9 nm band should be less affected by relative calibration issues, but the average of all lines should average out local fluctuations of the baseline and uncertainties on the g-factors.

The derived N${}_2^{+}$/CO⁺ column density ratio is 0.09 ± 0.03 for 22.9 December 2017 and 0.06 ± 0.02 for 22.8 February 2018, giving an average of 7 ± 2%. If we only take into account lines below 460 nm to limit biases with wavelength, the values become 0.08 ± 0.02 and 0.07 ± 0.03, respectively, yielding an average of 8 ± 2%. Since ionization efficiencies are similar for the two molecules (${\beta }_{i,{\text{N}}_2}=3.52\times {10}^{-7}$s⁻¹ and β_i,CO = 3.80 × 10⁻⁷s⁻¹ at 1 AU for quiet Sun; Huebner et al. 1992), this translates to a N₂ /CO production rate ratio of 8 ± 2%, consistent with the value of 6% determined by Cochran & McKay (2018).

3.4 Dust production

Our analysis on the dust production of the comet is solely based on optical data, as the detection of continuum emission in the infrared or in the radio has not been reported so far. The optical spectra present only a weak continuum (Fig. 9). Images show a “CO⁺” blue coma, and do not reveal any dust tail (Cochran & McKay 2018; Noël 2018), suggesting a low dust production. From coma photometry, F. Kugel and H. Boussier (Noël 2018) report Afρ values from 530 to 830 cm, with a mean value of 670 cm for the period from 8 February to 8 March, 2018 (at a nearly constant phase angle of 21°). The Afρ parameter introduced by A’Hearn et al. (1984) is proportional to the dust loss rate multiplied by the cross section of the dust particles within the aperture; it can be used as a proxy for the dust production. According to A’Hearn et al. (1995), a value of 1000 cm corresponds to a dust production rate of about 1 ton per second. For comet C/2016 R2, this would suggest that the dust production is about one order of magnitude less than the total gaseous production, in mass (Afρ∕Q_gas ≈ 0.1 cm kg⁻¹ vs. 2–5 for other comets at 3 AU (Table 6) and $Af\rho /Q_{{\text{H}}_2\text{O}}=$ 0.5 cm kg⁻¹ on average from A’Hearn et al. (1995) for comets closer to the sun). The conversion of the Afρ parameter into dust production rate is highly dependent on assumptions made on the dust size distribution and maximum size. If released dust consists predominantly of large grains, then the dust production rate could be much higher than the estimation given above. These large grains could form a dust trail which is not seen in the images, including those acquired in mid-December 2017 when the Earth crossed the plane of the orbit of the comet. As discussed in the following section, the Afρ value for C/2016 R2 is comparatively much lower than values measured in comets showing comparable gaseous activity.

4 Discussion

The derived CO production rate (10.6 × 10²⁸ molec s⁻¹, equivalent to ~5 tons per second) corresponds to a large outgassing rate, only surpassed by comet Hale–Bopp and by 17P/Holmes during its massive outburst (Table 6). This was unexpected given the small brightness of the comet (total visual magnitude m₁ around 10–11). Indeed, from the correlation between visual magnitudes and CO production (Biver 2001) we obtain a CO production rate 20 times lower. This indicates a dust-poor coma, consistent with the very low Afρ values (670 ± 110 cm) measured by F. Kugel and H. Boussier during this period. At the same heliocentric distance for comet Hale–Bopp the Afρ was 100 times higher (Weiler et al. 2003; at a similar phase angle) for a CO production only twice higher. We note that the visual magnitude is difficult to compare to that of other comets as the optical brightness is dominated by CO⁺ “tail” emission lines.

4.1 Molecular abundances

Pending the assessment of the abundance of other possible major molecules (CO₂, O₂), for which no ionic emission has been reported so far, the major species in the coma of comet C/2016 R2 (PanSTARRS) at 2.8 AU from the Sun are CO and N₂ with an abundance ratio on the order of 100:8. Observations and searchesfor NH₃ or NH₂ and hydrocarbons are needed to complete the inventory, but these species are not expected to be as abundant as the previous species. Water is not detected with an upper limit H₂O/CO < 0.1. Comparing with comets observed at similar heliocentric distances (Table 6), methanol is relatively abundant relative to water, but deficient compared to CO. Hydrogen cyanide and sulfur species are strongly depleted. Figure 15 compares the abundances of CO and HCN relative to methanol in comets observed at 2.3–3.3 AU from the Sun. Molecular production rates measured in these comets are listed in Table 6.

Fig. 15 CO/CH3OH and HCN/CH3OH ratios in comets observed between 2.3 and 3.3 AU from the Sun.

4.2 Upper limits on minor species

The wide frequency coverage (Table 1) has enabled us to look for many other species previously detected in comets (Biver et al. 2015). The most significant upper limits are provided in Table A.1. For species where several lines of similar expected intensities are sampled, we provide the combined upper limit (weighted rms) derived from all lines. In general the upper limits in comparison to the reference molecule (HCN for CN-species, CH₃OH for CHO-species and H₂S for S-species) are not very constraining (Biver et al. 2015). The marginal detection of HNCO in C/2016 R2 requires confirmation as it corresponds to a HNCO/CH₃OH ratio that is between three and ten times higher than in other comets and a HNCO/HCN two orders of magnitude higher. H₂ S is clearly under abundant (<6% relative to CH₃OH) in comparison to other comets observed at ~3 AU from the Sun (Table 6: H₂ S/CH₃OH = 30–70%).

We also looked for the CO isotopologues ¹³CO and C¹⁷O, but they were not detected. The best S/N is obtained considering the blue-shifted part of the lines, and gives ¹²C/¹³C > 54 (82 ± 42 at 2 σ), which is compatible with the terrestrial value (90).

4.3 Ions observed with IRAM 30 m

The rotational lines of CO⁺ were observed for the first time in comet Hale–Bopp at the CSO (Lis et al. 1997) and the IRAM 30 m (Henry et al. 2001). The CO⁺ lines at 236 GHz are marginally detected, red-shifted and broader than the CO line, as would be expected for ions (cf Hale–Bopp spectra of CO⁺ and HCO⁺; Lis et al. 1997). The two strongest CO⁺(N = 2–1) transitions at 236062.553 and 235789.641 MHz show a marginal 3–2 σ peak of FWHM = 1 km s⁻¹ close to the zero velocity in the comet frame and a broader component of 25 km s⁻¹ in width between −5 and +30 km s⁻¹ more clearly detected (5 and 4σ, respectively). This asymmetry towards larger red-shifted velocities is due to the acceleration of ions by the solar wind in the anti-solar direction. Figure 8 shows the combined spectrum (average of the two lines). Column densities cannot be derived in a straightforward manner as the properties of this profile show evidence of the acceleration of CO⁺ ions in the beam due to interaction with the solar wind. A rough estimate for an expansion velocity of 10 km s⁻¹ (mean Doppler shift of the lines), and a rotational population at 23–200 K is 0.4–1.7 × 10¹²cm⁻², but we have no precise idea of the rotational population of CO⁺ ions which are in a collisionless regime. We assumed that CO⁺ comes from the photo-ionization of CO (scale-length on the order of 12 × 10⁶ km) which yields optically thin lines in any case. The marginal narrow component would yield a column density about five times lower. However, this profile can be used to obtain more appropriate g-factors to analyse the optical spectra taken in a similar aperture: given the heliocentric velocity of the comet and the small phase angle of the observations, this means that the heliocentric velocity of the CO⁺ ions contributing the most of the signal is in the −10 to +25 km s⁻¹ range. Hence, we used the corresponding L/N values from Magnani and A’Hearn (1986) to interpret the optical spectra.

HCO⁺ is not detected. It has been detected in comets with weaker productions of CO, but at closer distances to the Sun where water production and protonation of CO in the coma are much more efficient (Milam et al. 2004).

4.4 Variation with time

Since the comet exhibited a unexpectedly large outgassing of CO, with some optical images showing rapid changes in the CO⁺ structures, we looked for possible evidence of short-term variations due to either the rotation of the nucleus or a transient outburst phase. Figure 16 shows the production rates from Table A.1 and their evolution over the two days of observations. On this short time scale, generally the variations were less than 20% and possibly more related to pointing and calibration uncertainties. Indeed, at high elevations (above 60–70°) the beam efficiency of IRAM 30 m degrades. We made some modelled corrections of the beam efficiency but could not precisely track its variation. This may account for ~10% variations. No significant variations of the CO(2–1) line Doppler shift (− 0.28 ± 0.02 km s⁻¹) that could be correlated with a variation in the production pattern – and independently of any calibration issue – are observed either. So, on a timescale of several hours, we do not see any significant variation of the activity.

Fig. 16 Evolution of the production rates (and upper limits) measured in comet C/2016 R2 (PanSTARRS) during the time interval spanning 23 and 24 January, 2018.

5 Summary and conclusion

We performed a multi-wavelength (optical, millimetre, centimetre) compositional study of comet C/2016 R2 (PanSTARRS). This comet has the following characteristics.

• A very large CO production, on the order of 10²⁹ molec s⁻¹ (i.e. 5 tons s⁻¹) at r_h = 2.8 AU from the Sun, only a factor of two below that of C/1995 O1 (Hale–Bopp) at same r_h.

• Very low dust production, based on the Afρ which is lower than values measured in comets of similar gaseous activity at ~3 AU from Sun by more than a factor of 15.

• Unusual optical spectrum showing essentially CO⁺ and N${}_2^{+}$ lines.

• Abundant N₂ in the coma (N₂/CO ~ 0.08), with N₂ dominating the nitrogen budget.

• A strong depletion of H₂O, CH₃OH, H₂CO, HCN, H₂S relative to CO (by more than one order of magnitude) compared with other comets observed at a similar heliocentric distance.

• A depletion of HCN and sulfur species relative to methanol.

The origin of the huge production of CO of C/2016 R2 could be related to the size of its nucleus. The nucleus of comet Hale–Bopp has a radius estimated to R_n ~ 37 km (Altenhoff et al. 1999), making this comet one of the largest ever observed. Distant comet 29P/Schwassmann-Wachmann 1, which also shows a large CO production on its circular 6-AU orbit (~ 4 × 10²⁸ molec s⁻¹), has a radius of ~23 km (Stansberry et al. 2004). Unfortunately, measurements of the nucleus size of C/2016 R2 have not yet been reported. Given the pure CO ice sublimation rate at 2.8 AU Z_CO(2.8) ~ 1.6 × 10²² molec m⁻² s⁻¹ (Sekanina 1991), the outgassing rate of C/2016 R2 could be supplied by a pure CO ice object of 3 km in diameter. However, detailed thermal and structural modelling is necessary to provide valuable constraints on the nucleus size from the CO production rate.

The depletion of H₂O, CH₃OH, H₂ CO, HCN, H₂S relative to CO, together with the low dust production, suggest that the large CO activity reflects a CO-rich ice composition and large ice/dust ratio within the nucleus material of C/2016 R2. Indeed, due to low latent heat, sublimation of CO ice can release smaller particles and larger aggregates than the outgassing of H₂ O ice (Gundlach et al. 2015). Therefore, if C/2016 R2 had a dust/ice ratio comparable to that of other comets, significant dust production should have been observed. Here, we assume that the thermal properties of the nucleus surface and subsurface layers of C/2016 R2, controlled by the heat conductivity and porosity of the material, are overall comparable to other comets. This seems a reasonable assumption. Indeed, the illumination-driven CO outgassing indicated by the CO line profile (Sect. 3.1.2) shows that the CO sublimation front is close to the surface, consistent with a low-porosity material. Other comets observed at r_h > 3 AU from the Sun (Hale–Bopp and 29P) display a similar CO outgassing pattern (e.g. Gunnarsson et al. 2003). In comet Hale–Bopp, which had a higher Afρ/gas ratio and higher abundance of H₂O and HCN at 2.8 AU (Table 6), icy grains were also found in the coma (Lellouch et al. 1998). We can speculate that following the low dust abundance, icy grains are also deficient in the coma of this comet, which could explain the low abundance of water and HCN, which could be significantly produced from the sublimation of icy grains in other comets like Hale–Bopp at such heliocentric distances.

The other particularity of comet C/2016 R2 is the large abundance of N₂, dominating other N-bearing species. Measured N₂/HCN is < 0.006 for Hale-Bopp (derived from Cochran & McKay 2018, using HCN/CO = 0.01), and ~0.2 in 67P at r_h = 3 AU (from Rubin et al. 2015; Le Roy et al. 2015). For C/2016 R2, N₂ /HCN = 2000, that is, four orders of magnitude higher. Since there is no hint of NH₂ lines in optical spectra, and other N-bearing species (Table A.1) searched for at IRAM are less abundant than N₂ by at least two orders of magnitude, we can conclude that contrary to other comets, most of the nitrogen escaping the nucleus is in the form of N₂.

Comet C/2016 R2 shows a remarkably high N₂/CO ratio (~0.08), comparable to the largest values measured in comets so far. A detailed comparison with other comets is given in Cochran & McKay (2018). Values are, for example, N₂/CO < 6 × 10⁻⁵ for comet Hale–Bopp (Cochran et al. 2000), 5.7 × 10⁻³ for 67P (Rubin et al. 2015), 0.01 for 29P (from a marginal detection of N${}_2^{+}$, Korsun et al. 2008), and 0.06 in C/2002 VQ₉₄ (LINEAR) (Korsun et al. 2014).

The strong nitrogen deficiency in some comets, both in the ice and refractory phases, was revealed during the space investigation of 1P/Halley, and confirmed with the Rosetta mission (Rubin et al. 2015; Fray et al. 2017). Its interpretation remains elusive. According to Owen and Bar-Nun (1995), the trapping of N₂ by amorphous water ice in the cooling solar nebula was inefficient within Neptune’s orbit, resulting in the formation of planetesimals deficient in N₂. Laboratory experiments show that the trapped N₂/CO ratio is depleted by a factor of ~20 at 24 K, with respect to the gas phase value, the depletion factor being strongly dependent on temperature (Bar-Nun et al. 2007). With a protosolar ratio N/C = 0.29 and assuming that all C and N is in the form of CO and N₂, the N₂ /CO ratio in the solar nebula gas phase is 0.15, resulting in a trapped N₂/CO of a few times 10⁻³ at ~25 K, relativelyconsistent with the value measured for 67P (Rubin et al. 2015). If this interpretation is correct, the very low N₂ /CO values measured in, for example, comet Hale–Bopp would indicate a formation in warmer regions of the solar nebula.

Considering the gas-phase species detected in comet C/2016 R2, the N/C ratio for this comet is ~ 0.15, that is, close to the solar value of 0.29 ± 0.10 (Lodders et al. 2009). To explain this property, a possibility is that C/2016 R2 agglomerated from grains formed at an extremely low temperature, favouring the trapping of high quantities of N₂ both as trapped gas and in condensed form. However, the high CO content of C/2016 R2 ices together with the low dust/ice ratio (assuming our extrapolation from coma to nucleus composition properties is correct) may suggest another scenario. Models examining the thermal evolution of the relatively large Kuiper Belt object show that due to radiogenic heating, the most internal layers reach high temperatures (Prialnik et al. 2008; Sarid & Prialnik 2009); released gases migrate towards colder regions where they refreeze. The more volatile ice refreeze closer to the cold surface than the less volatile, so that the pristine dust/ice mixture becomes enriched in volatile ices such as CO and N₂ (Prialnik et al. 2008; Sarid & Prialnik 2009). Dynamical studies of the transneptunian population argue for a rich collisional history in the Kuiper Belt (Morbidelli & Rickman 2015, and references therein). In the second scenario, comet C/2016 R2 would be a fragment of the disruptive collision of a large Kuiper Belt object, with properties representative of volatile-enriched layers.

Comet C/2016 R2 is representative of a family of comets that we observe only rarely each century. Besides C/1908 R1 (Morehouse) and C/1961 R1 (Humason) (de la Baume Pluvinel & Baldet 1911; Greenstein 1962), other candidates are comets 29P/Schwassmann-Wachmann 1 and C/2002 VQ₉₄ (LINEAR) (Korsun et al. 2008) which both showed optical spectra dominated by strong emissions of CO⁺ and N${}_2^{+}$, characteristics of abundant CO and N₂ production and high N₂/CO ratio. The diversity of the dust/gas ratios seen in these comets (as judged from the Afρ values, which indicate that 29P and C/2002 VQ₉₄ are dust-rich in comparison to the other comets in the sample) may favour the second scenario in which these comets are collisional fragments of differentiated transneptunian objects.

Acknowledgements

IRAM observations were conducted under the target of opportunity proposal D06-17 and we gratefully acknowledge the support from the IRAM director for awarding us discretionary time and the IRAM staff for their support and for scheduling the observations at short notice. Observations of comet C/2014 S2 (PanSTARRS) were in part made during the eighth IRAM summer school in September 2015, with contributions from P. Gratier, E. Garcia Garcia, A. Khudchenko, I. Kushniruk, T. Michiyama, A. Petriella, F. Ruppin, Y. Shoham and P. Torne. The data were reduced and analysed thanks to the use of the GILDAS, class software (http://www.iram.fr/IRAMFR/GILDAS). The Nançay Radio Observatory is the Unité scientifique de Nançay of the Observatoire de Paris, associated as Unité de service et de recherche (USR) No.704 to the French Centre national de la recherche scientifique (CNRS). The Nançay Observatory also gratefully acknowledges the financial support of the Conseil régional of the Région Centre in France. This research has been supported by the Programme national de planétologie de l’Institut des sciences de l’univers (INSU).

Note added in proof. CO was also detected with the 10 m telescope of the Arizona Radio Observatory with a limit on HCN, in agreement with the low abundance we observed (Wierzchos & Womack 2018).

Appendix A: Additional table

References

All Tables

All Figures

	Fig. 1 Small map of the CO(2–1) emission in the coma of comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 24.85 January 2018. The velocity and intensity scales of the spectra are given in the top right. The arrow points towards the projected direction of the Sun (phase angle of 17°).
In the text

	Fig. 2 Colour-coded map of the CO(2–1) line integrated intensity of comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 24.85 January 2018. The (0,0) offset position corresponds to the JPL ephemeris position, with a ~ 2′′ pointing uncertainty. The arrow points towards the projected direction of the Sun (phase angle of 17°).
In the text

Fig. 3 Average on-nucleus (offset 1.7′′) spectrum of the CO(2–1) line at 230.538 GHz obtained with the IRAM 30 m telescope on 24.80 January 2018. The vertical scale is the main beam brightness temperature and the horizontalscale is the Doppler velocity in the comet rest frame. The simulated profile with QCO = 5 × 1028 molec s−1 in a 0–60° cone at vexp = 0.56 km s−1 and QCO = 5.5 × 1028 molec s−1 in a 60–120° cone at vexp = 0.50 km s−1 is superimposed in green dotted line.

In the text

Fig. 4 Average off nucleus spectra of the CO(2–1) line obtained with the IRAM 30 m telescope on 23.78–24.85 January 2018. Left panel: 12.1′′ tail-ward (Δ RA > 0), right panel: 7.9′′ sun-ward. Scales are as in Fig. 3. The simulated profiles with QCO = 5 × 1028 molec s−1 in a 0–60° cone at vexp = 0.56 km s−1 and QCO = 5.5 × 1028 molec s−1 in a 60–120° cone at vexp = 0.50 km s−1 are superimposed in green dotted line.

In the text

	Fig. 5 Series of methanol lines around 242 GHz from comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 23.8 January 2018. The inferred rotational temperature is Trot (242GHz) = 18 ± 2 K. The vertical scale is the main beam brightness temperature and the horizontal scale is the Doppler velocity in the comet rest frame (lower axis) or the rest frequency (upper axis).
In the text

	Fig. 6 Average on-nucleus spectra of the ortho H2CO(312–211) line at 225.698 GHz (left panel) and para H2CO(303–202) line at 218.222 GHz (right panel) in comet C/2016 R2 (PanSTARRS) obtained with the IRAM 30 m telescope on 23.81 and 24.78 January 2018, respectively. Scales as in Fig. 3.
In the text

	Fig. 7 Average on-nucleus spectrum of the HCN (3–2) line at 265886.434 MHz obtained with the IRAM 30 m telescope between 23 and 24 January 2018. Scalesare as in Fig. 3. Relative intensities and positions of the hyperfine components relative to F = 3 − 2 line at 265886.434 MHz are indicated.
In the text

Fig. 8 Weighted average of the two strongest CO+ N = 2 − 1 lines observed in C/2016 R2 (PanSTARRS) on 24.8 January 2018. Each line has been divided by its statistical weight (9/15 and 5/15 for the F = 5∕2 − 3∕2 and F = 3∕2 − 1∕2 transitions, respectively) before averaging. Scales are as in Fig. 3. The line shows a narrow peak close to the zero velocity and a broader component mostly red-shifted indicative of antisunward acceleration of CO+ ions by the solar wind.

In the text

Fig. 9 Visible spectrum of comet C/2016 R2 (PanSTARRS) obtained by H. Boussier on 22.93 December 2018 UT with a 0.28 m telescope + LISA spectrometer and Atik 414L. The wavelengths covered are from 385 to 655 nm (from left to right). Vertical lines are atmospheric lines, mostly due to light pollution, while the cometary lines (mostly CO+ doublets) are less extended vertically. Scattering of the light of the mercury vapour lamps (Hg-I lines at 404.7, 435.8, 546.1, and 577+579 nm) and high pressure sodium lamps (568 nm line and broad emission around 589 nm) are the dominating sky emissions.

In the text

Fig. 10 Black line: visible spectrum of comet C/2016 R2 (PanSTARRS) obtained by H. Boussier on 22.92 December 2017 UT (40 min integration) from MPC station K21. The spectrum was extracted from the central 10 (binned) pixels rows. Green filled lines represent the simulated CO+ spectrum with respective average g-factors. Blue filled lines show the N${}_2^{+}$ spectrum fora N${}_2^{+}$/CO+ column density ratio of 0.06. The vertical intensity scale has been normalized. The dotted line is the sum of N${}_2^{+}$ and CO+ synthetic spectra. The signal towards the red end (>500 nm) of the spectrum is likely under-corrected because the LISA spectrometer is not optimised to cover the full wavelength range and the comet was slightly out of focus in that part of the spectrum. As a consequence, the signal of the lines beyond ~510 nm is lower than expected and some lines are not detected. This part of the spectrum is not used in the analysis.

In the text

Fig. 11 Black line: visible spectrum of comet C/2016 R2 (PanSTARRS) obtained by E. Bertrand on 22.80 February 2018 UT (1 h integration) from Vaux-sur-Mer (France). The spectrum was extracted from the central 50 pixels rows. Green filled lines show the simulated CO+ spectrum with respective average g-factors. Blue filled lines indicate the N${}_2^{+}$ spectrum fora N${}_2^{+}$/CO+ column density ratio of 0.08. The vertical intensity scale has been normalized. Dotted lines is the sum of N${}_2^{+}$ and CO+ synthetic spectra.

In the text

	Fig. 12 Rotational diagram of the CH3OH lines at 165 GHz. The green open dots correspond to the expected values for Tgas =23 K, corresponding to Trot = 20.5 K. The red line is the fit to the observations (Trot = 19.2 ± 3.9 K).
In the text

Fig. 13 Inversion of the CO(2–1) line profile from Fig. 3, assuming a constant expansion velocity of 0.56 km s−1. The production rate per steradian is derived as a function of the colatitude angle (ϕ = 0° on the Earthside, 180° on the opposite), assuming symmetry along the Comet-Earth line. The dots correspond to the values derived for each spectral channel and the connected line with error bars is the running average taking into account thermal broadening. The grey shaded region corresponds to the production rate profile used to determine the total production rates, with vexp =0.56 and 0.50 km s−1 in the dark and light grey regions, respectively (see text).

In the text

	Fig. 14 Inversion of the average of four CH3OH lines (J = 1–4) at 165 GHz. Scales are as in Fig. 13.
In the text

	Fig. 15 CO/CH3OH and HCN/CH3OH ratios in comets observed between 2.3 and 3.3 AU from the Sun.
In the text

	Fig. 16 Evolution of the production rates (and upper limits) measured in comet C/2016 R2 (PanSTARRS) during the time interval spanning 23 and 24 January, 2018.
In the text