# The N2 Production Rate in Comet C/2016 R2 (PanSTARRS)

S. E. Anderson,1 P. Rousselot,1 B. Noyelles,1 C. Opitom,2 E. Jehin,3 D.
Hutsemékers,3 and J. Manfroid3
1 Institut UTINAM UMR 6213 / CNRS, Univ. Bourgogne Franche-Comté, OSU THETA,
BP 1615, 25010 Besançon Cedex, France
2 Institute for Astronomy, Univ. of Edinburgh, Royal Observatory, Edinburgh
EH9 3HJ, UK
3 STAR Institute, Université de Liège, Allée du 6 Août 19c, 4000 Liège,
Belgium E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Observations of comet C/2016 R2 (PanSTARRS) have revealed exceptionally bright
emission bands of N${}_{2}^{+}$, the strongest ever observed in a comet
spectrum. Alternatively, it appears to be poor in CN compared to other comets,
and remarkably depleted in H2O. Here we quantify the N2 production rate from
N${}_{2}^{+}$ emission lines using the Haser model. We derived effective
parent and daughter scalelengths for N2 producing N${}_{2}^{+}$. This is the
first direct measurement of such parameters. Using a revised fluorescence
efficiency for N${}_{2}^{+}$, the resulting production rate of molecular
nitrogen is inferred to be Q(N2) $\sim$ 1 × 1028 molecules.s-1 on average for
11, 12 and 13 Feb. 2018, the highest for any known comet. Based on a CO
production rate of Q(CO) $\sim$ 1.1 × 1029 molecules.s-1, we find
Q(N2)/Q(CO)$\sim$0.09, which is consistent with the N${}_{2}^{+}$/CO+ ratio
derived from the observed intensities of N${}_{2}^{+}$ and CO+ emission lines.
We also measure significant variations in this production rate between our
three observing nights, with Q(N2) varying by plus or minus 20% according to
the average value.

###### keywords:

comets: general – comets:individual: C/2016 R2 (PanSTARRS)– molecular data

††pubyear: 2015††pagerange: The N2 Production Rate in Comet C/2016 R2 (PanSTARRS)–The N2 Production Rate in Comet C/2016 R2 (PanSTARRS) Table 1: Observing circumstances of comet C/2016 R2 (PanSTARRS) with VLT/UVES. Table from Opitom et al. (2019). Date |  $r_{h}$ (au) |  $\dot{r}_{h}$ (km.s-1) |  $\Delta$ (AU) |  $\dot{\Delta}$ (km.s-1) | Exposure time (s) | UVES Setup | UVES slit
---|---|---|---|---|---|---|---
2018-02-11T00:27:07.326 | 2.76 | -6.09 | 2.40 | 19.7 | 4800 | DIC1-390+580 | 0.44”$\times$8” - 0.44”$\times$12”
2018-02-13T00:46:23.196 | 2.76 | -5.97 | 2.43 | 19.9 | 4800 | DIC1-390+580 | 0.44”$\times$8” - 0.44”$\times$12”
2018-02-14T00:47:40.759 | 2.75 | -5.91 | 2.44 | 20.1 | 4800 | DIC1-390+580 | 0.44”$\times$8” - 0.44”$\times$12”

##  1 Introduction

Comets are among the most pristine relics of the formation of the Solar
System, having agglomerated from different icy grains and dust particles
leftover from the planetary formation process and undergone little alteration
since. As comets approach the Sun, the sublimation of their ices creates a
large gaseous coma, which, along with the interaction with solar radiation,
leads to various spectroscopic emissions. This allows us to investigate the
composition of cometary ices and determine the physico-chemical conditions of
their formation, thus providing further insight to the nature of the early
Solar System at the time and place these comets formed.

The long period comet C/2016 R2 (PanSTARRS), displayed atypical coma
morphology in optical images since its passage near perihelion ($\sim$ 3.0 au)
in December 2017, emitting strongly in blue optical wavelengths due to ion
emission dominating in the coma. Radio observations revealed that it was a CO-
rich comet (Wierzchos & Womack, 2018) and strongly depleted in water, with an
H2O/CO ratio of 0.0032 (McKay et al., 2019) with an upper limit of $<$ 0.1
(Biver et al., 2018). The spectrum was dominated by bands of CO+ as well as
N${}_{2}^{+}$, the latter of which was rarely seen in such abundance in comets
(Cochran & McKay, 2018; Opitom et al., 2019). It was also found to be both CN-
weak and dust-poor, with an $Af\rho$ of 500 cm (Opitom et al., 2019). This CO-
rich and water-poor composition, along with none of the usual neutrals seen in
most cometary spectra, makes C/2016 R2 a unique and intriguing specimen.

The apparent N2 deficiency in comets has long been a matter of great debate.
Despite Pluto and Triton – which also formed in the outer Solar System – both
exhibiting a N2-rich surface (Cruikshank et al., 1993; Owen et al., 1993;
Quirico et al., 1999; Merlin et al., 2018), very few ground-based facilities
have ever observed N${}_{2}^{+}$ in cometary spectra. This mainly concerns the
following comets: C/1908 R1 (Morehouse) (de La Baume Pluvinel & Baldet, 1911),
C/1961 R1 (Humason) (Greenstein, 1962), 1P/Halley (Wyckoff & Theobald, 1989;
Lutz et al., 1993), C/1987 P1 (Bradfield) (Lutz et al., 1993),
29P/Schwassmann-Wachmann 1 (Korsun et al., 2008; Ivanova et al., 2016; Ivanova
et al., 2018), and C/2002 VQ94 (LINEAR) (Korsun et al., 2008; Korsun et al.,
2014), and potentially C/2001 Q4 (NEAT) Feldman (2015). Only C/2002 VQ94
presents a N2/CO ratio as high as C/2016 R2.

