Free Access
Issue
A&A
Volume 512, March-April 2010
Article Number A69
Number of page(s) 17
Section Astronomical instrumentation
DOI https://doi.org/10.1051/0004-6361/200912107
Published online 02 April 2010
A&A 512, A69 (2010)

An analysis of outgassing pressure forces on the Rosetta orbiter using realistic 3D+t coma simulations

E. Mysen1 - A. V. Rodionov2 - J.-F. Crifo3

1 - Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway
2 - Central Research Institute of Machine Building, Pionyerskaya St., 4, Korolev, Moscow region 141070, Russia
3 - CNRS, LATMOS, BP 3, 91371 Verrières-le-Buisson Cedex, France

Received 18 March 2009 / Accepted 7 December 2009

Abstract
A model for the interaction between a multicomponent Maxwellian atmosphere and a spacecraft is described. Multidimensional, time-dependent gasdynamical simulations of the gas coma around the recently reconstructed aspherical rotating nucleus of comet 67P/C-G is used to analyze the outgassing pressure forces on the ESA spacecraft Rosetta. The forces were in general found to be directed significantly away from the cometocentric position vector of the spacecraft. It was also found that in a maximum outgassing scenario at comet rendezvous, the outgassing pressure force exceeds the gravitational attraction from the nucleus in the cometocentric direction of the Sun. Furthermore, the highly non-spherical pressure field was found to undergo very large changes as the nucleus rotated. Still, it was possible to represent the mean pressure field experienced by Rosetta by a fairly simple model, which can be used for the determination of the comet mass and the so-called oblateness coefficient c20 from the spacecraft Doppler signal. The oblateness coefficient represents a type of asphericity of the gravity field. The determination of the so-called triaxiality coefficient of the gravity field c22 may require using the true pressure field instead of the mean pressure field.

Key words: celestial mechanics - hydrodynamics - space vehicles - comets: individual: 67P/Churyumov-Gerasimenko

1 Introduction

In 2014, the ESA probe Rosetta will reach the short-period comet 67P/Churyumov-Gerasimenko (hereafter called 67P). The spacecraft will approach the comet in the comet's own orbit, and the corresponding low relative velocity will make it possible to insert the probe into a bound orbit around its nucleus.

At the time of the approval of the mission, it was believed that by placing the approach at a large heliocentric distance ( $r_h\simeq{}3$ AU) one would deal with an inactive nucleus. Now it is widely believed though that the formation of cometary comae (the nucleus ``activity'') is due partly to solar-driven sublimation of water ice at the nucleus surface, and partly to the diffusion of more volatile molecules (for instance carbon monoxide) from subsurface regions. Observations and physical simulations show that the former process is limited to $r_h \leq 2.5$ AU, whereas the latter may extend far out into the solar system (see the review of Prialnik et al. 2004). Instrumental limitations generally preclude observations of the distant activity (exceptions are e.g. comets Halley and Hale-Bopp, found active up to 10 AU and 30 AU, respectively), hence it can often happen that a nucleus which has substantial distant activity is improperly designated as ``inactive''. Similar limitations also often forbid the detection of molecules other than H2O even at small rh in small comets. This is the case for 67P, of which only the water activity has been observed - for $r_h \leq 2.5$ AU. If all comets are similar in origin, however, volatile molecule activity can also be expected at all distances (see for instance Capria et al. 2004). Indeed, recent observations have revealed large dust grain production at $\simeq $3 AU from 67P (Marco Fulle and Jessica Agarwal, private communication, 2009) thus, in view of the very large mass-to-area ratio of about 15 kg/m2 which is reached in a worst-case scenario, it is possible that the Rosetta orbiter will experience significant outgassing pressure forces at rendezvous.

A stated goal of the Rosetta Radio Science Investigations is the determination of the nucleus' gravity field from Doppler data (Pätzold et al. 2007). In other words, the gravity field of 67P can be determined by the study of how the comet's gravity alters the line-of-sight velocity of Rosetta in its orbit around the comet. The cometocentric gravitational acceleration of an orbiting spacecraft can be derived from the gravitational potential V, which for our purposes can be parameterized by the nucleus mass $m_{\rm c}$ (treated in Mysen & Aksnes 2008) and the gravity coefficients c20 and c22, to be defined later. In order to obtain these parameters, it is essential to have a physical model that connects them to the observable data. In the case of the Rosetta orbiter, this physical model must include the outgassing pressure force from the rapidly expanding coma or atmosphere around the comet's nucleus. Cometary dust poses a hazard to the spacecraft, but does not contribute significantly to the dynamics due to the low relative velocity between the dust particles and the spacecraft (see Mysen & Aksnes 2006, and references therein).

Below (1) we base our analysis on state-of-the-art computations of the near-nucleus coma for the first time and (2) we use parameters thoroughly compatible with the existing observations of 67P. On the other hand, to make the use of the computed pressure field easier (as far as possible) we investigate whether a simple deterministic model of acceptable precision (of the kind proposed by Scheeres et al. 2000) can exist for the pressure field.

1.1 Past celestial dynamics pressure force models

Several force models for the outgassing pressure effects on a cometary orbiter have been proposed previously. Models of this kind necessarily include (a) a gas outflow model to compute the gas parameters at any point from specific assumptions on the nucleus size, shape, rotation and gas production, and (b) algorithms to assess the force magnitude and direction from the local coma gas parameters. An analysis of the existing publications shows that in all cases these two components suffer from deficiencies. It is the purpose of the present paper to propose improved approaches for the 67P case.

As a first example, in Scheeres et al. (2000) an expansion in orthogonal spherical harmonic functions is fitted to the simulated outgassing pressure field $P(r,\theta,\phi)$ of a spherical nucleus, probably derived from some kind of fluid dynamics computation (unfortunately, the details are only referenced in unpublished work). As with the force $\vec{F}$, it is taken equal (after misprint correction) to $P(r,\theta,\phi) A_{S/C} \vec{r}/r$, where AS/C is the probe cross-section and $\vec{r}/r$ the cometocentric unit radial vector. That is, the pressure force is assumed to act mainly away from the nucleus, and the force's temporal variations as the nucleus rotates are not accounted for.

In the second example (Byram et al. 2007) the gas coma is modelled as a series of twenty discrete narrow conical jets attached to a rotating ellipsoidal nucleus, all having the same radial flow velocity V, but differing openings and gas mass fluxes $\dot{M}_i(z_\odot^i)$ (kg/s) assumed to vary with the (time-dependent) solar zenith angle at their base $z_\odot^i$. In support of this assumption it is stated that such a coma structure has been ``identified'' in the 2004 flyby images of the comet P/Wild 2. Let us observe here that these are images of the dust (not gas) coma, and that no supporting dusty gasdynamical computation is mentioned in support of this ``identification''. The aerodynamic force is taken to be

\begin{displaymath}\vec{F}(t) = A_{S/C} V \sum_i{ \chi_i(t) \dot{M}_i(z_\odot^i(t)) \left(\frac{r_i^N}{r_i}\right)^2 {\mathbf e_i(t)} },
\end{displaymath} (1)

where ri and riN are the distance to the apex of the jet i of the spacecraft and of the nucleus surface respectively, ei(t) is the (time-dependent) unit vector from this apex to the probe, and $\chi_i =1$ if the probe lies within the jet iand otherwise $\chi_i =0$ (this depends upon the orientation at time t and opening of the jet i).

1.2 Lessons from past gasdynamical model results

More than a decade of 2D, 3D and 3D+t gasdynamical simulations of comae by various authors using various computational techniques (see the review of Crifo 2004a) have clearly revealed how gas outflowing in a vacuum from a geometrically irregular, rotating solid builds up a huge atmosphere around it. There are four essential characteristics.

In the first place, the surface topography of the object plays a strong role (as common sense evidences) upon the outflow, attaching specific flow structures to many topographic details (such as concavities). This is due in part to shadowing effects and in part to flow obstacle effects. The effect is that the rotation of the central object renders the gas coma flow pattern fully time-dependent (quasi-periodic). Furthermore, since gravitational asphericities are in general associated with such topographic details, the effect may introduce a correlation between pressure force variations and gravity variations. Let us stress that this topography control occurs if the gas is emitted from the whole object or from several discrete areas. Thus, Crifo et al. (2002) showed that comet Halley's dust coma fine structures are equally well reproduced whether one postulates a uniform or a non-uniform gas production. In Crifo (2004) the Wild 2 synthetic dust coma image is compared to a dusty gasdynamical computation of the near-nucleus coma around a random shaped nucleus, suggesting immediately that the Wild 2 claimed (without any supporting modelling) ``active areas'' are nothing but the signatures of Wild 2 non-smooth surface. This will be checked in the future. By smoothing their model nuclei, both Sheeres et al. (2000) and Byram et al. (2007) eliminate this essential topography control.

Secondly, it is a fundamental property of inhomogeneous vacuum outflows that they are structured by a pattern of shocks - the effect is well documented e.g. in multiple-thrusters rocket motors. This effect (also called jet-jet interaction) renders the relation between source pattern and outflow structure non-intuitive. Thus, if the discrete gas sources of Byram et al. (2007) really existed, the jet-jet interactions would create a very disordered gas coma instead of their physically impossible independent rigid gas cones. In particular, even if all jets had initially the same velocity, a considerable velocity dispersion would appear (subsonic and supersonic areas). Furthermore, as the jet fluxes vary with nucleus orientation, the resulting ``very disordered'' structure and the velocity dispersion will change dramatically. Thus, while the Byram et al. (2007) coma changes during the rotation (as requested), the changes are over-simplified compared to those of a real coma. At this point it is instructive to note that these authors cite a paper by Crifo et al. (1994) in support of their approach, which precisely evidences the jet-jet interaction, whereby the impossibility of their model is demonstrated.

Third, because the illumination of the surface varies during its rotation, so does its temperature field, hence, even if there existed a permanent corotating surface gas flux pattern linked to the topography, still the initial gas temperature and velocity would change considerably with the solar zenith angle: a constant velocity in all jets during the rotation is a physical impossibility.

Finally, in view of the enormous difference in inertia between molecules and grains, coma dust patterns are non-trivial tracers of the surface gas emission. Hence, claims to visually infer gas sources from dust coma images without any supporting dusty gasdynamic computations must be firmly discarded. This is true for comet Wild 2 beyond any possible doubt, as stated above. To be sure, an inhomogeneity of surface gas production is possible in principle, but this can only be confirmed by coupled dust/gas dynamical simulations that reproduce the images.

Below we do not assume a strong inhomogeneity of the nucleus for two reasons. First, as discussed above, most advocations in the literature of strongly localized gas emission due to a strong nucleus inhomogeneity are merely opinions, not conclusions from physical model fits of the observations. Let us observe here that the flyby images of comet Tempel 1 evidence different areas of production for different molecules (Feaga et al. 2007); still the total coma density does not exhibit dramatic changes with position on the day side. Secondly, a significant inhomogeneity of gas production does appear even when assuming a homogeneous nucleus, due to the effect of solar illumination. Thus our computed coma is representative of any possible gas coma. However, it is true that an extreme instance where the total gas production would be localized in a few separated areas might lead to locally stronger pressures than computed here.

1.3 The problem of the force model

It is well known that a solid placed in a gas flow undergoes a force with components parallel (drag) and transverse (lift) to the flow. As described above, earlier works took the latter to be negligible on an intuitive basis, but a quantitative proof of the validity of this assumption is still missing. Another problem, probably specific to the large size of the Rosetta orbiter, is that this size is not necessarily smaller than the gas mean free path. For 67P at perihelion, simple estimates show that this is certainly not the case close to the nucleus. As the gas-spacecraft relative velocity is supersonic, it follows that in such a case, the gas parameters at the spacecraft differ from those of the unperturbed flow (an attached or a detached shock wave may appear). Below we mostly address the lift vs. drag problem for which we rederive the pressure field equations from first principles. As with the question of whether or not the gas mean free path is much larger than the size of the spacecraft, we here only delineate that region of the coma (at 3 AU) where this is expected to be the case.

