Outgassing induced acceleration of comet 67P/Churyumov-Gerasimenko

TL;DR

This study uses Rosetta mission data to analyze the non-gravitational acceleration of comet 67P, linking it to outgassing patterns and sublimation processes, and compares it with empirical models and in-situ water production measurements.

Contribution

It provides a detailed analysis of the non-gravitational acceleration of 67P using precise trajectory data, revealing the relationship with diurnal outgassing and surface illumination.

Findings

Non-gravitational acceleration correlates with subsolar latitude.

Marsden parameters partially explain the acceleration.

Orbit-based acceleration matches water production data.

Abstract

Cometary activity affects the orbital motion and rotation state due to sublimation induced forces. The availability of precise rotation-axis orientation and position data from the Rosetta mission allows one to accurately determine the outgassing of comet Churyumov-Gerasimenko/67P. We derive the observed non-gravitational acceleration of 67P directly from the Rosetta spacecraft trajectory. From the non-gravitational acceleration we recover the diurnal outgassing variations and study a possible delay of the sublimation response with respect to the peak solar illumination. This allows us to compare the non-gravitational acceleration of 67P with expectations based on empirical models and common assumptions about the sublimation process. We use an iterative orbit refinement and Fourier decomposition of the diurnal activity to derive the outgassing induced non-gravitational acceleration. The uncertainties of the data reduction are established by a sensitivity analysis of an ensemble of best-fit orbits for comet 67P. We find that the Marsden non-gravitational acceleration parameters reproduce part of the non-gravitational acceleration but need to be augmented by an analysis of the nucleus geometry and surface illumination to draw conclusions about the sublimation process on the surface. The non-gravitational acceleration follows closely the subsolar latitude (seasonal illumination), with a small lag angle with respect to local noon around perihelion. The observed minor changes of the rotation axis do not favour forced precession models for the non-gravitational acceleration. In contrast to the sublimation induced torques, the non-gravitational acceleration does not put strong constraints on localized active areas on the nucleus. We find a close agreement of the orbit deduced non-gravitational acceleration and the water production independently derived from Rosetta in-situ measurement.

Figures

Tables

Table 1. Table 1: Representative initial conditions for the orbit integration at JDB 2456897.71970 2456897.71970 2456897.71970 ( ≈ − 350 absent 350 \approx-350  days from perihelion), mean equator and equinox of the Earth J2000 frame. Additional 30 30 30 nearby initial conditions have been identified which differ on average less than 22 22 22  km over the ( − 400 , 400 ) 400 400 (-400,400)  days interval from the partly nonphysical ESOC orbit.

quantity	value	unit
$x$	$+1.332184491942$	au
$y$	$-2.770636585665$	au
$z$	$-1.613532002605$	au
$\dot{x}$	$0.0042333005604361959$	au/d
$\dot{y}$	$0.0074132927293021402$	au/d
$\dot{z}$	$0.0034828287065240359$	au/d

Table 2. Table 2: Non-gravitational acceleration parameters used by Horizons in conjunction with the standard parameters α = 0.1112620426 𝛼 0.1112620426 \alpha=0.1112620426 , k = 4.6142 𝑘 4.6142 k=4.6142 , m = 2.15 𝑚 2.15 m=2.15 , n = 5.093 𝑛 5.093 n=5.093 , r 0 = 2.808 subscript 𝑟 0 2.808 r_{0}=2.808 for g ​ ( r ) 𝑔 𝑟 g(r) in Eq. ( 3 ).

quantity	value	unit
$A_1$	$1.066669896245\times {10}^{-9}$	au/d2
$A_2$	$-3.689152188599\times {10}^{-11}$	au/d2
$A_3$	$2.483436092734\times {10}^{-10}$	au/d2
$\Delta t$	$35.07142$	d

Equations

$$\ddot{\text{@vec}r}=\text{@vec}{a}_{\mathrm{G}}+\text{@vec}{a}_{\mathrm{NG}}.$$

