Constraining models of activity on comet 67P/Churyumov-Gerasimenko with Rosetta trajectory, rotation, and water production measurements
Nicholas Attree, Laurent Jorda, Olivier Groussin, Stefano Mottola, Nick Thomas, Yann Brouet, Ekkehard Kührt, Martin Knapmeyer, Frank Preusker, Frank Scholten, Jorg Knollenberg, Stubbe Hviid, Paul Hartogh, Rafael Rodrigo
AAstronomy & Astrophysics manuscript no. preprint c (cid:13)
ESO 2019January 10, 2019
Constraining models of activity on comet67P/Churyumov-Gerasimenko with Rosetta trajectory,rotation, and water production measurements
N. Attree , , L. Jorda , O. Groussin , S. Mottola , N. Thomas , Y. Brouet , E. Kührt , M. Knapmeyer ,F. Preusker , F. Scholten , J. Knollenberg , S. Hviid , P. Hartogh , and R. Rodrigo Aix Marseille Univ, CNRS, CNES, Laboratoire d’Astrophysique de Marseille, Marseille, France Earth and Planetary Observation Centre, Faculty of Natural Sciences, University of Stirling, UK e-mail: [email protected] Physikalisches Institut, Universität Bern, Sidlerstrasse 5, 3012 Berne, Switzerland Deutsches Zentrum für Luft-und Raumfahrt (DLR), Institut für Planetenforschung, Rutherfordstraße 2, 12489 Berlin,Germany Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany International Space Science Institute, Hallerstrasse 6, 3012 Bern, SwitzerlandJanuary 10, 2019
ABSTRACT
Aims.
We use four observational data sets, mainly from the Rosetta mission, to constrain the activity pattern of thenucleus of comet 67P/Churyumov-Gerasimenko.
Methods.
We develop a numerical model that computes the production rate and non-gravitational acceleration of thenucleus of comet 67P as a function of time, taking into account its complex shape with a shape model reconstructedfrom OSIRIS imagery. We use this model to fit three observational data sets: the trajectory data from flight dynamics;the rotation state, as reconstructed from OSIRIS imagery; and the water production measurements from ROSINA, of67P. The two key parameters of our model, adjusted to fit the three data sets all together, are the activity pattern andthe momentum transfer efficiency (i.e., the so-called “ η parameter” of the non-gravitational forces). Results.
We find an activity pattern able to successfully reproduce the three data sets simultaneously. The fitted activitypattern exhibits two main features: a higher effective active fraction in two southern super-regions ( ∼ %) outsideperihelion compared to the northern ones ( < %), and a drastic rise of the effective active fraction of the southernregions ( ∼ − %) around perihelion. We interpret the time-varying southern effective active fraction by cyclicformation and removal of a dust mantle in these regions. Our analysis supports moderate values of the momentumtransfer coefficient η in the range . − . ; values η ≤ . or η ≥ . degrade significantly the fit to the three data sets.Our conclusions reinforce the idea that seasonal effects linked to the orientation of the spin axis play a key role in theformation and evolution of dust mantles, and in turn largely control the temporal variations of the gas flux. Key words. comets: general, comets: individual (Churyumov-Gerasimenko), planets and satellites: dynamical evolutionand stability
1. Introduction
The sublimation of ices when a comet is injected from itsreservoir to the inner solar system triggers the emissionof molecules. This outgassing produces in turn a reactionforce which can accelerate the comet nucleus in theopposite direction. The perturbing effect of cometaryactivity on the trajectory of comets has been establishedin the 1950’s in the pioneering work by Whipple (1950).At that time, it was clear that most comet trajectorieswere affected by a significant nongravitational acceleration(hereafter “NGA”) linked to their activity around perihe-lion (Marsden 1968). Shortly after, a theoretical modeldescribing the nongravitational force (hereafter “NGF”)produced by the sublimation of ices was establishedby Marsden et al. (1973). This model was based on asimple function describing the heliocentric dependence ofthe sublimation of water ice of a comet, combined with constant scaling parameters A , A and A describing theamplitude of the NGA along the three components (radial,transverse, normal) in the orbital frame of the comet. Themodel has been modified by Yeomans & Chodas (1989) toincorporate an asymmetric term used to describe the shiftof the maximum of activity with respect to perihelion. It isworth mentioning here that these simple models are still inuse nowadays to fit astrometric measurements, and henceto describe cometary orbits.More sophisticated models have been introducedsince then. Images acquired during the flyby of comet1P/Halley by Giotto showed narrow dust jets (Kelleret al. 1987), leading to the idea that the activity may beconfined to localised areas. This led to a new model ofnongravitational acceleration (Sekanina 1993) in which theoutgassing originates from discrete areas at the surfaceof a rotating nucleus. Szutowicz (2000) used Sekanina’s Article number, page 1 of 12 a r X i v : . [ a s t r o - ph . E P ] J a n &A proofs: manuscript no. preprint model to fit astrometric measurements of comet 43P/Wolf-Harrington obtained during nine perihelion passages. Shealso compared the modelled production rate with visuallightcurves used as a proxy for the activity of the comet(Szutowicz 2000). Recently, Maquet et al. (2012) revisitedSekanina’s approach with a model in which the activityis parametrised by the surface fraction of exposed waterice in “latitudinal bands” at the surface of an ellipsoidalnucleus. More accurate ground-based measurements aswell as space missions to comets offered the opportunity toincorporate new physical processes in models of the NGA.Rickman (1989) used the change of the orbital period ofcomet 1P/Halley caused by the NGA to extract its massand density. The detailed description of the local outflowvelocity incorporated a “local momentum transfer coeffi-cient” ζ l , originally called η in the improved descriptionof the solid-gas interface introduced by Crifo (1987). Thiscoefficient represents the fraction of the emitted gas’smomentum (dependent on its thermal velocity) whichneeds to be considered in the calculation of the momentumtransfer, and thus of the NGF (see Crifo 1987). Davidsson& Gutiérrez (2004) used 2D thermal modelling includingthermal inertia, self-shadowing, self-heating and an activitypattern to fit the NGA of comet 19P/Borrelly. They couldretrieve the mass of the nucleus and constrain the directionof the spin axis. The method was also applied to comets67P/Churyumov-Gerasimenko (Davidsson & Gutiérrez2005), 81P/Wild 2 (Davidsson & Gutiérrez 2006), and9P/Tempel 1 (Davidsson et al. 2007), all targets of spacemissions.The exploration of comet 1P/Halley in 1986 also yieldedthe discovery of the non-principal axis rotation of this comet(Samarasinha & A’Hearn 1991). Since then, the torque ofthe NGF was thus identified as the major effect responsi-ble for changes of the rotational parameters (Samarasinhaet al. 2004a, and references therein). Its modelling is re-quired to understand the observed change in the rotationalparameters of cometary nuclei, as well as the apparitionof non-principal axis rotations. Changes in the spin pe-riod of several comets have been detected from the anal-ysis of lightcurves (Mueller & Ferrin 1996; Gutiérrez et al.2003; Samarasinha et al. 2004b; Drahus & Waniak 2006;Knight et al. 2011; Bodewits et al. 2018). Predictions ofthe expected change in the direction of the spin axis andspin period for comets 81P/Wild 2 (Gutiérrez & Davidsson2007) and 67P/Churyumov-Gerasimenko (hereafter “67P”,Gutiérrez et al. 2005) have been made. Recently, a changein the spin period of comet 67P has been clearly identifiedbetween the 2009 and 2015 perihelion passages (Mottolaet al. 2014) based on early images acquired by the OSIRIScamera aboard Rosetta. Keller et al. (2015) showed thatthe spin period variation curve of 67P is controlled at firstorder by the bilobate shape of the nucleus.Thanks to its long journey accompanying comet 67P,Rosetta provides a unique chance to record measure-ments of most parameters involved in the modelling ofthe NGF. The mass of the comet has been retrieved byPätzold et al. (2016) from the radio science experimentaboard the spacecraft. The shape has been retrieved from astereophotogrammetric analysis of a subset of OSIRIS im-ages (Preusker et al. 2017), leading to an accurate knowl-edge of the moments of inertia. The activity pattern hasbeen constrained in the early phase of the mission by Marschall et al. (2016, 2017) from ROSINA measurementsin the form of “effective active fractions” associated to ge-ological regions. The total water production rate has beenconstrained from ROSINA measurements, complementedby ground-based observations (Hansen et al. 2016). Finally,the spin period has been monitored throughout the missionby ESA’s flight dynamics and OSIRIS teams before (Jordaet al. 2016) and after (Kramer et al. 2019) perihelion. Thelatter (Kramer et al. 2019) modelled the temporal evolu-tion of the rotational parameters, comparing it with themeasured change in spin period and spin axis orientation.The aim of this article is to try to reproduce three data setsderived by several Rosetta instruments (see section 2) witha model of the NGF (see section 3) in order to retrieve: (i) the local effective active fraction and its temporal vari-ations around perihelion, and (ii) a recommended value forthe momentum transfer coefficient η (see section 4). Theresults will be discussed in section 5, together with recom-mendations for the calculation of NGF of other comets.
2. Observational constraints
In this section, we describe the observational data, mainlyobtained by the Rosetta spacecraft, to which we will at-tempt to fit our NGF model.
The total water production rate is an important constraintfor any model of cometary activity, and a significant ef-fort has been made to measure it for 67P. As summarisedby Hansen et al. (2016), a number of different instrumentshave all been used and comparing and synthesising theirresults in non-trivial. Here, we use the ROSINA (RosettaOrbiter Spectrometer for Ion and Neutral Analysis) mea-surements, empirically corrected for spacecraft position, asour observed data points, O Q . Hansen et al. (2016) describehow estimating the uncertainty in these data, σ Q , is difficultso we use the bounds on the power-law, fitted with helio-centric distance to the inbound ROSINA data, as given intheir Table 2.We point out that ROSINA data are inferred from lo-cal measurements in the coma of 67P. Around perihelion,the spacecraft was located at a distance of
200 km fromthe nucleus, making it difficult to infer a total productionrate, whilst ground-based observations (see, for example,Bertaux 2015) have also suggested a variation in peak pro-duction between perihelion passages. Production rate esti-mates from Rosetta’s line-of-sight instruments MIRO andVIRTIS are also generally lower that the Hansen et al.(2016) results. These are important caveats to bear in mindwhen interpreting our results. Finally, the possibility thatsublimating icy grains are emitted by the nucleus is notconsidered in this article, and neither are other gas species,such as CO , CO and O . Other species represent < ofthe gas number density at perihelion (Hansen et al. 2016)and their production curves are not as well constrained aswater. Therefore, for this study we choose to focus on water,which is the primary driver of non-gravitational forces. Article number, page 2 of 12. Attree et al.: Activity models of comet 67P/Churyumov–Gerasimenko constrained by Rosetta
Outgassing produces a back-reaction force, and a resultingnon-gravitational acceleration (NGA), on cometary nucleiwhich affects their trajectory in a measurable way. For 67P,the nucleus trajectory has been reconstructed by the flightdynamics team of ESA, by a combination of radio-trackingof the Rosetta spacecraft from Earth and optical navigationof the comet relative to it, and is available in the formof NASA SPICE kernels (Acton 1996). Therefore, the 3Dposition of the comet in a heliocentric reference frame hasan accuracy greatest in the Earth-comet range direction( R , claimed accuracy of ∼
10 m ) and much less in theperpendicular (cross-track) directions (claimed accuracy of ∼
100 km ).Theoretically, the NGA resulting from outgassing couldbe directly extracted from the residuals between this mea-sured trajectory and a modelled gravitational orbit (takinginto account general relativity and the gravitational accel-erations of all major planets, Pluto and the most massiveasteroids). During the course of this work, however, it wasdiscovered (and later confirmed by ESAC; B. Grieger, per-sonal communication) that the reconstructed comet andspacecraft trajectories contain a series of discontinuities, atwhich the objects’ positions vary over hundreds of metresto several kilometres in an instantaneous time, making theabove method impossible. Within the orbital segments be-tween these ‘jumps’, there is no difference, to machine preci-sion, between our own orbital integrations (see Sect. 3.2 be-low) and the reconstructed trajectory, demonstrating thatit is a purely gravitational solution. The jumps occur at theboundaries between the integration segments that make upthe reconstructed orbit, and represent the offset in the ob-jects’ positions over the course of each segment due to theun-modelled non-gravitational acceleration. Unfortunately,it proved impossible to extract the NGAs directly from thejumps themselves as they contain, not only, the NGA ef-fect but also the typical uncertainty in the state vector atthe start of each segment, which is of a comparable mag-nitude. The jumps therefore have random magnitudes withtime (Fig. 1) and must be considered an additional sourceof noise in the uncertainty in the measured positions.Despite these issues, the reconstructed kernels remaina good description of the comet’s trajectory over orbitaltimescales, and within the limits of accuracy of the typicaljump size. Therefore it is still possible to compare thesemeasurements with a model of the orbit, including a ther-mal outgassing NGA as well as N-body gravitational inter-actions, in order to constrain said model. To do so, we usethe magnitude of the comet-to-Earth-centre range as ourobservable O R , since the jump sizes are smallest in this di-rection (Fig. 1). Considering the jumps to be a source ofrandom error, we conservatively estimate the uncertaintyin O R as σ R ≈ km. Back-reaction from outgassing not only produces a net ac-celeration on the nucleus but can also, depending upon itsshape, produce a net torque, altering its rotation state. Therotational parameters of 67P, including its spin rate overtime, ω ( t ) , has been measured as part of the reconstructionof its 3D shape from OSIRIS images (Jorda et al. 2016). −400 −200 0 200 400 Days since Perihelion −60−40−20020406080 J m p S i z e ( k m ) −3 −2 −1 0 1 2 3 J mp Size (km) N m b e r zyxR Fig. 1.
Discontinuities identified in the position of comet 67Pfrom the SPICE kernels, in ( x, y, z ) heliocentric J2000 coordi-nates and Earth-comet range ( R ). On the left: as a function oftime and on the right as a histogram of sizes. −400 −300 −200 −100 0 100 200 300 400 Days from Perihelion −0.50−0.250.000.250.500.751.00 T o r q e ( N m ) Fig. 2.
Observed torque, derived from the 67P rotation statefrom OSIRIS measurements, and a smoothed cubic spline fit tothe data. The grey region represents the RMS of the residualsbetween the two.
After verifying that the cross terms relating to the an-gular velocities along the first and second principal axes ofthe comet are negligible, changes in spin rate, ˙ ω z , can bedirectly related to the z component of the torque ( τ z , ina body-fixed frame) by Eq.(1), where I z = 1 . × kg m is the third (largest) moment of inertia derived fromthe shape model assuming a constant density of
538 kg m − (Preusker et al. 2017). τ z ≈ I z ˙ ω z , (1)Differentiating the observed ω ( t ) by time exacerbates mea-surement uncertainties so that the produced torque curvebecomes extremely noisy. For comparing with our simula-tions below, we therefore smooth the data by fitting it toa cubic spline, as shown in Fig. 2. Our fitted spline is thenused as the torque observable, O τ , with an assumed un-certainty equal to the root mean squared residuals of thederived minus smoothed data, σ τ = 575000 N . m (shown ingrey bounds in Fig. 2). Article number, page 3 of 12 &A proofs: manuscript no. preprint
3. Modelling
A thermal model is required to compute the temperature ofthe sublimating layer (assumed to be at the surface on thenucleus), from which we can derive the non-gravitationalforces acting on the nucleus. Our thermal model takes intoaccount solar insolation, surface thermal emission, subli-mation of water ice, projected shadows and self-heating.We use a decimated version of the 67P shape model called
SHAP7 (Preusker et al. 2017), with
124 938 facets. Thetemperature is computed for each facet of the shape model,36 times per rotation (i.e., every 1240 s), and 70 times (i.e.,every ∼ days) over the 2 years of the Rosetta mission,to ensure a good temporal coverage. At each time step,the distance to the Sun, the orientation of each facet rela-tive to the Sun, and the projected shadows, are computedusing the OASIS software (Jorda et al. 2010). Due to thelarge number of facets (>100 000) and time steps (>2500),heat conduction is neglected in the thermal model for nu-merical reasons. To test this assumption, we computed theproduction rate (Eq. 4) and acceleration (Eq. 5) of a spher-ical nucleus, at perihelion (where the torque is maximum),for two cases: a thermal inertia of 0 and 100 J/m /K/s . .The production rate and the acceleration are ∼ /K/s . (compared to anull thermal inertia), and the direction of the accelerationvector differs by less than 3 deg. Neglecting the thermalinertia therefore appears reasonable compared to the datauncertainties (e.g., production rates).The surface energy balance of the thermal model is givenby Eq. (2), for a facet with index i , where A b = 0 . is theBond albedo at 480 nm (Fornasier et al. 2015), F sun = 1370 W/m is the solar constant, z i is the zenithal angle, r h the heliocentric distance, SH i is the self-heating given byEq. (3), (cid:15) = 0 . is the assumed infrared emissivity, T i isthe surface temperature, f i is the fraction of water ice (inour case, either 0 for pure dust or 1 for pure ice, see below), α = 0 . accounts for the recondensation of water ice onthe surface (Crifo 1987), L = 2 . × J/Kg is the latentheat of sublimation of water ice at 200 K, and Z i is thewater sublimation rate given by Eq. (4). (1 − A b ) F sun cos z i r + SH i = (cid:15)σT i + f i (1 − α ) LZ i ( T i ) (2)In Eq. (3), the self-heating SH i is the sum of the infraredflux coming from all the n facets with index j that see facet i , where S j is the surface of facet j , θ j the angle between thenormal to facet j and the vector joining facets i and j , θ i theangle between the normal to facet i and the vector joiningfacets i and j , and d ij is the distance between facets i and j .This formalism is similar to that of Gutiérrez et al. (2001).To look for which facets are seeing each others, we usedan algorithm developed at Laboratoire d’Astrophysique deMarseille based on ray tracing and hierarchical search. Thecomputation of the viewing factors ( S j cos θ j cos θ i ) / ( πd ij ) is purely geometric and depends only on the shape model:it is therefore only performed once at the beginning. SH i = n (cid:88) j =0 (cid:15)σT j S j cos θ j cos θ i πd ij (3) In Eq. (4), the sublimation mass flow rate per facet is calcu-lated with the molar mass M = 0 . kg, the two constants A = 3 . × Pa and B = 6162 K for water (Fanale &Salvail 1984), and the gas constant R = 8 . J K − mol − . Z i = Ae − B/T ice (cid:114) M πRT ice (4)Finally, to compute the non-gravitational forces (Sect. 3.2),we need the sublimation rate (Eq. (4)) and the gas velocity(Eq. (6)), which both depend on a different temperature:the temperature of water ice T ice for Z i , and the temper-ature of dust T dust for v i . We therefore run our thermalmodel (Eqs. (2) to (4)) with two extreme cases: [1] with f i = 0 , which corresponds to a pure dust model, to com-pute T dust , and [2] with f i = 1 , which corresponds to a purewater ice model, to compute T ice . The reaction force vector per facet is then calculated basedon this mass flow rate, and the total acceleration is thesum over all facets divided by the comet mass ( M P =9 . × kg; Pätzold et al. 2016): F i = − ηx i Z i v i S i , a NG = (cid:80) i F i M P . (5)where η is the momentum transfer coefficient (Crifo1987; Rickman 1989), which we here assume to be con-stant across the comet. S i is the surface area of each facet(in the direction of its normal) and x i is its effective activefraction. Mass flow rate is calculated with Eq. (4) and thegas velocity is taken as the thermal velocity v i = (cid:114) RT dust πM , (6)assuming equilibrium with the surface grey body tempera-ture, i.e. that of the dust from run [1] of the thermal model.This is the upper limit that the gas can reasonably reach,meaning our non-gravitational force will be on the high endof estimation and our effective active fractions are lower lim-its. Calculated dust temperatures range between ∼ − K (ice temperatures are limited by the sublimation to ∼ K), leading to thermal velocities of ∼ − m s − .Water production and torque are likewise summed overthe surface Q = (cid:88) i x i Z i S i . τ NG = (cid:88) i τ i , (7)where torque per facet is the vector product of eachforce vector with its radius vector to the centre of mass( r i ). This can also be expressed as the magnitude of theforce multiplied by a “torque efficiency” (see Keller et al.2015): τ i = r i × F i = F i ( r i × ˆ S i ) . (8)The z component of the total net torque can then com-pared with the observations, using Eq. (1). It is advanta-geous to use the torque efficiency formalism since this vectoris in the body-fixed frame, and can be calculated once at Article number, page 4 of 12. Attree et al.: Activity models of comet 67P/Churyumov–Gerasimenko constrained by Rosetta +X+Z+Y -Y-Z-X
Fig. 3.
Torque efficiency in metres, a factor determined by thegeometry as defined in Eqn. 8. The comet’s rotation axis is inthe + z direction through the ‘neck’ region in the centre. the beginning of the simulation run, rather than being re-calculated each time. Torque efficiency is also a useful wayof visualising the effects of differing spatial distributionsin activity, which will be important later during the opti-misation. Mapping torque efficiency onto the shape model(Fig. 3) shows how local variations in topography combinewith large scale orientations of regions, varying the effectsof activity on the comet’s rotation across its surface (com-pare with Fig. 1 in Keller et al. 2015, which uses an older,incomplete shape model that is missing the southern hemi-sphere).Non-gravitational forces are calculated for each of the ×
36 = 2520 runs of the thermal model and the relevantquantities (force and torque vectors and summed water pro-duction) are averaged over a day. This produces smoothlytime-varying curves which can be inspected at any time ofinterest, using bilinear interpolation, which we refer to asour model solution. For comparison with the observed wa-ter production rate and torque, we can simply evaluate ourmodel at the time of each measurement, producing C Q and C τ , and directly compare.For the trajectory, however, a full N-body integrationmust be performed and the resulting position comparedat each time ( C R ). To do this we use the open-source REBOU N D code (Rein & Liu 2012), complete with fullgeneral relativistic corrections (Newhall et al. 1983) as im-plemented by the REBOU N Dx extension package . Allthe major planets are included, as well as Ceres, Pallas, http://rebound.readthedocs.io/en/latest/ http://reboundx.readthedocs.io/en/latest/index.html and Vesta, and are initialised with their positions in theJ2000 ecliptic coordinate system according to the sameephemerides used in the Rosetta trajectory reconstructions(NASA/JPL solar system solution DE405; Standish 1998).An additional particle representing 67P is initialised withits position given by SPICE. The system is then integratedforward in time, using the IAS15 integrator (Rein & Spiegel2015) and the standard equations of motion, with the ad-dition of a custom acceleration, a NG , for 67P, provided byour model. The modelled comet begins to diverge from themeasured positions and at each time of interest we directlycompare the magnitudes of the computed and measuredcomet-to-Earth-centre ranges, R (the most accurate partof the trajectory as described in sect. 2.2 above). In order to constrain the unknown parameters in ourmodel, such as the the surface active fraction, we performa bounded least-squares fit to the data using the doglegoptimisation routine (Voglis & Lagaris 2018) provided inthe scientific Python package. The optimisation proceeds,attempting to minimise the standard objective function
Obj = N (cid:88) j =0 (cid:18) O j − C j σ j (cid:19) , (9)with observed minus computed residuals, O − C , and ob-servation uncertainties, σ , at each time-step, j , up to thetotal N . The root mean squared residuals are then RM S = (cid:112) Obj/N .In our case we have three separate datasets to fit to ( R , Q and τ ), a multi-objective optimisation problem, and wetherefore use a linearly-scaled combination of the three togenerate a combined objective function. The term insidethe brackets in Eq. (9) then becomes the sum of the threeterms λ R (cid:18) O Rj − C Rj σ Rj (cid:19) , for < j ≤ N R ,λ Q (cid:18) O Qj − C Qj σ Qj (cid:19) , for N R < j ≤ N R + N Q ,λ τ (cid:18) O τj − C τj σ τj (cid:19) , for N R + N Q < j ≤ N, (10)respectively, where N = N R + N Q + N τ runs over thecombined number of points in all three datasets, and the λ scaling coefficients are variables, which themselves mustbe optimised in order to give the desired weighting to eachdataset. We set λ R = 1 and scale the other lambdas relativeto it, by bootstrapping from the residuals of a preliminaryrun, to give roughly equal weighting to all three datasets.Due to the very small relative errors for range ( σ R /R ∼ ± km / AU), this results in large values for the other lambdas( λ Q ∼ λ τ ∼ ) to give equal weighting.For each optimisation, we fix the coupling factor, η , ata constant value and parametrise the model in terms ofeffective active fraction, x , across the surface. The effectsof η on the best-fit model can then be studied indepen-dently. To begin with, we use the division of 67P’s surfaceinto five ‘super regions’, as performed by Marschall et al.(2016, 2017) in fitting ROSINA/COPS and OSIRIS data, Article number, page 5 of 12 &A proofs: manuscript no. preprint +Z -Z
Fig. 4.
Map of the 5 ‘super regions’ defined by Marschall et al.(2016, 2017) and used in our solution A optimisation. and start the optimisation with initial values from their re-sults. This divides the comet into: a southern hemisphereregion, an equatorial region, a region covering the base ofthe body and top of the head, Hathor, and Hapi, as shown inFig. 4 below. The fitting routine then proceeds to optimisethese five free parameters, subject to the lower and upperboundaries of zero and one, i.e. − active fraction.As described below, we also use more detailed parameteri-sations, including all 26 regions of Thomas et al. (2015) for26 free parameters.
4. Results
Before optimising our activity model we first test the N-body component by calculating the comet trajectory overthe course of the Rosetta mission with no NGAs applied,and with the classic NGA parametrisation based on themodel of Marsden et al. (1973). This model computes thethree components of a NG (radial, along-orbit and normal-to-orbit) based on a power-law with heliocentric distance,and three scaling parameters A , , found by a formal best-fit to the orbit for each comet. We use the A , , values forthe 2010 apparition of 67P from ground-observations givenby NASA/JPL Horizons ephemerides .Figure 5 shows a plot of the observed range plus theresiduals to both models. The RMS residuals are 1029 kmand 675 km, respectively, showing the order of magnitudeof the accuracy of the Marsden et al. (1973) model withground-based observations of roughly arcsecond accuracy.This stands as a baseline, against which we can check ourown activity model. Figures 6 –8 show the residuals to our best-fit solution usingthe five super regions defined by Marschall et al. (2016,2017), which we refer to as our solution A.The RMS range residuals, of 163 km (see Table 1 forall results), represent a significant improvement over the ∼ km of the purely gravitational solution and the ∼ km of the ground-based solution, demonstrating thesignificance of our NGA/N-body model. However, the waterproduction curve is clearly not a good fit, failing to repro-duce both the peak production rate as well as overestimat-ing the production far from perihelion. Likewise, the torquecurve is an extremely bad fit, failing totally to reproducethe expected positive torque peak (spin up) at perihelion. https://ssd.jpl.nasa.gov/?horizons R a n g e ( A U ) −400 −300 −200 −100 0 100 200 300 400 Days f om Pe ihelion R e s i d u a l s ( k m ) No NGAYoemans & Chodas
Fig. 5.
Observed Earth-comet range, R , and residuals for twomodels: a purely gravitational N-body solution with no NGA,and a ground-based solution based on the model of Marsdenet al. (1973) (see text for details). Differences between the ob-served and computed solutions, as well as the jumps describedin Sect 2.2, are invisible at the scale of the top plot. −300 −200 −100 0 100 200 Days from Perihelion R e s i d a l s ( k m ) Comp ted Gro nd-based
Fig. 6.
Observed minus computed range, R , for solution A (bluecurve). For comparison, the residuals to the ground-based solu-tion, using the Marsden et al. (1973) model as in Fig. 5, areshown in orange. Both curves diverge from zero most sharplyfollowing the maximum perturbation around perihelion, but oursolution is an improvement. This is confirmed by the RMS (normalised) residuals of 6.59and 5.13, respectively.Figure 9 shows the active fractions for this solution,mapped onto the shape model. Some artefacts of the regiondefinition can be seen, introducing spurious active fractionsat the borders between regions, but these facets representa small fraction of the total area and should not influencethe general results. A large difference between the effectiveactive fractions in the northern and southern hemispherescan be seen in Fig. 9, with the southern hemisphere show-ing active fractions of up to , while the north is onlya few percent active. This is supported by previous works(Marschall et al. 2016, 2017) based on the interpretation ofROSINA data. It is also consistent with the higher active
Article number, page 6 of 12. Attree et al.: Activity models of comet 67P/Churyumov–Gerasimenko constrained by Rosetta
Solution η λ R λ Q λ τ RMS R (km) RMS Q RMS τ RMS
Obj
No NGA 1029 1029Ground-based 675 675A 0.7 1 50 50 163 6.6 5.1 259B 0.7 1 50 50 187 5.2 2.1 235
C 0.7 1 50 50 46 4.2 0.57 125
Ca 0.5 1 50 50 115 5.5 1.32 176Cb 0.6 1 50 50 64 4.3 0.61 130Cc 0.8 1 50 50 62 5.5 1.52 168
Table 1.
Parameters and root-mean-squared residuals (in range, water production rate, non-gravitational torque, and the totalobjective function) for the best-fit models. ‘No NGA’ is an N-body only model computed in REBOUND, while ‘Ground-based’ hasan additional force given by the Marsden et al. (1973) model with parameters from NASA Horizons. A, B and C use the thermalmodel described here. The best model (C) is highlighted in bold font. O u t g a ss i n g R a t e ( s − ) ComputedCo ected ROSINA data−300 −200 −100 0 100 200
Days f om Pe ihelion R e s i d u a l s Fig. 7.
Observed and computed water production rates andresiduals for solution A. T o r q u e ( N m ) Days f om Pe ihelion R e s i d u a l s Fig. 8.
Observed and computed torques and residuals for solu-tion A. fraction (up to a factor of 2–3 in the southern hemispherecompared to the northern hemisphere) found by (Krameret al. 2019) from a thorough study of the evolution ofthe rotational parameters of 67P. The southern regions ofcomet 67P receive higher insolation during southern sum- +Z -Z
Fig. 9.
Mapped active fraction for solution A. mer which occurs near perihelion: the summer solstice hap-pens on Aug 15, 2015, only a couple of days after perihelion.An analysis of OSIRIS images of the Anhur/Bes south-ern regions (Fornasier et al. 2017) shows a high activityoriginating from these regions, combined with high relativeerosion rates deduced from the relatively higher numberof boulders found (Pajola et al. 2016). An examination ofthe production and force curves produced by the individualsuper regions showed that the southern hemisphere domi-nates production after equinox and around perihelion, asexpected. Taken as a whole, the southern hemisphere hasa negative torque efficiency (see Fig. 3), which leads to thedifficulty in jointly fitting both the production and torquecurves, seen in our solution A residuals. Splitting the south-ern hemisphere into more regions is therefore a promisingnext step.
Figures 10 –12 show the residuals to our optimisation withthe full 26 comet regions defined by Thomas et al. (2015),which we refer to as our solution B.Figure 12 shows a clear improvement in how well wefit the torque peak, with the model now showing a posi-tive peak of roughly the right magnitude, although still notmatching the shape. The improved RMS residual of 2.09backs this up. Conversely, however, the range residuals arenow slightly increased, relative to solution A, and the waterproduction curve is still not well fitted; peak production in
Article number, page 7 of 12 &A proofs: manuscript no. preprint −300 −200 −100 0 100 200
Days from Perihelion R e s i d a l s ( k m ) Comp ted Gro nd-based
Fig. 10.
Observed minus computed range for solution B (bluecurve), and the ground-based solution. O u t g a ss i n g R a t e ( s − ) ComputedCo ected ROSINA data−300 −200 −100 0 100 200
Days f om Pe ihelion R e s i d u a l s Fig. 11.
Observed and computed water production rates andresiduals for solution B. T o r q u e ( N m ) Days from Perihelion −2.50.02.55.0 R e s i d u a l s Fig. 12.
Observed and computed torques and residuals for so-lution B. +Z -Z
Fig. 13.
Mapped active fraction for solution B. the model is still too low, and does not fall off fast enoughwith heliocentric distance post-perihelion.The active fraction map of Fig. 13 shows the same trendfor high activity in the southern super-regions as beforebut now with slightly higher activity in Hapi, matchingMarschall et al. (2016). Active fraction in the southern re-gions can be seen to vary significantly and it is instructive tocompare this distribution to the map of torque efficiency.As can be seen from Fig. 3, positive torque efficiency isclustered in several small areas, and these are given highactivity by our optimisation. However, the optimisation islimited by the fact that the geographically defined regionscontaining these areas also contain areas of negative torqueefficiency. In other words, torque efficiency varies at a localscale and is not necessarily correlated with the regions usedin this parametrisation.To address this, we could further subdivide our regions,introducing more parameters, but no obvious best way to dothis presents itself. Instead we take the, somewhat simpli-fied, approach of reverting back to the 5 super region modelof solution A, and simply splitting the southern super re-gion by torque efficiency. This creates two non-contiguousand “spotty” super regions, which do not necessarily corre-spond to particular morphological or structural regions onthe cometary surface, but do provide a parametrised way ofoptimising the NGA model. Experiments with this methodshow significant improvements over models A and B, butstill have problems reproducing the production rate curve,which we address in the next subsection.
Since the above solutions with constant active fractions failto adequately reproduce all the data, we now consider time-varying solutions. Temporal variations of the effective activefraction is an obvious way to try to reconcile the modelledpost-perihelion slope of the water production rate with themeasured data.We begin with the 6 super regions (including a south-ern hemisphere split by torque efficiency) and implementa time-varying active fraction for both the southern hemi-sphere regions (since these are the most important aroundperihelion) while keeping the others constant. We first con-sidered a decline of the active fraction with time, with ahalf-Gaussian fall-off from the initial value, fitting for boththe active fractions and start and decay times of the Gaus-
Article number, page 8 of 12. Attree et al.: Activity models of comet 67P/Churyumov–Gerasimenko constrained by Rosetta −300 −200 −100 0 100 200
Days from Perihelion R e s i d a l s ( k m ) Comp ted Gro nd-based
Fig. 14.
Observed minus computed range for solution C (bluecurve), and the ground-based solution. sian. While this produced encouraging results, it is a some-what unphysical situation: the comet’s orbit is cyclical sothat the active fraction must “reset” back to the initial valuein time for the next perihelion. This may happen slowlyaround aphelion, in which case it will not affect the fit here,or it may occur during the time-period studied by Rosetta.To explore this latter case, we perform a final optimisation,allowing two active fractions for each of the two southernhemisphere regions as well as two start times, t and t ,and two decay (or growth) half-times, t and t , for atotal of twelve fitted parameters. The relative magnitudesof the two active fractions are unconstrained, but the twotimes, t and t , are constrained to lie either side of perihe-lion to ensure that they do not cross over. The half-timesare constrained to be larger than zero but unconstrained atthe top end.Finally, we also vary the momentum transfer efficiency( η parameter) around our nominal value η = 0 . . A full op-timisation of the twelve parameters is performed adopting η values of . , . and . (models Ca, Cb and Cc in Table 1).As shown in Table 1, the smallest multi-objective function(Eq. (3.3)) is achieved for η = 0 . , with a minimum valueof . While the η = 0 . and η = 0 . solutions correspondto a significantly higher multi-objective function (equal to and , respectively), the η = 0 . solution (equal to ) is only marginally larger.Figures 14 –16 show the residuals for our best fit modelin this case: model C, while Figs. 17 and 18 show themapped (peak) active fraction distribution and how itvaries with time. High activity is again favoured in thesouthern regions, with the optimisation now selecting veryhigh effective active fractions, of over , in order to fit thehigh peak production rates at perihelion (Fig. 15). Higheractive fractions in areas producing a positive torque allowthem to dominate the rest of the southern hemisphere, pro-ducing a net positive torque curve, which now matches verywell the observations (Fig. 16). Note however that the smalldrop off observed 125 days before perihelion is not repro-duced by the model.Figure 18 shows that the preferred solution has south-ern active fractions increase quickly between equinox andperihelion (with a half-time of ∼ days) to their high O u t g a ss i n g R a t e ( s − ) ComputedCo ected ROSINA data−300 −200 −100 0 100 200
Days f om Pe ihelion R e s i d u a l s Fig. 15.
Observed and computed water production rates andresiduals for solution C. T o r q u e ( N m ) Days f om Pe ihelion −101 R e s i d u a l s Fig. 16.
Observed and computed torques and residuals for so-lution C. peak values, before falling off to the ∼ level after peri-helion and during southern summer, with a decay half-timeof around 50 days. This reproduces the high productionrate at perihelion while matching the swift fall-off over the +Z -Z Fig. 17.
Mapped final active fraction for solution C.Article number, page 9 of 12 &A proofs: manuscript no. preprint −300 −200 −100 0 100 200
Days from Perihelion A c t i e F r a c t i o n South -South +Region1Region2HathorHapi
Fig. 18.
Active fraction with time for solution C. succeeding two hundred days. Some discrepancies with thedata still remain, for example the “knee” feature seen as theproduction ramps up around a hundred days before perihe-lion, but the result is now a much better match to the dataoverall.The greatly reduced range residual, of 46 km, is furtherevidence of the improved model, and confirms that the ma-jority of the NGA effect is concentrated near perihelion. Thefact that non-gravitational forces and production are bothseen to be strongly peaked near perihelion, over and abovewhat one would expect from geometric considerations, isconsistent with a time-varying active fraction. Small dif-ferences seen in the torque and production curves occurbefore equinox, when outgassing is controlled by the activefraction in the Northern hemisphere, which we do not varywith time. This suggests that minor improvements might bemade to the fit by focusing on the North, although thesewould be unlikely to affect the trajectory and peak produc-tion, both of which are dominated by activity at perihelion(i.e. in the South).The best-fit results allow us to calculate integratedquantitative parameters resulting from the cometary ac-tivity around perihelion. The total water ice mass lossamounts to . kg , corresponding to a mean erosionof ∼ over the entire nucleus surface, assuming a den-sity of
538 kg m − (Preusker et al. 2017). We emphasizethat this estimate does not take into account the dust massloss - predicted to be much larger depending on the dust-to-ice ratio - and the outgassing of more volatile minor speciesthroughout the orbit. The actual water ice erosion can becalculated by integrating the time-varying sublimation rateof each facet. We find a local erosion of . − . in thetwo southern super-regions and < . in the northernones.
5. Discussion
The most striking features of our best solution (C) are thedichotomy of the effective active fraction between the south-ern and northern hemispheres on the one hand, and thedrastic rise of the effective active fraction around perihe-lion in the southern regions on the other hand. The latter is required in our approach to explain the steep slope of theproduction rate. We propose the following qualitative ex-planation, based on the seasonal formation and disappear-ance of a dust mantle, for this cyclic increase of effectiveactive fraction in the southern regions around perihelion.This idea has been introduced a long time ago in the liter-ature. Among the pioneering works, Brin & Mendis (1979),Brin (1980) and Fanale & Salvail (1984) introduced the ideathat a dust mantle can form and be subsequently disruptedby the gas pressure if it remains thin enough. Using a one-dimensional thermal model, Rickman et al. (1990) showedhow the obliquity of the spin axis can influence the stabil-ity of the mantle. They show that a temporary mantle canform at intermediate and polar latitudes for nuclei withradii equal to and high obliquities (equal to ◦ intheir simulations). For nuclei with smaller radii (equal to ), temporary mantles only appear at perihelion dis-tances smaller than . u . . In a more recent work based ona 3D thermal model, De Sanctis et al. (2010) tried to repro-duce the thermal evolution of 67P. The role of the obliquityis emphasized as being critical, high obliquities favouringthe appearance of a stable dust mantle at equatorial lati-tudes. Other models (e.g., Kossacki & Szutowicz 2006) findon the contrary that a stable mantle with a non-uniformdust layer can explain the observed water production rates.In our explanation, the southern regions - includingthe southern polar cap - become progressively illuminatedand the activity starts to increase after the spring equinox( days before perihelion – 11 May 2015). This rise of ac-tivity allows dust present at the surface to be lifted off,decreasing the depth at which water ice is present belowthe dust layer. This is a runaway process as the increasedgas flux resulting from this reduced dust depth is able tolift off larger and larger dust particles from the surface. Anincreasing fraction of these particles eventually reaches ve-locities large enough to escape the nucleus gravity, or tobe redeposited in other nucleus (Northern) regions. Thismechanism leads to an increase of the effective active frac-tion. After the summer equinox ( days after perihelion –15 August 2015), the energy input starts to decrease, re-sulting in a reduced gas flux and surface temperature. Thelarge dust particles can no longer be lifted off from the sur-face and start to be redeposited locally. The apparition ofa dust mantle quenches the cometary activity as the wateris no longer at the surface of the comet; instead, it subli-mates through a deeper dust layer beneath the surface. Thegas diffusion through this dust mantle produces lower gasfluxes, decreasing the effective active fraction of the surface.The process continues for several months until the autumnequinox ( days after perihelion – 24 March 2016) trig-gers the southern autumn, followed by the long southernwinter around aphelion. At that time, the southern regionsare covered again by a dust layer that will be partiallyremoved at the next perihelion passage. It is not totallyexcluded also that the outgassing of more volatile speciesduring the northern summer at aphelion contributes to there-accumulation of material in the southern hemisphere.This scenario is partially supported by observationsfrom instruments aboard Rosetta. Lai et al. (2016) mod-elled the dynamics of dust grains around 67P taking intoaccount OSIRIS observations. They predict that dust par-ticles ejected from the southern hemisphere are redepositedin the northern hemisphere during the southern summer. Arecent geomorphological analysis of the surface (Birch et al. Article number, page 10 of 12. Attree et al.: Activity models of comet 67P/Churyumov–Gerasimenko constrained by Rosetta . around perihelion and . otherwise, agree well with pre-vious estimates for 67P (Lamy et al. 2007) and other comets(see, e.g. A’Hearn et al. 1995), but the large differences be-tween hemispheres and seasons highlights the limitations ofinterpreting cometary activity with a single number fromthe ground. η parameter The η parameter, also called “momentum transfer effi-ciency” in the literature represents the fraction of theMaxwellian thermal velocity of the gas, calculated fromthe nucleus surface temperature, contributing to the mo-mentum transfer. This parameter depends on the local icecontent and on the detailed structure of the porous mate-rial. Its value couples the water production curve with theeffect of the spin/orbit variations, allowing to compensateany systematic effect possibly present in the input data or inthe model. The value of the η parameter has been calculatedfrom gas-kinetic models of the Knudsen layer. Delsemme &Miller (1971) adopted a value η = 0 . , between the pureice plane surface value η = 1 / and higher values up to / predicted for porous media. In early interpretations ofthe NGA of comet 1P/Halley, Rickman (1986) and Sagdeevet al. (1988) used a value of . and . respectively basedon different assumptions. Crifo (1987) recommended a value η ≈ . based on revised gas-kinetic theoretical descriptionof the solid-gas interface, taking into account the recon-densation of water ice. Rickman (1989) used a correctedvalue η = 0 . − . based on the work of Crifo (1987),which seem to be in good agreement with gas velocity mea-surements in the coma of comet 1P/Halley (see Rickman1989, and references therein). However, the calculations donot consider intimate ice-dust mixtures and do not takeinto account the porosity of the surface. Skorov & Rickman(1999) introduced a correction factor of . to the valuesadopted by Rickman (1989) based on an analytical modelof the Knudsen layer above a porous dust mantle. This leadsto very high η values in the range − . .From our analysis of the data presented in this paper,our best fits are obtained for η in the range . − . , in goodagreement with the moderate values adopted in the liter-ature. They do not support more extreme values (around . and greater than . ), which degrade the overall fit of the three data sets. We stress however that we set the sur-face gas temperature to the dust temperature T dust in ourmodel (see Eq. (6)). This may lead to overestimating thegas temperature, which would in turn tend to underesti-mate the fitted value of η . Note, however, that the depen-dence of the gas velocity with the temperature is a squareroot in Eq. (6): a large and constant temperature devia-tion would be needed to significantly bias the value of η .The second point of concern is a possible overestimate ofthe water production curve due to sublimating icy grains inthe coma. This would once again require a larger value of η to compensate the smaller local sublimation rates. Alto-gether, even if our simulations point to η < . we cannotentirely rule out at this point higher values of this param-eter. Consolidated values of the surface water productionaround perihelion, as well as estimates of the gas velocityabove the surface, would help to reduce the uncertainties.
6. Conclusions and perspectives
From our work, the following conclusions could be reached:1. We succeeded in finding an activity pattern explainingsimultaneously the three following 67P data sets, ex-tracted mainly from the Rosetta mission: (1) the Earth-comet ranging data reconstructed by the flight dynam-ics team of ESA/ESOC, (2) the water production ratededuced from a mix of ROSINA and ground-based ob-servations (Hansen et al. 2016), and (3) the rate of spinperiod change deduced from the OSIRIS images. Theresiduals of the ranging data describing the effect of thenon-gravitational acceleration are reduced by an orderof magnitude compared to the ground-based solutionbased on the simple model of Marsden et al. (1973).2. The fitted activity pattern exhibits two main features:a higher effective active fraction in two southern super-regions ( ∼ %) outside perihelion compared to thenorthern ones ( < %), and a drastic rise of the effectiveactive fraction of the southern regions ( ∼ − %)around perihelion.3. In order to successfully fit the positive rate of spin pe-riod change, we need to split the southern super regioninto two entities, depending on the sign of the torque ef-ficiency. These two entities correspond to two relativelywell delimited areas in Anhur, Bes and Khepry, but cre-ates a patchy separation in Wosret.4. We interpret the time-varying southern effective activefractions by cyclic formation and removal of a dust man-tle in these regions. According to our interpretation, thedust mantle could be progressively removed when activ-ity rises after the southern spring equinox and formedagain when activity decreases towards the southern au-tumn equinox (and possibly around aphelion during thenorthern summer). Several observations performed dur-ing the Rosetta mission, such as dust transport fromSouth to North (Lai et al. 2016; Birch et al. 2017) andbluer colours observed near perihelion (Fornasier et al.2016), support this interpretation.5. If it is confirmed, this interpretation would stronglysupport post-Halley thermal modelling (Rickman et al.1990; De Sanctis et al. 2010) which predicted that sea-sonal effects linked to the orientation of the spin axisplay a key role in the formation and evolution of dustmantles, and in turn largely control the temporal vari-ations of the gas flux. Article number, page 11 of 12 &A proofs: manuscript no. preprint
6. Our analysis supports moderate values of the momen-tum transfer coefficient η in the range . − . . For moreextreme values of this coefficient ( ≤ . and ≥ . ), thefit of the three data sets is degraded. However, we willnot be able to rule out higher values of this parameterwithout consolidated water production measurementsand, to a lesser extent, without estimates of the near-surface thermal gas velocity.More work will be needed to better understand the ac-tivity of 67P using data collected during the Rosetta mis-sion. The solution found in this article through the pro-cess described in section 4 is non-unique on the one hand,and could probably be improved on the other hand. Im-provements may come from a better understanding of theranging data, resulting in the extraction of clean and accu-rate non-gravitational acceleration of the comet as a func-tion of time around perihelion. The change in the directionof the spin axis or angular velocity could also provide avaluable additional data set that could help to constrainthe activity pattern and reduce the number of solutions. Intheir interpretation of the temporal evolution of the rota-tional parameters, Kramer et al. (2019) did not introducea temporal variation of the activity around perihelion butconsidered instead a spatially heterogeneous surface with36 “patches” having different water-ice coverage. They alsodiscuss the possibility of a decreasing dust layer near perihe-lion increasing activity, and found a solution explaining thewater production curve of 67P. It would be interesting totest if this solution could also reproduce the NGA of comet67P described in this article. Finally, a better understand-ing of the production rate around perihelion, reconcilingthe different Rosetta instrument measurements and includ-ing species other than water, would be of benefit in fittingthe activity model. Acknowledgements.
This project has received funding from the Euro-pean Union’s Horizon 2020 research and innovation programme undergrant agreement no. 686709. This work was supported by the SwissState Secretariat for Education, Research and Innovation (SERI) un-der contract number 16.0008-2. The opinions expressed and argumentsemployed herein do not necessarily reflect the official view of theSwiss Government. We thank the authors of several python pack-age libraries which were used extensively here; these include
NumPy,SciPy, SpiceyPy, REBOUND and
REBOUNDx . We are grateful toL. Maquet for fruitful discussions in the course of this work. We alsothank the reviewer, Dennis Bodewits, for his comprehensive and usefulreview.
References