Comet 67P/Churyumov-Gerasimenko rotation changes derived from sublimation induced torques

TL;DR

This paper models how sublimation-driven torques from comet 67P/Churyumov-Gerasimenko cause observable changes in its rotation period and axis orientation, providing insights into surface activity and physical properties.

Contribution

It introduces an analytic model linking sublimation-induced torques to rotation changes, estimating surface activity distribution and physical properties of the comet.

Findings

Rotation changes constrain surface activity and sublimation rates.

Surface active fraction distribution aligns with observed rotation data.

Sublimation activity increases faster than insolation near perihelion.

Abstract

(Context) The change of the rotation period and the orientation of the rotation axis of comet 67P/Churyumov-Gerasimenko (67P/C-G) is deducible from images taken by the scientific imaging instruments on-board the Rosetta mission with high precision. Non gravitational forces are a natural explanation for these data. (Aims) We describe observed changes for the orientation of the rotation axis and the rotation period of 67P/C-G. For these changes we give an explanation 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 solve the equation of state for the angular momentum in the inertial and the body-fixed frames and provide an analytic theory of the rotation changes in terms of Fourier coefficients, generally applicable to periodically forced rigid body dynamics. (Results) The torque induced changes of the 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/C-G. The torque movement confirms that the sublimation increases faster than the insolation towards perihelion. The derived relatively uniform activity pattern is discussed in terms of related surface features.

Figures

Equations

$$\text{@vec}{x}^{\prime }(\text{@vec}x,t)=\text{@vec}r(t)+\text{@tens}R(t)\text{@vec}x,$$

$$\dot{\text{@tens}R}(t)\text{@vec}x=\text{@vec}\omega (t)\times \text{@tens}R(t)\text{@vec}x,$$

$$M(t)$$

$$\overline{\text{@vec}x}(t)$$

$$\overline{\text{@vec}{x}^{\prime }}(t)=\frac{1}{M}\int \mathrm{d}\text{@vec}{x}^{\prime }\text{@vec}{x}^{\prime }\rho (\text{@vec}{x}^{\prime },t)=\text{@vec}r+R\overline{\text{@vec}x}.$$

$$\text{@vec}V(\text{@vec}x,t)={\partial }_t\text{@vec}{x}^{\prime }=\dot{\text{@vec}r}+\text{@vec}\omega \times \text{@tens}R\text{@vec}x.$$

$$\text{@vec}P(t)$$

$$\text{@vec}L(t)$$

$$\text{@vec}{F}_{\mathrm{NG}}(t)$$

$$\text{@vec}{T}_{NG}(t)$$

$$\dot{m}(\text{@vec}x,t)=\sum _{\mathrm{gas}}{f}_{\mathrm{gas}}(\text{@vec}x)Z_{\mathrm{gas}}(\text{@vec}x,t),$$

$$\text{@vec}{F}_{\mathrm{G}}(t)$$

$$\text{@vec}{T}_{\mathrm{G}}(t)$$

$$\dot{\text{@vec}P}=\text{@vec}{F}_{\mathrm{NG}}+M\text{@vec}a,\dot{\text{@tens}R}()=\text{@vec}\omega \times \text{@tens}R(),\dot{\text{@vec}L}=\text{@vec}{T}_{\mathrm{NG}}+M\text{@tens}R\overline{\text{@vec}x}\times \text{@vec}a,$$

$$M\dot{\text{@vec}r}+M\text{@vec}\omega \times \text{@tens}R\overline{\text{@vec}x}=\text{@vec}P,\text{@tens}I\text{@vec}\omega -M\dot{\text{@vec}r}\times \text{@tens}R\overline{\text{@vec}x}=\text{@vec}L$$

$$\dot{\text{@vec}P}=\text{@vec}{F}_{\mathrm{NG}}+M\text{@vec}a,M\dot{\text{@vec}r}=\text{@vec}P$$

$$\dot{\text{@tens}R}()=\text{@vec}\omega \times \text{@tens}R(),\dot{\text{@vec}L}=\text{@vec}{T}_{\mathrm{NG}},\text{@tens}I\text{@vec}\omega =\text{@vec}L$$