$$\text{@vec}{a}_{\mathrm{NG},\text{Marsden}}=g(r^{\prime })(A_1\text{@vec}{e}_r+A_2\text{@vec}{e}_t+A_3\text{@vec}{e}_n),$$

$$g(r^{\prime })=\alpha {\left(\frac{r}{r_0}\right)}^{-m}{\left(1+{\left(\frac{r}{r_0}\right)}^n\right)}^{-k}.$$

$$\text{@vec}{r}_{\mathrm{obs}}(t)=\text{@vec}{r}_{\text{Marsden}}(t)+\Delta \text{@vec}{r}_{\text{fit}}(t)$$

$$\text{@vec}{a}_{\mathrm{NG},\mathrm{obs}}={\ddot{\mathbf{r}}}_{\mathrm{obs}}(t)-\text{@vec}{a}_{\mathrm{G}}({\mathbf{r}}_{\mathrm{obs}}(t)).$$

$$\text{@tens}{R}_{\mathrm{com}\to \mathrm{equ}}=\left(\frac{\text{@vec}h}{|\text{@vec}h|},-\frac{\text{@vec}\omega \times \text{@vec}r}{|\text{@vec}\omega \times \text{@vec}r|},\frac{\text{@vec}\omega }{|\text{@vec}\omega |}\right),$$

$$\text{@vec}{F}^{\mathrm{bf}}({\lambda }_{\text{sun}})=\left(\begin{array}{c}{D}_{x0}-D_{x1}\sin {\lambda }_{\text{sun}}-D_{x2}\cos {\lambda }_{\text{sun}}+\dots \\ D_{y0}-D_{y1}\sin {\lambda }_{\text{sun}}-D_{y2}\cos {\lambda }_{\text{sun}}+\dots \\ D_{z0}-D_{z1}\sin {\lambda }_{\text{sun}}-D_{z2}\cos {\lambda }_{\text{sun}}+\dots \end{array}\right).$$

$$⟨\text{@vec}{F}^{\mathrm{equ}}⟩=\text{@tens}{R}_{\mathrm{com}\to \mathrm{equ}}{\int }_0^{T_{\mathrm{rot}}}\frac{\mathrm{d}t}{T_{\mathrm{rot}}}\text{@tens}{R}_{\mathrm{z}}(-\omega t)\text{@vec}{F}^{\mathrm{bf}}(-\omega t)$$

$$\text{@vec}{a}_{\mathrm{NG},\mathrm{obs}}=\frac{⟨\text{@vec}{F}^{\mathrm{equ}}⟩}{M}=\frac{1}{M}\text{@tens}{R}_{\mathrm{com}\to \mathrm{equ}}\left(\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right)$$

$$\left(\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right)$$

$$\text{@tens}{R}_{\mathrm{com}\to \mathrm{equ}}=(\text{@vec}P,\text{@vec}Q,\text{@vec}R)\text{@tens}{R}_{\mathrm{z}}(\Phi )\text{@tens}{R}_{\mathrm{x}}(-I)\text{@tens}{R}_{\mathrm{z}}(\pi -{\theta }_0)$$

$$(\text{@vec}{e}_{\mathrm{r}},\text{@vec}{e}_{\mathrm{t}},\text{@vec}{e}_{\mathrm{n}})=(\text{@vec}P,\text{@vec}Q,\text{@vec}R)\text{@tens}{R}_{\mathrm{z}}(-\nu )$$

$$\left(\begin{array}{c}{U}_{\mathrm{r}}\\ U_{\mathrm{t}}\\ U_{\mathrm{n}}\end{array}\right)=(\text{@vec}{e}_{\mathrm{r}},\text{@vec}{e}_{\mathrm{t}},\text{@vec}{e}_{\mathrm{n}}{)}^{-1}\text{@tens}{R}_{\mathrm{com}\to \mathrm{equ}}\left(\begin{array}{c}\cos {\varphi }_D\cos {\lambda }_D\\ \cos {\varphi }_D\sin {\lambda }_D\\ \sin {\varphi }_D\end{array}\right).$$

