Issue 
A&A
Volume 630, October 2019
Rosetta mission full comet phase results



Article Number  A18  
Number of page(s)  11  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834415  
Published online  20 September 2019 
Constraining models of activity on comet 67P/ChuryumovGerasimenko with Rosetta trajectory, rotation, and water production measurements
^{1}
AixMarseille Université, CNRS, CNES, Laboratoire d’Astrophysique de Marseille, Marseille, France
^{2}
Earth and Planetary Observation Centre, Faculty of Natural Sciences, University of Stirling,
UK
email: n.o.attree@stir.ac.uk
^{3}
Physikalisches Institut, Universität Bern,
Sidlerstrasse 5, 3012 Berne, Switzerland
^{4}
Deutsches Zentrum für Luftund Raumfahrt (DLR), Institut für Planetenforschung,
Rutherfordstraße 2,
12489 Berlin, Germany
^{5}
MaxPlanckInstitut für Sonnensystemforschung,
JustusvonLiebigWeg 3,
37077 Göttingen, Germany
^{6}
International Space Science Institute,
Hallerstrasse 6,
3012 Bern, Switzerland
Received:
11
October
2018
Accepted:
8
January
2019
Aims. We use four observational data sets, mainly from the Rosetta mission, to constrain the activity pattern of the nucleus of comet 67P/ChuryumovGerasimenko (67P).
Methods. We developed a numerical model that computes the production rate and nongravitational acceleration of the nucleus of comet 67P as a function of time, taking into account its complex shape with a shape model reconstructed from OSIRIS imagery. We used 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 of 67P. The two key parameters of our model, adjusted to fit the three data sets all together, are the activity pattern and the momentum transfer efficiency (i.e., the socalled η parameter of the nongravitational forces).
Results. We find an activity pattern that can successfully reproduce the three data sets simultaneously. The fitted activity pattern exhibits two main features: a higher effective active fraction in two southern superregions (~10%) outside perihelion compared to the northern regions (<4%), and a drastic rise in effective active fraction of the southern regions (~25−35%) around perihelion. We interpret the timevarying southern effective active fraction by cyclic formation and removal of a dust mantle in these regions. Our analysis supports moderate values of the momentum transfer coefficient η in the range 0.6–0.7; values η ≤ 0.5 or η ≥ 0.8 significantly degrade 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 the formation and evolution of dust mantles, and in turn, they largely control the temporal variations of the gas flux.
Key words: planets and satellites: dynamical evolution and stability / comets: individual: 67P/ChuryumovGerasimenko / comets: general
© N. Attree et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The sublimation of ices when a comet is injected from its reservoir into the inner solar system triggers the emission of molecules. This outgassing in turn produces a reaction force that can accelerate the comet nucleus in the opposite direction. The perturbing effect of cometary activity on the trajectory of comets has been established in the 1950s in the pioneering work by Whipple (1950). At that time, it was clear that most comet trajectories were affected by a significant nongravitational acceleration (hereafter NGA) linked to their activity around perihelion (Marsden 1968). Shortly after, a theoretical model describing the nongravitational force (hereafter NGF) produced by the sublimation of ices was established by Marsden et al. (1973). This model was based on a simple function describing the heliocentric dependence of the sublimation of cometary water ice, combined with constant scaling parameters A_{1}, A_{2}, and A_{3} describing the amplitude of the NGA along the three components (radial, transverse, and normal) in the orbital frame of the comet. The model has been modified by Yeomans & Chodas (1989) to incorporate an asymmetric term used to describe the shift of the maximum of activity with respect to perihelion. It is worth mentioning here that these simple models are still in use today to fit astrometric measurements, and hence to describe cometary orbits.
More sophisticated models have been introduced since then. Images acquired during the flyby of comet 1P/Halley by the Giotto spacecraft showed narrow dust jets (Keller et al. 1987), leading to the idea that the activity may be confined to localized areas. This led to a new model of nongravitational acceleration (Sekanina 1993) in which the outgassing originates from discrete areas at the surface of a rotating nucleus. Szutowicz (2000) used Sekanina’s model to fit astrometric measurements of comet 43P/WolfHarrington that were obtained during nine perihelion passages. She also compared the modeled production rate with visual light curves that were used as a proxy for the activity of the comet (Szutowicz 2000). Recently, Maquet et al. (2012) revisited Sekanina’s approach with a model in which the activity is parameterized by the surface fraction of exposed water ice in “latitudinal bands” at the surface of an ellipsoidal nucleus. More accurate groundbased measurements as well as space missions to comets offered the opportunity to incorporate new physical processes in models of the NGA. Rickman (1989) used the change in orbital period of comet 1P/Halley caused by the NGA to extract its mass and density. Thedetailed description of the local outflow velocity incorporated a “local momentum transfer coefficient” ζ_{l}, originally called η in the improved description of the solid–gas interface introduced by Crifo (1987). This coefficient represents the fraction of the emitted gas momentum (dependent on its thermal velocity), which needs to be considered in the calculation of the momentum transfer, and thus of the NGF (see Crifo 1987). Davidsson & Gutiérrez (2004) used 2D thermal modeling including thermal inertia, selfshadowing, selfheating, and an activity pattern to fit the NGA of comet 19P/Borrelly. They were able to retrieve the mass of the nucleus and constrain the direction of the spin axis. The method was also applied to comets 67P/ChuryumovGerasimenko (Davidsson & Gutiérrez 2005), 81P/Wild 2(Davidsson & Gutiérrez 2006), and 9P/Tempel 1 (Davidsson et al. 2007), all targets of space missions.
The exploration of comet 1P/Halley in 1986 also yielded the discovery of the nonprincipal axis rotation of this comet (Samarasinha & A’Hearn 1991). Since then, the torque of the NGF was thus identified as the main effect responsible for changes in rotationalparameters (Samarasinha et al. 2004, and references therein). Modeling the torque is required to understand the observed change in the rotational parameters of cometary nuclei, as well as the apparition of rotations of the nonprincipal axis. Changes in the spin period of several comets have been detected from the analysis of light curves (Mueller & Ferrin 1996; Gutiérrez et al. 2003; Drahus & Waniak 2006; Knight et al. 2011; Bodewits et al. 2018). Predictions of the expected change in the direction of the spin axis and spin period for comets 81P/Wild 2 (Gutiérrez & Davidsson 2007) and 67P/ChuryumovGerasimenko (hereafter 67P, Gutiérrez et al. 2005) have been made. Recently, a change in the spin period of comet 67P has clearly been identified between the 2009 and 2015 perihelion passages (Mottola et al. 2014) based on early images acquired by the OSIRIS camera on board the Rosetta spacecraft. Keller et al. (2015) showed that the spin period variation curve of 67P is controlled at first order by the bilobate shape of the nucleus.
Through its long journey accompanying comet 67P, Rosetta provides a unique chance to record measurements of most parameters involved in the modeling of the NGF. The mass of the comet has been retrieved by Pätzold et al. (2016) from the radio science experiment on board the spacecraft. The shape has been retrieved from a stereophotogrammetric analysis of a subset of OSIRIS images (Preusker et al. 2017), leading to an accurate knowledge of the moments of inertia. The activity pattern has been constrained in the early phase of the mission by Marschall et al. (2016, 2017) from ROSINA measurements in the form of “effective active fractions” associated with geological regions. The total water production rate has been constrained from ROSINA measurements, complemented by groundbased observations (Hansen et al. 2016). Finally, the spin period has been monitored throughout the mission by ESA’s flight dynamics and OSIRIS teams before (Jorda et al. 2016) and after (Kramer et al. 2019) perihelion. The latter (Kramer et al. 2019) modeled the temporal evolution of the rotational parameters, comparing it with the measured change in spin period and spin axis orientation. The aim of this article is to try to reproduce three data sets derived by several Rosetta instruments (see Sect. 2) with a model of the NGF (see Sect. 3) in order to retrieve (i) the local effective active fraction and its temporal variations around perihelion, and (ii) a recommended value for the momentum transfer coefficient η (see Sect. 4). The results are discussed in Sect. 5, together with recommendations for the calculation of NGF of other comets.
2 Observational constraints
In this section, we describe the observational data, mainly obtained by the Rosetta spacecraft. We attempt to fit our NGF model to these data.
2.1 Water production rate
The total water production rate is an important constraint for any model of cometary activity, and a significant effort has been made to measure it for 67P. As summarized by Hansen et al. (2016), a number of different instruments have all been used, and comparing and synthesizing their results in nontrivial. Here, we used the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) measurements, empirically corrected for spacecraft position, as our observed data points, O_{Q}. Hansen et al. (2016) described that estimating the uncertainty in these data, σ_{Q}, is difficult, therefore we used the bounds on the power law, fitted with heliocentric distance to the inbound ROSINA data, as given in their Table 2.
We point out that ROSINA data are inferred from local measurements in the coma of 67P. Around perihelion, the spacecraft was located at a distance of 200 km from the nucleus, making it difficult to infer a total production rate, while groundbased observations (see, e.g., Bertaux 2015) have also suggested a variation in peak production between perihelion passages. Production rate estimates from the Rosetta lineofsight instruments MIRO and VIRTIS are also generally lower that the results of Hansen et al. (2016). These are important caveats to bear in mind when interpreting our results. Finally, the possibility that sublimating icy grains are emitted by the nucleus is not considered in this article, and neither are other gas species, such as CO_{2}, CO, and O_{2}. Other species represent < 10% of the gas number density at perihelion (Hansen et al. 2016), and their production curves are not as well constrained as water. Therefore, for this study we chose to focus on water, which is the primary driver of nongravitational forces.
2.2 Trajectory
Outgassing produces a backreaction force and a resulting NGA on cometary nuclei that affects their trajectory in a measurable way. For 67P, the nucleus trajectory has been reconstructed by the flight dynamics team of ESA by a combination of radiotracking of the Rosetta spacecraft from Earth and optical navigation of the comet relative to it, and it is available in the form of NASA SPICE kernels (Acton 1996). The 3D position of the comet in a heliocentric reference frame therefore has the greatest accuracy in the Earthcomet range direction (R, claimed accuracy of ~10 m) and a much lower accuracy in the perpendicular (crosstrack) directions (claimed accuracy of ~ 100 km).
Theoretically, the NGA resulting from outgassing could be directly extracted from the residuals between this measured trajectory and a modeled gravitational orbit (taking into account general relativity and the gravitational accelerations of all major planets, Pluto and the most massive asteroids). During the course of this work, however, it was discovered (and later confirmed by ESAC; B. Grieger, priv. comm.) that the reconstructed comet and spacecraft trajectories contain a series of discontinuities, at which the objects’ positions vary over hundreds of meters to several kilometers in an instantaneous time, making the above method impossible. Within the orbital segments between these “jumps”, there is no difference, to machine precision, between our own orbital integrations (see Sect. 3.2 below) and the reconstructed trajectory, demonstrating that it is a purely gravitational solution. The jumps occur at the boundaries between the integration segments that make up the reconstructed orbit, and represent the offset in the objects’ positions over the course of each segment due to the unmodeled nongravitational acceleration. Unfortunately, it proved impossible to extract the NGAs directly from the jumps themselves as they contain not only the NGA effect, but also the typical uncertainty in the state vector at the start of each segment, which is of a comparable magnitude. The jumps therefore have random magnitudes with time (Fig. 1) and must be considered an additional source of noise in the uncertainty in the measured positions.
Despite these issues, the reconstructed kernels remain a good description of the cometary trajectory over orbital timescales and within the limits of accuracy of the typical jump size. It is therefore still possible to compare these measurements with a model of the orbit, including a thermal outgassing NGA as well as Nbody gravitational interactions, in order to constrain said model. To do so, we used the magnitude of the comettoEarthcenter range as our observable O_{R}, since the jump sizes are smallest in this direction (Fig. 1). Considering the jumps to be a source of random error, we conservatively estimate the uncertainty in O_{R} as σ_{R} ≈ 1 km.
Fig. 1 Discontinuities identified in the position of comet 67P from the SPICE kernels, in (x, y, z) heliocentric J2000 coordinates, and Earthcomet range (R). Left panel: as a function of time. Right panel: as a histogram of sizes. 
2.3 Rotation
Backreaction from outgassing not only produces a net acceleration on the nucleus, but can also, depending upon its shape, produce a net torque, altering its rotation state. The rotational parameters of 67P, including its spin rate over time, ω(t), has been measured as part of the reconstruction of its 3D shape from OSIRIS images (Jorda et al. 2016).
After verifying that the cross terms relating to the angular velocities along the first and second principal axes of the comet are negligible, changes in spin rate, , can be directly related to the z component of the torque (τ_{z}, in a bodyfixed frame) by Eq. (1), where I_{z} = 1.899 × 10^{19} kg m^{2} is the third (largest) moment of inertia derived from the shape model assuming a constant density of 538 kg m^{−3} (Preusker et al. 2017), (1)
Differentiating the observed ω(t) by time exacerbates measurement uncertainties so that the produced torque curve becomes extremely noisy. To compare this with our simulations below, we therefore smoothed the data by fitting them to a cubic spline, as shown in Fig. 2. Our fitted spline was then used as the torque observable, O_{τ}, with an assumed uncertainty equal to the root mean squared residuals of the derived minus smoothed data, σ_{τ} = 575 000 N.m (shown in gray bounds in Fig. 2).
Fig. 2 Observed torque, derived from the 67P rotation state from OSIRIS measurements, and a smoothed cubic spline fit to the data. The gray region represents the RMS of the residuals between the two. 
3 Modeling
3.1 Thermal model
A thermal model is required to compute the temperature of the sublimating layer (assumed to be at the surface on the nucleus), from which we can derive the nongravitational forces acting on the nucleus. Our thermal model takes into account solar insolation, surface thermal emission, sublimation of water ice, projected shadows, and selfheating. We used a decimated version of the 67P shape model called SHAP7 (Preusker et al. 2017), with 124 938 facets. The temperature was computed for each facet of the shape model, 36 times per rotation (i.e., every 1240 s), and 70 times (i.e., every ~ 10 days) over the two 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 relative to the Sun, and the projected shadows were computed using the OASIS software (Jorda et al. 2010). Because of the large number of facets (>100 000) and time steps (>2500), heat conduction is neglected in the thermal model for numerical reasons. To test this assumption, we computed the production rate (Eq. (4)) and acceleration (Eq. (5)) of a spherical nucleus at perihelion (where the torque is maximum) for two cases: a thermal inertia of 0 and 100 J m^{−2} K^{−1} s^{−0.5}. The production rate and the acceleration are ~7% smaller for a thermal inertia of 100 J m^{−2} K^{−1} s^{−0.5} (compared to a null thermal inertia), and the direction of the acceleration vector differs by less than 3°. Neglecting the thermal inertia therefore appears reasonable compared to the data uncertainties (e.g., production rates).
The surface energy balance of the thermal model is given by Eq. (2) for a facet with index i, where A_{b} = 0.0119 is the Bond albedo at 480 nm (Fornasier et al. 2015), F_{sun} = 1370 W m^{−2} is the solar constant, z_{i} is the zenithal angle, r_{h} the heliocentric distance, SH_{i} is the selfheating given by Eq. (3), ϵ = 0.95 is the assumed infrared emissivity, T_{i} is the surface temperature, f_{i} is the fraction of water ice (in our case, either 0 for pure dust or 1 for pure ice, see below), α = 0.25 accounts for the recondensation of water ice on the surface (Crifo 1987), L = 2.66 × 10^{6} J Kg^{−1} is the latent heat of sublimation of water ice at 200 K, and Z_{i} is the water sublimation rate given by Eq. (4), (2)
In Eq. (3), the selfheating SH_{i} is the sum of the infrared flux 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 the normal to facet j and the vector joining facets i and j, θ_{i} the angle between the normal to facet i and the vector joining facets i and j, and is the distance between facets i and j. This formalism is similar to that of Gutiérrez et al. (2001). To determine which facets are seeing each others, we used an algorithm developed at Laboratoire d’Astrophysique de Marseillebased on ray tracing and hierarchical search. The computation of the viewing factors is purely geometric and depends only on the shape model: it was therefore only performed once at the beginning: (3)
In Eq. (4), the sublimation mass flow rate per facet is calculated with the molar mass M = 0.018 kg, the two constants A = 3.56 × 10^{12} Pa and B = 6162 K for water (Fanale & Salvail 1984), and the gas constant R = 8.3144598 J K^{−1} mol^{−1}, (4)
Finally, to compute the nongravitational 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 temperature of dust T_{dust} for v_{i}. We therefore ran our thermal model (Eqs. (2)–(4)) with two extreme cases: (1) with f_{i} = 0, which corresponds to a pure dust model, to compute T_{dust}, and (2) with f_{i} = 1, which corresponds to a pure water ice model, to compute T_{ice}.
3.2 Nongravitational force model
The reaction force vector per facet was then calculated based on this mass flow rate, and the total acceleration is the sum over all facets divided by the comet mass (M_{67P} = 9.982 × 10^{12} kg; Pätzold et al. 2016): (5)
where η is the momentum transfer coefficient (Crifo 1987; Rickman 1989), which we here assumed to be constant across the comet. S_{i} is the surface area of each facet (in the direction of its normal), and x_{i} is its effective active fraction. The mass flow rate is calculated with Eq. (4) and the gas velocity is taken as the thermal velocity (6)
assuming equilibrium with the surface graybody temperature, that is, that of the dust from run (1) of the thermal model. This is the upper limit that the gas can reasonably reach, meaning that our nongravitational force will be at the high end of estimation and our effective active fractions are lower limits. Calculated dust temperatures range between ~ 20 and 390 K (ice temperatures are limited by the sublimation to ~200 K), leading to thermal velocities of ~ 155−658 m s^{−1}.
Water production and torque were likewise summed over the surface (7)
where torque per facet is the vector product of each force vector with its radius vector to the center of mass (r_{i}). This can also be expressed as the magnitude of the force multiplied by a torque efficiency (see Keller et al. 2015): (8)
The z component of the total net torque can then be compared with the observations using Eq. (1). It is advantageous to use the torque efficiency formalism since this vector is in the bodyfixed frame and can be calculated once at the beginning of the simulation run, rather than being recalculated each time. Torque efficiency is also a useful way of visualizing the effects of differing spatial distributions in activity, which is important later during the optimisation. Mapping torque efficiency onto the shape model (Fig. 3) shows how local variations in topography combine with largescale orientations of regions, varying the effects of activity on the cometary rotation across its surface (compare with Fig. 1 in Keller et al. 2015, which uses an older, incomplete shape model that lacks the southern hemisphere).
Nongravitational forces were calculated for each of the 70 × 36 = 2520 runs of the thermal model, and the relevant quantities (force and torque vectors and summed water production) were averaged over a day. This produces smoothly timevarying curves that can be inspected at any time of interest using bilinear interpolation, which we refer to as our model solution. For comparison with the observed water production rate and torque, we evaluated our model at the time of each measurement, producing C_{Q} and C_{τ}, and compare this directly.
For the trajectory, however, a full Nbody integration must be performed and the resulting position compared at each time (C_{R}). To do this, we used the opensource REBOUND code^{1} (Rein & Liu 2012), complete with full general relativistic corrections (Newhall et al. 1983), as implemented by the REBOUNDx extension package^{2}. All the major planets are included, as well as Ceres, Pallas, and Vesta, and they are initialized with their positions in the J2000 ecliptic coordinate system according to the same ephemerides as were used in the Rosetta trajectory reconstructions (NASA/JPL solar system solution DE405; Standish 1998). An additional particle representing 67P was initialized with its position given by SPICE. The system was then integrated forward in time using the IAS15 integrator (Rein & Spiegel 2015) and the standard equations of motion, with the addition of a custom acceleration, a_{NG}, for 67P, provided by our model. The modeled comet begins to diverge from the measured positions and at each time of interest we directly compare the magnitudes of the computed and measured comettoEarthcenter ranges, R (the most accurate part of the trajectory as described in Sect. 2.2 above).
3.3 Optimization
In order to constrain the unknown parameters in our model, such as the surface active fraction, we performed a bounded leastsquares fit to the data using the dogleg optimization routine (Voglis & Lagaris 2018) provided in the scientific Python package. The optimization proceeds, attempting to minimize the standard objective function (9)
with observed minus computed residuals, OC, and observation uncertainties, σ, at each timestep, j, up to the total N. The root mean squared residuals are then .
In our case we have three separate data sets to fit to (R, Q and τ), a multiobjective optimisation problem, and we therefore used a linearly scaled combination of the three to generate a combined objective function. The term inside the brackets in Eq. (9) then becomes the sum of the three terms (10)
respectively, where N = N_{R} + N_{Q} + N_{τ} runs over the combined number of points in all three data sets, and the λ scaling coefficients are variables, which themselves must be optimized in order to give the desired weighting to each data set. We set λ_{R} = 1 and scaled the other lambdas relative to it by bootstrapping from the residuals of a preliminary run, to give roughly equal weighting to all three data sets. Because of the very small relative errors for range (σ_{R} ∕R ~±1 km ∕1 AU), this results in high values for the other lambdas (λ_{Q} ~ λ_{τ} ~ 50) to give equal weighting.
For each optimization, we fixed the coupling factor, η, at a constant value and parameterized the model in terms of effective active fraction, x, across the surface. The effects of η on the bestfit model can then be studied independently. To begin with, we used the division of the 67P surface into five “super regions”, as performed by Marschall et al. (2016, 2017) in fitting ROSINA/COPS and OSIRIS data, and started the optimization with initial values from their results. This divided the comet into a southern hemisphere region, an equatorial region, and a region covering thebase of the body and top of the head, Hathor, and Hapi, as shown in Fig. 4 below. The fitting routine then proceeded to optimize these five free parameters, subject to the lower and upper boundaries of zero and one, that is, 0–100% active fraction. As described below, we also used more detailed parameterizations, including all 26 regions of Thomas et al. (2015) for 26 free parameters.
Fig. 3 Torque efficiency in meters, a factor determined by the geometry as defined in Eq. (8). The cometary rotation axisis in the + z direction through the neck region in the center. 
Fig. 4 Map of the five “super regions” defined by Marschall et al. (2016, 2017) and used in our solution A optimization. 
4 Results
Before optimizing our activity model, we first tested the Nbody component by calculating the comet trajectory over the course of the Rosetta mission with no NGAs applied, and with the classic NGA parameterization based on the model of Marsden et al. (1973). This model computes the three components of a_{NG} (radial, alongorbit, and normaltoorbit) based on a power law with heliocentric distance, and three scaling parameters A_{1,2,3} found by a formal bestfit to the orbit for each comet. We used the A_{1,2,3} values for the 2010 apparition of 67P from groundobservations given by NASA/JPL Horizons ephemerides^{3}.
Figure 5 shows a plot of the observed range plus the residuals to both models. The RMS residuals are 1029 km and 675 km, respectively, showing the order of magnitude of the accuracy of the model of Marsden et al. (1973) with groundbased observations of roughly arcsecond accuracy. This stands as a baseline against which we can check our own activity model.
Fig. 5 Observed Earthcomet range, R, and residuals for two models: a purely gravitational Nbody solution with no NGA, and a groundbased solution based on the model of Marsden et al. (1973; see text for details). Differences between the observed and computed solutions, as well as the jumps described in Sect 2.2, are invisible at the scale of the top plot. 
4.1 Fiveparameter solution (A)
Figures 6 –8 show the residuals to our bestfit solution using the 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 for all results), represent a significant improvement over the ~ 1000 km of the purely gravitational solution and the ~600 km of the groundbased solution, demonstrating the significance of our NGA/Nbody model. However, the water production curve is clearly not a good fit, failing to reproduce both the peak production rate as well as overestimating the production far from perihelion. Likewise, the torque curve is an extremely bad fit, failing totally to reproduce the expected positive torque peak (spin up) at perihelion. This is confirmed by the RMS (normalized) residuals of 6.59 and 5.13, respectively.
Figure 9 shows the active fractions for this solution, mapped onto the shape model. Some artifacts of the region definition can be seen, introducing spurious active fractions at the borders between regions, but these facets represent a small fraction of the total area and should not influence the general results. A large difference between the effective active fractions in the northern and southern hemispheres can be seen in Fig. 9, with the southern hemisphere showing active fractions of up to 12%, while the north is only a few percent active. This is supported by previous works (Marschall et al. 2016, 2017) based on the interpretation of ROSINA data. It is also consistent with the higher active fraction (up to a factor of 2–3 in the southern hemisphere compared to the northern hemisphere) found by Kramer et al. (2019) from a thorough study of the evolution of the rotational parameters of 67P. The southern regions of comet 67P receive higher insolation during southern summer, which occurs near perihelion: the summer solstice occured on August 15, 2015, only some days after perihelion. An analysis of OSIRIS images of the Anhur/Bes southern regions (Fornasier et al. 2017) shows a high activity originating from these regions, combined with high relative erosion rates deduced from the relatively higher number of boulders found (Pajola et al. 2016). An examination of the production and force curves produced by the individual super regions showed that the southern hemisphere dominates production after equinox and around perihelion, as expected. Taken as a whole, the southern hemisphere has a negative torque efficiency (see Fig. 3), which leads to the difficulty in jointly fitting both the production and torque curves, seen in our solutionA residuals. Splitting the southern hemisphere into more regions is therefore a promising next step.
Fig. 6 Observed minus computed range, R, for solution A (blue curve). For comparison, the residuals to the groundbased solution, using the model of Marsden et al. (1973) as in Fig. 5, are shown in orange. Both curves diverge from zero most sharply following the maximum perturbation around perihelion, but our solution is an improvement. 
Fig. 7 Observed and computed water production rates and residuals for solution A. 
Fig. 8 Observed and computed torques and residuals for solution A. 
Parameters and rootmeansquared residuals (in range, water production rate, nongravitational torque, and the total objective function) for the bestfit models.
Fig. 9 Mapped active fraction for solution A. 
4.2 Twentysixparameter solution (B)
Figures 10–12 show the residuals to our optimization with the 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 the quality of the torque peak fit, with the model now showing a positive peak of roughly the correct magnitude, although it still does not match the shape. The improved RMS residual of 2.09 supports this. Conversely, however, the range residuals are now slightly increased relative to solution A, and the water production curve is still not well fit; peak production in the model is still too low, and does not fall off fast enough with heliocentric distance postperihelion.
The active fraction map of Fig. 13 shows the same trend for high activity in the southern super regions as before, but now with slightly higher activity in Hapi, matching Marschall et al. (2016). The active fraction in the southern regions can be seen to vary significantly, and it is instructive to compare this distribution to the map of torque efficiency. Figure 3 shows that positive torque efficiency is clustered in several small areas, and these are given high activity by our optimization. However, the optimization is limited by the fact that the geographically defined regions containing these areas also contain areas of negative torque efficiency. In other words, torque efficiency varies at a local scale and is not necessarily correlated with the regions used in this parameterization.
To address this, we could further subdivide our regions, introducing more parameters, but no obvious best way to do this presents itself. Instead we take the somewhat simplified approach of reverting to the five super region model of solution A, and simply splitting the southern super region by torque efficiency. This creates two noncontiguous and “spotty” super regions, which do not necessarily correspond to particular morphological or structural regions on the cometary surface, but do provide a parameterized way of optimizing the NGA model. Experiments with this method show significant improvements over models A and B, but still have problems in reproducing the production rate curve, which we address in the next subsection.
Fig. 10 Observed minus computed range for solution B (blue curve), and the groundbased solution. 
Fig. 11 Observed and computed water production rates and residuals for solution B. 
Fig. 12 Observed and computed torques and residuals for solution B. 
Fig. 13 Mapped active fraction for solution B. 
4.3 Timevarying solution (C)
Because the above solutions with constant active fractions failed to adequately reproduce all the data, we now consider timevarying solutions. Temporal variations of the effective active fraction is an obvious way to try to reconcile the modeled postperihelion slope of the water production rate with the measured data.
We begin with the six super regions (including a southern hemisphere split by torque efficiency) and implement a timevarying active fraction for both the southern hemisphere regions (because these are the most important around perihelion) while keeping the others constant. We first considered a decline in active fraction with time, with a halfGaussian falloff from the initial value, fitting for both the active fractions and start and decay times of the Gaussian. While this produced encouraging results, it is a somewhat unphysical situation: the cometary orbit is cyclical, so that the active fraction must be “reset” to the initial value in time for the next perihelion. This may occur slowly around 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 performed a final optimization and allowed two active fractions for each of the two southern hemisphere regions as well as two start times, t_{0} and t_{1}, and two decay (or growth) halftimes, and , for a total of 12 fitted parameters. The relative magnitudes of the two active fractions were unconstrained, but the two times t_{0} and t_{1} were constrained to lie either side of perihelion to ensure that they did not cross over. The halftimes were constrained to be larger than zero but unconstrained at the top end.
Finally, we also varied the momentum transfer efficiency (η parameter) around our nominal value η = 0.7. A full optimization of the 12 parameters was performed by adopting η values of 0.5, 0.6, and 0.8 (models Ca, Cb, and Cc in Table 1). As shown in Table 1, the smallest multiobjective function (Eq. (10)) is achieved for η = 0.7, with a minimum value of 124. While the η = 0.5 and η = 0.8 solutions correspond to a significantly higher multiobjective function (equal to 175 and 168, respectively), the η = 0.6 solution (equal to 130) is only marginally larger.
Figures 14–16 show the residuals for our bestfit model in this case: model C, while Figs. 17 and 18 show the mapped (peak) active fraction distribution and how it varies with time. High activity is again favored in the southern regions, with the optimization now selecting very high effective active fractions, of over 35%, in order to fit the high peak production rates at perihelion (Fig. 15). Higher active fractions in areas producing a positive torque allow them to dominate the rest of the southern hemisphere, producing a net positive torque curve, which now matches the observations very well (Fig. 16). We note, however, that the small dropoff observed 125 days before perihelion isnot reproduced by the model.
Figure 18 shows that in the preferred solution, the southern active fractions increase quickly between equinox and perihelion (with a halftime of ~ 25 days) to their high peak values, before falling off to the ~10% level after perihelion and during southern summer, with a decay halftime of around 50 days. This reproduces the high production rate at perihelion while matching the swift falloff over the succeeding 200 days. Some discrepancies with the data still remain, for example, the knee feature seen as the production ramps up around 100 days before perihelion, but the result is now a much better match to the data overall.
The greatly reduced range residual of 46 km is further evidence of the improved model, and confirms that the majority of the NGA effect is concentrated near perihelion. The fact that nongravitational forces and production are both seen to be strongly peaked near perihelion, over and above what would be expected from geometric considerations, is consistent with a timevarying active fraction. Small differences seen in the torque and production curves occur before equinox, when outgassing is controlled by the active fraction in the northern hemisphere, which we did not vary with time. This suggests that minor improvements might be made to the fit by focusing on the north, although these would be unlikely to affect the trajectory and peak production, both of which are dominated by activity at perihelion (i.e., in the south).
The bestfit results allow us to calculate integrated quantitative parameters resulting from the cometary activity around perihelion. The total waterice mass loss amounts to 4.5 × 10^{6} kg, corresponding to a mean erosion of ~ 9 cm over the entire nucleus surface, assuming a density of 538 kg m^{−3} (Preusker et al. 2017). We emphasize that this estimate does not take into account the dust mass loss – which is predicted to be much larger depending on the dusttoice ratio – and the outgassing of more volatile minor species throughout the orbit. The actual waterice erosion can be calculated by integrating the timevarying sublimation rate of each facet. We find a local erosion of 0.4– 1.4 m in the two southern super regions and <0.1 m in the northern ones.
Fig. 14 Observed minus computed range for solution C (blue curve), and the groundbased solution. 
Fig. 15 Observed and computed water production rates and residuals for solution C. 
Fig. 16 Observed and computed torques and residuals for solution C. 
Fig. 17 Mapped final active fraction for solution C. 
Fig. 18 Active fraction with time for solution C. 
5 Discussion
5.1 Temporal variation of the effective active fraction
The most striking features of our best solution (C) are the dichotomy of the effective active fraction between the southern and northern hemispheres on the one hand, and the drastic rise of the effective active fraction around perihelion in the southern regions on the other hand. The latter is required in our approach to explain the steep slope of the production rate. We propose the following qualitative explanation, based on the seasonal formation and disappearance of a dust mantle, for this cyclic increase in effective active fraction in the southern regions around perihelion. This idea has previously been introduced. Among the pioneering works, Brin & Mendis (1979), Brin (1980) and Fanale & Salvail (1984) introduced the idea that a dust mantle can form and be subsequently disrupted by the gas pressure if it remains thin enough. Using a 1D thermal model, Rickman et al. (1990) showed that the obliquity of the spin axis can influence the stability of the mantle. They showed that a temporary mantle can form at intermediate and polar latitudes for nuclei with radii equal to 5 km and high obliquities (equal to 90° in their simulations). For nuclei with smaller radii (equal to 1 km), temporary mantles only appear at perihelion distances smaller than 1 a.u. In a more recent work based on a 3D thermal model, De Sanctis et al. (2010) tried to reproduce the thermal evolution of 67P. The role of the obliquity was emphasized as being critical, high obliquities favoring the appearance of a stable dust mantle at equatorial latitudes. Other models (e.g., Kossacki & Szutowicz 2006) found in contrast that a stable mantle with a nonuniform dust layer could explain the observed water production rates.
In our explanation, the southern regions – including the southern polar cap – become progressively illuminated and the activity starts to increase after the spring equinox (93 days before perihelion, 11 May 2015). This rise in activity allows dust at the surface to be lifted off, decreasing the depth at which water ice is present below the dust layer. This is a runaway process as the increased gas flux resulting from this reduced dust depth is able to lift off increasingly larger dust particles from the surface. An increasing fraction of these particles eventually reaches velocities that are high enough to escape the nucleus gravity, or to be redeposited in other nucleus (northern) regions. This mechanism leads to an increase in the effective active fraction. After the summer equinox (3 days after perihelion, 15 August 2015), the energy input starts to decrease, resulting in a reduced gas flux and surface temperature. The large dust particles can no longer be lifted off from the surface and start to be redeposited locally. The apparition of a dust mantle quenches the cometary activity as the water is no longer at the surface of the comet; instead, it sublimates through a deeper dust layer beneath the surface. The gas diffusion through this dust mantle produces lower gas fluxes, decreasing the effective active fraction of the surface. The process continues for several months until the autumn equinox (224 days after perihelion, 24 March 2016) triggers the southern autumn, followed by the long southern winter around aphelion. At that time, the southern regions are covered again by a dust layer that will be partially removed at the next perihelion passage. It is not entirely excluded either that the outgassing of more volatile species during the northern summer at aphelion contributes to the reaccumulation of material in the southern hemisphere.
This scenario is partially supported by observations from instruments on board Rosetta. Lai et al. (2016) modeled thedynamics of dust grains around 67P taking into account OSIRIS observations. They predicted that dust particles ejected from the southern hemisphere are redeposited in the northern hemisphere during the southern summer. A recent geomorphological analysis of the surface (Birch et al. 2017) also supports the hypothesis of dust being ejected from the southern hemisphere and redeposited in the northern regions, as do several other works such as Thomas et al. (2015), Keller et al. (2017), and Hu et al. (2017). Finally, Fornasier et al. (2016) observed seasonal variations of the color of the nucleus, which becomes bluer near perihelion. This color change is attributed to the removal of the dust mantle covering regions with low gravitational slopes as the comet approaches perihelion.
It must be noted that our model interprets an active area as a literal fraction of the surface area that is covered with water ice and is outgassing. Since exposed water ice is rare on the surface, this is a simplification of the real process, which should involve gas flow through pores, dust layers, erosion, and so on, the details of which we consider beyond the scope of this paper. Nonetheless, an effective active fraction is a useful proxy with which to quantify activity, with the abovenoted spatial and temporal variations providing an insight into the changing activity patterns over the cometary seasons. Our areaweighted global average values of 6.4% around perihelion and 1.9% otherwise agree well with previous estimates for 67P (Lamy et al. 2007) and other comets (see, e.g., A’Hearn et al. 1995), but the large differences between hemispheres and seasons highlight the limitations of interpreting cometary activity with a single number from the ground.
5.2 Refined value for the η parameter
The η parameter, also called momentum transfer efficiency in the literature, represents the fraction of the Maxwellian thermal velocity of the gas, calculated from the nucleus surface temperature, that contributes to the momentum transfer. This parameter depends on the local ice content and on the detailed structure of the porous material. Its value couples the water production curve with the effect of the spin/orbit variations, allowing us to compensate for any systematic effect that might be present in the input data or in the model. The value of the η parameter has been calculated from gaskinetic models of the Knudsen layer. Delsemme & Miller (1971) adopted a value η = 0.6, which lies between the pure ice plane surface value η = 1∕2 and higher values of up to 2∕3 predicted for porous media. In early interpretations of the NGA of comet 1P/Halley, Rickman (1986) and Sagdeev et al. (1988) used values of 0.5 and 0.79, respectively, based on different assumptions. Crifo (1987) recommended a value η ≈ 0.5 based on a revised gaskinetic theoretical description of the solid–gas interface, taking into account the recondensation of water ice. Rickman (1989) used a corrected value η = 0.53–0.67 based on the work of Crifo (1987), which seem to agree well with gas velocity measurements in the coma of comet 1P/Halley (see Rickman 1989, and references therein). However, the calculations did not consider intimate icedust mixtures and did not take into account the porosity of the surface. Skorov & Rickman (1999) introduced a correction factor of 1.8 to the values adopted by Rickman (1989) based on an analytical model of the Knudsen layer above a porous dust mantle. This leads to very high η values in the range 1–1.2.
From our analysis of the data presented in this paper, our best fits are obtained for η in the range 0.6–0.7, in good agreement with the moderate values adopted in the literature. They do not support more extreme values (around 0.5 and greater than 0.8), which degrade the overall fit of the three data sets. We stress, however, that we set the surface gas temperature to the dust temperature T_{dust} in our model (see Eq. (6)). This may lead to an overestimate of the gas temperature, which in turn would tend to underestimate the fitted value of η. We note, however, that the dependence of the gas velocity on the temperature is a square root in Eq. (6): a large and constant temperature deviation is required to significantly bias the value of η. The second point of concern isa possible overestimation of the water production curve through sublimating icy grains in the coma. This would once again require a higher value of η to compensate for the lower local sublimation rates. Altogether, even though our simulations point to η < 0.8, we cannot entirely rule out higher values of this parameter at this point. Consolidated values of the surface water production around perihelion, as well as estimates of the gas velocity above the surface, would help to reduce the uncertainties.
6 Conclusions and perspectives
We reach the following conclusions based on our work:
 1.
