© The Authors 2025

1 Introduction

The sublimation of ice from a cometary nucleus transfers momentum from the nucleus to the surrounding space. The spatial integration over the momentum field on the surface results in nongravitational acceleration (NGA) and nongravitational torque (NGT) acting on the rigid body of the comet. The NGA leads to variations in osculating elements and in trajectories through the solar system, while the NGT changes the rotational angular momentum, and thus the rotation state of the nucleus. The overall magnitudes of the NGA and the NGT depend on the effective cross section of the nucleus illuminated by the sun.

The significance of NGAs and in particular of NGTs for the change in the rotation state is pointed out by Whipple (1950). Nongravitational torques are investigated; for example, for comet 2P/Encke by Whipple & Sekanina (1979), for comet 22P/Kopff by Sekanina (1984), for comet 1P/Halley by Peale & Lissauer (1989), Julian (1990), and for Halley-like comets by Samarasinha & Belton (1995). For rotational dynamics of general comets, reviews are outlined by Samarasinha et al. (2004), Yeomans et al. (2004), and Knight et al. (2024).

A rotating jet model on a spherical body used by Chesley & Yeomans (2004) features NGAs but NGTs are always zero. The torque formation and changes in rotation state vary with the underlying geometry; for example, oblate or prolate ellipsoids lead to different torque effects. Sitarski (1995) controls the evolution of the spin axis with an additional parameter in a linear precession model. Whipple & Sekanina (1979) and subsequent authors, for example Sekanina (1981), Królikowska (1998), and Królikowska et al. (2001), study the evolution of the rotation axis for a nonspherical nucleus under the assumption of a forced precession model. Jorda & Gutiérrez (2002) and Gutiérrez et al. (2002) investigate the sublimation activity of irregularly shaped cometary nuclei. Kramer et al. (2019) develop a Fourier series-based approach to parametrize the observed torque variations of comet 67P/Churyumov-Gerasimenko (67P/C-G) in a body-fixed coordinate system and compared the observations with different sublimation models. Attree et al. (2019, 2023, 2024) discuss the evolution of the torque of comet 67P/C-G based on the geomorphological properties defined in Marschall et al. (2017). Davidsson & Gutiérrez (2004) approximate the surface of an ellipsoid with a triangular mesh and regionally weighted activity. Maquet et al. (2012) weight the sublimation activity based on latitudinal bands on an ellipsoid surface.

The rotational dynamics of a cometary nucleus can be considered as the rotation of a rigid body, as is pointed out by Samarasinha et al. (2004). Early approaches such as that of Whipple & Sekanina (1979) are restricted to simple rotations, which are principal axis rotations. These are characterized as rotations around only one of the principal axes of the tensor of inertia. In this situation the directions for the rotational angular momentum and for the angular velocity are the same as that of the principal axis. A special case of a principal axis rotation is the state of lowest rotational energy for a prescribed angular momentum. It is the rotation around the axis of the largest moment of inertia. In general, a nucleus can be in a nonprincipal axis rotation where the vectors for rotational angular momentum and for angular velocity do not point into the same direction. As is discussed by Julian (1987), Samarasinha & Belton (1995), Gutiérrez et al. (2003), and Knight et al. (2024), nonprincipal axis rotations can be in short axis mode or in long axis mode. The lowest-energy state is the limiting case of a short axis mode. As is described by Jewitt (1997) and Frouard & Efroimsky (2018), in the presence of any damping for kinetic energy, a cometary nucleus in an excited state (i.e., one that is above lowest energy) will always tend to reach a lowest-energy state. The two comets 1 P and 103 P are regarded as undergoing nonprincipal axis rotation, with the possibility that other comets may also exhibit this behavior, as is proposed by Knight et al. (2024). Conversely, the majority of comets do not exhibit indications of strongly excited rotation states. In the case of comet 67P/C-G, Gutiérrez et al. (2016), Godard et al. (2017), and Kramer et al. (2019) report only a minor excitation from the lowest-energy state. This and the fact that most comets do not show excited rotations (Knight et al. 2024) leads us to the assumption that the lowest-energy state reflects the properties of most known comets in an appropriate manner. In this case, observational data for angular velocity can be used to assess torque directions.

Momentum is created by sublimation processes close to the surface of the cometary nucleus. Incoming solar irradiation heats up the upper surface layers and is mostly reradiated into space. The sublimation process consumes some fraction of the incoming energy, while the remaining energy propagates into subsurface layers, where other frozen volatiles might be present. Thermophysical models (TPMs) describe these processes at a point on the surface of the nucleus, as is discussed by Keller et al. (2015a) and Blum (2018). Complex TPMs consider the composition of the material on the surface and within the column of the subsurface domain, the temporal evolution of solar irradiation, and the energy fluxes within the column. The analyses of the Rosetta images of comet 67P/C-G by Kramer & Noack (2016), Gundlach et al. (2020), and Fulle et al. (2020) point to very low values of heat conductance, and thus sublimation of water occurs almost instantaneously with respect to the incoming radiation. Due to a higher sublimation rate at lower temperatures, CO₂ emission can occur even at night. Here, we constrict the discussion to TPMs that establish sublimation flux as an instantaneous response to the incoming irradiation. These models balance reradiation due to scattering or thermal activity and sublimation energy. Heat transfer into the cometary material is neglected. An increase in sublimation is related to an increase in irradiation (monotonous function). In this case, no heat flow to the subsurface layers or any other delaying process is considered. An example of such an approach is model A in Keller et al. (2015b), which is based on the energy balance for incoming radiation and outgoing thermal and sublimation fluxes on the surface.

The shape of the nucleus is of particular importance for rotational dynamics. In addition to the tensor of inertia, the formation of torques is strongly affected by the local geometry of the surface faces with respect to the rotational axis. The basic geometry-driven effects can already be studied for simplified ellipsoidal shapes. More complex body shapes further complicate the analysis, for instance, by introducing concavities. Concavities change the irradiation by self-shadowing and lead to significant variability in the regional torque directions. The basic effect of concavities can be studied for a shape with a simplified bilobed structure. Shapes that approximate the nucleus of comet 67P/C-G bring along further complexity, including surface roughness and boulders, which are not resolved in the typical grid sizes used for torque computations. We study the dependence of the torque on the mesh size using remeshed versions of the shape model of Preusker et al. (2017) for comet 67P/C-G. Also for other comets shape models have been derived. The combination of light curve inversion for convex bodies of Kaasalainen (2001) and radar imaging of Ostro et al. (2002) yield coarse resolution for surface geometry such as in Donaldson et al. (2023) for 162P/Siding Spring and in Magri et al. (2011) for Asteroid (8567) 1996 HW1. From spacecraft imagery the shapes of comets 1P, 9P, 19P, 67P, 81P, and 103P are constrained (see Snodgrass et al. 2022). Knight et al. (2024) classify four of these comets as having a bilobed shape.

Observational data for the rotation state of comets are mainly extracted from light curves and radar detections. The retrieved rotation periods are in the range between a few hours and several days. The rotation period can change over time as is reported by Belton et al. (2011) and Chesley et al. (2013) for 9P/Tempel 1. Even the tendency might change from increase to decrease as is reported by Farnham et al. (2021) for comet 46P/Wirtanen during one perihelion passage. Earth-bound observations constrain the orientation of the rotation axis by studying coma morphology. Measurements are typically associated with uncertainties on the order of 15^∘, see Knight et al. (2021). For comet 67P/C-G, orientations of the rotation axis are extracted by Jorda et al. (2016) based on the method of Gaskell et al. (2008). Models for the rotation state (see Kramer et al. 2019; Attree et al. 2019, 2023, 2024) might match complementing quantities, such as the rotation period or the production curve for global gas production.

Gas emanating from the surface carries along momentum determined by the molecular mass and velocity of the volatiles. Thermophysical models with a monotonous response of the sublimation to the solar irradiation together with the normal direction at each surface point establish the vector fields for forces and torques on the surface. The torque vector field can be reinterpreted as a weighted superposition of the vectorial torque efficiency. This quantity depends on the short axis vector and local properties of the shape geometry only. One scalar component of this torque efficiency that affects the rotation period is introduced by Keller et al. (2015b). Jewitt (1997) and Samarasinha & Mueller (2013) consider surface integrated torque terms to reflect cancelation effects. By integrating spatially the vector field for torque, the resulting torque vector influences the rigid body dynamics. This depends on the body geometry, on the heliocentric distance, and on the orientation state (the illuminated cross section) of the body relative to the sun. Already for a monotonous TPM, surface features and self-shadowing can shift the maximum sublimation activity away from the current subsolar point. Qualitatively, a resulting force caused by irradiation pushes the nucleus into the antisolar direction. For the torque and changes of the rotation state, such a simple relation does not persist since the shape of the nucleus and the orientation of the rotation axis have to be considered in detail. For a lowest-energy rotation, Kramer et al. (2019) parameterize diurnal torques by the subsolar longitude and derive a torque decomposition as a product of an orbital term and a shapedependent matrix. Subsolar latitude, together with solar distance, is associated with seasonal variations.

Although uniform activity on the complex shape of comet 67P/C-G explains observed changes for the rotation period, there is a mismatch between this activity and observations for the axis orientation. For this problem within a model for comet 67P/C-G, a solution is to prescribe a domain decomposition consisting of surface regions and to fit the corresponding activity weights (effective active fraction). These weights are only weakly constrained and it is difficult to obtain appropriate decompositions of the surface. For any cometary body geometry, the vectorial torque efficiency (with its local characterization of the geometry) offers the option to specify a domain decomposition consisting of eight regions. For lowest-energy rotations and independent of the shape, these regions have the significant property that all points of one region share the same direction for torque formation.