$$a_{\mathrm{NG},\mathrm{Marsden}}=\sqrt{A_1^2+A_2^2+A_3^2}g(r)\approx Q(r)v_{\mathrm{gas},\mathrm{av}}/M$$

$$Q$$

$$Q\ge {\int }_0^{T_{\mathrm{rot}}}\frac{\mathrm{d}t}{T_{\mathrm{rot}}{v}_{\mathrm{gas},\mathrm{av}}}\left|\text{@vec}{F}^{\mathrm{bf}}(-\omega t)\right|,v_{\mathrm{gas},\mathrm{av}}={\int }_{\text{surf}}\mathrm{d}S|\text{@vec}{v}_{\mathrm{gas}}(\text{@vec}{r}_{\mathrm{bf}},t)|.$$

$$Q\approx \frac{4\sqrt{\Delta }}{2\pi v_{\mathrm{gas},\mathrm{av}}}\mathrm{E}\left(\frac{-D_{x1}^2-D_{y1}^2-D_{z1}^2+D_{x2}^2+D_{y2}^2+D_{z2}^2}{\Delta }\right),$$

Full text

*Aims. * We derive the observed non-gravitational acceleration of 67P directly from the Rosetta spacecraft trajectory. From the non-gravitational acceleration we recover the diurnal outgassing variations and study a possible delay of the sublimation response with respect to the peak solar illumination. This allows us to compare the non-gravitational acceleration of 67P with expectations based on empirical models and common assumptions about the sublimation process.

*Methods. * We use an iterative orbit refinement and Fourier decomposition of the diurnal activity to derive the outgassing induced non-gravitational acceleration. The uncertainties of the data reduction are established by a sensitivity analysis of an ensemble of best-fit orbits for comet 67P.

*Results. * We find that the Marsden non-gravitational acceleration parameters reproduce part of the non-gravitational acceleration but need to be augmented by an analysis of the nucleus geometry and surface illumination to draw conclusions about the sublimation process on the surface. The non-gravitational acceleration follows closely the subsolar latitude (seasonal illumination), with a small lag angle with respect to local noon around perihelion. The observed minor changes of the rotation axis do not favor forced precession models for the non-gravitational acceleration.

*Conclusions. * In contrast to the sublimation induced torques, the non-gravitational acceleration does not put strong constraints on localized active areas on the nucleus. We find a close agreement of the orbit deduced non-gravitational acceleration and the water production independently derived from Rosetta in-situ measurement.

1 Introduction

Studying the non-gravitational acceleration of comets provides important insights into the sublimation of cometary ices. Earth-bound astrometry allows one to determine the orbital corrections arising from the sublimation of volatiles. Non-gravitational accelerations are often invoked to explain the orbital evolution of comets and even interstellar objects entering the solar system (Whipple (1950); Marsden et al. (1973); Micheli et al. (2018); Sekanina (2019)).

As pointed out by Yeomans et al. (2004) it is therefore of interest to compare widely used models with the precise data returned from spacecraft missions. For comet 67P Krolikowska (2003) provided an assessment of the Marsden parameters within the asymmetric outgassing model of Yeomans & Chodas (1989) before the arrival of Rosetta at the comet. The determination of accurate parameters from Earth bound telescopic observations requires to monitor the position of the comet at several apparitions.

The situation is different for the orbit of 67P as observed by Rosetta. Rosetta accompanied 67P for more than two years and provided measurements of the three-dimensional position vector with an accuracy better than 100 km (see Godard et al. (2017)). The discrepancy between Earth bound orbit prediction and the position where Rosetta encountered 67P $590$ days before perihelion requires to adjust the Marsden-type orbit by $2000$ km. From this it can be estimated that Rosetta provided an up-to twenty times higher accuracy compared to previous orbit determinations.

In conjunction with the precisely known rotation state of 67P we perform an attribution of the actual orbital changes to the sublimation activity across the nucleus. To account for the observed gas release of 67P requires to extend the non-gravitational acceleration models developed by Whipple, Marsden, and Sekanina. To this end we introduce a Fourier decomposition of the sublimation induced force and establish a general expression connecting the diurnal variations of the sublimation rate with the orbital evolution. This formalism simplifies the analysis by eliminating intermediate angles and emphasizes the lack of proportionality between total production rate and non-gravitational acceleration.

