Cold Traps of Hypervolatiles in the Protosolar Nebula at the Origin of the Peculiar Composition of Comet C/2016 R2 (PanSTARRS)

Recent observations of the long-period comet C/2016 R2 (PanSTARRS; hereafter R2) indicate an unusually high N2/CO abundance ratio, typically larger than ∼0.05, and at least 2–3 times higher than the one measured in 67P/Churyumov–Gerasimenko. Another striking compositional feature of this comet is its heavy depletion in H2O (H2O/CO ∼ 0.32%), compared to other comets. Here we investigate the formation circumstances of a generic comet whose composition reproduces these two key features. We first envisage the possibility that this comet agglomerated from clathrates, but we find that such a scenario does not explain the observed low water abundance. We then alternatively investigate the possibility that the building blocks of R2 agglomerated from grains and pebbles made of pure condensates via the use of a disk model describing the radial transport of volatiles. We show that N2/CO ratios reproducing the value estimated in this comet can be found in grains condensed in the vicinity of the CO and N2 ice lines. Moreover, high CO/H2O ratios (>100 times the initial gas-phase value) can be found in grains condensed in the vicinity of the CO ice line. If the building blocks of a comet assembled from such grains, they should present N2/CO and CO/H2O ratios consistent with the measurements made in R2’s coma. Our scenario indicates that R2 formed in a colder environment than the other comets that share more usual compositions. Our model also explains the unusual composition of the interstellar comet 2l/Borisov.