The manuscript describes our approach to decomposing a cometary surface into eight regions that determine the direction of the torque and the resulting evolution of the rotation state close to lowest energy. Sect. 2 introduces the dynamical equations and discusses the driving forces for rotational evolution. We consider instantaneous TPMs with a weighted activity on the cometary surface. Sect. 3 describes the formation of a torque field at each point on the cometary surface. The relation between the specified regions and the torque direction is examined for representative examples, namely, for an oblate ellipsoid, a nearprolate ellipsoid, a bilobed shape, and a shape model of comet 67P. These torques guide the evolution of the rotation state close to lowest energy in Sect. 4. Rotational data for comet 67P/C-G serve as constraints for effective active fractions with respect to our domain decomposition.

2 Dynamics of rotation

The rotation state of a rigid body can be described by the temporal evolution of a coordinate mapping x(ξ, t) from the body frame to the inertial frame at time t with a position vector ξ ∈ ℝ³ in the body frame. x is a linear mapping x(ξ, t) = Rξ represented by an orthogonal rotation matrix R(t) ∈ ℝ^{3 × 3}. Within the body frame, a tensor of inertia J_BF ∈ ℝ^{3 × 3} is prescribed. The corresponding tensor is J(t) = R J_BF R⁻¹ in the inertial frame. Based on the vectors for the angular velocity, ω(t) ∈ ℝ³, for the rotational angular momentum, L(t) ∈ ℝ³, and the prescribed torque, T(t) ∈ ℝ³, the equations for rotation in the inertial frame are $$\dot{R}()=\omega \times R(),\dot{L}=T,J\omega =L.$$(1)

Together with the initial conditions, R(t₀) = R₀ and L(t₀) = L₀, at time t₀, this equation poses an initial value problem that determines the temporal evolution of the rotation state, R(t) and L(t). From the solution of the differential equation, additional variables such as ω, x(ξ, t) are determined. We take as directions of the inertial frame the J2000 system, with the origin at the center of mass of the cometary body. Far away from further gravitationally interacting bodies, this translation does not lead to further torques or modifications of Eq. (1). Within this manuscript all results are based on numerical simulations for the full rotational dynamics.

The equation in the inertial frame has an equivalent representation in the body frame with the Euler equations for rotation, $$J_{\mathrm{B}\mathrm{F},ij}{\dot{\omega }}_{\mathrm{B}\mathrm{F},j}=T_{\mathrm{B}\mathrm{F},i}-{ϵ}_{ijk}{\omega }_{\mathrm{B}\mathrm{F},j}{J}_{\mathrm{B}\mathrm{F},kl}{\omega }_{\mathrm{B}\mathrm{F},l},$$(2)

for variables ω_BF = R⁻¹ ω and T_BF = R⁻¹ T.

2.1 Rigid body geometry

Within the body frame with the standard basis vectors, e₁, e₂, e₃ ∈ ℝ³, the spatial domain Ω ⊂ ℝ³ describes the geometrical shape (body shape) of the cometary nucleus. In the inertial frame at time t the location of the body is the spatial domain Ω_t = x(Ω, t). We consider three idealized body shapes, Ω, complemented by one additional shape, an approximation of the cometary nucleus of comet 67P/C-G. The first two shapes are ellipsoids and are oriented such that their symmetry axes coincide with e₁, e₂, e₃. The dimensions of the ellipsoids are characterized by the lengths of the semimajor axis a₁ ≥ a₂ ≥ a₃ with respect to e₁, e₂, e₃. The first shape is an oblate ellipsoid with a₁ = a₂ = 2 a₃ and is therefore axially symmetric around e₃. An oblate ellipsoid underlies the forced precession models discussed in Whipple & Sekanina (1979) and Sekanina (1984), together with further assumptions on the activity of regional jets on the surface. The second shape is a near-prolate ellipsoid characterized by 9a₁ = 18a₂ = 20a₃ without axial symmetry. As a consequence of three different axis lengths this shape is not axially symmetric. That is why we prefer the denotation near-prolate instead of prolate. For the numerical study here, we describe the boundaries of Ω for these two shapes using a surface, ∂Ω, consisting of a triangular mesh with N_E = 4000 elements. The third geometry is given by a bilobed shape and consists of two slightly oblate ellipsoids placed side by side without overlap. Each ellipsoid satisfies the relation 9a₁ = 9a₂ = 10a₃ such that the bilobed shape holds the same spatial dimensions as the near-prolate ellipsoid. This bilobed shape configuration could represent a contact binary and is similar to the shape of comet 67P/C-G. The base vector e₁ is located on the line through the centers of the spheres and e₂, e₃ span the perpendicular plane. The surface of the domain, Ω, for the bilobed shape consists of N_E = 4000 triangular elements. The last case is a shape model of comet 67P/C-G. The N_E = 3996 surface elements of the last shape are the same as those of Läuter et al. (2022), which are derived from the high-resolution (meter-scale) shape model of Preusker et al. (2017). All shapes share the same normalized volume |Ω| = 18.56 km³. Based on these choices, all shapes share the same spatial resolution of approximately Δ x ≈ 150 m.

The components of the tensor of inertia are given by $$J_{\mathrm{B}\mathrm{F},ij}=\underset{\mathrm{\Omega }}{\int }\rho ({\delta }_{ij}{\xi }^2-{\xi }_i{\xi }_j)\mathrm{d}\xi .$$(3)

For the two ellipsoids and for the bilobed shape within domain Ω, we assume the constant density ρ = 539 kg/m³. For shape 67P, the tensor of inertia J_BF is set to Eq. (A.1). For this case, the assumption of uniform density is suspended, in good agreement with the gravitational constraints of Laurent-Varin et al. (2024). For all shapes considered, e₃ is the principal axis of J_BF with the largest eigenvalue, J₃. The vector e₃ in the body frame corresponds to the principal axis, Re₃, in the inertial frame. In the case of ellipsoids the moments of inertia are ordered in inverse proportion to the lengths of the semimajor axes. The principal axis e₃ corresponds to the largest moment of inertia and to the shortest axis. Let us denote b₃ = Re₃ as the short axis in the inertial frame.

For a lowest-energy rotation state (represented by R and L), the directions of L, of ω, and of the short axis coincide. During the 2015 apparition Godard et al. (2017), Preusker et al. (2017), Jorda et al. (2016), and Gutiérrez et al. (2016) consider comet 67P/C-G to be in lowest-energy rotation. A free nutation motion for the angular velocity vector has not been observed within the accuracy range of a few tens of millidegrees. A lowest-energy rotation implies ω = |ω| b₃ for the angular velocity and L = J₃|ω| b₃ for the rotational angular momentum. In the following, the initial states, R₀ and L₀, are chosen to be at the lowest energy. With the comet becoming active, the NGTs introduce some perturbations leading to an excited state. A torque, T = T₃ b₃, directed along b₃, is able to change the vector length |ω| (and thus the rotation period), while the directions of b₃, ω, and L remain unchanged. A torque, T, that is orthogonal to b₃, namely with T · b₃ = 0, breaks the initial alignment of b₃ and ω in general. For linearized Euler equations in Eq. (2), an analysis for the angle γ between b₃ and ω yields the estimation tan γ ≈ |T|/J₃/ω². Assuming characteristic values for ice sublimation on comet 67P/C-G with |T| = 10⁹ Nm, J₃ = 10¹⁹ kg m², and |ω| = 10⁻⁴ s⁻¹ yields an upper bound |γ| ≈ 0.5^∘.

To analyze the orientation of vector quantities (e.g., torque and angular velocity) relative to the sun and to the short axis, we defined two different basis systems with the (orthogonal) vectors $$b_1(y)=\frac{(b_3\times y)\times b_3}{|(b_3\times y)\times b_3|},b_2(y)=\frac{b_3\times y}{|b_3\times y|},b_3=R{e}_3.$$(4)

We chose the vector y ∈ ℝ³ in two different variants that yield the local basis and the solar basis, respectively. As shorter denotation we also use B(y) = (b₁(y), b₂(y), b₃).

For the local basis at time t and a position vector x ∈ ∂Ω_t on the surface, we defined y = v(x, t). Fig. 1 shows the outward normal vector v on ∂Ω_t. At a fixed time t the local basis, B(v), depends on the surface position x and on the principal axis b₃. With advancing time t at a fixed position ξ ∈ ∂Ω, this basis in the body frame (namely R⁻¹ b₁(v), R⁻¹ b₂(v), R⁻¹ b₃) is constant. This means that the local basis follows the rotational movement in the inertial frame together with x = Rξ. The three vectors b₁, v, and b₃ are located in a shared plane. For the solar basis, B(r_s), in Eq. (4) we define y = r_s(t) with the vector r_s pointing toward the sun at time t. Fig. 2 shows this geometric setup. At a fixed time t the basis is constant along the surface, ∂Ω_t. With advancing time t, it does not follow the rotational movement in the inertial frame and it rotates in the body frame. The vector b₁ is related to the solar direction; that is, such that b₁, r_s, and b₃ are located in a shared plane. Except for the signs for a lowest-energy state, this solar basis was already applied by Kramer & Läuter (2019) and Sekanina (1981). Kramer et al. (2019) apply this basis to represent the diurnally averaged torque affecting the rotation state. The local and solar basis share the property that b₃ is aligned to the short axis so that the equatorial plane of the rotation is spanned by the vectors b₁ and b₂. The equatorial planes of both basis systems are the same. b₁(v), b₂(v) and b₁(r_s), b₂(r_s) are just two choices for a basis of this plane.