We succeeded in finding an activity pattern that simultaneously explains the three following 67P data sets, extracted mainly from the Rosetta mission: (1) the Earthcomet ranging data reconstructed by the flightdynamics team of ESA/ESOC, (2) the water production rate deduced from a mix of ROSINA and groundbased observations (Hansen et al. 2016), and (3) the rate of spinperiod change deduced from the OSIRIS images. The residuals of the ranging data describing the effect of the nongravitational acceleration are reduced by an order of magnitude compared to the groundbased solution based 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 (~10%) outside perihelion compared to the northern ones (<4%), and a drastic rise of the effective active fraction of the southern regions (~25–35%) around perihelion.
 3.
In order to successfully fit the positive rate of spin period change, we split the southern super region into two entities, depending on the sign of the torque efficiency. These two entities correspond to two relatively welldelimited areas in Anhur, Bes, and Khepry, but creates a patchy separation in Wosret.
 4.
We interpreted the timevarying southern effective active fractions by cyclic formation and removal of a dust mantle in these regions. According to our interpretation, the dust mantle could be progressively removed when activity rises after the southern spring equinox and formed again when activity decreases toward the southern autumn equinox (and possibly around aphelion during the northern summer). Several observations performed during the Rosetta mission, such as dust transport from south to north (Lai et al. 2016; Birch et al. 2017) and bluer colours observed near perihelion (Fornasier et al. 2016), support this interpretation.
 5.
