Adrienn Luspay-Kuti

Olivier Mousis

Francoise Pauzat

Ozge Ozgurel

Yves Ellinger

Jonathan I Lunine

Stephen A Fuselier

Kathleen E Mandt

Karlheinz J Trattner

Steven M Petrinec

Gerasimenko Adrienn Luspay-Kuti

Françoise Pauzat

Dual storage and release of molecular oxygen in comet 67P/Churyumov-Gerasimenko

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

One of the biggest surprises of the Rosetta mission was the detection of O2 in the coma of 67P/Churyumov-Gerasimenko in remarkably high abundances. The measured levels of O2 in the coma are generally assumed to reflect the overall abundance and chemical origin of cometary O2 in the nucleus. Along with its strong association with H2O and weak association with CO and CO2, these measurements led to the consensus that the source and release of cometary O2 are linked to H2O.

We analyzed ROSINA observations and found a previously unrecognized change in the correlations of O2 with H2O, CO2 and CO that contradict the prevailing notion that the release of O2 is linked to H2O at all times. These findings can be explained by the presence of two distinct reservoirs of O2: a pristine source in the deeper nucleus layers dating back prior to nucleus formation, and a water-trapped secondary reservoir formed during the thermal evolution of the nucleus. These results imply that O2 must have been incorporated into the nucleus in a solid and distinct phase during accretion in significantly lower abundances than previously assumed.

Main

The Rosetta Orbiter Spectrometer for Ion and Neutral Analysis/Double Focusing Mass Spectrometer (ROSINA/DFMS 1 ) detected O2 for the first time in a cometary environment [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF] . The detected levels of O2 in the coma of 67P/Churyumov-Gerasimenko (67P) were about an order of magnitude higher than predicted by models describing the chemistry in molecular clouds [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF][START_REF] Taquet | Multilayer modeling of porous grain surface chemistry-I. The GRAINOBLE model[END_REF] , and well over 2 orders of magnitude higher than the O2 observed in the gas phase in ρ Oph [START_REF] Larsson | Molecular oxygen in the ρ Ophiuchi cloud[END_REF][START_REF] Liseau | O 18 O and C 18 O observations of ρ Ophiuchi A[END_REF] .

The first measurements of O2 in the coma of 67P were taken over a six-month period (September 2014-March 2015) early in the mission [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF] . Over this time, the in situ O2 abundance relative to H2O was between 1 and 10%, which made it the fourth most abundant coma species following H2O, CO2 and CO. O2 had a persistently high abundance relative to H2O, the O2/H2O ratio remained fairly constant, and O2 showed a strong correlation with H2O, but it weakly correlated with CO and N2 -two molecules with similar volatility [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF] . These observations led to the initial conclusion that the origin and/or the release of O2 and H2O were linked, and that cometary O2 is likely primordial [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF] .

The O2-H2O correlation was found to extend through later parts of the comet's orbit [6][7][8] , further suggesting that O2 is controlled by H2O.

The high relative O2 abundance, the strong O2 -H2O correlation, and the weak O2 -CO correlation were key requirements that all proposed formation and production hypotheses needed to satisfy 9 . While the strong association of O2 with H2O provides a very important clue to the mystery of cometary O2, the details of when, where, and for how long this association is truly present make the interpretation less obvious. 67P experiences strong hemispheric differences in outgassing, largely due to the seasons, with the majority of CO2 coming from the southern hemisphere 6,8 . The minor species also showed strong latitudinal heterogeneity, with the less volatile species typically following H2O and the more volatile species following CO2 (refs. 10-14). However, the correlations between the major and minor volatiles appear to be very different when looking at shorter (several days) versus longer (months to entire orbital segments) time intervals 8,13 . This makes the interpretation of the observations potentially dependent on the location and the time interval considered.