Fig. 1 Local basis, B(v) = (b1(v), b2(v), b3), for the choice y = v(x, t) in Eq. (4) at time t and at position x ∈ ∂Ωt with the outward normal v.

Fig. 2 Solar basis, B(rs) = (b1(rs), b2(rs), b3), for the choice y = rs(t) in Eq. (4) at time t.

2.2 Domain decomposition of the surface

For each ξ ∈ ∂Ω on the surface, x = Rξ, v(x, t), and b_i(v) (the local basis vector) follow the same rotational movement in the inertial frame. Based on these quantities, we define regions on ∂Ω_t that move correspondingly in the inertial frame. Their counterparts in the body frame are time-independent, so these definitions are based on geometry properties of the shape in the body frame, ∂Ω.

At time t we define the southern and the northern hemisphere, $$A_{\text{SH}}=\{x\in \mathrm{\partial }{\mathrm{\Omega }}_t|\nu \cdot b_3<0\},A_{\text{NH}}=\{x\in \mathrm{\partial }{\mathrm{\Omega }}_t|0<\nu \cdot b_3\}.$$

Their specification is based on the sign of the b₃ component for the surface normal, v, at point x ∈ ∂Ω_t. For a lowest-energy rotation state, this projection controls the seasonal intensity of the solar irradiation on the surface, ∂Ω_t. For negative subsolar latitude SSL = arcsin (r_s · b₃/|r_s|), the point of maximum irradiation is located in A_SH. For positive SSL, this point is located in A_NH. This relation expresses the mapping between A_SH and A_NH on the one side and both hemispheric seasons on the other side.

This concept of defining regions based on projections to the local basis in Fig. 1 can be applied to the vector x × v. We define the vectorial torque efficiency by the three scalar components, $$(x\times \nu )\cdot b_i,$$(5)

for i = 1, 2, 3. Fig. 3 shows the maps of the vectorial torque efficiency for the near-prolate ellipsoid. Fig. 4 shows the b₂ component for the more complex geometries of the bilobed shape and shape 67P. In Sect. 3 we show that both figures represent components of the vectorial torque efficiency. Based on these fields for i, j, k = ± 1, we define eight regions, $$A_{i,j,k}=\{x\in \mathrm{\partial }{\mathrm{\Omega }}_t|\begin{array}{c}i(x\times \nu )\cdot b_1<0,\\ j(x\times \nu )\cdot b_2<0,\\ k(x\times \nu )\cdot b_3<0\end{array}\},$$(6)

and a complementing zero region, $$A_0=\{x\in \mathrm{\partial }{\mathrm{\Omega }}_t|\begin{array}{c}(x\times \nu )\cdot b_1=0\text{or}\\ (x\times \nu )\cdot b_2=0\text{or}\\ (x\times \nu )\cdot b_3=0\end{array}\}.$$

The analysis of signs in the three plots of Fig. 3 yields these regions in Fig. 5 for the near-prolate ellipsoid. In this case, the regions A_1,−1,−1, A_−1,1,−1, A_−1,−1,1, and A_1,1,1 are empty, which yields the decomposition ∂Ω_t = A₀ ∪ A_{−1,−1,−1} ∪ A_1,1,−1 ∪ A_1,−1,1 ∪ A_−1,1,1 into four regions and A₀. For the oblate ellipsoid, the rotational symmetry around b₃ yields (x × v) · b₁ = (x × v) · b₃ = 0 and finally ∂Ω_t = A₀. These two examples show that symmetries can lead to empty regions A_{i, j, k}.

In general, all eight surface regions, A_i,j,k, in Eq. (6) are nonempty, as shown in Figs. 6, 7 for the bilobed shape and for the shape 67P. Applying Eq. (A.2), each region is a subset of A_SH or A_NH at most, and the following decompositions hold. $$\begin{array}{rl}& A_{\text{SH}}=A_{-1,\pm 1,-1}\cup A_{1,\pm 1,1}\cup A_0,\\ & A_{\text{NH}}=A_{-1,\pm 1,1}\cup A_{1,\pm 1,-1}\cup A_0.\end{array}$$(7)

Both properties can be recognized in Figs. 5, 6, and 7. Looking straight from the south or from the north, only four of the eight regions are visible. For a lowest-energy rotation state, Eq. (7) expresses that each hemispheric season is related to the predominant irradiation of one of both subsets.

Based on ice sublimation, torque creation will be discussed in the next section. For that, regionally differing weights for activity will be considered. This is realized by the function f_EAF for the effective active fraction on the surface, ∂Ω_t. The uniform activity is represented by the constant function f_EAF = 1. Based on the decomposition of ∂Ω_t into the regions A₀ and A_{i, j, k}, the effective active fraction is defined by $$f_{\text{EAF}}=\sum _{i,j,k=\pm 1}{f}_{\mathrm{E}\mathrm{A}\mathrm{F},i,j,k}{\chi }_{A_{i,j,k}}.$$(8)

χ_Ai,j,k is the characteristic function on the region A_{i, j, k}. The weights f_EAF,i,j,k are open parameters for f_EAF, which need to be restricted by the observational data.

Fig. 3 For the near-prolate ellipsoid, the values (x × v) · bi(v) (vectorial torque efficiency) are shown based on longitude, λ, and latitude, ϕ, in the body frame. bi(v) is the local basis vector in Fig. 1.

Fig. 4 For the bilobed shape (top) and for shape 67P (bottom), the values (x × v) · b2(v) (second component of the vectorial torque efficiency) are shown based on longitude, λ, and latitude, ϕ, in the body frame. b2(v) is the local basis vector in Fig. 1.

2.3 Thermophysical models

The torque, T, in Eq. (1) modifies the rotational angular momentum, L, directly and, as a consequence, controls the evolution of the full rotation state. On the nucleus, Ω_t, solar irradiation causes ice sublimation close to the surface, which leads to gas expansion into the surrounding space. On a surface point, we denote as a TPM the mappings from solar irradiation, I, (and its temporal variation) to mass generation (due to sublimation), ṁ, and to gas velocity, v_gas. Thermophysical models consider complex dynamical process for material, energy and momentum on the surface and within the volume of the nucleus below. Interactions with subsurface layers can introduce thermal inertia, which causes changes in I to appear for ṁ and v_gas with a temporal delay.

For the local evaluation of solar irradiation, I, we take into account the shadowing of surface areas due to the concave shape of Ω_t. If there is no additional intersection between the line of sight from the surface position x ∈ ∂Ω_t to the sun at position r_s, the solar irradiation is given by $$I(x,r_{\text{s}})=I_{\text{s}}{\left(\frac{1\mathrm{a}\mathrm{u}}{|r_{\text{s}}|}\right)}^2\frac{\nu \cdot r_{\text{s}}}{|r_{\text{s}}|}$$(9)

with the solar constant I_s. In the following, we restrict ourselves to instantaneous TPMs without thermal inertia. This means that both functions ṁ(I) and v_gas(I) are functions of solar irradiation, I, only. Another assumption for our TPMs is that both functions ṁ(I) and v_gas(I) increase monotonously with respect to I.

For quantification, we consider two specific choices for instantaneous TPMs, model A and model α. Model A is based on the instantaneous energy balance of a irradiated surface covered by pure ice. Our model A is the TPM for water of Läuter et al. (2022). For less intense irradiation, I, of up to 50 W/m² the mass production is close to zero. For higher intensities in the range of up to 1000 W/m² this functional relation is close to a linear relation. The second TPM model α with the functions ṁ_α(I) and v_α(I) has been introduced by Kramer et al. (2019). This approach goes back to the observation during the 2015 perihelion passage of comet 67P/C-G, that the gas production increases more rapidly than predicted from model A. For an exponent 1 ≤ α we defined the TPM by $${\dot{m}}_{\alpha }(I)={\dot{m}}_0{\left(\frac{I}{I_0}\right)}^{\alpha },v_{\alpha }=450\mathrm{m}/\mathrm{s}.$$

The constant ṁ₀ = 3.47 × 10⁻⁴ kg/s/m² was chosen so that ṁ_α and model A yield the same sublimation rate at irradiation I₀ = 1000 W/m². Although any choice for α is possible, we restrict our consideration to the value α = 2.

Assuming the gravitational acceleration acting on the body is constant along the body domain, the torque caused by gravitation is zero. As a consequence, the torque vector, T(t), is restricted to nongravitational contributions at time t. Based on an instantaneous TPM, gas expansion starts from each point x on the surface, ∂Ω_t. The torque, T, and the vector field, t_NG, for the torque density on ∂Ω_t are defined by $$T(t)=\underset{\mathrm{\partial }{\mathrm{\Omega }}_t}{\int }{t}_{\text{NG}}\mathrm{d}\sigma ,t_{\text{NG}}(x,t)=-f_{\text{EAF}}{v}_{\text{gas}}\dot{m}x\times \nu .$$(11)

Because the velocity ω × x (related to the angular velocity, ω) is small compared to v_gas v, we have omitted this fraction of the velocity. The momentum transfer coefficient applied by Rickman (1989), Gutiérrez et al. (2005), and Attree et al. (2019) is also omitted. Because it is a constant we include it in f_EAF, which is still an open parameter. Eq. (10) expresses a relation between the intensity of solar irradiation and the strength of the torque. In particular, each of both hemispheric seasons is now related to torque activity in one of the two separate subsets in Eq. (7).

