© ESO, 2016

1. Introduction

Immediately after the wake-up of ROSETTA and the recommissioning of the scientific instruments, the OSIRIS cameras (Keller et al. 2007) started to perform systematic observations of comet 67P/Churyumov-Gerasimenko (67P). One of the goals of these observations was to accurately determine the rotational parameters of 67P and establish a solid starting point for detecting possible variations. By using data from March to June, 2014, Mottola et al. (2014) analyzed these early photometric observations and obtained that 67P was rotating at this epoch with a period of 12.4043 ± 0.0007 h.Near-nucleus operations began in early August 2014, with a comet-characterization phase followed by a mapping phase. One of the main goals during those two phases was to reconstruct detailed topographic models of the nucleus of 67P and to retrieve its global shape by using the hundreds of images that the OSIRIS instrument acquired.

So far, two different methods have been applied to retrieve high-resolution shape models of 67P from OSIRIS images: stereo-photogrammetry (SPG; Preusker et al. 2015) and stereo-photoclinometry (SPC; Jorda et al. 2016). The reconstruction of the shape of 67P involves a very accurate determination using stereo landmarks of the comet position and attitude in camera frame. Combined with the knowledge of the S/C position and attitude reconstructed by the Rosetta Mission Operations Centre (RMOC) from the radio ranging and star tracker data, this provides a very accurate determination of the comet position and attitude in an inertial frame, usually the EME2000 reference frame. Preusker et al. (2015) described the reconstruction of the so-called SHAP4S model using the SPG technique. They showed that the accuracy of the retrieved topographic model significantly improves when a precessing spin axis is introduced. By defining seven different image blocks, Preusker et al. (2015) obtained that the right ascension and declination (RA/Dec) of the spin pole that allow the most accurate shape model describe a circle around (RA, Dec) = (69.54^° ± 0.1^◦, 64.11^° ± 0.05^°), with a radius of 0.14^° ± 0.03^° and a possible period of 256.8 h ± 12 h (10.7 d ± 0.5 d), see Fig. 2 in Preusker et al. (2015). This was interpreted by the authors as the precession of the spin axis around the angular momentum vector. It is the first evidence that comet 67P might be rotating in complex mode, although with a low excitation given the small angle between the spin axis and the presumable angular momentum vector.

Fig. 1 SPC shape of 67P/Churyumov-Gerasimenko (SHAP5, version 1.2, Jorda et al. 2016) in the so-called Cheops frame (see Scholten et al. 2015). The spin axis is the z-axis.

The results obtained by Jorda et al. (2016) with the SPC method are consistent with those of Preusker et al. (2015). By applying the SPC method, Jorda et al. (2016) retrieved a spin axis that moves around (RA, Dec) = (69.57^◦, +64.01^◦). Jorda et al. (2016) obtained that the spin axis does not describe a circumference, but approximately fills an ellipse in an isotropic plot (Fig. 2). This difference might arise because the data from SPC method cover a time span muchlonger than those of the SPG method. When SPC and SPG results are compared, theputative coordinates of the angular momentum vector are also slightly different byabout 0.1^◦. This small difference, slightly larger than the associated errors, still needs to be understood (Jorda et al. 2016). Interestingly, Jorda et al. (2016) analyzed the spin axis orientation (the body frame z-axis coordinates) by means of the phase-dispersion minimization technique (Stellingwerf 1978) and obtained a confident periodicity of 276 h ± 12 h (11.5d ± 0.5d) from separately considering the RA and Dec coordinates. This periodicity can be considered to agree within the error bars with the early determination by Preusker et al. (2015). Additionally, among other characteristics of 67P, Jorda et al. (2016) determined the pre-perihelion rotational period and its evolution with time, obtaining that up to December 2014 it was 12.4041 h ± 0.0001 h, which refines the earlier determination of Mottola et al. (2014) that was based on light curves. This means that the SPG and SPC results both point to a situation in which the motion of the spin pole is perturbed by an additional periodic motion with a period of ∼270 h, much longer than the spin period of 12.40 h measured by Mottola et al. (2014). This seems to confirm the complex nature of the rotation of 67P, but an interpretation of these measurements is still lacking.

The goal of this paper is to interpret the periodicity associated with the spin axis orientation measured by ROSETTA taking into account the pre-perihelion rotation period measured by Mottola et al. (2014) that was later refined by Jorda et al. (2016). We restrict our analysis to the observations obtained until the end of December 2014 (250 days before perihelion), during which the rotation period remained constant (within the accuracy of the OSIRIS and flight dynamics measurements). After this, the rotation period started to vary (Jorda et al. 2016; Rosetta Mission Operations Reports), indicating that an outgassing-induced torque started to have detectable effects on the angular momentum of comet 67P. An activity-driven torque has been simulated (see, e.g., Bertaux 2015; Keller et al. 2015) to explain the spin period evolution of 67P. These simulations and calculations suggest that the sublimation-induced torque is negligible at large heliocentric distances and does not affect the rotational evolution. We therefore assume in the following that 67P was rotating under torque-free conditions before December 2014.

Finally, the moments of inertia can be derived from the geometric global shape of 67P assuming it has a homogeneous density. Considering the latest SHAP5 model (SPC method, Jorda et al. 2016), the principal inertia moments would follow the relationship I_x:I_y:I_z = 1.00:1.83:1.98, where I_x, I_y, and I_z are the inertia moments along the largest, intermediate, and shortest principal axes¹, respectively.

2. Data

Fig. 2 Coordinates in J2000 of the spin axis of 67P for the time interval between –376 and –250 days from perihelion. These coordinates are obtained during the SPC procedure employed to retrieve the shape of 67P (Jorda et al. 2016).

Two very different data sets of the spin axis orientation are available: those obtained with the SPG method on the one hand, and those obtained with the SPC on the other hand. The first data set includes seven determinations considering time intervals (blocks) of 12–13 h and covering a time span of about one month. The SPC method yielded 232 estimates of the spin axis orientation that cover a time span of approximately 125 days (from August 2014 to December 2014), each of them defined for time intervals of 10 h. This means that the coordinates do not correspond to the instantaneous position, but rather to its average location within a time interval of 10 h.

The image blocks (set of images allowing the retrieval of the average spin axis orientation for the covered time interval) in the two data sets are slightly different, therefore we considered only the data set derived from SPC to simplify the procedures. This was done with the understanding that data from SPG are compatible with those of SPC and that the addition of the seven data points from SPG to the 232 spin axis orientations from SPC will not change the conclusions of the study. The RA/Dec coordinates of the body frame z-axis as retrieved by means of the SPC method considered here are shown in Fig. 2. The typical formal angular error associated with RA and Dec coordinates is on the order of 0.03^◦ (Jorda et al. 2016). Because our intention is to perform statistical tests by comparing observational and theoretical distributions, we decided to remove seven clearly identifiable outliers from these coordinates and kept 225 coordinates. This data set is called observationally derived data.

In Fig. 2 we show that the spin axis orientation is not constant, but varies in time around a central position. Relying on the accuracy of the SPC method, we consider that this might be a clear indication that the body is slightly excited from the rotational point of view, with the spin axis moving around the constant angular momentum orientation. When estimated as the mean value of the spin axis coordinates, the RA and Dec of the angular momentum vector would be around 69.5^◦ and 64.04^◦, respectively. These coordinates are slightly different from those reported by Jorda et al. (2016) (and slightly closer to those of Preusker et al. 2015). The small difference arises because outliers were removed from our data set.

If RA and Dec coordinates, considered separately, are analyzed by means of different time-series techniques, a very significant peak at periodicities around 270 h appears. Figure 3 shows the Lomb² periodograms (Lomb 1976) of the RA and Dec data displayed in Fig. 2. The Lomb periodogram shows the highest significance of RA data for a periodicity of approximately 264 h, while Dec data have a significant associated periodicity at approximately 270 h, both compatible with the periodicity of 276 h obtained by Jorda et al. (2016) with the PDM technique. In addition to this periodicity, it is known that the light curve has a periodicity of 12.40 h (Mottola et al. 2014), which was interpreted as the spin period, assuming that 67P was rotating in pure spin mode at that time. In principle, the 12.40 h might well correspond to the modulus of full rotational vector, Ω, the rotational period. Determinations of the omega vector modulus obtained with the SPC method (Jorda et al. 2016) point out that this interpretation is probably correct. Nevertheless, this still needs some clarification and is discussed below. In addition, the meaning of the 270 h periodicity found in the spin axis RA/Dec coordinates needs to be determined, as well as the information that might be extracted from it. The 270 h periodicity is called RA/Dec periodicity throughout.

