© The Authors 2026

1 Introduction

Comets exhibit intense surface activity near perihelion, when volatile sublimation drives the release of gas and dust as jets. Beyond redistributing surface material, the resulting net forces and torques exert a measurable influence on the long-term orbital and rotational evolution of comets (Jia et al. 2017; Kokotanekova et al. 2017; Whipple 1950). Consequently, non-gravitational acceleration (NGA) modelling has long been a central component of cometary orbit determination, evolving from empirical formulations to more physically motivated and semi-empirical approaches (Marsden et al. 1973; Yeomans & Chodas 1989; Chesley & Yeomans 2005). Although most mainstream NGA models share a common physical basis in sublimation-driven outgassing and similar parameterisations of heliocentric distance dependence, they differ in how they represent the spatial distribution of activity on the nucleus, the temporal response of outgassing, and the incorporation of rotational or geometric effects (Sekanina 1981; Gutiérrez & Davidsson 2007). When embedded in orbit determination frameworks, these assumptions propagate into the estimation of non-gravitational forcing, such that reconstructed ephemerides can encode systematically different effective descriptions of outgassing-driven effects.

This issue is particularly relevant for comet 67P/ChuryumovGerasimenko (hereafter 67P). Following the comprehensive operational orbit determination by the European Space Operations Centre (ESOC) (Godard et al. 2015), three post-Rosetta ephemerides were independently reconstructed using different NGA models as well as different selections and weighting of observational constraints (Farnocchia et al. 2021; Kramer & Läuter 2019; Manghi et al. 2025). As a result, these ephemerides cannot be straightforwardly compared based solely on their post-fit observational residuals. Because such metrics are inherently conditioned on the adopted observation sets, weighting schemes, and NGA models, they may reflect distinct dynamical representations of non-gravitational forcing (Carpino et al. 2003; Desmars et al. 2015). These considerations motivate an ephemeris-based analysis in which the ephemerides themselves are treated as the primary source for inferring robust properties of the implied NGA. Hereafter, these ephemerides are referred to as K & L (2019), Farnocchia et al. (2021), and Bologna (2024), respectively. Recently, Attree et al. (2024a) adopted a forward modelling approach using NGA derived from the Bologna ephemeris to interpret water-driven outgassing activity through spatially variable momentum transfer efficiencies and effective active fractions of the cometary surface.

In this work, we adopted a reconstruction-oriented approach to investigate the physical consistency of the post-Rosetta orbital solutions. Rather than introducing a new outgassing model, we extracted the effective NGA directly from the ephemerides to decode the dynamical signatures embedded within them. Our analysis focuses on isolating the intrinsic geometry of the non-gravitational forces, specifically their alignment with the nucleus’s spin axis and surface features, to determine whether independent reconstructions yield a coherent picture of the cometary activity despite their methodological differences.

To provide the necessary spatial context for the directional analysis of non-gravitational forces, Figure. 1 presents the high-resolution shape model (Chen et al. 2023), visualising the complex geometry of the nucleus within the standard Cheops reference frame. The surface is coloured according to the Z-coordinate height, visually delineating the contrast between the northern (blue) and southern (red) hemispheres.

Building upon this geometric framework, we examined the effective latitude of the inferred acceleration (the angle between the net non-gravitational force vector and the comet equatorial plane), which condenses the three-dimensional acceleration direction into a measure of hemispheric dominance. This quantity enables a direct connection between reconstructed non-gravitational forcing and the well-established southern hemisphere dominance of surface activity on 67 P around and after perihelion. Within this framework, we assessed to what extent the different ephemerides capture this hemispheric asymmetry, and whether the inferred effective latitude evolution is compatible with independent constraints from Rosetta remote-sensing and in situ observations. Ultimately, this analysis highlights the utility of effective latitude not merely as a comparative metric, but as a physical link that anchors the reconstructed dynamics to the observed seasonal evolution of the cometary surface.

Fig. 1 Orthographic projections of the nucleus of comet 67 P, based on the shape model of Chen et al. (2023). The model is shown in the Cheops reference frame, with the Z-axis aligned with the rotation axis. Left: equatorial view showing the bilobate structure. Centre: south polar view (−Z), emphasising the southern hemisphere. Right: north polar view (+Z). The colour scale represents the Z-coordinate, distinguishing the southern (warm colours) and northern (cool colours) hemispheres.

2 NGA models and post-Rosetta ephemerides of comet 67P

We first reviewed the principal approaches that were used to model NGA in cometary orbit determination, establishing the necessary context for the comparative analysis presented in this work. Importantly, differences in modelling philosophy can lead to ephemerides that encode systematically different representations of non-gravitational forcing, even when constrained by similar observational data. This motivates a careful examination of both the modelling frameworks and the resulting reconstructed ephemerides.

2.1 Overview of NGA modelling approaches

(1) Empirical standard models. The earliest empirical formulations, established by Marsden et al. (1973), are grounded in Delsemme’s (1973) ‘dirty snowball’ hypothesis. This approach, known as the standard model (SM), correlates the NGA with the water ice sublimation flux. The acceleration vector a_NG,SM is expressed in the heliocentric radial-transverse-normal (RTN) frame: $${\mathit{a}}_{\mathrm{N}\mathrm{G},\mathrm{S}\mathrm{M}}=(A_1{\mathit{e}}_r+A_2{\mathit{e}}_t+A_3{\mathit{e}}_n)g(r),$$(1)