Instantaneous TPMs create torque without delay in time. That is why the solar direction of the vector r_S and with it the vector b₁(r_s) are associated with the preferred direction of forces on the surface. For rotational dynamics, these forces lead to a preferred torque direction in Eq. (10) associated with ± b₂(r_s) orthogonal to the solar direction. Any vector v ∈ ℝ³ can be represented by the three components v · b_i(r_s). With these components, longitude, λ_s, and latitude, ϕ_s, of v with respect to the solar basis are defined by λ_s = arctan (v · b₂(r_s)/v · b₁(r_s)) and ϕ_s = arcsin (v · b₃/|v|). The other way around, knowing the solar basis, both angles λ_s and ϕ_s determine the direction of v. The case longitude λ_s ∈{0^∘, 180^∘} corresponds to v · b₂ = 0. In this case, the solar direction r_s and v are closely related in the way that both projections to the equatorial plane coincide. The case λ_s ∈{−90^∘, 90^∘} corresponds to v · b₁ = 0 such that the equatorial projections of r_s and v are perpendicular. For latitude ϕ_s = 0^∘, v · b₃ = 0 holds and v is in the equatorial plane. Finally, ϕ_s ∈{−90^∘, 90^∘} gives v · b₁ = v · b₂ = 0 and the representation v = ±|v| b₃. For the special case of v = r_s, ϕ_s is the same as SSL.

Fig. 5 For the near-prolate ellipsoid decomposition of ∂Ωt into the four regions Ai, j, k in Eq. (6). Region A0 is colored in white. Top: longitude, λ, and latitude, ϕ, in the body frame. Bottom: 3D projection using the same color code.

Fig. 6 For the bilobed shape (top) and for the shape 67P (bottom), decomposition of ∂Ωt into the eight regions Ai,j,k in Eq. (6). Region A0 is colored in white. Longitude, λ, and latitude, ϕ, in the body frame.

Fig. 7 For the bilobed shape (top) and for the shape 67P (bottom), decomposition of ∂Ωt into the eight regions. The same data as in Fig. 6 but in 3D projection.

2.4 Data for the rotation state of 67 P/C−G

For a numerical solution of the initial value problem of Eq. (1) a set of parameters needs to be specified. These include the initial conditions, the tensor of inertia, and the mechanism for torque formation in Eq. (10). The solution of the inverse problem consists of the determination of the open parameters based on known observational data. Angular velocity, ω, was observed during the ESA mission of comet 67P/C-G during the apparition in 2015. This rotational information refers to the angular frequency, f = |ω|, of the rigid body, as well as to the directional angles ra(ω), dec(ω) of the rotation axis, ω/|ω|. Our data set is based on the measurements for ω discussed in Kramer et al. (2019). In the equatorial frame J2000 (EME2000) the angles are extracted based on the method of Gaskell et al. (2008) applied by Jorda et al. (2016) to comet 67P/C-G. Data sets data(ω) = (ra(ω), dec(ω), |ω|) are averaged along a moving time interval of 21 days and are considered between August 2014 and July 2016. Time sampling defines N_time = 58 time points t_i with i = 1, …, N_time. The sampling rate varies between a 7 day to a 33 day distance between two times with a higher sampling rate close to perihelion passage in August 2015. The uncertainty of the data at time t is estimated based on the spread of the data close to that time. This leads to uncertainties between 0.01^∘ and 0.11^∘ for dec(ω) and between 11 s and 60 s for rotation period 2π/f. The latter range is much higher compared to uncertainties on the order of 1 s for periods from Godard et al. (2017) and Mottola et al. (2014). We chose these higher weights to distribute the uncertainties uniformly in the direction and in the length of the vector, ω. We combined these preparations into the data set, $$(t_i,dat{a}_i,{\sigma }_i)\text{for}i=1,...,N_{\text{time}},$$(11)

of 3N_time data points data_i = data(ω_i) ∈ ℝ³ and uncertainties σ_i ∈ ℝ³ at time t_i.

3 Torque formation for lowest-energy rotation

The torque, T, in Eq. (10) depends on the geometric properties of the shape Ω_t, on the TPM, and on the weights of the effective active fraction f_EAF in Eq. (8). The sensitivity of T to Ω_t and f_EAF is discussed in the present section. Ω_t specifies the orientation of the surface regions and, along with it, the formation of the torque density on the surface. The uniform activity with f_EAF = 1 represents the dynamics of the complete shape. The combination of spatially weighted effective active fractions in Eq. (8) and complex shape geometries gives options to modify the rotational dynamics. Within this section we restrict the discussion to lowest-energy rotation states. As a consequence, two significant vectors share the same direction, the rotation axis, ω, and the principal axis, Re₃ (the same as b₃ in Eq. (4)).

The discussion of parameters acting on torque generation is based on the torque components, T · b_i(r_s), with respect to the solar basis in Fig. 2. Because torque, T(t), in Eq. (1) determines the derivative $\dot{L}$, the dynamical effect of T is represented by its integral or its average T̅ over time. Due to the rotation of point ξ ∈ ∂Ω with outward normal ν_ξ in the body frame, the surface point x(t) = R(t) ξ ∈ ∂Ω_t, the outward normal v(t) = R(t) v_ξ, the irradiation, I ∼ r_s · v(t), the gas velocity, v_gas, and the mass production, ṁ, change in time. Consequently, in general, the torque, T, is not constant in time. In particular this is true for the nearprolate ellipsoid and the bilobed shape holding fewer rotational symmetry than the oblate ellipsoid. At time t for T and t_NG, the temporal averages over one rotation period Δ t can be expressed by $$\overline{T}(t)=\underset{\mathrm{\partial }{\mathrm{\Omega }}_t}{\int }\overline{t_{\text{NG}}}\mathrm{d}\sigma ,\overline{t_{\text{NG}}}(x,t)=\frac{1}{\mathrm{\Delta }t}\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}{t}_{\text{NG}}(x(s),s)\mathrm{d}s.$$(12)

The average of t_NG is evaluated along the trajectory x(s) ∈ ∂Ω_s corresponding to ξ = R⁻¹(t) x and v_ξ = R⁻¹(t)v. Considering constant ω and r_s, Eq. (B.3) replaces the term $\overline{t_{\text{NG}}}(x,t)$, which gives $$T_i=b_i(r_{\text{s}})\cdot \overline{T}=-\underset{\mathrm{\partial }{\mathrm{\Omega }}_t}{\int }{f}_{\text{EAF}}{c}_i{b}_i(\nu )\cdot (x\times \nu )\mathrm{d}\sigma .$$(13)

As is described in Appendix B, the numbers c_i are Fourier coefficients for the source term v_gas ṁ representing the TPM. A similar idea has been introduced by Kramer et al. (2019) for the diurnal and seasonal sublimation cycles. There, three Fourier components are derived for each surface region that produce the torque movement. On the left-hand side of Eq. (13) are the components of torque, T̅, with respect to the solar basis, B(r_s), in Fig. 2. These represent directions toward and perpendicular to the sun as described in Sect. 2.3. On the right-hand side, T_i is constructed by components with respect to the local basis, B(v), in Fig. 1. Adjacent to the TPM terms, the integrand in Eq. (13) contains the vectorial torque efficiency, (x × v) · b_i(v), defined in Eq. (5).

Fig. 8 Direction of torque, T̅, for the near-prolate ellipsoid. Torque, T̅, in Eq. (12) is shown with model A and fEAF = χAi,j,k for region Ai,j,k in Eq. (6). Uniform activity fEAF = 1 is shown with dotted black lines for the oblate ellipsoid and for the near-prolate ellipsoid. λs and ϕs denote longitude and latitude for T̅ with respect to the solar basis in Fig. 2. The solar vector, rs, changes SSL from −45∘ to 45∘.

3.1 Idealized shape geometries

From Eq. (3) we know that the shape Ω determines the tensor of inertia J_BF together with the density function ρ. Complementing this, Ω_t influences torque formation in Eqs. (10), (12), and (13). In this section, we assume an instantaneous TPM and a constant solar position, r_s.

We discuss the oblate ellipsoid with uniform activity f_EAF = 1 on the surface. The rotational symmetry around b₃ yields different properties. The torque, T, is constant in time and finally T = T̅. Both projections (x × v) · b₁(v) and (x × v) · b₃ vanish everywhere on ∂Ω_t with the consequence T₁ = T₃ = 0. For southern points x ∈ A_SH, the feature 0 < (x × v) · b₂(v) holds. In the case of negative SSL (r_s · b₃ < 0), the southern region A_SH receives more irradiation from the sun so that v_gas ṁ has higher values than the northern region. Then the integral in Eq. (13) is dominated by the southern contribution such that T₂ < 0. The same argument leads to 0 < T₂ for positive SSL. Using Eq. (13) to evaluate T · b_i(r_s) yields ϕ_s = 0^∘ and directions of T that are limited to ± b₂(r_s). Fig. 8 shows λ_s, which flips from −90^∘ to 90^∘ when SSL changes from negative to positive.

The near-prolate ellipsoid is more complex, such that the torque, T(t), is a function of time with uniform activity f_EAF = 1. Fig. 3 shows the components (x × v) · b_i(v) for the nearprolate ellipsoid. We discuss the plots for the components b₁ and b₃. Along circles of constant latitude, ϕ_s, for the integrand in Eq. (13), the distribution is symmetric with respect to zero, such that the integrals are zero, T₁ = T₃ = 0, and finally T̅ · b₁(r_s) = T̅ · b₃ = 0. T₂ in Eq. (13) does not vanish because canceling for (x × v) · b₂(v) in Fig. 3 does not happen along these line integrals. This term shows the same signs for points on the hemispheres, A_SH and A_NH, as discussed for the oblate ellipsoid. Finally, the directions of T agree for both ellipsoids in Fig. 8.

