© ESO, 2012

1. Introduction

Back in the 19th century, cometary observers noted that the orbital period of comets was subject to significant changes from one orbit to another. Having predicted the return of the comet 2P/Encke a few hours too late, Encke (1826) proposed that the advance of the comet should have been due to a resisting medium causing a deceleration of the comet and therefore a decrease in its orbital period. As this resisting medium does not affect the planet motion, Bessel (1836) was the first to suggest that the orbital period changes stemmed from the ejection of material by a solid nucleus. When a comet approaches the Sun, the solar radiation is used to sublimate ices anisotropically, which triggers a time-dependent recoil force. Because observers noticed that the orbital period could also grow, Bessel (1836) proposed that the outgassing may act more on one side of perihelion.

To describe the trajectories of cometary nuclei more accurately, Marsden (1969) introduced the first model of the non-gravitational accelerations in the equations of motion of the comet. This model was later on improved by Marsden et al. (1973), who assumed that the sublimation of water ice acts symmetrically with respect to perihelion on a fast-rotating nucleus. The nucleus is also assumed to be spherical, and the whole surface is assumed to be isotropically outgassing. Marsden et al. (1973) introduced the dimensionless function g(r_h), which represents the variation in the sublimation rate as a function of the heliocentric distance r_h of the comet: (1)The non-gravitational acceleration vector is decomposed into a radial (A₁ g(r_h)), a transverse (A₂ g(r_h)), and a normal (A₃ g(r_h)) component in the orbital reference frame. The three free parameters A₁, A₂, and A₃ are solved in a least square fit of astrometrical observations, together with the six Keplerian osculating orbital elements. In this model, the parameter A₂ models the effect of the thermal inertia on a rotating nucleus, which introduces a lag angle between the direction of the Sun and that of the non-gravitational force (Marsden et al. 1973). This semi-empirical model soon became the standard to describe the non-gravitational acceleration. It has been used successfully since the seventies to significantly improve the accuracy of the comets ephemeris. This model has been used for an important study of the comet 1P/Halley by Rickman & Froeschle (1983a,b) with respect to the need of accurate ephemerides for the observations and the Giotto spacecraft flyby during the 1986 apparition. Then, they extended this study in a general framework to all the short-period comets in adding a thermal model interpretation of the cometary nuclei (Froeschle & Rickman 1986).

However, measurements of the dust and gas production rates around perihelion acquired during the last perihelion return of comet 1P/Halley confirmed that the ice sublimation can be asymmetric with respect to perihelion (Schleicher et al. 1998). To take the observed asymmetry into account, Yeomans & Chodas (1989) replaced the function g(r_h) by , where . The free parameter Δt allows shifting the maximum of the gas production rate with respect to the perihelion passage. This model is very useful for improving the ephemeris in the case of comets presenting a large asymmetry in their outgassing. For comet 6P/d’Arrest, the best rms residual (1.473 arcsec versus 1.857 arcsec with the Marsden et al. 1973 model) were obtained for Δt = 40 days (Yeomans & Chodas 1989). It should be noted that this model is commonly used by the JPL to calculate cometary ephemeris.

In parallel, electronic imaging data unambiguously revealed the presence of dust and gas structures in the coma of several comets. Dust jets and filaments also showed up in close-up images of the inner coma of 1P/Halley by the imaging systems aboard the Giotto and Vega space probes (Keller et al. 2004).

Since then, a succession of space mission have repeated these observations:

• Deep Space 1 on 19P/Borrelly (Soderblom et al. 2002);

• Stardust on 81P/Wild 2 (Sekanina et al. 2004);

• Deep Impact on 9P/Tempel 1 (A’Hearn et al. 2005);

• and EPOXI (Deep Impact extended mission) on 103P/Hartley 2 (A’Hearn et al. 2011).