where A₁, A₂, and A₃ are constant coefficients. The basis vectors are defined as follows: e_r is the radial unit vector pointing from the Sun to the comet, e_n is normal to the orbital plane and aligned with the orbital angular momentum, and e_t completes the right-handed triad in the orbital plane. The function g(r), which describes the heliocentric dependence of the water sublimation flux, is $$g(r)=\alpha {\left(\frac{r}{r_0}\right)}^{-m}{[1+{\left(\frac{r}{r_0}\right)}^n]}^{-k},$$(2)

where r is the heliocentric distance, r₀ is the scale distance, and α is a normalisation constant such that g(1 AU) = 1.

This empirical parameterisation was necessitated by the observational constraints of the 1970s. Since ground-based astrometry could not resolve the nucleus or determine its spin axis orientation, physical models linking outgassing jets to a nucleus-fixed frame were impractical for orbit determination at that time. Consequently, the SM relies on the heliocentric RTN frame, which effectively takes the specific surface activity into global orbital perturbations.

The SM assumes that cometary outgassing, and hence the NGA, reaches its maximum at perihelion. To account for the commonly observed asymmetry of comet activity with respect to perihelion, Yeomans & Chodas (1989) introduced a time-shift parameter Δt, by replacing the heliocentric distance dependence g(r) with g(r′), where r′ = r(t + Δt). This modification, known as the extended standard model (ESM), is now widely adopted in Jet Propulsion Laboratory (JPL) Horizons System comet orbit solutions.

While such empirical formulations have proven effective for orbit prediction, their parameters are not uniquely tied to the spatial or directional properties of surface activity, making it difficult to directly relate fitted coefficients to the geometry of outgassing on the nucleus. For this reason, ESM-type solutions are not used here as primary cases for nucleus-fixed geometric interpretation, which is instead reserved for the three post-Rosetta reconstructed ephemerides analysed below.

(2) Physically motivated models. Studies of discrete outgassing sources on cometary nuclei predate modern spacecraft missions. Early work by (Sekanina 1986, 1988a,b) introduced the concept of localised jets on a rotating nucleus to interpret observed coma morphology and to relate activity patterns to the geometry of solar illumination and nucleus rotation. Subsequent developments adapted the rotation-averaged effects of such jets for cometary orbit determination (Szutowicz 2000). These efforts culminated in what is commonly referred to as the rotating jet model (RJM), in which discrete active regions give rise to a rotation-averaged NGA acting on the cometary nucleus. Related jet-based concepts were also applied in the development of spacecraft dynamics and navigation models (Byram et al. 2007).

In orbit determination implementations of the RJM, the NGA is obtained by summing the contributions of multiple discrete sources and averaging the resulting instantaneous thrust over one rotation period (typically a few hours). A commonly used rotation-averaged form can be written as $${\overline{\mathit{a}}}_{\mathrm{N}\mathrm{G},\mathrm{P}}=g\left(r^{\mathrm{\prime }}\right)\sum _{i=1}^n{A}_{J_i}{J}_i,$$(3)

with the g(r′) function still being applied, A_Ji denotes the a priori intensity of the i-th active source, and J_i is the rotation-averaged thrust vector associated with the i-th jet, incorporating the effects of illumination geometry over one nucleus rotation. The vectors J_i are expressed in the nucleus-aligned SQP reference frame, the P-axis is aligned with the nucleus spin pole (pointing north), the S-axis is given by the projection of the Sun direction onto the equatorial plane, and the Q-axis completes the right-handed triad. This formulation highlights the key distinction from empirical models such as the ESM, the acceleration is explicitly linked to the geometry of discrete active regions and the nucleus rotation, although it still requires assumptions about jet locations, activity patterns, and parameterisations.

(3) Soft physics-constrained stochastic models. Due to the intrinsic complexity of cometary outgassing processes and the limited resolution of available observations, purely physics-driven NGA models often struggle to capture the full range of activity variability. In practice, both standard empirical formulations and physically motivated jet models can become under-constrained when observational coverage is sparse or when surface activity evolves rapidly.

The Bologna (2024) ephemeris introduces semi-empirical stochastic approaches to address these limitations. A key development among these is the piecewise linear stochastic model (LSM), which has been applied to the reconstruction of the comet 67P orbit. Rather than prescribing a fixed functional form for NGA, LSM approximates the effective acceleration as a sequence of locally linear segments defined over adaptive time intervals. These intervals are typically chosen to reflect the expected temporal variability of surface activity, with shorter segments near perihelion where outgassing intensifies.

A key feature of this class of models is the introduction of soft mathematical constraints that ensure the continuity and smoothness of the reconstructed acceleration over adjacent time intervals. Such constraints prevent non-physical discontinuities while maintaining sufficient flexibility to account for rapid changes in the inferred non-gravitational forces. As a result, they provide good performance for navigation, at the cost of partially, at least, disregarding the underlying physics. Artificial intelligence orbital models (Embaby et al. 2025), which are currently the subject of intensive research, belong to this family of stochastic models. Therefore, additional analyses are needed to extract physically meaningful information from the inferred NGA.

2.2 Overview of the three ephemerides

After ESOC released the first reconstructed post-Rosetta ephemeris for comet 67P (CORB_DV_257_03___T19_00345.BSP), three independent teams subsequently developed updated orbital ephemeris using different strategies. Figure 2 shows the differences in the trajectories of 67P relative to the solar system barycentre (SSB), decomposed into an SSB-centric RTN frame where the radial vector points from the SSB to the comet.