The properties of the bilobed shape are similar to those of the near-prolate ellipsoid for uniform activity f_EAF = 1. T(t) is a function of time. The two components T̅ · b₁(r_s) = T̅ · b₃ = 0 vanish. The corresponding plot of (x × v) · b₂(v) for the bilobed shape as in Fig. 3 suggests cancelation even for T₂. Shadowing effects must be considered so that the symmetry of v_gas ṁ is perturbed. The dotted line in Fig. 9 shows the resulting longitude, λ_s, of T̅. When comparing the magnitudes of T̅, the torque is smaller for the bilobed shape compared to the other shapes. This goes back to the directional signal caused by shadowing and is more sensitive compared to the geometrical effects of the other shapes.

For all components T_i in Eq. (13) the integrands on the right-hand side have nonzero contributions, which is shown in Fig. 3 for the near-prolate ellipsoid (without the term v_gas ṁ). As a consequence effective active fractions f_EAF can influence torque formation. We focus our approach on the superposition of regional contributions in Eq. (8). T_i(A_{i, j, k}) denotes the evaluation of T_i in Eq. (13) assuming the choice f_EAF = χ_Ai,j,k. Because region A_{i, j, k} collects surface points that share the same sign for the integration of T_l(A_i,j,k), the sign of T_l can be determined with this choice. The definition of the regions in Eq. (6) yields the properties $$\text{sgn}{T}_1(A_{i,j,k})=i,\text{sgn}{T}_2(A_{i,j,k})=j,\text{sgn}{T}_3(A_{i,j,k})=k.$$(14)

Based on Eq. (13) this extends to the signs of T̅(A_{i, j, k}) · b_1,2,3(r_s). For the bilobed shape with model A, Fig. 9 shows the longitude, λ_s, and latitude, ϕ_s, of T̅(A_{i, j, k}). The same figure for the bilobed shape with model α looks similar, leading to the same conclusion. Eq. (14) confirms the assignment of λ_s and ϕ_s to the angular quadrants. For example T₁(A_{−1,−1, ± 1}) and T₂(A_{−1,−1, ± 1}) both feature negative signs that correspond to the range −180^∘ < λ_s < −90^∘ for the longitude of the vector T̅(A_{−1,−1, ± 1}). Fig. 9 confirms this property. A closer look at λ_s and ϕ_s reveals that they are close to the middle of the expected ranges. For the example A_{−1,−1, ± 1}, λ_s is close to −135^∘. Deviations from the ideal situation go back to the concave shape of this example, leading to self-shadowing, which perturbs the perfect illumination conditions. The approximation in Eq. (13) is applicable for this case. For the oblate and near-prolate ellipsoid, this full control of λ_s and ϕ_s for T̅ cannot be expected. For the oblate ellipsoid, all regions A_{i, j, k} vanish and the torque direction is limited to values λ_s = ± 90^∘. For the near-prolate ellipsoid in Fig. 8 only four of the eight regions are present.

Fig. 9 Direction of torque, T̅, for the bilobed shape. Torque, T̅, in Eq. (12) is shown with model A and fEAF = χAi,j,k for region Ai, j, k in Eq. (6). Uniform activity fEAF = 1 is shown with dotted black lines. λs and ϕs denote longitude and latitude for T̅ with respect to the solar basis in Fig. 2. The solar vector, rs, changes SSL from −45∘ to 45∘.

Fig. 10 Direction of torque, T̅, for the shape 67P. Torque, T̅, in Eq. (12) is shown with model A and fEAF = χAi,j,k for region Ai,j,k in Eq. (6). Uniform activity fEAF = 1 is shown with dotted black lines. λs and ϕs denote longitude and latitude for T̅ with respect to the solar basis in Fig. 2. The solar vector, rs, changes SSL from −45∘ to 45∘.

3.2 Shape of comet 67P/C-G

The shape of comet 67P is more complex than any of the idealized shapes. However, the torque formation shares some properties with simplified models, while lack of symmetry brings about new features. According to the discussion in the last section, we consider an instantaneous TPM. No component T_i in Eq. (13) is zero and no integrand (x × v) · b_i(v) vanishes everywhere on ∂Ω_t.

From Figs. 8 and 9 for both ellipsoids and the bilobed shape, we know that the longitude, λ_s, for T̅ flips between −90^∘ and 90^∘ for changing SSL. For uniform activity f_EAF = 1 we investigate this property for shape 67P with Eq. (13). Fig. 4 shows the component (x × v) · b₂(v) on the surface, ∂Ω_t. Except for the far southern latitudes (small regions around ϕ = −60^∘), positive signs are predominant. For irradiation conditions coming from the north (positive SSL), sublimation controlled by v_gas ṁ is most active on A_NH and thus T̅ · b₂(r_s) ≈ T₂ < 0. The same argument holds for irradiation coming from mid latitudes even from the south. For longitude, λ_s, of T̅ this yields the expectation −180^∘ < λ_s < 0^∘ in a wide range.

For uniform activity, the symmetry of all shapes in the last section leads to vanishing values T̅ · b₃ ≈ T₃ = 0. This feature is different for shape 67P. Fig. 10 shows the latitude, ϕ_s, for T̅ as a falling function with respect to SSL. In SSL ≈ 10^∘, ϕ_s changes the sign from positive to negative and with it T̅ · b₃. The southern irradiation conditions (with SSL < 10^∘) related to the predominant sublimation on A_SH imply 0 < T̅ · b₃. The northern irradiation conditions (with SSL > 10^∘) imply T̅ · b₃ < 0. This observation corresponds to an increasing (decreasing) rotation period before (after) the first equinox of the 2015 apparition of comet 67P/C-G as described by Keller et al. (2015a). Equation (A.2) is a pointwise relation between the components b₁(v) and b₃ of x × v. This relation does not transfer to the complete integrals T₁ and T₃ (and not to T̅ · b₁ and T̅ · b₃) in Eq. (13). We consider this as an indication for T̅ · b₁ and T̅ · b₃, only for equally distributed surface activity. For southern irradiation conditions (and negative SSL) we expect sgn T̅ · b₁ = sgn T̅ · b₃. Due to 0 < T̅ · b₃ we obtain 0 < T̅ · b₁(r_s) for these conditions. The same argument for northern irradiation conditions (and positive SSL) yields the same expectation 0 < T̅ · b₁(r_s). For longitude, λ_s, of T̅ this gives the constraint −90^∘ < λ_s < 90^∘. Together with the analysis of the b₂ component in the last paragraph, we obtain the expectation −90^∘ < λ_s < 0^∘. Fig. 10 confirms this, showing that the course of λ_s is in the range between −90^∘ and −70^∘.

For shape 67P the eight regions A_{i, j, k} defined in Eq. (6) are shown in the Figs. 6 and 7. The orientations of T̅(A_{i, j, k}) are shown in Fig. 10 based on effective active fractions f_EAF = χ_Ai,j,k. The longitude, λ_s, and latitude, ϕ_s, of the different T̅ are within the ranges −135^∘ < λ_s < 135^∘ and −75^∘ < ϕ_s < 75^∘. This is similar to Fig. 9 for the bilobed shape. The appropriate weights for the superposition in Eq. (8) will be able to create torque formation for the prescribed values of λ_s and ϕ_s. Because the argumentation in Sect. 3.1 is independent of the specific shape, Eq. (14) is true for shape 67P as well.

4 Rotational angular momentum and velocity

We studied the evolution of the angular velocity, ω. Similarly to the discussion for T, the vector ω can be decomposed into components with respect to the solar basis, B(r_s), in Fig. 2. All initial states within the present section, R₀ and L₀, are assumed to be in lowest energy. In this case, the axes b₃, ω, and L (principal axis, rotation axis, and rotational angular momentum) are aligned to each other. The presence of torque, T, causes an excitation for the rotation state such that b₃, ω, and L are no longer parallel. For rotation states close to lowest energy, the torque formation discussed in Sect. 3 remains valid qualitatively. The longitude, λ_s, and the latitude, ϕ_s, of ω with respect to the solar basis are one option to express the alignment between ω and b₃. For the oblate ellipsoid, the near-prolate ellipsoid, the bilobed shape and the shape 67P, we considered model A with f_EAF = 1/8. For the shape 67P, we considered a second variant based on the model α with f_EAF = 1/4.

To study the evolution of the rotation state, we studied the formation of the torque, T, according to Eq. (10). The first setup considers a prescribed vector for the solar direction r_s, a constant heliocentric distance |r_s| = 2 au and a changing SSL from −45^∘ to 45^∘ during a time span of 800 days. For this setup, the latitude, ϕ_s, shows small deviations in the range 89.99^∘ ≤ ϕ_s ≤ 90^∘. This significantly improves the analysis of the linearized equations (γ ≈ 0.5^∘) in Sect. 2.1 and shows the close alignment of the vectors ω and b₃. The rotation remains close to a lowest-energy state. At the same time, the longitude, λ_s, for ω covers the entire angular range −180^∘ ≤ λ_s ≤ 180^∘ that describes the precessional movement of ω surrounding b₃.