$$\displaystyle\@tens{I}_{\rm hom}=\left(\!\!\!\begin{array}[]{ccc}9.55529\times 10^{18}&\!\!\!1.73767\times 10^{16}&\!\!\!2.24462\times 10^{17}\\ 1.73767\times 10^{16}&\!\!\!1.76369\times 10^{19}&\!\!\!-7.45958\times 10^{16}\\ 2.24462\times 10^{17}&\!\!\!-7.45958\times 10^{16}&\!\!\!1.89825\times 10^{19}\\ \end{array}\!\!\!\right)$$

$$\text{kg m2}.$$

$$\displaystyle\@tens{I}_{\rm bf}=\left(\!\!\!\begin{array}[]{ccc}9.3408457\times 10^{18}&\!\!\!5.6695663\times 10^{16}&\!\!\!0\\ 5.6695663\times 10^{16}&\!\!\!1.6562414\times 10^{19}&\!\!\!0\\ 0&0&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!1.8192083\times 10^{19}\end{array}\!\!\!\right)$$

$$\text{kg m2}.$$

$$(1-A)f_{\text{sun}}=ϵ\sigma T^4+Z(T)L_{\mathrm{ice}}$$

$$Z_{\mathrm{Hertz}-\mathrm{Knudsen}}(T)=2P(T)/(\pi v_{\mathrm{th}}).$$

$$v_{\mathrm{th}}(T)=\sqrt{8RT/(\pi {\mu }_{H_{2}O})}$$

$$\text{@vec}{T}_{\mathrm{bf}}({\lambda }_s,{\varphi }_s,r_h)=\sum _i^{N_{\mathrm{faces}}}{f}_i{A}_i{Z}_{H_{2}O,i}{v}_{\mathrm{th},i}\text{@vec}{r}_i\times {\widehat{\text{@vec}n}}_i.$$

$$\text{@vec}{T}_{F,\mathrm{bf}}^{(i)}(\alpha )=\text{@vec}{C}_0^{(i)}+\sum _{n=1}^m\text{@vec}{C}_n^{(i)}\sin (n\alpha )+\sum _{n=1}^m\text{@vec}{C}_{m+n}^{(i)}\cos (n\alpha )$$

$$\text{@vec}{C}_n^{\mathrm{int}}(t)=\text{interpolation}(\text{@vec}{C}_n^{(i)},\dots ,\text{@vec}{C}_n^{(N_{\mathrm{intervals}})}),$$

$$\text{@vec}{T}_{F,bf}(t,{\lambda }_s)$$

$$\{x(t),y(t),z(t)\}$$

$${\lambda }_s(t)$$

$$\text{@tens}R(t_S)=\left(\begin{array}{ccc}\frac{s_z{s}_x^2+s_y^2}{s_x^2+s_y^2}& \frac{s_x{s}_y(s_z-1)}{s_x^2+s_y^2}& s_x\\ \frac{s_x{s}_y(s_z-1)}{s_x^2+s_y^2}& \frac{s_x^2+s_y^2{s}_z}{s_x^2+s_y^2}& s_y\\ -s_x& -s_y& s_z\end{array}\right).$$

$$\text{@vec}{\omega }_{\mathrm{bf}}(t_S)$$

$$\text{@vec}L(t_S)$$

$$\{x_s^0(t),y_s^0(t),z_s^0(t)\}$$

$${\lambda }_s^0(t)$$

Full text

*Aims. *We describe observed changes for the orientation of the rotation axis and the rotation period of 67P/C-G. For these changes we give an explanation 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 solve the equation of state for the angular momentum in the inertial and the body-fixed frames and provide an analytic theory of the rotation changes in terms of Fourier coefficients, generally applicable to periodically forced rigid body dynamics.

*Results. *The torque induced changes of the 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/C-G. The torque movement confirms that the sublimation increases faster than the insolation towards perihelion. The derived relatively uniform activity pattern is discussed in terms of related surface features.