Overall, the ephemerides from Bologna (2024) and K & L (2019) exhibit relatively small mutual differences across all SSB-centric RTN components. In contrast, the Farnocchia et al. (2021) ephemeris shows substantially larger discrepancies with respect to the other two solutions. Among the three directions, the normal (N) component displays the largest differences in both position and velocity. The epoch of the maximum position difference does not occur at any clearly preferred time, whereas the maximum velocity difference is distinctly concentrated around perihelion.

An investigation into the origin of these discrepancies is conducted by comparing the reconstruction strategies, NGA modelling approaches, and observational constraints across the three ephemerides. More details are shown in Table A.1.

K & L (2019) did not perform a conventional orbit determination based on the processing of raw tracking or astrometric observations. Instead, their approach relied on post-processing the ESOC multi-arc ephemeris. Specifically, as detailed in their work, they fitted the residual deviations between ESM and the ESOC SPICE kernel positions using a combination of exponentials and polynomials up to the fourth order. Therefore, while their solution is empirical, it represents a distinct formulation rather than a pure Marsden SM. This strategy was primarily designed to rapidly characterise the global properties and temporal evolution of NGA, rather than to generate a fully self-consistent ephemeris rigorously anchored to the measurements. Because Kramer & Läuter (2019) presented an ensemble of trajectories based on varied initial conditions, the specific trajectory analysed in this work is a representative high-resolution member of this ensemble. The SPICE kernel used here corresponds to exemplary trajectory #1 and is now publicly available via Zenodo (Kramer & Läuter 2026).

Farnocchia et al. (2021) reconstructed the orbit of comet 67P over a full orbital revolution from approximately 2012–2018 by fitting a combination of ground-based optical astrometry and Rosetta-derived pseudo-range data. A key aspect of their approach was the selection of an appropriate NGA model through prediction-based diagnostics. Several candidate models, including a gravity-only baseline and gravity models augmented with either the empirical sublimation model or the RJM, were tested by fitting pre-Rosetta optical-only datasets and evaluating their predictive performance against Rosetta-era measurements that were excluded from the fitting intervals. Based on the stability and accuracy of these predictions, the RJM with a heliocentric-distance dependence g(r′) ∝ r⁻⁵ (Hansen et al. 2016) was selected and used in the final reconstruction. While this ephemeris provides a smooth orbit over the 2012–2018 interval, the perihelion passage itself lacks observational coverage, with the non-gravitational behaviour near peak activity inferred solely through the adopted model.

Bologna (2024) reconstructed the ephemeris of comet 67P during Rosetta’s proximity phase by processing high-precision radiometric tracking data, specifically two-way X-band range and delta-differential one-way ranging observations acquired by the European Space Agency (ESA) through its European Space Tracking (ESTRACK) network and by the National Aeronautics and Space Administration (NASA) through its Deep Space Network (DSN) between 2014 and 2016. The deltadifferential one-way ranging technique, which employs alternating observations of the spacecraft and nearby reference quasars to measure differential group delays, provided high-precision constraints on the spacecraft’s plane-of-sky position, thereby significantly reducing uncertainties in the comet’s ephemeris during proximity operations. Ground-based optical astrometry was incorporated as a complementary dataset, with measurements retrieved from the Minor Planet Center (MPC) that had not undergone dedicated high-precision reprocessing excluded near perihelion, and higher-accuracy observations retained following Gaia-based astrometric bias corrections. The NGA was modelled using LSM, in which outgassing-induced accelerations are treated as time-correlated stochastic processes subject to physical regularisation and heliocentric-distance scaling. In contrast to strategies where a fixed physical model is pre-selected based on predictive diagnostics, the LSM parameters were estimated simultaneously with the orbit by fitting the full observational dataset, allowing NGA to be rigorously anchored to high-precision radiometric data during the perihelion interval.

Fig. 2 Differences between the three independently reconstructed post-Rosetta ephemerides of comet 67P, expressed in the SSB-centred RTN frame. Left panels: position differences. Right panels: velocity differences. Vertical dashed lines mark the equinoxes, perihelion, and solstice. The top axis gives the heliocentric distance in AU.

2.3 Extraction of effective NGA

The effective NGA, ${\mathit{a}}_{\text{NG}}^{\text{eff}}$, was extracted from the reconstructed trajectories by subtracting the reference gravitational acceleration (Sun and planets) from the total acceleration. This process is defined as $${\mathit{a}}_{\mathrm{N}\mathrm{G}}^{\mathrm{e}\mathrm{f}\mathrm{f}}(t)={\mathit{a}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}(t)-{\mathit{a}}_{\mathrm{G}}^{\mathrm{r}\mathrm{e}\mathrm{f}}(t),$$(4)

where a_total is obtained by differentiating the Chebyshev interpolation polynomials in the BSP kernels. We adopted JPL DE440 as the reference planetary ephemeris, consistent with the Gaia Celestial Reference Frame 3 (CRF3) and informed by modern spacecraft tracking datasets.

Since the post-Rosetta reconstructed trajectories were published over a span of years, the source datasets rely on different planetary ephemeris standards (i.e. DE430 and DE440). To validate the use of a unified DE440 baseline for all datasets, we evaluated the sensitivity of the extraction to these version differences. As detailed in Table 1, the total acceleration discrepancy between DE430 and DE440 is approximately 10⁻¹⁸–10⁻¹⁵ au d⁻². This difference is more than five orders of magnitude smaller than the typical outgassing-driven signal near perihelion (~10⁻¹⁰ au d⁻²), confirming that the choice of reference ephemeris has a negligible impact on the extraction results.

