Dust-to-Gas and Refractory-to-Ice Mass Ratios of Comet 67P/Churyumov-Gerasimenko from Rosetta Observations

This chapter reviews the estimates of the dust-to-gas and refractory-to-ice mass ratios derived from Rosetta measurements in the lost materials and the nucleus of 67P/Churyumov-Gerasimenko, respectively. First, the measurements by Rosetta instruments are described, as well as relevant characteristics of 67P. The complex picture of the activity of 67P, with its extreme North-South seasonal asymmetry, is presented. Individual estimates of the dust-to-gas and refractory-to-ice mass ratios are then presented and compared, showing wide ranges of plausible values. Rosetta’s wealth of information suggests that estimates of the dust-to-gas mass ratio made in cometary comae at a single point in time may not be fully representative of the refractory-to-ice mass ratio within the cometary nuclei being observed.


Introduction
Comets are believed to be the least processed remnants of the early days of the outer Solar System because their low gravity prevented high internal pressure, and their birth place in the outer protosolar disk (with low mass densities and moderate Kepler velocities) only led to moderate modifications by collisions after their formation. Comets consist of a mixture of refractory materials in the classical sense (minerals), an organic "soup" that ranges from simple molecules more volatile than water ice, such as methane and methanol, to heavy and refractory molecules such as kerogen-like insoluble organic matter (IOM), and a number of frozen volatiles (ices) like H 2 O, CO 2 , CO, H 2 S, NH 3 , and multiple other minor species (Bockelée-Morvan et al. 2004;Mumma and Charnley 2011;. Comet nuclei have a very low density/high porosity (Weissman et al. 2004;Weissman and Lowry 2008;, and the few samples available on Earth that likely originate from comets (chondritic-porous interplanetary dust particles, ultra-carbonaceous Antarctica micro-meteorites), as well as grains collected in the coma of 81P/Wild by the Stardust mission, all exhibit a very fine-grained texture (Joswiak et al. 2017). These compositional and physical characteristics point to a rather slow and gentle accretion under cold conditions 100 K, i.e. at large heliocentric distances within the protosolar disk. The relative proportions of ices and refractory materials within comet nuclei (the refractory-to-ice mass ratio) is important to trace how and where these objects formed and evolved early on. So far it is unclear if ice and refractories are mixed on a microscopic scale, where ice forms an ice mantle around µm-sized dust grains as suggested for the condensation sequence in a molecular cloud or the interstellar medium, or if mixing takes place at a "macroscopic" scale, i.e. independent ice grains and refractory grains mixed together, most likely in the outer disk, or a combination of both. The refractory-to-ice ratio is one indirect clue to address this, in combination with the composition of outgassed materials and refractories. It may help discriminate between comet formation models (Weissman and Lowry 2008).
Comet formation processes are still a subject of intensive research and no fully selfconsistent model is available yet. Dust delivered to the protosolar nebula from the interstellar medium was sub-µm to µm-sized. To form comet parent bodies (planetesimals), cometary dust had to accrete over many orders in magnitude in size. Dust coagulation models are based on collisional properties derived from laboratory experiments, which have been recently summarized by Blum (2018). Most of these experiments have been performed using silicate dust aggregates, which bounce at collision speeds of v ∼ 1 cm/s and fragment at v ∼ 1 m/s. This would suggest that no silicate aggregates larger than mm-sized were able to grow via this process (Zsom et al. 2010). Experiments including water ice suggested that it is considerably stickier, as the icy aggregates only fragment at v ∼ 10 m/s (Aumatell and Wurm 2014;Gundlach and Blum 2014) consistent with earlier predictions of numerical models (Wada et al. 2011). However, more recent work suggests that the tensile strength of icy aggregates is lower than expected at low temperatures (Gundlach et al. 2018;Musiolik and Wurm 2019). Nevertheless, the potential impact of the refractory-to-ice ratio on the sticking properties of mixtures between ices and refractories has not been addressed by laboratory experiments to date.
Numerical experiments focused on global redistribution of dust in the protosolar disk showed that if the aggregate fragmentation threshold indeed depends on the refractory-to-ice ratio, a pile-up of solids inside of the water snow line could be triggered and planetesimals may preferentially form just outside of the snow line, where icy aggregates can grow and accumulate (Drążkowska and Alibert 2017;Schoonenberg and Ormel 2017). However, even if the icy dust sticks very efficiently, its direct growth to km-sized bodies seems unlikely, as the growing aggregates are removed from the outer parts of the disk by inward drift (Krijt et al. 2016;Homma and Nakamoto 2018). The most widely accepted scenario of planetesimal formation is a multifaceted process, in which the µm-sized dust first grows to mmto cm-sized aggregates that subsequently undergo so-called streaming instability, forming overdense filaments that collapse to km-sized, gravitationally bound objects (Drążkowska and Dullemond 2014;Johansen et al. 2014;Lorek et al. 2018). The refractory-to-ice ratio of comets formed in this scenario would reflect that of dust aggregates participating in the streaming instability, expected from one model in the range 3-9 (Lorek et al. 2016).
The refractory-to-ice ratio has a most direct impact on the thermal properties and activity of comet nuclei. Outgassing of cometary materials is driven by sublimation of ices (H 2 O and CO/CO 2 mainly), mostly from subsurface layers, which in turn accelerates particles ("dust") as they leave the nucleus. These "dust" particles may be only composed of refractory materials, or may also contain ice. In the remainder of the manuscript, we refer to ice-free particles as "dry dust". The relative proportions of ice and dust in the near-surface of the nucleus directly impact its thermal properties, as the filling of pores by water ice, either primordial or due to recondensation of sublimated ice from below, would increase the thermal conductivity of this environment, easing propagation of heat from the Sun to greater depths and mobilization of ices (Prialnik et al. 2004(Prialnik et al. , 2008. For comets observed to date, the dust-to-gas ratio that can be measured in the coma was assumed to be directly representative of the refractory-to-ice ratio within the nucleus. The dust-to-gas ratio in the coma of 1P/Halley was first determined as 1. 3-3 (McDonnell et al. 1991), although a re-evaluation of the dust size distribution (Levasseur-Regourd et al. 1999) suggests much larger dust-to-gas ratio values (Fulle et al. 2000). The dust-to-gas ratio of 9P/Tempel is greater than 1 (Küppers et al. 2005). However, these estimates were derived from observations acquired at a single time along their orbit. ESA's Rosetta mission (Glassmeier et al. 2007) escorted comet 67P/Churyumov-Gerasimenko (hereafter 67P) through 2014-2016, and was equipped with a comprehensive suite of instruments to survey the comet surface, its interior, its activity, and evolution over time. Thus, Rosetta provided the first opportunity to investigate dust-to-gas and refractory-to-ice mass ratios from up close, with multiple measurements, over an extended period of time. This review summarizes pertinent results obtained to date, their uncertainties, and potential ways to investigate this in the future.

General Characteristics
67P/Churyumov-Gerasimenko, named after the two astronomers who discovered it in 1969, is a Jupiter-Family Comet (JFC). As such, its elliptical orbit around the Sun (current orbital period 6.44 years, perihelion distance 1.24 AU, aphelion distance 5.68 AU) is controlled by gravitational interactions with the largest planet of our Solar System. Dynamical models retracing the history of 67P's orbit have shown that a close encounter with Jupiter in 1959 injected it into its current orbit, only 9 revolutions before the arrival of Rosetta (Krolikowska 2003;Groussin et al. 2007;Maquet 2015). However, these dynamical history computations also showed that the evolution of 67P's orbit is highly chaotic. This stochastic source of uncertainty prevents from determining whether 67P was injected into the inner Solar System for the first time after the 1959 encounter with Jupiter, or if it could have further evolved as a JFC in the inner Solar System before that, possibly with an even closer perihelion distance than now. Nevertheless, all dynamical models agree that 67P's orbit must have been controlled by Jupiter for thousands of years, hence suggesting that it is a relatively evolved comet.
Despite 67P having likely passed through the inner Solar System multiple times in its geologically recent past, the volatiles emitted by its nucleus are consistent with their original ices being formed from outer protosolar disk and perhaps presolar materials, see hereafter and the review by Hoppe et al. (2018). Indeed, 67P exhibits the highest D/H ratio measured to date in cometary water , suggesting significant incorporation of outer protosolar disk materials. The abundances of highly volatile molecules like N 2 (Rubin et al. 2015) and O 2 (Bieler et al. 2015a(Bieler et al. , 2015b, and noble gases like Ar , are consistent with formation at extremely low temperatures found in the Kuiper Belt and beyond. The presence of sulfur allotropes and other sulfur-bearing molecules ) is consistent with their formation via radiolytic processes in the presolar cloud and/or the outer protosolar disk, as shown in laboratory experiments (Mahjoub et al. 2017), prior to their incorporation in 67P. The isotopic composition of sulfur-bearing molecules ) furthermore suggests a significant contribution of materials extraneous to the solar nebula. Thus, 67P appears to still retain at present-day the primordial presolar/protosolar materials from which it formed.
The 67P nucleus, revealed by the Rosetta mission's Optical, Spectroscopic, Infrared Remote Imaging System, OSIRIS (Keller et al. 2007), has a striking bilobate structure  with a stable spin axis and strong seasonal characteristics. The high obliquity of 67P results in significant differences in insolation between the northern and southern "hemispheres" Jorda et al. 2016;). The southern regions experienced a relatively short summer around perihelion in August 2015 (with equinoxes on May 2015 and March 2016). Although the 67P spin orientation is stable, its rotation period changes over time owing to sublimation-induced torques (Kramer et al. 2018;Kramer and Läuter 2019). Hirabayashi et al. (2016) computed the stresses required to propagate a hundred-meters-long crack between the two lobes to the extent of splitting the nucleus, and found that the rotation period would need to reach values as low as ∼ 7 hrs. If the rotation period maintains its linear decrease with every orbit, it would take ∼ 15 orbits, almost 100 years, for 67P to reach the point where it may split into two (provided that no catastrophic events occur before that).
67P exhibits a low surface albedo ∼ 4% , is dominated by organic materials  and has an extremely high porosity Hérique et al. 2016;Pätzold et al. 2016Pätzold et al. , 2019. These compositional and structural characteristics result in a low thermal inertia Marshall et al. 2018), which promotes large temperature gradients over a few cm and extreme temperature swings following rapid changes in illumination. Little to no ice is directly detected on the surface, except: 1) in specific shadowed areas where ice can remain stable for extended periods of time (Filacchione et al. 2016a), 2) as surface frost associated to a diurnal cycle where ice is mobilized from the subsurface during the day and then is recondensed on the surface during the night , and 3) when the activity level increased enough towards perihelion for nucleus erosion to expose fresher materials .
The Philae lander was deployed on November 12, 2014, and bounced off its initial landing site Agilkia on a hard subsurface layer beneath ∼10-20 cm of loose regolith (Biele et al. 2015), which was imaged at high resolution by the ROLIS camera (Mottola et al. Schematic overview of measurements by Rosetta instruments that inform the parameters needed to derive the dust-to-gas and refractory-to-ice mass ratios of 67P 2015). At Abydos, Philae's final resting place on 67P, the SESAME package (Seidensticker et al. 2007) made observations compatible with a meter-scale layer of compact materials in the subsurface via SESAME-PP (Lethuillier et al. 2016) and SESAME-CASSE (Knapmeyer et al. 2018). The locally high thermal inertia of 85 ± 35 J m −2 K −1 s −0.5 ("high" by cometary standards only) measured around Abydos by the MUPUS-TM sensors (Spohn et al. 2015) suggests that this layer could be formed by sublimation and recondensation of water ice at depth within the pores, similarly to observations made in the KOmet-SImulationen (KOSI) laboratory experiments (Kochan et al. 1989;Grün et al. 1991).