Fig. 3 Lomb periodograms of the RA and Dec coordinates of the spin axis displayed in Fig. 2. A significant peak (with a significance higher than 99.99%) is located at a periodicity of approximately 270 h.This result is fully compatible with those found by Jorda et al. (2016) with the PDM technique. Data phased to this periodicity are shown in Fig. 12 of Jorda et al. (2016).

3. Some considerations from Euler equations

Misunderstandings may arise when the rotation of a body in space is described with different references frames. Thereference frame of the motion needs to be specified and the transformation between body frame and an inertial reference frame has to be defined. We here use the following terminology. When the motion of the body is described as the rotation around the axes of the body frame, the term “intrinsic” must be used to define the rotational motion. When the rotation is described as the motion of the inertial frame around a fixed body frame, the rotations are considered as “extrinsic”. We define a body frame with the orthogonal axes along the principal axes of inertia. The x-axis has the direction of the principal axis associated to the lowest moment of inertia (I_x). The y-axis has the direction of the principal axis of inertia associated to the intermediate inertia moment (I_y). Finally, thez-axis is the orthogonal axis parallel to the principal axis associated with the hightest inertia moment (I_z). To describe the rotation of this xyz body frame (and therefore that of the body, which is considered to be rigid) in space, a convenient set of Euler angles must be selected. The Euler angles define the transformation matrix between the body frame (xyz) and the inertial frame (XYZ). By considering the standard ZXZ convention, the transformation between the body frame xyz and the inertial frame XYZ, and therefore the motion of the body in space, is defined by the standardEuler angles φ, θ, and ψ.These three angles and the corresponding associated velocities define a rotation of the body frame with regard to the inertial frame and have specific meanings. In the terminology we use, the ψ angle represents the intrinsic rotation of the body, which is a spin around the principal axis z that is associated with the highest inertia moment. The θ angle describes a rotation around the line of nodes defined by the two aforementioned reference frames. This motion is called nutation and can be described as an oscillation of the body frame z-axis with regard to the inertial Z-axis. Finally, φ, the angle between the line the nodes and the inertial X-axis, defines a rotation around the inertial Z-axis. If the angular momentum orientation is defined along the inertial Z-axis, this latter rotational motion is generally described as the precession of the body frame z-axis (i.e., the axis defining the intrinsic rotation of the body) around the angular momentum orientation. These reference frames, angles, and their corresponding velocities, are shown in Fig. 4. We here describe the rotational motion of the body as follows:

• φ and describe the precession of the axis with the highest inertiamoment (z) around the angular momentumdirection (Z).

• θ and describe the nutation of the axis with the highest inertia moment (z) with regard to the angular momentum direction (Z).

• ψ and represent the intrinsic rotation, spinning, of the body around the axis with the highest inertia moment (z).

• Ω is the total rotation of the body, formed by the three previous motions together.

We therefore reserve the use of the term spin to refer to the intrinsic rotation described by the angle ψ, and spin axis or spin pole for the body frame z-axis. According to this, Fig. 2 shows the RA/Dec coordinates of the spin pole. The term rotation is generally used to describe the combination of all three described motions. Therefore, rotation axis refers to the instantaneous rotation vector driven by the spin, nutation, and precession motions together. The rotation is represented as Ω_XYZ or Ω_xyz, depending on the reference frame in which is expressed. We note that because Ω = | Ω_XYZ | = | Ω_xyz |, the rotational period P = 2π/ Ω does not need the reference frame to be specified.

Fig. 4 Body frame (xyz) and inertial frame (XYZ) with the Euler angles and their velocities in the standard ZXZ convention. This figure is adapted from Fig. A.1 from Samarasinha & Mueller (2015). Reproduced with permission.

To continue with notations, it is true throughout this study that I_x<I_y<I_z. If we call E the rotational energy, L the modulus of the angular momentum vector, the quantity C = L²/ 2E fulfills I_x<C<I_z in torque-free. We therefore may call the quantity E_l = 1−C/I_z the excitation level (E_l), whose magnitude is 0 when the body rotates at the lowest energy for its angular momentum (i.e., a simple rotation around z-axis). When I_x<C<I_y, the body rotates in long-axis mode (LAM). When I_y<C<I_z is fulfilled, the body rotates in short-axis mode (SAM; see, e.g., Belton 1990). 67P is presumably only slightly excited (given the small amplitude of the very likely complex motion), the rotational excitation must therefore be comparatively low (E_l → 0, C → I_z) and the nucleus probably rotates in SAM, close to its lowest rotational energy for its angular momentum. In SAM, the body spins or intrinsically rotates around the principal axis with the highest inertia moment and at the same time performs both a precession of this principal axis around the angular momentum vector and a nutation, oscillation, also of this principal axis. This means that the previously defined Euler angles are perfectly suitable to describe the possible rotational motion of 67P in SAM.

Under the assumption of torque-free motion, the body must rotate following certain constraints derived from the Euler equations, that is, the Euler angles and their velocities follow some regular behavior depending on inertia moments and excitation level. More details on the description of excited rotational states are given by Samarasinha & Mueller (2015).The authors included a complete and comprehensive compendium of the equations and their dependences governing the complex motion of a torque-free body in SAM mode. In the following, we summarize the more relevant equations and conclusions for the present study.

Thus, in torque-free and in SAM,

• ψevolves in a periodic way with a period P_ψ(1)where (2)Even though this is periodic, it does not mean that is constant.

• is given by (3)which is periodic with a period P_ψ/ 2. In the previous expression, snτ is the Jacobian elliptic function of the argument τdefined as (4)and t denotes time.

• φis generally aperiodic, except for some particular cases. Nevertheless, a time-averaged period for φ can be defined as (5)

• For slightly excited cases, the two above periods fulfill the condition (6)

• φ and ψ circulate in opposite senses. The component of the total rotational velocity vector, Ω_xyz, along the principal axis of inertia with the highest inertia moment is given by .

• θoscillates with a period P_ψ/ 2 between a minimum (7)and a maximum (8)

As two periodicities have been detected for 67P, the 12.40 h from Mottola et al. (2014) and the RA/Dec periodicity of about 270 h, the first and most intuitive interpretation is to consider that the body spins with a period of P_ψ = 12.40 h while the body frame z-axis precesses around the angular momentum vector with a precession period P_φ ≈ 270 h. Nevertheless, in the standard terminology, this description is not supported by the Euler equations, in which, considering the aforementioned constraints, P_ψ, the spin period, must always be longer than the precession period, P_φ, particularly P_ψ> 2P_φ. Therefore, the RA/Dec periodicity can only be identified with P_ψ, that is, with the intrinsic rotation period or spin, the rotation of the body around the principal axis with the highest inertia moment. Several combinations of these motions may result in the two observationally detected periodicities. The first problem is to correctly identify the meaning of the 12.40 h periodicity detected by Mottola et al. (2014) that was confirmed by Jorda et al. (2016). If we assume that 12.40 h is the rotational period, P = 2π/ Ω, and that the excitation is comparatively low, we can verify that the precession period, P_φ, must be shorter than P. Thus, if P = 12.40 h, P_φ< 12.40 h, and the 270 h periodicity can be identified with the spin period, P_ψ, or any other possible combination of periodicities. We call thiscase 1, and it is characterized by P = 12.40 h and P_φ< 12.40 h. An alternative to this possibility arises if the periodicity detected by Mottola et al. (2014) that was confirmed by Jorda et al. (2016) is not the rotation period. From our own experience and results from simulations of light curves of bodies rotating in complex mode, coinciding with those of Samarasinha & Mueller (2015), the most commonly detected periodicity in simulated light curves of bodies rotating in SAM corresponds to the precessional motion, that is, P_φ, not to the total rotational period P = 2π/ Ω. If that were our case, the periodicity detected by Mottola et al. (2014), 12.4 h, might be P_φ. In this circumstance, as previously, the 270 h might correspond to the spin period, P_ψ, or to any combination of motions. We call this case 2, and it is characterized by P_φ = 12.40 h and P> 12.40 h. After considering these two cases, it is still necessary to point out that from series of OSIRIS images that cover a full period of 12.40 h, the body seems to rotate with that period, meaning that the orientation of the body in space approximately returns to its starting orientation after the 12.40 h period. This clearly favors case 1, but we still kept case 2, which is discarded if P, depending on the value of P_ψ, is very different from 12.40 h.