The detection of N2 in comets using spectroscopic methods has been
challenging: as a diatomic, symmetrical molecule, N2 has no permanent dipole
moment, which results in the absence of pure rotational transitions. It emits
no radiation at millimeter-wavelengths, making the molecule invisible to
observations in that range. The electronic transitions are visible in the UV
through instruments placed outside of our atmosphere. The presence of N2 in
comet 67P/ Churyumov-Gerasimenko was not detected through any transition but
using the ROSINA mass spectrometer in-situ measurements (Rubin et al., 2015).
However, its daughter-species N${}_{2}^{+}$ is detectable in the optical range
through the bands of the first negative group (B${}^{2}\Sigma_{u}^{+}$ \-
X${}^{2}\Sigma_{g}^{+}$) with the (0,0) bandhead located at 3914 Å. Not only
have this ion’s spectral lines been observed in C/2016 R2, they are also the
brightest ever seen in a cometary spectrum (Cochran & McKay, 2018). The
quantity of N2 present is thus of significant importance. By measuring the
band intensity of the observed N${}_{2}^{+}$ in C/2016 R2’s spectrum, assuming
that solar resonance fluorescence is the only excitation source, the observed
emission fluxes have been used to calculate ionic ratios of N${}_{2}^{+}$/CO+
in the coma between $0.06\pm 0.01$ (Cochran & McKay, 2018; Opitom et al.,
2019) and $0.08\pm 0.01$ (Biver et al., 2018). This would be the same ratio
for N2/CO since ionization efficiencies of N2 and CO are similar
($\beta_{i,N_{2}}=3.52\times 10^{-7}s^{-1}$ and $\beta_{i,CO}=3.80\times
10^{-7}s^{-1}$ at 1 au for quiet Sun (Huebner et al., 1992)) as long as they
are fully photodissociated in the coma. This is larger than the best
measurement in comet 67P, with a N2/CO ratio of $\sim 2.87\times
10^{-2}$(Rubin et al., 2020), though these measurements were obtained much
closer to the nucleus.

Figure 1: Observational spectrum of C/2016 R2 (black) obtained at the edges of
the slit ($\sim 6500$km), overlapped with a CN model (red) and a N${}_{2}^{+}$
model (blue) from Rousselot et al. (2022). The $y$ axis (intensity scale) has
arbitrary units.

The origin of the unusual composition of C/2016 R2 is controversial. It may be
a fragment of a differentiated object as suggested by Biver et al. (2018),
similar to what was suggested for the CO-rich Interstellar comet 2I/Borisov
(Cordiner et al., 2020). There were also suggestions that the high CO content
in C/2016 R2 might be due to it being formed just beyond the CO snowline in
the protosolar nebula (Eistrup et al., 2019; Mousis et al., 2021).
Understanding the exact composition of this comet is essential to determine
the likelihood of each formation scenario, adding additional insight into the
composition of the early Solar System.

Using the N${}_{2}^{+}$ emission lines identified in the spectrum of C/2016
R2, we can then determine how much N${}_{2}^{+}$ is being produced through the
use of the Haser model. This model relates the intensity of the observation to
the number of molecules responsible for this emission. We can also determine
scalelengths related to the lifetime of the N2 molecules and N${}_{2}^{+}$
ions in the coma of C/2016 R2, which can be used for future comets presenting
high quantities of N${}_{2}^{+}$. This approach allows the determination of
how much N${}_{2}^{+}$, thus N2, is being produced by the active surface of
the comet’s nucleus.

Understanding the presence - or lack of - N2 in comets allows us to
investigate key features in the timeline of planetesimal formation. An object
such as C/2016 R2 is a unique opportunity to set a baseline for identifying
and exploring N2 in cometary spectra. We begin by describing the method of
obtaining the UVES (the high resolution spectrograph of the European Southern
Observatory Very Large Telescope) observations of C/2016 R2. Then, we outline
the Haser model and how we apply it to these observations in the Methods
section (Sec. 3). We start by applying the model to the molecule of CN, as it
is well known in the literature and provides a reliable benchmark for which to
test the quality of our observations and application of our methods. Finally,
we describe how N${}_{2}^{+}$ is identified in the spectra of cometary coma
and how we constrain its properties, before calculating its production rate.
We then discuss the implications of these results and conclude as to how this
can be applied to future N2-rich comets.

##  2 Observations

The observations of C/2016 R2 which were used in our analyses were obtained on
2018 February 11, 13, and 14 (one single exposure of 4800 s of integration per
night) with the Ultraviolet-Visual Echelle Spectrograph (UVES) (Dekker et al.,
2000) mounted on the 8.2 m UT2 telescope of the European Southern Observatory
Very Large Telescope (ESO VLT). All observations were made when C/2016 R2 was
near its perihelion distance of 2.6 au, at 2.76 and 2.75 au. We used a 0.44”
wide slit, providing a resolving power of R$\sim$80,000. The slit length was
8”, corresponding to about 14,500 km at the distance of the comet (geocentric
distance varying from 2.40 to 2.44 au). This is summarized in Table 1. No
N${}_{2}^{+}$ atmospheric lines were detected in the high resolution UVES
spectra.