Additionally, solar radiation pressure (SRP) is explicitly omitted from our reduction model. Due to the comet’s high mass (~9.982 × 10¹³ kg; Pätzold et al. 2016), the magnitude of SRP (~10⁻²³ au d⁻²) is significantly below the extracted non-gravitational signal.

Fig. 3 Magnitude of the extracted NGA for comet 67 P as a function of time from perihelion. Results are shown for the three reconstructed ephemerides: K & L (2019), Farnocchia et al. (2021), and Bologna (2024). The dashed curve represents the sublimation law g(r′), shifted by 23 days to align with the summer solstice. The shaded regions indicate the active window defined by $\left|a_{\text{NG}}^{\text{eff}}\right|\ge 0.1A_{\text{max}}$.

3 Comparative analysis of NGA

3.1 Magnitude and temporal evolution

The overall intensity of cometary outgassing is characterised by examining the effective NGA magnitude $\left|a_{\text{NG}}^{\text{eff}}\right|$ extracted at a daily cadence over a ±350-day window around perihelion. As illustrated in Figure. 3, all three reconstructed ephemerides exhibit a consistent global trend with a pronounced rise towards perihelion followed by a gradual decay. To quantify deviations from standard heliocentric-distance dependence, we compared the profiles against a modified sublimation reference curve g(r′) ∝ r⁻⁵ with a fixed time delay of 23 days, which aligns with the southern summer solstice.

While the gross behaviour is similar, examining the component-wise evolution in Figure. 4 reveals some systematic differences. The Bologna (2024) and K & L (2019) solutions exhibit remarkable consistency across the radial (a_R), transverse (a_T), and normal (a_N) directions. Both their curve trajectories and overall trends align closely, indicating a shared dynamical representation despite their fundamentally different reconstruction methodologies. Conversely, the Farnocchia et al. (2021) ephemeris presents distinct morphological differences. In the transverse (a_T) direction, its curve diverges noticeably, exhibiting a more pronounced pre-perihelion positive surge followed by a deep post-perihelion negative trough. Furthermore, its peak magnitude in the normal (a_N) direction is substantially amplified compared to the other two solutions. Such discrepancies suggest that the prediction-driven physical model (RJM) adopted by Farnocchia et al. (2021) yields a systematically different spatial distribution of the non-gravitational forcing near perihelion.

We defined the active window duration (T_0.1) as the total time interval during which the acceleration magnitude exceeds 10% of its peak value ($\left|a_{NG}^{eff}\right|\ge 0.1A_{max}$); it serves to quantify the time span of significant outgassing. Based on this metric, all three ephemerides suggest a compressed activity phase of 180–208 days, noticeably shorter than the ~229 days predicted by the delayed theoretical reference curve. Within the T_0.1 window, the profiles exhibit a slight tendency towards inbound-outbound asymmetry. We characterised this morphology using the asymmetry index, defined as $$Asym=\frac{W_R-W_L}{W_R+W_L},$$(5)

where W_L and W_R are the temporal full widths at half maximum (FWHM) before and after perihelion, respectively. A positive index indicates a ‘long-tail’ behaviour. As shown in Table A.2, all solutions yield positive values, confirming a more gradual decline of activity post-perihelion. This effect remains most pronounced in the K & L (2019) solution (Asym = 0.11), which maintains significant activity up to approximately 145 days post-perihelion. In contrast, the Farnocchia and Bologna solutions converge to a faster shutdown, reaching their 10% activity thresholds at around 118 days post-perihelion.

We captured the rates of activity onset and decline by fitting exponential functions of the form e^|t|/τ to the rising and decaying branches of the profiles. The resulting decay timescales (τ_decay) reveal a clear divergence during the post-perihelion decay phase. While the Farnocchia and Bologna ephemerides exhibit faster decline rates clustered around 55–59 days, the K & L (2019) solution points a significantly slower decay with a timescale of approximately 68 days. This extended post-perihelion activity is further underscored by its higher persistence ratio (τ_ratio = 1.35). Conversely, the pre-perihelion growth timescales (τ_rise) show tight agreement across all three ephemerides, ranging narrowly from approximately 50 to 52 days.

Overall, the intensity analysis reveals that the three ephemerides remain embedded within a common non-gravitational activity regime. Although small but measurable differences exist in the detailed morphology of the profiles, including the duration, asymmetry, and post-perihelion decay behaviour, these variations do not alter the broader conclusion that all solutions recover a similar acceleration scale and a consistent level of sustained activity around perihelion. This suggests that the dominant outgassing signature is robustly captured across the different reconstruction strategies, while the residual discrepancies are mainly confined to secondary temporal characteristics. Consequently, the trajectory-level differences discussed in Sect. 2.2 are unlikely to be controlled primarily by the absolute strength of the reconstructed force. Rather, they are more plausibly associated with differences in its directional partitioning and geometric evolution, which motivates the component-wise analysis in the following sections.

Fig. 4 Time evolution of the NGA components in the heliocentric RTN frame for the three post-Rosetta ephemerides. The radial (aR), transverse (aT), and normal (aN) components are shown as functions of time relative to perihelion.

3.2 Decomposition into heliocentric RTN and nucleus-aligned SQP frames