2 Interaction between coma and spacecraft

Let us consider a pure gas (i.e. containing identical molecules), and let us in addition assume that the gas density is large enough for the molecule velocity distribution to be well represented by a drifting Maxwellian, but small enough that the mean free path of the molecules is significantly larger than the spacecraft, an assumption which will be motivated and tested later. Adopting a cometocentric gas outflow velocity V on the order of 0.65  ${\rm km}~{\rm s}^{-1}$ (see Fig. 9) and a mass of the target comet of $m_{\rm c}=1\times{}10^{13}~{\rm kg}$ (Davidsson & Gutierrez 2005), the cometocentric velocity of the Rosetta orbiter at a close cometocentric distance $r=10~r_{\rm c}$ ( $r_{\rm c}\equiv{}2~{\rm km}$ being the mean radius of the nucleus according to Lamy et al. 2007) is $0.2~{\rm m}~{\rm s}^{-1}$. Therefore the distribution of relative velocity $\vec{U}$ between the orbiter and the gas is determined by $\vec{V}$, and the Maxwellian velocity distribution (Valorge 1995)

\begin{displaymath}f_M=\frac{\varrho}{m}~(2\pi~R~T)^{-3/2}~{\rm exp}\left(-\frac{u_\bot{}^2+u_\parallel^2+u_z^2}{2R~T}\right)
\end{displaymath} (2)

by setting

\begin{displaymath}
u_\bot{}=U_\bot{}+V\cos\vartheta,\quad{}u_\parallel=U_\parallel+V\sin\vartheta,\quad{}u_z=U_z,
\end{displaymath} (3)

where the subscripts $\parallel, \bot{}, z$ refer to the mutually perpendicular directions illustrated by Fig. 1 ($U_\bot{}$ is normal to the surface, $U_\parallel$ parallel to it and Uz completes the orthogonal right-handed system), $\varrho$ is the mass density of the coma, $\vartheta \in (~-\pi,~\pi~)$ the angle of attack (Fig. 1), T the coma temperature and $R\equiv{}k_B/m$ where kB is the Boltzmann's constant and m the mass of the gas molecule. Here we must stress that the molecular flow along the unit vector $\vec{u}_\parallel$ changes as $\vartheta\rightarrow{}-\vartheta$, leading to two distinctively different forces along the plate as the sign of the attack angle is altered. That is, the attack angle can take on both negative and positive values, corresponding to different physical situations. Note that the attack angle is counted positive from $\vec{n}$ towards $\vec{u}_\parallel$. The vector $\vec{u}_\parallel$ is parallel to the plate, lies in the plane spanned by $\vec{n}$ and $\vec{V}$, and must be defined in a way that Eq. (3) is valid. If we restrict our attention to two dimensions as in Fig. 1, we are free to define the direction $\vec{u}_\parallel$ as we like as long as it is parallel to the plate. Still, it should be noted that the unit vector's definition relative to the plate and the definition of attack angle remain fixed once chosen if Eq. (3) is to be valid. This means that $u_\parallel{}$ does not always point in the same direction, but sometimes in the opposite direction as compared to the projection of $\vec{V}$ along the plate.

To compute the pressure force on a spacecraft, it is convenient to represent it by a polyhedron and to sum the pressure forces exerted on each face of the polyhedron. Here we will further assume that this polyhedron is convex. That is to say, the whole half-space in front of each face is unobstructed by other faces of the spacecraft (thus also no face can receive molecules scattered from other faces). Then the pressure force due to any face is equal to the force exerted on one side of an isolated flat plate, computed with an appropriate area and orientation. Below, we will distinguish for convenience between ``front side plates'' and ``rear side plates''. That is, for a plate where both sides are exposed to gas pressure (as is the case in a Maxwellian gas), one is free to decide which side is ``front''. The other side then automatically becomes the ``rear'' side. It is important to note that once defined, this definition of ``front'' and ``rear'' remains fixed, independent of how the spacecraft is oriented with respect to the flow of the coma.

\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f01.eps}}
\end{figure} Figure 1:

A flat plate is seen from its side. The plate moves with a velocity $\vec{V}$ relative to the atmosphere. The unit vector $\vec{n}\equiv{}\vec{u}_\bot{}$ is normal to the surface. The unit vectors $\vec{u}_D$ and $\vec{u}_L$ define the drag and lift directions respectively. As for the attack angle, it is counted positive from $\vec{n}$ towards $\vec{u}_\parallel$. The vector $\vec{u}_\parallel$ is parallel to the plate, lies in the plane spanned by $\vec{n}$ and $\vec{V}$, and must be defined in a way that Eq. (3) is valid.

Open with DEXTER

2.1 Pressure force on a front side flat plate

Let us consider the elementary area ${\rm d}A$ of the front side of a flat plate moving with a velocity $\vec{V}$ relative to the atmosphere, with $\vec{n}\equiv{}\vec{u}_\bot{}$ as the unit vector normal to ${\rm d}A$. Let $\vec{u}_\parallel$ be a unit vector parallel to the plate, lying in the plane spanned by $\vec{V}$ and $\vec{n}$. The force ${\rm d}\vec{F}$ experienced by ${\rm d}A$ is due to the pressure of the impacting molecules and to the recoil from the re-emitted particles. We do not consider here the case where part of or all the impinging molecules would be condensed on the surface (for such a case, see e.g. Crifo 1995). We will also assume that the plate is perfectly smooth, in the sense that the flux of the reemitted molecules is symmetrical with respect to the incidence plane ( $\vec{u}_\parallel, \vec{n}$). Then, ${\rm d}\vec{F}$ has the form

\begin{displaymath}
{\rm d}\vec{F}={\rm d}F_\bot{}(\vartheta,T_w)~\vec{n}+{\rm d...
...T_w)~\vec{u}_\parallel{},\quad{}\vec{n}\equiv{}\vec{u}_\bot{},
\end{displaymath} (4)

where Tw is the plate temperature, and

\begin{displaymath}{\rm d}F_\bot{}={\rm d}F^i_\bot{}+{\rm d}F^{r}_\bot{},\quad{}{\rm d}F_\parallel={\rm d}F^i_\parallel-{\rm d}F^{r}_\parallel,
\end{displaymath} (5)

where the superscripts (i) and (r) denote the contributions from the incident molecules and from the re-emitted molecules respectively.

If we calculate the pressure exerted by the incoming gas molecules following the Maxwellian velocity distribution, we obtain for the force per unit area along $\vec{n}$

\begin{displaymath}
\frac{{\rm d}F^i_\bot{}}{{\rm d}A}=-\int_{-\infty}^{\infty}{...
...lel}\int_{-\infty}^{0}{\rm d}U_\bot{}~U_\bot{}~f_M~m~U_\bot{},
\end{displaymath} (6)

yielding
                           $\displaystyle \frac{{\rm d}F^i_\bot{}}{{\rm d}A}$ = $\displaystyle -\frac{\varrho{}~V^2}{2\sqrt{\pi}~s^2}~\{~s\cos\vartheta{}~\exp{(-s^2\cos^2\vartheta)}+$  
    $\displaystyle \sqrt{\pi}~\left(1/2+s^2\cos^2\vartheta\right)~[~1+{\rm erf}(s\cos\vartheta)~] \}.$ (7)

Here

\begin{displaymath}{\rm erf}(X)=\frac{2}{\sqrt{\pi}}\int_0^{X}{\rm d}x~\exp{\left(-x^2\right)},
\end{displaymath} (8)

and the pressure force from the incoming gas molecules is dependent on the flow Mach number M through the speed ratio s

\begin{displaymath}
s=\frac{V}{\sqrt{2~R~T}} \equiv M \sqrt{\gamma/2},
\end{displaymath} (9)

where $\gamma$ is the gas specific heat ratio. Likewise, we obtain for the component of the force per unit area parallel to the surface in the $\vec{u}_\parallel$ direction (Fig. 1)

\begin{displaymath}
\frac{{\rm d}F^i_\parallel{}}{{\rm d}A}=-\int_{-\infty}^{\in...
...int_{-\infty}^{0}{\rm d}U_\bot{}~U_\bot{}~f_M~m~U_\parallel{},
\end{displaymath} (10)

resulting in
                      $\displaystyle \frac{{\rm d}F^i_\parallel{}}{{\rm d}A}$ = $\displaystyle -\frac{\varrho{}~V^2~\sin\vartheta}{2\sqrt{\pi}~s}~\{~\exp{(-s^2\cos^2\vartheta)}$  
    $\displaystyle +\sqrt{\pi}~s~\cos\vartheta{}~[~1+{\rm erf}(s\cos\vartheta)~] \}.$ (11)

There is some uncertainty (Harrison & Swinerd 1996) related to what happens to the gas molecules after they strike a spacecraft surface, but in any case, mass is conserved. Therefore, assuming that there is no condensation of molecules at the surface, the expression for the number of gas molecules striking the surface per unit time will be useful:

\begin{displaymath}
\frac{{\rm d}N^i}{{\rm d}A}=-\int_{-\infty}^{\infty}{\rm d}U...
...d}U_{\parallel}\int_{-\infty}^{0}{\rm d}U_\bot{}~U_\bot{}~f_M,
\end{displaymath} (12)

which gives
                        $\displaystyle \frac{{\rm d}N^i}{{\rm d}A}$ = $\displaystyle \frac{\varrho{}}{m}\frac{V}{s}\frac{1}{2\sqrt{\pi}}~\{~\exp{(-s^2\cos^2\vartheta)}$  
    $\displaystyle +\sqrt{\pi}~s\cos\vartheta{}~[~1+{\rm erf}(s\cos\vartheta)~]~\}.$ (13)

The contribution from the molecules leaving the surface, ${\rm d}F^{r}$, is partly given by the momentum transfer coefficient normal to the surface

\begin{displaymath}\sigma_{\bot}\equiv{}\frac{{\rm d}F^i_\bot{}-{\rm d}F^{r}_\bot{}}{{\rm d}F^i_\bot{}-{\rm d}F^w_\bot{}},
\end{displaymath} (14)

in which ${\rm d}F^w_\bot{}$ is the force exerted by the molecules if they are completely thermalized by the spacecraft plate

\begin{displaymath}{\rm d}F^w_\bot{}=-{\rm d}N^i~P^w_\bot{},
\end{displaymath} (15)

where

\begin{displaymath}P^w_\bot{}=m~\sqrt{\frac{\pi}{2}}~\sqrt{R~T_w}=m~\frac{\sqrt{\pi}}{2}~\frac{V}{s}~\sqrt{\frac{T_w}{T}}
\end{displaymath} (16)

is the thermal momentum of a gas molecule leaving the plate with temperature Tw. The definition

\begin{displaymath}\sigma_{\parallel}\equiv{}\frac{{\rm d}F^i_\parallel-{\rm d}F...
...rallel},\quad{}{\rm d}F^w_\parallel={\rm d}N^i~P^w_\parallel=0
\end{displaymath} (17)

completes the parametrization where the contribution from both the incoming molecules as well as from those leaving the surface have been accounted for. Note that $P^w_\parallel=0$ since a molecule thermalized by a perfectly smooth plate does not have a momentum preference in this direction. For $\sigma_\parallel=\sigma_\bot=1$, the diffuse reflection approximation (Valorge 1995) is obtained.

Since the two directions of greatest interest here for the pressure force really are $\vec{n}$ and $\vec{u}_V\equiv{}\vec{V}/V$, we will use the force per unit area representation

\begin{displaymath}
\frac{{\rm d}\vec{F}}{{\rm d}A}=\varrho{}~V^2~\left(~C_{V}~\vec{u}_{V}+C_{N}~\vec{n}~\right),
\end{displaymath} (18)

derivable from Eq. (4) with the use of the relation

\begin{displaymath}
\vec{u}_V=\vec{n}~\cos\vartheta{}+\vec{u}_\parallel{}~\sin\vartheta{}.
\end{displaymath} (19)

As the pressure force coefficients will obey

