- Email: [email protected]
Halley
ICARUS 6 6 , 1 5 4 - 1 6 4
(1986)
A Model of Cometary Gas and Dust Production and Nongravitational Forces with Application to P/Halley F R A S E R P. F A N A L E AND JAMES R. S A L V A I L Planetar3, Geosciences Division, Hawaii Institute of Geophysics, University of Hawaii, Honohdu, Hawaii 96822 Received August 13, 1985: revised December 2, 1985 A program that c o m p u t e s gas and dust production rates and idealized nongravitational force c o m p o n e n t s has been developed and applied to the case of Comet Halley. We use a modified form of our earlier c o m e t model (F. P. Fanale and J. R. Salvail [(1984) h'arus 60,476-511] to which c o m a effects and a section on nongravitational forces have been added. The possibility of grain cohesion is also included. T h e s e models are used together with observations from 1910 and semiempirically derived data to investigate the effects of obliquity and thermal conductivity of the near surface nucleus on gas and dust production. The results indicate that the thermal conductivity of the nucleus is of the order of 105 ergs/cm-s-°K, which implies that the ice near the surface is in the crystalline form. A general method is presented for calculating the radii of cometary nuclei using theoretically derived and semiempirically derived nongravitational force c o m p o n e n t s . This m e t h o d is used to calculate possible radii for Comet Halley that depend on the model variation chosen. The m e t h o d used and the results presented herein should have greater significance and value when the observational data from Halley's current perihelion passage become available. ,c: 1986 Academic Press, Inc.
1. I N T R O D U C T I O N
In this paper we will use a modified form of our previous model (Fanale and Salvail, 1984), which now includes coma effects and a flux integration procedure, to calculate gas and dust production rates for Comet Halley for various values of obliquity and thermal conductivity. We will include an analytical procedure for calculating idealized nongravitational force components from first principles for a smooth, spherical, homogeneous nucleus. This procedure will yield the resultant nongravitational force, its components in the radial, transverse, and normal directions to the orbit plane, and the latitude and longitude where the resultant force acts, all as functions of heliocentric distance. Yeomans (1977) has already obtained P/Halley's nongravitational force components by fitting numerous observations of orbital position. We will develop a general method for calculating nuclear radii and apply it to the case of
P/Halley using Yeomans' semiempirical results and the results of our theoretical model. II. T H E G A S A N D D U S T P R O D U C T I O N M O D E L
Throughout this paper we will continually make reference to our previous paper (Fanale and Salvail 1984), which we will refer to hereafter as paper I. The model we have used is a modified version of the model used in paper I. We will not repeat the equations here but will present and discuss only the changes in the model. The model has been improved by including explicit calculations of the effects of the coma on the nucleus surface temperature. For this purpose we have used the simple radiative transfer model of Squyres et al. (1985). They assumed that a comet has a plane-parallel dusty atmosphere and that the coma has reached steady state. They further simplified the problem by assuming that all the dust has one characteristic radius and that the velocity of outflow is con154
0019-1035/86 $3.00 Copyright © 1986 by Academic Press. Inc. All rights of reproduction in any form reserved.
COMET MODEL
155
stant with distance from the nucleus. We will not rederive their results but will present only the final equations. The following equations in this section are from Squyres et al. (1985) unless otherwise indicated. The optical thickness of the coma is calculated from
the program. The coefficients of the second factor were computed by a least-squares procedure to be a = 1.247, b = 0.0233, c = - 2 . 0 4 5 × 10 4, and d = 6.52 × 10 -7. The radiation scattered over the nucleus by the dust coma is assumed to be uniform over the surface. The scattered radiation flux is given by
7" - 47rrvd"
Qdo'So Ims -- 64VdrR 2.
(1)
Qo is the dust particle production rate, # is the average extinction cross section, r is the radius of the nucleus, and vd is the particle terminal velocity. We calculated an average dust particle radius in a manner analogous to the derivation of Eq. (39) in paper I. The result is [1 - m][a~ 4 - m ) - a~in m) ]1/3. d = E[4 _ ml[a~c,_m, a~im, ]
(2)
ac is the critical particle radius, amin is the minimum radius, and m is the grain size distribution exponent. These are discussed in paper I. The program calculates a dust mass production rate based on the local gas fluxes and critical dust particle radii and the assumed dust/ice ratio. The dust mass production rate is translated into a dust particle production rate by dividing by the mass of an average-size spherical dust particle. The radius of such an average-size particle is calculated from Eq. (2), and the particle density is assumed. We used the results of our previous unpublished work on dust particle dynamics in cometary atmospheres to obtain an empirical relation for dust particle terminal velocity as a function of the temperature of the surface location from which the particle emanated. This relation is given as
vo =
[8kTc],/2 L 7"l"m J
(a + bTc + cT2c + dT3c).
(3)
The first factor on the right-hand side is just the thermal speed from kinetic theory. T~ is a characteristic average temperature that, if constant over the nucleus, would give the gas production rate computed by
(4)
Qd is the single scattering albedo, So is the solar constant at 1 AU, and R is the heliocentric distance. The infrared radiation flux emitted by the coma is calculated from
Itr = eT4(I - e-~i~).
(5)
e is the emissivity of the dust, assumed here to be 1.0, p is the S t e p h a n - B o l t z m a n n constant, and Td is the dust particle temperature given by
rd = [(1
Ad)S0] TM j ,
4eo.R2~
(6)
where Ad is the visible albedo of the dust. The infrared optical depth is calculated using Eq. (1) with particle cross sections appropriate for the infrared region of the spectrum. The program calculates surface and subsurface temperatures, gas fluxes, and dust fluxes for a grid of points on a rotating nucleus with an arbitrary axial orientation that varies along the orbit. Each grid point represents an area sector of the nucleus. The gas and dust fluxes from each gridpoint are multipled by the area of the sector, and the results are summed over the hour angle and then the latitude to yield the total gas and dust production rates for the comet at a given orbital position. The calculations are then repeated for different orbital positions. Preliminary information was obtained utilizing the program developed for paper I. Results for surface temperature, gas flux, dust flux, mantle thickness, and other pertinent data were obtained for the entire orbit of P/Halley. It was found that using a nu-
156
FANALE AND SALVAIL
clear radius of 2.5 km, the maximum mantle thickness is less than 50/xm. Such a thickness is not physically realistic and gives only an indication of the mass of residual dust. The thickness was calculated under the assumption that the dust layer is distributed uniformly like a fluid. However, larger grains contain a significant fraction of the mass of the nonvolatile component, although they are much fewer in number. Also, the smallest departure from absolute smoothness of the surface would cause uneven deposition of nonvolatile material, causing a large fraction of the surface to be exposed. In view of the above discussion, we have neglected any effect that such a thin discontinuous dust mantle might have in modifying the temperature of the ice surface or impeding the flow of gases. Therefore, in the present model we calculate the surface temperature of the ice for the case of a bare surface. This eliminates some of the complexity of the earlier model in calculating temperatures and reduces computer time. The thickness of liberated dust, the dust blowoff parameters, the dust flux, and the residual thickness of dust are calculated using the same equations as in the previous model. However, these items, except for the dust flux, are not of great interest in this paper, but are included in the program only for the purpose of calculating the dust production rate. An additional feature, not considered in paper I, is included in the calculation of the dust flux. In that paper we did not consider the possible cohesion of residual dust particles, assuming that there are no intergrain forces. Experiments simulating the emission of cometary gas and dust have been performed at the Jet Propulsion Laboratory (Fanale et al., 1983; Sutton et al., 1983). Results show that the sublimate liberated from the ice-dust matrix forms a fluffy residue with considerable cohesive strength, not unlike cotton candy. The cohesive force is thought to be electrostatic in nature, but the detailed physics of this phenomenon is still unknown. This fluffy resi-
due does not appear to adhere to the ice, only to cohere with itself. Therefore, we envision that when dust blowoff first occurs as a comet nears the Sun, the thin layer of fluffy residue is lifted from the nucleus as a series of patches. After this sheet of residue is lifted from the nucleus, additional grains liberated from the ice-dust matrix which are less than the critical size are instantly blown off by the escaping gases. The larger grains which remain may be too few in number to form the fluffy sublimate and instead, probably exist and eventually are blown off as isolated grains. In the computer model, when the lithostatic pressure of the thin layer of residue becomes less than the vapor pressure at the ice surface, the thin mantle is ejected, that is, when the actual thickness becomes less than the "critical thickness" given by Eq. (10) in paper I. This may result in a sudden surge in dust emission within a localized area. Within the range of heliocentric distances where residual dust cannot accumulate, the mechanism of dust blowoff is the same as discussed in paper I and others papers on the subject. While this change in the model is not too important in the case of Halley's comet we believe it is an important improvement in our earlier model in that it allows predicted cometary behavior to be bracketed. Calculation of global gas fluxes also makes possible the calculation of the nongravitational force components from first principles. The nongravitational force components for Comet Halley have already been accurately calculated from observations of P/Halley's orbital perturbations. Comparison of the results from the two methods could yield valuable information. Accordingly, we now develop the equations for the nongravitational force model. III. THE NONGRAVITATIONALFORCE MODEL Nongravitational forces cannot generally be accurately predicted by purely theoretical analysis because of various factors such
COMET MODEL as irregular surface topography, surface material inhomogeneity, and more or less nonspherical nuclei. Current methods of calculating nongravitational accelerations rely on many accurate observations of orbital position. H o w e v e r , if the various irregularities and inhomogeneities of cometary nuclei are ignored, idealized nongravitational forces can be calculated by theoretical means. In this section we derive equations for calculating the nongravitational force, the location where the resultant force acts in terms of latitude and longitude, and the radial, transverse, and normal components of the force. This analysis differs from those done previously by Marsden (1968, 1969), Marsden et al. (1973), and Yeomans (1977) in that it does not require observations of orbital position but is based on first principles of physics. This model is not intended to replace the aforementioned ones, since utilizing direct observations is obviously superior to utilizing a theoretical model which contains many unknown parameters whose values can only be guessed at. The value of a theoretical model is that the effects of various sensitive unknown parameters, such as obliquity and nuclear radius, can be evaluated quantitatively. This can be useful in determining reasonable bounds for these parameters. The nongravitational force is the resultant reaction force due to variable gas production o v e r the surface caused by uneven heating. The contribution to the nongravitational force at a given surface location is d F = qvodA,
(7)
where q is the gas flux, v0 is the radial efflux velocity at the surface, and d A is a differential area. The gas flux is given by [
q :/'v
m
t2~J
] 1/2
'
(8)
where Pv is the vapor pressure at the surface, m is the molecular mass, k is Boltzmann's constant, and T is the surface tern-
157
perature. The radial efflux velocity is the mean molecular speed multiplied by a constant determined by Delsemme and Miller (1971), [ 8 k r ] '~2
v0 = 0.6 L~--mmJ "
(9)
The vapor pressure can be expressed as Pv = A e B/r,
(10)
according to the Clausius-Clapeyron equations, where A and B are constants that depend on the gas. ff Eqs. (8), (9), and (10) are substituted into Eq. (7) along with the differential area in spherical coordinates and the result is integrated over the surface of the sphere, the nongravitational force is f
--
-"/7"
-Tr
e -Rff(°,*) cos OdOd4~, (11) where r, 0, and 4~ are nuclear radius, latitude, and hour angle, respectively. This is the magnitude of the nongravitational force. The location where this force acts is found in a way analogous to finding a centroid of some quantity. The latitude and hour angle of this point are found from 0c -
rr------F-
-~
e-B/r(°,e~tO cos OdOd4~, (12) 4,c = ~ - - - - p -
-~
e 8/r1°'~)4~ cos OdOdeb. (13) In the above equations T(O, 4~) is the variable surface temperature, which is assumed to be known at each point from the results of the previous section. It now remains to calculate the components of the nongravitational force in the radial, transverse, and normal directions, with respect to the comet's orbital plane. The latitude of the subsolar point, qJ, is calculated from sin 0 = sin I sin(fI + o0,
(14)
158
FANALE AND SALVAIL Using basic trigonometric identities for sin(/3 + I ) , cos fl, and the a b o v e expression for sin Ar, and simplifying with other basic identities yield the result
i
sin Ap --- sin(0 - 0c) cos I + cos(0 - 0c) sin ~bc sin I. V
(19)
The normal c o m p o n e n t of the force is
Ap tO Sort
Fn = F sin Ap = F[sin(to - 0~) cos I + cos(to - 0~) sin ~bc sin I].
(20)
An expression for the transverse component can be found in terms of trigonometry, but the derivation is more complicated than for the normal c o m p o n e n t . It is easier to FIG. 1. Geometry for calculation of nongravitational force components.
use
Ft = -+ N/F 2 - F~ - F~. where the obliquity, I, and the angle between the ascending node on the equator and the subsolar point at perihelion, 11, are specified, and the true anomaly, a, has been previously calculated for a given orbital position. Referring to Fig. 1, the angle b e t w e e n the line of action of F at 0~ and ~b~ and the radial direction is found using the spherical law of cosines as COS A r -- COS(t0 -
0c) c o s qSc.
(15)
Then the radial c o m p o n e n t of the force is Fr = F cos(to - 0c) cos 4~c.
(16)
It is convenient to consider next the normal c o m p o n e n t of force. Again referring to Fig. 1, the angle fl is e x p r e s s e d using the spherical law of sines as sin/3 =
sin(to - 0c) sin Ar sin(to - 0c) = k/1 - cos2(to - 0~) cos2~5c" (17)
An expression for A p , which is the angle b e t w e e n the line of action of F and the orbital plane, is found from the spherical law of sines as sin Ap = sin Ar sin(/3 + I).
(18)
(21)
In Eq. (21) the plus sign is used for prograde rotation or for obliquities less than 90 ° . The minus sign is used for retrograde rotation or for obliquities greater than 90 ° . The dust flux is a s s u m e d to have no effect on the nongravitational force since dust particles are ejected from the nucleus from zero velocity. Therefore, in the limit, mom e n t u m transfer takes place between dust particles and gas molecules after the particles are ejected and not with the nucleus. H o w e v e r , if dust particles are broken away, rather than being simply lifted, the energy of the gas flow would be decreased accordingly. This would cause the nongravitational force to be s o m e w h a t lower than given by this analysis. Theoretically determined nongravitational force c o m p o n e n t s , obtained from the a b o v e equations, can be used with empirical results for individual comets to estimate the sizes of c o m e t a r y nuclei. Marsden et al. (1973) used a semiempirical procedure to calculate nongravitational accelerations for several comets. In the equation of motion the magnitudes of the nongravitational acceleration c o m p o n e n t s are expressed as (Yeomans, 1977) ar = A I g ( R )
(22)
COMET MODEL
159
TABLEI
IV. R E S U L T S A N D D I S C U S S I O N
VARIABLE PARAMETERS SELECTED FOR RUNS
Results were obtained for various combinations of obliquity and thermal conductivity of the nucleus, as indicated in Table I. Values of important parameters are listed in Table II. Surface and subsurface temperatures, gas fluxes, and dust fluxes were calculated for a global grid of points for each 15° hour angle and 10° latitude and integrated to obtain gas and dust production rates. Nongravitational force components were obtained from the equations in the preceding section. These calculations were carried out over five complete orbits, starting and ending at aphelion, using Comet Halley's well-known orbital parameters. Calculations over five orbits allowed surface temperatures to converge to ~ 2 ° at
Model
Obliquity (deg)
Thermal conductivity
1
0
8 x 104
2 3 4
45 0 45
8 x 104 8 x lO3 8 x 103
and (23)
at = A 2 g ( R ) ,
fRj,j-, In the equation for g ( R ) , a , m , n, k , and R0 are constant parameters for a given comet. The parameters A1 and A2 for a given c o m e t are obtained from a leastsquares differential correction procedure. The theoretically determined nongravitational force components as derived earlier can be expressed in the form Fi
=
Ci re,
(25)
where Fi represents the ith nongravitational force component, and Ci is a parameter whose value is determined by the gas being sublimed and the array of surface temperatures. The corresponding acceleration is ai -
Fi m
-
Cir2 ]pn~rr 3
-
3Ci 4pnrrr"
(26)
On is the density of the cometary nucleus; Ci is obtained from Eq. (25). F~ is calculated independently from Eq. (16) if i = r or from Eqs. (16), (20), and (21) if i = t, and r is the assumed nuclear radius in the physical model. If a; is set equal to either ar or at in Eqs. (22) and (23), the actual radius of the cometary nucleus can be estimated from 3C,. r = - 4pnTrai
(27)
for i = r or i = t. This procedure will be used to estimate the radius o f Comet Halley in the next section.
T A B L E II VALUES OF IMPORTANT PARAMETERS a
Name Orbital and rotational parameters Semimajor axis Eccentricity Rotation period A r g u m e n t of subsolar meridian at perihelion Optical properties N u c l e u s albedo N u c l e u s emissivity Single scattering albedo of dust in the visible Single scattering albedo of d u s t in the infrared D u s t extinction c r o s s section in the visible Dust extinction cross section in the infrared Dust emissivity N u c l e u s properties Radius D u s t to ice ratio E x p o n e n t o f grain size distribution M a x i m u m grain size M i n i m u m grain size
Value
17.97 A U 0.9673 41.5 h b 340 °b 0.1 0.9 0.9' -0 C 6 × 10 -sc 2 x 10 -9c
1.0" 2.5 k m 0.5 -4.2 1.0 cm 0.5 x 10 4 cm
See Table I for obliquity and thermal conductivity. b Sekanina and L a r s o n (1984). c S q u y r e s et al. (1985).
160
FANALE AND SALVAIL
blc
810 9
8
CI 2
~7
--7
11
~s
~9
T~8
F, Fn
F,
$
$3
Z
g
2
J
6
I i
I
2
3
4
Heliocentric Dista~e (AU }
t
5
t
6
Heliocentric Distance (AU }
4t
i
1
2 3 4 5 Heliocentric Distance (AU)
FIG. 2. (a) Gas production rate, (b) dust production rate, and (c) nongravitational forces vs heliocentric distance for Model 2.
aphelion. Temperatures at smaller heliocentric distances were more accurate. Increased accuracy was not attempted due to long running times. The figures are numbered according to the model numbers shown in Table I. The results include the effects of attenuation, scattering, and thermal emission of radiation by dust in the coma. These effects cause the nucleus surface to be more isothermal. For example at 2 AU the maximum and minimum nucleus surface temperatures differ by about 100°K for the low-obliquity, high-thermal-conductivity case while at perihelion the difference is only about 25°K. The optical depth tends to lower slightly the dayside temperatures, while the radiation scattering and thermal emission by the dust significantly increase nightside temperatures. At perihelion the maximum surface temperature is lowered by about 2°K due to an optical depth of - 0 . 4 . H o w e v e r , the coma has only a small effect on the gas and dust production rates because the sublimation rates are governed essentially by the maximum temperatures due to the exponential dependence of sublimation rate on temperature. We will now focus, in turn, on the two
parameters chosen for investigation. For the obliquity we used 0 and 45 °. Two important characteristics distinguish the 0 and 45 ° obliquity plots. The first is that the 45 ° results show a small preperihelion asymmetry in the gas and dust production rates inside ~ 1 AU and an increasingly large postperihelion a s y m m e t r y beyond - 1 AU. The 0 ° results show only an increasing postperihelion asymmetry. The preperihelion asymmetry in the 45 ° results looks remarkably similar to the lightcurve produced by Yeomans (1977) from 1910 data except that the asymmetry is smaller. H o w e v e r , a consensus later developed among several observers that the dip in the lightcurve after perihelion was an artifact of the unusual observing conditions. This is discussed in detail by Newburn (1981). As a result the postperihelion portion of the curve inside 1 AU was modified to a consistent postperihelion a s y m m e t r y more like the 0 ° curves. The lightcurve produced by Yeomans (1977) and the postperihelion modification of it now being used for models at JPL are also shown in Newburn (1981). The second important characteristic is seen in the nongravitational force components. For the 45 ° cases Figs. 2c and 3c show that the
COMET MODEL
161
b,o
81C
9 GI 2
8
--7
I1
~7
~s
I
2 3 4 5 Heliocentric Distance(AU)
6
2
~6
I
5
O(
1 Heliocentric D~stor~e(AU)
2 3 4 5 Heliocentric Distonce(AU)
6
FIG. 3. (a) Gas productionrate, (b) dust productionrate, and (c) nongravitationalforces vs heliocentric distance for Model 4.
normal component of the nongravitational force is usually of the same order of magnitude as the transverse component. The normal components in both Figs. 2c and 3c appear to have discontinuities near perihelion. These are the positions where the normal component of the force reverses direction. The dotted lines represent the orbital positions where the normal component of the force is reversed, pointing south with respect to the orbital plane. The difference in the appearance of the two plots and, in particular, the range of heliocentric distances for which the normal component of the force is reversed are due only to the value used for the thermal conductivity. A larger thermal conductivity results in a larger angle between the subsolar point and the point where the resultant nongravitational force acts. This affects the angles 0c and 4~c in Eq. (20). The normal component of the force has the effect of changing the inclination of a comet's orbit. However, this effect would be small and is probably not able to be distinguished from the effects of planetary perturbations. Yeomans (1977) concluded that the acceleration component normal to the orbit plane has a negligible
effect on the orbital motion, a result that seems to be typical of other short period comets. This results from the fact that the normal component of the force has essentially no effect on the orbital velocity and, hence, the orbital period. For the 0 ° cases the normal component is always at least three orders of magnitude smaller than the transverse component. Ideally the normal component should be zero for 0° obliquity, but small errors in the computation of the temperatures result in finite but relatively insignificant normal component forces, which we have neglected in plotting. The asymmetry between the post- and preperihelion gas production rates, especially inside l AU, suggests a conclusion with regard to the thermal conductivity of the nucleus. Figures 4 and 5 show results of models that differ only in the value assumed for the thermal conductivity of the nucleus (see Table I). The results obtained using the higher thermal conductivity have a significantly higher asymmetry in gas production than those obtained using the lower thermal conductivity. Observations tend to support the results with the larger asymmetry in gas production. Fanale and Salvail
162
F A N A L E AND S A L V A I L
blo
alO
9 C12 F,
8 ~7
~7
11
Ft
$o d:
°t 1
I
i
2 3 4 5 HeLioce~ric Distar~e (AU)
6
I
I
6J
,
40
1
Heliocentric Distot~e (AU)
2 3 4 5 Heliocmtri(: Distance (AU)
,
6
FIG. 4. (a) Gas production rate, (b) dust production rate, and (c) nongravitational forces vs heliocentric distance for Model 1.
(1984) concluded that postperihelion asymmetries are largely caused by heat conduction into the nucleus. In this work we find that a higher thermal conductivity produces a larger postperihelion asymmetry. This fact has been discussed for some time, though here we do a more quantitative study. Previous physical models of come-
tary nuclei h a v e utilized thermal c o n d u c t i v ities ranging f r o m on the o r d e r o f 102ergs/ c m - s e c - ° K for u n c o m p a c t e d s n o w and dust to 105 e r g s / c m - s e c - ° K for crystalline w a t e r ice. By c o m p a r i n g the near perihelion a s y m m e t r i e s in gas p r o d u c t i o n for our m o d e l results with the near perihelion a s y m m e t r y in the o b s e r v e d light c u r v e
810
i
9 8
C
~-7
12 F, 11
F,
~g o
i
Z 2
5
1
o,
I
2 3 4 5 Heliocentric Distance (AU)
I
2
3
4
Heliocentric Distar~e (AU)
5
J 6
4
1
2 3 4 5 Heliocentric Distonce (AU)
FIG. 5. (a) Gas production rate, (b) dust production rate, and (c) nongravitational forces vs heliocentric distance for Model 3.
COMET MODEL 8
--
a- 6
--
A
O
I
5 0
1.0
I 2.0
Heliocentric Distance (Au)
FIG. 6. Gas production rate vs heliocentric distance for Newburn's (1981) model.
given by Yeomans (1977), we conclude that Comet Halley probably has a thermal conductivity near the higher end of the range given above. Therefore, Halley's nucleus is probably crystalline ice, at least near the surface. The results all show gas production rates lower than those of the Newburn model results shown in Fig. 6. Part of the discrepancy is due to the fact that Newburn arbitrarily assumed that 20% of the volatile c o m p o n e n t consists of volatiles other than water having an average molecular weight of 44 amu. We believe that a pure water model is just as realistic in view of the discussion given by Delsemme in Yeomans (1977), in which he concludes on the basis of statistics that "all short period comets lose all gases more volatile than water very fast, at least from the outer layers of the nucleus that are heated enough to participate in the sublimations." The remaining discrepancy could be the result of a number of assumptions, other than nuclear radius, that N e w b u r n (1981) made in his model. Finally, we will use the procedure outlined in the previous section to estimate the radius of Comet Halley. We will use the transverse force component, since Yeomans (1977) shows that the transverse component is determined far better than the radial component. In that paper A2 = 0.0159
163
x 10 -8 , ~ = 0.111262, r0 = 2.808 AU, m = 2.15, n = 5.093, and k -- 4.6142. We performed calculations for nine points in the orbit inside 2 AU before and after perihelion, for both cases having the high thermal conductivity. The values for g(R) were calculated from Eq. (24) for the appropriate value of R (AU). The values for at were obtained by first converting the parameter A2 to cgs units to get A2 = 3.1736 x 10 -7 cm/sec 2 and then multiplying by the corresponding value for g(R). The numbers for Ft were provided by our program, and they were divided by r 2 (r = 2.5 × 105 cm) to get the corresponding values for Ct. We must emphasize here that any realistic number could be used for the radius in the computer program, since Ct is independent of the assumed radius. The values for the computed nuclear radius were then calculated from Eq. (27), assuming that the bulk density of the nucleus, pn, is 1.0 g/cm 3. Using a dust to ice ratio of 0.5 and assuming that the ice and silicate are nonporous, the bulk nuclear density would be - 1 . 5 g/cm 3. A bulk density of 1.0 g/cm 3 is equivalent to assuming that both the ice and silicate have a porosity of - 0 . 6 5 . If an interested reader prefers using a different bulk density, the results for the radius can be easily converted. The computed values of the radius were integrated over time, and the result was divided by the elapsed time for the portion of the orbit chosen. This procedure yielded averaged values for the radius of 3.135 km for the 0 ° obliquity case and 1.6435 km for the 45 ° obliquity case. Our model and results have several shortcomings. The agreement between our model results and the facts indicated by observations was not the best. N o n e of the model variations satisfied all the observational constraints. Among the low-obliquity model variations, the low-thermal-conductivity case (Model 3) shows gas and dust production rates that are more consistent with the facts that Comet Halley was first seen in 1910 at - 3 . 4 AU (preperihelion) and last seen at - 5 . 2 AU (postperihelion)
164
FANALE AND SALVAIL
(Newburn, 1981). However, these model variations show the observed asymmetry near perihelion to be significantly smaller than indicated by observations. The highthermal-conductivity case (Model 1) produces a near perihelion asymmetry that agrees well with observations but has considerably lower gas and dust production rates beyond - 3 AU, especially preperihelion, than are indicated by observations. One possible explanation for the supposedly higher gas production rates beyond - 3 AU is that electrostatic forces produced by the solar wind cause small residual dust grains to be emitted from the nucleus in the absence of the usual gas production levels. This phenomenon is analyzed and discussed in detail by Flammer et al. (1985). This dust could then be interpreted by an observer as a higher gas production rate. The accuracy of the radius calculations depends not only upon the results of our model but upon Yeoman's semiempirical results. Although errors would arise due to the idealizations required for the theoretical analysis, we also believe that the semiempirical results are inaccurate. It is well known that comets generally have significant asymmetries in gas and dust production, and this is certainly true for Comet Halley. This strongly implies that nongravitational forces and accelerations must be asymmetric with respect to perihelion. Therefore, the dependence of the nongravitational acceleration on heliocentric distance must be different between the preperihelion and postperihelion legs of the orbit. This is not so for any of the semiempirical formulations currently available. Further refinements are obviously needed before accurate calculations can be made using this method. Although this work has provided results of questionable accuracy,
it has indicated the methodology by which physical and semiempirical models can be used to reveal characteristics of individual comets. REFERENCES DELSEMME, A. H., AND D. C. MILLER (1971). Physico-chemical phenomena in Comets - - III. Planet. Space Sei. 19, 1229-1257. FANALE, F. P., AND J. R. SALVAIL 0984). An idealized short period comet model: Surface insolation, H20 flux, dust flux and mantle evolution. Icarus 60, 476-511. FANALE, F. P., R. S. SAUNDERS, B. BANERDT, J. STEVENS, AND E. LANE (1983). Properties o f Sirnulated Comet Mantle, NASA TM-86246, pp. 56-58. FLAMMER, K. R., B. JACKSON, AND D. A. MENDIS
(1985). On the brightness variations of Comet Halley at large heliocentric distances. Preprint. MARSDEN, B. G. (1968). Comets and nongravitational forces. Astron. J. 73, 367-379. MARSDEN, B. G. (1969). Comets and nongravitational forces, I1. Astron. J. 74, 720-734. MARSDEN, B. G., Z. SCKANINA,AND D. K. YEOMANS (1973). Comets and nongravitational forces V. Astron. J. 78, 211-225.
NEWaURN, R. L. (1981). A semi-empirical photometric theory of cometary gas and dust production: Application to P/Halley's gas production rates. Comet Halley Dust and Gas Environment Workshop. ESA SP-174. NEWBURN, R. L., AND D. K. YEOMANS (1982). Halley's Comet. Annu. Rev. Earth Planet. Sci. 10. SEKAN1NA, Z., AND S. M. LARSON (1984). Coma Morphology and dust-emission pattern of periodic comet Halley. II. Nucleus spin vector and modeling of major dust features in 1910. Astron. J. 89, 1408-1447. SQUYRES, S. W., C. P. McKAY, AND R. T. REYNOLDS, (1985). Temperatures within comet nuclei.
Preprint. SUTTON, S. T., J. STEPHENS, R. S. SAUNDERS, AND F.
P. FANALE, (1983). Properties of residues from sublimed dirty ice: Implications for comets and Martian Polar deposits. Bull. Amer. Astron. Soc. 15, 847. WHIPPLE, F. L. (1980). Periodic Comet Halley. IAU Circ. 3459. YEOMANS, D. K. (1977). Comet Halley and nongravitational forces. In Comets, Asteroids" and Meteorites (A, R. Delsemme, Ed.), pp. 61-64. University of Toledo, Toledo.
(1986)
A Model of Cometary Gas and Dust Production and Nongravitational Forces with Application to P/Halley F R A S E R P. F A N A L E AND JAMES R. S A L V A I L Planetar3, Geosciences Division, Hawaii Institute of Geophysics, University of Hawaii, Honohdu, Hawaii 96822 Received August 13, 1985: revised December 2, 1985 A program that c o m p u t e s gas and dust production rates and idealized nongravitational force c o m p o n e n t s has been developed and applied to the case of Comet Halley. We use a modified form of our earlier c o m e t model (F. P. Fanale and J. R. Salvail [(1984) h'arus 60,476-511] to which c o m a effects and a section on nongravitational forces have been added. The possibility of grain cohesion is also included. T h e s e models are used together with observations from 1910 and semiempirically derived data to investigate the effects of obliquity and thermal conductivity of the near surface nucleus on gas and dust production. The results indicate that the thermal conductivity of the nucleus is of the order of 105 ergs/cm-s-°K, which implies that the ice near the surface is in the crystalline form. A general method is presented for calculating the radii of cometary nuclei using theoretically derived and semiempirically derived nongravitational force c o m p o n e n t s . This m e t h o d is used to calculate possible radii for Comet Halley that depend on the model variation chosen. The m e t h o d used and the results presented herein should have greater significance and value when the observational data from Halley's current perihelion passage become available. ,c: 1986 Academic Press, Inc.
1. I N T R O D U C T I O N
In this paper we will use a modified form of our previous model (Fanale and Salvail, 1984), which now includes coma effects and a flux integration procedure, to calculate gas and dust production rates for Comet Halley for various values of obliquity and thermal conductivity. We will include an analytical procedure for calculating idealized nongravitational force components from first principles for a smooth, spherical, homogeneous nucleus. This procedure will yield the resultant nongravitational force, its components in the radial, transverse, and normal directions to the orbit plane, and the latitude and longitude where the resultant force acts, all as functions of heliocentric distance. Yeomans (1977) has already obtained P/Halley's nongravitational force components by fitting numerous observations of orbital position. We will develop a general method for calculating nuclear radii and apply it to the case of
P/Halley using Yeomans' semiempirical results and the results of our theoretical model. II. T H E G A S A N D D U S T P R O D U C T I O N M O D E L
Throughout this paper we will continually make reference to our previous paper (Fanale and Salvail 1984), which we will refer to hereafter as paper I. The model we have used is a modified version of the model used in paper I. We will not repeat the equations here but will present and discuss only the changes in the model. The model has been improved by including explicit calculations of the effects of the coma on the nucleus surface temperature. For this purpose we have used the simple radiative transfer model of Squyres et al. (1985). They assumed that a comet has a plane-parallel dusty atmosphere and that the coma has reached steady state. They further simplified the problem by assuming that all the dust has one characteristic radius and that the velocity of outflow is con154
0019-1035/86 $3.00 Copyright © 1986 by Academic Press. Inc. All rights of reproduction in any form reserved.
COMET MODEL
155
stant with distance from the nucleus. We will not rederive their results but will present only the final equations. The following equations in this section are from Squyres et al. (1985) unless otherwise indicated. The optical thickness of the coma is calculated from
the program. The coefficients of the second factor were computed by a least-squares procedure to be a = 1.247, b = 0.0233, c = - 2 . 0 4 5 × 10 4, and d = 6.52 × 10 -7. The radiation scattered over the nucleus by the dust coma is assumed to be uniform over the surface. The scattered radiation flux is given by
7" - 47rrvd"
Qdo'So Ims -- 64VdrR 2.
(1)
Qo is the dust particle production rate, # is the average extinction cross section, r is the radius of the nucleus, and vd is the particle terminal velocity. We calculated an average dust particle radius in a manner analogous to the derivation of Eq. (39) in paper I. The result is [1 - m][a~ 4 - m ) - a~in m) ]1/3. d = E[4 _ ml[a~c,_m, a~im, ]
(2)
ac is the critical particle radius, amin is the minimum radius, and m is the grain size distribution exponent. These are discussed in paper I. The program calculates a dust mass production rate based on the local gas fluxes and critical dust particle radii and the assumed dust/ice ratio. The dust mass production rate is translated into a dust particle production rate by dividing by the mass of an average-size spherical dust particle. The radius of such an average-size particle is calculated from Eq. (2), and the particle density is assumed. We used the results of our previous unpublished work on dust particle dynamics in cometary atmospheres to obtain an empirical relation for dust particle terminal velocity as a function of the temperature of the surface location from which the particle emanated. This relation is given as
vo =
[8kTc],/2 L 7"l"m J
(a + bTc + cT2c + dT3c).
(3)
The first factor on the right-hand side is just the thermal speed from kinetic theory. T~ is a characteristic average temperature that, if constant over the nucleus, would give the gas production rate computed by
(4)
Qd is the single scattering albedo, So is the solar constant at 1 AU, and R is the heliocentric distance. The infrared radiation flux emitted by the coma is calculated from
Itr = eT4(I - e-~i~).
(5)
e is the emissivity of the dust, assumed here to be 1.0, p is the S t e p h a n - B o l t z m a n n constant, and Td is the dust particle temperature given by
rd = [(1
Ad)S0] TM j ,
4eo.R2~
(6)
where Ad is the visible albedo of the dust. The infrared optical depth is calculated using Eq. (1) with particle cross sections appropriate for the infrared region of the spectrum. The program calculates surface and subsurface temperatures, gas fluxes, and dust fluxes for a grid of points on a rotating nucleus with an arbitrary axial orientation that varies along the orbit. Each grid point represents an area sector of the nucleus. The gas and dust fluxes from each gridpoint are multipled by the area of the sector, and the results are summed over the hour angle and then the latitude to yield the total gas and dust production rates for the comet at a given orbital position. The calculations are then repeated for different orbital positions. Preliminary information was obtained utilizing the program developed for paper I. Results for surface temperature, gas flux, dust flux, mantle thickness, and other pertinent data were obtained for the entire orbit of P/Halley. It was found that using a nu-
156
FANALE AND SALVAIL
clear radius of 2.5 km, the maximum mantle thickness is less than 50/xm. Such a thickness is not physically realistic and gives only an indication of the mass of residual dust. The thickness was calculated under the assumption that the dust layer is distributed uniformly like a fluid. However, larger grains contain a significant fraction of the mass of the nonvolatile component, although they are much fewer in number. Also, the smallest departure from absolute smoothness of the surface would cause uneven deposition of nonvolatile material, causing a large fraction of the surface to be exposed. In view of the above discussion, we have neglected any effect that such a thin discontinuous dust mantle might have in modifying the temperature of the ice surface or impeding the flow of gases. Therefore, in the present model we calculate the surface temperature of the ice for the case of a bare surface. This eliminates some of the complexity of the earlier model in calculating temperatures and reduces computer time. The thickness of liberated dust, the dust blowoff parameters, the dust flux, and the residual thickness of dust are calculated using the same equations as in the previous model. However, these items, except for the dust flux, are not of great interest in this paper, but are included in the program only for the purpose of calculating the dust production rate. An additional feature, not considered in paper I, is included in the calculation of the dust flux. In that paper we did not consider the possible cohesion of residual dust particles, assuming that there are no intergrain forces. Experiments simulating the emission of cometary gas and dust have been performed at the Jet Propulsion Laboratory (Fanale et al., 1983; Sutton et al., 1983). Results show that the sublimate liberated from the ice-dust matrix forms a fluffy residue with considerable cohesive strength, not unlike cotton candy. The cohesive force is thought to be electrostatic in nature, but the detailed physics of this phenomenon is still unknown. This fluffy resi-
due does not appear to adhere to the ice, only to cohere with itself. Therefore, we envision that when dust blowoff first occurs as a comet nears the Sun, the thin layer of fluffy residue is lifted from the nucleus as a series of patches. After this sheet of residue is lifted from the nucleus, additional grains liberated from the ice-dust matrix which are less than the critical size are instantly blown off by the escaping gases. The larger grains which remain may be too few in number to form the fluffy sublimate and instead, probably exist and eventually are blown off as isolated grains. In the computer model, when the lithostatic pressure of the thin layer of residue becomes less than the vapor pressure at the ice surface, the thin mantle is ejected, that is, when the actual thickness becomes less than the "critical thickness" given by Eq. (10) in paper I. This may result in a sudden surge in dust emission within a localized area. Within the range of heliocentric distances where residual dust cannot accumulate, the mechanism of dust blowoff is the same as discussed in paper I and others papers on the subject. While this change in the model is not too important in the case of Halley's comet we believe it is an important improvement in our earlier model in that it allows predicted cometary behavior to be bracketed. Calculation of global gas fluxes also makes possible the calculation of the nongravitational force components from first principles. The nongravitational force components for Comet Halley have already been accurately calculated from observations of P/Halley's orbital perturbations. Comparison of the results from the two methods could yield valuable information. Accordingly, we now develop the equations for the nongravitational force model. III. THE NONGRAVITATIONALFORCE MODEL Nongravitational forces cannot generally be accurately predicted by purely theoretical analysis because of various factors such
COMET MODEL as irregular surface topography, surface material inhomogeneity, and more or less nonspherical nuclei. Current methods of calculating nongravitational accelerations rely on many accurate observations of orbital position. H o w e v e r , if the various irregularities and inhomogeneities of cometary nuclei are ignored, idealized nongravitational forces can be calculated by theoretical means. In this section we derive equations for calculating the nongravitational force, the location where the resultant force acts in terms of latitude and longitude, and the radial, transverse, and normal components of the force. This analysis differs from those done previously by Marsden (1968, 1969), Marsden et al. (1973), and Yeomans (1977) in that it does not require observations of orbital position but is based on first principles of physics. This model is not intended to replace the aforementioned ones, since utilizing direct observations is obviously superior to utilizing a theoretical model which contains many unknown parameters whose values can only be guessed at. The value of a theoretical model is that the effects of various sensitive unknown parameters, such as obliquity and nuclear radius, can be evaluated quantitatively. This can be useful in determining reasonable bounds for these parameters. The nongravitational force is the resultant reaction force due to variable gas production o v e r the surface caused by uneven heating. The contribution to the nongravitational force at a given surface location is d F = qvodA,
(7)
where q is the gas flux, v0 is the radial efflux velocity at the surface, and d A is a differential area. The gas flux is given by [
q :/'v
m
t2~J
] 1/2
'
(8)
where Pv is the vapor pressure at the surface, m is the molecular mass, k is Boltzmann's constant, and T is the surface tern-
157
perature. The radial efflux velocity is the mean molecular speed multiplied by a constant determined by Delsemme and Miller (1971), [ 8 k r ] '~2
v0 = 0.6 L~--mmJ "
(9)
The vapor pressure can be expressed as Pv = A e B/r,
(10)
according to the Clausius-Clapeyron equations, where A and B are constants that depend on the gas. ff Eqs. (8), (9), and (10) are substituted into Eq. (7) along with the differential area in spherical coordinates and the result is integrated over the surface of the sphere, the nongravitational force is f
--
-"/7"
-Tr
e -Rff(°,*) cos OdOd4~, (11) where r, 0, and 4~ are nuclear radius, latitude, and hour angle, respectively. This is the magnitude of the nongravitational force. The location where this force acts is found in a way analogous to finding a centroid of some quantity. The latitude and hour angle of this point are found from 0c -
rr------F-
-~
e-B/r(°,e~tO cos OdOd4~, (12) 4,c = ~ - - - - p -
-~
e 8/r1°'~)4~ cos OdOdeb. (13) In the above equations T(O, 4~) is the variable surface temperature, which is assumed to be known at each point from the results of the previous section. It now remains to calculate the components of the nongravitational force in the radial, transverse, and normal directions, with respect to the comet's orbital plane. The latitude of the subsolar point, qJ, is calculated from sin 0 = sin I sin(fI + o0,
(14)
158
FANALE AND SALVAIL Using basic trigonometric identities for sin(/3 + I ) , cos fl, and the a b o v e expression for sin Ar, and simplifying with other basic identities yield the result
i
sin Ap --- sin(0 - 0c) cos I + cos(0 - 0c) sin ~bc sin I. V
(19)
The normal c o m p o n e n t of the force is
Ap tO Sort
Fn = F sin Ap = F[sin(to - 0~) cos I + cos(to - 0~) sin ~bc sin I].
(20)
An expression for the transverse component can be found in terms of trigonometry, but the derivation is more complicated than for the normal c o m p o n e n t . It is easier to FIG. 1. Geometry for calculation of nongravitational force components.
use
Ft = -+ N/F 2 - F~ - F~. where the obliquity, I, and the angle between the ascending node on the equator and the subsolar point at perihelion, 11, are specified, and the true anomaly, a, has been previously calculated for a given orbital position. Referring to Fig. 1, the angle b e t w e e n the line of action of F at 0~ and ~b~ and the radial direction is found using the spherical law of cosines as COS A r -- COS(t0 -
0c) c o s qSc.
(15)
Then the radial c o m p o n e n t of the force is Fr = F cos(to - 0c) cos 4~c.
(16)
It is convenient to consider next the normal c o m p o n e n t of force. Again referring to Fig. 1, the angle fl is e x p r e s s e d using the spherical law of sines as sin/3 =
sin(to - 0c) sin Ar sin(to - 0c) = k/1 - cos2(to - 0~) cos2~5c" (17)
An expression for A p , which is the angle b e t w e e n the line of action of F and the orbital plane, is found from the spherical law of sines as sin Ap = sin Ar sin(/3 + I).
(18)
(21)
In Eq. (21) the plus sign is used for prograde rotation or for obliquities less than 90 ° . The minus sign is used for retrograde rotation or for obliquities greater than 90 ° . The dust flux is a s s u m e d to have no effect on the nongravitational force since dust particles are ejected from the nucleus from zero velocity. Therefore, in the limit, mom e n t u m transfer takes place between dust particles and gas molecules after the particles are ejected and not with the nucleus. H o w e v e r , if dust particles are broken away, rather than being simply lifted, the energy of the gas flow would be decreased accordingly. This would cause the nongravitational force to be s o m e w h a t lower than given by this analysis. Theoretically determined nongravitational force c o m p o n e n t s , obtained from the a b o v e equations, can be used with empirical results for individual comets to estimate the sizes of c o m e t a r y nuclei. Marsden et al. (1973) used a semiempirical procedure to calculate nongravitational accelerations for several comets. In the equation of motion the magnitudes of the nongravitational acceleration c o m p o n e n t s are expressed as (Yeomans, 1977) ar = A I g ( R )
(22)
COMET MODEL
159
TABLEI
IV. R E S U L T S A N D D I S C U S S I O N
VARIABLE PARAMETERS SELECTED FOR RUNS
Results were obtained for various combinations of obliquity and thermal conductivity of the nucleus, as indicated in Table I. Values of important parameters are listed in Table II. Surface and subsurface temperatures, gas fluxes, and dust fluxes were calculated for a global grid of points for each 15° hour angle and 10° latitude and integrated to obtain gas and dust production rates. Nongravitational force components were obtained from the equations in the preceding section. These calculations were carried out over five complete orbits, starting and ending at aphelion, using Comet Halley's well-known orbital parameters. Calculations over five orbits allowed surface temperatures to converge to ~ 2 ° at
Model
Obliquity (deg)
Thermal conductivity
1
0
8 x 104
2 3 4
45 0 45
8 x 104 8 x lO3 8 x 103
and (23)
at = A 2 g ( R ) ,
fRj,j-, In the equation for g ( R ) , a , m , n, k , and R0 are constant parameters for a given comet. The parameters A1 and A2 for a given c o m e t are obtained from a leastsquares differential correction procedure. The theoretically determined nongravitational force components as derived earlier can be expressed in the form Fi
=
Ci re,
(25)
where Fi represents the ith nongravitational force component, and Ci is a parameter whose value is determined by the gas being sublimed and the array of surface temperatures. The corresponding acceleration is ai -
Fi m
-
Cir2 ]pn~rr 3
-
3Ci 4pnrrr"
(26)
On is the density of the cometary nucleus; Ci is obtained from Eq. (25). F~ is calculated independently from Eq. (16) if i = r or from Eqs. (16), (20), and (21) if i = t, and r is the assumed nuclear radius in the physical model. If a; is set equal to either ar or at in Eqs. (22) and (23), the actual radius of the cometary nucleus can be estimated from 3C,. r = - 4pnTrai
(27)
for i = r or i = t. This procedure will be used to estimate the radius o f Comet Halley in the next section.
T A B L E II VALUES OF IMPORTANT PARAMETERS a
Name Orbital and rotational parameters Semimajor axis Eccentricity Rotation period A r g u m e n t of subsolar meridian at perihelion Optical properties N u c l e u s albedo N u c l e u s emissivity Single scattering albedo of dust in the visible Single scattering albedo of d u s t in the infrared D u s t extinction c r o s s section in the visible Dust extinction cross section in the infrared Dust emissivity N u c l e u s properties Radius D u s t to ice ratio E x p o n e n t o f grain size distribution M a x i m u m grain size M i n i m u m grain size
Value
17.97 A U 0.9673 41.5 h b 340 °b 0.1 0.9 0.9' -0 C 6 × 10 -sc 2 x 10 -9c
1.0" 2.5 k m 0.5 -4.2 1.0 cm 0.5 x 10 4 cm
See Table I for obliquity and thermal conductivity. b Sekanina and L a r s o n (1984). c S q u y r e s et al. (1985).
160
FANALE AND SALVAIL
blc
810 9
8
CI 2
~7
--7
11
~s
~9
T~8
F, Fn
F,
$
$3
Z
g
2
J
6
I i
I
2
3
4
Heliocentric Dista~e (AU }
t
5
t
6
Heliocentric Distance (AU }
4t
i
1
2 3 4 5 Heliocentric Distance (AU)
FIG. 2. (a) Gas production rate, (b) dust production rate, and (c) nongravitational forces vs heliocentric distance for Model 2.
aphelion. Temperatures at smaller heliocentric distances were more accurate. Increased accuracy was not attempted due to long running times. The figures are numbered according to the model numbers shown in Table I. The results include the effects of attenuation, scattering, and thermal emission of radiation by dust in the coma. These effects cause the nucleus surface to be more isothermal. For example at 2 AU the maximum and minimum nucleus surface temperatures differ by about 100°K for the low-obliquity, high-thermal-conductivity case while at perihelion the difference is only about 25°K. The optical depth tends to lower slightly the dayside temperatures, while the radiation scattering and thermal emission by the dust significantly increase nightside temperatures. At perihelion the maximum surface temperature is lowered by about 2°K due to an optical depth of - 0 . 4 . H o w e v e r , the coma has only a small effect on the gas and dust production rates because the sublimation rates are governed essentially by the maximum temperatures due to the exponential dependence of sublimation rate on temperature. We will now focus, in turn, on the two
parameters chosen for investigation. For the obliquity we used 0 and 45 °. Two important characteristics distinguish the 0 and 45 ° obliquity plots. The first is that the 45 ° results show a small preperihelion asymmetry in the gas and dust production rates inside ~ 1 AU and an increasingly large postperihelion a s y m m e t r y beyond - 1 AU. The 0 ° results show only an increasing postperihelion asymmetry. The preperihelion asymmetry in the 45 ° results looks remarkably similar to the lightcurve produced by Yeomans (1977) from 1910 data except that the asymmetry is smaller. H o w e v e r , a consensus later developed among several observers that the dip in the lightcurve after perihelion was an artifact of the unusual observing conditions. This is discussed in detail by Newburn (1981). As a result the postperihelion portion of the curve inside 1 AU was modified to a consistent postperihelion a s y m m e t r y more like the 0 ° curves. The lightcurve produced by Yeomans (1977) and the postperihelion modification of it now being used for models at JPL are also shown in Newburn (1981). The second important characteristic is seen in the nongravitational force components. For the 45 ° cases Figs. 2c and 3c show that the
COMET MODEL
161
b,o
81C
9 GI 2
8
--7
I1
~7
~s
I
2 3 4 5 Heliocentric Distance(AU)
6
2
~6
I
5
O(
1 Heliocentric D~stor~e(AU)
2 3 4 5 Heliocentric Distonce(AU)
6
FIG. 3. (a) Gas productionrate, (b) dust productionrate, and (c) nongravitationalforces vs heliocentric distance for Model 4.
normal component of the nongravitational force is usually of the same order of magnitude as the transverse component. The normal components in both Figs. 2c and 3c appear to have discontinuities near perihelion. These are the positions where the normal component of the force reverses direction. The dotted lines represent the orbital positions where the normal component of the force is reversed, pointing south with respect to the orbital plane. The difference in the appearance of the two plots and, in particular, the range of heliocentric distances for which the normal component of the force is reversed are due only to the value used for the thermal conductivity. A larger thermal conductivity results in a larger angle between the subsolar point and the point where the resultant nongravitational force acts. This affects the angles 0c and 4~c in Eq. (20). The normal component of the force has the effect of changing the inclination of a comet's orbit. However, this effect would be small and is probably not able to be distinguished from the effects of planetary perturbations. Yeomans (1977) concluded that the acceleration component normal to the orbit plane has a negligible
effect on the orbital motion, a result that seems to be typical of other short period comets. This results from the fact that the normal component of the force has essentially no effect on the orbital velocity and, hence, the orbital period. For the 0 ° cases the normal component is always at least three orders of magnitude smaller than the transverse component. Ideally the normal component should be zero for 0° obliquity, but small errors in the computation of the temperatures result in finite but relatively insignificant normal component forces, which we have neglected in plotting. The asymmetry between the post- and preperihelion gas production rates, especially inside l AU, suggests a conclusion with regard to the thermal conductivity of the nucleus. Figures 4 and 5 show results of models that differ only in the value assumed for the thermal conductivity of the nucleus (see Table I). The results obtained using the higher thermal conductivity have a significantly higher asymmetry in gas production than those obtained using the lower thermal conductivity. Observations tend to support the results with the larger asymmetry in gas production. Fanale and Salvail
162
F A N A L E AND S A L V A I L
blo
alO
9 C12 F,
8 ~7
~7
11
Ft
$o d:
°t 1
I
i
2 3 4 5 HeLioce~ric Distar~e (AU)
6
I
I
6J
,
40
1
Heliocentric Distot~e (AU)
2 3 4 5 Heliocmtri(: Distance (AU)
,
6
FIG. 4. (a) Gas production rate, (b) dust production rate, and (c) nongravitational forces vs heliocentric distance for Model 1.
(1984) concluded that postperihelion asymmetries are largely caused by heat conduction into the nucleus. In this work we find that a higher thermal conductivity produces a larger postperihelion asymmetry. This fact has been discussed for some time, though here we do a more quantitative study. Previous physical models of come-
tary nuclei h a v e utilized thermal c o n d u c t i v ities ranging f r o m on the o r d e r o f 102ergs/ c m - s e c - ° K for u n c o m p a c t e d s n o w and dust to 105 e r g s / c m - s e c - ° K for crystalline w a t e r ice. By c o m p a r i n g the near perihelion a s y m m e t r i e s in gas p r o d u c t i o n for our m o d e l results with the near perihelion a s y m m e t r y in the o b s e r v e d light c u r v e
810
i
9 8
C
~-7
12 F, 11
F,
~g o
i
Z 2
5
1
o,
I
2 3 4 5 Heliocentric Distance (AU)
I
2
3
4
Heliocentric Distar~e (AU)
5
J 6
4
1
2 3 4 5 Heliocentric Distonce (AU)
FIG. 5. (a) Gas production rate, (b) dust production rate, and (c) nongravitational forces vs heliocentric distance for Model 3.
COMET MODEL 8
--
a- 6
--
A
O
I
5 0
1.0
I 2.0
Heliocentric Distance (Au)
FIG. 6. Gas production rate vs heliocentric distance for Newburn's (1981) model.
given by Yeomans (1977), we conclude that Comet Halley probably has a thermal conductivity near the higher end of the range given above. Therefore, Halley's nucleus is probably crystalline ice, at least near the surface. The results all show gas production rates lower than those of the Newburn model results shown in Fig. 6. Part of the discrepancy is due to the fact that Newburn arbitrarily assumed that 20% of the volatile c o m p o n e n t consists of volatiles other than water having an average molecular weight of 44 amu. We believe that a pure water model is just as realistic in view of the discussion given by Delsemme in Yeomans (1977), in which he concludes on the basis of statistics that "all short period comets lose all gases more volatile than water very fast, at least from the outer layers of the nucleus that are heated enough to participate in the sublimations." The remaining discrepancy could be the result of a number of assumptions, other than nuclear radius, that N e w b u r n (1981) made in his model. Finally, we will use the procedure outlined in the previous section to estimate the radius of Comet Halley. We will use the transverse force component, since Yeomans (1977) shows that the transverse component is determined far better than the radial component. In that paper A2 = 0.0159
163
x 10 -8 , ~ = 0.111262, r0 = 2.808 AU, m = 2.15, n = 5.093, and k -- 4.6142. We performed calculations for nine points in the orbit inside 2 AU before and after perihelion, for both cases having the high thermal conductivity. The values for g(R) were calculated from Eq. (24) for the appropriate value of R (AU). The values for at were obtained by first converting the parameter A2 to cgs units to get A2 = 3.1736 x 10 -7 cm/sec 2 and then multiplying by the corresponding value for g(R). The numbers for Ft were provided by our program, and they were divided by r 2 (r = 2.5 × 105 cm) to get the corresponding values for Ct. We must emphasize here that any realistic number could be used for the radius in the computer program, since Ct is independent of the assumed radius. The values for the computed nuclear radius were then calculated from Eq. (27), assuming that the bulk density of the nucleus, pn, is 1.0 g/cm 3. Using a dust to ice ratio of 0.5 and assuming that the ice and silicate are nonporous, the bulk nuclear density would be - 1 . 5 g/cm 3. A bulk density of 1.0 g/cm 3 is equivalent to assuming that both the ice and silicate have a porosity of - 0 . 6 5 . If an interested reader prefers using a different bulk density, the results for the radius can be easily converted. The computed values of the radius were integrated over time, and the result was divided by the elapsed time for the portion of the orbit chosen. This procedure yielded averaged values for the radius of 3.135 km for the 0 ° obliquity case and 1.6435 km for the 45 ° obliquity case. Our model and results have several shortcomings. The agreement between our model results and the facts indicated by observations was not the best. N o n e of the model variations satisfied all the observational constraints. Among the low-obliquity model variations, the low-thermal-conductivity case (Model 3) shows gas and dust production rates that are more consistent with the facts that Comet Halley was first seen in 1910 at - 3 . 4 AU (preperihelion) and last seen at - 5 . 2 AU (postperihelion)
164
FANALE AND SALVAIL
(Newburn, 1981). However, these model variations show the observed asymmetry near perihelion to be significantly smaller than indicated by observations. The highthermal-conductivity case (Model 1) produces a near perihelion asymmetry that agrees well with observations but has considerably lower gas and dust production rates beyond - 3 AU, especially preperihelion, than are indicated by observations. One possible explanation for the supposedly higher gas production rates beyond - 3 AU is that electrostatic forces produced by the solar wind cause small residual dust grains to be emitted from the nucleus in the absence of the usual gas production levels. This phenomenon is analyzed and discussed in detail by Flammer et al. (1985). This dust could then be interpreted by an observer as a higher gas production rate. The accuracy of the radius calculations depends not only upon the results of our model but upon Yeoman's semiempirical results. Although errors would arise due to the idealizations required for the theoretical analysis, we also believe that the semiempirical results are inaccurate. It is well known that comets generally have significant asymmetries in gas and dust production, and this is certainly true for Comet Halley. This strongly implies that nongravitational forces and accelerations must be asymmetric with respect to perihelion. Therefore, the dependence of the nongravitational acceleration on heliocentric distance must be different between the preperihelion and postperihelion legs of the orbit. This is not so for any of the semiempirical formulations currently available. Further refinements are obviously needed before accurate calculations can be made using this method. Although this work has provided results of questionable accuracy,
it has indicated the methodology by which physical and semiempirical models can be used to reveal characteristics of individual comets. REFERENCES DELSEMME, A. H., AND D. C. MILLER (1971). Physico-chemical phenomena in Comets - - III. Planet. Space Sei. 19, 1229-1257. FANALE, F. P., AND J. R. SALVAIL 0984). An idealized short period comet model: Surface insolation, H20 flux, dust flux and mantle evolution. Icarus 60, 476-511. FANALE, F. P., R. S. SAUNDERS, B. BANERDT, J. STEVENS, AND E. LANE (1983). Properties o f Sirnulated Comet Mantle, NASA TM-86246, pp. 56-58. FLAMMER, K. R., B. JACKSON, AND D. A. MENDIS
(1985). On the brightness variations of Comet Halley at large heliocentric distances. Preprint. MARSDEN, B. G. (1968). Comets and nongravitational forces. Astron. J. 73, 367-379. MARSDEN, B. G. (1969). Comets and nongravitational forces, I1. Astron. J. 74, 720-734. MARSDEN, B. G., Z. SCKANINA,AND D. K. YEOMANS (1973). Comets and nongravitational forces V. Astron. J. 78, 211-225.
NEWaURN, R. L. (1981). A semi-empirical photometric theory of cometary gas and dust production: Application to P/Halley's gas production rates. Comet Halley Dust and Gas Environment Workshop. ESA SP-174. NEWBURN, R. L., AND D. K. YEOMANS (1982). Halley's Comet. Annu. Rev. Earth Planet. Sci. 10. SEKAN1NA, Z., AND S. M. LARSON (1984). Coma Morphology and dust-emission pattern of periodic comet Halley. II. Nucleus spin vector and modeling of major dust features in 1910. Astron. J. 89, 1408-1447. SQUYRES, S. W., C. P. McKAY, AND R. T. REYNOLDS, (1985). Temperatures within comet nuclei.
Preprint. SUTTON, S. T., J. STEPHENS, R. S. SAUNDERS, AND F.
P. FANALE, (1983). Properties of residues from sublimed dirty ice: Implications for comets and Martian Polar deposits. Bull. Amer. Astron. Soc. 15, 847. WHIPPLE, F. L. (1980). Periodic Comet Halley. IAU Circ. 3459. YEOMANS, D. K. (1977). Comet Halley and nongravitational forces. In Comets, Asteroids" and Meteorites (A, R. Delsemme, Ed.), pp. 61-64. University of Toledo, Toledo.