To investigate the directional structure of the extracted acceleration, we decompose the vectors into the heliocentric RTN frame (as defined in Sect. 2.1) in Figure. 4. While the magnitude evolution mimics the perihelion surge described in Sect. 3.1, the component-wise behaviour reveals a lack of geometric coherence. For all three solutions, the acceleration is distributed significantly across all three axes (R, T, N) without a stable, singular preferred direction. These features suggest that the heliocentric RTN frame, while natural for orbital dynamics, poorly represents the intrinsic geometry of the outgassing force, as it convolves the cometary activity with the evolving orbital aspect angle.

A strikingly more coherent picture emerges when the same vectors are projected into the nucleus-aligned SQP frame, which is aligned with the cometary spin axis, as shown in Figure. 5. In contrast to the heliocentric RTN scattering, the nucleus-aligned SQP decomposition uncovers a unifying structural property across all three ephemerides, the clear dominance of the P-component (polar direction) during the active phase. Across the solutions, the S and Q equatorial components remain relatively small, implying that the net thrust is primarily directed along the nucleus rotation axis. While the K & L (2019) ephemeris exhibits slightly more pronounced evolution in the S-components compared to the Farnocchia et al. (2021) and Bologna (2024) solutions, the P-component consistently captures the bulk of the acceleration magnitude for all models. This structural alignment validates the physical interpretation that the effective NGA is driven by sublimation from active regions that are effectively averaged or distributed such that the net force aligns with the spin axis.

The normalised P-contribution was introduced to quantify this directional preference: $$f_P(t)=\frac{\left|a_P\right|}{\parallel {\mathit{a}}_{\mathrm{S}\mathrm{Q}\mathrm{P}}\parallel },$$(6)

The f_P(t) measures the fraction of the total acceleration projected onto the spin axis. As shown in Figure 6, this metric confirms a highly stable force orientation across the different models. Throughout the active window, the K & L (2019) and Farnocchia et al. (2021) ephemerides yield high mean f_P values of 0.909 and 0.939, respectively. Similarly, the Bologna (2024) solution maintains a strong polar alignment with a mean f_P of 0.861. Rather than highlighting a fundamental divergence, this widespread agreement on P-axis dominance with all mean values remaining robustly high, demonstrates that the post-Rosetta ephemerides consistently encode a polar-aligned outgassing force, independent of the reconstruction strategy. Consequently, the nucleus-aligned SQP frame provides a physically interpretable framework that naturally decouples the intrinsic outgassing geometry from complex orbital kinematics.

Fig. 5 Same NGA vectors as in Figure. 4 but expressed in the nucleus-aligned SQP frame. The three panels show the S, Q, and P components, respectively.

Fig. 6 Temporal evolution of the normalised P-component contribution, fP(t), for the three ephemerides. The shaded regions denote the active window defined by the T0.1 threshold applied to the acceleration magnitude.

Fig. 7 Temporal evolution of the angles between the NGA vector and the basis vectors of the heliocentric RTN frame (left panels) and the nucleus-aligned SQP frame (right panels) for the three post-Rosetta ephemerides. Results are shown as functions of time relative to perihelion. Vertical dashed lines mark the equinoxes, perihelion, and summer solstice.

3.3 Directional properties of the acceleration

To further clarify the directional coherence of the extracted NGA and to compare the suitability of the two reference frames, we examined the time evolution of the angles between the NGA vector and the basis vectors of the heliocentric RTN and nucleus-aligned SQP frames. Figure 7 illustrates these angle variations as a function of time relative to perihelion.

In the heliocentric RTN frame, the angular behaviour during the perihelion phase exhibits a limited degree of coherence and lacks a clear geometric organisation. For all three ephemerides, the angles with respect to the R, T, and N directions exhibit substantial fluctuations near perihelion, with the T and N component showing particularly severe variations. This indicates a high geometric instability of the acceleration direction within this frame during the interval of strongest activity. Crucially, throughout this phase, the acceleration vector does not strictly align with any single cardinal axis (i.e. the angles do not approach 0° or 180° for any component). However, as observed in the component magnitudes (Figure 4), the R acceleration remains marginally dominant. This behaviour is physically expected, as outgassing is primarily driven by direct solar insolation on the sunward hemisphere, generating a strong antisunward reaction force. Nevertheless, despite this radial bias, the overall vector does not maintain a stable, singular alignment within the RTN frame. Consequently, while the heliocentric RTN frame reflects the temporal variability and the physically expected radial dominance of the NGA, it does not provide an optimal representation of its intrinsic directional structure during perihelion.

A markedly different and more physically coherent picture emerges in the nucleus-aligned SQP frame. For all three ephemerides, the angle between the NGA vector and the P axis decreases sharply around perihelion, reaching a distinct minimum during the active phase, while the corresponding angles with respect to the S and Q axes remain closer to 90°. This demonstrates that the force vector effectively collapses towards the polar axis, a structural alignment absent in the heliocentric RTN representation.

Moreover, prior to perihelion (at approximately −300 days), the P-axis angle in the Farnocchia et al. (2021) and Bologna (2024) solutions is approximately 140°–150°. This represents a reversal in the acceleration orientation that is consistent with the expected physical evolution of the active subsolar region. The variations of the S- and Q-axis angles are mainly concentrated around the equinox epochs, exhibiting a seasonally coherent pattern that is consistently reproduced by the ephemerides. Such systematic geometric features confirm that the nucleus-aligned SQP frame naturally decouples the complex orbital kinematics from the intrinsic outgassing geometry.