This study presents abundances (count rates) from ROSINA/DFMS 1 of the four most abundant species in the innermost coma (< 35 km) over significant time intervals at ~2.5 -3.7 au. We looked at short-and long-term correlations over the southern hemisphere only (sub-spacecraft latitudes between 25 and 90°S) 3-7 months before the pre-perihelion equinox (May 2015) and after the post-perihelion equinox (March 2016). The southern hemisphere was poorly illuminated at these times, which resulted in lower temperatures, hence lower H2O production rates over the southern hemisphere compared to the northern hemisphere 9-13, [START_REF] Biver | Long-term monitoring of the outgassing and composition of comet 67P/Churyumov-Gerasimenko with the Rosetta/MIRO instrument[END_REF][START_REF] Bockelée-Morvan | First observations of H2O and CO2 vapor in comet 67P/Churyumov-Gerasimenko made by VIRTIS onboard Rosetta[END_REF][START_REF] Migliorini | Water and carbon dioxide distribution in the 67P/Churyumov-Gerasimenko coma from VIRTIS-M infrared observations[END_REF] . That, along with the stronger erosion experienced by the southern hemisphere [START_REF] Keller | Seasonal mass transfer on the nucleus of comet 67P/Chuyumov-Gerasimenko[END_REF] provided a unique opportunity to see relationships otherwise masked by the strong H2O outgassing [11][12][13] .

Correlations with the major cometary species

Figs. 1 and2 show the correlations between the major cometary species before the preperihelion and after the post-perihelion equinox over the southern hemisphere. The strong O2-H2O correlation before the pre-perihelion equinox (Fig. 1) is consistent with previous results [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF]7 . At the same time, the Pearson correlation coefficients suggest a moderate correlation between O2-CO2 and O2-CO. These data are strongly heteroscedastic (nonequal variance), so the Pearson coefficient of correlation is meaningless and misleading as an indicator of a relationship between two species. Looking at shorter time intervals, there is a moderate correlation early on in the mission at ~3.3-3.2 au (September 2014) between O2 and CO2 as well as O2 and CO. That correlation weakens and disappears by ~3.17 au (mid-October of 2014), and the species remain uncorrelated/weakly correlated all the way to the end of our pre-perihelion equinox interval at ~2.4 au in January 2015 (Table 1).

O2 and H2O seem to be well-correlated again at first glance after the post-perihelion equinox (Fig. 2). However, when divided into shorter time intervals O2 and H2O only correlate for the earlier three periods, but show no correlation over the short intervals past ~3.3 au post-perihelion equinox (Table 1). This is key because before the pre-perihelion equinox, when the O2-H2O correlation was this strong over the whole interval (Fig. 1), the correlations over the shorter time periods were also very strong. Moreover, when O2 does not correlate with H2O, it strongly correlates with CO and CO2 (Table 1). While the count rates of all species have decreased by this time, H2O production decreases even more rapidly and drops below those of CO2 and CO in the southern hemisphere 8,13 . This also happens to be around the same time that the correlation between O2 and H2O begins to weaken. For O2-CO2 and O2-CO it is the exact opposite: they do not correlate up to ~3.3 au following the post-perihelion equinox, but they correlate very strongly beyond (Table 1). The correlations between O2, CO2 and CO and the lack of correlation between O2 and H2O are key because they show a previously unrecognized change in outgassing that comes with important consequences. This shift in correlations may happen because the measurements over the time interval after the post-perihelion equinox span different sublimation regimes. Over the time intervals, or distances covered by the yellow and dark blue dots (Fig. 2), the comet still experienced an appreciable H2O sublimation, albeit decreasing over time. After that point, the H2O sublimation began to turn off due to the low surface temperatures over the poorly illuminated southern hemisphere. That is precisely where the transition in correlations after the post-perihelion equinox occurs, and where the shorter time intervals indicated by the purple, green and light-blue markers in Fig. 2 begin to show no correlation between O2 and H2O. Thus, in Figure 2 we are mixing two different regimes of sublimation and outgassing. When mixing data for different time intervals driven by different physical processes, as is the case in Fig. 2, relationships revealing key information may be masked, and the linear regression may lead to inadvertently misconstrued interpretations. In order to get a true sense of the changes occurring between the two sublimation regimes, we honed in on this later time interval between ~3.3 -3.7 au, starting where the correlation between O2 and H2O declines and where H2O production is lower than it was at any time before the pre-perihelion equinox (Fig. 3). The interesting thing to notice here is that despite the lack of correlation between O2 and H2O for the shorter time periods (Table 1), there is still a moderate correlation for the whole interval (Fig. 3).

