A&A 399, 789-794 (2003)
DOI: 10.1051/0004-6361:20021832

The dust tail of Comet C/1999 T1 McNaught-Hartley

F. Moreno¹ - L. M. Lara¹ - J. J. López-Moreno¹ - O. Muñoz¹ - A. Molina^1,2 - J. Licandro³

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

Abstract
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 $Af\rho$ 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

1 Introduction

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).

2 Data analysis

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

The images were bias, dark, and flat corrected, and calibrated in flux. The spatial scale was 0.4325 arcsec/pixel, which corresponds to a resolution of about 408 km/pixel and 406 km/pixel on 2001 January 26 and 2001 February 5, respectively. The images were rotated to the (M,N)coordinate system (M being the prolonged radius vector). Each image was resampled as to have a total of $100\times70$ values of the surface light intensity. This procedure has the advantage of reducing the noise considerably, as a pixel in the resampled images corresponds to a square of $10\times10$ pixels in the original images. The processed images are shown in Fig. 1. While the image of February 5 is very clean regarding the presence of field stars, the image of January 26 is, unfortunately, heavily spotted with the presence of many field stars in the frame. These stars distort the cometary isophote field and make difficult the process of retrieving the modeled images. We will return to this point later.

3 The model

The theory of the Monte Carlo approach to retrieve the reconstructed images was developed by Fulle (1989), and will not be repeated here. For clarity, we will follow the same notation that appears in Fulle (1989).

As in Fulle (1989), we consider a sample of $N_{t} N_\mu N_{s}$ dust grains, where N_t is the number of samples in the time interval of dust ejection, $N_\mu$ is the number of samples in size and N_s is the number of grains of a given size uniformly distributed on a dust shell ejected at time t. We set $N_{t}= N_\mu =N_{s} =$ 200. The grain velocity is given by:

$$\vec{V}_{d} = \vec{V}_{\rm c}+v(t)(1-\mu)^{1/k}\vec{u}_q ~~~~~~~~~ q = 1... N_{s},$$ (1)

where $\vec{V}_{\rm c}$ is the comet nucleus velocity vector, $\vec{u}$ is the unit vector along the ejection direction and v(t) is the time-dependent ejection velocity, a function to be determined by a trial-and-error procedure. The parameter $(1-\mu)$ is defined as the ratio of the radiation force to that of gravity, and it is given for spherical particles by (Finson & Probstein 1968a; Fulle 1989):

$$(1-\mu) = C_{\rm pr}Q_{\rm pr}(\rho_{d} d)^{-1},$$ (2)