Rosetta Measurements Relevant to Dust-to-Gas (Lost Materials) and Refractory-to-Ice (Nucleus) Mass Ratios
Rosetta was the first mission to Rendezvous with a Jupiter Family Comet and escort it from 3.6 AU pre-perihelion to 3.6 AU post-perihelion. It carried 11 instruments on its orbiter and deployed a fully instrumented lander on the surface of 67P. It was thus ideally equipped to study the surface, interior, and activity of 67P and their evolution through the active portion of this comet's orbit (Fig. 1). The Radio Science Investigation, RSI (Pätzold et al. 2007), measured total mass and mass loss during the perihelion passage from radio tracking of the Rosetta spacecraft. OSIRIS (Keller et al. 2007) constrained the total volume of the nucleus, from which the nucleus density could be derived, and surveyed activity and mass transport processes. The Comet Nucleus Sounding Experiment by Radiowave Transmission, CONSERT (Kofman et al. 2007), used electromagnetic waves to sound the interior and derive dielectric properties of the nucleus. The Microwave Instrument on the Rosetta Orbiter, MIRO (Gulkis et al. 2007), the Visual InfraRed Thermal Imaging Spectrometer, VIRTIS (Coradini et al. 2007), and the Alice UV spectrograph (Stern et al. 2007) determined remotely the production rates of water and other volatiles, and their integration over the course of the mission. The Rosetta Orbiter Spectrometer for Ion and Neutral Analysis, ROSINA (Balsiger et al. 2007), monitored gas density and volatile abundances at the location of the spacecraft throughout the mission. The Grain Impact Analyzer and Dust Accumulator, GIADA ), sensed dust grains and measured their velocity, cross-section, and mass. The COmetary Secondary Ion Mass Analyzer, COSIMA (Kissel et al. 2007), collected dust particles, imaged them and measured their composition using secondary ion mass spectrometry. The Micro-Imaging Dust Analysis System for the Rosetta Mission, MIDAS (Riedler et al. 2007), used atomic force microscopy to image the dust particles to sub-micron scale.
All these complementary measurements, illustrated in Fig. 1, provide various pieces of information needed to determine the dust-to-gas δ DG mass ratio in the coma and the refractory-to-ice δ RI mass ratio within the 67P nucleus. In the subsequent sections, we present the results and discuss sources of uncertainties related to the measurements themselves (calibration over long-term mission, potential sources of instrument bias) and the models that need to be used to derive these parameters from the data. Furthermore, Rosetta found that 67P's activity involves nucleus-scale mass transfer processes, see next section. These processes contribute to materials being emitted and redistributed around the nucleus. All these issues may affect the determination of the δ DG and δ RI mass ratios if they were not adequately considered.

Mass Transport Processes on 67P
The irregular shape and very high obliquity of 55°of 67P, combined with its highly elliptical orbit, are responsible for an asymmetric distribution of illumination from the Sun. This results in extreme seasonal variations between the Northern and Southern regions, with the Northern regions being subject to a long but mild summer around aphelion, while the Southern regions conversely undergo a brief but extremely intense summer around perihelion, e.g. . Smooth plains observed mainly in the northern 67P hemisphere are evidence of dust redistribution on the nucleus, referred to as "airfall" (Thomas et al. 2015) and explained as a dominant mass transfer from South to North occurring mainly around perihelion .  investigated the complex shape of the nucleus and predicted the impact of insolation over the period of the mission, providing inferences of potential water production rates and surface erosion. Fougère et al. (2016a) used a Direct Simulation Monte Carlo (DSMC) model to invert the activity distribution from ROSINA data, shape model, and illumination, and found that H 2 O seemed to originate predominantly from the Northern regions, in particular Hapi, in the early phases of the mission (pre-perihelion).
As 67P approached perihelion, subsequent OSIRIS observations of the nucleus at visible wavelengths suggested a general increase in water ice content, likely associated with 67P shedding off its dust mantle and exposing fresher materials in particular in the Southern regions ). This inference is consistent with previous comparisons between OSIRIS spectrophotometry data and VIRTIS-M IR spectra of exposed water ice (Filacchione et al. 2016a). These observations confirm strong erosion in the South during the perihelion passage.
The analysis of ROSINA/DFMS ) and MIRO (Biver et al. 2019) data throughout the Rosetta mission shows how the source locations of gases in the coma varied with seasons and insolation. The long-term monitoring of dust jets and their source location over the course of the mission (Shi et al. 2016;Vincent et al. 2016a;Schmitt et al. 2017;Shi et al. 2018) confirm that these features are driven by local illumination, and further supports the derivation of gas sources locations and their evolution. Violent outbursts have also been observed Vincent et al. 2016b;), however Schematic evolution in activity of 67P along its orbit. Top: evolution in sub-solar latitude of a sphere with the same orbital parameters and spin pole orientation as 67P to illustrate the extreme seasonal variability. Modified after . Bottom: Images of 67P's activity at locations along the orbit representative of the three periods mentioned (Credit: ESA/NAVCAM/OSIRIS). Note that the re-deposition of dm-size particles may occur within the perihelion window (1) and up to several months after perihelion ) within period (2). Period (3) consists of the Northern summer around aphelion and early during 67P's next approach to the Sun for particles deposited in the Hapi region. The period for reactivation of these large particles deposited elsewhere would depend on local illumination and their ice content these features likely originate from other processes, such as thermal cracking, rupture of a pressurized gas reservoir, and/or large-scale collapse as observed by .
Hapi is also suggested as a location where sublimation from airfall particles (Thomas et al. 2015) emitted near perihelion would contribute significantly to the activity in the Northern regions around aphelion, during the northern summer . Large dm-size "chunks" of 67P materials could be emitted from the Southern regions near perihelion . These boulders could be transported and redeposited in Hapi and elsewhere, retain a significant portion of their ice content, and be later on activated once the northern regions exit polar night and are illuminated during the next approach towards the Sun . The enhanced backscattering ratio observed in dust particles in the inner coma of 67P during the post-perihelion phase by  is interpreted as consistent with this scenario.
Thus, 67P offers a picture of the complex processes playing a role in cometary activity, with a heterogeneity that is both spatial and temporal in nature. Figure 2 illustrates the evolution in subsolar latitude with heliocentric distance (in the approximation where 67P is a sphere, thus should not be taken literally), and evolution of activity and history of mass transport in this scenario. Potential sources of uncertainty and their impact on this scenario are discussed in Sect. 4.2.
3 Dust-to-Gas and Refractory-to-Ice Mass Ratios of Comet 67P