Introduction
Comets are supposed to be water-rich bodies that are relics of the formation of the solar system. Because they have undergone little alteration during the 4.6 billion years of the solar system evolution, investigating their composition provides indications of the physicochemical conditions that were at play during their formation in the protosolar nebula (PSN). An outstanding question remaining to be solved is the origin of the apparent N 2 deficiency observed in comets, while both Pluto and Triton, also formed in the outer solar system, harbor N 2 -rich surfaces and atmospheres (Lellouch et al. 2011;Cruikshank et al. 2015;Mandt et al. 2017). Decades of remote-sensing observations of comets suggest that they are depleted in N 2 (Cochran et al. 2000). Only upper limits of ∼10 −4 for + N 2 /CO + were derived in the coma of comets 122P/ 1995 S1 (de Vico) and C/1995 O1 (Hale-Bopp) (Cochran et al. 2000), while very few detections of + N 2 emission lines have been reported in other comets from ground-based facilities. This apparent N 2 depletion was interpreted as the result of the selective trapping of CO at the expense of N 2 in the building blocks of comets presumably agglomerated from clathrates (Iro et al. 2003), or as the result of their partial devolatilization due to radiogenic heating (Mousis et al. 2012). To a lesser extent, these important depletions are in agreement with the measurements of the N 2 /CO ratio acquired in comet 67P/Churyumov-Gerasimenko (hereafter 67P/C-G) by the ROSINA mass spectrometer aboard the Rosetta spacecraft very early in the mission in 2014 October, which suggest a value of (5.70 ± 0.66) × 10 −3 that is depleted by a factor of ∼25 as compared to the one derived from protosolar N and C abundances (Rubin et al. 2015). However, this N 2 /CO ratio was measured at a heliocentric distance beyond 3 au, far beyond perihelion. Rubin et al. (2020) later derived an N 2 /CO ratio of ∼2.87 × 10 −2 in 67P/C-G, a value obtained in 2015 May a few months before perihelion passage, and roughly 5 times larger than the previous one.
Contrasting with most of these previous measurements, recent optical spectra of the long-period comet C/2016 R2 (PanSTARRS; hereafter R2) performed at a heliocentric distance of ∼3 au showed that its spectrum was largely dominated by the emission band of CO + , but also surprisingly by the presence of the emission band of + N 2 (Cochran & McKay 2018). Subsequent observations confirmed an unusually high N 2 /CO ratio in R2, a value estimated to range between 0.06 ± 0.01 (Opitom et al. 2019) and 0.08 (Biver et al. 2018), and at least 2-3 times higher than the one measured closer to perihelion in 67P/C-G. Another striking compositional feature of this comet is its heavy depletion in H 2 O (H 2 O/CO ∼ 0.32%; McKay et al. 2019). These trends are at odds with our understanding of the thermodynamic evolution of cometary nuclei, which should be, in principle, depleted in the most volatile species instead of being enhanced, as R2 shows, after multiple solar passages. Since sublimation of supervolatile ices can occur possibly beyond 40 au for CO (Womack et al. 2017), one should expect CO outgassing, and then cometary activity over a significant fraction of R2ʼs orbit. Figure 1 represents a comparison of the volatile composition of R2 to the average comet made by McKay et al. (2019). The figure shows that all oxygen appears to be locked in CO and CO 2 , and not H 2 O, as is the case for typical comets.
In this paper, we aim at investigating the formation circumstances of a generic comet whose composition reproduces the two key features of R2ʼs coma, namely, its high N 2 /CO ratio and the observed significant water depletion, assuming that they reflect its primordial composition (McKay et al. 2019). We first investigate the possibility that this comet could have agglomerated from clathrates as has been suggested in the case of 67P/C-G (Luspay- Mousis et al. 2016Mousis et al. , 2018, but we find that such a scenario cannot explain the observed low water abundance. We then alternatively investigate the possibility that the building blocks of this comet agglomerated from grains and pebbles made of pure condensates crystallized in the vicinity of the CO ice line in the PSN. This scenario matches both observed features and suggests that this comet formed in a colder environment than the other comets sharing usual compositions.

Agglomeration from Clathrates?
For decades, the agglomeration of comets from clathrates has been regularly invoked to explain their compositional and thermodynamic properties (Lunine & Stevenson 1985;Klinger et al. 1986;Iro et al. 2003;Marboeuf et al. 2010Marboeuf et al. , 2011Marboeuf et al. , 2012. Recent measurements of 67P/C-Gʼs composition by the Rosina mass spectrometer aboard the Rosetta mission have also been interpreted in favor of the presence of clathrates in its interior Mousis et al. 2016Mousis et al. , 2018. This possibility is now investigated in the case of R2, which presents unusually high CO and N 2 abundances in its coma, compared to typical comets. To do so, we aim at comparing the plausible composition of multiple guest (MG) clathrates formed in the PSN to that of R2, assuming that the latter would have agglomerated from these icy structures. These MG clathrates are assumed to be formed from a gaseous mixture composed of N 2 and CO, assuming that these species are the two main Cand N-bearing volatiles in the gas phase of the disk, a hypothesis consistent with the thermochemical models (Lewis & Prinn 1980;Prinn & Fegley 1989;Mousis et al. 2002). Because CO has a much higher propensity for clathration than N 2 , the calculations utilized in this work are those performed in the framework of structure I clathrates, namely, the structure predicted for a CO-dominated clathrate (Mohammadi et al. 2005). Given the fact that half of protosolar carbon is expected to be used in organic compounds in the PSN (Pollack et al. 1994), only the remaining half is assumed to be in the form of CO. Assuming that 90% of protosolar nitrogen is in N 2 form (Lewis & Prinn 1980), we derive a N 2 /CO ratio of 2.66 × 10 −1 in the initial gas phase of the PSN, based on the protosolar abundances given by Lodders et al. (2009). In our model, clathrates form in the PSN as long as crystalline water is available. Here we assume that once all the water budget has been used for clathration, the remaining volatiles are not incorporated into solids.
Our model assumes the formation of an MG clathrate with an equilibrium pressure expressed as (Lunine & Stevenson 1985) where y i is the mole fraction of the component i in the fluid phase. The equilibrium pressure curves of each species are determined by fitting the available theoretical and laboratory data (Lunine & Stevenson 1985) with equations of the form , where P eq,i and T are the partial equilibrium pressure and temperature of the considered species i, respectively. The relative abundances of guest species incorporated in an MG clathrate formed at a given temperature and pressure from the PSN gas phase are calculated following the method described in Lunine & Stevenson (1985) and Mousis et al. (2010). This formalism has been used to interpret 67P/C-Gʼs ice structure and composition from Rosetta/ ROSINA observations of its coma (Mousis et al. , 2018. Figure 2 shows the influence of the abundance of crystalline water ice in the outer PSN on the N 2 /CO ratio in cometary grains at clathration temperatures of 25, 30, 40, and 50 K.  These temperatures are all above the formation temperature of N 2 and CO pure condensates. For water abundances below 1.3 × (O/H 2 ) e , the N 2 /CO ratios in icy grains correspond to those calculated in clathrates at specified formation temperatures. For example, the N 2 /CO ratios measured in Hale-Bopp, de Vico, and 67P/C-G are closely matched at 25, 30, and 50 K, respectively. At higher water abundances, CO is fully entrapped and significant amounts of N 2 become progressively trapped in clathrates, implying an increase of the global N 2 /CO in icy grains. A water abundance of ∼1.4 × (O/H 2 ) e in the PSN allows one to retrieve an N 2 /CO ratio in clathrates corresponding to the value measured in R2. Water abundances greater than ∼1.7 × (O/H 2 ) e allow the full trapping of the available N 2 in clathrates crystallized in the PSN, leading to an N 2 /CO ratio of ∼2.66 × 10 −1 in cometary grains, corresponding to the value calculated from protosolar N and C (with half C in CO and 90% N in N 2 ). These calculations illustrate the fact that, in principle, a slight increase of the water abundance in the PSN can lead to the full trapping of CO and enough N 2 to reproduce the N 2 /CO ratio in R2.
However, the amount of water required by clathrates formed in the PSN to match the N 2 /CO ratio measured in R2ʼs coma is inconsistent with its estimated low abundance. To match the N 2 /CO measured in R2, our clathrate model requires an H 2 O/CO ratio of 6.1. In contrast, the H 2 O/CO mixing ratio measured at 2.8 au in R2 is ∼0.32%, a value several orders of magnitude lower than those measured in 67P/C-G (∼200%-1900%) at a similar distance range (Gasc et al. 2017) or in many other comets (McKay et al. 2019). As a result, this mechanism fails at challenging to explain the H 2 O depletion in R2-like comets via this mechanism.

Agglomeration from Particles Formed in the Vicinity of CO and N 2 Ice Lines
Disk models including radial transport of volatiles can display significant enrichments of both volatile and refractory matter at its condensation location, implying potentially drastic changes of the local metallicity in the disk (Stevenson & Lunine 1988;Cyr et al. 1999;Ali-Dib et al. 2014;Mousis et al. 2019Mousis et al. , 2020Aguichine et al. 2020). Here we investigate this effect by considering the possibility that R2 formed in the vicinity of the condensation locations of CO and N 2 crystalline ices in the PSN.

Model
The volatile transport and distribution model used in our work is the one described in Aguichine et al. (2020) and Mousis et al. (2020), to which the reader is referred for details. Three energy sources are considered in our model, namely, viscous heating, irradiation from the current Sun, and ambient constant irradiation giving a background temperature of 10 K (see Equation (6) of Aguichine et al. 2020). Here the accretion rate onto the star declines over time, and the disk develops a transition radius between inward and outward gas flows. The location of this transition radius is moving outward with time. In a few words, our time-dependent PSN model is governed by the following differential equation (Lynden-Bell & Pringle 1974): This equation describes the time evolution of a viscous accretion disk of surface density Σ g of viscosity ν, assuming hydrostatic equilibrium in the z-direction. The viscosity is calculated in the framework of the α-formalism (Shakura & Sunyaev 1973) using the following method. For each distance r to the Sun, the diskʼs properties are calculated by solving the equation of energy balance between viscous heating and radiative transfer at the midplane level. This gives us ν, as well as the pressure and temperature profiles of the disk as a function of r. The evolution of the disk starts with an initial profile given by n S µ -r exp p g 2 ( ), with = p 3 2 for an early disk (Lynden-Bell & Pringle 1974). In our computations, the initial diskʼs mass is fixed to 0.1 M e . The computational box is set equal to 500 au, allowing 99% of the disk mass to be encapsulated within ∼100 au. Mass can only be lost by accretion onto the Sun or via outward diffusion, the boundary between these two regimes being defined by the location of R c , the centrifugal radius. R c is initially located at ∼6 au and expands up to ∼20 and 45 au after 0.1 and 1 Myr of evolution, respectively. The diskʼs initial mass accretion rate onto the Sun is set to 10 −7.6 M e yr −1 (Hartmann et al. 1998).
The size of dust particles used in our model is determined by a two-population algorithm derived from Birnstiel et al. (2012). This algorithm computes the representative size of particles through the estimate of the limiting Stokes number in various dynamical regimes. In our model, dust is initially present in the form of particles of sizes a 0 = 10 −7 m and grows through mutual collisions. This growth is limited by the maximum sizes imposed by fragmentation or by the drift velocity of the grains (see Aguichine et al. 2020 for details). The dust surface density is the sum over all surface densities of available solids at a given time and location, assuming a protosolar ice-to-rock ratio of ∼1.58 (Lodders et al. 2009) and a bulk density of 1000 kg m −3 .
We follow the approaches of Desch et al. (2017) and Drążkowska & Alibert (2017) for the dynamics of trace species in terms of motion and thermodynamics, respectively. We assume that the disk is uniformly filled with H 2 O, CO, and N 2 . The abundances of CO and N 2 are those derived in Section 2, and the abundance of H 2 O is that of the leftover oxygen. No chemistry is assumed to happen between the trace species. In our simulations, grains are considered as homogeneous mixtures of solid species in proportions determined from the aforementioned assumptions (relative abundance ratios, condensation, sublimation). Sublimation of grains occurs during their inward drift when partial pressures of trace species become lower than the corresponding vapor pressures. Once released, vapors diffuse both inward and outward. Because of the outward diffusion, vapors can recondense back in solid form following the rates defined by Drążkowska & Alibert (2017), and condensation occurs either until thermodynamic equilibrium is reached or until no more gas is available to condense. The position of the ice line of a given species is defined as the location where its gaseous abundance equals that of its coexisting solid (Lodders 2003;Öberg & Wordsworth 2019;Aguichine et al. 2020). This approach results in ice lines closer to the Sun than those computed via a simple comparison between the partial and saturation pressures.
The motion of dust and vapor is computed by integrating the 1D radial advection−diffusion equation derived from Birnstiel et al. (2012) and Desch et al. (2017) and detailed in Aguichine et al. (2020). The vapor pressures of trace species are taken from Fray & Schmitt (2009). Figure 3 represents the time evolution of the radial profiles of the N 2 /CO ratio relative to its assumed initial ratio (∼2.66 × 10 −1 ), defined by the enrichment factor f, in both solid and gas phases of the PSN, and for three values of the viscosity parameter α, namely, 5 × 10 −3 , 10 −3 , and 10 −4 . The adopted α values are well within the range of those typically used in models of protoplanetary disks (Hersant et al. 2001;Nelson et al. 2013;Simon et al. 2015). The figure shows that f is flat almost everywhere in the PSN, except in the vicinity of the CO and N 2 ice lines, which are located at very close distances from each other, i.e., less than a few tenths of an astronomical unit, in the 10-15 au region of the PSN, with the CO ice line situated a bit closer in in the disk.

Results
Depending on the choice of the α parameter, the N 2 /CO ratio increases to more than ∼10 in the gas phase in the area of the two ice lines. This peak corresponds to the supply of N 2 vapor when N 2 -rich dust drifts inward of the N 2 ice line, which is in excess compared to the CO vapor supplied via backward diffusion beyond the CO ice line. When significant, this peak is preceded by a decrease of the N 2 /CO, which corresponds to the supply of CO vapor when CO-rich dust drifts inward of the CO ice line. The N 2 /CO ratio in solid phase also experiences important depletions, which correspond to the location where N 2 essentially forms vapor while CO remains in solid phase. The N 2 /CO ratio in solid phase also reaches peaks up to f ∼ 10, depending on the α value, at the location of the N 2 ice line, corresponding to the formation of solid N 2 from N 2 vapor diffusing backward its corresponding ice line. The figure also shows that it is possible to form dust in a narrow region of the PSN, i.e., within 10-15 au, with N 2 /CO ratios matching the value estimated in R2. A comet agglomerated from these grains would present an N 2 /CO ratio consistent with the R2 value, in both its interior and its coma, because the two molecules present similar volatilities, as testified by the proximity of their ice lines in the PSN. Figure 4 represents the time evolution of the radial profiles of the CO/H 2 O ratio relative to its assumed initial ratio (∼2.96 × 10 −1 ), defined by the enrichment factor f, in both solid and gas phases of the PSN, and for three values of the viscosity parameter α, namely, 5 × 10 −3 , 10 −3 , and 10 −4 . The CO/H 2 O ratio in the gas phase increases steeply in the vicinity of the H 2 O ice line, which is located at much closer distances from the Sun than the CO and N 2 ice lines. This steep increase is due to the freezing of the H 2 O vapor when it diffuses outward through its ice line. The CO/H 2 O ratio is deeply depleted in the solid phase at distances inward of the CO ice line because this species remains in gaseous form. On the other hand, when approaching the CO ice line, the CO/H 2 O ratio presents values exceeding its initial value. A peak forms at this location, and its magnitude depends on the adopted viscosity parameter α of the disk. The CO/H 2 O enrichment factor f can typically reach a value of ∼10 with α = 10 −3 , while it peaks at ∼200 with α = 10 −4 . If the building blocks of a comet were assembled from grains formed at this location of the PSN, they should present a CO-rich and H 2 O-poor composition. Assuming f = 100, the CO/H 2 O ratio would be ∼30, a value exceeding by far those known in typical comets (Bockelée-Morvan & Biver 2017). We note that this value is still about 10 times lower than the one inferred in R2ʼs coma, but the latter can be seen as an overestimate of the bulk CO abundance in the comet. Indeed, CO is much more volatile than H 2 O, and its vapor comes from deeper layers than H 2 O vapor in the nucleus, potentially inducing an enriched CO/H 2 O ratio in the coma. This effect has been investigated in depth by Marboeuf & Schmitt (2014), who found that the abundance of volatile molecules (relative to H 2 O) released from the interiors of nuclei can vary by several orders of magnitude along the cometʼs orbital evolution.
Interestingly, Figures 3 and 4 show that the locations of the CO and N 2 ice lines do not significantly evolve with time and are only weakly affected by the value adopted for the viscosity parameter. The main reason for this is that viscous heating occurs mainly in the dense regions of the PSN. With our parameterization, the gas density quickly decreases after 10 au, implying that the temperature profile beyond this distance is mostly ruled by solar irradiation, here assumed to be constant.

Conclusions
In this work, we first show that a comet agglomerated from clathrates crystallized in the PSN could, in principle, display N 2 /CO ratios corresponding to the value measured in R2ʼs coma. This finding contrasts with the suggestion of Wierzchos & Womack (2018) to rule out the presence of clathrates in R2 Figure 3. From left to right: radial profiles of the N 2 /CO ratio relative to its initial abundance (defined by the enrichment factor f ) calculated as a function of time in the PSN for viscosity parameters α = 5 × 10 −3 , 10 −3 , and 10 −4 . Dashed and solid lines correspond to vapor and solid phases, respectively. The blue bar corresponds to the N 2 /CO ratio measured in R2.
only because its N 2 /CO ratio is higher than those measured in other comets. This statement is correct in the case of a limited budget of available crystalline water in the PSN. However, if water is abundant enough, as shown in this paper, both molecules are efficient clathrate formers. On the other hand, the amount of water required by clathrates formed in the PSN to match the N 2 /CO ratio measured in R2ʼs coma, namely, H 2 O/CO ∼ 6, is inconsistent with its estimated extremely low value (∼0.32%).
We have then used a disk model describing the radial transport of volatiles in the PSN, evolving in both gaseous and solid (pure condensates) phases. We show that N 2 /CO ratios reproducing the value estimated in R2 can be found in dust formed in the vicinity of the CO and N 2 ice lines, i.e., within the 10-15 au region of the PSN, depending on the adopted viscosity parameter α of the disk. Meanwhile, very high CO/H 2 O ratios, i.e., up to more than 200 times the value derived from a protosolar gaseous mixture, can be found in dust formed in the vicinity of the CO ice line (∼10 au), the extent of which also depends on the adopted viscosity parameter. If the building blocks of a comet assembled from grains formed close to the CO ice line in the PSN, they should present N 2 /CO and CO/H 2 O ratios consistent with the measurements made in R2ʼs coma.
The likelihood of finding R2-like comets versus H 2 O-rich comets in the PSN can be, to first order, assessed via comparing the masses of their respective reservoirs. To do so, we calculated the mass of H 2 O and CO ices contained in a PSN annulus in which CO was found to be more abundant than H 2 O (R2-like comets-reservoir 1), and in another annulus where the CO/H 2 O ratio was set to 2.96 × 10 −1 , i.e., the value derived from protosolar O and C abundances (average comets -reservoir 2). The time evolution of the reservoir 1/reservoir 2 mass ratio is represented in Figure 5 for α = 10 −3 and 10 −4 . Higher α values do not allow for the formation of reservoir 1 in the PSN. This is also shown by Figure 4 in the case α = 5 × 10 −3 , where the CO/H 2 O ratio in the solid phase never exceeds twice the initial ratio. Figure 5 illustrates the fact that the mass ratio between the two reservoirs strongly depends on the adopted α parameter and the epoch of the PSN evolution. In the case α = 10 −3 , the mass of R2-like comets is about 1% that of the average comets, if the grains from which they assembled condensed between 0.1 and 1 Myr in the PSN. In the case α = 10 −4 , the mass of R2-like comets varies between less than ∼1% and 100% that of the average comets, depending on the particular grain formation epoch in the 0.1-1 Myr time frame. We conclude that the likelihood of finding R2-like comets versus H 2 O-rich comets in the PSN strongly depends on the adopted disk parameters and thus remains difficult to predict a priori.
The heliocentric distances indicated for the different ice lines must be taken with caution and are model dependent. However, the agglomeration of R2 from dust condensed in the region of the CO ice line indicates that this comet formed at a greater heliocentric distance than an H 2 O-rich comet formed from clathrates because the formation temperatures of N 2 -and COrich clathrates are always higher than those of the pure condensates at PSN conditions. Our model also applies to other protoplanetary disks and can explain the unusual composition of the interstellar comet 2l/Borisov, which seems to contain substantially more CO than H 2 O gas in its coma (Bodewits et al. 2020).
Interestingly, physicochemical models of nitrogen chemistry in protostellar disks show that photodissociation of N 2 leads to production of HCN