These structures were interpreted as products by “active areas” (Sekanina 1993) or more recently also by concavities (Crifo et al. 1999) on the surface of a rotating nucleus. This led to the hypothesis that “seasonal effects” of the successive illumination of different areas with different properties on the nucleus’s surface can cause the production rate to vary during a perihelion passage, thus affecting the non-gravitational acceleration (Sekanina 1993). This led Sekanina (1993) to propose a “spotty model” of non-gravitational acceleration in which the perihelion asymmetry of the gas production rate was explained by the outgassing of discrete sources distributed at the surface of the nucleus. In his model, the maximum outgassing of a given active area took place when the local elevation of the Sun above the area was at its maximum. The non-gravitational acceleration is expressed as the sum of rotation-averaged contributions from all the active sources taking their location at the surface of a spherical nucleus into account. The other parameters involved in his model are the direction of the spin axis, the mass of the nucleus, the surface of the active sources, and their sublimation rate as a function of the heliocentric distance. The free parameters were then fitted until the non-gravitational acceleration got close to what is predicted by the standard model, in which the parameters A_i best-fit the astrometrical data. Sekanina (1993) noted a correlation between the sign of A₂, the locations of the active sources, and the asymmetry of the outgassing. He also pointed out that, unlike the standard model, a negative value of A₁ can be produced by an active area located close to a pole for a given orientation of the spin-axis.

The first attempt to directly introduce the rotation-averaged components of the non-gravitational acceleration of the spotty model into the equations of motion of a comet was made by Sitarski (1990) for comet 22P/Kopff. The parameters involved in the model are the cometocentric latitude β_j of the jth active regions, the constants A_j proportional to the active area, and the mass of the nucleus and the three angles η, I, and φ, which are respectively the lag angle of the outgassing maximum behind the subsolar meridian, the obliquity of the orbital plane to the nucleus equator, and the solar longitude at perihelion. These parameters, together with a set of osculating orbital elements, are directly fitted by minimizing a chi-square function that represents the distance between the modelled and the measured astrometrical positions.

The last and most sophisticated model was proposed by Davidsson & Gutiérrez (2004). They modelled the nucleus with a prolate ellipsoid, and the surface is divided into thousands of facets. Each facet has its own thermic equilibrium. They tested thousands of models where each facet is chosen randomly to be active or not. They fit their models on the observed water production rate. An important result of Davidsson & Gutiérrez (2004) is that they can reproduce the asymetric production rate curve. With the knowledge of the outgassing rate, they computed the non-gravitational force and therefore the mass of the nucleus.

In this article, we propose a new implementation of the “spotty” model like Sekanina (1993) into the equations of motion of a comet. Compared to the implementation of Sitarski (1990), our model has the following specific characteristics.

• We consider latitudinal bands on the surface of an ellipsoidalnucleus, and we fit one parameter per band, which is thepercentage of active surface. Our idea is to develop a model with afew parameters with a physical meaning.

• We model the water production rate from the energy balance equation.

The very low thermal inertia of cometary nuclei (e.g. Julian et al. 2000; Groussin et al. 2007) points toward a significant influence by the seasonal effects compared to that of the thermal inertia in the asymmetric behaviour of the production rates.

In Sect. 2, we give a thorough description of the model and of the algorithms used to retrieve the parameters. In Sect. 3, we apply our model to the astrometrical measurements of comet 19P/Borrelly and we compare the results with those of the standard model. We conclude in Sect. 4.

2. Description of the model

The spotty models are based on determining the latitude and active surface for each spot. Averaging the effects of this spot over one rotation is equivalent to a model with latitudinal bands that are more or less active. Therefore, we restrict our model to a few bands in order to simplify it.

2.1. Expression of the non-gravitational acceleration

The nucleus is modelled as a triaxial ellipsoid divided into latitudinal bands (Fig. 1). The latitude is defined as the angle between the normal to the surface of the ellipsoid and the plane perpendicular to the nucleus rotational axis. The thermal inertia of the nucleus is neglected and the gas velocity is considered to be proportional to the thermal gas velocity. These hypotheses allow seasonal effects.

The goal at the beginning is to compute the maximum non-gravitational force for each band if the whole band is active. After that it will be possible to fit a coefficient that represents the percentage of the active zone. For each point of the surface, we have to compute the local non-gravitational force in order to, at the end, integrate it on one nucleus rotation. We need to compute the gas production rate and the ejection velocity.

2.1.1. The surface temperature

The surface temperature is needed to compute the water sublimation rate for each surface element of the model. The following surface energy balance is solved by dichotomy to obtain the temperature of a point on the nucleus surface (2)with

• A_b: the Bond albedo (equal to the product of the geometric albedo by the phase integral);

• F_⊙: the solar constant (W   m^-2);

• r_h: the heliocentric distance (AU);

• z: the elevation of the Sun;

• η: the beaming factor introduced by Lebofsky & Spencer (1989);