4. Theoretical simulations

To study the behavior of the spin axis in an inertial frame, simulations with different inertia moments and excitation levels were performed. We calculated Euler angles and their velocities and transformed the instantaneous orientation of the spin axis to J2000 to compare it with the data shown in Fig. 2. Additionally, simulated RA/Dec coordinates were analyzed by means of the Lomb technique to determine which periodicities would be detected. Noise and averaging of theoretical data were considered to simulate actual circumstances as best possible. In the simulations, the real dependence of Euler angles with all involved variables, as described in Samarasinha & Mueller (2015), for instance, was considered. The only simplification was that to be able to perform the calculations in a systematic and easy way, we considered that the angular momentum modulus and the rotational energy can initially be defined as L = I_z·Ω and , respectively. This approach assumes that Ω ≈ Ω_z, neglecting the appropriate contribution to Ω_xyz of the components Ω_x and Ω_y. This approach is a reasonable sacrifice as starting point given the expected low excitation level, in which E_l ≈ 0, and C/I_z ≈ 1. Evaluated a posteriori, the highest excitation level necessary to interpret the observationally derived data is on the order of 10^-4. In case 1 we considered that the rotation velocity is Ω = 2π/P = 2π/ (12.4 ∗ 3600) rad s^-1, with P the periodicity detected by Mottola et al. (2014). For case 2, in which the precession period P_φ is identified with the periodicity detected by Mottola et al. (2014), the rotational period P, and therefore the rotational energy and the modulus of the angular momentum, depend on the meaning assigned to the RA/Dec periodicity. Thus, in case 2, P is explicitly discussed below.

Fig. 5 a) Bullets (•) show the spin axis orientation as determined with the SPC method (Fig. 2). Small dots (·) are the spin vector orientation for a simulated body with inertia moments equal to those derived from the shape of 67P assuming homogeneity, and an El = 1.9 × 10-6 under case 1 circumstances. These data are displayed with a time separation of 1 h. Diamonds (◆), to be compared with the bullet distribution, correspond to the orientation of the spin axis for the simulation performed at the time and date associated with each determination made with the SPC method. To locate the diamonds, neither noise nor averaging have been included. b) Lomb periodogram of the Dec-simulated coordinates (diamonds) shown in Fig. 5a. The RA coordinate periodogram is exactly the same as that of the Dec data, only small differences in spectral power (height of the peak) can be appreciated.

4.1. Simulating the homogeneous body

The orientation of the spin axis in J2000 was calculated for a body with inertia moments equal to those derived from the geometric shape of 67P assuminga homogeneous density distribution. Calculations were performed for different excitation levels E_l by visually checking if the RA/Dec distribution of theoretical data agreed with the observed distribution (displayed in Fig. 2). This case was also considered to study the effect of averaging on the spin axis RA/Dec distribution and on the detected periodicities. The simulation for case 1 with avalue of E_l = 1.9 × 10^-6 is shown in Fig. 5a. This figure shows that the instantaneous simulated spin axis orientation displays a pattern that may barely be considered similar to the observationally derived one. An important difference between observational and theoretical coordinates comes from the inner hole, which is clearly visible in the simulated data and missing in the actual spin axis orientation. This inner hole exists because the angle of nutation has a lowest value different from zero, mainly due to the excitation level. To decrease the size of this inner hole, the excitation level must be reduced. This would demand a modification of the inertia moments, which would imply some inhomogeneity, if we still wish to cover the observational range in RA and Dec coordinates.In addition to visual comparison, if a two-sided Kolmogorov-Smirnov test is applied to estimate compatibility between the observationally derived and the theoretical RA distributions, the probability of being compatible is lower than 0.05. This probability drops to 0.02 if Dec distributions are compared. Differences and poor compatibility between observational and theoretical coordinates are at the moment expected because simulated data are instantaneous values while the observationally derived coordinates are the result of an averaging within 10 h. Before we describe the effect of averaging on the RA/Dec distributions, we describe the detected periodicities when a time-series analysis software is applied to the RA/Dec distribution. This is to be compared with Fig. 3.

4.1.1. Meaning of the detected periodicities

Figure 5b shows the Lomb periodogram of the simulated Dec distribution shown in Fig. 5a that is considered ideal, meaning that neither noise nor averaging were included. The periodogram of the simulated RA coordinates would be exactly the same as that of simulated Dec coordinates, with only small differences in spectral power. This figure shows that two clear periodicities are found. They are at 9.67 h and at 17.3 h. From the Euler equations it is possible to estimate that, under case 1 circumstances and for the inertia moments and excitation level considered, this simulation has a P = 2π/ Ω = 12.4 h while the theoretical P_φ and P_ψ would be 9.66 h and 43.75 h, respectively.This means that while the rotational period would be 12.4 h, the body would be spinning with a mean period of 43.75 h around the third axis of inertia and that this axis, at the same time, would be precessing around the angular momentum vector with a mean precession period of 9.66 h. The periodogram shows that the most significant detected periodicity (working with ideal data) corresponds exactly to P_φ, the mean precession period (for the inertia moments and E_l considered, i.e., 9.66 h). The other significant periodicity corresponds to a combination P_c of P_φ and P_ψ, particularly to (9)which in our simulation is exactly 17.3 h. P_c is therefore obtained as a combination of the precession period and half the spin period, which is exactly the nutation period. These two significant periodicities, P_φ and P_c, are the only ones found in all simulations performed considering a wide range of inertia moments and excitation levels. None of the simulations showed a periodicity different from P_φ and/or P_c in the analysis of the simulated RA/Dec coordinates. In none of the simulations, P_ψ or 2π/ Ω were detected. Moreover, from the analysis of instantaneous simulated data, the most significant periodicity was always P_φ followed by P_c.

When comparing the periodogram of the observationally derived coordinates (Fig. 3) with the obtained periodogram from simulated data (Fig. 5b), we see that they are very different from each other. This indicates, again, that the homogeneous case seems to be incompatible withour measurements, at least when averaging is not considered and under case 1 circumstances.

Fig. 6 a) As in Fig. 5, bullets (•) show the spin axis orientation as determined withthe SPC method (Fig. 2). Small dots (·) are the spin vector orientation for a simulated body with inertia moments equal to those derived from the shape of 67P assuming homogeneity, and an El = 1.9 × 10-6 under case 1 assumptions , that is, exactly the same circumstances as those of the simulation shown in Fig. 5a. These data are displayed with a time separation of 1 h. Diamonds (◆), to be compared with the bullet distribution and also with Fig. 5a, correspond to the simulated spin axis orientation at the moments associated with each determination with the SPC method, but now calculated as the average value of the instantaneous coordinates within intervals of 10 h around these moments. b) Lomb periodogram of the simulated averaged Dec coordinates (diamonds) shown in Fig. 6a.

As for case 2 circumstances (P_φ = 12.40 h), we recall that it is not possible to estimate the rotational energy without attributing a meaning to the 270 h periodicity. Nevertheless, it is possible to proceed the other way around. Considering the inertia moments of the homogeneous body, we may search for the rotational energy corresponding to avalue of P_φ equal to 12.40 h. If this exercise were carried out, we would find that the body spins with a period P_ψ ≈ 56.16 h and precesses with a P_φ = 12.4 h, displaying a total rotational period of P ≈ 15.9 h. This situation would lead to a P_c of approximately 22.21 h. None of these periodicities can be related with the RA/Dec periodicity of 270h. Additionally, as previously, the periodogram of simulated data only shows the periodicities corresponding to P_φ and P_c, that is, 12.40 h and 22.21 h.