The ESO UVES pipeline along with custom routines was used to perform the data
reduction; extraction; cosmic rays removal; and correct for the Doppler shift
due to the relative velocity of the comet with respect to the Earth. Both the
archived master response curve or the response curve determined from a
standard star observed close to the science spectrum were used to calibrate
the spectra in absolute flux using either with no significant differences.
More details can be found in the UVES ESO pipeline manual. This data
processing produced a 2D spectrum for each night of 11, 13, and 14 Feb,
calibrated in wavelength and absolute flux units. The spatial extent reached
30 rows each corresponding to 0.25 arcsec on the slit with different
cometocentric distances centred on the nucleus. A full description of the
observations and further data reduction can be found in Opitom et al. (2019).

##  3 Methods

Figure 2: Line of sight and variables used for the integration of the Haser
profile. The origin of the coordinates is the location at which $R$ is
minimum.

A usual method for estimating the absolute gas production rates in a cometary
coma is to fit the observed flux to a Haser profile (Haser, 1957). This gives
an analytical solution to the column density of parent and daughter species in
the coma along the line of sight. Assuming the medium optically thin, the
intensity of the coma observed at a projected distance $\rho$ from the nucleus
would be proportional to the number of molecules along the line of sight that
produce it, and we can infer how much of the considered species is present in
the coma at that specific cometocentric projected distance. These parameters
are displayed in Fig. 2.

The parent molecules are disintegrated by photodissociation by UV photons
following the law $n=n_{0}\mathrm{e}^{-t/\tau_{p}}$, with $n_{0}$ the number
of parent molecules at $t$ = 0 and $\tau_{p}$ the average lifetime of a parent
molecule. The daughter-species produced from the photodissociation of the
parent-molecules are seen to be ejected radially with a velocity of $v_{d}$
from the nucleus. In the case of $N_{2}$ it is an ionization process but it
obeys this hypothesis in the Haser model.

The photodissociation (or photoionization) lifetime, $\tau_{p}$, is an
important factor in the behaviour of the parent and daughter species. We use
the relationship $l_{p,d}=\tau_{p,d}v_{p,d}$ which determines the scalelength
of the parent or daughter species, respectively. If we call $Q=n(R_{0})\cdot
4\pi R_{0}^{2}v_{p}$ the production rate (in molecules.$s^{-1}$) of a given
parent species and $n(R)$ the volume density of the parent species at distance
$R$ from the nucleus (of radius $R_{0}$), we then integrate over the line of
sight:

$N_{\text{parent}}(\rho)=\mathrm{e}^{\frac{R_{0}}{l_{p}}}\cdot\frac{Q}{2\pi
v_{p}\rho}\int_{0}^{+\pi/2}\mathrm{e}^{\frac{-\rho}{l_{p}\cos{\alpha}}}\mathrm{d}\alpha$
(1)

where $N_{\text{parent}}(\rho)$ is the column density of the parent species in
molecules/cm2, and $\rho$ the projected cometocentric distance. This provides
a fair assessment of the number of parent molecules at a given projected
distance from the nucleus of the comet. Because N${}_{2}^{+}$ is a daughter
species we need to take into account the photoionization of its parent
molecule at each successive distance in order to determine the quantity of
daughter species being produced by successive disintegrations of the parent
molecule.

One of the limitations of the Haser model for computing the profile of a ionic
species such as N${}_{2}^{+}$ is that it considers neutral products (like OH,
C2, CN…) not affected by electric forces due to the solar wind particles (or
any other ionic species present in the coma). This work being restricted,
however, to the inner coma (projected distance inferior to 8000 km), this
limits the effects of electric forces on the dynamics of radical species and
the physical hypotheses of the Haser model are still valid. Raghuram et al.
(2020) computed that the dominant process for creating N${}_{2}^{+}$ ions in
C/2016 R2 is the photoionization process due to solar UV photons, the
ionization process due to electrons inside the coma representing only a few
percents relative to photoionization. The validity of using the Haser model is
also confirmed by the quality of our fit (see Fig. 6) and the symetry in the
observed intensity profiles (see 3).

For the daughter species N${}_{2}^{+}$ we have:

$N(\rho)=\frac{Q}{2\pi
v_{d}\rho}\frac{l_{d}}{l_{p}-l_{d}}\int_{0}^{+\pi/2}\mathrm{e}^{\frac{1}{l_{p}}(R_{0}-\frac{\rho}{\cos(\alpha)})}-\mathrm{e}^{\frac{1}{l_{d}}(R_{0}-\frac{\rho}{\cos(\alpha)})}\mathrm{d}\alpha$
(2)

We then integrate N’($\rho$) along the slit in order to obtain
$N_{\text{tot}}$, the total number of daughter species present in the slit for
a production rate of 1 molecules.$s^{-1}$. We apply this value to the
relationship:

$N_{\text{tot}}Q=\frac{4\pi F_{\text{tot}}\Delta^{2}}{g}$ (3)

with $Q$ the production rate in molecules.s of the parent molecule$s^{-1}$,
and $\Delta$ the geocentric distance of C/2016 R2 at the moment of
observation. $F_{\text{tot}}$ represents the total flux observed through the
slit (in erg.s-1.cm-2) and $g$ is the fluorescence efficiency (expressed here
in erg.s-1.molecule-1). We can compute the production rate $Q$ with:

$Q=\frac{4\pi F_{\text{tot}}\Delta^{2}}{gN_{\text{tot}}}$ (4)

