© N. Attree et al. 2019

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₁, A₂, and A₃ 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/Wolf-Harrington 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 ground-based 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, self-shadowing, self-heating, 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/Churyumov-Gerasimenko (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 non-principal 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 non-principal 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/Churyumov-Gerasimenko (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 ground-based 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 non-trivial. 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 ground-based observations (see, e.g., Bertaux 2015) have also suggested a variation in peak production between perihelion passages. Production rate estimates from the Rosetta line-of-sight 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₂, CO, and O₂. 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 non-gravitational forces.

2.2 Trajectory

Outgassing produces a back-reaction 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 radio-tracking 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 Earth-comet range direction (R, claimed accuracy of ~10 m) and a much lower accuracy in the perpendicular (cross-track) 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 non-gravitational 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 N-body gravitational interactions, in order to constrain said model. To do so, we used the magnitude of the comet-to-Earth-center 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 Earth-comet range (R). Left panel: as a function of time. Right panel: as a histogram of sizes.

2.3 Rotation

Back-reaction 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 body-fixed frame) by Eq. (1), where I_z = 1.899 × 10¹⁹ kg m² is the third (largest) moment of inertia derived from the shape model assuming a constant density of 538 kg m⁻³ (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 non-gravitational forces acting on the nucleus. Our thermal model takes into account solar insolation, surface thermal emission, sublimation of water ice, projected shadows, and self-heating. 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⁻² K⁻¹ s^−0.5. The production rate and the acceleration are ~7% smaller for a thermal inertia of 100 J m⁻² K⁻¹ 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⁻² is the solar constant, z_i is the zenithal angle, r_h the heliocentric distance, SH_i is the self-heating 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⁶ J Kg⁻¹ 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 self-heating 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¹² Pa and B = 6162 K for water (Fanale & Salvail 1984), and the gas constant R = 8.3144598 J K⁻¹ mol⁻¹, (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 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 Non-gravitational 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¹² 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 non-gravitational 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⁻¹.

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 body-fixed 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 large-scale 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).

Non-gravitational 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 time-varying 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 N-body integration must be performed and the resulting position compared at each time (C_R). To do this, we used the open-source REBOUND code¹ (Rein & Liu 2012), complete with full general relativistic corrections (Newhall et al. 1983), as implemented by the REBOUNDx extension package². 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 comet-to-Earth-center 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 least-squares 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, O-C, and observation uncertainties, σ, at each time-step, 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 multi-objective 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 best-fit 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 N-body 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, along-orbit, and normal-to-orbit) based on a power law with heliocentric distance, and three scaling parameters A_1,2,3 found by a formal best-fit to the orbit for each comet. We used the A_1,2,3 values for the 2010 apparition of 67P from ground-observations given by NASA/JPL Horizons ephemerides³.

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 ground-based observations of roughly arcsecond accuracy. This stands as a baseline against which we can check our own activity model.

Fig. 5 Observed Earth-comet range, R, and residuals for two models: a purely gravitational N-body solution with no NGA, and a ground-based 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 Five-parameter solution (A)

Figures 6 –8 show the residuals to our best-fit 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 ground-based solution, demonstrating the significance of our NGA/N-body 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 ground-based 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.

Fig. 9 Mapped active fraction for solution A.

4.2 Twenty-six-parameter 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 post-perihelion.

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 non-contiguous 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 ground-based 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 Time-varying solution (C)

Because the above solutions with constant active fractions failed to adequately reproduce all the data, we now consider time-varying solutions. Temporal variations of the effective active fraction is an obvious way to try to reconcile the modeled post-perihelion 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 time-varying 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 half-Gaussian fall-off 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₀ and t₁, and two decay (or growth) half-times, and , for a total of 12 fitted parameters. The relative magnitudes of the two active fractions were unconstrained, but the two times t₀ and t₁ were constrained to lie either side of perihelion to ensure that they did not cross over. The half-times 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 multi-objective 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 multi-objective 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 best-fit 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 drop-off 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 half-time of ~ 25 days) to their high peak values, before falling off to the ~10% level after perihelion and during southern summer, with a decay half-time of around 50 days. This reproduces the high production rate at perihelion while matching the swift fall-off 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 non-gravitational 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 time-varying 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 best-fit results allow us to calculate integrated quantitative parameters resulting from the cometary activity around perihelion. The total water-ice mass loss amounts to 4.5 × 10⁶ kg, corresponding to a mean erosion of ~ 9 cm over the entire nucleus surface, assuming a density of 538 kg m⁻³ (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 dust-to-ice ratio – and the outgassing of more volatile minor species throughout the orbit. The actual water-ice erosion can be calculated by integrating the time-varying 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 ground-based 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 non-uniform 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 above-noted spatial and temporal variations providing an insight into the changing activity patterns over the cometary seasons. Our area-weighted 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 gas-kinetic 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 gas-kinetic 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 ice-dust 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 Earth-comet ranging data reconstructed by the flight-dynamics team of ESA/ESOC, (2) the water production rate deduced from a mix of ROSINA and ground-based observations (Hansen et al. 2016), and (3) the rate of spin-period change deduced from the OSIRIS images. The residuals of the ranging data describing the effect of the non-gravitational acceleration are reduced by an order of magnitude compared to the ground-based 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 well-delimited areas in Anhur, Bes, and Khepry, but creates a patchy separation in Wosret.

• 4.

  We interpreted the time-varying 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 post-Halley 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 near-surface 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 non-unique 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 non-gravitational 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 water-ice 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.0008-2. 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

All Tables

All Figures

	Fig. 1 Discontinuities identified in the position of comet 67P from the SPICE kernels, in (x, y, z) heliocentric J2000 coordinates, and Earth-comet 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 Earth-comet range, R, and residuals for two models: a purely gravitational N-body solution with no NGA, and a ground-based 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 ground-based 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 ground-based 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 ground-based 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