Taken together, this angular analysis provides a geometric confirmation of the component-wise results presented in the previous section. It indicates that during periods of intense outgassing, the NGA is predominantly aligned with the cometary polar axis. This physically meaningful organisation is naturally expressed in the nucleus-aligned SQP frame, providing the geometric basis for isolating the seasonal drivers of cometary activity, which are the focus of Sect. 4.

Fig. 8 Evolution of the daily maximum cosine of the solar zenith angle, max(0, cos z), for three representative latitudes on the nucleus of comet 67P: the equator (0°, top), northern mid-latitude (+45°, middle), and southern mid-latitude (−45°, bottom). Vertical lines denote the equinoxes, perihelion, and southern summer solstice.

4 Constraints on non-gravitational forcing from acceleration geometry

Building upon the geometric characterisation presented in Sect. 3, we next examined the latitudinal evolution of the extracted acceleration to connect the dynamical solutions with the physical state of the cometary nucleus.

4.1 The seasonal driver: Insolation geometry

To interpret the hemispheric migration of the NGA, it is essential to first consider the solar forcing acting on the nucleus. Figure 8 presents the daily maximum cosine of the solar zenith angle (cos z) for three representative latitudes on the nucleus. This metric illustrates the extreme seasonality of comet 67P driven by its high obliquity (52°), a well-established fundamental driver of the nucleus’s activity and evolutionary dichotomy (Keller et al. 2015; Groussin et al. 2025).

The northern hemisphere (Figure 8, middle panel) receives peak insolation well before perihelion but undergoes a rapid decline in energy input, transitioning into polar night shortly after the inbound equinox. Conversely, the southern hemisphere (Figure 8, bottom panel) experiences a sharp transition from darkness to intense, continuous illumination around perihelion. This coincides with the comet’s closest approach to the Sun, resulting in a maximised energy deposition that drives the vigorous ‘southern summer’ activity. This insolation asymmetry provides the fundamental physical mechanism for the polar alignment observed in the nucleus-aligned SQP decomposition. As the subsolar point migrates south, the active regions responsible for the bulk of the outgassing shift from the northern to the southern hemisphere, anchoring the net force vector to the rotation axis.

Fig. 9 Temporal evolution of the effective NGA latitude as a function of time relative to perihelion for the three post-Rosetta ephemerides. The latitude of the anti-subsolar point is shown for comparison as a purely geometric reference. Vertical dashed lines mark the equinoxes, perihelion, and summer solstice.

4.2 Evolution of effective latitude

We defined the effective latitude (λ_eff) of the acceleration vector to assess how accurately the ephemerides encode the seasonal driver. The quantity was computed in the standard nucleus-fixed reference frame (Jorda et al. 2016), where the Z-axis aligns with the cometary rotational pole. In this frame, λ_eff represents the angle between the NGA vector and the equatorial plane. Figure 9 compares the temporal evolution of this quantity extracted from the three ephemerides against the geometric latitude of the antisubsolar point.

As expected, the effective latitude undergoes a global transition from the northern to the southern hemisphere, broadly tracking the migration of the Sun. However, the precise timing and persistence of this migration reveal significant divergences among the reconstructed solutions, which can be directly evaluated against independent observational constraints (Table 2).

The first critical divergence appears during the inbound transition. The Bologna (2024) ephemeris places the equatorial crossing at approximately 1.69 au (−96.73 days), and the K & L (2019) solution follows closely at 1.61 au (−85.98 days). In contrast, the Farnocchia et al. (2021) ephemeris shifts the transition closer to perihelion (1.53 au). Independent multi-instrument analyses of Rosetta data indicate that the dominant H₂O production migrates from northern to southern latitudes between approximately 1.7 and 1.5 au (Hansen et al. 2016). Consequently, the earlier transitions of both the Bologna and K & L solutions appear more consistent with the observed onset of water-driven activity in the southern hemisphere, while the Farnocchia et al. (2021) ephemeris implies a dynamical response that is delayed relative to this observational baseline.

An even more pronounced discrepancy emerges during the outbound phase. Both the Farnocchia et al. (2021) and the K & L (2019) ephemerides suggest a return of the NGA vector towards the equator. While Farnocchia et al. (2021) crosses back rapidly before the outbound equinox, the K & L (2019) solution delays the drop until after 120 days post-perihelion but still ultimately undergoes a sharp decline, crossing the equator at approximately 200 days. In striking contrast, the Bologna (2024) solution remains persistently associated with high northern latitudes (λ_eff > 40°) well beyond the equinox and throughout the plotted window. The persistence is strongly supported by spatially resolved compositional data. At larger heliocentric distances, cometary mass loss is dominated by CO₂ outgassing. This persistence is also substantiated by the column density maps obtained with the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) Reflectron-type Time-Of-Flight (RTOF) instrument (Hoang et al. 2020; Groussin et al. 2025). The data reveal that H₂O emissions are still heavily focused in the southern hemisphere within the 2.5–3.1 au range. Subsequently, at distances of 3.4–3.9 au, although the dominant volatile source transitions from H₂O to CO₂, the outgassing activity continues to originate predominantly from the southern regions, consistent with the high northern latitudes of the effective acceleration vector.