where $C_{\rm pr} = 1.19\times 10^{-4}$ g cm^-2, $Q_{\rm pr}$ is the scattering efficiency for radiation pressure ( $Q_{\rm pr}\sim 1$ for large absorbing grains, see Burns et al. 1979), and $\rho_{d}$ is the mass density of the spherical grain of diameter d. The density will be assumed to be $\rho_{d} = 1$ g cm^-3. The dependence on size of the grain velocity is parameterized by the exponent 1/k. Models of the inner coma (Probstein 1968) give k=4, while dust tail analysis by Fulle gives a range of possible k depending on the comet analyzed, ranging from k=2 to k=6. In general, dust tail analysis seems to be quite insensitive to the assumptions regarding the size dependence of the dust velocity (e.g. Fulle 1989; Fulle et al. 1992). Another important parameter of the model is the degree of anisotropy of the ejection velocities. The anisotropy parameter $\omega$ of the ejection velocity is defined as the half width of the cone with axis toward the Sun, into which the particles are emitted. Cremonese & Fulle (1990) and Fulle et al. (1993) have found an important correlation between comet age and degree of anisotropy of the ejection velocities. Thus, dynamically new comets in the Oort sense, such as comets C/1987 VII Wilson or C/1990 V Austin, tend to show isotropic dust ejections, the opposite behavior to that shown by comet C/1988 V Liller, a long-period but dynamically old comet in which the observed dust tail is found to be inconsistent with isotropic dust ejections (Fulle et al. 1992). Comet C/1999 T1 is also a long-period comet but it is not a new comet in the Oort sense. Kenji Muraoka estimated the original semi-major axis of its orbit as a=897.3 AU from 357 observations, which is significantly less than 10⁴ AU, the estimated value required for a comet to be considered as dynamically new.

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 8$\times$10⁶ 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 $(1-\mu)$ value) are defined by the syndyne-synchrone network. The largest $(1-\mu)$ value at each dust grain ejection time is defined by the syndyne crossing the synchrone at the image edge, and the smallest $(1-\mu)$ value is simply the largest one divided by the number of samples in size ($N_\mu=200$). 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 $(1-\mu)$, F is the output vector, which contains the time-dependent $(1-\mu)$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, $r_{\rm h}=1.267$ AU), than on Dec. 28, 2000 ( $r_{\rm h}=1.198$ 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 $(1-\mu)$ 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.

4 Results and discussion

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 $\rho=1$ 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 ( $\omega =180^\circ$). The brightness level of the isophotes correspond to $5\times 10^{-10}$, $1\times 10^{-9}$, and $7\times 10^{-9}$. 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 - $\omega$ combinations for k=2, 4, and 6, and $\omega =45^\circ$ (highly anisotropic ejection), 90$^\circ$ (hemispherical ejection), and 180$^\circ$ (isotropic ejection), with a large number of trial velocity functions v(t). We quantified the goodness of the fits as given by $\chi^2 = \Sigma(I_{\rm obs}-I_{\rm model})/I_{\rm obs}$, where  $I_{\rm obs}$ and $I_{\rm model}$ 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, $\omega =45^\circ$, with little influence of the assumed k. The best fits for $\omega =45^\circ$ give $\chi^2=4.7\times10^{-7}$, while the best fits for isotropic or hemispherical ejections give about $\chi^2=1.5\times10^{-6}$, 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 $\chi^2$ for a given $\omega$ 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 ( $\omega =45^\circ$) for the image obtained on February 5, 2001. The brightness level of the isophotes correspond to $5\times 10^{-10}$, $1\times 10^{-9}$, and $7\times 10^{-9}$. 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 ( $\omega =45^\circ$) 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 $5\times 10^{-10}$, $1\times 10^{-9}$, and $7\times 10^{-9}$. 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 $\chi^2=5.6\times10^{-7}$, 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 $Af\rho$quantity (A'Hearn et al. 1984) from the two observations. The $Af\rho$ 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 ($\rho$). When the brightness of the coma is proportional to $\rho^{-1}$ 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 $Af\rho$ quantity. The quantity $Af\rho$ 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 $Af\rho$ 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 $2\sigma$ standard deviation.
Open with DEXTER

The dust mass loss rates were computed assuming an albedo of $Ap(\alpha)=0.04$, and should be rescaled for other values of the albedo. The dust loss rate shows a marked departure from the r_h^-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 $5.75\times10^{28}$ s^-1, as reported by Schleicher et al. (2001, IAUC 7558) on Dec. 28, 2000 ( $r_{\rm h}=1.198$ AU) and on Jan. 2, 2001, to $8.2\times10^{28}$ s^-1 on Jan. 13, 2001 ( $r_{\rm h}=1.267$ AU) (Mumma et al. 2001, IAUC 7578) or to $1.6\times10^{29}$ s^-1 also on Jan. 13, 2001, as reported by Conrad & Chaffee (2001, IAUC 7578). This gives a range of gas loss rates of $1.7\times 10^{6}$ g s^-1 to $4.8\times 10^{6}$ g s^-1. Our dust mass loss rates range from about $3\times 10^5$ to 10⁵ 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 $Af\rho$ parameter, an exactly opposite behavior to that of the dust loss rate would be found if we were considering the $Af\rho$ 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 $Af\rho$ is used as a measure of the dust production rate.

	Figure 8: Upper panel: time dependence of the quantity $Af\rho$ 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

Our results concerning the anisotropy of the grain ejection and the low values of the size distribution agree with those found generally by Fulle (e.g. Fulle 1989; Fulle 1992). The fact that the better fits correspond to highly anisotropic ejection velocities agrees with the fact that comet C/1999 T1 is not dynamically new, a result that also favors the previous results by Fulle (1992). In addition, the power index of the time dependent size distribution show generally quite small values, which also agrees with the fact that comet C/1999 T1 is not dynamically new. In fact, the power index of the time-averaged size distribution is $-3.8 \pm 0.1$. Notwithstanding this, more analysis of cometary dust tails are needed to confirm such behavior.

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 $1.7{-}2.7\times 10^{12}$ 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.

5 Conclusions

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 ( $-3.8 \pm 0.1$) also supports that conclusion.

The calculation of the time-dependent $Af\rho$ 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 $Af\rho$ quantity are almost anticorrelated at times t>+15 days since perihelion passage. Thus, the $Af\rho$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 $4{-}7\times 10^7$ 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.

Acknowledgements

This work was supported by contracts PNE-001/2000-C-01 and AYA-2001-1177.

References

 

• A'Hearn, M. F., Schleicher, D. G., Feldman, P. D., Millis, R. L., & Thompson, D. T. 1984, AJ, 89, 579 In the text NASA ADS  
• Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1 In the text NASA ADS  
• Cremonese, G., & Fulle, M. 1990, AJ, 100, 1285 In the text NASA ADS  
• Festou, M., Rickman, H., & West, R. M. 1993, A&AR, 4, 363 In the text  
• Finson, M. L., & Probstein, R. F. 1968a, ApJ, 154, 327 In the text NASA ADS  
• Finson, M. L., & Probstein, R .F. 1968b, ApJ, 154, 353 In the text NASA ADS  
• Fulle, M. 1989, A&A, 217, 283 In the text NASA ADS  
• Fulle, M. 2000, Icarus, 145, 239 In the text NASA ADS  
• Fulle, M., Cremonese, G., Jockers, K., & Rauer, H. 1992, A&A, 253, 615 In the text NASA ADS  
• Fulle, M., Bosio, S., Cremonese, G., et al. 1993, A&A, 272, 634 In the text NASA ADS  
• Fulle, M., Mikuz, H., & Bosio, S. 1997, A&A, 324, 1197 In the text NASA ADS  
• Fulle, M., Cremonese, G., & Böhm, C. 1998a, AJ, 116, 1470 In the text NASA ADS  
• Fulle, M., Mikuz, H., Nonino, M., & Bosio, S. 1998b, Icarus, 134, 235 In the text NASA ADS  
• Jockers, K. 1997, Earth, Moon & Planets, 79, 221 In the text NASA ADS  
• Lara, L. M., Licandro, J., & Tozzi, G. P. 2002, A&A, submitted In the text  
• McDonnell, J. A. M., Lamy, P. L., & Pankiewicz, G. S. 1991, in Comets in the Post-Halley Era, ed. R. L. Newburn, et al. (Kluwer Academic Publishers, London), 1043 In the text  
• Moreno, F., Muñoz, O., López-Moreno, J. J., Molina, A., & Ortiz, J. L. 2002, Icarus, 156, 474 In the text NASA ADS  
• Probstein, R. F. 1968, in Problems of Hydrodynamics and Continuum Mechanics (Soc. Industr. Appl. Math., Philadelphia), 568 In the text