Derivation of the Dust-to-Gas and Refractory-to-Ice Mass Ratios
The dust-to-gas mass ratio in the lost materials of 67P, δ DG , can be derived from the total mass lost by the nucleus (M Tot ) over the course of the Rosetta mission, (10.5 ± 3.4) × 10 9 kg (Pätzold et al. 2019) and the integrated amounts of volatiles that left the nucleus (M gas ) over the same timespan: Or, if considering measurements conducted at specific points in time t (from instruments onboard Rosetta or ground-or space-based observations), then in the absence of any fallback materials δ DG (t) would be related to the volatile and dust loss rates Q V and Q D at that time: The latter approach has been most commonly used to derive the δ DG (assumed constant) of comets to date from production rates at specific times, due to the observations available (Fulle et al. 2004;Sykes et al. 2004). However, the production rate approach is likely subject to large uncertainties at comets that exhibit strong seasonal variations in insolation, as δ DG could change over time.
In addition, one must note that most of these previous measurements at other comets considered that the gas phase was only composed of water vapor. Rosetta enabled detailed measurements of the composition of several additional volatiles in the coma. For an appropriate comparison with previous measurements, the dust-to-water δ H 2 O DG and dust-to-all-volatiles δ V DG ratios are treated separately hereafter.
In an ideal case where all dust particles in a given volume of the nucleus can be lifted by the sublimation of ices within that volume, and all these materials are lost, δ RI = δ DG = (M Tot − M gas )/M gas . However, deriving the refractory-to-ice mass ratio of the 67P nucleus, δ RI , is more difficult and prone to assumptions in the framework of a spatially and temporally varying activity. Rosetta showed that a substantial amount of materials emitted from the nucleus actually fall back onto the surface of 67P (Sect. 2.3).
The mass transfer processes observed by Rosetta at 67P strongly suggest a much more complex scenario, which is depicted in Fig. 3. We define the following sources of mass transfer (loss to space or fallback) around 67P: M gas PE is the mass of gas primarily emitted directly from the nucleus and lost; M (d+i) PE is the mass of dust (potentially still containing frozen volatiles) emitted directly then lost, which would after sublimation of those volatiles equate to a mass of refractory dust M d PE and a mass of gas in extended sources M ext PE ; FB is the mass of dust (and possibly ices) that falls back on the nucleus; M ext FB is the amount of gas emitted from large (∼dm) particles that fall back onto the nucleus but sublimate a portion of their internal ices during their journey in the inner coma; M gas RE is the mass of gas re-emitted subsequently from the sublimation of ices within the particles that have fallen back on the nucleus; M d RE is the mass of refractory materials subsequently emitted upon reactivation.

Fig. 3
Sketch of the quantities involved in mass transfer processes around 67P along its orbit. All these quantities contribute to the dust-to-gas δ DG and refractory-to-ice δ RI mass ratios determined in the lost materials and the nucleus, respectively. Depending on times in the mission and activity, the Rosetta spacecraft may be located inward or outward of the location where ice embedded in decimeter-size particles (coma boulders) would be completely sublimated, i.e. where dust particles (perhaps with ice) M (d+i) turn into dry dust M d and possible extended sources M ext , whether these particles are lost to space at primary emission (superscript "PE") or fall back on the nucleus and may be reactivated later (superscripts fallback "FB" and re-emitted "RE"). See text for details None of these individual parameters are directly accessible from measurements made by Rosetta. What can be assessed, albeit with rather large uncertainties, is the sum of all material lost during the mission from the comet: by considering the measured change in nucleus mass over the mission. In addition, the extent of fallback materials M (d+i) FB can be estimated from the observed transport of dm-size particles around perihelion ). The total mass loss and reasonable ranges of dust-to-gas and refractory-to-ice mass ratios suggest that fallback materials onto the nucleus may amount to multiple times the lost mass, see Sect. 3.2 (Pätzold et al. 2019). The total loss of volatiles (M gas from Eq. (1)) can also be determined by integrating the in-situ and remote sensing measurements made throughout the 2015 apparition of 67P, see Sect. 3.4.
Using production rates measured at discrete times in the mission to derive the refractoryto-ice mass ratios in the nucleus might be tempting, as suggested above for coma measurements, but this approach has the significant drawback that the respective mass contributions of fallback and re-emitted dust may take place at different times along the orbit of 67P. Thus, integrated mass losses in Eq. (3) cannot be simply replaced by production rates. However, as suggested in , one may do so near perihelion where most of the activity takes place, and then derive the amount of fallback materials to retrieve their refractory-toice mass fraction at that time.

Total Mass and Mass Loss Due to 2015 Perihelion Passage
67P's gravity field was measured from close distances in August-November 2014 by RSI, allowing to determine the total mass of the 67P nucleus: (9.982 ± 0.003) × 10 12 kg (Pätzold et al. 2016). Assuming a uniform nucleus, this mass and the 67P nucleus volume, first estimated by  and (Jorda et al. 2016), then refined by  to 18.56 ± 0.02 km 3 , yield a bulk density of 537.8 ± 0.6 kg m −3 (Pätzold et al. 2019).
Additional RSI measurements have been conducted in July-September 2016 (after the August 2015 perihelion passage), which allow to constrain the total mass lost by the 67P nucleus over the course of the Rosetta mission: (10.5 ± 3.4) × 10 9 kg (Pätzold et al. 2019), i.e. 0.1% of the nucleus mass. The bulk density of the nucleus post-perihelion, 537.3 ± 0.6 kg m −3 , remains within errors consistent with the value determined from pre-perihelion measurements.
Pätzold et al. (2019) discuss ranges of porosity and ice contents compatible with the bulk mass and density values from RSI data. For a range of compact dust material density from 2000 to 3500 kg/m3, the porosity varies between 65%-79% when the refractory-to-ice mass ratio δ RI for the nucleus lies in the range 3-7 (Pätzold et al. 2019). Values as low as 1 cannot be excluded but are not favored in their analysis. The nucleus is interpreted as a highly porous dusty body with low ice content.
By comparing the total mass loss with water production rate from ROSINA and MIRO (Sect. 3.4), Pätzold et al. (2019) suggested a dust-to-gas mass ratio of 0.5 to 1. From this dust-to-gas mass ratio and the range of refractory-to-ice mass ratios, they derived a range of plausible mass of fallback materials of 1.8 to 4.8 times the mass lost by 67P around the 2015 perihelion passage.

Dielectric Properties of the 67P Nucleus
CONSERT probed 67P's small lobe ("head") in the vicinity of Abydos up to a propagation length (along the line-of-sight Philae-Rosetta) of about 1 km and a penetration depth of about 100 m. The measured average dielectric permittivity ε = 1.27 ± 0.05 shows a highly porous nucleus .
The composition of the nucleus interior deduced from this low permittivity (Hérique et al. 2016) suggests an organic-rich refractory part of the comet, consisting of at least 75% volume fraction of organics (66% mass fraction) and of less than 25% silicates, which is consistent with surface observations from other Rosetta instruments (Sect. 2.1). This leads to a nucleus with an ice volume fraction ranging from 6% to 11%, a refractory volume fraction from 16% to 21% and a porosity from 73% to 76%. It is however worth noting a discrepancy between the organic matter abundance from CONSERT data and that measured in dust grains collected in the coma by COSIMA, which is around 45% by mass ), see Sect. 3.5. CONSERT measurements suggest that the refractory-to-ice mass ratio of 67P, δ RI , is larger than 3. This is an average value within the portion of the 67P nucleus probed by CON-SERT and cannot constrain local variation of the ratio (Hérique et al. 2019). It constitutes a lower limit, with all uncertainties included, and results from modelling. Large refractory-toice mass ratios are required to explain the low average permittivity measured by CONSERT: the constitutive materials must fit the 1.27 value and have a high enough density for small volume fractions to fit the 67P average density as derived from RSI. It is the case for some organics (ε ≈2, ρ ≈2), less for silicates (ε ≈5∼7, ρ ≈3.5), and definitively not the case for ices (ε ≈3.1, ρ ≈1 for H 2 O ice; ε ≈2.7, ρ ≈0.94 for a mixture of 71% H 2 O, 14% CO 2 and 15% CO in mass), as discussed in .

H 2 O
Water was directly observed throughout the mission by the MIRO, VIRTIS and ROSINA instruments. ROSINA measured local water number densities at the spacecraft, while the remote sensing instruments provided line-of-sight integrated information of the water distribution. The interpretation of these measurements is complex and the determination of water sources and total production rate requires extensive modelling (Sect. 4.1.3). We summarize in Table 1 estimates of the peak water production rate and total water mass loss obtained by these instruments. Hansen et al. (2016) analyzed ROSINA measurements from June 2014 to May 2016 using an empirical model of the water coma derived from the Direct Simulation Monte Carlo model (DMSC) of Fougere et al. (2016b), which was constrained by a restricted set of ROSINA data. Most recently, analyses of ROSINA data identified surface location of gas sources and derived production rates from August 2014 to September 2016, encompassing the entire comet close proximity operations phase of Rosetta Combi et al. 2020 Biver et al. (2019) uses a different approach. The ∼ 100 inner coma maps obtained from July 2014 to June 2016 were analyzed using a radiative transfer model, describing the outgassing pattern by a Gaussian shape. Both the line intensity and spectral shape were combined to derive the latitude of the water emitting region and total water production rate. Water production rates derived from VIRTIS observations are still sparse and published values rely on DSMC calculations using activity maps constrained by ROSINA data Fougère et al. 2016a;Fougere et al. 2016b).
Remote observations from spacecraft far from the comet also constrained water mass loss rates from 67P at various points in time. SOHO (SOlar Heliospheric Observatory) SWAN (Solar Wind ANisotropies) Lyman α observations were presented in Bertaux (2015), and by Shinnaka et al. (2017) using PROCYON/LAICA (PRoximate Object Close flYby with Optical Navigation/Lyman Alpha Imaging Camera) (Funase et al. 2015). Table 1 illustrates a disparity in water mass loss estimates, with results that appear incompatible when considering their error bars. These measurements are used to constrain the dust-to-water ratio in materials emitted by 67P in Sect. 3.4.3, after considering other volatiles, see next section. Further discussion of the approaches used to derive these values and their respective sources of uncertainties is provided in Sect. 4.1.3.