Diurnal re-freezing of sublimating volatiles from the subsurface onto the nucleus surface can be considered as a possible mechanism at play to explain the O2-CO-CO2 correlations. Diurnal refreezing and subsequent sublimation of H2O was invoked as an explanation of outbursts observed on comet 9P/Tempel 1 (ref. 19). Surface ice patches of H2O and CO2 detected on 67P were attributed to the refreezing of volatiles coming from deeper layers, providing a source of localized activity [START_REF] Filacchione | Exposed water ice on the nucleus of comet 67P/Churyumov-Gerasimenko[END_REF][START_REF] Sanctis | The diurnal cycle of water ice on comet 67P/Churyumov-Gerasimenko[END_REF] . However, the measured amounts of CO2 and CO at >3.3 au post-perihelion (Fig. 3) are simply too high to explain the observed strong correlations between O2, CO and CO2. The strong correlation with O2 extends all the way to the highest measured CO and CO2 counts and the observations suggest a continuous background activity instead of short-lived surges or minioutbursts [START_REF]Exposed water ice deposits on the surface of comet 9P/Tempel 1[END_REF] . The high observed fluxes of CO and CO2 indicate that these species come out together from the interior of the nucleus, suggesting that significant surface contribution to the observed correlations is unlikely.

Two distinct O2 nucleus reservoirs

Instead, our observations can be explained by the following. As the nucleus cools near the post-perihelion equinox following an intense southern hemisphere summer, the surface and near-surface temperatures drop until they become low enough that H2O sublimation declines precipitously. After the post-perihelion equinox when the H2O counts were well below the pre-perihelion equinox levels (Fig. 3), the hypervolatile species CO, CO2 and O2 still leave the porous nucleus from the deeper, colder, and more volatile-rich layers. The subsurface temperature is too low for H2O to sublimate at these depths. During this release, some of the O2 leaves the nucleus and enters the coma, but some of it gets trapped in the porous near-surface H2O ice depleted in the more volatile species (Fig. 4a). O2 can be (re-)trapped in porous H2O very easily because O2 attaches to the surrounding solid water with significant stabilization energy but without actually modifying the ice environment. A related scenario was considered for the radiolytic production of O2 from and incorporation into H2O ice grains in the ISM 23 . In that case, H2O molecules are removed by radiolysis of icy grains, creating voids for the trapping of O2. Here, we propose that macroscopic cavities pocket the H2O ice consistent with the high porosity of the nucleus 24 and as a result of the H2O ice in the upper layers being depleted in the more volatile species due to sublimation; the latter process leaving voids behind. Hence, O2 trapping in porous H2O ice as proposed here is a process that could be happening in the current-day nucleus of 67P (Methods).

Since H2O sublimation beyond ~3.3 au is low (Fig. 3), the trapped O2 may be bound to the H2O structure and unable to leave until the next onset of H2O sublimation. In that case, the O2 measured in the coma by ROSINA/DFMS at distances >3.3 au postperihelion is O2 released from depth together with CO and CO2. This O2 coming from the deeper nucleus layers may either leave the nucleus or keep accumulating in the nearsurface H2O ice while away from the Sun. This accumulation can happen as long as there is O2 release happening at depth and there are enough cavities in the near-surface H2O ice to trap it. This mechanism leads to two O2 reservoirs with distinct ice properties: a polar-ice dominated source with O2 trapped in H2O closer to the surface, and a non-polar (or weakly polar) source, where O2 is released from together with CO and CO2.

As 67P makes its way back toward the next pre-perihelion equinox, H2O sublimation in the southern hemisphere gradually turns on and the trapped O2 is released together with the sublimating H2O (Fig. 4b). The release of accumulated O2 from the secondary, watertrapped reservoir closer to the surface by the sublimation of H2O ice explains the surprisingly high O2 relative abundances measured by ROSINA/DFMS early in the mission [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF] . While O2 at depth continues to be released along with CO2 and CO, the accumulated O2 source in H2O ice is significantly stronger; hence the strong correlation between O2 and H2O up until the H2O sublimation begins to turn off.