2 Orbit changes by non-gravitational acceleration

The equation of motion for the position vector $\@vec{r}$ of the cometary nucleus is described by contributions of solar system gravitational acceleration $\@vec{a}_{\mathrm{G}}$ and the additional non-gravitational acceleration $\@vec{a}_{\mathrm{NG}}$ (see e.g. Yeomans et al. (2004))

[TABLE]

In general this equation holds in any inertial system but we considered it in the Earth’s equatorial system. The gravitational part $\@vec{a}_{\mathrm{G}}=\@vec{a}_{\mathrm{G}}(\@vec{r},t)$ is evaluated at the momentary position vector and takes into account the acceleration due to all solar system planets and the Earth Moon, Pluto, Vesta, and Ceres. Additional corrections due to relativistic effects are ignored here, since they play only a minor role for 67P. Relativistic corrections could be added for other objects moving faster around the sun. The non-gravitational part is affected by the sublimation of volatiles on the nucleus across the surface. For given initial conditions $\@vec{r}(t_{0})$, $\dot{\@vec{r}}(t_{0})$ and a suitable model of the non-gravitational acceleration, the orbit of a comet can be integrated with high precision.

2.1 Marsden-type orbits

To understand the origin of the non-gravitational acceleration $\@vec{a}_{\mathrm{NG}}$ it is helpful to consider different reference frames. Starting from the icy snowball model introduced by Whipple (1950, 1951) to account for sublimation processes, Marsden and co-authors developed a powerful parametrization for the non-gravitational acceleration in a series of papers Marsden (1968, 1969, 1970); Marsden & Sekanina (1971); Marsden (1972); Marsden et al. (1973). Marsden expressed $\@vec{a}_{\mathrm{NG}}$ with respect to three right-handed orthogonal unit vectors with $\@vec{e}_{\mathrm{r}}$ pointing from the sun to the nucleus, $\@vec{e}_{\mathrm{n}}$ directed along the orbital angular momentum perpendicular to the orbital plane, and $\@vec{e}_{\mathrm{t}}$ being orthogonal to both, $\@vec{e}_{\mathrm{r}}$ and $\@vec{e}_{\mathrm{n}}$. By comparing 14 cometary orbits Marsden et al. (1973) derived the following model for the non-gravitational acceleration

[TABLE]

with the constant parameters $A_{1}$, $A_{2}$, $A_{3}$ and the empirical activity function

[TABLE]

The solar distance $r^{\prime}=r(t-\Delta t)$ includes a time-shift $\Delta t$ introduced by Yeomans & Chodas (1989) to account for asymmetries of the activity with respect to the perihelion as studied by Sekanina (1988). As we discuss later, $g(r)$ is not directly proportional to the sublimation flux as originally stated by Marsden et al. (1973), Eq. (4). The Marsden parameters provide an excellent, albeit empirical description of the non-gravitational acceleration. The integration of a Marsden-type orbit proceeds by solving Eq. (1) with the non-gravitational acceleration Eq. (2).

2.2 Orbit determination based on spacecraft data

To move beyond the Marsden model requires more detailed data from spacecraft missions or radar observations about the magnitude and direction of the observed non-gravitational acceleration in space. For 67P the Rosetta mission provides a two years data set for both, the rotation-axis orientation (see Kramer et al. (2019)) and the orbital evolution, but required an iterative process to obtain the three-dimensional non-gravitational acceleration with high precision. We started with the the multi-arc orbital solution from the flight dynamics team at ESOC, which is available as SPICE kernel (CORB_DV_257_03___T19_00345.BSP) and is returned from the Horizons system as position for 67P during the Rosetta mission. In the following we refer to it as $\@vec{r}_{\text{ESOC}}(t)$. Before and after the Rosetta mission, Horizons returns the solution 67P/K154/2 based on Marsden parameters with a discontinuous jump by $2000$ km into the Rosetta period. Besides this discontinuity in Horizons, further discontinuities in the ESOC position vector $\@vec{r}_{\text{ESOC}}(t)$ during the Rosetta mission exist, which do not allow one to obtain the acceleration by a second derivative of the position vector (see the discussion by Attree et al. (2019)).