\begin{displaymath}
C_{V,N}=C_{V,N}(\cos\vartheta{}),
\end{displaymath} (20)

the attack angle need only be defined in the range $\vartheta\equiv{}(~0,~\pi~)$, and we can make the simple substitution $\cos\vartheta{}=\vec{u}_V\cdot{}\vec{n}$ in the equations which follow. Accordingly,
                          CV = $\displaystyle -\frac{\sigma_\parallel{}}{2\sqrt{\pi}~s}~\bigg\{~\exp{\left(-s^2\cos^2\vartheta\right)}$  
    $\displaystyle +\sqrt{\pi}~s\cos\vartheta{}~\left[~1+{\rm erf}(s\cos\vartheta)\right] \bigg\}.$ (21)

For the force component normal to the flat plate we get

\begin{displaymath}
C_N=C_{\Delta{}}+C_{\parallel}+C_{\bot},
\end{displaymath} (22)

where
                               $\displaystyle C_{\Delta{}}$ = $\displaystyle -\frac{1}{\sqrt{\pi}~s^2}~\{~s\cos\vartheta{}\exp{\left(-s^2\cos^2\vartheta\right)}$  
    $\displaystyle +\sqrt{\pi}~(1/2+s^2\cos^2\vartheta)~[~1+{\rm erf}(s\cos\vartheta) ]~\}$ (23)

with

\begin{displaymath}
C_{\parallel}=-C_V~\cos\vartheta,
\end{displaymath} (24)

and

\begin{displaymath}
C_{\bot}=-\frac{\sigma_\bot{}}{2}~C_\Delta{}+\frac{\sqrt{\pi...
...\bot}{\sigma_\parallel}~
\frac{1}{s} \sqrt{\frac{T_w}{T}}~C_V.
\end{displaymath} (25)

These expressions correspond to those of Schaaf & Chambré (1958) as given in Moore & Sowter (1991) if projected along the $\vec{n}$ and $\vec{u}_\parallel$ directions.

At this point we may note that the mass dependence of the force due to the inflowing molecules is contained in $\varrho$ and s, which means that if the incoming gas is a mixture of molecules, it will be in general sensitive to changes of the molecular composition with position in the coma.

A significant simplification occurs in the so-called hypersonic limit $ M\rightarrow{}\infty{}$. In this case, $s\rightarrow {}\infty $, whereby it can be shown that Eqs. (21) and (23) are reduced to

\begin{displaymath}
C_V^\infty{}=-\Theta{}(\cos\vartheta)\cos\vartheta
\end{displaymath} (26)

and

\begin{displaymath}
C_\Delta^\infty{}=-2~\Theta{}(\cos\vartheta)\cos^2\vartheta,
\end{displaymath} (27)

where $\Theta{}(X)$ is the step function which is zero for negative and one for positive arguments, X. Included as Fig. 2 is a plot of CV of Eq. (21) for some values of s (dashed) as well as its hypersonic representation Eq. (26) (solid).
\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f02.eps}}
\end{figure} Figure 2:

The pressure force coefficient CV of Eq. (21) is given for different values of s, dashed lines. From bottom to top we have s=1, s=3, and $s\rightarrow {}\infty $ (solid), Eq. (26).

Open with DEXTER
In Fig. 3, realizations of the force coefficient $C_\Delta $, which together with CV determines CN, are given in the same way. Note that the functions of Figs. 2 and 3 are symmetric with respect to $\vartheta=0$, and that we therefore may operate with a positive $\vartheta $ only, notably if we use the formulation Eq. (18).
\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f03.eps}}
\end{figure} Figure 3:

From bottom to top, $C_\Delta $ of Eq. (23) is given for s=1, s=3, and $s\rightarrow {}\infty $ (solid), Eq. (27).

Open with DEXTER
For the terms invoving exponentials and/or erf functions, the so-called hypersonic approximation $s\rightarrow {}\infty $ is acceptable for $s\ga{}3$, as shown by the preceding figures. However, the case of the remaining $ C_{\bot}$ may be different since according to Eq. (25) the condition for the hypersonic approximation to hold is here $s \gg \sqrt{T_w/T}$. It is well known that hypersonic motion in a rarefied atmosphere can lead to very high surface temperatures Tw. If we assume that the design of the Rosetta spacecraft maintains its surface around 300 K, then the hypersonic condition becomes $s \gg \sqrt{3}$ for an inner coma at 100-200 K. Note that the coma temperature drops off very rapidly with increasing distance to the nucleus and can be as low as 15 K already at a cometocentric distance of 20 km. At larger distances, the coma is heated again by solar radiation.

If the hypersonic approximation is met, four important implications appear. First, the front side part of the total force, Ff, is a function of the total mass density, whereby Eq. (18) can be used with $\varrho$ as the total mass density; thus, Ff is insensitive to the molecular composition. If this is not the case, the contributions of each kind of the molecules must be summed, using its fractional mass density and speed ratio (not anymore related to M - the mixture Mach number - by the Eq. (9)). Secondly, Ff is also independent from the coma temperature, which therefore does not need to be accurately assessed: it is only necessary to make sure that it is low enough. Third, it is also independent from the spacecraft surface temperature (which is anyway probably monitored in-flight). Fourth and most importantly, it is easy to see that Ff is in this case the same as the one obtained if the molecules were all at rest with respect to one another. Thus, the velocity dispersion playing no role, any velocity distribution - not necessarily a Maxwellian - yields the same result: the gas need not be in fluid regime. We return to this fundamental result below.

2.2 Pressure force on a rear side flat plate

A rear side plate has its surface normal vector in the opposite direction to that of a ``front side'' plate parallel to it, independently from the orientation of the spacecraft with respect to the coma flow.

As our previously derived formulas for ${\rm d}F(\vartheta,T_w)$ are applicable for all attack angles $\vartheta $, the force ${\rm d}\vec{F}^b$ due to the rear side is equal in magnitude and opposite in direction to that computed above for an attack angle $\vartheta_b = \vartheta_f + \pi$:

\begin{displaymath}
{\rm d}\vec{F}^b= - {\rm d}F_\bot{}(\vartheta_b,T_b)~\vec{n}...
...{\rm d}F_\parallel(\vartheta_b,T_b)~(\vec{u}_{\parallel{}})_f,
\end{displaymath} (28)

where $\vec{n}_f$ and $(\vec{u}_{\parallel})_f$ are the unit vectors orthogonal and parallel to the front surface (Fig. 1), Tb is the back side temperature, $\vartheta_f$ is the attack angle of the arbitrarily chosen front side, and the functions $ {\rm d}F_\bot{}$and ${\rm d}F_\parallel$ are the same as in Eq. (4).

It is clear that in the hypersonic limit the simplifications found for the front side of the plate, ${\rm d}\vec{F}^f$, also apply to ${\rm d}\vec{F}^b$.

2.3 Total force from both sides of a flat plate

The ``front side'' and ``rear side'' elements of a spacecraft do not necessarily have to be equivalent in temperature, nor does there have to be a corresponding rear plate parallel or of the same size to every front plate. For the present application however, our spacecraft shape model is constructed so it can be subdivided into pairs of front side and back side plates parallel to one another and with the same area. In this case, it is convenient to treat it as a collection of isolated plates (with both sides exposed to the gas). The total force acting on an isolated flat plate, $\sum{}{\rm d}\vec{F}$, is the sum of ${\rm d}\vec{F}^f$and ${\rm d}\vec{F}^b$ given by Eqs. (4) and (28) respectively:

                             $\displaystyle \sum{}{\rm d}\vec{F}$ = $\displaystyle \left[~{\rm d}F_\bot{}(\vartheta_f,T_f)-{\rm d}F_\bot{}(\vartheta_f+\pi,T_b)~\right]~\vec{n}_f$  
    $\displaystyle +\left[~{\rm d}F_\parallel{}(\vartheta_f,T_f)-{\rm d}F_\parallel{}(\vartheta_f+\pi,T_b)~\right]~(\vec{u}_\parallel)_f.$ (29)

Defining the pressure field coefficients $\sum{}C_{V,N}$ through

\begin{displaymath}
\frac{\sum{}{\rm d}\vec{F}}{{\rm d}A}=\varrho{}~V^2~\left[~\...
...right)~\vec{u}_{V}+\left(\sum{}C_{N}\right)~\vec{n}_f~\right],
\end{displaymath} (30)

it can easily be reasoned from Eqs. (19) and (29) that
                   $\displaystyle \sum{}C_N$ = $\displaystyle C_N(T_w=T_f,\vartheta=\vartheta_f)$  
    $\displaystyle -C_N(T_w=T_b,\vartheta=\vartheta_f+\pi),$ (31)

and
                   $\displaystyle \sum{}C_V$ = $\displaystyle C_V(T_w=T_f,\vartheta=\vartheta_f)$  
    $\displaystyle +C_V(T_w=T_b,\vartheta=\vartheta_f+\pi).$ (32)

For s > 3 we have the ``weak'' hypersonic expressions
                         $\displaystyle \sum{}C_V^\infty{}$ = $\displaystyle -\sigma_\parallel{}~\cos\vartheta_f~[~\Theta(\cos\vartheta_f)-\Theta(-\cos\vartheta_f)~]$  
  = $\displaystyle -\sigma_\parallel{}~\vert~\cos\vartheta_f~\vert,$ (33)

which is $\pi$-periodic as one must require for physical reasons. In the same way, one obtains with the use of Eqs. (27) and (31):

\begin{displaymath}\sum{}C_\Delta^\infty{}=-2~\cos\vartheta_f~\vert~\cos\vartheta_f~\vert,
\end{displaymath} (34)

\begin{displaymath}\sum{}C_\parallel^\infty{}=\sigma_\parallel{}~\cos\vartheta_f~\vert~\cos\vartheta_f~\vert
\end{displaymath} (35)

and
                         $\displaystyle \sum{}C_\bot{}$ = $\displaystyle \sigma_\bot{}~\cos\vartheta_f~\vert~\cos\vartheta_f~\vert$  
    $\displaystyle -\sigma_\bot{}~\frac{\sqrt{\pi{}}}{2}~\frac{\sqrt{2R}}{V}~\left[~\sqrt{T_f}\cos\vartheta_f\right.$  
    $\displaystyle + \left.(\sqrt{T_b}-\sqrt{T_f})~\Theta(-\cos\vartheta_f)\cos\vartheta_f~\right],$ (36)

where the last term in Eq. (36) has been retained because of the uncertainty in the value of the ratio Tf/T.

For the present application to Rosetta, we will assume that the thermal design of the spacecraft ensures approximate surface isothermicity, whereby $\sqrt{T_f}-\sqrt{T_b}\approx{}0$, hence

\begin{displaymath}
\sum C_N = (-2+\sigma_\parallel+\sigma_\bot)\cos\vartheta_f\...
...rtheta_f\vert -\sigma_\bot\frac{\Lambda_f}{V} \cos\vartheta_f,
\end{displaymath} (37)

an expression that in accordance with the physical situation changes sign if $\vartheta_f$ is increased by $\pi$, and in which the term

\begin{displaymath}
\Lambda_f = \sqrt{\frac{\pi{}~R~T_f}{2}} \equiv \frac{1}{2s}\sqrt{\pi\frac{T_f}{T}}
\end{displaymath} (38)

has been retained (``weak hypersonicity'') as Tf/T ratio's value is uncertain.

Based on Eqs. (30), (33) and (37) we propose as a first attempt to interpret the radiometric data from Rosetta that three estimation parameters for the description of the interaction between spacecraft and coma can be

\begin{displaymath}\sigma_\parallel,\quad{}\sigma_\bot{},\quad{}\sigma_\bot{}\Lambda_f.
\end{displaymath} (39)

Hopefully, applications will show that they are fairly constant. Otherwise the parameters above can be allowed to change over time to a certain extent. In Moore & Sowter (1991), where the last term of Eq. (37) $\sim{}V^{-1}$ was defined to be negligible, it was noted for instance that $\sigma_\bot$ itself might be a function of the attack angle (see also Moe & Moe 2005):