Since the dominant source of O2 measured in the coma of 67P is the water-trapped O2 reservoir, the relatively small correlation between O2 and H2O in Fig. 3 can be explained by the general decrease of H2O as sublimation ceases. When H2O sublimation stops, O2 release from the trapped reservoir stops with it. With the dominant source of O2 cut off, the O2 density in the coma also dramatically decreases with the decreasing H2O.

However, looking at the short time intervals individually, especially where the spread of H2O counts is large enough (purple dots in Fig. 3), it becomes clear that O2 does not truly follow H2O. Combining several uncorrelated time intervals into a longer period makes it appear as if the O2-H2O data were moderately correlated due to the decline of the dominant O2 source, when in fact the correlation in this case does not imply a true association related to the formation of comet 67P. This dominant water-trapped reservoir explains the key pre-equinox observations 2 : the strong correlation of O2 with H2O pre-equinox, the small scatter in the O2/H2O ratio, and the persistently high O2 abundance relative to H2O pre-perihelion. Consequently, the surprisingly high O2 abundance is not indicative of the true abundance of O2 overall in the nucleus, but rather of the amount of O2 accumulated in and then released from the secondary, trapped H2O ice reservoir. That means that trying to reproduce these high relative abundances with chemical models of the origin of O2 is unnecessary, as these values overestimate the true importance of O2 relative to the major species. At the same time, the O2 coming from depth beyond 3.3 au. after the post-perihelion equinox does not represent the bulk O2 content of the nucleus either, as some of that O2 gets trapped in the porous, volatile-depleted, near-surface H2O ice. However, the strong erosion in the southern hemisphere around perihelion likely exposed fresh, more volatile-rich layers 12,[START_REF] Filacchione | Exposed water ice on the nucleus of comet 67P/Churyumov-Gerasimenko[END_REF]25,26 , so one might argue that the relative O2 abundances measured farther than ~3.3 au may be used to derive a lower limit for the bulk nucleus O2 content.

CO and CO2 can also incorporate into to the porous H2O ice via inclusion with a similar, but somewhat higher propensity (Methods, Table 3). Both of these molecules are able to fill the holes in the overlying H2O ice, but the large difference between their bulk abundances must be considered when evaluating their outgassing. CO is the second most abundant species in the ISM after H2, and it is one of the most abundant cometary species. While our interpretation suggests that the large count rates of O2 observed by ROSINA/DFMS closer than 3.3 au can be attributed to the accumulation and subsequent release of water-trapped O2, the majority of CO comes from the deeper and colder layers instead, as the bulk abundance of CO in the nucleus is much higher than that of O2. Thus, CO sublimating from depth at any time will overpower the CO that may be trapped and released from the H2O ice in the shallower nucleus layers upon H2O sublimation. The difference between the bulk abundances and the dominant source of O2 and CO leads to their fundamentally different correlations with H2O.

Implications for origin