Fig. 11 Longitude, λs, and latitude, ϕs, for angular velocity, ω, with respect to the fixed solar basis in Fig. 2, based on rs(t0) and b3 = R(t0) e3 at initial time, t0 . T is evaluated with |rs| = 2 au and changing SSL from −45∘ to 45∘. Considered are the oblate ellipsoid, the near-prolate ellipsoid, the bilobed shape, the shape 67P with model A/fEAF = 1/8, and the shape 67P with model α/fEAF = 1/4.

4.1 Uniform activity

We consider a uniformly distributed surface ice coverage. In this case, the principal axis, b₃ = R(t) e₃, changes in time. To plot the evolution of ω in figures with respect to an inertial frame, we fix the solar basis of Fig. 2 based on the state at initial time t₀, namely with r_s(t₀) and b₃ = R(t₀) e₃. Fig. 11 shows λ_s and ϕ_s for ω with respect to this fixed solar basis. Due to lowest-energy state at time t₀, ω and b₃ start aligned, which is reflected by ϕ_s = 90^∘. The close alignment for b₃, ω, and L implies that Fig. 11 shows the directions for the principal axis, b₃(t), and for rotational angular momentum, L(t), up to the order of 10⁻² degree. The torques in Sect. 3 cause directional changes on the order of 5^∘ for b₃, ω, and L. Due to the smaller torque values for the bilobed shape, the changes in ω are smaller by an order. Qualitatively, the shapes evolve ω along the line λ_s = −90^∘. According to Sect. 2.1, the equatorial projections of r_s and ω are close to perpendicular. The torque orientations in Figs. 8 and 10 confirm important properties of the evolution of ω. At SSL = 0^∘ for both ellipsoids the jump of torque longitude from −90^∘ to 90^∘ is reflected by corresponding turning points for the ω trajectory in Fig. 11. For shape 67P the ω trajectories (model A and model α) do not contain a turning point. This observation agrees with the smooth curve for longitude, λ_s, in Fig. 10. Longitude, λ_s, is in the range between −90^∘ and −70^∘, which explains the deviation from λ_s = −90^∘ for ω in Fig. 11.

Fig. 12 shows the evolution of ω in time from 400 days before to 400 days after the perihelion passage of comet 67P/C-G. At initial time, the evolution starts with a fixed angular velocity, ω, and a lowest-energy state. Unlike Fig. 11, the solar vector, r_s, describes the trajectory of the sun during apparition 2015 relative to comet 67P/C-G. In the beginning, before the first equinox, SSL is positive and irradiation comes from the north. For both ellipsoids in Fig. 8, we know that the longitude, λ_s, for T̅ is at 90^∘. For shape 67P in Fig. 10, λ_s for T̅ is close to −75^∘. In Fig. 12 these opposite directions are reflected in the ω trajectories. At the two equinox times (95 days before and 224 days after perihelion), the SSL is zero. For the two ellipsoids, these points are related to directional switches between −90^∘ and 90^∘ for λ_s of T̅ in Fig. 8. For the corresponding ω trajectories this is related to turning points. This geometric phenomenon is known and reported; for example, from Whipple & Sekanina (1979), Sekanina (1984), and Gutiérrez & Davidsson (2007). Between the two equinox times, the SSL is negative. This setup, shown in Figs. 8 and 10, leads to torque directions close to λ_s = −90^∘ for all shapes. Figure 12 reflects this by showing similar slopes along all curves around the perihelion passage (marked by a point on the curve). For the bilobed shape, the relative small magnitude of the torque (discussed in Sect. 3) is represented by the small movement of ω. For both ellipsoids and the bilobed shape, the latitude, ϕ_s, of T̅ is always zero, which yields a constant rotation period, as shown in Fig. 12. For the shape 67P, the decrease in |ω| until day −80 (followed by an increase) is directly related to the positive SSL and the corresponding torque direction (with ϕ_s < 0) in Fig. 10.

Fig. 12 Angular velocity, ω. R0, L0 is in lowest-energy state. T is evaluated for solar vector, rs, describing the trajectory of the sun during the apparition 2015 relative to comet 67P/C-G. Considered are the oblate ellipsoid, the near-prolate ellipsoid, the bilobed shape, the shape 67P with model A/fEAF = 1/8, and the shape 67P with model α/fEAF = 1/4. The markers on each curve denote the ω direction at perihelion passage in August 2015. Top: right ascension (ra) and declination (dec) for ω in the equatorial J2000 (EME2000) reference frame. Bottom: angular frequency, |ω|.

4.2 Weighted activity

In Fig. 12 the values f_EAF = 1/8 for modal A and f_EAF = 1/4 for model α are chosen such that the angular frequency, |ω|, changes by the same order of magnitude as the observational data for comet 67P. The corresponding modeled changes for the ω direction (ra and dec) are on the order of 2^∘−3^∘, which is too high by the order of four to six compared to the observation. This property has been reported by Kramer et al. (2019) and Attree et al. (2024). The second observation for the shape 67P in Fig. 12 applies to the slope along the ω trajectory at perihelion (and during the complete pathway). Those curves seem to be rotated between models and observation. A uniform activity on the surface of the shape 67P with either models A or α is not able to bring into agreement the model results and the observations in 2015. This limitation can be resolved by the weighted activity in Eq. (8). With the assumption of a lowest-energy state and a fixed phase angle for R₀ at initial time, an initial vector ω₀ ∈ ℝ³ for angular velocity (three parameters) determines the complete initial condition, R₀ and L₀, for the solution of Eq. (1). Because the activity weights f_EAF,i,j,k in Eq. (8) are eight parameters, a model evolution has n_P = 11 open parameters p_i for i = 1, …, n_P. We fitted the model to the data in Eq. (11) and minimized the quadratic norm, $${\sigma }^2=\sum _{i=1}^{N_{\text{time}}}\sum _{j=1}^3{\left(\frac{{\text{data}}_j(\omega (t_i))-dat{a}_{i,j}}{{\sigma }_{i,j}}\right)}^2.$$(15)

Because the mapping from the parameter vector (p₁, …, p_nP) to the vector (ω(t₁), …, ω(t_Ntime)) is nonlinear, we solved the minimization problem for σ with a standard Gauss-Newton method. The solver attributes negative values to f_{EAF,−1,1,−1} and f_EAF,−1,1,1, so that we set both values to zero. These locations do not agree with low activity of Kramer et al. (2019). This agreement cannot be expected because the patch definitions of Kramer et al. (2019) mix contributions of several regions A_i,j,k from Eq. (6) and optimize the rotation axis under the additional constraint to keep the overall active surface fraction as homogeneous as possible.

We considered four setups with different solutions for the weighted surface activity. These were the bilobed shape and the shape 67P, each in combination with model A and model α. The two ellipsoids were not considered, because they do not have eight nonempty regions A_i,j,k. Fig. 13 shows the resulting activity weights f_EAF,i,j,k. Fig. 14 shows the associated evolution for ω. The norms σ in Eq. (15) are as good as 0.1^∘ for the directional fraction (ra and dec) and as 75 s for the fraction describing the rotation period. Modifications to model configurations produce fits with similar accuracy, σ, and allow us to draw the same conclusions. The first kind of modification is to change the number of elements and, along with it, the shape resolution that we discussed in Sect. 2.1. For the bilobed shape half and double the number of elements N_E change the value of σ in the range of up to 0.01 σ. The second kind refers to the bilobed shape with the modified value α = 1.5 that we introduced in Sect. 2.3. Similarly to the variations in Fig. 13, an agreement on the order of 25% is obtained only for f_EAF,1,1,1. For all remaining regions, the uncertainty is on the order of two or greater. The last kind of modification refers to variations in the sampling rate and uncertainty weights for the data set in Sect. 2.4. A variation in the data incorporating an equidistant time sampling of 7 days and constant uncertainties, σ_i, of 0.05^∘ for direction and 37 s for time period does not change the results.

For the four model setups, the solutions in Fig. 14 solve the mismatch between models and observations in Fig. 12. The general changes for the ω direction and for the angular frequency, |ω|, are qualitatively consistent with the data in Sect. 2.4. For all setups, the weights f_EAF,i,j,k in Fig. 13 are significantly smaller than for uniform activity, on the order of three for model A and on the order of five for model α. Because f_EAF is a unit free quantity, different values f_EAF,i,j,k for model A and model α go back to the normalization constants for both TPMs in Sect. 2.3. All values of the effective active fraction are bounded by f_EAF ≤ 0.085. In the regions A_{−1,−1,−1}, A_1,±1,1 ⊂ A_SH and A_1,±1,−1, A_−1,−1,1 ⊂ A_NH all setups agree to attribute nonzero activity. The uncertainties for the weights f_EAF,i,j,k suffer from significant limitations. The only region in which all the setups agree is A_1,1,1 ⊂ A_SH with weights in the range 0.06 ≤ f_EAF,1,1,1 ≤ 0.08. This is related to the fact that the data in Sect. 2.4 describe the dynamics for negative SSL (during the perihelion passage of comet 67P/C-G in 2015) predominantly and thus for activity on A_SH. In Fig. 6 for the shape 67P, A_1,1,1 is located in the regions (0^∘, 60^∘) × (−90^∘,−45^∘) and (150^∘, 210^∘) ×(−90^∘,−30^∘) for λ and ϕ. This fits well with the location of increased activity in Fig. 10 of Kramer et al. (2019). The same authors report a decrease in the effective active fraction in the regions (−45^∘, 15^∘) × (−30^∘, 0^∘) and (150^∘, 210^∘) × (−15^∘, 15^∘) for λ and ϕ. This statement correlates to reduced weights in Fig. 13 for the regions A_{−1,−1,−1} and A_1,−1,1 in Fig. 6. Restricted to model A, some agreement can be obtained between the two shapes on the regions A_{−1,−1,−1}, A_1,1,−1, and A_−1,−1,1. On the remaining regions the weights contain much higher uncertainties based on the consideration of all different setups.