4.1.2. Averaging effect

The observational RA/Dec coordinates derivedby the SPC method include an averaging within 10 h, which was not included in previous figures (Figs. 5a and 5b). Figures 6a and 6b show the effect of such averaging both in the RA/Dec pattern and in the Lomb periodogram. The first evident consequence is that the averaging significantly reduces the range of the RA and Dec coordinates. This is a logical result because the averaging interval (10 h) is larger than the precession period (9.66 h). Consequently, to make the simulation compatible with observations, an excitation level higher than the one we considered is necessary. As for the Lomb periodogram (Fig. 6b), while the same two periodicities are still clearly detected (9.67 h and 17.3 h), Figs. 5b and 6b are different. The averaging reduces the significance of the periodicity-affected (P_φ), while it increases that of P_c. It could be said that the averaging “transfers” significance from P_φ to P_c. Figure 6a clearly shows that the averaging makes the simulated RA/Dec pattern incompatible with the observational pattern. Keeping the inertia moments, the only variable available to make the two data sets more compatible would be E_l. It could be shown that an E_l of about 3 × 10^-5 is necessary when averaging is included in theoretical data to have a θ_max covering the full range of the observationally derived RA/Dec data. The problem is that when E_l is varied toincrease the value of θ_max, θ_min is also significantly increased. For the case of E_l = 3 × 10^-5, we verified that θ_min, when an averaging of 10 h is included in the theoretical data, is on the order of 0.2^◦ for the inertia moments of the homogeneous body, that is, a quantity clearly incompatible with the observationally derived RA/Dec data. Again, to make theoretical and observationally derived results compatible, it would be necessary to modify the inertia moments with regard to those of the homogenous body in addition to increasing E_l.

4.2. Slightly inhomogeneous case: the near prolate body

After we discarded that P_φ, the precession period, can be identified with the RA/Dec periodicity, the most simple possible interpretation of this observational periodicity compatible with the Euler equations is to consider that it corresponds to the period P_ψ. Equation (1) indicates that P_ψ depends on the inertia moments, the rotational energy, and the excitation level. It can be shown that the dependence with inertia moments is stronger than that with excitation level, for example, especially if the latter is close to 1. Therefore, by fixing the excitation level to reasonable values, it is possible to find the approximate relationship I_z/I_y that corresponds to avalue P_ψ of approximately 270 h, always assuming that L ≈ I_zΩ. Under case 1 circumstances, we recall that Ω = 0.0014075 rad s^-1. Under case 2 assumptions, if rad s^-1, and rad s^-1, it is possible to approach rad s^-1 if the excitation level is considered to be low.

Some examples of these ratios are presented in Table 1. It shows that even if it is possible to find an I_z/I_y ratio for any value of I_y, the body could not be homogeneous under the aforementioned interpretation. As an example, for the excitation level considered in Table 1 and an I_y of approximately 1.83 (the value estimated from homogeneity), I_z would have to be 1.835, lower than thevalue of 1.98 estimated when the body is considered to be homogeneous. We verified that for higher excitation levels the ratio I_z/I_y necessary to obtain a value of P_ψ equal to 270 h also increases, but hardly reaches 1% when I_y is that of the homogeneous body. Therefore, if P_ψ has to be around 270 h, regardless of the interpretation of thedetected 12.40 h periodicity (Mottola et al. 2014), the body must be slightly inhomogeneous and near-prolate (i.e., I_y ≈ I_z) from the standpoint of inertia moments.

Fig. 7 a) As previously, bullets (•) show the observationally derived spin axis (Fig. 2). Small dots (·) are the spin axis orientation for a simulated body with inertia moments 1:1.90:1.905 and an El = 2.9 × 10-6 (see text). The data are displayed with a time separation of 1 h. As previously, diamonds (◆), to be compared with the bullets distribution, correspond to the averaging within 10 h of the instantaneous theoretical values. b) Lomb periodogram of the Dec-simulated data (diamonds) shown in Fig. 7a. The x-axis here ranges up to 300 h to show that no significant periodicity is detected beyond Pc = 13 h. Given the wide range covered by the x-axis, the two significant periodicities Pφ = 11.85 h and Pc appear very close in the figure.

We also considered both the theoretical RA/Dec distribution and the Lomb periodograms of a body with the characteristics described above to compare them with the observationally derived ones. For illustrative purposes, a simulation under case 1 circumstances in which I_y = 1.90, that is, it is increased by 4% with regard to the homogeneous situation, and I_z = 1.0025·I_y = 1.905, that is, a reduction of 4% compared to the homogeneous case, is shown in Figs. 7a and 7b. The selected I_z and I_y, which can be considered reasonable in the sense that they are not very different from the actual ones assuming homogeneity, give the ratio necessary (see Table 1) for P_ψ at 268 h for an excitation level of 2.9 × 10^-6 (quantity selected to cover a range of RA/Dec coordinates similar to the observational one).

Figure 7a shows the theoretical RA/Dec distribution for this simulation considering an averaging of data within 10 h, compared to the observationally derived distribution. Both distributions, simulated and real, arevisually similar, although the simulated one shows a more regular pattern.A two-sided K-S test would show that theoretical and observational RA distributions are similar with a probability of approximately 0.5. If theoretical and observational Dec distributions are compared, the probability is much lower, of approximately 0.02. Regardless of the similarity of the RA/Dec distributions, Fig. 7b shows that the Lomb periodogram of the theoretical data is quite different from the observational one (Fig. 3). For the inertia moments and excitation level considered in this simulation, we calculated that if we assume that Ω = 0.00014075 rad s^-1, then P_φ = 11.85 h, P_ψ = 268 h, and P_c = 13.0 h. The periodogram shows that as previously (Figs. 5b and 6b), only two significant periodicities are detected, coinciding, again as previously, with P_φ = 11.85 h, and P_c = 13.0 h. The periodogram does not show any trace of periodicity at P_ψ = 270 h. This test was repeated for different inertia moments and excitation levels (covering the range in Table 1) and also for case 2 circumstances, obtaining always the same results: no simulation showed a significant periodicity at 270 h although P_ψ ≈ 270 h in all of them. The periodograms only showed the periodicities corresponding to P_φ and P_c.

Thus, in spite of the reasonable similarity between the simulated and observationally derived RA/Dec distributions (Fig. 7a), we find that the near prolate case is unlikely as a possible interpretation of the complex rotation of 67P. On the one hand, it is highly unlikely that such an irregular nucleus, with such a complex surface, might be near prolate from the standpoint of inertia moments. On the other hand, Lomb periodograms of theoretical data, showing just P_φ and P_c as detectable periodicities, do not support the interpretation that P_ψ is at approximately 270 h.

5. Interpretating the RA/Dec periodicity

Fig. 8 Relationship between Iz and Iy for a Pc of approximately 270 h for case 1 (continuous line) and case 2 (dashed line). Both Iy and Iz are normalized to Ix. The asterisk (∗) represents the homogeneous body.

A very likely interpretation of the complex rotation of 67P, supported by the Lomb periodograms of theoretical data, is that what has been detected in the periodograms of the observationally derived data (Fig. 3), that is, that the RA/Dec periodicity is actually P_c (Eq. (9)), with P_φ missing probably due to both averaging and errors (noise) in the determination of the RA/Dec coordinates.