Because of the association of O2 with CO and CO2, at least some and perhaps all of the comet's O2 is primordial, incorporated into the nucleus pre-agglomeration. Were the bulk of the O2 formed from H2O after nucleus agglomeration, i.e., via in situ reactions, we would see a much stronger correlated signal of co-sublimation of the two species. We cannot rule out, of course, that a minority of the O2 formed in this way. Recent experimental work indicates that high abundances of O2 can be formed from the photolysis of H2O:CO2 ices 27 . However, the fact that the isotopic ratios of O2 do not match the isotopic ratios of H2O in comet 67P (ref. 28) makes it less likely that O2 formed from H2O by radiolysis of ISM grains prior to their agglomeration into the nucleus. Instead, O2 in the comet nucleus largely or entirely originates from a solid phase independent of and distinct from H2O, and the O2 in comet 67P must have formed either during the PSN phase, or in the ISM via various chemical routes prior to nucleus agglomeration [29][30][31] . Once formed, O2 (and other volatiles) could have been adsorbed by amorphous grains in the ISM before being agglomerated by forming comets, if the PSN was cold enough to prevent the destabilization of amorphous ice. If the PSN was warm enough to destabilize the amorphous ice, then O2 was released in the gas phase of the disk. Subsequently, O2 either condensed or it was trapped in clathrates during the cooling of the disk. This implies that the primordial O2 may have been supplied either in the form of pure condensate, clathrates, or amorphous ice in the forming nucleus. Whatever the O2-bearing solid phase that has been incorporated into the cometary nucleus, its thermal evolution led to the progressive release of O2, and then to the formation of a secondary accumulated reservoir in the voids and pores of the overlying H2O matrix. This process could very well be happening on other comets too e.g., on 1P/Halley, where the presence of O2 was confirmed in comparable levels to 67P (ref. 32). A flyby, however, does not provide the same temporal and spatial resolution necessary for constraining the observations to specific times and locations with low H2O outgassing, as we did in this paper. Further analysis of O2 correlations with various minor volatile species (e.g., CH4, C2H6, N2) may narrow how O2 was supplied to the nucleus, and consequently, may unravel the origin of O2 in comet 67P.

Materials and Methods

ROSINA/DFMS Data Treatment

We used ROSINA/DFMS 1 mass spectra taken in high resolution mode for the analysis for masses 32 (O2), 18 (H2O), 28 (CO), and 44 (CO2). The L3 spectra have been corrected for the ageing of the Micro Channel Plate (MCP) by calibration measurements that were performed every few weeks throughout the mission.

For each individual mass spectrum, we performed a spectral fit to the species of interest by fitting a Gaussian to the pass peak. On some occasions, based on the shape of the peak two Gaussians were necessary to properly fit the shape (Supplementary Information). The area under the Gaussian(s) fitted to the mass peaks provides the abundance of a species at a given time. We then produced time series files containing the count rates (signal) of each species for each time a measurement was taken. We emphasize that we used count rates and not densities or concentrations.

The measurements have an associated ion statistical counting uncertainty given by √𝑛𝑛,

where n is the number of ions hitting the detector. Because the data were not corrected for instrument sensitivity, the counts in the figures and the slopes of the linear fits of the species do not represent absolute abundances 11 . Since we look at correlations between count rates and not absolute ratios or densities, any uncertainties in the absolute number stemming from using different gains for the different species is not important for the purpose of this study. Because the ion beam typically falls at the center of the MCP it causes a preferential degradation of these center pixels over time compared to the pixels along the edges of the MCP. The uncertainty in the individual pixel gains (i.e., mass peaks do not always fall on the same pixel) is important when determining concentrations, but it is less important when trying to get an idea about the correlation between two species. The individual pixel gain uncertainty was around 4% for the studied time intervals. We note that the major cometary species were well-measured throughout the mission; thus, the expected uncertainties are significantly lower than for other cometary species. Nevertheless, we applied a 10% uncertainty 2 in addition to the ion statistical counting uncertainty in order to quantify the possible variation of the Pearson correlation coefficients. The consistently small standard deviations around the mean correlation coefficients show that even with this more conservative uncertainty our determined correlations are robust.

