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



Article Number  A3  
Number of page(s)  11  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834349  
Published online  20 September 2019 
Comet 67P/ChuryumovGerasimenko rotation changes derived from sublimationinduced torques
^{1}
Zuse Institute Berlin, Supercomputing Division,
Takustr. 7,
14195 Berlin, Germany
email: kramer@zib.de
^{2}
Department of Physics, Harvard University,
17 Oxford St,
Cambridge,
MA 02138, USA
^{3}
Laboratoire d’Astrophysique de Marseille, UMR7326 CNRS/Université AixMarseille,
38 rue Frédéric JoliotCurie,
13388 Marseille Cedex 13, France
^{4}
Institute of Planetary Research, Deutsches Zentrum für Luft und Raumfahrt (DLR),
Rutherfordstrasse 2,
12489 Berlin, Germany
Received:
30
September
2018
Accepted:
7
January
2019
Context. The change in rotation period and the orientation of the rotation axis of comet 67P/ChuryumovGerasimenko (67P) can be deduced with high precision from images taken by the scientific imaging instruments on board the Rosetta mission. Nongravitational forces are a natural explanation for these data.
Aims. We describe observed changes in orientation of the rotation axis and the rotation period of 67P. We explain them based on a sublimation model with a best fit for the surface active fraction (model P). Torque effects of periodically changing gas emissions on the surface are considered.
Methods. We solved the equation of state for the angular momentum in the inertial and the bodyfixed frames and provide an analytic theory of the rotation changes in terms of Fourier coefficients, which are generally applicable to periodically forced rigidbody dynamics.
Results. The torqueinduced changes in rotation state constrain the physical properties of the surface, the sublimation rate, and the local active fraction of the surface.
Conclusions. We determine a distribution of the local surface active fraction in agreement with the rotation properties, period, and orientation of 67P. The torque movement confirms that the sublimation increases faster than the insolation toward perihelion. The derived relatively uniform activity pattern is discussed in terms of related surface features.
Key words: comets: general / comets: individual: 67P/ChuryumovGerasimenko / methods: analytical
© ESO 2019
1 Introduction
The Rosetta instruments have probed the gas and dust environment during almost the entire apparition of 67P/ChuryumovGerasimenko (67P) from 2014 to 2016. The observed changes in rotation period and axis orientation provide an independent measure of the sublimation activity of the nucleus by the induced torque. The rotation period of 67P was shortened by 21 min. The same amount has been reported for the previous apparition by Mottola et al. (2014). Keller et al. (2015a) proposed a sublimationdriven model and explained the changing rotation period in terms of a homogeneous ice model of the entire surface (Keller et al. 2015b) without studying the axis orientation. For many comets, changes in rotation axis and of the orbital elements due to activity are observed and predicted (see Whipple 1950; Jewitt 1997; Samarasinha et al. 2004; Mueller & Samarasinha 2018). Before the shape of 67P was known in detail, Gutiérrez et al. (2003) explored several scenarios for rotationaxis changes of small, irregularly shaped comets. Typical changes in rotation axis caused by sublimation forces range from 0.1 to several tens of degrees. For 67P, the observed change is at the lower end of this range (0.5°).
For model A proposed by Keller et al. (2015a), our analysis predicts a change greater by five times of the direction of the rotation axis. This change is furthermore in a different direction compared to the observations. In addition, model A predicts a slowerthanobserved increase in total gas production with decreasing heliocentric distance, while the analysis of the coma shows a faster increase in activity based on the pressure sensors (see Hansen et al. 2016; Kramer et al. 2017; Läuter et al. 2019). To overcome these discrepancies, we established a formalism to match models and observations in terms of a Fourier analysis of the gasinduced torque and derived a possible ice distribution on the surface that explains the rotationperiod changes, the movement of the direction of the angular momentum, and the increase in activity with heliocentric distance (model P).
2 Forced rigidbody dynamics
We reviewed the response of the rotation state of the nucleus to sublimation and other processes that alter the rotation and motion, see Thomson (1986). The comet is viewed as a moving and rotating rigid body. The cometary nucleus is defined within the threedimensional body frame by a prescribed bodyframe density ρ_{bf}(x, t) for each x in the body frame at time t to accommodate slow changes within the internal mass distribution. The body frame is linked to the inertial frame by the coordinate transformation at time t (1)
for each bodyframe point x. Here, r(t) denotes the center of figure, R(t) the orthogonal rotation matrix with the property (2)
and ω(t) the angular velocity. The mapping x′(x, t) reflects the movement of the center of figure with time t, but leaves the shape geometry (defined in the body frame) unchanged. The timedependent density ρ_{bf} (x, t) allows incorporating slow changes in the density and porosity of the comet. The density in the inertial frame ρ(x′, t) = ρ_{bf}(x, t) is obtained from the bodyframe density considering Eq. (1) and carries along the time dependence of the orbital and rotational movement. The comet mass M and the center of mass in the body frame are in general time dependent:
For the center of mass in the inertial frame, the relation holds that (5)
The time derivative of Eq. (1) yields the inertialframe velocity V of a fixed bodypoint (6)
To obtain the linear momentum P and angular momentum L of the whole nucleus in the inertial system, we integrate
with the tensor of inertia I(t) = RI_{bf}(t)R^{−1} in the inertial frame with respect to the center of figure r and the tensor of inertia I_{bf}(t) with respect to the bodyframe center 0. For the case of a timedependent body density, I_{bf}(t) needs to be computed for the bodyframe density at time t. The momentum changes are generated by the sum of nongravitational and gravitational forces and torques . Gas sublimation at point x on the surface leads to a mass loss ṁ and generates the nongravitational components
The gas velocity in the inertial frame consists of two components, one into normal direction on the nucleus surface in the inertial frame, and another due to the body rotation, u_{gas} denotes the thermal gas velocity from Eq. (22). The gas velocity in the body frame is given by with the outward surface normal at surface location x within the body frame. According to Jorda & Gutiérrez (2002), Eq. (11), the mass production ṁ for a mixture of gas species reads (11)
with the surface active fraction f_{gas} and the sublimation rate Z_{gas}. The mass loss ṁ changes the shape, reduces the total mass, and affects the tensor of inertia of the nucleus. The integrated mass loss of comet 67P during the 2015 apparition is estimated to be about 1/1000 of the total mass (see Godard et al. 2015, 2017), and therefore we neglect botheffects and assume a timeindependent mass and tensor of inertia. The gravitational components for force and torque yield
with the gravitational acceleration a due to other solar system bodies. For both volume integrals in Eqs. (12) and (13), a(x′) is assumed to be constant over the cometary body, which implies that tidal forces, for instance, are neglected.
Equations (2), (7)–(10), (12), and (13) result in a system of coupled algebraic and differential equations,
for the state variables r(t), P(t), R (t), and L(t). Equations (7) and (8) couple linear and rotational momenta through the center of mass in the body frame . If the density distribution ρ_{bf} satisfies , becomes thecenter of mass in the inertial frame and Eq. (14) decouple into two blocks. The first block, (16)
describes the translational movement for the state variables r(t), P(t), and the second block, (17)
describes the rotational dynamics for R(t) and L(t). The model for the changing rotation period by Keller et al. (2015a) is contained as a special case in Eq. (17). For this, we denote the eigenvector e of I_{bf} in the body frame with the largest moment of inertia I_{z} and assume that L is initially aligned with e′ = Re, that is, L(0) = L_{z}e′(0). Then ω(0) = L_{z}∕I_{z}e′(0) and consequently . Thus, e′ = R e is constant in time, and the angular velocity changes with .
For given initial conditions of all state variables, the system of Eqs. (14), (16), and (17) can be solved numerically. For accuracy, we used the LSODE package provided within Mathematica/FORTRAN (Radhakrishnan & Hindmarsh 1993). Following Shoemake (1985), the matrix–matrix operation R →ω × R was replaced by a quaternion multiplication for improved stability. There is a onetoone mapping between R and a quaternion q such that the matrix–matrix operation is substituted by q →q ⋅ (0, ω)∕2.
3 Rotation state of 67P
After the arrival of Rosetta at 67P in 2014, observations of the rotation by Preusker et al. (2015) and Godard et al. (2017) showed 67P in an excited state of rotation, albeit with the rotation axis close to the axis e′ with the largest moment of inertia (rotation state with minimum energy). This points to an alignment of the two axes, which is also compatible with observations (Jorda et al. 2016). The total mass was estimated to be 10^{13} kg by Godard et al. (2015). The assumption of a strictly homogeneous density leads to a tensor of inertia (18)
with respect to the center of mass. This tensor is incompatible with the observations, since then the axis e′ would be tilted by 2.9° with respect to the rotation axis (Preusker et al. 2017), and in addition, an offset of the center of mass with respect to the center of figure r would exist (Jorda et al. 2016). Inhomogeneities in the density have also been reported by Brouet et al. (2016) and Knapmeyer et al. (2018) from the CONSERT and SESAME/MUPUS Rosetta data. In the following, we assume a timeindependent, nonhomogeneous density distribution that aligns the rotation axis with the e′ under the constraint that the total mass iskept fixed at M = 10^{13} kg and . For clarity, we give a possible mass distribution with a resulting tensor of inertia (19)
This solution equals placing (1∕11)M on a thin ring centered at {−159.8, 275.5, −220.5} m in the plane − 0.4065601x + 0.01699493y + 0.9134659z = −131.7705 m with radius 1 km, while distributing (10∕11)M homogeneously throughout the entire nucleus. The solution is not unique and leads to an increased density in the big lobe, in agreement with Jorda et al. (2016). The small offdiagonal entries in Eq. (19) align the shape file with the body frame by a rotation of 0.4°. The changes in rotation state of the comet is determined by the relative change in angular momentum with respect to the initial state. Changing the total mass does not affect the resulting dynamics if the average surface active fraction is changed in the same proportion.
3.1 Sublimation
Changes in rotation state are primarily due to the torque induced by sublimation of ice. The total production of water and CO_{2} has been estimated by Hansen et al. (2016) and Läuter et al. (2019) from ROSINA COPS and DFMS data to be about 6.2 ± 2 × 10^{9} kg, corresponding to about 1∕1600 of the total mass of 67P (M = 10^{13} kg). The water production steeply increases with heliocentric distance around southern solstice, with exponents α ranging from − 6.5 up to − 7. The total gas production from the radiationdriven sublimation model A by Keller et al. (2015a) yields smaller exponents α = − 2.8. Our model for the rotation state only considers water emission from the surface for driving the torque evolution. The CO_{2} activity liberates decimetersized chunks (Keller et al. 2017) that contain additional water, which is seen as production but does not contribute to the torque and does not have the same diurnal signature as the surface. The CO_{2} contribution(about 1/7 of the water mass estimated from ROSINA/DFMS by Läuter et al. 2019) is not considered separately since the CO_{2} sources around perihelion coincide with the water regions (see Fougere et al. 2016a and Läuter et al. 2019)and drive the torque in a similar direction as water. In addition, the diurnal variation of CO_{2} is less pronounced than for water (see Filacchione et al. 2016) and thus has less influence on the axis orientation, as discussed in Sect. 4. At heliocentric distances larger than 3 au, and in particular on the outbound orbit, CO_{2} becomes thedominant species (Läuter et al. 2019) and does not follow the subsolar illumination. At these distances the rotation period of the comet has settled, and these times are beyond the scope of the present analysis. The rotation axis of 67P shows the largest movements around ± 100 days from perihelion, in agreement with the larger α exponents derived from the gas instruments. Keller et al. (2015a; model A) considered the gas production based on a shape model. On each surface element with a given Bond bolometric albedo A and solar irradiance at heliocentric distance r_{h} (au), solar constant S_{⊙} = 1361 W m^{−2}, the energy balance (20)
is solved for the sublimation rate Z, given by the HertzKnudsen relation (21)
The parameters are taken from Keller et al. (2015a), with emissivity ϵ = 0.9, latent heat of sublimation for water ice L_{ice} = 2.6 × 10^{6} J kg^{−1} (assumed to be constant), water vapor pressure P(T) = 3.56 × 10^{12} e^{−6141.667∕T} (kg m^{−1} s^{−2}), and thermal velocity of water molecules with molar mass (22)
The gas constant is denoted by R and the StefanBoltzmann constant by σ. The solution to Eqs. (20) and (21) in terms of the sublimation rate (s^{−1} m^{−2}) and mass of water molecule (kg), is shown in Fig. 1 for A = 0.01 (Keller et al. 2015a, model A). The changes observed by Hansen et al. (2016); Kramer et al. (2017); Läuter et al. (2019) point to a faster increase in total production with decreasing perihelion distance. To account for the observation requires assuming a sublimation rate that increases faster than linearly with illumination, as exemplified by the dashed red line in Fig. 1. A possible physical mechanism behind the increase in sublimation could be a decrease in dust layer when the comet approaches the Sun, leading to a steeper slope (transition of model C to model A, as described by Keller et al. 2015a). The required adjustment of the sublimation rate with heliocentric distance has also been suggested by Marsden et al. (1973) and is used in the DSMC coma models for 67P by Fougere et al. (2016b), Eq. (1).
3.2 Observed rotation axis changes
The (RA, Dec) orientation of the angular velocity ω of the comet nucleus was derived using the set of about 25 000 control points defined as the center of the maplets created in the stereophotoclinometry (SPC) method of Gaskell et al. (2008) applied to comet 67P by Jorda et al. (2016). The coordinates of the stereo control points measured on sequences of Rosetta/OSIRIS images (Keller et al. 2007) combined with startracker pointing measurements were used to determine the direction of the angular velocity vector in the equatorial J2000 (EME2000) reference frame during the Rosetta mission, see Fig. 2. The fluctuations in the resulting data set are caused by a possible nutation combined with the uncertainties in the determination of the direction of the ω vector. They led us to strongly smooth the data, as illustrated by the solid line in Fig. 2, which represents the timeaveraged motion of the rotation axis by fitting the data to a Gaussian function for the RA and to a hyperbolic tangent function for the Dec, in addition to a thirdorder polynomial. A similar data set has been retrieved by Godard et al. (2017).
Fig. 1
Radiationdriven sublimation (a) rate and (b) velocity as function of received irradiation. Solid blue lines show model A (surface active fraction = 1) and dashed red lines the effective sublimation curve with enhanced radiation response near perihelion. The effective sublimation curve leads to agreement with observations (see Figs. 6 and 7). 
4 Fourier theory of torqueinduced motion
The shape model from Preusker et al. (2017) was remeshed using the approximated centroidal Voronoi diagrams (ACVD) tool by Valette et al. (2008) into N_{faces} = 3996 triangular elements with area A_{i} and surface normal , chosen to be approximately of equal area. The results are robust against variations in shape model as long as N_{faces} ≫ 100. The torque was evaluated according to Eqs. (10) and (11). The sublimation rate Z_{gas} on each facet was evaluated for a given subsolar latitude ϕ_{s} and heliocentric distance r_{h} during one rotation period, starting with subsolar longitude λ_{s} = 0. We computed the total torque arising from the water sublimation curves, either the HertzKnudsen rate from Eq. (21), or alternatively, the effective sublimation curve in Fig. 1: (23)
The contributions of shadowed (including selfshadowed) faces were set to zero and the entire surface active fraction was first set to the same value. In the next iteration, the surface active fraction was considered to be a spatial fit parameter (but fixed in time) to be determined from the observed axis changes. The heliocentric orbit (Cartesian coordinates r_{h} (t)) was taken from the NASA Horizon system (Earth mean equator and equinox of reference epoch J2000). The initial orientation of the rotation axis in the equatorial frame was set to be in the lowest energy state (rotation axis and angular momentum aligned, pointing to RA α = 69.427°, Dec δ = 64.000° at t = −350 days). For 67P, the observed axis changes are small, and we tabulated the subsolar latitude ϕ_{s} and heliocentric distance r_{h} in N_{intervals} = 81 10day intervals and stored the bodyframe torque as function of λ_{s}. A typical diurnal torque evolution at perihelion is shown in Fig. 3. To gain more physical insight into the dynamics underlying the axis changes, we expanded the torque components into a Fourier cosine/sine series. The periodic argument α = [−π : π[ of the Fourier series is not time, but the subsolar longitude λ_{s} = α + π to accommodate changes in the rotation period. In the ith time interval, we extracted the first N = 2m + 1 Fourier coefficients , which yield the Fourier series representation of the torque (24)
The final parameterization of the body torque as function of subsolar longitude along the entire orbit is given by the time evolution of the Fourier coefficients (25)
The slowly changing subsolar latitude and heliocentric distance are implicitly contained through the time t argument. The time evolution of the coefficients is shown in Fig. 4 for N = 3. The first row in Fig. 4 displays the nondiurnal torque coefficients. Because the rotation axis is aligned with the z bodyaxis, the torque component directly affects the rotation period. The orientation of the rotation axis changes by the action of the diurnal components (Fig. 3). The diurnal components are not linearly independent, as discussed in Sect. 4.1.
4.1 Computation of the torque in the inertial system
First, we neglected the changes in orbital elements that are due to nongravitational momentum (Eq. (16)) and took the orbital evolution r_{h} (t) as fixed. For a given rotation matrix R(t) (body frame to equatorial inertialsystem) and position of the comet r_{h} (t), the subsolar longitude is given by
At initial time t_{S}, the rotation matrix R(t_{S}), which transfers the bodyframe zaxis to the rotation axis in the equatorial inertial frame, is given by (29)
Equation (17) is then integrated with the initial angular velocity and momentum set to
The Fourier decomposition provides additional insights into the axis changes. An important parameter is the angle around the zaxis to point the x–zplane toward the Sun (zero subsolar latitude). When the 0.5° tiltchange of the rotation axis is neglected, is given by
We obtain a good approximation of the angular momentum change during one rotation period T_{rot} = 2π∕ω by keeping fixed during this rotation and integrating the torque in the body frame T_{bf}, parameterized by the subsolar longitude and the Fourier components from Eq. (26), (34)
The angular momentum change along the entire orbit is then approximated by adding all contributions to the initial angular momentum. The “orbit” matrix does not affect the magnitude of the “shape” vector. All shown results were obtained without this approximation and use the full numerical solution of Eq. (17). Equation (34) was used to determine the physically relevant Fourier components, (35)
for analyzing the observations.
Fig. 2
Orientation changes in rotation axis. Black dots denote the reconstructed right ascension multiplied with cos 64° and declination, and the solid line shows a smooth fit to the data. 
Fig. 3
Calculated sublimationinduced torque at perihelion for one rotation period at perihelion as a function of subsolar longitude from model P using the effective sublimation curve from Fig. 1 with the bestfit solution from Fig. 10. The dashed lines show the Fourier representation (constant, cos, and sin terms) of the corresponding solid lines to parameterize the diurnal cycle. 
4.2 Extract observed torque from the rotationaxis evolution
Next, we considered the inverse problem of finding a plausible torque function in the cometary body frame as function of subsolar coordinates and solar distance. We inferred the torque in the body frame from the observation as function of time t under the assumption of an initial alignment of rotation axis and angular momentum, and with the tensor of inertia given in Eq. (19).
To parameterize the observed torque as function of subsolar longitude, we computed λ_{s} (t) from Eq. (27) at each observation time. Every 10 days, we determined the closest instance t_{i} of λ_{s} (t_{i}) = 0 and computed the Fourier coefficients to represent T_{bf}(λ_{s} = 0…2π). Only the three Fourier combinations from the components of Eq. (35) should be retrieved from the fit (Fig. 5) because the axis motion is not sensitive to the other Fourier components (see Eq. (34)).
Fig. 4
Torque evolution of a uniform active surface model with the effective sublimation curve from Fig. 1. The time evolution of the first three Fourier coefficients , , and is shown for (−300 : 300) days from perihelion for each Cartesian component of the bodyframe torque, Eq. (25). Units of 10^{6} kg m^{2} s^{−2}. 
Fig. 5
Time evolution of the three physical relevant combinations of Fourier coefficients (Eq. (35)) for (−300 :300) days from perihelion, units of 10^{6} kg m^{2} s^{−2}. (a) Model A, global uniform active surface 1∕12. (b) Model P (patches with effective sublimation curve). (c) Observed Fourier coefficients inferred from the rotationaxis movement and the tensor of inertia. 
5 Matching Fourier coefficients with the observed torque
The simplest sublimation model A results in a rotation axis movement shown in Fig. 6 with the green line. The evolution is rotated by 90° with respectto the observed torque movement (Fig. 6, red line) and leads to a largely increased axis tilt compared to observations. To explain the observations requires considering a spatially heterogeneous surface with varying waterice coverage (model P). We could show that an alternative explanation is a large thermal lag of several hours of the maximum sublimation with respectto the maximum irradiation caused by a dust layer of some millimeters thickness. However, this scenario is unlikely because the response of sublimation rates to radiation is almost immediate, as seen by short (< 1 h) delays of jet outbreaks (see Lai et al. 2016; Shi et al. 2016) and the good agreement of inner dust structures with illuminationdriven dust release (see Kramer & Noack 2016; Kramer et al. 2018). Measurements of VIRTIS and MIRO found that the thermal inertia is lower than 320 JK^{−1}m^{−2}s^{−0.5} when the error bars are included (see Marshall et al. 2018 for an overview). Thermal inertia this low is not able to provide the needed phase lag of several hours of the maximum sublimation with respect to the maximum irradiation, as shown by thermal simulations and measurements of the activity maxima compared to noon time, see Shi et al. (2016). We show that the ofteninvoked unrealistic thermal lag to explain the nongravitational forces acting on the cometary orbit can be at least partly replaced by the effects of a complex nucleus shape and its slightly nonuniform activity (see Davidsson & Gutierrez 2005; Sosa & Fernández 2009).
For the nonuniform case, we divided the surface into 36 equally spaced patches and computed their separate contributions to the torque using the Fourier method described before. Each patch provides a specific contribution to the Fourier components C_{I}, C_{II}, and C_{III} in Eq. (35) of the complete comet. For a uniform activity, the resulting extrema of the Fourier components are shown in Fig. 8 for each patch. To match the observed rotation state, a linear combination of the patch contributions must yield the observed values of C_{I}, C_{II}, and C_{III} in Fig. 5, indicated by the dashed lines in Fig. 8. The relative ratio of the three components for a single patch is a prescribed property of the sublimation curve. The largest difference of a singlepatch contribution to the observation is that for component C_{I} on patch 21. The activity of patch 21 has to be reduced, while patches 26–36 with opposite sign for C_{I} are candidates for an increased activity. Additional constraints on the activity arise from the simultaneous fitting of components C_{II} and C_{III}. To find the activity across all patches, we minimize the deviation of observed torque and observations every 20 days with respect to the L^{1} norm. Solving a nonlinear optimization problem does not always provide a unique solution. By changing the norm and using different time slices, we verified that the resulting solution with the identified depleted and enhanced surface active regions remains unaffected by the data selection. The fit leads to a closer alignment of observation and model A/patches for the axis movement (Fig. 6), but does not fix the exponent of the total production rate, Fig. 7, which remains at Q_{tot}(r) ~ r^{−2.8}. In contrast, observations from COPS/DFMS point to a larger exponent α approximately −6 to − 7. The change in sublimation with heliocentric distance is directly reflected by a small southern excursion of the rotation axis (300–100 days) before perihelion. The observations show that the sublimation activity increases nonlinearly with insolation, as discussed in Sect. 3.1. The effective sublimation curve in Fig. 1a, dashed line, yields the total production displayed in Fig. 7, with larger exponent α < −5 as measuredby several Rosetta instruments (see, e.g., Hansen et al. 2016; Kramer et al. 2017; Läuter et al. 2019) and modeled by Hu et al. (2017). The rotationaxis motion of this modified sublimation model is shown in Fig. 6 and agrees better with observations than the other considered scenarios.
Fig. 6
Rotationaxis movement. The red line shows the observation from Fig. 2, and the other linesrepresent different sublimation models. Model A with a globally constant surface active fraction (1∕12), model A/patches with bestfit adjustment of patches, and model P with effective sublimation curve and bestfit adjustment of patches. The gray inset shows a magnification of the curves. 
Fig. 7
Rotation period and total water production. The red line shows the observation, and the black lines the different sublimation models. Model A with globally constant surface active fraction (1∕12), model A/patches with bestfit adjustment of patches, and model P with effective sublimation curve and bestfit adjustment of patches. In all cases, the rotation period agrees reasonably well with observations. The total water production rate drops for model A and model A/patches scenarios with , while observations indicate . 
Fig. 8
Influence of the different surface areas on the torque evolution for a uniform active surface. Shown are the extrema of the Fourier torque components for each surface patch (see Fig. 10 for the patch boundaries). Model P seeks linear combinations of patches that in sum match the extrema derived from the observation, indicated by the dashed lines. 
6 Implication for the surface composition
The surface active fraction of the bestfit model is shown in Fig. 9 and as planar map in Fig. 10. The maps show the active fraction relative to the mean active fraction to highlight the differences to a uniformly active surface. The absolute value of the surface activity depends on the precise values of the cometary mass and the sublimation curve, while the relative distribution is not strongly affected. Patches with increased active water fraction are located in the southern hemisphere, which agrees with the increased southern relative activity shown in Fig. 6 by Fougere et al. (2016a) derived from Rosetta ROSINA/COPS/DFMS in situ gas densities. Toward the northern hemisphere, the zonal averaged active water fraction (side panel in Fig. 10) increases from the equatorial region, but stagnates on a lower value compared to the southern regions. Fougere et al. (2016a) retrieved a circular feature of increased relative activity around the north pole, which is not contained in our reconstruction. The direct use of measured gas densities from the ROSINA instruments to constrain the diurnal activity and the rotation state is limited because for operational reasons, Rosetta predominantly sampled gas in terminator illumination. Attree et al. (2019) discussed a fitting method for active fractions using three scalar properties: the Earthcomet distance, the total production rate, and the rotation period. They inferred a strong increase in active fraction on the southern hemisphere, limited to the time period of about 100 days around perihelion. In contrast, we worked with different observational data (all vector components of the rotation axis and the production rate) and derived an activity map that is constant in time.
Overall, the standard deviation from the homogeneous active surface (mean value 1.10) is 0.28, with the smallest activity on patch 21 (six times reduced active surface fraction). This confirms that the entire surface of 67P shows activity when it is insolated. A detailed correlation of our 36 patches with all geological regions cannot be expected because the resolution is just not good enough considering that the number of defined regions is now about twice as large; see Thomas et al. (2018). In Fig. 11 the surface activefraction regions fromFig. 10 have been draped onto the shape model of 67P (shown as color overlay). Thirty OSIRIS NAC images were mapped onto the shape to provide the morphological context. The images were acquired during the SHAP4S and SHAP5 mission phases for the northern hemisphere (September to October 2014), and SHAP7 and SHAP8 for the southern hemisphere (April to June 2016). Some image boundaries are visible in the mosaic because of the varying illumination condition during the mission phases. In general, the surface active fraction shows a northsouth trend, with the highest active fraction being in the roughconsolidated terrain of the southoriented regions (in particular around the southern neck regions, Fig. 11 − Z). The northern dustcovered regions, such as the Seth and Hapi region in the northern neck (Fig. 11 + Z), show intermediate levels of active fraction. This is compatible with the northern neck region being the most active in dust production during the early parts of the Rosetta mission. Some other features are seen: the activefraction map shows a dichotomy between the northern neck region of the big lobe (Seth) and the northern foot regions of the big lobe (see Fig. 12c). This dichotomy is not reflected in the surface morphology. Both sides of the big lobe show the same kind of smooth dustcovered terrain. It does, however, make sense from an insolation point of view. The northern neck is in polar night during the perihelion passage, while the foot of the big lobe is permanently illuminated throughout the comet year (Keller et al. 2015b). The volatiles in the northern neck are being replenished by seasonal mass transport on the comet (Keller et al. 2017). Mass transport on the foot of the comet will tend to accumulate in the Imhotep region, which is a gravitational low point on the comet. The Imhotep region (Fig. 11 − X) does indeed show comparable levels of active fraction to those of the northern neck. The Khonsu region (Fig. 12a) shows a southeast to northwest gradient in active fraction. The Khonsu region is a depression with a very rough terrain. Khonsu may be the result of an earlier fragmentation event that has caused parts of the surface to break off the nucleus. There is no significant morphological difference between one end of Khonsu and the other, and the integrated insolation is comparable. This may be a modeling artifact caused by the nonrandom choice of patch boundaries. The Wosret region (Fig. 12b) shows a surprisingly low level of active fraction. The Wosret region is the main part of the polar circle that receives permanent diurnal illumination during the perihelion passage of the comet. The region therefore has the highest potential level of activity of any region on the comet (highest integrated insolation). The morphology of the region is, however, quite different from the other southoriented regions on the comet. Wosret has a highly smooth but consolidated terrain toward the top of the small lobe and a much rougher consolidated terrain toward the southern neck. The areas of the activefraction map with lowest values are correlated with the smooth consolidated terrain. The rougher parts show a significantly higher active fraction. These activefraction values are more compatible with the levels found in the southern neck, which has a comparable terrain morphology. A possible explanation is that the smooth consolidated terrain is simply depleted of volatiles and will therefore exhibit no activity, regardless of the insolation. The smooth consolidated Wosret region could represents the final state terrain of cometary evolution.
Fig. 9
Surface active fraction f_{i} (Eq. (23)) relative to 1∕6 determined from the torque fit using the effective sublimation curve from Fig. 1a, dashed line, with the 36 patches shown in Fig. 10. The dashed lines indicate the mean value and the standard deviation. 
Fig. 10
Surface map showing the surface active fraction f_{i} (Eq. (23)) relative to 1∕6 corresponding to Fig. 9. The numbers indicate the patch label. Side panel: zonal averaged active water fraction. 
Fig. 11
Surface active fraction map projected onto the DLR SHAP7 shape model (Preusker et al. 2017). The shape has been textured using 30 OSIRIS NAC images acquired during the SHAP4S and SHAP5 mission phases for the norther hemisphere and the SHAP7 and SHAP8 mission phases for the southern hemisphere. The color overlay shows the active surface fraction from Fig. 10 with the view vector indicated by the basis vectors X, Y, Z in the body frame. 
Fig. 12
Zoom into the shape shown in Fig. 11. Part A: Khonsu region on the big lobe of the comet nucleus. Part B: Wosret region of the small lobe of the comet nucleus. Part C: northern edge of the big lobe. The color intensity of the Wosret (panel B) view has been decreased as compared to Fig. 11 to allow better visibility of the background image data. The wedgelike feature from the left side is an image artifact caused by an image acquired with a high (~ 90°) Sun incidence angle. 
7 Conclusions
We have presented a method to parameterize the observed rotationaxis movement in terms of a theory of Fourier coefficients. The sublimationinduced torques are encoded in three physically relevant combinations of the Fourier coefficients, which steer the rotation period changes and the rotationaxis movement. In particular, the rotation state of 67P is determined from the orbital evolution of the subsolar longitude and the specific shape. The increase in rotation period is caused by the diurnal average of the rotationaxisaligned torque (Fourier coefficient C_{III} = C_{0,z}), while the orientation change is caused by the diurnal torque cycle of the perpendicular components (Fourier coefficients C_{I} and C_{II}). Only by taking all three Fourier components together does a consistent fit result that constrains the local surface active fraction. From our analysis we conclude the following:

The sublimation model P contains a best fit for the surface active fraction to the observed rotation state, namely period and axis orientation.

The model includes a sublimation curve that increases much faster than linearly with insolation and reproduces the water production of 67P in Hansen et al. (2016) and Läuter et al. (2019).

A relatively small local variability (standard deviation 0.28) of the active surface fraction yields the required changes in rotation state.

Some area around Wosret on the small lobe seems to be less active, while the southern latitudes <−60° show an increased surface active fraction.
A further argument for a mostly uniform gas release comes from the observation of the dust structures in the inner coma modeled by Kramer & Noack (2016) and Kramer et al. (2018). The developed Fourier theory could be applied to other solar system bodies for which accurate measurements of the rotation axis motion and the shape are available.
Acknowledgements
The authors acknowledge the NorthGerman Supercomputing Alliance (HLRN) for providing computing time on the Cray XC40. This project has received funding from the European Unions Horizon 2020 research and innovation program under grant agreement No 686709 (MiARD).
References
 Attree, N., Jorda, L., Groussin, O., et al. 2019, A&A, 630, A18 (Rosetta 2 SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brouet, Y., LevasseurRegourd, A. C., Sabouroux, P., et al. 2016, MNRAS, 462, S89 [CrossRef] [Google Scholar]
 Davidsson, B., & Gutierrez, P. 2005, Icarus, 176, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Filacchione, G., Raponi, A., Capaccioni, F., et al. 2016, Science, 354, 1563 [NASA ADS] [CrossRef] [Google Scholar]
 Fougere, N., Altwegg, K., Berthelier, J. J., et al. 2016a, MNRAS, 462, S156 [CrossRef] [Google Scholar]
 Fougere, N., Altwegg, K., Berthelier, J.J., et al. 2016b, A&A, 588, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaskell, R. W., BarnouinJha, O. S., Scheeres, D. J., et al. 2008, Meteorit. Planet. Sci., 43, 1049 [Google Scholar]
 Godard, B., Budnik, F., Muñoz, P., Morley, T., & Janarthanan, V. 2015, Proceedings 25th International Symposium on Space Flight Dynamics–25th ISSFD, Munich, Germany [Google Scholar]
 Godard, B., Budnik, F., Bellei, G., & Morley, T. 2017, in International Symposium on Space Flight Dynamics26th ISSFD [Google Scholar]
 Gutiérrez, P. J., Jorda, L., Ortiz, J. L., & Rodrigo, R. 2003, A&A, 406, 1123 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Jewitt, D. 1997, Earth, Moon, Planets, 79, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Jorda, L., & Gutiérrez, P. 2002, in Cometary Science after HaleBopp, eds. H. Boehnhardt, M. Combi, M. R. Kidger, & R. Schulz (Dordrecht: Springer), 135 [Google Scholar]
 Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, H. U., Barbieri, C., Lamy, P., et al. 2007, Space Sci. Rev., 128, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, H. U., Mottola, S., Skorov, Y., & Jorda, L. 2015a, A&A, 579, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keller, H. U., Mottola, S., Davidsson, B., et al. 2015b, A&A, 583, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keller, H. U., Mottola, S., Hviid, S. F., et al. 2017, MNRAS, 371, 357 [Google Scholar]
 Knapmeyer, M., Fischer, H.H., Knollenberg, J., et al. 2018, Icarus, 310, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, T., & Noack, M. 2016, ApJ, 823, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, T., Läuter, M., Rubin, M., & Altwegg, K. 2017, MNRAS, 469, S20 [CrossRef] [Google Scholar]
 Kramer, T., Noack, M., Baum, D., Hege, H.C., & Heller, E. J. 2018, Adv. Phys. X, 3, 1404436 [Google Scholar]
 Lai, I.l., Ip, W.h., Su, C.c., et al. 2016, MNRAS, 462, S533 [CrossRef] [Google Scholar]
 Läuter, M., Kramer, T., Rubin, M., & Altwegg, K. 2019, MNRAS, 483, 852 [Google Scholar]
 Marsden, B. G., Sekanina, Z., & Yeomans, D. K. 1973, AJ, 78, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Marshall, D., Groussin, O., Vincent, J.B., et al. 2018, A&A, 616, A122 [NASA ADS] [CrossRef] [EDP Sciences] [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., & Samarasinha, N. H. 2018, AJ, 156, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Preusker, F., Scholten, F., Matz, K.D., et al. 2015, A&A, 583, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Preusker, F., Scholten, F., Matz, K.D., et al. 2017, A&A, 607, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Radhakrishnan, K., & Hindmarsh, A. C. 1993, Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations, NASA Reference Publication 1327, NASA [Google Scholar]
 Samarasinha, N. H., Mueller, B. E. A., Belton, M. J. S., & Jorda, L. 2004, in Comets II (Tucson, AZ: University of Arizona Press), 281 [Google Scholar]
 Shi, X., Hu, X., Sierks, H., et al. 2016, A&A, 586, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shoemake, K. 1985, ACM SIGGRAPH Computer Graphics, 19, 245 [CrossRef] [Google Scholar]
 Sosa, A., & Fernández, J. A. 2009, MNRAS, 393, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, N., El Maarry, M., Theologou, P., et al. 2018, Planet. Space Sci., 164, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Thomson, W. 1986, Introduction to Space Dynamics (New York, NY: Dover Publications) [Google Scholar]
 Valette, S., Chassery, J.M., & Prost, R. 2008, IEEE Trans. Vis. Comput. Graph., 14, 369 [CrossRef] [Google Scholar]
 Whipple, F. L. 1950, ApJ, 111, 375 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1
Radiationdriven sublimation (a) rate and (b) velocity as function of received irradiation. Solid blue lines show model A (surface active fraction = 1) and dashed red lines the effective sublimation curve with enhanced radiation response near perihelion. The effective sublimation curve leads to agreement with observations (see Figs. 6 and 7). 

In the text 
Fig. 2
Orientation changes in rotation axis. Black dots denote the reconstructed right ascension multiplied with cos 64° and declination, and the solid line shows a smooth fit to the data. 

In the text 
Fig. 3
Calculated sublimationinduced torque at perihelion for one rotation period at perihelion as a function of subsolar longitude from model P using the effective sublimation curve from Fig. 1 with the bestfit solution from Fig. 10. The dashed lines show the Fourier representation (constant, cos, and sin terms) of the corresponding solid lines to parameterize the diurnal cycle. 

In the text 
Fig. 4
Torque evolution of a uniform active surface model with the effective sublimation curve from Fig. 1. The time evolution of the first three Fourier coefficients , , and is shown for (−300 : 300) days from perihelion for each Cartesian component of the bodyframe torque, Eq. (25). Units of 10^{6} kg m^{2} s^{−2}. 

In the text 
Fig. 5
Time evolution of the three physical relevant combinations of Fourier coefficients (Eq. (35)) for (−300 :300) days from perihelion, units of 10^{6} kg m^{2} s^{−2}. (a) Model A, global uniform active surface 1∕12. (b) Model P (patches with effective sublimation curve). (c) Observed Fourier coefficients inferred from the rotationaxis movement and the tensor of inertia. 

In the text 
Fig. 6
Rotationaxis movement. The red line shows the observation from Fig. 2, and the other linesrepresent different sublimation models. Model A with a globally constant surface active fraction (1∕12), model A/patches with bestfit adjustment of patches, and model P with effective sublimation curve and bestfit adjustment of patches. The gray inset shows a magnification of the curves. 

In the text 
Fig. 7
Rotation period and total water production. The red line shows the observation, and the black lines the different sublimation models. Model A with globally constant surface active fraction (1∕12), model A/patches with bestfit adjustment of patches, and model P with effective sublimation curve and bestfit adjustment of patches. In all cases, the rotation period agrees reasonably well with observations. The total water production rate drops for model A and model A/patches scenarios with , while observations indicate . 

In the text 
Fig. 8
Influence of the different surface areas on the torque evolution for a uniform active surface. Shown are the extrema of the Fourier torque components for each surface patch (see Fig. 10 for the patch boundaries). Model P seeks linear combinations of patches that in sum match the extrema derived from the observation, indicated by the dashed lines. 

In the text 
Fig. 9
Surface active fraction f_{i} (Eq. (23)) relative to 1∕6 determined from the torque fit using the effective sublimation curve from Fig. 1a, dashed line, with the 36 patches shown in Fig. 10. The dashed lines indicate the mean value and the standard deviation. 

In the text 
Fig. 10
Surface map showing the surface active fraction f_{i} (Eq. (23)) relative to 1∕6 corresponding to Fig. 9. The numbers indicate the patch label. Side panel: zonal averaged active water fraction. 

In the text 
Fig. 11
Surface active fraction map projected onto the DLR SHAP7 shape model (Preusker et al. 2017). The shape has been textured using 30 OSIRIS NAC images acquired during the SHAP4S and SHAP5 mission phases for the norther hemisphere and the SHAP7 and SHAP8 mission phases for the southern hemisphere. The color overlay shows the active surface fraction from Fig. 10 with the view vector indicated by the basis vectors X, Y, Z in the body frame. 

In the text 
Fig. 12
Zoom into the shape shown in Fig. 11. Part A: Khonsu region on the big lobe of the comet nucleus. Part B: Wosret region of the small lobe of the comet nucleus. Part C: northern edge of the big lobe. The color intensity of the Wosret (panel B) view has been decreased as compared to Fig. 11 to allow better visibility of the background image data. The wedgelike feature from the left side is an image artifact caused by an image acquired with a high (~ 90°) Sun incidence angle. 

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.