A distinct line of evidence supporting the asymmetry comes from thermophysical modelling of the nucleus surface evolution, which independently corroborates the persistence of southern activity. Benseguane et al. (2022) demonstrate that the intense insolation received at and near perihelion drives a distinct pattern of differential erosion. Their simulations show that, integrated over ten orbital revolutions, the maximum surface retreat reaches approximately 80 m in the southern hemisphere (Guilbert-Lepoutre et al. 2023), compared to generally less than 30 m in northern regions. The north and south dichotomy, which is primarily controlled by seasonal insolation patterns common to Jupiter-family comets, provides a permanent geological record of the mass loss asymmetry. The persistent northern orientation of the Bologna (2024) acceleration vector is therefore physically consistent with the reaction force expected from a region undergoing such significant long-term mass wasting.

Furthermore, comparing the continuous evolution of λ_eff with the discrete fitted parameters of the Farnocchia et al. (2021) RJM provides a useful consistency check on how different non-gravitational parameterisations encode the implied outgassing geometry. In the Farnocchia et al. (2021) solution, the activity is represented by predefined sources at specific latitudes (approximately +50.7° and −60.6°, corresponding to thrust colatitudes of 39.3° and 150.6°). While this RJM framework can reproduce the macroscopic orbital perturbations by fixing a small number of persistent source regions, our λ_eff(t) analysis indicates a more continuous, seasonally modulated migration of the effective thrust direction. Here we emphasise that λ_eff is not a directly observable quantity, but a derived geometric diagnostic computed from the reconstructed NGA vector embedded in an orbit solution. In this sense, the discrete RJM jet latitudes should be interpreted as temporally averaged proxies that anchor the dominant hemisphere of activity, rather than as instantaneous locations of an actively migrating sublimation front. This comparison is particularly relevant for the orbit determination of unresolved comets: when only astrometry (and/or limited radiometric data) is available, fitting a discrete RJM can still provide a physically interpretable, model-based inference on hemispheric dominance, whereas the time evolution of a derived quantity such as λ_eff(t) offers a compact way to summarise how the inferred thrust direction varies seasonally within the adopted model class.

The persistent alignment of the NGA direction with the comet’s rotation axis has been previously pointed in the literature, even if it was not always stated in the same geometric terms adopted here. For example, Attree et al. (2023) show a rough alignment between the radial and polar directions around perihelion, and Attree et al. (2024a) find that their best-fit solutions maintain a systematically higher effective active fraction (EAF) for southern-facing terrains and a lower EAF for northern-facing terrains after perihelion. Importantly, their framework varies the EAF itself (rather than attributing the fit solely to momentum-transfer efficiency), and the resulting south-dominated post-perihelion activity is consistent with their water-production constraints and with the spatially resolved compositional maps of Hoang et al. (2020) and Groussin et al. (2025). Similarly, Kramer & Läuter (2019) compared inferred non-gravitational behaviour against the subsolar-latitude evolution as a geometric proxy. Within this context, our contribution is to make the directional statement explicit and compact. By extracting the effective latitude λ_eff, we provide a complementary and direct geometric confirmation of the persistent post-perihelion southern dominance implied by these earlier analyses.

Regarding the volatile drivers, it is well established that H₂O outgassing dominates over CO₂ until a heliocentric distance of about 3 AU (Läuter et al. 2020; Groussin et al. 2025). Therefore, the sustained northern-pointing NGA during the initial post-perihelion phase is primarily driven by southern water activity. The simultaneously high CO₂ flux observed in the same regions is plausibly linked to intense water-driven erosion that brings the sublimation front closer to the surface, potentially aided by CO₂-driven ejections, as discussed in related studies (e.g. Attree et al. 2024b, 2025).

At distances beyond 3 AU, where CO₂ overtakes H₂O as the dominant volatile, the Bologna (2024) acceleration vector maintains its northern alignment. This behaviour hints that CO₂ outgassing contributes to the residual NGA in this regime. However, we caution against over-interpreting this signature. This phase lies well outside the primary active window defined by the NGA magnitude, so the extracted acceleration becomes very small and any inferred direction correspondingly less reliable. Moreover, this persistent alignment is not strongly reproduced in the K&L (2019) or Farnocchia et al. (2021) solutions. Fully disentangling the dynamical contribution of CO₂ at large heliocentric distances would therefore benefit from dedicated thermophysical modelling coupled to orbit determination, and should be regarded as a topic for future work rather than a firm conclusion of the present analysis.

5 Discussion and conclusions

We investigated the geometry and temporal evolution of the effective NGA derived from three post-Rosetta ephemerides of comet 67P, utilising a uniform, ephemeris-based diagnostic framework. Although these ephemerides were reconstructed using diverse NGA formulations, strategies, and observational constraints, our analysis reveals several consistent qualitative features across all solutions.

Most notably, during periods of intense activity, the inferred acceleration consistently exhibits a pronounced alignment with the cometary polar axis. Furthermore, its overall directional geometry follows a coherent seasonal modulation governed by heliocentric geometry. These shared characteristics indicate that the core directional structure of the NGA is robustly encoded in the reconstructed ephemerides, largely independent of the specific modelling approach employed.

Despite these broad similarities, the detailed directional behaviour of the three solutions is not identical. Systematic differences in the evolution of the effective latitude are primarily associated with the Bologna (2024) ephemeris, which displays a persistent hemispheric asymmetry after perihelion. In contrast, the other ephemerides show a more rapid return towards near-equatorial effective latitudes following the perihelion passage. This divergence highlights the sensitivity of directional diagnostics to reconstruction choices, particularly during phases where observational coverage and modelling flexibility differ most significantly.