An accurate estimate of the non-gravitational acceleration is tied to finding the best possible initial condition during a time with negligible cometary activity. We performed an iterative orbit refinement to identify the initial conditions which minimize the error norm with respect to the ESOC data in the time period $(-400,-200)$ days from perihelion (see Table 1). During this initial search only the gravitational acceleration was considered, thus we evaluated the cometary orbits with Eq. (1) with the setup $\@vec{a}_{\mathrm{NG}}=0$.

In the next step we added Marsden-type non-gravitational accelerations with the parameters used by Horizons (see Table 2 and the parameter estimation by Krolikowska (2003)) and we solved Eqs. (1), (2) to obtain $\@vec{r}_{\text{Marsden}}(t)$. We extended the integration to the years 1959-2022 (limited by two Jupiter encounters) and verified that our initial conditions lead to orbital solutions within $4400$ km with respect to the Horizons solution. We compared our results with the ones obtained by submitting the corresponding osculating elements to the Advanced Horizons Asteroid & Comet SPK web interface Horizons (2019). For the given initial conditions our Marsden-type solution (Fig. 1) is considerably closer to $\@vec{r}_{\text{ESOC}}$ than the Marsden solution given by Attree et al. (2019). Attree et al. (2019) reported a Marsden solution with a difference in Earth bound range $|\@vec{r}_{\text{ESOC}}-\@vec{r}_{\text{Earth}}|-|\@vec{r}_{\text{Marsden}}-\@vec{r}_{\text{Earth}}|$ of $1500$ km $400$ days after perihelion, while our Marsden solution at that time differs only by $50$ km. This highlights the necessity to perform a throughout statistical ensemble analysis of initial conditions to validate the non-gravitational part of the orbital acceleration. In the final iteration we fitted the remaining difference vector $\Delta\@vec{r}(t)=\@vec{r}_{\text{ESOC}}-\@vec{r}_{\text{Marsden}}$ to a combination of exponentials and polynomials up to fourth order, which provided (upon adding it to the Marsden solution) a differentiable representation of the observed orbit

[TABLE]

(see Fig. 1). Only after these iterative steps the non-gravitational acceleration was obtained from

[TABLE]

The validation of the second order derivative of the initially noisy position vector required a careful analysis of errors introduced by the fit. We have repeated the entire analysis by applying the smoothing and differentiation filter introduced by Savitzky & Golay (1964) with the corrections by Steinier et al. (1972) to the larger difference vector $\@vec{r}_{\text{ESOC}}-\@vec{r}_{0}$, where $\@vec{r}_{0}$ includes only gravitational forces. The ten times larger distance increases the fit uncertainties but leads to qualitative similar results in the period $\pm 100$ days from perihelion. Remaining systematic errors are discussed in Sect. 4.

The resulting non-gravitational acceleration is shown in Fig. 2. In combination with the initial condition (Table 1), the retrieved $\@vec{a}_{\mathrm{NG},\mathrm{obs}}$ reproduced the multi-arc solution $\@vec{r}_{\text{ESOC}}$ for $\pm 400$ days around perihelion with an mean error of $21$ km, mainly caused by the nonphysical jumps (see black curve in Fig. 1). At most times the error is around 10 km ($\approx 5$ cometary radii). To test the sensitivity of the orbit with respect to the initial conditions, we investigated an ensemble of 1000 nearby initial conditions and verified that no further improvement is obtained. Out of the 1000 initial conditions we identified $31$ orbits with a mean error $<22$ km in the $\pm 400$ days interval around perihelion, which all originate from a phase space volume extending about $3$ km around the initial position listed in Table 1 and with $<10^{-5}$ variations in the velocities. In the following we show results for this set of $31$ solutions in the form of shaded bands to estimate the uncertainties of derived quantities from the orbits.