• ϵ: the nucleus infrared emissivity;

• σ: the Stefan-Boltzmann constant (J   K^-4   m^-2   s^-1);

• T: the surface temperature (K) depending on the time and the latitude of the point on the surface;

• L(T): the sublimation latent heat (J   mol^-1);

• : the Avogadro constant (mol^-1);

• Z(T): the surfacic sublimation rate (molecules   s^-1   m^-2).

The elevation of the Sun z is the angle between the sub-solar point and a point on the nucleus surface. It is defined as (3)with

• (θ_⊙,ϕ_⊙): the coordinates in the cometocentric frame of the sub-solar point;

• (θ,ϕ): the considered point on the nucleus.

The axis of the cometocentric frame is along the pincipal axis of inertia of the ellipsoidal nucleus.

2.1.2. Water sublimation rate

The water sublimation rate is the key value in computing non-gravitational forces because with the ejection velocity, it gives the quantity of motion. The sublimation latent heat of the water ice is interpolated as the following three-degree polynomial from the values of Washburn (1928)(4)The surfacic sublimation rate is given by (5)where M_H2O is the water molecule mass (kg), k the Boltzmann constant (J   K^-1), α a coefficient for the recondensation of water ice on the surface introduced by Crifo (1987) (we adopted Crifo’s value: α = 0.25), and P(T) is the water saturation vapor pressure (Pa) given by Fanale & Salvail (1984) as (6)with A = 3.56 × 10¹²  Pa and T₀ = 6141.6607  K.

We consider a surfacic sublimation rate that is equal in each latitudinal band point. Thus, the maximal water production rate of a band is the product of the rotation-averaged surfacic sublimation rate (over one rotation period P) of a point situated at the mid-latitude θ_j and the band total surface S_j. The water production rate of the comet Q_H2O(t) is given by (7)where N_b is the number of bands, and C_j are coefficients describing the activity of each surface band. (C_j = 0 if the band is not outgassing, C_j = 1 if the whole band is outgassing.)

2.1.3. Non-gravitational accelerations

From the water production rate, we derive the non-gravitational acceleration with the Newton second law of motion (8)where is the surfacic force of the jth band depending on the sublimation rate Z_j and the thermal gas velocity V_gj(9)where M_H2O is the water molecular mass and N_i is the surface normal.

We used the expression of the gas ejection velocity proposed by Crifo (1987): (10)with k the Boltzman constant (k = 1.380662 × 10^-23 J K^-1), and η a local parameter varying in relation with the position on the cometary nucleus. This parameter allows us to consider that the comet surface is not composed of pure ice and is not perfectly smooth. Crifo (1987) suggests choosing this parameter in the interval  [0.39,0.5] . We fix the value of this parameter to 0.45 for the application to comets. It should be noted that, from the expression of the non-gravitational acceleration, the determination of the nucleus mass is inversely proportional to this parameter.

This expression of the non-gravitational acceleration can be introduced into the equation of the comet’s motion after decomposition into three components in the equatorial heliocentric reference frame.

2.2. Equations of motion

We develop the equations of motion in an equatorial heliocentric reference frame (P_⊙,x,y,z). The variables are the Cartesian coordinates of the position (r_h) of the comet. The differential equations of the motion of the comet can be written as (11)with

• r_h: heliocentric position vector of the comet;

• G: Gauss gravitation constant;

• M_⊙: solar mass;

• m_i: mass of the ith perturbating body;

• r_i: heliocentric position vector of the ith body;

• R: acceleration due to relativistic effects;

• A_NG: non-gravitational acceleration (Sect. 2.1).

The acceleration due to relativistic effects is given by Beutler (2005)(12)where c is the light velocity. It is important to introduce this acceleration in the equation of motion because, for a comet passing close to the Sun, it can modify the value of the non-gravitational accelerations.

We fit the orbital motion on astrometric measurements and derive the parameters.

2.3. The determination of the nucleus mass

From Eq. (7), we can see that the water production rate is proportional to the coefficients C_i. Nevertheless, thanks to astrometric measurements, we can determine the ratio C_i/M_c. Thus, we can introduce these ratios into the water production rate in order to deduce the nucleus mass: (13)where Q_obs(t_j) is the water production rate observed, the maximal water production rate calculated for the ith band at the same time t_j, and N_sub and N_b are the number of observations and the number of bands, respectively.