Other Volatiles
The analysis of ROSINA/DFMS data by  retrieved mission-integrated total mass losses for CO 2 , CO, and O 2 , amounting altogether to (10.1 ± 3.7) × 10 8 kg. These results are comparable to relative abundances to H 2 O determined previously at fixed points in time (Bieler et al. 2015a(Bieler et al. , 2015b). The measured relative abundances to water of sulfur-bearing molecules, 5% , hydrocarbon molecules, 5% , CH 4 , NH 3 and CH 3 OH together ∼2%, amount to (5.8 ± 1.8) × 10 8 kg, where the uncertainty is primarily driven by that on water production. Adding the integrated masses of CO 2 , CO, and O 2 ) yields a mass loss of volatiles other than water of (15.9 ± 5.5) × 10 8 kg. Hence, the total mass loss of volatiles from in-situ measurements is (6.4 ± 1.9) × 10 9 kg.
The recent analysis of MIRO data throughout the 67P rendezvous phase (Biver et al. 2019) included all compounds accessible at the fixed-tuned frequencies of the MIRO Chirp Transform Spectrometer: H 2 O, CH 3 OH, NH 3 , and CO. In addition, VIRTIS-H data allowed to determine the relative abundances to water of CO 2 , CH 4 , and OCS around perihelion and their variations (Bockelée-Morvan et al. 2016). CH 4 and OCS have very low abundances relative to H 2 O, thus would not contribute significantly to the total volatile losses. Conversely, the CO 2 abundance relative to H 2 O of ∼ 15%, with a large distinction between pre-perihelion and post-perihelion, makes it the second most abundant volatile at 67P (Bockelée-Morvan et al. 2016). By combining the MIRO and VIRTIS measurements, and including relative abundances of additional species to water as done above, Biver et al.
(2019) report a total mass loss of volatiles from remote sensing and MIRO map data only of (4.2 ± 0.2) × 10 9 kg, and a possible range (4.0-5.8) × 10 9 kg when including MIRO nadir observations. However, greater uncertainties are associated with retrievals from nadir observations, thus the high end of the range is deemed less likely. Biver et al. (2019) also derived a dust-to-gas ratio for 67P of 1.5 ± 0.8.
Analysis of Alice remote sensing UV data yielded much higher O 2 abundances relative to water (Keeney et al. 2017) than the ROSINA/DFMS analysis ), up to 60% O 2 /H 2 O. A more recent analysis of Alice stellar occultation data conversely yielded O 2 /H 2 O abundances comparable to those measured by ROSINA (Keeney et al. 2019). The source of these differences is not understood yet, although high-energy electron impact not accounted for in the original analyses has been suggested.

Dust-to-Water and Dust-to-Gas Ratios δ DG
We determine the dust-to-all-volatiles mass ratio δ DG V by combining the volatiles mass losses from in-situ measurements of (6.4 ± 1.9) × 10 9 kg, and from remote sensing of (4.2 ± 0.2) × 10 9 kg (MIRO maps) with a possible range (4.0-5.8) × 10 9 kg (Biver et al. 2019), with the total nucleus mass loss of (10.5 ± 3.4) × 10 9 kg (Pätzold et al. 2019). From these values, we first derive the dry dust mass loss in the two cases, M d (in-situ) of (4.1 ± 3.9) × 10 9 kg and M d (remote sensing) of (6.3 ± 3.4) × 10 9 kg. Volatiles loss from remote sensing (MIRO maps) yields a δ DG V (remote sensing) of 1.5 for average values and a range that spans 0.7-2.3 when accounting for 1-sigma uncertainties. Including nadir observations extends the low end of the range and yields 0.2-2.3. Volatiles loss from in-situ measurements yields an average δ DG V (in-situ) of 0.64 and a range 0.01-1.28. We then determine the dust-to-water mass ratio δ DG H 2 O by combining the water mass losses from in-situ measurements of (4.8 ± 1.5) × 10 9 kg ) and remote sensing of (2.4 ± 0.2) × 10 9 kg (Biver et al. 2019) with the dust mass loss as calculated above for in-situ and remote sensing measurements. Water mass loss from in-situ yields a δ DG H 2 O (in-situ) of 0.85 for average values and a range that spans 0.01-1.71 accounting for 1-sigma uncertainties. Similarly, water mass loss from remote sensing yields a δ DG H 2 O (remote sensing) of 2.6 for average values and a range that spans 1.2-4.0 when accounting for 1-sigma uncertainties.