\begin{displaymath}\sigma_\bot=\sigma^{(0)}_\bot{}-\sigma^{(1)}_\bot{}\sec\vartheta_f.
\end{displaymath} (40)

Interestingly, this dependency of the normal accommodation coefficient will reintroduce a $\sum{}C_N$ with the same dependency on the attack angle as in Eq. (37).

As with Ff and Fb the total force is independent of the coma temperature in the hypersonic limit, and hence from the flow regime of the gas. However, for Rosetta it may be that this condition is only partly met as discussed precedingly, leaving in Eq. (37) the term dependent on $R\equiv{}k_B/m$, i.e., on the mass of the gas molecule. In this case we need to know $\varrho{}V^2$ for each different chemical species to calculate this component of the force. In navigation and radio science, where scientific quantities are found by inverting radiometric data from the spacecraft (Pätzold et al. 2007), only the properties of the chemical species which contributes the most can be found. That is, if we can find a proper parametrization of $\varrho{}V^2$, we do not necessarily have to operate with one such set of parameters for each chemical species. Of course, this simplification does not hold if the probe enters a region during the tracking interval where the dominant contributor to the pressure shifts from one species to another.

In this paper, where the molecules CO and ${\rm H}_2{\rm O}$ are considered and the CO pressure usually dominates, we will make the approximation

\begin{displaymath}\Lambda_f({\rm H}_2{\rm O})= \sqrt{28/18} ~\Lambda_f({\rm CO})\approx{}\Lambda_f({\rm CO}).
\end{displaymath} (41)

As a result, we only need a proper description of the sum $(\varrho_{\rm CO}+\varrho_{{\rm H}_2{\rm O}})V^2$, and not descriptions for each species (given that V is common for all kinds of molecules).

Although our coefficients in Eq. (30) are to a great extent determined by the defined parametrization, we will nevertheless conclude that there is little a priori evidence for the hypothesis of a pressure force on a cometary orbiter that acts mainly in the direction away from the comet nucleus in general, as assumed in Weeks (1995), Scheeres et al. (2000) and Byram et al. (2007). That is, $\sum{}C_N$ is zero for all attack angles $\vartheta_f$ only if there is full accommodation in all directions and the spacecraft temperature is zero (unless, of course, the possibilities $\sigma_{\parallel,\bot{}}>1$ also are considered). This represents the most important property of the presented simple model: it does not exclude forms of interaction between the gas molecules and the spacecraft.

2.4 Rosetta's effective pressure force area

Let us define $\vartheta_f$ as the angle between the relative velocity vector $\vec{u}_V$ and the unit vectors normal to the solar cell arrays with a total area $64~{\rm m}^2$. Also, let one surface of the bus with an area of $2.8\times{}2.0~{\rm m}^2$ point in the same direction as the solar cell arrays. Finally, assume that the bus area in the $\vec{u}_V$ direction is $2.1\times{}2.0~{\rm m}^2$ when $\vartheta_f=\pi/2$, i.e., when the solar cell arrays are not exposed to the pressure field. With the use of Eqs. (33) and (37) we are now in a position to plot the effective outgassing pressure force area of Rosetta in different directions. The directions used will be the drag direction $\vec{u}_D\equiv{}\vec{u}_V$ and the orthogonal lift direction $\vec{u}_L$ (Fig. 1), related to a vector $\vec{n}$ normal to a spacecraft surface through

\begin{displaymath}\vec{n}=\vec{u}_V~\cos\vartheta+\vec{u}_L~\sin\vartheta,\quad{}\vartheta\in{}(-\pi,\pi).
\end{displaymath} (42)

The contribution to the drag force from both the back side and front side of a spacecraft area ${\rm d}A_f$; with a front side surface normal $\vec{n}_f$ this is then

\begin{displaymath}
\frac{\sum{}{\rm d}\vec{F}_D}{{\rm d}A_f} = \varrho{}V^2 \ve...
...(\vartheta_f)+~\cos\vartheta_f~\sum{}C_N(\vartheta_f)~\right].
\end{displaymath} (43)

Likewise, for the lift component

\begin{displaymath}
\frac{\sum{}{\rm d}\vec{F}_L}{{\rm d}A_f} = \varrho{}V^2 \vec{u}_L~\left[~\sin\vartheta_f~\sum{}C_N(\vartheta_f)~\right].
\end{displaymath} (44)

Included as Fig. 4 is the effective area of Rosetta in the drag direction, based on Eq. (43). A typical value of $V=500~{\rm m}~{\rm s}^{-1}$ has been used, together with a spacecraft temperature of $T_f=200~{\rm K}$, assuming full accommodation normal to the plates, $\sigma_\bot=1$, and no accommodation of the gas molecules along the plates, $\sigma_\parallel=0$. The low mass of water molecules is chosen for the calculations to maximize the effect.
\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f04.eps}}
\end{figure} Figure 4:

The figure illustrates effective areas of Rosetta in the drag direction: total (solid), without bus (dotted), or with the spacecraft temperature set to $T_f=0~{\rm K}$ (densely dotted).

Open with DEXTER
Clearly, the effective area, amplified by the thermalization of gas molecules by the spacecraft walls, can be very large. The lift to drag ratio is plotted in Fig. 5. That is, one cannot argue that the force is directed mainly away from the comet nucleus in general based on the claim that the lift component is negligible compared to the drag component.
\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f05.eps}}
\end{figure} Figure 5:

The lift to drag ratio is shown above as a function of the attack angle $\vartheta $.

Open with DEXTER

3 Comet 67P model

The only reliable way to predict the behaviors of the important functions $\varrho{}V^2$ and V of the previous section is to base the predictions on state-of-the art computations of the 3D+t gas coma, making sure that all adopted parameters are fully compatible with the present information available on the 67P nucleus, summarized in Lamy et al. (2007). To be sure, considerable uncertainty persists on the properties of this nucleus. In particular, the CO production rate cannot be measured because its upper-limit lies below the detection limit of existing ground-based instruments, and the nucleus shape and perhaps rotation parameter are provisory. One can see the evolution in the best-fit nucleus parameters by comparing the present shape and coma properties (Figs. 7-9) with those given in Crifo et al. (2004b), which were based on 2002 data (let us remind that 67P has a period of about 6.5 years, hence new data are acquired every six to seven years). While there is an unavoidable uncertainty in the present parameters we can warrant that our assumptions are compatible with all existing observations by basing our parameter choice on Lamy et al. (2007), as opposed to assumptions which are in direct conflict with the observations, as those of a spherical or ellipsoidal nucleus shape. By using a 3D + t gasdynamical approach, we can guarantee the robustness of the physical model, as opposed to using algorithms which rely on untenable concepts, as a corotating rigid gas coma pattern.

The main model parameters are summarized in Table 1. We call them ``maximum model'' because they correspond to the upper limits of the total water and CO production rates. Evidently, the real conditions could as well be those of minimum production rates, which are an order of magnitude smaller - our scenario is in this sense a worst-case scenario from the point of view of pressure effect on the orbiter. However, as discussed below, the surface distribution of these fluxes also influences the pressure and at this point in time it is not known whether our present assumptions are worst-case assumptions or not.

Table 1:   Comet 67P model parameters.

3.1 Nucleus model

Using HST images of the comet, Lamy et al. (2007) have been able to extract an approximate size, shape, and rotation mode of this nucleus. Figure 6 shows that, as all imaged cometary nuclei, the nucleus of 67P appears to be small (km size) and very irregular in shape. We have adopted the shape of Lamy et al. (2007), and - consistently - their rotation model. The nucleus radii are distributed around $r_{\rm c}\equiv 2 $ km. The rotation is uniaxial and periodic with a period P=12.55 h. The constant angle between the cometocentric direction of the Sun and the axis of rotation is  $70^\circ{}$. Uniaxial rotation has been inferred from the presently existing lightcurves by Lamy et al. (2007), but a more realistic mode should not be entirely excluded.

\begin{figure}
\includegraphics*[scale=0.25,angle=0]{12107f06.eps}
\end{figure} Figure 6:

Present best-fit shape of the comet 67P nucleus. The largest diameter is about 5 km. The uncertainty is probably on the order of the magnitude of the topographic features. The corresponding best-fit monoaxial rotation is about the z-axis. The model will be periodically improved as new observational data are accumulated. After Lamy et al. (2007).

Open with DEXTER

The surface flux Z(CO) at any point with the solar zenith angle $z_\odot$ is given by:

\begin{displaymath}Z({\rm CO}) = Q_{{\rm CO}}\left(\frac{A_0}{{\cal A}_{\rm ext}...
...(1 -A_0) ~{\rm max}[\cos{z_\odot},0.]}{{\cal A}_\odot}\right), \end{displaymath} (45)

where ${\cal A}_{\rm ext}$ is the computed nucleus external area, ${\cal A}_\odot$ its cross-section seen from the Sun, A0 a dimensionless parameter, and $Q_{{\rm CO}}=10^{27}$ mol.s-1 is the desired total CO production rate. One sees that a fraction A0 of $Q_{{\rm CO}}$ is distributed uniformly over the surface, and a fraction (1 - A0) is distributed over the sunlit surface in proportion to its illumination. Thus, for the presently adopted value A0 = 0.2, the CO production has a pronounced day-to-night asymmetry. The water flux at any point is calculated by solving a sublimation energy budget equation at each time. In this equation, a heuristic term Fs represents the exchanges of heat between the nucleus interior and the surface by conduction. The surface is assumed to have at all points the same small amount of ice, representing an icy area fraction f = 0.07. The resulting total production rate $Q_{{\rm H_2O}}$varies during the rotation around $\sim{}6\times{}10^{26}~{\rm mol.}{\rm s}^{-1}$ at the heliocentric distance $r_h=3~{\rm AU}$.

As stated precedingly, other assumptions with respect to the surface distribution of the ice and of the CO emission could have been made, because at the present time there is no observational constraint on it. For instance, both could have been assumed localized in small areas of the nucleus, leading to a much more contrasted gas coma structure. It would certainly be interesting to investigate in the future the implications for the pressure force field.

3.2 Molecular outflow model

The vacuum outflow of the molecules emitted by the nucleus is computed using the gasdynamical approach described in Rodionov et al. (2002), except that the method of solutions here is more powerful than the one described there, as one true time-dependent solutions is acquired (in lieu of a succession of steady state solutions). Thus, the time dependent tridimensional Euler Equations governing the gas outflow are solved by the shock-capturing Godunov method, using periodic boundary conditions for the integration over time. The time mesh is 1/72 of a rotation period, hence the solution consists of a set of 72 3D files. A brief description of this new family of solutions is given in Rodionov and Crifo (2006). The adequacy of Eulerian Equations for the present parameter choice and the presently considered cometocentric distances will be discussed below.

4 Properties of the computed 3D+t gas coma