The nucleus mass and the corresponding error are determined as the mean value and the standard deviation of a Monte-Carlo distribution taking the precision on the observations and errors on the ratio C_i/M_c into account.

3. Application to comet 19P/Borrelly

The comet 19P/Borrelly was discovered in 1904 by the French astronomer Alphonse Borrelly at the Marseilles observatory (Borrelly 1905). With its short period (6.75 years) and low inclination (30°), it belongs to the Jupiter family. It is a well-documented comet observed at many of its returns. It was the target of Deep Space 1 mission in 2001 (Soderblom et al. 2002).

We also chose this comet for its stable axis (Schleicher et al. 2003) to avoid problem due to a moving axis.

3.1. Observational data and parameters of the model

This present investigation is based on two kinds of observations: astrometric positions and water production rates observations. Astrometric measurements were taken from the IAU Minor Planet Center. The observational material contains 2206 observations from June 12, 1994 to June 16, 2009. These observations cover three orbits of the comet but are not evenly distributed (Table 1). The number of observations was increased since 2001 perihelion because of the support to the space mission Deep Space 1.

To determine the mass of the nucleus, we used a compilation of 219 observations of the water production rates covering the same time period as astrometrical positions (see Maquet 2012, for a detailed table of these observations). This compilation consists of data sets coming from different methods of observations of water or of its photolysis products:

• direct observations of water from its submillimetric line at 557 GHz with the Odin satellite (Lecacheux et al. 2003) (2001 passage);

• ground-based observations of OH from its 18-cm lines, with the Nançay radio telescope (Bockelée-Morvan et al. 2004) (1994 and 2001 passages);

• ground-based observations of the OH radical in the near UV (Schleicher et al. 2003) (2001 passage);

• space observations of the Lyman-α line of hydrogen with the SOHO/SWAN instrument (Combi et al. 2011) (2001 passage);

• observations of the OI forbidden line in the visible (Fink 2009) (1994 passage).

The OH production rates were converted to water production rates by assuming Q_H2O = 1.1Q_OH (Huebner et al. 1992). The water production rate is observed to peak around 6 × 10²⁸ molecules/s about 20 days before perihelion (Fig. 5). We see that Lyman-α observations give significantly different results than the other methods.

In this work, we consider the nucleus of comet 19P/Borrelly as a triaxial prolate ellipsoid (a = 4.6 km, b = 1.8 km, c = 1.6 km). rotating in 26 h (Lamy et al. 1998).

3.2. The fit to astrometrical data

The goal of the fit is to optimize the parameters of the model (such as the initial conditions and the non-gravitational parameters divided by the mass of the nucleus) in order to minimize the differences between the observed positions of the comet and those delivered by our model at the same time t_n. These differences are noted (O–C)_n.

We first fit independently the non gravitational parameters and the initial position and velocity of the comet in a iterative process. Denoting ϕ_j as any Cartesian coordinates of the comet, the (O–C)_n depend on the constant of the model p_i as (14)or, in matrix notation (O–C)_n = HφX, as (15)where (16)To solve this kind of system by a least square method, we have to compute the partial derivates matrix φ. Each element of the matrix φ is obtained by numerical derivation as (17)with ϵ, a small variation on one of the model parameters p_i.

We fit the different parameters in an iterative process. In one iteration, the optimization is divided into two sub-fits because of the difference of order between magnitude of the parameters. We begin by the fit of the ratios C_i/M_c and then optimize the initial position and velocity of the comet. The errors on the parameters are determined through the variance-covariance matrix of all the parameters as (18)with χ² the sum of (O–C), N_ast the number of observations, N_par the number of parameters, and σ_ij the elements of the variance-covariance matrix.

3.3. Spin axis orientation determination

Up to now, the spin axis orientation has mainly been determined by observations of the coma morphology, which is dominated by a primary jet (Samarasinha et al. 2004). The orientation of the axis seems to be stable since Deep Space 1 measurements are compatible with the previous determination. This primary jet, shown as stable in orientation and morphology by Deep Space 1 mission, seems to be closely aligned with the rotation axis (Soderblom et al. 2002). Several authors measure the location of this primary jet from (α = 214° ± 4°, δ = −5° ± 4°) by Farnham & Cochran (2002) to (α = 233°, δ = −18°) by Thomas et al. (2001).