3 Transformation of the non-gravitational acceleration to the cometary body

Next we connected the observed non-gravitational acceleration to the cometary activity on the surface. We assume that the nucleus is not in tumbling rotation and that the rotation period and axis orientation are fixed during a single cometary rotation. For 67P this is an excellent approximation since the orientation over $800$ days changed only by $0.5^{\circ}$ and the comet rotation period $T_{\rm rot}$ decreased by $21$ minutes from $12.4$ h $300$ days before the 2015 perihelion (see Godard et al. (2017); Kramer et al. (2019)).

The Marsden non-gravitational acceleration given by Eq. (2) is restricted to the direction dictated by the time-independent linear combination $A_{1},A_{2},A_{3}$ of the co-moving basis. This rigid link ignores the physical properties of the nucleus, in particular the rotation axis orientation and rotation period encoded in the angular velocity vector $\@vec{\omega}$. For cometary activity driven by the solar illumination on the nucleus the $A_{1}$, $A_{2}$, $A_{3}$ components are no longer time independent. Sekanina (1967) studied the arising temporal variation of the Marsden parameters under the assumption of a fixed orientation of the rotation axis. Sekanina introduced a coordinate system that takes into account the obliquity of the comet equatorial plane with the orbital plane to study the illumination conditions of the subsolar point during the orbital motion. This approach was further extended by Whipple & Sekanina (1979), Sekanina (1984) and Sitarski (1990) to time-dependent Marsden parameters, including precession models with a changing rotation axis and associated oblateness of the nucleus. Before the rotation state of 67P was known, Krolikowska (2003) applied different models to 67P including forced precession solutions.

For 67P a detailed shape model is available from Preusker et al. (2017) and the changes of the rotation state are known (see Jorda et al. (2016); Kramer et al. (2019)). 67P showed a very repetitive diurnal pattern of gas and dust in the coma across the entire illuminated nucleus, indicating a very regular and periodically repeating outgassing Kramer & Noack (2016); Kramer et al. (2017, 2018); Läuter et al. (2019). Of particular interest is the surface activity with respect to the subsolar point in the body frame. For a given rotation vector $\@vec{\omega}$ and position vector $\@vec{r}$ of the nucleus the rotation matrix $\@tens{R}_{\mathrm{com\rightarrow equ}}$ is the transformation from the cometary equatorial frame (without nucleus rotation) to the Earth’s equatorial frame,

[TABLE]

with $\@vec{h}=-(\@vec{\omega}\times\@vec{r})\times\@vec{\omega}$. It puts the sun at a fixed subsolar longitude. This construction is very similar to Sekanina’s system (both share the basis vector $\@vec{\omega}/|\@vec{\omega}|$). Any arising force (observed or modeled) from cometary activity $\@vec{F}^{\mathrm{equ}}=\@tens{R}_{\mathrm{com\rightarrow equ}}R_{\rm z}(-\omega t)\@vec{F}^{\rm bf}$ is expanded in a Fourier series with respect to the subsolar longitude

[TABLE]

The Fourier coefficients $D$, in principle defined for each single rotation, are slowly varying functions with the orbital positions around the sun. This expression encompasses comets with few active regions (see Jewitt (1997) for a simple model of a rectangular shaped comet) as well as globally active ones. The Fourier representation facilitates the rotational averaging across one rotation period ($\lambda_{\sun}(t)=-\omega t$ for 67P) in the inertial system

[TABLE]

with the rotation matrix around the $z$ axis $\@tens{R}_{\mathrm{z}}$ in the notation of Montenbruck & Gill (2000). The final expression for the non-gravitational acceleration acting on the orbit of a comet with mass $M$ ($10^{13}$ kg for 67P) is given by inserting Eq. (6) into Eq. (7)

[TABLE]

with the three linear combinations remaining from the complete Fourier expansion

[TABLE]