The time-dependent, periodic but non-co-rotating structure of the coma can be seen on Fig. 7, which shows snapshots of the total gas density in the rotation equatorial plane at nine equispaced times during one rotation period. One notes a clear day-to-night asymmetry, due in part to our assumed (modest) solar dependence of the CO flux, and in part to the solar-driven water production. Also evident on the figure is the correlation between density structures and nucleus topography, and the resulting time-dependence of the gas density at any point.
\begin{figure}
\par\includegraphics[scale=0.25,angle=-90]{12107f71.eps}\hspace*{...
...\hspace*{4mm}
\includegraphics[scale=0.25,angle=-90]{12107f79.eps}
\end{figure} Figure 7:

Variation of the total mass density in the equatorial plane during the rotation. From left to right and from top to bottom: isocontours of Log10 (mass density, kg/m3) in the equatorial plane are shown every 1/9 of the rotation period (i.e., every 1.4 h.) The projected Sun direction is to the right. Notice the clockwise rotation of the nucleus cross section. On this and the two next figures, references and dates appear at the bottom left of each panel for safety of identification: the reader must overlook them. The dust contribution has not been included.

Open with DEXTER

To have a better feeling about the 3D structure of the coma and about its compositional variability, Fig. 8 shows the H2O and CO density distributions in the equatorial plane and in one meridional plane at rotation phase 0. It is clear that there is a very strong day-to-night change in the chemical composition due to the different physical origin of the molecules. Notice that the two patterns depend on one another, since under the conditions of interest here, the two species are closely coupled together by collisions and form a single fluid. As a consequence, if only water or only CO was produced, their coma patterns would differ from the present ones, because the distribution of surface gas flux would be different.

\begin{figure}
\par\includegraphics[width=17cm,clip]{12107f08.eps}
\end{figure} Figure 8:

Log10 of the H2O ( upper panels) and CO ( lower panels) density in two orthogonal planes at rotation phase 0. The projected Sun direction is to the right on the X-axis ( right panels) or at 200 to the X-axis ( left panels). The total density shown on Fig. 7 ( upper left panel) is the sum of the two present right panel densities.

Open with DEXTER

Figure 9 shows the distribution of the flow velocity and the flow Mach number in the same planes. There is not much day-to-night asymmetry in the Mach number, a general property of inviscid flows, but there is a very large difference in the flow velocity, due to the strong day-to-night surface temperature asymmetry. This temperature asymmetry is controlled in our model by the parameter Fs. Our present value for this parameter is probably a minimum, yielding a probable maximum day-to-night temperature asymmetry.

\begin{figure}
\par\includegraphics[width=16.5cm,clip]{12107f09.eps}
\vspace*{-2mm}
\end{figure} Figure 9:

Velocity ( upper panels) and Mach number ( lower panels) in two orthogonal planes at rotation phase 0. The projected Sun direction is to the right on the X-axis ( right panels) or at 200 to the X-axis ( left panels).

Open with DEXTER

Finally, the direction of the gas flow is also of interest. For instance, in Crifo et al. (1995) it was shown that the gas flow can be almost along the nucleus surface close to the nucleus i.e., highly non-radial. For the present coma simulations, we find that the angle between the direction of flow and the cometocentric radius vector is on average less than $7.5^\circ{}$ at $r=5~r_{\rm c}$. At distances r=10 and $25~r_{\rm c}$ the angle is on average less than $4.5^\circ{}$ and $2.5^\circ{}$, respectively.

In conclusion, we hope that these results (in all respects similar to those of previous simulations, see Crifo et al. 2004a) will convince the reader that there is no co-rotating feature of any kind in a realistic gas coma, nor any rotationally invariant parameter, in particular not the velocity or the temperature.

4.1 Validity of the present 3D+T Eulerian solution

A gas is said to be in fluid regime, transition regime, or free-molecular regime, according to whether its velocity distribution is a Maxwellian, a first-order-distorted Maxwellian, or a non-Maxwellian function (see e.g., Rodionov et al. 2002). If the hypersonic approximation holds, this distinction is immaterial to compute the pressure force from coma gas parameters. But this approximation does not hold in the immediate vicinity of the nucleus surface (see below), in which case the presently derived force formulas are valid only in fluid regime. On the other hand, the gas flow is governed by Euler equations only in fluid regime, hence the presently computed gas parameters are, strictly speaking, valid only in this case. It can be shown that the gas fluid regime can be inferred from the comparison between the characteristic flow scale ${\cal L}$ in the solution and the computed gas collisional mean free path (m.f.p.) $\Lambda_g$. Still, the precise boundaries between the different regimes vary with the problem considered, and can only be found from specific Direct Monte-Carlo Simulations of the flow (DSMC). A program of systematic investigations of this kind for cometary comae is in progress (see Zakharov et al. 2008, and references therein). The results obtained hitherto show that Euler equations provide correct density and velocity in the present conditions at least for $\Lambda_g / {\cal L} \leq 0.1$.

The definition and computation of $\Lambda_g$ in a gas mixture are complicated and will not be dealt with here. Instead, we compute a pseudo water free path, as if all the gas were pure water, and a pseudo CO free path, as if it were pure CO. For the calculation of the m.f.p. of the molecules, we use the so-called ``variable hard sphere'' (``VHS'') model (Rodionov et al. 2002). Here another difficulty appears: the evaluation of the m.f.p. at very low temperatures. In this case, the binary collisions should be treated using quantum mechanics instead of the VHS approximation. Since these conditions are never met at the laboratory, such computations are not available. Thus, to avoid the (unphysical) divergence of the VHS approximation at very low temperatures, we introduce an ad hoc lower-limit temperature $T_{\rm min}$. To estimate the inaccuracy of this, we have played with this lower-limit. Figure 10 shows the computed m.f.p. on the subsolar axis at rotation phase $0^\circ {}$for two values of $T_{\rm min}$. We see that the maximum uncertainty is of a factor of 5.

\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f10.eps}}
\end{figure} Figure 10:

``Pseudo'' mean free paths for water (solid lines) and CO (dashed lines) are shown, expressed in multiples of the characteristic orbiter size $D_{{\rm S/C}}=32~{\rm m}$ as a function of cometocentric distance expressed in multiples of $r_{\rm c}=2$ km along the subsolar axis. The time is the one we defined as the initial rotational phase of the nucleus. Two curves are plotted for each molecule, corresponding to the two minimum ad hoc temperatures $T_{{\rm min}}=15~{\rm K}$ ( lower) and $T_{{\rm min}}=50~{\rm K}$ ( upper).

Open with DEXTER
This being said, the worst case is found picking the largest value of the m.f.p. in Fig. 10 and comparing it to the characteristic gradient scale ${\cal L}$ in the flow. Taking the latter to be r/2 as in a spherical outflow, one computes easily that the (maximum) m.f.p. is smaller than ${\cal L}/10$ in the range of Fig. 10, which ensures the validity of the Euler solution at least to 50 $r_{\rm c}$.

5 Properties of the computed pressure field

5.1 Validity of the free-molecular spacecraft-gas interaction model

An essential assumption made in the preceding evaluation of the pressure force is that the spacecraft is submitted to a molecular flux given by the unperturbed coma gas flow velocity distribution. It can be shown that this requires that $\Lambda_g$ be larger than all dimensions of the spacecraft (the so-called free-molecular condition). If this is not the case, it acts as an obstacle to the flow, resulting in changes in the flow parameters close to it.

Here, the worst case occurs with the smallest values of $\Lambda_g$. This would lead to a restriction of our treatment to cometocentric distances larger than about $10~r_{\rm c}$. That is, below this altitude and given the supersonic relative spacecraft-gas motion, the gas-spacecraft interaction occurs via the formation of a shock wave (attached or detached) in front of the spacecraft. The latter is submitted to the pressure force of the shocked gas, hence its computation first requires that of the shock structure.

5.2 Range of validity of the hypersonic approximation

For the present conditions where CO is dominant, $\gamma$ is close to the vibrationally relaxed CO value 1.4, hence s = 0.84 M from Eq. (9), and the condition $s\ga{}3$ is equivalent to $M\ga{3.6}$. Figure 9 indicates that this condition is always met beyond 10 km ( $r \geq 5 r_{\rm c}$) and never met at r < 5 km ( $r \leq 2 r_{\rm c}$). Let us note that these distances do not depend much upon the gas density, but depend strongly upon the nucleus shape and size.

5.3 Spatial structure of the pressure field

Now we return to the important prefactor $\varrho{}V^2$ of Eq. (30). First one should note that for the unphysical case of a spherical nucleus with homogeneous sublimation properties across its surface, the prefactor is not dependent on the probe's longitude $\varphi $ in the cometocentric solar plane-of-sky (plane through the comet center, normal to the direction of the Sun). If this property also holds for the adopted nucleus of irregular shape, then the contours on which the prefactor is constant should be nearly horizontal in $\varphi $-$\theta $ space, where $\theta $ is the angle between the cometocentric direction of the spacecraft and the cometocentric direction of the Sun. We also want to test if the pressure field drops as the inverse square of the distance r from some center of the comet. If so, the function

\begin{displaymath}
\chi=\xi~\varrho{}~V^2~(r/r_{\rm c})^2,\quad{}\xi\equiv{}120~{\rm m^2}/1800~{\rm kg}
\end{displaymath} (46)

should be fairly independent of r. The prefactor $\xi$ is chosen to display the maximum drag acceleration according to Fig. 4, with almost all of Rosetta's fuel spent. Note that the physical maximum exposed area of Rosetta is around 70 m2, but that the effective area with respect to outgassing pressure forces has been nearly doubled to 120 m2 due to the non-zero spacecraft temperature.

In Fig. 11 the isocontours of $\chi $ are plotted in $\varphi $-$\theta $ space for the initial rotational phase of the nucleus, where $\varphi $ is measured in the solar plane-of-sky, relative to the plane which contains both the rotation axis and the cometocentric direction vector of the Sun. The different types of curves represent different cometocentric distances of the spacecraft: r=10 (dashed), 25 (solid) and $50~r_{\rm c}$ (densely dotted).

\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f11.eps}}
\end{figure} Figure 11:

The contours in $\varphi $-$\theta $ space on which the pressure field prefactor $\chi\equiv{}\xi~\varrho{}V^2~(r/r_{\rm c})^2$ is constant are shown for a nucleus rotation phase of $0^\circ {}$. Above, $\chi $ takes on the values 0.1 (in the anti-solar direction), 0.25, 0.5, 0.75, and with the peak $1~r_{\rm c}~{\rm h}^{-2}$ in the general direction of the Sun, $\theta \sim {}0-\pi /4$. The different types of curves represent different cometocentric distances: r=10 (dashed), 25 (solid) and $50~r_{\rm c}$ (densely dotted).

Open with DEXTER

From the figure we can infer three important properties of the pressure field. First of all, near $\theta \sim {}0$ the strength of the acceleration is $\sim$ $1~(r_{\rm c}/r)^2~r_{\rm c}~{\rm h}^{-2}$, which is comparable to the gravitational acceleration generated by the comet

\begin{displaymath}\ddot{\vec{r}}_{{\rm gravity}}=1.1~(r_{\rm c}/r)^2~r_{\rm c}~{\rm h}^{-2},\quad{}{\rm spherical}~~{\rm nucleus}
\end{displaymath} (47)

with a nominal nucleus mass of $m_{\rm c}=1\times{}10^{13}~{\rm kg}$(see for instance Mysen & Aksnes (2008) and references therein). That is, precautions need to be taken (with respect to array orientation, orbit) if the orbit is not to be continuously stabilized by thrusters. Secondly, the irregularity of the nucleus' shape distorts the isocontours considerably, i.e., the contours are not horizontal as one would expect from an unrealistic nucleus shape.

Third, the contour types (solid, dashed, densely dotted) representing different cometocentric distances overlap fairly well. In other words, the outgassing pressure field drops roughly as the inverse square distance from some center of the nucleus. From the figure one can see that this property begins to degrade as r decreases, but studies show that important structures are still preserved as close to the nucleus as $r\sim{}5~r_{\rm c}$. This means that we can approximate the pressure field, or more precisely $\chi=\chi(\varphi,\vartheta)$, by Fourier expansions in $\theta $ and $\varphi $, if not too close to the comet's nucleus ( $r\sim5{-}10~r_{\rm c}$). However, Fig. 11 demonstrates that such Fourier expansions will be of a relatively high order due to the large degree of variability of $\chi $. Also, the coefficients of the expansions will clearly exhibit large temporal variations as the nucleus rotates.

5.4 Variability

This is shown explicitly in Fig. 12 where isocontours of $\Delta{}\chi\equiv{}\chi_{{\rm max}}-\chi_{{\rm min}}$ are mapped over a rotation period.

\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f12.eps}}
\end{figure} Figure 12:

The isocontours of the amplitude $\Delta {}\chi $ over a rotation period are shown. Above, $\Delta {}\chi $ takes on the values 0.05 (dashed, $\theta \sim {}3\pi /4$) 0.1, 0.25, 0.5, 0.75 and $1~r_{\rm c}~{\rm h}^{-2}$ (dashed, $\theta \sim {}\pi /4$). The probe's cometocentric distance is set to $r=10~r_{\rm c}$ for convenience.

Open with DEXTER
These temporal variations with rotation are particularly important when the second degree and order gravity field of the nucleus are to be extrapolated from Doppler data (Pätzold et al. 2007), since both effects induce variations in the Doppler data with roughly the same frequencies. Notice the small amplitudes in the night-side of the coma in Fig. 12. This means that the anti-solar direction is to be preferred if outgassing effects on the extrapolation of the higher order gravity momenta are to be minimized. Figure 13 shows the $\chi $ variations in two regions where they are largest, demonstrating that this type of variation is in general non-trivial.
\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f13.eps}}
\end{figure} Figure 13:

The variation of $\chi $ with nucleus rotation for $(\varphi ,\theta )=(\pi ,\pi /4)$ (dashed), and for $(\varphi ,\theta )=(0,\pi /4)$ (solid). The probe's cometocentric distance has for convenience been set to $r=10~r_{\rm c}$.

Open with DEXTER

However, if the outgassing pressure perturbation of the spacecraft's Keplerian motion is sufficiently weakened by a minimized exposure of the probe's arrays, it follows from first-order perturbation theory that the mean motion of the orbit, and therefore also the Doppler observable, can be characterized by a time-averaged pressure field. This way, a systematic drift of the Doppler observable related to the non-periodic evolution of the spacecraft's cometocentric orbit can be accounted for over a limited time-scale. It should be stressed though that the Doppler signal related to the short-periodic (similar to rotation period of nucleus) variations of the orbit, induced by the outgassing pressure, are not accounted for in a mean field approximation. So the variations which are not accounted for can mask the Doppler signal related to the second degree and order gravity field of the nucleus. Still, to determine the nucleus oblateness (second degree, zero order gravity field), nucleus mass and position of the spacecraft relative to the nucleus, it is a proper description of the mean pressure field which is important, but only if it can be viewed as a weak perturbation.

5.5 Mean field

We will return to these topics later on in this paper, but before doing so, some properties of the mean outgassing pressure field $\overline {\chi }$ (strictly speaking not the pressure field, but the spacecraft's response to the pressure field) based on Eq. (46) will be quantified. The isocontours of the time-averaged field with equal weight on the 72 available rotational phases of the comet nucleus are shown in Fig. 14. The different types of curves represent the cometocentric distances r=10 (dashed), 25 (solid), and $50~r_{\rm c}$ (densely dotted). From plots like Fig. 11 we already know that $\overline {\chi }$ will be fairly independent of distance r, but Fig. 14 shows that this property is slightly improved.

\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f14.eps}}
\end{figure} Figure 14:

The curves on which the time averaged pressure field $\overline {\chi }$ is constant are shown for the values 0.075 ( $\theta \sim {}\pi $), 0.1, 0.25, 0.5, 0.75, 1, and $1.2~r_{\rm c}~{\rm h}^{-2}$ ( $\theta \sim {}0$). Again, the different types of curves represent different cometocentric distances: r=10 (dashed), 25 (solid) and $50~r_{\rm c}$ (densely dotted).

Open with DEXTER

As one would perhaps expect, this mean field is more well-behaved than the instantaneous ones of Fig. 11. In order to quantify this property, we have fitted $\overline {\chi }$ to expansions in spherical harmonic functions (Kaula 2000) of the order N

                      $\displaystyle \overline{\chi}^{(N)}(r)$ = $\displaystyle \sum_{n=0}^{N}\sum_{m=0}^{n}\gamma_{nm}~P_{nm}(\cos\theta)$  
    $\displaystyle \times{}\left(~\overline{a}_{nm}\cos{m\varphi}+\overline{b}_{nm}\sin{m\varphi}~\right)$ (48)

by the least-squares method. Above, $P_{nm}(\cos\theta)$ are the associated Legendre functions, and

\begin{displaymath}\gamma_{nm}=\sqrt{(n-m)!(2n+1)(2-\delta_{0,m})/(n+m)!}
\end{displaymath} (49)

is a normalization factor so that the values of $\overline {a}_{nm}$ and $\overline{b}_{nm}$ reflect the convergence of the series. The quantity

\begin{displaymath}
\overline{\Delta{}}(r_1,r_2)\equiv{}\frac{1}{100^{2}}\sum_{i...
...erline{\chi}_{ij}(r_2)}{\overline{\chi}_{ij}(r_2)}\right\vert,
\end{displaymath} (50)

defined on a $100\times{}100$ grid in $\varphi $-$\theta $ space where $\overline{\chi}_{ij}(r)\equiv{}\overline{\chi}(r,\varphi_i,\theta_j)$, measures the overall adequacy of the fit when r=r1=r2. For $r_1\ne{}r_2$, the quantity is a measure of how well the expansion Eq. (48) fitted to $\overline{\chi}(r_1)$ also works as an approximation for $\overline{\chi}(r_2)$. The maximum deviation on the grid

\begin{displaymath}
\Delta{}(r_1,r_2)\equiv{}{\rm MAX}_{i,j}\left\vert\frac{\ove...
...verline{\chi}_{ij}(r_2)}{\overline{\chi}_{ij}(r_2)}\right\vert
\end{displaymath} (51)

will also be of interest.

Accordingly, we see that if we define $\overline{b}_{nm}\equiv{}0$ due to the approximate symmetry of Fig. 14, a truncation at N=3 of Eq. (48) represents a good balance between keeping N low on one side and $\Delta{}$ and $\overline{\Delta}$ small on the other, Table 2.

Table 2:   Fit adequacy parameter $\overline {\Delta }(r_1,r_2)$ (and $\Delta (r_1,r_2)$) for an N=3 expansion Eq. (48) fitted to $\overline{\chi}(r_1=25~r_{\rm c}$).

The ten estimated parameters of $\overline{\chi}^{(3)}\approx{}\overline{\chi}(r_1=25~r_{\rm c})$ are given in Table 3.

Table 3:   Coefficients $\overline {a}_{nm}$ of an expansion in spherical harmonics fitted to $\overline{\chi}(r=25~r_{\rm c})$ of Eq. (46).

The fit also works well for other distances, see Table 2.

Of course, whether or not the approximation $\overline{b}_{nm}\equiv{}0$ is valid for a real nucleus is highly uncertain. For instance, the outgassing across the nucleus surface may be highly dependent on the different types of nucleus ices and the location of the nucleus sublimation fronts, and the mean field may therefore have other directions of preference than that of the Sun and the nucleus rotation axis.

In order to obtain the same quality of the fit for the instantaneous field of Fig. 11 as of Table 2, the order of the expansion in spherical harmonics must be carried out up to and including order N=10. This means that by adopting the mean field approach and thereby neglecting the outgassing induced short-period variations of the Doppler signal the number of parameters is reduced from more than 121 ( N=10) to $10~(N=3,\overline{b}_{nm}\equiv{}0)$.

5.6 The expansion velocity

From Eq. (37) we see that the drift velocity V of the outgassing pressure field is important if the incident gas molecules are thermalized by the spacecraft walls. From plots of the V in $\varphi $-$\theta $ space for different rotational phases one can see that the structures of the field are non-trivial, but this time the field's variability in the angles is significantly reduced. This is clearly illustrated in Fig. 15 where $\Delta{}V=V_{\rm max}-V_{{\rm min}}$ over the 72 rotational phases (one rotation period) are plotted at the distance $r=10~r_{\rm c}$ and should be compared to a global mean of $\overline{V}\sim{}650~{\rm m}~{\rm s}^{-1}$.

\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f15.eps}}
\end{figure} Figure 15:

The isocontours of the velocity amplitude over a rotation period are plotted above. Here the curves with the values $\Delta {}V=150$ (dashed, solar plane-of-sky), 100, 50, and $10~{\rm m}~{\rm s}^{-1}$ (dashed, anti-solar direction) are shown.

Open with DEXTER

Also, the simulations show that V is fairly independent of r beyond some distance. To quantify this property, we have again fitted simplified expansions like the one of Eq. (48) to the mean velocity field $\overline{V}$, calculated from the 72 available phases. Using the criteria Eqs. (50) and (51), and the simple expression

\begin{displaymath}
\overline{V}\approx{}[~635+142~\cos\theta~]~{\rm m}~{\rm s}^{-1},
\end{displaymath} (52)

we obtain the acceptable results of Table 4.

Table 4:   Fit parameter $\overline {\Delta }(r_1,r_2)$ (and $\Delta (r_1,r_2)$) for the expression Eq. (52), calculated as a least-squares approximation to $\overline{V}(r_1=25~r_{\rm c}$).

Evidently, the approximation Eq. (52) is a fair representation of the mean velocity field at other distances, but is degraded as the nucleus is approached.

6 Applications

6.1 Mean field model

The purpose of a mean field model like Table 2 is that it can be used to reduce the difference between observed and computed data during a comet rendezvous mission. If the differences, also called residuals, are similar to the observation noise, we can say that we have an adequate description of reality.

As for Rosetta, one essential data type for the radio science objectives, like the determination of the comet nucleus' gravity field, are Doppler data with an observation noise of $0.01~{\rm mm}~{\rm s}^{-1}$ for long integration times (Pätzold et al. 2007). Below we shall use the cometocentric velocity component $\dot{Z}$ of the spacecraft in the direction of the Sun as the Doppler observable, an acceptable approximation since the cometocentric direction of the Earth approximately coincides with that of the Sun at the early stages of the Rosetta mission (Mysen & Aksnes 2008). The unit vector $\vec{u}_V$, which describes the velocity direction of the spacecraft relative to the atmosphere is set to $-\vec{r}/r$.

In Fig. 16 the difference $\Delta{}\dot{Z}$ between ``real'' Doppler data, generated using the full coma simulations, and the mean field model Table 2, is shown as the solid line. Three important initial orbit elements are the cometocentric semi-major axis $a_0=20~r_{\rm c}$ (initial orbit period of 560 h), the eccentricity e0=0.2 and the orbit's inclination with respect to the cometocentric solar plane-of-sky I0=0.5. The choice of a small inclination value is motivated by stability reasons. So if the orbit is to remain bound for the duration of the studied time interval of $400~{\rm hr}$, then the surface of Rosetta exposed to outgassing has to be minimized, see Figs. 4 and 5. Still, the spacecraft's cometocentric orbit elements are seen to undergo significant changes over the time span shown in the plot. Radiation pressure from the Sun is very important, but is nevertheless not included since we choose to focus solely on the properties of the pressure forces due to outgassing from the comet nucleus.

\begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f16.eps}}
\end{figure} Figure 16:

The differences between ``real'' Doppler data from a cometary orbiter, generated using the full coma simulations, and model calculations are shown. For the dashed curve, the applied model is pure Kepler motion. The solid line is produced using the mean field of Table 2 as model. The dips in the lines are due to crossings between the ``real'' Doppler data and those simulated by the different models (Kepler, mean field).

Open with DEXTER

Although the mean field model reduces the residuals in comparison to a model which only includes Kepler motion as represented by the dashed line in Fig. 16, the resulting residuals are far from similar to one, as should be required for a proper description of reality of this maximum CO production scenario. When the $\overline {a}_{nm}$ coefficients of the mean field model are fitted to real Rosetta data in 2014, the residuals will become smaller than those of Fig. 16, provided that the prerequisites for the coma simulations and reality are the same. However, the solve for coefficients $\overline {a}_{nm}$ will then not represent the global field as those of Table 2, but the local pressure forces in parts of the coma which the probe has traversed.

The mean field approach is improved as the probe's distance to the comet increases.

6.2 Determination of gravity field

As previously mentioned, a stated goal of the Rosetta Radio Science Investigations is the determination of the nucleus' gravity field from Doppler data (Pätzold et al. 2007). The cometocentric gravitational acceleration $\ddot{\vec{r}}$ of an orbiting spacecraft can be obtained by taking the gradient of the gravitational potential

\begin{displaymath}
\ddot{\vec{r}}=-\nabla{}(V_0+V_1),\quad{}V_0=-\frac{\mu_{\rm c}}{r},
\end{displaymath} (53)

where $\mu_{\rm c}$ is the product of the constant of gravitation and the nucleus mass. The most important contributions from the part of the field which corotates with the nucleus are given by (Kaula 2000)

\begin{displaymath}
V_1=\frac{\mu_{\rm c}}{r}~\left(\frac{r_{\rm c}}{r}\right)^2...
...P_{20}(\sin\beta)+c_{22}P_{22}(\sin\beta)\cos2\lambda~\right].
\end{displaymath} (54)

Here, P2m are Legendre functions, and $\lambda$ and $\beta$ are the equatorial east longitude and north latitude, respectively, of the probe relative to the corotating nucleus. See for instance Mysen & Aksnes (2005) for more details. The parameters to be found from Doppler data are then $\mu_{\rm c}$ (treated in Mysen & Aksnes 2008), and the gravity coefficients c20 and c22.

6.2.1 Nucleus oblateness

The Doppler signals related to the two latter parameters are mainly experienced through variations in the probe's cometocentric velocity component in the direction of the observer (Mysen & Aksnes 2008), $\dot{Z}$. To investigate whether or not these gravity induced variations are drowned in the corresponding variations induced by outgassing pressure, we have mapped the difference

\begin{displaymath}\Delta{}\dot{Z}=\dot{Z}_{{\rm Kepler+out./obl.}}-\dot{Z}_{\rm Kepler},
\end{displaymath} (55)

where the first component is generated using Kepler acceleration and the full outgassing pressure model (solid curve of Fig. 
17), or Kepler acceleration and the oblateness parameter c20=0.088 (arbitrary choice) of Eq. (54) (dashed curve of Fig. 17). From the simulations we found that the the outgassing pressure prefactor $\varrho{}V^2$ had to be down-scaled one order of magnitude for a meaningful comparison to be made at realistic orbit sizes. In addition, the initial orbit inclination with respect to the solar plane-of-sky had to be set to the small value I0=0.2 to minimize the outgassing exposed surface of Rosetta over an orbit period.
\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f17.eps}}
\end{figure} Figure 17:

Doppler signals related to the outgassing pressure (solid) and the nucleus' c20 gravity field (dashed) are shown for an initial orbit semi-major axis of $a_0=5~r_{\rm c}$ and a pressure field that is scaled down one order of magnitude in comparison to the previously presented worst-case scenario.

Open with DEXTER
\begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f18.eps}}
\end{figure} Figure 18:

The oscillations of the Doppler signals related to the short-periodic variations of the outgassing pressure field (solid) and the c22 gravity field (dashed) are shown as the spacecraft travels through the night side of the coma, with an initial orbit semi-major axis $a_0=10~r_{\rm c}$. Again, the pressure fields used were scaled down one order of magnitude in comparison to the worst-case scenario presented earlier in this paper. Notice the short time scale in comparison to Fig. 17.

Open with DEXTER

If we make the conservative assumption that we have no a priori knowledge of the properties of the coma, i.e., that realistic coma simulations are not available, we can see that the signal related to c20 begins to become discernible at the cometocentric distance $r=5~r_{\rm c}$ of Fig. 17. Above this height, the component is drowned in the outgassing signal since the acceleration due to this gravity term drops as $\sim$r-4, while the outgassing induced acceleration diminishes as $\sim$r-2. That is, the amplitude of the periodic variations induced by a perturbing acceleration $\ddot{\vec{r}}$, with a frequency equal to the orbit frequency $n=\mu^{1/2}~a^{-3/2}$, can be approximated by the expression

\begin{displaymath}
\Delta{}\dot{Z}\sim{}\int{}{\rm d}t~\vert\ddot{\vec{r}}\vert...
...ec{r}}\vert\sim{}\frac{1}{n}~\Delta{}\vert\ddot{\vec{r}}\vert,
\end{displaymath} (56)

where f is the true anomaly of the spacecraft in its orbit, and $\Delta{}\vert\ddot{\vec{r}}\vert$ is the amplitude of the acceleration during an orbit period. Hence, one obtains

\begin{displaymath}\Delta{}\dot{Z}_{{\rm out.}}\sim{}\Delta{}\chi~r^{-1/2},\quad{}
\Delta{}\dot{Z}_{{\rm obl.}}\sim{}c_{20}~r^{-5/2},
\end{displaymath} (57)

where $\Delta {}\chi $ is the amplitude of the function Eq. (46) over an orbit revolution, corrected for probe surface exposure to outgassing. The increasing oscillation amplitudes of Fig. 17 are mainly caused by model differences and not by an increase in the spacecraft's distance to the comet.

6.2.2 Nucleus triaxiality

As for the $c_{22}\equiv{}-0.044$ (arbitrary choice) gravity coefficient related to the nucleus' triaxiality, we have studied the velocity differences

\begin{displaymath}\Delta{}\dot{Z}=\dot{Z}_{{\rm Kepler+out.(+triax.)}}-\dot{Z}_{\rm Kepler+mean.}
\end{displaymath} (58)

to see if one can discern the c22 induced variations. We compared the difference in the Doppler signal between a full outgassing model and the mean outgassing model (solid curve of Fig. 
18) to the result when the full model is extended to include the c22 component (dashed curve of Fig. 18). Again, the pressure field prefactor $\varrho{}V^2$ is down-scaled one order of magnitude in comparison to the coma simulations. The initial orbit semi-major axis is $a_0=10~r_{\rm c}$, with an orbit geometry in a way that the spacecraft travels through the night side of the comet nucleus, where the short-periodic variations of the pressure field are small, see Fig. 13. The arrays are still defined to point towards the Sun though. Thus with a proper array orientation, the outgassing perturbation can be reduced one order of magnitude. Accordingly, we see that the c22 oscillations are clearly discernible.

From a similar approach as that of Eq. (56), we derive the amplitude of the oscillations of Fig. 18 as

\begin{displaymath}\Delta{}\dot{Z}\sim{}\int{}{\rm d}t~\vert\ddot{\vec{r}}\vert\sim{}\frac{1}{\omega{}}~\Delta{}\vert\ddot{\vec{r}}\vert,
\end{displaymath} (59)

where $\omega=2\pi/P$ is the nucleus rotation frequency and $\Delta{}\vert\ddot{\vec{r}}\vert$ is the amplitude of the acceleration in question over the period P. As a result

\begin{displaymath}
\Delta{}\dot{Z}_{{\rm out.}}\sim{}\frac{\Delta{}\chi}{\omega...
...Delta{}\dot{Z}_{{\rm tri.}}\sim{}\frac{c_{22}}{\omega}~r^{-4},
\end{displaymath} (60)

where $\Delta {}\chi $ this time is the amplitude of the function Eq. (46) over a nucleus rotation period, again corrected for the surface exposure of the probe to outgassing.

It is important to note that the outgassing induced variations of the Doppler signal, with a period similar to the rotation period of the nucleus, and the estimated amplitude given by Eq. (60) cannot be removed from the data by the mean field approach of Eq. (48). According to Eq. (60), the amplitude of these outgassing induced Doppler signal variations are several orders of magnitude larger than the Doppler observation noise in the day side of the coma at cometocentric distance $10~r_{\rm c}$. Consequently the mean field model alone is not an adequate description of reality in general during the Rosetta mission for the presented maximum CO production scenario. However, as already shown in Fig. 18, at the same distance and in the night side of the coma, the outgassing induced variations of the Doppler signal are comparable to the Doppler noise, if the CO pressure field is down-scaled one order of magnitude in comparison to the maximum production scenario.

7 Discussion

In support of navigation and radio science during ESA's Rosetta mission, we have reviewed a simple model for the interaction between a spacecraft and a Maxwellian atmosphere and applied it to the Rosetta spacecraft with some simplifying assumptions with respect to its outer shape, molecular scattering properties, and temperature.

A state-of-the-art 3D+t gasdynamical simulation of the coma around the irregularly shaped nucleus of 67P was used to analyze the interaction. The simulation parameters are selected to yield upper-limit 3 AU total production rates of water - by solar-driven sublimation - and of CO (which dominates at this heliocentric distance). The production of CO was assumed to be moderately uneven over the surface, giving birth to a moderately contrasted pressure field, in comparison with what would have been obtained if the production had been assumed concentrated in a localized area of the nucleus. The spatial resolution of the model was three degrees in longitude and latitude, and the time resolution was 1/72 of the rotation period.

The present results have only been established for probe distances between 5 and 50 mean comet radii.

First of all, it was found that the probe distance to the nucleus must be larger than 5-10 mean comet radii if the mean free path of the molecules in the cometocentric direction of the Sun is to be much larger than the size of the spacecraft. Below this distance, one may need to include the spacecraft as part of the coma simulation to obtain reliable results.

It was also found that the pressure force may be dependent on the coma temperature in the night-side of the coma, and that beyond about 5-10 mean comet radii outwards, the gas flow direction was mainly directed away from the nucleus.

We found that the pressure force on Rosetta is amplified by a non-zero spacecraft temperature and that the force does not act mainly away from the comet nucleus in general. Three parameters were suggested to describe the interaction, assuming that one type of molecule is the main contributor to the pressure.

For probe orbit sizes above 5-10 mean comet radii, the pressure field was seen to drop as the inverse square of the distance to the comet nucleus. This separability of the pressure field dependency in probe direction and probe distance simplifies the search for representative pressure field models. However, the force field is significantly distorted due to the non-spherical shape of the comet nucleus. As a result, deterministic models which aim to represent the properties of the instantaneous coma will contain too many parameters and will therefore not be practical.

The average pressure field over a nucleus rotation period, on the other hand, was seen to be well represented by a relatively small number of parameters. However, for the ``nearly maximum'' gas production scenario studied, it was established that this mean field model alone was not able to describe reality to sufficient precision during the Rosetta mission.

In the general direction of the Sun, the maximum drag force on Rosetta was seen to be larger than the gravitational attraction from the nucleus at all probe distances.

As for the determination of the gravity field coefficients related to the nucleus oblateness and triaxiality through Doppler data, it was seen that it was easier to discern the latter from the ``noise'' induced by outgassing pressure forces on Rosetta. The reason for this is that the Doppler signature of the nucleus triaxiality is resolved over a shorter time interval compared to that of the nucleus oblateness. As a result, the nucleus triaxiality can be determined in a limited region of the coma where the variability of the outgassing pressure field, or ``noise'', is small. The Doppler signal related to the nucleus oblateness is, on the other hand, only resolved on a time scale similar to that of the spacecraft's orbit period. During this time scale, the probe will have experienced the global, and not only local, properties of the coma. If the nucleus rotation is not uniaxial as assumed above, the gravity field related to the nucleus oblateness could also be discernible on a time scale similar to the nucleus rotation period (Mysen et al. 2006).

As similar state-of-the-art coma simulations based on Rosetta observations will be available in 2014-2015, it should to a large extent be possible to remove the signal related to outgassing pressure forces from Doppler data. As a result, it should be easier to determine the gravity field coefficients of the nucleus.

It should be noted that in case the rotation of the target nucleus is not found to be uniaxial, the nucleus rotation will nevertheless be approximately periodic (also called quasi-periodic) on a limited time scale. Coma simulations can then be executed covering nucleus rotation phases over this approximate rotation period.

Acknowledgements
This work was possible due to supports from the Research Council of Norway, project 184686, and, in France, from the CNES ``Rosetta Croisière'' funding.

References

  1. Byram, S. M., Scheeres, D. J., & Combi, M. R. 2007, J. Guid. Cont. Dyn., 30, 5, 1445 [Google Scholar]
  2. Capria, M. T., Coradini, A., De Sanctis, M. C., & Fulle, M. 2004, in The New ROSETTA Targets, ed. L. Colangeli, E. M. Epifani, & P. Palumbo (Dordrecht: Kluwer), ASSL, 311, 177 [Google Scholar]
  3. Crifo, J. F. 1995, A&A, 300, 597 [NASA ADS] [Google Scholar]
  4. Crifo, J. F. 2006, Adv. Space Res., 38, 1911 [NASA ADS] [CrossRef] [Google Scholar]
  5. Crifo, J. F., Rodionov, A. V., Szegö, K., & Fulle, M. 2002, Earth, Moon and Planets, 90, 227 [Google Scholar]
  6. Crifo, J. F., Fulle, M., Kömle, N. I., & Szegö, K. 2004a, in Comets II, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson: The University of Arizona Press), 471 [Google Scholar]
  7. Crifo, J. F., Lukyanov, G. A., Zakharov, V. V., & Rodionov, A. V. 2004b, in The New ROSETTA Targets, ed. L. Colangeli, E. M. Epifani, & P. Palumbo (Dordrecht: Kluwer), ASSL, 311, 119 [Google Scholar]
  8. Davidsson, B. J. R., & Gutiérrez, P. J. 2005, Icarus, 176, 453 [NASA ADS] [CrossRef] [Google Scholar]
  9. Feaga, L. M., A'Hearn, M. F., Sunshine, J. M., Groussin, O., & Farnham, T. L. 2007, Icarus, 190, 345 [NASA ADS] [CrossRef] [Google Scholar]
  10. Harrison, I. K., & Swinerd, G. G. 1996, Planet. Space Sci., 44, 2, 171 [CrossRef] [Google Scholar]
  11. Kaula, W. M. 2000, Theory of Satellite Geodesy (New York: Dover Publications) [Google Scholar]
  12. Lamy, P. L., Toth, I., Davidsson, B. J. R., et al. 2007, Space Sci. Rev., 128, 23 [NASA ADS] [CrossRef] [Google Scholar]
  13. Moe, K., & Moe, M. M. 2005, Planet. Space Sci., 53, 793 [NASA ADS] [CrossRef] [Google Scholar]
  14. Moore, P., & Sowter, A. 1991, Planet. Space Sci., 39, 10, 1405 [Google Scholar]
  15. Mysen, E., & Aksnes, K. 2006, A&A, 455, 1143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Mysen, E., & Aksnes, K. 2008, A&A, 490, 447 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Mysen, E., Olsen, Ø, & Aksnes, K. 2006, Planet. Space Sci., 54, 750 [NASA ADS] [CrossRef] [Google Scholar]
  18. Pätzold, M., Häusler, B., Aksnes, K., et al. 2007, Space Sci. Rev., 128, 599 [NASA ADS] [CrossRef] [Google Scholar]
  19. Prialnik, D.,Benkhoff, J., & Podolak, M., in Comets II, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson: The University of Arizona Press), 359 [Google Scholar]
  20. Rodionov, A. V., & Crifo, J. F. 2006, Adv. Space Res., 38, 1923 [NASA ADS] [CrossRef] [Google Scholar]
  21. Rodionov, A. V., Crifo, J. F., Szegö, K., Lagerros, J., & Fulle, M. 2002, Planet. Space Sci., 50, 983 [NASA ADS] [CrossRef] [Google Scholar]
  22. Scheeres, D. J., Bhargava, S., & Enzian, A. 2000, A Navigation Model of the Continuous Outgassing Field Around a Comet, TMO Progress Report 42 [Google Scholar]
  23. Valorge, C. 1995, in Spaceflight Dynamics Part I, ed. J.-P. Carrou, CNES, 217 [Google Scholar]
  24. Weeks, C. J. 1995, J. Astronaut. Sci., 43, 3, 327 [Google Scholar]
  25. Zakharov, V. V., Rodionov, A. V., Lukyanov, G. A., & Crifo, J.-F. 2008, Icarus 194, 327 [CrossRef] [Google Scholar]