If it is confirmed, this interpretation would strongly support postHalley thermal modeling (Rickman et al. 1990; De Sanctis et al. 2010), which predicted that seasonal effects linked to the orientation of the spin axis play a key role in the formation and evolution of dust mantles, and in turn largely control the temporal variations of the gas flux.
 6.
Our analysis supports moderate values of the momentum transfer coefficient η in the range 0.6 − 0.7. For more extreme values of this coefficient (≤0.5 and ≥0.8), the fit of the three data sets is degraded. However, we cannot rule out higher values of this parameter without consolidated water production measurements, and to a lesser extent, without estimates of the nearsurface thermal gas velocity.
More work is required to better understand the activity of 67P using data collected during the Rosetta mission. The solution found in this article through the process described in Sect. 4 is nonunique on the one hand, and it might be improved on the other hand. Improvements may come from a better understanding of the ranging data, resulting in the extraction of clean and accurate nongravitational acceleration of the comet as a function of time around perihelion. The change in the direction of the spin axis or angular velocity might also provide a valuable additional data set that could help to constrain the activity pattern and reduce the number of solutions. In their interpretation of the temporal evolution of the rotational parameters, Kramer et al. (2019) did not introduce a temporal variation of the activity around perihelion, but instead considered a spatially heterogeneous surface with 36 patches with different waterice coverage. They also discussed the possibility of a decreasing dust layer near perihelion with increased activity, and found a solution that explained the water production curve of 67P. It would be interesting to test whether this solution could also reproduce the NGA of comet 67P we described here. Finally, a better understanding of the production rate around perihelion, reconciling the different Rosetta instrument measurements and including species other than water, would be beneficial in fitting the activity model.
Acknowledgements
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 686709. This work was supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number 16.00082. The opinions expressed and arguments employed herein do not necessarily reflect the official view of the Swiss Government. We thank the authors of several python package libraries that were used extensively here; these include NumPy, SciPy, SpiceyPy, REBOUND, and REBOUNDx. We are grateful to L. Maquet for fruitful discussions in the course of this work. We also thank the referee, Dennis Bodewits, for his comprehensive and useful review.
References
 Acton, C. H. 1996, Planet. Space Sci., 44, 65 [NASA ADS] [CrossRef] [Google Scholar]
 A’Hearn, M. F., Millis, R. C., Schleicher, D. O., Osip, D. J., & Birch, P. V. 1995, Icarus, 118, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Bertaux, J.L. 2015, A&A, 583, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birch, S. P. D., Tang, Y., Hayes, A. G., et al. 2017, MNRAS, 469, S50 [CrossRef] [Google Scholar]
 Bodewits, D., Farnham, T. L., Kelley, M. S. P., & Knight, M. M. 2018, Nature, 553, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Brin, G. D. 1980, ApJ, 237, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Brin, G. D., & Mendis, D. A. 1979, ApJ, 229, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Crifo, J. F. 1987, A&A, 187, 438 [NASA ADS] [Google Scholar]
 Davidsson, B. J. R., & Gutiérrez, P. J. 2004, Icarus, 168, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Davidsson, B. J. R., & Gutiérrez, P. J. 2005, Icarus, 176, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Davidsson, B. J. R., & Gutiérrez, P. J. 2006, Icarus, 180, 224 [NASA ADS] [CrossRef] [Google Scholar]
 Davidsson, B. J. R., Gutiérrez, P. J., & Rickman, H. 2007, Icarus, 191, 547 [NASA ADS] [CrossRef] [Google Scholar]
 De Sanctis, M. C., Lasue, J., Capria, M. T., et al. 2010, Icarus, 207, 341 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delsemme, A. H., & Miller, D. C. 1971, Planet. Space Sci., 19, 1229 [NASA ADS] [CrossRef] [Google Scholar]
 Drahus, M., & Waniak, W. 2006, Icarus, 185, 544 [NASA ADS] [CrossRef] [Google Scholar]
 Fanale, F. P., & Salvail, J. R. 1984, Icarus, 60, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Fornasier, S., Hasselmann, P. H., Barucci, M. A., et al. 2015, A&A, 583, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fornasier, S., Mottola, S., Keller, H. U., et al. 2016, Science, 354, 1566 [NASA ADS] [CrossRef] [Google Scholar]
 Fornasier, S., Feller, C., Lee, J.C., et al. 2017, MNRAS, 469, S93 [CrossRef] [Google Scholar]
 Gutiérrez, P. J., & Davidsson, B. J. R. 2007, Icarus, 191, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Gutiérrez, P. J., Ortiz, J. L., Rodrigo, R., & Ló pezMoreno, J. J. 2001, A&A, 374, 326 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gutiérrez, P. J., de León, J., Jorda, L., et al. 2003, A&A, 407, L37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gutiérrez, P. J., Jorda, L., Samarasinha, N. H., & Lamy, P. 2005, Planet. Space Sci., 53, 1135 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, K. C., Altwegg, K., Berthelier, J.J., et al. 2016, MNRAS, 462, S491 [Google Scholar]
 Hu, X., Shi, X., Sierks, H., et al. 2017, A&A, 604, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jorda, L., Spjuth, S., Keller, H. U., Lamy, P., & Llebaria, A. 2010, in Computational Imaging VIII, Proc. SPIE, 7533, 753311 [Google Scholar]
 Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, H. U., Delamere, W. A., Reitsema, H. J., Huebner, W. F., & Schmidt, H. U. 1987, A&A, 187, 807 [NASA ADS] [Google Scholar]
 Keller, H. U., Mottola, S., Skorov, Y., & Jorda, L. 2015, A&A, 579, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keller, H. U., Mottola, S., Hviid, S. F., et al. 2017, MNRAS, 469, S357 [CrossRef] [Google Scholar]
 Knight, M. M., Farnham, T. L., Schleicher, D. G., & Schwieterman, E. W. 2011, AJ, 141, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Kossacki, K. J., & Szutowicz, S. 2006, Planet. Space Sci., 54, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, T., Läuter, M., Hviid, S., et al. 2019, A&A, 630, A18, (Rosetta 2 SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lai, I.L., Ip, W.H., Su, C.C., et al. 2016, MNRAS, 462, S533 [CrossRef] [Google Scholar]
 Lamy, P. L., Toth, I., Davidsson, B. J. R., et al. 2007, Space Sci. Rev., 128, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Maquet, L., Colas, F., Jorda, L., & Crovisier, J. 2012, A&A, 548, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marschall, R., Su, C. C., Liao, Y., et al. 2016, A&A, 589, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marschall, R., Mottola, S., Su, C. C., et al. 2017, A&A, 605, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marsden, B. G. 1968, AJ, 73, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, B. G., Sekanina, Z., & Yeomans, D. K. 1973, AJ, 78, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Mottola, S., Lowry, S., Snodgrass, C., et al. 2014, A&A, 569, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mueller, B. E. A., & Ferrin, I. 1996, Icarus, 123, 463 [NASA ADS] [CrossRef] [Google Scholar]
 Newhall, X. X., Standish, E. M., & Williams, J. G. 1983, A&A, 125, 150 [NASA ADS] [Google Scholar]
 Pajola, M., Lucchetti, A., Vincent, J.B., et al. 2016, A&A, 592, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Preusker, F., Scholten, F., Matz, K.D., et al. 2017, A&A, 607, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Liu, S.F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H. 1986, in Comet Nucleus Sample Return Mission, ed. O. Melita, ESA SP, 249 [Google Scholar]
 Rickman, H. 1989, Adv. Space Res., 9, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H., Fernandez, J. A., & Gustafson, B. A. S. 1990, A&A, 237, 524 [NASA ADS] [Google Scholar]
 Sagdeev, R. Z., Elyasberg, P. E., & Moroz, V. I. 1988, Nature, 331, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Samarasinha, N. H., & A’Hearn, M. F. 1991, Icarus, 93, 194 [NASA ADS] [CrossRef] [Google Scholar]
 Samarasinha, N. H., Mueller, B. E. A., Belton, M. J. S., & Jorda, L. 2004, Comets II, eds. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson, AZ: University of Arizona Press), 281 [Google Scholar]
 Sekanina, Z. 1993, A&A, 277, 265 [NASA ADS] [Google Scholar]
 Skorov, Y. V., & Rickman, H. 1999, Planet. Space Sci., 47, 935 [NASA ADS] [CrossRef] [Google Scholar]
 Standish, E. 1998, IOM, 312.F98048 [Google Scholar]
 Szutowicz, S. 2000, A&A, 363, 323 [NASA ADS] [Google Scholar]
 Thomas, N., Sierks, H., Barbieri, C., et al. 2015, Science, 347, aaa0440 [CrossRef] [Google Scholar]
 Voglis, C., & Lagaris, I. E. 2018, A Rectangular Trust Region Dogleg Approach for Unconstrained and Bound Constrained Nonlinear Optimization, WSEAS International Conference on Applied Mathematics, Corfu, Greece, 2004. [Google Scholar]
 Whipple, F. L. 1950, ApJ, 111, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Yeomans, D. K., & Chodas, P. W. 1989, AJ, 98, 1083 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Parameters and rootmeansquared residuals (in range, water production rate, nongravitational torque, and the total objective function) for the bestfit models.