This formula corresponds to the simple case where all the parent species
photodissociate into the same daughter species. If this is not the case the
branching ratio (i.e. the fraction of daughter species produced compared to
all the destruction processes) must be taken into account. The production rate
given by the formula above must be divided by this branching ratio.

Figure 3: Intensity profiles for identified CN (left) and identified
N${}_{2}^{+}$ (right) for the nights of 11 (blue), 13 (green), and 14 (red)
Feb over the cometocentric distances. The average value of all three nights is
shown in black. The variation of the intensity of N${}_{2}^{+}$ over all three
nights is striking, compared to the consistency of CN. This is possibly due to
the rapid variation of the ion tail, seen to change from hour to hour.

To compute the production $Q$ we need, consequently, to compute the parameter
$N_{\text{tot}}$ and $F_{\text{tot}}$ the total flux inside the slit. The
first quantity implies the knowledge of the Haser parameters mentioned above,
i.e. $v_{d}$, $l_{p}$ and $l_{d}$.

For $v_{d}$ most authors use a value of 1 km.s-1 for this parameter, because
it is a reasonable estimate for comets at $\sim$1 au (based, e.g. on the
results from space probes). Cochran & Schleicher (1993) discuss the
variability of this parameter both with heliocentric distance and production
rate. The production rate seems to influence this parameter only for very
active comets, which is not the case for C/2016 R2. The heliocentric distance
influences significantly the gas expansion radial velocity inside the coma.
This effect cannot be neglected in our case because of the large heliocentric
distance of C/2016 R2 at the time of our observations (2.76 au). Cochran &
Schleicher (1993) conclude their discussion by using the law $v_{d}=850\times
r_{h}^{-0.5}$ with $v_{d}$ expressed in km.s-1. Such a law provides
$v_{d}=0.511$ km.s-1. It is why we decided to use 0.5 km.s-1, which is also
the value used by Raghuram et al. (2021). It is important to keep in mind,
nevertheless, that such a law is based on observations of comets with coma
dominated by water molecules. In the case of C/2016 R2 the coma is dominated
by CO, which can have an influence on this parameter. This influence being
difficult to quantify the above mentioned value seems the best approximation
that we can use. For our discussion on the scalelengths and for a better
comparison with other works, that consider scalelengths varying as $r_{h}^{2}$
(i.e. that consider a constant velocity with heliocentric distance and only
the variation of lifetime with heliocentric distance) we will use the same
scaling in this paper.

Because no work has ever been published so far for $l_{p}$ and $l_{d}$ values
in N${}_{2}^{+}$ in a cometary coma, our method consisted of computing this
profile along the slit and fitting it with different sets of parameters to
find the best values. The process is described in Sec. 4.1.These parameters
were then used to compute N${}_{\text{tot}}$. For the fluorescence
efficiencies $g$, we used recently revised values (Rousselot et al., 2022).
The true size of the radius of the nucleus is unknown; a reasonable
approximation must be used. The following results are based on calculations
using a 5 km radius as this is the value consistently used when applying the
Haser model to comets of unknown size.

From the 2D spectra obtained with UVES, we extracted row by row along the slit
length 30 1D spectra for each of the three different observing nights. For all
of them we subtracted individually the solar continuum scattered by the dust
grains, by adjusting a solar spectrum (Kurucz et al., 1984) convolved with the
instrument response function in the region without cometary emission lines.

These 1D spectra were then used to compute the total flux $F_{tot}$ both for
N${}_{2}^{+}$ and CN, by summing their different emission lines (we describe
this process in Sections 4.1 and 4.2). Each 1D spectrum corresponding to a
different cometocentric distance, it was possible to compute radial profiles
and to test different fits based on a Haser profile with different sets of
parameters. The seeing was taken into account for computing these profiles, by
convolving the theoretical profile with a Gaussian having a FWHM similar to
the seeing value recorded during the observations (i.e. 0.95”, 1.1”, and 0.94”
respectively for the nights February 11, 13, 14).

In order to properly estimate the production rate of N${}_{2}^{+}$, we first
measured the production rate of CN as a benchmark to ensure our methods are
sound. CN has long been studied in comets, and is well referenced in the
literature. We first fitted the Haser model on the CN (0,0)
B${}^{2}\Sigma^{+}$ \- X${}^{2}\Sigma^{+}$ band as a test to know if
scalelengths can properly be determined from data at hand.

##  4 Results and Discussion

###  4.1 CN production rate

Figure 4: Confidence level contours defined from the $\chi^{2}$ distribution showing the range of parent and daughter scalelengths that fit the observed data for CN molecule (left) and N${}_{2}^{+}$ ion (right). Values displayed are the result of the $\chi^{2}$ fit $\times 10^{-16}$. The darkest region indicates the best fit. Multiple pairs of parent and daughter scalelengths provide almost identical fits to the observed data. In the case of N${}_{2}^{+}$, as the scalelengths are of the same magnitude, numerical errors emerge when approaching a pair of identical scalelengths as it results in a division by zero. Figure 5: Our best fit of the convolved Haser model (purple) on the observed flux of CN (black) along with models from scalelengths obtained from the literature (Table 2). The observed CN flux is averaged from the observations of 11 Feb, 13 Feb, and 14 Feb (one profile for both sides of the slit from the nucleus), then averaged again over cometocentric distance (dashed line). Table 2: Parent and daughter scalelengths $l_{p}$ and $l_{d}$ at 1 au for the CN molecule found in the literature. This provides an indication of what to expect to find in Comet C/2016 R2. The values of the daughter scalelengths are not proportional to the parent scalelengths. Our fits with these scalelengths are shown in figure 2. Reference |  $l_{p}$ [km] |  $l_{d}$ [km]
---|---|---
This work | 1.3 $\times 10^{4}$ | 2.8 $\times 10^{5}$
A’Hearn et al. (1995) | 1.3 $\times 10^{4}$ | 2.2 $\times 10^{5}$
Cochran et al. (1992) | 1.4 $\times 10^{4}$ | 3.0 $\times 10^{5}$
Randall (1992) | 1.5 $\times 10^{4}$ | 1.9 $\times 10^{5}$
Wyckoff et al. (1998) | 1.6 $\times 10^{4}$ | 3.3 $\times 10^{5}$
Newburn & Spinrad (1989) | 1.8 $\times 10^{4}$ | 4.2 $\times 10^{5}$
Womack et al. (1994) | 2.5 $\times 10^{4}$ | 1.9 $\times 10^{5}$
Fink et al. (1991) | 2.8 $\times 10^{4}$ | 3.2 $\times 10^{5}$