The measurements used in this work are orders of magnitude higher than the background [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67P/Churyumov-Gerasimenko[END_REF] and well above the limiting sensitivity of ROSINA/DFMS. ROSINA/DFMS scanned through successive masses covering a mass range of 1-150 amu, with a maximum cadence of a high-resolution measurement at mass 32 every 45 minutes. In order to be able to make direct comparisons between species located on different masses, we restricted our data to measurements where the time difference between the measurements of two different species were less than 20 minutes. In addition, we discarded the first measurement following a data gap in the acquisition of longer than 2 hours (e.g., due to reaction wheel offloading), and only considered measurements when the instrument was pointing directly to the comet (nadir angle <10 deg). Lastly, we restricted our analysis to measurements that were taken within a distance of 40 km from the nucleus in order to exclude potential reactions between species in the coma. Reactions in the coma were minimal at heliocentric distances >~3 au during the time that a molecule propagates from the comet to Rosetta, and collisions were negligible within ~100 km (ref. 33).

Furthermore, the measurements show a 1/r 2 dependence as expected for a species that is not affected by any significant loss processes. In addition, the diurnal, longitudinal and latitudinal variations were clearly visible within this distance, indicating that the coma measurements by ROSINA/DFMS do reflect outgassing from the nucleus, and not that of a point source. This variation was visible to even larger spacecraft distances from the nucleus; thus, our approach weighing in the number of measurements vs. nucleus proximity was rather conservative. The above constraints on our data treatment then defined our pre-inbound and post-inbound equinox time intervals as: September 2014 through February 2015, and May 2016 through the end of mission in September 2016.

Statistical Analysis of the Derived Data: When looking at correlations for the shorter time intervals within the long (pre-and post-equinox) time periods, we selected the length of the short time periods making sure that:1. there was homoscedasticity in the data; i.e., the data showed fairly equal variability in the scatter plots; 2. no visibly discontinuous distribution of data was present; 3. multiple slopes were not visually present. These steps allowed us to rely on the derived Pearson correlation coefficients as indicators of a true link between two species, and ensured that we do not miss potential changes in outgassing. We have extensively tested the possibilities by splitting up the long time periods into several short ones and visually examining the scatter plots. We selected the data used in the short time intervals listed in Table 1 by ensuring all the above was true and the correlation was statistically significant at p < 0.01.

The Pearson correlation coefficient for each time interval was determined by performing

Monte Carlo simulations: we took 10,000 samples with a normal distribution within the boundaries of the uncertainties (discussed above) on each measurement point, and calculated the Pearson correlation coefficient 10,000 times. We report the mean of the resulting distribution of the correlation coefficients. Considering the overall small standard deviations, the correlations listed in Table 1 are robust.

Modeling gas trapping in solid H2O:

To determine whether O2 could be trapped and remain stable in the porous ice of the bulk nucleus, we modeled the interactions of H2O ice structures with O2 either adsorbed on the surface, or trapped in the voids or cracks of the H2O ice. We employed a method based on the first principle periodic density functional theory (DFT) to quantify these interactions between O2 and the icy bulk, using the Vienna ab initio simulation package 34-36 (VASP). VASP is well adapted to model bulk and surface ice structures interacting with atomic and molecular volatiles, and has been used many times for various planetary applications 30, 37-39 .

The stabilization energy (Estab) of O2 is defined as the difference between the energy of the whole system, [O2 interacting with the ice], and the sum of the energies of the isolated O2 and H2O ice:

Estab = E[ice+O2] -(Eice + EO2).
A negative value of Estab means repulsion, whereas a positive Estab indicates an effective trapping in that particular environment. In all cases of adsorption or insertion into the perfect crystalline ice or into the cavities found inside the bulk, the neighboring shells of water molecules are relaxed in order to allow O2 to fit into the void. At the same time, the outer shells are kept frozen to maintain a global bulk structure. In the cases when O2 inclusion is possible, that is when the hole is large enough (at least two neighboring H2O molecules extracted), the stabilization energy of O2 interacting with the ice is around 0.2-0.3 eV (Supplementary Table S1). While these values are lower than the ice cohesion energy (~0.55 eV), they stay in the order of that of an H2O dimer. Thus, O2 can effectively be stabilized in the ice cavities and pores, and the presence of O2 will not disturb the structure of the surrounding ice until the onset of H2O sublimation.

Taking into account the fact that several O2 molecules can be present in the same cavity, we consider the possibility that these O2 molecules could agglomerate, and eventually form dimers. According to high-level ab-initio calculations, the ground state of the O2 dimer is a singlet state contrary to the monomer that is a triplet 40 . To make space large enough for a dimer, a large cavity is needed. Therefore, we considered only those obtained following the loss of 4 H2O molecules. The order of magnitude of the values reported in Table 2 indicates that a couple of the O2 molecules can also easily be trapped in the same cavity within the ice either in the form of two O2 molecules in triplet state or as a proper O2 dimer in its ground singlet state. It can be noted that the dimer seems to be stabilized the most strongly of all; this might be a first step toward bigger structures, possibly a route to the construction of pure condensate. For comparison, some identical situations have been evaluated for the CO and CO2 molecules (see Supplementary Table S2). The strength of the interaction is systematically larger, but is within the same order of magnitude. It stays in the same order of magnitude. CO and CO2 stick slightly better than O2 but the difference seems too small to induce significantly different behaviors.

Data Availability: All data processed as described above are included in this published article and in a supplementary source file. The ROSINA/DFMS spectra used to derive count rates are available in the NASA Planetary Data System (PDS) and ESA Planetary Science Archive. AL-K conceptualized and led the study, performed the data analysis, interpreted the results and wrote the manuscript. OM contributed to the conceptualization and methodology of gas trapping, and the implications of the results for cometary origin. FP, YE and OO performed the chemical modeling of gas trapping in H2O ice and assisted in the interpretation. JIL contributed to the interpretation of the results and to the chemistry of ice trapping. SAF supported all steps of data analysis and assisted the correlation study. KEM contributed to the data analysis and interpretation. KJT and SMP developed the fitting routine, performed the spectral fits to the ROSINA/DFMS mass spectra and provided the time series files of the count rates of the analyzed species. All authors contributed to the discussion of the results and the interpretation, and the writing of the manuscript.

Competing interests:

The authors declare no competing interests. Table 1. Short-term correlations over the southern hemisphere before the pre-perihelion equinox (top half of the table above the line) and after the post-perihelion equinox (bottom half). Pearson correlation coefficients for shorter time intervals before the pre-perihelion equinox and after the postperihelion equinox. The coefficients were determined as the mean of 10,000 Monte Carlo simulations for each time interval, assuming a normal distribution within the measurement uncertainties (Methods). The standard deviations of the correlation coefficients are smaller than 0.02 for each time period. The dates in the table are color-coded to match the appropriate data in Figs. 1,

Fig. 4 .

 4 Fig. 4. O2 trapping in and release from the two distinct nucleus reservoirs (arbitrary

Correlation of O2 with H2O (left), CO2 (middle), and CO (right) prior to the pre- perihelion

  2 and 3. Time intervals after the postperihelion equinox are listed below the horizontal line in the middle of the table.

	Figure legends
	Fig. 1. Fig. 3. O2 correlations with the major cometary volatile species during low H2O
	production. The heliocentric distance for the corresponding time period was between 3.32-3.70
	au. The different colors correspond to different short time periods as shown in Table 1 (rows 10,
	Time interval 11, 12). The uncertainty on the correlation coefficients is < 0.01. Heliocentric distance (au) Average cometocentric distance (km) Sub-spacecraft latitude (°S) No. data points	O2-H2O	CO2 O2-	O2-CO
	1	09/27-10/12/2014 3.29-3.19	17.9	25.6 -56	128	0.94	0.51	0.69
	2	10/15-10/27/2014 3.17-3.10	9.9	25.1 -50.1	48	0.97	0.12	0.32

equinox. The data were taken between 3.29-2.43 au. The colored dots indicate the short time intervals listed in Table

1

(rows 1-6). The uncertainty on the correlation coefficients is < 0.01. Fig. 2. Correlations of O2 with H2O (left), CO2 (middle), and CO (right) after the postperihelion equinox. The data were taken between 2.91-3.70 au. The different colors designate different short time periods as indicated in

Table 1 (rows 7-12). The uncertainty on the correlation coefficients is < 0.01.

Acknowledgments

This work was supported by the National Aeronautics and Space Administration (NASA) through grant number 80NSSC18K1620 (AL-K, SAF, KEM, KJT and SMP), the Excellence Initiative of Aix-Marseille Université -A*Midex, a French "Investissements d'Avenir programme" AMX-21-IET-018 and CNES (OM), NASA grant NNX17AL71A (JIL), and CNRS-INSU national programs PCMI and PNP (FP, YE, OO).

Author contributions: