Issue 
A&A
Volume 548, December 2012



Article Number  A81  
Number of page(s)  8  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201220198  
Published online  27 November 2012 
CONGO, model of cometary nongravitational forces combining astrometric and production rate data
Application to comet 19P/Borrelly
^{1} IMCCE, Observatoire de Paris, CNRS, UPMC, Université Lille 1, 77 avenue DenfertRochereau, 75014 Paris, France
email: lucie.maquet@imcce.fr
^{2} LAM – Laboratoire d’Astrophysique de Marseille, Pôle de l’Étoile, Site de ChâteauGombert, 38 rue Frédéric JoliotCurie, 13388 Marseille Cedex 13, France
^{3} LESIA – Observatoire de Paris, CNRS, UPMC, Université ParisDiderot, 5 place Jules Janssen, 92195 Meudon, France
Received: 10 August 2012
Accepted: 12 October 2012
The gravitational orbit of a comet is affected by the sublimation of water molecules by the nucleus when the comet approaches perihelion. This outgassing triggers a nongravitational force (NGF) that significantly modifies the orbit of the comet. Up to now, modelling of this effect is mostly based on an empirical model defined in the early 70s that uses a simplified outgassing model. Attempts have been made to use advanced anisotropic thermal models both to calculate the NGF taking several observational constraints into account and to retrieve the nucleus’s mass and density, but (i) this approach is restricted to a handful of cometary nuclei that are sufficiently well known to allow this type of modelling, and (ii) the authors usually do not fit the astrometric measurements directly but rather nongravitational parameters calculated with the abovementioned empirical model. We present a new model for nongravitational forces with the aim of revisiting the problem of NGF calculation and nucleus density determination. Our model is closer to the nucleus outgassing physics with only a few free parameters. The amplitude of the perturbation depends on several parameters describing the comet activity that can be constrained by visible, infrared, and radio observations of the coma and the nucleus of the comet. It also depends on the nucleus mass, which can in turn be determined by modelling the effect of the NGF on the orbit of a comet. The method is based on the decomposition of the surface of the nucleus in several elements located at different latitudes. The contribution of each surface element to the overall NGF is fitted from the astrometric measurements, together with the density of the nucleus. This method is the only one available so far to estimate the density of cometary nuclei from groundbased observations. This method is tested on the wellknown comet 19P/Borrelly. The density found for these comet is between 150 and 600 kg m^{3}.
Key words: celestial mechanics / comets: individual: 19P/Borrelly
© 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 timedependent 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 nongravitational 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 fastrotating 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 nongravitational acceleration vector is decomposed into a radial (A_{1} g(r_{h})), a transverse (A_{2} g(r_{h})), and a normal (A_{3} g(r_{h})) component in the orbital reference frame. The three free parameters A_{1}, A_{2}, and A_{3} are solved in a least square fit of astrometrical observations, together with the six Keplerian osculating orbital elements. In this model, the parameter A_{2} 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 nongravitational force (Marsden et al. 1973). This semiempirical model soon became the standard to describe the nongravitational 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 shortperiod 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 closeup 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 nongravitational acceleration (Sekanina 1993). This led Sekanina (1993) to propose a “spotty model” of nongravitational 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 nongravitational acceleration is expressed as the sum of rotationaveraged 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 nongravitational acceleration got close to what is predicted by the standard model, in which the parameters A_{i} bestfit the astrometrical data. Sekanina (1993) noted a correlation between the sign of A_{2}, 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_{1} can be produced by an active area located close to a pole for a given orientation of the spinaxis.
The first attempt to directly introduce the rotationaveraged components of the nongravitational 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 chisquare 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 nongravitational 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 nongravitational 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.
Fig. 1 Geometrical view of a nucleus modelled with seven latitudinal strips. 
The goal at the beginning is to compute the maximum nongravitational 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 nongravitational 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 StefanBoltzmann 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 subsolar point and a point on the nucleus surface. It is defined as (3)with

(θ_{⊙},ϕ_{⊙}): the coordinates in the cometocentric frame of the subsolar 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 nongravitational 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 threedegree 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^{12} Pa and T_{0} = 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 rotationaveraged surfacic sublimation rate (over one rotation period P) of a point situated at the midlatitude θ_{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. Nongravitational accelerations
From the water production rate, we derive the nongravitational 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 nongravitational acceleration, the determination of the nucleus mass is inversely proportional to this parameter.
This expression of the nongravitational 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}: nongravitational acceleration (Sect. 2.1).
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 MonteCarlo 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 welldocumented 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.
Distribution of the astrometric observations from MPC.
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);