A rare feature of C/2016 R2 is the weakness of the CN (0,0)
B${}^{2}\Sigma^{+}$ \- X${}^{2}\Sigma^{+}$ emission lines at 3880 Å, usually
one of the strongest emission features observed in optical spectra of comets.
Here it is almost entirely dominated by the N${}_{2}^{+}$ of first negative
group (B${}^{2}\Sigma_{u}^{+}$ \- X${}^{2}\Sigma_{g}^{+}$) (0,0) bandhead at
3914 Å.

In order to identify the CN emission lines, we used a theoretical fluorescence
spectrum computed by interpolation at the right heliocentric velocity from a
spectrum calculated by Zucconi & Festou (1985). This model was convolved by
the instrument response function. The resulting synthetic CN model is shown in
figure 1.

The CN lines are then identified and summed along the pixels of the
spectrometer using the synthetic model as a mask. For each wavelength for
which the CN spectrum is non-null, we identified a potential match. We used
the high-resolution spectral atlas of cometary lines in 122P/de Vico (Cochran
& Cochran, 2002) in order to identify potential emission lines due to other
species than CN. In some cases, there is an overlap of several species, which
both could produce the observed intensity. In particular, CH (0-0) is intense
in C/2016 R2 (contrary to CN and C2), particularly around 3881 Å , around 3886
Å , and around 3892 Å . In the case where CN lines overlap with other species,
the intensity was not measured, in order to avoid contamination. As these CN
lines are remarkably weak compared to the intense overlap, these two lines
were eliminated with no effect to the CN $g$ factor or the total flux. The
selected lines were integrated in wavelength to compute their total flux
$F_{\text{tot}}$ so as to draw an intensity profile on both sides of the
comet’s nucleus, which is shown in Fig. 3.

We then averaged the flux of each line for the observations of the three
nights. The resulting averaged intensity profile is shown in black in Fig. 5.
The final intensity profile is a 1D profile starting from the nucleus of the
comet and sweeping outwards along the coma ending at a distance of $\rho\sim
6.5\times 10^{3}$ km. We first computed this profile separately for both side
of the slit (with respect to the nucleus, located in the middle of the slit)
and then averaged again these two 1D profiles into a single one. The total
flux measured for CN over the entire slit and averaged over the three nights
of observation was $2.1\times 10^{-15}$ erg.s-1.cm${-^{2}}$.

By using a $\chi^{2}$ test we were able to estimate the best fit of the Haser
model to the observed intensity profile and determine the scalelengths of both
the parent (HCN) and daughter (CN) species in the coma of C/2016 R2. The best
value for the scalelength of HCN was determined to be $l_{p}=1.3\times 10^{4}$
km while for CN radical it was determined to be $l_{d}=2.8\times 10^{5}$ km,
see Fig. 4 (values scaled to 1 au using an $r_{h}^{2}$ law). Our fit, as well
as other fits using scalelengths from the literature, are shown after the
convolution by the average seeing of the three nights in Fig. 5. Table 2 shows
that our scalelengths values are in good agreement with the literature. It
confirms that, despite a relatively small slit, we have a sufficient range of
cometocentric distances to model the intensity distribution of a daughter
product with a good accuracy. However, it should be noted that multiple pairs
of scalelengths could also fit the data.

The branching ratio for CN produced by HCN is nearly equal to one (Huebner et
al., 1992). The parent-species is most likely HCN as seen by the dust-poor
composition of C/2016 R2 (Opitom et al., 2019). Using the fluorescence factor
$g=3.52\times 10^{-2}$ photons.s-1.molecule-1 at 1 au provided by Schleicher
(2010), we estimate a production rate of $\sim 10^{25}$ mol.s-1 if we assume
the usual value of $v_{d}=1$ km.s-1). This latter values is about three times
the one published by Opitom et al. (2019), based on a similar value of
$v_{d}$. The authors of Opitom et al. (2019) have since reviewed their
computation about the overall CN flux measured in the slit and an error was
found. After correction, it agrees with our revised CN production rate. If we
assume $v_{d}=0.5~{}km.s^{-1}$, as it seems closer to the reality we obtain
Q(CN)$\sim 5\times 10^{24}$ mol.s-1.

###  4.2 N${}_{2}^{+}$ production rate