1 Introduction

Rosetta’s instruments have probed the gas and dust environment during almost the entire apparition of 67P/Churyumov-Gerasimenko (67P/C-G) from 2014-2016. The observed changes of the 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/C-G shortened by 21 min, the same number has been reported for the previous apparition by Mottola et al. (2014). Keller et al. (2015b) proposed a sublimation driven model and explain the changing rotation period in terms of a homogeneous ice model of the entire surface (Keller et al. (2015a)) without studying the axis orientation. For many comets, changes in the 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/C-G was known in detail, Gutiérrez et al. (2003) explored several scenarios for rotation-axis changes of small, irregularly shaped comets. Typical changes of the rotation axis caused by sublimation forces range from $0.1$ to several tens of degrees. For 67P/C-G the observed change is on the lower end of this range (0.5∘).

For the model A proposed by Keller et al. (2015b), our analysis predicts a five times larger change of the direction of the rotation axis, also along a different direction compared to the observations. In addition, model A predicts a slower than observed increase of the total gas production with decreasing heliocentric distance, while the analysis of the coma by Hansen et al. (2016); Kramer et al. (2017); Läuter et al. (2018) shows a faster increase of the activity based on the pressure sensors. To overcome these discrepancies, we establish a formalism to match models and observations in terms of a Fourier analysis of the gas induced torque and derive a possible ice distribution on the surface which explains the rotation-period changes, the movement of the direction of the angular momentum, and the increase of activity with heliocentric distance (model P).

2 Forced rigid body dynamics