One major difference between the evolutions of uniform (Fig. 12) and weighted (Fig. 14) activity is the rotated trajectory for the two ω directions. For uniform activity, Fig. 11 shows this trajectory (with respect to the solar basis) close to the line λ_s = −90^∘. For weighted activity, Fig. 15 shows the same trajectory close to the line λ_s = 0^∘. Especially for negative SSLs (related to the early part of the evolution in Fig. 15), weighted activity causes this rotation for the ω trajectory (and for the slope along the trajectory) by approximately 90^∘. Negative SSLs represent predominant illumination conditions during the perihelion passage of comet 67P/C-G in 2015. Due to rotation states close to lowest energy, this directional change applies to L as well. For weighted activity, the torque, T̅, points mainly in the λ_s = 0^∘ direction that expresses T̅ · b₁>0. Thus this directional evolution in Fig. 15 confirms the increased weights for region A_1,1,1 in Fig. 13.

Fig. 13 Activity weights fEAF,i,j,k for regions Ai,j,k in Eq. (8) for the bilobed shape and the shape 67P with model A and model α. The values of fEAF,i,j,k are constrained by the ω data in Sect. 2.4.

Fig. 14 Angular velocity, ω, similar to Fig. 12 for the bilobed shape and the shape 67P with model A and model α. Weighted activity from Fig. 13. Top: right ascension (ra) and declination (dec) for ω. Bottom: angular frequency, |ω|.

Fig. 15 Longitude, λs, and latitude, ϕs, for angular velocity, ω. The same projection with respect to the fixed solar basis and the same rs as in Fig. 11. Considered are the bilobed shape with model A, the shape 67P with model A, and the shape 67P with model α, and weighted activity fEAF in Fig. 13.

5 Conclusions

The vectorial torque efficiency introduced in Eq. (5) can be evaluated for any general cometary body. It provides a generalization of the scalar torque efficiency in Keller et al. (2015b) and Attree et al. (2019) that controls the rotational period alone. The signs of all three components of the vectorial torque efficiency lead to the specification of a domain decomposition into a maximum of eight distinct regions on the surface. In particular, this applies to the bilobed shape as well as to the body shape of comet 67P/C-G. Geometric symmetry for idealized cases can reduce the number of regions. The oblate ellipsoid features zero and the near-prolate ellipsoid features four nonempty regions. The domain decomposition is divided into two subsets containing four regions each. One subset is linked to the southern hemisphere, A_SH, and the other is linked to the northern hemisphere, A_NH.

For comet 67P/C-G in the time between 2014 and 2016, uniform activity leads to a mismatch between model and observations for directions of the rotation state. This mismatch is resolved with a best-fit solution for effective active fractions. This directional switch by approximately 90^∘ is illustrated in Fig. 14 for angular velocity and in Fig. 15 for torque. The fit quality for the data is on the order of 0.1^∘ for the ω direction and 75 s for the rotation period. Only one region A_1,1,1 is constrained by the data with values at 0.06 ≤ f_EAF,1,1,1 ≤ 0.08. It is located on the southern hemisphere, A_SH, which is the hemisphere that has a high solar irradiation (and negative SSL) around the perihelion passage of comet 67P/C-G in 2015. Thus, torque activity at that time is high for A_1,1,1 and according to the construction of A_1,1,1 the torque direction has three positive components with respect to the solar basis. In particular, this produces a positive b₃ component for the torque, an increase in the angular frequency, |ω|, and a decrease in the rotation period at that time. For the weights f_EAF,i,j,k of all other regions with (i, j, k) ≠ (1,1,1), the nonlinear fits hold high uncertainties. In particular, this affects all regions of the northern hemisphere, A_NH. Low solar irradiation around the perihelion passage results in low dynamical activity and weak signals for these locations. In addition, the discussion of uncertainties remains unchanged once the shape resolution of Δ x ≈ 150 m is moderately decreased or increased. We consider two types of TPM, model A and model α. The qualitative discussion does not differ for both model types and for the values α = 1.5 and α = 2. In general, the temporal sampling and the uncertainties prescribed for the ω data affect the fit. Our discussion does not change for a data set with modified sampling and uncertainty weights.

The same control for the rotation state is possible for each combination between the two shapes, the bilobed shape and the shape 67P, and the two TPMs, model A and model α. This ability is a characteristic of the domain decomposition. For lowest-energy states and sublimation following an instantaneous TPM, each region is related to torque formation pointing into one of the eight spatial quadrants with respect to the solar basis. Each region collects one specific combination of (signs for) the vectorial torque efficiency. Thus, the weights of the effective active fraction control all torque directions. Eight is the minimum number for weights (and regions), and it is obtained by our approach.

On the surface of comet 67P/C-G for a wide range of locations, the weights for effective active fraction are not well constrained by rotation data. This indetermination stems from the specifics of both, the geometry of the nucleus and the activity model. We have shown this comparing results for two shapes and two activity models. A better localization of the activity must consider complementary data. Attree et al. (2023) shows the option to add NGA to form a joint data set. Globally integrated gas production such as in Läuter et al. (2020) and locally resolved gas activity in Läuter et al. (2022) can also be included.

Any approach for a joint data set must address the difficulty of a relative weighting of the different data sets, which introduces a bias to the final result. With knowledge of the shape geometry, the localization of surface activity can be applied to other solar system objects that have weaker sources of torque activity. Nongravitational momenta have this potential on dark comets reported by Seligman et al. (2023) as well as the YORP effect on asteroids reviewed by Bottke et al. (2015).

Acknowledgements

The authors gratefully acknowledge the computing time made available to them on the high-performance computer Lise at the NHR center NHR@ZIB. This center is jointly supported by the federal ministry of education and research and the state governments participating in the NHR. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #547 (Understanding the Activity of Comets Through 67P’s Dynamics) with N. Attree acting as PI. M.L. acknowledges the support by Johannes Kepler Universität Linz Austria for a research visit in 2024. We thank the referee for the helpful comments, which improved the presentation of the manuscript.

Appendix A Geometric properties

For the shape model of comet 67P given by Preusker et al. (2017), Kramer et al. (2019) have derived the tensor of inertia. Based on a strictly homogeneous density distribution within the body their tensor suffers from an incompatibility of the local basis vector, e₃, and the observed angular velocity vector. To circumvent this issue they add inhomogeneities to the density and obtain an adapted tensor of inertia, $$\left(\begin{array}{c}9.3408457\times {10}^{18}5.6695663\times {10}^{16}0\\ 5.6695663\times {10}^{16}1.6562414\times {10}^{19}0\\ 001.8192083\times {10}^{19}\end{array}\right)\mathrm{k}\mathrm{g}{\mathrm{m}}^2,$$(A.1)

which is applied to Sect. 2.1.

At position x ∈ ∂Ω_t, the torque density, t_NG, in Eq. (10) has the direction of the term x × v. Two components with respect to the local basis, B(v), in Fig. 1 satisfy the relation $$(x\times \nu )\cdot b_1(\nu )=\frac{-\nu \cdot b_3}{\sqrt{1-(\nu \cdot b_3{)}^2}}(x\times \nu )\cdot b_3.$$

This transfers to the signs of components for t_NG. $$\begin{array}{rl}\text{sgn}(t_{\text{NG}}\cdot b_1(\nu ))& =\text{sgn}(t_{\text{NG}}\cdot b_3)\text{\hspace{0.17em}\hspace{0.17em}on\hspace{0.17em}\hspace{0.17em}}\mathrm{\partial }{\mathrm{\Omega }}_{\text{SH}},\\ \text{sgn}(t_{\text{NG}}\cdot b_1(\nu ))& =-\text{sgn}(t_{\text{NG}}\cdot b_3)\text{\hspace{0.17em}\hspace{0.17em}on\hspace{0.17em}\hspace{0.17em}}\mathrm{\partial }{\mathrm{\Omega }}_{\text{NH}}\end{array}$$(A.2)

Appendix B Fourier modes for torque

Fourier decompositions have been evaluated for NGA by Kramer & Läuter (2019) and for NGT by Kramer et al. (2019). Both approaches apply Fourier modes for vector quantities, which then simplify in the rotating frame. They share the same scalar field for mass generation. For the present section we apply the Fourier decomposition to this scalar field only. We consider specific assumptions for the rotation state. First the solution, R and L, of Eq. (1) describes a principal axis rotation around the short axis b₃ with constant |ω| and the rotation period Δ t = 2π/|ω|. With elementary matrices R_X and R_Z in Montenbruck & Gill (2000) in this situation, R can be represented with Euler angles by $$R(s)=R_{\text{z}}(-{\varphi }_0)R_{\text{x}}(-{\theta }_0)R_{\text{z}}(-{\psi }_0-|\omega |s)$$

at time s with appropriate angles ϕ₀, θ₀, ψ₀. In Eq. (10) torque, T, depends on the vector r_s pointing toward the sun. The second assumption is to consider a constant r_s. This leads to a constant solar basis, B_s = B(r_s), in Eq. (4). Because the orthogonal matrices B_s and R(s) agree in the third column, we can find a time t₀ to bring into agreement B_s = R(t₀). With λ(s) = |ω|(s−t₀), this leads to $B_{\text{s}}^T$ R(s) = R_z(−λ(s)).