The intensity profiles can be used to get a determination of the scalelengths
both for the parent molecule and the daughter species. N${}_{2}^{+}$ is the
daughter ion produced by the photoionization of N2, so we apply the Haser
model. We use the same parameters and methodology as those in our fit of the
Haser profile to the CN intensity profile as it was shown to provide good
results. From the synthetic fluorescence spectrum computed by Rousselot et al.
(2022) we determined the line wavelengths of the N${}_{2}^{+}$ (0,0) band and
the corresponding lines in the UVES spectra. Here, we limit the identification
process to the 3885.5 Å to 3915.0 Å interval in order to extract only the
lines of the (0,0) band and to exclude the CN emission lines region. The Haser
profile is then fit to our observations using the same methods as presented in
Sec. 4.1. The resulting intensity profile is shown in Fig. 3.

We explore this interval with the $\chi^{2}$ test and find new scalelengths of
$l_{p}=2.8\times 10^{6}$ km and $l_{d}=3.8\times 10^{6}$ km (scaled to 1 au).
This fit is shown in figure 6 after being convolved by the average seeing of
the three nights. These values are within the expected range estimated from
the rate coefficients. However, at this scale, multiple pairs of scalelengths
could be selected for N${}_{2}^{+}$ with an almost equally good fit. The
difference between the fit of different scalelength pairs is very small and
the effect on the slope is almost negligible. We understand our intensity
profile presents a certain degree of noise, and even the slightest alteration
to our selection criteria would give a different best fit profile.

These scalelengths can be used to compute the N${}_{2}^{+}$ production rate.
When we take into account the branching ratio, only a small fraction of the N2
molecules is ionized to N${}_{2}^{+}$, most of them being photodissociated.
Huebner et al. (1992) provide rate coefficients for the photodissociation and
ionization of N2. From these data we computed that the branching ratio for
creating N${}_{2}^{+}$ ion is 0.34 for a quiet Sun (and 0.36 for an active Sun
but with a total photodestruction rate being about twice). Wyckoff & Wehinger
(1976) also provide an estimate of both the photodissociation rate and
photoionization rate, significantly lower than the one of Huebner et al.
(1992) and a branching ratio of 0.42. They admit, like Huebner et al. (1992),
that there is a large uncertainty on their values. The photodissociation rate,
especially, is rather uncertain because both of uncertainties on the cross
sections and the fact that a predissociation line has a similar wavelength to
the solar Lyman $\gamma$ line (highly variable with solar activity and that
could be Doppler shifted). We use an average of these two branching ratios,
i.e. 0.38.

The fluorescence factor $g$ estimated by Rousselot et al. (2022) for the (0,0)
band in the 3885.5 Å to 3915.0 Å interval is $5.41\times 10^{-3}$
photon.s-1.ion-1 at $r_{h}=2.76$ au. Using this value of $g$, the measured
total integrated flux $F_{\text{tot}}=1.0\times 10^{-14}$ erg.s-1.cm-2, a
branching ratio of 0.38, and a velocity $v_{d}$ of 0.5 km.s-1 at $r_{h}$ (see
Section 3) we estimate a N2 production rate of $\sim 1.0\times 10^{28}$
molecules.s-1.

Figure 6: Best fit of the Haser model (blue) to the observed flux of
N${}_{2}^{+}$ (black). The observed N${}_{2}^{+}$ flux is averaged from the
observations of 11 Feb, 13 Feb, and 14 Feb (a single profile for both sides of
the slit centred on the nucleus), then averaged again over cometocentric
distance (dashed line).

Our result can be compared to other estimates based on the ratio with the CO
production rate. A first result has been calculated by Wierzchos & Womack
(2018) as $Q(N_{2})=(2.8\pm 0.4)\times 10^{27}$ mol.s-1 by determining the N2
column density and production rate using the N2/CO abundance ratio calculated
from optical spectra and their own CO results, purposefully chosen so as to be
as close as possible to McKay et al. (2019), with Q(N2) = (4.8 $\pm$ 1.1)
$\times 10^{27}$ mol.s-1. While not said explicitly in Biver et al. (2018),
their N2 production rate can be inferred from their CO production rate and
N2/CO radio to be in the order of $Q(N_{2})\sim 8.5\times 10^{27}$ mol.s-1.
These production rates were all derived directly from fluxes by using the
relative ratio of N${}_{2}^{+}$ /CO+, as:

$\frac{N_{2}^{+}}{CO^{+}}=\frac{g_{CO^{+}}}{g_{N_{2}^{+}}}\frac{I_{N_{2}^{+}}}{I_{CO^{+}}}$
(5)

where $I$ is the intensity of the bands. The relative ratio was determined
using a $g$ factor of $7\times 10^{-2}$ photons.s-1.ion-1 from Lutz et al.
(1993) scaled to 1 au. We estimate these values with the updated $g$ factor
from Rousselot et al. 2022, using $g=4.90\times 10^{-2}$ photons.s-1.ion-1
(for the whole band). Prior measurements of the N2 production rate become
Q(N2) = 4.6 $\times 10^{27}$ mol.s-1 for Wierzchos & Womack (2018), Q(N2) =
8.0 $\times 10^{27}$ mol.s-1 for McKay et al. (2019), and Q(N2) = 1.4 $\times
10^{28}$ mol.s-1 for Biver et al. (2018), providing the upper and lower values
of the expected production rate. Our N2 production rate, consequently, is
right within the adjusted production rates found using relative ratios. It is,
nevertheless, the first direct estimate of this production rate, independent
from any assumption about the CO production rate.