All Tables

Table 1:   Comet 67P model parameters.

Table 2:   Fit adequacy parameter $\overline {\Delta }(r_1,r_2)$ (and $\Delta (r_1,r_2)$) for an N=3 expansion Eq. (48) fitted to $\overline{\chi}(r_1=25~r_{\rm c}$).

Table 3:   Coefficients $\overline {a}_{nm}$ of an expansion in spherical harmonics fitted to $\overline{\chi}(r=25~r_{\rm c})$ of Eq. (46).

Table 4:   Fit parameter $\overline {\Delta }(r_1,r_2)$ (and $\Delta (r_1,r_2)$) for the expression Eq. (52), calculated as a least-squares approximation to $\overline{V}(r_1=25~r_{\rm c}$).

All Figures

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f01.eps}}
\end{figure} Figure 1:

A flat plate is seen from its side. The plate moves with a velocity $\vec{V}$ relative to the atmosphere. The unit vector $\vec{n}\equiv{}\vec{u}_\bot{}$ is normal to the surface. The unit vectors $\vec{u}_D$ and $\vec{u}_L$ define the drag and lift directions respectively. As for the attack angle, it is counted positive from $\vec{n}$ towards $\vec{u}_\parallel$. The vector $\vec{u}_\parallel$ is parallel to the plate, lies in the plane spanned by $\vec{n}$ and $\vec{V}$, and must be defined in a way that Eq. (3) is valid.

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f02.eps}}
\end{figure} Figure 2:

The pressure force coefficient CV of Eq. (21) is given for different values of s, dashed lines. From bottom to top we have s=1, s=3, and $s\rightarrow {}\infty $ (solid), Eq. (26).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f03.eps}}
\end{figure} Figure 3:

From bottom to top, $C_\Delta $ of Eq. (23) is given for s=1, s=3, and $s\rightarrow {}\infty $ (solid), Eq. (27).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f04.eps}}
\end{figure} Figure 4:

The figure illustrates effective areas of Rosetta in the drag direction: total (solid), without bus (dotted), or with the spacecraft temperature set to $T_f=0~{\rm K}$ (densely dotted).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f05.eps}}
\end{figure} Figure 5:

The lift to drag ratio is shown above as a function of the attack angle $\vartheta $.

Open with DEXTER
In the text

  \begin{figure}
\includegraphics*[scale=0.25,angle=0]{12107f06.eps}
\end{figure} Figure 6:

Present best-fit shape of the comet 67P nucleus. The largest diameter is about 5 km. The uncertainty is probably on the order of the magnitude of the topographic features. The corresponding best-fit monoaxial rotation is about the z-axis. The model will be periodically improved as new observational data are accumulated. After Lamy et al. (2007).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[scale=0.25,angle=-90]{12107f71.eps}\hspace*{...
...\hspace*{4mm}
\includegraphics[scale=0.25,angle=-90]{12107f79.eps}
\end{figure} Figure 7:

Variation of the total mass density in the equatorial plane during the rotation. From left to right and from top to bottom: isocontours of Log10 (mass density, kg/m3) in the equatorial plane are shown every 1/9 of the rotation period (i.e., every 1.4 h.) The projected Sun direction is to the right. Notice the clockwise rotation of the nucleus cross section. On this and the two next figures, references and dates appear at the bottom left of each panel for safety of identification: the reader must overlook them. The dust contribution has not been included.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=17cm,clip]{12107f08.eps}
\end{figure} Figure 8:

Log10 of the H2O ( upper panels) and CO ( lower panels) density in two orthogonal planes at rotation phase 0. The projected Sun direction is to the right on the X-axis ( right panels) or at 200 to the X-axis ( left panels). The total density shown on Fig. 7 ( upper left panel) is the sum of the two present right panel densities.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=16.5cm,clip]{12107f09.eps}
\vspace*{-2mm}
\end{figure} Figure 9:

Velocity ( upper panels) and Mach number ( lower panels) in two orthogonal planes at rotation phase 0. The projected Sun direction is to the right on the X-axis ( right panels) or at 200 to the X-axis ( left panels).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f10.eps}}
\end{figure} Figure 10:

``Pseudo'' mean free paths for water (solid lines) and CO (dashed lines) are shown, expressed in multiples of the characteristic orbiter size $D_{{\rm S/C}}=32~{\rm m}$ as a function of cometocentric distance expressed in multiples of $r_{\rm c}=2$ km along the subsolar axis. The time is the one we defined as the initial rotational phase of the nucleus. Two curves are plotted for each molecule, corresponding to the two minimum ad hoc temperatures $T_{{\rm min}}=15~{\rm K}$ ( lower) and $T_{{\rm min}}=50~{\rm K}$ ( upper).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f11.eps}}
\end{figure} Figure 11:

The contours in $\varphi $-$\theta $ space on which the pressure field prefactor $\chi\equiv{}\xi~\varrho{}V^2~(r/r_{\rm c})^2$ is constant are shown for a nucleus rotation phase of $0^\circ {}$. Above, $\chi $ takes on the values 0.1 (in the anti-solar direction), 0.25, 0.5, 0.75, and with the peak $1~r_{\rm c}~{\rm h}^{-2}$ in the general direction of the Sun, $\theta \sim {}0-\pi /4$. The different types of curves represent different cometocentric distances: r=10 (dashed), 25 (solid) and $50~r_{\rm c}$ (densely dotted).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f12.eps}}
\end{figure} Figure 12:

The isocontours of the amplitude $\Delta {}\chi $ over a rotation period are shown. Above, $\Delta {}\chi $ takes on the values 0.05 (dashed, $\theta \sim {}3\pi /4$) 0.1, 0.25, 0.5, 0.75 and $1~r_{\rm c}~{\rm h}^{-2}$ (dashed, $\theta \sim {}\pi /4$). The probe's cometocentric distance is set to $r=10~r_{\rm c}$ for convenience.

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f13.eps}}
\end{figure} Figure 13:

The variation of $\chi $ with nucleus rotation for $(\varphi ,\theta )=(\pi ,\pi /4)$ (dashed), and for $(\varphi ,\theta )=(0,\pi /4)$ (solid). The probe's cometocentric distance has for convenience been set to $r=10~r_{\rm c}$.

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f14.eps}}
\end{figure} Figure 14:

The curves on which the time averaged pressure field $\overline {\chi }$ is constant are shown for the values 0.075 ( $\theta \sim {}\pi $), 0.1, 0.25, 0.5, 0.75, 1, and $1.2~r_{\rm c}~{\rm h}^{-2}$ ( $\theta \sim {}0$). Again, the different types of curves represent different cometocentric distances: r=10 (dashed), 25 (solid) and $50~r_{\rm c}$ (densely dotted).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f15.eps}}
\end{figure} Figure 15:

The isocontours of the velocity amplitude over a rotation period are plotted above. Here the curves with the values $\Delta {}V=150$ (dashed, solar plane-of-sky), 100, 50, and $10~{\rm m}~{\rm s}^{-1}$ (dashed, anti-solar direction) are shown.

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{9cm}{!}{\includegraphics{12107f16.eps}}
\end{figure} Figure 16:

The differences between ``real'' Doppler data from a cometary orbiter, generated using the full coma simulations, and model calculations are shown. For the dashed curve, the applied model is pure Kepler motion. The solid line is produced using the mean field of Table 2 as model. The dips in the lines are due to crossings between the ``real'' Doppler data and those simulated by the different models (Kepler, mean field).

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f17.eps}}
\end{figure} Figure 17:

Doppler signals related to the outgassing pressure (solid) and the nucleus' c20 gravity field (dashed) are shown for an initial orbit semi-major axis of $a_0=5~r_{\rm c}$ and a pressure field that is scaled down one order of magnitude in comparison to the previously presented worst-case scenario.

Open with DEXTER
In the text

  \begin{figure}
\par\resizebox{8.8cm}{!}{\includegraphics{12107f18.eps}}
\end{figure} Figure 18:

The oscillations of the Doppler signals related to the short-periodic variations of the outgassing pressure field (solid) and the c22 gravity field (dashed) are shown as the spacecraft travels through the night side of the coma, with an initial orbit semi-major axis $a_0=10~r_{\rm c}$. Again, the pressure fields used were scaled down one order of magnitude in comparison to the worst-case scenario presented earlier in this paper. Notice the short time scale in comparison to Fig. 17.

Open with DEXTER
In the text


Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.