- Email: [email protected]
The Death of Comet Tabur 1996 Q1: The Tail without the Comet
ICARUS
134, 235–248 (1998) IS985943
ARTICLE NO.
The Death of Comet Tabur 1996 Q1: The Tail without the Comet M. Fulle Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy E-mail: [email protected]
H. Mikuzˇ Department of Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
and M. Nonino and S. Bosio Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy Received August 26, 1997; revised January 23, 1998
After a normal brightness increase, Comet Tabur 1996 Q1 showed a remarkable photometric behavior by rapidly fading in late October 1996. In this paper we analyze three CCD images of the remnant dust tail observed during the fading of the comet around perihelion and model them by means of the inverse dust tail model (M. Fulle, 1989, Astron. Astrophys. 217, 283–297). Assuming hemispherical sunward dust emission from the nucleus, satisfactory fits of the observed tail brightness distribution, turning axis and temporal fading allow us to conclude that only dust was observed, and contamination by gas and/or ions in the images is negligible. The model results include the temporal variation of the dust ejection velocity, the size distribution and dust mass loss rate. These values show a strong correlation during fading with strong drops consistent with the comet’s deactivation. In particular, the slow increase of the dust mass loss rate in September and its low absolute values allow us to exclude outbursts preceding fading and to exclude that the disappearance was due to a complete nucleus disruption. In this case, the nucleus mean radius should have been no more than 350 m (for a nucleus bulk density of 100 kg m23), which seems inconsistent with the observed water loss rate. A probable explanation of the comet fading is that the comet nucleus deactivation was due either to seasonal effects, putting all active areas in permanently night sides, or to the complete end of the whole nucleus surface activity (possibly due either to nucleus mantling or to the end of the ice reservoirs). 1998 Academic Press Key Words: comets; C/1996 Q1; dust tails.
1. INTRODUCTION
Comet Tabur 1996 Q1 was discovered by Vello Tabur on 19 August 1996, at 1.5 AU from the Sun and 1.3 AU
from the Earth (Green 1996a). Its parabolic orbit (with a perihelion at 0.84 AU from the Sun occurred on 3.5 November 1996; Marsden 1996a) was recognized to be very similar to that of long period Comet Liller 1988V (Jahn 1996). The C/1988 V orbital period of 2900 years makes probable that the splitting between the two objects occurred close to the preceding perihelion passage, despite their actual separation of 9 years (Marsden 1996b). During the 2 months following the discovery, the comet displayed a normal behavior, with a brightness increase up to mag 5 and the detection of X-rays from its coma (Dennerl et al. 1996). During the first half of October, the comet was characterized by Q(HCN) 5 2.5 1025 s21 (Womack and Suswal 1996) and by Q(OH) 5 4 1028 s21 (Crovisier, pers. commun.). However, in the second half of October, a strong drop of the brightness was observed: Q(OH) dropped below 1028 s21 and the coma magnitude dropped to 9.5 on 4 November, 1996, and to 10.5 on 10th (Green 1996b). Marsden (1996c) concluded that in December the comet was virtually lost. The morphology of the comet showed a similar behavior, with a strong change of the comet shape after the first half of October. During the 2 months following the discovery, the comet displayed a bright coma with a well developed ion tail, which, however, disappeared at the end of October. Then, the coma became diffuse, without any central condensation, and became typically spindle-shaped, with decreasing surface brightness, and with the axis turning with respect to the Sun direction. The color and the shape evolution of this feature pointed out that it was mainly composed of dust. Thus, it can be interpreted as a remnant dust tail, with the difference, with respect to an usual dust tail, that
235 0019-1035/98 $25.00 Copyright 1998 by Academic Press All rights of reproduction in any form reserved.
236
FULLE ET AL.
it was lacking the dust coma. Therefore, a lot of interesting questions about this feature are raised: is its brightness and shape evolution in time consistent with the dynamics of a dust tail? Is it possible to obtain a remnant dust tail with negligible dust production from the parent nucleus? Is the observed remnant dust tail built up by the continuous past activity of the parent nucleus, or is it built up by all the dust released by the possible disintegration of the whole nucleus? Is any dust tail model able to disentangle between these two extreme possibilities, thus inferring information about the evolution of the nucleus of C/1996 Q1? This paper is an attempt to answer to all these questions. In this paper we discuss the analysis of a homogeneous set of CCD images of the dust tail developed by C/1996 Q1 at the end of October. In order to infer the total dust mass in the dust tail, i.e., to infer whether it is consistent with the possible nucleus mass and also to infer the dust to gas ratio of this unusual object, particular attention is payed to the absolute calibration of the images by means of stars contained in the same frames regarding the dust tail. Then, the reduced dust tail images are modeled by means of the inverse tail model developed by Fulle (1989) and successfully tested on several dust tails (Fulle et al. 1994 and references therein). This model allows us to obtain the most probable fit (in the least square fit sense) of all the considered images and to obtain as output the time evolution of the dust loss rate, of the dust ejection velocity, and of the size distribution. The time evolution of these quantities and the absolute levels of the dust loss rate will allow us to provide the most probable answers to the questions raised by the fading of this unusual fragment comet. 2. DATA REDUCTION
The observations were performed at the Cˇrni Vrh Observatory (Slovenia). All the images were collected with the same instrument–CCD–filter combination and cover a period of 10 days during the fading of the comet (Table I). The telescope is a Schmidt–Cassegrain reflector of 0.36 m diameter and f/6.8 focal ratio. The CCD is a thinned, backilluminated Wright detector of 353 3 540 pixel2, thus covering a field of 17 3 11 arcmin2, with a scale of 1.88 arcsec pixel21. All the images were taken throughout a Cousins R filter, in order to limit the pollution due to possible ion and gas emissions. Each CCD frame was bias, dark, and flat corrected. During the observations, we preferred to avoid to take calibration images of standard stars at different zenithal angles and decided to use the stars available on the same comet frames to calibrate the corrected sky background, which was adopted as a reference level of the surface light intensities of the tail. Due to the small covered field, no standard stars were found on the available frames. However, a literature search
TABLE I Log of Observations No.
UT
r
D
l
a
Exp
K
96Q1R4 96Q1R5 96Q1R6
Oct 31.748 Nov 2.729 Nov 9.707
0.842 0.840 0.848
0.797 0.840 0.992
8.7 10.5 15.5
74.5 72.3 64.5
300 300 300
1.0 1.6 1.2
Note. No., image number; UT, observation time (start-exposure, 1996); r, D, Sun–comet and Earth–comet distances (AU), respectively; l, Earth cometocentric latitude on the comet orbital plane (degrees); a, phase angle (degrees); Exp, exposure time (s); K, correction factor of the absolute intensities.
pointed out that a secondary photometric standard star field, measured in order to monitor the dwarf nova TT Bootis, was placed a few arcmin outside image 96Q1R4 (Misselt 1996). Therefore, it was possible to find a single R plate of the digitized Palomar Sky Survey (with a scale of 1.66 arcsec pixel21) containing both the TT Bootis field and the comet field covered by image 96Q1R4. Then, we adopted the following procedure. We cut two subimages from the single digitized Palomar Sky Survey image: 96Q1C1, containing the comet field, and 96Q1C2, containing the TT Bootis field. We point out that, since these two subimages are cut from the same digitized plate, they are characterized by the same calibration zero magnitude, R0 5 RS 1 2.5 log10(IS 2 IB),
(1)
where RS is the R-magnitude of the TT Bootis standard star, and IS and IB are the integrated light intensity (in arbitrary units) computed over the same sky area AS centered on the star and on a near free star sky area, respectively. Then, we used the TT Bootis standard stars to compute R0 of image 96Q1C2, thus obtaining the same R0 for image 96Q1C1. Therefore, we were able to measure the magnitude of six stars in the comet field by means of the soobtained zero magnitude R0 . In order to avoid saturation of the brightest stars and undersampling of the faintest stars in the digitized Palomar Sky Survey, we considered stars of R-magnitudes between 14.5 and 15.5. The results of the calibration procedure are shown in Table II, where the Rmagnitudes and coordinates of the six tertiary standards are reported. At this point, it was possible to use these six tertiary standards to calibrate image 96Q1R4, containing the comet image to be calibrated. For each tertiary standard star, we measured the sky background R-magnitude arcsec22, RB 5 RS 1 2.5 log10
(IS 2 IB)AS . IB
(2)
237
THE DEATH OF C/1996 Q1
TABLE II Image Calibration No.
Star
96Q1C2 96Q1C2 96Q1C1 96Q1C1 96Q1C1 96Q1C1 96Q1C1 96Q1C1 96Q1R4 96Q1R4 96Q1R4 96Q1R4 96Q1R4 96Q1R4
TT Boo 1 TT Boo 2
RA
14h58.3m 14h58.2m 14h57.1m 14h57.1m 14h57.2m 14h57.3m 14h58.3m 14h58.2m 14h57.1m 14h57.1m 14h57.2m 14h57.3m
Dec.
RS
AS
IS
IB
R0
398579 398589 398599 398539 398539 408019 398579 398589 398599 398539 398539 408019
15.364 15.309 15.07 15.08 15.09 15.46 15.12 14.68 15.07 15.08 15.09 15.46 15.12 14.68
1100 1100 1100 1100 1100 1100 1100 1100 1400 1400 1400 1400 1400 1400
1498 1512 1553 1544 1545 1492 1520 1622 2272 2283 2245 2240 2258 2317
1349 1351 1355 1349 1351 1354 1331 1338 2190 2188 2166 2189 2189 2175
28.298 28.326 28.31 28.31 28.31 28.31 28.31 28.31
RB
19.37 19.54 19.36 19.25 19.24 19.58
Note. No., image number; Star, star (TT Boo standard stars from Misselt (1996)) of magnitude RS and coordinates RA and Dec. (2000.0); AS , sky area (arcsec2) over which the integrated intensities of the star IS and of the near sky-background IB were measured; R0 , zero magnitude; RB , surface magnitude of the sky background (R mag arcsec22). In image 96Q1C2, RS is the input and R0 is the output. In image 96Q1C1, R0 is the input and RS is the output. In image 96Q1R4, RS is the input and RB is the output.
The results shown in Table II point out that the image calibration error is not larger than 0.3 magnitudes, so that the absolute error in the model outputs should be lower than 30%, lower than the uncertainties due to the least square fit approach of the model (typically of the order of 50%). Unfortunately, images 96Q1R5 and 96Q1R6 were not covered by the same digitized Palomar Sky Survey plate, so that it was impossible to follow the same procedure to calibrate these images, too. Nevertheless, the Palomar Sky Survey allowed us to precisely find the north direction in all the three input images, so that they were precisely rotated along the sky-projected antisolar direction. In order to calibrate images 96Q1R5 and 96Q1R6, we adopted the following procedure, based on the concept that all the images must be fitted in the absolute levels by the same dust tail model. As a first step, we assumed for these two images the same RB value of image 96Q1R4. After the model fit of the tail images, it resulted that the model tail levels were systematically lower than the observed values. Therefore, to obtain the best fits of images 96Q1R5 and 96Q1R6, it was necessary to multiply the model tail brightness by a correction factor K, which reflects the difference between the RB values among the three input images. Therefore, the K value tells us which is the error introduced in the model outputs by assuming for all the input images the same RB 5 19.4 mag arcsec22. It turns out that the error is not larger than 60%. We point out that this error regards the absolute scaling of the dust mass loss rate only, which is affected by much larger uncer-
tainties due to the unknown real value of the dust albedo and dust scattering efficiency (see Discussion below). After calibration, the sky background level was subtracted from each image, which was resampled to the 30 3 30-pixel format required by the inverse tail model. The resampling was chosen in such a way to obtain about a ratio of 15 between the input and output sample numbers, usually the best ratio to obtain the faster convergence of inverse oversampled linear systems. This resampling procedure has also the advantage of strongly reducing the noise of the input image, because a pixel of the input image corresponds to a box of 3 3 9 pixels of the original image. Such input images are shown in Figs. 1–3 for selected isophote values, already corrected by the factor K provided by the model fit. In order to understand the schematic structure of the dust tail of C/1996 Q1, we plot in Fig. 4 the synchrone–syndyne network for all the three input images. This plot allows us to conclude that the way adopted to absolutely calibrate the images 96Q1R5 and 96Q1R6 is consistent with the time evolution of the dust tail. In fact, Fig. 4 points out that the number of particles inside each input image pixel is well conserved, since the dust radial motion in the tail among all the observations is negligible. Moreover, it allows us to conclude that the axis of the dust tail remains directed along the synchrones ejected between 50 and 70 days before perihelion: with respect to the Sun direction, the tail axis has the same turning rate of the synchrones. This fact both reinforces our assumption that we observe a dust tail and allows us to predict that probably at this time something may have
238
FULLE ET AL.
FIG. 1. Observed (continuous lines) and computed (dashed lines) isophotes of the dust tail image 96Q1R4 of Comet Tabur 1996 Q1. The magnitude of the innermost isophote is 21.0 mag arcsec22 in the R Cousins photometric system. The difference between isophotes is 0.75 mag (a log v(t, d) and w is the half width of the Sun-pointing dust ejection cone. The image vertical axis is the factor 2 in surface light intensity). u 5 log d antisolar direction.
occurred in the dust environment of the comet (this is the case; see Discussion below). The fading dust tail of Comet Tabur is extremely faint. From Table II and Figs. 1–3, we can conclude that the outermost isophote brightness is only 1% of the sky background. Therefore, to give physical meaning to this brightness, very careful data analysis was performed. The flat fielding procedure was very accurate. After each observing session, three flat fields were obtained, taking the image of a uniformly illuminated white screen. The flat field of each observing session was the mean of these three flats. Then, we subtracted from the raw images the zero and the dark image and then divided the result by the mean flat minus the zero image. After this flat fielding procedure, the resulting disomogeneities of the sky background were lower than the camera readout noise. The used CCD is a professional grade with EEV 02-06-1-206 sensor, thinned with 385 3 578 pixels, grade 1. The dark current is very low (0.06e pixel21 s21), resulting in a very uniform and constant zero image at the level of 51.7e pixel21. The linearity of the complete imaging system was tested on several
occasions on M67 standard photometric sequences, with a perfect linearity over the brightness range of the Tabur dust tail. Finally, the histogramic analysis of all images has confirmed that the sky background is perfectly constant within the readout noise. The data in Table II confirm this conclusion. The selected stars adopted to calibrate the 96Q1R4 image are well distributed around the tail, and the sky brightness variations (IB values) show a dispersion lower than 1%. Therefore, we can conclude that the photometric accuracy (i.e., the errors in brightness from a pixel to another in Figs. 1–3) is about 1% of the sky brightness. Obviously, the absolute photometric accuracy (measured in mag arcsec22) is much worse (RB values in Table II), due to the errors of the measurements of the star brightness. This good accuracy is fundamental in dust tail models, based on the brightness distribution over the image, and was favored by the location of the comet in the sky, with a zenithal distance not larger than 558 for image 96Q1R4 and 658 for 96Q1R6: our experience with the observing site has shown that at these elevations the sky brightness
THE DEATH OF C/1996 Q1
239
240
FULLE ET AL.
FIG. 2. Observed (continuous lines) and computed (dashed lines) isophotes of the dust tail image 96Q1R5 of Comet Tabur 1996 Q1. The magnitude of the innermost isophote is 21.4 mag arcsec22 in the R Cousins photometric system. The difference between isophotes is 0.75 mag (a log v(t, d) factor 2 in surface light intensity). u 5 and w is the half width of the Sun-pointing dust ejection cone. The image vertical axis is the log d antisolar direction.
is uniform over lengths of 10 arcmin, as is the case for the Tabur dust tail. 3. THE MODEL
In this section we briefly describe the inverse tail model adopted to fit the available images of the dust tail and to obtain the dust physical parameters most consistent with the images themselves. Details of the model can be found in Fulle (1989) and Fulle et al. (1992); it can be subdivided in two main sections: (i) the computation of the model dust tail and (ii) the fit of the observed tail by means of the model one, which provides the dust loss rate and the size distribution. The quality and the stability of the fit can be improved by changing the nonlinear free parameters of the first model section. Therefore, this is the most critical part of the approach, because it requires a trial and error procedure, which must cover all the possible values of the nonlinear parameters themselves to ensure their uniqueness. In order to limit the computation time within reasonable limits, it is therefore necessary to limit the number of
the nonlinear parameters to the lowest possible number compatible with a sufficiently realistic description of the dust properties. Several tests on tens of cometary dust tails have shown that the nonlinear parameters to which dust tail shapes and brightness distributions are most sensitive are the dust ejection anisotropy and the dust ejection velocity from the inner coma parameterized by v(t, d) 5 v(t, d0) 3 (d/d0)u, where t is the time of dust ejection from the inner coma, d is the dust diameter, d0 is a dust reference diameter (assumed here 1 mm), v(t, d0) describes the time evolution of the dust ejection, and u describes the size dependence (a power law) of the dust velocity. Therefore, the nonlinear parameters of the first model section are: (i) the time evolution of the dust ejection velocity v(t, d0), (ii) the power index u (which is constant for d . 10 em; Crifo 1991), and (iii) the dust ejection anisotropy parameter w, which is the half width of the Sun-pointing cone inside which the dust ejection is limited, thus selecting the possible directions of the dust ejection velocity vector. The dust ejection velocity vector allows us to compute
THE DEATH OF C/1996 Q1
241
242
FULLE ET AL.
FIG. 3. Observed (continuous lines) and computed (dashed lines) isophotes of the dust tail image 96Q1R6 of Comet Tabur 1996 Q1. The magnitude of the innermost isophote is 21.9 mag arcsec22 in the R Cousins photometric system. The difference between isophotes is 0.75 mag (a log v(t, d) and w is the half width of the Sun-pointing dust ejection cone. The image vertical axis is the factor 2 in surface light intensity). u 5 log d antisolar direction.
the rigorous keplerian orbit of each sample grain building up the model dust tail, computed by means of a Monte Carlo procedure involving about 107 particles. Dust dynamics depend on the ratio between solar radiation pressure and gravity forces, 1 2 e 5 CQ(rd)21, where C 5 1.19 3 1023 kg m22 is independent of the dust chemistry and physics (Burns et al. 1979). By adopting Q 5 1 (large absorbing grains) and r 5 103 kg m23, 1 2 e is converted to sizes: changes of these parameters imply a simple scaling of the model physical outputs, without changes of the dust dynamics and of the model inversion procedure. Since the 1 2 e variable is adopted, the model assumes the dust ejection velocity parametrized by v(t, 1 2 e) 5 v(t, 1 2 e0) [(1 2 e0)/(1 2 e)]u, where 1 2 e0 5 1.2 3 1023. The second model section concerns the automatic fit of the observed dust tail, which is performed by solving the oversampled linear system AF 5 I, where A is the kernel matrix containing the model dust tail, F is the output vector, and I is the data vector containing the surface brightness of all the input images. F depends on t and 1 2 e: when it is normalized versus 1 2 e, it provides the dust loss
rate and size distribution. A contains the surface density of the sampling particles in the model dust tail integrated over t and 1 2 e, so that it has dimension s m22. The integration time interval ranges back from the observation time of the first image to a starting time, which is determined as follows: the older the dust ejection time, the more diluted the ejected dust shell on the sky, the smaller its brightness contribution to the dust tail; the starting time is therefore the time before which the brightness contribution of the dust shell to the total tail brightness is negligible. The 1 2 e integration interval depends on the syndyne– synchrone network (Fig. 4) and therefore is time-dependent. Along each synchrone (i.e., for each dust ejection time), the largest 1 2 e value [namely (1 2 e)2(t)] is provided by the syndyne crossing the synchrone at the image edge. Syndynes are not circles (which would provide the same (1 2 e)2 for all times), but spirals, so that the older the synchrone, the smaller (1 2 e)2 . The smallest 1 2 e value [namely (1 2 e)1(t)] is simply given by (1 2 e)2 divided by the number of tail model size samples (usually 100). From the output F we compute the dust mass loss rate,
THE DEATH OF C/1996 Q1
243
244
FULLE ET AL.
FIG. 4. Synchrone–syndyne network for images 96Q1R4 (left), 96Q1R5 (center), and 96Q1R6 (right). The dotted line is the sky projection of the comet orbit, thus providing the real comet motion (toward the low right image corner). The syndynes (solid lines) are characterized by 1 2 e values of (from right to left, anticlockwise) 0.12, 0.04, 1.2 3 1022, 4 3 1023, 1.2 3 1023, 4 3 1024, and 1.2 3 1024, respectively. If Q 5 1 and r 5 103 kg m23, these 1 2 e values correspond to dust diameters of 0.01, 0.03, 0.1, 0.3, 1, 3, and 10 mm, respectively. The synchrones (dashed lines) are ejected at the times (from right to left, anticlockwise) 10, 30, 50, 70, 90, 110, 130, and 150 days before perihelion, respectively. The image vertical axis is the antisolar direction.
. 2fp2r 2CQ 100.4[R(2R(M,N)] M (t) 5 3Ap(a) I(M, N)
E
(12e)2(t)
(12e)1(t)
(3)
(1 2 e) F(t, 1 2 e) d(1 2 e), 21
where p is a radiant in arcsec, r is the observation Sun– comet distance in AU, R( is the Sun magnitude in the R passband, R(M, N) is the dust tail brightness in the R passband expressed as a function of the sky coordinates M and N, and Ap(a) is the albedo times the function of the phase angles a (Table I), with Ap(a) 5 Af for total diffusion over 4f of the incoming light. In Eq. (3), the surface brightness data I(M, N) used in input for the solving linear system are dimensionless, and F(t, 1 2 e) is the related solution vector F (m2 s21). Therefore, the ratio between 1020.4R(M,N) and I(M, N) is a constant dimensionless normalization factor of the input vector I. The time dependent integration limits (1 2 e)1 and (1 2 e)2 give the size range to which all the solutions are related. It follows that the dust loss rates provided by any tail model are surely lower limits of the real ones, because they are computed over a subset of the sizes of the really ejected grains. Equation (3) points out that the mass loss rate computed by means of tail models is independent of the poorly known dust bulk density r: it depends linearly on the scattering efficiency Q and inversely on the dust albedo Ap(a), assumed here Q 5 1 and Ap(a) 5 0.02 for 658 , a , 758. Changes of the albedo values imply a simple scaling of the mass loss rates values only.
The normalization of F(t, 1 2 e) versus 1 2 e allows us to obtain the 1 2 e distribution, from which we obtain the dust size distribution (Finson and Probstein 1968). Due to the adopted inverse procedure, the model provides the time evolution of the 1 2 e distribution sampled in 10 1 2 e values. In order to avoid confusion with 3D plots, we prefer to fit the obtained size distribution by means of a power law, so that the time evolution of the size distribution is described by the time evolution of its power index. This approximation, which is not required by the model, but is only adopted to show the model outputs in a less complex way, is consistent with all the available information on cometary dust size distribution. For instance, Fulle et al. (1993) have demonstrated that a power law size distribution is perfectly consistent with the fluences collected by the GIOTTO–DIDSY experiment during the Halley flyby (McDonnell et al. 1991). For completeness, we will show an example of 3D plot of the 1 2 e distribution. 4. RESULTS
The comparison between observed and model dust tails is shown in Figs. 1–3. We recall that the model fits not only the shown isophotes, but all the brightness values in the sampled pixels, so that it may happen that the fit is worse for some isophote than for others. Moreover, in order to evaluate the fit quality for a given combination of the u–w parameters, it is necessary to consider the fits
THE DEATH OF C/1996 Q1
245
FIG. 5. Dust environment of Comet Tabur 1996 Q1 for the combinations of the free parameters providing the best tail image fits. The continuous lines refer to u 5 2As and w 5 f, the dot-dashed lines to u 5 2As and w 5 f/2, the dashed lines to u 5 2Af and w 5 f/2, and the dotted lines to log v(t, d) u 5 2Ah and w 5 f/2, where u 5 and w is the half-width of the Sun-pointing dust ejection cone. For each u–w combination, we plot log d the dust velocity, the mass loss rate, the power index of the differential size distribution, and the time-dependent dust diameter interval to which all the physical quantities are related. The diameter interval derives from the model 1 2 e interval, which was converted assuming Q 5 1 and r 5 103 kg m23. For other values of these parameters, it should be scaled according to the relation 1 2 e 5 CQ(rd)21.
of all the three input images as a whole, in order to take into account the time evolution of the fading dust tail. It is apparent that all the fits characterized by w 5 f/4 overestimate the sunward extension of the tail, while the fits with w 5 f and u . 2As underestimate the left-hand side of the tail. Therefore, we conclude that the u–w combinations able to provide the best tail fits are only four, i.e., those with w 5 f/2 and any u and w 5 f and u 5 2As. In other words, the tail brightness distribution and its time evolution seem to constrain much more the dust ejection anisotropy (the best fits are those with hemispherical dust ejection) rather than the dust velocity size dependence. In Fig. 5 we plot the physical outputs related to the u–w combinations providing the four best tail fits, so that it is possible to extract from the results the time evolution common to the most significant fits. All the dust ejection velocities show a strong dependence on the Sun–comet
distance r (AU units), close to 20 r24 m s21 for grains of 1 2 e 5 1.2 3 1023 (corresponding to a diameter of 1 mm if Q 5 1 and r 5 103 kg m23), from t 5 2100 days (corresponding to r 5 2.0 AU) until t 5 225 days (corresponding to r 5 0.95 AU). After this time, the velocities show a strong decline, dropping to almost null values shortly before perihelion. The dust size distribution and mass loss rate show a similar behavior: a trend without strong jumps up to t 5 260 days and strong changes after, with a strong correlation among the solutions related to all u–w parameters. In particular, the mass loss rate begins a strong decline at t 5 230 days, corresponding to the velocity decline, dropping to values 103 times lower with respect to its peak shortly before perihelion. At the same times, the size distribution power index shows strong jumps, with a significant strong decline corresponding to the velocity and mass loss rate decline, dropping to 25,
246
FULLE ET AL.
FIG. 6. 1 2 e distribution for u 5 2As and w 5 f. The b parameter is 1 2 e 5 CQ(rd)21. Decreasing b distributions imply size distribution power indices larger than 24, whereas increasing b distributions imply size distribution power indices lower than 24. The 1 2 e distribution is set to zero where the model is unable to determine its values.
a value rarely observed in dust released by comets. The index oscillations after t 5 215 days may have no physical meaning, being related to a loss rate hundreds of times lower with respect to previous times, so that it may well be a spurious numerical output. We plot an example of the 1 2 e distribution in Fig. 6. We recall that a 1 2 e distribution constant versus 1 2 e implies a size distribution power index equal to 24 (Finson and Probstein 1968). In Fig. 5 the time-dependent size interval to which all the physical outputs are related is also shown. We recall that the model is forced (due to the syndyne–synchrone network) to consider a limited size range, a fact common (for other reasons) to all observation techniques (although this fundamental particular is often neglected in data analysis discussions). Therefore, the simplest way to explain the observed time evolution of the physical outputs seems to relate them to the decrease of the mean observed size: the dust size distribution power index might be strongly sizedependent, so that changes in time of the observed size interval would introduce spurious time dependences of the power index. The spurious decrease of the distribution power index would introduce a similar spurious drop of the mass loss rate, because the dust loss would become dominated by smaller and less massive grains. However, only a negligible fraction of the time changes of the observed outputs can be due to the decrease of the observed mean size, because most of the covered size range remains the same over all the considered times. In this case, a size dependence much larger than that observed in the physical output would be required to explain the observed time
changes. The fact that the 1 2 e values between 2 3 1024 and 6 3 1023 (corresponding to the diameters between 0.2 and 6 mm if Q 5 1 and r 5 103 kg m23) are observed during all the times considered by the model (only the sizes between 0.06 and 0.2 mm are observed after t 5 250 only, and only the size between 6 and 20 mm are observed before t 5 225 only) forces us to conclude that the observed time changes of the physical output are mostly real. The velocity increase for 2100 , t , 225 days is much steeper than that predicted by steady conditions (r21/2), in which the coma gas cooling is balanced by the solar heating (Delsemme 1982). Crifo and Rodionov (1997) have shown that it is possible to have steeper velocity increase thanks to seasonal effects. According to this scenario, the comet behavior was in some sense normal until 25 days before perihelion, after which something happened to drop the velocity to almost null values. This interpretation is supported by the photometric and gas loss rate time evolution up to the first half of October, which evidenced a steady increase of the comet activity. The available observations cannot provide information on short-term variations during such an increase. The most obvious interpretation of the sharp maximum of the dust ejection velocity at t 5 225 days, an outburst, is excluded by the time evolution of the mass loss rate, which shows a slow increase until t 5 240 days and a decrease thereafter. The striking strong correlation among the drops of the dust velocity, mass loss rate and size distribution power index for 225 , t , 210 days, seems to agree with what is expected in the case of a fast cometary activity decrease. A probable explanation of the dust velocity drop is a decrease of gas release from the nucleus surface, which implies a decrease of dust drag efficiency. Such a decrease of dust drag by gas has effects before on large grains than on small ones, so that a net effect is that the ejected size distribution becomes dominated by smaller and smaller grains: we exactly observe such a strong decrease of the size distribution power index. Due both to the decrease of dust release and to the fact that less and less massive grains are ejected, the mass loss rate shows a drop even larger than the velocity one, exactly what we observe (a factor 103 against 50). Moreover, the velocity drop should begin after the drop of the size distribution power index and, consequently, of the mass loss rate. In fact, even if the gas is able to drag less and less large grains, it drags them to the same terminal velocity. Only when the gas is unable to drag almost all the large grains, we can observe a drop of their velocity, which should be much faster than the corresponding index and mass drops. This is exactly what we observe for the velocity of 1-mm-sized grains. The strong correlation among all the physical outputs, and their good agreement with what is physically plausible, further
247
THE DEATH OF C/1996 Q1
points out the reliability of the time evolution of the physical outputs provided by the dust tail fading fit. 5. CONCLUSIONS
The results presented in the previous section allow us to select the most probable answers to the questions raised in the introduction about the fading of C/1996 Q1. Both the brightness distribution and turning axis of the observed dust tail changing in time are consistent with a dust tail model. This fact supports the assumption that most of the observed brightness is coming from dust, with probably negligible pollutions by gas and ions. The physical outputs provided by the fits of all the input images point out that the observed dust tail is a remnant dust tail, built up by the dust production that occurred before perihelion only, which becomes visible because it is no more masked by any bright coma produced during and after perihelion. Moreover, the features of the obtained physical outputs allow us to select the physical explanation most consistent with the time fading of the C/1996 Q1 dust tail. In fact, the time evolution of the dust mass loss rate allows us to answer the fundamental question raised in the Introduction: the obtained mass loss rate is inconsistent with a dust tail built up by the disintegration of the whole comet nucleus. In this case, we would wait for a strong peak of the mass loss rate in correspondence of the nucleus disintegration time. Over all, we would expect for much higher absolute values of the mass loss rate, which, on the contrary, remain at levels typical of short-period comets, with a very low dust to gas ratio (about 0.2 for an albedo times the phase function of 0.02). The total released dust mass during the perihelion approach is 1.2 3 109 kg, so that, according to the estimated dust to gas ratio, the total released mass would be 6 times larger, corresponding to a sphere of 250 m radius (for a nucleus bulk density of 100 kg m23), which can produce at most 7 3 1026 s21 water molecules from its surface, a loss rate 60 times lower than the observed one (Crovisier, pers. commun.). It seems improbable that a fractal nucleus surface also could increase the possible water loss rate of such a high factor. Even taking into account the uncertainties of the Q and Ap(a) parameters, it seems improbable that the total released dust mass would be 10 times larger, i.e., 1.2 3 1010 kg, with a dust to gas ratio equal to 2. In this case, the total released mass would be 1.5 times larger, corresponding to a sphere of 350 m radius (for a nucleus bulk density of 100 kg m23), which can produce at most 1.5 3 1027 s21 water molecules from its surface, a loss rate 30 times lower than the observed one (Crovisier, pers. commun.). We can conclude that the C/1996 Q1 nucleus mean radius was surely larger than 350 m. The strong drop following a very slow mass loss rate increase strongly favors the other interpretation: the nucleus was deactivated after
t 5 225 days, so that no new bright coma masked the normal dust tail released before perihelion. As was discussed in the previous section, the strong time correlation of the drops of all the physical outputs seems well consistent with what is expected in the case of a nucleus which is strongly decreasing the gas loss (usually called nucleus deactivation). However, the model outputs do not allow us to select the possible cause of such a deactivation. The strong dust velocity increase for 2100 , t , 225 days points out that the gas release from the nucleus surface may be dominated by seasonal effects (Crifo and Rodionov 1997), typical effects occurring in short-period comets. Also the changes in time of the size distribution power index before the drop starting at t 5 225 days might reflect changes of the active areas releasing the ejected dust. Finally, the same seasonal changes might have caused the nucleus deactivation, putting all the nucleus active areas in permanently night sides. However, the best tail fits are obtained for approximately hemispherical dust ejections, and not for strongly sunward anisotropic dust ejections, as would be expected for small active areas, which can be put in permanently night sides by the nucleus spin orientation. Hemispherical dust ejections are consistent with a nucleus activity dominated by seasons, as shown by the nucleus models of 46P/ Wirtanen, whose activity is probably dominated by seasons and characterized by an active area perhaps larger than half of the whole nucleus surface. However, the hemispherical dust ejection might point out that most of the Sunexposed nucleus was active during all the considered times, a fact in apparent contradiction with the sudden passage of all active areas in the nucleus night side. In other words, the results provided by the tail fit do not exclude a drop of the gas loss rate over all the surface independent of seasonal effects, e.g., due to dust mantling of the nucleus surface, or to exhaustion of the ice reservoir of C/1996 Q1. ACKNOWLEDGMENTS Standard stars were identified by means of the SIMBAD procedure, developed at the Observatoire Astronomique de Strasbourg. The referees M. F. A’Hearn and S. Larson suggested significant improvements to the previous version of the paper.
REFERENCES Burns, J. A., P. L. Lamy, and S. Soter 1979. Radiation forces on small particles in the Solar System. Icarus 40, 1–48. Crifo, J. F. 1991. Hydrodynamic models of the collisional coma. In Comets in the Post-Halley Era (R. L. Newburn, Jr., M. Neugebauer, and J. Rahe, Eds.), pp. 937–990. Kluwer, Dordrecht. Crifo, J. F. and A. V. Rodionov 1997. The dependence of the circumnuclear coma structure on the properties of the nucleus. I. Comparison between a homogeneous and an inhomogeneous spherical nucleus, with application to P/Wirtanen. Icarus 127, 319–353.
248
FULLE ET AL.
Delsemme, A. H. 1982. Chemical composition of cometary nuclei. In Comets (L. L. Wilkening, Ed.), pp. 85–130. Univ. of Arizona Press, Tucson. Dennerl, K., J. Englhauser, J. Truemper, and C. Lisse 1996. IAU Circ. 6495. Finson, M. L., and R. F. Probstein 1968. A theory of dust comets. I. Model and equations Astrophys. J. 154, 327–352. Fulle, M. 1989. Evaluation of cometary dust parameters from numerical simulations: Comparison with analytical approach and role of anisotropic emissions. Astron. Astrophys. 217, 283–297. Fulle, M., C. Bo¨hm, G. Mengoli, F. Muzzi, S. Orlandi, and G. Sette 1994. Current meteor production of Comet P/Swift–Tuttle. Astron. Astrophys. 292, 304–310. Fulle, M., L. Colangeli, V. Mennella, A. Rotundi, E. Bussoletti 1993. The sensitivity of the size distribution to the grain dynamics: Simulation of the dust flux measured by GIOTTO at P/Halley. Astron. Astrophys. 304, 622–630.
Fulle, M., G. Cremonese, K. Jockers, and H. Rauer 1992. The dust tail of Comet Liller 1988V. Astron. Astrophys. 253, 615–624. Green, D. W. E. 1996a. IAU Circ. 6455. Green, D. W. E. 1996b. IAU Circ. 6499, 6501, 6506, 6513. Jahn, J. 1996. IAU Circ. 6464. Marsden, B. G. 1996a. MPE Circ. 1996-R05. Marsden, B. G. 1996b. IAU Circ. 6464. Marsden, B. G. 1996c. IAU Circ. 6521. McDonnell, J. A. M., P. L. Lamy, and G. S. Pankiewicz 1991. Physical properties of cometary dust. In Comets in the Post-Halley Era (R. L. Newburn, Jr., M. Neugebauer, and J. Rahe, Eds.), pp. 937–990. Kluwer, Dordrecht. Misselt, K. A. 1996. Secondary photometric standards in selected northern dwarf-nova fields. Publ. Astron. Soc. Pacific 108, 146–165. Womack, M., and D. Suswal 1996. IAU Circ. 6485.
134, 235–248 (1998) IS985943
ARTICLE NO.
The Death of Comet Tabur 1996 Q1: The Tail without the Comet M. Fulle Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy E-mail: [email protected]
H. Mikuzˇ Department of Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
and M. Nonino and S. Bosio Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy Received August 26, 1997; revised January 23, 1998
After a normal brightness increase, Comet Tabur 1996 Q1 showed a remarkable photometric behavior by rapidly fading in late October 1996. In this paper we analyze three CCD images of the remnant dust tail observed during the fading of the comet around perihelion and model them by means of the inverse dust tail model (M. Fulle, 1989, Astron. Astrophys. 217, 283–297). Assuming hemispherical sunward dust emission from the nucleus, satisfactory fits of the observed tail brightness distribution, turning axis and temporal fading allow us to conclude that only dust was observed, and contamination by gas and/or ions in the images is negligible. The model results include the temporal variation of the dust ejection velocity, the size distribution and dust mass loss rate. These values show a strong correlation during fading with strong drops consistent with the comet’s deactivation. In particular, the slow increase of the dust mass loss rate in September and its low absolute values allow us to exclude outbursts preceding fading and to exclude that the disappearance was due to a complete nucleus disruption. In this case, the nucleus mean radius should have been no more than 350 m (for a nucleus bulk density of 100 kg m23), which seems inconsistent with the observed water loss rate. A probable explanation of the comet fading is that the comet nucleus deactivation was due either to seasonal effects, putting all active areas in permanently night sides, or to the complete end of the whole nucleus surface activity (possibly due either to nucleus mantling or to the end of the ice reservoirs). 1998 Academic Press Key Words: comets; C/1996 Q1; dust tails.
1. INTRODUCTION
Comet Tabur 1996 Q1 was discovered by Vello Tabur on 19 August 1996, at 1.5 AU from the Sun and 1.3 AU
from the Earth (Green 1996a). Its parabolic orbit (with a perihelion at 0.84 AU from the Sun occurred on 3.5 November 1996; Marsden 1996a) was recognized to be very similar to that of long period Comet Liller 1988V (Jahn 1996). The C/1988 V orbital period of 2900 years makes probable that the splitting between the two objects occurred close to the preceding perihelion passage, despite their actual separation of 9 years (Marsden 1996b). During the 2 months following the discovery, the comet displayed a normal behavior, with a brightness increase up to mag 5 and the detection of X-rays from its coma (Dennerl et al. 1996). During the first half of October, the comet was characterized by Q(HCN) 5 2.5 1025 s21 (Womack and Suswal 1996) and by Q(OH) 5 4 1028 s21 (Crovisier, pers. commun.). However, in the second half of October, a strong drop of the brightness was observed: Q(OH) dropped below 1028 s21 and the coma magnitude dropped to 9.5 on 4 November, 1996, and to 10.5 on 10th (Green 1996b). Marsden (1996c) concluded that in December the comet was virtually lost. The morphology of the comet showed a similar behavior, with a strong change of the comet shape after the first half of October. During the 2 months following the discovery, the comet displayed a bright coma with a well developed ion tail, which, however, disappeared at the end of October. Then, the coma became diffuse, without any central condensation, and became typically spindle-shaped, with decreasing surface brightness, and with the axis turning with respect to the Sun direction. The color and the shape evolution of this feature pointed out that it was mainly composed of dust. Thus, it can be interpreted as a remnant dust tail, with the difference, with respect to an usual dust tail, that
235 0019-1035/98 $25.00 Copyright 1998 by Academic Press All rights of reproduction in any form reserved.
236
FULLE ET AL.
it was lacking the dust coma. Therefore, a lot of interesting questions about this feature are raised: is its brightness and shape evolution in time consistent with the dynamics of a dust tail? Is it possible to obtain a remnant dust tail with negligible dust production from the parent nucleus? Is the observed remnant dust tail built up by the continuous past activity of the parent nucleus, or is it built up by all the dust released by the possible disintegration of the whole nucleus? Is any dust tail model able to disentangle between these two extreme possibilities, thus inferring information about the evolution of the nucleus of C/1996 Q1? This paper is an attempt to answer to all these questions. In this paper we discuss the analysis of a homogeneous set of CCD images of the dust tail developed by C/1996 Q1 at the end of October. In order to infer the total dust mass in the dust tail, i.e., to infer whether it is consistent with the possible nucleus mass and also to infer the dust to gas ratio of this unusual object, particular attention is payed to the absolute calibration of the images by means of stars contained in the same frames regarding the dust tail. Then, the reduced dust tail images are modeled by means of the inverse tail model developed by Fulle (1989) and successfully tested on several dust tails (Fulle et al. 1994 and references therein). This model allows us to obtain the most probable fit (in the least square fit sense) of all the considered images and to obtain as output the time evolution of the dust loss rate, of the dust ejection velocity, and of the size distribution. The time evolution of these quantities and the absolute levels of the dust loss rate will allow us to provide the most probable answers to the questions raised by the fading of this unusual fragment comet. 2. DATA REDUCTION
The observations were performed at the Cˇrni Vrh Observatory (Slovenia). All the images were collected with the same instrument–CCD–filter combination and cover a period of 10 days during the fading of the comet (Table I). The telescope is a Schmidt–Cassegrain reflector of 0.36 m diameter and f/6.8 focal ratio. The CCD is a thinned, backilluminated Wright detector of 353 3 540 pixel2, thus covering a field of 17 3 11 arcmin2, with a scale of 1.88 arcsec pixel21. All the images were taken throughout a Cousins R filter, in order to limit the pollution due to possible ion and gas emissions. Each CCD frame was bias, dark, and flat corrected. During the observations, we preferred to avoid to take calibration images of standard stars at different zenithal angles and decided to use the stars available on the same comet frames to calibrate the corrected sky background, which was adopted as a reference level of the surface light intensities of the tail. Due to the small covered field, no standard stars were found on the available frames. However, a literature search
TABLE I Log of Observations No.
UT
r
D
l
a
Exp
K
96Q1R4 96Q1R5 96Q1R6
Oct 31.748 Nov 2.729 Nov 9.707
0.842 0.840 0.848
0.797 0.840 0.992
8.7 10.5 15.5
74.5 72.3 64.5
300 300 300
1.0 1.6 1.2
Note. No., image number; UT, observation time (start-exposure, 1996); r, D, Sun–comet and Earth–comet distances (AU), respectively; l, Earth cometocentric latitude on the comet orbital plane (degrees); a, phase angle (degrees); Exp, exposure time (s); K, correction factor of the absolute intensities.
pointed out that a secondary photometric standard star field, measured in order to monitor the dwarf nova TT Bootis, was placed a few arcmin outside image 96Q1R4 (Misselt 1996). Therefore, it was possible to find a single R plate of the digitized Palomar Sky Survey (with a scale of 1.66 arcsec pixel21) containing both the TT Bootis field and the comet field covered by image 96Q1R4. Then, we adopted the following procedure. We cut two subimages from the single digitized Palomar Sky Survey image: 96Q1C1, containing the comet field, and 96Q1C2, containing the TT Bootis field. We point out that, since these two subimages are cut from the same digitized plate, they are characterized by the same calibration zero magnitude, R0 5 RS 1 2.5 log10(IS 2 IB),
(1)
where RS is the R-magnitude of the TT Bootis standard star, and IS and IB are the integrated light intensity (in arbitrary units) computed over the same sky area AS centered on the star and on a near free star sky area, respectively. Then, we used the TT Bootis standard stars to compute R0 of image 96Q1C2, thus obtaining the same R0 for image 96Q1C1. Therefore, we were able to measure the magnitude of six stars in the comet field by means of the soobtained zero magnitude R0 . In order to avoid saturation of the brightest stars and undersampling of the faintest stars in the digitized Palomar Sky Survey, we considered stars of R-magnitudes between 14.5 and 15.5. The results of the calibration procedure are shown in Table II, where the Rmagnitudes and coordinates of the six tertiary standards are reported. At this point, it was possible to use these six tertiary standards to calibrate image 96Q1R4, containing the comet image to be calibrated. For each tertiary standard star, we measured the sky background R-magnitude arcsec22, RB 5 RS 1 2.5 log10
(IS 2 IB)AS . IB
(2)
237
THE DEATH OF C/1996 Q1
TABLE II Image Calibration No.
Star
96Q1C2 96Q1C2 96Q1C1 96Q1C1 96Q1C1 96Q1C1 96Q1C1 96Q1C1 96Q1R4 96Q1R4 96Q1R4 96Q1R4 96Q1R4 96Q1R4
TT Boo 1 TT Boo 2
RA
14h58.3m 14h58.2m 14h57.1m 14h57.1m 14h57.2m 14h57.3m 14h58.3m 14h58.2m 14h57.1m 14h57.1m 14h57.2m 14h57.3m
Dec.
RS
AS
IS
IB
R0
398579 398589 398599 398539 398539 408019 398579 398589 398599 398539 398539 408019
15.364 15.309 15.07 15.08 15.09 15.46 15.12 14.68 15.07 15.08 15.09 15.46 15.12 14.68
1100 1100 1100 1100 1100 1100 1100 1100 1400 1400 1400 1400 1400 1400
1498 1512 1553 1544 1545 1492 1520 1622 2272 2283 2245 2240 2258 2317
1349 1351 1355 1349 1351 1354 1331 1338 2190 2188 2166 2189 2189 2175
28.298 28.326 28.31 28.31 28.31 28.31 28.31 28.31
RB
19.37 19.54 19.36 19.25 19.24 19.58
Note. No., image number; Star, star (TT Boo standard stars from Misselt (1996)) of magnitude RS and coordinates RA and Dec. (2000.0); AS , sky area (arcsec2) over which the integrated intensities of the star IS and of the near sky-background IB were measured; R0 , zero magnitude; RB , surface magnitude of the sky background (R mag arcsec22). In image 96Q1C2, RS is the input and R0 is the output. In image 96Q1C1, R0 is the input and RS is the output. In image 96Q1R4, RS is the input and RB is the output.
The results shown in Table II point out that the image calibration error is not larger than 0.3 magnitudes, so that the absolute error in the model outputs should be lower than 30%, lower than the uncertainties due to the least square fit approach of the model (typically of the order of 50%). Unfortunately, images 96Q1R5 and 96Q1R6 were not covered by the same digitized Palomar Sky Survey plate, so that it was impossible to follow the same procedure to calibrate these images, too. Nevertheless, the Palomar Sky Survey allowed us to precisely find the north direction in all the three input images, so that they were precisely rotated along the sky-projected antisolar direction. In order to calibrate images 96Q1R5 and 96Q1R6, we adopted the following procedure, based on the concept that all the images must be fitted in the absolute levels by the same dust tail model. As a first step, we assumed for these two images the same RB value of image 96Q1R4. After the model fit of the tail images, it resulted that the model tail levels were systematically lower than the observed values. Therefore, to obtain the best fits of images 96Q1R5 and 96Q1R6, it was necessary to multiply the model tail brightness by a correction factor K, which reflects the difference between the RB values among the three input images. Therefore, the K value tells us which is the error introduced in the model outputs by assuming for all the input images the same RB 5 19.4 mag arcsec22. It turns out that the error is not larger than 60%. We point out that this error regards the absolute scaling of the dust mass loss rate only, which is affected by much larger uncer-
tainties due to the unknown real value of the dust albedo and dust scattering efficiency (see Discussion below). After calibration, the sky background level was subtracted from each image, which was resampled to the 30 3 30-pixel format required by the inverse tail model. The resampling was chosen in such a way to obtain about a ratio of 15 between the input and output sample numbers, usually the best ratio to obtain the faster convergence of inverse oversampled linear systems. This resampling procedure has also the advantage of strongly reducing the noise of the input image, because a pixel of the input image corresponds to a box of 3 3 9 pixels of the original image. Such input images are shown in Figs. 1–3 for selected isophote values, already corrected by the factor K provided by the model fit. In order to understand the schematic structure of the dust tail of C/1996 Q1, we plot in Fig. 4 the synchrone–syndyne network for all the three input images. This plot allows us to conclude that the way adopted to absolutely calibrate the images 96Q1R5 and 96Q1R6 is consistent with the time evolution of the dust tail. In fact, Fig. 4 points out that the number of particles inside each input image pixel is well conserved, since the dust radial motion in the tail among all the observations is negligible. Moreover, it allows us to conclude that the axis of the dust tail remains directed along the synchrones ejected between 50 and 70 days before perihelion: with respect to the Sun direction, the tail axis has the same turning rate of the synchrones. This fact both reinforces our assumption that we observe a dust tail and allows us to predict that probably at this time something may have
238
FULLE ET AL.
FIG. 1. Observed (continuous lines) and computed (dashed lines) isophotes of the dust tail image 96Q1R4 of Comet Tabur 1996 Q1. The magnitude of the innermost isophote is 21.0 mag arcsec22 in the R Cousins photometric system. The difference between isophotes is 0.75 mag (a log v(t, d) and w is the half width of the Sun-pointing dust ejection cone. The image vertical axis is the factor 2 in surface light intensity). u 5 log d antisolar direction.
occurred in the dust environment of the comet (this is the case; see Discussion below). The fading dust tail of Comet Tabur is extremely faint. From Table II and Figs. 1–3, we can conclude that the outermost isophote brightness is only 1% of the sky background. Therefore, to give physical meaning to this brightness, very careful data analysis was performed. The flat fielding procedure was very accurate. After each observing session, three flat fields were obtained, taking the image of a uniformly illuminated white screen. The flat field of each observing session was the mean of these three flats. Then, we subtracted from the raw images the zero and the dark image and then divided the result by the mean flat minus the zero image. After this flat fielding procedure, the resulting disomogeneities of the sky background were lower than the camera readout noise. The used CCD is a professional grade with EEV 02-06-1-206 sensor, thinned with 385 3 578 pixels, grade 1. The dark current is very low (0.06e pixel21 s21), resulting in a very uniform and constant zero image at the level of 51.7e pixel21. The linearity of the complete imaging system was tested on several
occasions on M67 standard photometric sequences, with a perfect linearity over the brightness range of the Tabur dust tail. Finally, the histogramic analysis of all images has confirmed that the sky background is perfectly constant within the readout noise. The data in Table II confirm this conclusion. The selected stars adopted to calibrate the 96Q1R4 image are well distributed around the tail, and the sky brightness variations (IB values) show a dispersion lower than 1%. Therefore, we can conclude that the photometric accuracy (i.e., the errors in brightness from a pixel to another in Figs. 1–3) is about 1% of the sky brightness. Obviously, the absolute photometric accuracy (measured in mag arcsec22) is much worse (RB values in Table II), due to the errors of the measurements of the star brightness. This good accuracy is fundamental in dust tail models, based on the brightness distribution over the image, and was favored by the location of the comet in the sky, with a zenithal distance not larger than 558 for image 96Q1R4 and 658 for 96Q1R6: our experience with the observing site has shown that at these elevations the sky brightness
THE DEATH OF C/1996 Q1
239
240
FULLE ET AL.
FIG. 2. Observed (continuous lines) and computed (dashed lines) isophotes of the dust tail image 96Q1R5 of Comet Tabur 1996 Q1. The magnitude of the innermost isophote is 21.4 mag arcsec22 in the R Cousins photometric system. The difference between isophotes is 0.75 mag (a log v(t, d) factor 2 in surface light intensity). u 5 and w is the half width of the Sun-pointing dust ejection cone. The image vertical axis is the log d antisolar direction.
is uniform over lengths of 10 arcmin, as is the case for the Tabur dust tail. 3. THE MODEL
In this section we briefly describe the inverse tail model adopted to fit the available images of the dust tail and to obtain the dust physical parameters most consistent with the images themselves. Details of the model can be found in Fulle (1989) and Fulle et al. (1992); it can be subdivided in two main sections: (i) the computation of the model dust tail and (ii) the fit of the observed tail by means of the model one, which provides the dust loss rate and the size distribution. The quality and the stability of the fit can be improved by changing the nonlinear free parameters of the first model section. Therefore, this is the most critical part of the approach, because it requires a trial and error procedure, which must cover all the possible values of the nonlinear parameters themselves to ensure their uniqueness. In order to limit the computation time within reasonable limits, it is therefore necessary to limit the number of
the nonlinear parameters to the lowest possible number compatible with a sufficiently realistic description of the dust properties. Several tests on tens of cometary dust tails have shown that the nonlinear parameters to which dust tail shapes and brightness distributions are most sensitive are the dust ejection anisotropy and the dust ejection velocity from the inner coma parameterized by v(t, d) 5 v(t, d0) 3 (d/d0)u, where t is the time of dust ejection from the inner coma, d is the dust diameter, d0 is a dust reference diameter (assumed here 1 mm), v(t, d0) describes the time evolution of the dust ejection, and u describes the size dependence (a power law) of the dust velocity. Therefore, the nonlinear parameters of the first model section are: (i) the time evolution of the dust ejection velocity v(t, d0), (ii) the power index u (which is constant for d . 10 em; Crifo 1991), and (iii) the dust ejection anisotropy parameter w, which is the half width of the Sun-pointing cone inside which the dust ejection is limited, thus selecting the possible directions of the dust ejection velocity vector. The dust ejection velocity vector allows us to compute
THE DEATH OF C/1996 Q1
241
242
FULLE ET AL.
FIG. 3. Observed (continuous lines) and computed (dashed lines) isophotes of the dust tail image 96Q1R6 of Comet Tabur 1996 Q1. The magnitude of the innermost isophote is 21.9 mag arcsec22 in the R Cousins photometric system. The difference between isophotes is 0.75 mag (a log v(t, d) and w is the half width of the Sun-pointing dust ejection cone. The image vertical axis is the factor 2 in surface light intensity). u 5 log d antisolar direction.
the rigorous keplerian orbit of each sample grain building up the model dust tail, computed by means of a Monte Carlo procedure involving about 107 particles. Dust dynamics depend on the ratio between solar radiation pressure and gravity forces, 1 2 e 5 CQ(rd)21, where C 5 1.19 3 1023 kg m22 is independent of the dust chemistry and physics (Burns et al. 1979). By adopting Q 5 1 (large absorbing grains) and r 5 103 kg m23, 1 2 e is converted to sizes: changes of these parameters imply a simple scaling of the model physical outputs, without changes of the dust dynamics and of the model inversion procedure. Since the 1 2 e variable is adopted, the model assumes the dust ejection velocity parametrized by v(t, 1 2 e) 5 v(t, 1 2 e0) [(1 2 e0)/(1 2 e)]u, where 1 2 e0 5 1.2 3 1023. The second model section concerns the automatic fit of the observed dust tail, which is performed by solving the oversampled linear system AF 5 I, where A is the kernel matrix containing the model dust tail, F is the output vector, and I is the data vector containing the surface brightness of all the input images. F depends on t and 1 2 e: when it is normalized versus 1 2 e, it provides the dust loss
rate and size distribution. A contains the surface density of the sampling particles in the model dust tail integrated over t and 1 2 e, so that it has dimension s m22. The integration time interval ranges back from the observation time of the first image to a starting time, which is determined as follows: the older the dust ejection time, the more diluted the ejected dust shell on the sky, the smaller its brightness contribution to the dust tail; the starting time is therefore the time before which the brightness contribution of the dust shell to the total tail brightness is negligible. The 1 2 e integration interval depends on the syndyne– synchrone network (Fig. 4) and therefore is time-dependent. Along each synchrone (i.e., for each dust ejection time), the largest 1 2 e value [namely (1 2 e)2(t)] is provided by the syndyne crossing the synchrone at the image edge. Syndynes are not circles (which would provide the same (1 2 e)2 for all times), but spirals, so that the older the synchrone, the smaller (1 2 e)2 . The smallest 1 2 e value [namely (1 2 e)1(t)] is simply given by (1 2 e)2 divided by the number of tail model size samples (usually 100). From the output F we compute the dust mass loss rate,
THE DEATH OF C/1996 Q1
243
244
FULLE ET AL.
FIG. 4. Synchrone–syndyne network for images 96Q1R4 (left), 96Q1R5 (center), and 96Q1R6 (right). The dotted line is the sky projection of the comet orbit, thus providing the real comet motion (toward the low right image corner). The syndynes (solid lines) are characterized by 1 2 e values of (from right to left, anticlockwise) 0.12, 0.04, 1.2 3 1022, 4 3 1023, 1.2 3 1023, 4 3 1024, and 1.2 3 1024, respectively. If Q 5 1 and r 5 103 kg m23, these 1 2 e values correspond to dust diameters of 0.01, 0.03, 0.1, 0.3, 1, 3, and 10 mm, respectively. The synchrones (dashed lines) are ejected at the times (from right to left, anticlockwise) 10, 30, 50, 70, 90, 110, 130, and 150 days before perihelion, respectively. The image vertical axis is the antisolar direction.
. 2fp2r 2CQ 100.4[R(2R(M,N)] M (t) 5 3Ap(a) I(M, N)
E
(12e)2(t)
(12e)1(t)
(3)
(1 2 e) F(t, 1 2 e) d(1 2 e), 21
where p is a radiant in arcsec, r is the observation Sun– comet distance in AU, R( is the Sun magnitude in the R passband, R(M, N) is the dust tail brightness in the R passband expressed as a function of the sky coordinates M and N, and Ap(a) is the albedo times the function of the phase angles a (Table I), with Ap(a) 5 Af for total diffusion over 4f of the incoming light. In Eq. (3), the surface brightness data I(M, N) used in input for the solving linear system are dimensionless, and F(t, 1 2 e) is the related solution vector F (m2 s21). Therefore, the ratio between 1020.4R(M,N) and I(M, N) is a constant dimensionless normalization factor of the input vector I. The time dependent integration limits (1 2 e)1 and (1 2 e)2 give the size range to which all the solutions are related. It follows that the dust loss rates provided by any tail model are surely lower limits of the real ones, because they are computed over a subset of the sizes of the really ejected grains. Equation (3) points out that the mass loss rate computed by means of tail models is independent of the poorly known dust bulk density r: it depends linearly on the scattering efficiency Q and inversely on the dust albedo Ap(a), assumed here Q 5 1 and Ap(a) 5 0.02 for 658 , a , 758. Changes of the albedo values imply a simple scaling of the mass loss rates values only.
The normalization of F(t, 1 2 e) versus 1 2 e allows us to obtain the 1 2 e distribution, from which we obtain the dust size distribution (Finson and Probstein 1968). Due to the adopted inverse procedure, the model provides the time evolution of the 1 2 e distribution sampled in 10 1 2 e values. In order to avoid confusion with 3D plots, we prefer to fit the obtained size distribution by means of a power law, so that the time evolution of the size distribution is described by the time evolution of its power index. This approximation, which is not required by the model, but is only adopted to show the model outputs in a less complex way, is consistent with all the available information on cometary dust size distribution. For instance, Fulle et al. (1993) have demonstrated that a power law size distribution is perfectly consistent with the fluences collected by the GIOTTO–DIDSY experiment during the Halley flyby (McDonnell et al. 1991). For completeness, we will show an example of 3D plot of the 1 2 e distribution. 4. RESULTS
The comparison between observed and model dust tails is shown in Figs. 1–3. We recall that the model fits not only the shown isophotes, but all the brightness values in the sampled pixels, so that it may happen that the fit is worse for some isophote than for others. Moreover, in order to evaluate the fit quality for a given combination of the u–w parameters, it is necessary to consider the fits
THE DEATH OF C/1996 Q1
245
FIG. 5. Dust environment of Comet Tabur 1996 Q1 for the combinations of the free parameters providing the best tail image fits. The continuous lines refer to u 5 2As and w 5 f, the dot-dashed lines to u 5 2As and w 5 f/2, the dashed lines to u 5 2Af and w 5 f/2, and the dotted lines to log v(t, d) u 5 2Ah and w 5 f/2, where u 5 and w is the half-width of the Sun-pointing dust ejection cone. For each u–w combination, we plot log d the dust velocity, the mass loss rate, the power index of the differential size distribution, and the time-dependent dust diameter interval to which all the physical quantities are related. The diameter interval derives from the model 1 2 e interval, which was converted assuming Q 5 1 and r 5 103 kg m23. For other values of these parameters, it should be scaled according to the relation 1 2 e 5 CQ(rd)21.
of all the three input images as a whole, in order to take into account the time evolution of the fading dust tail. It is apparent that all the fits characterized by w 5 f/4 overestimate the sunward extension of the tail, while the fits with w 5 f and u . 2As underestimate the left-hand side of the tail. Therefore, we conclude that the u–w combinations able to provide the best tail fits are only four, i.e., those with w 5 f/2 and any u and w 5 f and u 5 2As. In other words, the tail brightness distribution and its time evolution seem to constrain much more the dust ejection anisotropy (the best fits are those with hemispherical dust ejection) rather than the dust velocity size dependence. In Fig. 5 we plot the physical outputs related to the u–w combinations providing the four best tail fits, so that it is possible to extract from the results the time evolution common to the most significant fits. All the dust ejection velocities show a strong dependence on the Sun–comet
distance r (AU units), close to 20 r24 m s21 for grains of 1 2 e 5 1.2 3 1023 (corresponding to a diameter of 1 mm if Q 5 1 and r 5 103 kg m23), from t 5 2100 days (corresponding to r 5 2.0 AU) until t 5 225 days (corresponding to r 5 0.95 AU). After this time, the velocities show a strong decline, dropping to almost null values shortly before perihelion. The dust size distribution and mass loss rate show a similar behavior: a trend without strong jumps up to t 5 260 days and strong changes after, with a strong correlation among the solutions related to all u–w parameters. In particular, the mass loss rate begins a strong decline at t 5 230 days, corresponding to the velocity decline, dropping to values 103 times lower with respect to its peak shortly before perihelion. At the same times, the size distribution power index shows strong jumps, with a significant strong decline corresponding to the velocity and mass loss rate decline, dropping to 25,
246
FULLE ET AL.
FIG. 6. 1 2 e distribution for u 5 2As and w 5 f. The b parameter is 1 2 e 5 CQ(rd)21. Decreasing b distributions imply size distribution power indices larger than 24, whereas increasing b distributions imply size distribution power indices lower than 24. The 1 2 e distribution is set to zero where the model is unable to determine its values.
a value rarely observed in dust released by comets. The index oscillations after t 5 215 days may have no physical meaning, being related to a loss rate hundreds of times lower with respect to previous times, so that it may well be a spurious numerical output. We plot an example of the 1 2 e distribution in Fig. 6. We recall that a 1 2 e distribution constant versus 1 2 e implies a size distribution power index equal to 24 (Finson and Probstein 1968). In Fig. 5 the time-dependent size interval to which all the physical outputs are related is also shown. We recall that the model is forced (due to the syndyne–synchrone network) to consider a limited size range, a fact common (for other reasons) to all observation techniques (although this fundamental particular is often neglected in data analysis discussions). Therefore, the simplest way to explain the observed time evolution of the physical outputs seems to relate them to the decrease of the mean observed size: the dust size distribution power index might be strongly sizedependent, so that changes in time of the observed size interval would introduce spurious time dependences of the power index. The spurious decrease of the distribution power index would introduce a similar spurious drop of the mass loss rate, because the dust loss would become dominated by smaller and less massive grains. However, only a negligible fraction of the time changes of the observed outputs can be due to the decrease of the observed mean size, because most of the covered size range remains the same over all the considered times. In this case, a size dependence much larger than that observed in the physical output would be required to explain the observed time
changes. The fact that the 1 2 e values between 2 3 1024 and 6 3 1023 (corresponding to the diameters between 0.2 and 6 mm if Q 5 1 and r 5 103 kg m23) are observed during all the times considered by the model (only the sizes between 0.06 and 0.2 mm are observed after t 5 250 only, and only the size between 6 and 20 mm are observed before t 5 225 only) forces us to conclude that the observed time changes of the physical output are mostly real. The velocity increase for 2100 , t , 225 days is much steeper than that predicted by steady conditions (r21/2), in which the coma gas cooling is balanced by the solar heating (Delsemme 1982). Crifo and Rodionov (1997) have shown that it is possible to have steeper velocity increase thanks to seasonal effects. According to this scenario, the comet behavior was in some sense normal until 25 days before perihelion, after which something happened to drop the velocity to almost null values. This interpretation is supported by the photometric and gas loss rate time evolution up to the first half of October, which evidenced a steady increase of the comet activity. The available observations cannot provide information on short-term variations during such an increase. The most obvious interpretation of the sharp maximum of the dust ejection velocity at t 5 225 days, an outburst, is excluded by the time evolution of the mass loss rate, which shows a slow increase until t 5 240 days and a decrease thereafter. The striking strong correlation among the drops of the dust velocity, mass loss rate and size distribution power index for 225 , t , 210 days, seems to agree with what is expected in the case of a fast cometary activity decrease. A probable explanation of the dust velocity drop is a decrease of gas release from the nucleus surface, which implies a decrease of dust drag efficiency. Such a decrease of dust drag by gas has effects before on large grains than on small ones, so that a net effect is that the ejected size distribution becomes dominated by smaller and smaller grains: we exactly observe such a strong decrease of the size distribution power index. Due both to the decrease of dust release and to the fact that less and less massive grains are ejected, the mass loss rate shows a drop even larger than the velocity one, exactly what we observe (a factor 103 against 50). Moreover, the velocity drop should begin after the drop of the size distribution power index and, consequently, of the mass loss rate. In fact, even if the gas is able to drag less and less large grains, it drags them to the same terminal velocity. Only when the gas is unable to drag almost all the large grains, we can observe a drop of their velocity, which should be much faster than the corresponding index and mass drops. This is exactly what we observe for the velocity of 1-mm-sized grains. The strong correlation among all the physical outputs, and their good agreement with what is physically plausible, further
247
THE DEATH OF C/1996 Q1
points out the reliability of the time evolution of the physical outputs provided by the dust tail fading fit. 5. CONCLUSIONS
The results presented in the previous section allow us to select the most probable answers to the questions raised in the introduction about the fading of C/1996 Q1. Both the brightness distribution and turning axis of the observed dust tail changing in time are consistent with a dust tail model. This fact supports the assumption that most of the observed brightness is coming from dust, with probably negligible pollutions by gas and ions. The physical outputs provided by the fits of all the input images point out that the observed dust tail is a remnant dust tail, built up by the dust production that occurred before perihelion only, which becomes visible because it is no more masked by any bright coma produced during and after perihelion. Moreover, the features of the obtained physical outputs allow us to select the physical explanation most consistent with the time fading of the C/1996 Q1 dust tail. In fact, the time evolution of the dust mass loss rate allows us to answer the fundamental question raised in the Introduction: the obtained mass loss rate is inconsistent with a dust tail built up by the disintegration of the whole comet nucleus. In this case, we would wait for a strong peak of the mass loss rate in correspondence of the nucleus disintegration time. Over all, we would expect for much higher absolute values of the mass loss rate, which, on the contrary, remain at levels typical of short-period comets, with a very low dust to gas ratio (about 0.2 for an albedo times the phase function of 0.02). The total released dust mass during the perihelion approach is 1.2 3 109 kg, so that, according to the estimated dust to gas ratio, the total released mass would be 6 times larger, corresponding to a sphere of 250 m radius (for a nucleus bulk density of 100 kg m23), which can produce at most 7 3 1026 s21 water molecules from its surface, a loss rate 60 times lower than the observed one (Crovisier, pers. commun.). It seems improbable that a fractal nucleus surface also could increase the possible water loss rate of such a high factor. Even taking into account the uncertainties of the Q and Ap(a) parameters, it seems improbable that the total released dust mass would be 10 times larger, i.e., 1.2 3 1010 kg, with a dust to gas ratio equal to 2. In this case, the total released mass would be 1.5 times larger, corresponding to a sphere of 350 m radius (for a nucleus bulk density of 100 kg m23), which can produce at most 1.5 3 1027 s21 water molecules from its surface, a loss rate 30 times lower than the observed one (Crovisier, pers. commun.). We can conclude that the C/1996 Q1 nucleus mean radius was surely larger than 350 m. The strong drop following a very slow mass loss rate increase strongly favors the other interpretation: the nucleus was deactivated after
t 5 225 days, so that no new bright coma masked the normal dust tail released before perihelion. As was discussed in the previous section, the strong time correlation of the drops of all the physical outputs seems well consistent with what is expected in the case of a nucleus which is strongly decreasing the gas loss (usually called nucleus deactivation). However, the model outputs do not allow us to select the possible cause of such a deactivation. The strong dust velocity increase for 2100 , t , 225 days points out that the gas release from the nucleus surface may be dominated by seasonal effects (Crifo and Rodionov 1997), typical effects occurring in short-period comets. Also the changes in time of the size distribution power index before the drop starting at t 5 225 days might reflect changes of the active areas releasing the ejected dust. Finally, the same seasonal changes might have caused the nucleus deactivation, putting all the nucleus active areas in permanently night sides. However, the best tail fits are obtained for approximately hemispherical dust ejections, and not for strongly sunward anisotropic dust ejections, as would be expected for small active areas, which can be put in permanently night sides by the nucleus spin orientation. Hemispherical dust ejections are consistent with a nucleus activity dominated by seasons, as shown by the nucleus models of 46P/ Wirtanen, whose activity is probably dominated by seasons and characterized by an active area perhaps larger than half of the whole nucleus surface. However, the hemispherical dust ejection might point out that most of the Sunexposed nucleus was active during all the considered times, a fact in apparent contradiction with the sudden passage of all active areas in the nucleus night side. In other words, the results provided by the tail fit do not exclude a drop of the gas loss rate over all the surface independent of seasonal effects, e.g., due to dust mantling of the nucleus surface, or to exhaustion of the ice reservoir of C/1996 Q1. ACKNOWLEDGMENTS Standard stars were identified by means of the SIMBAD procedure, developed at the Observatoire Astronomique de Strasbourg. The referees M. F. A’Hearn and S. Larson suggested significant improvements to the previous version of the paper.
REFERENCES Burns, J. A., P. L. Lamy, and S. Soter 1979. Radiation forces on small particles in the Solar System. Icarus 40, 1–48. Crifo, J. F. 1991. Hydrodynamic models of the collisional coma. In Comets in the Post-Halley Era (R. L. Newburn, Jr., M. Neugebauer, and J. Rahe, Eds.), pp. 937–990. Kluwer, Dordrecht. Crifo, J. F. and A. V. Rodionov 1997. The dependence of the circumnuclear coma structure on the properties of the nucleus. I. Comparison between a homogeneous and an inhomogeneous spherical nucleus, with application to P/Wirtanen. Icarus 127, 319–353.
248
FULLE ET AL.
Delsemme, A. H. 1982. Chemical composition of cometary nuclei. In Comets (L. L. Wilkening, Ed.), pp. 85–130. Univ. of Arizona Press, Tucson. Dennerl, K., J. Englhauser, J. Truemper, and C. Lisse 1996. IAU Circ. 6495. Finson, M. L., and R. F. Probstein 1968. A theory of dust comets. I. Model and equations Astrophys. J. 154, 327–352. Fulle, M. 1989. Evaluation of cometary dust parameters from numerical simulations: Comparison with analytical approach and role of anisotropic emissions. Astron. Astrophys. 217, 283–297. Fulle, M., C. Bo¨hm, G. Mengoli, F. Muzzi, S. Orlandi, and G. Sette 1994. Current meteor production of Comet P/Swift–Tuttle. Astron. Astrophys. 292, 304–310. Fulle, M., L. Colangeli, V. Mennella, A. Rotundi, E. Bussoletti 1993. The sensitivity of the size distribution to the grain dynamics: Simulation of the dust flux measured by GIOTTO at P/Halley. Astron. Astrophys. 304, 622–630.
Fulle, M., G. Cremonese, K. Jockers, and H. Rauer 1992. The dust tail of Comet Liller 1988V. Astron. Astrophys. 253, 615–624. Green, D. W. E. 1996a. IAU Circ. 6455. Green, D. W. E. 1996b. IAU Circ. 6499, 6501, 6506, 6513. Jahn, J. 1996. IAU Circ. 6464. Marsden, B. G. 1996a. MPE Circ. 1996-R05. Marsden, B. G. 1996b. IAU Circ. 6464. Marsden, B. G. 1996c. IAU Circ. 6521. McDonnell, J. A. M., P. L. Lamy, and G. S. Pankiewicz 1991. Physical properties of cometary dust. In Comets in the Post-Halley Era (R. L. Newburn, Jr., M. Neugebauer, and J. Rahe, Eds.), pp. 937–990. Kluwer, Dordrecht. Misselt, K. A. 1996. Secondary photometric standards in selected northern dwarf-nova fields. Publ. Astron. Soc. Pacific 108, 146–165. Womack, M., and D. Suswal 1996. IAU Circ. 6485.