Dust Particles Characteristics
Observations of dust particles by GIADA, COSIMA, and OSIRIS yield overlapping dust size distributions up to m-sizes (Fulle et al. 2016a(Fulle et al. , 2016b). The reader is referred to  for a synthesis of dust particle types, substructures, characteristics, and nomenclature.
GIADA measured directly the momentum and the speed of individual dust particles with particle sizes ranging from about 50 µm up to about 1 mm. Combining these two measurements, the mass of each particle was derived. The cross-section of each particle was derived from calibrated optical measurements , and together with the mass measurement, the particle density was derived Fulle et al. 2016aFulle et al. , 2016b). In addition, the GIADA microbalance system measured the integrated mass (1 reading every 300s) of particles <5 µm (Della Corte et al. 2019).
GIADA observations have revealed two distinct dust populations: a 'fluffy' one with size range 0.2-2.5 mm and of density <1 kg m −3 , and more compact particles 80-800 µm of density (1.9 ± 1.1) × 10 3 kg m −3 (Della Corte et al. 2015; . Following acquisition of data throughout the Rosetta comet operations at 67P, the bulk density of dust particles collected by GIADA was subsequently updated to 785 +520 −115 kg m −3 by , from which a volumetric nucleus composition of 54 ± 5% hydrocarbons, 22 ± 2% silicates, 4 ± 1% sulfides and 20 ± 8% ices was deduced. COSIMA measured directly the composition of dust particles via secondary ion mass spectrometry (SIMS) and found average relative proportions by weight of 45% for organic and 55% for inorganic matter, or 55% and 45% by volume, respectively, with about 30% variation between individual particles . Thus, the interpretations of GI-ADA data and the measurements made by COSIMA are compatible within their uncertainties.  report a detailed typology of dust particles, enabled by COSIMA/ COSISCOPE, identifying a variety of clusters or aggregates of particles in the range 10several 100 µm constituting 85% of their collected particles, while the remaining 15% are made up of >100 µm compact particles that did not fragment on impact with the collection targets. Some of these compact particles subsequently fragmented during SIMS analyses, owing to Lorentz forces caused by particle charging , while others were lifted and could not be analyzed.
A set of recent numerical simulations also suggested that part of the COSIMA collection is compatible with "fluffy" particle simulations, with either two populations with fractal dimensions <2 and about 2.6 and one impact velocity or one particle population with a fractal dimension of 2.6 and a spread of impact velocities on collection .
The particles collected by COSIMA were assumed ice-free, as the instrument was kept at about 10°C within the sample storage portion, and particles were analyzed days to weeks after their collection. The particles were illuminated in sequence with a red LED (640 nm) at phase angles between 72 and 84 degrees from two opposite directions. The particles showed a surprisingly diverse reflectance factor in the range 3-22% . The measured reflectance, in particular the comparison between left and right illumination, was explained by scattering centers inside the agglomerate volume and an optical mean free path in the 20-25 µm range. The required porosity depends on the size of the scattering centers and is estimated to be in the 50-90% range .
The aggregate nature of the particles from 67P extends below 10 µm size range according to measurements by MIDAS (Riedler et al. 2007;, and indicates a fractal-like structure in some particles . COSIMA and MIDAS dust observations suggest similarities with IDPs (Schulz et al. 2015;).

δ DG and δ RI from Dust Particles
GIADA data and the derived density of dust particles have enabled several determinations of dust-to-gas and refractory-to-ice mass ratios of 67P while considering various time periods over the course of the Rosetta mission. Early estimates of the dust-to-water δ DG H 2 O of 6 ± 2 and dust-to-all-volatiles δ DG V of 4 ± 2 were first derived from GIADA and MIRO measurements around 3.6 AU pre-perihelion in late 2014 . Updated dust and gas mass loss estimates from that time through perihelion suggested an average δ DG H 2 O of 6, and a minimum value of 5 and maximum of 10 or greater (Fulle et al. 2016a(Fulle et al. , 2016b. Considering GIADA data obtained throughout the comet operations phase of Rosetta,  then derived an average refractory-to-ice mass ratio δ RI of 7.5 and within the range of 3 or greater.
COSIMA measurements also enabled estimating the density of dust particles , from which the potential original ice content after collection and sublimation, assumed to be representative of the refractory-to-ice mass ratio within the nucleus, can be derived. The following analysis could not be made in that publication because the density of the 67P nucleus was only constrained later in the mission (Pätzold et al. 2016(Pätzold et al. , 2019.
In the COSIMA dust particle collection, all particles can be broken up by impact or electrostatic forces into elements, and can be considered as agglomerates of smaller elements shaped by hierarchical growth. The particles show sub-structures down to the instrument optical resolution limit of 14 µm/pixel. This indicates that the sub-structures, possibly again agglomerates, stick together with higher forces than the larger-scale agglomerates. To infer mechanical properties from the impact fragmentation,  used this model of agglomerates and inferred the strength limits for the fragmented and non-fragmented particles. The strength limits are determined by the binding force between sub-agglomerate structures or elements, and are around 1000 Pa for elements in the 10-40 µm size range.
A macroscopic porosity of 0.5 ± 0.1 is derived from observations of sub-structures (or elements of non-fragmented agglomerates) in the 60-300 µm size range. These elements are then assumed to be porous with the smallest solid unit of ∼ 0.2 µm diameter (Greenberg and Li 1999). This is formally described by their size-dependent filling factor ϕ ∝ r −0.4 , as implied by a hierarchical cascade on scales down to the solid grain. The fragmentation model constrains strength boundaries rather than void spaces.
Evaluating impact dynamics, size distribution of the fragments, and underlying hierarchical cascade of sub-structures and porosities, the density of the ice-free dust particles ρ D collected by COSIMA is in the range 100-400 kg m −3 ). This value is lower than the mean nucleus density, ρ 67P = 537.8 ± 0.6 kg m −3 (Pätzold et al. 2019).
Modeling the filling of the COSIMA dry dust particles with ice and deriving the putative refractory-to-ice mass ratio of these particles, assumed to be representative of the 67P nucleus, is relatively straightforward in that context. Assuming an average density of the monomers within dry dust particles ρ M , which are mixtures of silicates, sulfides, and organic materials as described above, the COSIMA-scale (microscopic) porosity of the dry particles is φ P = 1 − ρ P /ρ M , and the volume fraction (not to be confused with the filling factor) of compact dry dust monomers is f M = ρ P /ρ M . Given the density of compact water ice ρ I and the density of the nucleus ρ 67P , the volumetric fraction of porous ice f I that needs to be added to the pores to match the density of the 67P nucleus is: The refractory-to-ice ratio of the particles can then be derived in the scenario of ice filling porosity within dust particles. Indeed, the relative volumetric fractions of solid ice f I and of dust monomers f M that match the bulk density of the 67P nucleus are known: ρ 67P = f M ρ M +f I ρ I = ρ P +f I ρ I . So are their approximate densities, of ρ M ∼ 2500 kg m −3 and ρ I ∼ 920 kg m −3 . In this scenario where ice fills the porosity within dust particles, themselves consisting of aggregates of compact monomers up to the size of the nucleus (i.e. not considering any macroporosity), the refractory-to-ice mass ratio δ RI is: Adding ice to the ice-free porous dust particles, of density in the range 100-400 kg m −3 , results in a porosity of the icy particles in the range of 0.43 to 0.7, see Fig. 4. A direct transposition of these results to the 67P nucleus implies a refractoryto-ice mass ratio δ RI in the range 0.2 to 3, in stark contrast with the estimates derived with dust densities in the mm size range larger than the mean nucleus density (Fulle et al. 2016a(Fulle et al. , 2016b). The latter would then require further macroscopic porosity considerations to sizes much larger than the dust collected by COSIMA, in which case the ice content of the nucleus may not be reliably constrained from the ice-free dust particle density shown above.
The difference between the refractory-to-ice estimates from dust particles based on GI-ADA and COSIMA arises from the density of dust particles. GIADA data suggest dust particles denser than the nucleus itself (which then requires a substantial macroporosity), Fig. 4 Porosity of icy particles as derived from the ice-free particle density and the mean nucleus density (left), and refractory-to-ice mass ratio as function of the dry dust particle density and the mean nucleus density (right). The potential parameter ranges based on the COSIMA dust collection are indicated by the shaded areas while COSIMA data suggest dry dust particles less dense than the nucleus. In the latter case, ice is added until the density of dust particles prior to collection is equal to that of the 67P nucleus. Possible reasons for the discrepancy between GIADA and COSIMA dust density estimates are briefly mentioned in Levasseur-Regourd et al. (2018) and . Further discussion on potential collection biases and uncertainties associated to the density derivation principles can be found in Sect. 4.1.4.

Fallback Materials and Impact on Refractory-to-Ice Mass Ratio
Early estimates of dust mass loss were around 15 kg s −1 at 2.9 AU , rising by a few orders of magnitude during perihelion, where dust mass loss rates of (1.7 ± 0.9) × 10 3 kg s −1 (Fulle et al. 2016a(Fulle et al. , 2016b have been reported for boulders not all on radial trajectories. Subsequently, an estimate of up to 8.3 × 10 3 kg s −1 is found for boulders that are on radial trajectories ). These high dust production rates are such that the total nucleus mass loss measured by RSI corresponds to 71.5 and 14.6 days of activity at the Fulle et al. (2016aFulle et al. ( , 2016b and  rates, respectively. The above estimates neglect any ice or water vapor emission, thus are upper bounds on the duration of the activity peak that would be commensurate with these data. Given the commonly accepted overall activity scenario for 67P, see However, the abundances of the super volatile CO and the volatile CO 2 relative to water derived from the South pre-perihelion in May 2015 ), probably representative for the overall nucleus composition, are not higher than in the northern regions when Table 2 Summary of estimates of δ DG H 2 O , δ DG V , and δ RI for comet 67P. Values in parentheses are lower bound estimates that cannot be excluded but are not preferred in the original analyses. A "+" indicates that a higher upper bound cannot be excluded. "N /A" = not applicable. "inf." = results could be compatible with an infinitely high value, although a minimal fraction of ice is needed far from the Sun (Le Roy et al. 2015). This suggests that either ices more volatile than water have been retained in the fallback material, or the fallback layer in the north is not thick enough to dampen the release of CO and other volatiles from below. The high CO/H 2 O and CO 2 /H 2 O in the South early in the mission could then be explained by low water emission due to frigid polar night temperatures, as suggested from MIRO observations .
Fallback materials are likely dominated in mass by dm-size particles, observed in the inner coma of 67P ) and suggested to settle slowly during the postperihelion portion of the orbit ). This large size would allow them to preserve ices within their interior, and enable their reactivation later on ).  suggest that the inferred original composition of these "chunks" may be representative of the interior of the 67P nucleus and indicates a refractory-to-ice mass ratio of 3 or more. Table 2 summarizes published values of the dust-to-water, dust-to-gas, and refractory-toice(s) mass ratios, and adds new determinations for these parameters based on recent literature data considered in this article (see previous sections). It describes the methods used to determine these parameters and the period at which they have been conducted. The wide ranges of values of individual determinations, as well as the very restricted regions where independent determinations overlap, are quite striking. Figure 5 is a graphical representation of the data shown in Table 2, in an attempt to tease out further information from these estimates.

67P's Dust-to-Gas and Refractory-to-Ice Mass Ratio Estimates
The top portion of Fig. 5a presents dust-to-water mass ratios in the lost materials, while its bottom portion includes all volatiles. The total mass losses, integrated over the 2015 apparition of 67P, of H 2 O, other volatiles, and nucleus as a whole (Sect. 3.4), yield the lowest dust-to-water and dust-to-volatiles mass ratios. Figure 5b shows published values of the refractory-to-ice mass ratio within 67P, based on the analyses of GIADA and OSIRIS dust data , of RSI data (Pätzold et al. 2019), CONSERT data (Hérique et al. 2016), and the analysis of dry dust Summary of ranges of plausible values derived for the dust-to-gas δ DG (panel a) and the refractoryto-ice δ RI (panel b) mass ratios in 67P's lost materials and nucleus, respectively. In panel a) the mean value of dust-to-water or dust-to-all-volatiles mass ratios is indicated numerally, in addition to entire plausible ranges based on error bars of individual parameters. In panel b) the triangulated ends indicate possible lower values that are not favored in the interpretations. In both panels, the X-axis scale is limited at 15 for legibility. See text for details density from COSIMA data based on the morphology of collected particles ), see Sect. 3.5.2. As in the case of the dust-to-gas ratio, the individual estimates differ substantially.
The dust-to-all-volatiles mass ratio δ DG V is the mass ratio that would actually be expected to correspond most directly to the refractory-to-ices mass ratio within the nucleus of 67P, except if ice is largely incorporated in coma boulders that might act as extended sources of water. Indeed, it includes all the volatile species that contribute to both dielectric properties as determined by CONSERT and total mass loss around the 2015 perihelion passage measured by RSI.
Furthermore, an additional complication arises from the timing of the derived values. The 4 ± 2 value of δ DG V  was obtained early in the mission, and differs from the mission-integrated δ DG V determinations. The dust-to-gas ratio may have changed over the course of the mission (Fulle et al. 2016a(Fulle et al. , 2016b. And/or the activity of 67P and volatile-bearing fallback materials could result in different values/ranges for δ DG V and δ RI . Presently, there are two opposed interpretations of the refractory-to-ice mass ratio estimates shown in Table 2 and Fig. 5 within the community. One interpretation is that the uncertainties in the measurements themselves and in the steps needed to derive estimates of the refractory-to-ice mass ratio of Rosetta data analyzed to date are too large to draw a firm conclusion. Another interpretation is that a potential bias associated to the collection principle of COSIMA (see Sect. 4.1.4) justifies discarding those results and accepting a refractory-to-ice mass ratio of 3 or greater as consistent with all other Rosetta analyses to date. Such a minimum value would then be consistent with model predictions of 67P formation by gravitational collapse of pebbles driven by streaming instability in the protosolar disk (Lorek et al. 2016;Blum et al. 2017), albeit with the potential caveat that no predictions of refractory-to-ice mass ratios for alternate formation models have been published to date, to the best of our knowledge.

Uncertainties on Instrument Measurements and Derived Parameters
Each and every parameter involved in the estimation of δ DG and δ RI is derived from instrument measurements with their own individual statistical and systematic errors that may not be completely known. Furthermore, a significant a posteriori modeling effort is often required to derive the parameters needed to constrain the mass ratios. Therefore, uncertainties and approximations in the models used to analyze the data should be considered, as well as potential biases associated to the nature and origin of materials being measured. This translates into the wide ranges of plausible dust-to-gas and refractory-to-ice mass ratios reported in Table 2 and Fig. 5.

Total Mass and Mass Loss
The nucleus mass of 67P is amongst the best constrained values used in this work, because of the very low spacecraft velocity relative to 67P, high-precision X-band tracking, and extended stay at close distances that maximized spacecraft trajectory perturbations by the interactions of the gravity field of the comet with the spacecraft. The nucleus mass was determined both before and after the 2015 perihelion passage to a precision of 0.03% (Pätzold et al. 2016) and 0.01% (Pätzold et al. 2019), respectively. The mass loss due to the 2015 perihelion passage is 10.5 × 10 9 kg, or around 0.1% of the total nucleus mass (Pätzold et al. 2019), at an uncertainty of 32%. The total mass loss places constraints on the mass loss of volatiles determined by the other instruments. Fig. 6 Ternary diagram of the dust permittivity as a function of the volume percentages of ice, dust (i.e. carbon with silicates), and porosity, with the constraints on the average density (green lines), on the refractory to ice ratio (blues lines) and from CONSERT measurements (red lines). Refractory-Ice-Porosity volume fractions that meet the three constraints are located within the triangle defined by the upper green line, the right blue line, and the red line. Modified after Hérique et al. (2016)

Dielectric Properties
The composition estimation by CONSERT and, consequently, the δ RI of the nucleus, are based on the modeling of the effective permittivity from models of ice and dust fractions (Hérique et al. 2016). This model uses the following physical quantities as inputs.
Density: The modelled cometary material (ice, refractory and porosity) has to be compliant with both the average permittivity and the average density which is assumed to be ρ = 533 ± 6 kg m −3 , derived from mass measurements obtained by the Rosetta Radio Science Investigation (Pätzold et al. 2016) and the observed volume of 67P nucleus based on the 67P shape model . In the retrieval process, the lower limit of δ RI comes from the upper limit of the density (upper green line in Fig. 6).
Refractory-to-Ice ratio: The data retrieval approach developed in Hérique et al. (2016) takes a δ RI = 4 ± 2 (blue lines in Fig. 6) as input . Model calculations providing a good fit to the CONSERT data fall in the range of δ RI 3-6, which implies that for pure water ice, 3 is a lower limit. Values greater than 6 cannot be excluded as they have not been included in this analysis.
Ice fraction: The ice fraction considered consists of H 2 O, CO 2 and CO ices, with a maximum of 14% mass fraction of CO 2 (8% molecular abundance) and 15% CO (13% mol.), with permittivity ranging from 2.7 to 3.3 for density values of 931 to 1171 kg m −3 . Considering the ice fraction, the main uncertainties on the retrieved δ RI is from the CO 2 and CO contents: a larger fraction of CO 2 and CO, up to 30% each, would decrease the permittivity of the ice fraction, and consequently may be consistent with δRI values as low as ∼2.2.
The dielectric constant of such an ice-dust-porosity mixture remains compatible with that derived from CONSERT measurements.
The portion of the 67P nucleus probed by CONSERT appears absent of heterogeneities at the few-meters scale or more (Hérique et al. 2019). The measurement approach is not sensitive to local variability in δRI. CONSERT measurements only allow for estimating the average δRI value within the probed volume.
It must be noted that the analyses presented in this article assume that the whole 67P nucleus is homogeneous and that its fraction probed by CONSERT at the southern border of the Hatmehit region is truly representative of the nucleus. The composition of the nucleus surface and its interior appear homogeneous at large scales, however no Rosetta data can fully validate this assumption.

Uncertainties on Production Rates and Mass Loss Derivation
Deriving production rates and mass loss from measurements involves multiple modeling steps, typically with simplifying assumptions in order to keep the computation times manageable, and also calibration steps that will not be described here. Remote sensing instruments are sensitive to the molecular absorption or emission of photons at specific wavelengths (rotational modes of H 2 O, NH 3 , CH 3 OH, and CO for MIRO, and vibrational modes excited by fluorescence for VIRTIS) along their line of sight, thus the determination of production rates relies on both coma models and excitation models. ROSINA measured in-situ, at the location of the spacecraft, the total gas density directly using its COma Pressure Sensor (COPS), and the nature and relative proportions of coma constituents ionized within the instrument for detection and counting.
Coma models are used to predict the entire outgassing pattern from the nucleus and determine the local density, temperature, and velocity of the gas within the portions of the coma probed along the line of sight for remote sensing instruments (Biver et al. 2019), or at the spacecraft for ROSINA (Bieler et al. 2015a(Bieler et al. , 2015bKramer et al. 2017;. These models typically use as basis either the physical interactions between gas molecules, such as the Direct Simulation Monte Carlo approach (Fougère et al. 2016a;Zakharov et al. 2018) or a simplified collosionless coma as used to invert source location from ROSINA/COPS (Kramer et al. 2017) and ROSINA/DFMS data , or the geometric characteristics of the coma as observed (Biver et al. , 2019. For all approaches, the 67P coma is found to be extremely asymmetric, with water activity originating primarily from the illuminated portions of the nucleus at any point in time, and additionally with latitudinal variability for some species (Hässig et al. 2015). The asymmetry of the coma is three dimensional, and although 1D simulations can mimic some of this complexity by adding plume contributions with very good performance (Biver et al. 2019), the actual structure of the coma may add unforeseen uncertainties.
The analysis of remote sensing spectral data requires excitation models, as well as radiative transfer models. Combined, they enable the retrieval of the column density along the line of sight of the molecule(s) of interest. When lines are optically thick, the retrievals are highly dependent on assumptions on the gas local density, temperature and velocity field. These parameters have been evaluated at best when interpreting MIRO observations (Biver et al. 2019). As far as possible, optically thin lines have been used to infer column densities and production rates from MIRO and VIRTIS data Biver et al. 2019).
In-situ volatile data acquired by ROSINA is also dependent on coma models and laboratory measurements. Coma models are needed to extrapolate the density of the molecule(s) of interest throughout the coma from measurements made locally (Fougère et al. 2016a) to derive production rates. Laboratory measurements are used to constrain the fragmentation pattern and ionization efficiency of these molecules (and/or their fragments) in ROSINA/DFMS and ROSINA/RTOF, and thus their individual contribution to specific peaks observed within the mass spectra obtained. Note that the analysis of ROSINA/COPS data does not require such laboratory efforts, as the COPS sensor solely measures gas density.
Lastly, all the individual sets of measurements need to be integrated over time to derive total volatile mass losses (Biver et al. 2019;, as was done in the initial empirical estimate by Hansen et al. (2016). This step adds further uncertainty to the results despite the high cadence of mapping measurements by MIRO and the regular acquisition of ROSINA data that have been planned throughout the comet operations phase Vallat et al. 2017).

Extended Water Sources: Hidden Ice in Escaping Boulders?
Extended sources may be associated to ice-bearing grains that start sublimating in the coma once they are beyond Rosetta's position (within the first few tens -hundreds of km around the nucleus). Such sources of gas could explain some OSIRIS dust observations Gicquel et al. 2016). Extended sources of water are abundant at 103P/Hartley, a hyperactive comet (A' Hearn et al. 2011), and have been suggested on other comets such as 1P/Halley where distant sources of specific compounds (not water) were identified (Eberhardt 1999).
In the particular case of 67P, Odin (Frisk et al. 2003) conducted space-based observations of the 557 GHz rotational line of H 2 O at 67P, which is amongst the lines measured by MIRO, in November 2015. The line was not detected after 40 hrs of integration time and the 3-sigma noise yields an upper limit in production rate of <3.3 × 10 27 molecules/s (Biver et al. 2019). From the MIRO measurements at the same time, the water production rate was (2.2 ± 0.1) × 10 27 molecules/s (Biver et al. 2019). Furthermore, the peak water production rate derived from PROCYON Lyman α measurements (Shinnaka et al. 2017) is also about 50% higher than determined by MIRO, see Table 1.
Even though extended sources are not required to explain the measurements, due to the measurement uncertainties it cannot be excluded that up to ∼ 50% of the water production (and presumably other volatiles) could originate from coma boulders that are escaping and do not have time to fully release their volatiles while in the vicinity of the comet explored by Rosetta. Since these volatiles would be undetected by Rosetta, this contribution of volatiles would be part of the dry dust mass loss calculated in Sect. 3.4.3. As a consequence, the dustto-water δ DG H 2 O and dust-to-all-volatiles δ DG V mass ratios determined in the inner coma by Rosetta (Table 2 and Fig. 5) could be somewhat overestimated.

Production Rates and Mass Loss of Dust
Determining the production rate and mass loss of dust requires several elements: a good knowledge of the grain size distribution, composition and/or density, overall distribution around the nucleus to expand from local measurements to spatially-integrated production rates, and their evolution over time. The latter includes fallback materials and their uncertainties, which are discussed separately in Sect. 4.2.2.
GIADA and COSIMA used different approaches to constrain the density of dust particles. GIADA measured directly scattered light, velocity, and momentum of dust particles of a few µm to mm diameter. Particle mass and density was then derived from these measurements. COSIMA/COSISCOPE relied on optical observations to infer the properties of the particles. The density of particles collected by GIADA is on order of 800 kg m −3 , while the density of dust particles after fragmentation within COSIMA is significantly lower, in the range 100-400 kg m −3 .
COSIMA might have collected particles from one or two potential families, one with a fractal dimension of 2.5-3 and a spread in impact velocity, and the other <2 if all particles impacts with the same velocity . A large fraction of the dust entering COSIMA is reflected, owing to a mass transfer efficiency estimated to 1 to 10% depending on size and velocity of incoming particles (Ellerbroek et al. 2019). This potential collection bias of COSIMA based on a collection or detection area suggests that a fractionation of the dust populations might be applicable. The GIADA measurement principle, which catches dust particles as they fly into the instrument, is not affected by this type of bias.
A coma model is then needed to extrapolate to the scale of the entire coma. Dust is lifted by gas drag, therefore all uncertainties on outgassing of volatiles propagate to the dust lifting, while additional sources of uncertainties exist due to the size, shape, density of the dust, potential ice content that is not measurable optically nor by any in-situ method conducted near room temperature, and processes through which it is lifted (see Vincent et al. 2019, this book).
For instance, this approach was used by  to derive the total mass loss over 53 days around perihelion from parallax imaging of ∼1-30 cm dust particles. The latter observations were restricted to 6 km from the Rosetta spacecraft in the few occurrences where such particles were imaged, implying that these dm-scale dust particles were on escape trajectories (radial velocity >1 m/s, 200 km distance from the nucleus). The extrapolation of the individual measurements over the 53-day period of observation yields a greater integrated amount of ejected (and lost) materials than the total 67P mass loss determined by RSI. Two lines of interpretation can explain this discrepancy. A fraction of these dm-scale coma particles, outside of the terminator plane where Rosetta was located, could fall back on the nucleus . Or the sources of uncertainties associated to the derivation of volume, or subsequent extrapolation to particle mass then total mass, may be poorly constrained. These sources of uncertainties include distance/size from the parallax measurement , density (a very high value of 1000 kg m −3 was used), albedo, and statistical uncertainties. Note that the density does not affect the subsequent study by , and the assumed albedo of these particles is consistent with that measured by .

Thermal Processing, Erosion, and Water Outgassing
Often the dust-to-volatile ratio of a comet is measured in its coma, e.g.  and therein. But how representative is it of the refractory-to-ice ratio in the nucleus? In other words, what can we learn from coma measurements for the nucleus composition when accepting that: 1) The amount of hidden ice in the dm-scale escaping or falling back dust particles is widely unknown (Sect. 4.1.3). It is mostly inaccessible for remote sensing instruments and difficult to model because of the lack of detailed information on parameters needed. 2) The escaping materials may have been altered by thermal processing and may not represent the original bulk properties.   c Corresponding to a dust/ice mass ratio of 2 at 120 K (Winter and Saari 1969;Herman and Weissman 1987) (Groussin et al. 2019, this issue). Higher dust/ice ratios would yield lower specific heat values and hence higher conductivity, resulting in larger skin depth values

Thermal Properties
Thermal and mechanical models can evaluate both the loss of volatiles in the escaping materials and the consequences of activity for the refractory-ice-ratio at the surface of the nucleus. For modeling such heat transfer processes, one has to know the thermal and mechanical properties of the materials involved. Rosetta instruments (MIRO, VIRTIS, OSIRIS, RSI, MUPUS, and SESAME) provided some physical characteristics that are related to the thermal inertia: density, strength, diurnal period, and thermal conductivity (for a review, see Marshall et al. (2018) and , this issue). The thermal inertia was derived from Rosetta measurements with considerable error bars. It controls the temperature profile at the surface and the thermal skin depth that characterizes the penetration depth of the diurnal and seasonal heat waves, respectively, e.g. Prialnik et al. (2004). Table 3 offers simple calculations of the thermal conductivity of the nucleus material and the corresponding thermal skin depths for lower and upper limits of the most plausible range of thermal inertias (Groussin et al. 2019 this issue). The diurnal skin depth is extremely low, no more than a few cm, as shown from MIRO continuum observations , while the seasonal skin depth can be up to a few meters. Typically, the heat wave propagates to ∼2-3 skin depths, as the skin depth only corresponds to a 1/e attenuation of the amplitude of temperature variations as compared to the surface.

Nucleus Erosion
These estimates of thermal processing length scale can then be compared to the length scale of nucleus erosion. Neglecting uncertainties on reported values for simplicity, the total mass loss of 67P around the 2015 perihelion passage is 10.5 × 10 9 kg (Pätzold et al. 2019), the nucleus volume is 18.56 km 3 and its surface area is 51.7 km 2 ). If we assume the erosion to be the same everywhere on the surface, the thickness of eroded materials is 37.7 cm. Non-uniform erosion due to the extreme seasonal asymmetry may further increase the thickness of this layer in the southern regions. Since a considerable amount of material that is lifted off later falls back to the surface, a near-surface eroded/redistributed layer is likely on order of ∼ a meter thick, and perhaps more in the Southern regions (about 4 m suggested by . Figure 7 compares the skin depths with the thickness of the erosion layer for several values of thermal inertia and percentage of fallback materials. The higher the inertia and the lower the fallback are, the more seasonal processing of the nucleus affects the near-surface materials that are later on measured in the coma. Seasonal processing is of the same order as nucleus erosion. Therefore, even the dm-size particles found in the coma have likely been already processed to some extent, perhaps even prior to their ejection. This eventuality has been incorporated in the recent analysis of coma boulders ejection timescale and redeposition, and their implications for a possibly high (>3) refractory-to-ice mass ratio .

Refractory-to-Ice Mass Ratio, Dust Mantle Thickness, and Outgassing
As illustrated in Table 2 and Fig. 5, and previously suggested , the dustto-gas ratio measured in the lost materials could differ substantially from the refractory-toice ratio within the 67P nucleus. If indeed the dust-to-gas ratio of the escaping materials is continuously lower than the refractory-to-ice in the nucleus, refractories would become progressively enriched at the surface of the nucleus either by losing volatiles or by re-deposition of fallback materials processed in the coma. This would lead to the formation of a dry lag deposit that could dampen activity, at least in some regions. Then, counter processes such as outbursts excavating deep-seated fresh materials, or extreme heterogeneities in distribution of the erosion and re-deposition of fallback materials, or mass wasting processes, would be necessary to sustain the activity over geologic timescales.
Detailed thermophysical modeling efforts further investigated the joint influence of the refractory-to-ice mass ratio and the thickness of a dry surface mantle on predicted erosion and water outgassing. A study assuming mm-size pebbles, comprised of a top dehydrated single layer over a dust-ice mixture with a δ RI of 7.5, found that such a near-surface structure is compatible with multiple Rosetta and Philae measurements, and supports the interpretation that 67P formed through streaming instability (Blum et al. 2017). Model calculations with dust mantle thicknesses of 5 and 10 mm over a dust-ice mixture with δ RI of 10 and 100 suggest that the water activity of 67P early on, well before perihelion, can be matched with very high δ RI values (Hu et al. 2017a). In a follow-on study, the analysis of OSIRIS images for erosion rates in specific regions and dust cover re-deposition suggest that extremely high δ RI values up to 100 cannot be excluded, and that much lower δ RI as low as 4 may also be consistent with the results (Hu et al. 2017b). Computations using the complex Rome thermophysical model (Capria et al. 2017), which couples orbital evolution, starting composition, thermal properties, and mass balance for loss of volatiles and gas drag for loss of refractories, were run for two example locations that received intense insolation around the 2015 perihelion and δ RI values of 4 and 10, including in the Anhur region where VIRTIS-M IR data suggested the temporary exposure of carbon dioxide ice (Filacchione et al. 2016b). These calculations suggest that a multi-decimeter erosion could occur around perihelion and preserve almost unaltered materials at relatively shallow depths <1 m after erosion (Capria et al. 2017).
Two aspects of the studies listed above are noteworthy. First, the refractory-to-ice mass ratio is an input parameter to the models, thus cannot be fully constrained by them unless inversion schemes were used to estimate sensitivity to, and uncertainties on, input parameters (which these studies did not investigate). Second, these works relied on published δ RI values available at the time , all of which were of 3 and above.
A recent study  constrains the activity of 67P by using both the ROSINA data of the water production and the non-gravitational forces derived from the spin change. They found that a complex pattern of the ice surface fraction that changes with time and surface region explains the measurements best. Further studies are needed to investigate the sensitivity of thermophysical models to their input parameters, to determine whether they could be used to reliably constrain δ RI , and to investigate if low δ RI values commensurate with the mission-integrated δ DG can be excluded.

Contribution of Re-Activation of Deposited Materials to Activity
The evidence for North-South asymmetry in morphology and dust cover seems naturally related to seasonal variability in activity, as discussed in Sect. 2.3. It is an unexpected feature of 67P, with high obliquity of its spin axis together with its elliptic orbit and the specific positions of the nodes between its orbital and equatorial planes. There are however several aspects that need confirmation to validate the proposed scenario, and its implication of a possibly extremely high δ RI suggested by .
The main line of evidence for the scenario illustrated in Fig. 2 is a combination of factors: morphology consistent with fallback of material (originating from the Southern regions at perihelion) in the Northern regions, as well as predicted fallback in the North greater than mass loss of volatiles and dust associated to native activity, and lastly activity during the Northern summer originating presumably from Hapi where this fallback may dominate. However, this combination of factors could be at least in part circumstantial.
It must be noted that the coma models used to invert ROSINA DFMS data into activity maps (Fougère et al. 2016a) do not have unique solutions (in particular because each measurement is at one specific location within the coma at a time, and coverage can only be partial due to spacecraft trajectories). Zakharov et al. (2018) have shown that multiple other solutions to the activity distribution can yield the results obtained by ROSINA, some of these solutions having little to no H 2 O originating from the Hapi region itself. Furthermore, the surface activity maps recently derived from ROSINA/DFMS data  show a widely varying distribution of activity, with Hapi and potential fallback materials only partially contributing to the activity in the North.
In addition, Vincent et al. (2016a) have suggested that jets observed in the northern regions during their summer are mostly consistent with activity occurring along cliffs and walls, where mass wasting would expose fresher materials and their ice content get sublimated fairly rapidly to form dust jets. A similar scenario has also been put forth to explain the large Ashwan cliff collapse ). Nevertheless, some dust jets in the Northern regions clearly originate from the smooth areas within Hapi (Shi et al. 2018). Thus, the proposed re-activation of fallback material as a key contributor to the Northern activity is consistent with some observations, but may not be the only mechanism involved. Lastly, the longevity of different ices within dm-scale particles that would be transported from the South to the North is advocated for . Nevertheless, it remains to be fully evaluated, as function of the thermal and mechanical properties of these materials, the volatility of the ice species, how long their transport requires, how they rotate in the coma and are illuminated, how they lose their ices, and what conditions they experience following re-deposition.

Future Research Using Rosetta Data
Despite more than two years spent in orbit around 67P, the analysis of Rosetta data to date only allows to determine the dust-to-gas and refractory-to-ice mass ratios with very low precision, see Fig. 5. Further analysis of Rosetta data may help better constrain some of the parameters involved in these estimates. A few examples follow.
1) A quantification of the distribution of erosion and the possible extent of fallback materials from OSIRIS data could shed a new light on lateral transport processes on 67P. First erosion estimates have been conducted in selected areas from shadows cast by boulders and cliffs to estimate height, and from excavation of cavities (Fornasier et al. 2019;Hasselmann et al. 2019). In addition, high-resolution local Digital Terrain Models pre-and post-perihelion can help measure surface changes. The combination of Multi-resolution StereoPhotoClinometry by Deformation MSPCD (Capanna et al. 2015) with stereophotogrammetry shape modeling , can achieve much higher resolution than currently available. These new approaches will likely improve the quantification of the mass transport processes in the near future.
2) Improving our knowledge of the history and thermal properties of the dm-scale particles before ejection/during flight/after re-deposition, as well as the relative contributions of fallback materials and fresher materials exposed by mass wasting to the 67P activity, may provide additional constraints on the original δ RI of these materials.
3) An extensive set of thermophysical calculations that encompasses the entire surface of 67P, and uses inversion to constrain the dust mantle thickness and the refractory-to-ice ratio over extreme ranges of possible input values, would be a necessary step towards understanding how sensitive these models are to their input parameters. They could in particular test the hypothesis that δ DG and δ RI differ largely, whether at single points in time and/or over extended periods of time.

Potential Future Instruments and Missions
Another factor that likely contributes to the low precision of the dust-to-gas and refractoryto-ice mass ratio estimates is the lack of direct observations of the deep interior and shallow subsurface of the nucleus. Onboard Rosetta, CONSERT probed a limited part of 67P's small lobe, while RSI did not conclusively find density differences between the two lobes. Direct observation of materials from the deep interior would be crucial to characterize the composition of pristine material and then to access directly δ RI . The direct observation of the shallow subsurface is also crucial to better understand comet activity processes including re-deposition and re-condensation.
Radar is the most mature technique to probe both the pristine interior and the shallow subsurface. Two different instruments or channels can achieve these objectives: 1) A low frequency radar to probe the deep interior, which could be a monostatic radar like MARSIS onboard Mars Express (Picardi et al. 2005) or bistatic radar with a Lander like CONSERT onboard Rosetta (Kofman et al. 1998). 2) A higher frequency radar to probe the first few tens of meters of the nucleus with decimeter resolution, in order to detect layers and undergrounds structures like HFR developed in the frame of AIDA mission (Michel et al. 2016;Hérique et al. 2018).
For example, two monostatic radars have been proposed to instrument the Castalia mission concept to Main Belt comet 133P/Elst-Pizzaro, proposed to ESA for the recent M5 cosmic vision round (Snodgrass et al. 2018).
Future comet missions carrying such radar instruments could better understand the vertical structure of the first ten to tens of meters of the subsurface with high vertical resolution. This would help better constrain the composition of the nucleus and its gradients over the thickness most affected by cometary activity (Fig. 7), which is necessary to fully address the δ RI issue of 67P and other comets. In addition, better knowledge of the fine-scale structure of the most active layer of comet nuclei would likely enhance our understanding of cometary activity.
Lastly, the ideal way to determine the refractory-to-ice ratio within the nucleus of 67P and other comets would obviously be to analyze directly representative materials sampled at depth within the nucleus. Rosetta was intended to be a sample return mission at its inception (following the 1986 flyby of 1P/Halley by Giotto), but the concept was modified into bringing a laboratory to 67P because of the technological challenges associated with sampling, returning, and curating cryogenic samples.
The NRC Planetary Decadal survey report "Visions and Voyages" for the decade 2013-2022 advocated that a cryogenic comet sample return mission would have by far the highest scientific return for the Small Bodies community. That report also acknowledged the lack of maturity for the acquisition, return and curation of a cryogenic comet sample, then advocated for the present decade for a non-cryogenic comet surface sample return mission.
Cryogenic sampling and containment systems are being developed as part of other technology development efforts, for example as suggested for the Triple-F mission concept (Küppers et al. 2009). A cryogenic comet sample return mission equipped with instruments capable of assessing heterogeneities within the nucleus may be able to identify a most volatile-rich and porous, least processed region at depth. If the deep interior is homogeneous, a sample may be representative and inform the refractory-to-ice mass ratio. If the deep interios is heterogeneous at the scale of the sample, this may not enable a determination of the refractory-to-ice ratio of the nucleus to a precision better than Rosetta, however the sample would retain the textural information that holds the clues as to processes through which it formed in the outer protosolar disk. More than 30 years after the flyby of 1P/Halley by Giotto (Reinhard 1986), further efforts by space agencies around the world are needed to turn a cryogenic comet sample return mission concept into a reality.