groundbased observations of OH from its 18cm lines, with the Nançay radio telescope (BockeléeMorvan et al. 2004) (1994 and 2001 passages);
groundbased 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).
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 nongravitational 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}.
Fig. 2 Reducedχ^{2} map in relation with the spin axis orientation. 102 points have been calculated on one hemisphere only because the problem is symmetric. There are 25 points on the celestial equator and the number decreases like cosδ because of the degeneration in the pole. The other points have been interpolated. 
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}.
Parameters of the model found after fitting to the astrometrical observations and corresponding osculating elements for comet 19P/Borrelly.
We fit the different parameters in an iterative process. In one iteration, the optimization is divided into two subfits 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 variancecovariance matrix of all the parameters as (18)with χ^{2} the sum of (O–C), N_{ast} the number of observations, N_{par} the number of parameters, and σ_{ij} the elements of the variancecovariance 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 nongravitational 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χ^{2} map with different axis orientations (Fig. 2).
This map shows two large zones where the reducedχ^{2} 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χ^{2} 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 χ^{2}. But, as shown in Fig. 2, the axis position is within a flat χ^{2} 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/ChuryumovGerasimenko.
3.4. The CONGO model (COmetary NonGravitational 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_{1}/M_{C} represents the southern cap, C_{2}/M_{C} the equatorial band, and C_{3}/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 Radauintegretor (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 χ^{2}, 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 χ^{2} of the two models: 2.20 with the Marsden et al. (1973) model and 2.06 with the CONGO model.
Fig. 3 Postfit residuals in arcseconds in right ascension and declination. 
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^{13} ± 0.27 × 10^{13} 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^{10} m^{3}, see Lamy et al. 1998) and supports the hypothesis of a lowdensity object (see Table 3).
Fig. 4 Distribution of the mass estimations from the MonteCarlo simulation. 
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éeMorvan et al. 2004; Fink 2009; Lecacheux et al. 2003; Schleicher et al. 2003). This was also viewed on the comets HaleBopp (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_{2}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_{2}O production rate.The mass and the density found with this set of data are 33 × 10^{12} kg and 600 kg m^{3}. In another band, the estimations with OH are indirect measurements of the H_{2}O production rate, and they can be considered as the minimum estimation. The mass and the density found with these data are 15 × 10^{12} 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_{2}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_{2}O production rates. We note a factor of 2.5 difference between the method measuring H and the others measuring OH.
Fig. 5 The gas prodution rates obtained for ellipsoidal nucleus with a spin axis orientation (α,δ) = (192°,29°) compared with observations for the 2001 perihelion passage. 
Masses and densities (for a volume of 55 km^{3}) estimated in other studies compared to the present work.
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 postperihelion 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).
Percentage of active areas in the three bands.
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_{2}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_{2}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 lowdensity, porous object. Most of the uncertainty comes from the H_{2}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/ChuryumovGerasimenko (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/ChuryumovGerasimenko 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
 A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2005, Science, 310, 258 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396 [NASA ADS] [CrossRef] [Google Scholar]
 Bessel, F. W. 1836, Astron. Nachr., 13, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Beutler, G. 2005, Methods of celestial mechanics, Vol. I: Physical, mathematical, and numerical principles (Berlin: Springer) [Google Scholar]
 BockeléeMorvan, D., Biver, N., Colom, P., et al. 2004, Icarus, 167, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Borrelly, A. 1905, CRAS, 140, 104 [Google Scholar]
 Combi, M. R., Reinard, A. A., Bertaux, J.L., Quemerais, E., & Mäkinen, T. 2000, Icarus, 144, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Combi, M. R., Mäkinen, J. T. T., Bertaux, J.L., & Quemérais, E. 2005, Icarus, 177, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Combi, M. R., Lee, Y., Patel, T. S., et al. 2011, AJ, 141, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Coradini, A., Capaccioni, F., Drossart, P., et al. 2007, Space Sci. Rev., 128, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Crifo, J. F. 1987, A&A, 187, 438 [NASA ADS] [Google Scholar]
 Crifo, J. F., Rodionov, A. V., & BockeleeMorvan, D. 1999, Icarus, 138, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Davidsson, B. J. R., & Gutiérrez, P. J. 2004, Icarus, 168, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Encke, J. F. 1826, Berliner Astronomische jahrbuch fur das jahr, 1826, 124 [Google Scholar]
 Everhart, E. 1985, in Dynamics of comets: Their Origin and Evolution, 185 [Google Scholar]
 Fanale, F. P., & Salvail, J. R. 1984, Icarus, 60, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Farnham, T. L., & Cochran, A. L. 2002, Icarus, 160, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Fink, U. 2009, Icarus, 201, 311 [Google Scholar]
 Froeschle, C., & Rickman, H. 1986, A&A, 170, 145 [NASA ADS] [Google Scholar]
 Groussin, O., A’Hearn, M. F., Li, J.Y., et al. 2007, Icarus, 191, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Gulkis, S., Frerking, M., Crovisier, J., et al. 2007, Space Sci. Rev., 128, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Huebner, W. F., Keady, J. J., & Lyon, S. P. 1992, Ap&SS, 195, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Julian, W. H., Samarasinha, N. H., & Belton, M. J. S. 2000, Icarus, 144, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, H. U., Britt, D., Buratti, B. J., & Thomas, N. 2004, in Comets II, eds. M. C. Festou, H. U. Keller, & H. A. Weaver, 211 [Google Scholar]
 Lamy, P. L., Toth, I., & Weaver, H. A. 1998, A&A, 337, 945 [NASA ADS] [Google Scholar]
 Lebofsky, L. A., & Spencer, J. R. 1989, in Asteroids II, eds. R. P. Binzel, T. Gehrels, & M. S. Matthews, 128 [Google Scholar]
 Lecacheux, A., Biver, N., Crovisier, J., et al. 2003, A&A, 402, L55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maquet, L. 2012, Ph.D. Thesis, Institut de Mécanique Céleste et de Calculs d’Éphémérides [Google Scholar]
 Marsden, B. G. 1969, AJ, 74, 720 [Google Scholar]
 Marsden, B. G., Sekanina, Z., & Yeomans, D. K. 1973, AJ, 78, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Pätzold, M., Häusler, B., Aksnes, K., et al. 2007, Space Sci. Rev., 128, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H., & Froeschle, C. 1983a, in Cometary Exploration, ed.T. I. Gombosi, 3, 109 [Google Scholar]
 Rickman, H., & Froeschle, C. 1983b, in Cometary Exploration, ed. T. I. Gombosi, 1, 75 [Google Scholar]
 Rickman, H., Kamel, L., Festou, M. C., & Froeschle, C. 1987, in Diversity and Similarity of Comets, eds. E. J. Rolfe, & B. Battrick, ESA SP, 278, 471 [Google Scholar]
 Samarasinha, N. H., Mueller, B. E. A., Belton, M. J. S., & Jorda, L. 2004, Comets II, eds. Festou, M. C., Keller, H. U., & Weaver, H. A., 281 [Google Scholar]
 Schleicher, D. G., Millis, R. L., & Birch, P. V. 1998, Icarus, 132, 397 [NASA ADS] [CrossRef] [Google Scholar]
 Schleicher, D. G., Woodney, L. M., & Millis, R. L. 2003, Icarus, 162, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Sekanina, Z. 1993, AJ, 105, 702 [NASA ADS] [CrossRef] [Google Scholar]
 Sekanina, Z., Brownlee, D. E., Economou, T. E., Tuzzolino, A. J., & Green, S. F. 2004, Science, 304, 1769 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sitarski, G. 1990, Acta Astron., 40, 405 [NASA ADS] [Google Scholar]
 Soderblom, L. A., Becker, T. L., Bennett, G., et al. 2002, Science, 296, 1087 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sosa, A., & Fernández, J. A. 2009, MNRAS, 393, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, N., A’Hearn, M. F., Boice, D. C., et al. 2001, in BAAS 33, AAS/Division for Planetary Sciences Meeting Abstracts, 1074 [Google Scholar]
 Washburn, E. 1928, International Critical Tables, Vol. III (New York: McGrawHill) [Google Scholar]
 Yeomans, D. K., & Chodas, P. W. 1989, AJ, 98, 1083 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Parameters of the model found after fitting to the astrometrical observations and corresponding osculating elements for comet 19P/Borrelly.
Masses and densities (for a volume of 55 km^{3}) estimated in other studies compared to the present work.
All Figures
Fig. 1 Geometrical view of a nucleus modelled with seven latitudinal strips. 

In the text 
Fig. 2 Reducedχ^{2} map in relation with the spin axis orientation. 102 points have been calculated on one hemisphere only because the problem is symmetric. There are 25 points on the celestial equator and the number decreases like cosδ because of the degeneration in the pole. The other points have been interpolated. 

In the text 
Fig. 3 Postfit residuals in arcseconds in right ascension and declination. 

In the text 
Fig. 4 Distribution of the mass estimations from the MonteCarlo simulation. 

In the text 
Fig. 5 The gas prodution rates obtained for ellipsoidal nucleus with a spin axis orientation (α,δ) = (192°,29°) compared with observations for the 2001 perihelion passage. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.