It was previously mentioned that after the inertia moments are defined, the excitation level produces minor changes in the periodicities detected in the periodograms. Thus, it is possible to easily determine the I_z/I_y ratio, which provides a periodicity P_c at approximately 270 h for cases 1 and2. These ratios are shown in Fig. 8, always assuming that I_x is equal to 1, meaning that the inertia moments are always normalized to the smallest inertia moment. As a result of the complex dependence of P_c, this figure was built from a set of simulations by considering discrete values of I_y and performing a systematic search of the I_z defining a P_c with the desired value. A linear fit of the discrete values allows us to define the continuos relationship between I_z and I_y,(10)Strictly speaking, these relationships would only be valid for the E_l used to build it. Nevertheless, the effect on P_c of varying E_l within reasonable values (searching for compatibility with the range of the observationally derived RA/Dec coordinates) is always smaller than the error in the periodicity detected in the periodograms of the observationally derived RA and Dec coordinates (12 h, from Jorda et al. 2016). Considering these two relationships, forcases 1 and2 the differences are irrelevant as for the inertia moments. This accuracy in inertia moments would be very difficult to reach from any observational data given the uncertainty in cometary nature. Nevertheless, in spite of the similarity in inertia moments between cases 1 and2, the rotational characteristics of the body depend very much on the considered case. Recalling that we always considered that L ≈ I_zΩ, under case 1 assumptions, a body with an I_z given by the first row of Eq. (10) would have P ≈ 12.40 h, P_φ ≈ 6.35 h, P_ψ ≈ 13 h, and P_c ≈ 270 h, with small variations of these values depending on the excitation level (assumed to be small). This means that the body would be rotating with a period of approximately 12.40 h, with a complex rotational motion, spinning around its principal axis with the highest inertia moment (z) with a period of approximately 13 h, while this axis precesses around the angular momentum vector with a period of approximately 6.35 h. Under case 2 circumstances, a body with I_z given by the second row of Eq. (10) would have P ≈ 23.71 h, P_φ ≈ 12.40 h, P_ψ ≈ 26 h, and P_c ≈ 270 h. A body like this would be spinning around the third axis of inertia with a period of 26 h while precessing around the angular momentum with a period of 12.40 h, the rotational period being approximately 23.7 h. This fact practically rules out case 2. Under these circumstances, the body as seen from space would reach a similar orientation approximately every 24 h, which is incompatible with the images obtained with OSIRIS. Therefore, case 2 was not considered below, except for some last comments.

Under the present interpretation (that what was detected in the periodograms of observationally derived data is P_c), the previous relations would indicate that 67P cannot be homogeneous. Equation (10), for case 1, tell us that to have a period P_c around 270 h, I_z − I_y ≈ 0.96, a quantity much larger than the corresponding difference assuming that the body has a homogeneous density, I_z − I_y = 1.98−1.83 = 0.15.

In other words, while the homogeneous situation implies a ratio I_z/I_y = 1.08, if the RA/Dec periodicity is P_c, the ratio I_z/I_y would be higher, ranging, for example, between 1.4 for I_y = 2.5 and 1.6 for I_y = 1.5 for case 1 circumstances. Except for inhomogeneity, always assuming that our starting hypothesis are valid, nothing else can be said in principle about the internal mass distribution. Inertia moments I_y and I_z are estimated from the integrals of the mass with the distance across the xz and yx planes, respectively. This means that the two principal inertia moments share the effect of density along the x-axis, nearly parallel to the largest geometric axis of 67P. If the ratio I_z/I_y is higher than that of the homogeneous case, the density along the z-axis might be lower and/or that the density along the y-axis might be higher than in the homogeneous case. Interestingly, z-axis goes across the neck of 67P, a region certainly peculiar even in the already quite odd surface of 67P. In any case, we understand that the previous I_z − I_y difference can be considered quite extreme. The equivalent ellipsoid would tend to be oblate, with an intermediate-to-short axes ratio of approximately 7. This difference in inertia moments would therefore request a comparatively high increase of mass especially along the y-axis, which is practically impossible to reach with smooth and continuous variations of density across the body. For example, and in principle, the inertia moments ratio found in this study would not be compatible with a comparatively small difference in densities of the two lobes (Rickman et al. 2015; Jorda et al. 2016). One physical possibility to explain the estimated inertia moments could be a body with mass concentrated at the equatorial plane and at the surface, with a significant density variation. Another physical possibility to explain the inertia moments, again speculating with the internal structure of 67P, might be a body with comparatively small and randomly distributed holes in its interior, with mean sizes on the order of 100 m. This latter picture would point to the heterogeneity suggested by Vincent et al. (2015). These authors concluded from their analysis of the size and spatial distribution of pits that larger heterogeneities exist in the properties of the first few hundred meters below the current nucleus surface of 67P. The possibilities given above are just some examples among many others that need appropriate modeling to be confirmed or rejected. Without constraining the ratio of inertia moments, any attempt to find a possible explanation would require additional constraints from independent measurements and/or additional hypotheses on the internal structure of 67P in addition to the appropriate modeling.

Fig. 9 a) K-S probability as a function of Iy and El when the theoretical spin axis RA distribution is compared to the observationally derived one for case 1. b) The same as a), but in this case the theoretical and observationally derived Dec distributions are compared. In the two figures, the continuous white line, El,RA = El [ Iy ], indicates the excitation level for each Iy associated with the highest K-S probability when RA distributions are compared. The dashed white line, El,Dec = El [ Iy ], represents the excitation level for each Iy associated with the highest K-S probability when the Dec distributions are compared.

In addition to Vincent et al. (2015), at least certain heterogeneity in 67P has previously been suggested so far by Jorda et al. (2016), who detected an offset between the center of mass of the body shape assuming homogeneity and the origin of the Cheops frame. This offset suggests some heterogeneity. The modeling to interpret the inertia moments must also account for this offset. Luspay-Kuti et al. (2015) also pointed out that the compositional heterogeneity detected in the coma of 67P might suggest chemical heterogeneity in the southern hemisphere of the nucleus. Analysis of the measurements performed by the CONSERT instrument may shed some light on the internal mass distribution. Although it is still premature to draw definite conclusions because of the complexity of analyzing CONSERT data, last results points out some inhomogeneity, at least for the small lobe of 67P. Ciarletti et al. (2015) explained CONSERT measurements as due to a decrease with depth of the dielectric constant. This was interpreted by the authors as a possible increase of porosity with depth and/or a possible decrease of the dust-to-ice ratio. Both interpretations could point to an increase in the inertia moments in comparison with those of a homogeneous body and to a possible increase of the inertia moment ratio. During the revision process, Pätzold et al. (2016) presented their results on the gravitational field of 67P as determined from the Radio Science Instrument (RSI). According to the authors, these measurement seem to point to a homogeneous nucleus, ruling out large caverns in the nucleus interior. Nevertheless, the spatial scale of the ruled-out caverns is not clear yet, and we need to wait for the future controlled impact of the Rosetta orbiter with the surface nucleus in September 2016. At this moment, the complex gravity field could reveal details on the internal structure on the scale of hundred of meters.

The internal structure and its possible homogeneity are key characteristics for understanding cometary nuclei formation. Inhomogeneity has been suggested for several cometary nuclei as derived from different sources. For example, Kawakita et al. (1997) concluded that the presumable parent of comets C/1996 Q1 (Tabur) and C/1988A1 (Liller) was inhomogeneous because the gas-to-dust ratio measured in the two fragments was very different. Gibb et al. (2007) also suggested cometary heterogeneity from a compositional analysis of fragments of comet C/2001 A2 (LINEAR). There is other possible evidence of heterogeneity from other fragmenting comets, although, as pointed out by Dello Russo et al. (2007), no confident conclusion can be drawn. Chemical inhomogeneities related to possible formation mechanisms were also suggested for comets 8P/Tuttle (Bonev et al. 2008), 9P/Tempel 1 (Feaga et al. 2007), and 103P/Hartley 2 (A’Hearn et al. 2011), although evolutionary circumstances (mainly seasonal effects) cannot be ruled out (see, e.g., A’Hearn 2011).

As expected given the cometary diversity, there are also arguments favoring homogeneity in cometary nuclei. From high-quality spectra, Dello Russo et al. (2007) found that relative abundances measured from fragments B and C of comet 73P/Schwassmann-Wachmann 3 were remarkably similar, which suggests chemical homogeneity of the nucleus of 73P. We would like to stress again that discussing the heterogeneity in 67P is still very premature. More certain conclusions need from dedicated calculations and a more complete analysis of the internal structure of 67P by CONSERT and RSI data. If finally the suggested I_z − I_y is much higher than the actual one, it would be necessary to revisit our starting hypothesis, mainly rigid-body and free-torque conditions. At the end of the paper, we briefly comment on the first of these assumptions.