In this investigation, we try to determine the rotational axis through the dynamics of the body. Indeed, as seen in Sect. 2.1, the non-gravitational forces are directly related to the insolation of the nucleus, hence to its axis of rotation position. In our work, we consider a fixed axis position. We explored all the space with a 15° grid. This hypothesis was based on Samarasinha et al. (2004). Moreover, as noted in Samarasinha et al. (2004), the comet 19P/Borrelly is in an unexcited spin mode, and its axis slowly precesses by 5° − 10° per century, which is much longer than a revolution period (about 6.8 years). This last consideration allows us to perform this work now. To this end, we performed a astrometrical fit reduced-χ² map with different axis orientations (Fig. 2).

This map shows two large zones where the reduced-χ² is very low (about 2.06), which correspond to the optimal axis of rotation for a prograde or a retrograde nucleus. As noted in Sect. 2.1, the forces are averaged over one rotation of the nucleus and, as our model does not contain thermal inertia, it is impossible to distinguish between the two solutions for the spin axis orientation. These zones, projected onto a sphere, are describing two caps and the bad reduced-χ² are situated on the great circle of the celestial sphere. To find the optimal solutions, we fit the poles to correspond to this great circle. We found two possible poles (α = 190°, δ = 30°) and (α = 10°, δ = −30°). Error bars were impossible to compute with our “great circle” fit on bad χ². But, as shown in Fig. 2, the axis position is within a flat χ² zone, so the peak to valley error has an order of 40°. The previous determinations are far from our central determination but still within error bars. A possible explanation of these discrepanscies is that our model is still too simple. The thermal inertia must play a large role in polar axis determination. Nevertheless, our method works, and we hope an improved thermal model can be used after the Rosetta exploration of 67P/Churyumov-Gerasimenko.

3.4. The CONGO model (COmetary Non-Gravitational Orbits)

In the following, we adopted the spin axis orientation (192°, 29°). The nucleus is divided into three latitudinal bands: two polar caps situated between the latitudes  −90°,−45° and 90°,45° and an equatorial band between  −45°,45°. The ratio C₁/M_C represents the southern cap, C₂/M_C the equatorial band, and C₃/M_C the northern cap.

3.4.1. Astrometrical data fitting

As explained in Sect. 2, we first fit astrometrical data in order to determine parameters. This fit is sufficient to produce ephemeris. We carried out the fit on the astrometrical data in equatorial J2000 Cartesian coordinates. The equation of motion used are those presented in Sect. 2.2. We performed the integration using the Radau-integretor (Everhart 1985). The results of the fit are given in Table 2.

For our best fit (see Fig. 3), we obtain an rms of 1.43′′. To compute the reduced χ², we assumed an arbitrary error of 1′′ for all the astrometric data: (19)where Δ_obs is the accuracy of the observations. We performed a similar fit with the Marsden et al. (1973) model in the same conditions as for our model. We can compare the reduced χ² of the two models: 2.20 with the Marsden et al. (1973) model and 2.06 with the CONGO model.

3.4.2. Mass determination

To get the mass from Eq. (13), we used a Monte Carlo method on  and Q_obs From this distribution (see Fig. 4), we deduce that the optimal mass of the comet is 2.21 × 10¹³   ±   0.27 × 10¹³ kg. The error on the mass may seem very small, but it should be noted that the errors on the water production rate observations just take the S/N into account and not the uncertainty due to the Haser model and any other bias. In reality, we expect an error close to 50%. This estimation of the nucleus mass corresponds to a density of about 400 kg   m^-3. This result is quite comparable to the estimate 150   kg   m^-3 ≤ ρ ≤ 450   kg   m^-3 obtained by Davidsson & Gutiérrez (2004) (considering a nucleus volume equal to 5.5 × 10¹⁰  m³, see Lamy et al. 1998) and supports the hypothesis of a low-density object (see Table 3).

3.5. Discussion of results

To analyse our results, we can draw our theoretical water production rate curve on data. We can see that the fit is mainly made on the big amount of data coming from Combi et al. (2011) (Fig. 5). Nevertheless, there is a systematic offset by a factor of 2.5 between Combi et al. (2011) data and the other measurements (Bockelée-Morvan et al. 2004; Fink 2009; Lecacheux et al. 2003; Schleicher et al. 2003). This was also viewed on the comets Hale-Bopp (Combi et al. 2000) and Hyakutake (Combi et al. 2005) data. As this offset is quite constant, we decided to compute two separate fits, one on the Combi et al. (2011) data and the other on all the other data.