Conclusions
Rosetta data obtained at 67P enabled multiple estimates of the dust-to-gas mass ratio in the coma and the refractory-to-ice mass ratio within the nucleus. These estimates are overall of low precision, because of the uncertainties associated to individual measurements, the numerous models that need to be used to derive global production rates and mass loss from measurements made at one location and one point in time (hence large extrapolations both spatially and over time), and complications arising from the redistribution of materials around the nucleus and the poor knowledge of their composition and ice content.
The mission-integrated dust-to-gas ratio δ DG , based on total mass loss of the nucleus and the total mass loss of volatiles, is probably the best constrained parameter to date. A δ DG value of 2.3 or lower encompasses both remote sensing and in-situ determinations (considering the 1-sigma upper limit from the 1.5 ± 0.8 remote sensing determination, and although one cannot discard at this time the potential for systematic errors to affect this result).
However, it appears difficult to place firm lower or upper limits on the refractory-toice mass ratio. Opposing views are expressed to date on this topic within the community. Some argue that the data and analyses presently available have reached the level of maturity necessary to draw conclusions, while others caution that the numerous sources of uncertainties associated to the derivation of the refractory-to-ice mass ratio (retrievals, activity scenario, extent and composition of fallback materials, along with putative heterogeneities in composition within the nucleus and variability in comet activity through time) may not be understood and quantified well enough yet to support such conclusions.
It has been previously suggested that the dust-to-gas mass ratio measured in the coma may deviate significantly from the refractory-to-ice mass ratio in the nucleus, at least at a single point in time. One would expect that the integrated mass losses of dust and ice would help bound the refractory-to-ice mass ratio in the nucleus, but that is not the case at present.
Upcoming analyses of Rosetta data may shed a new light on the extent of fallback materials and erosion, and bring new quantified information on lateral transport processes. In turn, these new pieces of information may provide new constraints on the refractory-to-ice mass ratio within the 67P nucleus.
The most important lesson from Rosetta on the topic of dust-to-gas and refractory-to-ice mass ratios is probably that dust-to-gas mass ratio estimates from observations of cometary comae at limited points in time may not be representative of the refractory-to-ice mass ratio within the observed nuclei. This emphasizes the need for instrumentation tailored to investigating the interior and near-surface of comet nuclei and their contribution to cometary activity, and for future space missions to comets. vided the ESTRACK and DSN ground station networks during the data acquisition periods. This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 686709. This work was supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number 16.0008-2. The opinions expressed and arguments employed herein do not necessarily reflect the official view of the Swiss Government. This work was supported in part by the H2020 MiARD project. Part of this work has been conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract to NASA. MC acknowledges support from the US Rosetta Project, which manages the NASA contributions to Rosetta. Government sponsorship acknowledged.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.