A&A 399, 789-794 (2003)
F. Moreno1 - L. M. Lara1 - J. J. López-Moreno1 - O. Muñoz1 - A. Molina1,2 - J. Licandro3
1 - Instituto de Astrofísica de Andalucía, CSIC, PO Box 3004, 18080 Granada, Spain
2 - Departamento de Física Aplicada, Universidad de Granada, 18071 Granada, Spain
3 - Isaac Newton Group of Telescopes & Instituto de Astrofísica de Canarias, PO Box 321, 38700 Santa Cruz de La Palma, Tenerife, Spain
Received 9 October 2002 / Accepted 29 November 2002
A numerical method to compute cometary dust tail brightnesses has been developed. The method, analogous in many ways to the inverse Monte Carlo method proposed by Fulle (1989), has been applied to red continuum images of comet C/1999 T1 McNaught-Hartley obtained by Lara et al. (2002). The time and size dependence of the particle terminal velocities, the size distribution functions, and the dependence with time of the power index of the size distribution have been derived. The best fits to the observed isophote fields are obtained for highly anisotropic dust grain ejection velocities. The power indexes of the time-dependent dust size distributions show small values, particularly shortly after perihelion is reached, in which a minimum of -5.5 in the power index of the size distribution function is found, suggesting the apparition of small particles. Our results are in line with previous results by Fulle (1992) relating a systematic decrease of size distribution power index and dust ejection isotropy as the comet age increases. We also derived the dust mass loss rates and the quantity as a function of time, providing evidence that these two quantities do not have to be correlated if the size distribution and the dust ejection velocity are both time-dependent, in accordance with the results by Fulle (2000) for comet 46/P Wirtanen.
Key words: comets: general - comets: individual: C/1999 T1 McNaught-Hartley
The theory of the morphology of dust tails as a tool to study the dust environment of comets is based on an early mechanical theory which goes back to the nineteenth century with Bessel. The theory was subsequently improved by many others such as the Russian astronomer Feodor A. Bredikhin (see e.g. Festou et al. 1993). Using this mechanical theory, Finson & Probstein (1968a) presented a method to derive information on the size distribution, production rate, and size- and time-dependent velocities of the grains leaving the inner coma. The technique was successfully applied to comet Arend-Roland (Finson & Probstein 1968b). The method was, however, not free from simplifying assumptions, such as the isotropic dust ejection, the hypersonic approximation, and the linear expansion of the dust shells. Fulle (1989) proposed a new method based on the sampling of grains, computing rigorously their orbital motion, and avoiding the simplifying hypothesis above mentioned. The orbit of the dust grain is defined by the keplerian motion resulting from its initial velocity in the solar gravity field reduced by the pressure of solar radiation. Fulle has applied his method to many comets (e.g., Fulle 1989; Fulle et al. 1992; Fulle et al. 1998b).
We have recently envisaged a Monte Carlo method to compute energy fluxes impinging on cometary comae surrounded by dust, including the effects of polarization of light (Moreno et al. 2002). Our purpose is to have a tool to extract complementary information on the physical properties of cometary grains and the grain velocities in the coma vicinity with the interpretation of dust tail images, by building a numerical method similar to that by Fulle (1989). The model is described here, along with its application to a set of images of comet McNaught-Hartley 1999 T1 as obtained by Lara et al. (2002).
A detailed description of the observations and data reduction is given
by Lara et al. (2002). Briefly, the
two post-perihelion images
of comet C/1999 T1 analyzed here were obtained
on 2001 January 26, 06:34 UT, and 2001 February 5, 05:24 UT through a red
continuum filter on a CCD camera at the
Observatorio del Teide IAC80 telescope. The comet reached perihelion on
13.47 December 2000 (MPC 42106).
|Figure 1: Images of comet C/1999 T1 obtained on January 26 and February 5, 2001. The images are oriented to the (M,N) coordinate system and resampled as to have a resolution of 4080 km/pixel and 4060 km/pixel on Jan. 26 and Feb. 5, respectively.|
|Open with DEXTER|
As in Fulle (1989), we consider a
dust grains, where
Nt is the number of samples in the time interval of dust ejection,
is the number of samples in size and Ns is the number of grains
of a given size uniformly distributed on a dust shell ejected at
time t. We set
200. The grain velocity is given by:
Once derived the orbital elements of the dust grain according to its size and velocity, its position at a time t is derived, and then it is projected on the photographic coordinate system (M,N) (Finson & Probstein 1968a, 1968b). The Monte Carlo procedure involves 8106 particles. To obtain the modeled image densities we followed Fulle (e.g. Fulle 1989; Fulle et al. 1992; Fulle et al. 1998b). The size of the particles involved (or, equivalently, their value) are defined by the syndyne-synchrone network. The largest value at each dust grain ejection time is defined by the syndyne crossing the synchrone at the image edge, and the smallest value is simply the largest one divided by the number of samples in size (). The integration time ends at the observation time of the image and starts at a time where the contribution to the tail density is negligible The fit of the dust tail involves the inversion of the oversampled linear system AF=I, where A is the kernel matrix containing the model dust tail (i.e., the surface density of the sampling particles integrated over t and , F is the output vector, which contains the time-dependent distribution, and I is the observed surface brightness of the image in the selected (M,N)region. Since the kernel matrix A is often sparse and ill-conditioned, a regularizing condition is usually added to the system (Fulle 1989). The condition (in the least square fit sense) imposed by Fulle (1989) has the physical meaning of constraining the size distribution to be slowly time-dependent and the dust production rate to depend on the inverse square of the heliocentric distance. From the water production rates reported so far, this canonical law concerning gas production is strongly violated: the water production rate has been reported to be a factor of 2.8 higher on Jan. 13, 2001 (at heliocentric distance, AU), than on Dec. 28, 2000 ( AU) (see Schleicher et al. 2001, IAUC 7558, and Conrad & Chaffee 2001, IAUC 7578). Also, as stated by Jockers (1997), it is not clear what the regularization does to the model. For the present application, we have therefore preferred not to include any regularizing condition. From the vector F we obtained the and time-dependent distribution from which we obtained the size distribution (Finson & Probstein 1968a). The method should be applied to a sequence of calibrated images obtained during a few days apart. In this case, two images are available, although they are separated about 10 days in time. In addition, the image of January 26, 2001 is quite distorted by the presence of field stars. We have therefore preferred to apply the method only to the image of February 5 first, and give the maximum weight to the parameters derived with that single observation. After that, we modeled the image of January 26, and then, we compare the results obtained for the two observations.
Figure 2 shows the grain diameter range to which the
solutions for the images on the two dates aforementioned are related
(for an assumed density of g cm-3). Our approach was to
model first the image of February 5, 2001.
|Figure 2: The grain diameter to which the solutions of the inverse Monte Carlo procedure are related for the images of January 26, 2001 (filled circles), and February 5, 2001 (open circles).|
|Open with DEXTER|
|Figure 3: Upper panel: observed (solid contours) and reconstructed (dotted contours) isophote field for the February 5, 2001 image, for isotropic ejection velocity conditions ( ). The brightness level of the isophotes correspond to , , and . Each pixel corresponds to 4060 km. The lower panel shows the correlation between the observed and estimated fluxes.|
|Open with DEXTER|
We have tried all the possible k - combinations for k=2, 4, and 6, and (highly anisotropic ejection), 90 (hemispherical ejection), and 180 (isotropic ejection), with a large number of trial velocity functions v(t). We quantified the goodness of the fits as given by , where and are the observed and modeled images, the sum being extended to all the pixels (M,N) included in the region under analysis. The best fits correspond to those having the largest ejection anisotropy, , with little influence of the assumed k. The best fits for give , while the best fits for isotropic or hemispherical ejections give about , roughly a factor of 3 larger than those corresponding to the anisotropic case. Two examples of the reconstructed images are shown in Figs. 3 and 4. In Fig. 3 results are shown for the case of isotropic ejection velocities, and in Fig. 4 for highly anisotropic ejection velocities. It is clear how the fit improves and the noise in the reconstruction decreases in the anisotropic case. On the other hand, the size dependence of the ejection velocity is very weakly constrained, as we could find very similar values of for a given independently of the value of k, after modifying the velocity profile accordingly. The best fit for the time-dependent velocity function is shown in Fig. 5.
|Figure 4: Upper panel: observed (solid contours) and reconstructed (dotted contours) isophote field for highly anisotropic ejection conditions ( ) for the image obtained on February 5, 2001. The brightness level of the isophotes correspond to , , and . Each pixel corresponds to 4060 km. The lower panel shows the correlation between the observed and estimated fluxes.|
|Open with DEXTER|
|Figure 5: The time dependence of the dust ejection velocity as derived for the February 26, 2001 image.|
|Open with DEXTER|
|Figure 6: Upper panel: observed (solid contours) and reconstructed (dotted contours) isophote field for highly anisotropic ejection conditions ( ) for the image obtained on January 26, 2001. The regions of strong contamination due to field stars have been removed from the analysis. The brightness level of the isophotes correspond to , , and . Each pixel corresponds to 4080 km. The lower panel shows the correlation between the observed and estimated fluxes.|
|Open with DEXTER|
Once obtained the best fit for the February 5 image, we proceed to the application of the inverse Monte Carlo method to the January 26 image. Our approach was to consider, in principle, the same size and time-dependent velocity profile as for the best fitted image of February 5. With those input data, the resulted reconstruction for the January 26 image is shown in Fig. 6. The resulting fit gives , which indicates that the fit is nearly as good as that obtained for the February 5 image. This shows the consistency of the inverse Monte Carlo method in the sense that the same particle time- and size-dependent velocity function can give good fits for two images obtained 10 days apart. To ensure the robustness of the method, we should find, in addition, the same time-dependence of the power index of the derived size distribution functions. This condition may not be completely satisfied, as the grain sizes to which the solutions are related are not exactly the same for both dates (see Fig. 2), and the derived power index may well have a strong size dependence (e.g. McDonnell et al. 1991) in addition to the time dependence. Nevertheless, the agreement between both cases is good, as Fig. 7 shows. As can be seen, only the time range -10 to +58 days since perihelion is shown. For earlier times, the power index could not be estimated reliably. The wide dispersion in the power indexes may be related to the fact of not applying any regularizing constraint to the system AF=I, as that condition has the physical meaning of constraining the size distribution to be slowly time-dependent, as stated in the previous section. The power indexes show in fact a strong time dependence having a minimum close to -5.5 shortly after perihelion and increasing to values near -3 several days before and after perihelion is reached. Fulle (1989) also found a similar behavior in comets C/1962 III and C/1973 XII, in which there also appeared a population of small particles after perihelion, and in comet C/1996B2, in which there also appeared a strong drop in the size distribution power index about 20 days before perihelion (Fulle 1997), accompanied by a strong drop in the dust mass production rate. The time-dependent velocity function shows a maximum about 10 days after perihelion, and displays a similar behavior to that found by Fulle et al. (1992) for the dynamically old comet C/1988 V.
We have also calculated the time-dependent dust loss rates and the quantity (A'Hearn et al. 1984) from the two observations. The quantity, which has linear dimensions, is the product of the single scattering albedo of the dust particles (A) times the fraction of the observed field filled by the grains (f) times the radius of the diaphragm around the comet nucleus position (). When the brightness of the coma is proportional to this quantity is independent of the radius of the aperture used. These quantities were computed in our model according to the formulae given e.g. in Fulle (1998a), and are shown in Fig. 8. With the exception of the point near t=-17 days, the results derived from the two images are in good agreement in both the dust loss rates and the quantity. The quantity can be interpreted as the radius of the equivalent dust disk of unitary albedo which scatters the same light as the observed coma, and must not be directly related to a dust loss rate, as clearly demonstrated by Fulle (2000) in his analysis of comet 46/P Wirtanen. As can be seen in Fig. 8, the derived dust loss rates and the quantity are not correlated after perihelion, showing a clear anticorrelation at times t>+15 days since perihelion passage. The explanation resides in the time evolution of the size distribution function (Fulle 2000).
|Figure 7: Power index as a function of time for the size distribution functions derived for the image of January 26 (filled circles) and February 5 (open circles). The error bars represent standard deviation.|
|Open with DEXTER|
The dust mass loss rates were computed assuming an albedo of , and should be rescaled for other values of the albedo. The dust loss rate shows a marked departure from the rh-2law, a consequence of the rejection of the regularizing condition, but still it does show a maximum just before perihelion and a strong decline with heliocentric distance in both pre- and post-perihelion. Our dust loss rates can be compared with the gas loss rates derived from other observations. The water production rate has been found to increase with heliocentric distance from s-1, as reported by Schleicher et al. (2001, IAUC 7558) on Dec. 28, 2000 ( AU) and on Jan. 2, 2001, to s-1 on Jan. 13, 2001 ( AU) (Mumma et al. 2001, IAUC 7578) or to s-1 also on Jan. 13, 2001, as reported by Conrad & Chaffee (2001, IAUC 7578). This gives a range of gas loss rates of g s-1 to g s-1. Our dust mass loss rates range from about to 105 g s-1 in the same period, giving dust-to-gas ratios much smaller than 1. This result is different to the previous results by Fulle (e.g. Fulle 1997), who found dust-to-gas ratios in excess of 1 for all comets analyzed with his inverse dust tail model. Regarding the quite different time evolution of the mass loss rate and the parameter, an exactly opposite behavior to that of the dust loss rate would be found if we were considering the quantity as a measure of the dust loss rate, as its tendency is, as the water production rate, also to increase between Dec. 28, 2000, and Jan. 13, 2001. This illustrates the inconsistency in assessing variations of dust-to-gas ratios with heliocentric distance when the quantity is used as a measure of the dust production rate.
|Figure 8: Upper panel: time dependence of the quantity derived from the image obtained on January 26, 2001 (filled circles) and those derived from the image of February 5, 2001 (open circles). Lower panel: time dependence of the dust mass loss rate from the image of January 26 (filled circles) and those from the image of February 5 (open circles).|
|Open with DEXTER|
The Monte Carlo procedure allows to compute the individual orbit of each sampled grain, being possible to compute the total dust grain mass ejected that remains in bound orbits, and the minimum size of those particles contributing to the zodiacal dust cloud. Our calculations reveal that the total dust mass ejected is about g, from which 2-3% of that mass will remain in bound orbits in the Solar System. This has to be considered as a lower limit, since the solutions always refer to a limited size and time range. The minimum diameter of those grains injected into the Solar System is found to be 0.34 mm.
Numerical dust tail models can give important information concerning the dust component of comets such as the time-dependent size distribution, the ejection velocity of the grains as a function of time and its degree of spatial anisotropy, and the dust mass loss rate. We have built a Monte Carlo model on the same grounds of that described by Fulle (1989), and we have applied it to the long-period comet C/1999 T1. The model produces good fits to the tail isophotes only when the anisotropy of the ejected grains is very high, a result that seems to be common to dynamically old comets, as suggested by Fulle (1992). The low time-averaged power index of the dust size distribution found for C/1999 T1 ( ) also supports that conclusion.
The calculation of the time-dependent quantity and the dust mass loss rate from the numerical model allowed us to conclude that both parameters are not always correlated, in agreement with the analysis of Fulle (2000) of comet 46/P Wirtanen. Specifically, for comet C/1999 T1, the dust mass loss rate and the quantity are almost anticorrelated at times t>+15 days since perihelion passage. Thus, the parameter cannot be directly related to a dust loss rate, unless the size distribution and the dust ejection velocity are independent of time, physical conditions that hardly occur in practice.
The total dust mass injected into bound orbits in the Solar System is calculated from the Monte Carlo procedure to be about kg. This should be considered as a lower limit, due to the limited dust size range employed in the calculations. Also, the minimum size of those particles contributing to the zodiacal dust cloud is found to be 0.34 mm in diameter.
This work was supported by contracts PNE-001/2000-C-01 and AYA-2001-1177.