In our estimate of this production rate the main sources of uncertainties come
from the determination of $l_{p}$ and the radial expansion velocity of
daughter products $v_{d}$. We also assume that the solar wind has little
influence on the N${}_{2}^{+}$ column density in the region close to the
nucleus (the radial profile being then well described by the Haser model). Our
model provides a value for $l_{p}$ with the same order of magnitude as the
estimates of N2 destruction rates provided by Huebner et al. (1992) and
Wyckoff & Wehinger (1976), but the difference between these two estimates is
large (a factor of three). A lower value of $l_{p}$ would reduce the
N${}_{2}^{+}$ production rate roughly in the same proportion, leading to a
significantly different value of the one computed from the N${}_{2}^{+}$/CO+
ratio and CO production rate. The relationship from these two parameters is
based on a similar ionization rate for N2 and CO. It must also be kept in mind
that photodestruction of N2 is strongly dependant on solar activity. Our
estimate of $l_{p}$ corresponds to a low solar activity and could be
significantly different for comets observed during maxima of solar activity.
Errors also arise from our estimation of the total flux $F_{\text{tot}}$: Our
selection criteria on the emission lines can be more restrictive or
permissive, depending on whether we identify a match where our model is non-
null, or try to restrict noise by setting a positive lower limit. However,
this only has an effect of $\pm 0.1$ erg.s-1.cm-2.

Another important result from our work concerns the variability of N2
production rate between observation dates. The production rate given above
corresponds to the average of our three observing nights but, as shown Figure
3, the observed intensity significantly changes from one night to another. If
we compute N2 production rates separately for these three nights with the same
parameters mentioned above we find, respectively for 11, 12 and 13 Feb.:
$1.0\times 10^{28},8.0\times 10^{27}$, and $1.2\times 10^{28}$. Such rapid
changes in the N2 line intensities are not an instrumental effect because CN
profiles do not reveal any similar change
($F_{\text{tot}}=2.31,2.01,2.01\times 10^{-15}$ erg/s/cm2 for $Q$(CN) = 5.4,
4.7, 4.7 $\times 10^{24}$ mol.s-1 on Feb 11, 13, and 14, respectively). The
origin of such rapid changes is unknown but likely related to inhomogeneities
of the nucleus, thought it may be due to changes in the solar wind, as we do
not see such such rapid changes for neutral species.

##  5 Summary and Conclusions

We investigated the CO-dominated and water-poor comet C/2016 R2 (PanSTARRS)
using the Haser model in order to determine the N2 production rate. The
intensity of the N${}_{2}^{+}$ lines is truly incomparable to what has been
seen in comets so far, confirming the rarity of this type of comet. We
determined the following characteristics:

  * $-$

Parent (HCN) and daughter (CN) scalelengths of $l_{p}=1.3\times 10^{4}$ km and
$l_{d}=2.8\times 10^{5}$ km respectively (scaled to 1 au using $r_{h}^{2}$),
consistent with the literature;

  * $-$

A CN production rate of $Q$(CN)$\sim 5\times 10^{24}$ molecules.s-1, lower
than usually observed in cometary spectra at similar heliocentric distances
but consistent with other estimates for this comet;

  * $-$

Parent (N2) and daughter (N${}_{2}^{+}$) scalelengths of $l_{p}=2.8\times
10^{6}$ km and $l_{d}=3.8\times 10^{6}$ km respectively (for 1 au, when using
a scaling with $r_{h}^{2}$). ;

  * $-$

A N2 production rate of $Q$(N2) $\sim 10^{28}$ molecules.s-1, exceptionally
high, within the values estimated by other teams from the ratio
N${}_{2}^{+}$/CO+;

  * $-$

When compared to a CO production rate of Q(CO) $\sim$ 1.1 × 1029
molecules.s-1, we find a N2/CO ratio of 0.09, which is consistent with
observed intensity ratios. This is the highest of such ratios observed for any
comet so far;

  * $-$

Some large variations of N2 production rates over the course of the three
observing nights, the N2 production rate given above being the average of
these values. The production rate computed with the same parameters varies
between $8.0\times 10^{27}$ and $1.2\times 10^{28}$ molecules.s-1.

Table 3 summarises these results.

Table 3: Parent and daughter scalelengths calculated for CN and N${}_{2}^{+}$ (at 1 au, scaled with $R_{h}^{2}$) and the resulting production rates for their parent species. |  $l_{p}(\times 10^{3}$ km) |  $l_{d}(\times 10^{3}$ km) | Q (molecules.s-1)
---|---|---|---
CN | 13 | 280 | $5\times 10^{24}$
N${}_{2}^{+}$ | 2800 | 3800 | $1\times 10^{28}$

Large uncertainties still remain, and as it is so far the only comet of this
kind to ever be observed with modern equipment, we are lacking a point of
comparison. A detailed observation of C/2016 R2 at high heliocentric distances
could still be possible, as CO would continue to sublimate under 40 au, and
would provide further detail of the nucleus while inactive, in order to create
a detailed thermal and structural model of this peculiar comet.

## Data availability

No new data were generated or analysed in support of this research.