All Figures
Fig. 1 Discontinuities identified in the position of comet 67P from the SPICE kernels, in (x, y, z) heliocentric J2000 coordinates, and Earthcomet range (R). Left panel: as a function of time. Right panel: as a histogram of sizes. 

In the text 
Fig. 2 Observed torque, derived from the 67P rotation state from OSIRIS measurements, and a smoothed cubic spline fit to the data. The gray region represents the RMS of the residuals between the two. 

In the text 
Fig. 3 Torque efficiency in meters, a factor determined by the geometry as defined in Eq. (8). The cometary rotation axisis in the + z direction through the neck region in the center. 

In the text 
Fig. 4 Map of the five “super regions” defined by Marschall et al. (2016, 2017) and used in our solution A optimization. 

In the text 
Fig. 5 Observed Earthcomet range, R, and residuals for two models: a purely gravitational Nbody solution with no NGA, and a groundbased solution based on the model of Marsden et al. (1973; see text for details). Differences between the observed and computed solutions, as well as the jumps described in Sect 2.2, are invisible at the scale of the top plot. 

In the text 
Fig. 6 Observed minus computed range, R, for solution A (blue curve). For comparison, the residuals to the groundbased solution, using the model of Marsden et al. (1973) as in Fig. 5, are shown in orange. Both curves diverge from zero most sharply following the maximum perturbation around perihelion, but our solution is an improvement. 

In the text 
Fig. 7 Observed and computed water production rates and residuals for solution A. 

In the text 
Fig. 8 Observed and computed torques and residuals for solution A. 

In the text 
Fig. 9 Mapped active fraction for solution A. 

In the text 
Fig. 10 Observed minus computed range for solution B (blue curve), and the groundbased solution. 

In the text 
Fig. 11 Observed and computed water production rates and residuals for solution B. 

In the text 
Fig. 12 Observed and computed torques and residuals for solution B. 

In the text 
Fig. 13 Mapped active fraction for solution B. 

In the text 
Fig. 14 Observed minus computed range for solution C (blue curve), and the groundbased solution. 

In the text 
Fig. 15 Observed and computed water production rates and residuals for solution C. 

In the text 
Fig. 16 Observed and computed torques and residuals for solution C. 

In the text 
Fig. 17 Mapped final active fraction for solution C. 

In the text 
Fig. 18 Active fraction with time for solution C. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.