5.1. Is it possible to constrain the inertia moment ratio?

After we defined the relation between I_y and I_z to have P_c on the order of 270 h, we performed a full set of simulations for case 1, considering as free parameters both I_y and E_l, and fixing I_z according to Eq. (10). In the simulations, I_y ranged between 1.2 and 2.6, which is a reasonable interval, while a range between 3 × 10^-4 and 2 × 10^-6 was explored for E_l. This latter interval was defined after a trial-and-error procedure in which we searched for visual compatibility between observations and simulations. The results from these simulations were compared with the observationally derived data, in particular the RA/Dec pattern and the Lomb periodograms. Considering the periodograms, we verified that all simulations performed under case 1 circumstances have the most significant peak located at P_c ≈ 270 h followed by a second peak located at P_φ ≈ 6.35 h. As for the RA/Dec patterns, a two-sided Kolmogorov-Smirnov test was performed to statistically evaluate the significance of the similarity between the observationally derived and the theoretical distributions. The two-sided K-S test was performed separately over RA and Dec distributions.

Fig. 10 Simulation with the highest K-S probability from the comparison of RA coordinate distributions. a) As previously, bullets (•) show the observationally derived spin axis (Fig. 2). Diamonds (◆) correspond to the simulated spin axis orientation of a body with 1:1.90:2.84, and El = 6.7 × 10-5 when an averaging of 10 h and a uniformly distributed random noise with an amplitude of 0.03◦ have been considered. b) (Thick lines) Lomb periodograms of the RA and Dec simulated data (diamonds shown in Fig. 10a). The Lomb periodograms of the observationally derived coordinates (Fig. 3) have also been included (with thin lines).

Fig. 11 Simulation with the highest K-S probability from the comparison of Dec coordinate distributions. a) As previously, bullets (•) show the observationally derived spin axis (Fig. 2). Diamonds (◆) correspond to the simulated spin axis orientation of a body with 1:2.25:3.19, and El = 3.1 × 10-5 when an averaging of 10 h and a uniformly distributed random noise with an amplitude of 0.03◦ have been considered. b) (Thick lines) Lomb periodograms of the RA and Dec simulated data (diamonds shown in Fig. 11a). The Lomb periodograms of the observationally derived coordinates (Fig. 3) have also been included (with thin lines).

Fig. 12 K-S probability estimated as the mean of the RA and Dec probabilities for the two excitation curves, El = El [ Iy ], leading to the highest probability of the K-S tests performed on RA (•) and Dec (⋄) data separately.

The similarity between the RA theoretical distribution and the observationally derived RA pattern, understood as the probability that both patterns belong to the same distribution as determined from the two-sided K-S statistic as a function of both I_y and E_l, is shown in Fig. 9a for case 1. From the K-S statistics of the RA patterns, the region with I_y between 1.8 and 2.0 can be clearly identified with the highest probability. The best fit, identified as the highest probability, 0.94, is obtained for I_y = 1.90, I_z = 2.84 and E_l = 6.7 × 10^-5. This body would be rotating with P ≈ 12.4 h, P_φ = 6.34 h, P_ψ = 12.99 h, and P_c = 273 h. The simulated spin axis orientation and corresponding Lomb periodograms for these parameters are shown in Figs. 10a and 10b. Considering Fig. 10a, it may be accepted that theoretical spin axis coordinates have a distribution similar to the observationally derived one. For this distribution, in which a uniformly distributed noise with an amplitude of 0.03^◦ was included, the Lomb periodograms of the simulated RA and Dec coordinates show that the most significant peak is located at a periodicity of 275 h, identifiable with P_c, but with a spectral power significantly higher than the observational one. In addition, the periodicity corresponding to P_φ = 6.35 h can clearly be identified in simulated data, even though noise and an averaging significantly larger than P_φ were added. In principle, the differences between simulation and observations might arise because actual errors in observations are larger than formal errors. If this were true, adding a larger noise to the simulated data could erase the P_φ signature and also reduce the spectral power of P_c. Nevertheless, in this case, the K-S tests and the highest probability case would not be valid anymore because it would be necessary to at least modify E_l to adjust the range of the observationally derived RA/Dec.

The K-S probability that compares the observationally derived Dec distribution with the theoretical one as a function of I_y and E_l is shown in Fig. 9b for case 1. This figure shows that several regions with comparatively high probabilities can be identified, with a highest value of 0.82. Even though it is smaller than the value obtained from the RA comparison, this value could be considered to have a similar significance because it was estimated that the dispersion in the probability due to the observational error is on the order of 0.1. The K-S probabilities within that range may therefore be considered to have similar significance. From the Dec comparison, the highest probability is obtained for the combination I_y = 2.25, I_z = 3.19, and E_l = 3.1 × 10^-5. As previously, and as in all simulations under case 1 circumstances, this body would also be rotating with P ≈ 12.4 h, P_φ = 6.35 h, P_ψ = 13.0 h, and P_c = 271 h. The simulated spin axis orientation and corresponding Lomb periodograms for these parameters are shown in Figs. 11a and 11b. The comparison between simulated and observationally derived spin axis coordinates show that the discrepancies are noticeable here, especially in the RA range. Nevertheless, the Lomb periodogram still shows the most significant periodicity close to P_c, at 268 h, with very weak traces of P_φ. We note from simulations not shown that in addition to noise and averaging, the value of I_y also affects the significance (i.e., height of the peak in Lomb periodograms) of P_φ. For values higher than I_y ≈ 2.0, the peak corresponding to P_φ is hardly detectable. As for P_c, the spectral power of simulated data is much higher, as in the RA case, than that of the observationally derived ones. Again, a larger actual error in observations might conciliate both spin axis RA/Dec distribution and the spectral power of P_c, but it would be necessary to recalculate the K-S probability.

Bearing in mind that the dispersion introduced by the error in RA/Dec in the K-S probability is on the order of 0.1, we cannot favor any set of parameters. An additional statistical test was performed to determine whether it was possible to find a best common region for RA and Dec coordinates together. For each I_y, the inertia moment I_z was calculated by using Eq. (10), and first E_l was defined by using the continouos white line shown in Figs. 9a and b.This white line defines the highest probability as a function of I_y when RA coordinates are compared. After we defined the set of [I_y, I_z and E_l], the theoretically obtained RA and Dec coordinates were perturbed with random and uniform noise with an amplitude of 0.03^◦. One hundred perturbed RA and Dec patterns were built for each [I_y, I_z and E_l] combination. The K-S probabilities for all these patterns were calculated by comparing them with the observationally derived RA and Dec distributions, building the final K-S probability as ([K-S]_RA+[K-S]_Dec)/2. We then repeated the calculations by using the dashed white line shown in Figs. 9a and 9b, that is, the excitation level corresponding to the highest K-S probability for each I_y when simulated and observationally derived Dec coordinates were compared. The results, in terms of median and standard deviation of the 100 different random patterns, are shown in Fig. 12. This figure shows that when the RA and Dec K-S probability is averaged and random low-amplitude noise is added, there is no preferred combination of [I_y, I_z] for any of the two excitation curves that best fit the RA and Dec coordinates separately. Figure 12 shows that any I_y within the range 1.5–2.3 is virtually similar in terms of K-S probability when RA and Dec coordinates are considered together. Therefore, no [I_y, I_z and E_l] combination can be defined from the K-S tests.

5.2. Is there any trace of P_φ in the periodogram of observationally derived data?

Fig. 13 Lomb periodograms of the observationally derived RA (black) and Dec (Red) coordinates (Fig. 3) for periodicities lower than 14 h. Arrows indicate the position of Pφ for the two possibilities discussed in the text, 6.35 h and 12.40 h.