The Combi et al. (2011) data are based on direct observation of H. The problem is that H come from H₂O but also from many other molecules. This is taken into account in the model, but we can say that it can be considered as the maximum value of H₂O production rate.The mass and the density found with this set of data are 33  × 10¹² kg and 600 kg m^-3. In another band, the estimations with OH are indirect measurements of the H₂O production rate, and they can be considered as the minimum estimation. The mass and the density found with these data are 15  × 10¹² kg and 270 kg m^-3.

We computed the error on the mass, with the astrometrical fit and the errors given by authors of the H₂O production rate. We found a 15% error, which seems rather optimistic. In fact, the errors on the mass estimation are mainly due to H₂O production rates. We note a factor of 2.5 difference between the method measuring H and the others measuring OH.

3.5.1. Active zones

Knowing the value of the mass and the C_i/M_c parameters, we can now deduce the fraction of active area on the three bands of the considered nucleus. As for the mass and the density determination, we can calcule the two extreme values of the active surface interval by considering the Combi et al. (2011) observations or not. These results are reported in Table 4. The total active area percentage is between 14% and 31%. This estimation is compatible with the previous determination (8–18%) of Davidsson & Gutiérrez (2004).

We can note that the southern hemisphere is less active than the northern one. This result corresponds to the Deep Space 1 obervations of the northern hemisphere and to the post-perihelion observations by Farnham & Cochran (2002) and by Schleicher et al. (2003) when the southern hemisphere was in the direction of the Sun (Davidsson & Gutiérrez 2004).

Figure 5 shows the gas production rate compared to observations. We note that the peak of the curve is shifted by 20 days before perihelion. This is due to the ellipsoidal shape and the specific axis of rotation of the nucleus (Davidsson & Gutiérrez 2004). Indeed, before perihelion, the northern hemisphere is the most insolated part of the comet. Considering the locations of active areas on the nucleus presented above, it is easy to understand that the maximum outgassing occurs before perihelion.

4. Conclusions

This work shows that it is possible to use a simple realistic model based on the physics of H₂O outgassing for a comet ephemeris. The model will evolve in the future with a better thermal model. But even now, it can produce more accurate ephemerides than previous models (Marsden et al. 1973).

Moreover, just with astrometrical data, it is possible to estimate the spin axis orientation. This determination is still raw, mainly because of our simplistic thermal model. But, it can be applied to any comet, even a faint one.

The most important product of our model is determining of the comet nucleus mass by fitting the H₂O production rate. Using Deep Space 1 volume determination of the nucleus of 19P/Borrelly, we evaluated its density to be between 270 and 600 kg m^-3. It confirms the hypothesis of a low-density, porous object. Most of the uncertainty comes from the H₂O production rate estimation. An additional error is due to the rather arbitrary choice of the η parameter that governs the gas expansion velocity in Eq. (10).

Our estimation is compatible with previous determinations, especially with Davidsson & Gutiérrez (2004) who are using a more sophisticated thermic model.

Up to now, it has not been possible to measure cometary masses from spacecraft flybys, because these masses are too low to significantly affect the spacecraft trajectory. It is anticipated, however, that this will be possible with the Rosetta spacecraft orbiting 67P/Churyumov-Gerasimenko (Pätzold et al. 2007). It will be a good test for our model.

We can also measure the fraction of active surfaces. Our model is quite robust since it uses the same number of parameters as do Marsden et al. (1973), but it is closer to the physical reality. We can apply this method to all comets for which the size is known. If the size is unknown, we can just compute ephemerides of the comet.

Our model is a first attempt to bring together celestial mechanics and a realistic comet nucleus model. It can be used right now to compute ephemerides. Its physical parameters still need to be improved. In this resspect, the future in situ exploration of 67P/Churyumov-Gerasimenko by Rosetta will be helpful. Infrared sounding of the nucleus (Coradini et al. 2007) will lead to an accurate thermal model, and the radio spectroscopic observations (Gulkis et al. 2007) will directly measure the gas jet velocities.

The next step will also be to take the moving axis of rotation into account. With this, we really hope to extend the domain of validity of the orbits, and perhaps to have a better idea of the past of comets.

References

All Tables

All Figures