We review the response of the rotation state of the nucleus to sublimation and other processes which 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 three-dimensional body frame by a prescribed body-frame density $\rho_{\rm bf}(\@vec{x},t)$ for each $\@vec{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$

[TABLE]

for each body-frame point $\@vec{x}$. Here, $\@vec{r}(t)$ denotes the center of figure, $\@tens{R}(t)$ the orthogonal rotation matrix with the property

[TABLE]

and $\@vec{\omega}(t)$ the angular velocity. The mapping $\@vec{x}^{\prime}(\@vec{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 time-dependent density $\rho_{\rm bf}(\@vec{x},t)$ allows one to incorporate slow changes in the density and porosity of the comet. The density in the inertial frame $\rho(\@vec{x}^{\prime},t)=\rho_{\rm bf}(\@vec{x},t)$ is obtained from the body-frame 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 $\bar{\@vec{x}}$ are in general time-dependent

[TABLE]

For the center of mass in the inertial frame the relation holds

[TABLE]

The time derivative of Eq. (1) yields the inertial-frame velocity $\@vec{V}$ of a fixed body-point

[TABLE]

To obtain the linear momentum $\@vec{P}$ and angular momentum $\@vec{L}$ of the whole nucleus in the inertial system we integrate

[TABLE]

with the tensor of inertia $\@tens{I}(t)=\@tens{R}\@tens{I}_{\rm bf}(t)\@tens{R}^{-1}$ in the inertial frame with respect to the center of figure $\@vec{r}$ and the tensor of inertia $\@tens{I}_{\rm bf}(t)$ with respect to the body-frame center [math]. For the case of a time-dependent body density, $\@tens{I}_{\rm bf}(t)$ needs to be computed for the body-frame density at time $t$. The momentum changes are generated by the sum of non-gravitational and gravitational forces $\@vec{F}_{\rm NG}+\@vec{F}_{\rm G}=\dot{\@vec{P}}$ and torques $\@vec{T}_{\rm NG}+\@vec{T}_{\rm G}=\dot{\@vec{L}}$. Gas sublimation at point $\@vec{x}$ on the surface leads to a mass loss $\dot{m}$ and generates the non-gravitational components

[TABLE]

The gas velocity in the inertial frame $\@vec{v}_{\rm is}=u_{\rm gas}\hat{\@vec{n}}_{0}+\@vec{\omega}\times(\@vec{x}^{\prime}-\@vec{r})$ consists of two components, one into normal direction $\hat{\@vec{n}}_{0}(\@vec{x}^{\prime},t)$ on the nucleus’ surface in the inertial frame and one due to the body rotation, $u_{\rm gas}$ denotes the thermal gas velocity from Eq. (28). The gas velocity in the body frame is given by $\@vec{v}_{\rm bf}=\@tens{R}^{-1}\@vec{v}_{\rm is}=u_{\rm gas}\hat{\@vec{n}}+(\@tens{R}^{-1}\@vec{\omega})\times\@vec{x}$ with the outward surface normal $\hat{\@vec{n}}(\@vec{x},t)$ at surface location $\@vec{x}$ within the body frame. According to Jorda & Gutiérrez (2002), Eq. (11), the mass production $\dot{m}$ for a mixture of

gas species reads

[TABLE]

with the surface active fraction $f_{\rm gas}$ and the sublimation rate $Z_{\rm gas}$. The mass loss $\dot{m}$ changes the shape, reduces the total mass, and affects the tensor of inertia of the nucleus. The integrated mass loss of comet 67P/C-G during the 2015 apparition is estimated to be about 1/1000 of the total mass (see Godard et al. (2015) and Godard et al. (2017)) and therefore we neglect both effects, we assume time-independent mass and tensor of inertia.

The gravitational components for force and torque yield

[TABLE]

with the gravitational acceleration $\@vec{a}$ due to other solar system bodies. For both volume integrals in Eqs. (12), (13), $\@vec{a}(\@vec{x}^{\prime})$ is assumed to be constant over the cometary body, which implies that e.g. tidal forces are neglected.

Eqs. (2), (7), (8), (9), (10), (12), (13) result in a system of coupled algebraic and differential equations

[TABLE]

for the state variables $\@vec{r}(t)$, $\@vec{P}(t)$, $\@tens{R}(t)$, $\@vec{L}(t)$. Eqs. (7), (8) couple linear and rotational momenta through the center of mass in the body frame $\bar{\@vec{x}}(t)$. If the density distribution $\rho_{\rm bf}$ satisfies $\bar{\@vec{x}}=0$, $\@vec{r}=\overline{\@vec{x}^{\prime}}$ becomes the center of mass in the inertial frame and Eqs. (14) decouple into two blocks. The first block

[TABLE]

describes the translational movement for the state variables $\@vec{r}(t)$, $\@vec{P}(t)$, and the second one

[TABLE]

the rotational dynamics for $\@tens{R}(t)$ and $\@vec{L}(t)$. The model for the changing rotation period by Keller et al. (2015b) is contained as special case in Eq. (17). For that let us denote the eigenvector $\@vec{e}$ of $I_{\rm bf}$ in the body frame with the largest moment of inertia $I_{\rm z}$ and assume $\@vec{L}$ initially aligned with $\@vec{e}^{\prime}=\@tens{R}\@vec{e}$, that is $\@vec{L}(0)=L_{\rm z}\@vec{e}^{\prime}(0)$. Then $\@vec{\omega}(0)=L_{\rm z}/I_{\rm z}\@vec{e}^{\prime}(0)$ and consequently $\dot{\@tens{R}}\@vec{e}=0$. Thus, $\@vec{e}^{\prime}=\@tens{R}\@vec{e}$ is constant in time and the angular velocity changes with $\dot{\@vec{\omega}}\cdot\@vec{e}^{\prime}=\@vec{T}_{\rm NG}\cdot\@vec{e}^{\prime}/I_{\rm z}$.

For given initial conditions of all state variables, the system of Eqs. (14), (16), and (17) can be solved numerically. For accuracy, we use the LSODE package provided within Mathematica/FORTRAN (Radhakrishnan & Hindmarsh (1993)). Following Shoemake (1985), the matrix-matrix operation $\@tens{R}\rightarrow\@vec{\omega}\times\@tens{R}$ is replaced by a quaternion multiplication for improved stability. There is a one-to-one mapping between $\@tens{R}$ and a quaternion $\@vec{q}$ such that the matrix matrix operation is substituted by $\@vec{q}\rightarrow\@vec{q}\cdot(0,\@vec{\omega})/2$.

3 Rotation state of 67P/C-G

After the arrival of Rosetta at 67P/C-G in 2014, observations of the rotation by Preusker et al. (2015) and Godard et al. (2017) show 67P/C-G in an excited state of rotation, albeit with the rotation axis close to the axis $\@vec{e}^{\prime}$ with the largest moment of inertia (rotation state with minimum energy). This points to an alignment of both axes which is also compatible with observations (Jorda et al. (2016)). The total mass is estimated to be $10^{13}$ kg by Godard et al. (2015). The assumption of a strictly homogeneous density leads to a tensor of inertia

[TABLE]

with respect to the center of mass. This tensor is incompatible with the observations, since then the axis $\@vec{e}^{\prime}$ would be tilted by $2.9^{\circ}$ with respect to the rotation axis (Preusker et al. (2017)), and in addition an offset of the center of mass $\overline{\@vec{x}^{\prime}}$ with respect to the center of figure $\@vec{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 time-independent, non-homogeneous density distribution which aligns the rotation axis with the $\@vec{e}^{\prime}$ under the constraint that the total mass is kept fixed at $M=10^{13}$ kg and $\overline{\@vec{x}^{\prime}}=\@vec{r}$. For definiteness we give a possible mass distribution with resulting tensor of inertia

[TABLE]

This solution equals putting $(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 off-diagonal entries in Eq. (25) align the shape file with the body frame by a rotation of $0.4^{\circ}$. The changes of the rotation state of the comet is determined by the relative change of the 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 the 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); Läuter et al. (2018) from ROSINA COPS and DFMS data to be about $6.2\pm 2\times 10^{9}$ kg, corresponding to about $1/1600$ of the total mass of 67P/C-G ($M=10^{13}$ kg). The water production shows a steep increase with heliocentric distance $r_{h}^{\alpha}$ around Southern solstice, with exponents $\alpha$ ranging from $-6.5$ up to $-7$. The total gas production from the radiation driven sublimation model A by Keller et al. (2015b) yields smaller exponents $\alpha=-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 decimeter sized 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. (2018)) is not considered separately since the $CO_{2}$ sources around perihelion coincide with the water regions (Fougere et al. (2016b), Läuter et al. (2018)) 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 periodic torque components. At heliocentric distances larger than $3$ a.u., and in particular on the outbound orbit, $CO_{2}$ becomes the dominant species (Läuter et al. (2018)) and does not follow the subsolar illumination. At these distances the rotation period of the comet has settled and these times are outside the scope of the present analysis. The rotation axis of 67P/C-G shows largest movements around $\pm 100$ days from perihelion, in agreement with the larger $|\alpha|$ exponents derived from the gas instruments. Keller et al. (2015b) (model A) consider the gas production based on a shape model. On each surface element with a given Bond bolometric albedo $A$ and solar irradiance $f_{\sun}=S_{\sun}/r_{h}^{2}$ at heliocentric distance $r_{h}$ [a.u], solar constant $S_{\sun}=1361$ W m*-2*,

the energy balance

[TABLE]

is solved for the sublimation rate $Z$, given by the Hertz-Knudsen relation

[TABLE]

The parameters are taken from Keller et al. (2015b), with emissivity $\epsilon=0.9$, latent heat of sublimation for water ice $L_{\rm ice}=2.6\times 10^{6}$ J kg*-1* (assumed to be constant), water vapor pressure $P(T)=3.56\times 10^{12}\,{\rm e}^{-6141.667/T}$ [kg m*-1* s*-2*], and thermal velocity of water molecules with molar mass $\mu_{H_{2}O}$

[TABLE]

The gas constant is denoted by $R$, the Stefan-Boltzmann constant by $\sigma$. The solution to Eqs. (26,27) in terms of the sublimation rate $\dot{N}\equiv Z_{\rm Hertz-Knudsen}/m_{H_{2}O}$ [s*-1* m*-2*], mass of water molecule $m_{H_{2}O}$ [kg], is shown in Fig. 1 for $A=0.01$ (Keller et al. (2015b), model A). The observed changes by Hansen et al. (2016); Kramer et al. (2017); Läuter et al. (2018) point to a faster increase of the total production with decreasing perihelion distance.

To account for the observation requires to assume a sublimation rate which increases faster than linearly with illumination, as exemplified by the dashed red line in Fig. 1.

A possible physical mechanism behind the increase of the sublimation could be a decrease of the dust layer when the comet approaches the sun, leading to a steeper slope (model C to model A transition by Keller et al. (2015b)). 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/C-G by (Fougere et al. (2016a), Eq. (1)).

3.2 Observed rotation axis changes

The (RA,Dec) orientation of the angular velocity $\@vec{\omega}$ of the comet nucleus is derived using the set of about 25000 control points defined as the center of the maplets created in the stereophotoclinometry (SPC) method of Gaskell et al. (2008) applied to comet 67P/C-G 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 star tracker pointing measurements are 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 $\@vec{\omega}$ vector. They led us to strongly smooth the data as illustrated by the solid line in Fig. 2, which represents the time-averaged 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 third order polynomial. A similar data set has been retrieved by Godard et al. (2017).

4 Fourier theory of torque induced motion

The shape model from Preusker et al. (2017) is re-meshed using the Aproximated Centroidal Voronoi Diagrams (ACVD) tool by Valette et al. (2008) into $N_{\rm faces}=3996$ triangular elements with area $A_{i}$ and surface normal $\hat{\@vec{n}}_{i}$, chosen to be approximately of equal area. The results are robust against variations of the shape model as long as $N_{\rm faces}\gg 100$. The torque is evaluated according to Eqs. (10), (11). The sublimation rate $Z_{\rm gas}$ on each facet is evaluated for a given subsolar latitude $\phi_{s}$ and heliocentric distance $r_{h}$ during one rotation period, starting with subsolar longitude $\lambda_{s}=0$.

We compute the total torque arising from the water sublimation curves, either the Hertz-Knudsen rate from Eq. (27) or alternatively the effective sublimation curve in Fig. 1:

[TABLE]

The contributions of shadowed (including self-shadowed) faces are set to zero and all the surface active fraction are first set to the same value. In the next iteration, the surface active fraction is 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)$) is 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 is set to be in the lowest energy state (rotation axis and angular momentum aligned, pointing to RA $\alpha=69.427^{\circ}$, Dec $\delta=64.000^{\circ}$ at $t=-350$ days). For 67P/C-G the observed axis changes are small and we tabulate the subsolar latitude $\phi_{s}$ and heliocentric distance $|r_{h}|$ in $N_{\rm intervals}=81$ 10-days intervals and store the body-frame torque as function of $\lambda_{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 expand the torque components in a Fourier cosine/sine series. The periodic argument $\alpha=[-\pi:\pi[$ of the Fourier series is not time, but the subsolar longitude $\lambda_{s}=\alpha+\pi$ to accommodate changes in the rotation period. In the $i$th time interval we extract the first $N=2m+1$ Fourier coefficients $\@vec{C}_{n}^{(i)}=\{C^{(i)}_{n,x},C^{(i)}_{n,y},C^{(i)}_{n,z}\}$, which yield the Fourier series representation of the torque

[TABLE]

The final parametrization of the body torque as function of subsolar longitude along the entire orbit is given by the time evolution of the Fourier coefficients

[TABLE]

and rewriting

[TABLE]

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 non-diurnal torque coefficients. Since the rotation axis is aligned with the $z$ body-axis, the torque component $C^{\rm int}_{0,z}(t)$ directly affects the rotation period. The orientation changes of the rotation axis are caused by the diurnal $C^{\rm int}_{1,x},C^{\rm int}_{2,x},C^{\rm int}_{1,y},C^{\rm int}_{2,y}$ 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 neglect the changes of the orbital elements due to non-gravitational momentum (Eq. 16) and take the orbital evolution $\@vec{r}_{h}(t)$ as fixed. For a given rotation matrix $\@tens{R}(t)$ (body-frame to equatorial inertial-system) and position of the comet $\@vec{r}_{h}(t)$, the sub-solar longitude is given by

[TABLE]

At initial time $t_{S}$ the rotation matrix $\@tens{R}(t_{S})$, which transfers the body-frame $z$-axis to the rotation axis $\hat{\@vec{s}}=\{s_{x},s_{y},s_{z}\}$ in the equatorial inertial frame, is given by

[TABLE]

Eq. (17) is then integrated with the initial angular velocity and momentum set to

[TABLE]

The Fourier decomposition provides additional insights into the axis changes. An important parameter is the angle $\lambda_{s}^{0}$ around the $z$-axis to point the $x$-$z$-plane towards towards the sun (zero sub-solar latitude). Neglecting the $0.5^{\circ}$ tilt-change of the rotation axis, $\lambda_{s}^{0}$ is given by

[TABLE]

We obtain a good approximation of the angular momentum change $\boldsymbol{\Delta}\tilde{\@vec{L}}$ during one rotation period $T_{\rm rot}=2\pi/\omega$ by keeping $\lambda_{s}^{0}(t)$ fixed during this rotation and integrating the torque in the body frame $\@vec{T}_{\rm bf}$, parametrized by the subsolar longitude and the Fourier components from Eq. (32)

[TABLE]

The angular momentum change along the entire orbit is then apprximated by adding all $\@tens{R}(t_{S})\Delta\tilde{\@vec{L}}(t)$ contributions to the initial angular momentum. The “orbit” matrix does not affect the magnitude of the “shape” vector. All shown results are done without this approximation and use the full numerical solution of Eq. (17). Eq. (58) is used to determine the physically relevant Fourier components

[TABLE]

for analyzing the observations.

4.2 Extract observed torque from the rotation-axis evolution

Next, we consider the inverse problem of finding a plausible torque function in the cometary body frame as function of subsolar coordinates and solar distance. We infer 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. (25).

To parametrize the observed torque as function of subsolar longitude, we compute $\lambda_{s}(t)$ from Eq. (33) at each observation time. Every 10 days, we find the closest instance $t_{i}$ of $\lambda_{s}(t_{i})=0$ and compute the Fourier coefficients $C_{x,y,z}^{\rm obs}$ to represent $\@vec{T}_{\rm bf}(\lambda_{s}=0\ldots 2\pi)$. Only the three Fourier combinations from Eq. (4.1) components should be retrieved from the fit (Fig. 5), since the axis motion is not sensitive to the other Fourier components (see Eq. (58)).

5 Matching Fourier coefficients with the observed torque

The simplest sublimation model A results in a rotation axis movement shown in Fig. 6, green line. The evolution is rotated by 90∘ with respect to the observed torque movement (Fig. 6, red line) and leads to a largely increased axis tilt compared to observations. To explain the observations requires to consider a spatially heterogeneous surface with varying water-ice coverage (model P). We could show that an alternative explanation is a large thermal lag of several hours of the maximum sublimation with respect to the maximum irradiation caused by a dust layer of some millimeters thickness. However, this scenario is unlikely, since the response of sublimation rates to radiation is almost immediate as seen by short ($<1$ h) delays of jet outbreaks (Lai et al. (2016); Shi et al. (2016)) and the good agreement of inner dust structures with illumination driven dust release (Kramer & Noack (2016); Kramer et al. (2018)).

Measurements of VIRTIS and MIRO found that the thermal inertia is lower than $320$ JK*-1m-2s-0.5* when including the error bars (see Marshall et al. (2018) for an overview). Such a small thermal inertia is not able to provide the needed phase lag of several hours of the maximum sublimation with respect to the maximum irradiation as thermal simulations show and measurements of the activity maxima compared to noon time show Shi et al. (2016).

In fact we show that the often invoked unrealistic thermal lag to explain the non-gravitational forces acting on the cometary orbit can be at least partly replaced by the effects of a complex nucleus shape and its slightly non-uniform activity (Davidsson & Gutierrez (2005),Sosa & Fernández (2009)).

For the non-uniform case, we divide the surface in 36 equally spaced patches and compute 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}$, $C_{III}$ in Eq. (4.1) 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 single patch contribution to the observation is that for the 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 the $C_{II}$ and$C_{III}$ components. 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. Details of the data selection and the chosen norm influence the final fit result, but the general structure with the identified depleted and enhanced surface active regions remains unaffected.

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_{\rm tot}(r)\sim r^{-2.8}$. In contrast, observations from COPS/DFMS point to a larger exponent $\alpha\sim-6$ to $-7$. The change of 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 non-linearly with insolation, as discussed in Sect. 3.1. The effective sublimation curve in Fig. 1(a), dashed line yields the total production displayed in Fig. 7, with larger exponent $\alpha<-5$ as measured by several Rosetta instruments (see e.g.  Hansen et al. (2016); Kramer et al. (2017); Läuter et al. (2018)) and modeled by Hu et al. (2017). The rotation-axis motion of this modified sublimation model is shown in Fig. 6 and are in better agreement with observations than the other considered scenarios.

6 Implication for the surface composition

The surface active fraction of the best-fit 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 activity shown in Fig. 6 by Fougere et al. (2016b) derived from Rosetta ROSINA/COPS/DFMS in-situ gas densities. The direct use measured gas densities from the ROSINA instruments to constrain the diurnal activity and the rotation state is limited, since for operational reasons Rosetta predominantly sampled gas in terminator illumination.

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 all the surface of 67P/C-G shows activity whenever insolated.

A detailed correlation of our 36 patches with all geological regions cannot be expected since the resolution is just not good enough considering that the number of defined regions are now about twice as large Thomas et al. (2018). In Fig. 11 the surface active-fraction regions from Fig. 10 have been draped onto the shape model of 67P/C-G (shown as color overlay). 30 OSIRIS NAC images have been mapped onto the shape to provide the morphological context. The images have been 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 present during the mission phases. In general the surface active fraction shows a north south trend with the highest active fraction being in the rough consolidated terrain of the south oriented regions (In particular around the southern neck regions, Fig. 11 $-Z$). The northern dust covered regions like the Seth and Hapi region in the northern neck (Fig. 11 $+Z$) shows 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 active fraction map shows a dichotomy between 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 dust covered 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. (2015a)). 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 south-east to north-west 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 non random choice of patch boundaries. The Wosret region (Fig. 12b) shows a surprising low level of active fraction. The Wosret region is the major 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 south oriented regions on the comet. Wosret has a highly smooth but consolidated terrain towards the top of the small lobe and a much rougher consolidated terrain towards the southern neck. The areas of the active fraction map with lowest values are correlated with the smooth consolidated terrain. The rougher parts shows significantly higher active fraction. These active fraction values are more compatible with the levels found in the southern neck which has comparable terrain morphology. A possible explanation is that the smooth consolidated terrain is simply depleted of volatiles and will therefore exhibit no activity no matter the insolation. The smooth consolidated Wosret region could represents the final state terrain of cometary evolution.

7 Conclusions

We have presented a method to parametrize the observed rotation-axis movement in terms of a theory of Fourier coefficients. The sublimation induced torques are encoded in three physical relevant combinations of the Fourier coefficients, which steer the rotation period changes and the rotation axis movement. In particular, the rotation state of 67P/C-G is determined from the orbital evolution of the subsolar longitude and the specific shape. The increase of the rotation period is caused by the diurnal-average of the rotation axis aligned 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}$, $C_{II}$). Only by taking all three Fourier components together, a consistent fit results which constrains the local surface active fraction. From our analysis we conclude the following points:

• •

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/C-G in Hansen et al. (2016); Läuter et al. (2018).

• •

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

• •

Some area around Wosret on the small lobe seems to be less active, while the southern latitudes $<-60^{\circ}$ 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); 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 North-German Supercomputing Alliance (HLRN) for providing computing time on the Cray XC40.