The parameters $D_{1}$, $D_{2}$, $D_{3}$ are the force coefficients with respect to the cometary equatorial frame (represented by the transformation in Eq. (5)). The link from Eq. (5) to the commonly used Sekanina angles is established by expressing

[TABLE]

in terms of the orbital vectors $\@vec{P}$, $\@vec{Q}$, $\@vec{R}$, the obliquity $I$ of the orbit plane to the equator of the comet, the argument of the subsolar meridian at perihelion $\Phi$, and the longitude of the subsolar meridian from the ascending node of the orbit plane on the equator $\theta_{0}$. The Marsden basis vectors in Sect. 2.1 have the representation

[TABLE]

with the true anomaly $\nu$. The vector with longitude $\lambda_{D}$ and latitude $\phi_{D}$ in the cometary equatorial frame transforms to the components $U_{\mathrm{r}},U_{\mathrm{t}},U_{\mathrm{n}}$ for the Marsden basis by the relation

[TABLE]

Setting $\theta+\pi=\lambda_{D}$ and $\phi=-\phi_{D}$ yields Eq. (2) in Sekanina (1981). Assuming additionally $-\phi_{D}$ to be the latitude of the subsolar point this relation simplifies to the classical Eq. (4) in Sekanina (1981).

4 Physical properties from non-gravitational acceleration

The parameters $D_{1}$, $D_{2}$, $D_{3}$ and with it the non-gravitational acceleration arise from the diurnally averaged activity along the spin axis and the amplitudes in the equatorial plane of the comet. This is reflected in the observed data once decomposed in the cometary equatorial frame in terms of spherical coordinates, see Fig. 3. The longitude $\lambda_{D}$ and latitude $\phi_{D}$ denote the direction of the acceleration with respect to the subsolar point. Around perihelion the observed acceleration points towards the north in accordance with the subsolar latitude $\phi_{\rm sun}$ around $-50^{\circ}$, while at the two equinox crossings the acceleration lies in the equatorial plane. Up to a shift towards more southern latitudes the seasonal variation of the subsolar latitude is reflected in the observations. This provides a confirmation of the validity of Sekanina’s approach (see Sect. 3) with respect to the seasonal component. The determination of the diurnal lag angle with respect to the solar illumination shows larger uncertainties with no lag discernible up to perihelion. After perihelion the acceleration vector lags behind (in time) with respect to the momentary anti-solar direction up to $50^{\circ}$. This lag disappears around 140 days after perihelion, when the coma gets more and more CO2 dominated (see Läuter et al. (2019)). The lag cannot be explained by a forced precession model, since the rotation state changes of 67P are small. In principle a varying surface activity across the nucleus can lead to a lag angle if local surface normal directions do not point towards the sun (see Samarasinha et al. (1996) for an illustration). Davidsson & Gutiérrez (2004) studied the variations of the non-gravitational force with respect to varying activity patterns across an ellipsoidal nucleus (see their Fig. 6). They find small directional changes of the non-gravitational acceleration. Coma observations of 67P indicate a very repetitive gas and dust release across the entire illuminated surface, see Kramer et al. (2018). This is inline with the observation that sublimation models of 67P are not capable of reproducing the reported orbit $\@vec{r}_{\mathrm{ESOC}}$ within an Earth-bound range error $<10$ km. For instance, the model of Attree et al. (2019) results in a deviation of $140$ km with respect to the observed Earth bound range $150$ days after perihelion. Besides the shape of the nucleus, another possible cause of an diurnal lag angle is a delay between the maximum illumination and the highest gas release in some areas.

From the observation of sunset jets on 67P Shi et al. (2016) obtain thermal delays of peak surface and sublimation temperatures of about $1-2$ h (corresponding to a rotation of $29^{\circ}-58^{\circ}$) in the Ma’at region in April 2015. Their calculations are based on a thermal inertia of $50$ Wm*-2K-1s1/2*. At perihelion and at other locations on the surface, the thermal delays might differ from these values, or effects can cancel. Finally, the observed lag angle could be affected by an (unknown) systematic error of the ESOC orbit, which cannot be detected by the present analysis. A conclusive attribution of the observed lag angle to a physical process requires further modelling.