We showed that simulations indicate that P_φ, the precessional period, may be detectable in the periodograms even when an averaging period longer than P_φ is used when retrieving RA/Dec coordinates. The reason might be that the φ angle is not strictly periodic and that P_φ is defined just from the mean value of the instantaneous . Therefore, even when the averaging is large, some signatures of may remain in the space position of the spin axis. We therefore explored the periodogram of RA/Dec coordinates in the region of periods shorter than 14 h to determine whether there was any trace of the 6.35 h periodicity. The relevant region of the periodogram of the observationally derived RA/Dec coordinates is shown in Fig. 13. This figure shows that neither RA nor Dec periodograms show any signature at 6.35 h. Nevertheless, the RA coordinates periodogram displays a weak signature at 12.40 h that is also present in the Dec coordinates periodogram. If our interpretation that the highest peak in the periodogram is P_c with some potential signatures of P_φ is correct, this weak signature would point to case 2 circumstances. Nevertheless, as stated above, this would imply that 67P spins around the axis with the highest inertia moment with a period of 26 h while this axis precesses around the angular momentum vector with a period of 12.40 h, with P_c ≈ 270 h and the total rotational period is approximately 23.7 h. This solution would require reinterpreting the periodicity detected by Mottola et al. (2014) as the precession period, but more importantly, it would not be compatible with OSIRIS images. In principle, the signature at 12.40 h is really weak with a significance level of 3%, and therefore could be spurious. Nevertheless, its coincidence with the periodicity detected by Mottola et al. (2014) indicates that it might be real. In that case, its presence is unexpected, and based on the information obtained from theoretical periodograms, there is currently no satisfactory explanation, except if it is a coincidence.

5.3. Comments on the rigid-body assumption

The previous analysis and discussions were performed under the assumption that the body is considered rigid. In principle, this assumption is not necessarily true and it might be that non-rigidity is the cause of the motion of the spin axis described in Fig. 2, or at least of part of it. It is beyond the scope of this paper to study the sources of non-rigidity in depth as a possible cause of the RA/Dec periodicity, but we consider the potential effects of the most likely ones below.

A body rotating in complex mode, that is, on an excited level, dissipates rotational energy mostly as heat as a result of the friction of small and slow internal mass displacements. Because this is an internal adjustment, the angular momentum does not change, which means that the small displacements are translated into changes of the rotational velocity vector, inertia moments, and mass center. We know that the rotational period (2π/ Ω) or the precessional period (P_φ), depending on the interpretation, did not change within a very small error during the time covered by the observations used in this study. This would imply that inertia moments did not change and that therefore any motion of the spin axis that is due to some potential damping of the excitation energy is safely discarded. This agrees with the estimate of the timescale necessary to damp an excited rotational state by internal dissipation as derived by Burns & Safronov (1973) and applied to some comets by Jorda & Gutiérrez (2000). These latter authors showed that by using typical values for the parameters involved in the timescale, it may range from a few thousand of years (for large nuclei) to up to 10⁷ years (for small nuclei), which is considerably longer than the period of time covered by our observations. Other sources of non-rigidity such as erosion by activity and deposition of the heaviest particles falling back onto the surface that also induce changes in the inertia moments (and therefore in the rotational velocity) have equally very long timescales (see, e.g., Jorda & Gutiérrez 2000).

Related to non-rigidity, one possibility that immediately comes to mind is whether the comet nucleus may be suffering a kind of Chandler wobble. Before considering this possibility, we describe this motion including the Euler angles. The Chandler wobble is generally described as the torqueless precession of the Earth with a period of about 433 days. This means that the motion of the axis of rotation (Ω, North Pole at Earth) around the axis of symmetry of Earth (z) which is assumed to be an oblate body. It is important not to confuse this precession with the term precession as used throughout this study, meaning the motion of the z -axis around the angular momentum. This may lead to some confusion with the periodicities. In our case, as we described, the mean precession period would always be P_φ. If the motion of the rotation vector around the third axis of inertia is considered, the former precesses around the axis with the highest inertia moment (z) exactly with the period P_ψ, that is, our spin period. On the one hand, Earth, which is assumed to be an oblate body, does not fulfill the Euler equations and the period of Ω around the axis of symmetry is longer than expected from Euler equations for the known inertia moments. On the other hand, it has been calculated that the Chandler wobble would be damped in less than 100 years, unless some internal force were causing it. A clear explanation for the Chandler wobble is still lacking, but it might be related to fluctuating pressure at the bottom of the oceans, caused by temperature, salinity, and circulation changes (Gross 2000), or even to centrifugal deformation of the portion of the Earth’s mass contained in circulating fluids (Jenkins 2015). Regardless of the final explanation, it seems that Chandler wobble is due to the presence of liquid (or, in general, movable) mass that may deform and/or act on the rest of the body. Interestingly, the deformability of Earth may be estimated as (see, e.g., Jenkins 2015) (11)where ω_Chandler is the frequency associated with the Chandler wobble and ω_Euler is the frequency of the precession of the rotation axis around the axis with the highest inertia moment. According to the Euler equations for the inertia moments of Earth, ω_Euler is estimated to be 2π/P_ψ = 2π/ 306 days^-1. If the corresponding figures are used, k of Earth is estimated to be 0.293, which is fully consistent with the Love number k₂ computed form the magnitudes of terrestrial and oceanic tides (see, e.g., Stacey & Davis 2008).

Assuming that in our case Ω were precessing around the axis with the highest inertia moment and that for some reason it translates its periodicity into the detected RA/Dec periodicity, we could also estimate the deformability necessary by using the previous expression. Assuming homogeneity, that is, that the inertia moments as derived from the body shape if density is assumed to be homogeneous, the precession period of Ω around the axis with the highest inertia moment would be P_{ψ Euler} = 43.75 h for case 1. If our body were suffering some type of Chandler wobble with a P_Chandler = 270 h (i.e., the RA/Dec periodicity), then the deformability could be estimated as 0.838. This would represent an unreasonably high number if deformability, as for Earth, is identified with the Love number k₂.

The previous discussion does not totally discard that body deformation, or non-rigidity in general,is the cause ofthe RA/Dec periodicity detected from observations. An in-depth study is beyond the scope of this paper, which is based on the rigid-body assumption. We merely wished to illustrate that the most intuitive possibilities related to non-rigidity are ruled out.

6. Summary and conclusions

Spin axis orientation data as obtained during the SPC procedure to retrieve the shape of 67P show evidence that its nucleus was rotating in complex mode. This evidence is compatible with and confirms the early detection of possible complex rotation by Preusker et al. (2015). When the spin axis coordinates were studied by means of different time-series analysis techniques, the periodograms showed a very significant periodicity at approximately 270 h, different from the periodicity determined by Mottola et al. (2014). In this study, simulations of slightly excited rigid-bodies fulfilling Euler equations under free-torque conditions were performed to interpret the periodicity detected in the coordinates of the spin axis. The main results of these simulations are summarized in Table 2.

From simulations considering slightly excited bodies, we found that Lomb periodograms of the RA/Dec distributions only show two significant periodicities. In ideal data, the most significant periodicity corresponds to the mean precession period, P_φ. The second periodicity appearing in the periodograms always corresponds to P_c (Eq. (9)), a combination of the precession periodicity, P_φ, and the nutation periodicity, the latter being half the intrinsic rotation period, P_ψ.

When ideal simulated data were perturbed by averaging over a period on the order of P_φ and noise was added, the significance was transferred from P_φ to P_c, the latter becoming the most significant peak in the periodograms. The peak corresponding to P_φ may eventually disappear, leaving P_c as the single clearly detected periodicity.

When we assumed that the periodicity detected by Mottola is either the rotational period, P = 2π/ Ω, or the precessional period, P_φ, then P_ψ, the spin around the principal axis corresponding to the highest inertia moment, might in principle be associated with the periodicity detected in observationally derived RA/Dec coordinates (≈270 h). This possibility was ruled out for two reasons. On the one hand, it would require that 67P must be nearly prolate from the standpoint of the inertia moments, which is unlikely because of the highly irregular shape of its nucleus. On the other hand, no single simulation showed the periodicity corresponding to P_ψ in the periodograms. Lomb periodograms of simulated data in which the body is rotating with a P_ψ around 270 h are quite different from the periodogram of the observationally derived data.