The emergence of a clearer post-perihelion hemispheric asymmetry in the Bologna (2024) ephemeris should not be interpreted as evidence of shortcomings in more physically prescriptive NGA models. Rather, it reflects the combined influence of the continuous, high-precision observational constraints available during the mission, which directly constrained the perihelion phase, and a stochastic representation that allows greater freedom in capturing spatially asymmetric activity. Within the scope of the present analysis, it is not possible to unambiguously decouple the relative contributions of improved observational coverage versus increased model expressiveness. Nevertheless, this comparison demonstrates that directional quantities, such as effective latitude, are particularly effective in exposing reconstruction-dependent behaviours that remain hidden when analysing acceleration magnitude alone.

A key methodological outcome of this work is the demonstration that the nucleus-aligned SQP reference frame provides a superior representation of NGA geometry that is both compact and physically interpretable. By aligning one axis with the cometary spin direction, the nucleus-aligned SQP frame reduces the effective dimensionality of the problem, allowing directional diagnostics to directly reflect seasonal and hemispheric variations in activity. Compared with conventional heliocentric RTN decompositions, this representation concentrates physically meaningful directional information into fewer parameters, facilitating a clearer interpretation of the reconstructed non-gravitational forcing.

Finally, the comparative analysis shows that differences in the magnitude of NGA among post-Rosetta ephemerides are relatively modest and exhibit broadly similar temporal profiles. This suggests that the overall strength of outgassing is comparatively well constrained, whereas the distinctions among reconstructions are primarily encoded in the directional and hemispheric geometry of the inferred acceleration. Within the model-agnostic framework adopted here, directional diagnostics therefore provide a more sensitive and physically informative basis for comparing non-gravitational effects across heterogeneous ephemeris solutions. Taken together, these results demonstrate that the overall magnitude and polar alignment of the NGA are robust features encoded across all post-Rosetta solutions, whereas the post-perihelion evolution of the effective latitude, specifically the persistence of southern hemisphere activity, remains sensitive to the specifics of reconstruction strategies and data coverage.

Acknowledgements

We thank the Rosetta mission team for their hard work, which made this work possible. We are grateful to Tobias Kramer and Matthias Läuter for helpful correspondence regarding the ephemerides and for making the corresponding SPICE kernel publicly available via Zenodo (DOI: 10.5281/zenodo.19181638). This work is supported by the National Natural Science Foundation of China (grant No. 42241116) and the National Key Research and Development Program of China (grant No. 2022YFF0503202). J.Y. acknowledges support from the 2022 Project of Xinjiang Uygur Autonomous Region of China for the Heaven Lake Talent Program.

References

Appendix A Supplementary tables

All Tables

All Figures

Fig. 1 Orthographic projections of the nucleus of comet 67 P, based on the shape model of Chen et al. (2023). The model is shown in the Cheops reference frame, with the Z-axis aligned with the rotation axis. Left: equatorial view showing the bilobate structure. Centre: south polar view (−Z), emphasising the southern hemisphere. Right: north polar view (+Z). The colour scale represents the Z-coordinate, distinguishing the southern (warm colours) and northern (cool colours) hemispheres.

In the text

	Fig. 2 Differences between the three independently reconstructed post-Rosetta ephemerides of comet 67P, expressed in the SSB-centred RTN frame. Left panels: position differences. Right panels: velocity differences. Vertical dashed lines mark the equinoxes, perihelion, and solstice. The top axis gives the heliocentric distance in AU.
In the text

Fig. 3 Magnitude of the extracted NGA for comet 67 P as a function of time from perihelion. Results are shown for the three reconstructed ephemerides: K & L (2019), Farnocchia et al. (2021), and Bologna (2024). The dashed curve represents the sublimation law g(r′), shifted by 23 days to align with the summer solstice. The shaded regions indicate the active window defined by $\left|a_{\text{NG}}^{\text{eff}}\right|\ge 0.1A_{\text{max}}$.

In the text

	Fig. 4 Time evolution of the NGA components in the heliocentric RTN frame for the three post-Rosetta ephemerides. The radial (aR), transverse (aT), and normal (aN) components are shown as functions of time relative to perihelion.
In the text

	Fig. 5 Same NGA vectors as in Figure. 4 but expressed in the nucleus-aligned SQP frame. The three panels show the S, Q, and P components, respectively.
In the text

	Fig. 6 Temporal evolution of the normalised P-component contribution, fP(t), for the three ephemerides. The shaded regions denote the active window defined by the T0.1 threshold applied to the acceleration magnitude.
In the text

	Fig. 7 Temporal evolution of the angles between the NGA vector and the basis vectors of the heliocentric RTN frame (left panels) and the nucleus-aligned SQP frame (right panels) for the three post-Rosetta ephemerides. Results are shown as functions of time relative to perihelion. Vertical dashed lines mark the equinoxes, perihelion, and summer solstice.
In the text

	Fig. 8 Evolution of the daily maximum cosine of the solar zenith angle, max(0, cos z), for three representative latitudes on the nucleus of comet 67P: the equator (0°, top), northern mid-latitude (+45°, middle), and southern mid-latitude (−45°, bottom). Vertical lines denote the equinoxes, perihelion, and southern summer solstice.
In the text

	Fig. 9 Temporal evolution of the effective NGA latitude as a function of time relative to perihelion for the three post-Rosetta ephemerides. The latitude of the anti-subsolar point is shown for comparison as a purely geometric reference. Vertical dashed lines mark the equinoxes, perihelion, and summer solstice.
In the text