The initial motivation for the study of non-gravitational acceleration was to deduce the total production flux of a comet. The total sublimation flux $Q(r)$ is approximately given by the product of the average gas velocity $v_{\rm gas,av}$ with the magnitude $a_{\mathrm{NG}}$ of the non-gravitational acceleration divided by the cometary mass $M$

[TABLE]

The Fourier decomposition of the total force allows one to clarify the relation between sublimation flux and non-gravitational acceleration after factoring out the average gas velocity $v_{\rm av}$:

[TABLE]

Next, we applied the triangle inequality to obtain a lower bound for the sublimation flux

[TABLE]

Inserting Eq. (6) for the force components and neglecting faster oscillating terms allowed us to perform the integration analytically

[TABLE]

where $\Delta=D_{x0}^{2}+D_{y0}^{2}+D_{z0}^{2}+D_{x2}^{2}+D_{y2}^{2}+D_{z2}^{2}$, ${\rm E}(x)$ denotes the complete elliptic integral. The result shows that the magnitude of the non-gravitational acceleration contains additional force components besides the $D_{1}$, $D_{2}$, $D_{3}$ components, which lead to minor corrections.

An entirely independent determination of the gas production for 67P has been performed by Läuter et al. (2019) by analyzing the in-situ data of the Double Focusing Mass Spectrometer (DFMS) of the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA). The non-gravitational acceleration by Eqs. (13) shows a remarkable agreement to the ROSINA derived water production $Q_{\rm H2O}(r)$ of 67P from Läuter et al. (2019), Fig. 3, with gas velocity $v_{\rm gas,av}=480$ m/s.

5 Conclusions

The Rosetta mission to comet 67P provided the unique opportunity to retrieve the non-gravitational acceleration of a comet independently from Earth bound observations. In conjunction with the known rotation state the position data derived from Rosetta allowed us to verify commonly invoked assumptions about the non-gravitational acceleration. We have shown the sensitivity of the orbit reconstruction to the initial conditions and identified a phase-space volume 350 days before perihelion that leads to orbital solutions following the reported ESOC positions within a mean deviation of $22$ km. This close match allowed us to extract the three-dimensional non-gravitational acceleration and to relate it to the activity on the nucleus.

Using a Fourier series we have decomposed the non-gravitational acceleration into the averaged outgassing along the rotation axis and the amplitudes of the outgassing along the equatorial plane of the comet. A similar analysis and Fourier theory has been carried out by Kramer et al. (2019) for the rotation state of 67P. We provided error bounds for all derived quantities based on an extensive analysis of initial conditions for the orbital integration. We find up to perihelion no clear signal of a lag angle between illumination and force direction, while at later times deviations from the instantaneous illumination become apparent. The seasonal effect of the solar illumination on the non-gravitational acceleration is reflected by a strong correlation of the subsolar latitude and the sublimation force direction. The agreement of the observed non-gravitational acceleration with the seasonal illumination conditions demonstrates that the non-gravitational acceleration can be explained by the water-ice distribution on the entire nucleus. The derivation of the diurnal lag carries considerably larger errors at times later than 100 days after perihelion, but points to a shift toward larger delay times 50-100 days after perihelion. The non-gravitational acceleration alone is most indicative about active areas in terms of a zonal mean.

Due to the lack of longitudinal variations the detection of local active areas across the surface can not be expected. This is in contrast to the torque affecting the rotation axis orientation and rotation period: both react sensitively to local activity variations (Kramer et al. (2019)).

Finally, we have shown that solely based on the analysis of the Rosetta orbit a close agreement with the in-situ measurements of the gas coma and production of 67P exists. This connection allows one to relate Earth bound astrometry and production estimates with accurate in-situ measurements and models of cometary activity, a prerequisite to advance non-gravitational acceleration models.

Acknowledgements.

The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing computing time on the Cray XC40. We thank the anonymous referee for helpful comments.