According to the results obtained from the simulations, the significant peak found in the periodograms of the observationally derived data could be interpreted as the combination P_c of periodicities P_φ and P_ψ. If the 12.40 h periodicity detected by Mottola et al. (2014) is considered the rotational period, I_y and I_z must follow the relation I_z = 0.96 + 0.99I_y. For this relation between inertia moments, the third inertia axis of the body (spin axis) would be precessing around the angular momentum vector with a mean period P_φ ≈ 6.35 h and, at the same time, the body would be spinning around its third axis of inertia with a period P_ψ = 13.0 h, being the rotational period, P = 2π/ Ω, around 12.4 h.

Under these circumstances, a significant periodicity in the periodograms P_c ≈ 270 h of the spin axis RA/Dec simulated coordinates was found, which became highly significant when noise and an averaging of 10 h (as in observations) were considered.

The ratio I_z/I_y when the RA/Dec periodicity is identified with P_c, regardless of the meaning of the periodicity detected by Mottola et al. (2014), is significantly higher than the ratio I_z/I_y estimated assuming that 67P is homogeneous. If the RA/Dec periodicity around 270 h appearing in the periodograms is identified with P_c, 67P cannot have a homogenous density distribution.

To evaluate whether it is possible to constrain the inertia moments and excitation level, a systematic search of the probability of compatibility between simulated and actual RA/Dec patterns by means of two-sided K-S tests was performed. Even if it is possible to find very significant combinations of parameters [I_y,I_z,E_l] when RA and Dec coordinates are considered separately, K-S probabilities when RA and Dec data are considered together do not allow constraining the inertia moments and excitation level.

Acknowledgments

OSIRIS was built by a consortium of the Max-Planck- Institut für Sonnensystemforschung, Göttingen, Germany, CISAS University of Padova, Italy, the Laboratoire d’Astrophysique de Marseille, France, the Instituto de Astrofísica de Andalucía, CSIC, Granada, Spain, the Research and Scientific Support Department of the European Space Agency, Noordwijk,The Netherlands, the Instituto Nacional de Técnica Aeroespacial, Madrid, Spain, the Universidad Politécnica de Madrid, Spain, the Department of Physics and Astronomy of Uppsala University, Sweden, and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universität Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France(CNES), Italy(ASI), Spain(MEC), Sweden(SNSB), and the ESA Technical Directorate is gratefully acknowledged. P.J.G. and L.M.L. particularly acknowledge funding from projects AyA 2012-32237 and ESP 2014-54062-R.

References

All Tables

All Figures

	Fig. 1 SPC shape of 67P/Churyumov-Gerasimenko (SHAP5, version 1.2, Jorda et al. 2016) in the so-called Cheops frame (see Scholten et al. 2015). The spin axis is the z-axis.
In the text

	Fig. 2 Coordinates in J2000 of the spin axis of 67P for the time interval between –376 and –250 days from perihelion. These coordinates are obtained during the SPC procedure employed to retrieve the shape of 67P (Jorda et al. 2016).
In the text

	Fig. 3 Lomb periodograms of the RA and Dec coordinates of the spin axis displayed in Fig. 2. A significant peak (with a significance higher than 99.99%) is located at a periodicity of approximately 270 h.This result is fully compatible with those found by Jorda et al. (2016) with the PDM technique. Data phased to this periodicity are shown in Fig. 12 of Jorda et al. (2016).
In the text

	Fig. 4 Body frame (xyz) and inertial frame (XYZ) with the Euler angles and their velocities in the standard ZXZ convention. This figure is adapted from Fig. A.1 from Samarasinha & Mueller (2015). Reproduced with permission.
In the text

Fig. 5 a) Bullets (•) show the spin axis orientation as determined with the SPC method (Fig. 2). Small dots (·) are the spin vector orientation for a simulated body with inertia moments equal to those derived from the shape of 67P assuming homogeneity, and an El = 1.9 × 10-6 under case 1 circumstances. These data are displayed with a time separation of 1 h. Diamonds (◆), to be compared with the bullet distribution, correspond to the orientation of the spin axis for the simulation performed at the time and date associated with each determination made with the SPC method. To locate the diamonds, neither noise nor averaging have been included. b) Lomb periodogram of the Dec-simulated coordinates (diamonds) shown in Fig. 5a. The RA coordinate periodogram is exactly the same as that of the Dec data, only small differences in spectral power (height of the peak) can be appreciated.

In the text

Fig. 6 a) As in Fig. 5, bullets (•) show the spin axis orientation as determined withthe SPC method (Fig. 2). Small dots (·) are the spin vector orientation for a simulated body with inertia moments equal to those derived from the shape of 67P assuming homogeneity, and an El = 1.9 × 10-6 under case 1 assumptions , that is, exactly the same circumstances as those of the simulation shown in Fig. 5a. These data are displayed with a time separation of 1 h. Diamonds (◆), to be compared with the bullet distribution and also with Fig. 5a, correspond to the simulated spin axis orientation at the moments associated with each determination with the SPC method, but now calculated as the average value of the instantaneous coordinates within intervals of 10 h around these moments. b) Lomb periodogram of the simulated averaged Dec coordinates (diamonds) shown in Fig. 6a.

In the text

Fig. 7 a) As previously, bullets (•) show the observationally derived spin axis (Fig. 2). Small dots (·) are the spin axis orientation for a simulated body with inertia moments 1:1.90:1.905 and an El = 2.9 × 10-6 (see text). The data are displayed with a time separation of 1 h. As previously, diamonds (◆), to be compared with the bullets distribution, correspond to the averaging within 10 h of the instantaneous theoretical values. b) Lomb periodogram of the Dec-simulated data (diamonds) shown in Fig. 7a. The x-axis here ranges up to 300 h to show that no significant periodicity is detected beyond Pc = 13 h. Given the wide range covered by the x-axis, the two significant periodicities Pφ = 11.85 h and Pc appear very close in the figure.

In the text

	Fig. 8 Relationship between Iz and Iy for a Pc of approximately 270 h for case 1 (continuous line) and case 2 (dashed line). Both Iy and Iz are normalized to Ix. The asterisk (∗) represents the homogeneous body.
In the text

Fig. 9 a) K-S probability as a function of Iy and El when the theoretical spin axis RA distribution is compared to the observationally derived one for case 1. b) The same as a), but in this case the theoretical and observationally derived Dec distributions are compared. In the two figures, the continuous white line, El,RA = El [ Iy ], indicates the excitation level for each Iy associated with the highest K-S probability when RA distributions are compared. The dashed white line, El,Dec = El [ Iy ], represents the excitation level for each Iy associated with the highest K-S probability when the Dec distributions are compared.

In the text

Fig. 10 Simulation with the highest K-S probability from the comparison of RA coordinate distributions. a) As previously, bullets (•) show the observationally derived spin axis (Fig. 2). Diamonds (◆) correspond to the simulated spin axis orientation of a body with 1:1.90:2.84, and El = 6.7 × 10-5 when an averaging of 10 h and a uniformly distributed random noise with an amplitude of 0.03◦ have been considered. b) (Thick lines) Lomb periodograms of the RA and Dec simulated data (diamonds shown in Fig. 10a). The Lomb periodograms of the observationally derived coordinates (Fig. 3) have also been included (with thin lines).

In the text

Fig. 11 Simulation with the highest K-S probability from the comparison of Dec coordinate distributions. a) As previously, bullets (•) show the observationally derived spin axis (Fig. 2). Diamonds (◆) correspond to the simulated spin axis orientation of a body with 1:2.25:3.19, and El = 3.1 × 10-5 when an averaging of 10 h and a uniformly distributed random noise with an amplitude of 0.03◦ have been considered. b) (Thick lines) Lomb periodograms of the RA and Dec simulated data (diamonds shown in Fig. 11a). The Lomb periodograms of the observationally derived coordinates (Fig. 3) have also been included (with thin lines).

In the text

	Fig. 12 K-S probability estimated as the mean of the RA and Dec probabilities for the two excitation curves, El = El [ Iy ], leading to the highest probability of the K-S tests performed on RA (•) and Dec (⋄) data separately.
In the text

	Fig. 13 Lomb periodograms of the observationally derived RA (black) and Dec (Red) coordinates (Fig. 3) for periodicities lower than 14 h. Arrows indicate the position of Pφ for the two possibilities discussed in the text, 6.35 h and 12.40 h.
In the text