At position ξ ∈ ∂Ω with the outward normal ν_ξ in the body frame, the vector field for torque in Eq. (10) can be expressed by t_NG(x(s), s) = B_s t_BF($B_{\text{s}}^T$ x(s)) with x(s) = R(s) ξ and the torque in the body frame, $$t_{\text{BF}}(\xi )=-{\dot{m}}_{\text{BF}}({\nu }_{\xi })\xi \times {\nu }_{\xi },{\dot{m}}_{\text{BF}}({\nu }_{\xi })=v_{\text{gas}}\dot{m}(I(B_{\text{s}}{\nu }_{\xi },r_{\text{s}})).$$

For the time average in Eq. (12), we obtain $$\begin{array}{rl}\overline{t_{\text{NG}}}(x(t),t)& =\frac{B_{\text{s}}}{\mathrm{\Delta }t}\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}{t}_{\text{BF}}(B_{\text{s}}^{\mathrm{T}}x(s))ds\\ & =\frac{B_{\text{s}}}{2\pi }\underset{0}{\overset{2\pi }{\int }}{t}_{\text{BF}}(R_{\text{z}}(-\lambda )\xi )d\lambda .\end{array}$$(B.1)

The torque t_BF(R_z(−λ) ξ) can be separated into a product of the scalar term f(λ) = ṁ_BF(R_z(−λ) v_ξ) and the vectorial component −R_z(−λ)(ξ × ν_ξ). During one rotation the angle of maximum momentum is related to the angle $${\lambda }_{\text{M}}=\underset{0\le \lambda \le 2\pi }{\text{argmax}}f(\lambda ).$$

The value of λ_M is specific for ξ and we denote the corresponding time t_M holding the property λ(t_M) = λ_M. Because of monotonicity for ṁ together with irradiation, I, in Eq. (9), t_M can be expressed by argmax_s(r_s · R(s) v_ξ), equivalently. We observe two properties for this maximum. For the local basis in Eq. (4) the maximization property of t_M yields the relation B(v(t_M)) = B_s, which gives $$R_{\text{z}}(-{\lambda }_{\text{M}})(\xi \times {\nu }_{\xi })=B^{\mathrm{T}}(\nu (t))(x(t)\times \nu (t)).$$(B.2)

Secondly f(λ) is a symmetric function with respect to λ_M. That is why function f(λ) has vanishing coefficients related to sin functions and it has the Fourier coefficients $$c_n^{\prime }(\xi ,r_{\text{s}})=\frac{1}{\pi }\underset{0}{\overset{2\pi }{\int }}\cos (n(\lambda -{\lambda }_{\text{M}}))f(\lambda )d\lambda .$$

Following Chesley & Yeomans (2004) for spherical geometry, irradiation, I(B_s R_z(−λ) v_ξ, r_s), can be expressed representing the diurnal trend at the position ξ. In general f(λ) is nonlinear and it reflects the properties of the TPM in Sect. 2.3. The coefficients c_n^′ are the result of numerical integration resulting in $$\begin{array}{rl}& f(\lambda )=\frac{c_0^{\prime }}{2}+\sum _{n=1}^{\mathrm{\infty }}{c}_n^{\prime }\cos (n(\lambda -{\lambda }_{\text{M}})),\\ & \frac{1}{\pi }\underset{0}{\overset{2\pi }{\int }}f(\lambda )R_{\text{z}}(-\lambda )d\lambda =\text{diag}(c_1^{\prime },c_1^{\prime },c_0^{\prime })R_{\text{z}}(-{\lambda }_{\text{M}}).\end{array}$$

Together with Eq. (B.2) this can be used to evaluate Eq. (B.1), which gives $$b_i(r_{\text{s}})\cdot \overline{t_{\text{NG}}}(x,t)=-c_i{b}_i(\nu )\cdot (x\times \nu )$$(B.3)

with c₁ = c₂ = c₁^′/2 and c₃ = c₀^′/2. As a consequence the torque average, $\overline{t_{\text{NG}}}(x(t),t)$, is specific for the position ξ, and it is constant in time t.

For example at position ξ ∈ ∂Ω, the specific choice $$f(\lambda )=\{\begin{array}{ll}{\dot{m}}_0& \text{for}|\lambda -{\lambda }_M|\le \pi /2\\ 0& \text{otherwise}\end{array}$$

represents a constant activity along the rotating position x(t) with some constant ṁ₀ = ṁ₀(ξ). This idealized assumption yields the result c₁ = c₂ = ṁ₀/π, c₃ = ṁ₀/2.

References

All Figures

	Fig. 1 Local basis, B(v) = (b1(v), b2(v), b3), for the choice y = v(x, t) in Eq. (4) at time t and at position x ∈ ∂Ωt with the outward normal v.
In the text

	Fig. 2 Solar basis, B(rs) = (b1(rs), b2(rs), b3), for the choice y = rs(t) in Eq. (4) at time t.
In the text

	Fig. 3 For the near-prolate ellipsoid, the values (x × v) · bi(v) (vectorial torque efficiency) are shown based on longitude, λ, and latitude, ϕ, in the body frame. bi(v) is the local basis vector in Fig. 1.
In the text

	Fig. 4 For the bilobed shape (top) and for shape 67P (bottom), the values (x × v) · b2(v) (second component of the vectorial torque efficiency) are shown based on longitude, λ, and latitude, ϕ, in the body frame. b2(v) is the local basis vector in Fig. 1.
In the text

	Fig. 5 For the near-prolate ellipsoid decomposition of ∂Ωt into the four regions Ai, j, k in Eq. (6). Region A0 is colored in white. Top: longitude, λ, and latitude, ϕ, in the body frame. Bottom: 3D projection using the same color code.
In the text

	Fig. 6 For the bilobed shape (top) and for the shape 67P (bottom), decomposition of ∂Ωt into the eight regions Ai,j,k in Eq. (6). Region A0 is colored in white. Longitude, λ, and latitude, ϕ, in the body frame.
In the text

	Fig. 7 For the bilobed shape (top) and for the shape 67P (bottom), decomposition of ∂Ωt into the eight regions. The same data as in Fig. 6 but in 3D projection.
In the text

Fig. 8 Direction of torque, T̅, for the near-prolate ellipsoid. Torque, T̅, in Eq. (12) is shown with model A and fEAF = χAi,j,k for region Ai,j,k in Eq. (6). Uniform activity fEAF = 1 is shown with dotted black lines for the oblate ellipsoid and for the near-prolate ellipsoid. λs and ϕs denote longitude and latitude for T̅ with respect to the solar basis in Fig. 2. The solar vector, rs, changes SSL from −45∘ to 45∘.

In the text

	Fig. 9 Direction of torque, T̅, for the bilobed shape. Torque, T̅, in Eq. (12) is shown with model A and fEAF = χAi,j,k for region Ai, j, k in Eq. (6). Uniform activity fEAF = 1 is shown with dotted black lines. λs and ϕs denote longitude and latitude for T̅ with respect to the solar basis in Fig. 2. The solar vector, rs, changes SSL from −45∘ to 45∘.
In the text

	Fig. 10 Direction of torque, T̅, for the shape 67P. Torque, T̅, in Eq. (12) is shown with model A and fEAF = χAi,j,k for region Ai,j,k in Eq. (6). Uniform activity fEAF = 1 is shown with dotted black lines. λs and ϕs denote longitude and latitude for T̅ with respect to the solar basis in Fig. 2. The solar vector, rs, changes SSL from −45∘ to 45∘.
In the text

Fig. 11 Longitude, λs, and latitude, ϕs, for angular velocity, ω, with respect to the fixed solar basis in Fig. 2, based on rs(t0) and b3 = R(t0) e3 at initial time, t0 . T is evaluated with |rs| = 2 au and changing SSL from −45∘ to 45∘. Considered are the oblate ellipsoid, the near-prolate ellipsoid, the bilobed shape, the shape 67P with model A/fEAF = 1/8, and the shape 67P with model α/fEAF = 1/4.

In the text

Fig. 12 Angular velocity, ω. R0, L0 is in lowest-energy state. T is evaluated for solar vector, rs, describing the trajectory of the sun during the apparition 2015 relative to comet 67P/C-G. Considered are the oblate ellipsoid, the near-prolate ellipsoid, the bilobed shape, the shape 67P with model A/fEAF = 1/8, and the shape 67P with model α/fEAF = 1/4. The markers on each curve denote the ω direction at perihelion passage in August 2015. Top: right ascension (ra) and declination (dec) for ω in the equatorial J2000 (EME2000) reference frame. Bottom: angular frequency, |ω|.

In the text

	Fig. 13 Activity weights fEAF,i,j,k for regions Ai,j,k in Eq. (8) for the bilobed shape and the shape 67P with model A and model α. The values of fEAF,i,j,k are constrained by the ω data in Sect. 2.4.
In the text

	Fig. 14 Angular velocity, ω, similar to Fig. 12 for the bilobed shape and the shape 67P with model A and model α. Weighted activity from Fig. 13. Top: right ascension (ra) and declination (dec) for ω. Bottom: angular frequency, |ω|.
In the text

	Fig. 15 Longitude, λs, and latitude, ϕs, for angular velocity, ω. The same projection with respect to the fixed solar basis and the same rs as in Fig. 11. Considered are the bilobed shape with model A, the shape 67P with model A, and the shape 67P with model α, and weighted activity fEAF in Fig. 13.
In the text