## References

  * A’Hearn et al. (1995) A’Hearn M. F., Millis R. C., Schleicher D. O., Osip D. J., Birch P. V., 1995, Icarus, 118, 223
  * Biver et al. (2018) Biver N., et al., 2018, Astronomy & Astrophysics, 619, A127
  * Cochran & Cochran (2002) Cochran A. L., Cochran W. D., 2002, Icarus, 157, 297
  * Cochran & McKay (2018) Cochran A. L., McKay A. J., 2018, ApJ, 854, L10
  * Cochran & Schleicher (1993) Cochran A. L., Schleicher D. G., 1993, Icarus, 105, 235
  * Cochran et al. (1992) Cochran A. L., Barker E. S., Ramseyer T. F., Storrs A. D., 1992, Icarus, 98, 151
  * Cordiner et al. (2020) Cordiner M. A., et al., 2020, Nature Astronomy, 4, 861
  * Cruikshank et al. (1993) Cruikshank D. P., Roush T. L., Owen T. C., Geballe T. R., de Bergh C., Schmitt B., Brown R. H., Bartholomew M. J., 1993, Science, 261, 742
  * Dekker et al. (2000) Dekker H., D’Odorico S., Kaufer A., Delabre B., Kotzlowski H., 2000, in Optical and IR Telescope Instrumentation and Detectors. pp 534–545
  * Eistrup et al. (2019) Eistrup C., Walsh C., van Dishoeck E. F., 2019, A&A, 629, A84
  * Feldman (2015) Feldman P. D., 2015, ApJ, 812, 115
  * Fink et al. (1991) Fink U., Combi M. R., Disanti M. A., 1991, ApJ, 383, 356
  * Greenstein (1962) Greenstein J. L., 1962, ApJ, 136, 688
  * Haser (1957) Haser L., 1957, Bulletin de la Societe Royale des Sciences de Liege, 43, 740
  * Huebner et al. (1992) Huebner W. F., Keady J. J., Lyon S. P., 1992, Ap&SS, 195, 1
  * Ivanova et al. (2016) Ivanova O. V., Luk‘yanyk I. V., Kiselev N. N., Afanasiev V. L., Picazzio E., Cavichia O., de Almeida A. A., Andrievsky S. M., 2016, Planet. Space Sci., 121, 10
  * Ivanova et al. (2018) Ivanova O. V., Picazzio E., Luk’yanyk I. V., Cavichia O., Andrievsky S. M., 2018, PASS, 157, 34
  * Korsun et al. (2008) Korsun P. P., Ivanova O. V., Afanasiev V. L., 2008, Icarus, 198, 465
  * Korsun et al. (2014) Korsun P. P., Rousselot P., Kulyk I. V., Afanasiev V. L., Ivanova O. V., 2014, Icarus, 232, 88
  * Kurucz et al. (1984) Kurucz R. L., Furenlid I., Brault J., Testerman L., 1984, Solar flux atlas from 296 to 1300 nm. National Solar Observatory Atlas, Sunspot, New Mexico: National Solar Observatory
  * Lutz et al. (1993) Lutz B., Womack M., Wagner R., 1993, The Astrophysical Journal, 407, 402
  * McKay et al. (2019) McKay A. J., et al., 2019, AJ, 158, 128
  * Merlin et al. (2018) Merlin F., Lellouch E., Quirico E., Schmitt B., 2018, Icarus, 314, 274
  * Mousis et al. (2021) Mousis O., Aguichine A., Bouquet A., Lunine J. I., Danger G., Mandt K. E., Luspay-Kuti A., 2021, arXiv e-prints, p. arXiv:2103.01793
  * Newburn & Spinrad (1989) Newburn R. L., Spinrad H., 1989, AJ, 97, 552
  * Opitom et al. (2019) Opitom C., et al., 2019, Astronomy & Astrophysics, 624, A64
  * Owen et al. (1993) Owen T. C., et al., 1993, Science, 261, 745
  * Quirico et al. (1999) Quirico E., Douté S., Schmitt B., de Bergh C., Cruikshank D. P., Owen T. C., Geballe T. R., Roush T. L., 1999, Icarus, 139, 159
  * Raghuram et al. (2020) Raghuram S., Bhardwaj A., Hutsemékers D., Opitom C., Manfroid J., Jehin E., 2020, Monthly Notices of the Royal Astronomical Society, 501, 4035
  * Raghuram et al. (2021) Raghuram S., Bhardwaj A., Hutsemékers D., Opitom C., Manfroid J., Jehin E., 2021, MNRAS, 501, 4035
  * Randall (1992) Randall C., 1992, BAAS, 24, 3
  * Rousselot et al. (2022) Rousselot P., Anderson S., Alijah A., Noyelles B., Opitom C., Jehin E., Hutsemékers D., J. M., 2022, A&A, in press
  * Rubin et al. (2015) Rubin M., et al., 2015, Science
  * Rubin et al. (2020) Rubin M., Engrand C., Snodgrass C., Weissman P., Altwegg K., Busemann H., Morbidelli A., Mumma M., 2020, Space Sci.Rev., 216, 102
  * Schleicher (2010) Schleicher D. G., 2010, AJ, 140, 973
  * Wierzchos & Womack (2018) Wierzchos K., Womack M., 2018, The Astronomical Journal, 156, 34
  * Womack et al. (1994) Womack M., Lutz B. L., Wagner R. M., 1994, ApJ, 433, 886
  * Wyckoff & Theobald (1989) Wyckoff S., Theobald J., 1989, Advances in Space Research, 9, 157
  * Wyckoff & Wehinger (1976) Wyckoff S., Wehinger P. A., 1976, ApJ, 204, 604
  * Wyckoff et al. (1998) Wyckoff S., Heyd R. S., Fox R., 1998, in AAS/Division for Planetary Sciences Meeting Abstracts #30. p. 29.03
  * Zucconi & Festou (1985) Zucconi J. M., Festou M. C., 1985, A&A, 150, 180
  * de La Baume Pluvinel & Baldet (1911) de La Baume Pluvinel A., Baldet F., 1911, ApJ, 34, 89